for Surface Vector Partial Differential Equations
Von der Fakult¨at f¨ur Mathematik, Informatik und Naturwissenschaften
der RWTH Aachen University zur Erlangung des akademischen Grades eines Doktors der Naturwissenschaften genehmigte Dissertation
vorgelegt von
Thomas Jankuhn, M.Sc.
aus Berlin
Berichter: Univ.-Prof. Dr. Arnold Reusken Univ.-Prof. Dr. Maxim Olshanskii Tag der m¨undlichen Pr¨ufung: 08.06.2021
Diese Dissertation ist auf den Internetseiten der Universit¨atsbibliothek verf¨ugbar.
First of all, I would like to express my sincere gratitude to my advisor Prof. Dr. Arnold Reusken for his support and guidance during my pursuit of a doctoral degree. I also want to express my gratitude to my co-advisor Prof. Dr. Maxim Olshanskii for his efforts related to this thesis.
Furthermore, I would like to thank Prof. Dr. Christoph Lehrenfeld for providing continuous support on the Netgen/NGSolve code with his add-on package ngsxfem. Without his ngsxfem package I would still look for a bug in my DROPS code.
I would like to thank my colleagues from IGPM (Institut f¨ur Geometrie und Praktische Mathematik), especially my close colleagues from LNM (Lehrstuhl f¨ur Numerische Mathe- matik). A special thanks goes to Hauke Saß for the many fruitful discussions.
At last, I want to thank my friends and family for their continuous support.
In this thesis we develop and analyze efficient and highly accurate finite element meth- ods for the numerical simulation of the surface vector-Laplace equation and the surface (Navier-)Stokes equations on a stationary surface.
First, the incompressible surface Navier-Stokes equations on an evolving surface are de- rived. The derivation is based on fundamental continuum mechanical principles for a viscous material surface embedded in an ambient continuum medium. The resulting equations are formulated in terms of tangential differential operators in Cartesian coordinates, which makes the formulation more convenient for our numerical purposes. We use a directional splitting of the system into a coupled system of equations for the tangential density flow and the normal velocity of the surface in order to gain a further insight into both of these components. On an a priori given stationary surface the normal velocity vanishes and one is interested in the tangential density flow only. One can deduce the following three simplified models on a stationary surface: the surface vector-Laplace equation, the surface Stokes equations and the surface Navier-Stokes equations.
In surface flow problems we have the constraint that the flow must be tangential to the surface. It is not obvious how this constraint should be treated numerically. To better under- stand this issue we first consider discretizations of the surface vector-Laplace equation. We use the higher order parametric trace finite element method. Three different natural tech- niques for treating the tangential constraint are studied, namely a consistent penalty method, a simpler inconsistent penalty method and a Lagrange multiplier method. A complete anal- ysis of all three methods is presented that reveals how the discretization error bounds in the energy norm depend on relevant parameters such as the degree of the polynomials used for the approximation of the solution, the degree of the polynomials used for the approximation of the level set function that characterizes the approximation of the surface, the penalty parameter, the order of the normal vector approximation used in the penalty terms in both penalty methods and the degree of polynomials used for the approximation of the Lagrange multiplier. Furthermore, for the consistent penalty method we derive a new optimalL2-error bound. The results are confirmed and illustrated by numerical experiments. In a compar- ison of the three methods we conclude that the consistent penalty method has significant advantages over the other two methods.
For the surface Stokes equations an isoparametric Pk-Pk−1 trace Taylor-Hood finite ele- ment discretization (k≥2) is introduced. For the tangential constraint we employ the same consistent penalty technique used and analyzed for the surface vector-Laplace equation. With the help of a uniform inf-sup stability result for the Taylor-Hood pair well-posedness is shown and an optimal energy norm error bound is derived. Furthermore, we present a new error analysis that yields an optimal L2-error bound. Results of numerical experiments confirm the error analysis.
To discretize the surface Navier-Stokes equations on a stationary surface we apply an isoparametric Pk-Pk−1 trace Taylor-Hood finite element method (k ≥ 2) in space and a second order accurate semi-implicit BDF2 method in time. For the case that no outer forces
as a validation of numerical simulation results.
1 Introduction 1
1.1 Motivation . . . 1
1.2 Mathematical models . . . 2
1.3 Finite element discretization methods . . . 4
1.3.1 Numerical challenges . . . 5
1.3.2 Main contributions . . . 6
1.4 Outline of the thesis . . . 7
2 Basic notation and surface differential operators 11 3 The trace finite element method 15 3.1 Basic concepts of the trace finite element method . . . 16
3.2 The higher order parametric trace finite element method . . . 20
3.2.1 The parametric trace finite element space . . . 21
3.2.2 Notation and fundamental estimates . . . 25
3.2.3 The parametric TraceFEM . . . 30
3.2.4 Numerical experiments . . . 33
3.3 On the normal vector approximation . . . 34
3.3.1 Discretization error analysis . . . 36
3.4 Summary . . . 45
4 Modeling of incompressible fluid problems on embedded surfaces 47 4.1 Preliminaries . . . 49
4.2 Modeling of material surface flows . . . 51
4.3 Directional splitting of the surface Navier-Stokes equations . . . 53
4.4 Summary . . . 56
5 Well-posed variational formulations 57 5.1 Well-posed surface vector-Laplace variational formulations . . . 59
5.2 Well-posed surface Stokes variational formulations . . . 63
5.3 Summary . . . 66
6 Trace finite element methods for the surface vector-Laplace equation 67 6.1 Parametric trace finite element methods . . . 69
6.2 Well-posedness of discretizations . . . 72
6.3 Strang lemmas . . . 78
6.4 Approximation error bounds . . . 79
6.5 Consistency error analysis . . . 80
6.5.1 Geometry errors . . . 81
6.5.2 Discrete Korn’s type inequality . . . 85
6.5.3 Consistency error of the penalty methods . . . 86
6.5.4 Consistency error of the Lagrange multiplier method . . . 88
6.6 Energy norm error bounds . . . 90
6.7 L2-error bound . . . 92
6.8 Numerical experiments . . . 102
6.8.1 Results for the penalty methods . . . 103
6.8.2 Results for the Lagrange multiplier method . . . 108
6.9 Comparison of methods . . . 111
7 A trace finite element method for the surface Stokes equations 115 7.1 Isoparametric trace finite element method . . . 116
7.2 Well-posedness of discretization . . . 118
7.3 Error analysis . . . 119
7.3.1 Consistency error analysis . . . 120
7.3.2 Finite element error bounds . . . 125
7.4 Numerical Experiments . . . 128
7.4.1 Results on the unit sphere . . . 129
7.4.2 Results on a more complex geometry . . . 130
7.4.3 Results on a torus . . . 132
7.5 Summary . . . 134
8 A trace finite element method for the surface Navier-Stokes equations 135 8.1 Trace finite element method . . . 136
8.2 Numerical Experiments . . . 138
8.2.1 Numerical dissipation on the unit sphere . . . 140
8.2.2 Geometric interaction on a biconcave shape . . . 143
8.3 Summary . . . 144
9 Summary and Outlook 147 9.1 Summary . . . 147
9.2 Outlook . . . 148
Bibliography 153
1.1 Motivation
Fluid equations posed on surfaces have played an increasingly important role in many ap- plications in chemistry, biology and engineering. Surface fluids arise in the mathematical modeling of foams, emulsions, liquid crystals and biological membranes and often appear at phase boundaries in multiphase systems [Bre13, Pru17].
For example in engineering applications one is interested in the lifetime and stability of foams and emulsions [CCAE+13, Bik13, ZZQ+17, BCP99]. The characteristics of the inter- face separating the dispersed and continuous phases play a key role for the stability. These can be influenced by using surface-active molecules such as detergents, emulsifiers, lipids and proteins [Wil00, McC04]. Liquid crystals [DGP93, SS74], which are a state of matter inter- mediate between that of a crystalline solid and an isotropic liquid, can be found in different applications in technology and biology such as liquid-crystal displays or cell membranes. For example certain animal tissue, the so-called epithelium, can be modeled as an active nematic liquid crystal [SDN+17]. Experimental and numerical studies of such liquid crystals show a tight coupling between the topology of the surface, geometric properties of the surface, and topological defect dynamics on the surface [DSL00, FNVU+07, DSDLC12, NRV19]. Bi- ological membranes are interfaces composed of lipids and proteins [Eva18]. Such membranes exhibit a coupling between out-of-plane elastic membrane shape changes, resulting from pro- tein binding and unbinding reactions, and in-plane viscous lipid fluid flow and in-plane species diffusion making it a viscoelastic material surface [RDA13]. Examples are certain cell mem- branes such as the red blood cell membrane. Numerical simulations and modeling of biological membranes can be found in, e.g., [SOSM20, FHL+11, TSMA19, AD09, ESSW11, RAM+13].
The understanding of most of these interfacial/surface phenomena is still very poor. In particular, there is a lack of validated mathematical models that describe these processes ap- propriately. Research on modeling and simulation of interfacial flow and transport processes is a very active and strongly growing research field. Until now most research results in this field have been published in the engineering literature. There are only a few papers that have appeared in the numerical mathematics literature and that address the development and analysis of numerical simulation methods for these models.
In the mathematical modeling of interfacial flow and transport processes different scalar, vector or tensor valued surface partial differential equations arise. In particular, the viscous surface fluid flow of certain biological membranes, liquid crystals or foam and emulsions can be modeled using the (generalized) surface Navier-Stokes equations on evolving surfaces.
The importance of incorporating a surface viscosity into the interfacial constitutive law of the continuum mechanical model of two-phase flows was first discovered by Boussinesq [Bou13].
Later, Scriven [Scr60] extended this model to material interfaces of arbitrary curvature.
Our long-term goal is the numerical simulation of fluid flows on evolving surfaces (or interfaces). As shown in Chapter 4 the surface Navier-Stokes equations used to model incom- pressible surface fluid systems on evolving surfaces admit a natural splitting into a coupled
system of equations for the tangential fluid flow and the normal motion of the surface. Such a splitting yields a subproblem that is very similar to the surface (Navier-)Stokes equations on astationarysurface. Hence, the detailed study of finite element discretizations for the surface (Navier-)Stokes equations on astationary surface is an important step in the development of a robust and efficient finite element solver for the surface Navier-Stokes equations onevolving surfaces.
In this thesis we will restrict ourselves to the case of stationary surfaces. Our main goal is to develop and analyze efficient and highly accurate finite element methods for the numerical simulation of the surface (Navier-)Stokes equations on stationary surfaces.
Several difficulties arise in the discretization of the surface (Navier-)Stokes equations in addition to those well-known for the equations posed in Euclidean domains. For example one has to approximate the surface and the covariant derivatives. Another difficulty stems from the need to recover a tangential vector field on the surface. The latter difficulty does not occur in the discretization of scalar surface partial differential equations.
In the next section we give a brief overview of the surface Navier-Stokes equations and the simplifications thereof that we use in this thesis. In the remaining part of this chapter we discuss the finite element method of choice, explain the key challenges arising in the development and analysis of this finite element method and outline the main contributions of the thesis (Section 1.3). In Section 1.4 we give an outline of the thesis.
1.2 Mathematical models
In this section we give a brief overview of the mathematical model (cf. [JOR18]) and its sim- plifications that we will consider throughout the thesis. A detailed derivation and discussion can be found in Chapter 4.
There are different approaches to derive the surface Navier-Stokes equations on evolving surfaces. Since a long-term goal is to incorporate the surface Navier-Stokes equations into a two-phase flow problem, we present and use a model that is based on fundamental surface continuum mechanical principles and assumes the conservation laws of mass and momentum for an evolving surface embedded in an ambient continuum medium.
The formulation of fluid equations on surfaces involves the covariant derivative of a tan- gential vector field. The covariant derivative of a tangential vector field on general surfaces is defined via intrinsic variables of a local coordinate system on a surface. In this thesis we exploit the embedding of the two-dimensional surface inR3 and formulate the tangential differential operators in Cartesian coordinates. This notion in Cartesian coordinates does not involve local coordinates or parametrizations, which makes the formulation more convenient for our numerical purposes. For the definition of the surface differential operators we refer to Chapter 2.
Let Γ(t), t ∈ [0,∞), be a smooth, closed, evolving material surface embedded in R3. We denote by u = u(x, t) the smooth velocity field of the density flow on Γ(t). The four- dimensional space-time manifold defined by the evolution of Γ(t) is denoted by
S := [
t∈[0,∞)
(Γ(t)× {t}), S ⊂R4.
Based on the conservation laws of mass and momentum for a viscous material surface embed- ded in an ambient continuum medium one can derive the surface Navier-Stokes equations.
Surface Navier-Stokes equations
For a given force vector field f: S →R3 determine the velocity flow field u:S → R3 and the surface pressurep:S →Rsuch that
ρu˙ =−∇Γp+ 2µdivΓ(Es(u)) +f+pκn onS,
divΓu= 0 on S.
Here, ˙u:= ∂u∂t + (∇u)udenotes the material derivative, Es(u) := 12(∇Γu+ (∇Γu)T) is the surface rate-of-strain tensor, n is the outward pointing unit normal on Γ(t), κ denotes the mean curvature,ρdenotes the density andµis the interface shear viscosity. The bulk medium interacts with the fluidic surface through the area forcesf. Together with the equations ˙ρ= 0 and VΓ = u·n, where VΓ is the normal velocity of Γ(t), and suitable initial conditions the surface Navier-Stokes equations form a closed system of six equations for six unknownsu,p, ρ, andVΓ. These surface Navier-Stokes equations determineu=Pu+ (u·n)n=uT +uNn, which describes the dynamics both of the normal velocity uN of the manifold and of the tangential density flowuT over Γ(t).
In this thesis we only consider finite element methods for three simplified models of the surface Navier-Stokes equations. All three simplified models are posed on a stationary surface.
On an a priori given stationary surface one is interested in the tangential density flow only.
Using a directional splitting of the system above we obtain a coupled system equations for the tangential density flow and the normal velocity of the surface. If we restrict ourselves to a stationary surface then by assumption we have u·n= 0. The model above reduces to the incompressible surface Navier-Stokes equations for the tangential velocityuT on a stationary surface Γ.
Surface Navier-Stokes equations on a stationary surface
For a given tangential force vector field fT: Γ×[0,∞)→R3 determine the tangential velocity uT: Γ×[0,∞)→R3 and the surface pressurep: Γ×[0,∞)→Rsuch that
ρ ∂uT
∂t + (uT · ∇Γ)uT
=−∇Γp+ 2µPdivΓEs(uT) +fT on Γ×[0,∞),
divΓuT = 0 on Γ×[0,∞).
Together with a suitable initial condition these equations form a closed system.
A simplification of these equations is given by the surface Stokes equations. To derive these equations we additionally assume that the viscous surface forces dominate and therefore omit the nonlinear term (uT · ∇Γ)uT. We also restrict ourselves to the equilibrium flow problem, i.e., ∂u∂tT = 0, setµ= 12 and introduce an additional mass term.
Surface Stokes equations
For a given tangential force vector field fT: Γ → R3, a source term g: Γ → R and α ∈ [0,1] determine the tangential velocity uT: Γ → R3 and the surface pressure p: Γ→Rsuch that
−PdivΓ(Es(uT)) +αuT +∇Γp=fT on Γ, divΓuT =g on Γ.
We add the additional zero order termαuT on the left-hand side to avoid technical details related to the kernel of the differential operator uT 7→ −PdivΓ(Es(uT)). To better under- stand some phenomena that are related to the approximation of the differential operator for the velocity part of the surface Stokes problem we introduce the surface vector-Laplace equation.
Surface vector-Laplace equation
For a given tangential force vector field fT: Γ → R3 and α ∈ [0,1] determine the tangential vector field uT: Γ→R3 such that
−PdivΓ(Es(uT)) +αuT =fT on Γ.
In particular, we use this problem to better understand how to numerically treat the constraint that the solution must be tangential.
We consider these three problems posed on stationary surfaces in reverse order to the one presented here. In Chapter 5 we treat well-posed variational formulations for the surface vector-Laplace equation and surface Stokes equations. Based on these formulations we con- sider and analyze finite element discretizations in Chapters 6 and 7. Results of numerical experiments of a finite element method for the surface Navier-Stokes equations are presented in Chapter 8.
1.3 Finite element discretization methods
The study of finite element methods for partial differential equations on general surfaces can be traced back to the pioneering work of Dziuk in [Dzi88]. In that paper, a finite element method for thescalar Laplace-Beltrami equation on a stationary surface Γ is considered. The surface is approximated by a regular family Γhof consistent triangulations where the vertices of the triangulation are assumed to lie on the surface Γ. The finite element space is defined on Γh and consists of scalar functions that are continuous on Γh and linear on each triangle of the triangulation. This so-called surface finite element method (SFEM) was extended by Demlow to higher order finite elements in [Dem09] and by Elliot et al. [DE07, DE13a, EV15]
to evolving surfaces Γ(t). For the evolving surface finite element method (ESFEM) the initial surface Γ(0) is approximated by a triangulation as described above with the nodes lying on the exact surface. The vertices are advected by a bulk velocity field w along the space-time manifold so that the vertices of the approximated surface Γh(t) remain on Γ(t) for all times t. Hence, this method is based on a Lagrangian tracking of the surface. The finite element space consists of scalar functions that are continuous on Γh(t) and linear on each triangle of the triangulation Γh(t) for fixedt. If a surface undergoes strong deformations or topological changes or it is defined implicitly, e.g., as the zero level of some level set function, then numerical methods based on a Lagrangian tracking of the surface have certain disadvantages.
For an overview of the surface finite element method and other finite element techniques for surface partial differential equations we refer to [DE13b].
In this thesis we consider the trace finite element method (TraceFEM) also known as CutFEM. TraceFEM belongs to the class of Eulerianunfitted finite element methods, where the method avoids a triangulation of the surface or any other fitting of a mesh to the surface.
Unfitted finite element methods originated from interface problems, where two phases are separated by the phase boundary Γ. Here, one uses a standard finite element space defined
on the bulk triangulation of the two phases and then “cuts” the functions from this space at the interface Γ, which is equivalent to taking the restrictions of these functions to one of the phases. An overview of such unfitted finite element methods for interface problems can be found in [LR17b]. The trace finite element method uses the traces of the finite element functions from the bulk triangulation on the interface or surface. The first trace finite element method was introduced in [ORG09]. The main ideas are as follows. We define a surface independent triangulation of a fixed polygonal bulk domain Ω ⊂ R3, such that Γ ⊂ Ω. On the bulk mesh of Ω we define the finite element space Vh of continuous and piecewise linear functions. For the surface representation the level set technique is used, i.e., the surface Γ is implicitly given by the zero level of a level set functionφ. The approximated surface Γh is defined by the zero level of the Lagrangian interpolant (or some other piecewise linear approximation) φh ∈Vh of φ. The traces on Γh of the finite element functions inVh are used to define the method. TraceFEM was extended to higher order finite elements in [GR16, GLR18] and to evolving surfaces in [ORX14a, OR14, Gra14, OX17, LOX18]. For an overview of the trace finite element method for (evolving) surfaces we refer to [OR17].
The TraceFEM approaches for evolving surfaces differ in the way they handle the numerical integration in time but all approaches are Eulerian methods since they use a fixed surface independent bulk triangulation (in space). In comparison to surface finite element methods the trace finite element methods, or unfitted finite element methods in general, have certain attractive properties concerning the flexibility (no remeshing). Moreover, results of numerical experiments show that the combination of the level set technique with TraceFEM yields a very robust discretization method that is capable of handling the case of surfaces undergoing topological changes without any modifications and without stability restrictions on mesh or time step sizes [GOR14, OX17]. This is a unique property of TraceFEM in comparison to SFEM.
On stationary surfaces the SFEM seems to be the better choice in terms of complexity [BJP+21]. The choice of the geometrically unfitted TraceFEM over the SFEM is motivated by the long-term goal to consider fluid equations on evolving surfaces, including cases, where a parametrization of the surface is not explicitly available and the surface may undergo large deformations or even topological changes. Furthermore, we may already need the volume mesh for the discretization of a flow problem in the bulk, where the interface Γ(t) is transported by the bulk velocity.
1.3.1 Numerical challenges
We briefly summarize the key issues in the development and analysis of the higher order trace finite element methods for the surface vector-Laplace equation, surface Stokes equations and surface Navier-Stokes equations.
Tangential flow constraint. In surface flow problems we have the constraint that the flow must be tangential to the surface. It is not obvious how this constraint should be treated numerically. We cannot enforce this constraint strongly, since this may lead to a so-called
“locking effect” (cf. Remark 5.7). Instead we allow for a velocity solution that is fully three- dimensional and use a penalty or Lagrange multiplier approach to control the component of the velocity field that is normal to the surface.
TraceFEM for vector PDEs. In recent years, several papers on the analysis of TraceFEM
for scalar surface PDEs have appeared. The analysis of TraceFEM, or finite element methods in general, for surfacevector PDEs is still at an early stage. In particular, for the vector case we have to derive bounds for the consistency error caused by the geometry approximation.
For the TraceFEM such estimates have not been available in the literature.
Sufficiently accurate geometry approximation. The original TraceFEM, which uses a piecewise linear approximation of the surface, is not suitable for higher order methods due to the dominating geometry error of second order. A major issue with higher order methods in the context of unfitted finite element methods in general is the realization of numerical integration on domains which are represented implicitly by a level set function. In [Leh16, GLR18] an isoparametric mapping is introduced and analyzed in the context of the Laplace-Beltrami equation. This mapping yields a higher order accurate approximation of the surface and allows for numerical integration on the piecewise planar surface approxima- tion.
Stabilization. From the analysis of TraceFEM for scalar surface PDEs it is known that one needs to apply an appropriate stabilization to control instabilities caused by “small cuts” and to improve the algebraic properties of the resulting algebraic systems. In addition, for our purposes the stabilization must be suitable for higher order finite element approximations.
The normal gradient volume stabilization introduced and analyzed in [GR16] meets these criteria.
Stable finite element pair. The Pk-Pk−1 continuous Taylor-Hood element (k ≥ 2) is a popular inf-sup stable finite element pair for the discretization of incompressible fluid flow problems in the Euclidean space. It is not obvious whether the continuous Taylor-Hood elements defined in the bulk space are uniformly inf-sup stable for TraceFEM. In [ORZ21, JORZ20] it is shown that a particular stabilization technique for the pressure variable is essential for the central inf-sup stability result to hold.
1.3.2 Main contributions
Below we outline the main contributions of this thesis. Part of the results presented in the thesis are based on the joint publications [JOR18, GJOR18, JR21, JR19, JORZ20]:
Modeling of fluid problems. Based on fundamental continuum mechanical prin- ciples such as the conservation laws of mass and momentum for a viscous material surface embedded in an ambient continuum medium we derive the incompressible sur- face Navier-Stokes equations on an evolving surface. With a directional splitting of the system into a coupled system of equations for the tangential density flow and the nor- mal velocity of the surface we give a further insight into the interplay between normal and tangential velocity. The resulting equations are formulated in terms of tangen- tial differential operators in Cartesian coordinates, which makes the formulation more convenient for our numerical purposes.
Surface vector-Laplace equation. We introduce three different well-posed vari- ational formulations for the surface vector-Laplace problem, namely an inconsistent penalty approach, a consistent penalty approach and a Lagrange multiplier approach, to handle the tangential constraint weakly. Based on these formulations we propose
three higher order parametric trace finite element discretizations for the surface vector- Laplace equation. For all three methods well-posedness is shown and and a complete discretization error analysis is presented. The analysis reveals, how the discretization error bounds in the energy norm depend on relevant parameters such as the degree of the polynomials used for the approximation of the solution, the degree of the poly- nomials used for the approximation of the level set function that characterizes the approximation of the surface and the order of the normal vector approximation used in the penalty terms in both penalty methods. The results are confirmed and further illustrated by numerical experiments.
Surface Stokes equations. Based on the consistent penalty approach to handle the tangential constraint that was studied for the surface vector-Laplace problem we present an isoparametric Pk-Pk−1 trace Taylor-Hood finite element discretization (k ≥ 2) for the surface Stokes equations. Using a uniform inf-sup stability result for the Taylor- Hood pair well-posedness is shown and an optimal energy norm error bound is derived.
Results of numerical experiments confirm the error analysis.
Further important contributions not published in the literature are the following:
L2-error analysis. For the consistent penalty trace finite element discretization of the surface vector-Laplace equation and for the presented isoparametricPk-Pk−1 trace Taylor-Hood finite element discretization (k ≥ 2) for the surface Stokes equations we derive an optimalL2-error bound in each case.
Surface Navier-Stokes equations. We apply an isoparametric trace Taylor-Hood finite element method in space and a second order accurate semi-implicit BDF2 method in time to discretize the surface Navier-Stokes equations on a stationary surface. For the case that no outer forces are present, asymptotic properties of the solution in the continuous case are derived and used as a validation of numerical simulation results.
Normal vector approximation. We present an error analysis and results of numeri- cal experiments for a discretization of the Laplace-Beltrami equation that uses theexact normal vector to discretize the Laplace-Beltrami operator. The error analysis shows that this choice of the normal vector approximation introduces an additional geometric error that leads to a suboptimal order of convergence. This shows that the normal vector approximation in the surface differential operators should be chosen carefully.
1.4 Outline of the thesis
The thesis is structured as follows.
In Chapter 2 we present some basic notation and introduce the surface differential operators. Furthermore, we define the Sobolev spaces on the surface Γ.
InChapter 3we give an overview of the concepts of the trace finite element method us- ing the Laplace-Beltrami equation as a model problem. We discuss the basic concepts of TraceFEM such as the definition of the finite element spaces and the need for stabilizing the bilinear forms. Furthermore, we introduce the higher order parametric trace finite element method and present results of the analysis and numerical experiments for the
Laplace-Beltrami equation. These results are known from the literature and are used in Chapters 6–8 for the definition and analysis of the trace finite element discretizations for the surface vector-Lapalce equation, the surface Stokes equations and the surface Navier-Stokes equations on a stationary surface. We present a new analysis and results of numerical experiments for a discretization of the Laplace-Beltrami equation using the normal vector of the exact surface to discretize the Laplace-Beltrami operator.
InChapter 4we derive and discuss the surface Navier-Stokes equations that are used as the mathematical model throughout this thesis. Three simplifications of this model are used as a basis for the variational formulations and trace finite element methods in Chapters 5–8. We introduce the assumptions of our model that are based on the conservation laws of mass and momentum. With these assumptions we derive the surface Navier-Stokes equations for a viscous material surface. In a next step we apply a directional splitting of the system to obtain a coupled system of equations for the tangential density flow uT and the normal velocity of the surface uN. Finally, we consider two simplified models of the surface Navier-Stokes equations. One of them is the surface Navier-Stokes equations on a stationary surface. The derivation that we present in this chapter is taken from [JOR18].
In Chapter 5 we introduce three well-posed variational formulations for both the vector-Laplace equation and the surface Stokes equations. These formulations are suit- able for Galerkin finite element methods and are the basis for the trace finite element methods introduced and analyzed in Chapters 6–8. The difference of these three meth- ods lies in the way how they deal with the tangential constraint, which is a crucial point in the development and analysis of finite element methods for the surface vector- Laplace equation and the surface Stokes equations. For both problems we introduce an inconsistent penalty formulation, a consistent penalty formulation and a Lagrange multiplier formulation to deal with the tangential constraint.
In Chapter 6 we introduce and analyze three different higher order parametric trace finite element methods for the surface vector-Laplace equation that are based on the inconsistent penalty formulation, the consistent penalty formulation and the Lagrange multiplier formulation introduced in Chapter 5. We show the well-posedness of these methods and present a complete discretization error analysis in the energy norm based on a Strang lemma, which bounds the discretization error in terms of an approximation error and a consistency error. Most of these results can also be found in [JR21]. For the trace finite element method that is based on the consistent penalty formulation we present an optimal L2-error analysis. Such an L2-error analysis is not known in the literature. We present results of numerical experiments that support our theoretical results and also give further insights into these methods. We conclude this chapter with a comparison of the three methods presented and analyzed in this chapter.
InChapter 7we consider a higher order isoparametric trace Taylor-Hood finite element discretizaton for the surface Stokes equations. The method uses the same technique for the tangential constraint as the consistent penalty method introduced and analyzed in Chapter 6 for the surface vector-Laplace equation. We present well-posedness results of the method and an optimal discretization error analysis in the energy norm. Most of these results are published in [JORZ20]. We present an error analysis that yields an optimal L2-error bound. Such a bound is not known in the literature. We provide
results of numerical experiments that support the theoretical results and that illustrate new phenomena that do not occur in Euclidean domains.
InChapter 8we apply a higher order isoparametric trace Taylor-Hood finite element method in space and a second order accurate BDF2 method from [OY19] in time for the discretizaton of the surface Navier-Stokes equations on a stationary surface. The trace Taylor-Hood finite element method that is used for the space discretization was introduced and analyzed in Chapter 7 for the surface Stokes equations. For the case that no outer forces are present we derive a sharp estimate for the kinetic energy decay and asymptotic properties of the solution of the surface Navier-Stokes equations. In two numerical experiments we illustrate the numerical dissipation of the method and the interplay of curvature and vortices in flow on the surface.
InChapter 9we summarize the main results and give an outlook on topics that have not been addressed, but that we consider to be of interest for future research.
In this chapter we recall some basics of tangential calculus for two-dimensional (evolving) manifolds embedded inR3. For further reading on manifolds we refer to the standard litera- ture, e.g., [War83, AMR12]. For a more elaborate introduction into the notion of tangential calculus on manifolds in the context of finite elements we refer to [DE13b].
We consider a two-dimensional closed, connected, orientable, smooth evolving manifold Γ = Γ(t)⊂R3 fort≥0. Note that a closed manifold is compact and without boundary. We comment on the smoothness properties of Γ(t) in Remark 2.2. For everyt≥0 we define an open tubular neighborhood of Γ(t) by
Uδ(t) :=
x∈R3 | |d(x, t)|< δ ,
with δ > 0 and d(·, t) the signed distance function to Γ(t), which we take negative in the interior of Γ(t). On Uδ(t) we define n(x, t) :=∇d(x, t), the outward pointing unit normal on Γ(t),H(x, t) :=∇2d(x, t), the Weingarten map,P(x, t) :=I−n(x, t)n(x, t)T, the orthogonal projection onto the tangential plane and p(x, t) := x−d(x, t)n(x, t), the closest point pro- jection. The Gauss curvature, i.e., the product of the two principal curvatures, is denoted by K and the mean curvature is denoted by κ:= tr(H). We assume δ >0 to be sufficiently small such that for everyt≥0 the decomposition
x=p(x, t) +d(x, t)n(x, t) (2.1) is unique for all x ∈ Uδ(t). For functions v: Γ(t) → Rm, m = 1,2, . . ., we define the componentwise constant normal extensionve:=v◦pto the neighborhoodUδ(t). For vector functions w:R3 → R3 we denote by∇w(x) the Jacobian at x ∈ R3, (∇w(x))i,j = ∂w∂xi
j(x), 1 ≤ i, j ≤ 3. Note that for scalar functions f:R3 → R the gradient ∇f is still the usual column vector with (∇f)i=∂if,i= 1,2,3.
Remark 2.1. The unit normal n is extended constantly in normal direction by definition, i.e., ne =n. We further note that the Weingarten mappingH, which is the derivative of a function that is extended constantly in normal direction, however, is not extended constantly in normal direction.
Surface differential operators
We consider Γ := Γ(t) for some fixed t ≥ 0. By Ck(Γ) we denote the space of k-times continuously differentiable functions on Γ that is defined in the usual way using local charts (cf. [DE13b]). We use the notation Ck(Γ) := Ck(Γ)3 and Ck(Γ)3×3 := Ck(Γ)3×3 for the componentwise defined vector and tensor spaces. We define the tangential derivative∇Γg of a scalar functiong∈C1(Γ) and the covariant derivative∇Γvof a vector functionv∈C1(Γ):
∇Γg(x) :=P(x)∇g(x), x∈Γ, ∇Γv(x) :=P(x)∇v(x)P(x), x∈Γ, (2.2)
where g and v are smooth extensions to the neighborhood Uδ. Note that the tangential and covariant derivatives are independent of the chosen extension, since they only depend on the values of g and v on Γ (see [DE13b, Lemma 2.4]). In general, there is no “natural extension” of functions defined on Γ. Thus, we will often use the constant normal extension for functions defined on Γ. In particular, for the constant normal extension we have on Γ the identities∇ge =∇(g◦p) =P∇ge =∇Γge and ∇ve=∇(v◦p) = ∇veP forg∈C1(Γ) and v∈C1(Γ). Note that the definition of a covariant derivative in (2.2) coincides with the covariant derivative from differential geometry for tangential vector fieldsv∈C1(Γ).
For scalar functions f, g∈C1(Γ) and vector functionsu,v∈C1(Γ) we have the following product rules:
∇Γ(f g) =g∇Γf +f∇Γg,
∇Γ(u·v) = (∇Γu)Tv+ (∇Γv)Tu, if Pu=u, Pv=v,
∇Γ(fu) =f∇Γu+Pu(∇Γf)T.
The surface divergence operator for vector-valued functions u ∈ C1(Γ) and tensor-valued functions A∈C1(Γ)3×3 are defined by
divΓv:= tr(∇Γv) = tr(P(∇v)P) = tr(P(∇v)) = tr((∇v)P), divΓA:= divΓ(ATe1), divΓ(ATe2), divΓ(ATe3)T
,
withei being theith unit vector inR3. Forv∈C1(Γ) andφ∈C1(Γ) we define the surface curl operators (cf. [Reu20b])
curlΓu:= divΓ(u×n) = (∇Γ×u)·n, curlΓφ:=n× ∇Γφ.
The Laplace-Beltrami operator applied to a scalar functionf ∈C2(Γ) is given by
∆Γf := divΓ(∇Γf).
Remark 2.2. In the following chapters we will always assume that for any given t≥0 the surface Γ(t) is at leastC2-smooth. Together with the fact that the manifold is embedded in R3, this allows us to formulate the surface differential operators on Γ(t) in terms of differential operators in the Euclidean space R3 with respect to the standard basis inR3. For example, for aC2-manifold, a normal extension off ∈Ck(Γ),k= 0,1, is a Ck-smooth function in the neighborhood Uδ(t).
Sobolev spaces
For a bounded subset S ⊂ R3 we denote by Lp(S) the standard Lp-spaces of (equivalence classes of) measurable functions f:S → R (with respect to the Lebesgue measure), with finite norm given by
kfkLp(S):=
Z
S|f|pdx 1p
, forp∈[1,∞), kfkL∞(S):= inf
N ⊂S
|N |=0
sup
S\N|f|, (2.3)
is a Hilbert space with the scalar product h·,·iL2(S). For the componentwise defined vector Lebesgue spaces Lp(S) := Lp(S)3 and Lp(S)3×3 := Lp(S)3×3 we use the same notation for the scalar product and norm. We also use the space
L20(S) :={f ∈L2(S)| Z
S
fdx= 0}.
Similarly, we define the Lp-spaces with respect to the one- and two-dimensional Hausdorff measures. In this case we write dl for the one-dimensional Hausdorff measure and ds for the two-dimensional Hausdorff measure. Moreover, we adopt the notation dsh for the two-dimensional Hausdorff measure if we integrate over a discrete surface. For the two- dimensional Hausdorff measure of a setS ⊂R3 we use meas2(S).
For a Lebesgue-measurable subset Ω⊂R3, which is usually either open or closed and has non-empty interior, we denote by Wk,p(Ω) the standard Sobolev spaces equipped with the norm
kfkpWk,p(Ω):= X
|α|≤k
kDαfkpLp(Ω),
where α∈ N3 is a multi-index with |α|=α1+α2+α3 and Dα := ∂xα1∂|α|
1 ∂xα22∂xα33. For p = 2 we have that Hk(Ω) := Wk,2(Ω) is a Hilbert space. For the componentwise defined vector and tensor Sobolev spaces Wk,p(Ω) :=Wk,p(Ω)3 and Wk,p(Ω)3×3 :=Wk,p(Ω)3×3 we use the same notation for the scalar product and norm.
For further reading on the standard Lebesgue and Sobolev spaces we refer the reader to [BS07, AF03, Bre11, Eva10]. In particular, in [BS07] one can find the definition of Sobolev spaces that also includes closed subsets.
The definition of the Sobolev spaces on Γ follows along [DE13b]. We first define the weak derivative on Γ. The definition of the weak derivative is based on the formula for integration by parts on Γ for smooth functions. For the i-th component of the tangential differential operator∇Γ=P∇we use the notation∇i:=P3
k=1Pik∂k, where Pik is the (i, k)-entry ofP.
Definition 2.3. For a smooth, connected, orientable manifold Γ (possibly with boundary) we say that a functionf ∈L1(Γ) has the weak derivativevi=∇if ∈L1(Γ),i= 1,2,3, if for everyϕ∈C1(Γ) with compact support{x∈Γ|ϕ(x)6= 0} ⊂Γ, we have
Z
Γ
f∇iϕds=− Z
Γ
ϕvids+ Z
Γ
f ϕκnids.
With the weak derivative on Γ we define the Sobolev spaces on Γ.
Definition 2.4. For 1≤p≤ ∞ and k∈N the Sobolev spaces on Γ are recursively defined by
Wk,p(Γ) :={f ∈Wk−1,p(Γ)| ∇i1∇i2· · · ∇ikf ∈Lp(Γ), i1, . . . , ik∈ {1,2,3}}, (2.4)
whereW0,p(Γ) =Lp(Γ). The corresponding norm is defined by kfkpWk,p(Γ):=
k
X
l=0
X
i1,...,il
∈{1,2,3}
k∇i1· · · ∇ilfkpLp(Γ), forp∈[1,∞),
kfkWk,∞(Γ):= max
l=0,1,...,k max
i1,...,il
∈{1,2,3}
k∇i1· · · ∇ilfkL∞(Γ).
The spaceHk(Γ) :=Wk,2(Γ) is a Hilbert space. ByH∗k(Γ) we denote the space of functions inv ∈Hk(Γ) withR
Γvds= 0. The vector and tensor Sobolev spaces Wk,p(Γ) := Wk,p(Γ)3 and Wk,p(Γ)3×3 :=Wk,p(Γ)3×3 are defined componentwise. For the norms and scalar prod- uct on these Sobolev spaces we use the same notation as for the scalar-valued spaces. In particular, the surface Sobolev space of weakly differentiable vector-valued functions H1(Γ) is equipped with the norm
kuk2H1(Γ):=
Z
Γ
kuk22+k∇uek22ds. (2.5) Note thatk∇uek2 =k(∇ue)Tk2 and on Γ we have
(∇ue)T =P ∇ue1,∇ue2,∇ue3
=
3
X
i=1
∇ΓuieTi . (2.6) Hence, the norm is a natural extension to vector-valued functions of the usual scalar H1(Γ)- norm. We note that for the previous definitions we only need Γ to be C2-smooth. In this case the mean curvature κ of Γ is well defined.
Remark 2.5. The definition of the weak derivatives and, therefore, the definition of the Sobolev spacesWk,p(Γ) also applies to the case of smooth surfaces with Lipschitz boundaries, i.e., in particular for curved triangles or quadrilaterals. This will be needed in the next chapter to define a Sobolev space on the discrete surface, which is only locally smooth.
Remark 2.6. In the norm on H1(Γ) we only consider tangential derivatives of the compo- nents of u. This is different from the covariant derivative defined in (2.2), where we have an additional projection on the left-hand side. For tangential vector functions u ∈ H1(Γ) one can show (see [HLL20]) that the norm (2.5) is equivalent to
kuk2H1
tan(Γ):=
Z
Γkuk22+k∇Γuk22ds.
In this chapter we introduce the trace finite element method using the Laplace-Beltrami problem as a scalar model problem for an elliptic equation. Throughout this chapter we assume that Γ is a smooth, closed, connected and orientable manifold embedded in R3. The surface Γ is given by the zero level of a smooth level set function φ: Uδ → R, i.e., Γ = {x ∈ R3 | φ(x) = 0}. We start by recalling the well-known weak formulation of the Laplace-Beltrami problem:
Continuous Problem 3.1
For a given f ∈L20(Γ) determine u∈H∗1(Γ) such that Z
Γ∇Γu· ∇Γvds= Z
Γ
f vds for all v∈H1(Γ).
The solution u exists, is unique and satisfies u ∈ H2(Γ) with kukH2(Γ) ≤ ckfkL2(Γ) and some constantcindependent off (see [DE13b]). Iff ∈Hs(Γ),s >0, then the unique solution u satisfies the smoothness property u ∈Hs+2(Γ) and kukHs+2(Γ) ≤c(s)kfkHs(Γ) holds with some constant c(s) > 0 (cf. [War83, Fundamental Inequality 6.29]). The well-posedness follows from the Poincar´e inequality (cf. [DE13b]):
kvkL2(Γ)≤ck∇ΓvkL2(Γ) for allv ∈H∗1(Γ). (3.1) Note that since we introduce several different continuous and discrete formulations in the remainder, we use the term “Continuous Problem” for continuous variational problems and
“Discrete Problem” for discrete variational problems.
In Section 3.1 we present the basic concepts of the trace finite element method. We discuss the background and surface triangulation, the definition of the finite element spaces and the necessity for a stabilization. In this section we only consider continuous piecewise linear finite elements and a surface approximation where the surface is represented as the zero level of the continuous piecewise linear finite element interpolation of the level set functionφ. This is a special case of the higher order parametric trace finite element method presented in Section 3.2. In Section 3.2 we introduce the higher order parametric trace finite element spaces, list the fundamental inequalities that are needed for the analysis of TraceFEM and summarize the most important results of the error analysis for the parametric trace finite element method of the Laplace-Beltrami problem. We conclude this section with some numerical examples. The following Section 3.3 is devoted to the choice of the normal vector approximation to discretize the Laplace-Beltrami operator. We answer the question whether the normal vector of the exact surface is a better choice to discretize the Laplace-Beltrami operator in comparison to the usual normal of the discrete surface.
3.1 Basic concepts of the trace finite element method
In this section we describe the basic concepts of the trace finite element method using the example of the Laplace-Beltrami problem. We present two different discrete formulations, namely an unstabilized formulation and a stabilized formulation using the so-called normal gradient stabilization. The unstabilized method is the original trace finite element method introduced in [ORG09] for piecewise linear finite elements. The case of higher order finite elements was analyzed in [Reu15]. The normal gradient stabilized method was introduced later in [GLR18] for a higher order isoparametric TraceFEM (presented in Section 3.2) and in [BHLM18] for linear finite elements on embedded manifolds of arbitrary codimensions. Other stabilization methods for linear finite elements, however, were introduced earlier. Those will be discussed at the end of this section.
Let Ω ⊂ R3 be a polygonal domain that contains Γ. For the discretizations of the prob- lem we use the finite element spaces induced by the surface independent family of (closed) tetrahedral triangulations {Th}h>0 of the outer domain Ω. For ease of presentation we as- sume quasi-uniformity of the triangulation. Byh we denote the characteristic mesh size, i.e., h∼hT := diam(T), T ∈ Th. Here and from now on x ∼y means that there exist constants c1, c2 > 0, which are independent of h, such that c1x ≤y ≤ c2x. First we need a suitable approximation Γh of Γ. In this section we restrict to the case that Γh is the zero level of the continuous piecewise linear finite element nodal interpolationφh of φon Th, i.e.,
Γh:={x∈Ω|φh(x) = 0}.
Note that we do not need any knowledge of the signed distance functiond. Using basic finite element approximation theory, e.g., [EG13], one gets dist(Γh,Γ) ≤ ch2 with some constant c >0.
Remark 3.2. This choice of Γh is a special case of the higher order parametric trace finite element method (see Section 3.2). In [ORG09] and [Reu15] no specific construction of Γh was considered. Instead, the Lipschitz surface is only assumed to satisfy certain accuracy conditions concerning the distance and the normal vector.
Let ThΓ be the set of all tetrahedra T ∈ Th that have a nonempty intersection with Γh. The union of all tetrahedraT ∈ ThΓ is denoted by ΩΓh. Forx∈T ∈ ThΓ we define the normal approximation
nh(x) := ∇φh(x) k∇φh(x)k2
.
Forx∈Γh this normal approximation coincides with the unit normal on Γh (in the direction of φh >0).
Remark 3.3. Every ΓT := Γh∩T, T ∈ ThΓ is either a triangle or a quadrilateral. This surface triangulation of Γh is not necessarily regular, i.e., the elements ΓT may have very small internal angles and the size of the neighboring elements can vary strongly, cf. Figure 3.1. We further note that Γh is not a “surface triangulation” of Γ in the usual sense, i.e., a O(h2) approximation of Γ consisting of regular triangles with nodes lying on Γ, cf. Figure 3.1.
On the domain ΩΓh we define the space of continuous piecewise linear finite elements Vh :={vh∈C(ΩΓh)|vh|T ∈ P1(T) for all T ∈ ThΓ},
Figure 3.1: Example of a background meshTh and induced surface mesh in 3D (left) and 2D (right).
whereP1(T) is the space of polynomials of degree one onT. Using the spaceVhwe introduce the following trace spaces on Γh :
VhΓ :={vh|Γh |vh∈Vh}, Vh,0Γ :=n
vh ∈VhΓ
Z
Γh
vhdsh = 0o .
For the discrete problems we need a suitable extension of the data f to Γh denoted by fh. We assume that R
Γhfhdsh = 0 holds. Specific choices for fh are discussed in Remark 3.21.
We also introduce the following space on which the bilinear forms introduced below are well-defined:
Vreg,h:=n
v∈H1(ΩΓh)|v|Γh ∈H1(Γh)o
⊃Vh. (3.2)
For the definition of H1(Γh) we refer to Definition 3.11. We introduce the bilinear form ah(u, v) :=
Z
Γh
∇Γhu· ∇Γhvdsh, u, v∈Vreg,h, where
Ph(x) :=I−nh(x)nh(x)T, x∈ΩΓh, ∇Γhu(x) :=Ph(x)∇u(x), x∈Γh. The unstabilized discrete problem is as follows:
Discrete Problem 3.4 Find uh ∈Vh,0Γ such that
ah(uh, vh) = (fh, vh)L2(Γh) for all vh∈VhΓ. This discrete problem is well-posed (cf. [GLR18]).
Remark 3.5. Problem 3.4 has a unique solution in the trace spaceVhΓ. In the implementation of the method we use the bulk spaceVh. For everyvh ∈VhΓthere may be differentwh,w˜h ∈Vh withvh =wh|Γh = ˜wh|Γh. This is related to the fact that the traces of the outer finite element nodal basis functions only form a frame (in general not a basis) of the trace spaceVhΓ. Hence, there may be different representations of the same unique solution (on Γh) and the resulting stiffness matrix may be singular. The chosen solver determines a (non-unique) representation of the unique solution (on Γh).
Remark 3.6. In [OR10] the condition number of the diagonally scaled mass and stiffness matrix were investigated. Note that the mass matrix Mh has a one-dimensional kernel (the vector representation of the discrete level set function φh) and the stiffness matrixAh has a two-dimensional kernel (the vector representation of the discrete level set function φh and a constant function). The authors consider the effective condition number, which is the largest eigenvalue divided by the smallest non-zero eigenvalue. Certain numerical examples show an effective spectral condition number of O(h−2) for Mh and Ah for the three-dimensional case, and O(h−2) for Ah and O(h−3) for Mh for the two-dimensional case. Under certain assumptions on the distribution of the nodes near the surface they present an analysis for the bounds (up to an additional logarithmic term |lnh|) for the two-dimensional case. However, these assumptions are in practice not always fulfilled. In [GR11, Remark 13.2.7] a two- dimensional example was constructed where arbitrary small cuts of the surface Γh through a background triangle lead to an arbitrarily high effective condition number of the diagonally scaled stiffness matrix. This example can easily be extended to the three-dimensional case.
Hence, in certain cases the mass and stiffness matrices may be very ill-conditioned.
The discrete problem 3.4 has two disadvantages. We may have different representations of the same unique solution (see Remark 3.5) and the stiffness matrix may be very ill- conditioned, depending on how the surface cuts through the outer triangulation (see Remark 3.6). Both difficulties can be eliminated with a suitable stabilization. Forρ >0 we introduce the normal gradient volume stabilization:
sh(u, v) :=ρ Z
ΩΓh
(nh· ∇u)(nh· ∇v) dx, u, v∈Vreg,h. (3.3) This stabilization method was first introduced in [GLR18] for higher order finite elements and in [BHLM18] for linear finite elements. We discuss the choice of ρ in Remark 3.22. We define the stabilized bilinear form
Ah(u, v) :=ah(u, v) +sh(u, v), u, v∈Vreg,h.
In the stabilized problem we look for a solution in the outer finite element space Vh instead of the trace space Vh,0Γ . Hence, we define
Vh,0 :=n
vh ∈Vh
Z
Γh
vhdsh = 0o .
