arXiv:1712.02520v2 [cs.CE] 25 Mar 2018

Higher-order surface FEM for incompressible
Navier-Stokes flows on manifolds
T.P. Fries
March 28, 2018

Institute of Structural Analysis
Graz University of Technology
Lessingstr. 25/II, 8010 Graz, Austria
www.ifb.tugraz.at
fries@tugraz.at

Abstract
Stationary and instationary Stokes and Navier-Stokes flows are considered on
two-dimensional manifolds, i.e., on curved surfaces in three dimensions. The higherorder surface FEM is used for the approximation of the geometry, velocities, pressure,
and Lagrange multiplier to enforce tangential velocities. Individual element orders
are employed for these various fields. Stream-line upwind stabilization is employed
for flows at high Reynolds numbers. Applications are presented which extend classical benchmark test cases from flat domains to general manifolds. Highly accurate
solutions are obtained and higher-order convergence rates are confirmed.
Keywords: Stokes, Navier-Stokes, higher-order FEM, surface FEM, surface PDEs,
manifold

2

CONTENTS

Contents
1 Introduction

4

2 Preliminaries

6

2.1

Surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

2.2

Surface operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

2.2.1

The tangential projector . . . . . . . . . . . . . . . . . . . . . . . .

7

2.2.2

Surface gradient of scalar quantities . . . . . . . . . . . . . . . . . .

7

2.2.3

Surface gradient of vector quantities

. . . . . . . . . . . . . . . . .

8

2.2.4

Divergence operators and divergence theorem . . . . . . . . . . . .

8

3 Governing equations
3.1

3.2

3.3

9

Flow models in strong form . . . . . . . . . . . . . . . . . . . . . . . . . .

9

3.1.1

Stationary Stokes flow . . . . . . . . . . . . . . . . . . . . . . . . .

9

3.1.2

Stationary Navier-Stokes flow . . . . . . . . . . . . . . . . . . . . .

11

3.1.3

Instationary Navier-Stokes flow . . . . . . . . . . . . . . . . . . . .

11

Flow models in weak form . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

3.2.1

Stationary Stokes flow . . . . . . . . . . . . . . . . . . . . . . . . .

12

3.2.2

Stationary Navier-Stokes flow . . . . . . . . . . . . . . . . . . . . .

13

3.2.3

Instationary Navier-Stokes flow . . . . . . . . . . . . . . . . . . . .

13

Surface FEM for flows on manifolds . . . . . . . . . . . . . . . . . . . . . .

14

3.3.1

Surface meshes . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14

3.3.2

Surface FEM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

3.3.3

Stationary Stokes flow . . . . . . . . . . . . . . . . . . . . . . . . .

16

3.3.4

Stationary Navier-Stokes flow . . . . . . . . . . . . . . . . . . . . .

17

3.3.5

Instationary Navier-Stokes flow . . . . . . . . . . . . . . . . . . . .

18

CONTENTS

4 Numerical results

3
19

4.1

Stokes flow on an axisymmetric surface . . . . . . . . . . . . . . . . . . . .

20

4.2

Driven cavity flows on manifolds . . . . . . . . . . . . . . . . . . . . . . . .

25

4.3

Flows on zero-level sets . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

4.4

Cylinder flows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

5 Conclusions

36

6 Acknowledgements

39

References

39

4

1

Introduction

Introduction

The solution of boundary value problems on curved surfaces has many practical applications in mathematics, physics, and engineering. For example, there are transport processes
on interfaces, e.g., in foams, biomembranes and bubble surfaces [16, 25, 45], or structurerelated phenomena such as in membranes and shells [2, 7]. Herein, Stokes and incompressible Navier-Stokes flows on curved, two-dimensional manifolds are considered. The
governing equations for flows on moving surfaces are discussed in [3, 33] based on fundamental surface continuum mechanics and conservation laws and in [35], an energetic approach
is presented. Earlier works in a similar context may be traced back to [15, 26, 43, 47].
For an excellent overview, the reader is refered to [33]. The references given above often
focus on mathematical properties such as the existence and uniqueness of the solutions or
stabilitity analyses. Applications are often two-phase flows where the fluid field in the bulk
and on the moving interface are coupled. However, it is also worthwhile to consider the
situation for fixed manifolds, e.g., related to meterology and oceanography where the flows
take place on (part of) a sphere. Special geometries such as hyperbolic planes and spheres
are discussed in [6, 34, 36].
Herein, the focus is on the approximation of stationary and instationary (Navier-)Stokes
flows on fixed manifolds based on the surface finite element method as outlined in [11, 13,
14]. The governing equations resemble the three-dimensional (Navier-)Stokes equations
where the classical gradient and divergence operators are replaced by their tangential counterparts derived from tangential differential calculus [33]. The equations are formulated in
the classical stress-divergence form, contrasted to the approach in [37]. An additional constraint is required to enforce that the velocities remain in the tangent space of the manifold;
it is labelled “tangential velocity constraint”. The models are first given in strong form and
are then transformed to the weak form to enable a numerical solution based on the surface
FEM. Finite element spaces of different orders are employed for the approximation of the
geometry and of the involved physical fields, i.e., the velocities, pressure and the Lagrange
multiplier field required to enforce the tangential velocity constraint. It is found that the
balance of these element orders is critical for the accuracy and conditioning of the system
of equations. In particular, the well-known Babuška-Brezzi condition applies [1, 4, 17] as
both, the incompressibility constraint and the tangential velocity constraint are enforced
using Lagrange multipliers. For the case of the instationary Navier-Stokes equations, the
Crank-Nicolson time stepping scheme is employed for the semi-discrete sytem of equations

5
resulting from using the surface FEM in space. Surface FEM based on linear elements is
used in the recent work [41], where the penalty method is employed to enforce tangential velocities and a projection method rather than a monolithic approach is suggested to
solve for the different physical fields. Alternatives for the surface FEM are the TraceFEM
[8, 22, 40] and CutFEM [27, 28], where the basis functions are generated from a background
mesh in the bulk surrounding the manifold of interest.
Using the FEM for the Navier-Stokes flows at large Reynolds numbers requires stabilization. Herein, the streamline-upwind Petrov-Galerkin (SUPG) approach is used [5, 49].
Alternatively, other variants such as the Galerkin least squares stabilization [32] and variational multiscale approaches [29, 23] may also be employed. Stabilization for advectiondiffusion applications on manifolds are considered in [38].
The numerical results show that higher-order convergence rates are achieved provided that
the finite element spaces are properly chosen. Also the conditioning of the system of
equations depends on the element orders employed for the approximation of the individual
physical fields. The presented results are based on well-known benchmark test cases in
two dimensions such as driven cavity flows and cylinder flows with vertex shedding which,
herein, are extended to curved surfaces. Due to the higher-order elements, the results are
highly accurate and may serve as future benchmarks in the context of (Navier-)Stokes flows
on manifolds. Most test cases are carried out on parametrized surfaces, however, also the
situation of flows on zero-isosurfaces is covered herein.
To the best of our knowledge, this is the first time, where (i) general higher-order surface FEM is used for the (in)stationary (Navier-)Stokes equations on manifolds including stabilization, (ii) numerical convergence studies are presented confirming higher-order
convergence rates, and (iii) benchmark test cases are proposed and solutions presented.
Furthermore the notation employed is closely related to the typical engineering literature
and aims to provide a bridge from the mathematical to the engineering community.
The paper is organized as follows: In Section 2, some requirements and properties of
surfaces are described and tangential differential operators are defined based on [10, 14].
Section 3 covers the governing equations for (i) Stokes flow, (ii) stationary, and (iii) instationary Navier-Stokes flows on two-dimensional manifolds. They are given in strong form,
weak form, and discretized weak form according to the surface FEM. Numerical results
are presented in Section 4. Convergence studies are performed for a test case for which an
analytic solution is available and it is shown that higher-order convergence rates can be

6

Preliminaries

(a)

(b)

(c)

(d)

Figure 1: Some examples of (a, b) compact manifolds with ∂Γ = ∅ and (c, d) manifolds
with boundary ∂Γ.

achieved. For the other test cases where no analytic solutions are available it is confirmed
that in the flat two-dimensional case, well-known reference solutions are reproduced. Various meshes with different orders and resolutions have been employed to obtain highly
accurate results on curved surfaces. Finally, a summary and outlook are given in Section
5.

2
2.1

Preliminaries
Surfaces

The task is to solve a boundary value problem (BVP) on an arbitrary surface Γ in three
dimensions. Let the surface be fixed in space over time, possibly curved, sufficiently smooth,
orientable, connected (so that there is only one surface), and feature a finite area. There
is a unit normal vector nΓ ∈ R3 on Γ. The surface may be compact, i.e., without a
boundary, ∂Γ = ∅, see Figs. 1(a) and (b) for examples. Otherwise, it may be bounded by
∂Γ as shown in Figs. 1(c) and (d). Then, associated with ∂Γ, there is a tangential vector
t∂Γ pointing in direction of ∂Γ and a co-normal vector n∂Γ = nΓ × t∂Γ pointing “outwards”
and being normal to ∂Γ and tangent to Γ. The surface may be given in parametrized form
or implied, e.g., based on the level-set method; both situations are considered herein. For
the equivalence of these two cases and more mathematical details, see, e.g., [14].

2.2

7

Surface operators

2.2
2.2.1

Surface operators
The tangential projector

On the manifold Γ, the tangential projector P (x) ∈ R3×3 is defined by the normal vector
as
P (x) = I − nΓ (x) ⊗ nΓ (x) .
Some important properties are: (i) P · nΓ = 0, (ii) P = PT , and (iii) P · P = P.

2.2.2

Surface gradient of scalar quantities

The tangential gradient operator ∇Γ of a differentiable scalar function u : Γ → R on the
manifold is given by
∇Γ u (x) = P (x) · ∇ũ (x) ,

x ∈ Γ,

(2.1)

where ∇ is the standard gradient operator, and ũ is a smooth extension of u in a neighborhood U of the manifold Γ. Of course, ũ may also be some given function (rather
than an arbitrary extension) in global coordinates, i.e., ũ (x) : R3 → R. For the case of
parametrized surfaces defined by the map x (r) : R2 → R3 , and a given scalar function
u (r) : R2 → R, the tangential gradient may be determined without explicitly computing
an extension ũ using
∇Γ u (x (r)) = J (r) · G−1 (r) · ∇r u (r) ,
(2.2)
with J = ∂x/∂r being the (3 × 2)-Jacobi matrix and G = JT · J being the metric tensor
(first fundamental form). Equation (2.2) shall be used later in the context of the FEM
to determine tangential gradients of shape functions. It is noteworthy that ∇Γ u is in the
tangent space of Γ and, thus, P · ∇Γ u = ∇Γ u and ∇Γ u · nΓ = 0. The components of the
tangential gradient are denoted by


T

∇Γ u (x) = ∂xΓ u, ∂yΓ u, ∂zΓ u

,

representing the first-order partial derivatives on Γ. Second-order partial derivatives may
be denoted by


Heij (u (x)) = ∂xΓi xj u (x) = ∂xΓi ∂xΓj u (x) ,

8

Preliminaries

where Heij (u (x)) is the tangential Hessian matrix. In the context of manifolds, this matrix
is not symmetric [10], i.e., for mixed second derivatives ∂xΓi xj u 6= ∂xΓj xi u for i 6= j.
2.2.3

Surface gradient of vector quantities

Next, operators for vector quantitites u (x) : Γ → R3 are considered. The “directional
gradient” of u is the tensor of tangential derivatives and defined as





dir
∇dir
Γ u (x) = ∇Γ 

u (x)
v (x)
w (x)










=





∂xΓ u ∂yΓ u ∂zΓ u
∂xΓ v ∂yΓ v ∂zΓ v
∂xΓ w ∂yΓ w ∂zΓ w







= ∇ũ · P.

In contrast, the covariant derivatives are
dir
∇cov
Γ u (x) = P · ∇Γ u (x) = P · ∇ũ · P.

One has to carefully distinguish these two different gradient operators. It is noted that
∇cov
Γ u appears frequently in the modeling of physical phenomena on manifolds, i.e., in the
governing equations. On the other hand, ∇dir
Γ u is often used in straightforward extensions
of identities such as product rules and divergence theorems. For example, we have for a
scalar function f (x) and vector functions u (x), v (x)
dir
∇dir
Γ (f · u) = ∇Γ f ⊗ u + f · ∇Γ u,
T
dir
T
v T · ∇dir
Γ u = −u · ∇Γ v + ∇Γ (u · v) ,

however, the relations are less straightforward for the covariant counterparts ∇cov
Γ (f · u)
T
cov
and v · ∇Γ u, respectively. Later on, in the context of FEM implementations, it proves
useful to transform covariant derivatives systematically to directional ones. This allows
the computation of directional derivatives of FE shape functions with respect to x ∈ R3
independent of the integration of the weak form of the governing equations.

2.2.4

Divergence operators and divergence theorem

The divergence of a vector function u (x) : Γ → R3 is given as




cov
divΓ u (x) = tr ∇dir
Γ u = tr (∇Γ u) =: ∇Γ · u.

9
For a tensor function A (x) : Γ → R3×3 , there holds





divΓ A (x) = 

divΓ (A11 , A12 , A13 )
divΓ (A21 , A22 , A23 )
divΓ (A31 , A32 , A33 )







=: ∇Γ · A.

It may be shown that divΓ P = −κ · nΓ with κ = tr (H) being the mean curvature and
H = ∇cov
Γ nΓ being the second fundamental form.
The following divergence theorem on manifolds is later needed for deriving the weak forms
[9, 10],
Z
Γ

u · divΓ A dA = −

Z
Γ

∇dir
Γ u

: A dA +



Z
Γ

κ · u · A · nΓ dA +

Z
∂Γ

u · A · n∂Γ ds,

(2.3)



dir
T
where ∇dir
. For tangential tensor functions with A = P · A · P,
Γ u : A = tr ∇Γ u · A
the term involving the curvature κ vanishes because then A · nΓ = 0. In this case, we also
cov
have ∇dir
Γ u : A = ∇Γ u : A.

3

Governing equations

In the following, we consider (i) stationary Stokes flow, (ii) stationary Navier-Stokes flow,
and (iii) instationary Navier-Stokes flow on fixed manifolds. The governing equations are
first given in strong and weak forms. The surface FEM is then applied for the discretization
of the weak forms. As mentioned above, these models are also considered, e.g., in [3, 33, 35]
among others.

3.1
3.1.1

Flow models in strong form
Stationary Stokes flow

Starting point is stationary Stokes flow on a manifold. Let u (x) ∈ C 2 (Γ) be the threedimensional velocity field on the surface Γ, p (x) ∈ C 1 (Γ) a pressure field, and ft (x)
a tangential body force, e.g., with unit N/m2 . The governing field equations (in stress-

10

Governing equations

divergence-form [12]) to be fulfilled ∀x ∈ Γ are
−P · divΓ σ (u, p) = ft ,

(3.1)

divΓ u = 0,

(3.2)

u · nΓ = 0.

(3.3)

Equation (3.1) expands to three momentum equations, equation (3.2) is the incompressibility constraint and equation (3.3) represents the tangential velocity constraint that restricts
the velocities to the tangent space of Γ. Two different strain tensors are introduced,

T
1
dir
· ∇dir
u
+
∇
u
,
ε (u) =
Γ
Γ
2

1  cov
T
εcov (u) =
· ∇Γ u + (∇cov
u)
,
Γ
2
dir





(3.4)
(3.5)

which are related to each other as εcov (u) = P · εdir (u) · P. The stress tensor is then
defined as
σ (u, p) = −p · P + 2µ · εcov (u)
where µ ∈ R+ is the (constant) dynamic viscosity. It is easily shown that
−P · divΓ σ (u, p) = ∇Γ p − 2µP · divΓ εcov (u)
Suppose there exists a boundary ∂Γ of the manifold that consists of two non-overlapping
parts, the Dirichlet boundary, ∂ΓD , and the Neumann boundary, ∂ΓN . The corresponding
boundary conditions are given as
u (x)
= û (x)
σ (x) · n∂Γ (x) = t̂ (x)

on ∂ΓD ,
on ∂ΓN ,

(3.6)

where the prescribed velocities û and tractions t̂ are in the tangent space of Γ, i.e., û·nΓ =
t̂ · nΓ = 0.
Note that, in general, there are no explicit boundary conditions needed for the pressure p.
In cases where no Neumann boundary is present, i.e., ∂ΓN = ∅ and ∂ΓD = ∂Γ, the pressure
is defined up to a constant [12, 24]. This includes compact manifolds where ∂Γ = ∅. In
such situations, the pressure may be prescribed at a given point on Γ or it is imposed by
R
a constraint in the form of Γ p dA = 0.

3.1

11

Flow models in strong form

Vorticity on manifolds. The vorticity ω is a physical quantity frequently computed in
flow problems. In the context of manifolds, we shall define
ω = ∇cov
Γ × u.

(3.7)

Note that ω is co-linear to the normal vector nΓ , hence, P · ω = 0. Therefore, it is useful
to determine the signed magnitude of ω, that is, the scalar function
ω ? (x) = ω · nΓ = ± kωk

∀x ∈ Γ.

(3.8)




This scalar quantity may also be obtained using directional derivatives, i.e., ω ? = ∇dir
Γ ×u ·
nΓ .

3.1.2

Stationary Navier-Stokes flow

For stationary Navier-Stokes flow, a non-linear advection term is added to equation (3.1)
resulting into
% · (u · ∇cov
(3.9)
Γ ) u − P · divΓ σ (x) = ft (x) ,
cov
where % ∈ R+ is the (constant) fluid density with unit kg/m2 and (u · ∇cov
Γ ) u := (∇Γ u) · u.
It is quite common to express the body force in the form ft (x) = % · g t (x) where g t
may consider gravity as g t = P · [0, 0, −9.81]T m/s2 for instance. The remaining equations
(3.2) and (3.3) and the boundary conditions (3.6) remain unchanged. The solution of the
non-linear governing equations can be obtained iteratively based on the Newton-Raphson
method or other fixed-point iterations such as Picard iterations. Because the advection operator is not self-adjoint, well-known stability issues may arise for large Reynolds numbers
in a numerical context.

3.1.3

Instationary Navier-Stokes flow

For instationary Navier-Stokes flow, the momentum equation (3.1) changes to
% · (∂t u (x, t) + (u · ∇cov
Γ ) u − g t (x, t)) − P · divΓ σ (x, t) = 0.

(3.10)

The functions representing the physical fields live in space (on Γ) and time, i.e., in the
time interval τ = [0, T ]. Therefore, Eqs. (3.10), (3.2), and (3.3) have to be solved in the

12

Governing equations

space-time domain Γ × τ . Herein, we restrict ourselves to spatially fixed manifolds Γ.
The boundary conditions (3.6) also extend in time dimension, hence, there are prescribed
velocities û (x, t) along ∂ΓD × τ and tractions t̂ (x, t) along ∂ΓN × τ . Furthermore, an
initial condition is needed,
u (x, 0) = u0 (x) , with divΓ u0 = 0 and u0 · nΓ = 0 ∀x ∈ Γ at t = 0.

3.2

(3.11)

Flow models in weak form

The following trial and test function spaces are introduced,
Su =

n

Vu =

n

o

u ∈ H1 (Γ)3 , u = û on ∂ΓD ,

(3.12)
o

wu ∈ H1 (Γ)3 , wu = 0 on ∂ΓD ,

(3.13)

Sp = Vp = L2 (Γ) ,

(3.14)

Sλ = Vλ = L2 (Γ) .

(3.15)

As mentioned previously, if no Neumann boundary exists, i.e., ∂ΓN = ∅, the pressure is
defined up to a constant and one may replace Sp by
n

Sp0 = p ∈ L2 (Γ) ,
3.2.1

Z

o

p dA = 0 .

Γ

Stationary Stokes flow

The weak form of the Stokes problem becomes: Given viscosity µ ∈ R+ , body force f (x)
in Γ, and traction t̂ (x) on ∂ΓN , find the velocity field u (x) ∈ Su , pressure field p (x) ∈ Sp ,
and Lagrange multiplier field λ (x) ∈ Sλ such that for all test functions (wu , wp , wλ ) ∈
Vu × Vp × Vλ , there holds in Γ
Z
Γ

∇dir
Γ w u : σ (u, p) dA +

Z
Γ

λ · (wu · nΓ ) dA =

Z

Γ

wu · f dA +

Z

wu · t̂ ds, (3.16)

∂ΓN

wp · divΓ u dA = 0,

(3.17)

wλ · (u · nΓ ) dA = 0.

(3.18)

Z Γ
Γ

Z

In order to obtain Eq. (3.16), the divergence theorem (2.3) was applied to − Γ wu ·
divΓ σ dA where the curvature term vanishes due to σ · nΓ = 0. Using the definition
R

3.2

13

Flow models in weak form

of the stress tensor, we get
Z
Γ

∇dir
Γ wu

: σ (u, p) dA = −

Z
Γ

∇dir
Γ wu

Z

: (p · P) dA + 2µ ·

Γ

cov
∇dir
(u) dA
Γ wu : ε

The following relations are easily derived:




dir
∇dir
Γ w u : (p · P) = p · tr ∇Γ w u · P

= p · divΓ wu




cov
cov
(u)
(u) = tr ∇dir
∇dir
Γ wu · ε
Γ wu : ε

= tr (εcov (wu ) · εcov (u))


dir
= tr P · ∇dir
(u) · P
Γ wu · ε



(3.19)

It is readily verified that solutions of the strong form also fulfill the weak form from above.
This is obvious for Eqs. (3.17) and (3.18) due to Eqs. (3.2) and (3.3), respectively. For
the momentum equations, it is noted that (3.16) is fufilled for −divΓ σ (u, p) + λ · nΓ = f .
Restricting this to the tangential space by multiplication with the projector P yields the
strong form of the momentum equations (3.1) because P · nΓ = 0. It is thus also seen
that the Lagrange multiplier field λ may be physically interpreted as a force in normal
direction.

3.2.2

Stationary Navier-Stokes flow

The weak form of the stationary Navier-Stokes equations is similar to the Stokes problem
from above, however, Eq. (3.16) is replaced by
%·

Z
Γ

wu · (u · ∇cov
Γ ) u dA +

Z
Γ

∇dir
Γ w u : σ (u, p) dA +
=

Z
Γ

Z
Γ

λ · (wu · nΓ ) dA

wu · f dA +

Z

wu · t̂ ds

∂ΓN

where the added advection term is readily identified.

3.2.3

Instationary Navier-Stokes flow

The weak form of the instationary Navier-Stokes problem is: Given density % ∈ R+ , viscosity µ ∈ R+ , body force % · g (x, t) in Γ × τ , traction t̂ (x, t) on ∂ΓN × τ , and initial condition

14

Governing equations

u0 (x) on Γ at t = 0 according to (3.11), find the velocity field u (x, t) ∈ L2 (τ ; Su ), pressure
field p (x, t) ∈ L2 (τ ; Sp ), and Lagrange multiplier field λ (x, t) ∈ L2 (τ ; Sλ ) such that for
all test functions (wu , wp , wλ ) ∈ Vu × Vp × Vλ , there holds in Γ × τ
%·

Z
Γ

+

Z
Γ

wu · (∂t u + (u · ∇cov
Γ ) u − g) dA

∇Γ wdir
u : σ (u, p) dA +

Z
Γ

λ · (wu · nΓ ) dA =

Z

3.3
3.3.1

wu · t̂ ds,

∂ΓN

wp · divΓ u dA = 0,

(3.21)

wλ · (u · nΓ ) dA = 0.

(3.22)

Z Γ
Γ

(3.20)
Z

Surface FEM for flows on manifolds
Surface meshes

Assume that a suitable surface mesh composed by higher-order triangular or quadrilateral
Lagrange elements of order q may be generated with desired element sizes and all nodes
on Γ. Well-known, necessary requirements of meshes such as the shape regularity of the
elements and bounds on inner angles, are fulfilled. The shape of each (physical) element
in the mesh results from a map of the corresponding reference element with nq nodes,


x (r) =






x (r, s)
y (r, s)
z (r, s)







=

nq
X

Niq (r) xi .

(3.23)

i=1

Niq (r) are classical Lagrangean shape functions of order q in reference coordinates r ∈ R2
and xi ∈ Γ are the nodal coordinates. The resulting mesh is an approximation Γhq ∈ C 0 of
the exact surface Γ. Clearly, Γhq is defined parametrically through the map (3.23) even if
the original Γ was implicitly given, e.g., by the zero-isosurface of a level-set function. See
[18, 19, 20] for the automatic generation of higher-order meshes on zero-isosurfaces. The
discrete unit normal vector is
∂r x × ∂s x
nhΓ =
k∂r x × ∂s xk
and is not smooth across element edges due to the C 0 -continuity of the surface mesh. The
discrete tangent and co-normal vectors th∂Γ and nh∂Γ are easily obtained along the element
edges on the boundary of Γhq . The definitions of the surface operators from Section 2.2

3.3

15

Surface FEM for flows on manifolds

readily extend to the case of a discrete manifold Γhq and are not repeated here.

3.3.2

Surface FEM

We use higher-order surface FEM as detailed, e.g., in [11, 14] for the discretization of the
weak forms from above. Finite element spaces of different orders are involved. As mentioned before, suitable surface meshes of order q may be generated defining approximations
Γhq ∈ C 0 of Γ. Let there be a “geometry mesh” of order q = kgeom with the sole purpose to
approximate the geometry of the manifold Γh = Γhkgeom and define the element maps (3.23).
In particular, this mesh is not used to imply a finite element space for the approximation
of the weak forms.
Next, a finite element space of order k is generated on Γh for which it is assumed that there
is a second mesh of order k. The two meshes feature the same element types and number
of elements with identical coordinates at the corners, however, the total number of nodes
differs due to the individual orders. It is emphasized that the coordinates of the nodes
in the k-th order mesh are, in fact, never needed and it is only the connectivity which is
required to set up the finite element space.
Associated to triangular or quadrilateral elements in the k-th order mesh, there is a fixed
n
o
set of local basis functions Nik (r) defined in a reference element with i = 1, . . . , nk
and nk being the number of nodes per element. Classical Lagrange basis functions with
Nik (r j ) = δij are used herein. Based on the map (3.23) which is completely determined
n
o
by the geometry mesh, one may generate Nik (x (r)) for all x ∈ Γh and tangential
derivatives ∇Γ Nik (x (r)) are determined based on Eq. (2.2). This is only an iso-parametric
map when k = kgeom . Summing up the element contributions for nodes belonging to several
n
o
elements, this generates a set of global, C 0 -continuous basis functions Mik (x (r)) in Γh
with i = 1, . . . , nknodes and nknodes being the number of nodes of the k-th order surface mesh.
n
o
Note that to generate the nodal basis Mik (x (r)) , only the coordinates of the geometry
mesh are needed, however, not from the k-th order mesh. A general finite element space
of order k is now defined by
Qhk



h

= u ∈ C0



Γhkgeom



nknodes
h

,u =

X
i=1

Mik







(x (r)) · ui , ui ∈ R ⊂ H1 Γhkgeom .

16

Governing equations

Based on this, the following discrete trial and test function spaces are defined,
Suh



=


Vuh =

h

u ∈

i3
Qhku

h

h

whu ∈ Qhku

h

h

, u = û on

i3

∂ΓhD



,

(3.24)

, whu = 0 on ∂ΓhD ,

(3.25)



Sph = Vph = Qhkp ,

(3.26)

Sλh = Vλh = Qhkλ .

(3.27)

Although shape functions for the pressure and the Lagrange multiplier for enforcing the
tangential velocity constraint may be discontinuous, we restrict ourselves to classical C 0 continuous approximations. Note that individual orders ku , kp , and kλ are associated to the
approximations of velocities uh , pressure ph , and Lagrange multiplier field λh , respectively.
Analogous to the continuous case, one may impose that the functions in Sph have to fulfill
R h
Γ p dA = 0 if no Neumann boundary is present.

3.3.3

Stationary Stokes flow

The discrete weak form of the Stokes problem reads: Given viscosity µ ∈ R+ , body force
h
f h (x) in Γh , and traction t̂ (x) on ∂ΓhN , find the velocity field uh (x) ∈ Suh , pressure field
ph (x) ∈ Sph , and Lagrange multiplier field λh (x) ∈ Sλh such that for all test functions


whu , wph , wλh ∈ Vuh × Vph × Vλh , there holds in Γh
Z
Γ





h
h h
∇dir
dA +
Γ wu : σ u , p

Z



Γ

Z
Γ

Z
Γ



λh · whu · nhΓ dA =

Z
Γ

whu · f h dA +

Z

wph · divΓ uh dA = 0,


whu · tˆh ds,
(3.28)

∂ΓN

(3.29)



wλh · uh · nhΓ dA = 0.

(3.30)

The usual element assembly yields a linear system of equations in the form






 





K G L   u   f
 


· p = 0
GT 0 0 
 


LT 0 0
λ
0




,


(3.31)

with [u, p, λ]T = [u, v, w, p, λ]T being the sought nodal values of the velocity components,
pressure, and Lagrange multiplier. For the implementation, it is interesting to compare the

3.3

17

Surface FEM for flows on manifolds

system (3.31) with the system obtained for a classical three-dimensional Stokes problem,


 







K3D G3D   u   f 

·
=
.
GT
0
p
0
3D

(3.32)

Assume a function which generates K3D and G3D based on three-dimensional FE shape
functions (including classical partial derivatives with respect to x) evaluated at given integration points in 3D. The same function may be used for generating K and G provided that (i) the integration points are restricted to Γh with proper weights, (ii) the
classical partial derivatives in ∇ are replaced by the tangential derivatives as in ∇dir
Γ ,
and (iii) the contribution to K at the current integration point, K (xi ), is projected as
K (xi ) = P (xi ) · K3D (xi ) · P (xi ) which is due to Eq. (3.19). The same shall later hold for
the advection matrix C (u) in the Navier-Stokes equations.
As expected in the context of the Lagrange multiplier method, the matrix in Eq. (3.31)
has a saddle-point structure and is typical for a mixed FEM. The well-known BabuškaBrezzi condition [1, 4, 17] must be fulfilled to obtain useful solutions for all involved fields.
This may be achieved by adjusting the orders of the approximation spaces for the different
fields and is further detailed in the numerical results. It is noted that stabilization may
be employed to circumvent the Babuška-Brezzi condition rather than to fulfill it, see, e.g.,
[17, 30, 31] which is, however, beyond the scope of this work.

3.3.4

Stationary Navier-Stokes flow

The discrete weak form of the stationary Navier-Stokes problem reads: Given density
h
% ∈ R+ , viscosity µ ∈ R+ , body force % · g h (x) in Γh , and traction t̂ (x) on ∂ΓhN , find
the velocity field uh (x) ∈ Suh , pressure field ph (x) ∈ Sph , and Lagrange multiplier field


λh (x) ∈ Sλh such that for all test functions whu , wph , wλh ∈ Vuh × Vph × Vλh , there holds in
Γh
%·
−
+

Z

ZΓ

whu ·





h

whu · t̂ ds +

∂ΓN
nel Z
X

e=1 Γe



uh · ∇cov
uh − g h dA +
Γ

τSUPG



Z
Γ

Z

wph · divΓ uh dA +


 h

uh · ∇cov
whu · % ·
Γ



ZΓ





h
h h
∇dir
dA +
Γ wu : σ u , p

Γ



Z



Γ



λh · whu · nhΓ dA



wλh · uh · nhΓ dA






uh · ∇cov
uh − g h − divΓ σ uh , ph
Γ

i

= 0.

18

Governing equations

The equations related to the different field equations were added up for brevity. The last
row adds a stabilization term which is needed to obtain stable solutions for flows at high
Reynolds numbers [12, 24]. In particular, the streamline upwind Petrov-Galerkin (SUPG)
method is used for the stabilization. Different definitions of the stabilization parameter
τSUPG are found [44, 49, 48] and

τSUPG =





2
∆t

2

2 kue k
+
he

!2

4µ
+
he 2



2 −1/2


is used herein with element-averaged velocity ue , element length he and ∆t → ∞ for
the stationary case. When stabilization is not necessary because no oscillations occur,
τSUPG = 0. Note that in the stabilization term, second-order derivatives appear (only in
the element interiors). The definition of tangential second-order derivatives is given, e.g.,
in [10].
Element assembly results in a non-linear system of equations of the form






 







K? + C (u) G? L   u   f 
 



·  p  =  0 ,
GT
0 0 
 



LT
0 0
λ
0

(3.33)

which is no longer symmetric (partly) due to the advection matrix C (u). The distinguishing feature of K? and G? (compared to K and G of the Stokes problem) are the added
SUPG-stabilization terms. The issues related to mixed FEMs and the Babuška-Brezzi
condition remain relevant.

3.3.5

Instationary Navier-Stokes flow

The discrete weak form of the instationary Navier-Stokes problem is: Given density % ∈
h
R+ , viscosity µ ∈ R+ , body force % · g h (x, t) in Γh × τ , traction t̂ (x, t) on ∂ΓhN × τ ,
and initial condition uh0 (x) on Γh at t = 0 according to (3.11), find the velocity field




uh (x, t) ∈ L2 τ ; Suh , pressure field ph (x, t) ∈ L2 τ ; Sph , and Lagrange multiplier field








λh (x, t) ∈ L2 τ ; Sλh such that for all test functions whu , wph , wλh ∈ Vuh × Vph × Vλh , there

19
holds in Γh × τ
%·
−
+

Z

ZΓ



whu

e=1 Γe

h

· ∂t u + u ·

whu

∂ΓN
nel Z
X



h

h

· t̂ ds +

τSUPG

Z
Γ



wph

∇cov
Γ



h

u −g

h



Γ
h

· divΓ u dA +
 h



dA +

Z

Z
Γ

wλh





h
∇dir
Γ wu
h

· u ·

nhΓ







h

h

: σ u ,p



dA +

Z



Γ



λh · whu · nhΓ dA

dA






uh · ∇cov
whu · % · ∂t uh + uh · ∇cov
uh − g h − divΓ σ uh , ph
Γ
Γ

i

= 0.

This yields a system of non-linear semi-discrete equations for t ∈ τ
M · u̇ (t) + (K? + C (u)) · u (t) + G? · p (t) + L · λ (t) = f (t) ,
GT · u (t) = 0,
LT · u (t) = 0,
with initial condition u (0). This system may be advanced in time by using finite difference
schemes and the Crank-Nicolson method is employed herein.

4

Numerical results

The following error measures are computed in the convergence studies. When analytic
(exact) velocity and pressure fields, uex and pex , are known, the velocity error is determined
by
s
εu =

3
X

Z 

i=1

Γ

uhi (x) − ui,ex (x)

2

dA

(4.1)

and the pressure error calculated as
εp =

sZ
Γ

(ph (x) − pex (x))2 dA.

(4.2)

When analytic solutions are not available, it is useful to evaluate the error of the FE
approximations in the strong form of the momentum or continuity equations, integrated
over the domain. For the example of stationary Stokes flow, the corresponding residual
errors are defined as
εmom =

sZ

Γ





P · divΓ σ uh , ph + f h

2

dA

(4.3)

20

Numerical results

and
εcont =

sZ

Γ

divΓ uh

2

dA.

(4.4)

This can be easily extended to the case of Navier-Stokes flows where the advection term is
added to the integrand in (4.3). Also the error in the tangential velocity constraint from
Eq. (3.3) may be computed in a similar manner. The evaluation of the error εmom involves
second-order derivatives and convergence can only be expected for higher-order elements
and sufficiently smooth solutions.

4.1

Stokes flow on an axisymmetric surface

A test case is developed for which analytic solutions are available. An axisymmetric surface
with height L = 5 and radius
r (z) = 1 + 1/5 · sin (1 + 3 · z) ,

z ∈ [0, L] ,

is generated as illustrated in Fig. 2(a). Let r0 = r(0) and r00 =
one may also define Γ based on the map x (a) : R2 → R3 ,



x (a) = 



cos a · r(b)
sin a · r(b)
b

dr(0)
.
dz

In parametrized form,







with a ∈ [0, 2π] , b ∈ [0, 5] .

The lower boundary at z = 0 is the Dirichlet boundary ∂ΓD , where the inflow in co-normal
direction of the manifold is prescribed as


û (x) =

u?
ku? k





and u? = 

r00 · cos θ
r00 · sin θ
1




,


with angle θ given by tan θ = y/x. The upper boundary at z = L is the outflow boundary where zero-tractions are applied as Neumann boundary conditions. The density and
viscosity are set to % = 1 and µ = 0.01, respectively.

4.1

21

Stokes flow on an axisymmetric surface

(a) domain

(b) quadmesh,
nz = 6

(c) quadmesh,
nz = 10

(d) quadmesh,
nz = 20

(e)
trimesh,
nz = 6

Figure 2: Manifold for the axisymmetric test case and meshes with different resolutions.
In (a), the blue arrows indicate the direction of the flow at the inflow (bottom) and outflow
(top) boundaries.

The mass flow on the lower boundary is
Q0 =

Z
∂ΓD

û (x) · n∂Γ ds = 2π · r0

and due to mass conservation, the mass flow along the height follows as Q (z) = 2π · r0/r(z).
As the flow field is expected to be axisymmetric for this test case, and the tangential
velocity constraint applies, one may compute the velocity components as






uex (x)
vex (x)
wex (x)









r0

=

r

r·

1+




 2 · 

dr
dz

dr
dz
dr
dz

1

· x/r
· y/r




.


See Fig. 3 for a graphical representation. It is noted that the mass flow Q (z), velocity
magnitude kuk, and the vertical velocity component w are only functions of z, that is,
they do not vary in x- and y-directions.
Finite element approximations are carried out on various meshes composed by triangular
or quadrilateral Lagrange elements of different orders. For the convergence studies, meshes
with nz = {4, 6, 10, 14, 20, 30, 40, 60} elements over the height are chosen; the number of
elements in circumferential direction is nθ = round (2πr0/L · nz ). The meshes are perturbed,
as illustrated in Figs. 2(b) to (d), to avoid perfectly axisymmetric meshes which, otherwise,

22

Numerical results

5

height z

4
3
2
1
|u|
w
r
0
0 0.5 1 1.5
velocities [m/s] / radius [m]

(a) u (x)

(b) v (x)

(c) w (x)

(d) ku (x)k

(e) w (z), ku (z)k

Figure 3: Exact solution for the axisymmetric test case.
could have improved the convergence rates for this special case.
The individual element orders used for the convergence studies are indicated by a 4-tuple
{kgeom , ku , kp , kλ }. To be precise, this tuple summarizes the employed orders for the geometry, kgeom , the velocities, ku , the pressure, kp , and the Lagrange multiplier for enforcing
the tangential velocity constraint, kλ . For each tuple, meshes with different resolutions
(given by nz and nθ ) are considered and errors calculated, each time resulting in one curve
in the convergence plots as indicated in the legends.
Systematic studies of different combinations of element orders showed that equal-order
approximations for the velocity and pressure, i.e., kp = ku do not converge satisfactory (or
at all), which is well-known from the standard context of the incompressible Navier-Stokes
equations in 2D and 3D due to the Babuška-Brezzi condition. For the studies outlined in
this paper, we shall choose kp = ku − 1 which is a popular choice for FEM approximations
of classical incompressible flows and known as Taylor-Hood elements [46].
For the first study, we use 2 ≤ ku ≤ 5, kp = ku − 1, and kλ = ku − 1 which, lateron,
becomes the recommended standard setting. Convergence plots for εu and εp are given in
Fig. 4. The thick solid lines are for kgeom = ku + 1. It is noteworthy that for quadrilateral
elements, setting kgeom = ku leads to almost identical results as seen from the thin dashed
lines in Figs. 4(c) and (d). This does not necessarily hold for triangular elements, see
Figs. 4(a) and (b), where the convergence may drop by one order when setting kgeom = ku
rather than kgeom = ku + 1. This is later confirmed for the errors εmom and εcont in Fig. 6.
Therefore, we recommend to choose the geometry one order higher than ku which is done
in the remainder of this work. Another reason is that the normal vector n is present in

4.1

23

Stokes flow on an axisymmetric surface

0

10

−1

10
−1

10

−2

10

10

2

−2

1

εu10

−3

3

3,
4,
5,
6,
2,
3,
4,
5,

1
−4

10

−5

4
1

10

−6

10

2

−3

5
1

0.025

2,
3,
4,
5,
2,
3,
4,
5,

1,
2,
3,
4,
1,
2,
3,
4,

0.1

1
2
3
4
1
2
3
4

10
εp
−4
10

1
3

−5

10

−6

3,
4,
5,
6,
2,
3,
4,
5,

1
4
1

10

5
1

0.025

0.25

2,
3,
4,
5,
2,
3,
4,
5,

1,
2,
3,
4,
1,
2,
3,
4,

0.1

1
2
3
4
1
2
3
4
0.25

h

h

(a) tri, εu

(b) tri, εp

0

10

−1

10
−1

10

−2

10

εu −3
10

−2

1
3

−4

−6

10

3,
4,
5,
6,
2,
3,
4,
5,

4
1

10

2

−3

1

10

−5

10

2

5
1

0.025

0.1
h

(c) quad, εu

2,
3,
4,
5,
2,
3,
4,
5,

1,
2,
3,
4,
1,
2,
3,
4,

1
2
3
4
1
2
3
4
0.25

10
εp
−4
10

1
3

3,
4,
5,
6,
2,
3,
4,
5,

1

−5

10

4

−6

1

−7

1

10
10

5

0.025

0.1

2,
3,
4,
5,
2,
3,
4,
5,

1,
2,
3,
4,
1,
2,
3,
4,

1
2
3
4
1
2
3
4
0.25

h

(d) quad, εp

Figure 4: Convergence results in εu and εp for the axisymmetric test case, (a) and (b) for
triangular elements, (c) and (d) for quadrilateral elements. The legends decode the orders
{kgeom , ku , kp , kλ } of the meshes.

24

Numerical results

0

1012

10

-1

1011

10

-2

10

2
1

εu10-3

3

10-4

1

10-5

1

4

5

10-6
10-7

3,
4,
5,
4,
5,
3,
4,
5,

1

0.025

0.1
h

(a) quad, εmom

2,
3,
4,
3,
4,
2,
3,
4,

1,
2,
3,
2,
3,
1,
2,
3,

1
1
1
2
3
2
3
4
0.25

1
6

1010

κ 10

3,
4,
5,
4,
5,
3,
4,
5,

9

108

1
2

107

2,
3,
4,
3,
4,
2,
3,
4,

1,
2,
3,
2,
3,
1,
2,
3,

1
1
1
2
3
2
3
4

106
105
0.025

0.1

0.25

h

(b) quad, κ

Figure 5: Influence of the order of the Lagrange multiplier field for enforcing the tangential velocity constraint: (a) convergence results in εu and (b) conditioning κ for the
axisymmetric test case.
the governing equations and is computed based on the Jacobi matrix, i.e., first derivatives
of the element mappings of order kgeom are involved.
It is important to note in Fig. 4 that the convergence rates in the pressure are optimal,
mp = kp +1, however, in the velocities one order sub-optimal, mu = ku . We have traced this
back to the influence of the order kλ of the Lagrange multiplier field. This is demonstrated
in Fig. 5 where (a) shows the error εu and (b) the condition number κ of the corresponding
system of equations (obtained with Matlab’s condest-function). As before, kgeom = ku + 1
and kp = ku − 1. Fig. 4 shows that setting kλ = 1 yields convergence rates mu = 2
independent of the other orders (black lines). Setting kλ = ku yields optimal convergence
rates mu = ku + 1 for the velocities (red lines), however, there is a dramatic influence on
the conditioning which scales with κ ∼ O (h−6 ) in this case rather than with O (h−2 ) for
all choices where kλ < ku . Therefore, we set kλ = ku − 1 in the following and accept the
sub-optimal convergence in the velocities.
Next, the error is observed in the strong form of the momentum and continuity equations,
εmom and εcont , see Eqs. (4.3) and (4.4). Results for kp = kλ = ku − 1 are depicted in Fig. 6
for triangular and quadrilateral elements. Again, the thick lines refer to kgeom = ku + 1 and
the thin dashed lines to kgeom = ku . As mentioned above for the L2 -errors in the velocities
and pressure, this makes a difference (of one order) for triangular elements, however, not
for quadrilateral elements. When using kgeom = ku + 1 on the safe side, the convergence

4.2

25

Driven cavity flows on manifolds

100
10-1
10-1
10-2

1
1

εmom
10

-3

10

-4

2

1

4

3

10

-5

3,
4,
5,
6,
2,
3,
4,
5,

2,
3,
4,
5,
2,
3,
4,
5,

1,
2,
3,
4,
1,
2,
3,
4,

1
2
3
4
1
2
3
4

10-2

εcont
10-3

1
3
1

10-4
10-5

5

4
1

1

1

0.025

0.1

0.025

10

0.1

10

3,
4,
5,
6,
2,
3,
4,
5,

1
2
1

-4

4

3
1

0.25

0

10-1

1

10

1
2
3
4
1
2
3
4

(b) tri, εcont

10-1

-3

1,
2,
3,
4,
1,
2,
3,
4,

h

(a) tri, εmom

εmom

2,
3,
4,
5,
2,
3,
4,
5,

1

0.25

h

10-2

3,
4,
5,
6,
2,
3,
4,
5,

2

2,
3,
4,
5,
2,
3,
4,
5,

1,
2,
3,
4,
1,
2,
3,
4,

1
2
3
4
1
2
3
4

10

-2

εcont
-3
10
10

-4

10

-5

1
3
1

0.1

0.25

h

(c) quad, εmom

5

4
1

1

0.025

3,
4,
5,
6,
2,
3,
4,
5,

2

2,
3,
4,
5,
2,
3,
4,
5,

1,
2,
3,
4,
1,
2,
3,
4,

1
2
3
4
1
2
3
4

1

0.025

0.1

0.25

h

(d) quad, εcont

Figure 6: Convergence results in εmom and εcont for the axisymmetric test case, (a) and
(b) for triangular elements, (c) and (d) for quadrilateral elements. The legends decode the
orders {kgeom , ku , kp , kλ } of the meshes.
rate in εmom is mmom = ku − 1 as expected due to the presence of second-order derivatives
of u in the momentum equations. The expected convergence rate in εcont is mcont = ku
due to the presence of first order derivatives of u in the continuity equation.

4.2

Driven cavity flows on manifolds

The stationary Navier-Stokes model is considered in this example. Starting point is the
driven cavity for the case of a flat 2D domain as depicted in Fig. 7(a). This case has
well-documented reference solutions for a variety of Reynolds numbers [21]. There, a flow
inside a quadratic domain Ω2D = (0, 1) × (0, 1) with no-slip boundary conditions on the

26

Numerical results

(b) α = 0.4

(a) flat

(c) β = 0.4

Figure 7: Different manifolds for the driven cavity test case, (a) flat, (b) map A with
α = 0.4, and (c) map B with β = 0.4.

left, right and lower wall develops under a shear flow of u = 1.0 and v = 0.0 applied on the
upper boundary until a stationary solution is reached. The Reynolds number is computed
as Re = % · u · L/µ.
Herein, the situation is extended to curved surfaces in 3D by deforming the flat 2D domain
Ω2D in z-direction using functions z(x, y). In particular, two different maps A and B are
used,




map A :

z(x, y) = α · −1 + 8x + 2y − 8x2 · (1 − y) ,

map B :

z(x, y) = β · (1 − y) · sin ((2x − 1) π) · cos ((2y − 1) π) ,

where α and β scale the height in z-direction, see Figs. 7(b) and (c) for examples. The
advantage is that for α = 0 and β = 0, the flat situation is recovered and the reference
solutions in [21] are relevant. We have confirmed that these solutions are recovered with
great accuracy also for any rigid body tranformation of Ω2D into three dimensions. The
density is chosen as % = 1 and two different viscosities of µ = 0.01 and µ = 0.001 leading
to Reynolds numbers of Re = 100 and Re = 1000 for the flat case, respectively. Solutions
for the velocity magnitude and pressure field for some example manifolds are displayed in
Fig. 8.
The meshes feature quadrilateral elements of different orders and are refined towards the
boundaries to capture the resulting boundary layers. See Fig. 9 for the meshes in Ω2D which

4.2

27

Driven cavity flows on manifolds

(a) α = 0.4, ku(x)k

(b) α = 0.4, p(x)

(c) β = 0.4, ku(x)k

(d) β = 0.4, p(x)

Figure 8: Velocity and pressure fields for the driven cavity test case with µ = 0.001 for
map A with α = 0.4 and map B with β = 0.4.

(a) 10 × 10 elements

(b) 20 × 20 elements

(c) 50 × 50 elements

Figure 9: Different meshes for the driven cavity test case in top view.

28

Numerical results

are mapped to 3D according to map A and B from above for various scaling coefficients
α and β. The number of elements per dimension is n = {10, 20, 30, 50, 70, 100}. For the
numerical studies, kp = kλ = ku − 1 and kgeom = ku + 1 is used as recommended above.
Just as for the reference solutions in [21], the results are presented as velocity profiles
along the horizontal and vertical centerlines in Ω2D . Fig. 10 shows the profiles for the
velocity component u along the vertical centerline and v along the horizontal centerline for
the two maps with different scaling factors α and β, respectively. The crosses indicating
the reference solution from [21] are only relevant for the flat case where α = β = 0.
The results for the velocity component w along the two centerlines are given in Fig. 11.
These results have the quality of benchmark solutions and have been obtained with ku = 4
and 100 elements per dimensions. The convergence of other element orders and mesh
resolutions towards these profiles has been confirmed, and a small selection is shown in
Fig. 12. Without stabilization, the typical oscillations are seen for this rather high Reynolds
number for coarse meshes with low order. As no analytical solutions for the velocities
and pressure are available, it is impossible to provide convergence results in εu and εp .
Furthermore, the singular pressure in the upper left and right corners lead to singularities
in the derivatives of other physical fields. Thus, it cannot be expected that (optimal)
convergence in εmom and εcont is achieved.

4.3

Flows on zero-level sets

The next test case shows the potential to solve flows on zero-level sets with the proposed
models. Stationary Stokes and Navier-Stokes flows are considered. The scalar function
φ (x) : R3 → R is based on [14] and defined as
φ (x) =




2

x2 + y 2 − 4
2

x2 − 1



2





+ x2 + z 2 − 4
2

+ y2 − 1



2

+ y2 + z2 − 4
2

+ z2 − 1

+

− 15.

The zero-isosurface of φ implies the compact manifold of interest, Γ = {x : φ (x) = 0}
and is depicted in Fig. 13. In a first step, meshes with linear triangular elements are
generated using distmesh [39]. A scaling parameter h may be chosen which defines an
average element length. In a second step, higher-order elements are mapped to this linear
surface mesh and their element nodes are “lifted” [14] such that they are on the manifold
Γ. Thereby, a higher-order accurate representation Γh is obtained.

4.3

29

Flows on zero-level sets

1

1

0.75

0.75

0.5

0.5
α=0
α = 0.2
α = 0.4
α = 0.6
α = 0.8
α = 1.0
reference

0.25
0
0

0.25

0.5

0.75

1

0.25
0
0

(a) map A, µ = 0.01

0.25

0.5

0.75

1

(b) map B, µ = 0.01

1

1

0.75

0.75

0.5

0.5
α=0
α = 0.2
α = 0.4
α = 0.6
α = 0.8
α = 1.0
reference

0.25
0
0

β=0
β = 0.1
β = 0.2
β = 0.3
β = 0.4
reference

0.25

0.5

0.75

1

(c) map A, µ = 0.001

β=0
β = 0.1
β = 0.2
β = 0.3
β = 0.4
reference

0.25
0
0

0.25

0.5

0.75

1

(d) map B, µ = 0.001

Figure 10: Velocity profiles for the driven cavity test case for different α and β. The vertical
profiles show the velocity component u (x), the horizontal profiles v (x). The scaling factor
of the velocities is 0.5.

30

Numerical results

1

1

0.75

0.75

0.5

0.5
α=0
α = 0.2
α = 0.4
α = 0.6
α = 0.8
α = 1.0

0.25

0
0

0.25

0.5

0.75

1

0.25

0
0

(a) map A, µ = 0.01
1

0.75

0.75

0.5

0.5

0
0

α=0
α = 0.2
α = 0.4
α = 0.6
α = 0.8
α = 1.0
0.25

0.5

0.75

0.5

0.75

=0
= 0.1
= 0.2
= 0.3
= 0.4

β
β
β
β
β

=0
= 0.1
= 0.2
= 0.3
= 0.4

1

(b) map B, µ = 0.01

1

0.25

0.25

β
β
β
β
β

1

0.25

0
0

(c) map A, µ = 0.001

0.25

0.5

0.75

1

(d) map B, µ = 0.001

Figure 11: Velocity profiles for the driven cavity test case for different α and β. The
horizontal and vertical profiles show the velocity component w (x) scaled by the factor 1.

1

1

0.75

0.75

0.5

0.5

10 × 10 elems, ku = 2
10 × 10 elems, ku = 3
10 × 10 elems, ku = 4
100 × 100 elems, ku = 4

0.25

0
0

0.25

0.5

0.75

1

(a) flat

20 × 20 elems, ku = 2
10 × 10 elems, ku = 3
10 × 10 elems, ku = 4
100 × 100 elems, ku = 4

0.25

0
0

0.25

0.5

0.75

1

(b) map A, α = 1.0

Figure 12: Velocity profiles for the driven cavity test case following Fig. 10. Results of
coarse meshes with 10 × 10 elements are compared to the high-accuracy results from above
(with 100 × 100 elements with ku = 4), (a) flat manifold, (b) manifold according to map
A with α = 1.0.

4.3

31

Flows on zero-level sets

(b) mesh with h = 0.3

(a) domain

(c) mesh with h = 0.1

Figure 13: Manifold for the zero-isosurface test case and meshes with different resolutions.

As there are no boundaries present, an accelaration field in z-direction drives the flow.
That is, on the right hand side, g = P · [0, 0, gz ]T where gz is determined by

gz (x) =




exp


0

z2
− 2
2σ0

!

with σ0 = 0.15 for x < 0 and y < 0,
else,

and visualized in Fig. 14(a). It is virtually non-zero only for the left front “pillar” of the
domain. The density is % = 1 and the viscosity is µ = 0.05. For the case of stationary
Navier-Stokes flow, the corresponding velocity magnitude, pressure fields and vorticity ω ?
according to Eq. (3.8) are seen in Figs. 14(b) to (d), respectively.
In the numerical studies, 2 ≤ ku ≤ 5, kp = kλ = ku − 1, and kgeom = ku + 1 are used. As
there is no analytical solution available, convergence results are only shown in εmom and
εcont in Fig. 15. Higher-order rates are clearly achieved. In order to make the solution
more quantitative, the velocity profiles for w (x) in the horizontal xy-plane (at z = 0) are
shown in Fig. 16. The four closed black lines represent the intersection of the plane with
the vertical “pillars” of the zero-isosurface. Fig. 16(a) shows w(x) as a third dimension,
and (b) shows the same result where w(x) is plotted in normal direction of the plane-pillar
intersections with a scaling factor of 0.4. A clear convergence to these profiles was observed
when using meshes with different resolutions and orders.

32

Numerical results

(a) gz (x)

(b) ku(x)k

(c) p(x)

(d) ω ? (x)

Figure 14: Physical fields of the zero-isosurface test case: (a) acceleration gz , (b) velocity
magnitude kuk, (c) pressure p, (d) vorticity ω ? .

101

10-1
1

εmom
10

εcont
10-2

1
0

10-1
0.03

ku
ku
ku
ku

ku
ku
ku
ku

=2
=3
=4
=5
2
1

2

=2
=3
=4
=5

3

1
4

3
1

0.1
h

(a) εmom

1

10-3
1

1

0.3

0.03

5

4

0.1
h

1

0.3

(b) εcont

Figure 15: Convergence results in (a) εmom and (b) εcont for the zero-isosurface test case.

4.4

33

Cylinder flows

2

y

1
Zero-level set
Stokes
Nav.-Stokes

0
−1
−2
−2

(a) w(x)

−1

0
x

1

2

(b) w(x)

Figure 16: Profiles for the vertical velocity w (x) in the plane with z = 0, the scaling factor
in (b) is 0.4.

4.4

Cylinder flows

As an example for the instationary Navier-Stokes equations, the following test case is based
on a channel flow around a cylinder according to [42]. The geometry is first described in
2D, labelled Ω2D , and later on mapped to obtain curved surfaces in 3D. In 2D, the cylinder
with a diameter of 0.1 is placed slightly unsymmetrically in y-direction of the channel in
[0, 2.20]×[0, 0.41], see Fig. 17(a). No-slip boundary conditions are applied on the upper and
lower wall and on the cylinder surface. A quadratic velocity profile for u, with umax = 1.5,
and v = 0 is applied at the inflow on the left side of the domain. At the outflow, tractionfree boundary conditions are used. The density and viscosity are prescribed as % = 1.0 and
µ = 0.001. This results in a Reynolds number of Re = % · um · L/µ = 100 when taking the
cylinder diameter as a length scale L and the average inflow velocity um = 1.0 at the inflow.
At this Reynolds number, periodic flow patterns known as the Kármán vortex street are
observed behind the cylinder. Reference solutions are given for the lift and drag coefficients
cL and cD of the cylinder [42] and the current implementation confirms these numbers for
the flat case (i.e., in 2D or when the flat 2D domain is transformed by a rigid body motion
to 3D). The reference Strouhal number St = D/ (um T ), with the diameter D = 0.1 of the
cylinder, and the time T for 2 periods of the curve of cD , is given as 0.295 ≤ St ≤ 0.305,
resulting in a frequency of about f = 3.33 1/s.
The 2D domain is mapped to three dimensions using two different maps. Assume that the
coordinates of the 2D domain Ω2D , as seen in Fig. 17(a), are given in coordinates (a, b).

34

Numerical results

(a) flat

(b) map A

(c) map B

Figure 17: Manifolds for the cylinder flow test case.
Map A, x (a) : R2 → R3 , is defined as
π·a
· (b + 0.35),
2.2


π·a
· (b + 0.35),
y(a) = sin
2.2





x(a) = cos

q

q



z(x(a), y(a)) = 2 + 1/2 x2 + y 2 − sin 3 x2 + y 2 .
For map B, we first define an intermediate mapping r (a) : R2 → R3 applying some twist
to the domain,
r (a) = a,
s (a) = − (1 + q (a)) · (b − 0.205) · cos (π/6 (1 − 25/11 · a)) ,
t (a) = − (1 + q (a)) · (b − 0.205) · sin (π/6 (1 − 25/11 · a)) .
with q(a) = −0.2/2.42 · a2 + 0.44/2.42 · a. This is further mapped by x (r) : R3 → R3
defined as
x (r) = cos (50/198 · π · r) · (s + 6/5) ,
y (r) = sin (50/198 · π · r) · (s + 6/5) ,
z (r) = t + 1/5 sin (3r) .
The resulting curved manifolds according to map A and B are visualized in Figs. 17(b)
and (c), respectively. Note that also the inflow velocities are mapped accordingly based on
the Jacobians of the respective mappings to ensure that they are in the tangent space at

4.4

35

Cylinder flows

(a) ku(x)k

(c) ω ? (x)

(b) p(x)

Figure 18: Physical fields for the cylinder flow test case according to map A: (a) velocity
magnitude kuk, (b) pressure p, (c) vorticity ω ? .

(a) ku(x)k

(c) ω ? (x)

(b) p(x)

Figure 19: Physical fields for the cylinder flow test case according to map B: (a) velocity
magnitude kuk, (b) pressure p, (c) vorticity ω ? .
∂ΓD .
The initial condition on the manifolds is u0 (x) = 0. The observed time interval is τ = [0, 6]
and the inflow velocities are ramped by a cubic function in time,

R (t) =



−2 · (t/t? )3

1

+ 3 · (t/t? )2

for t ≤ t? ,
else,

with t? = 0.96. That is, after t? , the full velocity profile is active at the inflow. Figs. 18
and 19 show the velocity magnitude, pressure field, and vorticity ω ? at time t = 6 for the
two mappings. The expected vortex shedding can be clearly seen.
Two different meshes with 972 and 1920 elements each are used which are refined at
the no-slip boundaries to resolve the boundary layers. They are visualized for Ω2D in
Fig. 20 and mapped to the manifolds accordingly. We use element orders of kgeom = 4,
ku = 3, kp = 2 and kλ = 2 in the numerical studies shown here. Higher orders achieved
virtually indistinguishable results for the quantities shown below. It is also noted that the

36

Conclusions

(a) coarse, 972 elements

(b) fine, 1920 elements

Figure 20: Coarse and fine mesh for the cylinder flow test case; these meshes are later
mapped according to map A or B.
Crank Nicolson method used for the time discretization is only second-order accurate. For
the time discretization nstep = {150, 300, 600, 1200, 2400, 4800} time steps are used. To
make the results more quantitative, the stresses at the cylinder wall are summed up to
obtain a force resultant F (t) = kF (t)k in 3D. This is the equivalent of the lift and drag
coefficients for the flat 2D case. Furthermore, the pressure difference between the front
and back position of the cylinder (in Ω2D , mapped to three dimensions) is computed, i.e.,
∆p (t) = pfront (t) − pback (t).
The results for map A are shown in Fig. 21 for the different number of time steps. It can
be seen that after about 2s, the expected vortex shedding is almost established. After 3s,
the resulting oscillations remain virtually unchanged. The time interval [5.2, 6] is shown
in more detail in Figs. 21(b) and (d) for F (t) and ∆p (t), respectively. The convergence
with increasing number of time steps is clearly demonstrated. Fig. 22 shows the results in
the same style for map B; the same conclusions may be drawn. The spatial convergence
is investigated in Fig. 23 where it is found that the coarse and fine mesh employed here
obtain very similar results for the chosen element orders. The frequency of the oscillations
for map A is fA = 2.191 1/s and for map B is fB = 3.078 1/s; for the flat case the frequency
is f = 3.33 1/s.

5

Conclusions

The surface FEM with higher-order elements is applied to solve Stokes and Navier-Stokes
flows on (fixed) manifolds. For the governing equations, the classical gradient and divergence operators are replaced by their tangential counterparts. An additional constraint is
needed to ensure that the velocities are in the tangent space of the manifold. Stabilization
is required for the case of Navier-Stokes flows at large Reynolds numbers and the standard
streamline-upwind Petrov-Galerkin (SUPG) approach is used herein.

37

0.25
0.24

0.15
0.1
0.05
0
0

1

2

3
time t

4

150
300
600
1200
2400
4800
5
6

force resultant F

force resultant F

0.2
0.23

0.22

0.21

0.2
5.2

(a) F (t), t ∈ (0, 6)

150
300
600
1200
2400
4800
5.4

5.6
time t

5.8

6

(b) F (t), t ∈ (5.2, 6)
2.7

2
1.5
1
0.5
0
0

1

2

3
4
time t

(c) ∆p(t), t ∈ (0, 6)

150
300
600
1200
2400
4800
5
6

pressure difference ∆p

pressure difference ∆p

2.5
2.6
2.5
2.4
2.3
2.2
2.1
5.2

150
300
600
1200
2400
4800
5.4

5.6
time t

5.8

6

(d) ∆p(t), t ∈ (5.2, 6)

Figure 21: Force resultant F (t) and pressure difference ∆p(t) obtained on the fine mesh of
the cylinder flow test case according to map A.

38

0.2
0.18
0.16
0.14
0.12
0.1
0.08
0.06
0.04
0.02
0
0

0.195

1

2

3
time t

150
300
600
1200
2400
4800
5
6

4

force resultant F

force resultant F

Conclusions

0.19
0.185
0.18
0.175
5.2

(a) F (t), t ∈ (0, 6)

5.8

6

2.85

2.5
2
1.5
1
150
300
600
1200
2400
4800
5
6

0.5
1

2

3
4
time t

pressure difference ∆p

pressure difference ∆p

5.6
time t

(b) F (t), t ∈ (5.2, 6)

3

0
0

150
300
600
1200
2400
4800
5.4

2.8

2.75

2.7

2.65
5.2

(c) ∆p(t), t ∈ (0, 6)

150
300
600
1200
2400
4800
5.4

5.6
time t

5.8

6

(d) ∆p(t), t ∈ (5.2, 6)

Figure 22: Force resultant F (t) and pressure difference ∆p(t) obtained on the fine mesh of
the cylinder flow test case according to map B.

0.24
force resultant F

force resultant F

0.195
0.23

0.22

0.21

0.2
5.2

0.19
0.185
0.18

coarse
fine
5.4

5.6
time t

5.8

(a) map A, F (t), t ∈ (5.2, 6)

6

0.175
5.2

coarse
fine
5.4

5.6
time t

5.8

6

(b) map B, F (t), t ∈ (5.2, 6)

Figure 23: Force resultant F (t) for the cylinder flow test case according to map A and B
for the fine and coarse meshes with 4800 time steps.

39
For the discretization, the surface FEM is employed with quadrilateral or triangular elements. Element spaces of different orders are used for (i) the geometric approximation
of the manifold, kgeom , (ii) the approximation of the velocity fields, ku , (iii) the pressure
field, kp , and (iv) the Lagrange multiplier field for the enforcement of the tangential velocity constraint, kλ . The choice of these orders affects the properties of the resulting FEM
in terms of conditioning, accuracy, and stability. Particularly useful combinations for a
chosen order ku are kgeom = ku + 1 and kp = kλ = ku − 1. Some benchmark test cases
for flows on manifolds are proposed and higher-order convergence rates are achieved. The
notation used in this work is closely related to the engineering literature for the FEM in
fluid mechanics. Implementational matters are outlined.
There is a large potential for future research related to this work: One may investigate
different stabilization methods such as Galerkin least-squares stabilization and variational
multiscale methods. Stabilization may also be useful to circumvent the Babuška-Brezzi
condition and enable equal-order shape functions for the velocities and pressure. The
tangential velocity constraint may be more efficiently enforced based on penalty methods
or other Lagrange multiplier approaches such as the Uzawa method. We believe that flows
on manifolds have a strong potential for fundamental research in mathematics, physics,
and engineering.

6

Acknowledgements

The fruitful discussions with Dr. Sven Groß, Dr. Thomas Rüberg and Prof. Günther Of
are gratefully acknowledged.

References
[1] Babuška, I.: Error-bounds for finite element method. Numer. Math., 16, 322 – 333,
1971.
[2] Blaauwendraad, J.; Hoefakker, J.H.: Structural Shell Analysis, Vol. 200, Solid Mechanics and Its Applications. Springer, Berlin, 2014.
[3] Bothe, D.; Prüss, J.: On the Two-Phase Navier-Stokes Equations with BoussinesqScriven Surface Fluid. J. math. fluid mech., 12, 133 – 150, 2010.

40

REFERENCES

[4] Brezzi, F.: On the existence, uniqueness and approximation of saddle-point problems
arising from Lagrange multipliers. RAIRO Anal. Numér., R-2, 129 – 151, 1974.
[5] Brooks, A.N.; Hughes, T.J.R.: Streamline upwind/Petrov-Galerkin formulations for
convection dominated flows with particular emphasis on the incompressible NavierStokes equations. Comp. Methods Appl. Mech. Engrg., 32, 199 – 259, 1982.
[6] Chan, C.H.; Yoneda, T.: On the stationary Navier-Stokes flow with isotropic streamlines in all latitudes on a sphere or a 2D hyperbolic space. Dynamics of Partial
Differential Equations, 10, 209 – 254, 2013.
[7] Chapelle, D.; Bathe, K.J.: The Finite Element Analysis of Shells – Fundamentals.
Computational Fluid and Solid Mechanics. Springer, Berlin, 2011.
[8] Deckelnick, K.; Elliott, C.M.; Ranner, T.: Unfitted finite element methods using bulk
meshes for surface partial differential equations. SIAM J. Numer. Anal., 52, 2137 –
2162, 2014.
[9] Delfour, M.C.; Zolésio, J.P.: Tangential Differential Equations for Dynamical ThinShallow Shells. Journal of Differential Equations, 128, 125 – 167, 1996.
[10] Delfour, M.C.; Zolésio, J.P.: Shapes and geometries—Metrics, Analysis, Differential
Calculus, and Optimization. SIAM, Philadelphia, PA, 2011.
[11] Demlow, A.: Higher-order finite element methods and pointwise error estimates for
elliptic problems on surfaces. SIAM J. Numer. Anal., 47, 805 – 827, 2009.
[12] Donea, J.; Huerta, A.: Finite Element Methods for Flow Problems. John Wiley &
Sons, Chichester, 2003.
[13] Dziuk, G.: Finite Elements for the Beltrami operator on arbitrary surfaces, Chapter 6,
142 – 155. Springer Berlin Heidelberg, Berlin, Heidelberg, 1988.
[14] Dziuk, G.; Elliott, C.M.: Finite element methods for surface PDEs. Acta Numerica,
22, 289 – 396, 2013.
[15] Ebin, D.G.; Marsden, J.: Groups of diffeomorphisms and the motion of an incompressible fluid. Ann. of Math. (2), 92, 102 – 163, 1970.

REFERENCES

41

[16] Edwards, D.A.; Brenner, H.; Wasan, D.T.: Interfacial Transport Processes and Rheology. Butterworth-Heinemann, Oxford, 1991.
[17] Franca, L.P.; Hughes, T.J.R.: Two classes of mixed finite element methods. Comp.
Methods Appl. Mech. Engrg., 69, 89 – 129, 1988.
[18] Fries, T.P.; Omerović, S.: Higher-order accurate integration of implicit geometries.
Internat. J. Numer. Methods Engrg., 106, 323 – 371, 2016.
[19] Fries, T.P.; Omerović, S.; Schöllhammer, D.; Steidl, J.: Higher-order meshing of
implicit geometries—part I: Integration and interpolation in cut elements. Comp.
Methods Appl. Mech. Engrg., 313, 759 – 784, 2017.
[20] Fries, T.P.; Schöllhammer, D.: Higher-order meshing of implicit geometries—part II:
Approximations on manifolds. Comp. Methods Appl. Mech. Engrg., 326, 270 – 297,
2017.
[21] Ghia, U.; Ghia, K.N.; Shin, C.T.: High-Re solutions for incompressible flow using the
Navier-Stokes equations and a multi-grid method. J. Comput. Phys., 48, 387 – 411,
1982.
[22] Grande, J.; Reusken, A.: A Higher Order Finite Element Method for Partial Differential Equations on Surfaces. SIAM J. Numer. Anal., 54, 388 – 414, 2016.
[23] Gravemeier, V.: The variational multiscale method for laminar and turbulent flow.
Archives of Computational Methods in Engineering, 13, 249, 2006.
[24] Gresho, P.M.; Sani, R.L.: Incompressible Flow and the Finite Element Method, Vol.
1+2. John Wiley & Sons, Chichester, 2000.
[25] Gross, S.; Reusken, A.: Numerical Methods for Two-phase Incompressible Flows,
Vol. 40, Springer Series in Computational Mathematics. Springer, Berlin, 2011.
[26] Gurtin, M.E.; Murdoch, I.A.: A continuum theory of elastic material surfaces. Archive
for Rational Mechanics and Analysis, 57, 1975.
[27] Hansbo, A.; Hansbo, P.: A finite element method for the simulation of strong and
weak discontinuities in solid mechanics. Comp. Methods Appl. Mech. Engrg., 193,
3523 – 3540, 2004.

42

REFERENCES

[28] Hansbo, P.; Larson, M.G.: Finite element modeling of a linear membrane shell problem
using tangential differential calculus. Comp. Methods Appl. Mech. Engrg., 270, 1 –
14, 2014.
[29] Hughes, T.J.R.; Feijóo, G.R.; Mazzei, L.; Quincy, J.P.: The variational multiscale
method—a paradigm for computational mechanics. Comp. Methods Appl. Mech. Engrg., 166, 3 – 24, 1998.
[30] Hughes, T.J.R.; Franca, L.P.: A new finite element formulation for computational
fluid dynamics: VII. The Stokes problem with various well-posed boundary conditions:
symmetric formulations that converge for all velocity/pressure spaces. Comp. Methods
Appl. Mech. Engrg., 65, 85 – 96, 1987.
[31] Hughes, T.J.R.; Franca, L.P.; Balestra, M.: A new finite element formulation for
computational fluid dynamics: V. Circumventing the Babuška-Brezzi condition: a
stable Petrov-Galerkin formulation of the Stokes problem accommodating equal-order
interpolations. Comp. Methods Appl. Mech. Engrg., 59, 85 – 99, 1986.
[32] Hughes, T.J.R.; Franca, L.P.; Hulbert, G.M.: A new finite element formulation
for computational fluid dynamics: VIII. The Galerkin/Least-squares method for
advective-diffusive equations. Comp. Methods Appl. Mech. Engrg., 73, 173 – 189,
1989.
[33] Jankuhn, T.; Olshanskii, M.A.; Reusken, A.: Incompressible fluid problems on embedded surfaces: Modeling and variational formulations. arXiv:1702.02989, 2017.
[34] Khesin, B.; Misiołek, G.: Euler and Navier-Stokes equations on the hyperbolic plane.
Proceedings of the National Academy of Sciences (PNAS), 109, 18324 – 18326, 2012.
[35] Koba, H.; Liu, C.; Giga, Y.: Energetic variational approaches for incompressible fluid
systems on an evolving surface. Quart. Appl. Math., 75, 359 – 389, 2017.
[36] Kobayashi, M.H.: On the Navier-Stokes equations on manifolds with curvature. J.
Eng. Math., 60, 55 – 68, 2008.
[37] Nitschke, I.; Voigt, A.; Wensch, J.: A finite element approach to incompressible twophase flow on manifolds. Journal of Fluid Mechanics, 708, 418 – 438, 2012.

REFERENCES

43

[38] Olshanskii, M.A.; Reusken, A.; Xu, X.: A stabilized finite element method for
advection-diffusion equations on surfaces. IMA Journal of Numerical Analysis, 34(2),
732 – 758, 2014.
[39] Persson, P.O.; Strang, G.: A Simple Mesh Generator in MATLAB. SIAM Review,
46, 329 – 345, 2004.
[40] Reusken, A.: Analysis of trace finite element methods for surface partial differential
equations. IMA J Numer Anal, 35, 1568 – 1590, 2014.
[41] Reuther, S.; Voigt, A.: Solving the incompressible surface Navier-Stokes equation by
surface finite elements. Physics of Fluids, 30, 012107, 2018.
[42] Schäfer, M.; Turek, S.: Benchmark Computations of Laminar Flow around a Cylinder. In Flow Simulation with High-Performance Computers II. (Hirschel, E.H., Ed.),
Vieweg Verlag, Braunschweig, 1996.
[43] Scriven, L.: Dynamics of a fluid interface equation of motion for Newtonian surface
fluids. Chemical Engineering Science, 12, 98 – 108, 1960.
[44] Shakib, F.; Hughes, T.J.R.; Johan, Z.: A new finite element formulation for computational fluid dynamics: X. The compressible Euler and Navier-Stokes equations. Comp.
Methods Appl. Mech. Engrg., 89, 141 – 219, 1991.
[45] Slattery, J.C.; Sagis, L.; Oh, E.S.: Interfacial transport phenomena, Vol. 2. Springer,
Berlin, 2007.
[46] Taylor, C.; Hood, P.: A Numerical Solution of the Navier-Stokes Equations Using the
Finite Element Technique. Computers & Fluids, 1, 73 – 100, 1973.
[47] Temam, R.: Infitite dimensional dynamical systems in mechanics and physics. Springer, Berlin, 1988.
[48] Tezduyar, T.; Sathe, S.: Stabilization Parameters in SUPG and PSPG Formulations.
J. Comput. Appl. Math., 4, 71 – 88, 2003.
[49] Tezduyar, T.E.; Osawa, Y.: Finite element stabilization parameters computed from
element matrices and vectors. Comp. Methods Appl. Mech. Engrg., 190, 411 – 430,
2000.

