Delft Aerospace Computational Science 0000 0000 0000 0000 0000 0000 0000 0000 0000 1111 1111 1111 1111 1111 1111 1111 1111 1111
000000000
000000000
000000000
000000000
000000000
111111111
111111111
111111111
111111111
111111111
ENGINEERING MECHANICSREPORT
November 2006Discontinuous Galerkin methods with
solenoidal elements for the Stokes problem
T.K. Heemskerk, E.H. van Brummelen and K.G. van der Zee
All rights reserved. This work may not be translated or copied in whole or in part without the written permission of the author.
TUD/LR/EM
Kluyverweg 1, 2629HS Delft, The Netherlands
P.O. Box 5058, 2600 GB Delft, The Netherlands
Phone: +31 15 278 5460
Fax: +31 15 261 1465
Stokes problem
T.K. Heemskerk, E.H. van Brummelen and K.G. van der Zee Delft University of Technology, Faculty of Aerospace Engineering
P.O. Box 5058, 2600 GB Delft, The Netherlands
November 30, 2006
ABSTRACT
The essential difficulty in the numerical solution of the incompressible Navier-Stokes (NS) equations is the coupling between the pressure and the velocity. The coupling enforces a constraint on the relation between the pressure and the velocity space. This constraint can be studied using a simplified form of the NS equations, viz. the Stokes equations. The present work examines the applicability of discontinuous Galerkin (DG) methods to the Stokes equations. Moreover, the application of DG-methods with locally solenoidal velocity fields is investigated. One of the main advantages of using DG-methods is their suitability to hp-adaptivity, as they can easily handle complex geometries and approximations that have polynomials of different degrees in different elements. We have restricted ourselves to the two most prominent DG methods for the Laplace operator, viz. the asymmetric Baumann-Oden (BO) method and the Symmetric Interior Penalty Galerkin (SIPG) method.
This work presents stability analysis for various elements with diverse finite-element spaces, together with analytically and numerically determined error convergence rates. Also a measure for the computational effort is given for each element. Emphasis is put on varying combinations of polynomial orders for the velocity and pressure spaces, showing the possibilities for adaptivity. Using the results of these tests, we see that applicability of DG-methods for the Stokes problem is very promising, especially with respect to stability, implementation of non-homogeneous boundary conditions and adaptivity.
Taking into account the error convergence, computational cost and the adaptivity properties, the solenoidal DG-element with discontinuous pressure on the skeleton is the most promising element. For this element the SIPG-method is the best method, because it is more stable and convergence is somewhat better than the BO-method, unless a penalty parameter is undesirable. In that case one should opt for the BO-method.
Keywords and Phrases: Stokes problem, finite element methods, solenoidal velocity field.
1. Introduction
The essential difficulty in the numerical solution of the incompressible Navier-Stokes (NS) equations is the coupling between the pressure and the velocity. The coupling enforces a constraint on the relation between the pressure and the velocity space, the so-called compatibility condition, see, for instance [6, 16, 20]. This is in contrast to compressible flow, where the continuity equation may be used as a transport equation for density and the energy equation as a transport equation for temperature. The pressure is then obtained from the density and the temperature using the equation of state. For incompressible flow the velocity is divergence free, that is, solenoidal. This is enforced by the pressure, which can therefore be seen as a Lagrange multiplier. However, it also remains a physical entity responsible for guiding the flow, although it is not a thermodynamic variable, since there is no equation of state.
To study the incompressible NS-equations various techniques can be used. However, to show how these techniques cope with the coupling between pressure and velocity, we can omit the non-linear part to get equations, which do possess the coupling but are easier to solve. These are the Stokes equations. In this report we will consider two general finite element techniques for the solution of the Stokes equations:
• Continuous Galkerin (CG) methods.
In this method the velocity is approximated by means of continuous shape functions. The standard way of solving the Stokes equations is based on conforming subspaces of the infinite-dimensional spaces that furnish the setting of the problem.
• Discontinuous Galerkin (DG) methods.
These are relatively new methods in which a solution is approximated by shape functions that are discontinuous over element edges. The inter-element continuity is then weakly enforced. The reduced continuity requirement endows these methods with two essential advantages. First, DG methods are well suited for adaptive-refinement strategies. Second, the reduced continu-ity requirement facilitates the construction of bases with particular properties, e.g., piecewise solenoidal bases.
Our main focus will be on the DG formulations of the Stokes problem. Our objective is to assess their advantages and disadvantages relative to the CG methods, especially with respect to the compatibility condition and its influence on the convergence of the solution.
For both the CG and the DG methods we can constrain the shape functions to solenoidal ones, for which the pressure term cancels out or becomes limited to the element edges, thereby bypassing or relaxing the compatibility condition. These methods have been neglected for a long time since it is hard to obtain a global solenoidal test space. However, with the rise in DG finite element formulations, interest in these methods has grown, see, for instance [5, 8, 9, 10, 11].
The outline of this work is as follows; In Chapter 2 the Stokes equations are derived, together with the weak formulations. Chapter 3 describes the theory of well-posedness of weak formulations, which contains the criteria to derive uniqueness and stability of a solution. Here the compatibility condition will be derived. In Chapter 4 CG approximations of the Stokes problem are given, together with compatibility criteria and convergence rates. Chapter 5 describes the theory of DG methods and its application to the Stokes problem, with an emphasis on DG methods with solenoidal velocity spaces. Analytical bounds are derived for the convergence of the solution, together with the criteria for which they hold. In Chapter 6 the analytically derived stability parameters and convergence bounds are verified numerically. Also the relative costs of the different finite-element methods are compared. Using these results, in Chapter 7 conclusions will be drawn on the suitability of DG methods for approximation of the solution of the Stokes problem. Also recommendations for future research will be given.
2. The Stokes equations
This work is concerned with finite-element methods for the Stokes equations for incompressible viscous fluid flow. The incompressible Stokes equations serve as a prototype for the incompressible Navier– Stokes equations. This section presents the Stokes equations. Moreover, two distinct variational formulations will be presented, viz., the mixed formulation and the constrained formulation. These variational formulations provide the basis for the finite-element methods in Chapters 4 and 5.
2.1 Derivation of the Stokes equations
equations. In this work we restrict ourselves to the steady Navier–Stokes equations:
−ν∆u + (u · ∇)u + ∇p = f in Ω, (2.1a)
divu = 0 in Ω, (2.1b)
u= g on ΓD, (2.1c)
−pn + νn · ∇u = h on ΓN, (2.1d)
where u = (u1, ..., un) : Ω → RN represents the velocity field, p : Ω → R the kinematic pressure, f= (f1, ..., fn) : Ω → RN a prescribed body force, g = (g1, ..., gn) : ΓD → RN a prescribed velocity field on the boundary ΓD and h = (h1, ..., hn) : ΓN → RN a prescribed traction on the boundary ΓN. The constant ν > 0 designates the viscosity. In the particular case that Γ = ΓD and, accordingly, ΓN =∅, g must satisfy the complementarity condition:
I Γ
n· g ds = 0.
The second term in the left member of equations (2.1a), the so-called convective term, renders the Navier–Stokes equations nonlinear. This nonlinearity forms a fundamental complication. Because the nonlinear term is inessential for our study of velocity/pressure compatibility for finite-element methods, we shall primarily restrict our considerations to the so-called Stokes equations, which are obtained from (2.1) by omitting the term (u· ∇)u:
−ν∆u + ∇p = f in Ω, (2.2a)
div u = 0 in Ω, (2.2b)
u= g on ΓD, (2.2c)
−pn + νn · ∇u = h on ΓN. (2.2d)
2.2 Classical weak formulations of the Stokes problem
2.2.1 Functional settingLet Ω be an open bounded domain in RN with Lipschitz boundary Γ. Let Hm(Ω) denote the Sobolev space of functions from Ω into R that reside in L2(Ω) together with their generalized partial derivatives up to and including order m; see also [14]. Equipped with the inner product
(u, v)m= X |α|≤m
Z Ω
Dαu Dαv dx, (2.3)
the space Hm(Ω) is a Hilbert space, i.e., Hm(Ω) is an inner product space that is complete with respect to the norm induced by the inner product; see, for instance, [4]. In equation (2.3), Dα = Dα1
1 Dα
2
2 · · · Dα
N
N denotes the generalized partial derivative and|α| = α1+ α2+· · · + αN. The norm and semi norm for u ∈ Hm(Ω) are denoted kuk
m and |u|m, respectively. The Sobolev spaces are
extended to vector-valued functions in the usual manner: [Hm(Ω)]N represents the space of functions from Ω into RN that are component wise in Hm(Ω). The space [Hm(Ω)]N is a Hilbert space under the inner product
(u, v)m= N X i=1
(ui, vi)m. (2.4)
For the L2inner product (·, ·)0 we also use the condensed notation (·, ·).
The dual space X′ of a Hilbert space X is the aggregate of all continuous linear functionals on X. The norm on X′ is defined byk · kX′ = supv
∈Xkvkh·,viX, where h·, ·i denotes the duality pairing between
2.2.2 Mixed formulation
We consider the Stokes equations (2.2) with homogeneous Dirichlet boundary conditions, i.e., Γ = ΓD and g = 0. To obtain a weak formulation, equations (2.2a)–(2.2b) are multiplied with test functions v∈ V and q ∈ Q and integrated over the domain Ω, where V and Q represent spaces of admissible test functions, to be specified later:
Z Ω−ν∆u · v dx + Z Ω∇p · v dx = Z Ω f· v dx ∀v ∈ V, (2.5a) Z Ω q∇ · u dx = 0 ∀q ∈ Q. (2.5b)
In accordance with the specification of Dirichlet boundary conditions, we require that functions in V vanish at the boundary Γ. Integration by parts then yields:
ν Z Ω∇u : ∇v dx − Z Ω p∇ · v dx = Z Ω f· v dx ∀v ∈ V, (2.6a) Z Ω q∇ · u dx = 0 ∀q ∈ Q. (2.6b)
This leads to the following weak mixed formulation for the Stokes problem:
Find (u, p)∈ V × Q such that :
a(u, v) + b(v, p) = (f, v) ∀ v ∈ V,
b(u, q) = 0 ∀ q ∈ Q,
(2.7)
where
a(u, v) = ν(∇u, ∇v), (2.8a)
b(v, p) =−(p, ∇ · v). (2.8b)
Formulation (2.7) constitutes a saddle-point or mixed problem, 2.2.3 Constrained formulation
Using the kernel space Vdiv0 ={v ∈ V : div v = 0} of the divergence operator, the mixed
formu-lation (2.7) can be transformed into a velocity-only form. The elements of Vdiv0 are referred to as
solenoidal vector fields. The kernel space enables the reformulation of equation (2.7) into the reduced variational problem:
Find u∈ Vdiv0 such that :
a(u, v) = (f, v) ∀v ∈ Vdiv0,
(2.9) The pressure can be calculated as a post-processing step. Equivalence between problems (2.7) and (2.9) follows from the following proposition [4]:
Proposition 2.1. Let u be a function in [H1
0(Ω)]N. The following statements are equivalent: (i) There is p in L2
0(Ω) such that the pair (u, p) solves (2.7). (ii) u solves (2.9).
The proof is outside the scope of this thesis.
3. Theory of well-posedness for weak formulations
3.1 Primal variational form
To obtain a definition for well-posedness, we consider the following abstract problem:
Find u∈ V such that :
a(u, v) =hf, vi ∀ v ∈ V, (3.1)
where
• V is a Hilbert space equipped with norm k · kV.
• a(·, ·) is a continuous bilinear form on V × V , i.e. a : V × V → R with a(u, v) ≤ γakukVkvkV. • f is a continuous linear form on V , i.e. f ∈ V′ = (V → R).
Well-posedness is understood as in [4]:
Definition 3.1 (Well-posedness). Problem (3.1) is said to be well-posed if it admits one and only one solution and if the following a priori estimate holds:
∃c > 0, ∀f ∈ V′, kukV ≤ ckfkV′ (3.2)
Well-posedness for problem (3.1) can be determined using the Lax-Milgram theorem or the gen-eralized Lax-Milgram theorem; see, for instance, [4, 18, 19]:
Theorem 3.2 (Lax-Milgram lemma). Let V be a Hilbert space, let a ∈ V × V → R be continuous and let f∈ V′. Assume that the bilinear form a is coercive, i.e.,
∃α > 0, ∀u ∈ V, a(u, u) ≥ αkuk2V, (3.3) then problem (3.1) is well-posed with a priori estimate
∀f ∈ V′, kukV ≤ 1
αkfkV′ (3.4)
Theorem 3.2 provides sufficient conditions for well-posedness. However, these conditions are not necessary. Necessary and sufficient conditions are provided by the premises of the generalized Lax-Milgram (or BNB [4]) condition:
Theorem 3.3 (Generalized Lax-Milgram lemma).
Let U be a Banach space and let V be a reflexive Banach space, and let a∈ U × V → R and let f ∈ V′. If and only if there exists an α > 0 such that
inf u∈U\{0}v∈V \{0}sup a(u, v) kukUkvkV ≥ α, (3.5a) sup u∈U\{0} a(u, v)≥ 0 ∀v ∈ V \ {0}, (3.5b)
then problem (3.1) is well-posed with a priori estimate ∀f ∈ V′, kuk
U≤ 1
αkfkV′. (3.6)
see [4].
3.2 Well-posedness of the mixed formulation
The constrained formulation of the Stokes problem conforms to (3.1). Hence, its well-posedness can be verified on the basis of Theorems 3.2 or 3.3. For the mixed formulation, further elaboration is required. Returning to the mixed formulation of the Stokes problem (2.7), we see that it is of the form:
Find (u, p)∈ V × Q such that :
a(u, v) + b(v, p) =hf, vi ∀ v ∈ V,
b(u, q) =hg, qi ∀ q ∈ Q.
(3.7)
where
• V and Q represent Hilbert spaces with corresponding norms k · kV andk · kQ.
• a(·, ·) and b(·, ·) are continuous bilinear forms on V × V and V × Q , i.e. a : V × V → R and b : V× Q → R, with a(u, v) ≤ γakukVkvkV and b(v, q)≤ γbkvkVkqkQ.
• f is a continuous linear form on V, i.e. f ∈ V′, and g is a continuous linear form on Q, i.e. g∈ Q′.
The bilinear form b defines an operator B : V→ Q. The kernel of this operator is defined as: ker B ={u ∈ V : b(u, q) = 0 ∀ q ∈ Q} .
Using this definition, necessary and sufficient conditions for well-posedness of (3.7) are provided by the following theorem:
Theorem 3.4 (Well-posedness mixed formulation). Using the framework above assume that:
inf
u∈ker B\{0}v∈ker B\{0}sup
a(u, v) kukVkvkV ≥ α > 0, (3.8a) sup u∈ker B\{0} a(u, v)≥ 0 ∀v ∈ ker B \ {0}, (3.8b)
and assume that the following compatibility condition, also known as the LBB condition, between the spaces V and Q holds:
inf
q∈Q\{0}v∈V\{0}sup
b(v, q)
kvkVkqkQ ≥ β > 0.
(3.9) Then, for every f ∈ V′ and g∈ Q′, problem (3.7) is well-posed and the following a-priori estimates hold kukV≤ 1 α kfkV′+α + γa β kgkQ′ (3.10) kpkQ≤ 1 β 1 +γa α kfkV′ +γa(α + γa) αβ kgkQ′ ; (3.11)
see, for instance, [4, 16, 18]:
Equations (3.7) can be recast into a primal form by introducing the product space W = V× Q and the operator c((u, p), (v, q)) = a(u, v) + b(v, p) + b(u, q):
Find (u, p)∈ W such that :
3.3 Well-posedness of a weak formulation in approximation spaces
Theorems 3.2, 3.3 and 3.4 are also suitable for establishing the well-posedness of finite-element ap-proximations of the Stokes problem, simply by replacing the infinite-dimensional Hilbert spaces in the above formulations by the finite-element approximation spaces. In Chapters 4 and 5, analytical results will be presented on the well-posedness of continuous and discontinuous Galerkin finite-element discretizations, respectively. Based on these results, bounds can be derived for the convergence of the solution in the approximation spaces.
It is to be noted that if the test and trial spaces in Theorems 3.3 and 3.4 are finite dimensional, then the second requirements, (3.5b) and (3.8b), are implied by the first part of the inf-sup condition, (3.5a) and (3.8a), respectively; see [4] Proposition 2.21 (iii).
4. Continuous Galerkin approximations
In this Chapter, we consider the Stokes problem in continuous finite-dimensional spaces. First, the mixed form of the Stokes problem (2.7) is examined in continuous approximation spaces. Next, we consider the constrained formulation of the Stokes problem (2.9) on a continuous solenoidal approxi-mation space.
4.1 Mixed CG Stokes formulation
4.1.1 Approximate finite-element formulationsWe partition the domain Ω into a number of triangular elements Ωe using standard triangula-tion techniques. The solutriangula-tion of the continuous saddle point problem (2.7) can be approximated by means of finite elements using finite-dimensional subspaces Vk1
h ⊂ V and Q
k2
h ⊂ Q of piece-wise polynomial functions, where Vk1
h = vh∈ V : vh|Ωe ∈ {Pk1(Ωe)}
N
∀Ωe∈ Th(Ω) and Qkh2 = {qh∈ Q : qh|Ωe ∈ Pk2(Ωe) ∀Ωe∈ Th(Ω)}. Here Pk(Ωe) denotes the space of polynomials of degree
≤ k on Ωeand h is a measure of the size of the elements. The spaces V and Q are identified as [H01]N and L2 0 respectively , where H01(Ω) =v ∈ H1(Ω) : v = 0 on ΓD , L20(Ω) = q∈ L2(Ω) : Z Ω q dx = 0 , together with the norms
k · kV=k · k1, k · kQ=k · k0.
In Appendix B.1 well-posedness of the Stokes problem in these spaces is shown for the continuum problem.
Using the finite-element approximation spaces Vk1
h and Q
k2
h , the discrete mixed problem becomes: Find (uh, ph)∈ Vk1 h × Q k2 h such that : a(uh, vh) + b(vh, ph) = (f, vh) ∀ vh∈ Vhk1, b(uh, qh) = 0 ∀ qh∈ Qkh2, (4.1) where a(u, v) = ν(∇u, ∇v), b(v, p) =−(p, ∇ · v).
Existence of a unique discrete solution can now be established on the basis of Theorem 3.4 with the spaces V, Q replaced by Vk1
h , Q
k2
000 000 111 111 000 000 111 111 000 000 111 111 000 000 111 111 000 000 111 111 00 00 11 11 00 00 11 11 00 00 11 11 00 00 11 11 000 000 111 111 000 000 111 111 000 000 111 111 000 000 111 111 000 000 000 111 111 111 000 000 111 111 000 000 111 111 00 00 11 11 00 00 11 11 000 000 111 111 Velocity
Velocity and continuous pressure Discontinuous pressure
Crouzeix-Raviart Taylor-Hood
Figure 1: 2D conforming triangular elements
4.1.2 Well-posedness of the finite-element formulation
Well-posedness of (4.1) generally depends on the finite-element type and triangulation. We consider two common types of elements, elements with continuous pressures and elements with discontinuous pressures, with best known members the Taylor-Hood element and the Crouzeix-Raviart elements, respectively. These elements are illustrated in Figure 1.
The two mentioned elements comply with the criteria set in Theorem 3.4; see, for instance, [4, 6, 18], although for Taylor-Hood there is a minimum number of internal edges, equal to space dimension N . It is to be mentioned, however, that many conforming-space elements are not divergence stable, e.g., all equal order (k1 = k2) continuous-velocity/continuous-pressure elements1. Certain of these finite elements can be stabilized by restricting the class of allowable grids, but we will not further explore this here.
Of all possible continuous velocity space elements we will examine the well-posedness of the con-tinuous pressure, concon-tinuous velocity elements, like the Taylor-Hood element. As said before, the Taylor-Hood element (with continuous second-order velocity space and continuous first-order pressure space) satisfies the divergence stability criterion, as well as all other Taylor-Hood type elements, hav-ing a continuous pressure space one polynomial degree lower than the velocity space, as long as the number of internal edges is equal to or exceeds the space dimension N . We will check this criterion numerically and compare the error convergence with other solution methods. The Crouzeix-Raviart element is mentioned for completeness. It will not be further considered here.
To examine well-posedness, compliance with Theorem 3.4 is checked for the spaces defined in the previous section. Upper bounds for the continuity constants can be derived in the following manner:
|a(uh, vh)| = ν R Ω∇uh∇vhdx , ≤ νqR Ω(∇uh)2dx q R Ω(∇vh)2dx, (Schwartz inequality), ≤ νkuhkVkvhkV, (kvkV=kvk1≥ |v|1), (4.2) and |b(v, p)| = |R Ωp∇ · v dx|, ≤qR Ωp2dx q R Ω(∇ · v) 2 dx, (Schwartz inequality), ≤ kpkQkvkV, (k(∇ · v)k0≤ kvkVand Q = L20). (4.3)
The Schwartz inequality can be found in, for instance, ([19], Appendix A). The above exposition reveals that γa = ν and γb = 1, which are indeed independent of h.
To show coercivity of a(., .) use a(u, u) = ν
Z
Ω∇u : ∇udx = ν|u| 2
H1≥ νCkuk2V ∀u ∈ [H1]N, (4.4)
where C is independent of mesh size h. The final inequality follows from the first Poincar´e inequality (Theorem A.25 from [18]). Hence, the coercivity constant is α = νC. Now since coercivity transfers to subspaces and because Vk1
h is a subspace of [H01]N, the inf-sup condition (3.8) is satisfied with constant αh independent of h.
A derivation of a lower bound for the inf-sup coefficient βhis much more intricate, and we will not further pursue it here; see, for instance, [1, 4, 18]. Inf-sup constants for particular choices of the mesh and the polynomial degree are presented in the numerical experiments section.
4.1.3 Convergence
To determine the error convergence and compare it with other elements, mixed or velocity-only, we are primarily concerned with the convergence of the velocity. Assuming well-posedness, Theorem 3.4 and the continuity conditions can be used to establish the following error estimate [4]:
ku − uhkV≤ 1 + γa αh 1 + γb βh inf vh∈Vhku − v hkV + γb αhqhinf∈Qhkp − q hkQ, (4.5)
where the second term can be dropped if ker(Bh)⊂ker(B). Equation (4.5) then reduces to ku − uhkV≤ 1 + γa αh 1 + γb βh inf vh∈Vhku − v hkV. (4.6)
To obtain a bound for the error between the finite-element approximation and the exact solution, it is necessary to know the accuracy of the interpolation. The approximation properties of Vk1
h and
Qk2
h are determined by standard approximation estimates (see [15]): inf
vh∈Vk1h
ku − vhks1 ≤ C0h
ρ1
kukr1 and inf
qh∈Qk2h kp − qhks2≤ C1h ρ2 kpkr2, (4.7) where u∈ [Hr1]N, p ∈ Hr2, ρ1= min(r1
− s1, k1+ 1− s1), ρ2= min(r2− s2, k2+ 1− s2) and C0and
C1are constants independent of h.
Combination of the error estimate (4.5) and the bound on the interpolation error (4.7) yields: ku − uhkV≤ 1 + γa αh 1 + γb βh C1hρ1 kukr1 + γb αhC2h ρ2 kpkr2, (4.8) where u∈ [Hr1]N, p ∈ Hr2, ρ1= min(r1
− 1, k1) and ρ2= min(r2, k2+ 1). Sincek · kV =k · k1 and k · kQ=k · k0we have s = 1 in the leftmost inequality in (4.7) and s = 0 in the rightmost.
Equations (4.2)–(4.4) convey that the continuity constants and the inf-sup coefficient αhare inde-pendent of h. If the inf-sup coefficient βhis also independent of h, the following convergence behaviour is implied: ku − uhkV≤ C1hρ1kukr1+ C2h ρ2 kpkr2, (4.9) where u ∈ [Hr1]N, p ∈ Hr2, ρ1 = min(r1
4.2 Constrained CG Stokes formulation
4.2.1 Approximate finite-element formulationAs in Section 4.1, we partition the domain Ω into triangular elements Ωe. A finite-element approx-imation of the continuous constrained problem (2.9) is obtained by replacing the space Vdiv0 by a
finite-dimensional, continuous and piecewise polynomial subspace Vk
div0:h ⊂ Vdiv0, where V
k
div0:h =
vh∈ Vdiv0 : vh|Ωe∈ {Pk(Ωe)}
N
∀Ωe∈ Th(Ω) . Here Pk(Ωe) denotes the space of polynomials of
de-gree≤ k on Ωeand h is a measure of the size of the elements. For more information on how to obtain Vkdiv0:h, see Appendix C. The space Vdiv0 is identified as
Vdiv0(Ω) =v ∈ [H
1
0(Ω)]N : div v = 0 in Ω , equipped with the norm
k · kV=k · k1.
Using the above finite-dimensional solenoidal spaces, the discrete variational problem becomes:
Find (uh)∈ Vk
div0:h such that:
a(uh, vh) = (f, vh) ∀ vh∈ Vkdiv0:h.
(4.10) The weak formulation now involves the velocity only. Therefore, this formulation is devoid of a compatibility condition.
4.2.2 Well-posedness of the finite-element formulation
Stability of problem (4.10) can be proved by establishing strong coercivity of the bilinear form on Vdiv0,
because coercivity transfers to subspaces. To show coercivity of the bilinear form in the constrained CG Stokes problem, we use (4.4), which shows that the inf-sup constant α is equal to νC with C independent of the mesh size h. The space Vk
div0:h being a subspace of [H
1]N, it immediately follows that the approximate problem (4.10) is well posed.
4.2.3 Convergence
To obtain a bound for the error between the finite-element approximation and the actual solution, it is again necessary to know the accuracy of the interpolation. The approximation properties of Vk
div0:h are given in [1]: inf vh∈Vk div0:h ku − vhk1≤ Chskuks+1, (4.11)
where 0 ≤ s ≤ k and C is independent of the mesh size. Using Theorem 3.2 and the continuity condition the following a-priori estimate can be derived [4]:
ku − uhkV≤ γa αh 12 inf vh∈Vku − v hkV (4.12)
From the interpolation estimate it follows that:
ku − uhkV≤ γa αh 12 Chs kuks+1, (4.13)
where 0≤ s ≤ k and C is independent of the mesh size. The continuity parameter γa is independent of h according to (4.2) and so is the coercivity coefficient αh. This leads to the following error estimate:
ku − uhkV≤ chskuks+1, (4.14)
As said before, the great difficulty is to obtain the space Vk
div0:h. We obtain this space by
em-ploying the method described in Appendix C. For small values of k, few elements and/or various grid configurations, this space can be empty, for instance for all k = 1 on the grid in Section 6.1. Then, the space of approximate solutions is obviously empty. The stability and convergence of this formulation are verified numerically in Chapter 6. The pressure can be recovered using the method in Appendix D.
4.3 Applying non-homogeneous Dirichlet boundary conditions to CG problem
In the previous sections we have investigated the Stokes problem with homogeneous Dirichlet boundary conditions in CG variational form. In this section we consider methods for enforcing non-homogeneous Dirichlet conditions.
There are two fundamental mechanisms to enforce Dirichlet boundary conditions. The first is strong enforcement or lifting (not to be confused with lifting operators for DG formulations) and the second is weak enforcement. For continuous Galerkin methods, strong enforcement of Dirichlet boundary conditions is standard practice. Here a demonstration of strong enforcement is given.
The inhomogeneous Stokes equations (2.2) with Dirichlet boundary conditions can be made ho-mogeneous in the following way [17]: Let w∈ [H1(Ω)]N be such that w = g on ∂Ω and div w = 0 in Ω. Set z = u− w in (2.2) to obtain the following homogeneous Stokes problem for z:
−ν∆z + ∇p = f − ν∆w in Ω, (4.15a)
div z = 0 in Ω, (4.15b)
z= 0 on ∂Ω. (4.15c)
A method for numerically obtaining a function w is presented in Section C.3 as a part of the process of determining the solenoidal shape-functions.
In solving the homogeneous equations we only consider approximations that are zero on Dirichlet boundaries. This reduces the number of degrees of freedom in comparison to weak enforcement of Dirichlet boundary conditions. Weak enforcement is accomplished by introduction of a boundary term into the variational statement. In the finite-element approximation, this translates into additional terms in elements adjacent to the Dirichlet boundary. Unlike strong enforcement, weak enforcement does not remove degrees of freedom associated with boundary nodes from the algebraic problem. Another dissimilarity is that a weak enforcement of the boundary conditions does not require the boundary condition to be met exactly. Instead, they are consistently approximated. The spatial resolution and order of convergence should be the same as for the solution in the interior, if applied properly [1].
5. Discontinuous Galerkin approximations
In this section, the Stokes problem in discontinuous finite-element spaces will be considered. First, a short introduction to DG-methods is given. Next, weak formulations and well-posedness criteria are derived for mixed, mixed solenoidal and fully constrained DG formulations.
5.1 DG introduction
In their quest to unify the existing DG-methods for second-order elliptic problems, Arnold et al.[2] identify two main ideas behind the development of the DG-methods. These are the methods inspired by the original IP methods and the methods inspired by the finite-volume techniques for hyperbolic problems. The first group has originally been developed in the primal-formulation form, while the second group has been developed in flux form. As a consequence, not all formulations developed in flux form can be transformed into a primal formulation without having to resort to an artificial construct called the lifting operator. Obtaining this lifting operator is in itself a variational problem, thus adding to the computational cost of the method. Still it can be rewarding to use these methods since they can be shown to yield favorable convergence properties. However, methods with lifting operators will not be considered here. For more information, see [8].
Various finite-dimensional spaces can be defined for the discontinuous Galerkin approximation of the Stokes equations. Here a division is made between methods with standard velocity DG spaces and methods using broken divergence free velocity spaces, either in mixed form with a pressure term or decoupled with a jump penalty.
5.2 Notations for broken spaces
To derive the discontinuous Galerkin formulation, we integrate over a sub-domain Ke of Ω, instead of over the complete domain Ω as is done for the continuous Galerkin methods. Let us first define notations and conventions for broken spaces, cf [14]. Let P = {Ph(Ω)}h>0 be a family of regular partitions of domain Ω⊂ RN into N (
Ph) sub-domains Ke, such that forPh∈ P,
¯ Ω = N (Ph) [ e=1 ¯
Ke, and Ke∩ Kf =∅ for e 6= f.
Let us denote by ξI
h the set of inter-element edges pertaining to Ph, by ξDh the set of Dirichlet boundary edges, and by ξN
h the set of Neumann boundary edges. On every edge in ξhI, we define ¯
n= ne on (∂Ke∩ ∂Kf)⊂ ξIh for indices e and f such that e < f , where e and f are the element numbers. Furthermore, we define on every ξI
h the jump [[q]] as (qe− qf) and the average{q} as 1
2(qe+ qf). On boundary edges this is set as [[q]] ={q} = qe. Γ is the set of element edges where every element has its own edge, so all internal edges are counted twice, which is in contrast to ξ. We write (p, φ)Ph = ΣKe∈Ph(p, φ)Ke, where (p, φ)Ke represents the L
2 inner product on the element Ke. If Fh is a subset of ξh we use the notation (f, φ)Fh = Σe∈Fh(f, φ)e for the L
2 product on all the edges in Fh.
On the spacePhwe define the normskvkm,Ph =
P K∈Phkvkm,K andkvkξ = P e∈ξkh 1 2 evke, where
heis the length of edge e.
5.3 Mixed DG Stokes formulation
5.3.1 Approximate finite-element formulationsIn this section we derive the weak formulation for the Stokes problem on broken velocity spaces and continuous pressure space. We define the following finite-dimensional spaces onPh:
Vk1 h =vh∈ [L2(Ω)]N : vh|Ke ∈ [H 1(Ke)]N ∩ {Pk1(Ke) N } ∀ Ke∈ Ph(Ω) , Qk2 h = qh∈ H1(Ω) : qh|Ω∈ Pk2(Ω) : Z Ω q dx = 0 ⇐⇒ ΓN =∅ . Here Pl(D) denotes the space of polynomials of degree≤ l on D.
To obtain a weak formulation we start by multiplying equations (2.2) by test functions vh∈ Vk1
h and qh ∈ Qk2
defined on element Ke: ν Z Ke ∇u : ∇vhdx− ν Z ∂Ke vh∂u ∂nds + Z Ke vh· ∇p dx = Z Ke f· vhdx, (5.1a) Z Ke qh∇ · u dx = 0. (5.1b)
Summation over all sub-domains Ke and application of Green’s formula on the second equation then yields: ν (∇u, ∇vh)Ph − ν vh,∂u ∂n Γ + (vh,∇p)Ph= (f, vh)Ph, (5.2a) (u,∇qh)Ph− (qh, u· n)Γ= 0, (5.2b)
where it should be noted that there is no connectivity in velocity between the elements.
On account of the continuity of the velocity field of the actual solution, the functionals in (5.2) can be augmented with terms that vanish for continuous velocity fields u. Such terms are referred to as consistency-preserving augmentations. Such terms do not interfere with the consistency of the formulation. However, they can improve the stability by enhancing the coercivity; see, for instance [19]. For a particular class of these stabilizing terms, equations (5.2) can be cast into the following weak formulation. Find (uh, ph)∈ Vk1 h × Q k2 h such that : ν(∇uh,∇vh)Ph± ν([[uh]],{∇vh· ¯n})ξ −ν({∇uh· ¯n}, [[vh]])ξ+ νσh([[uh]], [[vh]])ξ+ (vh,∇ph)Ph = (f, vh)Ph± ν(g, ∇vh· n)ξD+ ν σ h(g, vh)ξD ∀ vh∈ Vk1 h , (uh,∇qh)Ph= (qh, g· n)ξ D ∀ qh∈ Qk2 h . (5.3)
Two options for the values of ± and σ are shown in Table 1. The name by which the choice of parameters is often referred to, is also given.
We can write the weak formulation as: Find (uh, ph)∈ Vk1 h × Q k2 h such that : a(uh, vh) + b(vh, ph) = (f, vh)Ph± ν(g, ∇vh· n)ξD + ν σ h(g, vh)ξD ∀ vh∈ V k1 h , b(uh, qh) = (qh, g· n)ξD ∀ qh∈ Qk2 h , (5.4) where
a(uh, vh) =ν(∇uh,∇vh)Ph± ν([[uh]],{∇vh· ¯n})ξ− ν({∇uh· ¯n}, [[vh]])ξ
+ νσ
h([[uh]], [[vh]])ξ, b(vh, ph) = (vh,∇ph)Ph.
We will restrict ourselves to the two most prominent DG methods for the Laplace operator, viz. the asymmetric Baumann-Oden (BO) method and the Symmetric Interior Penalty Galerkin (SIPG) method for the given velocity and pressure spaces.
5.3.2 Well-posedness of the formulation in finite dimensions
To ensure existence and uniqueness of a solution we have to show that the formulation complies with the premises of Theorem 3.4. To do so, we first have to specify the norm for the space Vk1
Method ± σ
SIPG − σ
Baumann-Oden + 0
Table 1: The± sign and parameters for several DG formulations
Baumann and Oden DG method
The following bilinear form is the Baumann Oden discontinuous Galerkin functional for the Laplace operator [14]:
a(uh, vh) =ν(∇uh,∇vh)Ph+ ν([[uh]],{∇vh· ¯n})ξ− ν({∇uh· ¯n}, [[vh]])ξ, (5.5)
where ν is incorporated for consistency with equation 5.4.
Boundedness. For the boundedness of the bilinear form, the following derivation can be made, using the Schwartz inequality:
|a(uh, vh)| ≤ + ν|(∇uh,∇vh)Ph| + ν|([[uh]],{∇vh· ¯n})ξ| + ν|({∇uh· ¯n}, [[vh]])ξ| , ≤ + ν X K∈Ph s Z K (∇uh)2dx Z K (∇vh)2dx + νX e∈ξ s Z e [[uh]]2ds Z e{∇v h· ¯n}2ds + s Z e [[vh]]2ds Z e{∇u h· ¯n}2ds ! .
The discrete Schwartz inequality then gives:
|a(uh, vh)| ≤ + ν s X K∈Ph Z K (∇uh)2dx s X K∈Ph Z K (∇vh)2dx + ν s X e∈ξ 1 he Z e [[uh]]2ds s X e∈ξ he Z e{∇v h· ¯n}2ds + ν s X e∈ξ he Z e{∇u h· ¯n}2ds s X e∈ξ 1 he Z e [[vh]]2ds.
The Schwarz and discrete Schwartz inequalities can be found in, for example, [19]. If we apply the norm below on the velocity space:
kvhk2Vbo = X K∈Ph Z K (∇vh)2dx +X e∈ξ 1 he Z e [[vh]]2ds +X e∈ξ he Z e{∇v h· ¯n}2ds =k∇vhk2Ph+kh −1 e [[vh]]k2ξ+k{∇vh· n}k2ξ, (5.6) we obtain the following bound for the bilinear form:
|a(uh, vh)| ≤ νkuhkVbokvhkVbo. (5.7)
Coercivity. On account of the asymmetry of the BO formulation, well-posedness of the formulation relies on weak coercivity, rather than coercivity. The exposition in [19] reveals that the BO formulation is not coercive on the broken continuum space. However, in [13, 14] it is shown that the variational formulation is stable for finite-element approximations, provided that the polynomial degree of the approximation space exceeds one. For such spaces, the inf-sup parameter is strictly positive, although it does depend on the shape of the elements and the polynomial order. We expect the same behaviour for the BO formulation of the Stokes problem.
Symmetric Interior Penalty Galerkin
The following bilinear operator is the Symmetric Interior Penalty Galerkin (SIPG) operator for elliptic equations [3]: a(uh, vh) =ν(∇uh,∇vh)P − ν([[uh]],{∇vh· ¯n})ξ− ν({∇uh· ¯n}, [[vh]])ξ + νX e∈ξ σe he Z e [[uh]][[vh]]ds. (5.8)
An appropriate norm on the velocity space for the SIPG formulation is: kvhk2Vsipg= X K∈Ph Z K (∇vh)2dx +X e∈ξ σe he Z e [[vh]]2ds =k∇vhk2Ph+k √σe he [[vh]]k 2 ξ, (5.9)
where ν is incorporated for consistency with equation 5.4. In [3] it is proven that for two-dimensional Poisson problems continuity and the inf-sup condition (with parameters γa and αh:sipg independent of h) are satisfied for
σe> 6(k1)(k1+ 1) cot(θPh), (5.10)
where θPh is the smallest angle of the triangles in the triangulation. It is important to note that
the derivation in [3] holds for a primal elliptic formulation and that we want to satisfy the inf-sup condition for the mixed formulation. Therefore, the value of the penalty σe can be different from the one given above for the Poisson problem.
Stability ofb(vh, qh)
In this section the continuity and the inf-sup constant of the following bilinear functional is investi-gated:
b(vh, ph) = (vh,∇ph)Ph. (5.11)
Boundedness. The functional in (5.11) can be bounded as follows: |b(vh, ph)| ≤ X K∈Ph Z ∂K phvh· nds − Z K ph∇ · vhdx ≤ X e∈ξ Z e ph[[vh· ¯n]]ds − X K∈Ph Z K ph∇ · vhdx ≤X e∈ξ s he Z e p2 hds s 1 he Z e [[vh· n]]2ds + X K∈Ph s Z K p2 hdx s Z K (∇ · vh)2dx ≤ s X e∈ξ he Z e p2 hds s X e∈ξ 1 he Z e [[vh· n]]2ds + s X K∈Ph Z K p2 hdx s X K∈Ph Z K (∇ · vh)2dx ≤CkphkQkvhkV, (5.12)
with C independent of h. The first, third and fourth inequalities follow from the divergence theorem, the Schwartz inequality and the discrete Schwartz inequality, respectively.
The final step in (5.12) should hold for both the norm associated with the Baumann Oden method (5.6) and the norm of the SIPG method (5.9). Because pressure is defined up to a constant, the following norm is chosen:
kphk2Q\R= min
c∈R kph− ck
2
ξ+kph− ck2Ph . (5.13)
With the above choice of the normk · kQ, the operator b(vh, ph) is bounded for both the Baumann Oden and the SIPG method with their corresponding normsk · kV.
Inf-sup constant. A precise derivation of the inf-sup constant for equation (5.11) is intricate, and we will not further pursue it here. However, let us note that the space of discontinuous velocity fields encompasses the space of continuous velocity fields used in Section 4.1. Loosely speaking, the compatibility problem can be associated to the fact that the pressure space is so large that the velocity field is over-constrained. Therefore, it is anticipated that extension of the velocity space from continuous to discontinuous functions will in have a positive effect on the inf-sup constant.
5.3.3 Convergence
The approximation properties of Vk1
h will be examined by means of standard local-approximation estimates (see [14]). Let u ∈ [Hr(Ke)]N; there exists a constant C depending on r and on the the shape of Kebut independent of u and he such that for any 0≤ s ≤ r the following estimate holds:
inf
vh∈Vhk1
ku − vhks,Ke ≤ Ch
ρ
ekukr, (5.14)
where ρ = min(k1+ 1− s, r − s). For the space Qk2
h the same approximation properties as given in the rightmost inequality in (4.7) hold:
inf
qh∈Qk2h
kp − qhks≤ Chρkpkr,
where p∈ Hr2, ρ = min(k2+ 1
− s, r − s) and C is a constant independent of h. According to [13], we have the following trace inequality:
kvk2
Using this inequality, element-wise approximations for parts of the norms can be made: kh−1e [[v]]k2ξ ≤ C X K∈Ph h−1K kvkK h−1K kvkK+kvk1,K , k{∇v · ¯n}k2ξ ≤ C X K∈Ph hKk∇vkK h−1K k∇vkK+k∇vk1,K , ≤ C X K∈Ph kvk1,K(kvk1,K+ hKkvk2,K) ,
Now for the BO method the following holds for the velocity:
ku − uhk2Vbo≤ C ku − uhk21,Ph+ X K∈Ph h−1K ku − uhkK h−1K ku − uhkK+ku − uhk1,K + X K∈Ph ku − uhk1,K(ku − uhk1,K+ hKku − uhk2,K) ! ,
and for the SIPG method:
ku − uhk2Vsipg ≤ C ku − uhk21,Ph+ X K∈Ph h−1K ku − uhkK hK−1ku − uhkK+ku − uhk1,K ! .
For the pressure space we obtain kp − phk2Q≤ C X K∈Ph hKkp − phk2∂K+kp − phk2K ≤ C X K∈Ph kp − phk2K+ hKkp − phkKkp − phk1,K .
In the above equations, C represents various constants, possibly distinct from one instantiation to the next. In all instances, however, C is independent of the element size hK. If all elements are equally sized, hK = h.
Assuming that the formulation satisfies Theorem 3.4, we can again use equation (4.5). Now using the approximation properties described in equations (4.7) and (5.14), we get the following estimate for the convergence of the error:
ku − uhkV≤ 1 + γa αh 1 + γb βh inf vh∈Vhku − v hkV + γb αhqhinf∈Qhkp − q hkQ, (5.16) ≤C 1 + γa αh 1 + γb βh hρ1 X K∈Ph kuk2r1,K !12 + γb αhh ρ2 X K∈Ph kpk2r2,K !12 , (5.17)
5.4 Mixed solenoidal DG Stokes formulation
5.4.1 Approximate finite-element space formulationsIn this section we use broken velocity spaces to obtain an element-wise divergence-free velocity field. The pressure can be taken continuous or discontinuous. We define the following finite-dimensional spaces onPh: Vk1 h =vh∈ L 2(Ω)N : vh |Ke ∈ [H 1(Ke)]N ∩ {Pk1(Ke) N : div vh= 0 } ∀ Ke∈ Ph(Ω) . (5.18) For the pressure space we use one of the following two spaces:
Qk2 h = qh∈ L2(ξ) : qh|ξ ∈ Pk2(ξ) : Z ξ qhds = 0 , (5.19) or Qk2 h = qh∈ L2(ξ) : qh|e∈ {Pk2(e)} ∀ e ∈ ξ : Z ξ qhds = 0 . (5.20)
The difference between the two spaces is that for the first one the pressure is continuous on the skeleton and for the second the pressure is only edgewise continuous.
Equations (5.4) are modified to account for the element wise divergence-free velocity by applying Green’s formula to b(vh, ph): b(vh, ph) = (vh,∇ph)Ph, = X K∈Ph Z K vh∇phdx, = X K∈Ph Z K∇ · p hvhdx− X K∈Ph Z K ph∇ · vhdx, = X K∈Ph Z ∂K phvh· n ds, =X e∈ξ Z e ph[[vh· ¯n]] ds,
where the divergence term is dropped because of the solenoidal velocity space. Hence, the weak problem becomes: Find (uh, ph)∈ Vk1 h × Q k2 h such that : a(uh, vh) + b(vh, ph) = (f, vh)Ph± ν(g, ∇vh· n)ξD + ν σ h(g, vh)ξD ∀ vh∈ Vk1 h , b(uh, qh) = (qh, g· n)ξD ∀ qh∈ Qk2 h , (5.21) where
a(uh, vh) =ν(∇uh,∇vh)Ph± ν([[uh]],{∇vh· ¯n})ξ− ν({∇uh· ¯n}, [[vh]])ξ
+ νσ
h([[uh]], [[vh]])ξ, b(vh, ph) =(ph, [[vh· ¯n]])ξ.
5.4.2 Well-posedness of the formulation in finite dimensions
For existence and uniqueness of a solution we again have to show that the formulation complies with Theorem 3.4. Since the bilinear form a(uh, vh) is the same as in Section 5.3 for both the BO method and the SIPG method, the same norms as derived above will be used for space Vk1
h . Hence, only the stability of b(vh, qh) remains to be investigated.
Stability ofb(vh, qh)
In this section the continuity and the inf-sup constant of the following bilinear form are investigated: b(vh, ph) = (ph, [[vh· ¯n]])ξ = X e∈ξ Z e ph[[vh· ¯n]]ds, (5.22)
where the inf-sup constant is specified by β in equation (3.9).
Boundedness. Using the Schwartz and discrete Schwartz inequalities the functional in equa-tion (5.22) can be bounded as follows:
|b(vh, ph)| = X e∈ξ Z e ph[[vh· ¯n]]ds ≤X e∈ξ s he Z e p2 hds s 1 he Z e [[vh· ¯n]]2ds ≤ s X e∈ξ he Z e p2 hds s X e∈ξ 1 he Z e [[vh· ¯n]]2ds ≤ s X e∈ξ he Z e p2 hds s X e∈ξ 1 he Z e [[vh]]2ds ≤ kphkQkvhkV (5.23)
This derivation is valid for both the norm of the BO method (5.6) and the norm of the SIPG method (5.9), since for the SIPG method σe is always greater than one. Recalling that the pres-sure is specified up to a constant, the norm that bounds the prespres-sure term ph is:
kphkQ\R= min c∈R s X e∈ξ he Z e (ph− c)2ds = min c∈Rkph− ckξ, (5.24) so that b(vh, ph) is bounded for both the BO and the SIPG method with their corresponding norms in V. This norm is the same for either of the spaces defined in equations (5.19) and (5.20).
Inf-sup constant. According to equation (3.9) the inf-sup condition for b(vh, qh) is
inf qh∈Q\R sup vh∈V b(vh, qh) kvhkVkqhkQ ≥ β > 0, which can be rewritten as
sup
vh∈V
b(vh, qh)
kvhkV ≥ βkqhk
Q, (5.25)
To prove the inf-sup condition for the functional b(vh, qh), we will use a result from [8], derived for spaces defined in equations (5.18) and (5.20).
Theorem 5.1. For any qh∈ Qk
h∃ vh∈ Vkh such that X K∈Ph kvhk20,K+ h2Kk∇vhk20,K≤ C|||qh|||2Q, (5.26) with |||·|||2 Q = P K∈Ph|∂K| −1kq
hk20,∂K/R and there is a constant C independent of the mesh size such
that
b(vh, qh)≥ C|||qh|||2Q ∀vh∈ Vkh, (5.27) Although this Theorem has been derived for the space Qk
h in equation (5.20), it is evident from the derivation of Theorem 5.1 in [8] that it also holds for the space Qk
h in equation (5.19).
Also given in [8] is the multiplicative trace inequality, which asserts that for each edge e of an element K∈ Phit holds that:
h−1e kvhk20,e≤ Ckvhk0,K h−2K kvhk0,K+ h−1K k∇vhk0,K ,
≤ Ch−2K kvhk20,K+ Ck∇vhk20,K. (5.28) For the SIPG method the relation between the normkvhkVas in equation (5.9) and the norm|||qh|||Q can be derived using equations (5.26) and (5.28):
kvhk2Vsipg= X K∈Ph Z K (∇vh)2dx +X e∈ξ σe he Z e [[vh]]2ds, ≤ X K∈Ph k∇vhk20,K+ X K∈Ph C hK Z ∂K v2hds, ≤ C X K∈Ph 1 h2 Kkv hk20,K+k∇vhk20,K , ≤ h2C min|||q h|||2Q. (5.29)
For the Baumann Oden method the derivation is more complex due to the occurrence of the term P
e∈ξhe
R
e{∇vh· ¯n}2ds in the norm. However, according to [13], the following inequality holds: X e∈ξ he Z e{∇v h· ¯n}2ds≤ C X K∈Ph Z K (∇vh)2dx ∀vh∈ Vkh1, (5.30)
where the constant C depends on the polynomial degree k, but not on the mesh-size h. Using this inequality, a bound similar to that in equation (5.29) can be derived.
According to [8] the following holds:
with C independent of mesh size and hmin= minK∈PhhK. Let us note that this bound is not sharp,
and that the inf-sup constant of b(vh, qh) can in fact be independent of h. Anticipating the results in Section 6.3, we remark that the numerical experiments indeed indicate that the inf-sup constant of b(vh, qh) is independent of h.
5.4.3 Convergence
The approximation properties of Vk1
h given in equation (5.18) will be estimated using the following lo-cal approximation estimates, see Theorem C.1. Let u∈ [Sr]N, where [Sr]N =u ∈ [Hr(Ke)]N : div u = 0 ; there exists a constant C depending on r and on the shape of Kebut independent of u and hesuch that for any 0≤ s ≤ r the following holds:
inf vh∈Vhk1 ku − vhks,Ke ≤ Ch ρ ekuhkr, (5.33) where ρ = min(k1+ 1− s, r − s).
Using the same inequalities as in Section 5.3.3 for the BO and the SIPG method, we obtain the following inequality for the BO method:
ku − vhk2Vbo ≤ C ku − vhk21,Ph+ X K∈Ph h−1ku − vhkK h−1ku − vhkK+ku − vhk1,K + X K∈Ph ku − vhk1,K(ku − vhk1,K+ hku − vhk2,K) ! ,
and for the SIPG method:
ku − vhk2Vsipg ≤ C ku − vhk21,Ph+ X K∈Ph h−1ku − vhkK h−1ku − vhkK+ku − vhk1,K ! ,
with C distinct but independent of h in both equations.
We have defined two pressure spaces. However, the elements of both are bounded in the norm in (5.24). Assuming the actual pressure to be continuous on the skeleton, both spaces yield the same standard approximation estimates
inf
qh∈Qk2h
kp − qhks,e≤ Chρkpkr,e (5.34)
where p∈ Hr(e), ρ = min(k2+ 1
− s, r − s) and C is independent of h.
Using the approximation properties in equations (5.33) and (5.34), we get the following estimate for the convergence of the error in either the BO or the SIPG norm:
ku − uhkV≤ 1 + γa αh 1 + γb βh inf vh∈Vhku − v hkV + γb αhqhinf∈Qhkp − q hkQ, (5.35) ≤C 1 + γa αh 1 + γb βh hρ1 X K∈Ph kuk2r1,K !12 + γb αhh ρ2 X e∈ξ kpk2 r2,e 1 2 , (5.36)
5.5 Decoupled solenoidal DG Stokes problem
5.5.1 Approximate finite-element formulationsIn this section we again use broken velocity spaces to obtain an element-wise divergence-free velocity field. We use the same velocity space as in the previous section:
Vkdiv0 =vh∈ L
2(Ω)N : vh
|Ke∈ [H
1(Ke)]N
∩ {Pk(Ke)N : div vh= 0} ∀ Ke∈ Ph(Ω) . (5.37)
To circumvent the compatibility condition, we should reformulate equation (5.21) such that it depends on uh and vh only. However, to achieve this, we have to eliminate the (p, [[v· ¯n]])ξ term. Let us for the moment assume that this term can be ignored. Then the primal formulation becomes
Find u∈ Vk
div0 such that :
a(u, v) = (f, v)Ph± ν(g, ∇v · n)ξD + ν σ h(g, v)ξD ∀ v ∈ Vk1 div0, (5.38) where
a(u, v) =ν(∇u, ∇v)Ph± ν([[u]], {∇v · ¯n})ξ− ν({∇u · ¯n}, [[v]])ξ
+ νσ
h([[u]], [[v]])ξ.
Again the various options for the values of ± and σ are shown in Table 1. It is to be mentioned that formulation (5.38) is inconsistent, independent of the choice of the sign of± and the parameter σ, on account of the omission of (p, [[v· ¯n]])ξ. However, adding a term that assigns a very large penalty to the normal components of the jumps in velocity in principle enables equation (5.38) to be used to approximate the velocity in the Stokes problem. Both the BO and the SIPG-method can be endowed with a super penalty on the normal component of the jumps. The penalty term that we add to a(u, v) is:
νθ
h([[u· ¯n]], [[v · ¯n]])ξ, (5.39)
which adds the following term to the right-hand side of equation (5.38): νθ
h(g· n, v · n)ξD. (5.40)
Clearly, formulation (5.38) is devoid of a compatibility condition, as the formulation only involves the velocity. In the numerical experiments, we will set θ = 106. This renders the discontinuity in the normal components of the velocity field over the element edges sufficiently small to be ignored for most choices of polynomial order and mesh size, without degenerating the continuity of the bilinear form (and, accordingly, the condition number of the system of algebraic equations to be solved) too much; see Section 6.5.
Let us remark that a formulation similar to the above is described in [7], where an SIPG method with a super penalty on the normal component of the jumps is used.
5.5.2 Well-posedness of the formulation in finite dimensions Baumann Oden method with extra jump penalty
The weak formulation of the Baumann-Oden method with the extra jump penalty is: (
Find u∈ Vk
div0 such that :
a(u, v) = (f, v)Ph+ ν(g,∇v · n)ξD+ ν θ h(g· n, v · n)ξD ∀ v ∈ V k1 div0, (5.41) where
a(u, v) =ν(∇u, ∇v)Ph+ ν([[u]],{∇v · ¯n})ξ− ν({∇u · ¯n}, [[v]])ξ
+ νθ
The norm that bounds a(·, ·) is the same as the BO norm: kvk2
Vbo =k∇vk2Ph+kh
−1
e [[v]]k2ξ+k{∇v · ¯n}k2ξ, (5.42) because the extra penalty term h−1([[v· ¯n]], [[v · ¯n]])
ξ can be bounded bykh−1e [[v]]k2ξ. For large θ, the continuity constant will be significantly larger than for the velocity part of the original BO bilinear form. However, the continuity constant is still independent of mesh size h if θ is independent of h.
To check the coercivity of the operator we make the following derivation: a(vh, vh) = ν|vh|21,Ph+ νσk[[h
−1
e v· ¯n]]kξ. (5.43)
Using equation (5.30), it follows that thatkvhk2Vbo ≤ Ck∇vhk 2 Ph+kh −1 e [[vh]]k2ξ. If we now assume that σkh−1 e [[vh· ¯n]]kξ ≥ kh−1e [[vh]]kξ then: a(vh, vh)≥ ν|vh|21,Ph+ νk[[vh]]k 2 ξ ≥ Ckvhk2Vbo (5.44)
with C independent of h and the problem is well posed. An indication of the validity of the assumption will be given by the numerical results in Section 6.3.
Symmetric Interior Penalty Galerkin method with extra jump penalty The weak formulation of the SIPG method with extra jump penalty is:
(
Find u∈ Vk
div0 such that :
a(u, v) = (f, v)Ph− ν(g, ∇v · n)ξD+ ν σ h(g, v)ξD+ ν θ h(g· n, v · n)ξD ∀ v ∈ Vk1 div0, (5.45) where
a(u, v) =ν(∇u, ∇v)Ph− ν([[u]], {∇v · ¯n})ξ− ν({∇u · ¯n}, [[v]])ξ
+ νσ
h([[u]], [[v]])ξ+ ν θ
h([[u· ¯n]], [[v · ¯n]])ξ. The norm that bounds a(·, ·) is the same as the SIPG norm:
kvhk2Vsipg=k∇vhk2Ph+k
√σe h [[vh]]k
2
ξ. (5.46)
because the extra penalty term h−1([[v· ¯n]], [[v · ¯n]])
ξ can be bounded by k
√σ
he[[v]]k
2
ξ. As for the BO method, the continuity constant will increase dramatically for a(u, v) through to the influence of θ. However, it remains independent of h.
In Section 5.3.2 compliance with the inf-sup condition was established for the SIPG method, for a certain minimum penalty parameter σ. In that case the SIPG method is coercive. The additional jump penalty will have a positive effect on the inf-sup parameter, because of its positivity. Hence, we will have a well-posed problem for this method.
5.5.3 Convergence
Well-posedness of the decoupled solenoidal DG-formulation (5.38) does not immediately imply con-vergence, on account of the consistency error. The consistency error is defined as:
Rh(u, p) = sup vh∈V |f(vh)− a(u,vh)| kvhkV = sup vh∈V |(p, vh· ¯n)ξ| kvhkV (5.47) Now according to [4] the following holds:
0 1 0
1
x
y
Figure 2: Illustration of the mesh layout for h = 16
which shows the influence of increased coercivity on the error-convergence. Using the same approxi-mation estimates and norm inequalities as in Section 5.4.3, we get:
ku − uhkV≤C 1 + γa αh hρ X K∈Ph kuk2r,K !12 + 1 αhRh(u, p), (5.49)
where the constant C is independent of h and where ρ = min(k, r− 1). Again, we will evaluate this behaviour in Chapter 6, where the consistency term Rh(u, p) is omitted to facilitate comparissons with other elements. This is a premise that is probably incorrect.
The pressure field can be derived as a post-processing step, as explained in Appendix D.
6. Numerical experiments
In this section we conduct numerical experiments, to verify the theoretical results and to establish the inf-sup parameters for those cases in which an analytical derivation proved too complicated. We will compare the analytical and numerical results for the various formulations and comment on their similarities and differences.
6.1 Mesh layout
The calculations were performed on sequences of regular meshes, composed of isosceles right triangles. The mesh size h is taken as the length of the equal sides. For the domain Ω, we take the unit square. Figure 2 shows the mesh for h = 1
6.
The number of elements is equal to 2
h2. We opted for a regular grid to ensure that we have a
6.2 Finite-elements spaces
In the previous sections various finite elements for the Stokes equations have been described and evaluated analytically. These elements shown in Figure 3. The following elements can be identified by their velocity and pressure spaces:
• Element 1. Continuous velocity and pressure spaces with homogeneous Dirichlet velocity bound-ary conditions: Vk1 0 ≡v ∈ [H1(Ω)]N : v = 0 on ΓD : v|K ∈ {Pk1(K)} N ∀K ∈ T (Ω) ||v||V=||v||H1 Qk2 0 ≡ q∈ L2(Ω) : Z Ω q dx = 0 : q|K ∈ Pk2(K)∀K ∈ T (Ω) ||q||Q=||q||L2\R
Corresponding weak formulation is given in equation (4.1).
• Element 2. Continuous solenoidal velocity space with homogeneous Dirichlet boundary condi-tions: Vk1 div0:0≡v ∈ [H 1(Ω)]N : v = 0 on ΓD : div v = 0 in Ω : v|K ∈ {Pk1(K)} N ∀K ∈ T (Ω) ||v||V=||v||H1
Corresponding weak formulation is given in equation (4.10).
• Element 3. Discontinuous velocity space and continuous pressure space: Vk1 =v ∈ L2(Ω)N : v |Ke∈ {Pk1(Ke) N } ∀ Ke∈ Ph(Ω) ||v||V =||v||Vbo or||v||Vsipg Qk2 0 = q∈ H1(Ω) : q|K ∈ Pk2(K)∀K ∈ Ph(Ω) : Z Ω q dx = 0 ||q||Q =||q||L2 \R
Corresponding weak formulation is given in equation (5.4).
• Element 4. Discontinuous solenoidal velocity space and continuous pressure space on skeleton: Vk1 div0 =v ∈ L 2(Ω)N : v |K ∈ {Pk1(K) N : div v = 0 } ∀ K ∈ Ph(Ω) ||v||V=||v||Vbo or ||v||Vsipg Qk2 0 = q∈ L2(ξh)∪ C0(ξh) : q|e∈ Pk2(e)∀ e ∈ ξh: Z ξ q ds = 0 ||q||Q= min c∈Rkph− ckξ
Corresponding weak formulation is given in equation (5.21).
• Element 5. Discontinuous solenoidal velocity space and discontinuous pressure space on skeleton: Vk1 div0=v ∈ L 2(Ω)N : v |Ke∈ {Pk1(K) N : div v = 0 } ∀ K ∈ Ph(Ω) ||v||V=||v||Vbo or||v||Vsipg Qk2 0 = q∈ L2(ξh) : q |e∈ Pk2(e)∀ e ∈ ξh: Z ξ q ds = 0 ||q||Q= min c∈Rkph− ckξ
• Element 6. Discontinuous solenoidal velocity space. Vk1 div0=v ∈ L 2(Ω)N : v |Ke∈ {Pk1(K) N : div v = 0 } ∀ K ∈ Ph(Ω) ||v||V=||v||Vbo or||v||Vsipg
Corresponding weak formulation is given in equation (5.38).
000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 111111111111 111111111111 111111111111 111111111111 111111111111 111111111111 111111111111 111111111111 111111111111 111111111111 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 111111111111 111111111111 111111111111 111111111111 111111111111 111111111111 111111111111 111111111111 111111111111 111111111111 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 111111111111 111111111111 111111111111 111111111111 111111111111 111111111111 111111111111 111111111111 111111111111 111111111111 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 000000000000 111111111111 111111111111 111111111111 111111111111 111111111111 111111111111 111111111111 111111111111 111111111111 111111111111 000 000 000 111 111 111 000 000 111 111 000 000 111 111 000 000 111 111 000 000 111 111 000 000 111 111 000 000 111 111 000 000 111 111 000 000 111 111 000 000 111 111 000 000 111 111 0000 0000 0000 1111 1111 1111 0000 0000 0000 1111 1111 1111 0000 0000 0000 1111 1111 1111 0000 0000 0000 1111 1111 1111
Element 1 Element 2 Element 3
Element 4 Element 5 Element 6
Continuous velocity Continuous pressure
Continuous pressure on skeleton Discontinuous velocity
Discontinuous pressure on skeleton Divergence free basis function for velocity
Figure 3: 2D triangular elements
To determine the stability of the elements we determine the inf-sup constant, which is constant α from the first equation of Theorem 3.3. For mixed problems we determine the inf-sup constant α of the aggregated problem using equation (3.12). As said in Section 3.3, although the inf-sup condition exists of two equations we only have to satisfy the first of them, since the second follows from the first as we take a standard Galerkin approximation (solution space is test space). Table 2 summarizes the information about the various elements above. Here, div0 indicates solenoidal test and solution spaces.
6.3 Numerical determination of the inf-sup constant
Element Weak-form FEM type Pressure term div0 1 (4.1) CG Yes No 2 (4.10) CG No Yes 3 (5.4) DG Yes No 4 (5.21) DG Yes Yes 5 (5.21) DG Yes Yes 6 (5.38) DG No Yes
Table 2: List of elements and corresponding weak formulations of the Stokes equations
6.3.1 Numerical calculation method
To numerically determine the value of the inf-sup constant we use the following approach [14]. To calculate the following inf-sup constant:
inf
u∈V\{0}v sup
∈V\{0}
a(u, v) ||u||V||v||V ≥ α,
we first take the norm for which a(u, v) is continuous, here || · ||V. We then construct the matrix corresponding to the inner product, Wij = (ψi, ψj)V, where{ψi} forms a basis of the finite-element space. Being a positive definite matrix, W admits a Cholesky factorization according to W = UT
V·UV.
The norm of a function in V is then related to the vector norm of its coefficients ¯vas follows: ||v||2V= ¯v
TWVv¯ = ¯vTUT
VUVv¯ =||UVv¯||2 ¯v∈ RN,
and the following holds max v∈V |a(u, v)| ||v||V = max ¯ v∈RN |(¯v, A¯u)| ||UV¯v|| = maxw¯∈RN |( ¯w, U−TV A¯u)| || ¯w|| =||U −T V A¯u||,
where in the last step we find the maximum by taking ¯w= U−TV A¯u. Now using this result the inf-sup constant becomes: αh= min u∈Vmaxv∈V |a(u, v)| ||u||V||v||V = min ¯ u∈RN ||U−TV A¯u|| ||UVu¯|| , = min ˜ u∈RN ||U−TV AU−1V u˜|| ||˜u|| , = min ˜ u∈RN, ||˜u||=1|| ˜A˜ u||, = min ˜ u∈RN ,||˜u||=1 q ˜ uTA˜TA˜˜u,
where ˜u= UVu¯ and ˜A= U−TV AU−1V . Using the Courant-Fischer theorem the last equation finally becomes αh= q λmin( ˜ATA),˜ = q σ2 min( ˜A), = σmin(U−TV AU−1V ), (6.1)
6.3.2 Numerically determined inf-sup constants
With the method described in the previous section, the inf-sup constants of the finite elements are derived. The mesh size h−1is varied between one and six as is the polynomial interpolation degree of the velocity field. For elements with pressure fields the polynomial interpolation degree is taken equal to that of the velocity field. The values of the determined inf-sup constants can be seen in Figures 4 to 13. −0.8 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 −19.5 −19 −18.5 −18 −17.5 −17 −16.5 −16 −15.5 k= 1 k= 2 k= 3 k= 4 k= 5 k= 6 Element 1 log10(h) lo g10 (α )
Figure 4: Inf-sup coefficient for element 1: As can be seen, α ≈ 0. Hence, the element is unstable for k2= k1. This was expected, since the Taylor-Hood type elements are only stable for k2≤ k1− 1. −0.8 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 −9 −8 −7 −6 −5 −4 −3 −2 −1x 10 −3 k= 1 k= 2 k= 3 k= 4 k= 5 k= 6 log10(h) lo g10 (α ) Element 2
Figure 5: Inf-sup coefficient for element 2: As said before, for some of the combinations of k and h there is no continuous solenoidal ve-locity field for the chosen field. The values of α for these combinations are not derived, for instance for k = 1. For other values of k, the inf-sup constant α seems to converge to a constant value slightly smaller than one. For those values the element is stable.
For elements one, three, four and five also tests have been performed with k2 = k1 − 1 and k2= k1+ 1. Results of those tests are summed up in Table 3, where k1 varies between one and six for k2= k1− 1 and between two and six for k2= k1+ 1.
Element k2= k1− 1 k2= k1+ 1
1 bounded for all k1and h≤ 0.5 not bounded 3 BO-DG bounded for all k1 bounded for k1≥ 2 3 SIPG-DG bounded for all k1 bounded for all k1
4 BO-DG bounded for all k1 bounded for all even order k1 4 SIPG-DG bounded for all k1 bounded for all even order k1 5 BO-DG bounded for all k1 not bounded
5 SIPG-DG bounded for all k1 not bounded
Table 3: Examination of the boundedness of inf-sup coefficient α for k2= k1− 1 and k2= k1+ 1
−0.8 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 −20 −18 −16 −14 −12 −10 −8 −6 −4 −2 0 k= 1 k= 2 k= 3 k= 4 k= 5 k= 6 Element 3 log10(h) lo g10 (α )
Figure 6: Inf-sup coefficient for element 3 with BO-DG: α ≈ 0 for k1 = 1, but for all other k1, α > 0. Hence, this element is stable with BO-DG and k2= k1 as long as k1≥ 2.
−0.8 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 −3 −2.8 −2.6 −2.4 −2.2 −2 −1.8 −1.6 k= 1 k= 2 k= 3 k= 4 k= 5 k= 6 Element 3 log10(h) lo g10 (α )
Figure 7: Inf-sup coefficient for element 3 with SIPG-DG: α > 0 for all combinations of k2= k1 and α is bounded from below for all k1. Hence, the element is stable.
element when the number of internal edges is lower than the space dimension N . Taking the pressure space one degree higher than the velocity space is possible for some elements, this can be especially useful for grid adaptive strategies.
−0.8 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 −18 −16 −14 −12 −10 −8 −6 −4 −2 0 k= 1 k= 2 k= 3 k= 4 k= 5 k= 6 Element 4 log10(h) lo g10 (α )
Figure 8: Inf-sup coefficient for element 4 with BO-DG: α ≈ 0 for k1 = 1, but for all other k1, α > 0. Hence, this element is stable with BO-DG and k2= k1 as long as k1≥ 2.
−0.8 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 −3.4 −3.2 −3 −2.8 −2.6 −2.4 −2.2 −2 −1.8 −1.6 k= 1 k= 2 k= 3 k= 4 k= 5 k= 6 Element 4 log10(h) lo g1 0 (α )
Figure 9: Inf-sup coefficient for element 4 with SIPG-DG: α > 0 for all combinations of k2= k1 and α is bounded from below for all k1. Hence, the element is stable.
−0.8 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 −2.6 −2.4 −2.2 −2 −1.8 −1.6 −1.4 −1.2 k= 1 k= 2 k= 3 k= 4 k= 5 k= 6 Element 5 log10(h) lo g10 (α )
Figure 10: Inf-sup coefficient for element 5 with BO-DG: α > 0 for all k ≥ 1 and seems to converge to a constant value with respect to h for k1 ≥ 2. Interestingly also for k1 = 1 α > 0, however α is not bounded from below.
−0.8 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 −3.6 −3.4 −3.2 −3 −2.8 −2.6 −2.4 −2.2 −2 −1.8 k= 1 k= 2 k= 3 k= 4 k= 5 k= 6 Element 5 log10(h) lo g1 0 (α )
−0.8 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 −1.6 −1.5 −1.4 −1.3 −1.2 −1.1 −1 −0.9 −0.8 −0.7 k= 1 k= 2 k= 3 k= 4 k= 5 k= 6 Element 6 log10(h) lo g10 (α )
Figure 12: Inf-sup coefficient for element 6 with BO-DG: α is again not bounded from below for k1 = 1, but for all other k1, α > 0. Hence, this element is stable for k1≥ 2.
−0.8 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 −0.2 −0.18 −0.16 −0.14 −0.12 −0.1 −0.08 −0.06 −0.04 −0.02 k= 1 k= 2 k= 3 k= 4 k= 5 k= 6 Element 6 log10(h) lo g10 (α )
