# Application of a minimal compatible element to incompressible and nearly incompressible continuum mechanics

In this note we will explore some applications of the recently constructed piecewise affine, H^1-conforming element that fits in a discrete de Rham complex [Christiansen and Hu, Generalized finite element systems for smooth differential forms and Stokes' problem. Numer. Math. 140 (2018)]. In particular we show how the element leads to locking free methods for incompressible elasticity and viscosity robust methods for the Brinkman model.

## Authors

• 27 publications
• 1 publication
• 11 publications
• ### Exact sequences on Worsey-Farin Splits

We construct several smooth finite element spaces defined on three–dimen...
08/12/2020 ∙ by Johnny Guzmàn, et al. ∙ 0

• ### Local Finite Element Approximation of Sobolev Differential Forms

We address fundamental aspects in the approximation theory of vector-val...
11/01/2020 ∙ by Evan S. Gawlik, et al. ∙ 0

• ### A lowest-degree strictly conservative finite element scheme for incompressible Stokes problem on general triangulations

In this paper, we propose a finite element pair for incompressible Stoke...
08/24/2021 ∙ by Wenjia Liu, et al. ∙ 0

• ### Helicity-conservative finite element discretization for MHD systems

We construct finite element methods for the magnetohydrodynamics (MHD) s...
07/15/2020 ∙ by Kaibo Hu, et al. ∙ 0

• ### Preconditioning immersed isogeometric finite element methods with application to flow problems

Immersed finite element methods generally suffer from conditioning probl...
08/11/2017 ∙ by Frits de Prenter, et al. ∙ 0

• ### Fortin Operator for the Taylor-Hood Element

We design a Fortin operator for the lowest-order Taylor-Hood element in ...
04/28/2021 ∙ by Lars Diening, et al. ∙ 0

• ### An iterative generalized Golub-Kahan algorithm for problems in structural mechanics

This paper studies the Craig variant of the Golub-Kahan bidiagonalizatio...
08/23/2018 ∙ by Mario Arioli, 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

It is well known that standard finite element methods are not in general well-suited for the approximation of nearly incompressible elasticity or incompressible flow problems. Indeed, in particular low order approximation spaces often suffer from locking in the incompressible limit [2]. They may typically also exhibit instability when Darcy flow is considered if the element was designed for Stokes’ problem [16]. These problems can be alleviated using stabilization [4, 3], but such stabilizing terms, although weakly consistent to the right order, may upset local conservation of e.g. mass, momentum, and introduce an additional layer of complexity to the computational method and its analysis. Recently some new results on -conforming piecewise polynomial approximation spaces compatible with the de Rham complex have been published [11, 12, 15, 13, 8, 6]. Such elements are interesting, since they provide a simple tool for the robust approximation of models in mechanics where a divergence constraint is present. Herein we will focus on the piecewise affine element derived in the last reference. The advantage of this approach is that it offers a simple low order locking free element in arbitrary space dimensions. Observe that for the Scott-Vogelius element the polynomial order of the spaces typically depends on the number of dimensions [8]

. We discuss how this element can be implemented in engineering practice and show the basic, robust, error estimates that may be obtained for linear elasticity and incompressible flow. In this paper we will consider two different models, linear elasticity and the Brinkman model for porous media flow. The idea is to show the locking free property of the element on the elasticity model and then illustrate how the element seamlessly can change between the Stokes’ equations modelling free flow and Darcy’s equations modelling porous media flow, while remaining

-conforming. The two models are introduced in section 2. The construction of the element is discussed in section 3 and the finite element discretizations of the model problems and their analysis are the topics of section 4 and 5. In section 6 we discuss how boundary conditions may be imposed weakly using Nitsche’s method, without sacrifying the good properties of the element. Finally section 7 gives some numerical illustrations to the theory.

## 2 Model problems: linear elasticity and the Brinkman model

We will consider two model problems with solutions in , initially assuming homogeneous Dirichlet boundary conditions. Let , denote a convex polyhedral domain with boundary . The first model problem is linear elasticity. Here we wish to find , where , such that

 −∇⋅σ(u)=f, in Ω, (2.1)

where , with

the symmetric part of the gradient tensor,

the identity matrix and

the Lamé coefficients and . This system can be written on weak form: find such that

 aE(u,v)=l(v), for all v∈V0,

where

 aE(w,v):=∫Ωσ(u):∇sv dx, (2.2)

where the tensor product is defined by and

 l(v):=∫Ωf⋅v dx. (2.3)

It is well-known that the problem (2.1) admits a unique weak solution in the space through application of Lax-Milgram’s lemma, and that the following regularity holds [10],

 ∥u∥H2(Ω)+λ∥∇⋅u∥H1(Ω)≤CR∥f∥Ω for μ∈[μ1,μ2] and λ∈(0,∞). (2.4)

The second model problem is the Brinkman problem where we look for a velocity-pressure couple , where denotes the set of square integrable functions with mean zero, such that

 −μΔu+σu+∇p=f in Ω∇⋅u=g in Ω. (2.5)

Here , , is the viscosity coefficient and a possibly space dependent coefficient modelling friction due to the porous medium. Observe that if we recover the Darcy model for porous media flow and if we obtain the classical Stokes’ system for creeping incompressible flow.

The corresponding weak formulation reads: find such that:

 AB[(u,p),(v,q)]=l(v), for all % (v,q)∈V0×Q.

Here the bilinear forms are given by

 AB[(u,p),(v,q)]:=aB(u,v)−b(p,v)+b(q,u) (2.6)

with

 aB(w,v):=∫Ωμ∇w:∇v+σw⋅v dx,
 b(q,v):=∫Ωq∇⋅v dx

and

 lB(v,q):=∫Ωf⋅v dx+∫Ωgq dx. (2.7)

By the surjectivity of the divergence operator we may write where . Unique existence of the part of the solution is ensured through the application of the Lax-Milgram lemma in the space , where

 Hdiv0:={v∈V:∇⋅v=0}.

A unique pressure is then guaranteed by the Ladyzhenskaya-Babuska-Brezzi condition [2].

## 3 The finite element space

Let denote a conforming, shape regular tesselation of into simplices . We denote the set of faces of the simplices in by and the subset of faces that lie on the boundary by . We let denote the space of functions in that are constant on each element,

 Xh:={x∈L2(Ω):x|T∈P0(T);∀T∈Th}.

The -projection on , is defined by for all . satisfies the stability for all and the approximation error estimate

 ∥π0v−v∥Ω≤Ch|v|H1(Ω),∀v∈H1(Ω).

We also introduce the -projection of the trace of a function

 ~π0:L2(∂Ω)↦∂Xh

where

 ∂Xh:={x∈L2(∂Ω):x|F∈P0(F);∀F∈Fb}

where is the set of faces in such that . We let denote the space of vectorial piecewise affine functions on ,

 Wh:={v∈[H1(Ω)]d:v|T∈[P1(T)]d;∀T∈Th}

and define . It is well known that the space is not robust for nearly incompressible elasticity and that velocity-pressure space is unstable for incompressible flow problems. To rectify this we will enrich the space with vectorial bubbles on the faces, following the design in [6], that allows us to remain conforming in , resulting in an extended space, that we will denote . The detailed construction of this space is the topic of the next section. We then apply in the finite element method for the system of compressible elasticity and for the Brinkman system. For the space with built in homogeneous Dirichlet boundary conditions we write . Observe that by construction all functions satisfy .

### 3.1 Construction of the finite element space Vh

The finite element space is constructed by decomposing every simplex in subelements. On these subelements face bubbles are constructed, similar to the face bubbles used in the Bernardi-Raugel element [1]

, but in this case they are constructed using piecewise affine elements. Using the subgrid degrees of freedom similar degrees of freedoms as in the Bernardi-Raugel element are designed as well. The upshot here is that the piecewise affine basis functions are designed so that the divergence restricted to each simplex in the original tesselation is constant. The pressure space then consists of one constant pressure degree of freedom per (macro) simplex, allowing for exact imposition of the divergence free condition. Although the numerical examples in this work are restricted to the two-dimensional case we below for completeness also give a detailed description of the construction in three space dimensions.

We first treat the 2D case for which our numerical examples are implemented and then describe how this extends to the three dimensional case. Consider a triangular element twice subdivided. We call the triangle type I, the first subdivision type II, and the second subdivision type III, cf. Fig 1. The first subdivision is created by joining the centroid of triangle I with its corner nodes. The second subdivision splits each triangle II by the line joining the centroid of triangle I with the centroid of is neighbouring type I triangle sharing the edge to be split. On the boundary we have a free choice of how to split the edge; we here choose to split the edge along the line in the direction of the normal to the boundary. On triangles of type I the approximation is piecewise linear with two velocity degrees of freedom in each corner node. On triangles of type III we add a hierarchical “bubble” approximation in the following way. To the node on the exterior edge

of each triangle of type I is assigned a unit vector

along the line of the split into type III triangles, see Figure 2. The unknown in the corresponding edge node is the vector where is a hierarchical scalar unknown. The centroid–to–centroid nature of the split then ensures continuity of the discrete solution. In the centroid node the bubble has two velocity components determined a priori by setting the divergence equal (with ) on the triangles sharing node and the triangles of type II not being split by . The divergence is set by

 d:=∫Eνi⋅nEds.

The hierarchical bubble is then piecewise linear on these type II triangles and the type III triangles sharing node . Thus, each edge on triangle I has its own unique hierarchical bubble and the total approximation is the sum of the linear function on type I and the three (vector-valued) bubbles.

A closed form for the velocities defining the bubble associated with an edge can be computed beforehand. With the location of the corner, center, and edge nodes according to Fig. 2, with the area of triangle , we find

 um=D(xm−xo) (3.1)

where

 D:=xr(ym−yl)+xm(yl−yr)+xl(yr−ym)2A|xi−xm|.

This gives equal divergence on all subtriangles.

### 3.2 The construction of Vh in three space dimensions

The construction in 3D is analogous to the one in 2D: any given tetrahedron is decomposed using the Worsey Farin (WF) split [17], defined as follows. An inpoint is chosen for the tetrahedron, typically (but not necessarily) the center of the inscribed sphere. As inpoint on the (triangular) faces, one chooses (crucially) the point on the line joining the inpoints on the two neighboring tetrahedra. The faces are then split in three subfaces by joining the inpoint to its vertices. The tetrahedron is split in 12 small tetrahedra, three for each face, based on a subface and with summit at the inpoint of the tetrahedron.

The finite element space on the tetrahedron can then be described as the the space of continuous vectorfields on the WF split which are divergence free, to which one adds one vector field with constant divergence on , namely . As shown in [6] this space has dimension 16. It contains the vectorfields on (dimension 12), and four bubbles attached to faces (dimension 4). As degrees of freedom one may use vertex values and integrals of normal components on faces.

A face bubble can be defined explicitely for a face , as follows. We let be the normalized vector parallel to the line joining the inpoints of the two neighboring tetrahedra of . The vectorfield on has value at vertices of , at the inpoint of the face , and at inpoints of the other faces. At the inpoint of we determine the vector by the condition that the divergence of the vector field is the same on all the small tetrahedra of the WF split and satisfies Stokes’ theorem on the three that are based on .

### 3.3 The Fortin interpolant

For every there exists such that in the vertices of type I simplices, where

denotes the Clément interpolant, and for all

 ∫Fπhu⋅nF ds=∫Fu⋅nF ds.

Note that the interpolant satisfies the approximation error estimate

 ∥πhu−u∥Ω≤C1h|u|H1(Ω),h∥∇(πhu−u)∥Ω+∥πhu−u∥Ω≤C2h2|u|H2(Ω). (3.2)

The proof of the existence of is identical to that of the interpolant for the Bernardi-Raugel element [1]. Note that for functions such that there holds that .

It follows from this construction that for all and for all , using the divergence theorem we have

 ∫T∇⋅πhuqh dx=∫∂T(πhu⋅n∂T)qh ds=∫∂T(u⋅n∂T)qh ds=∫T∇⋅uqh dx=∫Tπ0∇⋅uqh dx.

A consequence of the existence of the Fortin interpolant is the existence of a non-trivial subspace such that

 Vdiv(v):={vh∈Vh:∇⋅vh=π0∇⋅v}.

As a consequence, for every there exists

 ζq∈V0h such that ∇⋅ζq=qh and ∥ζq∥H1(Ω)≤C0∥qh∥Ω. (3.3)

To see the note that by the surjectivity of the divergence operator from to for every there exists such that and and if we now consider we see that and we conclude that may be chosen in directly.

## 4 Finite element discretization of the model problems

We consider the finite element spaces that were defined in the previous section. The finite element discretization of the problem (2.1) then takes the form: find such that

 aE(uh,vh)=l(vh), for % all vh∈V0h, (4.1)

where and are defined by (2.2) and (2.3). The finite element method for the problem (2.5) on the other hand takes the form find such that

 AB[(uh,ph),(vh,qh)]=lB(vh,qh), for all (vh,qh)∈V0h×Qh. (4.2)

Both the problem (4.1) and (4.2) admit a unique solution by the same arguments as for the continuous problem. This is also a consequence of the stability estimates that we derive in the next section.

## 5 Stability and error analysis

We introduce two triple norms. First for the elasticity system,

 |||vh|||2E:=2∥μ12∇svh∥2Ω+∥λ12∇⋅vh∥2Ω. (5.1)

Observe that by Korn’s inequality and Poincaré’s inequality the -seminorm is a norm on . Then for the incompressible model we have the triple norm,

 |||vh,yh|||2B:=∥μ12∇vh∥2Ω+∥σ12vh∥2Ω+∥∇⋅vh∥2Ω+∥(μ+σ)−12yh∥2Ω. (5.2)

For the problem (4.1) Korn’s inequality leads to the coercivity, there exists such that for all

 αE|||vh|||2E≤aE(vh,vh). (5.3)

For the problem (4.2) we need to prove an inf-sup condition for stability.

###### Proposition 5.1

(inf-sup stability for the Brinkman problem) There exists such that for all there holds

 αB|||vh,yh|||B≤supwh,qh∈(V0h∖0)×(Qh∖0)AB[(vh,yh),(wh,qh)]|||wh,qh|||B.

Proof. First we take and to obtain

 ∥μ12∇vh∥2Ω+∥σ12vh∥2Ω=AB[(vh,yh),(wh,qh)].

Then we chose , where is defined by (3.3) so that

 (μ+σ)−1∥yh∥2Ω=AB[(vh,yh),(wh,0)]−(μ∇vh,wh)Ω−(σvh,wh)Ω.

Observing now that

 (μ∇vh,∇wh)Ω≤∥μ12∇vh∥Ωμ12(μ+σ)−1C0∥yh∥Ω

and

 (σvh,wh)Ω≤∥σ12vh∥Ωσ12(μ+σ)−1C0∥yh∥Ω

it follows that

 12(μ+σ)−1∥yh∥2Ω≤AB[(vh,yh),(wh,0)]−C20(∥μ12∇vh∥2Ω+∥σ12vh∥2Ω).

Taking and we conclude that

 min(12,12C0)|||vh,yh|||2B≤12∥μ12∇vh∥2Ω+12∥σ−12vh∥2Ω+12C0(μ+σ)−1∥yh∥2Ω≤AB[(vh,yh),(wh,qh)]

To finish the proof note that

 |||wh,qh|||B≤|||vh,yh|||B+|||(2C0)−1(μ+σ)−1ζy,0|||B≤|||vh,yh|||B+(2C0)−1μ12(μ+σ)−1C0∥yh∥Ω+∥∇⋅vh∥Ω≤C|||vh,yh|||B.

Using the stability estimates we may now prove error estimates for the approximations of (4.1) and (4.2).

###### Proposition 5.2

Let be the solution of (2.1) and the solution of (4.1) then

 ∥μ12∇(u−uh)∥Ω+∥λ12(π0∇⋅u−∇⋅uh)∥Ω≤Cinfvh∈Vdiv(u)∥μ12∇(u−vh)∥Ω

and

 ∥μ12∇(u−uh)∥Ω+∥λ12(∇⋅u−∇⋅uh)∥Ω≤Ch(μ12∥u∥H2(Ω)+λ12∥∇⋅u∥H1(Ω))≤CEh∥f∥Ω.

where is independent of .

Proof. Let , with . Note that by adding and subtracting and using the triangle inequality and Korn’s inequality we have

 ∥μ12∇(u−uh)∥Ω+∥λ12(π0∇⋅u−∇⋅uh)∥Ω≤∥μ12∇(u−vh)∥Ω+|||eh|||E.

For the second term we apply the coercivity (5.3), followed by Galerkin orthogonality

 aE(u−uh,wh)=0 for all wh∈V0h

to obtain

 αE|||eh|||2E≤aE(eh,eh)=aE(u−vh,eh).

Noting that

 (λ∇⋅(u−vh),∇⋅eh)Ω=(λ(∇⋅u−π0∇⋅u),∇⋅eh)Ω=0 (5.4)

we may write

 αE|||eh|||2E≤(2μ∇s(u−vh),∇seh)Ω≤2∥μ12∇(u−vh)∥Ω|||eh|||E, (5.5)

which proves the first claim.

The second claim is immediate, taking and using the approximation properties of , (3.2) and the regularity bound (2.4). To show that the constant is independent of observe that .

###### Proposition 5.3

Let be the solution to (2.5), with , and the solution to (4.2). Then there holds

 |||u−uh,π0p−ph|||B≤Cinfvh∈Vgdiv|||u−vh,0|||B

where and

 |||u−uh,π0p−ph|||B≤C(hμ12|u|H2(Ω)+min(C1hσ12|u|H1(Ω),C2h2σ12|u|H2(Ω)).

Proof. We introduce, as before, discrete errors , with and . Using the triangle inequality we see that

 |||u−uh,0|||B≤|||u−vh,0|||B+|||eh,ηh|||B.

For the second term in the right hand side we apply the stability of Proposition 5.1 to obtain

 |||eh,ηh|||B≤supwh,qh∈(V0h∖0)×(Qh∖0)AB[(eh,ηh),(wh,qh)]|||wh,qh|||B.

using Galerkin orthogonality

we have

 AB[(eh,ηh),(wh,qh)]=AB[(u−vh,p−π0p),(wh,qh)]. (5.6)

Observe that by construction we have

 b(qh,u−vh)=0 and b(p−π0p,wh)=0.

The only remaining term in the right hand side of (5.6) is bounded using the Cauchy–Schwarz inequality,

 aB(u−πhu,wh)≤|||u−vh,0|||B|||wh,qh|||B.

This proves the first claim and the second follows as before taking and using the approximation properties of the Fortin interpolant (3.2).
Since we have imposed the boundary conditions strongly above we can not take in the Brinkman model corresponding to the case of the Darcy equations. In order to make this limit feasible we will now discuss weak imposition of boundary conditions using Nitsche’s method.

## 6 Weakly imposed boundary conditions, Nitsche’s method

Here we will discuss how to impose non-penetration conditions on the space as one wishes to do in the case of zero-traction boundary conditions in elasticity and how to relax the no-slip condition when for the Brinkman model. Therefore we here propose Nitsche methods for the imposition of boundary conditions that preserve the locking free character for elasticity and are robust in the limit of pure porous media flow for the Brinkman model.

### 6.1 Zero traction conditions for linear elasticity

Consider first the elasticity problem (2.1), with the boundary decomposed in where and each consists of a set of entire polyhedral faces. We assume that

 tu=gD on ∂ΩD and u⋅n=gN on ∂Ω and t(σ(u)n)=0 on ∂ΩN. (6.1)

Here the tangential projection is defined by . The Nitsche formulation then takes the form: Find such that

 AE,h(uh,vh)=L(vh) (6.2)

with

 AE,h(uh,vh):=aE(uh,vh)−c(uh,vh)−c(vh,uh)+s(uh,vh)

and

 L(vh)=l(vh)+lc(vh)

where

 c(uh,vh):=(n⋅(σ(uh)n),vh⋅n)∂Ω+(t⋅(σ(uh)n),tvh)∂ΩD
 s(uh,vh):=(γ/h(μ+λ~π0) uh⋅n,vh⋅n)∂Ω+(γμ/h tuh,tvh)∂ΩD

and

 lc(vh)=(gN,γ/h(μ+λ~π0) vh⋅n−n⋅(σ(vh)n))∂Ω+(gT,γμ/h tvh−t⋅(σ(vh)n))∂ΩD.

Observe that the projection in the boundary penalty of the normal component is necessary to avoid locking.

We define the stabilization semi-norm by

 |vh|s:=s(vh,vh)12

and the following augmented energy norm defined on

 |||vh|||2E,h:=|||vh|||2E+|vh|2s.

We recall that is a norm by Korn’s inequality and Poincaré’s inequality. We recall the trace inequalities

 ∥v∥∂T≤CT(h−12∥v∥T+h12∥∇v∥T)∀T and v∈H1(T) (6.3)

and

 ∥vh∥∂T≤CTh−12∥vh∥T∀T and vh∈Vh. (6.4)

Using these inequalities it is straightforward to prove the following approximation estimate in the norm and a bound on the form .

###### Lemma 6.1

The following approximation inequality holds

 |||u−πhu|||E,h≤Ch(μ12|u|H2(Ω)+λ12|∇⋅u|H1(Ω)). (6.5)

Proof. The inequality

 |||u−πhu|||E≤Ch(μ12|u|H2(Ω)+λ12|∇⋅u|H1(Ω)).

is immediate by the commuting property and approximation properties of the Fortin interpolant. Considering the stabilization part we see that using (6.3) on each boundary face followed by the approximation (3.2),

 (μ/h)12∥(u−πhu)⋅n∥∂Ω≤Chμ12|u|H2(Ω).

Using the definition of we see that and therefore

 (λ/h)12∥~π0(u−πhu)⋅n∥∂Ω=0.

This last property is necessary to prove that the method is locking free.

 (6.6)
###### Lemma 6.2

For there holds

 c(uh,uh)≤ϵ|||uh|||2E+ϵ−1C2Tγ−1|uh|2s. (6.7)

Proof. This proof follows the ideas of [14], we include it here for completeness. First we note that

 c(uh,uh)=(2μn⋅∇suhn+λ∇⋅uh,uh⋅n)∂Ω+(2μt⋅∇suhn,uh⋅t)∂ΩD.

Since for , there holds

 (λ∇⋅uh,uh⋅n)∂Ω=(λ∇⋅uh,~π0uh⋅n)∂Ω.

Applying the Cauchy–Schwarz inequality followed by the trace inequality (6.4) we see that for all ,

 (2μn⋅∇suhn,uh⋅n)∂Ω≤2CT∥μ12∇suh∥Ω∥μ12h−12uh⋅n∥∂Ω≤ϵ∥μ12∇suh∥2Ω+C2Tϵ−1∥μ12h−12uh⋅n∥2∂Ω
 (2t⋅μ∇suhn,uh⋅t)∂ΩD≤2CT∥μ12∇suh∥Ω∥μ12h−12uh⋅t∥∂ΩD≤ϵ∥μ12∇su