Numerical Algorithms for Divergence-Free Velocity Applications

Preprint

Article

Numerical Algorithms for Divergence-Free Velocity Applications

Altmetrics

Downloads

107

Views

Comments

A peer-reviewed article of this preprint also exists.

This version is not peer-reviewed

Submitted:

15 July 2024

Posted:

16 July 2024

You are already at the latest version

Alerts

Abstract

This work focuses on the well-known issue of mass conservation in the context of the finite element technique for computational fluid dynamic simulations. Specifically, non-conventional finite element families for solving Navier-Stokes equations are investigated to address the mathematical constraint of incompressible flows. Raviart-Thomas finite elements are employed for the achievement of a discrete free-divergence velocity. In particular, the proposed algorithm projects the velocity field into the discrete free-divergence space by using the lowest-order Raviart-Thomas element. This decomposition is applied in the context of the projection method, a numerical algorithm employed for solving Navier-Stokes equations. Numerical examples validate the approach’s effectiveness, considering different types of computational grids. Additionally, the presented paper considers an interface advection problem using marker approximation, in the context of multiphase flow simulations. Numerical tests, equipped with an analytical velocity field for the surface advection, are presented to compare exact and non-exact divergence-free velocity interpolation.

Keywords:

Subject: Computer Science and Mathematics - Applied Mathematics

1. Introduction

In the context of incompressible flow simulations, Computational Fluid Dynamics codes play a pivotal role. In particular, divergence-free fields are crucial for mass conservation in numerical simulations of engineering applications (multiphase flows, porous-media flows, etc.). Typically, commercial codes employ finite element (FEM) and finite volume (FVM) methods to approximate solutions of the physics problem, both widely applicable across various engineering domains. While FEM and FVM methods share popularity and similar computational costs, the finite volume method is usually preferred for fluid dynamic simulations due to concerns about mass conservation related to the finite element method. Nevertheless, achieving the desired divergence-free field across the discrete domain remains challenging for both methods.

This study addresses the challenge of obtaining a finite element approximation of the solution for the Navier-Stokes system, specifically aiming for a divergence-free velocity field over the discrete domain. This is achieved through the utilization of Raviart-Thomas basis functions [1].

When dealing with fluid dynamics simulations, the mathematical analysis of the specific partial differential equations is fundamental for obtaining reliable numerical results [2]. Focusing on the incompressible Navier-Stokes equations, the literature of the past fifty years has increasingly emphasized the use of mixed finite elements [3,4,5,6], where different types of numerical discretizations are employed to represent different variables. Indeed, this discretization approach has proven to yield convergent numerical schemes, theoretical convergence rates, and other advantageous properties.

Naturally, in the framework of finite element discretization, the flexibility of mixed methods arises from the relaxation of the divergence constraint [7]. However, this relaxation entails a cost, which becomes evident when considering a standard error estimate for the considered equation. To provide a brief example, when dealing with error estimates for the velocity field in the Navier-Stokes system, non-divergence mixed finite elements present a connection between the discretized velocity and the continuous pressure [8]. This mathematical connection became a numerical drawback for certain finite element families that can lead to a loss of order of convergence.

In describing divergence-free discretization, various approaches have been proposed in the literature, such as the Scott-Vogelius element and the discontinuous Galerkin method [9,10,11,12,13,14,15]. Despite these efforts, the mathematical problem, commonly referred to as poor mass conservation [16], remains a subject of ongoing interest. Indeed, some stabilization techniques have been proposed to overcome this issue, such as the grad-div stabilization [17] or a transformation on the continuous pressure [18].

Moreover, some remarks can be pointed out regarding the equations described by the incompressible Navier-Stokes system. The existence of the solution depends on the divergence operator and specific surjectivity properties are required, such as the well-known inf-sup condition, which is a sufficient condition to find a unique solution for a saddle point problem. In addition, the problem must also preserve the invariance property, ensuring that a change in the external body force by adding a gradient field alters only the solution of the pressure and not the velocity.

As pointed out in [19], a lack of

L^{2} -

orthogonality between discretely-divergence-free vector and irrotational fields, such as the pressure gradient, may generate a poor momentum balance, which translates into poor mass conservation.

In this work, we aim to exploit a finite element family designed to address the numerical challenges and maintain a divergence-free velocity field. Expressly, within the framework of mixed problems, the Raviart-Thomas finite element family assumes an important role, suitable for resolving partial differential equations subject to divergence constraints.

This problem is examined in the context of coupled and split pressure-velocity formulations of the momentum equation. Indeed, the solution of the coupled incompressible Navier-Stokes system is known for its high computational effort, leading to the development of various numerical algorithms for treating the velocity and the pressure separately [20] and therefore decreasing significantly the computational burden. We recall that the coupling between velocity and pressure is due to the incompressibility constraint, resulting in saddle-point matrices in discrete form.

To address these challenges, Chorin and Temam first introduced the projection method in the framework of finite element methods, splitting the Navier-Stokes system into two distinct steps: one for resolving the velocity field and another for the pressure field [21,22]. This approach has proven to reduce the computational effort and encouraged the development of various projection methods [23].

The requirement for a divergence-free velocity gains significance in specific fluid dynamic simulations, especially in scenarios involving multiphase flows where mass conservation is crucial for reliable numerical outcomes. This mathematical constraint is addressed by leveraging a divergence-free representation of the velocity field in multiphase flows. Specifically, the multiphase problem considered in this work involves the surface advection of a single phase using marker technique approximation.

The paper is structured as follows: in the first section, we outline the mathematical framework with a description of the problem related to the divergence-free constraint. After that, we provide an overview of finite element discretization, with specific attention to the Raviart-Thomas finite elements within quadrilateral and hexahedral elements. Subsequently, a focus on a multiphase flow advection problem is reported considering the velocity orthogonal decomposition by using Raviart-Thomas basis functions. Finally, the algorithm proposed for the resolution of the Navier-Stokes system by projection technique is presented considering divergence-free finite elements for the projection step.

2. Mathematical Framework

In the following, the main basic properties regarding the numerical solution of the Navier-Stokes equations will be introduced. In particular, the mathematical framework will be described, with a major effort on the divergence-free constraint.

2.1. Notation

To present and describe the Navier-Stokes equation mathematically, the typical weak formulation requires foundational concepts from variational calculus, such as appropriate Sobolev spaces [24]. Let

Ω

be an open bounded subset of

R^{d}

where the boundary

Γ = \partial Ω

is locally represented by a Lipschitz function. This contextual description leads us to introduce

L^{2} (Ω) : = \{f : \int_{Ω} | f^{2} | d Ω < \infty\},

(1)

the space comprising square-integrable functions with inner product and norm as

(f, g) : = \int_{Ω} f g d Ω

and

{| | f | |}_{0}^{2} : = (f, f)

, respectively. Also, we denote

L_{0}^{2} (Ω)

the space of square-integrable functions with vanishing mean. The Sobolev spaces for integers

k \geq 0

, are defined as follows

\begin{matrix} H^{k} (Ω) & : = \{f \in L^{2} (Ω) : \partial^{α} f \in L^{2} (Ω), for | α | \leq k\}, \end{matrix}

(2)

where the non-negative integer indices

α = (α_{1}, α_{2}, \dots, α_{d})

denote the order of the partial derivatives, i.e.

| α | : = \sum_{k = 1}^{d} α_{k}

, and

\partial^{α} f = \partial_{x_{1}}^{α_{1}} \partial_{x_{2}}^{α_{2}} \dots \partial_{x_{d}}^{α_{d}} f

. The respective norms and semi-norms are defined as

\begin{matrix} {| | f | |}_{k} & : = {(\sum_{| α | \leq k} | | \partial^{α} f {| |}_{0}^{2})}^{\frac{1}{2}}, & {| f |}_{k} & : = {(\sum_{| α | = k} | | \partial^{α} f {| |}_{0}^{2})}^{\frac{1}{2}} . \end{matrix}

(3)

Additionally, by introducing the fractional-order Sobolev space on the boundary

Γ

, we get

H^{\frac{1}{2}} (Γ) : = \{f \in L^{2} (Γ) : {| f |}_{\frac{1}{2}, Γ} < \infty\},

(4)

with the corresponding norm and semi-norm given by

\begin{matrix} {| f |}_{\frac{1}{2}} & : = {(\int_{Γ} \int_{Γ} \frac{{| f (s) - f (t) |}^{2}}{{| s - t |}^{2}} d s d t)}^{\frac{1}{2}}, & {| | f | |}_{\frac{1}{2}} & : = {({| | f | |}_{0, Γ}^{2} + {| f |}_{\frac{1}{2}, Γ}^{2})}^{\frac{1}{2}} . \end{matrix}

(5)

For vector-valued functions such as

u = (u_{1}, u_{2}, \dots, u_{d})

, we can consider the following norms on

{(H^{k} (Ω))}^{d}

\begin{matrix} {| | u | |}_{k} & : = (| | u_{1} {| |}_{k}^{2} + | | u_{2} {| |}_{k}^{2} + \dots + | | u_{d} {| |}_{k}^{2})^{\frac{1}{2}}, \end{matrix}

(6)

while for

p \in H^{1} (Ω)

we can define on

Γ

a surjective trace map with continuous lifting

γ : H^{1} (Ω) \to H^{\frac{1}{2}} (Γ)

such that

{p |}_{Γ} = γ p

for all

p \in H^{1} (Ω)

. The space

H^{\frac{1}{2}} (Γ)

is equipped with the norm

{∥ p ∥}_{\frac{1}{2}, Γ} = {inf {∥ w ∥}_{1, Ω} | w \in H^{1} (Ω), γ w = p}

. We denote by

H^{- \frac{1}{2}} (Γ)

the dual space of

H^{\frac{1}{2}} (Γ)

with the norm

∥ p^{*} ∥_{- \frac{1}{2}, Γ} = sup {〈 p, p^{*} 〉 / {∥ p ∥}_{\frac{1}{2}, Γ} | p \in H^{\frac{1}{2}} (Γ)}

, where the bracket

〈 \cdot, \cdot 〉

denotes duality between the two spaces.

In addition, for a vector field

v \in L^{p} (Ω)

(p \geq 1)

the functional

〈 \nabla \cdot v, ψ 〉 = - \int_{Ω} \nabla ψ \cdot v d x \forall ψ \in C_{0}^{\infty} (Ω)

(7)

defines the distributional divergence in

Ω

. Moreover, if there exists a function

ρ \in L_{loc}^{1}

such that

〈 \nabla \cdot v, ψ 〉 = \int_{Ω} ψ ρ d x \forall ψ \in C_{0}^{\infty} (Ω)

(8)

then

ρ

is called weak divergence in

L_{loc}^{1}

, which is the space of locally integrable function. With the term ’locally’ we mean that the function is integrable on every compact subset of its definition domain. If

ρ = 0

, the vector is called a divergence-free vector, namely a divergence-free vector on

Ω

is a vector that is orthogonal to all gradient fields with compact support on

Ω

With the previously defined functional spaces, we can now introduce the space used for the mixed formulation of second-order elliptic problems. Hence, considering a vector-valued function, we define the space

H (div, Ω)

H (div, Ω) = \{v \in L^{2} (Ω) : \nabla \cdot v \in L^{2} (Ω)\},

(9)

that is a Hilbert space equipped with the norm

{| | v | |}_{H} (div) : = {({| | v | |}_{0}^{2} + | | \nabla \cdot {v | |}_{0}^{2})}^{\frac{1}{2}} .

(10)

It is important to note that

H^{1} (Ω)

is continuously embedded in

H (div, Ω)

. On this space it is then possible to define a surjective map called normal trace

γ_{n} : H (div, Ω) \to H^{- 1 / 2} (Γ)

such that

γ_{n} v = v \cdot n_{Γ}

for all

v \in H (div, Ω)

and

n_{Γ}

the normal unit vector to

Γ

. We recall the divergence theorem

\begin{matrix} \int_{Ω} \nabla \cdot v p d x = \int_{Γ} v \cdot n p d x - \int_{Ω} v \cdot \nabla p d x \end{matrix}

(11)

for all

v \in H (div, Ω)

and

p \in H^{1} (Ω)

. Also, we remark that

H^{1 / 2} (Γ)

does not contain characteristic functions that are zero on a proper subset

Γ_{0} \subset Γ

or, in other words, a function in

H^{1 / 2} (Γ_{0})

cannot be extended by zero outside

Γ_{0}

to a function in

H^{1 / 2} (Γ)

. This should be taken into account in (11) when the duality is used to characterize the trace function

γ_{n} v \in H^{- 1 / 2} (Γ)

. For details on Sobolev and divergence spaces, the reader can refer to [24,25].

2.2. The Mixed-Method-Robustness Property

We briefly recall now the main mathematical issues related to the numerical solution of the Navier-Stokes system. Let

Ω

be a subset of

R^{d}

, with

d \in {2, 3}

. We illustrate the mixed-method-robustness property in the case for the Stokes for simplicity but this can be applied to all mixed-method systems.

The Stokes equations can be stated in strong form as

\begin{matrix} - ν Δ u + \nabla p & = f, \\ - \nabla \cdot u & = g, \\ {u |}_{Γ} & = 0 . \end{matrix}

(12)

Naturally, the incompressible formulation holds when

g = 0

. The variational formulation of (12), for the state

(u, p) \in H_{0}^{1} (Ω) \times L_{0}^{2} (Ω)

, can be written as

\begin{matrix} \int_{Ω} ν \nabla u^{T} : \nabla w d x + \int_{Ω} \nabla p \cdot w d x & = \int_{Ω} f \cdot w d x, \forall w \in H_{0}^{1} (Ω) \\ - \int_{Ω} ψ \nabla \cdot u d x & = \int_{Ω} ψ g d x \forall ψ \in L_{0}^{2} (Ω), \\ {u |}_{Γ} & = 0 . \end{matrix}

(13)

Two fundamental properties of the mixed variable system (13) can be immediately observed: the inf-sup condition and the pressure invariance property. Since we aim to find a unique solution for (12), the inf-sup condition must be satisfied. Specifically, this condition reads as

inf_{q \in L_{0}^{2} (Ω) ∖ {0}} sup_{u \in H_{0}^{1} (Ω) ∖ {0}} \frac{(\nabla \cdot u, q)}{{| | \nabla u | |}_{L^{2} (Ω)} {| | q | |}_{L^{2} (Ω)}} \geq β > 0 .

(14)

Differently, the constraint

- \nabla \cdot u = g

does not hold.

Furthermore, the invariance property must be preserved, ensuring that the velocity is not affected by a gradient change in the external body force, namely

f \to f + \nabla ψ \Rightarrow (u, p) \to (u, p + ψ) .

(15)

In fact, an additional force that can be represented purely as the gradient

\nabla ψ

must be balanced by a pressure gradient when no-slip boundary conditions do not affect the pressure.

The violation of (14) or (15) in any discretizations has severe well-known consequences. Many finite element spaces satisfy (14) but not (15).

Consider now a finite element method on a pair of finite element spaces with piecewise polynomials. Let

X_{h} \subset H_{0}^{1} (Ω)

and

Y_{h} \subset L_{0}^{2} (Ω)

denote a pair of conforming spaces for a partition

F_{h}

Ω

parametrized by the characteristic length h. The Stokes system, for the state

(u_{h}, p_{h}) \in (X_{h}, Y_{h})

becomes

\begin{matrix} a (u_{h}, w_{h}) + b (w_{h}, p_{h}) & = (f, w_{h}) & \forall w_{h} \in X_{h} (Ω), \\ b (u_{h}, ψ_{h}) & = (g, ψ_{h}) & \forall ψ_{h} \in Y_{h} (Ω), \end{matrix}

(16)

where

a (u, w) = ν (\nabla u, \nabla w)

and

b (w, p) = - (\nabla \cdot w, p)

. For this system of equations an a priori error estimate reads as

| | u - u_{h} {| |}_{1, h} \leq C_{1} h^{k} {| u |}_{k + 1} + \frac{C_{2}}{ν} h^{k} {| p |}_{k},

(17)

where standard Lagrangian polynomial finite elements are considered. The order k is chosen for the velocity discretization, while an order

k - 1

is set for the discrete pressure [3,4]. The term

\frac{C_{2}}{ν} h^{k} {| p |}_{k}

in (17) represents the drawback when dealing with non-divergence mixed finite elements, showing the connection between the discrete velocity

u_{h}

and the continuous pressure p. The Stokes problem is said to have the mixed-method-robustness property when the pressure cannot influence the velocity approximation. Similarly, we say that an approximation of the Stokes problem does not have mixed-method-robustness property when a change in pressure has a large impact on the velocity field. One can refer to [8] for details on the mixed-method-robustness property for Stokes and Navier-Stokes problems.

In order to have a finite element approximation with the mixed-method-robustness property, we should construct an approximation that takes into account properties (14) and (15). We recall that [8]

Proposition 1.

Let

X_{h} \subset H_{0}^{1} (Ω)

and

Y_{h} \subset L_{0}^{2} (Ω)

be a couple of finite element spaces that satisfy the inf-sup stability property (14) with

\nabla \cdot X_{h} \subset Y_{h}

and let

Π_{Y_{h}} q \in Y_{h}

be the

L^{2}

projection of q. Then, for the solution

(u_{h}, p_{h}) \in (X_{h}, Y_{h})

of (16) we have:

i) the velocity error is pressure independent (mixed-method-robustness property holds)

ii) changing

f \to f + \nabla ψ \Rightarrow (u_{h}, p_{h}) \to (u_{h}, p_{h} + Π_{Y_{h}} ψ)

(the invariance property holds).

Standard finite element couples, such as the Taylor-Hood approximations, satisfy (14) but not (15). The couple consisting of functions in

X_{h} \subset H (d i v, Ω)

, such as linear Raviart-Thomas functions, and its

\nabla \cdot X_{h} \subset Y_{h}

piecewise constant finite element may satisfy such a requirement, and it can be considered to have the mixed-method-robustness property.

3. $H (div, Ω)$ Function Space Approximations

This section introduces the function spaces required in the following for the construction of the divergence-free approximation fields.

3.1. Partition of Domain and $H (div, Ω)$ Functions

In order to construct finite element approximations, fundamental insights into domain partitioning are presented for

H (div, Ω)

space. Establishing a finite element approximation within a domain

Ω

demands maintaining specific continuity properties at the interfaces between its constituent elements.

Consider a domain

Ω

partitioned into subdomains, such that

Ω = \cup_{r = 1}^{m} K_{r}

, where the generic element

K_{r}

can take the form of a triangle or a quadrilateral in the case of a bi-dimensional domain (tetrahedron or hexahedron for three dimensions). We consider only compatible meshes, so no hanging node is present. Additionally, following the literature notation, the mesh size is represented by the index h, which is also used for the maximum diameter of the element. Following [25], we now reformulate the continuous functional spaces previously described in the context of a partitioned domain. Starting with a scalar field q, we define

X (Ω)

X (Ω) : = \{q \in L^{2} {(Ω) : q |}_{K_{i}} \in H^{1} (K_{i}), \forall i \in N, i \leq m\} = \prod_{r} H^{1} (K_{r}),

(18)

equipped with the norm

{| | q | |}_{X (Ω)}^{2} : = \sum_{r} {| | q | |}_{1, K_{r}}^{2}

. Regarding vector-valued fields, we define also the space

Y (Ω)

\begin{matrix} Y (Ω) : = \{v \in L^{2} {(Ω) : v |}_{K_{i}} \in H (div, K_{i}), \forall i \in N, i \leq m\} = \prod_{r} H (div, K_{r}), \end{matrix}

(19)

with the norm

{| | v | |}_{Y (Ω)}^{2} : = \sum_{r} {| | v | |}_{div, K_{r}}^{2}

Consider now a partition of the domain boundary such as

Γ = D \cup N | D \cap N = \emptyset

, where the boundary

Γ

has been divided considering different types of boundary conditions, i.e. D for Dirichlet and N for Neumann. Hence, we define

H_{0, D}^{1} (Ω) : = \{q \in H^{1} {(Ω) : q |}_{D} = 0\},

(20)

where naturally we have

H_{0, D}^{1} (Ω) = H_{0}^{1} (Ω)

D = Γ

and

H_{0, D}^{1} (Ω) = H^{1} (Ω)

D = \emptyset

. Likewise, we have

H_{0, N} (div, Ω) : = \{v \in H (div, Ω) : 〈 v \cdot n, q 〉 = 0, \forall q \in H_{0, D}^{1} (Ω)\},

(21)

with

H_{0} (div, Ω) = H_{0, N} (div, Ω)

N = Γ

. Lastly, we recall another subspace of

H (div, Ω)

defined as

H (div 0, Ω) : = \{v \in H (div, Ω) : \nabla \cdot v = 0\} .

(22)

Considering now a function

q \in H^{1} (Ω)

and a function

v \in H (div, Ω)

, denoting the vector normal to

Γ_{r} = \partial K_{r}

n_{r}

, we have that

\sum_{r} {〈 v \cdot n_{r}, q 〉}_{Γ_{r}} = {〈 v \cdot n, q 〉}_{Γ},

(23)

where the operator

〈 \cdot, \cdot 〉

represents the duality product between

H^{\frac{1}{2}} (Γ_{r})

and

H^{- \frac{1}{2}} (Γ_{r})

. Therefore, inside each element, we can apply the Green formula obtaining

{〈 v \cdot n, q 〉}_{Γ} = \sum_{r} \{\int_{K_{r}} \nabla \cdot v q d x + \int_{K_{r}} v \cdot \nabla q d x\} .

(24)

As described in [25] we have the following proposition:

Proposition 2.

The spaces

H_{0, D}^{1} (Ω)

and

H_{0, N} (div, Ω)

can be expressed as

H_{0, D}^{1} (Ω) = \{q \in X (Ω) : \sum_{r} 〈 v \cdot n_{r}, q 〉 = 0, \forall v \in H_{0, N} (div, Ω)\},

(25)

and

H_{0, N} (div, Ω) = \{v \in Y (Ω) : \sum_{r} 〈 v \cdot n_{r}, q 〉 = 0, \forall q \in H_{0, D}^{1} (Ω)\},

(26)

respectively.

We can conclude that a function

v \in Y (Ω)

is in

H (div, Ω)

if and only if the normal trace is continuous at the interface. In order to preserve the normal trace continuity at the interfaces we need to use Piola’s transformation as the pullback in our finite element approximation.

3.2. Piola’s Transformation

In the framework of finite element approximation, the coordinate changes play a fundamental role, allowing the use of a reference element in which the computation is performed.

Consider a reference domain

\hat{K} \subset R^{n}

. Let

\partial \hat{K}

be its boundary and denote by

\hat{n}

the outward-oriented normal. Additionally, since we are dealing with integral relations, we define

d \hat{x}

the Lebesgue measure on

\hat{K}

and with

d \hat{σ}

the surface measure on

\partial \hat{K}

. We can introduce a smooth mapping

F : R^{n} \to R^{n}

, where smooth implies at least

C^{1}

continuity. Consequently, the mapping, or pullback, between the element and the reference element is defined by

K = F (\hat{K})

. We say that the element K is the image of

\hat{K}

under the diffeomorphism

F

We indicate with

DF (\hat{x})

the Jacobian matrix of the transformation, assuming that

DF (\hat{x})

is invertible for any

\hat{x}

and that

F

is globally invertible on the element K. Therefore, the following relationship holds:

{DF}^{- 1} (\hat{x}) = {(DF (\hat{x}))}^{- 1}

. When the transformation

F (\hat{x})

is linear, i.e.

F (\hat{x}) = x_{0} + B \hat{x}

, the map is said to be affine. In addition, the Jacobian matrix is constant,

DF (\hat{x}) = B

. Finally, we define the determinant of the matrix as

J (\hat{x}) : = | det DF (\hat{x}) |

, for

\hat{x} \in \hat{K}

One of the main features of functions in

H (div, Ω)

is the use of normal components on the element faces as degrees of freedom. In fact, remembering Proposition 2, we understand how a function belonging to

H_{0, N} (div, Ω)

is defined as the summation over every element such that the dot product between the fluxes of the function, i.e.

v \cdot n_{Γ}

, and the test function is equal to zero. On the other hand, the transformation

F

previously described does not preserve the normal components, and besides, we are not able to map

H (div, \hat{K})

into

H (div, K)

. For this reason, it is necessary to introduce the Piola’s transformation, which is a contravariant transformation. For a vector-valued function

\hat{q} \in H^{1} (\hat{K})

we consider the mapping

G (\hat{q}) (x) : = \frac{1}{J (\hat{x})} DF (\hat{x}) \hat{q} (\hat{x}), with x = F (\hat{x}) .

(27)

Since we have an invariance of the trace matrix with a change of variable, the following relation holds:

\nabla \cdot q = J^{- 1} \nabla \cdot \hat{q}

. Following the idea presented in [26,27], we report other important lemmas.

Lemma 1.

Consider the transformations

v = F (\hat{v})

and

q = G (\hat{q})

. The following integral relations hold

\begin{matrix} \int_{K} q \cdot \nabla v d x & = \int_{\hat{K}} \hat{q} \cdot \nabla \hat{v} d \hat{x} & \forall \hat{v} & \in L^{2} (\hat{K}), \end{matrix}

(28)

\begin{matrix} \int_{K} v \nabla \cdot q d x & = \int_{\hat{K}} \hat{v} \nabla \cdot \hat{q} d \hat{x} & \forall \hat{v} & \in L^{2} (\hat{K}), \end{matrix}

(29)

\begin{matrix} \int_{\partial K} q \cdot n v d σ & = \int_{\partial \hat{K}} \hat{q} \cdot \hat{n} \hat{v} d \hat{σ} & \forall \hat{v} & \in L^{2} (\partial \hat{K}) . \end{matrix}

(30)

Note that, from the last relation in Lemma 1, with the new transformation

G

the normal trace in

H^{- \frac{1}{2}}

is preserved.

With the definition of the transformation

F

we can also define the normal vector and the tangent vector on

\partial K

. Considering the unit normal

\hat{n}

and the unit tangent vector

\hat{t}

\partial \hat{K}

, we have that

n (x) = {[DF]}^{- T} \cdot \hat{n} (\hat{x}) / | | {[DF]}^{- T} \cdot \hat{n} (\hat{x}) | |

and

t (x) = DF \cdot \hat{t} (\hat{x}) / | | DF \cdot \hat{t} (\hat{x}) | |

3.3. Lowest-order Raviart-Thomas finite element approximation

In this section, we aim to describe a suitable criterion to understand if a finite element space is a subspace of

H (div, Ω)

. Therefore, considering the shape-regular triangulation

T_{h}

of the domain

Ω

, we define with

ε_{h}

the set of edges or faces, considering respectively a bi-dimensional or a three-dimensional domain. Moreover, we define the set of boundary edges/faces with

ε_{h}^{B} \subset ε_{h}

, namely the variable

e \in ε_{h}^{B}

e \cap Γ \neq 0

, and the set of interior edges/faces with

ε_{h}^{I} : = ε_{h} ∖ ε_{h}^{B}

. An important lemma regarding the normal component of finite element subspaces of

H (div, Ω)

is now reported, the proof can be found in [8].

Lemma 2.

Let

T_{h}

the partitioned domain Ω and consider a space of piecewise polynomials function

w_{h}

. Since

w_{h} \subset H (div, Ω)

, we have continuity across interelement boundaries

e \in ϵ_{h}^{I}

of the normal components.

Naturally, the space

w_{h}

represents in the context of the Stokes equation the finite element space for the velocity field. Furthermore, this lemma does not hold necessarily for the tangential components. Among the finite element spaces satisfying the previous lemma, i.e. conforming subspaces of

H (div, Ω)

, we can introduce the Raviart-Thomas space of order

k \geq 0

(

{RT}_{k}

) [28,29,30,31]. Given an element K of the partitioned domain

Ω

, the local Raviart-Thomas space of order

k \geq 0

is defined as

{RT}_{k} (K) = P_{k} {(K)}^{n} + x P_{k} (K) .

(31)

Therefore, we define the spaces as

\begin{matrix} {RT}_{k} & : = \{w_{h} \in H_{0} (div, Ω) : w_{h} |_{K} \in {RT}_{k} (K) \forall K \in T_{h}\}, \end{matrix}

(32)

where, we have used the space

H_{0} (div, Ω) = \{{v \in H (div, Ω) : v \cdot n |}_{Γ} = 0\}

. Moreover,

P_{k} = P_{k}^{n}

(

P_{k}

) represents the space of globally continuous vector (scalar)-valued piecewise polynomials of degree not exceeding k.

4. Divergence-Free Field for Two-Phase Flow

In this section, the first physical problem is discussed, i.e. the mass conservation issue in the context of multiphase flow problems.

4.1. Problem Model

A key point for two-phase applications is the conservation of mass. In order to introduce a typical two-phase problem, we consider an incompressible, viscous two-phase flow over an open bounded domain

Ω

with continuous boundary

Γ = \partial Ω

, consisting of two subdomains

Ω_{1} \subset Ω

and

Ω_{2} \subset Ω

, representing the Phase 1 and Phase 2, respectively. We can also define the indicator function

χ

χ (x, t)

that is one over

Ω_{1}

and zero otherwise. We may assume a single velocity

u

and pressure p field for both phases, the fluid motion to obtain the one-fluid model equations governing the multiphase flow motion, which read

\begin{matrix} ρ (\frac{\partial u}{\partial t} + u \cdot \nabla u) = - \nabla p + \nabla \cdot [μ (\nabla u + \nabla u^{T})] + f, \end{matrix}

(33)

\begin{matrix} \nabla \cdot u = 0, \end{matrix}

(34)

\begin{matrix} \frac{\partial χ}{\partial t} + u \cdot \nabla χ = 0, \end{matrix}

(35)

where (33) represents the conservation of linear momentum, (34) the conservation of volume (the phases are non-miscible and incompressible), and (35) the advection of the indicator function

χ

f

is the volumetric or the surface tension force. The two considered fluids, the reference phase 1 and the secondary phase 2 have both constant densities

ρ_{1}

and

ρ_{2}

and viscosities

μ_{1}

and

μ_{2}

, within each subdomain

Ω_{1}

and

Ω_{2}

, respectively. Therefore, we set

ρ = χ ρ_{1} + (1 - χ) ρ_{2}

, and

μ = χ μ_{1} + (1 - χ) μ_{2}

Suppose that we use Taylor-Hood finite element approximation for the velocity-pressure field

(u_{h}^{*}, p_{h}) \in X_{h} \times W_{h} \subset H^{1} (Ω) \times L_{0}^{2} (Ω)

. Once the Navier-Stokes equation is solved with this finite element pair, the discrete solution

u_{h}^{*}

is not a divergence-free function at each point and a projection into

H (div 0, Ω)

is very important for the conservation of the phase indicator. Thus, if we have

u^{*} \in H^{1}

, the objective is to find the velocity

u \in H (div 0, Ω)

by minimizing

\begin{matrix} F (u) & = \frac{1}{2} \int_{Ω} {(u - u^{*})}^{2} d Ω, \\ \nabla \cdot u & = 0, \end{matrix}

(36)

over the linear function subspace

H (div, Ω)

. This is equivalent to project

u^{*}

into

H (div 0, Ω)

by using the Helmholtz-Hodge decomposition (HDD) considering the Raviart-Thomas representation for

u

[8]. We recall now a lemma for the HDD.

Lemma 3.

Let Ω be a connected domain and

u^{*}

be in

L^{2} (Ω)

, there exists a vector field

u \in H (div 0, Ω)

and a scalar function

p \in H^{1} (Ω)

with

\begin{matrix} i) & u^{*} = u + \nabla p \end{matrix}

(37)

\begin{matrix} i i) & \nabla \cdot u = 0 \end{matrix}

(38)

\begin{matrix} i i i) & (u, \nabla φ) = 0, \forall φ \in H^{1} (Ω) . \end{matrix}

(39)

As a result of the minimization problem, we obtain the standard weak formulation

\begin{matrix} \int_{Ω} u \cdot v d x - \int_{Ω} k \nabla \cdot v d x & = \int_{Ω} u^{*} \cdot v d x & \forall v \in H (div, Ω), \\ \int_{Ω} q \nabla \cdot u d x & = 0 & \forall q \in L^{2} (Ω), \end{matrix}

(40)

where k represents the lagrange multiplier associated with the incompressibility constraint. It is important to remark that iii) in Lemma 3 is always satisfied when

u \in H (div 0, Ω)

with correct boundary conditions.

4.2. Approximation over Cartesian Cells

In order to define velocity fluxes and divergence operators in the interior of a domain one needs to define an internal domain structure. Every approximation, preserving free-divergence properties, needs an appropriate definition over its substructure domain. Let

Ω

be a polyhedron with mesh

T_{h}

that consists of N subdomains

K_{k}

(

k = 1, . . N

). Each element is assumed in the form of a polyhedron with four and six faces for two and three-dimensional geometry, respectively. Its maximum diameter h is obtained from the reference element

\hat{K}

through the mapping

F_{k} : \hat{K} \to K_{k}

. As shown in Figure 1 on the right, the reference element

\hat{K}

is assumed to be a square (two-dimensional geometry) or a cube (three-dimensional geometry). We recall the Jacobian matrix

{DF}_{k}

associated to each element and its determinant

J_{k}

. The mesh is a shape regular and non-degenerate mesh by assuming standard hypotheses. For details, the reader can refer to [32].

R^{2}

the partition of

Ω

consists of quadrilateral elements

K_{k}

of maximum diameter h obtained from the reference element

\hat{K} = {[- 1, 1]}^{2}

in agreement with standard piece-wise linear finite elements. As shown in Figure 1 on the left, we denote with

x_{j}^{k}

the vertex of the quadrilateral elements, with

x_{f M}^{k}

the middle point of the face f and

n_{\partial K_{k}^{f}}^{k}

its unit normal for

f = 1, \dots, 4

. In this case, each point

\hat{x} = ({\hat{x}}_{1}, {\hat{x}}_{2}) \in \hat{K} \subset R^{2}

sets a point

x = (x_{1}, x_{2}) \in K_{k} \subset R^{2}

through the mapping

F_{k} : \hat{K} \to K_{k} \subset R^{2}

defined by

x = F_{k} (\hat{x}) = \sum_{j = 1}^{4} x_{j}^{k} ϕ_{j} (\hat{x}) = x_{0}^{k} + {\hat{x}}_{1} t_{1, 0}^{k} + {\hat{x}}_{2} t_{2, 0}^{k} + {\hat{x}}_{1} {\hat{x}}_{2} d_{1, 2}^{k},

(41)

for all

\hat{x} \in \hat{K}

. The

ϕ_{j}

are the standard linear basis functions for QUAD4 finite element. The vectors

t_{1, 0}^{k} = (x_{2 M}^{k} - x_{M}^{k}) / 2

and

t_{0, 2}^{k} = (x_{3 M}^{k} - x_{M}^{k}) / 2

define the two middle point lines, while

x_{M}^{k}

is the quadrilateral middle point and

d_{1, 2}^{k} = (x_{3}^{k} + x_{1}^{k} - x_{2}^{k} - x_{4}^{k}) / 4

. For

d_{1, 2}^{k} = 0

the quadrilateral is a parallelogram. The Jacobian matrix

{DF}_{k}^{T}

becomes

{DF}_{k}^{T} (\hat{x}) = [t_{1, 0}^{k} + {\hat{x}}_{2} d_{1, 2}^{k}, t_{2, 0}^{k} + {\hat{x}}_{1} d_{1, 2}^{k}] .

From the previous definitions, the global Raviart-Thomas finite element spaces can be introduced. In this case, the mesh-size h has to be defined as

h = {max}_{K \in T_{h}} h_{K}

. Moreover, we assume that the hypothesis of regularity condition is satisfied in any element

K \in T_{h}

. Therefore, the global space can be defined as

{RT}_{0} (T_{h}) : = {v \in H (div, Ω) : v |_{K} \in {RT}_{0} (K) \forall K \in T_{h}},

(42)

where

{RT}_{0} (K)

is defined as the space spanned by the basis functions. For two-dimensional quadrilateral

\hat{K}

we set

v = \sum_{i = 1}^{4} p_{f} {\hat{b}}_{f}

with

\begin{matrix} {\hat{b}}_{1} & = \frac{{\hat{x}}_{2} - 1}{4} {\hat{i}}_{2}, & {\hat{b}}_{2} & = \frac{{\hat{x}}_{1} + 1}{4} {\hat{i}}_{1}, & {\hat{b}}_{3} & = \frac{{\hat{x}}_{2} + 1}{4} {\hat{i}}_{2}, & {\hat{b}}_{4} & = \frac{{\hat{x}}_{1} - 1}{4} {\hat{i}}_{1}, \end{matrix}

(43)

where

{\hat{i}}_{1} = {(1, 0)}^{T}

and

{\hat{i}}_{2} = {(0, 1)}^{T}

. The vectors are oriented in agreement with the outer direction. In

{RT}_{0} (K)

, applying the divergence theorem, we know that the continuity of the normal component across the edges elements is a necessary condition for a piecewise polynomial vector-valued function to be in

H (div, Ω)

In order to complete the description of the finite element subspaces required for the approximation of the mixed problem, it is important to introduce the approximating subspace of the scalar variable. Note that (40) does not have any derivatives of the pressure and the continuity requirement does not hold. Given these considerations, we introduce the standard space of piecewise polynomial functions with degree

k = 0

, not necessarily continuous, as

P_{k} (T_{h}) : = {q \in L^{2} (Ω) : q |_{K} \in P_{k} (K) \forall K \in T_{h}},

(44)

where the superscript d considers the possibility of discontinuous functions. With these spaces, the following identity holds:

\nabla \cdot {RT}_{0} = P_{0}

Having introduced the finite element subspace for the approximation of the elliptic mixed problem (40), we introduce the finite element discretization for the solution pair

(u_{h}, p_{h}) \in {RT}_{0} \times P_{k}

\begin{matrix} \int_{Ω} u_{h} \cdot v d x - \int_{Ω} p_{h} \nabla \cdot v d x & = \int_{Ω} u_{h}^{*} \cdot v d x & \forall v \in {RT}_{0} \\ \int_{Ω} q \nabla \cdot u_{h} d x & = \int_{Ω} f q d x & \forall q \in P_{0} . \end{matrix}

(45)

4.3. Numerical Tests

Now, we consider numerical results related to the multiphase advection problem under a discrete divergence-free field, with markers and interface tracking. Since we exploit the algorithm already described and validated in [33], our attention is focused only on the velocity interpolation of the marker position to satisfy the orthogonal decomposition reported in Lemma 3.

Flow fields characterized by uniform translations and rotations are usually employed to move and advect the cell markers. In fact, simple advection movements are designed to displace fluid bodies within a domain and verify the conservation of surface shape and volume. For this reason, standard tests from the literature have been considered as benchmark results to evaluate the goodness of interface tracking algorithm [34,35,36,37].

All the tests presented are two-dimensional problems, with a velocity field designed so that the resulting vorticity is not uniform in the domain. Consequently, significant distortions occur in the fluid body interface, making the maintenance of the interface not straightforward. Finally, cosinusoidal time-dependence analytical velocity functions are used to return to the initial configuration at the end of the period [38].

The

L_{1}

error norms are considered to compare the results [34]. The relative mass error

E_{m} (t_{1})

is defined to compare the total volume of a phase, specifically the reference volume, at the initial time

t_{0}

and the subsequent time

t_{1}

E_{m} (t_{1}) = \frac{| \sum_{i = 1}^{N_{e l}} A_{i} C_{i} (t_{1}) - A_{i} C_{i} (t_{0}) |}{\sum_{i = 1}^{N_{e l}} A_{i} C_{i} (t_{0})} .

(46)

In this case, the color function value at the cell i at time t is represented by

C_{i} (t)

A_{i}

represents the area of the cell i meanwhile the total number of the cells is

N_{e l}

. Another error, termed geometrical error

E_{g} (t_{1})

, is introduced as

E_{g} (t_{1}) = \sum_{i = 1}^{N_{e l}} A_{i} | C_{i} (t_{1}) - C_{i} (t_{0}) | .

(47)

The objective is to verify whether the final shape aligns with the initial configuration. Therefore, we have computed another type of error, considering a circular geometry as the initial shape. Therefore, given the center

(x_{c}, y_{c})

and the radius R, the distance between a marker m with position

(x_{m}, y_{m})

and the center is computed and compared with the radius R

E_{a l} = \sum_{m = 1}^{N_{m}} |\sqrt{{(x_{m} - x_{c})}^{2} + {(y_{m} - y_{c})}^{2}} - R| s_{m} .

(48)

The quantity

s_{m}

is the arc length. In all computations presented, the

C F L = u Δ t / h

has been considered, with u representing the velocity component along the x-axis,

Δ t

the time step simulation and h the grid spacing. Each variable presented above can be computed by considering the algorithm implemented in [39], whose description is omitted since it is beyond the scope of this work.

Furthermore, an initial attempt to incorporate Raviart-Thomas basis functions into the surface marker reconstruction framework has been investigated. The primary goal is to compare the advection test using two different representations of the velocity field.

Indeed, the numerical approximation of the velocity is naturally stored through finite element discretization, such as the classical 9 points of a biquadratic quadrilateral element. In the same way, the analytical velocity is stored based on the specific data structure arising from mesh discretization. The exact values of that field are found only in the standard degrees of freedom of the mesh elements. Notably, the library for marker reconstruction is designed to handle both kinematic and dynamic two-phase flow simulations, accommodating scenarios where the velocity field is either imposed or fully solved.

In situations where we need to compute the velocity at points beyond the nodes of the elements (such as marker coordinates), interpolation is used. This is essential for marker advection, where the Runge-Kutta method requires varying velocity values at different positions. This section compares two types of finite element interpolation techniques for determining the velocity field at marker locations. The standard approach involves Lagrangian interpolation, while the alternative method employs Raviart-Thomas interpolation.

Denoting

u_{m r k}

as the velocity field of the i-th marker in the cell with coordinates

x_{p}

, the comparison is performed using the following formulations

\begin{matrix} u_{m r k} & = \sum_{i = 1}^{n_{d o f}} u_{i} φ_{i} (x_{p}), & u_{m r k} & = \sum_{i = 1}^{n_{f a c e s}} b_{i} (x_{p}) p_{i}, \end{matrix}

(49)

where

n_{d o f}

represents the biquadratic nodes for a quadrilateral element, and

n_{f a c e s}

is the number of faces of the same element. The functions

φ_{i}

and

b_{i}

represent the biquadratic Lagrangian basis functions and the Raviart-Thomas basis functions that act on the fluxes faces

p_{i}

, respectively. Since the computational domain is a classical two-dimensional Cartesian mesh, any convergence issues related to using the lowest-order Raviart-Thomas finite elements do not affect this context.

The reason behind the use of the Raviart-Thomas basis function for the velocity field interpolation is motivated by two closely connected aspects:

Mass Conservation in Multiphase Flow: in multiphase flow simulations, maintaining mass conservation is crucial. Since the density of both phases is assumed constant, it is imperative to satisfy the mass conservation equation. The divergence-free velocity constraint should be managed through appropriate discretization techniques. The use of a divergence-free representation of the numerical velocity, facilitated by Raviart-Thomas finite element discretization, addresses this constraint effectively;
Divergence-Free Analytical Velocity Field: The analytical velocity fields implemented for advection tests are sinusoidal functions derived from stream functions, which are, by definition, divergence-free. Therefore, approximating these velocity fields using $H_{0} (div)$ basis functions provides a more realistic representation at each physical point. This approach helps avoid approximation errors that may arise when using standard finite element Lagrangian interpolation.

From a computational perspective, only the function responsible for the Runge-Kutta advection has been modified to allow both interpolations. The application has been developed by integrating two finite element libraries: the FEMuSTTU library [39], where routines for surface marker reconstruction are implemented, and the ProXPDE library [40], which provides the Raviart-Thomas finite element interpolation.

The approach for testing the Raviart-Thomas basis functions for velocity field interpolation involves a single bubble advection test. The analytical velocity field for this simulation is given by

\begin{matrix} u & = 2 {sin}^{2} (π x) sin (π y) cos (π y) cos (\frac{π t}{T}), \\ v & = - 2 sin (π x) cos (π x) {sin}^{2} (π y) cos (\frac{π t}{T}) . \end{matrix}

The domain considered is represented by

Ω = [- L / 2, L / 2] \times [- H / 2, H / 2]

, with

H = L = 1

, and the initial circular geometry is located at

(0, 0.25)

with a radius

R = 0.15

. The period T for this test is set to

4 s

, implying that the maximum stretch of the bubble is reached at

t = 2 s

. The purpose is to compare the performance of standard Lagrangian interpolation with Raviart-Thomas interpolation for the velocity field with surface marker reconstruction algorithms. The simulation evaluates how well each interpolation scheme maintains the accuracy and conservation properties of the marker positions during advection.

In Table 1, the error values for the Lagrangian-type interpolation, denoted with the superscript Q, are presented based on the number of mesh elements

N_{e l}

, ranging from a

32 \times 32

grid to a

128 \times 128

grid. Every computed error shows a decrease by increasing the mesh refinement, showing the same behavior as the previous tests. We notice that the convergence error rates for

p_{m}

and

p_{g}

reach a value close to 2, while

p_{a l}

seems to have a slightly larger value. We recall that the

E_{a l}

error provides a geometrical measure of the difference between the initial configuration and the corresponding final shape. For this reason, an order of convergence of

p_{a l}

between 2 and 3 indicates a good behavior of the algorithm.

The same error analysis has been performed for the

{RT}_{0}

velocity interpolation, and the results are reported in Table 2. Again, convergence errors decrease with the grid refinement. Regarding the rates of convergence, similar values have been obtained for the

{RT}_{0}

, except for

p_{a l}

which seems to be slightly lower than the one computed with the lagrangian interpolation.

A comparison between the two techniques shows that the

Q_{2}

Lagrangian interpolation provides better error results. The order of magnitude is almost lower by one compared to the corresponding

{RT}_{0}

interpolation on the same grid. However, we can note that these errors consider the entire simulation and the final values of the color function. As such, these errors, as defined, do not serve as a direct indicator of the approximation quality of the velocity.

The numerical algorithm incorporates several functions for surface marker reconstruction. Therefore, while the velocity interpolations of the marker may be similar inside the cell, a velocity field variation can lead to changes in the final position of the marker with the Runge-Kutta advection scheme.

In Figure 2, the final positions of the markers are reported for different grid refinements obtained with the standard Lagrangian

Q_{2}

interpolation. Some discrepancies from the analytical circular geometry are present in the marker positions only for the coarse grid (

32 \times 32

). As the grid refinement increases, the markers reach an end configuration in better agreement with the initial shape.

The same comments can be made for the

{RT}_{0}

interpolation, which is shown in Figure 3. The final marker positions match the initial circular geometry, except for the

32 \times 32

grid resolution. It is important to note that this aspect is not merely a graphical consideration; it is represented by the

E_{a l}

error, which measures the discrepancies from the initial shape. For both interpolation techniques, we observe orders of magnitude ranging from

10^{- 3}

up to

10^{- 5}

, confirming the algorithm’s behavior in regaining the initial circumference.

As previously described, it is crucial to verify the shape configuration when the maximum deformation is reached. For this reason, we now present the bubble configuration at half of the period T, i.e., at

t = 2

seconds, to check if the closed surface is preserved. Despite the thin filaments in the bubble tail, we observe that the algorithm can reconstruct the surface and rebuild the marker positions. This situation is depicted for both techniques, in Figure 4 for the

Q_{2}

interpolation and in Figure 5 for the

{RT}_{0}

interpolation, across the three tested grids.

5. Projection Method for Navier-Stokes Equation

This section recalls the framework of the projection methods used for solving the Navier-Stokes equations in incompressible flow simulations, also known in the literature as the fractional step method.

5.1. Model and Algorithm

This method broadly falls into three schemes: pressure-correction, velocity-correction, and consistent splitting methods. In particular, some steps present in the projection algorithm are essentially the same as the orthogonal decomposition of the velocity. Therefore, since the decomposition is treated by using the mixed finite element tools presented in Section 3.3, we establish a connection to fully resolve the Navier-Stokes system of equations. For an exhaustive review of this method, the interested reader can consult the work of Guermond et al. [23].

The projection method originates from the incompressibility constraint, which represents the main motivation of this work. In particular, one of the main issues related to the numerical solution of momentum and mass conservation equations is the coupling between velocity and pressure, imposed by the condition of zero divergence of velocity. The solution of the coupled velocity-pressure system can be expensive from a computational point of view, and this has led to the development of different numerical algorithms for the treatment of the split system for velocity and the pressure fields. In order to overcome these problems, many authors introduced the projection method by which the Navier-Stokes system is subdivided into two separate steps, one for the velocity and one for the pressure.

The first idea that appeared in literature is related to the work of Chorin and Temam [21,22], where the time-dependent solution of incompressible viscous flow has been proposed. Specifically, the idea is to build a sequence of decoupled elliptic equations for the velocity and the pressure and solve them at each time step. Indeed, this numerical method leads to an efficient simulation and reduces the computational effort.

Hereafter, the algorithm employed for the resolution of the Navier-Stokes system in the multigrid finite element library FEMuS [41] is described, considering the standard incremental pressure correction scheme. Thus, we recall the Navier-Stokes system for a generic incompressible fluid flow, with constant physical properties and a smooth force

f

\begin{matrix} \frac{\partial u}{\partial t} + (u \cdot \nabla) u & = - \nabla p + ν Δ u + f, \\ \nabla \cdot u & = 0, \end{matrix}

(50)

where homogeneous Dirichlet boundary conditions are taken into account.

In order to define the time-dependent solutions, we can consider the time step

Δ t > 0

and set

t^{k} = k Δ t

for

0 \leq k \leq K = [T / Δ t]

, where T is the upper bound of the time interval considered,

[0, T]

. To describe the adopted pressure-velocity split algorithm, the simple Euler scheme for the time discretization is used. Specifically, a fictitious intermediate velocity field

u^{n *}

is introduced leading to

\frac{u^{n *} - u^{n - 1 *}}{Δ t} + \frac{u^{n} - u^{n - 1}}{Δ t} + (u^{n} \cdot \nabla) u^{n} = - \nabla p^{n} + ν Δ u^{n} + f .

(51)

After that, following a standard approach, the equation can be subdivided into two different equations by introducing the incremental pressure

δ p

, and considering

p^{n} = p^{n - 1} + δ p^{n}

\begin{matrix} \frac{u^{n *} - u^{n - 1}}{Δ t} & = - (u^{n} \cdot \nabla) u^{n} - \nabla p^{n - 1} + ν Δ u^{n} + f, \end{matrix}

(52)

\begin{matrix} \frac{u^{n} - u^{n *}}{Δ t} & = - \nabla δ p^{n} . \end{matrix}

(53)

The first equation considers the viscous effect solving the velocity

u^{n *}

, which has a divergence different from zero, while the second equation represents the

L^{2}

projection between the two velocity fields, where the

\nabla \cdot u = 0

constraint is understood.

In this work, a different velocity-pressure split is proposed and analyzed with the objective of exploiting the orthogonal decomposition of the velocity field [42]. Specifically, the first equation, which solves the intermediate velocity

u^{n *}

remains the same, while the pressure equation is not transformed by applying the divergence operator as in standard split algorithms. Thus, for the ’pressure equation’ we solve directly the coupled system

\{\begin{matrix} u^{n} - u^{n *} & = - Δ t \nabla δ p^{n} \\ \nabla \cdot u^{n} & = 0 . \end{matrix}

(54)

In this way, we avoid the resolution of the Poisson equation for the incremental pressure, which would affect the imposition of the well-known non-physical homogenous Neumann boundary condition for

δ p^{n}

Therefore, we apply the orthogonal decomposition to (54), by using the properties of the Raviart-Thomas finite element family. In particular, (54) is solved by considering the same optimal minimization problem presented for the multiphase flow problem (36). The objective is to find a field

u

which is the closest velocity to

u^{*}

, under the constraint of the divergence equal to zero.

Differently, for the velocity

u^{*}

the solution space can be found considering the standard solution of the Navier-Stokes equations. In particular, consider the Taylor-Hood finite element space

X_{h} \subset H^{1} (Ω)

. Thus, if we have

u_{h}^{*} \in X_{h}

, the objective is to find the velocity

u_{h} \in {RT}_{0} \subset H (div, Ω)

by minimizing

\begin{matrix} F (u_{h}) & = \frac{1}{2} \int_{Ω} {(u_{h} - u_{h}^{*})}^{2} d Ω, \\ \nabla \cdot u_{h} & = 0, \end{matrix}

(55)

over the linear function subspace

{RT}_{0} \subseteq H (div, Ω)

Naturally, for

u_{h}^{*}

, the velocity approximation is defined as

u_{h}^{*} = \sum_{j} u_{j h}^{*} φ_{j} (x)

, where

u_{j h}^{*}

is the velocity field at the points j and

φ_{j} (x)

indicates the Lagrangian quadratic polynomial basis functions. Instead, for

u_{h}

, the approximation velocity reads as

u_{h} = \sum_{f} p_{f} N_{f} (x)

, where

p_{f}

represents the fluxes through the faces and

n_{f} (x)

corresponds to the Raviart-Thomas vector basis functions.

In order to derive the final system to be solved, we can consider now a piecewise constant pressure, i.e.

k_{h} = δ p_{h} \in S_{h} \subset L^{2} (Ω)

. Therefore, the system to be solved is described with the following two equations

\begin{matrix} \int_{Ω} δ k_{h} \nabla \cdot u_{h}^{n} d Ω = 0 & \forall δ k_{h} \in S_{h} \subset L^{2} (Ω), \\ \int_{Ω} (u_{h}^{n} - u_{h}^{* n}) δ u_{h}^{n} d Ω + \int_{Ω} k_{h} \nabla \cdot δ u_{h}^{n} d Ω = 0 & \forall δ u_{h}^{n} \in {RT}_{0} (Ω), \end{matrix}

(56)

where

k_{h}

can be interpreted as the discrete Lagrange multiplier of the divergence equation.

It is worth noting that (56) is only a step of the split algorithm. Indeed, the elliptic equation for the velocity

u_{h}^{* n}

in (52) is still solved by considering standard Lagrangian finite elements, in particular the Taylor-Hood type in our case.

We recall that by using the lowest-order

RT

element the velocity is determined considering only the fluxes through the element edges, i.e. four (six) degrees of freedom for a bidimensional (three-dimensional) domain. Naturally, in this discussion, only quadrilateral and hexahedral elements have been considered. Hence, the local matrix is equipped only with five (seven) rows that represent the four (six) fluxes through the faces and the central value for the pressure field, instead of using, for example, the nine nodes of a biquadratic quadrilateral with Lagrangian basis functions.

5.2. Numerical Tests

In this section, we aim to test the employment of the Raviart-Thomas finite element family in the framework of the projection method for the resolution of the Navier-Stokes equation. Specifically, (56) has been adopted where we recall that the orthogonal decomposition of the velocity field is considered as the second step of the split algorithm, i.e. the resolution of the pressure equation.

Three cases have been investigated, considering a bidimensional and a three-dimensional geometry: a channel, both in two and three dimensions and a bidimensional cavity. For these tests, the same comparison has been performed, considering three different algorithms for the numerical solution. The first one consists of a standard coupled algorithm for the resolution of the velocity and pressure, and by using classical lagrangian finite elements

P_{2} - P_{1}

for the field discretization. The second technique is based on the same finite element family but employs a standard projection algorithm for the Navier-Stokes system. Lastly, the third method is characterized by the employment of the Raviart-Thomas finite element family for the resolution of the pressure equation in the context of the split technique.

The first analyzed case focuses on the resolution of the Navier-Stokes equation characterized by a low Reynolds number inside a bidimensional channel, i.e. a laminar flow. It is well known that the solution of this kind of configuration is expressed by the Poiseuille profile for the streamwise component of the velocity. Indeed, we expect to obtain the classical parabolic profile.

In Figure 6 the investigated mesh discretizations of the channel geometry are reported. Note that, on the right, is reported a channel discretization where non-affine elements have been employed. On the other hand, the multigrid refinement produces elements that converge to a parallelogram shape, which is an affine element since the opposite edges are parallel, avoiding any convergence issue.

Considering the boundary conditions, at the inlet section

Γ_{i}

a fixed velocity has been imposed, a standard no-slip boundary condition has been imposed on the wall side

Γ_{w}

, while an outlet-type condition has been imposed at the outlet section

Γ_{0}

to fix the pressure value. For this kind of numerical simulation, it is not possible to have an analytical solution of the velocity field. In fact, even though a Poiseuille flow type is searched, if we consider an inlet boundary condition the analytical solution for the Navier-Stokes system does not exist. For this reason, the velocity error norm with respect to a reference solution does not have a numerical significance. On the other hand, in order to compare qualitatively the solutions with the three types of numerical discretizations the

L^{2}

norm of the velocity field over the entire domain is reported, for different grid refinements. In particular, for both regular and irregular meshes, we compute the velocity norm as reported in Table 3. The different methods are denoted with

P^{c}

P^{s}

respectively, for the coupled and split technique with Taylor-Hood finite elements, while with

{RT}_{0}

the Raviart-Thomas finite elements.

We can notice that with the increasing of the grid refinement, the velocity norm values tend to the same value for every method, ensuring reliable numerical results for this test. Consider that, even if the second mesh is characterized by non-affine quadrilateral elements, the numerical solution is sought considering a multigrid approach. Therefore, the initial irregular mesh tends to asymptotically affine quadrilateral elements, allowing the use of the standard Raviart-Thomas family of order 0. The last column of Table 3 confirms that, with this geometry, this spatial discretization is able to find the right numerical solution.

In addition, considering the framework of a laminar flow inside a channel, the main variable of interest is the streamwise component of the velocity field, denoted with v. For this reason, in Figure 7 the velocity profile of v is reported for the three different algorithms, as a function of the x coordinates. The plot has been done, considering a fixed y equal to 1. Since the numerical solutions between the two different grids are very similar, we report only the case with regular quadrilateral elements.

Considering Figure 7, the type of algorithm to solve the Navier-Stokes system is denoted with the subscript, c for the coupled system and s for the split one. With the superscript is denoted the type of finite element family employed for the velocity-pressure discretization,

P

for standard

P_{2} - P_{1}

lagrangian elements and with

RT

the Raviart-Thomas elements. The solid line represents the reference numerical result, that is, the numerical solution obtained with a coupled algorithm and classical Taylor-Hood lagrangian elements. With the markers, the numerical solutions of the split system are reported: the triangular markers represent the case where

P_{2} - P_{1}

elements have been employed in every equation, while the circular markers represent the solution obtained by using an

RT

approximation for the resolution of the pressure equation. We can notice a good agreement of the velocity profiles obtained by employing the three different algorithms.

In the second test, we want to verify the convergence error rate, by considering a cavity configuration, i.e. a rectangular closed geometry. Specifically, the domain is described with the same channel of the previous cases, i.e. the bidimensional channel described in Figure 6, where the elements employed for the spatial discretization are standard regular quadrilateral elements. Regarding the velocity field, the boundary condition imposed on every edge is a homogeneous Dirichlet boundary condition in order to have

u = 0

on the walls.

In order to compute the

L^{2}

norm of the velocity error, the steady exact Navier-Stokes solution has been imposed on the right-hand side of the equation. In fact, given a generic operator A which represents the left-hand side terms of the Navier-Stokes equation, we aim to solve

A u = A u^{*}

(57)

for a specified

u^{*}

, which represents the desired solution. Specifically, the exact solution for the velocity components reads as

\begin{matrix} u^{*} = [\begin{matrix} c u^{*} \\ v^{*} \end{matrix}] = [\begin{matrix} c \frac{π}{2} {sin}^{2} (π x) sin (π \frac{y}{2}) cos (π \frac{y}{2}) \\ - π {sin}^{2} (π \frac{y}{2}) sin (π x) cos (π x) \end{matrix}] . \end{matrix}

(58)

In Table 4 the velocity error norm and the corresponding order of convergence for different levels of refinement are reported, for the

{RT}_{0}

finite element approximation. The order of convergence p is reported in the last column, computed as

p = l n (ε_{l - 1} / ε_{l}) / l n (2)

. We can notice a linear trend for the velocity error norm, as expected by theoretical error estimates for the Raviart-Thomas velocity approximation.

The same test has also been computed considering the standard Taylor-Hood lagrangian basis function, both with a coupled and a split algorithm. The results are reported in Table 5. We can notice a good convergence trend for the velocity error for both methods, even though the order of convergence p is not reported. In fact, despite these parameters appearing to be equal to 3 for both simulations, we recall that the velocity error should be considered together with the pressure error norm. For coupled systems or for lagrangian-type basis functions, that are not pointwise divergence-free, we know that an error in the pressure field produces an effect also the error in the velocity.

An interesting result has been found considering the behavior of the velocity divergence. As expected, for the coupled system solved by using

{RT}_{0}

approximation, the velocity divergence reaches values close to machine precision, indicating that this technique is equipped with an exact zero divergence in every point. Considering the other two methods, the results have been reported in Table 6.

It is worth noting that the velocity divergence error is different from zero as expected. On the other hand, these values seem to converge with a quadratic order.

The last test represents the extension of the bidimensional Poiseuille flow to a three-dimensional configuration, i.e. a regular parallelepiped. The considered domain is represented in Figure 8, where the characteristic lengths have the same values (

L_{x} = 1, L_{y} = 1

and

L_{z} = 4

). Moreover, regular hexahedral elements have been employed for the domain discretization.

The flow configuration follows the same path of the bidimensional case, and thus the fluid enters at the inlet section

Γ_{i}

, and exits through the outlet section

Γ_{o}

. On the remaining boundaries, the standard no-slip boundary condition has been imposed.

Regarding the reliability of the numerical solutions, the same comments of the two-dimensional channel case can also be drawn for the case of a three-dimensional channel. Therefore, since it is not possible to compare the numerical solution with an analytical field, also in this case the

L^{2}

-norm of the velocity field has been computed for the three methods as an indicator of the solution goodness. For this reason, in Table 7 the velocity norm values have been reported considering three levels of refinements, for the three techniques. The notation is the same as the two-dimensional channel test.

The

L^{2}

-norm of the velocity field tends to the same value for each different algorithm, confirming the good behavior of the numerical solution, as already shown in the bidimensional case.

Also in this case, the same qualitative comparison between three types of algorithms has been performed on the numerical solution. Therefore, in Figure 9 the streamwise velocity component w has been reported, where the employed notation and symbols are the same as the bidimensional channel. The w component is represented as a function of the x coordinate, with a plot performed with a fixed y coordinate equal to

0.5

Figure 9 shows good accuracy for the computed numerical solution. Indeed, the velocity profile w obtained with a Raviart-Thomas approximation of the pressure equation in a split system (56) has a similar trend to the velocity profile obtained with a standard coupled (50) and split algorithm (52).

6. Conclusions

This work focuses on the challenges posed by mass conservation in incompressible flow simulations, employing Raviart-Thomas basis functions. The main goal was to achieve a divergence-free velocity field, crucial for accurate numerical solutions in many applications, ranging from multiphase to porous-media flows. An overview of finite element discretization has been presented, emphasizing the use of Raviart-Thomas finite elements built over quadrilateral and hexahedral elements. Specifically, this work has taken into account only the lowest-order Raviart-Thomas elements.

A new velocity-pressure projection method has been presented to obtain a divergence-free velocity field for the resolution of the Navier-Stokes system. Numerical results concerning the orthogonal velocity decomposition and the employment of the projection method have been shown. Convergence rates for velocity and pressure error have been proven to agree with the theoretical convergence rates. The methodology guarantees robustness against spillover of the divergence constraint into the velocity error, the so-called pressure invariance property.

Multiphase flow simulations have been exploited to demonstrate the methodology’s performance, emphasizing the significance of maintaining a divergence-free velocity in scenarios where mass conservation is critical, such as interface tracking problems. An interesting comparison has been performed, exploiting specific characteristics of the Lagrangian and Raviart-Thomas finite elements.

In summary, this work has contributed to the understanding of maintaining divergence-free fields in incompressible flow simulations. The Raviart-Thomas finite element family emerged as a valuable tool in addressing this challenge, with applications ranging from the projection method for the Navier-Stokes equations to the advection of interfaces in multiphase flows.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflicts of interest.

References

Raviart, P.A.; Thomas, J.M. Introduction à l’analyse numérique des équations aux dérivées partielles; Masson, 1983.
Linke, A. On the role of the Helmholtz decomposition in mixed methods for incompressible flows and a new variational crime. Computer methods in applied mechanics and engineering 2014, 268, 782–800. [Google Scholar] [CrossRef]
Girault, V.; Raviart, P.A. Finite element methods for Navier-Stokes equations: theory and algorithms. NASA STI/Recon Technical Report A 1986, 87, 52227. [Google Scholar]
Fortin, M.; Brezzi, F. Mixed and hybrid finite element methods; Vol. 2, New York: Springer-Verlag, 1991.
Elman, H.; Silvester, D.; Wathen, A. Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics; Oxford university press, 2014.
Layton, W. Introduction to the numerical analysis of incompressible viscous flows; SIAM, 2008.
Gunzburger, M. The inf-sup condition in mixed finite element methods with application to the Stokes system. Collected lectures on the preservation of stability under discretization 2002, pp. 93–121.
John, V.; Linke, A.; Merdon, C.; Neilan, M.; Rebholz, L.G. On the divergence constraint in mixed finite element methods for incompressible flows. SIAM review 2017, 59, 492–544. [Google Scholar] [CrossRef]
Vogelius, M. A right-inverse for the divergence operator in spaces of piecewise polynomials: Application to the p-version of the finite element method. Numerische Mathematik 1983, 41, 19–37. [Google Scholar] [CrossRef]
Qin, J. On the convergence of some low order mixed finite elements for incompressible fluids; The Pennsylvania State University, 1994.
Cockburn, B.; Kanschat, G.; Schötzau, D. A note on discontinuous Galerkin divergence-free solutions of the Navier–Stokes equations. Journal of Scientific Computing 2007, 31, 61–73. [Google Scholar] [CrossRef]
Evans, J.; Hughes, T. Isogeometric divergence-conforming B-splines for the unsteady Navier–Stokes equations. Journal of Computational Physics 2013, 241, 141–167. [Google Scholar] [CrossRef]
Falk, R.; Neilan, M. Stokes complexes and the construction of stable finite elements with pointwise mass conservation. SIAM Journal on Numerical Analysis 2013, 51, 1308–1326. [Google Scholar] [CrossRef]
Wang, J.; Ye, X. New finite element methods in computational fluid dynamics by H (div) elements. SIAM Journal on Numerical Analysis 2007, 45, 1269–1286. [Google Scholar] [CrossRef]
Zhang, S. A new family of stable mixed finite elements for the 3D Stokes equations. Mathematics of computation 2005, 74, 543–554. [Google Scholar] [CrossRef]
Linke, A. Collision in a cross-shaped domain: A steady 2d Navier–Stokes example demonstrating the importance of mass conservation in CFD. Computer Methods in Applied Mechanics and Engineering 2009, 198, 3278–3286. [Google Scholar] [CrossRef]
Olshanskii, M.; Reusken, A. Grad-div stabilization for Stokes equations. Mathematics of Computation 2004, 73, 1699–1718. [Google Scholar] [CrossRef]
Turek, S.; Ouazzi, A.; Hron, J. On pressure separation algorithms (PSepA) for improving the accuracy of incompressible flow simulations. International Journal for Numerical Methods in Fluids 2009, 59, 387–403. [Google Scholar] [CrossRef]
Linke, A. A divergence-free velocity reconstruction for incompressible flows. Comptes Rendus Mathematique 2012, 350, 837–840. [Google Scholar] [CrossRef]
Turek, S. On discrete projection methods for the incompressible Navier-Stokes equations: An algorithmical approach. Computer Methods in Applied Mechanics and Engineering 1997, 143, 271–288. [Google Scholar] [CrossRef]
Chorin, A.J. Numerical solution of the Navier-Stokes equations. Mathematics of computation 1968, 22, 745–762. [Google Scholar] [CrossRef]
Temam, R. Navier-Stokes equations: theory and numerical analysis; Vol. 343, American Mathematical Soc., 2001.
Guermond, J.L.; Minev, P.; Shen, J. An overview of projection methods for incompressible flows. Computer methods in applied mechanics and engineering 2006, 195, 6011–6045. [Google Scholar] [CrossRef]
Adams, R.A.; Fournier, J.J. Sobolev spaces; Elsevier, 2003.
Boffi, D.; Brezzi, F.; Fortin, M.; others. Mixed finite element methods and applications; Vol. 44, Springer, 2013.
Thomas, J.M. Sur l’analyse numérique des méthodes d’éléments finis hybrides et mixtes. PhD thesis, Université Pierre et Marie Curie, 1977.
Raviart, P.A.; Thomas, J.M. A mixed finite element method for second order elliptic problems. Mathematical Aspects of the Finite Element Method, Rome, 1975. [Google Scholar]
Durán, R. Mixed finite element methods. Mixed finite elements, compatibility conditions, and applications 2008, pp. 1–44.
Brezzi, F.; Douglas, J.; Durán, R.; Fortin, M. Mixed finite elements for second order elliptic problems in three variables. Numerische Mathematik 1987, 51, 237–250. [Google Scholar] [CrossRef]
Brezzi, F.; Fortin, M.; Marini, L.D.; others. Efficient rectangular mixed finite elements in two and three space variables. ESAIM: Mathematical Modelling and Numerical Analysis 1987, 21, 581–604. [Google Scholar] [CrossRef]
Brezzi, F.; Douglas, J.J.; Marini, D. Recent results on mixed finite element methods for second order elliptic problems. Vistas in applied mathematics. Numerical analysis, atmospheric sciences, immunology 1986, pp. 25–43.
Girault, V.; Raviart, P.A. Finite element methods for Navier-Stokes equations: theory and algorithms; Vol. 5, Springer Science & Business Media, 2012.
Capodaglio, G.; Aulisa, E. A particle tracking algorithm for parallel finite element applications. Computers & Fluids 2017, 159, 338–355. [Google Scholar] [CrossRef]
Aulisa, E.; Manservisi, S.; Scardovelli, R. A mixed markers and volume-of-fluid method for the reconstruction and advection of interfaces in two-phase and free-boundary flows. Journal of Computational Physics 2003, 188, 611–639. [Google Scholar] [CrossRef]
Aulisa, E.; Manservisi, S.; Scardovelli, R. A surface marker algorithm coupled to an area-preserving marker redistribution method for three-dimensional interface tracking. Journal of Computational Physics 2004, 197, 555–584. [Google Scholar] [CrossRef]
Lafaurie, B.; Nardone, C.; Scardovelli, R.; Zaleski, S.; Zanetti, G. Modelling merging and fragmentation in multiphase flows with SURFER. Journal of computational physics 1994, 113, 134–147. [Google Scholar] [CrossRef]
Harvie, D.; Fletcher, D. A new volume of fluid advection algorithm: the stream scheme. Journal of Computational Physics 2000, 162, 1–32. [Google Scholar] [CrossRef]
Leveque, R.J. High-resolution conservative algorithms for advection in incompressible flow. SIAM Journal on Numerical Analysis 1996, 33, 627–665. [Google Scholar] [CrossRef]
Aulisa, E.; Bná, S.; Bornia, G. FEMuS. https://github.com/eaulisa/MyFEMuS.git, 2014.
Cervone, A. ProXPDE. https://github.com/capitalaslash/proxpde.git, 2015.
Barbi, G.; Bornia, G.; Cerroni, D.; Cervone, A.; Chierici, A.; Chirco, L.; Viá, R.; Giovacchini, V.; Manservisi, S.; Scardovelli, R.; others. FEMuS-Platform: A numerical platform for multiscale and multiphysics code coupling. 9th International Conference on Computational Methods for Coupled Problems in Science and Engineering, COUPLED PROBLEMS 2021. International Center for Numerical Methods in Engineering, 2021, pp. 1–12.
Barbi, G.; Cervone, A.; Chierici, A.; Giovacchini, V.; Manservisi, S.; Scardovelli, R.; Sirotti, L. A New Projection Method for Navier-Stokes Equations by using Raviart-Thomas Finite Element. World Congress in Computational Mechanics and ECCOMAS Congress. Scipedia SL, 2022, pp. 1–12.

Figure 1. On the left a quadrilateral element

K_{k}

and on the right the square canonical element

\hat{K}

Figure 1. On the left a quadrilateral element

K_{k}

and on the right the square canonical element

\hat{K}

Figure 2. Single bubble test with

Q_{2}

velocity interpolation: comparison of the final interface position at

t = 4 s

(circular marker) with the initial circular geometry (dashed black line) for the meshes with

32 \times 32

(top left),

64 \times 64

(top right) and

128 \times 128

(bottom) cells.

Figure 2. Single bubble test with

Q_{2}

velocity interpolation: comparison of the final interface position at

t = 4 s

(circular marker) with the initial circular geometry (dashed black line) for the meshes with

32 \times 32

(top left),

64 \times 64

(top right) and

128 \times 128

(bottom) cells.

Figure 3. Single bubble velocity with

{RT}_{0}

velocity interpolation: comparison of the final interface position at

t = 4 s

(circular marker) with the initial circular geometry (dashed black line) for the meshes with

32 \times 32

(top left),

64 \times 64

(top right) and

128 \times 128

(bottom) cells.

Figure 3. Single bubble velocity with

{RT}_{0}

velocity interpolation: comparison of the final interface position at

t = 4 s

(circular marker) with the initial circular geometry (dashed black line) for the meshes with

32 \times 32

(top left),

64 \times 64

(top right) and

128 \times 128

(bottom) cells.

Figure 4. Single bubble test with

Q_{2}

velocity interpolation: interface position at maximum deformation (

t = 2 s

) for the meshes with

32 \times 32

(top left),

64 \times 64

(top right) and

128 \times 128

(bottom) cells.

Figure 4. Single bubble test with

Q_{2}

velocity interpolation: interface position at maximum deformation (

t = 2 s

) for the meshes with

32 \times 32

(top left),

64 \times 64

(top right) and

128 \times 128

(bottom) cells.

Figure 5. Single bubble test with

{RT}_{0}

velocity interpolation: interface position at maximum deformation (

t = 2 s

) for the meshes with

32 \times 32

(top left),

64 \times 64

(top right) and

128 \times 128

(bottom) cells.

Figure 5. Single bubble test with

{RT}_{0}

velocity interpolation: interface position at maximum deformation (

t = 2 s

) for the meshes with

32 \times 32

(top left),

64 \times 64

(top right) and

128 \times 128

(bottom) cells.

Figure 6. Two-dimensional channel flow: geometry (left), regular (center) and irregular coarsest mesh (right) for the channel flow.

Figure 7. Two-dimensional channel flow (regular mesh): streamwise velocity component v for coupled algorithm, split algorithm, and split algorithm with Raviart-Thomas approximation.

Figure 8. Three-dimensional channel flow geometry.

Figure 9. Three-dimensional channel flow: streamwise velocity component w for coupled algorithm, split algorithm, and split algorithm with Raviart-Thomas approximation.

Table 1. Values of the

E_{m}

E_{g}

and

E_{a l}

errors for different grids, with respective rates of convergence for the

Q_{2}

velocity interpolation.

Table 1. Values of the

E_{m}

E_{g}

and

E_{a l}

errors for different grids, with respective rates of convergence for the

Q_{2}

velocity interpolation.

$N_{e l}$	$E_{m}^{Q}$	$E_{g}^{Q}$	$E_{a l}^{Q}$	$p_{m}$	$p_{g}$	$p_{a l}$
$32 \times 32$	$1.64 \cdot 10^{- 2}$	$1.00 \cdot 10^{- 3}$	$9.66 \cdot 10^{- 4}$	-	-	-
$64 \times 64$	$9.88 \cdot 10^{- 4}$	$6.69 \cdot 10^{- 5}$	$8.52 \cdot 10^{- 5}$	$4.05$	$3.91$	$3.50$
$128 \times 128$	$2.35 \cdot 10^{- 4}$	$1.65 \cdot 10^{- 5}$	$1.36 \cdot 10^{- 5}$	$2.07$	$2.02$	$2.64$

Table 2. Values of the

E_{m}

E_{g}

and

E_{a l}

errors for different grids, with respective rates of convergence for

{RT}_{0}

velocity interpolation.

Table 2. Values of the

E_{m}

E_{g}

and

E_{a l}

errors for different grids, with respective rates of convergence for

{RT}_{0}

velocity interpolation.

$N_{e l}$	$E_{m}^{R T}$	$E_{g}^{R T}$	$E_{a l}^{R T}$	$p_{m}$	$p_{g}$	$p_{a l}$
$32 \times 32$	$6.43 \cdot 10^{- 2}$	$3.94 \cdot 10^{- 3}$	$2.24 \cdot 10^{- 3}$	-	-	-
$64 \times 64$	$1.57 \cdot 10^{- 2}$	$1.07 \cdot 10^{- 3}$	$5.79 \cdot 10^{- 4}$	$2.03$	$1.89$	$1.95$
$128 \times 128$	$3.73 \cdot 10^{- 3}$	$2.59 \cdot 10^{- 4}$	$1.42 \cdot 10^{- 4}$	$2.08$	$2.04$	$2.03$

Table 3. Two-dimensional channel test:

L^{2}

-norm of the velocity field for different levels l of grid refinement with

n_{e l}

number of elements, for regular and irregular mesh.

Table 3. Two-dimensional channel test:

L^{2}

-norm of the velocity field for different levels l of grid refinement with

n_{e l}

number of elements, for regular and irregular mesh.

l	$n_{e l}$	Regular mesh			Irregular mesh
l	$n_{e l}$	$P^{c}$	$P^{s}$	${RT}_{0}$	$P^{c}$	$P^{s}$	${RT}_{0}$
1	6.40E+01	1.551	1.408	1.515	1.516	1.408	1.526
2	2.56E+02	1.544	1.470	1.530	1.525	1.470	1.531
3	1.02E+03	1.539	1.502	1.532	1.530	1.501	1.532
4	4.10E+03	1.536	1.517	1.532	1.532	1.517	1.532
5	1.64E+04	1.535	1.525	1.532	1.533	1.525	1.531
6	6.55E+04	1.534	1.529	1.530	1.533	1.529	1.528

Table 4. Two-dimensional cavity test with projection method: velocity error norm and convergence rate for different levels l of grid refinement and corresponding number of elements

n_{e l}

, for

{RT}_{0}

finite element approximation.

Table 4. Two-dimensional cavity test with projection method: velocity error norm and convergence rate for different levels l of grid refinement and corresponding number of elements

n_{e l}

, for

{RT}_{0}

finite element approximation.

l	$n_{e l}$	${RT}_{0}$
l	$n_{e l}$	$∥ u^{*} - u_{h} ∥_{0}$	$\frac{ε_{l - 1}}{ε_{l}}$	p
1	6.40E+01	9.59E-01	-	-
2	2.56E+02	4.86E-01	1.976	0.997
3	1.02E+03	2.44E-01	1.994	0.997
4	4.10E+03	1.22E-01	1.998	1.002
5	1.64E+04	6.10E-02	2.000	0.998
6	6.55E+04	3.05E-02	2.000	0.999
7	2.62E+05	1.52E-02	1.999	1.000

Table 5. Two-dimensional cavity test with Taylor-Hood finite element approximation: velocity error norm for different levels l of grid refinement and corresponding number of elements

n_{e l}

, for coupled and split algorithm.

Table 5. Two-dimensional cavity test with Taylor-Hood finite element approximation: velocity error norm for different levels l of grid refinement and corresponding number of elements

n_{e l}

, for coupled and split algorithm.

l	$n_{e l}$	Coupled $P_{2} - P_{1}$		Split $P_{2} - P_{1}$
l	$n_{e l}$	$∥ u^{*} - u_{h} ∥_{0}$	$\frac{ε_{l - 1}}{ε_{l}}$	$∥ u^{*} - u_{h} ∥_{0}$	$\frac{ε_{l - 1}}{ε_{l}}$
1	6.40E+01	4.60E-02	-	4.56E-02	-
2	2.56E+02	5.81E-03	7.925	5.77E-03	7.905
3	1.02E+03	7.25E-04	8.018	7.23E-04	7.982
4	4.10E+03	9.05E-05	8.007	9.06E-05	7.986
5	1.64E+04	1.13E-05	8.002	1.13E-05	8.011
6	6.55E+04	1.43E-06	7.929	1.41E-06	7.995

Table 6. Two-dimensional cavity test with Taylor-Hood finite element approximation: velocity divergence error norm for different levels l of grid refinement and corresponding number of elements

n_{e l}

, for coupled and split algorithm.

Table 6. Two-dimensional cavity test with Taylor-Hood finite element approximation: velocity divergence error norm for different levels l of grid refinement and corresponding number of elements

n_{e l}

, for coupled and split algorithm.

l	$n_{e l}$	Coupled $P_{2} - P_{1}$		Split $P_{2} - P_{1}$
l	$n_{e l}$	$∥ \nabla \cdot u^{*} - \nabla \cdot u_{h} ∥_{0}$	$\frac{ε_{l - 1}}{ε_{l}}$	$∥ \nabla \cdot u^{*} - \nabla \cdot u_{h} ∥_{0}$	$\frac{ε_{l - 1}}{ε_{l}}$
1	6.40E+01	4.33E-01	-	4.42E-01	-
2	2.56E+02	1.13E-01	3.839	1.13E-01	3.92
3	1.02E+03	2.83E-02	3.978	2.83E-02	3.98
4	4.10E+03	7.09E-03	3.996	7.09E-03	4.00
5	1.64E+04	1.77E-03	3.999	1.77E-03	4.00
6	6.55E+04	4.44E-04	3.990	4.43E-04	4.00

Table 7. Three-dimensional channel test: velocity error norm for different levels l of grid refinement and the corresponding number of elements

n_{e l}

Table 7. Three-dimensional channel test: velocity error norm for different levels l of grid refinement and the corresponding number of elements

n_{e l}

l	$n_{e l}$	$P^{c}$	$P^{s}$	${RT}_{0}$
1	2.08E+02	2.147	1.962	2.301
2	4.10E+03	2.236	2.139	2.312
3	3.28E+04	2.280	2.231	2.323

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer

Numerical Algorithms for Divergence-Free Velocity Applications

Abstract

1. Introduction

2. Mathematical Framework

2.1. Notation

2.2. The Mixed-Method-Robustness Property

3. H ( div , Ω ) Function Space Approximations

3.1. Partition of Domain and H ( div , Ω ) Functions

3.2. Piola’s Transformation

3.3. Lowest-order Raviart-Thomas finite element approximation

4. Divergence-Free Field for Two-Phase Flow

4.1. Problem Model

4.2. Approximation over Cartesian Cells

4.3. Numerical Tests

5. Projection Method for Navier-Stokes Equation

5.1. Model and Algorithm

5.2. Numerical Tests

6. Conclusions

Funding

Conflicts of Interest

References

MDPI Initiatives

Important Links

Subscribe

3. $H (div, Ω)$ Function Space Approximations

3.1. Partition of Domain and $H (div, Ω)$ Functions