# A locally mass conserving quadratic velocity, linear pressure element

By supplementing the pressure space for the Taylor-Hood element a triangular element that satisfies continuity over each element is produced. Making a novel extension of the patch argument to prove stability, this element is shown to be globally stable and give optimal rates of convergence on a wide range of triangular grids. This theoretical result is extended in the discussion given in the appendix, showing how optimal convergence rates can be obtained on all grids. Two examples are presented, one illustrating the convergence rates and the other illustrating difficulties with the Taylor-Hood element which are overcome by the element presented here.

There are no comments yet.

## Authors

• 1 publication
• 1 publication
• ### Convergence rates for the numerical approximation of the 2D stochastic Navier-Stokes equations

We study stochastic Navier-Stokes equations in two dimensions with respe...
06/27/2019 ∙ by Dominic Breit, et al. ∙ 0

• ### Velocity-vorticity-pressure formulation for the Oseen problem with variable viscosity

We propose and analyse an augmented mixed finite element method for the ...
02/10/2021 ∙ by Veronica Anaya, et al. ∙ 0

• ### Convergence of a finite element method for degenerate two-phase flow in porous media

A finite element method with mass-lumping and flux upwinding, is formula...
01/24/2020 ∙ by Girault Vivette, et al. ∙ 0

• ### Optimal quadratic element on rectangular grids for H^1 problems

In this paper, a piecewise quadratic finite element method on rectangula...
03/03/2019 ∙ by Huilan Zeng, et al. ∙ 0

• ### Flux-mortar mixed finite element methods on non-matching grids

We investigate a mortar technique for mixed finite element approximation...
08/21/2020 ∙ by Wietse M. Boon, et al. ∙ 0

• ### Convergence analysis of oversampled collocation boundary element methods in 2D

Collocation boundary element methods for integral equations are easier t...
03/31/2021 ∙ by Georg Maierhofer, et al. ∙ 0

• ### Reconstructed Discontinuous Approximation to Stokes Equation in A Sequential Least Squares Formulation

We propose a new least squares finite element method to solve the Stokes...
03/04/2020 ∙ by Ruo Li, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

A popular triangular element for solving two-dimensional flows was introduced by Hood & Taylor [9]. It has the serious physical drawback that continuity is only satisfied over the whole domain of the problem and not over each element. A consequence of this phenomenon and the poor approximation that can result is given by Tidd, Thatcher & Kaye [11] and a further example of poor results is illustrated in section 6. Tidd et al. show that by supplementing the continuous linear pressure by functions constant in each element, not only is continuity satisfied locally but also the quality of the solution is greatly improved, at least for the particular problem that they were considering. This idea of supplementing the pressure space had previously been suggested by Gresho et al. [5] and is also discussed by Griffiths [7]. In this paper a stability analysis is presented which shows that this element is stable on a wide range of triangular grids.111Added in 2020: The stability result established here has since been generalised by Boffi et al. (Journal of Scientific Computing, 52:383–400, 2012) to cover Hood–Taylor triangular and tetrahedral meshes that are augmented by adding piecewise constant functions to the pressure space of continuous piecewise polynomials of degree ( in 2D and in 3D) under the restriction that every element has at least one vertex in the interior of the domain. Details of the formulation of the continuous and discrete equations for Stokes flow and the implications of the analysis of Stokes flow for Navier–Stokes flow will not be given here. Only sufficient detail to introduce the notation is presented; for further details see Girault & Raviart [4].

The Stokes equations, in weak form, may be written

 find (u,q)∈(H10(Ω))2×L20(Ω) such that ν(∇u,∇v)Ω−(divv,q)Ω =(f,v)Ω∀v∈(H10(Ω))2, (1) (divu,p)Ω =0∀p∈L20(Ω), (2)

where

1. ,

2. ,

3. .

The discrete analogue of (1) and (2) in the finite element subspaces and is given by

 find (uh,q)h∈Vh×Ph such that ν(∇uh,∇vh)Ω−(divvh,qh)Ω =(f,vh)Ω∀vh∈Vh, (3) (div\,uh,ph)Ω =0∀ph∈Ph. (4)

The essential result for stability and convergence is the discrete ‘inf–sup’ condition

 infp∈Phsupv∈Vh{(divv,p)Ω|v|1,Ω∥p∥L2(Ω)}≥β>0 (5)

with independent of . In general, such a condition can only be satisfied if all triangular elements satisfy a regularity condition of the form222Added in 2020: This condition may be relaxed. Certain approximation methods (including Taylor–Hood) are known to be inf–sup stable on highly stretched grids.

 h△≤σρ△ (6)

where and the parameters and are respectively the diameter of element and the diameter of the largest circle inside .

For the Taylor–Hood element on a triangulation of , and are given by

 Ph ={p∈L20(Ω)∩C(¯¯¯¯Ω)∣p∈P1(△) ∀△∈Th} (7) Vh ={v∈(H10(Ω))2∣vi∈P2(△)%fori=1,2 and ∀△∈Th}, (8)

where is the set of all polynomials of degree less than or equal to in . The fact that the Taylor–Hood element satisfies the ‘inf–sup’ condition (5) was originally shown by Bercovier & Pironneau [1], although the results can now be proved quite easily using the patch ideas of Boland & Nicolaides [2] or equivalently Stenberg [10].

For the element discussed here, the space is the same and given by (8) but is defined by

 Ph={p∣p=p0+p1, p0∈L20(Ω), p0∈P0(△), p1∈L20(Ω)∩C(¯¯¯¯Ω), p1∈P1(△),∀△∈Th}. (9)

## 2 Stability of patches of elements

Proving the discrete ‘inf–sup’ condition became relatively easy only after the ideas of locally stable patches of elements were developed. The continuous linear plus constant pressure element does not fit neatly into the local analysis of Stenberg [10] nor Boland & Nicolaides [2] because on any patch there are always two distinct ways of producing a constant function in the pressure space (i.e., constant at the vertex nodes and zero at the centroids or zero at the vertices and constant at the centroids). The analysis presented here is in the spirit of the approach by Stenberg [10].333Added in 2020: In retrospect, a simpler and more elegant way of establishing stability would be to employ the construction used in the analysis of the Taylor–Hood element in Girault & Raviart [4, pp.176–180] together with the overlapping patch framework developed by Stenberg in his follow-up paper (Mathematics of Computation, 54:495–508, 1990).

The notation is used for the class of patches of elements topologically equivalent to the patch of elements . By a patch of elements it is understood that it is a union of elements, each of which has at least one side in common with another element of the patch. For the precise definition, see Stenberg [10]. By way of illustration we note that the patches 1.2, 1.3, 1.4 in figure 1 are topologically equivalent but they are not equivalent to the patch 1.1. Indeed all patches of 3 elements are topologically equivalent to either 1.1 or 1.2. The only constraint on is that all elements must satisfy the regularity constraint (6) for some value of .

To prove patch stability for a class of patches , we first consider a typical patch . For this patch we define

 (i) VhM ={v∈(H10(M))2∣vi∈P2(△) for% i=1,2 and ∀△∈M}, (10) (ii) PhM ={p∣p=p0+p1, p0∈P0(△) ∀△∈M, p1∈C(¯¯¯¯¯¯M), p1∈P1(△) ∀△∈M}. (11)

The first step in proving patch stability is to find a subspace of that satisfies the condition

 infp∈PhMsupv∈VhM{(divv,p)M|v|1,M∥p∥L2(M)}≥βM>0. (12)

The condition (12) is equivalent to the condition

 p∈RhM,(divv,p)M=0∀v∈VhM⟹p=0. (13)

This condition can be established by looking at the null space of the matrix defined by

 v∈VhM,p∈PhM,(divv,p)M=v–TBp–, (14)

where

is the vector of nodal coordinates determining the finite element function

and is the vector of nodal coordinates determining . We will then show that the constraints to produce from annihilate this null space. Such a process will establish (12) for the particular patch and with the particular definition of .

Both of the three-element patches 1.1 and 1.2 are important in subsequent sections of this paper and we will define an and show that the constraints annihilate the null space in both cases.

### 2.1 The patch 1.1

A typical patch of this type, type 1, is illustrated444Added in 2020: The hand-drawn figures are reproduced here exactly as in the original report. in figure 2. For this patch we define

 RhM={p∣p∈PhM, p0∈L20(M), p1∈L20(M)}. (15)

The patch has three elements and we have

 v–T ={(V1)1,(V1)2,(V2)1,(V2)2,(V3)1,(V3)2,(V4)1,(V4)2,}, p–T ={P1,P2,…,P7},

with an matrix. We denote by and the area of and respectively and by the two-dimensional vector in the usual values

 (b(i)1)1=y(i)2−y(i)3,(b(i)1)2=x(i)3−x(i)2, etc.

with the local nodes of each element illustrated in figure 2. The matrix is given by

 16⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣−b(1)2b(2)1+b(1)1−b(2)3b(2)2+b(1)3−4b(1)2−4b(2)30−b(3)3−b(2)2b(3)1+b(2)1b(3)2+b(2)30−4b(2)2−4b(3)3b(1)1+b(3)1−b(1)3−b(3)2b(1)2+b(3)3−4b(1)30−4b(3)20000−b(1)1−b(2)1−b(3)1⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦ (16)

and has null space spanned by

 ST1=(1,1,1,1,0,0,0),ST2=(0,0,0,0,1,1,1). (17)

We denote by the matrix, the first column of which is and the second is . The constraints that give from are

1. ,

2. ,

which constrain the vector by

1. ,

2. ,

which we write as

 Hp–=0–, (18)

where is a matrix. The only vector in the null space of that satisfies both constraints is the zero vector if the matrix

 CM=HS (19)

is of full column rank. Here is given by

 [0|M|3|M|0]. (20)

Clearly (20) is nonsingular, therefore for any patch of type 1 with given by (15), the inequality (12) is satisfied provided is nonempty. In section 3 we will construct functions belonging to .

### 2.2 The patch 1.2

A typical patch of this type, type 2, is illustrated in figure 3 with a four-dimensional vector and eight dimensional.

The matrix is the matrix

 16[−b(1)3b(2)1+b(1)1−b(2)20b(2)3+b(1)2−4b(1)3−4b(2)200−b(2)3b(3)1+b(2)1−b(3)2b(3)3+b(2)20−4b(2)3−4b(3)2] (21)

which has a null space spanned by

 ST1=(1,1,1,1,1,0,0,0),ST2=(0,0,0,0,0,1,1,1),ST3=(4,0,0,0,0,−1,0,0),ST4=(0,0,0,4,0,0,0,−1). (22)

Here we note that the matrix is an matrix, the th column of which is . It is not possible for the two constraints to produce given by (15) to give a matrix of full column rank (because with these constraints is a matrix). Thus we need to apply further constraints on to annihilate the null space. We define by

 RhM={p∣p∈PhM, p0∈L20(M), p1∈L20(M),p∈L20(△(i))∀△(i) with two sides on the boundary of M}. (23)

Thus the constraints are

1. ,

2. ,

3. ,

4. ,

giving us the constraint equation of the form (18). Thus, here,

 CM=HS=⎡⎢ ⎢ ⎢ ⎢⎣0|M|−|△(1)|−|△(3)|3|M|04|△(1)|4|△(3)|33103301⎤⎥ ⎥ ⎥ ⎥⎦. (24)

Clearly, (24) is nonsingular, therefore for any patch of the type 2 with given by (23), the inequality (12) is satisfied provided is nonempty. In section 4 we will construct functions belonging to for this patch.

### 2.3 Stability over classes of topologically equivalent patches

In this section we shall use the important result in Stenberg [10], who observed that, given a class of topologically equivalent patches which satisfy (12) for every , then the value of

 βM=infp∈RhMsupv∈VhM{(divv,p)M|v|1,M∥p∥L2(M)} (25)

is independent of a transformation of the form

 ~x––=A(x––−α––), (26)

where and are fixed. This, together with the regularity constraint (6) allows us to represent the whole range of values of for as a function over a compact set , every point of which represents a patch (or indeed many patches) belonging to . Thus there exists such that

 β=minR(βM)=minM∈εσM(βM) (27)

for which

 infp∈RhMsupv∈VhM{(divv,p)M|v|1,M∥p∥L2(M)}≥β>0 (28)

for every with independent of and and depending only on the topology of and . Thus here, for a given , we have two classes and of elements topologically equivalent to the patches of type 1, and 2 respectively, and the inequality (28) holds for these two classes with different values of . We will actually use an equivalent form of (28), namely for each , there exists such that

 (divvM,p)M≥β∥p∥2L2(M),|vM|1,M≤C∥p∥L2(M), (29)

or, for ,

 (div(γvM),p)M≥βγ∥p∥2L2(M),|(γvM)|1,M≤Cγ∥p∥L2(M), (30)

with and depending only on the class of patches and on .

## 3 Global stability of grids made up of patches of the type 1

In order to use the results of the previous section we construct an operator by

 Π(1)Mp =1|M|∫MpdΩ=k0+k1, (31) wherek0=1|M|∫Mp0dΩ,k1=1|M|∫Mp1dΩ. (32)

Thus, since

 ∫M(p0−k0)dΩ=0,∫M(p1−k1)dΩ=0

then using the argument from Stenberg [10] we can establish the following theorem and corollary.

###### Theorem 3.1.

For every ,

1. , defined by (15),

2. .

###### Corollary (Corollary to Theorem 3.1).

For each , there exists such that

 (divvM,p)M≥β1∥p−Π(1)Mp∥2L2(M),|vM|1,M≤C1∥p−Π(1)Mp∥L2(M),

where and are independent of and but are dependent on .

The proof of global stability on a grid made up of patches of the type 1 now follows that given by Stenberg [10] or Boland & Nicolaides [2]. Moreover, we have the optimal rates of convergence to the solution of the Stokes equation555Added in 2020: Assuming additional smoothness (that is, regularity) of the target solution., namely

 |u−uh|1,Ω+∥q−qh∥L2(Ω)≤Ch2(|u|3,Ω+∥q∥2,Ω) (33)

where is the numerical solution in .

## 4 Further results on the patch of the type 2

Before we establish global stability and optimal rates of convergence, there are further results required for a patch of this type. Let be a patch of the type 2, we construct an operator on the set by

 Π(2)Mp={k(2)in △(2),4k(2)(1−3L(i))−3k(i)(1−4L(i))% in △(1) and △(3), (34)

where

1. is defined in figure 3,

2. is the areal coordinate in , and , equal to 1 at the intersection of the two sides on the boundary of ,666Added in 2020: Thus with the numbering shown in figure 3, we have in and in .

3. (35)

For every then

1. (36)

2. (37)

###### Proof.

(A) Clearly , thus we have to show that satisfies the constraints from to . We write

 [Π(2)Mp]0 ={(1−α)k(2)in △(2)(1−α)k(2)+3(k(2)−k(i))in △(i),i=1 and 3, (38) [Π(2)Mp]1 ={αk(2)in △(2)αk(2)−12(k(2)−k(i))L(i)in △(i),i=1 and 3, (39)

with

 α=∫Mp1dΩ+4(k(2)−k(1))|△(1)|+4(k(2)−k(3))|△(3)|k(2)|M| (40)

if . (If then the definition of for is independent of .)

We now observe that

 ∫Mp0−[Π(2)Mp]0dΩ=0, ∫Mp1−[Π(2)Mp]1dΩ=0, ∫△(1)p−Π(2)MpdΩ=0, ∫△(3)p−Π(2)MpdΩ=0.

Thus we have established part (A).

(B) Taking the left-hand side of (37),

 (divv,Π(2)Mp)M=(divv,k(2))M+3(k(2)−k(1))(divv,(1−4L(1)))△(1)+3(k(2)−k(3))(divv,(1−4L(3)))△(3). (41)

The first term in (41) is zero because . If is in then , thus we see immediately777Added in 2020: This can be verified by direction computation; whenever . that the second term in (41) is zero and so therefore is the third. ∎

###### Corollary (Corollary to Theorem 4.1).

For each there exists such that

 (divvM,p)M ≥β2∥p−μM∥2L2(△(2)), (42) |vM|1,M ≤C2∥p−μM∥L2(△(2)), (43)

with

 μM=1|△(2)|∫△(2)pdΩ, (44)

where and depend only on the regularity constant .

###### Proof.

By Theorem 4.1 and inequalities (30), for each there exists such that, for ,

 (div(γ~vM),p)M=(div(γ~vM),p−Π(2)Mp)≥β2γ∥∥p−Π(2)Mp∥∥2L2(M), |γ~vM|1,M≤c2γ∥∥p−Π(2)Mp∥∥L2(M).

The result follows by choosing

 (45)

noting that , and choosing . ∎

## 5 Global stability of grids made up of patches of the type 2

We assume that is a polygonal region which has been triangulated into triangular elements . We further assume that each element sits inside an extended patch of the type 2 as element . Firstly we note that these patches overlap. Secondly we note that this assumption does exclude some grids of triangles, namely those which either

1. contain triangular elements with two sides on the boundary, or

2. contain a patch of the form 1.1 with that patch having a side on the boundary.

We can see this by looking at figure 4, elements to all sit as element of a patch of the type 2 but elements and do not.

We consider such that is constant in each , Girault & Raviart [3] show that for each there exists such that

 (div~v,μ)Ω=∥μ∥2L2(Ω),|~v|1,Ω≤~c∥μ∥L2(Ω). (46)

Let , we define such that in the element is the constant value

 μ=1|E(i)|∫E(i)pdΩin E(i) % for i=1,…,N; (47)

thus by the corollary to Theorem 4.1, for each there exists such that

 (divvMi,p)Mi ≥β2∥p−μ∥2L2(E(i)), (48) |vMi|1,Mi ≤C2∥p−μ∥L2(E(i)). (49)

We note that and do not depend on the particular patch and depend only on the regularity constant . Thus, for each there exists such that

 (div~v,p)Ω ≥~β∥p−μ∥2L2(Ω), (50) |~v|1,Ω ≤~c∥p−μ∥L2(Ω), (51)

where

 (a) ~v =N∑i=1~vi, (b) ~vi ={vMiin Mi,0elsewhere in Ω, (52) (c) ~C =C2,~β=β2.

Following Stenberg [10] the inequalities (46) and (50) establish that for each then

 v=~v+(2~β1+~C2)~v (53)

satisfies

 (divv,p)Ω ≥(~β1+~C2)∥p∥2L2(Ω), (54) |v|1,Ω ≤(~C+2~β~C1+~C2)∥p∥L2(Ω). (55)

Thus we have global stability and optimal rates of convergence on all grids that satisfy the regularity constraint (6) and restrictions on the triangles mentioned at the beginning of this section.

It is only the former of these restrictions, namely that we must triangulate ‘into the corners’ that is an essential restriction. By including patches of the form 1.1 in the above argument then the second restriction can be removed but now and . Thus we have established optimal convergence rates on all grids of triangles provided the grid has been triangulated into the corners. Further discussion of this topic when the grid has not been triangulated into the corners is given in the appendix.

## 6 Numerical examples

In this section we shall consider two numerical examples. The first is a simple test problem to demonstrate that optimal convergence rates are achieved and the second is an example where the Taylor–Hood element gives poor results but the linear plus constant pressure element, which we shall call the LC element888Added in 2020: The mixed approximation method is referred to in the book by Elman et al. [Finite Elements and Fast Iterative Solvers, Oxford University Press, 2014), and as the enhanced Hood–Taylor scheme in the book by Boffi et al. [Mixed Finite Element Methods and Applications, Springer, 2013]., gives relatively good results.

### 6.1 Testing rates of convergence

The first test problem is one proposed by Griffiths & Mitchell [6]. It is an enclosed flow problem (namely a Stokes flow) in the unit square with solution

 vx=−20xy3,vy=5y4−5x4,p=−60x2y+20y3+5. (56)

Typical grids for this test problem are illustrated in figure 5 and the solutions (i.e. norms

of the errors) are presented in table 1 with the results for the LC element compared with the Taylor–Hood element and the Raviart bubble element; see Girault & Raviart [4]. It can be seen that the error for the LC and Taylor–Hood elements are comparable and both are much smaller than the Raviart bubble element for this problem. Moreover, for all three elements, the optimal convergence rates are observed.

### 6.2 Illustrating difficulties with the Taylor–Hood element

The solution of a non-Newtonian fluid in the volume of revolution of the region illustrated in figure 6, with the boundary conditions given in the figure, can be reduced to solving the following set of equations

 1R(∂2Vr∂r2+∂2Vr∂z2+1r∂Vr∂r−Vrr2)−∂p∂r =−NnV2θr+Ne[(∂Vθ∂r−Vθr)2+(∂Vθ∂z)2], (57a) 1R(∂2Vθ∂r2+∂2Vθ∂z2+1r∂Vθ∂r−Vθr2)