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 . They may typically also exhibit instability when Darcy
flow is considered if the element was designed for Stokes’ problem .
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  . 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
. 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
where the tensor product is defined by and
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
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:
Here the bilinear forms are given by
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
A unique pressure is then guaranteed by the Ladyzhenskaya-Babuska-Brezzi condition .
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,
The -projection on , is defined by for all . satisfies the stability for all and the approximation error estimate
We also introduce the -projection of the trace of a function
where is the set of faces in such that . We let denote the space of vectorial piecewise affine functions on ,
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 , 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
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  , 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.
, 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
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
of each triangle of type I is assigned a unit vectoralong 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
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
This gives equal divergence on all subtriangles.
3.2 The construction of 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 , 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  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
denotes the Clément interpolant, and for all
Note that the interpolant satisfies the approximation error estimate
The proof of the existence of is identical to that of the interpolant for the Bernardi-Raugel element . 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
A consequence of the existence of the Fortin interpolant is the existence of a non-trivial subspace such that
As a consequence, for every there exists
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
5 Stability and error analysis
We introduce two triple norms. First for the elasticity system,
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,
For the problem (4.1) Korn’s inequality leads to the coercivity, there exists such that for all
For the problem (4.2) we need to prove an inf-sup condition for stability.
(inf-sup stability for the Brinkman problem) There exists such that for all there holds
Proof. First we take and to obtain
Then we chose , where is defined by (3.3) so that
Observing now that
it follows that
Taking and we conclude that
To finish the proof note that
Proof. Let , with . Note that by adding and subtracting and using the triangle inequality and Korn’s inequality we have
For the second term we apply the coercivity (5.3), followed by Galerkin orthogonality
we may write
which proves the first claim.
Proof. We introduce, as before, discrete errors , with and . Using the triangle inequality we see that
For the second term in the right hand side we apply the stability of Proposition 5.1 to obtain
using Galerkin orthogonality
Observe that by construction we have
The only remaining term in the right hand side of (5.6) is bounded using the Cauchy–Schwarz inequality,
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
Here the tangential projection is defined by . The Nitsche formulation then takes the form: Find such that
Observe that the projection in the boundary penalty of the normal component is necessary to avoid locking.
We define the stabilization semi-norm by
and the following augmented energy norm defined on
We recall that is a norm by Korn’s inequality and Poincaré’s inequality. We recall the trace inequalities
Using these inequalities it is straightforward to prove the following approximation estimate in the norm and a bound on the form .
The following approximation inequality holds
Proof. The inequality
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),
Using the definition of we see that and therefore
This last property is necessary to prove that the method is locking free.
For there holds