1 Introduction

# An HDG method for linear elasticity with strong symmetric stresses

## Abstract.

This paper presents a new hybridizable discontinuous Galerkin (HDG) method for linear elasticity on general polyhedral meshes, based on a strong symmetric stress formulation. The key feature of this new HDG method is the use of a special form of the numerical trace of the stresses, which makes the error analysis different from the projection-based error analyzes used for most other HDG methods. For arbitrary polyhedral elements, we approximate the stress by using polynomials of degree and the displacement by using polynomials of degree . In contrast, to approximate the numerical trace of the displacement on the faces, we use polynomials of degree only. This allows for a very efficient implementation of the method, since the numerical trace of the displacement is the only globally-coupled unknown, but does not degrade the convergence properties of the method. Indeed, we prove optimal orders of convergence for both the stresses and displacements on the elements. In the almost incompressible case, we show the error of the stress is also optimal in the standard norm. These optimal results are possible thanks to a special superconvergence property of the numerical traces of the displacement, and thanks to the use of a crucial elementwise Korn’s inequality. Several numerical results are presented to support our theoretical findings in the end.

###### Key words and phrases:
hybridizable; discontinuous Galerkin; superconvergence; linear elasticity
65N30, 65L12

## 1. Introduction

In this paper, we introduce a new hybridizable discontinuous Galerkin (HDG) method for the system of linear elasticity

 Aσ––−ϵ–(u) =0 in Ω⊂R3, (1.1a) ∇⋅σ–– =f in Ω, (1.1b) u =g on ∂Ω. (1.1c)

Here, the displacement is denoted by the vector field . The strain tensor is represented by . The stress tensor is represented by , where denotes the set of all symmetric matrices in . The compliance tensor is assumed to be a bounded, symmetric, positive definite tensor over . The body force lies in , the displacement of the boundary is a function in and is a polyhedral domain.

In general, there are two approaches to design mixed finite element methods for linear elasticity. The first approach is to enforce the symmetry of the stress tensor weakly ([4, 5, 11, 17, 25, 28, 32, 33, 36]). In this category, is included the HDG method considered in [22]. The other approach is to exactly enforce the symmetry of the approximate stresses. The methods considered in [21, 1, 2, 3, 7, 8, 26, 30, 35, 37, 38] belong to the second category, and so does the contribution of this paper. In general, the methods in the first category are easier to implement. On the other hand, the methods in the second category preserve the balance of angular momentum strongly and have less degrees of freedom. Next, we compare our HDG method with several methods of the second category.

In [21], an LDG method using strongly symmetric stresses (for isotropic linear elasticity) was introduced and proved to yield convergence properties that remain unchanged when the material becomes incompressible; simplexes and polynomial approximations or degree in all variables were used. However, as all LDG methods for second-order elliptic problems, although the displacement converges with order , the strain and pressure converge sub-optimally with order . Also, the method cannot be hybridized. Stress finite elements satisfying both strong symmetry and -conformity are introduced in [1, 2]. The main drawback of these methods is that they have too many degrees of freedom of stress elements and hybridization is not available for them (see detailed description in [28]). In [3, 7, 8, 26, 30, 35, 37, 38], non-conforming methods using symmetric stress elements are introduced. But, methods in [3, 7, 8, 30, 37, 38] use low order finite element spaces only (most of them are restricted to rectangular or cubical meshes except [3, 7]). In [26], a family of simplicial elements (one for each ) are developed in both two and three dimensions. (The degrees of freedom of were studied in [26] and then used to design the projection operator in [27]). However, the convergence rate of stress is suboptimal. The first HDG method for linear and nonlinear elasticity was introduced in [34, 35]; see also the related HDG method proposed in [39]. These methods also use simplexes and polynomial approximations of degree in all variables. For general polyhedral elements, this method was recently analyzed in [23] where it was shown that the method converges optimally in the displacement with order , but with the suboptimal order of for the pressure and the stress. For , these orders of convergence were numerically shown to be sharp for triangular elements. In this paper, we prove that by enriching the local stress space to be polynomials of degree no more than , and by using a modified numerical trace, we are able to obtain optimal order of convergence for all unknowns. In addition, this analysis is valid for general polyhedral meshes. To the best of our knowledge, this is so far the only result which has optimal accuracy with general polyhedral triangulations for linear elasticity problems.

Like many hybrid methods, our HDG method provides approximation to stress and displacement in each element and trace of displacement along interfaces of meshes. In general, the corresponding finite element spaces are , which are defined to be

 V––h= {v––∈L––2(Ω): v––|K∈V––(K) ∀K∈Th}, Wh= {ω∈L2(Ω): ω|K∈W(K) ∀K∈Th}, Mh= {μ∈L2(Eh): μ|F∈M(F) ∀F∈Eh}.

Here denotes a triangulation of the domain and is the set of all faces of all elements . The spaces are called the local spaces which are defined on each element/face. In Table 1 we list several choices of local spaces for different methods. In this paper, our choice of the local spaces is defined as:

 V––(K)=P––k(S,K),W(K)=Pk+1(K),M(F)=Pk(F).

Here, the space of vector-valued functions defined on whose entries are polynomials of total degree is denoted by (). Similarly, denotes the space of symmetric-valued functions defined on whose entries are polynomials of total degree . In addition, our method allows to be any conforming polyhedral triangulation of .

Note the fact that the only globally-coupled degrees of freedom are those of the numerical trace of displacement along , renders the method efficiently implementable. However, the fact that the polynomial degree of the approximate numerical traces of the displacement is one less than that of the approximate displacement inside the elements, might cause a degradation in the approximation properties of the displacement. However, this unpleasant situation is avoided altogether by taking a special form of the numerical trace of the stresses inspired on the choice taken in [29] in the framework of diffusion problems. This choice allows for a special superconvergence of part of the numerical traces of the stresses which, in turn, guarantees that, for , the -order of convergence for the stress is and that of the displacement . So, we obtain optimal convergence for both stress and displacement for general polyhedral elements. Let us mention that the approach of error analysis of our HDG method is different from the traditional projection-based error analysis in [19, 20, 22] in three aspects. First, here, we use simple -projections, not the numerical trace-tailored projections typically used for the analysis of other HDG methods. Second, we take the stabilization parameter to be of order instead of of order one. And finally, we use an elementwise Korn’s inequality (Lemma 4.1) to deal with the symmetry of the stresses.

We notice that mixed methods in [17, 25] and HDG methods in [22] also achieve optimal convergence for stress and superconvergence for displacement by post processing. However, there are two disadvantages regarding of implementation. First, these methods enforce the stress symmetry weakly, which means that they have a much larger space for the stress. In additon, these methods usually need to add matrix bubble functions ( in [17]) into their stress elements in order to obtain optimal approximations. In fact, the construction of such bubbles on general polyhedral elements is still an open problem. In contrast, our method avoids using matrix bubble functions but only use simple polynomial space of degree . In Table 1, we compare methods which use for approximating trace of displacement on . There, is a post-processed numerical solution of displacement.

The remainder of this paper is organized as follows. In Section , we introduce our HDG method and present our a priori error estimates. In Section , we give a characterization of the HDG method and show the global matrix is symmetric and positive definite. In Section , we give elementwise Korn’s inequality in Lemma 4.1, then provide a detailed proof of the a priori error estimates. In Section , we present several numerical examples in order to illustrate and test our method.

## 2. Main results

In this section we first present the method in details and then show the main results for the error estimates.

### 2.1. The HDG formulation with strong symmetry

Let us begin by introducing some notations and conventions. We adapt to our setting the notation used in [20]. Let denote a conforming triangulation of made of shape-regular polyhedral elements . We recall that , and denotes the set of all faces of all elements. We denote by the set of all faces of the element . We also use the standard notation to denote scalar, vector and tensor spaces. Thus, if denotes a space of scalar-valued functions defined on , the corresponding space of vector-valued functions is and the corresponding space of matrix-valued functions is . Finally, denotes the symmetric subspace of .

The methods we consider seek an approximation to the exact solution in the finite dimensional space given by

 V––h= {v––∈L––2(S,Ω): v––|K∈P––k(S,K) ∀K∈Th}, (2.1a) Wh= {ω∈L2(Ω): ω|K∈Pk+1(K) ∀K∈Th}, (2.1b) Mh= {μ∈L2(Eh): μ|F∈Pk(F) ∀F∈Eh}. (2.1c)

Here denotes the standard space of polynomials of degree no more than on . Here we require .

The numerical approximation can now be defined as the solution of the following system:

 (Aσ––h,v––)Th+(uh,∇⋅v––)Th−⟨ˆuh,v––n⟩∂Th =0, (2.2a) (σ––h,∇ω)Th−⟨ˆσ––hn,ω⟩∂Th =−(f,ω)Th, (2.2b) ⟨ˆσ––hn,μ⟩∂Th∖∂Ω =0, (2.2c) ⟨ˆuh,μ⟩∂Ω =⟨g,μ⟩∂Ω, (2.2d) for all (v––,ω,μ)∈–Vh×Wh×Mh, where ˆσ––hn=σ––hn−τ(PMuh−ˆuh)on ∂Th. (2.2e)

In fact, in Christoph Lehrenfeld’s thesis, the author defines the numerical flux in this way for diffusion problems (see Remark in [29]). This method was then analyzed for diffusion recently in [31]. Here, denotes the standard -orthogonal projection from onto . We write , , and where denotes the integral of over . Similarly, we write and , where denotes the integral of over .

The parameter in (2.2e) is called the stabilization parameter. In this paper, we assume it is a fixed positive number on all faces. It is worth to mention that the numerical trace (2.2e) is defined slightly different from the usual HDG setting, see [20]. Namely, in the definition, we use instead of . Indeed, this is a crucial modification in order to get error estimate. An intuitive explanation is that we want to preserve the strong continuity of the flux across the interfaces. Without the projection , by (2.2c) the normal component of is only weakly continuous across the interfaces.

### 2.2. A priori error estimates

To state our main result, we need to introduce some notations. We define

 ∥v––∥L2(A,Ω)=√(Av––,v––)Ω,∀v––∈L––2(S,Ω).

We use to denote the usual norm and semi-norm on the Sobolev space . We discard the first index if . A differential operator with a sub-index means it is defined on each element . Similarly, the norm is the discrete norm defined as . Finally, we need an elliptic regularity assumption stated as follows. Let be the solution of the adjoint problem:

 Aψ––−ϵ–(ϕ) =0 in Ω, (2.3a) ∇⋅ψ–– =eu in Ω, (2.3b) ϕ =0 on ∂Ω. (2.3c)

We assume the solution has the following elliptic regularity property:

 ∥ψ––∥1,Ω+∥ϕ∥2,Ω≤Creg∥eu∥Ω, (2.4)

The assumption holds in the case of planar elasticity with scalar coefficients on a convex domain, see [9].

We are now ready to state our main result.

###### Theorem 2.1.

If the meshes are quasi-uniform and , then we have

 ∥σ––−σ––h∥L2(A,Ω)≤Chs(∥u∥s+1,Ω+∥σ––∥s,Ω), (2.5)

for all . Moreover, if the elliptic regularity property (2.4) holds, then we have

 ∥u−uh∥Ω≤Chs+1(∥u∥s+1,Ω+∥σ––∥s,Ω), (2.6)

for all . Here the constant depends on the upper bound of compliance tensor but it is independent of the mesh size .

This result shows that the numerical errors for both unknowns are optimal. In addition, since the only globally-coupled unknown, , stays in , the order of convergence for the displacement remains optimal only because of a key superconvergence property, see the remark right after Corollary 4.2. In addition, we restrict our result on quasi-uniform meshes to make the proof simple and clear. This result holds for shape-regular meshes also.

### 2.3. Numerical approximation for nearly incompressible materials

Here, we consider the numerical approximation of stress for isotropic nearly incompressible materials.

We define isotropic materials to be those whose compliance tensor satisfying the following Assumption 2.1.

###### Assumption 2.1.
 Aτ–– =PDτ––D+PTtr(τ––)3I3 (2.7) where τ––D =τ––−tr(τ––)3I3,

for any in , and and are two positive constants. An isotropic material is nearly incompressible if is close to zero.

###### Theorem 2.2.

If the material is isotropic (whose compliance tensor satisfies Assumption 2.1), is positive, the boundary data , the meshes are quasi-uniform and , then we have

 ∥σ––−σ––h∥L2(Ω)≤Chs(∥u∥s+1,Ω+∥σ––∥s,Ω), (2.8)

for all . Here, the constant is independent of .

This result shows that the HDG method (2.2) is locking-free for nearly incompressible materials. We emphasize that the convergence rate of stress for nearly incompressible materials is one order higher than [5, 26] with the same finite element space for numerical trace of displacement.

## 3. A characterization of the HDG method

In this section we show how to eliminate elementwise the unknowns and from the equations (2.2) and rewrite the original system solely in terms of the unknown , see also [35]. Via this elimination, we do not have to deal with the large indefinite linear system generated by (2.2), but with the inversion of a sparser symmetric positive definite matrix of remarkably smaller size.

### 3.1. The local problems

The result on the above mentioned elimination can be described using additional “local” operators defined as follows:

On each element , for any , we denote to be the unique solution of the local problem:

 (AQ––λ,v––)K+(Uλ,∇⋅v––)K =⟨λ,v––⋅n⟩∂K, (3.1a) −(∇⋅Q––λ,ω)K+⟨τPMUλ,ω⟩∂K =⟨τλ,ω⟩∂K, (3.1b)

for all .

On each element , we also denote to be the unique solution of the local problem:

 (AQ––Sf,v––)K+(USf,∇⋅v––)K =0, (3.2a) −(∇⋅Q––Sf,ω)K+⟨τPMUSf,ω⟩∂K =−(f,ω)K, (3.2b)

for all .

It is easy to show the two local problems are well-posted. In addition, due to the linearity of the global system (2.2),the numerical solution satisfies

 σ––h=Q––ˆuh+Q––Sf,uh=Uˆuh+USf. (3.3)

### 3.2. The global problem

For the sake of simplicity, we assume the boundary data . Then, the HDG method (2.2) is to find satisfying

 (Aσ––h,v––)Th+(uh,∇⋅v––)Th−⟨ˆuh,v––n⟩∂Th =0, (3.4a) −(∇⋅σ––h,ω)Th+⟨τ(PMuh−ˆuh),ω⟩∂Th =−(f,ω)Th, (3.4b) ⟨σ––hn−τ(PMuh−ˆuh),μ⟩∂Th∖∂Ω =0, (3.4c)

for all , where .

Combining (3.4c) with (3.3), we have that for all ,

 ⟨(Q––ˆuh)n−τ(PMUˆuh−ˆuh),μ⟩∂Th=⟨(Q––Sf)n−τPMUSf,μ⟩∂Th. (3.5)

Up to now we can see that we only need to solve the reduced global linear system (3.5) first, then recover by (3.3) element by element. Next we show that the global system (3.5) is in fact symmetric positive definite.

### 3.3. A characterization of the approximate solution

The above results suggest the following characterization of the numerical solution of the HDG method.

###### Theorem 3.1.

The numerical solution of the HDG method (2.2) satisfies

 σ––h=Q––ˆuh+Q––Sf,uh=Uˆuh+USf.

If we assume the boundary data , then is the solution of

 ah(ˆuh,μ)=⟨(Q––Sf)n−τPMUSf,μ⟩∂Th,∀μ∈M0h, (3.6)

where

 ah(ˆuh,μ)=(AQ––ˆuh,Q––μ)Th+⟨τ(PMUˆuh−ˆuh),PMUμ−μ⟩∂Th.

In addition, the bilinear operator is positive definite.

###### Proof.

In order to show (3.6) is true, we only need to show that for all , then

 ah(λ,μ)=⟨(Q––λ)n−τ(PMUλ−λ),μ⟩∂Th.

According to (3.1), we have

 (AQ––m,v––)Th+(Um,∇⋅v––)Th=⟨m,–v⋅n⟩∂Th, (3.7a) (∇⋅Q––m,ω)Th=⟨τ(PMUm−m),ω⟩∂Th, (3.7b)

for all , . Then, we have

 ⟨(Q––λ)n−τ(PMUλ−λ),μ⟩∂Th = ⟨μ,(Q––λ)n⟩∂Th−⟨τ(PMUλ−λ),μ⟩∂Th = (AQ––μ,–Qλ)Th+(Uμ,∇⋅Q––λ)Th−⟨τ(PMUλ−λ),μ⟩∂Th by (???) = Extra open brace or missing close brace = (AQ––μ,–Qλ)Th+⟨τ(PMUλ−λ),Uμ−μ⟩∂Th by (???) = (AQ––μ,–Qλ)Th+⟨τ(PMUλ−λ),PMUμ−μ⟩∂Th = ah(λ,μ).

So, we can conclude that (3.6) holds. We end the proof by showing the bilinear operator is positive definite.

If for some , from the previous result we have

 Q––λ=0,PMUλ−λ|∂Th=0.

We apply integration by parts on (3.1a), we have

 ⟨ϵ–(Uλ),–v⟩∂K=0,∀v––∈V––(K).

This implies that for all . So, for any , there are such that . Since , we have . Combining this result with the fact that and , we can conclude that and .

Finally, let us consider two adjacent element with the interface . In addition, we assume that on , can be expressed as

 Uλ=ai×x+bi,i=1,2.

We claim that and . This fact can be shown by considering the continuity of the function on the interface . We omit the detailed proof since it only involves elementary linear algebra.

From this result we conclude that there exist such that in . By the fact that , we can conclude that , hence . This completes the proof. ∎

###### Remark 3.2.

In Theorem 3.1, we assume the boundary data . Actually, if is not zero, we can still obtain the same linear system as in Theorem 3.1 by the same treatment of boundary data in [16].

## 4. Error Analysis

In this section we provide detailed proofs for our a priori error estimates - Theorem 2.1 and Theorem 2.2. We use elementwise Korn’s inequality (Lemma 4.1), which is novel and crucial in error analysis. We use to denote the standard -orthogonal projection onto , respectively. In addition, we denote

 eσ–––=ΠV–––––σ−σ––h,eu=ΠWu−uh,eˆu=PMu−ˆuh,

In the analysis, we are going to use the following classical results:

 ∥u−ΠWu∥Ω ≤Chs∥u∥s,Ω 0≤s≤k+2, (4.1a) ∥σ––−ΠV––––σ––∥Ω ≤Cht∥σ––∥t,Ω 0≤t≤k+1, (4.1b) ∥u−PMu∥Eh ≤Chs−12∥u∥s,Ω, 1≤s≤k+1, (4.1c) ∥u−ΠWu∥∂K ≤Chs−12∥u∥s,K, 1≤s≤k+2, (4.1d) ∥σ––n−ΠV––––σ––n∥∂K ≤Cht−12∥σ––∥t,K, 1≤t≤k+1, (4.1e) ∥v∥∂K ≤Ch−12∥v∥K, ∀v∈Ps(K), (4.1f) ∥σ––n−PM(σ––n)∥∂K ≤Cht−12∥σ––∥t,K, 1≤t≤k+1. (4.1g)

The above results are due to standard approximation theory of polynomials, trace inequality.

Let denote the discrete symmetric gradient operator, such that for any , . It is well known (see Theorem in [14]) the kernel of the operator is:

 kerϵ–h=Υh:={Λ∈L2(Ω),Λ|K=B––Kx+bK,B––K∈A–––,bK∈R3,K∈Th}.

Here, denotes the set of all anti-symmetric matrices in .

In the analysis, we need the following elementwise Korn’s inequality:

###### Lemma 4.1.

Let be a generic element with size and . Then for any function , we have

 infΛ∈Υ(K)∥∇(v+Λ)∥K≤C∥ϵ–(v)∥K,

Here is independent of the size . In addition, if is a tetrahedron, the above inequality holds for any .

###### Proof.

Let denote the reference tetrahedron element and . The mapping from to is where is a non-singular matrix and .

We define , which is the pull back of on , by

 A––−⊤Kˆv(ˆx)=v(x)∀ˆx∈ˆK.

So, we have

 ∇v(x)=∇(A––−⊤Kˆv)(x)=A––−⊤K(∇ˆv)(x).

The last equality above is due to the fact that every component of is constant. It is easy to see that

 (∇ˆv)(x)=ˆ∇ˆv(ˆx)A––−1K.

So, we have

 A––−⊤Kˆ∇ˆv(ˆx)A––−1K=∇v(x).

By taking the symmetric part of both sides of the above equation, we have

 A––−⊤Kˆϵ–ˆv(ˆx)A––−1K=ϵ–v(x). (4.2)

According to Theorem in [14], the following inequality holds:

 infˆΛ∈Υ(ˆK)∥ˆv+ˆΛ∥1,ˆK≤C∥ˆϵ–(ˆv)∥0,ˆK.

So, there is with and , such that

 ∥ˆ∇(ˆv+ˆΛ)∥0,ˆK≤C∥ˆϵ–(ˆv)∥0,ˆK. (4.3)

We define

 Λ(x)=A––−⊤KˆΛ(ˆx)∀x∈K.

It is easy to see that

 ∇Λ=A––−⊤Kˆ∇ˆΛA––−1K=A––−⊤KB––ˆKA––−1K∈A–––.

So, . Then, by standard scaling argument with (4.2, 4.3) and the shape regularity of the meshes, we can conclude that the proof for arbitrary tetrahedron element is complete.

Now, we consider the case of arbitrary shape regular element , which can be hexahedron, prism or pyramid. Let . It is well known that for any ,

 ∂j(∂kvi)=∂j(ϵik(v))+∂k(ϵij(v))−∂i(ϵjk(v)).

Here, . Consequently, we have

 ∥∇(∂jvi−∂ivj)∥0,K≤C∥∇ϵ–(v)∥0,K≤Ch−1K∥ϵ–(v)∥