`Mathematical Problems in EngineeringVolume 2010, Article ID 759547, 20 pageshttp://dx.doi.org/10.1155/2010/759547`
Research Article

## On the Local Discontinuous Galerkin Method for Linear Elasticity

1Department of Mathematics, Shanghai Jiao Tong University, Shanghai 200240, China
2Division of Computational Science, E-Institute of Shanghai Universities and Scientific Computing Key Laboratory of Shanghai Universities, Shanghai Normal University, China

Received 25 February 2010; Accepted 22 May 2010

Copyright © 2010 Yuncheng Chen et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Following Castillo et al. (2000) and Cockburn (2003), a general framework of constructing discontinuous Galerkin (DG) methods is developed for solving the linear elasticity problem. The numerical traces are determined in view of a discrete stability identity, leading to a class of stable DG methods. A particular method, called the LDG method for linear elasticity, is studied in depth, which can be viewed as an extension of the LDG method discussed by Castillo et al. (2000) and Cockburn (2003). The error bounds in -norm, -norm, and a certain broken energy norm are obtained. Some numerical results are provided to confirm the convergence theory established.

#### 1. Introduction

This paper is focused on systematically studying discontinuous Galerkin (DG) methods for the linear elasticity problem. Since the DG method was first introduced in 1970s, these methods have been applied for solving a variety of mathematical-physical problems including linear and nonlinear hyperbolic problems, Navier-Stokes equations, convection-dominated diffusion problems, and so on. The DG method may be viewed as high-order extensions of the classical finite volume method. Since no inter-element continuity is imposed, such methods can be defined on very general meshes including nonconforming meshes. Moreover, polynomials of arbitrary degree can be used on each element, making these methods suitable for -adaptivity. We refer to [1] for an excellent historical survey. As a generalization of the DG method proposed in [2] for the solution of the compressible Navier-Stokes equations, some local discontinuous Galerkin (LDG) methods were introduced and analyzed in [3, 4] for time-dependent convection-diffusion systems and second-order elliptic problems, respectively. After that, a unified theory was developed for conducting error analysis for DG methods of elliptic problems (cf. [5]). And a new framework was proposed in [6] for designing and analyzing DG methods, where stabilization mechanism frequently used in DG methods was also investigated. In [7], an important framework was proposed for constructing stable DG methods.

On the other hand, linear elasticity discusses how solid objects deform and become internally stressed on prescribed loading conditions, and is encountered extensively in structural analysis and engineering design. By use of usual displacement-based finite-element methods for the elasticity equation, we are able to numerically determine the displacement field directly. However, in many engineering applications, the stress field is a quantity of more interest. For this, we may apply mixed finite-element methods to solve the linear elasticity system, from which the aforementioned two physical quantities can be obtained simultaneously. Over the past four decades, there have been many efforts along this line. But due to the symmetry constraint on the stress tensor, it is extremely difficult to construct stable stress-displacement finite elements. In two spatial dimensions, the first stable finite element with polynomial shape functions is presented in [8]. For the lowest-order element, the discrete stress space is composed of certain piecewise cubic polynomials, while the discrete displacement space consists of piecewise linear polynomials. In three dimensions, a piecewise quadratic stress space is constructed with degrees of freedom on each tetrahedron (cf. [9]). Another approach in this direction is to use composite elements (macroelements), in which the discrete displacement space consists of piecewise polynomials with respect to one triangulation of the domain, while the discrete stress space consists of piecewise polynomials with respect to a more refined triangulation (cf. [1013]). It is mentioned that for solving the previous problem, several mixed elements with weakly imposed symmetry have also been developed (cf. [1416]).

Regarding the complexity of mixed elements given above, the discontinuous Galerkin method is naturally a suitable alternative for numerically solving linear elasticity problems. To the best of our knowledge, a local DG method and an interior penalty DG method for linear elasticity are presented in [17] and [18, 19], respectively. A mixed DG method is given in [20], which one may find is covered by our general formulation below. Moreover, a mixed formulation is also extended to the case of nonsymmetric stress tensors (cf. [21]).

In this paper, we are going to look for new DG methods for the linear elasticity problem. Following [4, 7], we build up a framework to construct our DG methods. Then a discrete stability identity is established, from which we derive feasible choices of numerical traces and get a class of stable DG methods for linear elasticity. With a parameter taken to be zero, namely, , we get an LDG method. Following [5] and using some technical arguments, we get optimal-order error estimates for the LDG method proposed in a certain broken energy norm, -norm, and -norm, respectively. It should be emphasized that in our formulation the symmetry of the stress tensor is preserved automatically.

The rest of this paper is organized as follows. The basic framework of DG methods and the determination of numerical traces are presented in Section 2. An error analysis for the LDG method is given in Section 3. And in Section 4, some numerical results are included to confirm our theoretical convergence orders.

#### 2. The DG Method for Linear Elasticity

##### 2.1. Basic Framework for the DG Method

Assume that is a bounded polygon or polyhedron. Let be the stress, let be the displacement, and let be the body force. Denote by the linearized strain tensor with , tr the trace operator, and the divergence operator. Then, the mathematical model of linear elasticity reads (cf. [22]) where is the fourth-order compliance tensor defined by where is the Kronecker tensor, and the positive constants and are the Lamé coefficients of the elastic material under consideration.

For defining our DG methods for problem (2.1), we introduce some notations first of all. For any Banach space , the subspace of symmetric matrix-valued function is denoted by . Given a bounded domain and a nonnegative integer , let be the usual Sobolev space consisting of all functions in whose weak derivatives up to degree are also -integrable (cf. [23]). The corresponding norm is denoted by . Let be the closure of with respect to the norm . Moreover, we simply write for .

Let be a regular family of triangulations of (cf. [23, 24]); and . Let be the union of all faces (edges) of the triangulation and the union of all interior faces (edges) of the triangulation . For any , is set to be the diameter of . Based on the triangulation , let The corresponding finite element spaces are given by where, for each , and are two finite-dimensional spaces of polynomials in containing and , respectively, with . Here, for a nonnegative integer , stands for the set of all polynomials in with the total degree no more than .

To guarantee uniqueness of the solution to the LDG method to be proposed, we always assume that For a function with for all , let be the usual broken -type norm of defined by If is a vector-valued or matrix-valued function, the corresponding term is defined in a similar manner. For a vector or a matrix , denote by the quantity or . Here the symbol stands for the double dot product operation of matrices. Throughout this paper, we use the notation “” to indicate “”, where is a generic positive constant independent of and other parameters, which may take different values at different appearances.

Let and be two adjacent elements of . Let be an arbitrary point of the set , and let and be the corresponding outward unit normals at that point. For a vector-valued function smooth inside each element , let us denote by the trace of on from the interior of . Then we define averages and jumps at as follows: If is on an edge/face lying on the boundary , the above terms are defined by where is the unit outward normal vector on . We define a matrix-valued jump of a vector as follows: where denotes the matrix whose th component is for two vectors and .

Now, we are ready to introduce a framework to derive DG methods for problem (2.1). Following [4, 7], we first derive a variational formulation for problem (2.1). Taking a double dot product with a matrix-valued function on both sides of the first equation of (2.1) and integrating by parts over , we have Multiplying the second equation of (2.1) by a vector-valued function and integrating by parts over yields

Motivated by the above two identities, we may define our DG method as follows. Find an approximate solution such that for all and all . Note that any function with the hat superscript is only defined over all faces of the triangulation , which is called a numerical trace. Since is symmetric, it is natural to choose as a symmetric matrix-valued function. Moreover, we only consider the case where the numerical traces are single valued over all faces (conservation).

##### 2.2. Numerical Traces and the LDG Method

We begin by producing a stability identity for the continuous problem (2.1), a crucial step in constructing feasible numerical traces to get a stable DG method from (2.12)-(2.13). For this, taking a double dot product with on both sides of the first equation of (2.1) and then integrating over , we have Multiplying the second equation of (2.1) by and then integrating over again, we find from the homogeneous boundary condition of that Now adding these two equations gives This is the stability identity corresponding to the continuous problem (2.1).

Next, we mimic the above derivation to get a discrete analogue of the stability identity (2.16) for the DG method (2.12)-(2.13). Taking in (2.12) and summing over all , we have Similarly, summing up (2.13) over all with , we come to Adding the last two equations gives where

Lemma 2.1. Assume and are both symmetric, and the numerical traces and are single valued over all . Then

Proof. We only prove the first identity, and the other two ones can be obtained in the similar manners. The left-hand side of (2.21) can be rewritten as On the other hand, Hence, the identity (2.21) is a direct consequence of (2.24)-(2.25).

By Lemma 2.1 and using integration by parts twice, we see that Thus, if , we take and if , we take where and are two nonnegative continuous functions on . When , the corresponding method is called the LDG method for linear elasticity, viewed as a generalization of the LDG method for second-order elliptic problems in [4, 7].

For the above choice of numerical traces, we have by some direct manipulation that Then, the combination of (2.19) and (2.29) allows a discrete stability estimate: which is essential in constructing a reliable DG method (cf. [7]).

Let us further show the unique solvability of problem (2.12)-(2.13) with the numerical traces given by (2.27) and (2.28), whenever and the finite element spaces and satisfy condition (2.5). In fact, it suffices to verify that this DG method only has zero solution when . By setting , the stability identity (2.19) gives which implies that in , on , and on , provided that . Therefore, from (2.12) the definition of (cf. (2.27) and (2.28)), and Lemma 2.1, we know that which, due to (2.5), implies that . By Korn's inequality (2.34) on the discontinuous finite element space, given below, it is easy to see that . Then with on and , we conclude that in , as required.

Remark 2.2. The standard Korn's inequality states that (cf. [25]) The following Korn's inequality on the discontinuous finite-element space is given in [26]:

As in [4, 7], the DG method (2.12)-(2.13) with the numerical traces (2.27)-(2.28) can also be written in a mixed formulation (cf. [27]). After some direct manipulation, the approximate solution can be characterized as the unique solution of the following variational problem. Find such that for all , where

Remark 2.3. If we take and , the mixed formulation (2.35)-(2.36) will be reduced to the one adopted in [20].

#### 3. Error Analysis for the LDG Method

In this section, we provide an error analysis for the DG method (2.35)-(2.36) in the case , which is named as the LDG method in the last section. For this purpose, we are going to derive a primal formulation for the approximate method. First we introduce a global lifting operator defined by Moreover, for each , we introduce a local lifting operator defined by Then it is easy to check that is only nonzero in the triangles with as one edge, and there holds the identity Now we express in terms of . From the first equation in Lemma 2.1 and (2.35), it follows that for all . Then by (2.5), we get that is, Substituting from the last equation into (2.36), we get a primal formulation for the LDG method as follows. Find such that where

Next, we consider the consistency of the method (3.7). Assume that , where is some nonnegative integer, and is the exact solution of (2.1). It is easy to see that on . Therefore, for any , we have From (3.9), the definition of lifting operator (see (3.1)), and the fact that , we further find that where Noting that in and using integration by parts and the same technique for deriving (2.21), we deduce from (2.1) that Hence, we know that the LDG method (3.7) is not consistent with respect to the bilinear form , but it admits the following identity:

It is worth noting that the LDG method does not contain any parameter which needs to be quantified a priori. In what follows, we choose on each with having a uniform positive bound from above and below.

Next, we present a useful estimate for the local lifting operator . For this, let and define the mesh-dependent energy norm for by

Lemma 3.1. For any and , one has

Proof. Since implies that , it suffices to verify the result for . Taking and in (3.2), we get On the other hand, by the scaling arguments, the trace theorem, and the inverse inequality, Therefore, (3.15) is a direct consequence of (3.16) and (3.17).

Lemma 3.2 (boundedness). For any , it holds that

Proof. According to the Cauchy-Schwarz inequality and the basic properties of and , we have which, together with (3.15), yields as required.

Lemma 3.3 (stability). For any , it holds that

Proof. Using the Cauchy-Schwarz inequality and (3.15), we get Here , , and is the constant in (3.15). Therefore (3.21) is true if we choose such that , that is, .

Now, we are in a position to give error analysis for the LDG method (3.7). The main idea of our derivation is based on the framework on error analysis of DG methods for second-order elliptic problems (cf. [5]). Let be an -projection operator from onto the finite-element space . Let denote the usual -orthogonal projection operator onto . For simplicity, we still write (resp., ) for (resp., ). Using the standard scaling arguments (cf. [24]), we can easily obtain error estimates for the operators and .

Lemma 3.4. Let , with . Then for all , one has

Theorem 3.5. Let , where is a nonnegative integer and is the solution of (2.1). Let be the solution of (3.7). Then one has

Proof. From the stability (3.21), the identity (3.13), the boundedness (3.18), and Lemmas 3.1 and 3.4, it follows that that is, Therefore, using the triangle inequality and Lemma 3.4, we further have

Theorem 3.6. Let , where is a nonnegative integer and is the solution of (2.1). Let be the solution of (3.7). Then one has

Proof. By the triangle inequality and Korn's inequality over (cf. (2.34)), we know that The desired result then follows from the above estimate combined with Lemma 3.4 and Theorem 3.5.

Theorem 3.7. Suppose that is a convex bounded polygon or polyhedron. Let , where is a nonnegative integer and is the solution of (2.1). Let be the solution of (3.7). Then one has

Proof. The proof relies on the usual duality argument. Let be the solution of the auxiliary problem Formally, (3.31) is problem (2.1) with replaced by . Since is a convex bounded polygon or polyhedron, we have an -regularity estimate (cf. [28])
Therefore, applying an elementwise integration by parts to the second equation of (3.31) with a test function and using the definitions of (cf. (3.1)) and , the first equation of (3.31), and the technique to derive (2.21), we get and with the identity (3.13), we further have Therefore, it follows from the boundedness (3.18), the regularity (3.31), and Lemmas 3.1 and 3.4 that This, along with Theorem 3.5, immediately leads to

We would like to end this section with some results on the optimality of our estimates derived. If we choose and to be and , respectively, it is easy to check that the first condition of (2.5) implies . Therefore, we can obtain the following result from the previous theorems directly.

Corollary 3.8. Let , where is a nonnegative integer and is the solution of (2.1). Let be the solution of (3.7) with and for all . Then one has the following optimal error estimates: In addition, if is a convex bounded polygon or polyhedron,

#### 4. Numerical Results

In this section, we report a numerical example to illustrate the theoretical results. Let , , , and It is not difficult to check that the exact solution of (2.1) is We use a quasiuniform triangulation over . For any , we take and where . On each edge , is set to be . The numerical results of the LDG method for a few choices of and are shown in Table 1 and Figures 1, 2, and 3. We observe that, when , the numerical convergence rates of , , and are , , and , respectively. However, there is no convergence when . These phenomena agree with theoretical results in Theorems 3.53.7. When , accuracies of the numerical results are nearly the same as those for and . It goes without saying that it is more convenient to simulate for than for .

Table 1: Numerical errors of LDG method for several choices of and . The theoretical convergence rates of , , and are , , and , respectively.
Figure 1: in the - scale for several choices of and .
Figure 2: in the - scale for several choices of and .
Figure 3: in the - scale for several choices of and .

#### Acknowledgments

The authors thank the referee for the valuable comments and suggestions, which greatly improved the presentation of the manuscript. The research of the second author is supported partially by the NNSFC (10371076) and E-Institutes of Shanghai Municipal Education Commission (E03004).

#### References

1. B. Cockburn, G. E. Karniadakis, and C.-W. Shu, Eds., Discontinuous Galerkin Methods. Theory, Computation and Applications, vol. 11 of Lecture Notes in Computational Science and Engineering, Springer, Berlin, Germany, 2000.
2. F. Bassi and S. Rebay, “A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations,” Journal of Computational Physics, vol. 131, no. 2, pp. 267–279, 1997.
3. B. Cockburn and C.-W. Shu, “The local discontinuous Galerkin method for time-dependent convection-diffusion systems,” SIAM Journal on Numerical Analysis, vol. 35, no. 6, pp. 2440–2463, 1998.
4. P. Castillo, B. Cockburn, I. Perugia, and D. Schötzau, “An a priori error analysis of the local discontinuous Galerkin method for elliptic problems,” SIAM Journal on Numerical Analysis, vol. 38, no. 5, pp. 1676–1706, 2000.
5. D. N. Arnold, F. Brezzi, B. Cockburn, and L. D. Marini, “Unified analysis of discontinuous Galerkin methods for elliptic problems,” SIAM Journal on Numerical Analysis, vol. 39, no. 5, pp. 1749–1779, 2001/02.
6. F. Brezzi, B. Cockburn, L. D. Marini, and E. Süli, “Stabilization mechanisms in discontinuous Galerkin finite element methods,” Computer Methods in Applied Mechanics and Engineering, vol. 195, no. 25-28, pp. 3293–3310, 2006.
7. B. Cockburn, “Discontinuous Galerkin methods,” Zeitschrift für Angewandte Mathematik und Mechanik, vol. 83, no. 11, pp. 731–754, 2003.
8. D. N. Arnold and R. Winther, “Mixed finite elements for elasticity,” Numerische Mathematik, vol. 92, no. 3, pp. 401–419, 2002.
9. S. Adams and B. Cockburn, “A mixed finite element method for elasticity in three dimensions,” Journal of Scientific Computing, vol. 25, no. 3, pp. 515–521, 2005.
10. B. M. Fraeijs de Veubeke, “Displacement and equilibrium models the finite element method,” in Stress Analysis, O. C. Zienkiewicz and G. S. Holister, Eds., pp. 145–197, Wiley, New York, NY, USA, 1965.
11. V. B. Watwood Jr. and B. J. Hartz, “An equilibrium stress field model for finite element solution of two-dimensional elastostatic problems,” International Journal of Solids and Structures, vol. 4, pp. 857–873, 1968.
12. C. Johnson and B. Mercier, “Some equilibrium finite element methods for two-dimensional elasticity problems,” Numerische Mathematik, vol. 30, no. 1, pp. 103–116, 1978.
13. D. N. Arnold, J. Douglas Jr., and C. P. Gupta, “A family of higher order mixed finite element methods for plane elasticity,” Numerische Mathematik, vol. 45, no. 1, pp. 1–22, 1984.
14. D. N. Arnold, F. Brezzi, and J. Douglas, Jr., “PEERS: a new mixed finite element for plane elasticity,” Japan Journal of Applied Mathematics, vol. 1, no. 2, pp. 347–367, 1984.
15. D. N. Arnold, R. S. Falk, and R. Winther, “Mixed finite element methods for linear elasticity with weakly imposed symmetry,” Mathematics of Computation, vol. 76, no. 260, pp. 1699–1723, 2007.
16. W. Qiu and L. Demkowicz, “Mixed $hb$-finite element method for linear elasticity with weakly imposed symmetry,” Computer Methods in Applied Mechanics and Engineering, vol. 198, no. 47-48, pp. 3682–3701, 2009.
17. A. Lew, P. Neff, D. Sulsky, and M. Ortiz, “Optimal BV estimates for a discontinuous Galerkin method for linear elasticity,” Applied Mathematics Research eXpress, vol. 2004, no. 3, pp. 73–106, 2004.
18. P. Hansbo and M. G. Larson, “Discontinuous Galerkin methods for incompressible and nearly incompressible elasticity by Nitsche's method,” Computer Methods in Applied Mechanics and Engineering, vol. 191, no. 17-18, pp. 1895–1908, 2002.
19. P. Hansbo and M. G. Larson, “Discontinuous Galerkin and the Crouzeix-Raviart element: application to elasticity,” M2AN. Mathematical Modelling and Numerical Analysis, vol. 37, no. 1, pp. 63–72, 2003.
20. Z. Cai and X. Ye, “A mixed nonconforming finite element for linear elasticity,” Numerical Methods for Partial Differential Equations, vol. 21, no. 6, pp. 1043–1051, 2005.
21. R. Bustinza, “A note on the local discontinuous Galerkin method for linear problems in elasticity,” Scientia Series A, vol. 13, pp. 72–83, 2006.
22. K. Feng and Z. Shi, Mathematical Theory of Elastic Structures, Springer, Berlin, Germany, 1995.
23. S. C. Brenner and L. R. Scott, The Mathematical Theory of Finite Element Methods, vol. 15 of Texts in Applied Mathematics, Springer, New York, NY, USA, 3rd edition, 2008.
24. P. G. Ciarlet, The Finite Element Method for Elliptic Problems, vol. 4 of Studies in Mathematics and Its Applications, North-Holland, Amsterdam, The Netherlands, 1978.
25. G. Fichera, “Existence theorems in elasticity,” in Handbuch der Physik, Band VI, pp. 347–389, Springer, Heidelberg, Germany, 1972.
26. S. C. Brenner, “Korn's inequalities for piecewise ${H}^{1}$ vector fields,” Mathematics of Computation, vol. 73, no. 247, pp. 1067–1087, 2004.
27. F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, vol. 15 of Springer Series in Computational Mathematics, Springer, New York, NY, USA, 1991.
28. P. Grisvard, Singularities in Boundary Value Problems, vol. 22 of Recherches en Mathématiques Appliquées, Masson, Paris, France; Springer, Berlin, Germany, 1992.