An HDG method for linear elasticity with strong symmetric stresses
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
2000 Mathematics Subject Classification:65N30, 65L12
In this paper, we introduce a new hybridizable discontinuous Galerkin (HDG) method for the system of linear elasticity
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 . 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 , 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 ). 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 , a family of simplicial elements (one for each ) are developed in both two and three dimensions. (The degrees of freedom of were studied in  and then used to design the projection operator in ). 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 . These methods also use simplexes and polynomial approximations of degree in all variables. For general polyhedral elements, this method was recently analyzed in  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
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:
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  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  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 ) 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 . 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
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:
|for all , where|
In fact, in Christoph Lehrenfeld’s thesis, the author defines the numerical flux in this way for diffusion problems (see Remark in ). This method was then analyzed for diffusion recently in . 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 . 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
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:
We assume the solution has the following elliptic regularity property:
The assumption holds in the case of planar elasticity with scalar coefficients on a convex domain, see .
We are now ready to state our main result.
If the meshes are quasi-uniform and , then we have
for all . Moreover, if the elliptic regularity property (2.4) holds, then we have
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.
for any in , and and are two positive constants. An isotropic material is nearly incompressible if is close to zero.
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
for all . Here, the constant is independent of .
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 . 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:
for all .
On each element , we also denote to be the unique solution of the local problem:
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
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
for all , where .
3.3. A characterization of the approximate solution
The above results suggest the following characterization of the numerical solution of the HDG method.
The numerical solution of the HDG method (2.2) satisfies
If we assume the boundary data , then is the solution of
In addition, the bilinear operator is positive definite.
In order to show (3.6) is true, we only need to show that for all , then
According to (3.1), we have
for all , . Then, we have
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
We apply integration by parts on (3.1a), we have
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
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. ∎
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
In the analysis, we are going to use the following classical results:
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 ) the kernel of the operator is:
Here, denotes the set of all anti-symmetric matrices in .
In the analysis, we need the following elementwise Korn’s inequality:
Let be a generic element with size and . Then for any function , we have
Here is independent of the size . In addition, if is a tetrahedron, the above inequality holds for any .
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
So, we have
The last equality above is due to the fact that every component of is constant. It is easy to see that
So, we have
By taking the symmetric part of both sides of the above equation, we have
According to Theorem in , the following inequality holds:
So, there is with and , such that
It is easy to see that
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 ,
Here, . Consequently, we have