### Computational Techniques on Two Phase Flow

### Authors: Charitini Stavropoulou, Erik Blom, Gunnar Dufwa Supervisors: Gunilla Kreiss, Murtazo Nazarov

### Project in Computational Science: Report February 8, 2021

## PROJECT REPORT

### Contents

1 Introduction 2

1.1 Two - Phase Flow . . . 2

1.2 Level Set Method . . . 2

1.3 Normal and Curvature . . . 3

2 Overview of the Literature and Principles 3 2.1 Level-set equation . . . 3

2.2 Principle of the artificial - viscocity - based compression method . . . 4

2.3 Incompressible Two-phase Flow . . . 4

2.3.1 Incompressible Flow . . . 4

2.3.2 Model of Two-Phase flow . . . 5

2.4 Scalar Conservation Laws . . . 5

2.5 A residual-based artificial viscocity method (RV) . . . 6

2.6 A residual-based artificial viscocity compression method (Antidiffusion) . . . 7

2.7 Runge -Kutta Method . . . 7

2.8 Incremental Pressure Correction Scheme (IPCS) . . . 8

2.9 Normal and Curvature . . . 8

2.9.1 Elliptic Projection . . . 9

2.9.2 Error Method for Curvature . . . 9

3 Results 10 3.1 Residual-based artificial viscosity method (RV) . . . 10

3.1.1 Smooth Initial Data . . . 10

3.1.2 Non Smooth Initial Data - The Level set Function . . . 11

3.1.3 Artificial Viscocity over Time . . . 13

3.2 A residual-based artificial viscocity compression method - Anti-diffusion . . . 13

3.3 Normal Vector . . . 14

3.4 Curvature . . . 15

3.5 Incompressible Navier stokes . . . 17

3.5.1 Channel Flow - Constant Inflow Profile . . . 18

3.5.2 Channel Flow - Variable Density - Constant Inflow profile . . . 18

3.5.3 Channel Flow - non Constant Inflow profile . . . 18

3.5.4 Channel Flow - Obstacle - non Constant Inflow profile . . . 18

3.5.5 Flow-past a cylinder . . . 19

4 Discussion 21

5 Conclusions 24

Abstract

In this work, we study a conservative level-set method [1] for two-phase flow. In parti-
cular, it is investigated how this level-set approach performs in combination with a residual-
based second order artificial viscosity stabilization technique for finite element approxima-
tions. It is proved that the convergence rate of the finite element approximation is close
to 0.5, measured in the L^{2}-norm. We use a stabilization method for the evaluation of
functionals, fundamental to two-phase flow - such as the curvature and the normal of the
interface. In addition, an a-posteriori error-estimate for the computational curvature is
proposed and used, however, subsequent results on the curvature convergence are inconclu-
sive. The simulations are performed with a known velocity field as well as coupled with the
incompressible Navier-Stokes equations for various settings that involve a channel flow.

### 1 Introduction

### 1.1 Two - Phase Flow

The concept of multiphase flow entails a rich field of study with a wide range of applications in industry processing, aquaculture and biomedicine to name only a few examples [2], [3]. A flow is considered multiphase when it consists of multiple components in distinct phases [4].

Multiphase flows can be categorized by the thermodynamic phases of the components involved, by the degree of mixing between components and by the topology of the flow. The richness and complexity in a complete description of multiphase flows together with its variety of applications has demanded the development of new simplifications and accurate numerical methods [2]. A complete description of multiphase flows (in particular immiscible flows where interfaces between fluids are well-defined) would have to include the effects of surface tension in the fluid mechanical model, which in turn requires accurate knowledge about the curvature of the interfaces between fluid components. It is, however, difficult to accurately evaluate the curvature from a numerical representation of an interface, and it is possible that a proper method for evaluating the curvature will depend upon the particular interface representation in consideration. This article considers simplified two-phase flow systems to analyse the accuracy of the interface curvature using a conservative level set method [1] to represent the interface.

### 1.2 Level Set Method

The level set method is a method for an implicit description of moving fronts. The basic idea is that the front location is given as a constant level set of an auxiliary field defined over the domain of interest. The level set method is based on the concept of implicit surfaces. To define an implicit surface, we construct a function φ(x, y, z) such that when

φ(x, y, z) = 0.5 (1)

then the point with coordinates (x, y, z) lies on the surface. Notice that φ is of more general form in order to describe a two-dimensional surface like the air-water interface. Specifically, φ is a three-dimensional function that describes a volume that contains the desired surface.

When φ > 0.5, the point (x, y, z) lies outside the surface while when φ < 0.5, the point (x, y, z) lies inside the surface. Only when φ = 0.5 the point (x, y, z) lies on the surface. Eq.(1) may appear to be abstract at first, however, it offers a simple means for describing a two-dimensional surface. For example, we can determine whether a point (x, y, z) is inside or outside a fluid body by simply checking the value of the implicit surface function φ [5].

The traditional level set method uses a signed distance function and φ = 0 represents the interface. It is useful for handling complex topological dynamics but suffers from a lack of mass conservation. In [1] a novel conservative level set technique is presented using a new compression strategy. Conservative level set methods have been the focus of recent research.

Techniques for enforcing mass conservation in level set methods differ depending on the context of the method. The assembly of the so-called antidiffusion term is of main interest in [1]. The proposed construction is straight-forward to implement. In addition, the proposed method does not require solving a re-initialization or interface reconstruction subproblem (see, e.g. Fedkiw et al. [6], Olsson and Kreiss [7] ).

### 1.3 Normal and Curvature

The implementation of an anti-diffusion construction in [1] raises questions upon the computa- tion of functionals. The task is to investigate how accurate one can compute the normal and the curvature of the interface in a setting when the interface is advected by a known velocity field.

If an implicit surface function φ is considered, provided that φ is differentiable, then the com- putation of the unit normal vector n and the curvature κ at any point of φ is straightforward.

The gradient of φ determines the normal direction, thus the unit normal is given by:

n = − ∇φ

||∇φ||_{l}2

, (2)

where the l^{2}-norm denotes the local L^{2}-norm such that the length of all n is unity. Furthermore,
the mean curvatures of the interface can be related to the divergence of the unit normal vector,
as:

κ = ∇ · n. (3)

### 2 Overview of the Literature and Principles

### 2.1 Level-set equation

The level set method typically refers to the transport of a smooth distance function whose constant isosurface is the interface to be tracked. The level set method allows to simulate flow of two immiscible fluids separated by a moving interface. If the advection velocity field is incompressible, the transport problem can be recast in conservative form which allows for the use of smoothed Heaviside level set functions, conservative discretizations, and anti-diffusion techniques for enhancing mass conservation. The basic equation for the level set method is an advection equation:

∂_{t}φ + u · ∇φ = 0. (4)

Eq. (4) is often referred to as "the level set equation", u is the flux term. The interface, I, is represented by the curve φ = 0.5. In the transition layer near the interface, φ changes smoothly form 0 to 1. Region , Ω1, where φ < 0.5 is filled with fluid 1 of density ρ1and dynamic viscocity µ1. Region , Ω2, where φ > 0.5 is filled with fluid 2 of ρ2 and µ2.

Ω1= {x|φ(x) < 0.5},
Ω_{2}= {x|φ(x) > 0.5},
I = {x|φ(x) = 0.5}.

### 2.2 Principle of the artificial - viscocity - based compression method

Artificial compression methods were proposed in [1] for use in the numerical solution of con- servation laws, of the form (4). Consider a one-dimensional, possibly stabilized, conservation equation

∂_{t}φ + ∂_{x}f (φ) = ∂_{x}(µ · ∂_{x}φ), φ(x, 0) =

(φ_{L} x < 0,

φR x > 0, (5)

with convex flux f and artificial viscocity µ > 0. The artificial compression method proceeds by augmenting (5) with an additional flux term g(φ) to produce,

∂_{t}φ + ∂_{x}(f (φ) + g(φ)) = ∂_{x}(µ · ∂_{x}φ). (6)
The function g(φ) is called an “artificial compression function”. The classical artificial compres-
sion approach to solving (6) is a two-stage method. First, (5) is solved followed by a sub-timestep
problem. An anti-diffusive formulation which can be solved in one stage is enabled by the careful
selection of the artificial compression function g(φ). In [1] g(φ) is chosen as :

g(φ) = cµ[(φ − φL)(φR− φ)^{+}] ∂xφ

|∂xφ|, (7)

where c ≥ 0 is a non-negative parameter. It can now be seen that the use of the positive part
[.]^{+} in (7) ensures that a numerical approximation to g(φ) avoids overshoots and undershoots.

Note that since µ scales like a wave speed times a length scale, and g(φ) must scale like a wave
speed times φ the parameter c must scale like c ∼ ccomp· `^{−1}where ccompis a universal constant
and ` is some length scale, ` = h, where h is the mesh size. Eq. (6) can be written as:

∂tφ + ∂x f (φ) − µ

1 − ccomp

[(φ − φL)(φR− φ)]^{+}
h|∂xφ|

^{+}

· ∂xφ

!

= 0. (8)

Hence, artificial compression can be achieved for any artificial-viscosity-based approximation method by redefining the viscosity appropriately. The generalization of (8) in higher dimension is:

∂_{t}φ + ∇ · f (φ) − µ

1 − c_{comp}[(φ − φ_{L})(φ_{R}− φ)]^{+}
h||∇φ||L^{2}

^{+}

∇φ

!

= 0. (9)

Eq (9) is the motivation of the method proposed in paper [1].

### 2.3 Incompressible Two-phase Flow

2.3.1 Incompressible Flow

The motion of fluids where the denisty is considered constant with the velocity of the fluid can be modeled as an incompressible flow. In fluid mechanics incompressible flow is described by the incompressible Navier-Stokes equation,

(ρ(∂tu + u · ∇u) = ∇ · σ(u, p) + b,

∇ · u = 0, (10)

where u is the velocity of the fluid. The right hand side of (10) represents the forces acting on the fluid. These forces includes the body force b and the symmetric Cauchy stress tensor σ(u, p):

σ(u, p) = 2µ(u) − pI, (11)

where µ is the dynamic viscosity, p is the pressure, I is the identity matrix. The (u) is the symmetric gradient,

(u) =1

2(∇u + ∇^{T}u). (12)

2.3.2 Model of Two-Phase flow

In the case of a two-phase flow, surface tension at the interface between fluids needs to be considered. The surface tension, Tη, can be included in (10), cf. eg. [8]. The two different fluids are represented by a level-set function φ, which is advected by the velocity u. Resulting in the system of equations describing two-phase flow,

ρ(∂_{t}u + u · ∇u) = ∇ · σ(u, p) + g + T_{η},

∇ · u = 0,

∂_{t}φ + u · ∇φ = 0,

(13)

where the surface tension is defined as

T_{η}= ηκnδ_{s}, (14)

considering a constant surface tension coefficient η. The functionals κ and n are the curvature and normal correspondingly, defined in Section 1.3. Since the surface tension only acts on the interface, a Dirac-delta function δsis included.

In model (13) the density ρ and the dynamic viscocity µ are defined as, (ρ = φ · ρ1+ (1 − φ)ρ2,

µ = φ · µ1+ (1 − φ)µ2, (15)

such that the density and viscosity is spatially defined where the corresponding fluid is located by the level-set function φ.

### 2.4 Scalar Conservation Laws

Let Ω be a fixed domain in R^{2}, with boundary ∂Ω over a time interval [0, T ] with initial time
zero and the final time T. We are interested in solving the following time-dependent scalar
conservation laws:

∂tφ + ∇ · f (φ) = 0, (x, t) ∈ Ω × (0, T ],
φ(x, 0) = φ_{0}(x), x ∈ Ω,

φ(x, t) = 0, ∀x ∈ ∂Ω.

(16)

Here φ represents the unknown variable, f ∈ C^{1}(R, R^{2}) is the flux term and φ0 ∈ L^{∞}(R^{2}) is
given initial data. Let 0 = t0 < t1 < ... < tN = T be a sequence of discrete time steps with
associated time intervals In = (tn−1, tn]of length kn= tn− tn−1 , n = 1, 2, ..., N

Variational Formulation

To derive a variational formulation of equation (16) we multiply (16) by a test function v, which is
assumed to vanish on the boundary ∂Ω,we choose the homogeneous Dirichlet Boundary condition
and integrate using Green’s formula (i.e., integration by parts). Further, introducing the space
V_{0}, defined by V0= {v : kvk_{L}2+ k∇vk_{L}2< ∞, v|_{∂Ω}= 0}. The weak form reads:

Find φ(x, t) ∈ V0such that:

(∂φ

∂t, v) + (∇ · f (φ), v) = 0, ∀v ∈ V0, (17)
where the norm (·, ·) is the L^{2}-norm over the entire domain Ω unless a subscript specifies another
domain.

Finite Element Approximation

Now, let Th be a triangulation of Ω and let Vh be the space of continuous piece-wise linear
functions on Th, Vh = {v : v ∈ C^{0}(Th), v|K ∈ P1(K), K ∈ Th}. Also, to satisfy the strong
boundary conditions, let Vh,0 ⊂ Vh be the subspace Vh,0 = {v ∈ Vh : v|∂Ω = 0}. With this
choice of approximation space the finite element method form of the problem is,

Find φh(x, t) ∈ Vh,0 such that:

(∂φh

∂t , v) + (∇ · f (φh), v) = 0 ∀v ∈ Vh,0. (18) Linear Advection Equation - A rotating disk

Let us consider the problem (16) in the unit disk Ω = {x : x^{2}1+x^{2}_{2}≤ 1}, and rewrite the flux term
in the non-conservative form, i.e., ∇·f(φ) := f^{0}(φ) · ∇φ. Assume that f^{0}(φ) := 2π(−x_{2}, x_{1})and
in addition the homogeneous Dirichlet boundary condition is applied on the entire boundary:

φ(x, t) = 0, ∀x ∈ ∂Ω.

Now using the non-conservative form of the flux term f^{0}(φ)the Finite element Approximation
reads:

Find φh(x, t) ∈ Vh,0 such that:

(∂φh

∂t , v) + (f^{0}(φh) · ∇φh, v) = 0 ∀v ∈ Vh,0. (19)

### 2.5 A residual-based artificial viscocity method (RV)

The numerical solution to the linear advection equation of the form (19) is susceptible to un- physical Gibbs oscillations that form in the proximity of strong gradients. The oscillations are obvious when we use a non smooth function as an initial condition to the finite element approx- imation, of the form (19). A residual-based artificial viscosity method or the RV method is a mesh-dependent approach to stabilize the finite element approximation (19). This modification adds stability by the introduction of an elliptic term which depends on the residual. The vis- cosity is added to parts of the solution where the residual is large, which usually indicates the regions where the solution is non-smooth.

For the given problem we modify the Galerkin FEM formulation in the following way.

Find φh(x, t) ∈ Vh,0 such that:

(∂φh

∂t , v) + (f^{0}(φh) · ∇φh, v) + (h∇φ, ∇v) = 0 ∀v ∈ Vh,0, (20)
where h:= h(tn)is the artificial viscosity term computed for each cell K ⊂ Th. It is based on
the residual

R(φh) = 1 kn

(φn− φn−1) + f^{0}(φn) · ∇φn,

and computed in the following way (see, e.g [9]),

n,K = min CvelhKβK, CRVh^{2}_{K} kRRVk_{∞,K}
φn− ¯φn

_{∞,Ω}

!

, (21)

where Cvel and CRV are stabilization parameters, hK is the length of the shortest edge of K, φ¯n is the space average of the solution over the domain and βK is the local element wave speed computed as presented in Equation (22).

βK ≡ kf^{0}(φn)k_{L}∞(K)= max

Ni∈K

[f_{1}^{0}(φn)^{2}+ f_{2}^{0}(φn)^{2}]^{1}^{2}(Ni)

. (22)

### 2.6 A residual-based artificial viscocity compression method (Antid- iffusion)

A residual-based artificial viscosity compression method is a mesh-dependent approach to sta- bilise the finite element approximation (19), while preserving the sharpness of the solution through an antidiffusion term (see, e.g. Section 2.2). For the given problem we modify the Galerkin FEM in the following way.

Find φh(x, t) ∈ Vh,0 such that

(∂φh

∂t , v) + (f^{0}(φh) · ∇φh, v) + (hΘ(φh)∇φh, ∇v) = 0 ∀v ∈ Vh,0,
Θ(φh) =

1 − ccomp

[(φh− φL)(φR− φh)]^{+}
h||∇φh||_{L}2

+

.

(23)

Here ccomp= 1, h is the mesh size and his defined in every triangle of the mesh as in Section 2.5.

### 2.7 Runge -Kutta Method

The explicit Runge-Kutta 4 (RK4) time-stepping scheme is used to circumvent difficulties arising from the non-linear term in the anti-diffusion method (23). The finite element scheme (23) is rearranged to isolate the time derivative

(∂φh

∂t , v) = −(f^{0}· ∇φh, v) − (hΘ∇φh, ∇v) ∀v ∈ Vh,0, (24)
where the advection velocity f^{0}(φh)and anti-diffusion term Θ(φh)is stated without dependence
on φhto emphasize that they do not change during the RK4 procedure (even though they might
change between time steps). Similarly, h does not change during the RK4 procedure. In the
Navier-Stokes implementation, the advection velocity is evaluated first at each time step and
used as constant in the subsequent level-set advection scheme.

The above equation is on a form where the right hand side at a given time step tn is only a function of the solution and the test function

(w, v) = rn(φh, v) ∀v ∈ Vh,0, (25)

where the subscript n denotes that a given time step is considered, and where w is the time
derivative of the finite element solution. To get the solution φ^{n}_{h}from φ^{n−1}_{h} , the RK4 weights are
obtained as

(w0, v) = rn(φ^{n−1}_{h} , v), (26)

(w_{1}, v) = r_{n}(φ^{n−1}_{h} + 0.5k_{n}w_{0}, v), (27)
(w2, v) = rn(φ^{n−1}_{h} + 0.5knw1, v), (28)
(w_{3}, v) = r_{n}(φ^{n−1}_{h} + k_{n}w_{2}, v), (29)
again, using the same f^{0}(φh), Θ(φh), h(tn), obtained at time step tn for each RK4 stage. The
final solution is obtained using the weights as

φ^{n}_{h}= φ^{n−1}_{h} +kn

6 (w0+ 2w1+ 2w2+ w3). (30)

### 2.8 Incremental Pressure Correction Scheme (IPCS)

The velocity and pressure from the incompressible Navier-Stokes equation (10) is numerically
computed using an efficient yet accurate IPCS. IPCS divides the computation of (10) into three
steps. Discretized in time, a tentative velocity (u^{∗}) is in the first step computed by

ρu^{∗}− u^{n}

k + ρu^{n}· ∇u^{n} = ∇ · σ(u^{n+}^{1}^{2}, p^{n}) + f^{n+1}, (31)
where k is the step in time.

The next step is to update the pressure as

∆p^{n+1}= ∆p^{n}−∇ · u^{∗}

k . (32)

The updated velocity is finally calculated by
u^{n+1}− u^{∗}

k + ∇p^{n+1}− ∇p^{n}= 0. (33)

IPCS is discretized using a finite element approximation as proposed by M. G. Larson et al. in [10].

### 2.9 Normal and Curvature

The finite element approximation of the normal (2) and curvature (3) of φ - nh, κh respectively - are obtained across the entire computational domain. Both nh and κh are obtained using the so-called elliptic projection to attain sufficient smoothness, in particular in a region of interest near the interface. A method for evaluating the curvature error on the interface is proposed in Section 2.9.2.

2.9.1 Elliptic Projection

The elliptic projection of a function is obtained by solving the following finite element approxi- mation: find uein Vh such that

(ue, v) + cellh^{2}(∇ue, ∇v) = (f, v), ∀v ∈ Vh, (34)
where f is the function in Vh which is projected, h is the global mesh size, cell is the elliptic
projection parameter. Here Vh is a space consisting of either scalar or vector-valued functions,
depending on the function f. Note that the elliptic projection is the standard finite element
projection in the limit as h goes to zero or if cell = 0. No boundary conditions are imposed in
the elliptic projection of the normal and curvature. The normal nh is obtained using elliptic
projection with f = ∇φh in an appropriate vector function space. The elliptic projection of the
gradient is normalized as,

nh= − ue

||u_{e}||_{L}2+ δ, (35)

where the additional δ-term in the denominator is much smaller than one which has the effect of neglecting the gradients of φ with magnitudes close to zero.

Subsequently, the curvature κh is evaluated using f = ∇ · nh in an appropriate scalar function space. Different elliptic projection parameters cell are considered for the evaluation of normal and curvature with regards to optimizing the error and convergence of the curvature. The goal is to find a set of elliptic projection parameters that at least results in a convergent curvature so that the conservative level set method can be implemented to include effects of, for example, surface tension.

2.9.2 Error Method for Curvature

The approach taken to analyse the error and convergence of the curvature in this report requires
some clarification as there is not only one valid method of approach. The error between analytical
curvature and numerical curvature after one rotation around the disk is evaluated as an L^{2}-norm
along the interface. The analytical solution to the interface curvature (the values of the curvature
at the interface) is known in the case of the rotating disk, e.g for an initial interface as a circle
with radius r0 then the curvature is κanalytic= 1/r_{0} along the interface. The curvature values
along the exact interface - the interface points as defined analytically - are evaluated by the
procedure explained above for the final state of the finite element approximation (23) after
one complete rotation around the disk. Let the curvature values from the final state after one
rotation be κ1. Then the L^{2}-error of the curvature is the line integral

e_{κ}=
s

Z

I

(κ_{analytic}− κ_{1})^{2}ds, (36)

along the interface I. The curvature error provides a measure of convergence when several meshes of different size are used. In the section Results, a circular initial interface with r0 = 0.25is used. The line integral is then simplified considering a coordinate transformation from cartesian to polar coordinates, with

e_{κ}=
s

Z 2π 0

r_{0}(κ_{analytic}(θ) − κ_{1}(θ))^{2}dθ, (37)
where the curvatures are considered as functions of the polar angle θ. The integral is evaluated
numerically using 600 equally distributed points along the interface circle.

### 3 Results

The Finite Element coding and simulation was done using the FENICS project [11].

### 3.1 Residual-based artificial viscosity method (RV)

In this section the linear advection equation is tested, onto the rotating disk (16) using a second order Residual-based artificial Viscosity method (see e.g, Section 2.5). The RV stabilization parameters are set to Cvel = 0.25 and CRV = 1. The test was done using two different initial conditions: ’smooth’ and ’non-smooth’ discontinuous initial data.

3.1.1 Smooth Initial Data

The smooth initial data is presented on Figure 1 and is described by

φ0(x, 0) = 1 2

1 − tanh (x_{1}− x^{0}_{1})^{2}+ (x_{2}− x^{0}_{2})^{2}

r^{2}_{0} − 1

, (38)

where (x^{0}1, x^{0}_{2}) = (0.3, 0)and r0 = 0.25. Figure 2 presents the solution at time T = 1 obtained
using the RV method derived previously (see Section 2.5). Although stabilization is not necessary
for the advection of the smooth initial data, it is included here to see that the stabilization
implementation has the appropriate convergence rate. The solution is plot for the finest mesh
size tested hmax= 0.033.

(a) (b)

Figure 1: ’Smooth’ initial data for the rotating disk problem using mesh size hmax= 0.033.

The convergence rate

In table 1 the error in the L^{2}-norm of the approximation was examined for varying parameter
hmax. Figure 3 presents a loglog-plot of hmax versus the L^{2}-norm of the errors. The slope
of the line determines the convergence rate α = 1.81. The function h^{α}max was plotted for
comparison.

(a) (b)

Figure 2: Solution to the rotating disk problem after one complete rotation (T = 1) using the ’smooth’ initial data, Residual Viscosity stabilization with mesh size hmax= 0.033.

Table 1: L^{2}-error and convergence results for the rotating disk problem, using

’smooth’ initial data, varying mesh size and Residual Viscocity (RV) stabilization.

hmax L^{2}-error convergence (q)

0.033 0.0043 1.81

0.065 0.016 1.70

0.137 0.057 -

Figure 3: A loglog -plot of hmaxversus the L^{2}-norm of the solution errors for the
rotating disk problem, using the ’smooth’ initial data. The red line is the measured
error and the blue line is h^{α} using α = 1.81.

3.1.2 Non Smooth Initial Data - The Level set Function

In this section the level set equation (4) is advected onto the rotating disk. The non-smooth discontinuous initial data from Figure 4 is described by the following formula

φ0(x, 0) =

(1 if (x1− x^{0}_{1})^{2}+ (x2+ x^{0}_{2})^{2}≤ r_{0}^{2},

0 otherwise, (39)

where (x^{0}1, x^{0}_{2}) = (0.3, 0) and r0= 0.25. Figure 5 presents the solution at time T = 1 obtained
using the RV method. The solution is plotted for the finest mesh size tested hmax = 0.033.

(a) (b)

Figure 4: ’Non-smooth’ initial data for the rotating disk problem using mesh size hmax= 0.033.

(a) (b)

Figure 5: Solution to the rotating disk problem after one complete rotation (T = 1) using the ’non-smooth’ initial data, Residual Viscosity stabilization with mesh size, hmax= 0.033.

The convergence rate

In Table 2 the error in the L^{2}-norm of the approximation is examined for varying parameter
hmax. Figure 6 presents a loglog-plot of hmaxversus the L^{2}-norm of the errors. The slope of the
line determines the convergence rate α = 0.46. The function h^{α}max is plotted for comparison.

Table 2: L^{2}-error and convergence results for the rotating disk problem using

’non-smooth’ initial data, varying mesh size and RV stabilization.

hmax L^{2}-error convergence (q)

0.033 0.1093 0.46

0.065 0.1518 0.43

0.137 0.2092 -

Figure 6: A loglog -plot of hmax versus the L^{2}-norm of the errors for the rotating
disk problem, using the ’non-smooth’ initial data. The red line is the measured
error and the blue line is h^{α} using α = 0.46.

3.1.3 Artificial Viscocity over Time

The artificial viscosity h was plot over time. Figure 7 presents a plot of the largest value of artificial viscosity found in each time step. On Figure 7b the artificial Viscocity for the smooth data from Figure 7a is plot.

(a) Non Smooth and Smooth Data (b) Smooth Data

Figure 7: Artificial Viscocity over time.

For the smooth solution, the residual stays small though out the simulation, which makes the artificial viscosity take small values. For the non-smooth solution, the residual is large at the beginning, hence larger values of h. The RV method smooths out the regions in which the residual is large, thus n decreases with time.

### 3.2 A residual-based artificial viscocity compression method - Anti- diffusion

In this section the motivated compression technique [1], (see e.g, Section 2.6) is implemented for the level set equation (4). Preliminary experiments showed that good solutions were obtained

using the RV stabilization parameters set to Cvel= 0.5 and CRV = 100and these values were used in the following analysis. Initial data are set as in Figure 4 and Eq. (39). Figure 11 presents the solution at time T = 1 obtained using the Antidiffusion method. The solution is plot for the finest mesh size tested hmax= 0.033.

(a) (b)

Figure 8: Solution to the rotating disk problem after one complete rotation (T

= 1) using the ’non-smooth’ initial data, RV stabilization and Antidiffusionn with mesh size hmax= 0.033.

Figure 9: Solution to the rotating disk problem after one complete rotation (T

= 1) using the ’non-smooth’ initial data, RV stabilization and Antidiffusion with mesh size hmax= 0.033. Here using a solid colormap.

The convergence rate

In Table 3 we plot the error in the L^{2}-norm of the approximation for varying parameter hmax.
Figure 10 presents a loglog-plot of hmax versus the L^{2}-norm of the errors. The slope of the line
determines the convergence rate α = 0.58. The function h^{α}maxis plotted for comparison.

### 3.3 Normal Vector

In this section the normal vector of the solution from the level set equation (4) is computed.

Figure 11 shows the normal vector field obtained at T = 0 and T = 1 using the Antidiffusion method during the advection using for cell= 1(see e.g, Sections 2.6 and 3.2).

Table 3: L^{2}-error and convergence results for the rotating disk problem using

’non-smooth’ initial data, varying mesh size, RV stabilization and Antidiffusion.

h_{max} L^{2}-error convergence (q)

0.033 0.0844 0.58

0.065 0.1162 0.68

0.137 0.1913 -

Figure 10: A loglog -plot of hmax versus the L^{2}-norm of the solution errors for
the rotating disk problem using the Antidiffusion stabilization. The red line is the
measured error and the blue line is h^{α}using α = 0.58.

(a) T = 0 (b) T = 1

Figure 11: Normal vector field for the initial state of the rotating disk and the
state after one complete rotation where Antidiffusion is used. Mesh size hmax =
0.033, elliptic smoothing coefficient cell= 1 and normalization term = 10^{−4}.

### 3.4 Curvature

In this section we present the curvature of the interface for T =1, and for cell = 1 (see, e.g Section 2.9). The interface forms a circle of radius r = 0.25 (39), as such the analytic curvature

is constant and is equal to κ = 4. The curvature of the advected interface at T = 1 is shown in
Figures 12 and 13, for three mesh sizes, on the same figure the analytic curvature is also plot for
comparison. The convergence rate of the curvature and the error in the L^{2}-norm is computed
for several mesh sizes, Table 4 using cell = 1and Table 5 using cell = 20. Figure 14 compares
the curvature contour lines for the mesh sizes used in Table 5 with the analytical contour line.

The contour lines do not appear to converge towards the exact analytical line. However, the contour lines for larger cell are more circular than for lower values.

Figure 12: Level 4 curvature contour lines from the rotating disk problem. An- alytical curvature line in white and numerical curvature contour line in grey after one complete rotation. Mesh size hmax= 0.137 and elliptic smoothing cell= 1.

(a) Mesh size hmax= 0.065 (b) Mesh size hmax= 0.033

Figure 13: Contour lines for curvature equal to 4 from the rotating disk problem.

Analytical curvature line in white and numerical curvature contour line in beige after one complete rotation (T=1). Two mesh sizes and elliptic smoothing cell= 1.

Figure 14: Contour lines for curvature equal to 4 from the rotating disk problem after one complete revolution (T=1). The mesh sizes correspond to the sizes used in table 5 with hmax= 0.023 in dark blue, hmax= 0.038 in cyan, hmax= 0.046 in light blue and the analytical contour line in white. The curvature values for the coarsest mesh, hmax= 0.091, are all smaller than 4 and the corresponding contour lines do not exist.

Table 4: L^{2}-error and curvature convergence results for the rotating disk problem
using Antidiffusion with cell = 1 in the elliptic smoothing of the normal and
curvature.

hmax L^{2}-error convergence (q)
0.033 0.6127 negative
0.065 0.4639 near zero

0.137 0.4660 -

Table 5: L^{2}-error and curvature convergence results the rotating disk problem
using Antidiffusion with cell = 20 in the elliptic smoothing of the normal and
curvature.

hmax L^{2}-error convergence (q)

0.023 0.45 1.2

0.038 0.0688 3.97

0.046 0.43 2.46

0.091 2.31 -

### 3.5 Incompressible Navier stokes

In this section we present numerical solutions to the two-phase flow system defined in equation
(13). Where the advection equation (4) is solved according to Section 2.7. The incompressible
Navier–Stokes equation (10) is solved with the IPCS given in Section 2.8. Body forces (g) and
the surface tension (Tη) is excluded from the system. The computational domain is a channel
placed on c1 = (0, 0), c2 = (4, 0), c3 = (0, 1), c4 = (4, 1). The boundaries of the fluid-flow are
given by: u|x_{1}=0 = u0(x2), u|x_{2}=0= u|x_{2}=1= 0 and p|x_{1}=4= 0. The density is ρ1= ρ2= 1, if

other is not given.

3.5.1 Channel Flow - Constant Inflow Profile

The inflow profile is set to u0(y) = (1, 0). The interface (φ(x, t) = ^{1}_{2}) between fluid 1 and fluid
2 is represented by the white curve in Figure 15. The magnitude of the velocity field from the
incompressible Navier-Stokes is represented as a colormap. The level-set function is initially set
to be a cylinder with radius r = ^{1}_{4} positioned at c5= (1,^{1}_{2}).

(a) T = 0 (b) T = 1

Figure 15: Magnitude of the fluid velocity and the interface of the initial setting and the final solution. Mesh size is hmax = 0.012 and the viscosity of both fluids are µ = 0.1.

3.5.2 Channel Flow - Variable Density - Constant Inflow profile

We consider a two-phase flow with identical setting as previous Section (3.5.1), but with a variable density. The density of fluid 1 is ρ1 = 1and of fluid 2 it is ρ2= 5, implemented as in equation (15). Streamlines are included, see Figure 16, in order to denote the more complicated flow profile. The system is shown at four distinct time steps of interest.

3.5.3 Channel Flow - non Constant Inflow profile

The inflow profile of the velocity field in (10) is set as inflow: u0(y) = (sin(y), 0). We plot the advected surface at T = 1 , Figure 19b. We plot the according velocity field both at T = 0, Figure 18a and at T = 1, Figures 18b and 19a. The unit cylinder is placed on c5 = (0.5, 0.5) Figure 17.

3.5.4 Channel Flow - Obstacle - non Constant Inflow profile

The inflow profile of the velocity field in (10) is set as inflow: u0(y) = (sin(y), 0). We plot the advected surface at T = 1 , Figure 21b. We plot the according velocity field both at T = 0, Figure 20a and at T = 1, Figures 20b and 21a.

We plot the contour line of the interface for several meshes to investigate whether the resulted
plot converges to a certain shape, Figure 22. The dark blue line corresponds to the finer mesh,
hmax = 0.033, the light blue line corresponds to hmax = 0.065and the white line corresponds
to hmax = 0.137. In order to measure mass conservation of the fluid inside the channel, the
total area of the advected surface at T = 1, denoted as A1, is computed for three different
meshes. The area A1 is then compared with the analytic area of the initial shape, denoted as
Aanalytic, Figure 17b. The analytic area of the initial shape is constant Aanalytic= π ∗ r^{2}, where
r = 0.25is the radius of the circle. The number e = ||Aanalytic− A1||_{L}^{2} is computed for a

(a) T = 0 (b) T = 0.03

(c) T = 0.25 (d) T = 1

Figure 16: Magnitude of the fluid velocity and the interface of four different T.

Black lines are streamlines of the velocity field. Mesh size is hmax= 0.012 and the viscosity of both fluids are µ = 0.1.

(a) Channel (b) Channel with obstacle

Figure 17: Initial condition for the NS setting involving flow inside a channel with and without an obstacle, the mesh size is , hmax= 0.033, and the viscocity of the fluid is µ = 0.1.

fine, a medium and a coarse mesh. As the mesh gets finer the area at T = 1, A1approaches the analytic area Aanalytic , while e approaches zero, Table 6.

3.5.5 Flow-past a cylinder

In this section we use a velocity field we obtained from Section 3.4 of the FEniCS project documentation [11], Figure 24. The initial condition is on Figure 23. The radius of the cylinder is set r = 0.2 and is placed at c0= (0.7, 0.2). The example was tested for two different viscocity values µ = 0.001, Figure 26 and µ = 0.5, Figure 25.

The final example, flow-past a cylinder, produces an interesting velocity field, two fluids with

(a) T = 0 (b) T = 1

Figure 18: The velocity field is plotted at T = 0 and T = 1, it is obvious , from Figure 18b that the higher velocity values are obtained in the middle of the channel at T = 1.

(a) velocity field (b) Solution

Figure 19: The advected surface is plotted at T = 1, for a mesh size of hmax= 0.033, and fluid viscocity µ = 0.1. In Figure 19a the corresponding vector velocity field is also plotted at T = 1.

Table 6: Analytic area of the initial cylinder Aanalytic and area of the advected
shape A1, the number e = ||Aanalytic− A1||_{L}2 is computed for varying mesh
size.

hmax Aanalytic A1 e = ||Aanalytic− A1||L^{2}

0.033 0.1963 0.1947 0.0015

0.065 0.1963 0.1925 0.0037

0.137 0.1963 0.1593 0.036

different viscocity values were chosen to test the behaviour of the advected cylinder. In Figure 25 the fluid with the higher viscocity value µ = 0.5 produces a promising result, as the final plot is smooth and the area of the shape is preserved. Whereas on Figure 26 when a fluid with very low viscocity is tested the resulted plot is broken apart into two pieces and some oscillations are appearing, mainly the on the white-faded tail on Figure 26b.

(a) T = 0 (b) T = 1

Figure 20: The Velocity Field of the channel flow involving an obstacle is plot at T = 0 and T = 1, it is shown on figure 20b that the higher velocity values are obtained in the middle of the channel at T = 1.

(a) Velocity Field (b) Solution

Figure 21: The advected surface is plotted at T = 1 for a mesh size of hmax = 0.033 and a fluid of viscocity µ = 0.1. The vector velocity field is plot on figure 21a.

### 4 Discussion

In this project we study a two-phase flow problem, more precisely a particular level-set method:

the conservative level-set approach in [1]. The interface between the immiscible fluids is obtained
by the finite element method using a particular compression term (Antidiffusion). The level set
equation (4) is advected and the L^{2}-norm of the error is measured after a full rotation around
the unit disk. The solution obtained by the Antidiffusion method is sharper around the edges
compared with the solution obtained by the RV method. Moreover the convergence rate at T

= 1 using the Antidiffusion method is greater compared with the convergence rate obtained by the RV method.

The computation of the normal vector and the curvature of the interface after one full rotation are of main interest. In order to compute these functionals, a particular form of the finite element projection is used, including an elliptic term. The elliptic term is governed by the cell

Figure 22: Contour lines of the advected interface for varying mesh sizes, T = 1, fine mesh (dark blue), medium mesh (light blue), coarse mesh(white)

Figure 23: The initial condition on the setting flow-past a cylinder, for a mesh size of hmax= 0.033

(a) T = 0 (b) T = 1

Figure 24: The Velocity Field at T = 0 and T = 1 for the setting flow-past a cylinder, the higher velocity values at T = 1 are obtained in the middle of the channel.

coefficient and the quadratic h^{2}, where h is the mesh size. The curvature error is measured
comparing a perfect circle with radius r = 0.25 against a contour line of value 4, resulted from
the computational curvature. The curvature error , measured in the L^{2}-norm increases as we

(a) Velocity field (b) solution

Figure 25: The advected surface at T = 1, for a mesh size hmax= 0.033, and a viscocity value of µ = 0.5. The according vector velocity field is also plotted.

(a) Velocity field (b) solution

Figure 26: the advected surface at T = 1, for a mesh size hmax = 0.033, and a viscocity value of µ = 0.001. The according vector velocity field is also plotted.

refine the mesh for values of elliptic smoothing cell< 14, Table 4 showed example results using
c_{ell} = 1, and no definite convergence value could be obtained. For large values of smoothing,
the convergence values varied greatly and after a certain point in the mesh refinement, the error
increased. Table 5 shows an example result using cell = 20.No definite convergence values could
be obtained. It is, however, possible that further mesh refinement would yield more consistent
convergence results, or that a different method of evaluating the curvature would yield conclusive
convergence results in the range of mesh sizes considered above.

The work is extended in order to include computations involving the incompressible Navier- Stokes (INS) equations. In order to implement the advection of the interface inside a velocity field computed from the solution of the INS, three different settings were used involving a channel flow. Moreover inside the INS setting two different velocity profiles were used involving a constant and a non constant inflow. In the first test the interface deforms as expected, considering the flow profile (Figure 15). In the variable density case the less dense fluid is forced around the denser until they are both accelerated to approximately the same stable velocity.

Variable density/viscosity systems would arguably be a more complicated to model without the sharp level-set function, produced by the antidiffusion. More complex geometries and inflow profiles also show satisfactory dynamical behavior of the interface. To test the convergence of

the INS solution the contour line of the advected interface was plot for three different meshes.

The resulted plots converged to a similar shape. The area of each fluid was preserved during the flow and the preservation was shown to be convergent.

The final INS experiment involved flow-past a cylinder, this setting was tested for two different viscocity values, µ = 0.5 and µ = 0.001. For the higher viscocity value µ = 0.5 the advected shape maintained its total area and the plot at T = 1 is smooth. For the lower viscocity value µ = 0.001 the plot at T = 1 is not smooth and the final shape is broken into two pieces. In order to resolve the turbulence-issue that appears when using a fluid with low viscocity, the INS setting should be enhanced with computations that involve the surface tension. However, including effects of surface tension requires an accurate knowledge about the interface curvature and the results on curvature convergence were inconclusive. Thus, the Antidiffusion technique together with the method for evaluating the curvature used in this report can not yet be used for the purpose of including surface tension effects in fluid simulations.

### 5 Conclusions

In this project we study a two-phase flow problem. The interface between the two immiscible fluids is obtained by the finite element method using a conservative level-set approach. The main interest is the computation of the normal vector and the curvature of the advected interface.

The following main conclusions are made.

• A residual-based artificial viscosity (RV) h acts as a good stabilization technique for
the discontinuous level-set equation. The convergence rate using RV, measured in the
L^{2}-norm, is 0.46

• The Antidiffusion technique produces a sharper solution compared with the RV method.

Also the convergence rate using Antidiffusion, measured in the L^{2}-norm, is 0.58 and thus
greater than the rate obtained by the RV method.

• Elliptic projection of the normal and curvature, on to the finite element space, yields a better estimate.

• The results on convergence rate for the interface curvature error are inconclusive. More precisely, the convergence varies as the mesh is refined.

• Two-phase flow simulations convey a positionally convergent behavior of the interface and the area is well preserved.

• In order to model more complicated flows, such as turbulent flows, surface-tension needs to be included. This would demand better approximations of the normal and the curvature.

As a final remark the RV and Antidiffusion technique for the level-set method appear to yield good results in tracking of the interface. However, the obtained results on the curvature con- vergence are inconclusive and further testing or developing different methods for evaluating the curvature and curvature error is recommended for future work on the topic.

### References

[1] Jean-Luc Guermond, Manuel Quezada de Luna, and Travis Thompson. An conservative anti-diffusion technique for the level set method. Journal of Computational and Applied Mathematics, 321:448 – 468, 2017.

[2] Markku Kataja, Antti Koponen, V.-M Luukkainen, M. Manninen, Ulla Ojaniemi, J. Po- ranen, Juha Salmela, Pasi Selenius, H. Eloranta, Jouko Halttunen, Markus Honkanen, T. Pärssinen, Pentti Saarenrinne, S. Zhou, M. Engblom, K. Fagerudd, M.F. Geant, A. Her- manson, Sirpa Kallio, and E. Rehn. Multiphase flows in process industry: Promoni. VTT Tiedotteita - Valtion Teknillinen Tutkimuskeskus, pages 3–177, 01 2005.

[3] J. Yao and M. Takei. Application of process tomography to multiphase flow measurement in industrial and biomedical fields: A review. IEEE Sensors Journal, 17(24):8196–8205, 2017.

[4] Christopher E. Brennen. Fundamentals of Multiphase Flow. Cambridge University Press, 2005.

[5] Nikolaos D. Katopodes. Chapter 13 - level set method. In Nikolaos D. Katopodes, editor, Free-Surface Flow, pages 804 – 828. Butterworth-Heinemann, 2019.

[6] Ronald P. Fedkiw, Guillermo Sapiro, and Chi-Wang Shu. Shock capturing, level sets, and pde based methods in computer vision and image processing: a review of osher’s contributions. Journal of Computational Physics, 185(2):309 – 341, 2003.

[7] Elin Olsson and Gunilla Kreiss. A conservative level set method for two phase flow. Journal of Computational Physics, 210(1):225 – 246, 2005.

[8] Stéphane Popinet. Numerical Models of Surface Tension. Annual Review of Fluid Mechan- ics, 50:49 – 75, January 2018.

[9] Murtazo Nazarov and Johan Hoffman. Residual-based artificial viscosity for simulation of turbulent compressible flow using adaptive finite element methods. International Journal for Numerical Methods in Fluids, 71:339–357, 01 2013.

[10] K. Selim, A. Logg, and M.G. Larson. An adaptive finite element splitting method for the incompressible navier–stokes equations. Computer Methods in Applied Mechanics and Engineering, 209-212:54 – 65, 2012.

[11] Martin S. Alnæs, Jan Blechta, Johan Hake, August Johansson, Benjamin Kehlet, Anders Logg, Chris Richardson, Johannes Ring, Marie E. Rognes, and Garth N. Wells. The fenics project version 1.5. Archive of Numerical Software, 3(100), 2015.