Mixed finite element methods for linear elasticity with weakly imposed symmetry

By Douglas N. Arnold, Richard S. Falk, and Ragnar Winther

Abstract

In this paper, we construct new finite element methods for the approximation of the equations of linear elasticity in three space dimensions that produce direct approximations to both stresses and displacements. The methods are based on a modified form of the Hellinger–Reissner variational principle that only weakly imposes the symmetry condition on the stresses. Although this approach has been previously used by a number of authors, a key new ingredient here is a constructive derivation of the elasticity complex starting from the de Rham complex. By mimicking this construction in the discrete case, we derive new mixed finite elements for elasticity in a systematic manner from known discretizations of the de Rham complex. These elements appear to be simpler than the ones previously derived. For example, we construct stable discretizations which use only piecewise linear elements to approximate the stress field and piecewise constant functions to approximate the displacement field.

1. Introduction

The equations of linear elasticity can be written as a system of equations of the form

Here the unknowns and denote the stress and displacement fields engendered by a body force acting on a linearly elastic body which occupies a region . Then takes values in the space of symmetric matrices and takes values in . The differential operator is the symmetric part of the gradient, the div operator is applied row-wise to a matrix, and the compliance tensor is a bounded and symmetric, uniformly positive definite operator reflecting the properties of the body. If the body is clamped on the boundary of , then the proper boundary condition for the system Equation 1.1 is on . For simplicity, this boundary condition will be assumed here. The issues that arise when other boundary conditions are assumed (e.g., the case of pure traction boundary conditions ) are discussed in Reference 9.

The pair can alternatively be characterized as the unique critical point of the Hellinger–Reissner functional

The critical point is sought among all , the space of square-integrable symmetric matrix fields with square-integrable divergence, and all , the space of square-integrable vector fields. Equivalently, is the unique solution to the following weak formulation of the system Equation 1.1:

A mixed finite element method determines an approximate stress field and an approximate displacement field as the critical point of over where and are suitable piecewise polynomial subspaces. Equivalently, the pair is determined by the weak formulation Equation 1.3, with the test space restricted to . As is well known, the subspaces and cannot be chosen arbitrarily. To ensure that a unique critical point exists and that it provides a good approximation of the true solution, they must satisfy the stability conditions from Brezzi’s theory of mixed methods Reference 12Reference 13.

Despite four decades of effort, no stable simple mixed finite element spaces for elasticity have been constructed. For the corresponding problem in two space dimensions, stable finite elements were presented in Reference 10. For the lowest order element, the space is composed of piecewise cubic functions, with 24 degrees of freedom per triangle, while the space consists of piecewise linear functions. Another approach which has been discussed in the two-dimensional case is the use of composite elements, in which consists of piecewise polynomials with respect to one triangulation of the domain, while consists of piecewise polynomials with respect to a different, more refined, triangulation Reference 5Reference 21Reference 23Reference 31. In three dimensions, a partial analogue of the element in Reference 10 has been proposed and shown to be stable in Reference 1. This element uses piecewise quartic stresses with 162 degrees of freedom per tetrahedron, and piecewise linear displacements.

Because of the lack of suitable mixed elasticity elements, several authors have resorted to the use of Lagrangian functionals which are modifications of the Hellinger–Reissner functional given above Reference 2Reference 4Reference 6Reference 27Reference 28Reference 29Reference 30, in which the symmetry of the stress tensor is enforced only weakly or abandoned altogether. In order to discuss these methods, we consider the compliance tensor as a symmetric and positive definite operator mapping into , where is the space of matrices. In the isotropic case, for example, the mapping has the form

where are positive scalar coefficients, the Lamé coefficients. A modification of the variational principle discussed above is obtained if we consider the extended Hellinger–Reissner functional

over the space , where denotes the space of skew symmetric matrices. We note that the symmetry condition for the space of matrix fields is now enforced through the introduction of the Lagrange multiplier, . A critical point of the functional is characterized as the unique solution of the system

Clearly, if is a solution of this system, then is symmetric, i.e., , and therefore the pair solves the corresponding system Equation 1.3. On the other hand, if solves Equation 1.3, then and, if we set to the skew-symmetric part of , then solves Equation 1.5. In this respect, the two systems Equation 1.3 and Equation 1.5 are equivalent. However, the extended system Equation 1.5 leads to new possibilities for discretization. Assume that we choose finite element spaces and consider a discrete system corresponding to Equation 1.5. If is a discrete solution, then will not necessarily inherit the symmetry property of . Instead, will satisfy the weak symmetry condition

Therefore, these solutions in general will not correspond to solutions of the discrete system obtained from Equation 1.3.

Discretizations based on the system Equation 1.5 will be referred to as mixed finite element methods with weakly imposed symmetry. Such discretizations were already introduced by Fraejis de Veubeke in Reference 21 and further developed in Reference 4. In particular, the so-called PEERS element proposed in Reference 4 for the corresponding problem in two space dimensions used a combination of piecewise linear functions and cubic bubble functions, with respect to a triangulation of the domain, to approximate the stress , piecewise constants to approximate the displacements, and continuous piecewise linear functions to approximate the Lagrange multiplier . Prior to the PEERS paper, Amara and Thomas Reference 2 developed methods with weakly imposed symmetry using a dual hybrid approach. The lowest order method they discussed approximates the stresses with quadratic polynomials plus bubble functions and the multiplier by discontinuous constant or linear polynomials. The displacements are approximated on boundary edges by linear functions. Generalizations of the idea of weakly imposed symmetry to other triangular elements, rectangular elements, and three space dimensions were developed in Reference 28, Reference 29, Reference 30 and Reference 24. In Reference 29, a family of elements is developed in both two and three dimensions. The lowest order element in the family uses quadratics plus the curls of quartic bubble functions in two dimensions or quintic bubble functions in three dimensions to approximate the stresses, discontinuous linears to approximate the displacements, and discontinuous quadratics to approximate the multiplier. In addition, a lower order method is introduced that approximates the stress by piecewise linear functions augmented by the curls of cubic bubble functions plus a cubic bubble times the gradient of local rigid motions. The multiplier is approximated by discontinuous piecewise linear functions and the displacement by local rigid motions. Morley Reference 24 extends PEERS to a family of triangular elements, to rectangular elements, and to three dimensions. In addition, the multiplier is approximated by nonconforming rather than continuous piecewise polynomials.

There is a close connection between mixed finite elements for linear elasticity and discretization of an associated differential complex, the elasticity complex, which will be introduced in §3 below. In fact, the importance of this complex was already recognized in Reference 10, where mixed methods for elasticity in two space dimensions were discussed. The new ingredient here is that we utilize a constructive derivation of the elasticity complex starting from the de Rham complex. This construction is described in Eastwood Reference 18 and is based on the the Bernstein–Gelfand–Gelfand resolution; cf. Reference 11 and also Reference 14. By mimicking the construction in the discrete case, we are able to derive new mixed finite elements for elasticity in a systematic manner from known discretizations of the de Rham complex. As a result, we can construct new elements in both two and three space dimensions which are significantly simpler than those derived previously. For example, we will construct stable discretizations of the system Equation 1.5 which only use piecewise linear and piecewise constant functions, as illustrated in the figure below. For simplicity, the entire discussion of the present paper will be given in the three-dimensional case. A detailed discussion in two space dimensions can be found in Reference 8. Besides the methods discussed here, we note that by slightly generalizing the approach of this paper, one can also analyze some of the previously known methods mentioned above that are also based on the weak symmetry formulation (see Reference 19 for details).

Figure 1.

Elements for the stress, displacement, and multiplier in the lowest order case in two dimensions and three dimensions.

Graphic without alt text

An alternative approach to construct finite element methods for linear elasticity is to consider a pure displacement formulation. Since the coefficient in Equation 1.1 is invertible, the stress can be eliminated using the first equation in Equation 1.1, the stress-strain relation. This leads to the second order equation

for the displacement . A weak solution of this equation can be characterized as the global minimizer of the energy functional

over the Sobolev space . Here denotes the space of all square integrable vector fields on , with square integrable derivatives, and which vanish on the boundary . A finite element approach based on this formulation, where we seek a minimum over a finite element subspace of is standard and discussed in textbooks, (e.g., Reference 16). However, for more general models, arising, for example, in viscoelasticity and plasticity (cf. Reference 15), the stress–strain relation is not local and an elimination of the stress is impossible. For such models, a pure displacement model is excluded, and a mixed approach seems to be an obvious alternative. The construction of stable mixed elements for linear elasticity is an important step in the construction of mixed methods for these more complicated models. Another advantage of the mixed approach is that we automatically obtain schemes which are uniformly stable in the incompressible limit, i.e., as the Lamé parameter tends to infinity. Since this behavior of mixed methods is well known, we will not focus further on this property here. A more detailed discussion in this direction can, for example, be found in Reference 5.

An outline of the paper is as follows. In §2, we describe the notation to be used, state our main result, and provide some preliminary discussion on the relation between stability of mixed finite element methods and discrete exact complexes. In §3, we present two complexes related to the two mixed formulations of elasticity given by Equation 1.3 and Equation 1.5. In §4, we introduce the framework of differential forms and show how the elasticity complex can be derived from the de Rham complex. In §5, we derive discrete analogues of the elasticity complex beginning from discrete analogues of the de Rham complex and identify the required properties of the discrete spaces necessary for this construction. This procedure is our basic design principle. In §6, we apply the construction of the preceding section to specific discrete analogues of the de Rham complex to obtain a family of discrete elasticity complexes. In §7 we use this family to construct stable finite element schemes for the approximation of the mixed formulation of the equations of elasticity with weakly imposed symmetry. Finally, in §8, we show how a slightly more complicated procedure leads to a simplified elasticity element.

2. Notation, statement of main results, and preliminaries

We begin with some basic notation and hypotheses. We continue to denote by the space of -vectors, by the space of real matrices, and by and the subspaces of symmetric and skew symmetric matrices, respectively. The operators and denote the symmetric and skew symmetric parts, respectively. Note that an element of the space can be identified with its axial vector in given by the map :

i.e., for any vectors and .

We assume that is a domain in with boundary . We shall use the standard function spaces, like the Lebesgue space and the Sobolev space . For vector-valued functions, we include the range space in the notation following a semicolon, so denotes the space of square integrable functions mapping into a normed vector space . The space denotes the subspace of (vector-valued) functions in whose divergence belongs to . Similarly, denotes the subspace of (matrix-valued) functions in whose divergence (by rows) belongs to .

Assuming that is an inner product space, then has a natural norm and inner product, which will be denoted by and , respectively. For a Sobolev space , we denote the norm by and for , the norm is denoted by . The space denotes the space of polynomial functions on of total degree . Usually we abbreviate this to just .

In this paper we shall consider mixed finite element approximations derived from Equation 1.5. These schemes take the form:

Find such that

where now , and .

Following the general theory of mixed finite element methods (cf. Reference 12Reference 13) the stability of the saddle–point system Equation 2.1 is ensured by the following conditions:

(A)

whenever satisfies , and ,

(A)

for all nonzero , there exists nonzero with ,

where and are positive constants independent of .

The main result of this paper, given in Theorem 7.1, is to construct a new family of stable finite element spaces , , that satisfy the stability conditions (A) and (A). We shall show that for , the choices of the Nédélec second family of H(div) elements of degree for (cf. Reference 26) and of discontinuous piecewise polynomials of degree for and provide a stable finite element approximation. In contrast to the previous work described in the introduction, no stabilizing bubble functions are needed; nor is interelement continuity imposed on the multiplier. In §8 we also discuss a somewhat simpler lowest order element () in which the local stress space is a strict subspace of the full space of linear matrix fields.

Our approach to the construction of stable mixed elements for elasticity is motivated by the success in developing stable mixed elements for steady heat conduction (i.e., the Poisson problem) based on discretizations of the de Rham complex. We recall (see, e.g., Reference 7) that there is a close connection between the construction of such elements and discretizations of the de Rham complex

More specifically, a key to the construction and analysis of stable mixed elements is a commuting diagram of the form

Here, the spaces and are the finite element spaces used to discretize the flux and temperature fields, respectively. The spaces and are additional finite element spaces, which can be found for all well-known stable element choices. The bottom row of the diagram is a discrete de Rham complex, which is exact when the de Rham complex is (i.e., when the domain is contractible). The vertical operators are projections determined by the natural degrees of freedom of the finite element spaces. As pointed out in Reference 7, there are many such discretizations of the de Rham complex.

A diagram analogous to Equation 2.3, but with the de Rham complex replaced by the elasticity complex defined just below, will be crucial to our construction of stable mixed elements for elasticity. Discretization of the elasticity complex also gives insight into the difficulties of constructing finite element approximations of the mixed formulation of elasticity with strongly imposed symmetry; cf. Reference 8.

3. The elasticity complex

We now proceed to a description of two elasticity complexes, corresponding to strongly or weakly imposed symmetry of the stress tensor. In the case of strongly imposed symmetry, relevant to the mixed elasticity system Equation 1.3, the characterization of the divergence-free symmetric matrix fields will be needed. In order to give such a characterization, define to be the differential operator defined by taking of each row of the matrix. Then define a second order differential operator by

It is easy to check that and that . In other words,

is a complex. Here the dependence of the domain is suppressed, i.e., , and denotes the six-dimensional space of infinitesimal rigid motions on , i.e., functions of the form with and . In fact, when is contractible, then Equation 3.2 is an exact sequence, a fact which will follow from the discussion below. The complex Equation 3.2 will be referred to as the elasticity complex.

A natural approach to the construction of stable mixed finite elements for elasticity would be to extend the complex Equation 3.2 to a complete commuting diagram of the form Equation 2.3, where Equation 3.2 is the top row and the bottom row is a discrete analogue. However, due to the pointwise symmetry requirement on the discrete stresses, this construction requires piecewise polynomials of high order. For the corresponding problem in two space dimensions, such a complex was proposed in Reference 10 with a piecewise cubic stress space; cf. also Reference 8. An analogous complex was derived in the three-dimensional case in Reference 3. It uses a piecewise quartic space, with 162 degrees of freedom on each tetrahedron for the stresses.

We consider the formulation based on weakly imposed symmetry of the stress tensor, i.e., the mixed system Equation 1.5. Then the relevant complex is, instead of Equation 3.2,

Here,

and denotes the extension of the operator defined on by Equation 3.1 such that on . We remark that may be written

where is the algebraic operator

with the identity matrix. Indeed, if is symmetric, then is trace free, and therefore the definition Equation 3.4 reduces to Equation 3.1 on . On the other hand, if is skew with axial vector , then , and so .

Observe that there is a close connection between Equation 3.2 and Equation 3.3. In fact, Equation 3.2 can be derived from Equation 3.3 by performing a projection step. To see this, consider the diagram

where the projection operators are defined by