# A DG method for a stress formulation of the elasticity eigenproblem with strongly imposed symmetry

We introduce a pure–stress formulation of the elasticity eigenvalue problem with mixed boundary conditions. We propose an H(div)-based discontinuous Galerkin method that imposes strongly the symmetry of the stress for the discretization of the eigenproblem. Under appropriate assumptions on the mesh and the degree of polynomial approximation, we demonstrate the spectral correctness of the discrete scheme and derive optimal rates of convergence for eigenvalues and eigenfunctions. Finally, we provide numerical examples in two and three dimensions.

## Authors

• 5 publications
04/07/2022

### A new DG method for a pure–stress formulation of the Brinkman problem with strong symmetry

A strongly symmetric stress approximation is proposed for the Brinkman e...
03/03/2022

### Symmetric mixed discontinuous Galerkin methods for linear viscoelasticity

We propose and rigorously analyse semi- and fully discrete discontinuous...
03/23/2022

### Spectral analysis of a mixed method for linear elasticity

The purpose of this paper is to analyze a mixed method for linear elasti...
04/22/2020

### Approximation of the Axisymmetric Elasticity Equations with Weak Symmetry

In this article we consider the linear elasticity problem in an axisymme...
04/17/2022

11/12/2019

### Tutorial on Hybridizable Discontinuous Galerkin (HDG) Formulation for Incompressible Flow Problems

A hybridizable discontinuous Galerkin (HDG) formulation of the linearize...
11/10/2020

### Locking free staggered DG method for the Biot system of poroelasticity on general polygonal meshes

In this paper we propose and analyze a staggered discontinuous Galerkin ...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The finite element determination of the vibration characteristics (natural frequencies and mode shapes) of elastic bodies is of great interest in structural mechanics. For example, the knowledge of the eigenfrequencies keeps the forced oscillations safe from resonance regimes, and the eigenmodes can be used to expand the solution of transient elastodynamic problems in a Fourier series. We approach this topic from the perspective of the mixed formulation derived from the Hellinger-Reissner variational principle. Namely, we are interested in variational formulations in which the Cauchy stress tensor prevails as the main unknown. In addition to the fact that accurate approximations of the stress are of paramount importance in many applications, it is well known that mixed formulations are immune to locking in the case of nearly incompressible materials.

In recent years, the theory of Descloux–Nassif–Rappaz [11, 12] for non-compact operators has been successfully applied to the mixed finite element analysis of eigenvalue problems in elasticity [26, 23, 25]. The same approach allowed to deal with mixed formulations of the Stokes eigenproblem formulated in terms of a pseudo-stress [27, 22] or the Cauchy stress tensor [25]

. The symmetry requirement for the stress tensor, which reflects the conservation of angular momentum, is a specific feature of the Hellinger-Reissner variational principle. The imposition of this restriction in association with H(div)-conformity gives rise to conforming Galerkin methods with a very large number of degrees of freedom, and which are difficult to implement

[2, 20]. A common practice to overcome this drawback consists in enforcing the symmetry constraint variationally through a Lagrange multiplier. In this context, [26, 25] validated the use of the weakly symmetric mixed finite elements [5, 4, 9, 18] for the stress formulation of the elasticity eigenproblem.

Motivated by the ability of DG methods to handle efficiently -adaptive strategies and to facilitate the implementation of high order methods, an H(div)-based interior penalty version of [26] (that retains the weak imposition of the symmetry) has been introduced in [23]. Nevertheless, on account of [3, 19, 33], it is known that relaxing H(div)-conformity by using non-conforming or DG approximations for the elasticity source problem allows the incorporation of the symmetry constraint in the energy space at a reasonable computational cost. To our knowledge, the eigenvalue numerical analysis of such non-conforming/DG mixed methods is not yet available. In this work, our main issue is to determine whether a strong imposition of the symmetry constraint in the scheme introduced in [23] provides a correct eigenvalue approximation.

The resulting DG-method approximates the stress by symmetric tensors with piecewise polynomial entries of degree , in 2D and 3D. We note that, the stress/displacement DG formulation introduced in [33] for the elasticity source problem relays on the same discrete space for the stress. However, the displacement field is not present as an independent variable in our DG formulation because it is eliminated via the momentum balance equation. The same equation can be used to post-process the displacement at the discrete level. We prove that the inf-sup stability of the Scott-Vogelius element [31] for the Stokes problem (see Assumption 2 below) is a sufficient condition for the spectral correctness of our DG method. We also obtain optimal error estimates for eigenvalues and eigenfunctions in an adequate DG norm.

We finally highlight that, unlike [26, 23], our analysis does not rely on any extra Sobolev regularity of an auxiliary elasticity source problem. Hence, our analysis remains valid for eigenproblems posed in general domains, with mixed boundary conditions and with minimal requirements on material coefficients.

Outline. The contents of this paper have been organized in the following manner. The remainder of this section contains notational conventions and definitions of Sobolev spaces. Section 2 presents the pure–stress formulation of the elasticity eigenproblem and provides a characterization of its spectrum. Preliminary definitions and auxiliary tools related with H(div)-based discontinuous Galerkin methods are collected in Section 3. The definition of the mixed DG method (with strong symmetry of the stress) is detailed in Section 4, where we also introduce a couple of operators that are useful in our analysis. The spectral correctness of the DG scheme is treated in Section 5

, together with the deduction of optimal error estimates for eigenvalues and eigenspaces. Several numerical results are presented in Section

6, confirming the expected rates of convergence for different parameter sets including the nearly incompressible regime.

Notations and Sobolev spaces. We denote the space of real matrices of order by and let be the subspace of symmetric matrices, where stands for the transpose of . The component-wise inner product of two matrices is defined by . We also introduce and denote by the identity in . Along this paper we convene to apply all differential operators row-wise. Hence, given a tensorial function

and a vector field

, we set the divergence , the gradient , and the linearized strain tensor as

 (divσ)i:=∑j∂jσij,(∇u)ij:=∂jui,%andε(u):=12(∇u+(∇u)t).

Let be a polyhedral Lipschitz domain of , with boundary . For , stands for the usual Hilbertian Sobolev space of functions with domain and values in E, where is either , or . In the case we simply write . The norm of is denoted and the corresponding semi-norm , indistinctly for . We use the convention and let be the inner product in , for , namely,

 (u,v):=∫Ωu⋅v∀u,v∈L2(Ω,Rd),(σ,τ):=∫Ωσ:τ∀σ,τ∈L2(Ω,M).

We consider the space of tensors satisfying , and denote the corresponding norm , where is either or . Let be the outward unit normal vector to . Let be a sufficiently regular symmetric tensor, Green’s formula

 (τ,ε(v))+(divτ,v)=∫∂Ωτn⋅vv∈H1(Ω,Rd), (1.1)

can be used to extend the normal trace operator to a linear continuous mapping , where is the dual of .

Throughout this paper, we shall use the letter to denote a generic positive constant independent of the mesh size , that may stand for different values at its different occurrences. Moreover, given any positive expressions and depending on , the notation means that .

## 2 A stress formulation of the elasticity eigenproblem

Our aim is to determine the natural frequencies of an elastic structure with mass density and occupying a polyhedral Lipschitz domain () . This amounts to solve the eigenproblem

 divσ+ω2ϱ(x)u =0in Ω, (2.1) A(x)σ =ε(u)in Ω, (2.2)

where is the displacement field and is the Cauchy stress tensor. The symmetric and positive-definite 4-order tensor involved in the linear material law (2.2) is known as the compliance tensor. We assume that there exist such that

 a−ζ:ζ≤A(x)ζ:ζ≤a+ζ:ζ∀ζ∈S,a.e. % in Ω.

We also suppose that there exists a polygonal/polyhedral disjoint partition of such that for all and let and .

We impose the boundary condition on a subset of positive surface measure and let the structure free of stress on . Here, we opt for combining the equilibrium equation (2.1) with the constitutive law (2.2) to eliminate the displacement field and impose as a primary variable. This procedure leads to the following eigensystem: find eigenmodes and eigenfrequencies such that,

 −ε(1ϱdivσ)=ω2Aσ in Ω,1ϱdivσ=0 on ΓD,σn=0 on ΓN, (2.3)

where stands for the exterior unit normal vector on .

In the following, we write for the space endowed with the -weighted inner product and denote the corresponding norm . The eigenfunctions will be sought in the closed subspace of defined by

where holds for the duality pairing between and . We introduce the symmetric and positive semidefinite bilinear form given by

 c(σ,τ):=(1ϱdivσ,divτ)

and endow with the Hilbertian inner product . We denote the corresponding norm .

Testing the first equation of (2.3) with and applying Green’s formula (1.1) we deduce, after a shift argument, the following pure–stress variational formulation of the eigenproblem: find and such that

 a(σ,τ)=κ(σ,τ)A,∀τ∈X. (2.4)

We introduce the source operator corresponding to the variational eigenproblem (2.4); which is defined for any by

 a(~Tf,τ)=(f,τ)A,∀τ∈X. (2.5)

Obviously, is linear and bounded, actually it holds,

 ∥∥~Tf∥∥X≤∥f∥A∀f∈H. (2.6)

We denote the -Sobolev space with incorporated Dirichlet boundary conditions on either or by

 H1⋆(Ω,Rd):={v∈H1(Ω,Rd); v|Γ⋆=0},⋆∈{D,N}.

It is important to notice that testing (2.5) with a tensor whose entries are indefinitely differentiable and compactly supported in proves that . Hence, by virtue of Korn’s inequality, and it follows readily from Green’s formula (1.1) that vanishes on . In other words, and there exists such that

 ∥∥1ϱdiv(~Tf)∥∥1,Ω≤C∥f∥A,∀f∈H. (2.7)

The operator is relevant in our analysis because its eigenvalues and those of problem (2.4) are reciprocal to each other and the corresponding eigenfunctions are the same. A full description of the spectrum of will then solve problem (2.4).

We consider the direct sum decomposition into closed subspaces

 K:={τ∈X; divτ=0 in Ω}andK⊥:={σ∈X; (σ,τ)A=0, ∀τ∈K},

which are orthogonal with respect to both and . It is clear that is an eigenvalue of (2.4) with associated eigenspace . Consequently, as is not a finite-dimensional subspace of , is not a compact operator.

###### Lemma 2.1.

The orthogonal projection in onto is characterized, for any , by where and is the unique solution of

 (A−1ε(˜u),ε(v))=−(divσ,v),∀v∈H1D(Ω,Rd). (2.8)
###### Proof.

We first point out that Korn’s inequality provides the stability estimate

 ∥~u∥1,Ω≤C∥divσ∥0,Ω. (2.9)

By definition, (2.9) also ensures that . Moreover, by construction, which ensures that is bounded. Moreover, it is clear that and . It remains to show that the range of coincides with . To this end, we notice that, for any ,

 (Pσ,τ)A=(ε(˜u),τ)=(∇˜u,τ)=0,∀τ∈K,

which proves that . The reciprocal inclusion is a consequence of , where we used that , and the result follows.

###### Lemma 2.2.

The inclusions and are compact.

###### Proof.

Let be a weakly convergent sequence in . The continuiuty of implies that the sequence is also weakly convergent in . By definition, , where solves (2.8) with right-hand side . It follows from (2.9) that is bounded in and the compactness of the embedding implies that admits a subsequence (denoted again ) that converges strongly in . Next, we deduce from Green’s identity

 (˜σp−˜σq,˜σp−˜σq)A=(ε(˜up−˜uq),˜σp−˜σq)=−(˜up−˜uq,div(˜σp−˜σq)),

that is a Cauchy sequence in , which implies that the embedding is compact.

Finally, it follows from (2.7) that

 T(X)∩P(X)⊂{σ∈P(X); 1ϱdivσ∈H1(Ω,Rd)},

and the compactness of the embedding is a consequence of the fact that the inclusion is compact. ∎

We point out that is symmetric with respect to , which implies that is -invariant. Consequently, it holds and Lemma 2.2 implies that the -symmetric and positive definite operator is compact. Therefore, we have the following characterization of the spectrum of .

###### Theorem 2.1.

The spectrum of is given by , where is a sequence of finite-multiplicity eigenvalues of that converges to 0. The ascent of each of these eigenvalues is and the corresponding eigenfunctions lie in . Moreover, is an infinite-multiplicity eigenvalue of with associated eigenspace and is not an eigenvalue.

## 3 Definitions and auxiliary results

We consider a sequence of shape-regular simplicial meshes that subdivide the domain into simplices of diameter . The parameter represents the mesh size of . We assume that is aligned with the partition and that is a shape-regular mesh of for all and all .

For all , we consider the broken Sobolev space

 Hs(∪jΩj):={v∈L2(Ω); v|Ωj∈Hs(Ωj), ∀j=1,…,J}

corresponding to the partition . Its vectorial and tensorial versions are denoted and , respectively. Likewise, the broken Sobolev space with respect to the subdivision of into is

 Hs(Th,E):={v∈L2(Ω,E):v|K∈Hs(K,E)∀K∈Th},for E∈{R,Rd,M}.

For each and the components and represent the restrictions and . When no confusion arises, the restrictions of these functions will be written without any subscript.

Hereafter, given an integer and a domain , denotes the space of polynomials of degree at most on . We introduce the space

 Pm(Th):={v∈L2(Ω): v|K∈Pm(K), ∀K∈Th}

of piecewise polynomial functions relatively to . We also consider the space of functions with values in and entries in , where is either , or .

Let us introduce now notations related to DG approximations of -type spaces. We say that a closed subset is an interior edge/face if has a positive -dimensional measure and if there are distinct elements and such that . A closed subset is a boundary edge/face if there exists such that is an edge/face of and . We consider the set of interior edges/faces, the set of boundary edges/faces and let be the set of edges/faces composing the boundary of . We assume that the boundary mesh is compatible with the partition in the sense that, if and then and . We denote

 Fh:=F0h∪F∂h% andF∗h:=F0h∪FNh,

and for all . Obviously, in the case we have that .

We will need the space given on the skeletons of the triangulations by . Its vector valued version is denoted . Here again, the components of coincide with the restrictions . We endow with the inner product

 (u,v)F∗h:=∑F∈F∗h∫FuF⋅vF∀u,v∈L2(F∗h,Rd),

and denote the corresponding norm . From now on, is the piecewise constant function defined by for all with denoting the diameter of edge/face . By virtue of our hypotheses on and on the triangulation , we may consider that is an element of and denote for all . We introduce defined by if and if .

Given and , with , we define averages and jumps by

 {v}F:=(vK+vK′)/2and⟦τ⟧F:=τKnK+τK′nK′∀F∈F(K)∩F(K′),

with the conventions

 {v}F:=vKand⟦τ⟧F:=τKnK∀F∈F(K),F∈F∂h,

where is the outward unit normal vector to .

For any , we let , with . Given , we define by for all and endow with the norm

If it happens that with , we also introduce

It is important to notice that for all .

The following discrete trace inequality is useful in our analysis.

###### Lemma 3.1.

There exists a constant independent of and such that

 (3.1)
###### Proof.

By definition of , for any , it holds

Applying in the last inequality the well-known estimate (see for example [13])

 h12K∥ϕ∥0,∂K≤Ctr∥ϕ∥0,K∀ϕ∈Pk(K), (3.2)

where is independent of , we obtain the result. ∎

For all and for a large enough given parameter , we consider the symmetric bilinear form

 ch(σ,τ):=c(σ,τ)+a(ϱ−1Fh−1F⟦σ⟧,⟦τ⟧)F∗h−({1ϱdivhσ},⟦τ⟧)F∗h−({1ϱdivhτ},⟦σ⟧)F∗h

and let

 ah(σ,τ):=(Aσ,τ)+ch(σ,τ).

For all satisfying with , a straightforward application of the Cauchy-Schwarz inequality gives

Moreover, if we take in the last estimate , we deduce from Lemma 3.1 that,

 (3.3)

with .

The bilinear form and the DG-norm are designed in such a way that the coercivity of the bilinear form on can be achieved with a stability parameter that is independent of the material coefficients, as shown in the following result.

###### Proposition 3.1.

There exists a constant , independent of and , such that if , then

 ch(τ,τ)≥12(∥∥∥ϱ−12divhτ∥∥∥20,Ω+∥∥∥δ−12Fh−12F⟦τ⟧∥∥∥20,F∗h),∀τ∈XDGh. (3.4)
###### Proof.

By definition, we have

 ch(τ,τ)=∥∥∥ϱ−12divhτ∥∥∥20,Ω+a∥∥∥ϱ−12Fh−12F⟦τ⟧∥∥∥20,F∗h−2({ϱ−1divhτ},⟦τ⟧)F∗h (3.5)

Using the Cauchy-Schwarz inequality, Young’s inequality together with the discrete trace inequality (3.1) we obtain the estimate

 2∣∣({ϱ−1divhτ},⟦τ⟧)F0h∣∣≤2∥∥∥γ12Fh12F{ϱ−1divhτ}∥∥∥0,F∗h∥∥∥γ−12Fh−12F⟦τ⟧∥∥∥0,F∗h≤2Ctr∥∥∥ϱ−12divτ∥∥∥0,Ω∥∥∥γ−12Fh−12F⟦τ⟧∥∥∥0,F∗h≤12∥∥∥ϱ−12divτ∥∥∥20,Ω+2C2tr∥∥∥γ−12Fh−12F⟦τ⟧∥∥∥20,F∗h. (3.6)

Combining (3.6) and (3.5) gives the result with . ∎

## 4 The pure–stress DG scheme

We are now in a position to introduce the following mixed DG discretization of (2.4): Find and such that

 ah(σh,τ)=κh(σh,τ)A,∀τ∈XDGh. (4.1)
###### Remark 4.1.

We are only considering via (4.1) the symmetric interior penalty DG method (SIP) because the non-symmetric DG versions, known in the literatura as NIP and IIP, have sup-optimal rates of convergence for the eigenvalues [1, 22].

In all what follows, we make the following stability assumption.

###### Assumption 1.

The parameter is greater than or equal to : .

Under this assumption, Proposition permits us to guaranty the well–posedness of the discrete source operator given, for any , by

 ah(~Thf,τh)=(f,τh)A∀τh∈XDGh. (4.2)

Actually, is uniformly bounded, namely,

 (4.3)

Similarly to the continuous case, we observe that is a solution of problem (4.1) if and only if , is an eigenpair of , i.e., . Moreover, it is clear that is an eigenvalue common to (4.1) and with corresponding eigenspace

 Kh:={τh∈XDGh; ch(τh,τh)=0}. (4.4)

The following result establishes a Céa estimate for the DG approximation (4.2) of (2.5).

###### Theorem 4.1.

Under Assumption 1, for all , it holds

 ∣∣∣∣∣∣(~T−~Th)f∣∣∣∣∣∣≤(1+2M)infτh∈XDGh∣∣∣∣∣∣~Tf−τh∣∣∣∣∣∣∗, (4.5)

with as in (3.3).

###### Proof.

We already know from (2.7) that . Hence, using the integration by parts (1.1) elementwise in

 a(~Tf,τ)=(f,τ)A,∀τ∈X

gives

 (~Tf−f,τ)A=(ε(u),τ)=−(u,divhτ)+({u},⟦τ⟧)F∗h∀τ∈XDGh.

Substituting back into the last expression we get

 (~Tf,τ)A+(1ϱdiv(~Tf),divhτ)−({1ϱdiv(~Tf)},⟦τ⟧)F∗h=(f,τ)A∀τ∈XDGh.

Combining the last identity with (4.2) yields the following consistency property

 ah((~T−~Th)f,τh)=0∀τh∈XDGh,∀f∈H. (4.6)

Now, by virtue of (3.3), (3.4) and (4.6), it holds

 ≤ah(~Thf−τh,~Thf−τh)=ah(~Tf−τh,~Thf−τh) ≤M∣∣∣∣∣∣~Tf−τh∣∣∣∣∣∣∗∣∣∣∣∣∣~Thf−τh∣∣∣∣∣∣∀τh∈XDGh,

and the result follows from the triangle inequality.

### 4.1 The operator Jsh

For technical reasons, we want to consider here the -conforming finite element space given by . The goal of this section is to prove that, under certain conditions on the mesh and on the polynomial degree , it holds

 infτh∈Xch∥σ−τh∥X⟶0,when