NIETHAMMER et al

*Matthias Niethammer, Mathematical Modeling and Analysis, Technische Universität Darmstadt, Peter-Grünberg-Straße 10, 64287 Darmstadt, Germany.

# An unstructured geometrical un-split VOF method for viscoelastic two-phase
flows

Matthias Niethammer  Muhammad Hassan Asghar  Tomislav Maric  Dieter Bothe
Mathematical Modeling and Analysis, Technische Universität Darmstadt,
Darmstadt, Germany  Simulation of reactive Thermo-Fluid Systems, Technische
Universität Darmstadt, Darmstadt, Germany<EMAIL_ADDRESS>
(26 April 2016; 6 June 2016; 6 June 2016)

###### Abstract

[Summary] We present an unstructured geometrical Volume-of-Fluid (VOF) method
for handling two-phase flows with a viscoelastic liquid phase. The
viscoelastic behavior is modeled via generic rate-type constitutive equations.
A one-field formulation is employed, which results from conditional volume
averaging of the local instantaneous bulk equations and interface jump
conditions. The method builds on the ‘plicRDF-isoAdvector’ geometrical VOF
solver that is extended and combined with the modular framework ‘DeboRheo’ for
viscoelastic CFD. A piecewise-linear geometrical interface reconstruction
technique on general unstructured meshes is employed for discretizing the
viscoelastic stresses across the fluid interface in a numerically consistent
and highly accurate way. Because of the numerical challenges posed by the high
Weissenberg number problem, we employ an appropriate stabilization approach to
the constitutive equation of the viscoelastic phase to increase the robustness
of the method at higher fluid elasticity. DeboRheo facilitates a flexible
combination of different rheological models with appropriate stabilization
methods to address the high Weissenberg number problem. We discuss the
theoretical formulation and implementation of the proposed method and
demonstrate its effectiveness using numerical examples of viscoelastic flows.
The results highlight the method’s ability to accurately capture the behavior
of viscoelastic fluids in various applications. The proposed method holds
promise for furthering our understanding and predictive capabilities of
viscoelastic flows in various industrial and natural processes.

###### keywords:

viscoelastic two-phase flows, direct numerical simulation (DNS), geometric
Volume-of-Fluid (VOF) method, unstructured Finite Volume method, high
Weissenberg number stabilization

††articletype: Research Article

##  1 Introduction

Two-phase flows with at least one viscoelastic phase are encountered in a wide
range of industrial and natural processes, including pharmaceutical
manufacturing, plastics technology, and food processing. One example is
laminar and dispersive mixing operations, such as industrial emulsification
processes, where well-defined shear and elongational flows subject large drops
to deformation and fragmentation, resulting in the formation of smaller drops
1, 2. A thorough understanding and accurate simulation of such viscoelastic
two-phase flows are crucial for adjusting product quality, ensuring safety
properties, and optimizing industrial mixing processes.

A fundamental issue in emulsion technology is the study of a single drop’s
deformation under shear 3, 4, 5, 6, as seen in shear mixing processes where
large drops deform and fragment. This phenomenon is critical in predicting
drop size distributions and rheological properties of emulsions 7. These
properties hold significant implications on various industrial applications,
such as the manufacturing of pharmaceuticals, where they profoundly impact the
performance and stability of the final products. Similarly, in plastics
technology, the properties of paints and varnishes are crucially determined by
the deformation and fragmentation behavior of droplets in the emulsion.

In Direct numerical simulations (DNS) of viscoelastic two-phase flows major
challenges arise from the abrupt change of physical properties and rheological
behaviors of the two phases at the fluid interface, as well as the proper
incorporation of viscoelastic constitutive equations into the computational
framework. Additionally, the presence of high Weissenberg numbers, associated
with the elastic nature of the viscoelastic phase, introduces additional
numerical complexities 8, 9 that demand specific numerical treatment 10.
Furthermore, the deformation of a drop in simple shear flow is characterized
by strongly deformed fluid interfaces and potential drop breakup, which can
impact the stability and computational efficiency of DNS. Under these
conditions, the Volume-of-Fluid (VOF) 11 method emerges as an appropriate
numerical modeling approach that has been previously employed with success in
addressing viscoelastic drop deformation in both 2D 12 and 3D shear flows 13,
14, 15, 16, 17. Moreover, a generic closure for the viscoelastic constitutive
equations, supporting various rheological models and stabilization techniques
for high Weissenberg numbers, was systematically derived in 18 through
conditional volume averaging 19, 20, 21. This framework has demonstrated the
ability to provide high-fidelity predictions of viscoelastic two-phase flow
phenomena 22 and serves as the foundation for the developments addressed in
the present study. The current work extends the VOF framework for complex
fluids ‘DeboRheo’18, 23 with the ‘plicRDF-isoAdvector’ geometrical VOF method
24, 25 from the TwoPhaseFlow library 26, 27. The result of this integration is
a new unstructured, dimensionally un-split VOF approach specifically designed
for viscoelastic two-phase flows which is introduced and validated below.

The VOF method is an interface capturing method widely employed for direct
numerical simulations of two-phase flows, because of its inherent global and
local volume conservation and ability to handle strong deformations,
coalescence and breakup of the fluid interface. Additionally, the VOF method
can be applied to unstructured domain discretization and provides a
straightforward parallel computation model. In the VOF method, the sharp
interface between phases is captured by a phase indicator function. The
governing equations are written down as a one-fluid formulation, utilizing a
single set of equations for both phases throughout the entire domain and
incorporating appropriate source terms to account for the fluid interface 28.

Two types of VOF methods can be categorized based on the specific advection
scheme utilized to propagate the sharp interface:

Algebraic VOF methods 29, 30, 31 advect the volumetric phase fraction
algebraically through numerical schemes by solving a scalar transport
equation. Over the last decades, algebraic VOF methods have seen widespread
application on general unstructured grids due to their flexibility and high
level of serial and parallel computational efficiency. Despite these
advantages, the advection schemes of algebraic VOF methods struggle to keep
the interface sharp and free of oscillations. The abrupt change of the volume
fraction field at the interface causes large interpolation errors and,
subsequently, artificial diffusion.

Geometric VOF methods32 utilize geometrical operations to approximate the
solution of the volume fraction advection equation, rather than algebraic
calculations. Contemporary geometrical VOF methods reconstruct the interface
as a piecewise-planar surface within multi-material cells. Subsequently, the
reconstructed interface is advected by computing the phase-specific volumes
that are fluxed across the boundary of the finite volume. Un-split geometrical
VOF methods32 facilitate the use of unstructured meshes to manage complex
geometrical domains and improve the overall solution accuracy, which motivated
the extension of the ‘DeboRheo’ computational framework with the ‘plicRDF-
isoAdvector’ geometrical VOF method 24.

The objective of this paper is threefold: (i) to derive an extension of the
geometrical VOF for viscoelastic fluids; (ii) to demonstrate the numerical
stability and robustness of the method at high Weissenberg number for
technically relevant geometries; (iii) to study the deformation of a single
drop in shear flow using transient simulations in both 2D and 3D, both for
scenarios involving a viscoelastic drop in a Newtonian matrix as well as a
Newtonian drop in a viscoelastic matrix.

##  2 Mathematical model

The extended Volume-of-Fluid (VOF) equations for viscoelastic fluids can be
derived rigorously from the local instantaneous conservation equations in both
bulk phases $\Omega^{\pm}(t)$ and the interfacial jump conditions on the
interface $\Sigma(t)$ by Conditional Volume Averaging (CVA). Readers
interested in a comprehensive derivation of the extended VOF equations can
refer to Niethammer (2019) 18. The outcome of this derivation is a one-field
formulation, where a single set of equations is valid for both phases across
the entire domain. The one-field formulation is summarized below under the
assumptions of laminar, isothermal flows with constant mass densities in the
bulk phases.

A physical domain, represented by $\Omega$, is cut by a two-dimensional smooth
surface, $\Sigma(t)$, into two distinct subdomains (bulk phases)
$\Omega^{\pm}(t)$, as depicted in Figure 1. These sub-domains are filled with
different types of incompressible fluids. The interface unit normal pointing
towards the phase $\Omega^{-}(t)$ is denoted as $\mathbf{n}_{\Sigma}$. The
phase-indicator function,

$\chi(\textbf{x},t)=\begin{cases}1&\textbf{x}\in\Omega^{+}(t)\\\ 0&\textbf{x}\
\textnormal{otherwise},\end{cases}$ (1)

is utilized to distinguish the sub-domains $\Omega^{\pm}$. The sub-domain
$\Omega^{+}$ can therefore be written as

$\Omega^{+}(t):=\\{\textbf{x}\in\Omega:\chi(\textbf{x},t)=1\\}.$ (2)

$\Omega$
---
$\Omega^{+}(t)$
---
$\Omega^{-}(t)$
---
$\Sigma(t)$
---
$\textbf{n}_{\Sigma}$
---
$\partial\Omega$
---
$\chi(t,\textbf{x})=1$
---
$\chi(t,\textbf{x})=0$
---
Figure 1: Schematic diagram of a twophase domain $\Omega$.

The phase-indicator function describes the distribution of the bulk phases in
$\Omega$. Although not defined on the interface $\Sigma(t)$, its left and
right-sided limits exist there.

In CVA, the phase-indicator function is applied to the local instantaneous
bulk equations for mass and momentum prior to performing volume averaging over
the selected volume. This technique requires the definition of volume averaged
quantities, which are introduced below. The volumetric phase fraction, or
volume fraction, is defined as the volume average of the phase indicator
function

$\alpha(\textbf{x},t):=\frac{1}{|V|}\int_{V}\chi(\textbf{x}+\textbf{y},t)\,d\textbf{y}\,.$
(3)

In this definition the centroid of the averaging volume is indicated by the
position vector x, while the integration over the volume $V$ is performed in
terms of the components of the relative position vector y. The volume fraction
describes the ratio of the volume $\Omega^{+}(t)$ occupied by the phase $+$ in
the volume $(\textbf{x}+V)$ at time $t$ to the total volume of $V$ and is
constrained by $0\leq\alpha\leq 1$. Alternatively, the volume fraction can be
defined as the average

$\alpha(\textbf{x},t):=\frac{1}{|V|}\int_{(\textbf{x}+V)\,\cap\,\Omega^{+}(t)}1\,d\textbf{y}\,=\frac{\left|(\textbf{x}+V)\,\cap\,\Omega^{+}(t)\right|}{|V|}$
(4)

over the phase-specific volume $(\textbf{x}+V)\,\cap\,\Omega^{+}(t)$. From (4)
it follows for a scalar or tensor variable $\mathbf{\Phi}(\textbf{x},t)$ that

$\frac{1}{|V|}\int_{V}\chi(\textbf{x}+\textbf{y},t)\,\mathbf{\Phi}(\textbf{x}+\textbf{y},t)\,d\textbf{y}\,=\frac{|(\textbf{x}+V)\,\cap\,\Omega^{+}(t)|}{|V|}\frac{1}{|(\textbf{x}+V)\,\cap\,\Omega^{+}(t)|}\int_{(\textbf{x}+V)\,\cap\,\Omega^{+}(t)}\mathbf{\Phi}(\textbf{y},t)\,d\textbf{y}\,=\alpha\,\overline{\mathbf{\Phi}(\textbf{x},t)}^{+}\,,$
(5)

where $\overline{\mathbf{\Phi}(\textbf{x},t)}^{+}$ denotes the intrinsic phase
average over the volume $(\textbf{x}+V)\,\cap\,\Omega^{+}(t)$. This phase
average can be defined for both phases

$\overline{\mathbf{\Phi}(\textbf{x},t)}^{\pm}=\frac{1}{|(\textbf{x}+V)\,\cap\,\Omega^{\pm}(t)|}\int_{(\textbf{x}+V)\,\cap\,\Omega^{\pm}(t)}\mathbf{\Phi}(\textbf{y},t)\,d\textbf{y}\,.$
(6)

With these definitions, the volume averaged mass and momentum conservation
equations can be written as

$\displaystyle\nabla{\scriptscriptstyle\mathchoice{\mskip 3.0mu\lower
0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip 3.0mu\lower
0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}{\mathbf{u}}$ $\displaystyle=0,$ (7)
$\displaystyle\partial_{t}{\rho}{\mathbf{u}}+({\mathbf{u}}{\scriptscriptstyle\mathchoice{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\nabla){\rho}{\mathbf{u}}$
$\displaystyle=-{\nabla p}+{\nabla{\scriptscriptstyle\mathchoice{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}{\boldsymbol{\tau}}}+{\rho}\mathbf{b}+\mathbf{f}_{\sigma},$
(8)

where $p$ is the pressure, $\mathbf{b}$ represents external body forces and
$\mathbf{f}_{\sigma}$ denotes the surface tension forces. Note that the one-
field representation (7) and (8) is formulated in terms of mixture quantities,
namely, the mixture density

$\rho=\alpha\rho^{+}+(1-\alpha)\rho^{-}\,,$ (9)

the mixture velocity

$\mathbf{u}=\frac{\alpha\rho^{+}\overline{\mathbf{u}}^{+}\,+\,(1-\alpha)\rho^{-}\overline{\mathbf{u}}^{-}}{\rho}\,,$
(10)

and the mixture extra stress

${\boldsymbol{\tau}}={\alpha\overline{\boldsymbol{\tau}}^{+}\,+\,(1-\alpha)\overline{\boldsymbol{\tau}}^{-}}\,.$
(11)

Note also that the term $f_{\sigma}$ in (8) originally is a surface averaged
surface tension force as this results from the CVA process together with the
jump conditions at the interface. But in the numerical treatment we use the
CSF model33 for the surface tension force which corresponds to
$\mathbf{f}_{\sigma}=\sigma\kappa\nabla\alpha$, where $\sigma$ denotes the
surface tension and $\kappa$ represents the curvature.

Additional closures are required for the stress term on the right-hand side of
the momentum equation (8).For the extra stress ${\boldsymbol{\tau}}$, we
employ the solvent-polymer stress splitting, i.e.

$\displaystyle{\boldsymbol{\tau}}={\boldsymbol{\tau}_{s}}+{\boldsymbol{\tau}_{p}},$
(12)

where ${\boldsymbol{\tau}_{s}}$ is the solvent stress and
${\boldsymbol{\tau}_{p}}$ corresponds to the polymer stress. The viscous
solvent stress obeys Newton’s constitutive relation, given by

$\displaystyle{\boldsymbol{\tau}_{s}}$
$\displaystyle=2{{\eta}_{s}}\mathbf{D}={{\eta}_{s}}\big{(}\nabla{{\mathbf{u}}}+\nabla{{\mathbf{u}}}\\!^{\top}\big{)},$
(13)

with

${{\eta}_{s}}=\alpha{\eta}^{+}_{s}+(1-\alpha){\eta}^{-}_{s}\,,$ (14)

denoting the mixture solvent viscosity. The equations (13) and (14) constitute
a complete closure for the solvent stress contribution in (12). Nonetheless,
considering that the stress appears in a stress divergence term in the
momentum equation (8), additional manipulations can enhance the numerical
evaluation of this specific term. Incorporating (13) and (7) into (8), the
solvent stress divergence can be reformulated according to

$\displaystyle{\nabla{\scriptscriptstyle\mathchoice{\mskip 3.0mu\lower
0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip 3.0mu\lower
0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\boldsymbol{\tau}_{s}}$
$\displaystyle=\nabla{\scriptscriptstyle\mathchoice{\mskip 3.0mu\lower
0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip 3.0mu\lower
0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\left({{\eta}_{s}}\nabla{{\mathbf{u}}}\right)+\nabla{\scriptscriptstyle\mathchoice{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\left({{\eta}_{s}}(\nabla{{\mathbf{u}}})\\!^{\top}\right)$
(15) $\displaystyle=\nabla{\scriptscriptstyle\mathchoice{\mskip 3.0mu\lower
0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip 3.0mu\lower
0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\left({{\eta}_{s}}\nabla{{\mathbf{u}}}\right)+\nabla{{\mathbf{u}}}{\scriptscriptstyle\mathchoice{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\nabla{{\eta}_{s}},$ (16)

where the elliptic (Laplacian-type) term in (16) is discretized implicitly
within the present work’s finite volume framework. Furthermore, it is
important to note that the discretization of (16) requires an evaluation of
the mixture viscosity ${{\eta}_{s}}$ at cell faces using the phase indicator
approximated by the geometrical VOF method, and it is necessary to ensure
numerical stability for phases with strongly different densities 34. This
scheme enables the computation of (16) from the face-reconstructed volume
fractions instead of the standard interpolation of cell-centered values to
faces. Incorporating the PLIC reconstruction scheme in the evaluation of the
mixture solvent viscosity at the cell-faces ensures a compact computational
stencil, reducing potential sources of interpolation errors. This numerical
procedure is described in section 4.

The polymer stress $\boldsymbol{\tau}_{p}$ in (8) can be represented as a
linear function of the polymer conformation tensor $\mathbf{C}$ 35, 36. This
relationship can be expressed as

$\displaystyle\boldsymbol{\tau}_{p}=G\left(\mathbf{C}-\mathbf{I}\right),$ (17)

where the mixture modulus $G$ is defined as $G=\alpha
G^{{{+}}}+\left(1-\alpha\right)G^{{{-}}}$ with the phase-specific constants
$G^{\pm}=\frac{\eta_{p}^{\pm}}{\lambda^{\pm}}$ and $\mathbf{I}$ is the unit
tensor. Here, ${\eta}_{p}$ is the polymer viscosity and $\lambda$ is the
relaxation time. The mixture conformation tensor $\mathbf{C}$ in (17) is
defined as

$\mathbf{C}=\frac{\alpha
G^{+}\overline{\mathbf{C}}^{+}\,+\,(1-\alpha)G^{-}\overline{\mathbf{C}}^{-}}{G}\,.$
(18)

Note that the conformation tensor is symmetric and positive definite by
definition. We exploit the diagonalizability of the conformation tensor, along
with its three real and positive eigenvalues, as a key element in the
numerical stabilization approach for high Weissenberg numbers. For this
reason, the mathematical model is formulated with the conformation tensor
rather than the polymer stress.

To describe the behavior of the conformation tensor, the introduction of
rheological constitutive equations is required. In this study, the focus is
placed on rate-type partial differential equations37, 35, 36. Up to this
point, the equations presented are applicable to scenarios with two
viscoelastic phases. However, for the polymer stress, a closure is provided
specifically for the case in which only one phase exhibits viscoelasticity.
Under this condition, (17) reduces to

$\displaystyle\boldsymbol{\tau}_{p}=\alpha\,G^{+}\left(\overline{\mathbf{C}}^{+}\,-\mathbf{I}\right).$
(19)

Then, a generic constitutive equation for $\overline{\mathbf{C}}^{+}$ is given
by18

$\displaystyle\partial_{t}\overline{\mathbf{C}}^{+}+\mathbf{u}{\scriptscriptstyle\mathchoice{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\nabla\overline{\mathbf{C}}^{+}-\alpha\,(\nabla\mathbf{u})\\!^{\top}{\scriptscriptstyle\mathchoice{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\overline{\mathbf{C}}^{+}-\alpha\,\overline{\mathbf{C}}^{+}{\scriptscriptstyle\mathchoice{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\nabla\mathbf{u}=\frac{\alpha}{\lambda^{+}}\boldsymbol{\mathcal{P}}(\overline{\mathbf{C}}^{+}).$
(20)

The function $\boldsymbol{\mathcal{P}}(\overline{\mathbf{C}}^{+})$ in (20) is
a real analytic and isotropic tensor function of the conformation tensor. As a
consequence, we can express
$\boldsymbol{\mathcal{P}}(\overline{\mathbf{C}}^{+})$, using the well-known
tensor representation theorems 38, 39 as

$\boldsymbol{\mathcal{P}}(\overline{\mathbf{C}}^{+})=g_{0}\mathbf{I}+g_{1}\overline{\mathbf{C}}^{+}+g_{2}{\overline{\mathbf{C}}^{+}}^{2},$
(21)

where the coefficients $g_{0}$, $g_{1}$, and $g_{2}$ in (21) are isotropic
scalar functions of the three principal invariants of the conformation tensor.
Thus, the specific expressions for $g_{0}$, $g_{1}$, and $g_{2}$ in (21)
determine the viscoelastic behavior of the material. Table 1 lists some
suitable constitutive expressions for the coefficients used in this work.

Table 1: Scalar coefficients for the generic constitutive equation (20) and (21). constitutive model |  $g_{0}$ |  $g_{1}$ |  $g_{2}$ |
---|---|---|---|---
Maxwell/Oldroyd-B | $1$ | $-1$ | $0$ |
Giesekus | $1-\alpha$ | $2\alpha-1$ | $-\alpha$ |

Equations (19) to (21) provide a closure for the polymer stress term in the
VOF equations, specifically when phase $\Omega^{+}$ is known to be
viscoelastic, while phase $\Omega^{-}$ is considered purely viscous. The
challenge of finding an appropriate closure when both phases are viscoelastic
remains, to the best of the authors’ knowledge, an unresolved issue in the
field of CVA and is not within the scope of the present study.

##  3 Stabilization

To address the high Weissenberg number problem (HWNP) 8, 9, numerical
stabilization techniques are necessary to enhance the numerical stability of
the method at moderate and high levels of fluid elasticity. Fattal and
Kupferman 40 demonstrated that using a logarithmic transformation for the
conformation tensor equation can circumvent the HWNP. Balci et al. 41 proposed
a square root representation for the conformation tensor, which eliminates the
need for diagonalization. Afonso et al. 42 generalized these approaches by
introducing generic kernel transformation functions. The present work employs
a unified mathematical and numerical stabilization framework proposed by
Niethammer et al.10, 23 and extended to two-phase flows within an algebraic
VOF framework 18. This section provides a brief summary of the general aspects
of stabilization framework in the context of the geometric VOF method.

The generic constitutive equation for the conformation tensor (20) can be
transformed through a change of variable to obtain an evolution equation for a
real-analytic tensor function
${\boldsymbol{\mathcal{F}}(\overline{\mathbf{C}}^{+})}$. Assuming that
${\boldsymbol{\mathcal{F}}(\overline{\mathbf{C}}^{+})}$ is an isotropic
tensor-valued function of the conformation tensor, then the volume averaged
change-of-variable representation is given by

$\partial_{t}{\boldsymbol{\mathcal{F}}(\overline{\mathbf{C}}^{+})}+(\mathbf{u}{\scriptscriptstyle\mathchoice{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\nabla){\boldsymbol{\mathcal{F}}(\overline{\mathbf{C}}^{+})}=\alpha\left({{2\mathbf{B}{\scriptscriptstyle\mathchoice{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\overline{\mathbf{\Upsilon}}^{+}{\scriptscriptstyle\mathchoice{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\overline{\mathbf{C}}^{+}}+\mathbf{\Omega}{\scriptscriptstyle\mathchoice{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}{\boldsymbol{\mathcal{F}}(\overline{\mathbf{C}}^{+})}-{\boldsymbol{\mathcal{F}}(\overline{\mathbf{C}}^{+})}{\scriptscriptstyle\mathchoice{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\mathbf{\Omega}+\frac{1}{\lambda}\overline{\mathbf{\Upsilon}}^{+}{\scriptscriptstyle\mathchoice{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\boldsymbol{\mathcal{P}}(\overline{\mathbf{C}}^{+})}\right).$
(22)

A comprehensive derivation of the representation (22) can be found in 18. In
(22), the deformation terms present in the convective derivative are
decomposed into the first three terms on the right-hand side, employing the
tensors $\mathbf{B}$ and $\mathbf{\Omega}$. This local decomposition was first
proposed in 40 and has been applied to the volume-averaged VOF equations in
18. Such a decomposition is generally necessary for transforming the
constitutive equation for the conformation tensor (20) into an evolution
equation for the real-analytic tensor function
${\boldsymbol{\mathcal{F}}(\overline{\mathbf{C}}^{+})}$, by a change of
variable. The tensors $\mathbf{B}$ and $\mathbf{\Omega}$ in (22) are defined
as mixture quantities that are dependent on the gradient of the mixture
velocity. Their functional relationship with the velocity gradient is
presented in Appendix A.

The generic constitutive equation (22) reduces to the classical logarithm or
root conformation representations (LCR or RCR for short), depending on the
choice of ${\boldsymbol{\mathcal{F}}(\overline{\mathbf{C}}^{+})}$ and
$\overline{\mathbf{\Upsilon}}^{+}$. Table 2

Table 2: Function ${\boldsymbol{\mathcal{F}}(\overline{\mathbf{C}}^{+})}$ and $\mathbf{\Upsilon}$ terms in (22) for the conformation tensor representation, the $k$-th root of $\mathbf{C}$ representation, and the logarithm of $\mathbf{C}$ to base $a$ representation. | $\boldsymbol{\mathcal{F}}(\mathbf{C})$ | $\mathbf{\Upsilon}$
---|---|---
conformation tensor | $\mathbf{C}$ | $\mathbf{I}$
$k$-th root of $\mathbf{C}$ | $\mathbf{R}=\mathbf{Q}{\scriptscriptstyle\mathchoice{\mskip 3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip 3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower 0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower 0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\mathbf{D}^{\frac{1}{k}}{\scriptscriptstyle\mathchoice{\mskip 3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip 3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower 0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower 0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\mathbf{Q}\\!^{\top}$ | $\frac{1}{k}\mathbf{R}^{1-k}$
logarithm of $\mathbf{C}$ to base $a$ | $\mathbf{S}=\mathbf{Q}{\scriptscriptstyle\mathchoice{\mskip 3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip 3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower 0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower 0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}{\log}_{a}{(\mathbf{D})}{\scriptscriptstyle\mathchoice{\mskip 3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip 3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower 0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower 0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\mathbf{Q}\\!^{\top}$ | $\frac{1}{\operatorname{ln}(a)}a^{-\mathbf{S}}$

shows the specific terms corresponding to different conformation
representations. In the present work, we choose the RCR, which reduces (22) to
a more concise form according to:

$\displaystyle\frac{\mathrm{d}\mathbf{R}}{\mathrm{d}t}$
$\displaystyle=\frac{2}{k}\mathbf{B}{\scriptscriptstyle\mathchoice{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\mathbf{R}+\mathbf{\Omega}{\scriptscriptstyle\mathchoice{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}{\mathbf{R}}-{\mathbf{R}}{\scriptscriptstyle\mathchoice{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\mathbf{\Omega}+\frac{1}{k\lambda}\left(g_{0}\mathbf{R}^{1-k}+g_{1}\mathbf{R}+g_{2}\mathbf{R}^{1+k}\right).$
(23)

For the scope of the present study, the square root formulation with $k=2$
turned out to be stable across a range of Deborah numbers from $0$ to $16$.

##  4 Numerical method

We utilize the unstructured geometrical Volume-of-Fluid method 32, because its
discretization of the fluid interface consistently matches with the
discretization of eq. 19.

###  4.1 The geometrical Volume-of-Fluid method

The flux-based geometrical VOF method updates volume fractions defined in eq.
3 by integrating the flux of the phase $\Omega^{+}(t)$ through the boundary
$\partial\Omega_{c}$ of the fixed control volume $\Omega_{c}$ over the time
step $[t^{n},t^{n+1}]$, namely

$\alpha_{c}(t^{n+1})=\alpha_{c}(t^{n})-\frac{1}{|\Omega_{c}|}\int_{t^{n}}^{t^{n+1}}\int_{\partial\Omega_{c}}\chi(\textbf{x},t)\mathbf{v}\cdot\mathbf{n}\leavevmode\nobreak\
dS\leavevmode\nobreak\ dt.$ (24)

The reader is directed to 32 for a derivation of eq. 24. Let us also note that
we have chosen to use $\Omega_{c}$ for the shifted averaging volume to have a
closer connection to the Finite Volume notation. Moreover, as we only use a
single centered value $\alpha_{c}$ for each cell as opposed to a full field of
the averaged quantity, there is no spatial variable left in $\alpha_{c}$ in
(24) and below. The boundary $\partial\Omega_{c}$ of the cell $\Omega_{c}$ in
eq. 24 is a union of surfaces (faces) bounded by line segments (edges), namely

$\partial\Omega_{c}:=\bigcup_{f\in F_{c}}S_{f}.$ (25)

Using eq. 25, eq. 24 is reformulated according to

$\alpha^{n+1}_{c}=\alpha^{n}_{c}-\frac{1}{|\Omega_{c}|}\sum_{f\in
F_{c}}\int_{t^{n}}^{t^{n+1}}\int_{S_{f}}\chi(\textbf{x},t)\mathbf{v}\cdot\mathbf{n}\leavevmode\nobreak\
dS\leavevmode\nobreak\ dt.$ (26)

Equation 26 is an exact identity provided that $\partial\Omega=\cup_{c\in
C}(\partial\Omega\bigcap\partial\Omega_{c})$, and because eq. 26 does not
involve any approximations. It is important to note that eq. 26 utilizes the
sharp phase indicator $\chi(\mathbf{x},t)$ to obtain
$\alpha^{n+1}_{c}:=\alpha_{c}(t^{n+1})$. Using $\chi(\mathbf{x},t)$ is
advantageous for discretizing the equations for two-phase viscoelastic flows
because it is also used for the conditional volume averaging in the derivation
of the single-field formulation of two-phase Navier-Stokes equations for
viscoelastic flows 18.

The second-order accurate Unstructured Finite Volume Method (UFVM)43 utilizes
second-order accurate centroid averages

$\phi_{D}=\frac{1}{|D|}\int_{D}\phi(\mathbf{x})dx=\phi(\mathbf{x}_{c})+\frac{1}{|D|}\int_{D}\nabla\nabla\phi:(\mathbf{x}-\mathbf{x}_{c})\otimes(\mathbf{x}-\mathbf{x}_{c})\,d\textbf{x}$
(27)

where $D$ can be the control volume $\Omega_{c}$ or one of its faces $S_{f}$,
in which case the integral is understood as a surface integral, i.e. $dx$ is
the surface measure $dS$, and $x_{c}$ is their centroid. Therefore, velocity
and the face-area-normal vector in eq. 26 are averaged, resulting in

$\displaystyle\alpha^{n+1}_{c}$
$\displaystyle=\alpha^{n}_{c}-\frac{1}{|\Omega_{c}|}\sum_{f\in
F_{c}}\int_{t^{n}}^{t^{n+1}}\mathbf{v}_{f}\cdot\frac{\mathbf{S}_{f}}{|\mathbf{S}_{f}|}\int_{S_{f}}\chi(\textbf{x},t)\leavevmode\nobreak\
dS\leavevmode\nobreak\ dt+e_{\alpha_{h}}(h^{2})$ (28)
$\displaystyle=\alpha_{c}(t^{n})-\frac{1}{|\Omega_{c}|}\sum_{f\in
F_{c}}\int_{t^{n}}^{t^{n+1}}F_{f}(t)|A_{f}(t)|\leavevmode\nobreak\
dt+e_{\alpha_{h}}(h^{2})$

with $|A_{f}|:=\int_{S_{f}}\chi(\textbf{x},t)\leavevmode\nobreak\ dS$ defining
the phase-specific face area, and
$F_{f}:=\int_{S_{f}}\mathbf{v}_{f}\cdot\mathbf{n}dS$ the volumetric flux. The
integral on the r.h.s. of eq. 28 represents the volume of the phase
$\Omega^{+}$ (_phase-specific volume_) fluxed through the face $S_{f}$ over
the time interval $[t^{n},t^{n+1}]$.

Each flux-based unstructured geometrical VOF method computes the fluxed phase-
specific volume differently32, we use specifically the plicRDF-isoAdvector
method 24 because of its computational efficiency.

###  4.2 The plicRDF-isoAdvector Volume-of-Fluid method

The plicRDF-isoAdvector method 24 approximates the phase indicator function
with $\tilde{\chi}(\mathbf{x},t)\approx\chi(\mathbf{x},t)$ as

$\tilde{\chi}_{c}(\mathbf{x},t)=\begin{cases}1,\quad P_{c}(\mathbf{x},t)\leq
0,\\\ 0,\quad P_{c}(\mathbf{x},t)>0,\end{cases}$ (29)

with

$P_{c}(\mathbf{x},t)=a_{c,0}(t)+a_{c,1}(t)x+a_{c,2}(t)y+a_{c,3}(t)z,$ (30)

defined in such a way that, given $\mathbf{n}_{c}=[a_{c,1},a_{c,2},a_{c,3}]$,
the value $a_{c,0}$ at some time $t$ is found that satisfies

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

the so-called interface positioning problem that ensures volume conservation
44, 45 of the phase-indicator function approximation. The plicRDF-isoAdvector
24 relies on the reconstructed distance function to iteratively improve the
initial estimate of $\mathbf{n}_{c}$. Cell-specific $P_{c}(\mathbf{x},t)$
makes the approximation of the fluid interface piecewise linear.

The plicRDF-isoAdvector further approximates eq. 28 into

$\displaystyle\alpha^{n+1}_{c}$
$\displaystyle\approx\alpha_{c}(t^{n})-\frac{1}{|\Omega_{c}|}\sum_{f\in
F_{c}}0.5(F_{f}^{n}+F_{f}^{n+1})\int_{t^{n}}^{t^{n+1}}|A_{f}(t)|\leavevmode\nobreak\
dt,$ (32)

computing $\int_{t^{n}}^{t^{n+1}}|A_{f}(t)|\leavevmode\nobreak\ dt$
geometrically, by approximating $|A_{f}|^{n+1}$ by translating the plane
$P_{c}(\mathbf{x},t^{n})$ along $n_{c}^{n}$ with the displacement given by
$\mathbf{v}_{c}\Delta t$, with the time step $\Delta t:=t^{n+1}-t^{n}$, as
schematically shown in fig. 2(a). This dimensional reduction in computing the
fluxed phase-specific volume via
$\int_{t^{n}}^{t^{n+1}}|A_{f}(t)|\leavevmode\nobreak\ dt$ instead of
reconstructing the 3D fluxed phase-specific volume46, 32 ensures a high-degree
of computational efficiency for the plicRDF-isoAdvector method.

((a)) Upwind $A_{f}^{n+1}$ computed in eq. 32 during $\alpha_{c}^{n+1}$ update
with $\mathbf{n}_{c}^{n}$.

((b)) $A_{f}^{n+1}$ reconstructed using $\mathbf{n}_{c}^{n+1}$ and
$\alpha_{c}^{n+1}$.

Figure 2: Different ways to compute $A_{f}^{n+1}$ for $\alpha_{f}$ weighting
of the polymer stress $\tau_{p}$ in the pressure Poisson equation 33.

###  4.3 Solution algorithm

The key to simulating a viscoelastic fluid phase with the geometrical VOF
method is eq. 19 and its use in the Poisson equation for the pressure in the
segregated solution algorithm of eqs. 8, 7 and 28 and eq. 22 with a
representation of the conformation tensor selected from table 2. Without
diving into a detailed derivation of segregated solution algorithms 47, 34, we
focus on the discrete pressure Poisson equation, specifically, on the
divergence of the polymer stress component $\nabla\cdot{\tau}_{p}$ from eq. 19
that appears on the right-hand side of the pressure Poisson equation,

$\sum_{f\in F_{c}}\left(\frac{1}{a_{c}}\right)_{f}(\nabla
p)_{f}\cdot\mathbf{S}_{f}=\sum_{f\in
F_{c}}\left(\frac{1}{a_{c}}\right)_{f}\alpha_{f}(\tau_{p})_{f}\cdot\mathbf{S}_{f}+\dots$
(33)

The consistency between $\alpha$ and $\alpha_{c}$ when the finite volume
$\Omega_{c}$ is used for conditional volume averaging makes $\alpha_{f}$ in
eq. 33 defined as

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

numerically consistent with the geometrical VOF method. Although the plicRDF-
isoAdvector method relies only on $|A_{f}|$ to update $\alpha_{c}^{n+1}$,
other geometrical VOF methods that compute $|A_{f}|$ as a part of the
approximation of the fluxed phase-specific volume in eq. 28 are equivalently
numerically consistent with eq. 19.

The segregated solution algorithm is outlined in algorithm 1 and it ensures
eqs. 8, 7, 28 and 22 are satisfied at $t^{n+1}$, therefore, we require
$\alpha_{f}^{n+1}=\frac{|A_{f}^{n+1}|}{|S_{f}|}$ in eq. 33. While uniquely
defined for each $S_{f}$, the obvious choice of the upwind $|A_{f}^{n+1}|$
from the plicRDF-isoAdvector eq. 32 is not the right choice (cf. fig. 2(a)).
Namely, when the plicRDF-isoAdvector method computes $|A_{f}^{n+1}|$ in eq.
32, it does so by translating the interface plane given by eq. 30 in the
direction of the interface-normal at the _old time step_ , namely
$\mathbf{n}_{c}^{n}$, thus computing $A_{f}^{n+1}$ as a first-order accurate
estimate in time, which is acceptable because in eq. 32 the integration over
$\Delta t$ recovers second-order accuracy in time. However, upwind
$A_{f}^{n+1}$ results from the old interface normal $\mathbf{n}_{c}^{n}$; only
after $\alpha_{c}^{n+1}$ has been updated using eq. 32, and the interface has
been again reconstructed at $t^{n+1}$ by optimizing $\mathbf{n}_{c}^{n+1}$ and
positioning the interface by finding $a_{c,0}$, do the area fractions
$\alpha_{f}^{n+1}$ correspond to the phase indicator approximation
$\tilde{\chi}_{c}(\mathbf{x},t)$ at $t^{n+1}$. Since the solution algorithm
must satisfy equation eq. 33 at $t^{n+1}$, $|A_{f}^{n+1}|$ should be defined
by $\chi_{c,k}(\mathbf{x},t^{n+1})$, where $c,k$ denote indices of cells that
share $S_{f}$. We therefore advect $\alpha_{c}^{n+1}$ using the plicRDF-
isoAdvector scheme to solve eq. 32, then reconstruct the interface to obtain
$\tilde{\chi}_{c}(\mathbf{x},t^{n+1})$, as shown schematically in fig. 2(b).
The issue in computing the reconstructed $|A_{f}^{n+1}|$ is the piecewise
planarity of $\tilde{\chi}$, resulting in two candidate values from cells
$\Omega_{k,c}$, namely

$|A_{k,c}(t)|=\int_{S_{f}}\tilde{\chi}_{k,c}(\mathbf{x},t)\,dS,$ (35)

and we solve the issue of uniqueness of $A_{f}^{n+1}$ by averaging

$|A_{f}^{n+1}|=0.5(|A_{k}|^{n+1}+|A_{c}|^{n+1}),$ (36)

ensuring sufficient mesh resolution for the plicRDF-isoAdvector scheme to
reduce the discontinuity of $\tilde{\chi}$ on $\partial\Omega_{c}$.

Algorithm 1 Solution Algorithm

1:for $o=1$ to $N_{\text{outer\\_iterations}}$ do

2: Solve the conformation tensor equation (22) for $(\tau_{p})^{o}_{c}$.

3: Assemble (solve) the momentum equation (8) for $\mathbf{v}_{c}^{o}$.

4: Solve the volume fraction equation (32).

5: Reconstruct the interface to obtain $\tilde{\chi}_{c}^{o}$.

6: Compute reconstructed area fractions $\alpha_{f}^{o}$.

7: Assemble the viscoelastic stress term (19) for the pressure equation (33).

8: for $i=1$ to $N_{\text{inner\\_iterations}}$ do

9: Solve the pressure equation (33).

10: Update the cell-centered velocity $\mathbf{v}_{c}^{i}$ and flux
$F_{f}^{i}$.

11: end for

12:end for

##  5 Results

###  5.1 2D drop deformation in simple shear flow

A two-dimensional single drop deforming under simple shear flow is studied as
an initial test of the method. This flow was studied numerically by several
researchers, see for example 48, 12, 49, for different combinations of the
droplet’s material properties and the surrounding fluid matrix. Here, we are
interested in the configurations drop/matrix: Newtonian/Newtonian (NN),
Newtonian/viscoelastic (NV) and viscoelastic/Newtonian (VN).

####  5.1.1 Problem Description

A two-dimensional circular drop (fluid 1) is immersed in the center of a
rectangular domain, as shown in Figure 3.

Figure 3: Schematic diagram of the two-dimensional drop deformation in shear
flow.

The channel height equals $H=8r$, where $r$ is the initial drop radius. The
horizontal channel length is $L=20r$. Periodic boundary conditions are imposed
at the two outer channel boundaries. At the channel walls, Diriclet boundary
conditions are prescribed with a constant value for the velocity. The (t)op
and (b)ottom walls move with a speed $U_{t}=-U_{b}$ to generate a shear flow.
For the stress and the generic constitutive variable, a zero gradient boundary
condition is defined for the channel walls. The simulation is initialized by a
fully developed shear flow velocity field, while the other fields are
initially zero.

The total viscosities of the fluid $i=1,2$ are defined as
$\eta_{i}={\eta}_{s}{{}_{i}}+{\eta}_{p}{{}_{i}}$, where ${\eta}_{s}{{}_{i}}$
is the solvent dynamic viscosity and ${\eta}_{p}{{}_{i}}$ the polymer
viscosity. The relaxation time is denoted as $\lambda_{i}$, the density
$\rho_{i}$ and the surface tension is $\sigma$. The constant shear rate reads
$\dot{\gamma}=(U_{t}-U_{b})/H$. With these definitions, the system is governed
by the following dimensionless numbers: Reynolds number $\operatorname{Re}$,
Deborah number $\operatorname{De}$, capillary number $\operatorname{Ca}$ and
viscosity ratio $\Pi$ of drop (fluid 1) to matrix (fluid 2) defined by

$\displaystyle\operatorname{Re}=\frac{\rho_{2}\dot{\gamma}r^{2}}{\eta_{2}},\quad\operatorname{De}=\lambda\dot{\gamma},\quad\operatorname{Ca}=\frac{\eta_{2}\dot{\gamma}r}{\sigma},\quad\Pi=\frac{\eta_{1}}{\eta_{2}}.$
(37)

For the present study, the parameters read $\operatorname{Re}=0.0003$, $\Pi=1$
and for the respective viscoelastic phase
${\eta}_{s}{{}_{i}}={\eta}_{p}{{}_{i}}$. The Deborah number and the Capillary
number are varied. The drop deformation is studied by means of the deformation
parameter3

$D=\frac{r_{\textnormal{max}}-r_{\textnormal{min}}}{r_{\textnormal{max}}+r_{\textnormal{min}}},$
(38)

where $r_{\textnormal{max}}$ and $r_{\textnormal{min}}$ are the largest and
shortest distances from the drop center to its boundary, respectively. The
orientation angle $\theta$ is defined as the angle between the principal axes
of the deformed drop and the flow direction.

####  5.1.2 Results and discussion

A mesh convergence study is carried out to verify the algorithm and determine
the required mesh cells per initial drop radius to achieve mesh-independent
results. Figure 4 displays the deformation for various mesh resolutions over
time. The coarsest mesh configuration has $20$ cells per drop radius, totaling
$64000$ cells. In contrast, the finest mesh consists of $50$ cells per drop
radius, summing up to $400000$ cells. The VN system appears to have a more
evident mesh sensitivity compared to the NV system. However, convergence
towards the finest resolution is achieved for both systems. A steady-state
solution is reached by the time $t\dot{\gamma}=8$. Using the densest mesh as
the reference solution, the deformation error on the coarsest configuration,
with $20$ cells per drop radius, remains under $1\%$ in the steady-state. This
can be seen from Table 3, which presents the deformation parameter across
different mesh resolutions. During the transient deformation phase, the
maximum error exceeds that of the steady-state. For subsequent analyses, a
mesh resolution of $50$ cells per drop radius is selected to ensure highest
accuracy during both transient and steady-state deformations.

Figure 4: Drop deformation as a function of time with varying mesh refinement for the VN and NV systems at $\operatorname{Ca}=0.24$ and $\operatorname{De}=0.4$.  Table 3: Comparison of drop deformation $D$ and orientation angle $\theta$ at time $t=8\leavevmode\nobreak\ \dot{\gamma}^{-1}$. NV and VN systems for $\operatorname{Ca}=0.24$ and $\operatorname{De}=0.4$ on different mesh resolutions. | cells per radius $=20$ |  | cells per radius $=35$ |  | cells per radius $=50$
---|---|---|---|---|---
| $D$ | $\theta$ |  | $D$ | $\theta$ |  | $D$ | $\theta$
NV | 0.2418 | 25.36 |  | 0.2408 | 25.14 |  | 0.241 | 25.97
VN | 0.2444 | 33.3 |  | 0.2473 | 31.21 |  | 0.2448 | 32.10

A further mesh convergence study is performed using two different curvature
models: RDF and height function, both of which are available in the
TwoPhaseFlow OpenFOAM project 27. More details on these models can be found in
26. Figure 5 shows a comparison of the drop deformation as a function of time,
using these curvature models across various mesh resolutions. The findings
suggest that the RDF curvature model offers superior convergence and accuracy
compared to the height function. As a result, the RDF curvature model is
selected for use in this research. Recall that the height function approach is
highly accurate for Cartesian meshes, its application to a general
unstructured mesh, as used here, significantly complicates the implementation.
In its current state the RDF curvature model is slightly more accurate, but we
intend to enhance the accuracy of the height function implementation in future
developments.

Figure 5: Drop deformation as a function of time for the VN system at
$\operatorname{Ca}=0.6$ and $\operatorname{De}=0.4$ with the height function
(hf) and RDF curvature models.

The results from the finest mesh with 50 cells per drop radius are validated
against simulation data from the literature in Table 4. This comparison covers
the Newtonian drop in a Newtonian matrix as well as the viscoelastic NV and VN
scenarios for $\operatorname{Ca}=0.24$ and $\operatorname{De}=0.4$. The data
refer to the time step at $t=3\,\dot{\gamma}$, at which a steady-state is not
yet reached completely in all configurations, which can be seen from the
comparison with Table 3. Notably, there are more substantial discrepancies
among all sources at $t=3\,\dot{\gamma}$ for the NV and VN systems than for
the Newtonian case. The results of the present study indicate a reduced drop
deformation in the viscoelastic cases compared to the Newtonian case,
consistent with existing literature. The smallest deformation parameter $D$ is
observed in the NV combination, confirming previous research findings.
Moreover, the values of the deformation parameter closely correspond with
those reported in 49 across all the systems: NN, NV, and VN, albeit being
slightly smaller in our case. The orientation angle $\theta$ lies within the
range reported in references 12 and 49 for both Newtonian and viscoelastic
cases. Consistent with prior research, the orientation angle in the VN system
is larger than that in the NN system, while the NV system exhibits the
smallest value for $\theta$.

Table 4: Comparison of drop deformation $D$ and orientation angle $\theta$ at time $t=3\leavevmode\nobreak\ \dot{\gamma}^{-1}$. NN, NV and VN systems for $\operatorname{Ca}=0.24$ and $\operatorname{De}=0.4$. | This work |  | 12 |  | 49
---|---|---|---|---|---
| $D$ | $\theta$ |  | $D$ | $\theta$ |  | $D$ | $\theta$
NN | 0.2559 | 31.63 |  | 0.2878 | 32.26 |  | 0.2674 | 28.90
NV | 0.2441 | 26.88 |  | 0.2656 | 28.15 |  | 0.2550 | 25.38
VN | 0.247 | 32.03 |  | 0.2799 | 32.53 |  | 0.2582 | 30.32

Figure 6 displays the transient deformation for varying $\operatorname{Ca}$
and $\operatorname{De}$ values.

((a)) VN,$\;\operatorname{Ca}=0.24,\,0.6$.

((b)) NV, $\operatorname{Ca}=0.24,\,0.6$.

((c)) VN,$\;\operatorname{Ca}=0.24,\,0.6$.

((d)) NV,$\;\operatorname{Ca}=0.24,\,0.6$.

Figure 6: Drop deformation parameter $D$ and orientation angle $\theta$ as a
function of time for the NN, VN, and NV systems at $\operatorname{Ca}=0.24$
and $\operatorname{Ca}=0.6$ and varying Deborah numbers.

At lower Capillary and Deborah numbers, a steady-state is achieved more
rapidly in terms of physical time. With increasing $\operatorname{De}$ and
$\operatorname{Ca}$, the drop exhibits prolonged dynamic changes. Moreover, at
$\operatorname{Ca}=0.6$, there is a significant increase in deformation across
all scenarios compared to $\operatorname{Ca}=0.24$. Comparing Figures 6(a) and
6(b) indicates that an increase in $\operatorname{De}$ has a more pronounced
effect on the NV system, where the matrix fluid exhibits viscoelastic
properties. Specifically, within the range of $\operatorname{De}=4$ to $16$,
the NV system shows a significantly more pronounced deformation behavior, with
the steady-state solution not being reached even after a simulation time of
$t\dot{\gamma}=80$. This can be attributed to the matrix fluid’s inherent
ability to accumulate and release elastic stress over time, observable even at
considerable distances from the deforming fluid interface.

Figures 6(c) and 6(d) display the transient orientation angle $\theta$. As the
Capillary numbers increase, there is a notable decrease in the size of the
orientation angle. In the VN system, an increase of the Deborah number up to
$\operatorname{De}=2$ results in increased $\theta$ values. At Deborah numbers
above this threshold, the orientation angle shows a slight reduction, yet the
values remain substantially higher than those observed in the Newtonian drop.
Conversely, in the NV system, the orientation angle decreases as the Deborah
numbers increase. The effect of varying Deborah numbers is more pronounced,
compared to the VN system. Specifically, at $\operatorname{Ca}=0.6$ and
$\operatorname{De}=16$, the orientation angle approaches almost zero.

Further analysis reveals a distinct pattern in the steady-state deformation,
as shown in Figure 7: there is a decline within the range of
$\operatorname{De}=0$ to $\operatorname{De}=1$, followed by an increase for
higher Deborah numbers, specifically when $\operatorname{De}>1$. This behavior
aligns with previous observations made by Yue 50 and Aggarwal and Sarkar 51
for a Deborah number range of up to $\operatorname{De}\leq 2$. The present
study confirms these earlier findings and demonstrates a sustained increase of
the deformation parameter in the range $1<\operatorname{De}\leq 16$. The
impact of the Deborah number on steady-state deformation is considerably more
pronounced in the NV system. As shown in Figure 7, the deformation surpasses
the Newtonian reference state when $\operatorname{De}>4$. Hence, depending on
the particular value of the Deborah number, viscoelasticity can either
suppress or enhance droplet deformation relative to a Newtonian drop in a
Newtonian matrix. We emphasize that no steady-state is reached within
$t=80\leavevmode\nobreak\ \dot{\gamma}^{-1}$ for $\operatorname{De}=16$; the
plot shows the final simulation values in that case. For future research, it
would be valuable to explore whether a critical Deborah number exists beyond
which a drop can no longer maintain a stationary shape. Previously, a critical
Capillary number $\operatorname{Ca}_{c}$ was identified at around $0.875$ 52
beyond which periodic drop oscillations occur. However, higher critical
Capillary numbers around $1.0$ have been reported 53. Our simulations suggest
that both the Deborah number and the Capillary number play significant roles.
It appears that at lower Capillary numbers, oscillatory behavior may emerge
when the Deborah number is sufficiently high.

((a)) VN,$\;\operatorname{Ca}=0.24,\,0.6$.

((b)) NV,$\;\operatorname{Ca}=0.24,\,0.6$.

Figure 7: Drop deformation as a function of Deborah numbers for the NN, VN,
and NV systems at $\operatorname{Ca}=0.24$ and $\operatorname{Ca}=0.6$. The
deformation parameter D value corresponds to the latest D values of Figure 6.

The drop shapes are visualized in Figure 8 at $t=10\,\dot{\gamma}^{-1}$ and
for the latest times up to $t=80\,\dot{\gamma}^{-1}$ in Figure 9. Figure 8
serves as a reference for comparing with literature results where simulations
were halted prior to reaching a steady-state 12. For the lower Capillary
number, the drops have a more ellipsoidal shape. At $\operatorname{Ca}=0.6$,
the stretching is significantly higher, which corresponds to the results of
the deformation parameter reported above. For the Newtonian matrix fluid, the
drop continues to deform and change its orientation, as evident when comparing
Figures 8(b) and (d) to Figures 9(b) and (d). In particular, with an increase
in the Deborah number, the drops tend to align more closely with the direction
of the flow, resulting in lower orientation angles. At $\operatorname{Ca}=0.6$
and $\operatorname{De}=16$, the drop shape differs significantly from an
ellipsoidal shape. A wave-like structure is formed on the upper and lower
sides of the interface. We assume that this is related to local extrema of the
first normal stress difference observed near the interface, which have an
impact on the flow field.

((a)) VN,$\;\operatorname{Ca}=0.24$.

((b)) NV,$\;\operatorname{Ca}=0.24$.

((c)) VN,$\;\operatorname{Ca}=0.6$.

((d)) NV,$\;\operatorname{Ca}=0.6$.

Figure 8: Drop shapes at $t=10\,\dot{\gamma}^{-1}$.

((a)) VN,$\;\operatorname{Ca}=0.24$.

((b)) NV,$\;\operatorname{Ca}=0.24$.

((c)) VN,$\;\operatorname{Ca}=0.6$.

((d)) NV,$\;\operatorname{Ca}=0.6$.

Figure 9: Drop shapes at the latest time corresponding to Figure 6.

Next, we examine the effect of different rheological models. In this context,
simulations are carried out using the Giesekus model, which introduces an
additional nonlinear term compared to the Oldroyd-B model. The mobility
parameter $\alpha$ in the Giesekus model, an empirical constant, determines
the degree of nonlinear features to be considered. When $\alpha$ is set to
zero, the Oldroyd-B model is recovered. A slightly increased $\alpha$ value of
$0.03$ already has a significant impact on the flow behavior, which can be
seen in Figure 10. While the deformation remains nearly unchanged in the VN
system, a significant deviation occurs in the NV system. Figure 10(b) shows
that the deformation decreases with a mobility parameter $\alpha=0.03$. This
effect becomes more pronounced as the Deborah number increases. In contrast to
the Oldroyd-B model, the Giesekus model shows a progressive reduction in the
deformation parameter as the Deborah number increases. Additionally, the size
and duration of the droplet dynamics are significantly smaller. The Giesekus
matrix fluid at $\operatorname{De}=16$ approaches a steady-state deformation
within the interval $t\dot{\gamma}=40$. For the Oldroyd-B model this time
interval is insufficient and it remains unclear whether a steady-state
solution exists.

((a)) VN.

((b)) NV.

Figure 10: Comparison of rheological models: Oldroyd-B (dotted) and the
Giesekus model with $\alpha=0.03$ (dashed) for droplet deformation as a
function of time for NN, VN, and NV systems at $\operatorname{Ca}=0.24$ and
varying Deborah numbers.

Figure 11 shows the corresponding drop orientation angles $\theta$ for the
Giesekus and the Oldroyd-B models. In the VN system, the orientation angles
are observed to be similar for both models. However, as shown in Figure 11(b),
the orientation angles in the Giesekus matrix fluid are consistently higher
compared to those in the Oldroyd-B fluid. This implies that in the Giesekus
fluid, the drop’s alignment with the flow direction is less pronounced.
Notably, the minimum $\theta$ is observed at a Deborah number of $4$. Beyond
this, with increasing $\operatorname{De}$, the orientation angle tends to rise
in the Giesekus fluid.

((a)) VN.

((b)) NV.

Figure 11: Comparison of rheological models: Oldroyd-B (dotted) and the
Giesekus model (dashed) for droplet deformation angle as a function of time
for NN, VN, and NV systems at $\operatorname{Ca}=0.24$ and varying Deborah
numbers.

Figure 12 illustrates the effect of further increasing the mobility parameter.
Raising $\alpha$ from $0.03$ to $0.1$ results in marginally increased
deformations in the VN system and substantially reduced deformations in the NV
system. In addition, with the increase of the Giesekus parameter, the drop
dynamics are further reduced in the NV system.

((a)) VN.

((b)) NV.

Figure 12: Drop deformation as a function of time for NN, VN and NV at
$\operatorname{Ca}=0.24$. Comparison of the Giesekus fluid for the model
parameters $\alpha=0.03$ (dotted) and $\alpha=0.1$ (dashed) at varying Deborah
numbers.

The similarities between the Oldroyd-B and Giesekus models in the VN system
are confirmed by the two-dimensional plot of the first normal stress
difference. Figure 13 shows the contours of the first normal stress difference
for $\operatorname{Ca}=0.24$ and $\operatorname{De}=0.4$ in the VN system. The
first normal stress difference is computed from the polymer stress components
as $N_{1}=\tau_{p,xx}-\tau_{p,yy}$. In the VN case, the distribution is nearly
identical. The first normal stress difference, $N_{1}$, has its maximum at the
poles and its lowest values near the equator. This $N_{1}$ distribution
indicates that with increasing $\operatorname{De}$ polymeric tensile stresses
at the poles impact the stretching of the drop, as described in 50. Figure 14
shows the corresponding contour maps for the NV cases. In the viscoelastic
matrix fluid, the first normal stress difference in the Giesekus fluid is
slightly smaller than in the Oldroyd-B fluid. The local extrema are positioned
similarly to those in the Oldroyd-B fluid; however, their spatial extent is
more confined. A maximum in $N_{1}$ is observed near the poles of the drop in
both cases.

((a)) VN, Oldroyd-B model.

((b)) VN, Giesekus model.

Figure 13: First normal stress difference, N1 contours for VN system for
$\operatorname{Ca}=0.24$, $\operatorname{De}=0.4$ for Oldroyd-B and Giesekus
model.

((a)) NV, Oldroyd-B model.

((b)) NV, Giesekus model.

Figure 14: First normal stress difference, N1 contours for NV system for
$\operatorname{Ca}=0.24$, $\operatorname{De}=0.4$ for Oldroyd-B and Giesekus
model.

###  5.2 3D drop deformation in simple shear flow

Three-dimensional numerical studies analyzing the deformation of viscoelastic
drops in simple shear have been extensively researched in the literature.
Khismatullin et al. 13 employed the Volume-of-Fluid method with a parabolic
interface reconstruction (VOF-PROST), utilizing Oldroyd-B and Giesekus fluids
for the viscoelastic phase. Aggarwal and Sarkar54, 51 applied a front-tracking
Finite-Difference method, leveraging the Oldroyd-B constitutive equation for
the viscoelastic fluid. Verhulst et al. 15 adopted the VOF-CSF method and used
both the Oldroyd-B and Giesekus models to characterize viscoelasticity.
Afkhami et al. 14 applied a Finite-Difference VOF method paired with a
piecewise linear interface reconstruction, using the Oldroyd-B constitutive
model for the viscoelastic phase. Mukherjee and Sarkar 55 utilized the front-
tracking method, as described by Aggarwal and Sarkar 54, to conduct in-depth
examinations of drop deformation with an Oldroyd-B viscoelastic phase.
Verhulst et al. 15, 56 employed the Volume-of-Fluid method, as described by
Khismatullin et al. 13, for comprehensive studies of the steady-state and
transient drop deformation, using the Oldroyd-B, Giesekus, Ellis, and multi-
mode Giesekus constitutive equations. The aforementioned studies have been
referenced in more recent publications for the purposes of validation and
comparison 57.

In this section, we demonstrate the capability of our geometric VOF method to
simulate three-dimensional drop deformations under shear flows. We explore
parameter variations of both the Capillary and Deborah numbers. Specifically,
we examine 3D drop deformations employing the Oldroyd-B constitutive equations
for Deborah numbers extending up to $\operatorname{De}=16$. This extends
beyond the scope of most of the previous studies referenced and also
underlines the efficacy of the high Deborah number stabilization implemented
in our study.

####  5.2.1 Problem description

The setup for the three-dimensional case is based on the aforementioned two-
dimensional study. The channel domain has the height $H=8r$, the length
$L=10r$ and the depth $Z=4r$. Initially, a spherical drop of radius $r$ is
placed at the center of this domain. The (t)op and (b)ottom walls move with a
constant velocity of $U_{t}$ = $-U_{b}$. Zero gradient boundary conditions for
the pressure, stress, and generic constitutive variable are defined at the
channel walls. Periodic boundary conditions are applied to the left and right
sides of the channel. The simulation is initialized by a fully developed shear
flow velocity field, while the other fields are initially zero.

The 3D simulation is governed by the dimensionless numbers defined in (37).
For all simulations, a Reynolds number of $Re=0.0003$ and a dimensionless
viscosity ratio of $\Pi=1$ are used. In the viscoelastic phase, the solvent
viscosity is set equal to the polymer viscosity. The analysis considers the
drop deformation parameter, given by (38), along with the orientation angle
$\theta$, representing the alignment of the principal axes of the deformed
drop with the flow direction.

####  5.2.2 Results and discussion

In Figure 15, a mesh convergence study is presented in 3D, examining drop
resolutions ranging from $12$ cells per radius to $32$ cells per radius. The
curve corresponding to a resolution of $12$ cells per radius is evidently
underresolved. However, the curves for the higher resolutions closely align,
indicating convergence towards a mesh-independent solution. The difference
between the resolutions of the two finest meshes is marginal. Based on this
observation, we consider the simulations with $32$ cells per radius to be
sufficiently resolved, and this mesh resolution is adopted for the 3D drop
deformation study. Note that the chosen 3D resolution of 32 cells per radius
is smaller than the resolution used in the 2D simulations. This decision
represents a compromise between accuracy and computational resources for the
parameter variation.

Figure 15: Drop deformation in 3D as a function of time with varying mesh
refinement for the VN system at $\operatorname{Ca}=0.24$ and
$\operatorname{De}=0.4$.

Figure 16 displays the transient deformation of a drop in 3D under varying
Capillary numbers: $\operatorname{Ca}=0.24$, $\operatorname{Ca}=0.4$, and
$\operatorname{Ca}=0.6$. For each Capillary number, the Deborah number is
varied from $\operatorname{De}=0$ to $\operatorname{De}=16$. The Oldroyd-B
fluid is used in the viscoelastic phase in all simulations. The VN system
results, shown in Figure 16(a), demonstrate that an increase in the Capillary
number leads to stronger drop deformation. Significant drop dynamics are
observed at $\operatorname{Ca}=0.6$ when $\operatorname{De}>1$, indicating
oscillative deformation patterns in $D$. It appears that the oscillatory
behavior may be damped up to at least $\operatorname{De}=8$; however, the
simulation time frame of $t\dot{\gamma}=50$ is insufficient to conclusively
determine the damping behavior. At $\operatorname{Ca}=0.4$, damped
oscillations in $D$ are observed on a much smaller scale, and within the given
simulation period, a steady-state deformation parameter is achieved for all
Deborah numbers. When compared to the 2D results shown in Figure 6, the 3D
simulations exhibit larger drop deformations for all Capillary numbers.
However, the overall qualitative behavior of the deformation remains very
similar. Beginning from the Newtonian reference state, an increase in the
Deborah number initially results in smaller deformation parameters, reaching a
minimum around $\operatorname{De}=2$. As the Deborah number continues to rise,
the deformation increases accordingly in the range up to
$\operatorname{De}=16$; the $D$ values remain smaller than the Newtonian
reference state within this range. Figure 16(b) presents the results for the
NV system, where again deformations in 3D are observed to be larger than those
in their 2D counterparts. Notably, in this system, the Deborah number has a
more substantial impact on drop deformation compared to the VN system.
Relative to the Newtonian reference state, the drop deformation parameter
decreases when $\operatorname{De}$ is increased to $0.4$, marking the point of
minimum deformation in 3D, which differs from the 2D case. As the Deborah
number is further increased, the deformation parameter begins to rise. For
$\operatorname{De}>1$, the deformation exceeds the Newtonian reference state.
The drop dynamics in the NV system appear to be less periodic compared to the
VN system. We do not observe any damped oscillations, as seen in the VN
system. The influence of the Deborah number on the drop dynamics is
substantial, even at the lowest Capillary number, $\operatorname{Ca}=0.24$. To
thoroughly study the drop dynamics at higher Deborah numbers, simulation times
greater than $t\dot{\gamma}=50$ would be necessary. Furthermore, in such
cases, it would be imperative to expand the size of the computational domain
to consider the extensive stretching of the droplet over longer time
intervals. These simulation setups require increased computational resources
and go beyond the scope of the present work. Our simulations were terminated
when the droplets approached the outer boundaries of the computational domain.
The simulations with high Deborah numbers did not reach their maximum droplet
deformation within the simulation time due to these constraints. For the
highest Capillary number $\operatorname{Ca}=0.6$, the drop deformation is so
intense that a larger domain size becomes necessary, even for the Newtonian
drop.

((a)) VN.

((b)) NV.

Figure 16: Three-dimensional droplet deformation as a function of time for NN,
VN, and NV systems for varying Capillary numbers and Deborah numbers. Some of
the simulations were terminated before $t\dot{\gamma}=50$ because: (i) at
$\operatorname{Ca}=0.6$ the drops are strongly deformed and the computational
domain’s length constraint with periodic boundaries prevents running the
simulations any further without risk of drop coalescence, or (ii) a steady-
state deformation was reached.

Figure 17 shows the transient orientation angles of the drops in 3D
simulations. In the VN system, the orientation angles in 3D are slightly
smaller than their 2D counterparts, suggesting a more pronounced alignment
with the flow direction in three dimensions. This discrepancy between 2D and
3D orientations becomes more pronounced as the Capillary number increases.
Specifically, at a Capillary number of $\operatorname{Ca}=0.6$, the
orientation angles in 3D are about $10$ degrees smaller than those in 2D. The
periodic droplet oscillations, as observed in the deformation parameter, are
also visible in the transient orientation angles. Therefore, these
oscillations affect both the deformation and orientation of the drop.
Regarding the variation of the Deborah number, a similar qualitative behavior
to that in 2D drops is observed: relative to the Newtonian reference state,
viscoelastic drops tend to exhibit stronger resistance to aligning with the
flow direction, resulting in larger orientation angles. Starting from the
Newtonian reference state, an increase in the Deborah number results in larger
orientation angles with a maximum $\theta$ at about $\operatorname{De}=2$. As
the Deborah number is further increased, the orientation angles begin to
decrease, indicating a better alignment with the flow direction. Figure 17(b)
displays the orientation angles in the NV system. With respect to the
Newtonian reference matrix fluid an increase in the Deborah number leads to
reduced drop orientation angles, meaning an improved alignment with the flow
direction. The influence of the Deborah number on orientation angles is
considerably more pronounced in the NV system than in the VN setup. At
$\operatorname{De}=1$, the orientation angles are approximately $10$ degrees
smaller than those in the Newtonian matrix fluid in the
$\operatorname{Ca}=0.24$ and $\operatorname{Ca}=0.4$ cases. A further increase
of the Deborah number results in $\theta$ values below $10$ degrees,
indicating a very strong impact of the flow in the viscoelastic matrix on the
drop orientation. At $\operatorname{Ca}=0.6$, the resistance of the drop is so
small that the orientation angles approach values below $10$ degrees for all
Deborah numbers.

((a)) VN.

((b)) NV.

Figure 17: Three-dimensional droplet orientation angle as a function of time
for NN, VN, and NV systems for varying Capillary numbers and Deborah numbers.

Figure 18 depicts the contours of the drop shapes, captured in a vertical
cross-section through the center of the three-dimensional drops at a fixed
time, $t=10\,\dot{\gamma}^{-1}$. At small Deborah and Capillary numbers the
drops exhibit nearly elliptical shapes. With increasing Capillary numbers, the
drops change into more elongated shapes. At $\operatorname{Ca}=0.6$, we
observe the characteristic dumbbell shapes 12, 58. Note that the drop shapes
at $t=10\,\dot{\gamma}^{-1}$ have not reached a steady-state. Actually, for
certain configurations a steady-state may not exist.

((a)) VN, $\operatorname{Ca}=0.24$.

((b)) NV, $\operatorname{Ca}=0.24$.

((c)) VN, $\operatorname{Ca}=0.4$.

((d)) NV, $\operatorname{Ca}=0.4$.

((e)) VN, $Ca=0.6$.

((f)) NV, $Ca=0.6$.

Figure 18: Contours of the drop shapes in a vertical cutting plane through the
center of the drop in 3D at $t=10\,\dot{\gamma}^{-1}$.

##  6 Summary and conclusion

A new geometric Volume-of-Fluid framework has been presented for simulating
viscoelastic fluids on 3D general unstructured finite volume meshes, building
on the ‘plicRDF-isoAdvector’ method in OpenFOAM and the modular platform
‘DeboRheo’ for computational rheology. This framework employs a change-of-
variable formulation in the constitutive equations for viscoelastic fluids,
which significantly enhances stability at moderate and high Deborah numbers. A
comprehensive mathematical modeling approach for a broad range of rheological
constitutive equations has been established within the VOF framework,
utilizing conditional volume averaging and appropriate closure assumptions.
Additionally, improvements in the reconstruction of viscosities at cell-faces
have been implemented to reduce local discretization errors and increase the
accuracy of the simulations. An iterative segregated solution algorithm is
proposed to solve the coupled system of equations.

The VOF framework has been extensively validated for transient single droplet
deformation in shear flows in 2D and 3D. The Oldroyd-B and Giesekus models are
considered, with the latter being evaluated across a range of mobility
parameters. Our results match well with previous numerical studies in 2D 12,
49 and in 3D 13, 54, 51, 55, 15, 56. The present study extends the
investigation to Deborah numbers up to $\operatorname{De}=16$, which, to the
best of our knowledge, is beyond the range explored in previous numerical
research. No stability issues were observed with the square root conformation
representation within this Deborah number range. Key observations from this
study include how various combinations of Capillary and Deborah numbers
correlate with the droplet deformation and orientation, as well as potential
periodic droplet oscillations. At higher Deborah numbers, damped droplet
oscillations are observed at lower Capillary numbers than those reported in
previous numerical works, where limited ranges of Deborah numbers were
considered. Despite the stability of the simulations, investigating higher
Deborah numbers requires considerably more computational resources, as both
the necessary domain size and the time scale for detailed droplet dynamics
analysis increase substantially. Nevertheless, it would be valuable in future
work to conduct a more systematic examination of critical combinations of
Capillary and Deborah numbers beyond which steady-state deformation parameters
and orientation angles can no longer be sustained.

## Acknowledgement

This work was partly supported from the German Research Foundation (DFG) –
Project-ID 265191195 – SFB 1194. This support is gratefully acknowledged.
Calculations for this research were conducted on the Lichtenberg high
performance computer of the Technische Universität Darmstadt.

##  Appendix A Definitions of the tensors $\mathbf{B}$ and $\mathbf{\Omega}$

By employing the diagonalization
$\mathbf{C}=\mathbf{Q}{\scriptscriptstyle\mathchoice{\mskip 3.0mu\lower
0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip 3.0mu\lower
0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\mathbf{D}{\scriptscriptstyle\mathchoice{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\mathbf{Q}\\!^{\top}$, the tensor
$\mathbf{B}$ is computed as
$\mathbf{B}=\mathbf{Q}{\scriptscriptstyle\mathchoice{\mskip 3.0mu\lower
0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip 3.0mu\lower
0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\tilde{\mathbf{B}}{\scriptscriptstyle\mathchoice{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\mathbf{Q}\\!^{\top}$. In this
formulation, the elements of the diagonal tensor $\tilde{\mathbf{B}}$ are
given by the tensor
$\mathbf{L}=\nabla\mathbf{u}\\!^{\top}=\mathbf{Q}{\scriptscriptstyle\mathchoice{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\tilde{\mathbf{L}}{\scriptscriptstyle\mathchoice{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\mathbf{Q}\\!^{\top}$ as
$\tilde{b}_{ii}=\tilde{l}_{ii}$. The tensor $\mathbf{\Omega}$ is computed as
$\mathbf{\Omega}=\mathbf{Q}{\scriptscriptstyle\mathchoice{\mskip 3.0mu\lower
0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip 3.0mu\lower
0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\tilde{\mathbf{\Omega}}{\scriptscriptstyle\mathchoice{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\mskip
3.0mu\lower 0.43054pt\hbox{\scalebox{3.0}{$\cdot$}}\mskip 3.0mu}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}{\lower
0.64583pt\hbox{\scalebox{1.5}{$\cdot$}}}}\mathbf{Q}\\!^{\top}$, where
$\tilde{\mathbf{\Omega}}$ has zero diagonal entries $\tilde{\omega}_{ii}=0$
and its off-diagonal elements are given by

$\tilde{\omega}_{{{ij},\;{i\neq j}}}=\frac{d_{ii}{\tilde{l}_{{ij},\;{i\neq
j}}}+d_{jj}{\tilde{l}_{{ji},\;{j\neq i}}}}{d_{jj}-d_{ii}},\ i,j=1,2,3.$ (39)

A comprehensive description on the local decomposition of the deformation
terms in the convective derivative can be found in 10.

##  Appendix B Benchmark data for drop deformation and orientation

Table 5: Comparison of (3D) drop deformation $D$ and orientation angle $\theta$ at time $t=3\leavevmode\nobreak\ \dot{\gamma}^{-1}$, $t=10\leavevmode\nobreak\ \dot{\gamma}^{-1}$ and $t_{f}$ (the final instance of simulations, corresponding to the results presented in section 5). NV systems for $\operatorname{Ca}=0.24$, $\operatorname{Ca}=0.4$ and $\operatorname{Ca}=0.6$ and varying Deborah numbers. $(^{*})$ The deformation values do not correspond to the equilibrium state values. | $t=3\leavevmode\nobreak\ \dot{\gamma}^{-1}$ |  | $t=10\leavevmode\nobreak\ \dot{\gamma}^{-1}$ |  | $t_{f}$
---|---|---|---|---|---
| $D$ | $\theta$ |  | $D$ | $\theta$ |  | $D$ | $\theta$
NV , $Ca=0.24$ , $De=0$ | 0.2805 | 29.60 |  | 0.2741 | 29.93 |  | 0.2746 | 29.95
NV , $Ca=0.24$ , $De=0.4$ | 0.2637 | 25.62 |  | 0.2649 | 24.16 |  | 0.2663 | 25.57
NV , $Ca=0.24$ , $De=1$ | 0.2529 | 21.53 |  | 0.2726 | 16.71 |  | 0.2714 | 18.03
NV , $Ca=0.24$ , $De=2$ | 0.2249 | 20.90 |  | 0.2844 | 7.43 |  | 0.2904 | 8.47
NV , $Ca=0.24$ , $De=4$ | 0.1897 | 23.59 |  | 0.2696 | 5.02 |  | 0.3383 | 4.86
NV∗ , $Ca=0.24$ , $De=8$ | 0.1633 | 27.90 |  | 0.2354 | 2.64 |  | 0.3755 | 4.49
NV∗ , $Ca=0.24$ , $De=16$ | 0.1481 | 28.51 |  | 0.1691 | 16.43 |  | 0.4551 | 3.29
NV , $Ca=0.4$ , $De=0$ | 0.4543 | 22.68 |  | 0.5088 | 18.40 |  | 0.5081 | 18.42
NV , $Ca=0.4$ , $De=0.4$ | 0.4298 | 21.14 |  | 0.4769 | 16.25 |  | 0.4731 | 16.35
NV , $Ca=0.4$ , $De=1$ | 0.4083 | 19.40 |  | 0.5079 | 9.68 |  | 0.5100 | 9.75
NV∗ , $Ca=0.4$ , $De=2$ | 0.3669 | 18.28 |  | 0.5096 | 7.79 |  | 0.5625 | 6.06
NV∗ , $Ca=0.4$ , $De=4$ | 0.3146 | 19.79 |  | 0.4582 | 5.68 |  | 0.6010 | 4.18
NV∗ , $Ca=0.4$ , $De=8$ | 0.2728 | 22.49 |  | 0.3889 | 4.17 |  | 0.618 | 2.57
NV∗ , $Ca=0.4$ , $De=16$ | 0.2466 | 23.43 |  | 0.3050 | 2.43 |  | 0.6618 | 2.31
NV∗ , $Ca=0.6$ , $De=0$ | 0.572 | 20.55 |  | 0.7771 | 10.25 |  | 0.8512 | 6.51
NV∗ , $Ca=0.6$ , $De=0.4$ | 0.5420 | 19.62 |  | 0.7561 | 9.40 |  | 0.8572 | 4.70
NV∗ , $Ca=0.6$ , $De=1$ | 0.5176 | 18.36 |  | 0.7714 | 7.36 |  | 0.856 | 4.68
NV∗ , $Ca=0.6$ , $De=2$ | 0.4741 | 18.68 |  | 0.7695 | 6.83 |  | 0.8500 | 7.73
NV∗ , $Ca=0.6$ , $De=4$ | 0.4188 | 19.35 |  | 0.7140 | 6.15 |  | 0.8508 | 3.96
NV∗ , $Ca=0.6$ , $De=8$ | 0.3725 | 20.82 |  | 0.6028 | 5.04 |  | 0.8442 | 2.97
NV∗ , $Ca=0.6$ , $De=16$ | 0.3433 | 21.77 |  | 0.4766 | 4.62 |  | 0.8345 | 2.28
Table 6: Comparison of (3D) drop deformation $D$ and orientation angle $\theta$ at time $t=3\leavevmode\nobreak\ \dot{\gamma}^{-1}$, $t=10\leavevmode\nobreak\ \dot{\gamma}^{-1}$ and $t_{f}$ (the final instance of simulations, corresponding to the results presented in section 5). VN systems for $\operatorname{Ca}=0.24$, $\operatorname{Ca}=0.4$ and $\operatorname{Ca}=0.6$ and varying Deborah numbers. $(^{*})$ The deformation values do not correspond to the equilibrium state values. | $t=3\leavevmode\nobreak\ \dot{\gamma}^{-1}$ |  | $t=10\leavevmode\nobreak\ \dot{\gamma}^{-1}$ |  | $t_{f}$
---|---|---|---|---|---
| $D$ | $\theta$ |  | $D$ | $\theta$ |  | $D$ | $\theta$
VN , $Ca=0.24$ , $De=0$ | 0.2805 | 29.60 |  | 0.2741 | 29.93 |  | 0.2746 | 29.95
VN , $Ca=0.24$ , $De=0.4$ | 0.2578 | 32.43 |  | 0.2549 | 31.03 |  | 0.2548 | 31.01
VN , $Ca=0.24$ , $De=1$ | 0.2550 | 34.27 |  | 0.2473 | 32.99 |  | 0.2464 | 32.68
VN , $Ca=0.24$ , $De=2$ | 0.2577 | 34.08 |  | 0.2473 | 33.00 |  | 0.2467 | 33.02
VN , $Ca=0.24$ , $De=4$ | 0.2604 | 33.24 |  | 0.2516 | 34.34 |  | 0.250 | 34.37
VN , $Ca=0.24$ , $De=8$ | 0.2633 | 33.68 |  | 0.2556 | 34.15 |  | 0.2512 | 34.42
VN , $Ca=0.24$ , $De=16$ | 0.2652 | 33.56 |  | 0.2588 | 33.99 |  | 0.2519 | 34.44
VN , $Ca=0.4$ , $De=0$ | 0.4543 | 22.68 |  | 0.5088 | 18.40 |  | 0.5081 | 18.42
VN , $Ca=0.4$ , $De=0.4$ | 0.4385 | 24.62 |  | 0.4491 | 21.74 |  | 0.4478 | 21.79
VN , $Ca=0.4$ , $De=1$ | 0.4388 | 25.84 |  | 0.4198 | 24.16 |  | 0.4205 | 24.08
VN , $Ca=0.4$ , $De=2$ | 0.4465 | 26.48 |  | 0.4201 | 26.41 |  | 0.4235 | 25.14
VN , $Ca=0.4$ , $De=4$ | 0.4554 | 26.07 |  | 0.4414 | 25.38 |  | 0.4447 | 25.25
VN , $Ca=0.4$ , $De=8$ | 0.4622 | 25.70 |  | 0.4663 | 24.25 |  | 0.4624 | 24.45
VN , $Ca=0.4$ , $De=16$ | 0.4662 | 25.48 |  | 0.4878 | 23.31 |  | 0.4878 | 23.94
VN∗ , $Ca=0.6$ , $De=0$ | 0.572 | 20.55 |  | 0.7771 | 10.25 |  | 0.8512 | 6.51
VN∗ , $Ca=0.6$ , $De=0.4$ | 0.5516 | 21.57 |  | 0.7374 | 11.20 |  | 0.8577 | 5.42
VN , $Ca=0.6$ , $De=1$ | 0.5550 | 22.42 |  | 0.6700 | 13.08 |  | 0.8030 | 7.62
VN , $Ca=0.6$ , $De=2$ | 0.5670 | 22.78 |  | 0.6984 | 13.70 |  | 0.6276 | 15.43
VN∗ , $Ca=0.6$ , $De=4$ | 0.5803 | 22.99 |  | 0.7375 | 12.70 |  | 0.6748 | 14.32
VN∗ , $Ca=0.6$ , $De=8$ | 0.5904 | 22.45 |  | 0.7807 | 11.49 |  | 0.7527 | 9.83
VN∗ , $Ca=0.6$ , $De=16$ | 0.6018 | 21.83 |  | 0.8108 | 10.77 |  | 0.8492 | 7.85

## References

  * 1 Fischer P, Erni P. Emulsion drops in external flow fields – The role of liquid interfaces. Curr. Opin. Colloid Interface Sci. 2007; 12(4-5): 196–205.  doi: 10.1016/j.cocis.2007.07.014
  * 2 Windhab E, Dressler M, Feigl K, Fischer P, Megias-Alguacil D. Emulsion processing – from single drop deformation to design of complex processes and products. Chem. Eng. Sci. 2005; 60(8-9): 2101–2113.  doi: 10.1016/j.ces.2004.12.003
  * 3 Taylor G. The formation of emulsions in definable fields of flow. Proc. R. Soc. London, Ser. A. 1934; 146(858): 501–523.
  * 4 Guido S, Greco F. Dynamics of a liquid drop in a flowing immiscible liquid. In: British Society of Rheology, Aberystwyth, United Kingdom. 2004 (pp. 99–142).
  * 5 Van Puyvelde P, Vananroye A, Cardinaels R, Moldenaers P. Review on morphology development of immiscible blends in confined shear flow. Polymer 2008; 49(25): 5363–5372.  doi: 10.1016/j.polymer.2008.08.055
  * 6 Guido S. Shear-induced droplet deformation: Effects of confined geometry and viscoelasticity. Curr. Opin. Colloid Interface Sci. 2011; 16(1): 61–70.  doi: 10.1016/j.cocis.2010.12.001
  * 7 Tucker III CL, Moldenaers P. Microstructural evolution in polymer blends. Annu. Rev. Fluid Mech. 2002; 34(1): 177–210.  doi: 10.1146/annurev.fluid.34.082301.144051
  * 8 Joseph D, Renardy M, Saut J. Hyperbolicity and change of type in the flow of viscoelastic fluids. Arch. Ration. Mech. Anal. 1985; 87(3): 213–251.  doi: 10.1007/BF00250725
  * 9 Keunings R. On the High Weissenberg Number Problem. J. Non-Newton. Fluid Mech. 1986; 20: 209–226.  doi: 10.1016/0377-0257(86)80022-2
  * 10 Niethammer M, Marschall H, Kunkelmann C, Bothe D. A numerical stabilization framework for viscoelastic fluid flow using the finite volume method on general unstructured meshes. Int. J. Numer. Methods Fluids 2018; 86(2): 131–166.  doi: 10.1002/fld.4411
  * 11 Hirt C, Nichols B. Volume of fluid (VOF) method for the dynamics of free boundaries. J. Comput. Phys. 1981; 39(1): 201–225.  doi: 10.1016/0021-9991(81)90145-5
  * 12 Chinyoka T, Renardy Y, Renardy M, Khismatullin D. Two-dimensional study of drop deformation under simple shear for Oldroyd-B liquids. J. Non-Newtonian Fluid Mech. 2005; 130(1): 45–56.  doi: 10.1016/j.jnnfm.2005.07.005
  * 13 Khismatullin D, Renardy Y, Renardy M. Development and implementation of VOF-PROST for 3D viscoelastic liquid–liquid simulations. J. Non-Newtonian Fluid Mech. 2006; 140(1-3): 120–131.  doi: 10.1016/j.jnnfm.2006.02.013
  * 14 Afkhami S, Yue P, Renardy Y. A comparison of viscoelastic stress wakes for two-dimensional and three-dimensional Newtonian drop deformations in a viscoelastic matrix under shear. Phys. Fluids 2009; 21(7): 1–8.  doi: 10.1063/1.3182830
  * 15 Verhulst K, Cardinaels R, Moldenaers P, Renardy Y, Afkhami S. Influence of viscoelasticity on drop deformation and orientation in shear flow. Part 1. Stationary states. J. Non-Newtonian Fluid Mech. 2009; 156(1-2): 29–43.  doi: 10.1016/j.jnnfm.2008.06.007
  * 16 Figueiredo R, Oishi C, Afonso A, Tasso I, Cuminato J. A two-phase solver for complex fluids: Studies of the Weissenberg effect. Int. J. Multiphase Flow 2016; 84: 98–115.  doi: 10.1016/j.ijmultiphaseflow.2016.04.014
  * 17 López-Herrera J, Popinet S, Castrejón-Pita A. An adaptive solver for viscoelastic incompressible two-phase problems applied to the study of the splashing of weakly viscoelastic droplets. J. Non-Newtonian Fluid Mech. 2019; 264(November 2018): 144–158.  doi: 10.1016/j.jnnfm.2018.10.012
  * 18 Niethammer M, Brenn G, Marschall H, Bothe D. An extended volume of fluid method and its application to single bubbles rising in a viscoelastic liquid. J. Comput. Phys. 2019; 387: 326–355.  doi: 10.1016/j.jcp.2019.02.021
  * 19 Anderson T, Jackson R. Fluid mechanical description of fluidized beds. Equations of motion. Ind. Eng. Chem. Fund. 1967; 6(4): 527–539.  doi: 10.1021/i160024a007
  * 20 Slattery J. Flow of viscoelastic fluids through porous media. AIChE J. 1967; 13(6): 1066–1071.
  * 21 Whitaker S. Diffusion and dispersion in porous media. AIChE J. 1967; 13(3): 420–427.  doi: 10.1002/aic.690130308
  * 22 Bothe D, Niethammer M, Pilz C, Brenn G. On the molecular mechanism behind the bubble rise velocity jump discontinuity in viscoelastic liquids. J. Non-Newton. Fluid Mech. 2022; 302: 104748.  doi: 10.1016/j.jnnfm.2022.104748
  * 23 Niethammer M, Marschall H, Bothe D. Robust Direct Numerical Simulation of Viscoelastic Flows. Chemie Ing. Tech. 2019; 91(4): 522–528.  doi: 10.1002/cite.201800108
  * 24 Scheufler H, Roenby J. Accurate and efficient surface reconstruction from volume fraction data on general meshes. J. Comput. Phys. 2019; 383: 1–23.  doi: 10.1016/j.jcp.2019.01.009
  * 25 Roenby J, Bredmose H, Jasak H. A computational method for sharp interface advection. R. Soc. Open Sci. 2016; 3(11): 160405.  doi: 10.1098/rsos.160405
  * 26 Scheufler H, Roenby J. TwoPhaseFlow: An OpenFOAM based framework for development of two phase flow solvers. 2021.  doi: 10.48550/arXiv.2103.00870
  * 27 Scheufler H. TwoPhaseFlow. https://github.com/DLR-RY/TwoPhaseFlow; 2022\.
  * 28 Tryggvason G, Esmaeeli A, Al-Rawahi N. Direct numerical simulations of flows with phase change. Comput. & Structures 2005; 83(6-7): 445–453.  doi: 10.1016/j.compstruc.2004.05.021
  * 29 Zalesak S. Fully multidimensional flux-corrected transport algorithms for fluids. J. Comput. Phys. 1979; 31(3): 335–362.
  * 30 Muzaferija S. A two-fluid Navier-Stokes solver to simulate water entry. In: National Academy Press. ; 1999: 638–651.
  * 31 Ubbink O, Issa R. A Method for Capturing Sharp Fluid Interfaces on Arbitrary Meshes. J. Comput. Phys. 1999; 153(1): 26–50.
  * 32 Marić T, Kothe DB, Bothe D. Unstructured un-split geometrical Volume-of-Fluid methods – A review. J. Comput. Phys. 2020; 420: 109695\.  doi: 10.1016/j.jcp.2020.109695
  * 33 Brackbill J, Kothe D, Zemach C. A continuum method for modeling surface tension. J. Comput. Phys. 1992; 100(2): 335–354.  doi: https://doi.org/10.1016/0021-9991(92)90240-Y
  * 34 Liu J, Tolle T, Bothe D, Marić T. An unstructured finite-volume level set/front tracking method for two-phase flows with large density-ratios. J. Comput. Phys. 2023; 493: 112426.
  * 35 Bird R, Armstrong R, Hassager O. Dynamics of Polymeric Liquids. 1.  John Wiley & Sons, New York .  1977\.
  * 36 Larson R. Constitutive Equations for Polymer Melts and Solutions (Butterworths Series in Chemical Engineering).  Butterworths, Boston .  1988\.
  * 37 Oldroyd J. On the Formulation of Rheological Equations of State. Proc. R. Soc. Lond. A 1950; 200(1063): 523–541.  doi: 10.1098/rspa.1950.0035
  * 38 Rivlin R, Ericksen J. Stress-Deformation Relations for Isotropic Materials. J. Rat. Mech. Anal. 1955; 4(2): 323–425.  doi: 10.1007/978-1-4612-2416-7_61
  * 39 Rivlin R. Further Remarks on the Stress-Deformation Relations for Isotropic Materials. J. Rat. Mech. Anal. 1955; 4(5): 681–702.  doi: 10.1007/978-1-4612-2416-7_62
  * 40 Fattal R, Kupferman R. Constitutive laws for the matrix-logarithm of the conformation tensor. J. Non-Newton. Fluid Mech. 2004; 123(2-3): 281–285.  doi: 10.1016/j.jnnfm.2004.08.008
  * 41 Balci N, Thomases B, Renardy M, Doering C. Symmetric factorization of the conformation tensor in viscoelastic fluid models. J. Non-Newton. Fluid Mech. 2011; 166(11): 546–553.  doi: 10.1016/j.jnnfm.2011.02.008
  * 42 Afonso A, Pinho F, Alves M. The kernel-conformation constitutive laws. J. Non-Newton. Fluid Mech. 2012; 167-168: 30-37.  doi: 10.1016/j.jnnfm.2011.09.008
  * 43 Jasak H. Error analysis and estimation for the finite volume method with applications to fluid flows. PhD thesis. Imperial College London (University of London), ; 1996.
  * 44 Marić T. Iterative Volume-of-Fluid interface positioning in general polyhedrons with Consecutive Cubic Spline interpolation. Journal of Computational Physics: X 2021; 11: 100093.  doi: doi.org/10.1016/j.jcpx.2021.100093
  * 45 Kromer J, Bothe D. Face-based Volume-of-Fluid interface positioning in arbitrary polyhedra. J. Comput. Phys. 2022; 449: 110776.  doi: 10.1016/j.jcp.2021.110776
  * 46 Marić T, Marschall H, Bothe D. An enhanced un-split face-vertex flux-based VoF method. J. Comput. Phys. 2018; 371: 967–993.  doi: 10.1016/j.jcp.2018.03.048
  * 47 Tolle T, Bothe D, Marić T. SAAMPLE: A segregated accuracy-driven algorithm for multiphase pressure-linked equations. Comput. Fluids 2020; 200: 104450.  doi: 10.1016/j.compfluid.2020.104450
  * 48 Pillapakkam S, Singh P. A Level-Set Method for Computing Solutions to Viscoelastic Two-Phase Flow. J. Comput. Phys. 2001; 174(2): 552–578.  doi: 10.1006/jcph.2001.6927
  * 49 Chung C, Hulsen M, Kim JM, Ahn K, Lee S. Numerical study on the effect of viscoelasticity on drop deformation in simple shear and 5:1:5 planar contraction/expansion microchannel. J. Non-Newtonian Fluid Mech. 2008; 155(1-2): 80–93.  doi: 10.1016/j.jnnfm.2008.06.002
  * 50 Yue P, Feng J, Liu C, Shen J. Viscoelastic effects on drop deformation in steady shear. J. Fluid Mech. 2005; 540: 427–437.  doi: 10.1017/S0022112005006166
  * 51 Aggarwal N, Sarkar K. Effects of matrix viscoelasticity on viscous and viscoelastic drop deformation in a shear flow. J. Fluid Mech. 2008; 601: 63–84.  doi: 10.1017/S0022112008000451
  * 52 Zhou H, Pozrikidis C. The flow of suspensions in channels: Single files of drops. Phys. Fluids A, Fluid Dyn. 1993; 5(2): 311–324.  doi: 10.1063/1.858893
  * 53 Chung C, Hulsen MA, Ahn KH, Lee SJ. Viscoelastic drop deformation in simple shear flow investigated by the front tracking method. In: . I of WIT Transactions on Engineering Sciences, Vol 56. WIT Press; 2007: 215–224
  * 54 Aggarwal N, Sarkar K. Deformation and breakup of a viscoelastic drop in a Newtonian matrix under steady shear. J. Fluid Mech. 2007; 584: 1–21.  doi: 10.1017/S0022112007006210
  * 55 Mukherjee S, Sarkar K. Effects of viscosity ratio on deformation of a viscoelastic drop in a Newtonian matrix under steady shear. J. Non-Newtonian Fluid Mech. 2009; 160(2-3): 104–112.  doi: 10.1016/j.jnnfm.2009.03.007
  * 56 Verhulst K, Cardinaels R, Moldenaers P, Afkhami S, Renardy Y. Influence of viscoelasticity on drop deformation and orientation in shear flow. Part 2. Dynamics. J. Non-Newtonian Fluid Mech. 2009; 156(1-2): 44–57.  doi: 10.1016/j.jnnfm.2008.10.003
  * 57 Luo Z, Bai B. Dynamics of capsules enclosing viscoelastic fluid in simple shear flow. J. Fluid Mech. 2018; 840: 656–687.  doi: 10.1017/jfm.2018.88
  * 58 Sheth KS, Pozrikidis C. Effects of inertia on the deformation of liquid drops in simple shear flow. Comput. Fluids 1995; 24(2): 101–119.  doi: 10.1016/0045-7930(94)00025-T
  * 59 Cristini V, Guido S, Alfani A, Bławzdziewicz J, Loewenberg M. Drop breakup and fragment size distribution in shear flow. J. Rheol. 2003; 47(5): 1283–1298.  doi: 10.1122/1.1603240
  * 60 Guido S, Simeone M, Greco F. Effects of matrix viscoelasticity on drop deformation in dilute polymer blends under slow shear flow. Polymer 2003; 44(2): 467–471.  doi: 10.1016/S0032-3861(02)00763-2
  * 61 Muzaferija S. Computation of free surface flows using interface-tracking and interface-capturing methods. in Mahrenholtz O., Markiewicz M. (eds.), Nonlinear water wave interaction, Computational Mechanics Publications, Southampton 1998\.
  * 62 Popinet S. Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries. J. Comput. Phys. 2003; 190(2): 572–600.  doi: 10.1016/S0021-9991(03)00298-5
  * 63 Renardy Y, Renardy M. PROST: A Parabolic Reconstruction of Surface Tension for the Volume-of-Fluid Method. J. Comput. Phys. 2002; 183(2): 400–421.  doi: 10.1006/jcph.2002.7190
  * 64 Scardovelli R, Zaleski S. Analytical Relations Connecting Linear Interfaces and Volume Fractions in Rectangular Grids. J. Comput. Phys. 2000; 164(1): 228–237.  doi: 10.1006/jcph.2000.6567
  * 65 Scardovelli R, Zaleski S. Interface reconstruction with least-square fit and split Eulerian-Lagrangian advection. Int. J. Numer. Methods Fluids 2003; 41(3): 251–274.  doi: 10.1002/fld.431
  * 66 Slattery J. Advanced Transport Phenomena.  Cambridge University Press .  1999\.
  * 67 Whitaker S. The method of volume averaging.  Kluwer Academic Publishers, Dordrecht, The Netherlands .  1999\.