# twoPhaseInterTrackFoam: an OpenFOAM module for Arbitrary Lagrangian/Eulerian
Interface Tracking with Surfactants and Subgrid-Scale Modeling

Moritz Schwarzmeier<EMAIL_ADDRESS>Suraj Raju<EMAIL_ADDRESS>darmstadt.de Željko Tuković<EMAIL_ADDRESS>Mathis Fricke
<EMAIL_ADDRESS>Dieter Bothe<EMAIL_ADDRESS>Tomislav
Marić<EMAIL_ADDRESS>
###### Abstract

We provide an implementation of the unstructured Finite-Volume Arbitrary
Lagrangian / Eulerian (ALE) Interface-Tracking method for simulating
incompressible, immiscible two-phase flows as an OpenFOAM module. In addition
to interface-tracking capabilities that include tracking of two fluid phases,
an implementation of a Subgrid-Scale (SGS) modeling framework for increased
accuracy when simulating sharp boundary layers is enclosed. The SGS modeling
framework simplifies embedding subgrid-scale profiles into the unstructured
Finite Volume discretization. Our design of the SGS model library
significantly simplifies adding new SGS models and applying SGS modeling to
Partial Differential Equations (PDEs) in OpenFOAM.

###### keywords:

finite volume , interface tracking , ALE , unstructured mesh , subgrid-scale

[addrTUDa]organization=Mathematical Modeling and Analysis Institute,
Mathematics department, TU Darmstadt, country=Germany

[addrUniZ]organization=Faculty of Mechanical Engineering and Naval
Architecture, University of Zagreb, country=Croatia

PROGRAM SUMMARY

Program Title: twoPhaseInterTrackFoam
CPC Library link to program files: (to be added by Technical Editor)
Developer’s repository link: https://gitlab.com/interface-
tracking/twophaseintertrackfoamrelease
Code Ocean capsule: (to be added by Technical Editor)
Licensing provisions: GPLv3
Programming language: C++
Nature of problem:
Two-phase flow problems involving surface-active agents (surfactants),
variable surface tension force and very sharp boundary layers.
Solution method:
An OpenFOAM implementation of the Arbitrary Lagrangian / Eulerian Interface
Tracking method.

##  1 Introduction

The unstructured Finite Volume Arbitrary Lagrangian / Eulerian Interface
Tracking method in OpenFOAM (ALE-IT) [31, 29] is a highly accurate method that
tracks fluid phases as deforming solution domains separated by domain
boundaries, making it possible to accurate account for jump conditions at
moving domain boundaries shared by two fluid phases, as well as to solve
Partial Differential Equations (PDEs) on moving domain boundaries. The ALE-IT
method has been successfully applied to challenging two-phase flow problems,
involving, e.g., two-phase flows with soluble surfactant [9, 10, 22] and two-
phase flows with interfacial mass transfer [33, 23].

The Subgrid-Scale (SGS) modeling [1, 4, 34] employs an analytical solution to
enhance the computation of passive scalar transport in interfacial boundary
layers, saving multiple local refinement levels of Finite Volume mesh. This
can accordingly speed up simulations, enable them in the first place or save
(computing) resources, while delivering results with high accuracy.

With this publication, we provide an open-source implementation of our ALE-IT
method as an OpenFOAM module. The version used to generate results from this
manuscript is v1.0, available in the GitLab repository [25] and archived as
[30]. Results reported in this paper are archived as [26]. The following
sections cover the details of the mathematical model, equation discretization,
and subgrid-scale modeling implemented in the version 1.0 of our ALE-IT
OpenFOAM module, as well as the results for some showcases.

##  2 Mathematical model

###  2.1 Arbitrary Lagrangian/Eulerian (Sharp) Interface (Tracking) Model

Gas phase (g), $\rho_{\rm g}$, $\mu_{\rm g}$Liquid phase (l), $\rho_{\rm l}$,
$\mu_{\rm l}$Mesh L$\mathbf{g}$Mesh G$S_{\rm g}$$S_{\rm l}$

Interface

$\sigma$$V_{\rm g}$$V_{\rm l}$$\Sigma(t)$$\mathbf{n}_{\Sigma}$

Figure 1: Definition of the solution domain for the ALE-IT method.

We model isothermal flow of two immiscible incompressible fluids,
schematically shown in fig. 1, e.g. as a liquid and a gas phase separated by a
sharp interface $\Sigma(t)$. Appropriate boundary conditions are then enforced
at the fluid interface $\Sigma(t)$, see [31]. We only briefly outline the
model, the reader is additionally referred to [15] for more details.

Isothermal flow of an incompressible Newtonian fluid inside an arbitrary
moving control volume $V$ bounded by a closed surface $S:=\partial V$, i.e.
$S$ is the boundary of the finite volume, moving with the boundary velocity
$\mathbf{v}_{s}$ is governed by the mass and linear momentum conservation
laws:

$\oint\limits_{S}\rho\mathbf{n}\cdot\mathbf{v}\,\,dS=0,$ (1)
$\frac{\partial}{\partial
t}\int\limits_{V}\rho\mathbf{v}\,\,dV+\oint\limits_{S}\mathbf{n}\cdot\rho(\mathbf{v}-\mathbf{v}_{s})\mathbf{v}\,\,dS=\oint\limits_{S}\mathbf{n}\cdot\boldsymbol{\tau}\,\,dS-\int\limits_{V}\nabla
p\,\,dV,$ (2)

where $\mathbf{n}$ is the outward pointing unit normal on $S$, $\rho$ is the
fluid density, $\mathbf{v}$ is the fluid velocity, $\mathbf{v}_{s}$ is the
velocity of surface $S$, and $\boldsymbol{\tau}$ is the viscous stress tensor
which together with the thermodynamic pressure $p^{\prime}$ makes the Cauchy
stress tensor $\mathbf{T}=-p^{\prime}\mathbf{I}+\boldsymbol{\tau}$. The
dynamic pressure in Eq. (2) is defined as

$p=p^{\prime}-\rho\mathbf{g}\cdot\mathbf{r},$ (3)

where $\mathbf{g}$ is the gravitational acceleration vector and $\mathbf{r}$
is the position vector. Furthermore, the viscous stress tensor of the
Newtonian fluid is given by

$\boldsymbol{\tau}=\eta(\nabla\mathbf{v}+\nabla\mathbf{v}^{T}),$ (4)

where $\eta$ is the dynamic viscosity.

The relationship between the rate of change of the volume $V$ and the velocity
$\mathbf{v}_{s}$ is defined by the geometric (space) conservation law (GCL,
see [27, 8]):

$\frac{\partial}{\partial
t}\int\limits_{V}\,dV-\oint\limits_{S}\mathbf{n}\cdot\mathbf{v}_{s}\,\,dS=0.$
(5)

The formulation of the above mathematical model is a well-known ALE
formulation.

###  2.2 Interface coupling conditions

If fluid phases are immiscible, fluid flow eqs. 1 and 2 can be applied for
each phase separately, while on the interface the proper boundary conditions
must be used. The relation between fluid velocities on the two sides of the
interface is determined by the kinematic condition, which states that the
velocity must be continuous across the interface:

$\mathbf{v}_{\rm l}=\mathbf{v}_{\rm g},$ (6)

where $\mathbf{v}_{\rm l}$ and $\mathbf{v}_{\rm g}$ are the fluid velocities
at the two sides (e.g. liquid and gas) of the fluid interface.

The dynamic condition follows from the momentum conservation law and states
that forces acting on the fluid at the interface are in equilibrium:

$(p_{\rm g}-p_{\rm l})\mathbf{n}_{\Sigma}-(\boldsymbol{\tau}_{\rm
g}-\boldsymbol{\tau}_{\rm
l})\cdot\mathbf{n}_{\Sigma}=\sigma\kappa\mathbf{n}_{\Sigma}+\nabla\\!_{s}\sigma-(\rho_{\rm
g}-\rho_{\rm l})\,(\mathbf{g}\cdot\mathbf{r})\,\mathbf{n}_{\Sigma},$ (7)

where $\mathbf{n}_{\Sigma}$ is the unit normal vector on the interface, which
points from the liquid phase to the gas phase,
$\kappa=-\nabla\\!_{s}{\scriptscriptstyle\bullet}\mathbf{n}_{\Sigma}$ is twice
the mean curvature of the interface, $\sigma$ is the surface tension
coefficient and $\nabla\\!_{s}\sigma$ is the gradient of the surface tension
coefficient, where
$\nabla\\!_{s}=\nabla(\mathbf{I}-\mathbf{n}_{\Sigma}\mathbf{n}_{\Sigma})$ is
the surface gradient operator. The last term on the right-hand side of eq. 7
appears due to transformation of the thermodynamic pressure to the dynamic
pressure at the interface. The pressure jump across the interface is obtained
by taking the normal component of the force balance (eq. 7):

$p_{\rm g}-p_{\rm l}=\sigma\kappa-2\left(\eta_{\rm g}-\eta_{\rm
l}\right)\nabla\\!_{s}\cdot\mathbf{v}-(\rho_{\rm g}-\rho_{\rm
l})\,(\mathbf{g}\cdot\mathbf{r}),$ (8)

The second term on the right-hand side of eq. 8 represents the jump of the
normal viscous force across the interface, expressed through the surface
divergence of the interface velocity. For example, the normal viscous force at
the interface can be expressed as (see [6]):

$(\mathbf{n}_{\Sigma}\mathbf{n}_{\Sigma}):\boldsymbol{\tau}=2\eta\,(\mathbf{n}_{\Sigma}\mathbf{n}_{\Sigma}):\nabla\mathbf{v}=-2\eta\nabla\\!_{s}\cdot\mathbf{v},$
(9)

where the following identity, being valid for an incompressible fluid flow
($\nabla\cdot\mathbf{v}=0$), is used:

$\nabla\cdot\mathbf{v}=\nabla\\!_{s}\cdot\mathbf{v}+(\mathbf{n}_{\Sigma}\mathbf{n}_{\Sigma}):\nabla\mathbf{v}.$
(10)

By taking the tangential component of the force balance (eq. 7) one obtains a
relation between the normal derivative of the tangential velocity on the two
sides of the interface:

$\eta_{\rm g}\left(\mathbf{n}_{\Sigma}\cdot\nabla\mathbf{v}_{\rm
t}\right)_{\rm g}-\eta_{\rm
l}\left(\mathbf{n}_{\Sigma}\cdot\nabla\mathbf{v}_{\rm t}\right)_{\rm
l}=-\nabla\\!_{s}\sigma-\left(\eta_{\rm g}-\eta_{\rm
l}\right)\left(\nabla\\!_{s}v_{\rm n}\right),$ (11)

where it is taken into account that the tangential component of the viscous
force at the interface can be expressed as follows:

$(\mathbf{n}_{\Sigma}\cdot\boldsymbol{\tau})\cdot(\mathbf{I}-\mathbf{n}_{\Sigma}\mathbf{n}_{\Sigma})=\eta\left(\mathbf{n}_{\Sigma}\cdot\nabla\mathbf{v}_{\rm
t}+\nabla\\!_{s}v_{\rm n}\right),$ (12)

where $\mathbf{v}_{\rm
t}=(\mathbf{I}-\mathbf{n}_{\Sigma}\mathbf{n}_{\Sigma})\cdot\mathbf{v}$ is the
tangential velocity component and $v_{\rm
n}=\mathbf{n}_{\Sigma}\cdot\mathbf{v}$ is the normal velocity component at the
interface. A non-zero gradient $\nabla\\!_{s}\sigma$ of the surface tension
coefficient can occur for example due to a nonuniform distribution of
surfactant at the interface or due to the presence of a temperature gradient.

###  2.3 Surfactant transport

If impact of surfactant on the flow field is to be analyzed, it is necessary
to extend the mathematical model by an equation which governs the transport of
surfactant in the liquid bulk and an equation which governs the transport of
surfactant along the interface. The bulk surfactant transport equation in ALE
form read as follows:

$\frac{\partial}{\partial
t}\int_{V}c\,dV+\oint_{S}\mathbf{n}\cdot(\mathbf{v}-\mathbf{v}_{s})\,c\,dS-\oint_{S}\mathbf{n}\cdot(D\nabla
c)\,dS=0,$ (13)

where $c$ is the volume-averaged bulk surfactant concentration and $D$ is the
bulk surfactant diffusion concentration.

The governing equation for surfactant transport along a surface patch $S$
attached to the fluid interface ($S\subset\Sigma(t)$) and bounded by a closed
curve $\partial S$ which can move along the interface reads:

$\frac{\partial}{\partial t}\int_{S}\Gamma\,dS+\oint_{\partial
S}\mathbf{m}\cdot(\mathbf{v}-\mathbf{b})\,\Gamma\,dL-\oint_{\partial
S}\mathbf{m}\cdot(D_{\Gamma}\nabla\\!_{s}\Gamma)\,dL=\int_{S}s_{\Gamma}\,dS,$
(14)

where $\Gamma$ is the surfactant concentration at the interface,
$\Gamma_{\infty}$ is the saturated surfactant concentration under the given
thermodynamic conditions, $\mathbf{m}$ is the outward pointing unit bi-normal
on $\partial S$, $\mathbf{b}$ is the velocity at which the curve $\partial S$
moves along the interface, $L$ is the arc length measured along $\partial S$,
$D_{s}$ is the diffusion coefficient of the surfactant along the interface and
$s_{\Gamma}$ is the source/sink of surfactants per unit area due to adsorption
and desorption.

For the sake of simplicity, we focus in this paper on the prominent Langmuir’s
kinetics law [16, 5, 22] for the transfer of surfactant between the bulk and
the interface due to adsorption and desorption:

$s_{\Gamma}=k_{\rm a}\left[c_{\rm
I}(\Gamma_{\infty}-\Gamma)-\beta\Gamma\right],$ (15)

where $k_{\rm a}$ and $\beta$ are the parameters of the adsorption and
desorption kinetics, respectively, and $c_{\rm I}$ is the value of the bulk
surfactant concentration at the interface. The implementation allows for a
straightforward extension to cope with other sorption models. The normal
derivative of bulk surfactant concentration at the interface is given by
following expression:

$(\mathbf{n}\cdot\nabla c)_{\rm\Sigma}=-\frac{s_{\Gamma}}{D}.$ (16)

In order to account for the surfactant bulk-interface mass transfer balance
(eq. 16), the normal derivative of the surfactant concentration at the
interface needs to be calculated with sufficient accuracy. Due to convection
along the interface and typically small surfactant diffusivities, extremely
thin boundary layers appear and make this calculation a challenging step. At
this point, use of subgrid-scale modeling for improved flux calculation as
described in section 4 is beneficial.

The surface tension is related to the surfactant concentration on the
interface and given for the Langmuir sorption model with its adsorption
isotherm by the following equation of state

$\sigma=\sigma_{0}+\mathcal{R}T\Gamma_{\infty}\ln{\left(1-\frac{\Gamma}{\Gamma_{\infty}}\right)},$
(17)

where $\sigma_{0}$ is the surface tension of a clean interface, $\mathcal{R}$
is the universal gas constant and $T$ is the temperature [16, 5, 22]. Equation
17 can also be written as

$\sigma=\sigma_{0}\left[1+\beta_{s}\ln{\left(1-\frac{\Gamma}{\Gamma_{\infty}}\right)}\right],$
(18)

where $\beta_{s}=\frac{\mathcal{R}T\Gamma_{\infty}}{\sigma_{0}}$ is the
elasticity number. For the test case described in section 6.3 we added a
linear equation of state:

$\sigma=\sigma_{0}\left(1+\beta_{s}\frac{\Gamma}{\Gamma_{\infty}}\right).$
(19)

##  3 Finite Volume and Finite Area discretisation and equation-system
solution

###  3.1 Discretisation of the computational domain

The cell-centered Finite Volume Method (FVM, [31]) is used to discretise the
bulk fluid models (volumetric transport equations) and the face-centered
Finite Area Method (FAM, [31]) is used to discretise the surfactant transport
model, a transport equation on the fluid interface (moving surface). The
finite volume/area discretisation is based on the integral form of the
conservation equation over a fixed or moving control volume/area. The
discretisation procedure is divided into two parts: discretisation of the
computational domain and the equation discretisation. Here we provide a brief
description of both discretisation methods together with a description of the
interface tracking method, while more details can be found in [31, 29, 28].

The time interval is split into a finite number of time-steps $\Delta t$ and
the equations are solved in a time-marching manner. The computational space is
divided into a finite number of polyhedral control volumes (CV) or cells
bounded by polygons. The cells do not overlap and fill the spatial domain
completely. Figure 2 shows a polyhedral control volume $V_{P}$ around the
computational point $P$ located in its centroid, face $f$ with area
$S_{\\!f}$, face unit normal vector $\mathbf{n}_{f}$ and the centroid $N$ of a
neighboring CV sharing the face $f$. The geometry of the CV is fully
determined by the position of its vertices.

$f$$N$$\mathbf{d}_{f}$$P$$V_{P}$$\mathbf{n}_{f}$$z$$y$$x$$\mathbf{r}_{P}$$S_{f}$

Figure 2: Polyhedral control volume (cell) and the notation used (adopted from
[31])

.

In case the solution domain is changing in time, e.g., when the fluid
interface deforms, the ﬁnite volume mesh is adjusted to the time-varying
position or shape of the solution domain boundaries using a deforming mesh
approach. The internal CV-vertices are moved based on the prescribed motion of
the boundary vertices, while the topology of the mesh stays unchanged. In the
implemented moving mesh interface tracking solver, a finite-volume automatic
mesh motion solver is applied, where the Laplace equation is solved for mesh-
point displacements with variable diffusion coefficient, using the polyhedral
unstructured cell-centered FVM.

The polygonal cell-faces which coincide with the interface constitute the
finite area mesh on which the surfactant transport equation is discretised.
Figure 3 shows a sample polygonal control area $S_{\\!P}$ around the
computational point $P$ located in its centroid, the edge $e$, the edge length
$L_{e}$, the edge unit bi-normal vector $\mathbf{m}_{e}$ and the centroid $N$
of the neighbouring control area sharing the edge $e$. The bi-normal
$\mathbf{m}_{e}$ is perpendicular to the edge normal $\mathbf{n}_{e}$ and to
the edge vector $\mathbf{e}$.

$P$$e$$\mathbf{m}_{e}$$\mathbf{n}_{P}$$j$$N$$S_{P}$$\mathbf{e}$$\mathbf{d}_{e}$$L_{e}$$\mathbf{n}_{e}$$i$$\mathbf{n}_{N}$

Figure 3: Polygonal control area.

###  3.2 Discretisation of volumetric equations

The second-order collocated FV discretisation of an integral conservation
equation transforms the surface integrals into sums of face integrals and
approximates them and the volume integrals to the second order accuracy by
using the mid-point rule. The spatially discretised form of the momentum eq. 2
for the moving control volume $V_{P}$ reads:

$\begin{split}&\frac{\partial\left(\mathbf{v}_{P}V_{P}\right)}{\partial
t}+\sum_{f}(\dot{V}_{f}-\dot{V}_{s,f})\,\mathbf{v}_{f}\\\
&=\sum_{f}\nu_{f}\,[\mathbf{n}\cdot(\nabla\mathbf{v})]_{f}\,S_{f}+(\nabla
p)_{P}V_{P},\end{split}$ (20)

where the subscripts P and f represent the cell-center and face-center values.
The volume flow rate due to fluid flow through the cell-face,
$\dot{V}^{n+1}_{f}=(\mathbf{n}\cdot\mathbf{v}S)^{n+1}_{f}$ must satisfy the
discretised mass conservation law, where the superscript n+1 represents the
new time-step, while the volume flux due to grid motion,
$\dot{V}^{n+1}_{s,f}$, must satisfy the discretised GCL, eq. 5. The evaluation
of the mesh-face volume fluxes is described in [11, 29].

The face-center values of all dependent variables are calculated using linear
interpolation of the neighbouring cell-center values. The exception is the
face value of the dependent variable in the convection term [$\mathbf{v}_{f}$
in eq. 20] which must be calculated using some of the bounded convection
discretisation schemes available in OpenFOAM.

The face-normal derivative of the velocity
$[\mathbf{n}\cdot\nabla\mathbf{v}]_{f}$ is discretised using the linear scheme
with non-orthogonal and skewness correction (see fig. 2):

$(\mathbf{n}\cdot\nabla\mathbf{v})_{f}=\frac{\mathbf{v}_{N}-\mathbf{v}_{P}}{d_{{\rm
n},f}}+\frac{\mathbf{k}_{N}\cdot(\nabla\mathbf{v})_{N}-\mathbf{k}_{P}\cdot(\nabla\mathbf{v})_{P}}{d_{{\rm
n},f}}\,,$ (21)

where $d_{{\rm n},f}=\mathbf{n}_{f}\cdot\mathbf{d}_{f}$,
$\mathbf{k}_{P}=(\mathbf{I}-\mathbf{n}_{f}\mathbf{n}_{f})\cdot(\mathbf{r}_{f}-\mathbf{r}_{P})$
and
$\mathbf{k}_{N}=(\mathbf{I}-\mathbf{n}_{f}\mathbf{n}_{f})\cdot(\mathbf{r}_{N}-\mathbf{r}_{f})$.
The first term on the right-hand side of eq. 21 is treated implicitly, while
the correction term is treated explicitly.

The cell-face volumetric fluid flux in the non-linear convection term is
treated explicitly after the discretisation, i.e. its values from the previous
iteration are used. The cell-center gradient of the dependent variables used
for the calculation of the explicit source terms and the non-orthogonal
correction in eq. 21, is obtained by using the least-squares fit [7]. This
method produces a second-order accurate gradient irrespective of the mesh
quality.

The temporal discretisation of all equations is performed by using the second-
order accurate implicit three time level scheme [11] referred to as the
backward scheme in OpenFOAM. All terms of eq. 20 are evaluated at the new time
instance $t^{n+1}=t^{n}+\Delta t$ and the temporal derivative is discretised
by using two old-time levels:

$\left(\\!\frac{\partial\mathbf{v}}{\partial
t}\\!\right)^{\\!{n+1}}=\frac{3\mathbf{v}^{n+1}-4\mathbf{v}^{n}+\mathbf{v}^{n-1}}{2\Delta
t},$ (22)

where $\mathbf{v}^{n+1}=\mathbf{v}(t^{n}+\Delta t)$,
$\mathbf{v}^{n}=\mathbf{v}(t^{n})$ and
$\mathbf{v}^{n-1}=\mathbf{v}(t^{n}-\Delta t)$. One can also use other temporal
discretisation schemes available in OpenFOAM.

When eqs. 21 and 22 are substituted into eq. 20 and the convection
discretisation scheme is applied, the discretised form of the momentum eq. 2
can be written in the form of a linear algebraic equation, which for cell $P$
reads:

$a_{P}\,\mathbf{v}^{n+1}_{P}+\sum_{N}a_{N}\,\mathbf{v}^{n+1}_{N}=\mathbf{r}_{P}+(\nabla
p)_{P},$ (23)

where the diagonal coefficient $a_{P}$, the off-diagonal coefficients $a_{N}$
and the source term $\mathbf{r}_{P}$ can be found in [29].

The mathematical model of fluid flow is solved using a segregated solution
procedure where the discretised momentum equation is solved decoupled from the
discretised pressure equation. The discretised pressure Poisson equation,
obtained by combining the discretised momentum and continuity equations,
reads:

$\sum_{f}\left(\frac{1}{a}\right)_{f}(\mathbf{n}\cdot\nabla
p)^{n+1}_{f}\,S^{n+1}_{f}=\sum_{f}\mathbf{n}^{n+1}_{f}\cdot\left(\frac{\mathbf{H}}{a}\right)_{f}S^{n+1}_{f},$
(24)

where the $(1/a)_{f}$ and $(\mathbf{H}/a)_{f}$ terms are calculated by using
the temporally consistent Rhie-Chow interpolation procedure proposed in [29].

The face normal derivative of the pressure at the left hand side of eq. 24 is
discretised using eq. 21 applied on the pressure field. The absolute volume
fluid flux $\dot{V}^{n+1}_{\\!f}$ through the cell face $f$ is calculated as
follows:

$\dot{V}^{n+1}_{f}=\mathbf{n}^{n+1}_{f}\cdot\left[\left(\frac{\mathbf{H}}{a}\right)_{f}-\left(\frac{1}{a}\right)_{f}(\nabla
p)^{n+1}_{f}\right]S^{n+1}_{\\!f}.$ (25)

Volume flux calculated in this manner will satisfy the discretised continuity
equation if the pressure field satisfies the pressure eq. 24.

###  3.3 Discretisation of surface equations

Applying a second-order collocated finite area method, the surfactant
transport eq. 14 can be discretised on the moving control area $S_{\\!P}$ (see
fig. 3) as follows:

$\frac{\partial\left(\Gamma_{\\!P}\,S_{\\!P}\right)}{\partial
t}+\sum_{e}(\mathbf{m}\cdot\mathbf{v})_{e}L_{e}\,\Gamma_{e}=\sum_{e}D_{\Gamma,e}(\mathbf{m}\cdot\nabla\\!_{s}\Gamma)_{e}L_{e}+(s_{\Gamma})_{P}\,S_{\\!P},$
(26)

where the subscripts P and e represent the face-center and edge-center values,
and it is assumed that finite-area points move in normal direction which means
that $(\mathbf{m}\cdot\mathbf{b})_{e}=0$.

The edge-center tangential velocity $\mathbf{v}_{\rm t}$ is calculated using
following linear interpolation formula:

$\mathbf{v}_{{\rm t},e}=(\mathbf{T}_{e})^{\rm
T}\cdot\left[e_{x}\mathbf{T}_{P}\cdot\mathbf{v}_{{\rm
t},P}+(1-e_{x})\mathbf{T}_{N}\cdot\mathbf{V}_{{\rm t},N}\right].$ (27)

where $e_{x}$ is the interpolation factor which is calculated as the ratio of
the geodetic distances $\overline{eN}$ and $\overline{PeN}$ (see fig. 3):

$e_{x}=\frac{\overline{eN}}{\overline{PeN}},$ (28)

and $\mathbf{T}_{P}$, $\mathbf{T}_{N}$ and $\mathbf{T}_{e}$ are the tensors of
transformation from the global Cartesian coordinate system to the
$(\mathbf{t}^{\prime}_{e},\mathbf{t}_{e},\mathbf{n}_{e})$ edge-based local
coordinate system defined in fig. 3. The convection term in eq. 26 is
discretised using the linear upwind discretisation scheme by taking into
account geodetic distances between the neighbouring control area centers. The
diffusion term is discretised using the central differencing scheme with non-
orthogonal correction [14]:

$(\mathbf{m}\cdot\nabla\\!_{s}\psi)_{e}=\mathop{\mathchoice{\vtop{\halign{#\cr$\hfil\displaystyle\left|\mathbf{\Delta}_{e}\right|\frac{\psi_{N}-\psi_{P}}{L_{PN}}\hfil$\crcr\kern
2.0pt\nointerlineskip\cr\@mathmeasure\displaystyle{\left|\mathbf{\Delta}_{e}\right|
\frac{\psi_{N}-\psi_{P}}{L_{PN}}}\@mathmeasure\displaystyle{\upbrace}\@mathmeasure\displaystyle{\upbraceg}\@mathmeasure\displaystyle{\upbracegg}\@mathmeasure\displaystyle{\upbraceggg}\@mathmeasure\displaystyle{\upbracegggg}$\displaystyle\bracelu\leaders\hbox{$\bracemid$}\hfill\bracemu\leaders\hbox{$\bracemid$}\hfill\braceru$\crcr}}}{\vtop{\halign{#\cr$\hfil\textstyle\left|\mathbf{\Delta}_{e}\right|\frac{\psi_{N}-\psi_{P}}{L_{PN}}\hfil$\crcr\kern
2.0pt\nointerlineskip\cr\@mathmeasure\textstyle{\left|\mathbf{\Delta}_{e}\right|
\frac{\psi_{N}-\psi_{P}}{L_{PN}}}\@mathmeasure\textstyle{\upbrace}\@mathmeasure\textstyle{\upbraceg}\@mathmeasure\textstyle{\upbracegg}\@mathmeasure\textstyle{\upbraceggg}\@mathmeasure\textstyle{\upbracegggg}$\textstyle\bracelu\leaders\hbox{$\bracemid$}\hfill\bracemu\leaders\hbox{$\bracemid$}\hfill\braceru$\crcr}}}{\vtop{\halign{#\cr$\hfil\scriptstyle\left|\mathbf{\Delta}_{e}\right|\frac{\psi_{N}-\psi_{P}}{L_{PN}}\hfil$\crcr\kern
2.0pt\nointerlineskip\cr\@mathmeasure\scriptstyle{\left|\mathbf{\Delta}_{e}\right|
\frac{\psi_{N}-\psi_{P}}{L_{PN}}}\@mathmeasure\scriptstyle{\upbrace}\@mathmeasure\scriptstyle{\upbraceg}\@mathmeasure\scriptstyle{\upbracegg}\@mathmeasure\scriptstyle{\upbraceggg}\@mathmeasure\scriptstyle{\upbracegggg}$\scriptstyle\bracelu\leaders\hbox{$\bracemid$}\hfill\bracemu\leaders\hbox{$\bracemid$}\hfill\braceru$\crcr}}}{\vtop{\halign{#\cr$\hfil\scriptscriptstyle\left|\mathbf{\Delta}_{e}\right|\frac{\psi_{N}-\psi_{P}}{L_{PN}}\hfil$\crcr\kern
2.0pt\nointerlineskip\cr\@mathmeasure\scriptscriptstyle{\left|\mathbf{\Delta}_{e}\right|
\frac{\psi_{N}-\psi_{P}}{L_{PN}}}\@mathmeasure\scriptscriptstyle{\upbrace}\@mathmeasure\scriptscriptstyle{\upbraceg}\@mathmeasure\scriptscriptstyle{\upbracegg}\@mathmeasure\scriptscriptstyle{\upbraceggg}\@mathmeasure\scriptscriptstyle{\upbracegggg}$\scriptscriptstyle\bracelu\leaders\hbox{$\bracemid$}\hfill\bracemu\leaders\hbox{$\bracemid$}\hfill\braceru$\crcr}}}}\limits_{\textrm{Orthogonal
contribution}}+\mathop{\mathchoice{\vtop{\halign{#\cr$\hfil\displaystyle\mathbf{k}_{e}\cdot(\nabla\\!_{s}\psi)_{e},\hfil$\crcr\kern
2.0pt\nointerlineskip\cr\@mathmeasure\displaystyle{\mathbf{k}_{e}\cdot($\nabla\\!_{s}$\psi)_{e},}\@mathmeasure\displaystyle{\upbrace}\@mathmeasure\displaystyle{\upbraceg}\@mathmeasure\displaystyle{\upbracegg}\@mathmeasure\displaystyle{\upbraceggg}\@mathmeasure\displaystyle{\upbracegggg}$\displaystyle\bracelu\leaders\hbox{$\bracemid$}\hfill\bracemu\leaders\hbox{$\bracemid$}\hfill\braceru$\crcr}}}{\vtop{\halign{#\cr$\hfil\textstyle\mathbf{k}_{e}\cdot(\nabla\\!_{s}\psi)_{e},\hfil$\crcr\kern
2.0pt\nointerlineskip\cr\@mathmeasure\textstyle{\mathbf{k}_{e}\cdot($\nabla\\!_{s}$\psi)_{e},}\@mathmeasure\textstyle{\upbrace}\@mathmeasure\textstyle{\upbraceg}\@mathmeasure\textstyle{\upbracegg}\@mathmeasure\textstyle{\upbraceggg}\@mathmeasure\textstyle{\upbracegggg}$\textstyle\bracelu\leaders\hbox{$\bracemid$}\hfill\bracemu\leaders\hbox{$\bracemid$}\hfill\braceru$\crcr}}}{\vtop{\halign{#\cr$\hfil\scriptstyle\mathbf{k}_{e}\cdot(\nabla\\!_{s}\psi)_{e},\hfil$\crcr\kern
2.0pt\nointerlineskip\cr\@mathmeasure\scriptstyle{\mathbf{k}_{e}\cdot($\nabla\\!_{s}$\psi)_{e},}\@mathmeasure\scriptstyle{\upbrace}\@mathmeasure\scriptstyle{\upbraceg}\@mathmeasure\scriptstyle{\upbracegg}\@mathmeasure\scriptstyle{\upbraceggg}\@mathmeasure\scriptstyle{\upbracegggg}$\scriptstyle\bracelu\leaders\hbox{$\bracemid$}\hfill\bracemu\leaders\hbox{$\bracemid$}\hfill\braceru$\crcr}}}{\vtop{\halign{#\cr$\hfil\scriptscriptstyle\mathbf{k}_{e}\cdot(\nabla\\!_{s}\psi)_{e},\hfil$\crcr\kern
2.0pt\nointerlineskip\cr\@mathmeasure\scriptscriptstyle{\mathbf{k}_{e}\cdot($\nabla\\!_{s}$\psi)_{e},}\@mathmeasure\scriptscriptstyle{\upbrace}\@mathmeasure\scriptscriptstyle{\upbraceg}\@mathmeasure\scriptscriptstyle{\upbracegg}\@mathmeasure\scriptscriptstyle{\upbraceggg}\@mathmeasure\scriptscriptstyle{\upbracegggg}$\scriptscriptstyle\bracelu\leaders\hbox{$\bracemid$}\hfill\bracemu\leaders\hbox{$\bracemid$}\hfill\braceru$\crcr}}}}\limits_{\textrm{Non-
orthogonal correction}}$ (29)

where $\mathbf{k}_{e}=\mathbf{m}_{e}-\mathbf{\Delta}_{e}$,
$\mathbf{\Delta}_{e}=\frac{\mathbf{t}_{e}}{\mathbf{t}_{e}\cdot\mathbf{m}_{e}}$,
$L_{PN}$ is the geodetic distance $\overline{LeP}$ and $\mathbf{t}_{e}$ is the
unit tangential vector to the geodetic line $\overline{PeN}$ at the point $e$
(see fig. 4).

$\mathbf{n}_{P}$$\mathbf{n}_{N}$$S_{{\rm
l}f}$$N$$\mathbf{t}^{\prime}_{P}$$\mathbf{t}_{P}$$\mathbf{t}_{e}$$\mathbf{t}_{N}$$\mathbf{t}^{\prime}_{N}$$\mathbf{t}^{\prime}_{e}$$e$${\rm
l}f\equiv P$$\mathbf{n}_{e}$

Figure 4: Edge-based local orthogonal coordinate system whose axis are aligned
with the orthogonal unit vectors $\mathbf{n}$, $\mathbf{t}$ and
$\mathbf{t}^{\prime}$, where vector $\mathbf{t}$ is tangential to the geodetic
line $\overline{PeN}$

###  3.4 Interface tracking procedure

The numerical modelling of the two-phase fluid flow with a sharp interface is
performed using a moving mesh interface tracking procedure. The computational
mesh consists of two separate parts, where each of the meshes covers only one
of the considered two fluid phases, see fig. 1. The two meshes are in contact
over two geometrically equal surfaces, $S_{\\!\rm l}$ and $S_{\\!\rm g}$, at
the boundary between the phases, i.e. the interface. Surface $S_{\\!\rm l}$
represents the liquid side of the interface, and surface $S_{\\!\rm g}$
represent the gas side of the interface. Each surface is defined by a set of
boundary faces, see fig. 5, where each face ${\rm l}f$ on the surface
$S_{\\!\rm l}$ has a corresponding geometrically equal face ${\rm g}f$ on the
surface $S_{\\!\rm g}$. Matching of the two meshes at the interface is assumed
in order to make the explanation of the interface tracking method clearer, and
is not required in general.

$\mathbf{n}_{f}$${\rm l}f$${\rm g}f$Phase gPhase l

Figure 5: Representation of the interface with the mesh boundary faces.

Coupling of flow equations between fluid phases is performed by applying
adequate boundary conditions at the boundary faces which define the side $\rm
l$ and $\rm g$ of the interface. At the side $\rm l$ the pressure $p_{\rm l}$
and normal velocity derivative $(\mathbf{n}\cdot\nabla\mathbf{v})_{\rm l}$ are
specified, while at the side $\rm g$ the velocity $\mathbf{v}_{\rm g}$ is
specified and normal derivative of dynamic pressure is set to zero,
$(\mathbf{n}\cdot\nabla p)_{\rm g}=0$. The boundary conditions are calculated
using the kinematic and dynamic conditions as follows:

  1. 1.

The value of dynamic pressure specified on the face ${\rm l}f$ (see fig. 5) is
calculated from the dynamic pressure at the face ${\rm g}f$ using eq. 8 as
follows:

$\begin{split}p_{{\rm l}f}&=p_{{\rm g}f}-(\rho_{\rm l}-\rho_{\rm
g})\,\mathbf{g}\cdot\mathbf{r}_{{\rm l}f}\\\ &-(\sigma\kappa)_{{\rm
l}f}-2\left(\eta_{\rm l}-\eta_{\rm
g}\right)(\nabla\\!_{s}\cdot\mathbf{v})_{{\rm l}f},\end{split}$ (30)

where $p_{{\rm g}f}$ is the dynamic pressure at the face ${\rm g}f$ calculated
by extrapolation from the fluid $\rm g$, and $\mathbf{r}_{{\rm l}f}$ is the
position vector of the face center ${\rm l}f$. The surface divergence of the
velocity vector ${(\nabla\\!_{s}\cdot\mathbf{v})_{{\rm l}f}}$ at ${\rm l}f$ is
calculated using the surface Gauss integral theorem [32]. The procedure for
calculating the surface force $(\sigma\kappa)_{{\rm l}f}$ is described later
in this section.

  2. 2.

The normal velocity derivative specified at ${\rm l}f$ is calculated from the
normal velocity gradient at ${\rm g}f$ using eq. 11, as follows:

$\begin{split}\left(\mathbf{n}\cdot\nabla\mathbf{v}\right)_{{\rm
l}f}&=\frac{\eta_{\rm g}}{\eta_{\rm l}}\,(\mathbf{I}-\mathbf{n}_{{\rm
l}f}\mathbf{n}_{{\rm l}f})\cdot\left[\mathbf{n}_{{\rm
l}f}\cdot\left(\nabla\mathbf{v}\right)_{{\rm g}f}\right]\\\
&+\frac{1}{\eta_{\rm l}}\,(\nabla\\!_{s}\sigma)_{{\rm l}f}-\mathbf{n}_{{\rm
l}f}\,(\nabla\\!_{s}\cdot\mathbf{v})_{{\rm l}f}\\\ &+\frac{(\eta_{\rm
g}-\eta_{\rm l})}{\eta_{\rm l}}\left(\nabla\\!_{s}\mathbf{v}_{\rm
n}\right)_{{\rm l}f},\end{split}$ (31)

where $\mathbf{n}_{{\rm l}f}=-\mathbf{n}_{{\rm g}f}$ is the unit normal of the
face ${\rm l}f$. Equation 31 is derived using the identity
$\mathbf{n}\cdot\nabla\mathbf{v}+\mathbf{n}\,(\nabla\\!_{s}\cdot\mathbf{v})=(\mathbf{I}-\mathbf{n}\mathbf{n})\cdot(\mathbf{n}\cdot\nabla\mathbf{v})$.
The surface gradient of the normal velocity component
$(\nabla\\!_{s}\mathbf{v}_{\rm n})_{{\rm l}f}$ at ${\rm l}f$ is calculated
using the surface Gauss integral theorem. Calculating the tangential surface
force $(\nabla\\!_{s}\sigma)_{{\rm l}f}$ is described later in this section.

  3. 3.

According to the kinematic condition (6), the tangential velocity component
specified on ${\rm g}f$ is transferred from ${\rm l}f$:

$(\mathbf{v}_{\rm t})_{{\rm g}f}=(\mathbf{v}_{\rm t})_{{\rm l}f}.$ (32)

The normal velocity component specified on ${\rm g}f$ is calculated from the
condition of zero net mass flux, $(\dot{V}_{{\rm g}f}-\dot{V}_{{\rm g}f})=0$,
i.e.

$(\mathbf{v}_{\rm n})_{{\rm g}f}=-\frac{\dot{V}_{{\rm g}f}}{S_{{\rm
g}f}}\mathbf{n}_{{\rm l}f},$ (33)

where $\dot{V}_{{\rm g}f}$ is the volume flux of the face ${\rm g}f$. Since
the displacement of mesh points on the side $\rm g$ is equal to the
displacement on the side $\rm l$ of the interface, $\mathbf{u}_{{\rm
g}i}=\mathbf{u}_{{\l}i}$, the same is valid for the volume fluxes of the faces
${\rm l}f$ and ${\rm g}f$:

$\dot{V}_{{\rm g}f}=\dot{V}_{{\rm l}f}.$ (34)

${\rm l}i$${{\rm l}f}$$\mathbf{f}_{{\rm l}i}$${\rm l}c$interface pointscontrol
points$\mathbf{n}_{{\rm l}i}$$\mathbf{f}_{{\rm l}c}$$\mathbf{n}_{{\rm l}f}$

Figure 6: Definition of the interface using control points.

The specification of interface boundary conditions using the above procedure
is done at the beginning of each outer iteration. In general, at the end of an
outer iteration, the net volume flux through the boundary side $\rm l$ is
different from zero, i.e.

$[\dot{V}^{\rm p}_{{\rm l}f}-(\dot{V}_{s})^{\rm p}_{{\rm l}f}]\neq 0,$ (35)

where $\dot{V}^{\rm p}_{{\rm l}f}$ is the fluid volume flux through the face
${\rm l}f$ and $(\dot{V}_{s})^{\rm p}_{{\rm l}f}$ is the cell-face volume flux
through the same face, both obtained in the previous outer iteration. In order
to correct the net volume flux, the interface points must be moved to
accomplish the following volume flux corrections:

$\dot{V}^{\prime}_{{\rm l}f}=[\dot{V}^{\rm p}_{{\rm l}f}-(\dot{V}_{s})^{\rm
p}_{{\rm l}f}].$ (36)

The displacement of the interface points is calculated based on the procedure
proposed by [21], where a control point ${\rm l}c$ is attached to the each
face ${\rm l}f$ at the side $\rm l$ of the interface as is shown in fig. 6.
The corrected position of the interface points is calculated using the
following procedure:

  1. 1.

Calculate the volume $\delta V^{\prime}_{{\rm l}f}$ which face ${\rm l}f$
sweeps on the way from the current to the corrected position in order to
cancel the net mass flux:

$\delta V^{\prime}_{{\rm l}f}=\frac{2}{3}\,\dot{V}^{\prime}_{{\rm l}f}\,\Delta
t,$ (37)

where $\dot{V}^{\prime}_{{\rm l}f}$ is the volume flux correction for face
${\rm l}f$, eq. 36.

  2. 2.

Using the above, the displacement of control points in the direction
$\mathbf{f}_{{\rm l}c}$ is:

$h^{\prime}_{{\rm l}c}=\frac{\delta V^{\prime}_{{\rm l}f}}{S^{\rm p}_{{\rm
l}f}\,\mathbf{n}^{\rm p}_{{\rm l}f}\cdot\mathbf{f}_{{\rm l}c}},$ (38)

where $S^{\rm p}_{{\rm l}f}$ and $\mathbf{n}^{\rm p}_{{\rm l}f}$ are the area
and unit normal of the face ${\rm l}f$ in the previous iteration and
$\mathbf{f}_{{\rm l}c}$ is the control point displacement direction. The new
corrected positions of the control points are calculated according to the
following expression:

$\mathbf{r}^{n}_{{\rm l}c}=\mathbf{r}^{\rm p}_{{\rm l}c}+h^{\prime}_{{\rm
l}c}\,\mathbf{f}_{{\rm l}c},$ (39)

where $\mathbf{r}^{\rm p}_{{\rm l}c}$ is the position vector of the control
point ${\rm l}c$ before the correction.

  3. 3.

The new position of the interface mesh point ${\rm l}i$ is obtained by
projection to the plane which is laid over the corresponding control points
using the least square method.

###  3.5 Calculation of interface curvature and surface tension

Regardless of the approach used for tracking the interface between the phases
in a multiphase fluid flow, the implementation of surface tension is always
demanding.

Let us assume the interface is discretised with an unstructured surface mesh
consisting of arbitrary polygonal control areas. The surface tension force
acting on the control area $S_{\\!\rm lf}$ (see fig. 7) can be expressed by
the following equation:

$\mathbf{F}_{{\rm l}f}^{\sigma}=\oint\limits_{\partial S_{{\rm
l}f}}\mathbf{m}\sigma\,\,dL=\sum_{e}\int\limits_{L_{e}}\mathbf{m}\sigma\,\,dL=\sum_{e}\sigma_{e}\mathbf{m}_{e}L_{e},$
(40)

where $\sigma_{e}$ and $\mathbf{m}_{e}$ are the surface tension and bi-normal
unit vector at the center of the edge $e$ and $L_{e}$ is the length of the
edge $e$. If the total surface tension force for each control area in the mesh
is calculated using eq. 40, then the total surface tension force for a closed
surface will be exactly zero if the unit bi-normals $\mathbf{m}_{e}$ for two
control areas sharing the edge $e$ are parallel and have opposite direction.

$\mathbf{m}_{e}$$e$$\hat{\mathbf{e}}$$\mathbf{n}_{e}$$\mathbf{n}_{i}$$\mathbf{n}_{j}$$i$$j$${\rm
l}f$$L_{e}$$S_{{\rm l}f}$$\mathbf{n}_{{\rm l}f}$

Figure 7: Control area $S_{{\rm l}f}$ at the interface.

It remains to decompose the surface tension force $\mathbf{F}_{{\rm
l}f}^{\sigma}$ into the tangential component $(\nabla\\!_{s}\sigma)_{{\rm
l}f}$ used in eq. 31, and the normal component $(\kappa\sigma)_{{\rm
l}f}\mathbf{n}_{{\rm l}f}$ used in eq. 30. Using the surface Gauss integral
theorem, the surface tension force acting on the control area $S_{\\!{\rm
l}f}$ can be expressed by the following equation:

$\mathbf{F}_{{\rm l}f}^{\sigma}=\int_{S_{{\rm
l}f}}\nabla\\!_{s}\sigma\,\,dS+\int_{S_{{\rm
l}f}}\kappa\sigma\mathbf{n}\,\,dS.$ (41)

When the right hand side of eq. 41 is approximated by the mid-point role and
the result of the discretisation is equalised with the right hand side of eq.
40, the following expression is obtained:

$(\nabla\\!_{s}\sigma)_{{\rm l}f}+(\kappa\sigma)_{{\rm l}f}\,\mathbf{n}_{{\rm
l}f}=\frac{1}{S_{{\rm l}f}}\sum_{e}\sigma_{e}\mathbf{m}_{e}L_{e}.$ (42)

Hence, the tangential component of the surface tension force acting on the
control area $S_{\\!{\rm l}f}$ is equal to the tangential component of the
right hand side of eq. 42:

$(\nabla\\!_{s}\sigma)_{{\rm l}f}=\frac{1}{S_{{\rm
l}f}}\,(\mathbf{I}-\mathbf{n}_{{\rm l}f}\mathbf{n}_{{\rm
l}f})\cdot\sum_{e}\sigma_{e}\mathbf{m}_{e}L_{e}$ (43)

and its normal component is equal to the respective normal component:

$(\kappa\sigma)_{{\rm l}f}\mathbf{n}_{{\rm l}f}=\frac{1}{S_{{\rm
l}f}}\,(\mathbf{n}_{{\rm l}f}\mathbf{n}_{{\rm
l}f})\cdot\sum_{e}\sigma_{e}\mathbf{m}_{e}L_{e}.$ (44)

If the surface tension coefficient is constant ($\sigma=\textrm{const.}$), eq.
43 will give a tangential component of the surface tension force equal to zero
if the normal unit vector of the control area $S_{\\!\rm{\rm l}f}$ satisfies
the following equation:

$\kappa_{{\rm l}f}\mathbf{n}_{{\rm l}f}=\frac{1}{S_{\\!{\rm
l}f}}\sum_{e}\mathbf{m}_{e}L_{e}$ (45)

or, if $\kappa_{{\rm l}f}\neq 0$:

$\mathbf{n}_{{\rm
l}f}=\frac{\sum\limits_{e}\mathbf{m}_{e}L_{e}}{|\sum\limits_{e}\mathbf{m}_{e}L_{e}|}.$
(46)

With eqs. 43, 44 and 45 we shall formulate a procedure for the calculation of
the surface tension force which ensures that the total surface tension force
on a closed surface will be exactly zero. Unfortunately, the fulfilment of
this condition is not sufficient for successful application of surface tension
forces in the calculation. Specifically, unphysical fluid flow near the
interface arises due to local (rather than global) inaccuracy in the
calculation of surface tension forces.

From eq. 40 one can see that the accuracy of surface tension force calculation
depends on the accuracy of calculation of the bi-normal unit vector
$\mathbf{m}_{e}$ which is calculated using the following expression:

$\mathbf{m}_{e}=\hat{\mathbf{e}}\times\frac{\mathbf{n}_{i}+\mathbf{n}_{j}}{2},$
(47)

where $\hat{\mathbf{e}}$ is the unit vector parallel with edge $e$ and
$\mathbf{n}_{i}$ and $\mathbf{n}_{j}$ are the interface normal unit vectors in
points $i$ and $j$ (see fig. 7). Using eqs. 45 and 47 one obtains the exact
value of curvature of the control area $S_{\\!{\rm l}f}$ if the points of the
control area lie on the surface of the sphere.

###  3.6 Solution procedure

Based on the described interface tracking method, one can now define the
solution procedure for the Navier-Stokes system on a moving mesh, which may be
used for simulating two-phase fluid flow using a moving mesh interface
tracking method. The procedure consists of the following steps:

  1. 1.

For the new time instance $t=t^{n+1}$, initialize the values of all dependent
variables with the corresponding values from the previous time instance;

  2. 2.

Define the displacement directions for the interface mesh points and the
control points;

  3. 3.

Start of outer iteration loop:

    1. (a)

Update pressure and velocity boundary conditions at the interface;

    2. (b)

Assemble and solve the discretised momentum equation eq. 20 on the mesh with
the current shape of the interface. The pressure field, face mass fluxes and
volume fluxes are used from the previous (outer) iteration;

    3. (c)

The velocity field obtained in the previous step is used for the assembly of
the discretised pressure eq. 24. After the pressure equation is solved, new
absolute mass fluxes through the cell faces are calculated. The net mass flux
through an interface is generally different from zero;

    4. (d)

In order to compensate the net mass flux obtained in the previous step, the
interface displacement is calculated using the interface points displacement
procedure defined in the previous section;

    5. (e)

The interface points displacement is used as a boundary condition for the
solution of the mesh motion problem. After mesh displacement, the new face
volume fluxes are calculated using the current points positions and the
position from the previous time instance;

    6. (f)

Convergence is checked and if the residual levels and the net mass flux
through the interface do not satisfy the prescribed accuracy, the procedure is
returned to step (a).

  4. 4.

If the final time instance is not reached, return to step 1.

##  4 Subgrid-Scale Model

###  4.1 Motivation

The simulation of realistic gas-liquid systems is still a huge challenge, e.g.
when one would like to simulate an industrial scale bubble column reactor. The
nature of this application is multi-scale. One of the smallest scales
occurring is that of concentration boundary layers of the rising bubbles at
realistic, i.e. high, Schmidt numbers. With conventional FV-methods these
scales have to be resolved, what is still unfeasible for an application as the
one mentioned above. This motivated the first paper to discuss SGS modeling
[1]. In SGS modeling, an analytical function is used to describe the profile
of a scalar field inside its boundary layer and to use this SGS information
for improved flux computation.

The method has been developed further since and has been applied in the
context of Volume of Fluid [34, 13, 4] as well as in interface tracking
frameworks [22, 23, 33]. An overview on the state of subgrid-scale modeling
was also given by [35]. Recently, the approach was taken up in [12], where
local curvature and tangential convection effects have been included.

###  4.2 Analytical model

Figure 8: The simplified model for species transfer across a bubble surface.
[34] adapted from [4]

The SGS model we implement is based upon analytical solution of a simplified
problem, depicted for a bubble in fig. 8. The idea behind the SGS modeling is
to zoom in on the interface, so that the curvature gets negligibly small,
allowing the analytical model therefore to assume a (locally) flat boundary.
This assumption is adequate for scenarios with a concentration boundary layer
thickness small in comparison to the radius of curvature of the given
boundary.

At the interface $\Sigma$ of the model problem (see fig. 8), a constant and
homogeneous concentration $c_{\Sigma}$ is prescribed, which can be computed
from Henry’s law under the sensible assumption of a well-mixed gas-side
concentration, as is pointed out by [4].

The velocity in the substitute model problem is parallel to the interface. No
velocity gradient normal to the boundary exists in the model problem
(${\partial v}/{\partial x}=0$). There could be a non-zero gradient in the
boundary-parallel direction, in which case the boundary layer thickness
$\delta$ in eq. 49 would change accordingly. There is no velocity component
normal to the boundary.

In the application within ALE-IT, the boundary geometry and the velocity field
have to be resolved. In a co-moving reference frame, the interface
neighbouring cell’s velocity should be nearly parallel to the interface and
with the assumption of zero interface-normal velocity gradient the next cells
velocity should be nearly identical to it. When a thin scalar boundary layer
exists, the SGS can model the scalar fluxes normal to the boundary without
resolving the scalar boundary layer itself.

Further assumptions of the substitute problem are a given, constant value
(e.g. zero concentration) upstream ($y<0$) of the domain and far away from the
boundary ($x\to\infty$) as well as negligible diffusion in streamwise
direction. We also assume quasi-steadiness and the transported scalar being
passively transported.

Applying these assumptions to the advection-diffusion eq. 13 with the sorption
term ${s^{\Sigma}+[\\![\boldsymbol{j}\cdot\boldsymbol{n}_{\Sigma}]\\!]=0\text{
on }\Sigma(t)}$ and Fick’s law $\boldsymbol{j}=-D\nabla c$ and $\nu$ being the
boundary-parallel velocity, lead to the expression

$\nu\frac{\partial c}{\partial y}=D\frac{\partial^{2}c}{\partial
x^{2}}\quad\text{for $x>0$ and $y>0$}$ (48)

with the boundary conditions

$c(x,y=0)=c_{\infty}\text{,}\quad
c(x\to\infty,y>0)=c_{\infty}\quad\text{and}\quad c(x=0,y>0)=c_{|\Sigma}.$

The analytical solution to this problem is

$c(x,y)=c_{|\Sigma}+(c_{\infty}-c_{|\Sigma})\cdot{\rm
erf}\left(\frac{x}{\delta(y)}\right)\qquad\text{with
}\delta(y)=\sqrt{4Dy/\nu}.$ (49)

The boundary layer thickness $\delta$ is considered a free parameter, which in
IT will be fitted for each interface attached cell individually as described
in section 4.3.2. The fitting selects the locally best approximation from the
family of given error-functions in eq. 49.

From the analytical solution for the fitted parameter, we can approximate a
more accurate boundary-normal gradient of the scalar concentration at the
interface as well as at the first cell faces normal to the boundary, $x_{f}$
in fig. 9. These gradients are used to scale the convective and diffusive
scalar fluxes in the Finite Volume discretization, keeping the discretization
of convective and diffusive operators implicit. Scaling is only applied in
interface normal direction and described in detail in sections 4.3.3 and
4.3.4.

Without the SGS, the maximally second-order unstructured Finite Volume
convective and diffusive fluxes are underestimated at the interface and
overestimated at the first cell faces normal to the boundary. Figure 9
visualizes the linear gradient calculation in red and the actual gradients at
the interface and the first cell faces normal to the boundary with the dashed
blue lines.

Figure 9: The concentration profile and linear/corrected gradient next to a
boundary. Adapted from [34].

Note that in the algorithm the concentration in the first cell is not altered
by the SGS. Instead the gradients are used to scale the diffusivities at the
respective locations, as described in sections 4.3.3 and 4.3.4.

Details of the SGS model can be found in [4, 34, 23]. It has been applied to
VoF by [34], where the implementation becomes more complicated.

###  4.3 Algorithm

####  4.3.1 Equation discretization in finite volume solvers

The species transport eq. 13 is discretised using the unstructured Finite
Volume Method, e.g. using the backward scheme as described above for Navier-
Stokes equations,

$\frac{3c_{P}^{n+1}V_{P}^{n+1}-4c_{P}^{n}V_{P}^{n}+c_{P}^{n-1}V_{P}^{n-1}}{\Delta
t}+\sum_{f}\Phi_{f}c_{f}^{n+1}=\sum_{f}D_{f}(\nabla
c)_{f}^{n+1}\cdot\boldsymbol{S}_{f},$ (50)

where $V_{P}$ is a control volume/cell,
$\Phi_{f}=\boldsymbol{S}_{f}\cdot(\mathbf{v}-\boldsymbol{w})_{f}$ is the face
flux and the superscripts ${}^{n+1},{}^{n}\text{ and }{}^{n-1}$ symbolize the
new, current and preceding timestep. The formula is written for a constant
time step size. The subscript ${}_{P}\text{ and }{}_{f}$ represent the value
at the cell centers, respectively the face centers.

For the SGS model to take effect in the computation, we will modify the
diffusivity coefficient and the face fluxes in the equation above.

####  4.3.2 Computation of model parameter boundary layer thickness

The model parameter boundary layer thickness $\delta$ is computed by
iteratively comparing the amount of the scalar in the first cell adjacent to
the boundary in the simulation, $\overline{\eta}_{c}$, with the amount that
would be contained in said cell for a given $\delta$, $\eta_{SGS}$. In fig. 9
both light brown areas have to be the same size or the area enclosed by the
real concentration profile and the cell centered FV-value has to be the same.
Mathematically this is being expressed in this expression:

$\displaystyle\overline{\eta}_{c}=\frac{\overline{c}-c_{|\Sigma}}{c_{\infty}-c_{|\Sigma}}\overset{!}{=}\frac{1}{V}\int_{V}\eta(x/\delta)\,dV=\eta_{SGS}$
(51) with
$\displaystyle\eta(x,y)=\frac{c(x,y)-c_{|\Sigma}}{c_{\infty}-c_{|\Sigma}}={\rm
erf}(x/\delta(y)).$ (52)

Using the derivative of $\eta_{SGS}$ with respect to $\delta$, as described in
[23], we use a Newton-bisection method until eq. 51 is fulfilled within a
given tolerance. For details like the initial value and the algorithm in
detail the reader is referred to the previously cited article.

For both, the advective and diffusive flux, the same idea is utilized: the
scaling of the fluxes according to the boundary layer thickness and therefore
the theoretical values given by the SGS.

####  4.3.3 Computation of diffusive fluxes

With $D$ being the molecular diffusivity, $S$ being the boundary (cell) area
and $\partial_{n}c$ being the surface normal gradient of the scalar c, the
diffusive flux $F^{D}$ is given by

$F^{D}=-D\cdot S\cdot\partial_{n}c.$ (53)

In the following, we will use the superscript SGS to indicate the value is
obtained from or modified by the SGS model. The SGS model solves the
analytical model for the actual surface normal gradient of $c$, which will be
$(\partial_{n}c)^{SGS}$. However, to stay implicit and increase numerical
stability, the diffusion coefficient is scaled. Therefore, we state that the
flux must be equal if computed with the actual gradient or a scaled diffusion
coefficient, meaning

$F^{D}=-D\cdot S\cdot(\partial_{n}c)^{SGS}\overset{!}{=}D^{SGS}\cdot
S\cdot\partial_{n}c$ (54a) which can be easily re-arranged into
$D^{SGS}=D\frac{\partial_{n}c}{(\partial_{n}c)^{SGS}},$ (54b)

where $D^{SGS}$ is a scaled diffusivity value that leads to the same flux, as
a calculation with $(\partial_{n}c)^{SGS}$ would. This scaling is applied to
the cell faces at the boundary as well as on the first cell faces normal to
the boundary. For both respective locations the diffusivity parameter needs to
be modified at the corresponding position, $x_{\Sigma}$ and $x_{f}$ in fig. 9.
There are exceptions, when the correction is not employed, the details of this
can be found in [23].

####  4.3.4 Computation of advective fluxes

In similar manner to the diffusion, the advective flux is given by

$F^{c}_{\\!f}:=c_{\\!f}^{SGS}F_{\\!f}=c_{\\!f}F_{\\!f}^{SGS}\Phi,$ (55)

where $F^{c}_{\\!f}$ is the volumetric flux of the concentration $c$ at the
face $f$, and $F_{\\!f}:=\mathbf{v}_{\\!f}\cdot\mathbf{S}_{f}$ is the
volumetric flux. Since no convection of the species normal to the interface
can take place, the scaling for advective fluxes is only applied to the first
cell faces normal to the boundary. Equation 55 analogously to our procedure
for the diffusive flux, leads to

$\Phi^{SGS}=\frac{c^{SGS}}{c}\Phi.$ (56)

Again, an exception handling takes places, of which the details can be found
in [23].

##  5 Software design

###  5.1 Software design of the SGS model library

The software design is closely related to the model as described in section 4.

Figure 10: Class diagram of the Subgrid-Scale model, also showing a potential
new model

The first benefit of our implementation is a strict adherence to the OpenFOAM
design principles, which enable easy library extension with new models
alternative to the SGS [23] and their Run-Time Selection (RTS) [17] using
configuration files. The implementation is depicted in fig. 10, already
showing how a new model newSubGridScaleModel would be implemented.

The second benefit of our software design is the possibility to apply SGS
models to any Partial Differential Equation (PDE) discretized in OpenFOAM,
without requiring extensive modification of existing solvers. Our SGS model
library can be applied to any part of the domain boundary with a minimal
requirement of a prismatic boundary-adjacent mesh with at least three layers
of cells in the boundary layer.

Our design of the SGS model implements an OpenFOAM library with a class
hierarchy and RTS [17]. A base class which does nothing is present, so an
option remains of not using SGS modeling. The SGS model is implemented in a
derived class, with virtual functions implementing the functionality from
[23]. All functionality needed by the SGS is either accessed from core
OpenFOAM libraries or is contained within the SGS model implementation, which
ensures the modularity of the SGS model. An alternative SGS model will require
similar addressing and geometric mesh information in the boundary layer and
the implementation provided in the original SGS model can be re-used.

To enable the user of the SGS a good understanding of the effects of the SGS
model and to make accessible insights for further developments of the SGS
model, we implemented the possibility to output various SGS model parameters
into OpenFOAM volScalarFields. These include the model parameter $\delta$ or
an indication of the exception handling, what represents the reason the SGS
model might not be employed. These can easily be inspected visually, an
example is given in fig. 16, where the diffusion coefficient after the SGS
scaling is visualized. One can see, that on the inflow side, the scaling is
strong and of opposite effect for the boundary faces versus the faces between
adjacent cell and the next outer cell layer. On the outflow side in contrast,
there is little effect to be seen. The third layer represents the un-modified
diffusion coefficient. Visualizing multiple operating figures of the SGS
enables insights into and deepens the understanding of the simulation case and
the effect of the SGS to it.

An exemplified code for the modification of the PDE in an OpenFOAM solver for
applying an SGS model from our model library is shown in listing
LABEL:lst:minSolver:

Listing 1: A minimal example of the application of the SGS model to a PDE
solver in OpenFOAM.

⬇

#include "baseSubGridScaleModel.H"

int main(int argc, char *argv[])

{

...

autoPtr<baseSubGridScaleModel> sgsPtr =

baseSubGridScaleModel::New(transportProperties.subDict("SGS"));

while (simple.loop())

{

Info<< "Time␣=␣" << runTime.timeName() << nl << endl;

while (simple.correctNonOrthogonal())

{

// Compute SGS data needed for scaling.

sgsPtr->update(psi);

fvScalarMatrix psiEqn

(

fvm::ddt(psi)

+ fvm::div(sgsPtr->scaledFlux(phi, psi), psi)

- fvm::laplacian(sgsPtr->scaledDiffCoeff(D, psi), psi)

==

fvOptions(psi)

);

psiEqn.relax();

fvOptions.constrain(psiEqn);

psiEqn.solve();

fvOptions.correct(psi);

}

runTime.write();

}

}

The SGS requires some values as user input, that are read from the
configuration file transportProperties in OpenFOAM, with the relevant SGS sub-
dictionary shown in listing LABEL:lst:transProp.

Listing 2: transportProperties sub-dictionary for the SGS model

⬇

SGS

{

// Choose the SGS model: "SubGridScale" or "inactive"

type SubGridScale;

// patch to apply SGS to

patchName freeSurface;

// far-field concentration

cInfinity cInfinity [0 -3 0 0 1 0 0] 0;

// (optional:) write parameters for visualization

// of SGS parameters

visualizeParameters true;

// (optional:) verbosity level for info statement output (0-3)

infoLvl 1;

}

The SGS model can be activated with the keyword-pair type SubGridScale for a
mesh boundary patch specified by the keyword patchName. Additional required
user input is the far-field concentration, defined with the keyword cInfinity.
Optional parameters are, whether to output additional parameters from the SGS
model as an OpenFOAM volScalarField and the level of verbosity from 0 to 3 to
print to the solver log. 0 corresponds to a quiet mode and 3 prints very
extensive information of almost all processes and results within the SGS
algorithm and calculations. As these settings affect simulation performance
and disk consumption, the first defaults to false while the later affects log
readability and defaults to 1, which gives limited log-output. These
visualizing capabilities as well the rapid info statement level enable
insights to the user, but they also help with future model development.

###  5.2 Automatic testing

We have set up our test cases in a Continuous Integration (CI) pipeline. The
pipeline is depicted in figs. 11 and 12. In our approach we follow the
practices presented in [18, 19]. This means the tests are run and evaluated
automatically in the CI with Jupyter Notebooks. They also produce the results
and graphs included in this work. A job in a succeeding stage to the
downstream pipeline has been omitted from fig. 12 for the sake of readability.
Said job collects and bundles all those artifacts from the two previous
stages.

Figure 11: Upstream CI pipeline Figure 12: Downstream CI pipeline

In [18, 19] CI is focused on tests and their automatic evaluation. We have
added tests to the CI pipeline, that check certain additional aspects of the
repository, like e.g. formatting and the formal correctness of the
CITATION.cff-file. We also build an image within the CI with OpenFOAM and
other needed packages. This job can not be seen in fig. 11, since it is only
run once a new OpenFAOM version is available. It is located before the trigger
job for the downstream pipelines.

We further expanded the CI to be able to test the code for multiple OpenFOAM
versions, which at the moment is being done for the latest three versions
OpenFOAM v2312, v2306 and v2212. In feature-branches only the latest version
is tested to save resources.

Another feature of the CI pipeline is the addition of a test, that checks the
OpenFOAM log-files for warnings and errors, where tutorials added to the
repository are automatically detected and run for a single timestep. This
almost ensures, that there are no errors in the set-up of test cases and would
also uncover some of errors in the code, that can only be detected, when
OpenFOAM is run. The automatisation of this test in finding and running all
simulation cases is especially valuable. Since running all of the simulations
for only one timestep needs a very limited amount of time it can be checked,
that changes in the code don’t corrupt any case before commiting and pushing
to the central repository. The automated finding of cases also prevents
forgetting to add a simulation case to the CI (to this testing method) and the
automated execution and error finding gently forces all participants of the
repository to set up their cases to run flawlessly and fully automatised, e.g.
they need to include everything needed to run a case in the Allrun script.
This ensures all the cases provided with the code can be run without
debugging.

Another topic not covered in the previously mentioned references is that for
repetitive tasks, like reading dictionary entries, reading files or
calculating error norms needed for the evaluation a set of python functions is
provided. This set of functions is also tested automatically with the testing
framework pytest as part of the CI pipeline and precedes the more
computationally expensive OpenFOAM-specific jobs, as can be seen in fig. 11.

##  6 Results

###  6.1 Flat plate test case

The smallest test for model development and code verification is what we call
the flat plate test case. It implements the assumptions made in section 4.2
and was already implemented by e.g. [23]. The test consists of a flat plate
with a prescribed constant concentration of
$1\text{\,}\mathrm{m}\mathrm{o}\mathrm{l}\mathrm{/}\mathrm{m}^{3}$. Parallel
to this plate a flow is prescribed, that has no interaction with the wall,
i.e. there is no velocity boundary layer or velocity gradient. Initially and
in the inflow, the concentration is set to zero, the length of the plate is
$0.005\text{\,}\mathrm{m}$ and the fluid velocity is
$0.1\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$. The timestep of
$0.0002\text{\,}\mathrm{s}$ is used for a physical duration of
$0.1\text{\,}\mathrm{s}$, which is double the time the fluid needs to be
transported through the domain. Figure 13 visualizes the setting.

Figure 13: Set-up of the flat plate test. Figure taken from [23].

With the possibilities described in fig. 10, it can be verified that
fundamental parts of the algorithm as described in section 4.3 are implemented
correctly, like the boundary layer thickness $\delta$ and the diffusion
correction. Beneficially, there is an analytical solution to this model for
comparison. For these reasons the test enables insights into the capabilities
of the SGS model and the implementation thereof. This simple test case does
circumvent most of the error-handling and also convective flux correction must
not have an effect on the outcome.

Figure 14: Results for the parameter study over of the flat plate test case

In fig. 14 the local Sherwood numbers of a parameter study along the boundary
surface are displayed. The parameter study has four levels of mesh refinement,
where the cell size halves between each execution from
$40\text{\,}\mu\mathrm{m}5\text{\,}\mu\mathrm{m}$, four different diffusivity
levels spanning four orders of magnitude from
$5\text{\times}{10}^{-8}\text{\,}\mathrm{m}^{2}\mathrm{/}\mathrm{s}5\text{\times}{10}^{-11}\text{\,}\mathrm{m}^{2}\mathrm{/}\mathrm{s}$
and every of these setups is run with and without active SGS model. Taking the
length of the plate as reference length $L$, the Peclet number is within the
interval $[10^{4},10^{7}]$.

It is apparent, how without the SGS model, the Sherwood number cannot be
computed in cases with a boundary layer that is small in comparison to the
cell size, i.e. when the cell size is big and the diffusivity is small. In
contrast to the linear interpolation with the SGS model the Sherwood number is
captured reasonably well in all cases. Note in fig. 14, that at some locations
the Sherwood number with the application of the SGS model is actually smaller
than without it.

###  6.2 Satapathy-Smith axisymmetric bubble or droplet

The second test case for the SGS model mimics a moving droplet or bubble. It
uses the axisymmetric solution of velocity field from Satapathy and Smith in
[24] for low Reynolds numbers. The computational domain is a wedge with
prescribed velocity. This test case was also implemented by and we take the
reference solution from [23]. This test case introduces additional complexity
compared to section 6.1, since there is flow non-parallel to the interface,
therefore violating some of the assumptions made when deriving the SGS
mathematical model, as described in section 4.2. From the visualization
capabilities we can see in the rear part of the ”bubble” the exception
handling partially taking place, as there the concentration rises, not
conforming to the assumption of a thin boundary layer (with parallel flow)
made in the derivation of the SGS. As in section 6.1 the inflow and initial
concentration is set to
$0\text{\,}\mathrm{m}\mathrm{o}\mathrm{l}\mathrm{/}\mathrm{m}^{3}$ and the
interface is set to have a constant concentration of
$1\text{\,}\mathrm{m}\mathrm{o}\mathrm{l}\mathrm{/}\mathrm{m}^{3}$. The
simulated physical duration is $0.6\text{\,}\mathrm{s}$ and the results shown
here represent the final state.

To see the best results with the SGS, it is important to set the interpolation
schemes to limitedLinear phi 1.0 for the concentration c and to linear for the
interpolation of the scaled diffusion coefficient as well as for the flux of U
in the fvSchemes-file.

We run a parameter study consisting of four different diffusivities from
$1\text{\times}{10}^{-8}\text{\,}\mathrm{m}^{2}\mathrm{/}\mathrm{s}1\text{\times}{10}^{-11}\text{\,}\mathrm{m}^{2}\mathrm{/}\mathrm{s}$,
four different mesh resolutions with the meshing parameter $\rm
N=[62,124,248,496]$ as the number of faces along the half-circled interface.
Every study is run with and without employing the SGS and with a constant
timestep of $0.0002\text{\,}\mathrm{s}$. The bubble diameter of
$2\text{\,}\mathrm{m}\mathrm{m}$ serves as reference length for the Sherwood
number. It can be observed from fig. 15 that without deploying the SGS, the
local Sherwood number is almost constant and except for the rear part of the
bubble too low in case of small diffusivities and coarse meshes. When the mesh
resolution is fine and the diffusivity is high, the boundary layer is not as
thin (in comparison to the mesh cell size), and therefore the results of
linear modeling and SGS modeling are similar.

Figure 15: Local Sherwodd numbers of a parameter study for Schmidt-numbers
from $1\text{\times}{10}^{4}\text{\,}1\text{\times}{10}^{7}\text{\,}$

In fig. 16 the diffusion coefficient after scaling with the SGS is displayed.
In the cell layer adjacent to the interface, the diffusion coefficient at the
interface is visible, whilst in the next outer cell layer the scaled diffusion
coefficient at the faces between interface-adjacent cells and that same cell
is visualized. The third row of cells represents the uncorrected diffusion
coefficient for reference.

The biggest correction is applied in the front part of the bubble, where the
flow impinges onto the interface and the boundary layer is very thin,
therefore the derivative in surface normal direction of the concentration is
very big. In the rear part of the bubble the boundary layer has grown and the
main direction of the flow is no longer parallel to the surface, making the
surface concentration gradient in interface-normal direction much smaller and
reasonable good represented by linear interpolation. Therefore less correction
is needed. To understand what is happening, we would like to point the reader
to fig. 9 again.

Figure 16: Visualization of corrected diffusion coefficients for the coarsest
mesh with a diffusivity of
$1\text{\times}{10}^{-9}\text{\,}\mathrm{m}^{2}\mathrm{/}\mathrm{s}$

###  6.3 Marangoni migration of an air bubble in silicone oil

The velocity of a bubble/droplet within an unbounded fluid subjected to
uniform temperature gradient ($\partial T/\partial z$) and zero gravity field
($g=0$) is defined by following expression [36]:

$v^{T}_{\rm
b,YBG}=\frac{2}{\left(2+\dfrac{\lambda_{i}}{\lambda_{o}}\right)\left(2+3\dfrac{\mu_{i}}{\mu_{o}}\right)}\frac{R}{\mu_{o}}\frac{\partial\sigma}{\partial
T}\frac{\partial T}{\partial z},$ (57)

where $R$ is the radius of the bubble, $T$ is the fluid temperature, $\sigma$
is the surface tension coefficient, $\mu_{i}$ is the dynamic viscosity of the
internal fluid, $\mu_{o}$ is the dynamic viscosity of the outside fluid and
$\lambda$ is the thermal conductivity of the fluid.

In the case the Marangoni convection is caused by a specified fixed gradient
of surfactant concentration along the interface, the bubble migration velocity
can be expressed based on eq. 57 as follows, see [20]:

$v^{\Gamma}_{\rm
b,YBG}=\frac{2}{3\left(2+3\dfrac{\mu_{i}}{\mu_{o}}\right)}\frac{R}{\mu_{o}}\frac{\partial\sigma}{\partial\Gamma}\frac{\partial\Gamma}{\partial
z},$ (58)

where $\Gamma$ is the concentration of surfactants at the interface.

A numerical simulation is performed for an air bubble of radius
$r_{b}=1~{}{\rm mm}$ moving through the silicone oil. Relevant properties of
both phases are listed in table 1. The bubble is assumed to be of spherical
shape and rigid and is fixed during the simulation. The spatial computational
domain consist of a volume of the spherical bubble filled with air and volume
around the bubble bounded by the spherical surface of radius $r=20r_{\rm b}$
and filled with silicone oil. The spatial domain is discretised by the
unstructured mesh consisting of $979^{\prime}520$ hexahedral cells.

The flow inside and outside of the bubble is driven by Marangoni forces due to
the nonuniform distribution of surfactant concentration along the interface,
which is fixed during the simulation. The distribution of surfactant along the
bubble surface is defined by following linear function:

$\frac{\Gamma}{\Gamma_{\infty}}=\frac{z+2r_{\rm b}}{L},$ (59)

where $\Gamma_{\infty}$ is the saturated surfactant concentration and $L$ is
the length scale, whose value is set to $525\text{\,}\mathrm{m}\mathrm{m}$.
The equation of state is also linear and reads as follows:

$\sigma=\sigma_{\rm s}\left(1-\beta_{\rm
s}\frac{\Gamma}{\Gamma_{\infty}}\right),$ (60)

where $\beta_{\rm s}=\frac{\mathcal{R}T\Gamma_{\infty}}{\sigma_{\rm s}}$ is
the elasticity number which amounts to $2$ in this study.

The unsteady computation is performed with time-step size $\Delta
t=$0.01\text{\,}\mathrm{s}$$ using the first order accurate Euler temporal
discretisation scheme, until steady state is reached. The gradient of the
surface tension coefficient as well as the gradients of volume fields is
calculated using the least square method. The diffusion terms in the momentum
and pressure equations are discretised using the skew-corrected central
differencing scheme, while the convective term in the momentum equation is
discretised using the linear-upwind scheme. At the outer side of the spatial
domain the velocity and normal derivative of pressure are set to zero and at
the interface corresponding kinematic and dynamics conditions are enforced.

The bubble rise velocity is calculated based on numerically determined volume
flow rates through the faces at the interface as follows:

$v_{\rm b}=\frac{\sum_{f}{\rm pos}(\phi_{\rm i,f})\phi_{\rm i,f}}{\sum_{f}{\rm
pos}(\mathbf{k}\cdot\mathbf{S}_{\rm i,f})(\mathbf{k}\cdot\mathbf{S}_{\rm
i,f})},$ (61)

where $\phi_{\rm i,f}$ is the volume flow rate of the fluid through the face
$f$ at the interface, $\mathbf{S}_{\rm i,f}$ is the area vector of the face
$f$ at the interface and $\mathbf{k}$ is the unit direction vector of bubble
motion. The bubble velocity calculated using eq. 61 based on simulation
results amounts to
$v_{b}=$0.001\,325\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$$, while the
corresponding analytical (eq. 58) gives the migration velocity as
$v^{\Gamma}_{\rm
b,YBG}=$0.001\,329\,5\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$$. The velocity
field inside and around the bubble along the plane $y=0$ is shown in fig. 17.

| $\rho,{\rm kg/m^{3}}$ | $\mu,{\rm Pas}$ | $\sigma,{\rm N/m}$
---|---|---|---
Air | $1.8\times 10^{-5}$ | $1.2$ | $0.02$
Oil | $955$ | $0.191$ | $0.02$
Table 1: Properties of air and silicone oil. Figure 17: Velocity field inside
and around the bubble along the plane $y=0$ for the bubble marangoni migration
case.

###  6.4 Transport with the interface

We test the framework’s ability of adequate transport of a species along an
interface with a rotating droplet with a prescribed concentration. The droplet
is represented in 2-D. The test setup follows the test done by [3, 2]. We
prescribe the initial concentration distribution and the velocity. For the
evaluation we compute one rotation, what equals one second with our prescribed
velocity referring to an angular velocity $\omega=1/s$.

The initial distribution is given by

$\Gamma(\theta,\tau=0)=\frac{\cos{\theta}+1}{2}\Gamma_{0},$ (62)

where $\Gamma$ is the surface concentration, $\tau$ is the time, $\theta$ is
the angle and $\Gamma_{0}$ is the maximum initial concentration. After a time
$\tau$ the expected distribution from an analytical approach is given by

$\Gamma(\theta,\tau)=\frac{e^{-D_{I}\tau}\cos(\theta-\omega\tau)+1}{2}\Gamma_{0},$
(63)

where $D_{I}$ is the diffusion in interface tangential direction.

Figure 18: Results from the rotating droplet parameter study

For proper evaluation a mesh refinement study with a parameter variation of
the diffusivity is conducted. The results are visualized in fig. 18. We use a
meshing parameter $N=[10,20,30]$, which correspond to $40$, $80$ and $160$
surface cells along the interface. The timestep is set to
$0.0005\text{\,}\mathrm{s}$ regardless of the mesh size. The surface
diffusivity is varied between
$0\text{\,}\mathrm{m}^{2}\mathrm{/}\mathrm{s}0.1\text{\,}\mathrm{m}^{2}\mathrm{/}\mathrm{s}$.

(a) Error norms

(b) Rate of convergences

Figure 19: Errors of the rotating droplet parameter study

Visually it is hard to see a difference between the discretizations in fig.
18. In fig. 19 the error norms and their rate of convergence are shown. The
infinity norm and L1-norm are close to $1$ and the L2-norm is around $1.5$ for
the resolutions analyzed. With the low error norms and their respective rate
of convergence we deduct, that our framework is well able to handle transport
along an interface.

###  6.5 Expanding fluid interface area

To test the frameworks conservation of surfactant mass for an expanding fluid
interface area, we implement a test case, that mimics an expanding droplet
like it was done by [3, 2]. The initial distribution of species is described
and there is no diffusion happening. The initial distribution is the same as
in eq. 62. The droplet is represented by a wedge that is inflated with a
prescribed velocity and thereby the surface is stretched. The analytical
solution to this setup is given by

$\Gamma(\theta,r(\tau))=\frac{\cos{\theta}+1}{2}\Gamma_{0}\frac{r_{0}^{2}}{r(\tau)^{2}},$
(64)

where $r(\tau)=\sqrt[3]{3V_{0}(\tau+\tau_{0})}$, where $\tau_{0}=1/3s$ and $r$
is the radius, respectively $r_{0}$ is the initial radius.

Figure 20 shows, how the species concentration per area at the interface has
decreased after $2.5\text{\,}\mathrm{s}$ of bubble growth. The concentration
precisely matches the expected result for all mesh resolutions tested. We
conclude, that our framework is fit to simulate phenomena involving surface
expansion/shrinking.

Figure 20: Surface concentration after $0.5\text{\,}\mathrm{s}$ for different
mesh resolutions

##  7 Conclusions

We provide a thoroughly automatically tested implementation of the
unstructured Finite-Volume Arbitrary Lagrangian / Eulerian (ALE) Interface-
Tracking (IT) method for simulating incompressible, immiscible two-phase flows
as an OpenFOAM module.

In conclusion, we have successfully developed a Subgrid-Scale (SGS) model
hierarchy within an Arbitrary Lagrangian-Eulerian (ALE) Interface Tracking
OpenFOAM module, with configurable output parameters for in-depth analysis.
The SGS model’s ability to accurately capture thin boundary layers and steep
gradients was demonstrated, alongside the critical need for flux correction
not only at the interfaces but also across adjacent cell boundaries.
Furthermore, the Interface Tracking module’s proficiency in accurately
simulating surfactants on moving fluid interfaces interfaces and their impact
on surface tension and Marangoni flow was demonstrated.

##  8 Acknowledgments

Moritz Schwarzmeier and Tomislav Marić would like to thank the Federal
Government and the Heads of Government of the Länder, as well as the Joint
Science Conference (GWK), for their funding and support within the framework
of the NFDI4Ing consortium. Funded by the German Research Foundation (DFG) -
project number 442146713.

Funded by the German Research Foundation (DFG) – Project-ID 265191195 – SFB
1194.

Part of this work was funded by the Hessian Ministry of Higher Education,
Research, Science and the Arts - cluster project Clean Circles.

## References

  * Alke et al. [2010] A. Alke, D. Bothe, M. Kröger, B. Weigand, Dominik Weirich, and H. Weking.  Direct numerical simulation of high schmidt number mass transfer from air bubbles rising in liquids using the volume-of-fluid-method.  _Ercoftac Bulletin_ , 82:5–10, 2010.
  * Antritter et al. [2024] T. Antritter, T. Josyula, T. Marić, D. Bothe, P. Hachmann, B. Buck, T. Gambaryan-Roisman, and P. Stephan.  A two-field formulation for surfactant transport within the algebraic volume of fluid method.  _Computers & Fluids_, page 106231, March 2024.  ISSN 0045-7930.  doi: 10.1016/j.compfluid.2024.106231.
  * Antritter [2022] Thomas Antritter.  _Numerical Simulation of Coupled Wetting and Transport Phenomena in Inkjet Printing_.  PhD thesis, Technische Universität Darmstadt, Darmstadt, 2022.  doi: 10.26083/tuprints-00021326.
  * Bothe and Fleckenstein [2013] Dieter Bothe and Stefan Fleckenstein.  A volume-of-fluid-based method for mass transfer processes at fluid particles.  _Chemical Engineering Science_ , 101:283–302, 2013.  ISSN 0009-2509.  doi: 10.1016/j.ces.2013.05.029.
  * Brenner [2013] H. Brenner.  _Interfacial Transport Processes and Rheology_.  Butterworth-Heinemann series in chemical engineering. Elsevier Science, 2013.  ISBN 9781483292274.
  * Chen et al. [2000] K. P. Chen, W. Saric, and H. A. Stone.  On the deviatoric normal stress on a slip surface.  _Physics of fluids A_ , 12(12):3280–3281, 2000\.
  * Demirdžić and Muzaferija [1995] I. Demirdžić and S. Muzaferija.  Numerical method for coupled fluid flow, heat transfer and stress analysis using unstructured moving meshes with cells of arbitrary topology.  _Computer methods in applied mechanics and engineering_ , 125:235–255, 1995.
  * Demirdžić and Perić [1988] I. Demirdžić and M. Perić.  Space conservation law in finite volume calculations of fluid flow.  _International journal for numerical methods in fluids_ , 8:1037–1050, 1988.
  * Dieter-Kissling et al. [2014] Kathrin Dieter-Kissling, Mohsen Karbaschi, Holger Marschaldenl, Aliyar Javadi, Reinhard Miller, and Dieter Bothe.  On the applicability of drop profile analysis tensiometry at high flow rates using an interface tracking method.  _Colloids and Surfaces A: Physicochemical and Engineering Aspects_ , 441:837–845, 2014.  doi: 10.1016/j.colsurfa.2012.10.047.
  * Dieter-Kissling et al. [2015] Kathrin Dieter-Kissling, Holger Marschall, and Dieter Bothe.  Direct numerical simulation of droplet formation processes under the influence of soluble surfactant mixtures.  _Computers & Fluids_, 113:93–105, 2015.  doi: 10.1016/j.compfluid.2015.01.017.
  * Ferziger and Perić [1995] J. H. Ferziger and M. Perić.  _Computational methods for fluid dynamics_.  Springer Verlag, Berlin-New York, 1995.
  * Grosso et al. [2024] Mathis Grosso, Guillaume Bois, and Adrien Toutant.  Thermal boundary layer modelling for heat flux prediction of bubbles at saturation: A priori analysis based on fully-resolved simulations.  _International Journal of Heat and Mass Transfer_ , 222:124980, May 2024.  ISSN 00179310.  doi: 10.1016/j.ijheatmasstransfer.2023.124980.
  * Gründing et al. [2016] Dirk Gründing, Stefan Fleckenstein, and Dieter Bothe.  A subgrid-scale model for reactive concentration boundary layers for 3D mass transfer simulations with deformable fluid interfaces.  _International Journal of Heat and Mass Transfer_ , 101:476–487, October 2016.  ISSN 0017-9310.  doi: 10.1016/j.ijheatmasstransfer.2016.04.119.
  * Jasak [1996] H. Jasak.  _Error analysis and estimation for finite volume method with applications to fluid flows_.  PhD thesis, Imperial College, University of London, 1996.
  * John C. Slattery et al. [2007] John C. Slattery, Leonard Sagis, and Eun-Suok Oh.  _Interfacial Transport Phenomena_.  Springer US, Boston, MA, 2007.  ISBN 978-0-387-38438-2.  doi: 10.1007/978-0-387-38442-9.
  * Kralchevsky et al. [2002] Peter A. Kralchevsky, Krassimir D. Danov, and Nikolai D. Denkov.  Chemical Physics of Colloid Systems and Interfaces.  In K. S. Birdi, editor, _Handbook of Surface and Colloid Chemistry_ , page 208. CRC Press, New York, second expanded and updated edition edition, 2002.  ISBN 978-0-429-12509-6.
  * Marić et al. [2021] Tomislav Marić, Jens Höpken, and Kyle G. Mooney.  _The OpenFOAM technology primer_.  March 2021.  doi: 10.5281/zenodo.4630596.  Version OpenFOAM-v2012.1.
  * Marić et al. [2022] Tomislav Marić, Dennis Gläser, Jan-Patrick Lehr, Ioannis Papagiannidis, Benjamin Lambie, Christian Bischof, and Dieter Bothe.  A Research Software Engineering Workflow for Computational Science and Engineering.  August 2022.  doi: 10.48550/arXiv.2208.07460.
  * Marić et al. [2023] Tomislav Marić, Dennis Gläser, Jan-Patrick Lehr, Ioannis Papagiannidis, Benjamin Lambie, Christian Bischof, and Dieter Bothe.  A pragmatic workflow for research software engineering in computational science.  2023\.  doi: 10.48550/ARXIV.2310.00960.
  * Muradoglu and Tryggvason [2008] Metin Muradoglu and Gretar Tryggvason.  A front-tracking method for computation of interfacial flows with soluble surfactants.  _Journal of Computational Physics_ , 227(4):2238–2262, February 2008.  doi: 10.1016/j.jcp.2007.10.003.
  * Muzaferija and Perić [1997] S. Muzaferija and M. Perić.  Computation of free-surface flows using the finite-volume method and moving grids.  _Numerical heat transfer, Part B_ , 32:369–384, 1997.
  * Pesci et al. [2017] Chiara Pesci, Holger Marschall, Talmira Kairaliyeva, Vamseekrishna Ulaganathan, Reinhard Miller, and Dieter Bothe.  Experimental and Computational Analysis of Fluid Interfaces Influenced by Soluble Surfactant.  In Dieter Bothe and Arnold Reusken, editors, _Transport Processes at Fluidic Interfaces_ , pages 395–444. Springer International Publishing, Cham, 2017.  ISBN 978-3-319-56602-3.  doi: 10.1007/978-3-319-56602-3˙15.
  * Pesci et al. [2018] Chiara Pesci, Andre Weiner, Holger Marschall, and Dieter Bothe.  Computational analysis of single rising bubbles influenced by soluble surfactant.  _Journal of Fluid Mechanics_ , 856:709–763, 2018.  ISSN 0022-1120, 1469-7645.  doi: 10.1017/jfm.2018.723.
  * Satapathy and Smith [1961] R. Satapathy and W. Smith.  The motion of single immiscible drops through a liquid.  _Journal of Fluid Mechanics_ , 10(4):561, 1961\.  ISSN 0022-1120, 1469-7645.  doi: 10.1017/S0022112061000366.
  * Schwarzmeier et al. [2024] M. Schwarzmeier, S. Raju, Z. Tuković, and T. Marić.  Twophaseintertrackfoam source code repository.  GitLab, 2024.  URL https://gitlab.com/interface-tracking/twophaseintertrackfoamrelease/-/tree/v1.0?ref_type=tags.  Accessed: 2024-03-24.
  * Schwarzmeier et al. [2024-03-24] Moritz Schwarzmeier, Suraj Raju, Zeljko Tuković, and Tomislav Marić.  twophaseintertrackfoam - data & results v1.0, 2024-03-24.
  * Thomas and Lombard [1979] P. D. Thomas and C. K. Lombard.  Geometric conservation law and its application to flow computations on movining grids.  _AIAA Journal_ , 17:1030–1037, 1979.
  * Tuković and Jasak [2008] Željko Tuković and Hrvoje Jasak.  Simulation of free-rising bubble with soluble surfactant using moving mesh finite volume/area method.  In _Proceedings of 6th International Conference on CFD in Oil & Gas, Metallurgical and Process Industries, no. CFD08-072_, 2008.
  * Tuković et al. [2018] Željko Tuković, Milovan Perić, and Hrvoje Jasak.  Consistent second-order time-accurate non-iterative PISO-algorithm.  _Computers & Fluids_, 166:78–85, April 2018.  doi: 10.1016/j.compfluid.2018.01.041.
  * Tuković et al. [2024-03-24] Zeljko Tuković, Moritz Schwarzmeier, Suraj Raju, and Tomislav Marić.  twophaseintertrackfoam - code, 2024-03-24.  doi: 10.48328/tudatalib-1395.
  * Željko Tuković and Jasak [2012] Željko Tuković and Hrvoje Jasak.  A moving mesh finite volume interface tracking method for surface tension dominated interfacial fluid flow.  _Computers & Fluids_, 55:70–84, February 2012.  doi: 10.1016/j.compfluid.2011.11.003.
  * Weatherburn [1972] C. E. Weatherburn.  _Differential geometry in three dimension_.  Cambridge University Press, London, 1972.
  * Weber et al. [2017] Paul S Weber, Holger Marschall, and Dieter Bothe.  Highly accurate two-phase species transfer based on ale interface tracking.  _International Journal of Heat and Mass Transfer_ , 104:759–773, 2017.
  * Weiner and Bothe [2017] Andre Weiner and Dieter Bothe.  Advanced subgrid-scale modeling for convection-dominated species transport at fluid interfaces with application to mass transfer from rising bubbles.  _Journal of Computational Physics_ , 347:261–289, 2017\.  ISSN 00219991.  doi: 10.1016/j.jcp.2017.06.040.
  * Weiner et al. [2021] Andre Weiner, Dirk Gründing, and Dieter Bothe.  Computing Mass Transfer at Deformable Bubbles for High Schmidt Numbers.  _Chemie Ingenieur Technik_ , 93(1-2):81–90, 2021\.  ISSN 1522-2640.  doi: 10.1002/cite.202000214.
  * Young et al. [1959] N. O. Young, J. S. Goldstein, and M. J. Block.  The motion of bubbles in a vertical temperature gradient.  _Journal of Fluid Mechanics_ , 6(03):350, October 1959.  doi: 10.1017/s0022112059000684.