# A benchmark for surface-tension-driven incompressible two-phase flows

Anja Lippert<EMAIL_ADDRESS>Tobias Tolle<EMAIL_ADDRESS>Aaron Dörr<EMAIL_ADDRESS>Robert Bosch GmbH, Corperate Research,
Robert-Bosch-Campus 1, 71272 Renningen, Germany  Tomislav Maric<EMAIL_ADDRESS>darmstadt.de  Mathematical Modeling and Analysis, Technische Universität
Darmstadt, Alarich-Weiss-Str. 10, 64287 Darmstadt, Germany

###### Abstract

The Volume-of-Fluid (VoF) method for simulating incompressible two-phase flows
is widespread in academic and commercial simulation software because of its
many advantages: a high degree of volume conservation, applicability to
unstructured domain discretization (relevant for engineering applications),
straightforward parallel implementation with the domain-decomposition and
message-passing approach (important for large-scale simulations), and
intrinsic handling of strong deformations and topological changes of the fluid
interface. However, stable and accurate handling of small-scale capillary
flows (dominated by surface tension forces) is still challenging for VoF
methods. With many different VoF methods making their way into commercial and
open-source software, it becomes increasingly important to compare them
quantitatively and directly. For this purpose, we propose a set of simulation
benchmarks and use them to directly compare VoF methods available in OpenFOAM,
Basilisk, and Ansys Fluent. We use Jupyter notebooks to document and process
the benchmark results, making a direct comparison of our results with other
two-phase simulation methods very straightforward. The publicly available
input data, secondary benchmark data, and post-processing Jupyter notebooks
can be re-used by any two-phase flow simulation method that discretizes two-
phase Navier-Stokes equations in a one-fluid formulation, which can save a
significant amount of person-hours.

###### keywords:

volume-of-fluid (VoF), benchmark , two-phase flow , capillary flows

[listing]style=Plaintop

##  1 Introduction

Immiscible multiphase flows on small scales become increasingly important not
only in academic research but also in industrial applications, e.g. in lab-on-
a-chip applications or for ensuring media tightness for corrosion prevention
in electrical components. Industrial multiphase flow processes generally occur
in geometrically complex flow domains making predictive and accurate numerical
methods in such domains essential. The Volume-of-Fluid (VoF) method is widely
used because it allows for geometrically complex flow domains, large
deformations and breakup/merging of the fluid interface, while ensuring a high
degree of volume conservation.

The VoF method dates back to 1981 [1] and since then a multitude of different
variants of the VoF method have emerged, distinguishable by different methods
used to approximate the fluid interface and advect it. Geometrical VoF
methods, although originally found to be computationally expensive, are
emerging nowadays in both efficient and accurate versions and finding their
use in commercial and open-source Computational Fluid Dynamics (CFD) software.
Since the geometrical VoF methods are substantially more accurate than
algebraic methods, they are the focus of this manuscript. For the un-split,
geometric transport on unstructured meshes alone numerous different methods
exist (cf. [2] for a recent review), showing very active ongoing research in
this field. Another crucial and actively researched aspect of the VoF method
is a force-balanced, accurate, and numerically stable discretization of the
surface tension force, reviewed recently in [3].

Comparing VoF methods directly and efficiently becomes increasingly relevant,
with many different VoF method implementations becoming available in open-
source and commercial software. Direct method comparison based on verification
and validation tests [4] is also relevant for developing a specific VoF method
- in that case, a comparison between an improved version and an existing
version is necessary. To facilitate VoF method comparison for small-scale
engineering problems, we adapt canonical verification and validation cases to
include fluid combinations and dimensions relevant to industrial applications.
This study captures the status quo of state-of-the-art VoF method
implementations and sheds light on open challenges. The automated post-
processing that relies on open-source software and a specification of an open
secondary-data format provides a basis for continuous benchmarking [5] for
two-phase flow simulation methods.

Six different VoF method implementations are benchmarked for advection and
hydrodynamic accuracy, namely: the interFoam and interIsoFoam OpenFOAM [6]
solvers, the interFlow solver from the OpenFoam sub-project TwoPhaseFlow [7,
8], as well as the Ansys Fluent [9] and Basilisk [10] software. All test cases
are resolved with three increasing mesh refinements to test for mesh
convergence, and the hydrodynamic cases are simulated with three different
realistic fluid/fluid pairings.

To enable future comparisons of these and other methods, the benchmarks
(results, setups and postprocessing scripts) are made publicly available. The
snapshot of the benchmark repository used for this publication is published at
the TUDatalib data repository at TUDarmstadt [11], and the research repository
is available at Bosch Research GitHub [12].

##  2 The Volume-of Fluid Method

###  2.1 Two-phase Navier-Stokes equations in the one-field formulation

$\Omega^{-}(t)$
---
$\Omega_{c}$
---
$\chi(\mathbf{x},t)=1$
---
$\chi(\mathbf{x},t)=0$
---
$\partial\Omega$
---
$\boldsymbol{n}_{\Sigma}$
---
$\Omega^{+}(t)$
---
$\Sigma(t)$
---
Figure 1: Two phase flows domain: the ellipsoidal heavy phase flow is
surrounded by the light phase flow.

The Volume of Fluid (VoF) method utilizes a one-field formulation of the
Navier-Stokes Equations (NSE) for two-phase flows.

The solution domain (cf. fig. 1) $\Omega\subseteq\mathbb{R}^{d},d=1,2,3$, is
separated into two evolving subdomains $\Omega^{\pm}(t)\subseteq\Omega$ each
occupied with a fluid phase ”$+$” or ”$-$” with respective constant densities
$\rho^{\pm}$ and viscosities $\mu^{\pm}$, separated by the fluid interface
$\Sigma(t):=\partial\Omega^{-}(t),\Sigma(t)\in\mathbb{R}^{d-1}$. The phase
indicator

$\displaystyle\chi\left(\mathbf{x},t\right)=\begin{cases}1&\text{if
}x\in\Omega^{-}(t),\\\ 0&\text{if }x\in\Omega^{+}(t)\end{cases}$ (1)

indicates fluid phases $\Omega^{\pm}(t)$. The phase-indicator defines the
single-field density and viscosity

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

with $\rho^{\pm},\mu^{\pm}$ as constant densities and viscosities of fluid
phases occupying $\Omega^{\pm}$, respectively.

Assuming the incompressibility of both fluid phases results in

$\nabla\cdot\mathbf{v}=0$ (3)

and eq. 2, together with the surface-tension force responsible for the
momentum transfer at $\Sigma(t)$, results in the single-field formulation of
the two-phase momentum conservation equation

$\partial_{t}(\rho\mathbf{v})+\nabla\cdot(\rho\mathbf{v}\mathbf{v})-\nabla\cdot(\mu\nabla\mathbf{v})=-\nabla
p+\rho\mathbf{g}+\mathbf{f}_{\Sigma}+\nabla\cdot(\mu(\nabla\mathbf{v})^{T}),$
(4)

where

$\mathbf{f}_{\Sigma}:=\sigma\kappa_{\sigma}\mathbf{n}_{\Sigma}\delta_{\Sigma}$
(5)

is the surface tension force with the surface tension coefficient $\sigma$,
mean curvature of the fluid interface $\Sigma(t)$ $\kappa_{\Sigma}$,
$\mathbf{n}_{\Sigma}$ the interface-normal (cf. fig. 1) and the interface
Dirac’s distribution $\delta_{\Sigma}$ that activates $\mathbf{f}_{\Sigma}$
only at $\Sigma(t)$. More information on the mathematical model is available
in [13].

For arbitrary initial and boundary conditions the exact solution to one-field
formulation of two-phase Navier-Stokes equations given by eqs. 3 and 4 cannot
be obtained. The sharpness of the phase indicator 1 further complicates
numerical solutions of eqs. 3 and 4 by algebraic discretization methods.

Equivalent to the unstructured Finite Volume Method [14, 15, 16, 17], the VoF
method decomposes $\Omega$ into a set of non-overlapping fixed-in-time control
volumes $\Omega_{c}$

$\Omega=\bigcup_{c\in
C}\Omega_{c},\quad\Omega_{c}\cap\Omega_{k}\in\mathbb{R}^{e},e<d,\forall c,k\in
C,c\neq k.$ (6)

For a fixed-in-time control volume $\Omega_{c}$, the volume fraction
$\alpha_{c}(t)$ is defined as

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

and it represents the fractional volume of the phase $\Omega^{-}(t)$ occupying
$\Omega_{c}$ at time $t$. Incompressible flow condition 3 is equivalent to
volume conservation in each fluid phase $\Omega^{\pm}$. Expressing the rate of
change of any of the phase-specific volumes

$\displaystyle\Omega_{c}^{-}(t):=\Omega^{-}(t)\cap\Omega_{c}=\int_{\Omega_{c}}\chi(\mathbf{x},t)\,dV,$
(8)
$\displaystyle\Omega_{c}^{+}(t):=\Omega^{+}(t)\cap\Omega_{c}=\int_{\Omega_{c}}(1-\chi(\mathbf{x},t))\,dV,$

say, $\Omega_{c}^{-}(t)$, using the Reynolds transport theorem (RTT), and
dividing by $|\Omega_{c}|$ leads to

$\partial_{t}\alpha_{c}(t)=-\frac{1}{|\Omega_{c}|}\int_{\partial\Omega_{c}}\chi(\mathbf{x},t)\mathbf{v}\cdot\mathbf{n}\,dS,$
(9)

an exact equation if

$\bigcup_{c\in C}(\partial\Omega_{c}\cap\partial\Omega)=\partial\Omega,$ (10)

which is frequently violated, e.g., if $\partial\Omega$ is non-linear and it
is discretized by boundaries of control-volumes that are linear.

###  2.2 Discretization

Since the goal of establishing a benchmark is to establish expectations in
terms of accuracy on different methods irrespective of their differences, we
only briefly cover the differences between VOF methods used here in direct
comparison. The reader is pointed to other sources of literature for further
details.

####  2.2.1 Volume fraction advection

The VoF method was developed into two different categories over the years, the
algebraic and geometric VoF methods, both start with integrating eq. 9 over a
time step $[t^{n},t^{n+1}]$, and a fixed polyhedral control volume, whose
boundary is bounded by possibly non-convex and non-planar faces
$\partial\Omega_{c}=\bigcup_{f\in F_{c}}S_{f}$, leading 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{x},t)\mathbf{v}\cdot\mathbf{n}\,dS\,dt,$
(11)

using the notation $\phi_{c}^{n}\equiv\phi_{c}(t=t^{n})$. The volume fraction
advection equation 11 is still exact.

##### Algebraic VoF methods

Unstructured algebraic VOF methods use

$\alpha_{f}:=\frac{1}{|S_{f}|}\int_{S_{f}}\chi(\mathbf{x},t)\,dS$ (12)

in eq. 11, and face-centered FV velocity approximation $\mathbf{v}_{f}$ to
obtain

$\displaystyle\alpha_{c}^{n+1}$
$\displaystyle=\alpha_{c}^{n}-\frac{1}{|\Omega_{c}|}\sum_{f\in
F_{c}}\int_{t^{n}}^{t^{n+1}}\mathbf{v}_{f}(t)\cdot\frac{\mathbf{S}_{f}}{|S_{f}|}\int_{S_{f}}\chi(\mathbf{x},t)\,dS\,dt,$
(13) $\displaystyle=\alpha_{c}^{n}-\frac{1}{|\Omega_{c}|}\sum_{f\in
F_{c}}\int_{t^{n}}^{t^{n+1}}F_{f}(t)\alpha_{f}(t)dt.$

The temporal integration quadrature used the phase-specific volumetric flux
determines the temporal accuracy. For example, the trapezoidal quadrature
results in a Crank-Nicolson discretization of eq. 9,

$\alpha_{c}^{n+1}-\alpha_{c}^{n}+0.5\frac{\Delta
t}{|\Omega_{c}|}(F_{f}^{n+1}\alpha_{f}^{n+1}+F_{f}^{n}\alpha_{f}^{n})=0.$ (14)

All algebraic methods rely on algebraic approximations of
$\alpha_{f}=\alpha_{f}(\alpha_{Of},\alpha_{Nf})$, which is less accurate than
the geometric approach. Still, algebraic approximation of $\alpha_{f}$ is
computationally much faster than the geometric integration of eq. 9 discussed
below. Furthermore, algebraic methods allow for an implicit discretization,
shown e.g. in eq. 14, leading to increased numerical stability [18]. Algebraic
methods suffer from dispersion errors. The phase-specific volumetric fluxes
$F_{f}\alpha_{f}$ must be limited to suppress these errors [19]. Some methods
extend eq. 9 with an anti-diffusion (compression) term in the attempt to
maintain sharpness of $\\{\alpha_{c}\\}_{c\in C}$ [19], others recover a
hyperbolic tangent in the interface-normal direction from
$\\{\alpha_{c}\\}_{c}$ [20].

##### Geometrical VoF methods

Geometrical VoF methods are called geometrical, because they discretize eq. 9
by geometrically approximating $\chi(\mathbf{x},t)$ on the r.h.s. of eq. 9,
and integrate eq. 9 in time using geometrical calculations. Unstructured
geometrical VoF methods also compute the r.h.s integral of eq. 9 over a finite
volume whose boundary $\partial_{c}$ consists of possibly non-planar faces
(surfaces), that are, however, bounded with edges (mesh edges)
$\partial_{c}\Omega_{c}:=\bigcup_{f\in F_{c}}S_{f}$.

We avoid repeating the details of approximate solutions of eq. 11 by different
VoF methods used in this benchmark, the reader is directed to [2] for a recent
review of unstructured geometrical VoF methods, [21, 22, 7] for the OpenFOAM
solvers implementing plicRDF-isoAdvector schemes, [23] for Basilisk, and [24]
for Ansys Fluent.

It suffices to know that every geometrical VoF method integrates the r.h.s. of
eq. 11 differently, but all methods result with a discrete equation of the
form

$\alpha_{c}^{n+1}=\alpha_{c}^{n}-\frac{1}{|\Omega_{c}|}\sum_{f\in
F_{c}}V_{f}^{\alpha},$ (15)

where $V_{f}^{\alpha}$ is the phase-specific volume fluxed through the face
$S_{f}\in\partial\Omega_{c}$ over the time step $[t^{n},t^{n+1}]$.

Volume conservation, monotonicity, stability, boundedness-preserving property
($\alpha_{c}\in[0,1]$ by definition 7), serial and parallel computational
efficiency (numerical consistency in parallel), absolute accuracy and
convergence order, ability to handle geometrically complex solution domains
$\Omega$, are all determined by the way a VoF method approximates
$V_{f}^{\alpha}$ [2]. Methods benchmarked here all approximate
$\chi(\mathbf{x},t)$ using a piecewise-linear approximation of
$\Sigma_{c}(t):=\Sigma\cap\Omega_{c}$, the part of the fluid interface in a
cell $\Omega_{c}$. This Piecewise Linear Interface Calculation, PLIC, (cf. [2]
for details) approximates $\Sigma_{c}(t)$ as a plane in each cell $\Omega_{c}$
intersected by the fluid interface $\Sigma(t)$.

Errors in $V_{f}^{\alpha}$ discretization, even if they are very small, may
lead to very large errors in $\alpha_{c}^{n+1}$. This becomes evident from eq.
15 because $|\Omega_{c}|\propto h^{d}$ ($d$ is the spatial dimension) divides
any residual geometrical error in $V_{f}^{\alpha}$. This error amplification
can lead to catastrophic failure in the approximation of the surface tension
force $\mathbf{f}_{\sigma}$ and is still challenging in some cases for
structured geometrical VoF methods, and even more so for unstructured VoF
methods.

####  2.2.2 Surface tension force

All VoF methods benchmarked in section 3 model the surface tension force using
eq. 5. For small-scale problems, surface tension force plays a crucial role,
making its correct discretization necessary to achieve accurate simulation
results. The surface tension force is discretized as a volumetric force
density in the vicinity of the interface using the Continuum Surface Force
(CSF) model [25]

$\displaystyle\mathbf{f}_{\Sigma,f}=\sigma\kappa_{\Sigma_{f}}(\nabla\alpha)_{f},$
(16)

averaged over $S_{f}$ at the face-centroid, where $(\nabla\alpha)_{f}$ denotes
the discrete gradient operator, and
$(\nabla\alpha)_{f}\approx\mathbf{n}_{\Sigma}\delta_{\Sigma}$ in eq. 5.

The main difference between different VoF methods lies in the approximation of
the mean curvature $\kappa$. The simplest VoF curvature model is

$\displaystyle\kappa_{\Sigma,f}=-\left[\nabla\cdot\left(\frac{\nabla\alpha}{|\nabla\alpha|}\right)\right]_{f}\,.$
(17)

which is obviously problematic to evaluate on finer meshes, since
$|\nabla\alpha|\to\delta_{\Sigma}\mathbf{n}_{\Sigma}$ when $h\to 0$, $h$ being
the discretization length. For structured, hexahedral meshes, Torrey et al.
[26] proposed a height-function based curvature approximation

$\displaystyle\kappa=\frac{h^{\prime\prime}}{\left(1+h^{\prime
2}\right)^{3/2}}\quad\text{(2D formulation)},$ (18)

which is used in Basilisk [10] and its predecessor Gerris [23]. While well-
suited for Cartesian meshes, this approach poses considerable challenges for
unstructured, polyhedral meshes, in terms of complexity and computational
costs [27].

Another approach for curvature approximation is to locally fit a paraboloid to
an interface cell and its neighbour cells. Curvature is then analytically
calculated using the paraboloid parameterization. Fitting is formulated as a
least squares problem to best match e.g. volume fractions [28] or centroids of
PLIC polygons [7].

The last method used by the plicRDF-isoAdvector method [22, 7] is to compute
the mean curvature from a Reconstructed Distance Function, RDF, [29]. In this
approach, a signed distance $\phi$ is reconstructed geometrically from the
PLIC elements. From $\phi$, the curvature is calculated analogously to eq. 17
by

$\kappa=-\nabla\cdot\left(\frac{\nabla\phi}{|\nabla\phi|}\right).$ (19)

###  2.3 Benchmarked methods

The six different solvers listed in table 1 are considered for the test cases
in section 3. Table 1 shows their main attributes, ranging from licensing to
discretization of phase advection and surface tension approximation. The
solvers are chosen based on their common usage in either industry or academia
or both. For all solvers we use the recommended settings and do not perform
any use-case-specific adaptation which might taint the results.

Solver |  Publisher |  license |  mesh |  $\alpha$-transport |  curvature appr.
---|---|---|---|---|---
interFoam |  OpenCFD Ltd |  Gnu GPL v3 |  unstructured |  algebraic VoF |  derivative (eq. 17)
interIsoFoam |  OpenCFD Ltd |  Gnu GPL v3 |  unstructured |  geom. VoF |  derivative (eq. 17)
TwoPhaseFlow |  - |  Gnu GPL v3 |  unstructured |  geom. VoF |  parabolic fit
TwoPhaseFlow |  - |  Gnu GPL v3 |  unstructured |  geom. VoF |  reconstr. distance funct. (eq. 19)
Fluent |  Ansys Inc. |  commercial |  unstructured |  geom. VoF |  derivative
Basilisk |  - |  Gnu GPL v3 |  structured |  geom. VoF |  height funct. (eq. 18)
Table 1: Overview of benchmarked solvers and some relevant properties.

All present solvers offer the additional capability of adaptive mesh
refinement, however, it was not applied in this study. In addition, except for
Basilisk, all tools allow polyhedral meshes. While beneficial for complex
geometries as present in industrial applications, this additional layer of
numerical complexity is not considered here. Hence, only Cartesian meshes are
used in the following section 3. Furthermore, regarding mesh quality Cartesian
meshes are optimal for the Finite Volume discretization used by the
benchmarked solvers [30]. For OpenFOAM-related cases OpenFOAM v2112 [6] is
used with `blockMesh` for grid generation. The OpenFoam sub-project
TwoPhaseFlow also uses `blockMesh`. Fluent was used with Ansys Meshing [31] to
generate meshes while the geometries are constructed and exported as STEP file
from either Gmsh [32] in 2D or FreeCAD [33] in 3D. Basilisk cases are build up
with the Basilisk library [10], the successor of Gerris [23].

##  3 Benchmark cases

The following cases serve to quantify each VoF-implementation separately as
well as compare the solvers described in Sec. 2.3.

###  3.1 Convection cases

The accuracy of the convective transport of the volume fraction field $\alpha$
is compared for three cases, namely 2D diagonal transport of a sphere, a
sphere in a 2D vortex flow and a sphere in a 3D vortex flow. All flows are
reversed halfway though the simulation so that an easy comparison between the
volume fraction field at the beginning and end is possible via subtraction.
The $L_{1}$ shape error is defined as

$e_{\text{shape}}=\sum_{i}^{N_{\text{cells}}}V_{i}\left|\alpha_{i}^{\text{end}}-\alpha_{i}^{\text{start}}\right|\,.$
(20)

Note that for the two-dimensional test cases $V_{i}$ refers to cell areas and
only in three dimensions to the cell volumes. The time discretization of
$\Delta t=1\cdot 10^{-5}$ s results in a CFL number below unity for all cases.
In addition, the two TwoPhaseFlow solvers are not tested separately for
convection since they build on interIsoFoam.

####  3.1.1 2D sphere in diagonal flow

To analyze grid-orientation dependence of transport, a disc of radius $R=0.25$
m is initialized at $(x_{0},y_{0})=(0.5,0.5)$ m and transported diagonally on
a Cartesian mesh in a rectangular domain of 3×2 m2 with the background
velocity $\mathbf{U}=\left(1,0.5\right)$ m/s. The flow is reversed after 2 s.
The domain is resolved with 90×60 cells and twice refined to capture the grid
dependence.

(a) L1-shape error for 2D sphere in reversed diagonal flow and different
resolutions

(b) Isocontour of drop with interFoam for res. 120 at: top $t=0$ s, middle
$t=2$ s, bottom $t=4$ s

(c) Isocontour of drop with interIsoFoam for res. 120 at: top $t=0$ s, middle
$t=2$ s, bottom $t=4$ s

Figure 2: Error visualization for 2D sphere in reversed diagonal flow for
different resolutions.

Figure 2(a) shows the absolute shape error after 4 s of transport. For all
codes the error overall reduces with increasing mesh resolution. interIsoFoam
and Basilisk exhibit the smallest error as expected due to their geometrical
transport. Surprisingly Fluent shows an higher error than expected, since a
geometric transport is chosen here as well. The errors span with 1 to 2 orders
of magnitude a much broader band than expected for such widely applied tools.
Visualizing the actual drop shape at the reversal point and at the end in
figs. 2(b) and 2(c) shows the actual ”meaning” behind the shape error, showing
a completely deformed interface contour. Determining any shape depending
forces, for instance surface tension, will result in a large error.

####  3.1.2 2D sphere in vortex flow

To test the accuracy of convection against a complex flow in two-dimensions, a
disc of radius $R=0.15$ m is initialized at $(x_{0},y_{0})=(0.5,0.75)$ m in a
domain of 1×1 m2. The velocity is prescribed as

$\displaystyle u$ $\displaystyle=\cos\left(\pi
t/t_{end}\right)\left(-2\sin^{2}(\pi x)\sin(\pi y)\cos(\pi y)\right)$ (21)
$\displaystyle v$ $\displaystyle=\cos\left(\pi
t/t_{end}\right)\left(2\sin^{2}(\pi y)\sin(\pi x)\cos(\pi x)\right)$ (22)

with an end time of $t_{end}=3$ s. The domain is resolved by 32×32 cells and
twice refined.

(a) $L_{1}$ shape error for 2D sphere in reversed vortex flow and different
resolutions.

(b) Isocontour of drop with interFoam for res. 64 at: top $t=0$ s, middle
$t=1.5$ s, bottom $t=3$ s

(c) Isocontour of drop with interIsoFoam for res. 64 at: top $t=0$ s, middle
$t=1.5$ s, bottom $t=3$ s

Figure 3: Error visualization for 2D sphere in reversed vortex flow for
different resolutions.

Figure 3(a) shows the absolute shape error. We see a similar error
distribution for different resolutions and solvers as before. Basilisk and
interIsoFoam exhibit the smallest error, while interFoam with its algebraic
transport accumulates the largest errors. The visualization of this error is
shown on the right in figs. 3(b) and 3(c). The strong contour deformation
exposes the reason for the shape error. In this case, interFoam achieves lower
$L_{1}$ errors as for the diagonal flow. This can be explained by the fact
that for the vortex flow the velocity aligns with normals of at least some
cell faces. This is favorable for the algebraic transport. However, for the
diagonal flow, there is no such alignment. The improved performance of Fluent
for the vortex flow is difficult to explain without further details on the
interface transport.

####  3.1.3 3D sphere in vortex flow

Extending the test of the previous paragraph to three dimensions, we
initialize a sphere of radius $R=0.15$ m at
$(x_{0},y_{0},z_{0})=(0.35,0.35,0.35)$ m in a domain of 1×1×1 m3. The velocity
is given by

$\displaystyle u$ $\displaystyle=2\sin^{2}(\pi x)\sin(2\pi y)\sin(2\pi
z)\cos\left(\pi t/t_{end}\right)$ (23) $\displaystyle v$
$\displaystyle=-\sin^{2}(\pi y)\sin(2\pi x)\sin(2\pi z)\cos\left(\pi
t/t_{end}\right)$ (24) $\displaystyle w$ $\displaystyle=-\sin^{2}(\pi
z)\sin(2\pi x)\sin(2\pi y)\cos\left(\pi t/t_{end}\right)\,.$ (25)

The resolution starts at 64×64×64 with two levels of refinement and an end
time of $t_{end}=3$ s.

(a) $L_{1}$ shape error for 3D sphere in reversed vortex flow and different
resolutions.

(b) Isocontour of drop with interFoam for res. 256 at: top $t=0$ s, middle
$t=1.5$ s, bottom $t=3$ s.

(c) Isocontour of drop with interIsoFoam for res. 256 at: top $t=0$ s, middle
$t=1.5$ s, bottom $t=3$ s.

Figure 4: $L_{1}$ shape error for 3D sphere in reversed vortex flow for
different resolutions.

Figure 4(a) shows that for a similar sized domain and drop, the error
distribution between the solvers and resolutions stays consistent. Only Fluent
appears to improve in three dimensions but still the error remains between the
algebraic advection solver interFoam and the other geometric transport solvers
Basilisk and interIsoFoam. Furthermore, the errors for all resolutions move
closer together. Considering that the accuracy of algebraic transport depends
on the alignment between flow direction and cell face normals, this is not
surprising as this alignment depends on cell topology not on spatial
resolution. However, looking at figs. 4(b) and 4(c), the deviation from the
initial sphere shape at the end of the simulation is clearly visible with
interFoam, even for the finest resolution of 256 cells.

###  3.2 Hydrodynamic cases

The hydrodynamic cases quantify the numerical accuracy of the physical
interaction between convection, diffusion and balancing surface tension
effects with pressure difference. Three fluid pairings are considered. Their
fluid properties can be found in table 2. These span a density ratio of 1.6 to
838.8 and a kinematic viscosity ratio of 0.06 to 15.3. A total of four cases
is investigated, some of them in 2D and 3D, namely stationary droplet,
translating droplet, oscillating droplet and oscillating wave.

@20 °C |  density $\rho$
( kg/m3) |  kin. viscosity $\nu$
( m2/s) |  dyn. viscosity $\mu$
( Pa s) |  surface tension $\sigma$
( N/m)
---|---|---|---|---
water[34, 35] |  998.2 |  1e-6 |  9.982e-4 |  72.74e-3[36]
air[35] |  1.19 |  1.53e-5 |  1.8207e-5
oil (Ravenol CLP 220 [37]) |  888 |  2.71e-4111extrapolated value |  0.240648 |  32.9e-3222value for NACA reference oil from [38]
air[35] |  1.19 |  1.53e-5 |  1.8207e-5
oil (Novec 7500 [39]) |  1614 |  0.77e-6 |  1.24278e-3 |  49.5e-3[40]
water[34, 35] |  998.2 |  1e-6 |  9.982e-4
Table 2: Fluid pairings under consideration for hydrodynamic test cases.

The time step for each case is restricted by the CFL number as well as the
capillary time step constraint

$\Delta t_{\sigma}=\sqrt{\frac{(\rho^{-}+\rho^{+})\Delta
x^{3}}{2\pi\sigma}}\,,$ (26)

as found in [41] for the respective resolutions $\Delta x$.

For the stationary and translating droplet test cases (both 2D and 3D), the
following error metrics considering velocity are examined:

$\displaystyle L_{1}(\mathbf{v})=$
$\displaystyle\frac{1}{N_{\text{cells}}}\sum_{i}^{N_{\text{cells}}}|\mathbf{v}_{i}-\mathbf{v}_{\text{ref}}|,$
(27) $\displaystyle L_{\infty}(\mathbf{v})=$
$\displaystyle\max(|\mathbf{v}_{i}-\mathbf{v}_{\text{ref}}|),$ (28)

with $\mathbf{v}_{\text{ref}}$ denoting the uniform background velocity, being
zero for the stationary droplet.

###  3.3 Volume fraction initialization

Figure 5: Influence of different initialization algorithms for the 2D
oscillating wave for three different resolutions, from left to right: 32, 64,
128. The solver interIsoFoam has been used for these computations with the
water/air fluid pairing.

Before continuing the discussion on hydrodynamic benchmark results, we want to
highlight the influence of initialization accuracy here. In the previous
section 3.1 this was less of an issue since the flow fields are all reversed
and the initial volume fraction field also serves as the reference. So while
the initial drop shape may differ between solvers, the error-measurement only
compares their drop shape to itself. For the hydrodynamic cases, this is not
valid anymore. Observe the solution quality difference, especially for coarse
resolutions, in fig. 5 which shows two options available in OpenFOAM. With the
”inside-outside” approach, the volume fraction value of a cell is determined
solely by evaluating eq. 1 for its centroid $\mathbf{x}_{c}$. This means
$\alpha_{c}$ assumes values of either 1 or 0, depending on whether
$\mathbf{x}_{c}\in\Omega^{+}$ or $\mathbf{x}_{c}\in\Omega^{-}$. Consequently,
this approach introduces rather large errors in the volume fraction. The more
accurate approach, as used in setAlphaField, determines the actual overlap
between initialized shape contour and grid cells and by doing so calculates
the volume fraction. The error decreases with increasing resolution, but must
be kept in mind when comparing different solvers. OpenFOAM, TwoPhaseFlow and
Basilisk all offer similar initialization algorithms. For Fluent, there are
different options to define the phase distribution (here, custom field
functions are used). However, phase initialization is based on the ”inside-
out” method only, albeit with a default ”smoothing” step controlled by a
relaxation parameter. The smoothing step adds one or more diffusion steps to
the starting form of the ”inside-out” volume fraction field [24]. Throughout
this work, the default of one smoothing step with a factor of 0.25 is used.

####  3.3.1 2D stationary droplet

(a) water/air, res: 64

(b) Ravenol/air, res: 64

(c) Novec/water, res: 64

Figure 6: Mean velocity error for 2D stationary droplet for different fluid
pairings.

(a) water/air, res: 64

(b) Ravenol/air, res: 64

(c) Novec/water, res: 64

Figure 7: Max velocity error for 2D stationary droplet for different fluid
pairings.

(a) water/air

(b) Ravenol/air

(c) Novec/water

Figure 8: Mean (eq. 27, dark color) and maximum (eq. 28, light color) velocity
errors over time $t>\frac{1}{2}t_{\text{max}}$ from the 2D stationary droplet
summerized as violin plots. For each fluid pairing, the results of different
solvers and resolutions are displayed.

To analyze so-called parasitic currents in 2D, a drop with $R=10^{-3}$ m is
initialized at the center of a Cartesian mesh of 32×32 cells in a rectangular
domain of 0.01×0.01 m2. The end time is set to $t_{\text{end}}=0.1$ s. Two
additional levels of refinement are investigated. The error evaluation is done
using eqs. 27 and 28. Plotting the mean and max error over time for the middle
resolution for all fluid pairings reveals the differing behavior of solvers
over time. In general, Basilisk and interFlow, using plicRDF for curvature
(abbr. interFlow-RDF), produce the smallest error for all fluid pairings. As
expected, the error magnitude decreases across-the-board for higher viscous
fluids, see Ravenol/air vs. water/air. The fluid pairing Novec/water with a
density and viscosity ratio of close to one reduces the error compared to
water/air by about one order of magnitude. Furthermore, the max error (fig. 7)
exceeds the mean error (fig. 6) by at least one order of magnitude. In
addition, the plots reveal that only for Basilisk and interFlow-RDF the error
reduces over time. None of the others are unstable, but reach some kind of
plateau state.

To capture the essence of the simulations over all resolutions, fluid pairings
and their time-dependent changes, we introduce in fig. 8 violin plots which
include for each solver and resolution the mean (darker color) and max error
(lighter color) with their respective max and min values as well as their
time-variant distribution. For these distributions, only the values from the
second half of simulated time, $t>\frac{1}{2}t_{\text{max}}$, are considered
to reduce the influence of initialization, shown in section 3.3. This overview
reveals a multitude of things, the most obvious is a missing trend of a
reduction in error with increased resolution, meaning the error for a
resolution 128×128 might be higher than for 64×64. This stands in strong
contrast to the previous findings regarding the convection error in section
3.1.

####  3.3.2 3D stationary droplet

Extending the above described setup to three dimension in a cubic domain of
0.01×0.01×0.01 m3 and a droplet radius of $R=10^{-3}$ m, the deviation from
the initialized drop at rest is studied. The drop is initialized at the center
of the domain. Errors are once more evaluated via eqs. 27 and 28. The end time
is set to $t_{\text{end}}=0.1$ s with the coarsest resolution being 32×32×32
followed by two refinement steps. Figure 9 shows again the missing correlation
of error magnitude and resolution. For instance, considering the data of
Fluent, for water/air the mean and max error increase for increased
resolution. In a Ravenol/air setup, however, the mean error decreases while
the max error still increases. For interIsoFoam, the trend is completely
reversed, decreasing errors for water-air, but increasing errors for Ravenol-
air. In general, an increasing error for finer refinement appears counter-
intuitive but can be explained by the less diffused interface and, hence,
higher more localized error. This primarily hits the geometrical solvers who
base their curvature on the second derivative of the $\alpha$-field, such as
Fluent and interIsoFoam. Basilisk’s poor performance for Ravenol-air for
128×128×128 cannot be explained by this, but should be noted here, as it
stands in harsh contrast to the 2D case in section 3.3.1. Overall, the errors
are close in magnitude, with Basilisk still outperforming the unstructured
solvers, but less distinct than in two dimensions (see fig. 8 for reference).
In general, compared to the 2D case, the mean and max errors differ more
strongly, covering at least one order of magnitude, often even two.

(a) water/air

(b) Ravenol/air

(c) Novec/water

Figure 9: Mean (eq. 27, dark color) and maximum (eq. 28, light color) velocity
errors over time $t>\frac{1}{2}t_{\text{max}}$ from the 3D stationary droplet
summerized as violin plots. For each fluid pairing, the results of different
solvers and resolutions are displayed.

####  3.3.3 2D translating droplet

The following case upgrades the 2D stationary droplet case, see section 3.3.1,
with a background velocity field of $\mathbf{U}=(1.0,0)$ m/s to evaluate their
error superposed with the advection errors investigated in section 3.1. For
this purpose, we initialize a drop with $R=10^{-3}$ m at $(2R,5R)$ in a
rectangular domain of 0.02×0.01 m2 and simulate up to an end time of
$t_{\text{end}}=0.015$ s. The domain is resolved with 64×32 and subsequently
twice refined for this study.

(a) water/air

(b) Ravenol/air

(c) Novec/water

Figure 10: Mean (eq. 27, dark color) and maximum (eq. 28, light color)
velocity errors over time $t>\frac{1}{2}t_{\text{max}}$ from the 2D
translating droplet summerized as violin plots. For each fluid pairing, the
results of different solvers and resolutions are displayed.

Figure 10 summarizes the different solver performances for varied resolution
and fluid pairings. The evaluation is done using eq. 27 and eq. 28 with
$\mathbf{v}_{\text{ref}}=\mathbf{U}$. Focusing first on the right hand side of
fig. 10, the OpenFOAM based solvers interFlow, interFoam and interIsoFoam,
perform more similarly compared to Fluent and Basilisk with their max and mean
error being of a comparable order. Only interFoam and interIsoFoam exhibit
outliers for the max error when it comes to the fluid pairing of Novec and
water, even though the density ratio is nearly one here. For all fluid-
pairings and resolutions, Basilisk shows the best results. Fluent stays
between Basilisk and the OpenFOAM solvers. This is surprising based on the 2D
convection and 2D stationary droplet case, where both errors exceeded for
instance interFlow-RDF by one order of magnitude. Due to the closed source and
commercial nature of Fluent, we can only assume that the errors either cancel
each other for this specific case or some additional unknown feature is at
work.

####  3.3.4 3D translating droplet

Switching to a cuboid domain of 0.02×0.01×0.01 m3, the drop is initialized
with $R=10^{-3}$ m at $(2R,5R,5R)$. The translation velocity is as in the 2D
case chosen as $\mathbf{U}=(1.0,0,0)$ m/s while translating to an end time of
$t_{\text{end}}=0.015$ s. The domain is resolved with 64×32×32 and twice
refined to study the influence of grid resolution.

(a) water/air

(b) Ravenol/air

(c) Novec/water

Figure 11: Mean (eq. 27, dark color) and maximum (eq. 28, light color)
velocity errors over time $t>t_{\text{max}}$ from the 3D translating droplet
summerized as violin plots. For each fluid pairing, the results of different
solvers and resolutions are displayed.

Comparing first the mean errors for different fluid pairings in 3D as visible
in fig. 11 with the 2D errors in fig. 10: The errors stay in a similar range,
even slightly improving in three dimensions. Especially Fluent and interFlow-
RDF improve. Across the board, the solvers perform worst for water/air as
expected for low viscosities, highest density ratio and largest surface
tension coefficient. It is also worth mentioning that while the mean error
decreases for most solvers with increasing resolution, the maximum error often
increases. This effect was less noticeable for the stationary droplet 3D, see
fig. 9 and must stem from the additional phase transport. Different
phase/interface positions for each time step might result in disadvantageous
normal calculations for the geometrical advection schemes which introduce
additional errors. However, even while max errors increase, the mean errors
stay in a similar range to the 3D stationary droplet or partly improve.

####  3.3.5 2D oscillating wave

The capillary wave with an initial interface position
$h(x,0)=a\cos{(2\pi/L(x-L/2))}+L/2$ with $a=5.5\cdot 10^{-5}$ m and
$L=2.6\cdot 10^{-3}$ m is initialized in domain of L×L. Gravitational effects
are not considered. The domain is resolved 32×32 and refined twice. The
considered physical time span is such that 5 oscillations are tracked. The
interface height at the center of the domain is tracked and compared to the
analytical solution by Prosperetti (eqs. (24) and (25) in [42]). Since the
high viscosity of Ravenol causes overdamping, see table 3 (frequency close to
zero), only the fluid pairings water/air and Novec/water are discussed. To
summarize the temporal evolution of the difference between the numerically
predicted interface height and the asymptotic reference solution, the time-
averaged relative amplitude error

$\bar{\epsilon}_{a,rel}=\frac{\Delta
t}{\tau}\sum_{i}\frac{|a_{i,\text{num}}-a_{i,\text{ref}}|}{a_{0}}$ (29)

is evaluated. Here, $\tau$ denotes the total physical time simulated, the
subscript $i$ refers to the $i$-th time step, $a_{i,\text{num}}$ and
$a_{i,\text{ref}}$ refer to the numerically approximated interface height and
the the interface height according to the reference at time step $i$ and
$a_{0}=5.5\cdot 10^{-5}$ m denotes the initial amplitude.

fluid pairing | decay rate $\delta$ [1/s] | frequency $f$ [1/s]
---|---|---
water/air | 11.46 | 161.13
Novec/water | 18.90 | 79.53
Ravenol/air | 181.05 | 0.001
Table 3: Parameters of the asymptotic solution $h(0,t)-L/2=a\exp(-\delta
t)\cos(2\pi ft)$ to the capillary wave problem for $t\to\infty$ by Prosperetti
(eqs. (24) and (25) in [42]). The remaining parameters are given by
$L=2.6\cdot 10^{-3}$ m and $a=5.5\cdot 10^{-5}$ m, respectively.

(a) water/air, res: 128

(b) Novec/water, res: 128

Figure 12: Tracking the film height at center over time for oscillating wave.
The reference is given by the analytical solution by [42].

(a) water/air

(b) Novec/water

Figure 13: Time averaged relative difference between simulated wave height and
analytical solution for a 2D oscillating wave.

Figure 12 displays the temporal evolution of the film height above the center
for water/air and Novec/water at the finest resolution. Deviations are visible
for interIsoFoam and Fluent. With interIsoFoam as solver, the deacy of the
amplitude is slightly underestimated for water/air, a behaviour that is more
pronounced for the Novec/water fluid pairing. Results from Fluent show a
longer oscillation period and a faster amplitude decay than given by the
reference solution for both fluid pairings. The behaviour of Fluent can be
explained, at least partially, by its rather inaccurate initialization
approach. As shown in section 3.3, initialization errors impact the computed
solution since they represent deviations from the exact interface. All other
solvers coincide with the reference solution for both fluid pairings. These
observations are also reflected in the errors evaluated according to eq. 29,
displayed in fig. 13. Again, interIsoFoam and Fluent display the largest
errors for both fluid pairings, with the exception of water/air at the
intermediate resolution. For that configuration, interIsoFoam yields the
lowest error among the benchmarked solvers. Given that the finest resolution
underestimates the decay of the amplitude and the coarse resolution tends to
overestimate the decay, interIsoFoam probably obtains the best agreement with
the reference for the intermediate resolution. Besides interIsoFoam, all
solvers show converging errors. Comparing fluid pairings, the errors are in
the same order of magnitude for fixed resolution and solver. A remark on the
choice of boundary condition for the volume fraction field when using the
interFlow solver: with OpenFOAM’s `zeroGradient` boundary condition, in
principle a valid boundary condition for this case which has been used for
interFoam and interIsoFoam, the results showed no resemblance with the
reference solution. However, when the `constantAlphaContactAngle` boundary
condition, provided by TwoPhaseFlow [7], with a conatct angle of 90 degree,
accuracy improved distinctly as shown in figs. 12 and 13.

####  3.3.6 3D oscillating droplet

As a test case in three dimensions with deforming interface an oscillating
droplet is chosen. The droplet is initialized as a rotational ellipsoid

$\frac{x^{2}}{a_{0}^{2}}+\frac{y^{2}+z^{2}}{b_{0}^{2}}=1$ (30)

with semi-axes $a_{0}=5.25\cdot 10^{-4}$ m and $b_{0}=4.8795\cdot 10^{-4}$ m,
respectively. The droplet’s volume is equivalent to that of a spherical
droplet of radius $R=0.5\cdot 10^{-3}$ m. The droplet is placed at the center
of a cubic domain of edge length $l=5\cdot 10^{-3}$ m and the domain’s center
coincides with the origin. The domain is resolved with 50 cells in each
direction and twice refined. Slip is prescribed as velocity boundary condition
and a Neumann boundary condition for the pressure. Both fields are set to zero
initially. A total physical time of $t=0.02$ s is simulated with a fixed time
step determined by eq. 26. For a droplet that is only slightly perturbed from
a spherical shape, an analytical solution for the eigenshape at eigenfrequency
$\omega_{n}$ has been given by Lamb [43]. Rather than performing a surface
harmonics expansion of the droplet and evaluating it over time to find the
interface evolution, the simplified approach by Shin and Juric [44] is
followed. Thus, in order to compute the oscillation of the longest semi axis
$a(t)$, the amplitude of the lowest-order eigenshape ($n=2$) is set to the
initial semi axis length $a_{0}$ minus the mean droplet radius $R$. This leads
to the reference solution

$a(t)=R+\left(a_{0}-R\right)\exp\\!{\left(-\frac{(n-1)(2n+1)\nu}{R^{2}}t\right)}\cos(\omega_{n}t)~{}\text{for}~{}n=2,$
(31)

where

$\omega_{n}^{2}=\frac{n(n+1)(n-1)(n+2)\sigma}{\left((n+1)\rho^{-}+n\rho^{+}\right)R_{0}^{3}}.$
(32)

Since eq. 31, by assumption of a viscous droplet surrounded by a fluid of
comparatively negligible kinematic viscosity, is only valid for the fluid
paring water/air, the kinematic viscosity $\nu$ is set to the value for water
(cf. table 2), i.e. $\nu=\nu^{-}$.

Figure 14 shows the temporal evolution of the initially longest semi-axis for
water/air and Novec/water fluid pairings. Ravenol/air is not reported here as
the high viscosity of the oil prevents oscillatory behaviour similar to the
oscillating wave. Note also that the reference solution is only plotted in
fig. 14a and not b since Novec/water does not fulfill the criterion of a
negligible ambient fluid viscosity. With water/air as fluid pairing, all
solvers match the oscillation period quite well. After $t=0.015$ s, small
offsets become apparent, for interIsoFoam and interFlow-RDF offsets are
noticeable at earlier times already. Based on their previous behavior in other
test scenarios, this is not explainable. Figure 15 highlights this discrepancy
and even shows a increasing deviation with increasing resolution. For Fluent,
the influence of volume fraction initialization is visible for around one and
a half oscillation periods. Afterwards its performance stabilizes and it
performs rather well. This finding agrees with the significant improvement of
Fluent’s time averaged relative amplitude error with increasing resolution as
visible in fig. 15, again pointing towards the drawback of inaccurate
initialization. In general, all solvers overestimate the damping, probably due
to realistic fluid pairing where the outer viscosity differs from zero.
Regarding the fluid pairing Novec/water only the relative solver performance
can be judged. Again, Fluent struggles with initialization which introduces
additional harmonics. The dominant oscillating frequency lies in a similar
range for all solvers, with interFoam and interIsoFoam as outliers. InterFlow-
RDF overshoots massively towards one side. Assuming a good performance of
Basilisk based on the previous results, only interFlow-fP shows a comparable
behavior.

(a) Water/air, res: 100

(b) Novec/water, res: 100

Figure 14: Tracking the initially longest semi-axis in 3D over time for the
oscillating droplet. Figure 15: Time averaged relative difference between
simulated initially longest semi-axis and the analytical solution for 3D
oscillating drop for water/air.

##  4 Conclusions

We provide a basis for continuously benchmarking numerical methods for
surface-tension-driven incompressible two-phase flows that discretize the
single-field formulation of two-phase Navier-Stokes equations. The post-
processing implemented in Jupyter notebooks can be re-used without
modification by methods developed in other software to compare with our
benchmark results directly. Providing the benchmark suite as a reusable open-
source project can save significant research time.

Overall, the benchmark results demonstrate a need for significant improvements
in the approximation of surface tension forces, especially for the
unstructured VoF methods available in OpenFOAM and Fluent. The structured VoF
method in Basilisk delivers the best results with the height-function
curvature approximation, but not for all fluid pairings, and with higher
accuracy in 2D compared to 3D.

This study also highlights the importance of an accurate field initialization
to enable an acceptable solver performance.

In future work, we will re-use the benchmark suite to continuously compare new
approximation methods for surface tension forces with other methods. We invite
the multiphase community to participate in this direct comparison, following
the instructions in the linked code repository [12].

##  5 Acknowledgments

The last author acknowledges the funding received from the German Research
Foundation (DFG)- Project-ID 265191195 - SFB 1194.

## References

  * Hirt and Nichols [1981] C. W. Hirt, B. D. Nichols,  Volume of fluid (vof) method for the dynamics of free boundaries,  Journal of Computational Physics 39 (1981) 201–225.
  * 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.
  * Popinet [2018] S. Popinet,  Numerical models of surface tension,  Annual Review of Fluid Mechanics 50 (2018) 49–75.
  * Oberkampf and Roy [2010] W. L. Oberkampf, C. J. Roy, Verification and validation in scientific computing, Cambridge University Press, 2010.
  * Grambow et al. [2019] M. Grambow, F. Lehmann, D. Bermbach,  Continuous Benchmarking: Using System Benchmarking in Build Pipelines,  in: 2019 IEEE International Conference on Cloud Engineering (IC2E), IEEE, Prague, Czech Republic, 2019, pp. 241–246. URL: https://ieeexplore.ieee.org/document/8790186/. doi:10.1109/IC2E.2019.00039.
  * OpenCFD Ldt. [2022] OpenCFD Ldt., OpenFOAM Repository, https://develop.openfoam.com/Development/openfoam/-/tree/OpenFOAM-v2112, 2022\. Git tag OpenFOAM-v2112.
  * Scheufler and Roenby [2021] H. Scheufler, J. Roenby,  TwoPhaseFlow: An OpenFOAM based framework for development of two phase flow solvers,  arXiv preprint arXiv:2103.00870 (2021).
  * Two [2022] TwoPhaseFlow Repository, version of2112, https://github.com/DLR-RY/TwoPhaseFlow/tree/of2112, 2022\. Git commit: cbaf314577987db34a9e41bb4bdca511c40d080.
  * ANSYS, Inc. [2020] ANSYS, Inc., Ansys Fluent User Guide 2020 R1, 2020.
  * S. Popinet and collaborators [2015] S. Popinet and collaborators, Basilisk, http://basilisk.fr, since 2015. Accessed: 2022-10-13.
  * Lippert et al. [1 28] A. Lippert, T. Tolle, A. Dörr, T. Maric, A benchmark for surface-tension-driven incompressible two-phase flows - benchmark data, 2022-11-28. URL: https://tudatalib.ulb.tu-darmstadt.de/handle/tudatalib/3630. doi:10.48328/tudatalib-989.
  * Robert Bosch GmbH [2022] Robert Bosch GmbH, BoschResearch/sepMutliphaseFoam, git tree/publications/ST-VoF-benchmark, https://github.com/boschresearch/sepMultiphaseFoam/tree/publications/ST-VoF-benchmark, 2022\. Created: 2022-11-20.
  * Tryggvason et al. [2011] G. Tryggvason, R. Scardovelli, S. Zaleski, Direct numerical simulations of gas–liquid multiphase flows, Cambridge university press, 2011.
  * Jasak [1996] H. Jasak, Error analysis and estimation for the finite volume method with applications to fluid flows., Ph.D. thesis, Imperial College London (University of London), 1996.
  * Hirsch [2007] C. Hirsch, Numerical computation of internal and external flows: The fundamentals of computational fluid dynamics, Elsevier, 2007.
  * Maric et al. [2014] T. Maric, J. Hopken, K. Mooney, The OpenFOAM technology primer, Sourceflux, 2014.
  * Moukalled et al. [2016] F. Moukalled, L. Mangani, M. Darwish, The finite volume method, Springer, 2016.
  * Denner et al. [2022] F. Denner, F. Evrard, B. van Wachem,  Breaching the capillary time-step constraint using a coupled vof method with implicit surface tension,  Journal of Computational Physics 459 (2022) 111128.
  * Deshpande et al. [2012] S. Deshpande, L. Anumolu, M. Trujillo,  Evaluating the performance of the two-phase flow solver interfoam,  Computational Science & Discovery 5 (2012). doi:10.1088/1749-4699/5/1/014016.
  * Chen et al. [2022] D. Chen, B. Xie, F. Xiao,  Revisit to the thinc/qq scheme: Recent progress to improve accuracy and robustness,  International Journal for Numerical Methods in Fluids 94 (2022) 719–755.
  * Roenby et al. [2016] J. Roenby, H. Bredmose, H. Jasak,  A Computational Method for Sharp Interface Advection,  Royal Society Open Science 3 (2016). URL: http://arxiv.org/abs/1601.05392. doi:10.1098/rsos.160405.
  * 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://www.sciencedirect.com/science/article/pii/S0021999119300269. doi:10.1016/j.jcp.2019.01.009, arXiv: 1801.05382 Publisher: Academic Press.
  * Popinet [2009] S. Popinet,  An accurate adaptive solver for surface-tension-driven interfacial flows,  Journal of Computational Physics 228 (2009) 5838–5866. doi:doi.org/10.1016/j.jcp.2009.04.042.
  * ANSYS, Inc. [2020] ANSYS, Inc., Ansys Fluent Theory Guide 2020 R1, 2020.
  * 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.
  * Torrey et al. [1985] M. D. Torrey, L. D. Cloutman, R. C. Mjolsness, C. W. Hirt, NASA-VOF2D: a computer program for incompressible flows with free surfaces, Technical Report, Los Alamos National Lab., NM (USA), 1985. URL: https://www.osti.gov/biblio/5934123.
  * Ivey and Moin [2015] C. B. Ivey, P. Moin,  Accurate interface normal and curvature estimates on three-dimensional unstructured non-convex polyhedral meshes,  Journal of Computational Physics 300 (2015) 365–386. doi:doi.org/10.1016/j.jcp.2015.07.055.
  * Jibben et al. [2019] Z. Jibben, N. Carlson, M. Francois,  A paraboloid fitting technique for calculating curvature from piecewise-linear interface reconstructions on 3d unstructured meshes,  Computers & Mathematics with Applications 78 (2019) 643–653. doi:10.1016/j.camwa.2018.09.009, proceedings of the Eight International Conference on Numerical Methods for Multi-Material Fluid Flows (MULTIMAT 2017).
  * Cummins et al. [2005] S. J. Cummins, M. M. Francois, D. B. Kothe,  Estimating curvature from volume fractions,  Computers & Structures 83 (2005) 425–434. doi:10.1016/j.compstruc.2004.08.017, frontier of Multi-Phase Flow Analysis and Fluid-Structure.
  * Juretić [2005] F. Juretić, Error analysis in finite volume CFD, Ph.D. thesis, Imperial College London (University of London), 2005\.
  * ANSYS, Inc. [2020] ANSYS, Inc., Ansys Meshing User Guide 2020 R1, 2020.
  * Christophe Geuzaine and Jean-Francois Remacle [2020] Christophe Geuzaine and Jean-Francois Remacle, Gmsh, 2020. URL: http://http://gmsh.info/.
  * Riegel et al. [2022] J. Riegel, W. Mayer, Y. van Havre, FreeCAD (Version 0.19), 2001-2022. URL: http://www.freecadweb.org.
  * Lemmon et al. [2022] E. W. Lemmon, I. H. Bell, M. L. Huber, M. O. McLinden, NIST Chemistry WebBook, NIST Standard Reference Database Number 69, National Institute of Standards and Technology, Gaithersburg MD, 2022, p. 20899\. URL: https://webbook.nist.gov.
  * VDI Gesellschaft [2010] VDI Gesellschaft (Ed.), VDI Heat Atlas, VDI-Buch, Springer Berlin Heidelberg, 2010\. URL: https://books.google.de/books?id=6benoAEACAAJ.
  * International Association for the Properties of Water and Steam [2014] International Association for the Properties of Water and Steam, IAPWS R1-76(2014) Revised release on surface tension of ordinary water substance, 2014. URL: http://www.iapws.org/relguide/Surf-H2O-2014.pdf.
  * Ravensberger Schmierstoffvertrieb GmbH [2022] Ravensberger Schmierstoffvertrieb GmbH, RAVENOL Getriebeöl CLP 220, 2022. URL: https://www.ravenol.de/storage/app/media/product-pdf/Tds_1332109_en.pdf.
  * Ross [1950] S. Ross, Variation With Temperature of Surface Tension of Lubricating Oils, Technical Report, United States National Advisory Committee for Aeronautics, 1950\. URL: https://apps.dtic.mil/dtic/tr/fulltext/u2/a801538.pdf.
  * 3M [2022] 3M, $\text{3M}^{\text{TM}}\text{ }\text{Novec}^{\text{TM}}$ 7500 Engineered Fluid, 2022. URL: https://multimedia.3m.com/mws/media/65496O/3m-novec-7500-engineered-fluid.pdf.
  * Brosseau et al. [2014] Q. Brosseau, J. Vrignon, J.-C. Baret,  Microfluidic dynamic interfacial tensiometry (µDIT),  Soft Matter 10 (2014) 3066–3076. URL: http://dx.doi.org/10.1039/C3SM52543K. doi:10.1039/C3SM52543K.
  * 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://www.sciencedirect.com/science/article/pii/S002199911500025X. doi:https://doi.org/10.1016/j.jcp.2015.01.021.
  * Prosperetti [1981] A. Prosperetti,  Motion of two superposed viscous fluids,  The Physics of Fluids 24 (1981) 1217–1223. URL: https://aip.scitation.org/doi/abs/10.1063/1.863522. doi:10.1063/1.863522. arXiv:https://aip.scitation.org/doi/pdf/10.1063/1.863522.
  * Lamb [1975] H. Lamb, Hydrodynamics, 6th ed., Cambridge University Press, 1975\.
  * Shin and Juric [2002] S. Shin, D. Juric,  Modeling three-dimensional multiphase flow using a level contour reconstruction method for front tracking without connectivity,  Journal of Computational Physics 180 (2002) 427–470.