A mixed finite element method for nearly incompressible elasticity and Stokes equations using primal and dual meshes with quadrilateral and hexahedral grids

# A mixed finite element method for nearly incompressible elasticity and Stokes equations using primal and dual meshes with quadrilateral and hexahedral grids

Bishnu P. Lamichhane School of Mathematical and Physical Sciences, University of Newcastle, Callaghan, NSW 2308, Bishnu.Lamichhane@anu.edu.au
###### Abstract

We consider a mixed finite element method for approximating the solution of nearly incompressible elasticity and Stokes equations. The finite element method is based on quadrilateral and hexahedral triangulation using primal and dual meshes. We use the standard bilinear and trilinear finite element space enriched with element-wise defined bubble functions with respect to the primal mesh for the displacement or velocity, whereas the pressure space is discretised by using a piecewise constant finite element space with respect to the dual mesh.

Key words. mixed finite elements, nearly incompressible elasticity, primal and dual meshes, Stokes equations, inf-sup condition

AMS subject classifications. 65N30, 65N15, 74B10

## 1 Introduction

Although there are many mixed finite element methods for nearly incompressible elasticity and Stokes equations leading to an optimal convergence, the search for simple, efficient and optimal finite element schemes is still an active area of research. In this article we present a mixed finite element method for nearly incompressible elasticity and Stokes equations using quadrilateral and hexahedral meshes. The displacement or velocity field is discretised by using the standard bilinear or trilinear finite element space enriched with element-wise defined bubble functions, whereas the pressure space is discretised by the piecewise constant finite element space based on a dual mesh. Such a finite element space for the simplicial mesh is presented in [12], where the inf-sup condition is proved by using the fact that the mini finite element [1] satisfies the inf-sup condition. Note that the mini finite element [1] consists of the linear finite element space enriched with element-wise defined bubble functions for the displacement or velocity and the linear finite element space for the pressure space. The enrichment of the displacement or velocity field increases one vector degree of freedom per element. A main hindrance to extend this approach to the case of quadrilateral and hexahedral meshes is that the displacement or velocity space should be enriched by more than a single bubble function to obtain the inf-sup condition [2]. In this article we show that a similar discretisation scheme can be applied to quadrilateral and hexahedral meshes. We prove that if the pressure space is discretised by using the piecewise constant function space with respect to the dual mesh, it is sufficient to enrich the standard bilinear and trilinear finite element space with a single bubble function per element.

## 2 The boundary value problem of linear elasticity

We introduce the boundary value problem of linear elasticity in this section. In particular, we present the standard weak formulation and a mixed formulation of a linear elastic problem. We consider a homogeneous isotropic linear elastic material body occupying a bounded domain , , with Lipschitz boundary . For a prescribed body force , the governing equilibrium equation in reads

 −div\boldmath{σ}=\boldmath{f}, \hb@xt@.01(2.1)

where is the symmetric Cauchy stress tensor. The stress tensor is defined as a function of the displacement by the Saint-Venant Kirchhoff constitutive law

 \boldmath{σ}=12C(∇% \boldmath{u}+[∇\boldmath{u}]t), \hb@xt@.01(2.2)

where is the fourth-order elasticity tensor. The action of the elasticity tensor on a tensor is defined as

 \boldmath{σ}=C\boldmath{d}:=λ(tr\boldmath{d})1+2μ\boldmath{d}. \hb@xt@.01(2.3)

Here, is the identity tensor, and and are the Lamé parameters, which are constant in view of the assumption of a homogeneous body, and they are assumed to be positive. For simplicity of exposition we assume that the displacement or velocity satisfies homogeneous Dirichlet boundary condition

 \boldmath{u}=0onΓ. \hb@xt@.01(2.4)

However, the approach works also for mixed boundary conditions.

Standard weak formulation.
Let be the set of square-integrable functions defined on , where the inner product and norm on this space is denoted by and , respectively. The Sobolev space is defined in terms of the space as

 H1(Ω)={u∈L2(Ω),∇u∈[L2(Ω)]d},

and , where a function in vanishes on the boundary in the sense of traces. The space is the subset of defined as

 L20(Ω)={p∈L20(Ω):∫Ωpd% \boldmath{x}=0}.

To write the weak or variational formulation of the boundary value problem, we introduce the space of displacement or velocity with inner product and norm defined in the standard way; that is, , with the norm being induced by this inner product.

We define the bilinear form and the linear functional by

 A:\boldmath{V}×\boldmath{V}→\twelvemsb R,A(\boldmath{u},\boldmath{v}):=∫ΩC\boldmath{ε}(\boldmath{u}):\boldmath{ε}(\boldmath{v}) d\boldmath{% x},ℓ:V→\twelvemsb R,ℓ(\boldmath{v}):=∫Ω\boldmath{f}⋅\boldmath{v} d\boldmath{x}.

Then the standard weak form of linear elasticity problem is as follows: given , find that satisfies

 A(\boldmath{u},\boldmath{v})=ℓ(\boldmath{v}) ,\boldmath{v}∈\boldmath{V}. \hb@xt@.01(2.5)

The assumptions on guarantee that is symmetric, continuous, and -elliptic. Hence by using standard arguments it can be shown that (2.5) has a unique solution . Furthermore, if the the domain is convex with polygonal or polyhedral boundary, , and there exists a constant independent of such that [3, 15, 9]

 ∥\boldmath{u}∥2+λ∥div% \boldmath{u}∥1≤C∥\boldmath{f}∥0. \hb@xt@.01(2.6)

Mixed formulation. There are many mixed formulation for the linear elasticity problem. The simplest one is given by introducing pressure as an extra variable, which leads to penalized Stokes equations. Defining , a mixed variational formulation of linear elastic problem (2.5) is given by: find such that

 a(\boldmath{u},\boldmath{v})+b(\boldmath{v},p)=ℓ(\boldmath{v}),\boldmath{v}∈\boldmath{V},b(\boldmath{v},q)−1λc(p,q)=0,q∈L20(Ω), \hb@xt@.01(2.7)

where

 a(\boldmath{u},\boldmath{v}):=2μ∫Ω\boldmath{ε}(\boldmath{u}):\boldmath{ε}(\boldmath{v})d\boldmath{x},b(% \boldmath{v},q):=∫Ωdiv\boldmath% {v}qd\boldmath{x},andc(p,q):=∫Ωpqd\boldmath{x}.

Using a standard theory of mixed finite elements [4], the existence and uniqueness of the solution of the problem (2.7) can be shown. Of particular interest for us is the incompressible limit, which corresponds to . We note that the saddle point equations (2.7) reduce to Stokes equations of fluid flow when . In this case, denotes the kinematic viscosity, and represents the velocity of the fluid.

## 3 Finite element discretizations

We consider a quasi-uniform triangulation - called the primal mesh - of the polygonal or polyhedral domain , where consists of convex quadrilaterals or hexahedras. The finite element meshes are defined by maps from a reference square or reference cube .

For nonnegative integer , we let denote the space of polynomials in two or three variables of total degree less than or equal to , and the space of polynomials in two variables of total degree less than or equal to in each variable. A typical element is generated by an iso-parametric map from the reference element , in which is defined using the basis functions corresponding to . It is clear that if , then is in general not a polynomial on the quadrilateral . However, in the following we assume that the map is affine for all .

The finite element space of displacements is taken to be the space of continuous functions whose restrictions to an element are obtained by maps of bilinear or trilinear functions from the reference element:

 Sh:={vh∈H10(Ω), vh|K=^vh∘F−1K,^\boldmath{v}h∈Q1(^K),  K∈Th} . \hb@xt@.01(3.1)

Let be the number of vertices in , and the set of all vertices in be denoted by

 Nh:={\boldmath{x}1,\boldmath{x}2,⋯,\boldmath{x}N}.

A dual mesh is introduced based on the primal mesh so that the elements of are called control volumes. Each control volume element is associated with a vertex . For simplicity we explain the construction of the dual mesh for quadrilateral meshes. The idea can be extended to hexahedral meshes in the standard way as the extension from triangular meshes to tetrahedral meshes in [6, 12]. For a vertex let and be the set of elements and edges touching , respectively. Then the volume element corresponding to the vertex is the polygonal domain joining centroids of all elements in and centroids of all edges in . If is a boundary vertex, all the boundary edges touching will also form the boundary of . The set of all volume elements in will form a non-overlapping decomposition of the polygonal domain . We refer to [6, 12] for a similar construction in simplicial case. A dual mesh for a quadrilateral grid is shown in Figure 3.1.

In the following we use a generic constant which will take different values at different places but will always be positive and independent of the mesh-size . We call the control volume mesh regular or quasi-uniform if there exists a positive constant such that

 Chd≤|Vi|≤hd,Vi∈T∗h,

where is the maximum diameter of all elements . It can be shown that, if is locally regular, i.e., there is a constant such that

 ChdT≤|T|≤hdT,T∈Th

with for all elements , then this dual mesh is also locally regular. The dual volume element space to discretize the pressure is now defined by

 S∗h:={p∈L20(Ω):p|V∈P0(V),V∈T∗h}.

Now any element and can be written as and , where are the standard nodal basis functions associated with the vertex , and are the characteristic functions of the volume . Let with on and , where is the centroid of , be a bubble function corresponding to the element . Let , where is the standard linear or trilinear basis function corresponding to the reference element associated with the origin. Defining the space of bubble functions

 \boldmath{B}h:={\boldmath{b}h∈[C0(Ω)]d:\boldmath{b}h|T=cT∇ϕTbT,cT∈R,T∈Th}, \hb@xt@.01(3.2)

we introduce our finite element space for displacement or velocity as .

Then, the finite element approximation of (2.7) is defined as a solution to the following problem: find such that

 a(\boldmath{u}h,\boldmath{v}h)+b(\boldmath{v}h,ph)=ℓ(\boldmath{v}h),% \boldmath{v}h∈\boldmath{V}h,b(\boldmath{u}h,qh)−1λc(ph,qh)=0,qh∈S∗h. \hb@xt@.01(3.3)

To establish a priori estimates for the discretization errors, we consider the saddle point formulation (3.3) of the elasticity problem and apply the theory of mixed finite elements. The continuity of the bilinear form on , of on and of on is straightforward. By using the Korn’s inequality, it is standard that the ellipticity of the bilinear form holds on . It remains to show that the uniform inf-sup condition holds for the bilinear form on . That means, for any , there exists a constant independent of the mesh-size such that

 sup\boldmath{v}h∈\boldmath{V}hb(\boldmath{v}h,qh)∥\boldmath{v}h∥1≥β∥qh∥0. \hb@xt@.01(3.4)

We note that we have proved optimal a priori error estimate for the finite element scheme (3.3) in [12] for simpilical meshes using the stability of the mini finite element [1], where the standard linear finite element space is enriched by an element-wise defined bubble function per element. For quadrilateral and hexahedral meshes the stability is not attained by enriching the standard bilinear or trilinear finite element space by a single element-wise defined bubble function, see [2]. Therefore, we recourse to another method to prove the stability of the scheme here. We prove the inf-sup condition (3.4) using a domain decomposition technique as in [14].

In the following we assume that we have a decomposition of in disjoint subdomains , where each subdomain consists of four quadrilaterals or eight hexahedra touching the vertex , see the right picture of Figure 3.1 for the quadrilateral case. Let be the restriction of the finite element space to the set satisfying the homogeneous Dirichlet boundary condition on the boundary of . First we observe that the necessary condition for the patch test is satisfied as the velocity space has degrees of freedom in two dimensions and in three dimensions, whereas the pressure space on has degrees of freedom in two dimensions and in three dimensions after excluding the constant functions on each . The proof of the following lemma can be obtained by a direct computation on one .

###### Lemma 3.1

The dimension of the space

 Bi={qh∈S∗h:b(\boldmath{v}h,qh)=0,%\boldmath$v$h∈\boldmath{V}ih}

is one.

Proof. We outline the proof in two dimensions. Each bubble function yields two equations leading to eight equations. Due to the symmetry we get linearly independent seven equations. We get additional two equations by using the vertex basis function associated with the vertex . However, only one of them is linearly independent to the previous seven equations as constant functions are in . Thus we have eight linearly independent equations leading to the fact that contains only constant functions. Note that the factor in the definition of the space of bubble functions in (3.2) is used to get that there are eight linearly independent equations. If we use the standard bubble function for discretising all components of the displacement, there will be only seven linearly independent equations.

Thus we have the following lemma.

###### Lemma 3.2

There exists a constant independent of the mesh-size such that

 sup\boldmath{v}h∈\boldmath{V}ihb(\boldmath{v}h,qh)∥\boldmath{v}h∥1≥C√∫Siq2hd\boldmath{x},qh∈L20(Si)∩S∗h,i=1,⋯,M. \hb@xt@.01(3.5)

This lemma is used to prove the inf-sup condition.

###### Theorem 3.3

There exists independent of the mesh-size such that

 sup\boldmath{v}h∈\boldmath{V}hb(\boldmath{v}h,qh)∥\boldmath{v}h∥1≥β∥qh∥0,qh∈S∗h.

Proof. Let

 qih=1|Si|∫Siqhd\boldmath{x},and~qh=M∑i=1qihχSi,

where is the characteristic function of the set . Then for

 sup\boldmath{v}h∈\boldmath{V}hb(\boldmath{v}h,qh)∥\boldmath{v}h∥1=sup\boldmath{v}h∈\boldmath{V}hb(\boldmath{% v}h,qh−~qh)∥\boldmath{v}h∥1+sup%\boldmath$v$h∈\boldmath{V}hb(\boldmath{v}h,~qh)∥\boldmath{v}h∥1

Note that

 sup\boldmath{v}h∈\boldmath{V}hb(\boldmath{v}h,qh−~qh)∥\boldmath{v}h∥1≥M∑i=1sup\boldmath{v}h∈\boldmath{V% }ihb(\boldmath{v}h,qh−~qh)∥% \boldmath{v}h∥1.

Since is supported only on for , we get

 sup\boldmath{v}h∈\boldmath{V}hb(\boldmath{v}h,qh−~qh)∥\boldmath{v}h∥1≥C∥qh−~qh∥0

from Lemma 3.5, and since is a piecewise constant function associated with one level coarser mesh, we have

 sup\boldmath{v}h∈\boldmath{V}hb(\boldmath{v}h,~qh)∥\boldmath{v}h∥1≥C∥~qh∥0.

Note that the final result follows from the fact that

 ∥qh−~qh∥0+∥~qh∥0≥C∥qh∥0.

The immediate consequence of the above discussion is the well-posedness of the discrete problem (3.3). From the theory of saddle point problem, see, e.g., [4], we have the following theorem.

###### Theorem 3.4

The discrete problem (3.3) has exactly one solution , which is uniformly stable with respect to the data , and there exists a constant independent of Lamé parameter such that

 ∥\boldmath{u}h∥1+∥ph∥0≤C∥% \boldmath{f}∥0.

The convergence theory is provided by an abstract result about the approximation of saddle point problems, see [4].

###### Theorem 3.5

Assume that and be the solutions of problems (2.7) and (3.3), respectively. Then, we have the following error estimate uniform with respect to :

 ∥\boldmath{u}−\boldmath{u}h∥1+∥p−ph∥0≤C(inf\boldmath{v}h∈\boldmath{V}h∥\boldmath{u}−\boldmath{v}h∥1+infqh∈S∗h∥p−qh∥0). \hb@xt@.01(3.6)

Since the space contains the space of piece-wise linear polynomials and contains the space of piece-wise constant functions with respect to the dual mesh , Theorem 3.5 yields the linear convergence of the discrete solution with respect to the energy norm for the displacement and with respect to the -norm for the pressure.

We note that the displacement field on the primal mesh and the stress field in the dual mesh in the Hellinger-Reissner problem of finding such that

 ∫ΩC−1\boldmath{σ}h:\boldmath{τ}h d\boldmath{x}−∫Ω% \boldmath{ε}(\boldmath{u}h):\boldmath{τ}h d\boldmath{x}=0,\boldmath{τ}h∈% \boldmath{S}∗h,∫Ω\boldmath{ε}(\boldmath{v}h):% \boldmath{σ}h d\boldmath{x}=ℓ(\boldmath{v}h) d\boldmath{x},\boldmath{v}h∈\boldmath{V% }h,

where , we arrive at the node-based uniform strain elements [7, 5]. However, the formulation is not stable [11].

Now we briefly describe how a displacement-based formulation is achieved for a nearly incompressible elasticity problem. From the second equation of (3.3), we can write with

 pi=λ|Vi|∫Vi∇⋅\boldmath{u}hd\boldmath{x}.

Hence, after condensing out the pressure from the formulation, we arrive at a problem of finding so that

 a(\boldmath{u}h,\boldmath{v}h)+N∑i=1λ|Vi|(∫Vi∇⋅\boldmath{u}hd\boldmath{x})(∫Vi∇⋅\boldmath{v}hd\boldmath{x})=ℓ(\boldmath{v}h),% \boldmath{v}h∈\boldmath{V}h.

If we look at the algebraic formulation of the finite element scheme for a nearly incompressible elasticity problem, we have

 (ABTB−1λC)(% \boldmath{u}hph)=(\boldmath{f}h0), \hb@xt@.01(3.7)

where and are matrices associated with the bilinear forms , and , respectively, and is the discrete vector associated with the linear form . The important thing to note is that the matrix is diagonal and the degrees of freedom corresponding to the pressure can easily be condensed out from the system leading to a positive definite formulation.

## 4 Numerical Results

### Example 1: Cook’s membrane problem

Our first test example is the popular benchmark problem known as Cook’s membrane problem [13, 10, 8]. Let be the quadrilateral connecting four points

 {(0,0),(48,44),(48,60),(0,44)}.

The left boundary of is clamped, and the right one is subjected to an in-plane shearing load of 100N along the -direction, as shown in the left picture of Figure 4.1. The material properties are taken to be and , so that a nearly incompressible response is obtained.

We have presented the vertical tip displacements at the point T computed using the mixed formulation and the standard displacement formulation are presented in the right picture of Figure 4.1, for different levels of uniform refinement, where the computation is started with the initial triangulation shown in the left picture of Figure 4.1. As can be seen from the right picture of Figure 4.1, the standard displacement approach exhibits extreme locking whereas the new mixed formulation shows rapid convergence.

### Example 2: Rectangular beam

The second example is concerned with a linear elastic beam of rectangular size subjected to a couple at one end, as shown in Figure 4.2. Along the edge , the horizontal displacement and vertical surface traction are zero. At the point , the vertical displacement is also zero. The exact solution is given by

 u(x,y)=2f(1−ν2)Elx(l2−y), and   v(x,y)=f(1−ν2)El[x2+ν1−νy(y−l)].

We set , , , , and . We have shown the setting of the problem in Figure 4.2, and the discretization errors with respect to the number of elements are presented in Figure 4.3. As can be seen from Figure 4.3, the standard approach locks completely, whereas we get very good numerical approximations with our new mixed formulation.

## 5 Conclusion

We have presented a finite element approach based on primal and dual meshes using quadrilateral and hexahedral meshes to approximate the solution of nearly incompressible elasticity or Stokes equations. Working with the space of bilinear or trilinear finite elements enriched with bubble functions for the displacement or velocity field we have proved the uniform inf-sup condition. As we have an orthogonal basis for the piecewise constant finite element space on the dual mesh, we can statically condense out the pressure variable from the system leading to a displacement-based formulation. The resulting displacement-based formulation is symmetric and positive-definite.

## References

• [1] D.N. Arnold, F. Brezzi, and M. Fortin. A stable finite element for the Stokes equations. Calcolo, 21:337–344, 1984.
• [2] W. Bai. A quadrilateral ’mini’ finite element for the Stokes problem. Computer Methods in Applied Mechanics and Engineering, 143:41–47, 1997.
• [3] S.C. Brenner and L. Sung. Linear finite element methods for planar linear elasticity. Mathematics of Computation, 59:321–338, 1992.
• [4] F. Brezzi and M. Fortin. Mixed and hybrid finite element methods. Springer–Verlag, New York, 1991.
• [5] C.R. Dohrmann, M.W. Heinstein, J. Jung, S.W. Key, and W.R. Witkowski. Node-based uniform strain elements for three-node triangular and four-node tetrahedral meshes. International Journal for Numerical Methods in Engineering, 47:1549–1568, 2000.
• [6] R. E. Ewing, T. Lin, and Y. Lin. On the accuracy of the finite volume element method based on piecewise linear polynomials. SIAM Journal on Numerical Analysis, 39:1865–1888, 2002.
• [7] D.P. Flanagan and T. Belytschko. A uniform strain hexahedron and quadrilateral with orthogonal hourglass control. International Journal for Numerical Methods in Engineering, 17:679–706, 1981.
• [8] E.P. Kasper and R.L. Taylor. A mixed-enhanced strain method. Part I: geometrically linear problems. Computers and Structures, 75:237–250, 2000.
• [9] V.A. Kozlov, V.G. Maz’ya, and J. Rossmann. Spectral Problems Associated with Corner Singularities of Solutions to Elliptic Equations. Mathematical Surveys and Monographs 85. American Mathematical Society, Providence, RI, 2001.
• [10] M. Küssner and B.D. Reddy. The equivalent parallelogram and parallelepiped, and their application to stabilized finite elements in two and three dimensions. Computer Methods in Applied Mechanics and Engineering, 190:1967–1983, 2001.
• [11] B.P. Lamichhane. From the Hu–Washizu formulation to the average nodal strain formulation. Computer Methods in Applied Mechanics and Engineering, 198:3957–3961, 2009.
• [12] B.P. Lamichhane. Inf-sup stable finite element pairs based on dual meshes and bases for nearly incompressible elasticity. IMA Journal of Numerical Analysis, 29:404–420, 2009.
• [13] J.C. Simo and M.S. Rifai. A class of assumed strain method and the methods of incompatible modes. International Journal for Numerical Methods in Engineering, 29:1595–1638, 1990.
• [14] R. Stenberg. A technique for analysing finite element methods for viscous incompressible flow. nternational Journal for Numerical Methods in Fluids, 11:935–948, 1990.
• [15] M. Vogelius. An analysis of the p-version of the finite element method for nearly incompressible materials. Uniformly valid, optimal error estimates. Numerische Mathematik, 41:39–53, 1983.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters