# Inconsistencies in Unstructured Geometric Volume-of-Fluid Methods for Two-
Phase Flows with High Density Ratios

Jun Liu<EMAIL_ADDRESS>Tobias Tolle<EMAIL_ADDRESS>Davide Zuzio<EMAIL_ADDRESS>Jean-Luc Estivalezes  Jean-
<EMAIL_ADDRESS>Santiago Marquez Damian<EMAIL_ADDRESS>Tomislav Marić<EMAIL_ADDRESS>Mathematical Modeling and Analysis,
Technische Universität Darmstadt, Germany  Bosch Research, Robert Bosch GmbH;
research conducted at MMA, Technische Universität Darmstadt  The French
Aerospace Lab ONERA, France  Centro de Investigaciones en Mecánica
Computacional (CIMEC), UNL/CONICET, Argentina

###### Abstract

Geometric flux-based Volume-of-Fluid (VOF) methods [1] are widely considered
consistent in handling two-phase flows with high density ratios. However,
although the conservation of mass and momentum is consistent for two-phase
incompressible single-field Navier-Stokes equations without phase-change [2],
discretization may easily introduce small inconsistencies that result in very
large errors or catastrophic failure. We apply the consistency conditions
derived for the $\rho$LENT unstructured Level Set / Front Tracking method [2]
to flux-based geometric VOF methods [1], and implement our discretization into
the plicRDF-isoAdvector geometrical VOF method [3]. We find that computing the
mass flux by scaling the geometrically computed fluxed phase-specific volume
destroys the equivalence between the scaled volume fraction equation and the
mass conservation equation, depending on the choice for the temporal and
convective term discretization schemes. We propose two solutions. First, based
on the analysis of discretization errors, we suggest a consistent combination
of the temporal discretization scheme and the interpolation scheme for the
momentum convection term. Second, similar to [2], we solve an auxiliary mass
conservation equation with a geometrical calculation of the face-centered
density. We prove the equivalence between these two approaches mathematically
and verify and validate their numerical stability for density ratios of
[$1$,$10^{6}$] and viscosity ratios of [$10^{2}$,$10^{5}$].

###### keywords:

volume-of-fluid , unstructured , finite volume , high density ratios

[listing]style=Plaintop

##  1 Introduction

Numerical simulation methods still face challenges when dealing with
incompressible two-phase flows that involve fluid phases of significantly
different densities (high density ratios). We have provided a detailed review
of two-phase flow simulation methods that address the challenges in handling
high density ratios recently in [2], and here we only list more current
literature contributions.

Huang et al. [4] proposed the mixed Upwind/Central WENO scheme on a staggered
structured grid, an extension of the conventional WENO scheme [5] on
collocated grids, to spatially discretize the nonlinear convective term in
conservative form. The Upwind WENO scheme is used to evaluate velocity at cell
faces, while 3 different forms of Central WENO scheme are applied to evaluate
mass flux in x-/y-directions as well as density at cell corners. In addition
to the mentioned spatial discretization scheme, Huang et al. [4] proposed also
a semi-implict projection scheme to decouple pressure and velocity in momentum
transport equation. A backward difference scheme combined with a special
treatment for the viscous term is introduced to discretize the momentum
equation without the pressure gradient. The viscous term containing the
intermediate velocity is split into two parts, one of which comprises the
constant arithmetic mean of two phases’ viscosity and the intermediate
velocity, another one comprises the updated viscosity and the explicit
velocity (referring to [6, Apendix A]). Numerous 2D cases with density ratio
$[1,1000]$ are tested, whose results are in good agreement with analytical and
experimental results. The interface capturing method in [4] is selected to
phase-field. The authors indicated that the mixed Upwind/Central WENO scheme
is able to be coupled with any interface capturing/tracking method.

Xie et al. [7] introduced a consistent and balanced-force model with level-set
and volume of fluid function (CBLSVOF) on polyhedral unstructured grids. A
hybrid algebraic/geometric VOF method based on THINC/QQ [8] is developed to
capture interface, while the LS function constructed from volume fraction is
exploited to improve the accuracy of interface curvature evaluation. To
achieve the consistency, the convective term of both the volume-fraction
equation and momentum equation in the conservative form are discretized in the
same manner, wherein the velocity and volume fraction are reconstructed by
using the identical high-order reconstruction scheme based on a quadratic
polynomial function. The balanced-force formulation proposed firstly by
Francois et al. [9] is adapted in this work to eliminate the parasitic current
at the interface. The interface curvature estimated by the continuous signed-
distance function instead of the abruptly changing volume fraction function
offers computational simplicity and numerical stability. A 3D single bubble,
and two bubbles coalescence cases with the density ratio of magnitude $10^{3}$
are tested. The results agree well with previous experimental results.

Desmons and Coquerelle [10] proposed a generalized approach known as the HOMP
(High-Order Momentum Preserving) method. A high-order temporal and spatial
scheme, i.e., RK2 and WENO5,3 [11], are used to discretize an auxiliary
advection equation of characteristic function $\chi$ that is compatible with
the mass equation and independent of the underlying interface-representation
function. The momentum equation is then discretized using the same high-order
schemes as before, and the density is deduced from the characteristic function
$\chi$, which is computed from the auxiliary advection equation rather than
from the interface transport step. During the discretization, both the
auxiliary advection equation and the momentum equation remain in the
conservative form. The authors combined several interface representation
methods, i.e., Volume Of Fluid, Level Set Method, and Moment Of Fluid, with
HOMP and then tested them. The results from various selected validation cases,
including water and air, demonstrate good agreement with the literature
results. It is worth noting that no theoretical explanation is offered as to
why the auxiliary advection equation is required, rather than just maintaining
the consistency between the interface transport equation and the momentum
equation. Alternatively, the higher-order WENO scheme requires very large
stencils of variable width when used with the finite-volume method, making its
parallel implementation inefficient using the domain-decomposition / message-
passing parallel programming model.

In this manuscript, we apply the consistency requirements imposed by the
single-field Navier-Stokes equations derived in [2] to flux-based geometric
Volume-of-Fluid methods [1], using the isoAdvector method to verify and
validate our findings [3, 12, 13]. The consistency requirements from [2]
indicate that geometrical flux-based VOF methods are inherently consistent in
handling high density ratios. On the modeling level, an exact solution of the
volume fraction equation is equivalent to solving the mass conservation
equation. On the discrete level, flux-based VOF methods consistently compute
the mass flux needed in the implicitly discretized single-field momentum
convection term, by scaling the fluxed phase-specific volume. However, we show
that inconsistencies easily arise through the choice of an inconsistent
combination of the temporal integration scheme and the interpolation scheme
for the two-phase momentum conservation term, leading to significant errors
for small density ratios and catastrophic failures for large density ratios.

In the following sections, we briefly review the single-field formulation of
incompressible two-phase Navier-Stokes equations and analyze their
discretization using the collocated unstructured Finite Volume method [14]. In
section section 3.1, we show that a consistent combination of the temporal
discretization scheme and the interpolation scheme for the momentum convection
term is necessary for ensuring stable solutions with high density ratios. In
section sections 3.2 and 3.3, we introduce an auxiliary mass conservation
equation with a geometrical calculation of the face-centered density to the
geometrical VOF method isoAdvector[3, 12, 13]. We prove the equivalence
between these two approaches and verify and validate their numerical stability
for density ratios of [$1$,$10^{6}$] and viscosity ratios of
[$10^{2}$,$10^{5}$] in the results section section 4.

##  2 Mathemathical model

Figure 1: The domain $\Omega$, split by the fluid interface $\Sigma(t)$ into
two sub-domains $\Omega^{\pm}$.

Figure 1 illustrates an incompressible two-phase flow system without phase
change. The flow domain $\Omega\subset\mathbb{R}^{3}$ with a boundary
$\partial\Omega$ is partitioned into two subdomains $\Omega^{+}(t)$ and
$\Omega^{-}(t)$, occupied by two different fluid phases. The boundary between
$\Omega^{\pm}(t)$, referred to as the fluid interface, is denoted by
$\Sigma(t)$ and has the normal $\mathbf{n}_{\Sigma}$ pointing, say, outwards
of the domain $\Omega^{-}(t)$.

A phase indicator function

$\chi(\mathbf{x},t):=\begin{cases}1,&\mathbf{x}\in\Omega^{-}(t),\\\
0,&\mathbf{x}\in\Omega^{+}(t),\end{cases}$ (1)

indicates $\Omega^{\pm}(t)$ and formulates single-field density and dynamic
viscosity as

$\rho(\mathbf{x},t)=(\rho^{-}-\rho^{+})\chi(\mathbf{x},t)+\rho^{+},$ (2)
$\mu(\mathbf{x},t)=(\mu^{-}-\mu^{+})\chi(\mathbf{x},t)+\mu^{+},$ (3)

where $\rho^{\pm}$ and $\mu^{\pm}$ are the constant densities and dynamic
viscosities in $\Omega^{\pm}(t)$. The one-fluid density and dyamic viscosity
are used in the Navier-Stokes equations without phase change, namely

$\displaystyle\nabla\cdot\mathbf{v}$ $\displaystyle=0,$ (4)
$\displaystyle\partial_{t}(\rho\mathbf{v})+\nabla\cdot(\rho\mathbf{v}\otimes\mathbf{v})$
$\displaystyle=-\nabla
p-(\mathbf{g}\cdot\mathbf{x})\nabla\rho+\nabla\cdot(\mu\left(\nabla\mathbf{v}+(\nabla\mathbf{v})^{T}\right))+\mathbf{f}_{\Sigma}.$
(5)

where $p$ is the modified pressure defined as the subtraction of hydrostatic
pressure from the total pressure, i.e., $p=P-\rho\mathbf{g}\cdot\mathbf{x}$.
The surface tension force $\mathbf{f}_{\Sigma}$ is exerted only on the
interface, as

$\mathbf{f}_{\Sigma}=\sigma\kappa\mathbf{n}_{\Sigma}\delta_{\Sigma},$ (6)

where $\sigma$ is the constant surface tension coefficient, $\kappa$ is twice
the local mean curvature of the interface $\Sigma(t)$, and $\delta_{\Sigma}$
is the interface Dirac distribution.

This singe-field formulation of incompressible two-phase Navier-Stokes
Equations (NSE) without phase change is especially relevant for engineering
applications because it provides a solid modeling basis for two-phase flows
with fluid interfaces that can arbitrarily deform, breakup and merge, provided
that the method responsible for advecting the phase indicator does not impose
its own restrictions regarding fluid interface deformation and topological
changes. Single-field NSE also embed strict consistency requirements [2]: an
equivalence between the conservation of mass and the phase-indicator (volume)
conservation. These requirements translate on the discrete level into
requirements for the computation of the mass flux in the discretized two-phase
momentum convection term in eq. 5, analyzed in the following section.

##  3 Numerical methodology

###  3.1 Single-field mass conservation and volume fraction conservation

The discretization of the single-field two-phase momentum convection term from
eq. 5 is reformulated by

$\int_{\Omega_{c}}\nabla\cdot(\rho\mathbf{v}\otimes\mathbf{v})\,dV=\int_{\partial\Omega_{c}}(\rho\mathbf{v}\otimes\mathbf{v})\cdot\mathbf{n}\,ds=\sum_{f\in
F_{c}}\rho_{f}F_{f}\mathbf{v}_{f}+\mathbf{e}_{\rho\mathbf{v}}(h^{2}),$ (7)

for the domain $\Omega$ discretized by $|C|$ non-overlapping control volumes
$\Omega:=\cup_{c\in C}\Omega_{c}$. Each finite volume $\Omega_{c}$ is bounded
by a number ($|F_{c}|$) of non-planar surfaces $S_{f}$, i.e.
$\partial\Omega_{c}=\cup_{f\in F_{c}}S_{f}$, which is the surface-area-normal
vector oriented outwards of $\Omega_{c}$. Linearizing the convective term for
the velocity $\mathbf{v}$ and using the second-order accurate -
$\mathbf{e}_{\rho\mathbf{v}}(h^{2})$ is the momentum equation discretization
error \- collocated Unstructured Finite Volume Method (UFVM) [15, 14, 16]
introduces the face-centered volumetric flux

$F_{f}:=\mathbf{v}_{f}\cdot\mathbf{S}_{f},$ (8)

where $\mathbf{v}_{f}$ is the second-order accurate face-centered velocity
average. The collocated UFVM denotes with $f$ second-order accurate face-
centered area-averaged quantities, and with $c$ second-order accurate cell-
centered volume-averaged quantities.

To compute $\rho_{f}F_{f}$, we first consider the volume fraction equation.
The volume fraction is defined using the phase-indicator eq. 1 as

$\alpha_{c}(t):=\frac{1}{|\Omega_{c}|}\int_{\Omega_{c}}\chi(\mathbf{x},t)\,dV.$
(9)

The volume fraction advection equation [1]

$\frac{d}{dt}\int_{\Omega_{c}}\chi\,dV=|\Omega_{c}|\frac{d}{dt}\alpha_{c}(t)=-\int_{\partial\Omega_{c}}\chi\mathbf{v}\cdot\mathbf{n}\,dS,$
(10)

is equivalent to the advection equation for the phase $\Omega^{+}(t)$
indicated by $1-\chi(\mathbf{x},t)$, namely

$\frac{d}{dt}\int_{\Omega_{c}}(1-\chi)\,dV=-|\Omega_{c}|\frac{d}{dt}\alpha_{c}(t)=-\int_{\partial\Omega_{c}}(1-\chi)\mathbf{v}\cdot\mathbf{n}\,dS,$
(11)

since $\frac{d}{dt}\int_{\Omega_{c}}\,dV=0$, as $\Omega_{c}\neq\Omega_{c}(t)$,
and
$\int_{\partial\Omega_{c}}\mathbf{v}\cdot\mathbf{n}\,dS=\int_{\Omega_{c}}\nabla\cdot\mathbf{v}\,dV=0$
for incompressible two-phase flows without phase change.

The mass conservation with the single-field density eq. 2 in a fixed control
volume $\Omega_{c}$ states

$\frac{d}{dt}\int_{\Omega_{c}}\rho\,dV=-\int_{\partial\Omega_{c}}\rho\mathbf{v}\cdot\mathbf{n}\,dS.$
(12)

Inserting eq. 2 into the r.h.s. of eq. 12 and integrating over the time step
$[t^{n},t^{n+1}]$ results in

$\displaystyle=\rho_{c}^{n}+\frac{1}{|\Omega_{c}|}\left[-\rho^{-}\int_{t^{n}}^{t^{n+1}}\sum_{f\in
F_{c}}\int_{S_{f}}\chi\mathbf{v}\cdot\mathbf{n}\,dS\,dt-\rho^{+}\int_{t^{n}}^{t^{n+1}}\sum_{f\in
F_{c}}\int_{S_{f}}(1-\chi)\mathbf{v}\cdot\mathbf{n}\,dS\,dt,\right]$ (13)

Inserting the eqs. 10 and 11 in the eq. 13 to replace the sums of surface
integrals results in

$\rho_{c}^{n+1}=\rho_{c}^{n}+\frac{\rho^{-}}{|\Omega_{c}|}|\Omega_{c}|\int_{t^{n}}^{t^{n+1}}\frac{d}{dt}\alpha_{c}(t)\,dt+\frac{\rho^{+}}{|\Omega_{c}|}|\Omega_{c}|\int_{t^{n}}^{t^{n+1}}-\frac{d}{dt}\alpha_{c}(t)\,dt,$
(14)

leading finally to

$\rho_{c}^{n+1}=\rho_{c}^{n}+(\rho^{-}-\rho^{+})(\alpha_{c}^{n+1}-\alpha_{c}^{n}).$
(15)

Equation 15 shows that solving eq. 12 for $\rho_{c}^{n+1}$ is equivalent to
adding the solution of the volume fraction equation scaled with
$(\rho^{-}-\rho^{+})$ to $\rho_{c}^{n}$.

Alternatively, integrating the single-field density model (eq. 2) over the
control volume $\Omega_{c}$ gives

$\int_{\Omega_{c}}\rho\,dV=\int_{\Omega_{c}}\rho^{-}\chi+\rho^{+}(1-\chi)\,dV.$
(16)

Dividing eq. 16 by $|\Omega_{c}|$ using eq. 9 results in the discrete single-
field density model,

$\rho_{c}(t)=\rho^{-}\alpha_{c}(t)+\rho^{+}(1-\alpha_{c}(t)),$ (17)

which, evaluated at $t^{n+1}$ and $t^{n}$ and subsequently subtracted, results
in

$\displaystyle\rho_{c}^{n+1}-\rho_{c}^{n}$
$\displaystyle=\rho^{-}\alpha_{c}^{n+1}+\rho^{+}(1-\alpha_{c}^{n+1})-\rho^{-}\alpha_{c}^{n}-\rho^{+}(1-\alpha_{c}^{n})$
$\displaystyle=(\rho^{-}-\rho^{+})(\alpha_{c}^{n+1}-\alpha_{c}^{n}),$

which is eq. 15: solving the volume fraction advection equation eq. 10 and
using the single-field density model to compute the cell centered density by
eq. 2 is equivalent to solving the mass conservation equation (eq. 12).

Note that all equations eq. 9 \- eq. 17 are exact, as
$\partial\Omega_{c}:=\bigcup_{f\in F_{c}}S_{f}$, where $S_{f}$ are non-linear
surfaces that bound the control volume $\Omega_{c}$, and reformulating the
integration in time is exact by the fundamental theorem of calculus.

Ghods and Herrmann [17] introduce the auxiliary mass conservation equation as
a means for ensuring the consistency of the two-phase momentum convection, and
in [2] we provide the theoretical reasoning for the auxiliary mass
conservation equation. Contrary to Ghods and Herrmann [17], eq. 15 and eq. 16
both demonstrate, that the solution of the volume fraction equation in the
context of the VOF method [1] is exactly equivalent to the solution of the
mass conservation equation, rendering an auxiliary mass conservation equation
unnecessary.

However, the following question arises: if flux-based algebraic/geometric VOF
methods inherently ensure numerical stability for the two-phase momentum
convection with high density ratios, where do the numerical inconsistencies
reported throughout the literature [18, 19, 20, 21, 22, 23] come from?

Although eqs. 15 and 16 show the inherent consistency of VOF methods in the
mathematical model, the discrete computation of $\alpha_{c}^{n+1}$, the
approximation of the mass flux $\rho_{f}F_{f}$, the choice of the temporal
scheme and the flux limiting scheme, can potentially cause inconsistencies.

We begin with the discrete computation of $\alpha_{c}^{n+1}$ and the
approximation of $\rho_{f}F_{f}$ by scaling the fluxed phase-specific volume.
For a second-order accurate flux-based VOF method, the approximations applied
to the temporal and surface integrals for the fluxed phase-specific volumes
$V_{f}^{\alpha}$ when solving eq. 10 lead to

$\alpha_{c}^{n+1}=\alpha_{c}^{n}-\frac{1}{|\Omega_{c}|}\sum_{f\in
F_{c}}\int_{t^{n}}^{t^{n+1}}\int_{S_{f}}\chi\mathbf{v}\cdot\mathbf{n}\,dS\,dt=\alpha_{c}^{n}-\frac{1}{|\Omega_{c}|}\sum_{f\in
F_{c}}|V_{f}^{\alpha}|_{s}+e_{\alpha_{t}}(\Delta
t^{p})+e_{\alpha_{h}}(h^{2}),$ (18)

with $e_{\alpha_{t}}(\Delta t^{p})$ and $e_{\alpha_{h}}(h^{2})$ as temporal
and spatial volume fraction discretization errors. Note that

$|V_{f}^{\alpha}|_{s}:=sgn(F_{f})|V_{f}^{\alpha}|$ (19)

is a signed magnitude of a phase-specific volume $V_{f}^{\alpha}$ [1], whose
sign is determined by the volumetric flux $\rho_{f}F_{f}$.

The phase-specific volume $|V_{f}^{\alpha}|_{s}$ can be used to approximate
the mass flux $\rho_{f}F_{f}$ using eq. 15 and eq. 18. Reordering eq. 18,
results in

$\alpha_{c}^{n+1}-\alpha_{c}^{n}=\frac{1}{|\Omega_{c}|}\sum_{f\in
F_{c}}|V_{f}^{\alpha}|_{s}+e_{\alpha_{t}}(\Delta
t^{p})+e_{(}\alpha_{h})(h^{2})$ (20)

Inserting eq. 20 into eq. 15 results in

$\rho_{c}^{n+1}-\rho_{c}^{n}=\frac{\rho^{-}-\rho^{+}}{|\Omega_{c}|}\sum_{f\in
F_{c}}|V_{f}^{\alpha}|_{s}+(\rho^{-}-\rho^{+})[e_{\alpha_{t}}(\Delta
t^{p})+e_{\alpha_{h}}(h^{2})].$ (21)

Equivalently integrating eq. 12 over $[t^{n},t^{n+1}]$ results in

$\rho_{c}^{n+1}-\rho_{c}^{n}=\frac{1}{|\Omega_{c}|}\sum_{f\in
F_{c}}\int_{t^{n}}^{t^{n+1}}\rho_{f}F_{f}\,dt=\frac{1}{|\Omega_{c}|}\sum_{f\in
F_{c}}|M_{f}|_{s}+e_{\rho_{t}}(\Delta t^{s})+e_{\rho_{h}}(h^{2}),$ (22)

with $e_{\rho_{t}}(t^{s}),e_{\rho_{h}}(h^{2})$ as the temporal and spatial
discretization errors of the mass conservation equation. The right-hand sides
of eqs. 21 and 22 express the mass fluxed through $\partial\Omega_{c}$ over
$[t^{n},t^{n+1}]$, i.e. $|M_{f}|_{s}$, as the phase-specific volume
$|V_{f}^{\alpha}|_{s}$ scaled by the density difference, thus connecting the
fluxed mass with the fluxed phase specific volume.

Consistency of the volume fraction advection and the mass conservation on the
discrete level requires the equivalence of eqs. 21 and 22. Equations 21 and 22
will be exactly the same, only if their errors on the r.h.s. cancel out. Error
cancellation is impossible if equations eqs. 21 and 22 are using different
numerical schemes for $|V_{f}^{\alpha}|_{s}$ and $|M_{f}|_{s}$, which strongly
connects the fluxed mass $|M_{f}|_{s}$ with the fluxed phase-specific volume.
In other words, if we use a specific VOF method for $|V_{f}^{\alpha}|_{s}$, we
should compute the fluxed mass from $|V_{f}^{\alpha}|_{s}$. This is
hypothetical, of course, since there is no need to actually solve two
equations that are equivalent. However, we do need the mass flux for the
discretized momentum equation (eq. 5). If the fluxed mass must be computed
from the scaled $|V_{f}^{\alpha}|_{s}$ to ensure the consistency of the volume
fraction advection and mass conservation, it follows that the mass flux
consistent with mass conservation must also be computed from
$|V_{f}^{\alpha}|_{s}$, using eq. 21.

However, geometrical VOF methods use temporal integration schemes and
geometrical algorithms to compute $|V_{f}^{\alpha}|_{s}$ \- a time-integrated
quantity - that are very unlike the algebraic discretization schemes applied
on eq. 5. Precisely this difference is a source of very strong instabilities
and catastrophic failures for high density ratios.

It helps to consider a concrete example. Consistency of the mass flux can be
shown when the mass flux is computed by scaling the fluxed phase-specific
volume in eq. 21 with $(\rho^{-}-\rho^{+})$, in a simplified first-order
geometrical VOF method, which uses the ”Euler” temporal integration (rectangle
quadrature) of $|V_{f}^{\alpha}|_{s}$, i.e.

$\displaystyle|V_{f}^{\alpha}|_{s}^{Euler}$
$\displaystyle=\int_{t^{n}}^{t^{n+1}}\int_{S_{f}}\chi\mathbf{v}\cdot\mathbf{n}\,dS\,dt=\frac{F_{f}^{n}}{|S_{f}|}\int_{t^{n}}^{t^{n+1}}\int_{S_{f}}\chi\,dS\,dt+e_{\alpha_{t}}(\Delta
t^{2})+e_{\alpha_{h}}(h^{2})$ (23)
$\displaystyle=F_{f}^{n}\int_{t^{n}}^{t^{n+1}}\alpha_{f}(t)\,dt+e_{\alpha_{t}}(\Delta
t^{2})+e_{\alpha_{h}}(h^{2})$ $\displaystyle=F_{f}^{n}\alpha_{f}^{n}\Delta
t+e_{\alpha_{t}}(\Delta t^{2})+e_{\alpha_{h}}(h^{2}),$

where we define the fraction of the wetted face area $A_{f}(t)$ as
$\alpha_{f}:=A_{f}(t)/S_{f}$, and $N_{t}e_{\alpha_{t}}(\Delta
t^{2})=\frac{T}{\Delta t}e_{\alpha_{t}}(\Delta t^{2})\propto\Delta t$ results
in a first-order temporal quadrature error over the simulated physical time
$T$. We emphasize that eq. 23 is a simplified scheme used here only to discuss
the mass-flux consistency condition on the discrete level, it is not a
practically usable scheme for advecting $\alpha_{c}$, because it is
significantly less accurate than modern geometrical schemes [1]. Multiplying
$|V_{f}^{\alpha}|_{s}^{Euler}$ from eq. 23 with $(\rho^{-}-\rho^{+})$ to
obtain $|M_{f}|_{s}$ ensures the consistency of eqs. 21 and 22, namely

$\displaystyle\frac{(\rho^{-}-\rho^{+})|V_{f}^{\alpha}|_{s}^{Euler}}{\Delta
t}$
$\displaystyle=(\rho^{-}-\rho^{+})F_{f}^{n}\alpha_{f}^{n}+(\rho^{-}-\rho^{+})(e_{\alpha_{t}}(\Delta
t)+\frac{1}{\Delta t}e_{\alpha_{h}}(h^{2})),$ (24)
$\displaystyle=(\rho_{f}^{n}F_{f}^{n})^{Euler}+(\rho^{-}-\rho^{+})(e_{\alpha_{t}}(\Delta
t)+\frac{1}{\Delta t}e_{\alpha_{h}}(h^{2})).$

Note that the temporal accuracy lost by dividing by $\Delta t$ is recovered
when the mass flux is integrated over (multiplied with) $\Delta t$, from the
temporal term in the momentum conservation equation. The next step, therefore,
uses the mass flux $(\rho_{f}^{n}F_{f}^{n})^{Euler}$ from eq. 24 in the
discretized convective term from eq. 5. Since the unstructured Finite Volume
Method linearizes the volumetric flux when discretizing the convective term,
an Euler implicit discretization of eq. 5 leads to the contribution from the
discretized momentum convection term in the form or

$\frac{\Delta
t}{|\Omega_{c}|}\left(\int_{\Omega_{c}}\nabla\cdot(\rho\mathbf{v}\otimes\mathbf{v})\,dV\right)_{t^{n+1}}\approx\frac{\Delta
t}{|\Omega_{c}|}\sum_{f\in F_{c}}\rho_{f}^{o}F_{f}^{o}\mathbf{v_{f}}^{n+1},$
(25)

with the factor $\frac{\Delta t}{|\Omega_{c}|}$ resulting from the finite-
difference approximation of the cell-centered average of the temporal
derivative term, i.e. $(\partial_{t}\rho\mathbf{v})_{c}$, recovering the
first-order temporal accuracy in eq. 24 for the implicit Euler temporal
discretization scheme physical time $T$. Note that eq. 5 is solved
(discretized) iteratively in a segregated solution algorithm (e.g., [24, 2]),
so the linearized mass flux is denoted with the outer iteration index $1\leq
o\leq N_{o}$. If the solution algorithm converges,
$\rho_{f}^{o}F_{f}^{o}=\rho_{f}^{n+1}F_{f}^{n+1}$. The Euler implicit scheme
in eq. 25 allows $\rho_{f}^{o}F_{f}^{o}$ to be approximated with the
consistent mass flux $(\rho_{f}^{o}F_{f}^{o})^{Euler}$ from eq. 24, which
ensures the consistency between mass conservation and volume fraction
advection. Namely, in a converged solution, the Euler implicit discretization
requires $\rho_{f}^{o}F_{f}^{o}=\rho_{f}^{n+1}F_{f}^{n+1}$, and we can only
obtain a single mass flux $(\rho_{f}F_{f})^{Euler}$ by scaling
$|V_{f}^{\alpha}|^{Euler}$ in eq. 23, a mass flux averaged over $\Delta t$.
Therefore, increasing $N_{o}$ brings the mass flux in eq. 25 closer to the one
expected by the implicit Euler scheme; however, it does not impact
consistency.

The consistency of the mass conservation and volume fraction advection is
crucial in the momentum conservation equation, as any deviation from mass
conservation in momentum conservation increases the the source-term
$\sum_{f\in
F_{c}}\left(\frac{1}{a_{c}}\right)_{f}\left[\mathbf{H}(F^{o},\rho^{o},\mathbf{v}^{i-1})\right]_{f}\cdot\mathbf{S}_{f}$
in the pressure Poisson equation, with $i$ denoting the inner iterations of
the pressure equation (cf. [24, 2] for details). In other words, errors in the
mass flux artificially accelerate or decelerate the fluid as the pressure
equation tries to ensure volume (mass) conservation. Examining the errors in
eq. 24, we see that an error in the mass flux scaled from
$|V_{f}^{\alpha}|_{s}$ will be multiplied by $(\rho^{-}-\rho^{+})$: small
errors in the volumetric flux (fluxed phase-specific volume) are scaled with
the density difference, and this leads to large errors in the velocity field
and catastrophic failures for large density differences.

We further consider an alternative temporal integration scheme for eq. 5, say,
Crank-Nicolson scheme, resulting in contributions from the convective term in
the form of

$0.5\frac{\Delta t}{|\Omega_{c}|}\left(\sum_{f\in
F_{c}}\rho_{f}^{n}F_{f}^{n}\mathbf{v_{f}}^{n}+\sum_{f\in
F_{c}}\rho_{f}^{o}F_{f}^{o}\mathbf{v_{f}}^{n+1}\right)$ (26)

In this case, it is impossible to ensure consistency. The Crank-Nicolson
scheme strictly requires $\rho_{f}^{o}F_{f}^{o}=\rho_{f}^{n+1}F_{f}^{n+1}$,
and we obtain $(\rho_{f}F_{f})^{Euler}$ from eq. 24 as a single average
quantity over $\Delta t$, so we do not have two mass fluxes for eq. 26, making
eq. 24 inconsistent with the Crank-Nicolson scheme eq. 26 already in the first
time step.

The main takeway point is that computing the mass flux by scaling the fluxed
phase-specific volume over $\Delta t$ limits the temporal discretization to
schemes that utilize a single mass flux term within $\Delta t$ \- e.g., Euler
explicit or implicit, or 2nd-order backward implicit schemes. Any other
temporal discretization scheme (e.g., Runge-Kutta) that utilises a linear
combination of different mass fluxes within $\Delta t$ are inconsistent with
the mass flux scaling given by eq. 24.

Additionally, any modification of the scaled mass flux causes inconsistencies.
Concretely, limiting the mass flux in the discretized eq. 5 causes a hidden
inconsistency between mass and volume fraction conservation. The inconsistency
is hidden because we are not solving eq. 22, we are however using a scaled
flux, e.g., in eq. 24, resulting in the integration of a mass flux (along with
momentum) over $\Delta t$ that is different from the consistent mass flux from
eq. 24.

Figure 2: Geometric upwinding for $|V_{f}^{\alpha}|^{isoAdvector}_{s}$ in
plicRDF-isoAdvector [12].

We now turn our attention from the simplified VOF scheme eq. 23 to a more
complex, geometric isoAdvector scheme [3], that computes

$\displaystyle|V_{f}^{\alpha}|_{s}^{isoAdvector}$
$\displaystyle=\frac{0.5(F_{f}^{n}+F_{f}^{o})}{|\mathbf{S}_{f}|}\int_{t^{n}}^{t^{o}}\int_{S_{f}}\chi(\mathbf{x},t)\,dS\,dt+e_{\alpha_{t}}(\Delta
t^{2})+e_{\alpha_{h}}(h^{2})$ (27)
$\displaystyle=\frac{0.5(F_{f}^{n}+F_{f}^{o})}{|\mathbf{S}_{f}|}\int_{t^{n}}^{t^{o}}A_{f}(t)\,dt+e_{\alpha_{t}}(\Delta
t^{2})+e_{\alpha_{h}}(h^{2})$
$\displaystyle=0.5(F_{f}^{n}+F_{f}^{o})\int_{t^{n}}^{t^{o}}\alpha_{f}(t)\,dt+e_{\alpha_{t}}(\Delta
t^{2})+e_{\alpha_{h}}(h^{2})$

with $\int_{t^{n}}^{t^{o}}\alpha_{f}(t)dt$ computed geometrically by
displacing the piecewise-linear interface from the upwind cell, using the
upwind-cell velocity and the first-order accurate displacement approximation
$\mathbf{v}_{c}\Delta t$, schematically shown in fig. 2. The Euler temporal
integration of the displacement $\mathbf{v}_{c}\Delta t$ for the evaluation of
$\int_{t^{n}}^{t^{o}}\alpha_{f}(t)dt$ on the r.h.s. of eq. 27 makes it
possible to consistently compute the mass flux by scaling the phase-specific
fluxed volume

$(\rho_{f}^{o}F_{f}^{o})^{isoAdvector}=\dfrac{(\rho^{-}-\rho^{+})|V_{f}^{\alpha}|^{isoAdvector}_{s}}{\Delta
t},$ (28)

equivalently to eq. 24. Even though $\int_{t^{n}}^{t^{o}}\alpha_{f}(t)dt$ is
evaluated using exact geometric integration (cf. [12] for details), the
$A_{f}^{o}$, used as the end-point of the geometric integration is
approximated by first-order displacement $\mathbf{v}\Delta t$, that makes the
average displacement velocity of the fluid interface constant over $\Delta t$,
and, consequently, the mass flux from eq. 28 consistent.

The same consistency would be ensured if $V_{f}^{\alpha}$ would be constructed
using a flux-based VOF method that maps $S_{f}$ using the reverse flow-map
(e.g., [25]), if the approximation of the flow-map is first-order accurate. In
other words, if $V_{f}^{\alpha}$ is constructed geometrically using
displacements given by $\mathbf{v}(\mathbf{x},t^{n})\Delta_{t}$, regardless of
the geometrical approximation, $V_{f}^{\alpha}$ is constructed using
displacements constant over $\Delta t$, so dividing the volume
$V_{f}^{\alpha}$ with $\Delta t$ results in the consistent volumetric flux,
constant over $\Delta t$.

Trying to compute the consistent mass flux by scaling the phase-specific
fluxed volume for the Crank-Nicolson scheme,

$\sum_{f\in F_{c}}0.5\Delta
t(\rho_{f}^{o}F_{f}^{o}+\rho_{f}^{n}F_{f}^{n})=(\rho^{-}-\rho^{+})\sum_{f\in
F_{c}}|V_{f}^{\alpha}|_{s},$ (29)

involves $\rho_{f}F_{f}$ at $t^{n}$ and $t^{o}$ in the same equation, making
it impossible to express either one without the other, using only the known
phase-specific volume $|V_{f}^{\alpha}|$, fluxed over $\Delta t$ before the
solution of the momentum equation. This confirms the difficulty in evaluating
mass flux from time-integrated volumetric fluxes pointed out by Arrufat et al.
[21]. Furthermore, in order to integrate $|V_{f}^{\alpha}|$ using the Crank-
Nicolson scheme (i.e. trapezoidal quadrature), the position of the interface
at $t^{o}$ must be known, namely $\Sigma^{o}$. The interface $\Sigma^{o}$ can
be reconstructed from $\alpha_{c}^{o}$, which is advected using the forward
extrapolation in time (eqs. 27 and 2, [3]). However, this kind of quadrature
does not exactly correspond to the Crank-Nicolson scheme used in the momentum
equation, because the advection of the interface $\Sigma^{o}$ uses old-time
velocities $\mathbf{v}_{c}^{n}$, and Crank-Nicolson in the momentum equation
uses the new, implicit velocity $\mathbf{v}_{c}^{n+1}$.

If, however, the mass flux is not computed by scaling $|V_{f}^{\alpha}|$,
$\rho_{f}F_{f}$ is uniquely determined by $\chi(\mathbf{x},t)$ [2]. Namely, at
any time $t$ omitted here for brevity, from eq. 2, we get

$\int_{S_{f}}\rho\mathbf{v}\cdot\mathbf{n}\,dS=:(\rho_{f}F_{f})^{\rho}\approx\frac{F_{f}}{|S_{f}|}\int_{S_{f}}[(\rho^{-}-\rho^{+})\chi(\mathbf{x})+\rho^{+}]\,dS=(\rho^{-}-\rho^{+})\alpha_{f}F_{f}+\rho^{+}F_{f}$
(30)

with superscript $\rho$ in $(\rho_{f}F_{f})^{\rho}$ denoting the mass flux
estimated directly from the model for $\rho$ eq. 2, i.e. from $\rho^{\pm}$ and
the phase indicator $\chi$ (eq. 1). The phase-specific volumetric flux
$\alpha_{f}F_{f}$ from eq. 30 is, however, inconsistent with the equivalent
volumetric flux used for advecting volume fractions eq. 18, which can be
discretized with a wide range of different flux-based geometric VOF methods.
Therefore, if the mass flux $(\rho_{f}F_{f})^{\rho}$ from eq. 30 is used in
the discretized momentum eq. 5, cell-centered density in the discrete momentum
equation $\rho_{c}$ cannot be evaluated from $\alpha_{c}$, because of the
phase-specific volumetric flux inconsistency. Instead, the auxiliary mass
conservation equation introduced in [17] must be used to advect the cell-
centered density, that is only used for the $p-\mathbf{v}$ coupling, and
discarded at the end of a time step.

###  3.2 Collocated segregated solution algorithm with the auxiliary density
equation

We base our consistent solution algorithms for the unstructured Volume-of-
Fluid methods on the plicRDF-isoAdvector method [12] and implement the
auxiliary density equation solution into  the segregated algorithm (i.e.
solver) ”interIsoRhoFoam”, which is summarized by algorithm 1.

In [2], we show on the level of the mathematical model why solving an
auxiliary mass conservation equation plays a key role in reducing numerical
inconsistency caused by high density ratios.

We solve the mass conservation equation in the outer loop of the segregated
algorithm algorithm 1 in the following discrete form

$\rho^{o+1}_{c}=\rho^{o}_{c}+\frac{\Delta
t}{|V_{\Omega_{c}}|}\sum_{f}\rho^{o}_{f}F^{o}_{f}.$ (31)

Equation 31 is the auxiliary mass conservation (density equation). It is
solved after updating the volume fraction in the outer loop to
$\alpha_{c}^{o}$ by utilizing eq. 18 with any flux-based VOF method. Interface
reconstruction computes $\tilde{\chi}_{c}(\mathbf{x},t^{o})$ in every finite
volume $\Omega_{c}$ intersected by the fluid interface. The piecewise-linear
interface approximation $\tilde{\chi}_{c}(\mathbf{x},t^{o})$, together with
eq. 2, provides $\rho_{f}^{o}$ for eq. 31. Since $F_{f}^{o}$ is also
available, eq. 31 can be explicitly evaluated. To compute the mass flux, we
utilize the consistency relationship between the phase indicator and the face-
centered density $\rho_{f}$, derived in [2] for the Level Set / Front Tracking
method. It can be applied on any two-phase flow simulation method that is
using a phase indicator $\chi$. The surface integral of mass flux at time step
$o$ is expressed as

$\int_{\partial\Omega_{c}}\rho^{o}\mathbf{v}^{o}\cdot\mathbf{n}dS=\int_{\partial\Omega_{c}}[\rho^{-}\chi^{o}+\rho^{+}(1-\chi^{o})]\mathbf{v}^{o}\cdot\mathbf{n}\,dS,$
(32)

and discretized as

$\displaystyle\sum_{f\in F_{c}}\rho_{f}^{o}F_{f}^{o}$
$\displaystyle=\sum_{f\in
F_{c}}\left[\rho^{-}\alpha_{f}^{o}+\rho^{+}(1-\alpha_{f}^{o})\right]F_{f}^{o},$
(33)

where

$\alpha_{f}^{o}=\frac{1}{|S_{f}|}\int_{S_{f}}\chi(\mathbf{x},t^{o})\,dS=\frac{|A_{f}(t^{o})|}{|S_{f}|}.$
(34)

is used to define the face-centered density

$\rho_{f}^{o}=\rho^{-}\alpha_{f}^{o}+\rho^{+}(1-\alpha_{f}^{o})$ (35)

in the mass flux. The momentum equation is discretized as

$\rho^{o}_{c}\mathbf{v}^{o}-\rho^{n}_{c}\mathbf{v}^{n}+\frac{\Delta
t}{|\Omega_{c}|}\sum_{f}\rho^{o}_{f}F^{o}_{f}\mathbf{v}_{f}^{o}=\frac{\Delta
t}{|\Omega_{c}|}\mathbf{M},$ (36)

and solved for $\mathbf{v}_{f}^{o}$ in the outer loop of the $p-\mathbf{v}$
solution algorithm. The source term $\mathbf{M}$ is a shorthand term that
contains all the remaining contributions from the discretizaiton, used here
for brevity. The volume fractions $\alpha_{c}^{o}$ are already solved for
using eq. 18 at the beginning of the ”$o$” outer iteration, and used to
approximate the phase indicator $\tilde{\chi}_{c}(\mathbf{x},t^{o})$ for the
mass flux $\rho_{f}^{o}F_{f}^{o}$ using eq. 33, used in the same way in the
auxiliary density equation eq. 31, solved for $\rho_{c}^{o}$, as
$\rho_{f}^{o}F_{f}^{o}$ in eq. 36, without using flux limiters. The values
computed from the last outer loop are regarded as the new values at the time
step. At last, we need to restore the density field with respect to the volume
fraction to maintain consistency between them. The solution algorithm
interIsoRhoFoam is summarized by algorithm 1.

Algorithm 1 The solution algorithm interIsoRhoFoam.

1:while $t\leq t_{end}$ do

2: $t^{n+1}=t^{n}+\Delta t$

3: for $o=0$; $o<N_{outer}$; ++$o$ do

4: Solve volume fraction for $\alpha_{c}^{o}$ $\triangleright$ Equation 18

5: Reconstruct the phase indicator $\tilde{\chi}_{c}(\mathbf{x},t^{o})$

6: Compute $\rho_{f}^{o}$ from $\tilde{\chi}_{f}(\mathbf{x},t^{o})$
$\triangleright$ Equation 35

7: Compute the mass flux
$\rho_{f}^{o}F_{f}^{o}:=\rho_{f}^{o}(\mathbf{v}_{f}^{o-1}\cdot\mathbf{S}_{f})$

8: Solve the density for $\rho_{c}^{o}$ with $\rho_{f}^{o}F_{f}^{o}$
$\triangleright$ Equation 31

9: Discretize momentum equation (eq. 5) with $\rho_{c}^{o}$ and
$\rho_{f}^{o}F_{f}^{o}$.

10: for $i=0$; $i<N_{inner}$; ++$i$ do

11: Solve the pressure equation for $p_{c}^{i}$. $\triangleright$ Cf. [24, 2].

12: Compute $F_{f}^{i}$ and $\mathbf{v}_{c}^{i}$ from $p_{c}^{i}$.
$\triangleright$ Cf. [24, 2].

13: end for

14: end for

15: Restore $\rho_{c}^{n+1}$ consistent with $\alpha_{c}^{n+1}$, i.e.
$\rho_{c}^{n+1}=(\rho^{-}-\rho^{+})\alpha_{c}^{n+1}+\rho^{+}$ $\triangleright$
Equation 2.

16:end while

###  3.3 Phase-specific face area calculation

$\tilde{\Omega}_{D}^{+}(t^{o})$
---
$\tilde{\Omega}_{D}^{-}(t^{o})$
---
$\tilde{\Sigma}^{n+1}(t^{o})$
---
$\tilde{\Omega}_{U}^{-}(t^{o})$
---
$\tilde{\Omega}_{U}^{+}(t^{o})$
---
(a) Upwind and downwind cells intersected with piecewise-linear VOF
interfaces.

$(\alpha_{f})_{U}(t^{o})$
---
$(\alpha_{f})_{D}(t^{o})$
---
(b) A face $S_{f}$ intersected with downwind (D) and upwind (U) piecewise-
linear VOF interfaces.

Figure 3: Interface reconstructed as $\tilde{\Sigma}^{o}_{U,D}$ in upwind
($U$) and downwind ($D$) cells. Green polygons are interface polygons
$\tilde{\Sigma}^{o}_{U,D}\cap\Omega_{U,D}$. Blue lines are intersection line
segments $\tilde{\Sigma}^{o}_{U,D}\cap S_{f}$. Red points are intersection
points $\tilde{\Sigma}^{o}_{U,D}\cap\partial S_{f}$.

The integral eq. 34 leaves room for alternative discretizations, and therefore
requires attention. Flux-based geometrical Volume-of-Fluid methods advect
volume fractions using geometrical upwind advection schemes, transporting
$\tilde{\Omega}^{-}(t)\approx\Omega^{-}(t)$ geometrically from the upwind
cell, denoted with $U$ in fig. 3(a), to the downwind cell, denoted with $D$ in
fig. 3(a). Flux-based geometrical VOF methods therefore already provide the
intersection points and the intersection line segments denoted as red dots and
the blue line segments in fig. 3, are available, simplifying the area fraction
calculation in eq. 34. Since the geometrical VOF methods approximate the phase
indicator as piece-wise continuous, with a jump discontinuity not only across
the fluid interface $\Sigma(t)$ but also across the finite volume boundary
$\partial\Omega_{c}$, the line segments forming the intersection of the piece-
wise continuous interface and a face $S_{f}$ it intersects (see fig. 3(a)), do
not overlap. Figure 3(b) shows schematically the two intersection line
segments on a cell face. The intersection line segment of the interface and
the upwind cell are selected as the submerged area border to calculate the
area fraction. The evaluation of $\alpha_{f}$ and, from it,
$\rho_{f},\mu_{f}$, is summarized by algorithm 2.

Algorithm 2 Sub-algorithm for calculating $\alpha_{f},\rho_{f},\mu_{f}$.

1:Initialize $\rho_{f}$, $\alpha_{f}$ with upwind cell-centered values
$\rho_{U}$, $\alpha_{U}$.

2:$\alpha_{f}=\frac{1}{|S_{f}|}\int_{S_{f}}\tilde{\chi}_{U}(\mathbf{x})\,dS=\frac{|\tilde{\Omega}^{-}_{U}\cap
S_{f}|}{|S_{f}|}$ $\triangleright$ Equation 34 with upwind $\chi$.

3:$\rho_{f}=(\rho^{-}-\rho^{+})\alpha_{f}+\rho^{+}$. $\triangleright$ Equation
2.

4:$\mu_{f}=(\rho^{-}\nu^{-}-\rho^{+}\nu^{+})\alpha_{f}+\rho^{+}\nu^{+}$.
$\triangleright$ Equation 3.

It is important to note that $\mu_{f}$ is geometrically evaluated from a
geometrical $\alpha_{f}$, and not interpolated, as we found that interpolation
leads to large errors in simulations with large differences in dynamic
viscosity.

###  3.4 Consistency of VOF methods for Two-Phase Flows with High Density
Ratios

Here we summarize our findings that lead to an equation discretization with
flux-based geometrical VOF methods that remains consistent for very high
density ratios:

  * 1.

Computing the mass flux $\rho_{f}F_{f}$ by scaling the fluxed phase-specific
volume $|V_{f}^{\alpha}|_{s}$ with $(\rho^{-}-\rho^{+})/\Delta t$ approximates
only one constant average mass flux over $\Delta t$.

  * 2.

The constant mass flux $\rho_{f}F_{f}$ scaled from $|V_{f}^{\alpha}|_{s}$ with
$\Delta t$ , as an average value over $\Delta t$, can only be consistently
used in first-order schemes, i.e. use a single mass flux value over $\Delta
t$.

  * 3.

Computing the mass flux $\rho_{f}F_{f}$ from the density model eq. 2
disconnects mass conservation from volume fraction advection, which uses
geometrical integration, i.e. a geometrically integrated phase-specific
volumetric flux. This requires a solution of an additional (auxiliary) density
equation ([17, 2]) for cell-centered density, discarded at the end of the time
step.

  * 4.

Upwinding geometric VOF methods that use the Euler temporal scheme to
approximate point displacements are equivalent to using $\rho_{f}F_{f}$ from
eq. 30 and solving the density equation for $\rho_{c}^{o}$, used further in
$p-\mathbf{v}$ coupling to obtain a divergence-free cell-centered velocity. In
other words, combining Euler temporal integration scheme with upwind scheme
for the momentum equation guarantees numerical consistency for any flux-based
VOF method, which uses temporally first-order accurate displacements in its
geomterical integration of the fluxed phase-specific volume.

These points we verify and validate in the following section.

##  4 Verification and validation

Data archives of the implementation of the interIsoRhoFoam algorithm, input
data, post-processing software and secondary data are publicly available [26,
27]. The method is actively developed in a publicly available git repository
[28].

###  4.1 Time step size

The time step is limited by the CFL condition in the explicit plicRDF-
isoAdvector method [12],

$\Delta t_{CFL}\leq\dfrac{CFLh}{|\mathbf{v}|},$ (37)

where $h$ is the discretizaiton length, and we use $CFL=0.2$ from [12].
Another restriction for the time step size considers the propagation of
capillary waves on fluid interfaces,

$\Delta t_{cw}\leq\sqrt{\dfrac{(\rho^{+}+\rho^{-})h^{3}}{2\pi\sigma}}.$ (38)

This time step constraint was introduced first by Brackbill et al. [29], and
revised by Denner and van Wachem [30]. We restrict the time step using the
relation from Tolle et al. [24], namely

$\Delta t=argmin_{\Delta t_{CFL},\Delta t_{cw}}\,(k_{1}\Delta t_{cw},\
k_{2}\Delta t_{CFL})$ (39)

where $k_{1}$ and $k_{2}$ are scaling factors. We set $k_{1}=1$, $k_{2}=0.5$
as the default value.

###  4.2 Translating droplet in ambient flow

A canonical test case, originally introduced by Bussmann et al. [31] in 2D,
involves a moving droplet in quiescent ambient flow. Zuzio et al. [19]
extended the 2D case to 3D using a density ratio of $10^{6}$. Following the
setup in [19], the droplet of radius $R=0.15$ has the initial velocity of
$(0,0,10)$. To avoid the perturbation caused by sudden acceleration of still
ambient flow, the velocity field of ambient flow is initialized slightly
differently from [19], where the ambient flow velocity is set to zero. Here,
the initial velocity field in the ambient flow domain is defined using
potential theory from [32, sec. 10.3.1]

$\displaystyle v_{x}$
$\displaystyle=\frac{3|\mathbf{v}_{d}|r_{0}^{3}d_{z}d_{x}}{2|\mathbf{r}|^{5}},$
(40) $\displaystyle v_{y}$
$\displaystyle=\frac{3|\mathbf{v}_{d}|r_{0}^{3}d_{z}d_{y}}{2|\mathbf{r}|^{5}},$
(41) $\displaystyle v_{z}$
$\displaystyle=|\mathbf{v}_{\infty}|-\frac{|\mathbf{v}_{d}|r_{0}^{3}(|\mathbf{r}|^{2}-3.0d_{z}d_{z})}{2|\mathbf{r}|^{5}},$
(42)

where $\mathbf{v}_{a}^{0}(\mathbf{x})=(v_{x},v_{y},v_{z})$ indicates velocity
of points $\mathbf{x}$ outside the droplet, $|\mathbf{v}_{d}|=10$,
$|\mathbf{v}_{\infty}|=0$ denote magnitude of droplet velocity
$\mathbf{v}_{d}$ and undisturbed infinite velocity $\mathbf{v}_{\infty}$
correspondingly, $r_{0}$ is radius of droplet. It is assumed that
$\mathbf{r}_{p}$ and $\mathbf{r}_{d}$ are the respective positions of the
point outside the droplet and the droplet center, and we define a relative
position $\mathbf{r}=\mathbf{r}_{p}-\mathbf{r}_{d}=(d_{x},d_{y},d_{z})$. The
droplet with initial center location $(0.5\ 0.5\ 0.5)$ translates to a
distance of $L=1$ and the simulation time is then $t_{end}=0.1$. The
computational domain has dimensions $L_{z}=L_{x}=L_{y}=1$. The periodic
boundary condition is applied to all boundary patches. During testing with the
periodic boundary condition, we found that the plicRDF-isoAdvector method [12]
implemented in [33] has an inconsistency at the periodic boundary, which we
fixed as described in A. The mesh setup from [19] is also followed: the mesh
resolution is in a range of $N\in(32,\ 48,\ 64)$ per unit side-length of the
computational domain, resulting in $\approx(10,\ 15,\ 20)$ mesh cells per
droplet diameter. Surface tension force is neglected in this case, so only
$\Delta t_{CFL}$ from eq. 39 is taken into account. The viscosity and
gravitational forces are neglected to highlight the numerically consistent
behavior of the mass and momentum convection. We adopted three error norms to
evaluate the results quantitatively, for mass, momentum, and sphericity:

$\displaystyle E_{mass}$
$\displaystyle=\frac{M(t)-M(0)}{M(0)}=\frac{\sum_{k}m_{k}(t)-\sum_{k}m_{k}(0)}{\sum_{k}m_{k}(0)}=\frac{\sum_{k}\rho_{k}(t)V_{k}}{\sum_{k}\rho_{k}(0)V_{k}},$
(43) $\displaystyle E_{mom}$
$\displaystyle=\frac{|\sum_{k}m_{k}(t)\mathbf{v}_{k}(t)|-|\sum_{k}m_{k}(0)\mathbf{v}_{k}(0)|}{|\sum_{k}m_{k}(0)\mathbf{v}_{k}(0)|},$
(44) $\displaystyle E_{sph}$
$\displaystyle=\Psi_{d}(t)-\Psi_{d}(0),\quad\Psi_{d}(t)=\frac{\pi^{1/3}(6V_{d}(t))^{2/3}}{A_{d}(t)}$
(45)

where the subscript $k$ indicates that the value is extracted from the cell
$\Omega_{k}$ and the $m_{k}$, $V_{k}$ denote the mass and volume of
$\Omega_{k}$. In eq. 45, $V_{d}$, and $A_{d}$ represent the volume and surface
area of the droplet, respectively. These two ratios of eqs. 43 and 44
represent the time evolution of the normalized error of the global sum of the
heavy phase mass and momentum. In this case, $E_{mass}$ is expected to be near
zero because there is no source and dissipation for both mass, and flux-based
geometrical VOF method have a very high degree of local volume (mass)
conservation.

In the absence of force terms on both sides of momentum transport equation,
i.e, eq. 5, assuming the periodic boundary condition is applied to all
boundary patches, momentum conservation dictates that the deviation $E_{mom}$
should also theoretically remain at zero over time. Bussmann et al. [31]
proposed that a droplet, characterized by a large density ratio ($10^{6}$),
should undergo translation without deformation in an ambient flow, much like a
solid sphere moving through a void. This conclusion has been widely accepted
and corroborated by several publications (e.g., [34, 35, 17, 36, 37]).
However, these studies qualitatively assessed droplet deformation based on
visual representations of droplet shape. In this work, we employ the
sphericity $\Psi_{d}$ as a quantitative measure. Additionally, we use
$E_{sph}$ to characterize the deviation of deformation from the initial
droplet shape.

####  4.2.1 Achieving consistency with scheme selection

As discussed in section 3.4, when the first-order accurate Euler and Gauss
upwind scheme are employed to discretize momentum conservation eq. 5, the mass
flux scaled from $|V_{f}^{\alpha}|_{s}$ by $\Delta t$ in eq. 28 will be
consistent. Since the choice of discretization schemes ensures consistency of
the discretization, there is no need to modify the implementation of the
numerical method. We have verified the analysis from section 3 for the
”interIsoFoam” solver and compared it with the ”interIsoRhoFoam” solution
algorithm 1 that implements the auxiliary density equation. The normalized
mass error eq. 43 is near machine epsilon as depicted in fig. 4, for both
interIsoFoam and interIsoRhoFoam, showing excellent conservation of mass
during the simulation. The momentum conservation eq. 44 is equally very high
for both solution algorithms, as depicted in fig. 5. We therefore verify the
consistency of the ”interIsoFoam” with Euler and Gauss upwind scheme, with
”interIsoRhoFoam”, that solves the auxiliary density equation.

Note that this verification case is extremely challenging, since it is an
inviscid case - there is no viscous force available in this case to dampen the
errors resulting from inconsistent two-phase mass and momentum transport.

Next, we demonstrate inconsistencies leading to large errors and often to
catastrophic failure when more than one mass flux is used over $\Delta t$ when
discretizing eq. 5, or the mass flux is limited in the discretized eq. 5.

(a) interIsoFoam

(b) interIsoRhoFoam

Figure 4: Temporal evolution of normalized mass conservation error with Euler
and Gauss upwind.

(a) interIsoFoam

(b) interIsoRhoFoam

Figure 5: Temporal evolution of normalized momentum conservation error with
Euler and Gauss upwind.

####  4.2.2 Comparison of different schemes

Time schemes Gaussian convection schemes Order of convection accuracy Category
of convection Boundedness of convection Mass flux consistency Crank-
Nicolson[38] upwind first-order NVD/TVD Bounded no Euler upwind first-order
NVD/TVD Bounded yes cubic[39] second-order non-NVD/TVD Unbounded no
limitedLinearV first-/second-order NVD/TVD Unbounded no linear second-order
non-NVD/TVD Unbounded no LUST second-order non-NVD/TVD Unbounded no MUSCL[40]
second-order NVD/TVD Unbounded no QUICK[41] second-order NVD/TVD Unbounded no
SuperBee[42] second-order NVD/TVD Unbounded no vanLeer[43] second-order
NVD/TVD Unbounded no

Table 1: The combinations of different time and convection schemes used to
test the effect of numerical consistency.

Using the Crank-Nicolson scheme to discretize eq. 5 reveals that the temporal
scheme involving implicit mass flux $\rho_{f}F_{f}$ results in a mismatch
between mass convection and volume fraction convection scaled by
$(\rho^{-}-\rho^{+})$, due to the fact that the Navier–Stokes equation using a
segregated method is solved iteratively within a time step, and the
interface’s advection velocity thus cannot be updated simultaneously with
$\mathbf{v}^{n+1}$ from the previous $p-v$ coupling iteration. The
inconsistency is amplified by the density-ratio. Combinations of schemes
listed in table 1, are tested to verify their effect on the mass flux
inconsistency.

Figure 6: Temporal evolution of normalized mass conservation error with
different schemes: interIsoFoam, $N=64$. Figure 7: Temporal evolution of
normalized momentum conservation error with different schemes: interIsoFoam,
$N=64$.

Figures 6 and 7 represent the temporal evolution of $E_{mass}$, and $E_{mom}$.
It is evident from fig. 6 that mass conservation is not maintained when
utilizing the Crank-Nicolson time discretization scheme, and divergence
schemes cubic and LUST. These simulations terminate at an early stage with
catastrophic failure. However, for cases that can run until the final time,
the magnitude of mass errors, as shown in the zoomed subfigure of fig. 6, is
on the order of $10^{-17}$, indicating strict mass conservation. This
observation aligns with the inherent characteristic of the volume of fluid
method, which is known for its mass conservation property. The vertical lines
from the results of cubic, QUICK and Crank-Nicolson can be observed in fig. 7.
Some combinations deliver stable momentum errors. A closer examination of the
stable cases in the zoomed view of fig. 7 confirms the use of the Euler
temporal scheme.

Figure 8: Temporal evolution of sphericity error with different schemes:
interIsoFoam, $N=64$.

(a) Crank-Nicolson $+$ upwind

(b) Euler $+$ LUST

(c) Euler $+$ limitLinearV

(d) Euler $+$ upwind

Figure 9: Final shape of the droplet calculated by interIsoFoam with different
schemes: $N=64$.

The temporal evolution of sphericity error is depicted in fig. 8 where the
steepness of the curves indicates the extent of droplet deformation from its
initial shape. The results obtained from four unstable scheme combinations
show significantly larger deviations in sphericitiy. It is notable that
although the schemes Euler$+$limitedLinearV, linear and vanLeer, demonstrate
comparable performances compared with the consistent-on-paper Euler$+$upwind
in terms of mass and momentum conservation in figs. 6 and 7, sphericity errors
calculated using these schemes deviate more from the expected value. Moreover,
the zoomed view reveals that, apart from Euler$+$upwind, all other stable
scheme combinations display varying degrees of divergence. This suggests that
if the simulations were to continue for a longer duration, these cases would
likely crash. Figure 9 presents the final shapes of the droplets simulated
using interIsoFoam with two stable scheme combinations and two unstable
combinations. For the unstable combinations, a common occurrence is observed:
the original droplets disintegrate into small, irregularly scattered pieces. A
comparison between the final shapes obtained using the inconsistent but in
this verification case, for these parameters, still stable scheme combination
Euler $+$ limitedLinearV (fig. 9(c)) and the consistent scheme combination
Euler $+$ upwind (fig. 9(d)) corroborates the findings from the sphericity
errors in fig. 8: an evident ”crown” with irregular bumps forms on the top
part of the droplet in fig. 9(c), leading to a reduction in sphericity. On the
contrary, in fig. 9(d), only a slight shrinkage occurs in the neck region of
the droplet, while the top surface remains smooth.

Table 2 summarizes these findings, providing information on early termination
of scheme combinations that are inconsistent, as well as mass, momentum and
sphericity errors for different mesh resolutions.

interIsoRhoFoam interIsoFoam Scheme combination Resolution End time Mass error
Momentum error Sphericity Error End time Mass error Momentum error Sphericity
Error Euler+upwind 32 0.1 0.0 0.0 -1.525e-3 0.1 0.0 0.0 -1.525e-3 48 0.1 0.0
0.0 -9.492e-4 0.1 0.0 0.0 -9.492e-4 64 0.1 0.0 0.0 -1.147e-3 0.1 0.0 0.0
-1.147e-3 Euler+cubic 32 0.006 5.164e15 5.350e45 4.642e8 0.036 1.016e-4
3.055e-4 -5.6182e-2 48 0.002 9.664e17 6.739e58 1.6335e10 0.009 0.0 3.578e-5
-5.372e-3 64 0.001 0.2515 1.258e21 -0.9522 0.1 7.119e-6 49.984 -6.155e-2
Euler+limitedLinearV 32 0.053 0.0 1.451e-05 -6.864e-4 0.078 0.0 0.0 -1.268e-3
48 0.008 0.0 7.156e-06 -8.807e-5 0.012 0.0 0.0 -6.165e-4 64 0.002 0.0 0.0
4.746e-4 0.1 0.0 0.0 4.894e-3 Euler+linear 32 0.055 0.0 7.257e-06 -5.625e-4
0.079 0.0 0.0 -2.259e-3 48 0.008 0.0 1.431e-05 -6.850e-5 0.012 0.0 0.0
-5.871e-4 64 0.002 0.0 0.0 4.647e-4 0.1 0.0 0.0 -1.054e-2 Euler+LUST 32 0.009
235.22 7.324e+21 0.2461 0.047 7.257e-6 -1.451e-5 -5.388e-2 48 0.005 -1.723e-2
2.827e+24 -0.9615 0.1 0.0 2.082e-6 -3.797e-3 64 0.002 0.0 9.254e-05 5.042e-4
0.044 1.423e-5 -1.423e-5 -0.1345 Euler+MUSCL 32 0.053 0.0 5.080e-05 -8.676e-4
0.075 0.0 1.267e-15 -1.630e-3 48 0.008 0.0 1.431e-05 -4.208e-4 0.010 0.0
1.353e-14 -4.306e-4 64 0.002 0.0 1.864e-15 5.042e-4 0.1 0.0 1.882e-12
-6.308e-3 Euler+QUICK 32 0.026 0.0 3.332e7 -0.1908 0.072 0.0 1.793e-12
-1.830e-3 48 0.007 0.0 1.431e-05 -5.871e-4 0.010 0.0 3.124e-15 -4.110e-4 64
0.002 0.0 2.279e-15 5.141e-4 0.057 0.0 -0.1327 -6.605e-3 Euler+SuperBee 32
0.052 0.0 3.628e-05 -4.004e-4 0.074 0.0 2.534e-14 -1.592e-3 48 0.007 0.0
1.431e-05 -5.774e-4 0.010 0.0 3.124e-15 -2.348e-4 64 0.002 0.0 0.0 5.162e-4
0.1 0.0 2.019e-12 -7.890e-3 Euler+vanLeer 32 0.053 0.0 5.080e-05 -9.916e-4
0.074 0.0 8.659e-15 -9.153e-4 48 0.008 0.0 1.431e-05 -1.761e-4 0.010 0.0
2.082e-16 -2.740e-4 64 0.002 0.0 1.450e-15 4.845e-4 0.1 0.0 3.895e-14
-5.804e-3 CrankNicolson+upwind 32 0.004 0.0 7.853e-4 8.104e-4 0.075 0.0
-7.257e-06 -2.498e-3 48 0.001 0.0 6.440e-05 4.501e-4 0.010 0.0 0.0 -3.229e-4
64 0.066 6.001e-3 1533.82 -0.1537 0.066 6.001e-3 1533.82 -0.1537

Table 2: The terminated time and the final normalized errors of
interIsoRhoFoam and interIsoFoam with different scheme combinations and mesh
resolutions for translating droplet in ambient flow case.

###  4.3 Translating sub-millimeter droplet with realistic physical properties

materials/properties ($25\text{\,}\mathrm{\SIUnitSymbolCelsius}$) density
($\mathrm{kg}\text{/}{\mathrm{m}}^{3}$) kinematic viscosity
(${\mathrm{m}}^{2}\text{/}\mathrm{s}$) surface tension
($\mathrm{N}\text{/}\mathrm{m}$) density ratio air $1.1839$ $---$ $---$ $---$
[44] mercury $13.5336\times 10^{3}$ $---$ $---$ (in air) $11431.37$(in air)
[44]

Table 3: Realistic fluid properties of the mercury droplet/air ambient pair.

The realistic physical properties of the mercury droplet and air ambient pair
are selected for the verification case of the translating sub-millimeter
droplet. The density ratio displayed in table 3 is around $10^{4}$. The rest
setups are same as in [2]. A spherical droplet of radius $R=$
$0.25\text{\,}\mathrm{mm}$ translates a distance of three diameters with
velocity $0.01\text{\,}\mathrm{m}\text{/}\mathrm{s}$ in $z$-direction of the
rectangular solution domain ($L_{x}=L_{y}=5R,L_{z}=10R$). The ambient flow has
the same velocity, that is $\mathbf{v}_{a}=(0,0,0.01)$. Three resolutions are
tested in this case: $N\in(16,\ 32,\ 64)$. The initial centroid position of
the droplet is $(2.5R,2.5R,2R)$. Surface tension and viscous forces are not
considered in this case. Since the droplet translates with the ambient flow
and there is no sink or source for the droplet moving, the velocity field
should keep unchanged.

The error norm $L_{\infty}$ is employed to measure the maximal deviation
between the numerical velocity and the analytical one among all cells, i.e.,

$L_{\infty}(\mathbf{v})=\max_{i}\left(\frac{\|\mathbf{v}_{i}-\mathbf{v}_{\infty}\|}{\|\mathbf{v}_{\infty}\|}\right),$
(46)

where $\mathbf{v}_{i}$ denotes the velocity of the cell $i$, and the
analytical velocity value is $\mathbf{v}_{\infty}=\mathbf{v}_{a}=(0,0,0.01)$.
The expected exact value of $L_{\infty}$ is $0$.

(a) interIsoFoam

(b) interIsoRhoFoam

Figure 10: Temporal evolution of the velocity error norm
$L_{\infty}(\mathbf{v})$ with pure advection: Euler, Gauss upwind, density
ratio: $10^{4}$, mesh resolution: $N\in(16,\ 32,\ 64)$. Figure 11: Temporal
evolution of the velocity error norm $L_{\infty}(\mathbf{v})$ with pure
advection - combining $10$ schemes, density ratio is $10^{4}$, mesh resolution
is $N=64$. Only Euler and upwind(ing) schemes remain consistent and stable.

The fig. 10 presents the temporal evolution of $L_{\infty}(\mathbf{v})$. The
same magnitude of $L_{\infty}(\mathbf{v})$ calculated through both solvers
reveals the equivalence between the volume fraction and mass conservation
equation under the schemes’ combination of Euler+upwind. The equivalence
bridges the volume fraction convection and the momentum convection. Thus,
consistency is assured in interIsoFoam when deploying the schemes’
combination, even without solving an auxiliary mass equation. A notable
outcome from fig. 10 is the value of the final converged
$L_{\infty}(\mathbf{v})$, which is at the magnitude of $1\times 10^{-11}$ and
almost reaches the machine tolerance.

We also conducted tests on all schemes listed in table 1 for this particular
case setup. The corresponding results are presented in fig. 11. The results
from Crank-Nicolson temporal scheme show significant numerical instability at
the initial stage of the simulation. In contrast, the errors obtained from all
other schemes remain stable throughout the entire simulation. Notably, there
is a substantial variation in accuracy among these schemes. The velocity
errors from Euler $+$ QUICK and upwind remain at magnitudes around $10^{-11}$,
while the errors from the other schemes initially increase and stabilize at
values that are $10^{7}$ to $10^{10}$ times larger.

###  4.4 Mixing layer

$x$
---
$y$
---
$\rho^{+}$
---
$\rho^{-}$
---
$1.5\text{\,}\mathrm{mm}$
---
$-1.5\text{\,}\mathrm{mm}$
---
$0\text{\,}\mathrm{mm}$
---
$3\text{\,}\mathrm{mm}$
---
Figure 12: 2D mixing layer

In this case, a 2D mixing layer case is tested. The 2D computational domain as
depicted in fig. 12 has the same length $L=$3\text{\,}\mathrm{mm}$$ in both
$x-$ and $y-$direction. The liquid with high density
$\rho^{-}=$1000\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}$$ flows
in the middle region
$$-0.15\text{\,}\mathrm{mm}$<y<$0.15\text{\,}\mathrm{mm}$$ of the square
computational domain with relatively low initial velocity
$v_{x}^{-}=$2\text{\,}\mathrm{m}\text{/}\mathrm{s}$$, while the gas with the
density $\rho^{+}=$1\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}$$
flows on both sides of the liquid area with very high initial velocity
$v_{x}^{+}=$30\text{\,}\mathrm{m}\text{/}\mathrm{s}$$. A spatial velocity
perturbation is initialized in the internal field and has the distribution

$v_{y}(y)=0.01v^{-}_{x}\sin{2\pi\frac{x}{L}}\exp{-(\frac{2y}{h})^{2}},$

where $h$ indicates the thickness of the liquid region, i.e.
$0.3\text{\,}\mathrm{mm}$. The simulations are tested with a resolution of
$N_{x}\times N_{y}\times N_{z}=128\times 128\times 1$. The duration of the
simulation is set to $0.003\text{\,}\mathrm{s}$, allowing for sufficient
number of time steps for the inconsistencies to develop. To highlight the
dominant impacts of the convection, the surface tension force, gravity, and
viscosity are excluded from the simulations. The periodic condition is
employed for all boundaries. We compared our results with the results computed
by the ONERA DYJEAT codes [45, 46, 47, 48, 19], which upholds consistent mass-
momentum transport through solving the temporary density equations together
with momentum equations on staggered meshes.

Considering the equivalence, we anticipate similarities between the results
obtained from interIsoFoam(Euler+upwind) and interIsoRhoFoam. Figure 13
illustrates the final state of the computational domain. Both solvers
successfully complete the simulation until the designated final time of
$t_{\text{end}}=$0.003\text{\,}\mathrm{s}$$, without any simulation breakdown.
The vorticity and velocity fields shown in fig. 13(a) and fig. 13(b)
respectively display very similar distributions. Furthermore, the contour
lines of the iso-value $\alpha=0.5$, represented by the black lines in the
sub-figures, also demonstrate very similar shapes. The consistent scheme in
DYJEAT generates similar results in terms of momentum conservation; however,
the liquid column ruptured into irregular liquid parts at the final state, as
shown in fig. 13(c), which manifests slightly stronger instabilities of the
results. The fig. 14 provides a quantitative comparison among multiple schemes
in table 1 and with the DYJEAT code, focusing on the temporal evolution of the
normalized momentum error evaluated using eq. 44. From the plot, it is evident
that only three cases using interIsoFoam remain stable, namely, as expected
using Euler $+$ upwind, but also Euler $+$ LUST and Euler $+$ vanLeer. The
zoomed-in view in fig. 14 highlights the accuracy of the results for these
four stable cases. The errors calculated from DYJEAT are larger than errors
from stable cases using interIsoFoam, which are around $4\%$. As shown in the
detail in fig. 14, the errors from consistent Euler $+$ are minimal, with
respect to all other combinations of schemes.

(a) interIsoFoam (Euler+upwind)

(b) interIsoRhoFoam (Euler+upwind)

(c) DYJEAT

Figure 13: The state of the mixing layer at the end of simulation time
$t=0.003s$: interIsoFoam and interIsoRhoFoam with Euler and Gauss upwind and
DYJEAT codes, density ratio: $10^{3}$, resolution: $N=128$. Figure 14: Time
evolution of normalized momentum error of mixing layer with different schemes
and DYJEAT codes, density ratio: $10^{3}$, resolution: $N=128$.

###  4.5 Validation of a single rising bubble

In the present study, we investigate the performance of the proposed method by
applying it to the simulation of a single bubble rising in a quiescent viscous
liquid. To validate our approach, we adopt the configuration presented by
Anjos et al. [49], who simplified the rising bubble experiments originally
conducted by Bhaga and Weber [50]. In their work, Anjos et al. [49] select
three distinct viscosities for comparative analysis. Our focus, in particular,
lies on the cases, characterized by larger liquid viscosities. The specific
cases correspond to a Morton number
$Mo=g\nu_{l}^{4}/\rho_{l}\sigma^{3}=(848,41.1)$, where $g$ represents the
gravitational acceleration, and $\nu_{l}$, $\rho_{l}$, and $\sigma$ denote the
viscosity, density of the ambient liquid, and surface tension coefficient,
respectively. Three resolutions are tested: $N\in(64,96,128)$.

For our simulations, we initialize the air bubble with an idealized spherical
shape, possessing a diameter of $D=$2.61\text{\,}\mathrm{cm}$$. The air
properties are defined by a viscosity of
$1.78\times${10}^{-5}\text{\,}\mathrm{kg}\text{/}\text{(}\mathrm{m}\text{\,}\mathrm{s}\text{)}$$
and a density of $1.225\text{\,}\mathrm{kg}\text{/}{\mathrm{m}}^{3}$, while
the liquid properties encompass viscosities of
$(2.73,1.28)\,$\mathrm{kg}\text{/}\text{(}\mathrm{m}\text{\,}\mathrm{s}\text{)}$$
and a density of $1350\text{\,}\mathrm{kg}\text{/}{\mathrm{m}}^{3}$.
Furthermore, the surface tension between the air bubble and the liquid is
$0.078\text{\,}\mathrm{N}\text{\,}\mathrm{m}$. The computational domain is
defined as $(-4D,-4D,-2D)\times(4D,4D,6D)$, where the positions of the space
diagonal vertices of the computational domain are delineated, with the initial
position of the bubble set at the origin, $(0,0,0)$. A set of dimensionless
normalized variables is introduced as follows:

$\mathbf{w}=\frac{\mathbf{v}}{\sqrt{gD}},\ t=\sqrt{\frac{g}{D}}\tau,$ (47)

where $\tau$ indicates the physical time in seconds.

(a) interIsoFoam: Mo$=41.4$

(b) interIsoFoam: Mo$=848$

(c) interIsoRhoFoam: Mo$=41.4$

(d) interIsoRhoFoam: Mo$=848$

Figure 15: Temporal evolution of rising velocity using interIsoFoam and
interIsoFoam: Euler $+$ upwind, $\rho^{-}/\rho^{+}\approx 10^{3}$.

(a) interIsoFoam: Mo$=41.4$

(b) interIsoFoam: Mo$=848$

(c) interIsoRhoFoam: Mo$=41.4$

(d) interIsoRhoFoam: Mo$=848$

Figure 16: Comparisons of final shapes of rising bubbles using interIsoFoam
and interIsoRhoFoam with the experimental visualization from Bhaga and Weber
[50] (reprinted with permission): Euler $+$ upwind, $\rho^{-}/\rho^{+}\approx
10^{3}$, $N=128$, $t=6$.

As illustrated in fig. 15, the utilization of the Euler and upwind schemes
ensures the preservation of equivalence between volume fraction and mass
advection equation. Consequently, the results obtained from interIsoFoam and
interIsoRhoFoam, when employing the identical Morton number, display a
substantial level of similarity. When considering Mo $=41.1$, as demonstrated
in fig. 15(a) and fig. 15(c), the velocities acquired from both solvers
demonstrate an initial increase until approximately $t=1$, followed by a
subsequent decline leading to a stable state. Notably, the differences among
the results become evident when dealing with the cases involving Mo $=41.1$.
Specifically, in instances where a coarse mesh ($N=64$) is employed, the
velocity decline is more pronounced, resulting in a smaller final stable
velocity compared to the values presented in prior works [50, 51, 49].
Conversely, higher resolutions yield stable velocities that agree well with
the experimental data documented by Bhaga and Weber [50]. On the other hand,
for cases with Mo $=848$, characterized by a larger viscosity, the impact of
resolution is less conspicuous, as depicted in fig. 15(b) and fig. 15(d). The
final velocities attained from different resolutions converge similarly to
values that fall between the results obtained in the simulation conducted by
Hua and Lou [51] and the experimental study conducted by Bhaga and Weber [50].

Furthermore, the profiles obtained by slicing the surfaces of the droplets,
passing through their centers, are compared with the experimental profiles
from Bhaga and Weber [50], as illustrated in fig. 16. For both solvers, the
droplets with varying viscosities exhibit final shapes that closely align with
the experimental visualizations, validating the hypothesis of the equivalence
between Euler+upwind discretization of eq. 5 using the scaled mass flux
$|V_{f}^{\alpha}|_{s}$ from eq. 28, and the solution of the auxiliary density
equation using algorithm 1.

###  4.6 Liquid jet in high speed gaseous cross-flow

Different from the parallel velocities of the mixing layers in section 4.4,
the liquid flows with a lower velocity is perpendicular to the velocity of the
gaseous phase in this case, which is called the injection of a liquid jet in a
gaseous cross-flow (LJCF) and is common in many engineering applications. The
geometry and the physical properties are configured by referring to Zuzio et
al. [19]. The rectangular computational domain
$\Omega:[-0.01,0,-0.01]\times[0.03,0.02,0.01]\,$\mathrm{m}$$ has two inlets.
The gas flows in with a velocity
$\mathbf{v}^{+}=[65,0,0]\,$\mathrm{m}\text{/}\mathrm{s}$$ from the left
boundary $x_{min}$. Thuillet [52] revealed the impact of the liquid inlet
velocity profile on the jet trajectory. He simulated the jet with an uniform
liquid inlet velocity profile, and with a velocity profile calculated through
simulating the injector. The jet trajectory results from the case with
calculated velocity profile showed better agreement with the experiment. We
followed the calculated liquid injected velocity profile from [52], which is

$\mathbf{v}^{-}_{y}=-21.434\left(\frac{r}{d}\right)^{3}+15.512\left(\frac{r}{d}\right)^{2}+8.6504$

, where $\mathbf{v}^{-}_{y}$ is the $y-$component of the liquid inlet velocity
$\mathbf{v}^{-}$, $r$ indicates the distance to the nozzle center, and $d$ is
the diameter of the nozzle. The $x$, $z-$component of $\mathbf{v}^{-}$ are set
to zero, whereas $d=$0.002\text{\,}\mathrm{m}$$. The nozzle’s center locates
at $[0,0,0]$ in the bottom boundary patch $y_{min}$. To save the computation
resource, the uniform Cartesian mesh with a moderate resolution of
$[N_{x},N_{y},N_{z}]=[128,64,64]$ has been adopted. Figure 17 depicts the flow
domain. The physical properties are
$\rho^{-}=$1000\text{\,}\mathrm{kg}\text{/}{\mathrm{m}}^{3}$$,
$\rho^{+}=$1.225\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}$$,
$\mu^{-}=1.0\times${10}^{-3}\text{\,}\mathrm{kg}\text{/}\text{(}\mathrm{m}\text{\,}\mathrm{s}\text{)}$$,
$\mu^{+}=1.78\times${10}^{-5}\text{\,}\mathrm{kg}\text{/}\text{(}\mathrm{m}\text{\,}\mathrm{s}\text{)}$$,
$\sigma=7.2\times${10}^{-2}\text{\,}\mathrm{N}\text{\,}\mathrm{m}$$,
$g=$9.81\text{\,}\mathrm{m}\text{/}{\mathrm{s}}^{2}$$.

$-0.01$
---
$0.03$
---
$-0.01$
---
$0.01$
---
$0$
---
$0.02$
---
$d$
---
$x$
---
$z$
---
$y$
---
$\rho^{-}$
---
$\rho^{+}$
---
Figure 17: Liquid in a cross flow.

Figure 19 shows the final state of the injected liquid at
$t=$7.1\text{\,}\mathrm{ms}$$ using interIsoFoam and interIsoRhoFoam with
Euler and Gauss upwind regarding the iso-value of the reconstructed distance
function $RDF=0$, as well as using DYJEAT [19] with a high resolution of
$[N_{x},N_{y},N_{z}]=[1024,512,512]$, whose shape is rendered by the iso-value
$0.5$ of the volume fraction. Figure 19(a) and Figure 19(b) illustrate the
ruptured liquid jet from the $y$-side view. The droplets’ distributions of
interIsoFoam and interIsoRhoFoam display many comparable characteristics. We
can observe two strips of droplets and a strip of bag-like liquid structure.
The outer liquid segregates into two $yz$-plane symmetric strips of droplets
at an early stage, i.e., at a low penetration height, as shown in the right
subfigure of fig. 19(a) and fig. 19(b). These droplets translate with the gas
along the $x$-direction and spread spanwise in the $y$-direction. The
remaining center liquid has a wavelike detachment and forms the bag-like
structure in the middle of the strips. Figure 19(c) demonstrates the results
from DYJEAT codes with a higher resolution. Similar liquid distributions can
be observed: the liquid in the center zone of the nozzle propagates like a
wave and breaks up into some large packets at a higher position, whereas the
liquid in the periphery of the nozzle zone rips at a low penetration height.

This validation case is a candidate for a benchmark case for validating two-
phase flow numerical methods that consistently handle high density ratios,
because the experimental form of the jet can be accurately reproduced on
coarser mesh resolutions. Using a coarser mesh resolution interIsoFoam and
interIsoRhoFoam of course do not capture the small structures such as liquid
streaks, sacs and droplets, as shown in DYJEAT’s results. However, the solvers
accurately predict the jet curve, which can be used as a quantifiable argument
for validity of a consistent method against experimental data.

The fig. 18 displays the final bent shape of liquid jet simulated by three
solvers and their comparison with the experimental observation made by ONERA
[53, 54]. The same case with two different resolutions are tested. We put the
liquid jets’ shape results in the same parallel view to compare them with each
other and also with the experimental results marked by the red line. The blue
translucent liquid jet represents the results from the case with a higher
resolution $N_{h}=[254,128,128]$, while the gray liquid jet comes from the
above low-resolution results. It can be seen from each figs. 18(a), 18(b) and
18(c) that more small droplets and complex structures can be captured when
deploying the higher mesh resolution. Despite the different resolutions, there
exists a very minor difference between the liquid jets with regard to the bent
shape. The windward surfaces of the two liquid jets in each sub-figure almost
attach to each other, which highlights that this case is insensitive to the
mesh resolution. As to the comparison with the experimental trajectories, both
the jets’ bent surface in figs. 18(a) and 18(b) show good correspondence to
the experimental shape, i.e. the red line at the low penetration height
$<$9\text{\,}\mathrm{mm}$$. The jets reattach to the red line in the upper-
right zone. The maximal deviation between the simulated jets and the
experiment shape is around $1\text{\,}\mathrm{m}\mathrm{m}$, which is $5\%$ of
the jet height. A more obvious discrepancy between the jets and the experiment
is shown in fig. 18(c). The liquid jet bent less than the experiment in the
high-speed flow after the given time, which results in a wrong prediction of
the impingement position.

(a) interIsoFoam

(b) interIsoRhoFoam

(c) DYJEAT

Figure 18: The instantaneous liquid jet shape at final time $t=7.1ms$ with
different resolutions (the blue translucent jet: $N_{h}=[256,128,128]$; the
gray jet: $N_{l}=[128,64,64]$) and its comparison with the experimental
results(the red line). This case makes it possible to evaluate performance on
coarser meshes as resolving finer structures does not impact the jet
trajectory.

(a) interIsoFoam: $t=7.1ms$

(b) interIsoRhoFoam: $t=7.1ms$

(c) DYJEAT: $t=7.07ms$

Figure 19: The shape of the injected liquid with interIsoFoam and
interIsoRhoFoam (Euler and Gauss upwind, density ratio: $816$, CFL number:
$CFL=0.2$, resolution: $N_{l}=[128,64,64]$), and with DYJEAT (resolution:
$N=[1024,512,512]$ [19])

##  5 Conclusions

We analyze the equivalence between the mass conservation and volume fraction
advection equation in the context of flux-based VOF methods. When an average
mass flux is computed on the discrete level by scaling the fluxed phase-
specific volume with the time step, we find that the equivalence is possible
only when the first-order Euler temporal discretization is used for momentum
conservation, no flux limiting is applied, and the flux-based VOF scheme uses
first-order quadrature for the integration of the fluxed phase-specific
volume. For any flux-based VOF method, the integration of the fluxed phase-
specific volume $V_{f}^{\alpha}$ therefore plays a key role in the stability
of two-phase flow simulations with high density ratios, when the mass flux is
estimated from it. Altering the mass flux $\rho_{f}F_{f}$ in the discrete
momentum conservation equation by applying flux limiting schemes or blending
schemes causes errors that are proportional to the density difference, and
lead either to large errors in the interface shape and topology, or
catastrophic failure. In other words, upwinding the two-phase momentum in a
flux-based VOF method is consistent to the geometrical upwinding of the volume
fraction by the flux-based VOF method.

We successfully apply the $\rho$LENT method for high-density ratio flows to
the isoAdvector-plicRDF Volume-of-Fluid method [12]. The adaption of the
$\rho$LENT method in the context of the geometrical Volume-of-Fluid method is
straightforward, requiring only the geometrical calculation of the upwind area
fraction $\alpha_{f}^{o}$ from the available geometric VOF interface (phase-
indicator) approximation.

We demonstrate the equivalence on the discrete level between the mass flux
scaled from the phase-specific volume and the solution of an auxiliary density
equation. We confirm this hypothesis by means of verification with challenging
inviscid cases, as well as by validation against experiments.

##  6 Acknowledgments

The last author acknowledges funding by the German Research Foundation (DFG) –
Project-ID 265191195 – SFB 1194. Calculations for this research were conducted
on the Lichtenberg high performance computer of the TU Darmstadt.

##  Appendix A Correct cyclic boundary condition for the plicRDF-isoAdvector
method

$\mathbf{x}_{i}^{own}$
---
$\mathbf{x}_{c}^{nei}$
---
$\mathbf{x}_{i}^{nei}$
---
$\mathbf{x}_{f}^{nei}$
---
$\mathbf{x}^{own}_{f}$
---
$\mathbf{v}$
---
owner-patch
---
neighbor-patch
---
RDF
---
RDF
---
Figure 20: The ghost cells layer to correct signed distance: the green line (
) presents fixed RDF from cell centroid $\mathbf{x}_{c}^{nei}$ to interface in
the cyclic neighbor cell, while the red line ( ) depicts the original unfixed
RDF.

A cyclic boundary condition (BC) treats two boundary patches as if they were
physically connected, with their respective cell layers placed next to each
other. For the geometrical VOF method such as the plicRDF-isoAdvector, the
cyclic boundary condition impacts the interface reconstruction and the volume
fraction advection. To achieve this, the cyclic BC performs calculations on
the so-called owner-patch, and then reflects the result to the so-called
neighbor-patch, as shown schematically in fig. 20 for a geometric VOF method
in two dimensions.

For the advection discretization at cyclic BC, we noticed that the fluxed
phase-specific volumes $V_{f}^{\alpha}$ field is not properly adjusted for the
cyclic patches.

In the first step, the cyclic BC initializes $V_{f}^{\alpha}$ at every face of
both cyclic patches using the upwind scheme, i.e.,

$V_{f}^{\alpha,init}=F_{f}\alpha_{U}\Delta t,$ (48)

where $F_{f}:=\mathbf{v}_{f}\cdot\mathbf{S}_{f}$ is the volumetric flux at the
centroid of the cell-face $S_{f}$, $\alpha_{U}$ is the volume fraction of the
upwind cell $U$ w.r.t the cell face $S_{f}$, and $\Delta t$ is the time step.

In the second step, the cyclic BC computes the geometric $V_{f}^{\alpha}$ at
the faces that belong to interface cells on both patches

$V_{f}^{\alpha}=\int_{t^{n}}^{t^{n+1}}\int_{S_{f}}\chi\mathbf{v}\cdot\mathbf{n}\,dS\,dt=\int_{t^{n}}^{t^{n+1}}\frac{F_{f}(t)}{|\mathbf{S}_{f}|}A_{f}(t)\,dt$
(49)

The existing cyclic BC [55] does not consider cyclic boundary conditions. The
discretization in cell layers adjacent to two cyclic boundary patches should
handle the cell layers as if they are placed next to each other as shown in
fig. 20. A phase-specific volume fluxed out of the domain on a face that
belongs to the cyclic owner-patch, should be fluxed into the corresponding
face in the cyclic-neighbor patch, as described in algorithm 3.

Algorithm 3 The modified fluxed phase-specific volumes $V_{f}^{\alpha}$ update
method in interIsoFoam.

1:Initialize $V_{f}^{\alpha,owner}$ using the upwind scheme. $\triangleright$
eq. 48

2:for all boundary patches do

3: if boundary patch is cyclic then

4: for all cyclic-patch faces $f\in[1,|P_{cyclic}|]$ do

5: if $F_{f}>0$ then

6: Geometrically compute the outflow $V_{f}^{\alpha}$.

7: else if $F_{f}<0$ then $\triangleright$ The inflow $F_{f}<0$ here is
outflow $F_{f}>0$ of the neighbor.

8: $V_{f}^{\alpha}=-V_{f}^{\alpha,neigbor}$

9: end if

10: end for

11: end if

12:end for

Aside from the consistent calculation of $V_{f}^{\alpha}$ in the advection,
the cyclic BC allso affects the geometric interface reconstruction. The
plicRDF reconstruction [12] uses the reconstructed distance function (RDF) and
its gradient for improving discrete interface-normal vectors. When considering
the cyclic boundary condition, the distance calculation should be treated
carefully, especially in the case where two cells sharing a vertex are also
attached to two cyclic patches. This issue is illustrated in fig. 20, where
$\mathbf{x}_{i}^{own}$ denotes a center of a VOF interface polygon located in
a corresponding interface cell $\Omega_{i}^{own}$ that belongs to the cell
layer attached to the owner patch of the cyclic BC. In the plicRDF
implementation in OpenFOAM-v2306 [55], the cyclic BC falsely uses VOF
interface polygon centers $\mathbf{x}_{i}^{own}$ to compute signed distances
at centers $\mathbf{x}_{c}^{nei}$ in the cell layer adjacent to the cyclic
neighbor-patch. A correct implementation of the cyclic BC requires positions
$\mathbf{x}_{i}^{nei}$, as shown in fig. 20. The VOF interface centroids from
the cell layer adjacent to the cyclic owner-patch $\mathbf{x}_{i}^{own}$ and
the face centers of the cyclic owner-patch can be used to compute
$\mathbf{x}_{i}^{nei}$ for every $\mathbf{x}_{i}^{own}$ and facilitate a
correct cyclic (periodic) computation of the Reconstructed Distance Function
(RDF) in the plicRDF-isoAdvector method.

We define a transformation of the position of the interface centers using

$\mathbf{x}_{i}^{nei}:=\mathbf{x}_{i}^{own}+(\mathbf{x}_{f}^{nei}-\mathbf{x}_{f}^{own}),$
(50)

where $(\mathbf{x}_{f}^{nei}-\mathbf{x}_{f}^{own})$ is the difference
(displacement, or transformation) vector between the face centers of cyclic BC
owner and neighbor patch face centers.

Note that there is no need to tranform the interface normals $\mathbf{n}_{i}$,
because they are orientation and not position vectors. The transformed PLIC
polygon centroids together with the cyclic BC neighbor-patch interface normals
$\mathbf{n}^{nei}_{i}=\mathbf{n}^{own}_{i}$ build a cyclic ghost-data layer
that is shown schematically by dashed cells in fig. 20.

Signed distances in the cell layer adjacent to the cyclic owner-patch at cell
centers $\mathbf{x}_{i}^{own}$, are thus computed from the transformed PLIC
interface information from the cyclic patch-owner data.

If the cyclic-patch-adjacent cell $\Omega_{c}$ does contain its own PLIC
interface with the PLIC centroid, then the signed distance to this PLIC
interface is used, if it is closer to $\mathbf{x}_{i}$ of that cell.

We verify our discretization of the cyclic BC in the plicRDF-isoAdvector
method using a constant flow case shown in fig. 21, where both the droplet and
the ambient flow have the same initial velocity $\mathbf{v}=(v_{x},0)$. The
left and right boundaries are set as cyclic (periodic). It is noteworthy that
the momentum equation eq. 5 is not solved in this test case, ensuring exactly
constant velocity and pressure over time, aiming to only to test the periodic
(cyclic) interface reconstruction and advection. We test two interface
reconstruction methods, e.g. isoAlpha and plicRDF from [12].

Figure 22 shows the droplet reaching the right cyclic boundary. With the
erroneous calculation of the fluxed phase-specific volume $V_{f}^{\alpha}$ at
the cyclic boundaries, as shown in fig. 22(a), the interfaces appear in the
patch neighbor cells incorrectly. The fig. 22(b) show the accurate result of
our modification with a single interface cell. The modification of
$V_{f}^{\alpha}$ impacts both isoAlpha and plicRDF methods, while adapting the
displacement is crucial only for the plicRDF reconstruction. As shown in fig.
23(a), without modifying the displacement vector, some liquid remains in the
cell layer adjacent to the neighbor-patch after the droplet crosses the cyclic
boundary, and reaches a location far from both cyclic boundaries. With
applying our modification from fig. 23(b), the droplet retains its initial
form. The fixed cyclic boundary condition is available in [26].

Figure 21: A droplet moves with the ambient constant flow.

(a) Unfixed

(b) Fixed

Figure 22: The alpha field and interfaces reconstructed by iso-Alpha method
reach to the right cyclic boundary; blue region: the ambient flow; red region:
the liquid droplet; white line segments: the PLIC interfaces.

(a) Unfixed

(b) Fixed

Figure 23: The alpha field and interfaces reconstructed by plic-RDF method
cross the right cyclic boundary; blue region: the ambient flow; red region:
the liquid droplet; white line segments: the PLIC interfaces.

## References

  * Marić et al. [2020] T. Marić, D. B. Kothe, D. Bothe,  Unstructured un-split geometrical volume-of-fluid methods–a review,  Journal of Computational Physics 420 (2020) 109695. URL: https://doi.org/10.1016/j.jcp.2020.109695.
  * Liu et al. [2023] J. Liu, T. Tolle, D. Bothe, T. Marić,  An unstructured finite-volume level set/front tracking method for two-phase flows with large density-ratios,  Journal of Computational Physics 493 (2023) 112426. URL: https://doi.org/10.1016/j.jcp.2023.112426.
  * Roenby et al. [2016] J. Roenby, H. Bredmose, H. Jasak,  A computational method for sharp interface advection,  Royal Society Open Science 3 (2016). URL: https://arxiv.org/abs/1601.05392. doi:10.1098/rsos.160405. arXiv:1601.05392.
  * Huang et al. [2019] Z. Huang, G. Lin, A. M. Ardekani,  A mixed upwind/central WENO scheme for incompressible two-phase flows,  Journal of Computational Physics 387 (2019) 455–480. URL: https://doi.org/10.1016/j.jcp.2019.02.043. doi:10.1016/j.jcp.2019.02.043.
  * Jiang and Shu [1996] G.-S. Jiang, C.-W. Shu,  Efficient implementation of weighted eno schemes,  Journal of computational physics 126 (1996) 202–228. URL: https://doi.org/10.1006/jcph.1996.0130. doi:10.1006/jcph.1996.0130.
  * Dodd and Ferrante [2014] M. S. Dodd, A. Ferrante,  A fast pressure-correction method for incompressible two-fluid flows,  Journal of Computational Physics 273 (2014) 416–434. URL: https://doi.org/10.1016/j.jcp.2014.05.024. doi:10.1016/j.jcp.2014.05.024.
  * Xie et al. [2020] B. Xie, P. Jin, Y. Du, S. Liao,  A consistent and balanced-force model for incompressible multiphase flows on polyhedral unstructured grids,  International Journal of Multiphase Flow 122 (2020) 103125. URL: https://www.sciencedirect.com/science/article/pii/S0301932219304501. doi:https://doi.org/10.1016/j.ijmultiphaseflow.2019.103125.
  * Xie and Xiao [2017] B. Xie, F. Xiao,  Toward efficient and accurate interface capturing on arbitrary hybrid unstructured grids: The thinc method with quadratic surface representation and gaussian quadrature,  Journal of Computational Physics 349 (2017) 415–440. URL: https://www.sciencedirect.com/science/article/pii/S0021999117305995. doi:https://doi.org/10.1016/j.jcp.2017.08.028.
  * Francois et al. [2006] M. M. Francois, S. J. Cummins, E. D. Dendy, D. B. Kothe, J. M. Sicilian, M. W. Williams,  A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework,  Journal of Computational Physics 213 (2006) 141–173. URL: https://www.sciencedirect.com/science/article/pii/S0021999105003748. doi:https://doi.org/10.1016/j.jcp.2005.08.004.
  * Desmons and Coquerelle [2021] F. Desmons, M. Coquerelle,  A generalized high-order momentum preserving (homp) method in the one-fluid model for incompressible two phase flows with high density ratio,  Journal of Computational Physics 437 (2021) 110322. URL: https://www.sciencedirect.com/science/article/pii/S0021999121002175. doi:https://doi.org/10.1016/j.jcp.2021.110322.
  * Wang and Spiteri [2007] R. Wang, R. J. Spiteri,  Linear instability of the fifth-order weno method,  SIAM Journal on Numerical Analysis 45 (2007) 1871–1901. URL: https://epubs.siam.org/doi/abs/10.1137/050637868. doi:10.1137/050637868.
  * Scheufler and Roenby [2019] H. Scheufler, J. Roenby,  Accurate and efficient surface reconstruction from volume fraction data on general meshes,  Journal of computational physics 383 (2019) 1–23. URL: https://doi.org/10.1016/j.jcp.2019.01.009.
  * Scheufler and Roenby [2023] H. Scheufler, J. Roenby,  Twophaseflow: A framework for developing two phase flow solvers in openfoam,  OpenFOAM® Journal 3 (2023) 200–224. URL: https://doi.org/10.51560/ofj.v3.80.
  * Marić et al. [2021] T. Marić, J. Höpken, K. G. Mooney, The OpenFOAM Technology Primer, Zenodo, 2021. URL: https://doi.org/10.5281/zenodo.4630596. doi:10.5281/zenodo.4630596.
  * Jasak [1996] H. Jasak,  Error analysis and estimation for the finite volume method with applications to fluid flows. (1996). URL: https://spiral.imperial.ac.uk/handle/10044/1/8335.
  * Moukalled et al. [2016] F. Moukalled, L. Mangani, M. Darwish, et al., The finite volume method in computational fluid dynamics:An Advanced Introduction with OpenFOAM® and Matlab®, volume 6, Springer, 2016, pp. 443–450. URL: https://link.springer.com/book/10.1007/978-3-319-16874-6.
  * Ghods and Herrmann [2013] S. Ghods, M. Herrmann,  A consistent rescaled momentum transport method for simulating large density ratio incompressible multiphase flows using level set methods,  Physica Scripta 88 (2013). URL: https://iopscience.iop.org/article/10.1088/0031-8949/2013/T155/014050. doi:10.1088/0031-8949/2013/T155/014050.
  * Fuster et al. [2018] D. Fuster, T. Arrufat, M. Crialesi-Esposito, Y. Ling, L. Malan, S. Pal, R. Scardovelli, G. Tryggvason, S. Zaleski,  A momentum-conserving, consistent, Volume-of-Fluid method for incompressible flow on staggered grids,  Computers & Fluids (2018). URL: http://arxiv.org/abs/1811.12327. arXiv:1811.12327.
  * Zuzio et al. [2020] D. Zuzio, A. Orazzo, J. L. Estivalèzes, I. Lagrange,  A new efficient momentum preserving Level-Set/VOF method for high density and momentum ratio incompressible two-phase flows,  Journal of Computational Physics 410 (2020) 109342. URL: https://doi.org/10.1016/j.jcp.2020.109342. doi:10.1016/j.jcp.2020.109342.
  * Pal et al. [2021] S. Pal, D. Fuster, S. Zaleski, A novel momentum-conserving, mass-momentum consistent method for interfacial flows involving large density contrasts, 2021. URL: https://arxiv.org/abs/2101.04142. doi:10.48550/ARXIV.2101.04142.
  * Arrufat et al. [2021] T. Arrufat, M. Crialesi-Esposito, D. Fuster, Y. Ling, L. Malan, S. Pal, R. Scardovelli, G. Tryggvason, S. Zaleski,  A mass-momentum consistent, Volume-of-Fluid method for incompressible flow on staggered grids,  Comput. Fluids 215 (2021) 104785. URL: https://doi.org/10.1016/j.compfluid.2020.104785. doi:10.1016/j.compfluid.2020.104785.
  * Liu et al. [2022] C. Liu, R. Gao, C. Hu,  A consistent mass–momentum flux computation method for the simulation of plunging jet,  Physics of Fluids 34 (2022) 032114. URL: https://doi.org/10.1063/5.0084894. doi:10.1063/5.0084894.
  * Jin et al. [2022] Q. Jin, D. Hudson, W. G. Price,  A combined volume of fluid and immersed boundary method for modeling of two-phase flows with high density ratio,  Journal of Fluids Engineering 144 (2022). URL: https://doi.org/10.1115/1.4052242.
  * Tolle et al. [2020] T. Tolle, D. Bothe, T. Marić,  Saample: A segregated accuracy-driven algorithm for multiphase pressure-linked equations,  Computers & Fluids 200 (2020) 104450. URL: https://doi.org/10.1016/j.compfluid.2020.104450.
  * Marić et al. [2018] T. Marić, H. Marschall, D. Bothe,  An enhanced un-split face-vertex flux-based vof method,  Journal of computational physics 371 (2018) 967–993. URL: https://doi.org/10.1016/j.jcp.2018.03.048.
  * Maric and Liu [2023a] T. Maric, J. Liu, Inconsistencies in unstructured volume-of-fluid methods for two-phase flows with high density ratios - code, 2023a. URL: https://tudatalib.ulb.tu-darmstadt.de/handle/tudatalib/4018. doi:10.48328/tudatalib-1276.
  * Maric and Liu [2023b] T. Maric, J. Liu, Inconsistencies in unstructured volume-of-fluid methods for two-phase flows with high density ratios - data, 2023b. URL: https://tudatalib.ulb.tu-darmstadt.de/handle/tudatalib/4019. doi:10.48328/tudatalib-1277.
  * Marić and Liu [2023] T. Marić, J. Liu, Inconsistencies in unstructured volume-of-fluid methods for two-phase flows with high density ratios - code repository, 2023. URL: https://github.com/tmaric/TwoPhaseFlow/tree/feature/density-ratio.
  * Brackbill et al. [1992] J. U. Brackbill, D. B. Kothe, C. Zemach,  A continuum method for modeling surface tension,  Journal of computational physics 100 (1992) 335–354. URL: https://doi.org/10.1016/0021-9991(92)90240-Y.
  * Denner and van Wachem [2015] F. Denner, B. G. van Wachem,  Numerical time-step restrictions as a result of capillary waves,  Journal of Computational Physics 285 (2015) 24–40. URL: https://doi.org/10.1016/j.jcp.2015.01.021.
  * Bussmann et al. [2002] M. Bussmann, D. B. Kothe, J. M. Sicilian,  Modeling high density ratio incompressible interfacial flows,  American Society of Mechanical Engineers, Fluids Engineering Division (Publication) FED 257 (2002) 707–713. URL: https://doi.org/10.1115/FEDSM2002-31125. doi:10.1115/FEDSM2002-31125.
  * Spurk and Aksel [2020] J. H. Spurk, N. Aksel, Potential Flows, Springer International Publishing, Cham, 2020, pp. 345–442. URL: https://doi.org/10.1007/978-3-030-30259-7_10. doi:10.1007/978-3-030-30259-7_10.
  * Scheufler and Roenby [2021] H. Scheufler, J. Roenby,  Twophaseflow: An openfoam based framework for development of two phase flow solvers,  arXiv preprint (2021). URL: https://doi.org/10.48550/arXiv.2103.00870.
  * Desjardins and Moureau [2010] O. Desjardins, V. Moureau,  Methods for multiphase flows with high density ratio,  Center for Turbulence Research Proceedings of the Summer Program (2010) 313–322. URL: https://web.stanford.edu/group/ctr/Summer/SP10/6_02_desjardins.pdf.
  * Raessi and Pitsch [2012] M. Raessi, H. Pitsch,  Consistent mass and momentum transport for simulating incompressible interfacial flows with large density ratios using the level set method,  Computers and Fluids 63 (2012) 70–81. URL: http://dx.doi.org/10.1016/j.compfluid.2012.04.002. doi:10.1016/j.compfluid.2012.04.002.
  * Chenadec and Pitsch [2013] V. L. Chenadec, H. Pitsch,  A 3D unsplit Forward/Backward Volume-of-Fluid approach and coupling to the level set method,  J. Comput. Phys. 233 (2013) 10–33. URL: http://dx.doi.org/10.1016/j.jcp.2012.07.019. doi:10.1016/j.jcp.2012.07.019.
  * Vaudor et al. [2014] G. Vaudor, A. Berlemont, T. Ménard, M. Doring,  A consistent mass and momentum flux computation method using rudman-type technique with a clsvof solver,  in: Fluids Engineering Division Summer Meeting, volume 46230, American Society of Mechanical Engineers, 2014, p. V01CT23A012. URL: https://doi.org/10.1115/FEDSM2014-21802.
  * Crank and Nicolson [1947] J. Crank, P. Nicolson,  A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type,  Mathematical Proceedings of the Cambridge Philosophical Society 43 (1947) 50–67. doi:https://doi.org/10.1017/S0305004100023197.
  * Rubin and Khosla [1976] S. Rubin, P. Khosla,  Higher-order numerical solutions using cubic splines,  AIAA Journal 14 (1976) 851–858. doi:https://doi.org/10.2514/3.61427.
  * van Leer [1979] B. van Leer,  Towards the ultimate conservative difference scheme. v. a second-order sequel to godunov’s method,  Journal of Computational Physics 32 (1979) 101–136. URL: https://www.sciencedirect.com/science/article/pii/0021999179901451. doi:https://doi.org/10.1016/0021-9991(79)90145-1.
  * Leonard [1979] B. Leonard,  A stable and accurate convective modelling procedure based on quadratic upstream interpolation,  Computer Methods in Applied Mechanics and Engineering 19 (1979) 59–98. URL: https://www.sciencedirect.com/science/article/pii/0045782579900343. doi:https://doi.org/10.1016/0045-7825(79)90034-3.
  * Roe [1986] P. L. Roe,  Characteristic-based schemes for the euler equations,  Annual review of fluid mechanics 18 (1986) 337–365. doi:https://doi.org/10.1146/annurev.fl.18.010186.002005.
  * van Leer [1974] B. van Leer,  Towards the ultimate conservative difference scheme. ii. monotonicity and conservation combined in a second-order scheme,  Journal of Computational Physics 14 (1974) 361–370. URL: https://www.sciencedirect.com/science/article/pii/0021999174900199. doi:https://doi.org/10.1016/0021-9991(74)90019-9.
  * Adamson et al. [1997] A. W. Adamson, A. P. Gast, et al., Physical chemistry of surfaces, 6th Edition, volume 150, New York: Wiley, 1997. URL: https://www.wiley.com/en-us/Physical+Chemistry+of+Surfaces%2C+6th+Edition-p-9780471148739.
  * Couderc [2007] F. Couderc, Développement d’un code de calcul pour la simulation d’écoulements de fluides non miscibles. Application à la désintégration assistée d’un jet liquide par un courant gazeux., Ph.D. thesis, Ecole nationale superieure de l’aeronautique et de l’espace, 2007. URL: https://theses.hal.science/tel-00143709/document.
  * Zuzio and Estivalezes [2011] D. Zuzio, J. Estivalezes,  An efficient block parallel amr method for two phase interfacial flow simulations,  Computers & Fluids 44 (2011) 339–357. URL: https://www.sciencedirect.com/science/article/pii/S0045793011000429. doi:https://doi.org/10.1016/j.compfluid.2011.01.035.
  * Zuzio et al. [2018] D. Zuzio, J.-L. Estivalèzes, B. DiPierro,  An improved multiscale eulerian–lagrangian method for simulation of atomization process,  Computers & Fluids 176 (2018) 285–301. URL: https://www.sciencedirect.com/science/article/pii/S0045793016304017. doi:https://doi.org/10.1016/j.compfluid.2016.12.018.
  * Xavier et al. [2020] T. Xavier, D. Zuzio, M. Averseng, J.-L. Estivalezes,  Toward direct numerical simulation of high speed droplet impact,  Meccanica 55 (2020) 387–401. URL: https://link.springer.com/article/10.1007/s11012-019-00980-x.
  * Anjos et al. [2014] G. Anjos, N. Borhani, N. Mangiavacchi, J. Thome,  A 3d moving mesh finite element method for two-phase flows,  Journal of Computational Physics 270 (2014) 366–377. URL: https://www.sciencedirect.com/science/article/pii/S0021999114002551. doi:https://doi.org/10.1016/j.jcp.2014.03.067.
  * Bhaga and Weber [1981] D. Bhaga, M. E. Weber,  Bubbles in viscous liquids: shapes, wakes and velocities,  Journal of Fluid Mechanics 105 (1981) 61–85. doi:10.1017/S002211208100311X.
  * Hua and Lou [2007] J. Hua, J. Lou,  Numerical simulation of bubble rising in viscous liquid,  Journal of Computational Physics 222 (2007) 769–795. URL: https://linkinghub.elsevier.com/retrieve/pii/S0021999106003949. doi:10.1016/j.jcp.2006.08.008.
  * Thuillet [2018] S. Thuillet, Simulation multi-échelle de l’atomisation d’un jet liquide sous l’effet d’un écoulement gazeux transverse en présence d’une perturbation acoustique, Theses, UNIVERSITE DE TOULOUSE, 2018. URL: https://hal.archives-ouvertes.fr/tel-02023561.
  * Desclaux et al. [2020] A. Desclaux, S. Thuillet, D. Zuzio, J.-M. Senoner, D. Sebbane, V. Bodoc, P. Gajan,  Experimental and Numerical Characterization of a Liquid Jet Injected into Air Crossflow with Acoustic Forcing,  Flow, Turbulence and Combustion 105 (2020) 1087–1117. URL: https://hal.archives-ouvertes.fr/hal-03104296. doi:10.1007/s10494-020-00126-0.
  * Bodoc et al. [2020] V. Bodoc, A. Desclaux, P. Gajan, F. Simon, G. Illac,  Characterization of Confined Liquid Jet Injected into Oscillating Air Crossflow,  Flow, Turbulence and Combustion 104 (2020) 1–18. URL: https://hal.archives-ouvertes.fr/hal-02780211. doi:10.1007/s10494-019-00037-9.
  * Ltd [2023] O. Ltd, Openfoam-v2306, 2023. URL: https://develop.openfoam.com/Development/openfoam/-/tree/OpenFOAM-v2306.