PROJECTREPORT ComputationalTechniquesonTwoPhaseFlow

(1)

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

(2)

In this work, we study a conservative level-set method [1] for two-phase flow. In particular, it is investigated how this level-set approach performs in combination with a residual- based second order artificial viscosity stabilization technique for finite element approximations. It is proved that the convergence rate of the finite element approximation is close to 0.5, measured in the L²-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 inconclusive. 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].

(4)

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 computation 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 computation 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 = − ∇φ

||∇φ||_l2

, (2)

where the l²-norm denotes the local L²-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}, Ω₂= {x|φ(x) > 0.5}, I = {x|φ(x) = 0.5}.

(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 conservation laws, of the form (4). Consider a one-dimensional, possibly stabilized, conservation equation

∂_tφ + ∂_xf (φ) = ∂_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 compression 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· `⁻¹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²

⁺

∇φ

!

= 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)

(6)

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

(u) =1

2(∇u + ∇^Tu). (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,







ρ(∂_tu + 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², 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) = φ₀(x), x ∈ Ω,

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

(16)

Here φ represents the unknown variable, f ∈ C¹(R, R²) is the flux term and φ0 ∈ L^∞(R²) 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₀, defined by V0= {v : kvk_L2+ k∇vk_L2< ∞, v|_∂Ω= 0}. The weak form reads:

(7)

Find φ(x, t) ∈ V0such that:

(∂φ

∂t, v) + (∇ · f (φ), v) = 0, ∀v ∈ V0, (17) where the norm (·, ·) is the L²-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⁰(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²1+x²₂≤ 1}, and rewrite the flux term in the non-conservative form, i.e., ∇·f(φ) := f⁰(φ) · ∇φ. Assume that f⁰(φ) := 2π(−x₂, x₁)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⁰(φ)the Finite element Approximation reads:

(∂φh

∂t , v) + (f⁰(φ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 approximation, 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 viscosity 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.

(∂φh

∂t , v) + (f⁰(φ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⁰(φn) · ∇φn,

(8)

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

n,K = min CvelhKβK, CRVh²_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⁰(φn)k_L∞(K)= max

Ni∈K

[f₁⁰(φn)²+ f₂⁰(φn)²]¹²(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⁰(φh) · ∇φh, v) + (hΘ(φh)∇φh, ∇v) = 0 ∀v ∈ Vh,0, Θ(φh) =

1 − ccomp

[(φh− φL)(φR− φh)]⁺ h||∇φh||_L2

+

.

(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⁰· ∇φh, v) − (hΘ∇φh, ∇v) ∀v ∈ Vh,0, (24) where the advection velocity f⁰(φ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)

(9)

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 φⁿ_hfrom φⁿ⁻¹_h , the RK4 weights are obtained as

(w0, v) = rn(φⁿ⁻¹_h , v), (26)

(w₁, v) = r_n(φⁿ⁻¹_h + 0.5k_nw₀, v), (27) (w2, v) = rn(φⁿ⁻¹_h + 0.5knw1, v), (28) (w₃, v) = r_n(φⁿ⁻¹_h + k_nw₂, v), (29) again, using the same f⁰(φh), Θ(φh), h(tn), obtained at time step tn for each RK4 stage. The final solution is obtained using the weights as

φⁿ_h= φⁿ⁻¹_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ⁿ

k + ρuⁿ· ∇uⁿ = ∇ · σ(uⁿ⁺¹², pⁿ) + fⁿ⁺¹, (31) where k is the step in time.

The next step is to update the pressure as

∆pⁿ⁺¹= ∆pⁿ−∇ · u^∗

k . (32)

The updated velocity is finally calculated by uⁿ⁺¹− u^∗

k + ∇pⁿ⁺¹− ∇pⁿ= 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.

(10)

2.9.1 Elliptic Projection

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

(ue, v) + cellh²(∇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||_L2+ δ, (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²-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₀ 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²-error of the curvature is the line integral

e_κ= s

Z

I

(κ_analytic− κ₁)²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

(11)

e_κ= s

Z 2π 0

r₀(κ_analytic(θ) − κ₁(θ))²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₁− x⁰₁)²+ (x₂− x⁰₂)²

r²₀ − 1

, (38)

where (x⁰1, x⁰₂) = (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²-norm of the approximation was examined for varying parameter hmax. Figure 3 presents a loglog-plot of hmax versus the L²-norm of the errors. The slope of the line determines the convergence rate α = 1.81. The function h^αmax was plotted for comparison.

(12)

(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²-error and convergence results for the rotating disk problem, using

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

hmax L²-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²-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⁰₁)²+ (x2+ x⁰₂)²≤ r₀²,

0 otherwise, (39)

(13)

where (x⁰1, x⁰₂) = (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.

In Table 2 the error in the L²-norm of the approximation is examined for varying parameter hmax. Figure 6 presents a loglog-plot of hmaxversus the L²-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²-error and convergence results for the rotating disk problem using

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

0.033 0.1093 0.46

0.065 0.1518 0.43

0.137 0.2092 -

(14)

Figure 6: A loglog -plot of hmax versus the L²-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

(15)

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.

In Table 3 we plot the error in the L²-norm of the approximation for varying parameter hmax. Figure 10 presents a loglog-plot of hmax versus the L²-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).

(16)

Table 3: L²-error and convergence results for the rotating disk problem using

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

h_max L²-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²-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⁻⁴.

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

(17)

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²-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.

(18)

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²-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²-error convergence (q) 0.033 0.6127 negative 0.065 0.4639 near zero

0.137 0.4660 -

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

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₁=0 = u0(x2), u|x₂=0= u|x₂=1= 0 and p|x₁=4= 0. The density is ρ1= ρ2= 1, if

(19)

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) = ¹₂) 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 = ¹₄ positioned at c5= (1,¹₂).

(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², where r = 0.25is the radius of the circle. The number e = ||Aanalytic− A1||_L² is computed for a

(20)

(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

(21)

(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||_L2 is computed for varying mesh size.

hmax Aanalytic A1 e = ||Aanalytic− A1||L²

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.

(22)

(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²-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

(23)

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², 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²-norm increases as we

(24)

(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

(25)

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²-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²-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 convergence 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.

(26)

[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.