Guaranteed Lower Bounds for the Elastic Eigenvalues by Using the Nonconforming Crouzeix–Raviart Finite Element

: This paper uses a locking-free nonconforming Crouzeix–Raviart ﬁnite element to solve the planar linear elastic eigenvalue problem with homogeneous pure displacement boundary condition. Making full use of the Poincaré inequality, we obtain the guaranteed lower bounds for eigenvalues. Besides, we also use the nonconforming Crouzeix–Raviart ﬁnite element to the planar linear elastic eigenvalue problem with the pure traction boundary condition, and obtain the guaranteed lower eigenvalue bounds. Finally, numerical experiments with MATLAB program are reported.


Introduction
The linear elasticity discusses how solid objects deform and become internally stressed under prescribed loading conditions, and is widely used in structural analysis and engineering design. It has been well-known that using the finite element methods for the elasticity equations/eigenvalue problems, the displacement field can be determined numerically. There have been many studies on the finite element methods for the elastic equations/eigenvalue problems (e.g., see ). For the elastic equations/eigenvalue problems with pure displacement boundary condition, [1][2][3] studied the conforming finite element methods. As we all know, when λ → ∞ (or Poisson ratio ν → 1 2 ), i.e., the elastic material is nearly incompressible, the locking phenomenon will occur by using the conforming finite elements to solve the equations/eigenvalue problems. In order to overcome the phenomenon of locking, various numerical methods for the linear elasticity equations have been developed. For instance, Brenner et al. [4,5] applied nonconforming Crouzeix-Raviart finite element (CR element). Based on the nonconforming CR element Hansbo [6] proposed a discontinuous Galerkin method, the discontinuous Galerkin method is closely related to the classical nonconforming CR element, which is obtained when one of the stabilizing parameters tends to infinity. Lee et al. [7] proposed a nonconforming Galerkin method based on triangular and quadrilateral elements. Botti et al. [8] constructed a low-order nonconforming approximation method. Rui [9] constructed a finite difference method on staggered grids. Zhang et al. [10] developed the nonconforming virtual element method. Chen et al. [11] presented a nonconforming triangular element and a nonconforming rectangular element. For the elasticity equations with interfaces, Lee [12] adopted the immersed finite element method based on the nonconforming CR element. Jo [13] introduced a low order finite element method for three dimensional elasticity problems. In recent
Let the displacement vector u(x) = (u 1 (x), u 2 (x)) T , and the displacement gradient tensor ∇u be defined by ∇u = ∂ x u 1 ∂ y u 1 ∂ x u 2 ∂ y u 2 .
Consider the elastic eigenvalue problem with homogeneous pure displacement boundary condition: where ρ(x) is the mass density. Without loss of generality, we assume that ρ ≡ 1 throughout this paper. The weak formulation of (1) is: where and A : B = tr(AB T ) is the Frobenius inner product of matrices A and B, and we use A 2 L ∼ 2 (Ω) = Ω A : A dx to denote the L 2 -norm of matrix A. From the following Korn s inequality (see [4], Corollary 11.2.25) ( we know that the bilinear form a(·, ·) is H 1 0 (Ω)-coercive. We use a(·, ·) and · a = a(·, ·) as an inner product and norm on H 1 0 (Ω), respectively, use b(·, ·) and · L 2 (Ω) = b(·, ·) as an inner product and norm on L 2 (Ω), respectively.
Let π h := {κ} be a regular triangular mesh of Ω, and the mesh diameter h = max is the diameter of element κ. Let the set of all interior edges of π h as ε i h , the set of the edges on the boundary as For the nonconforming CR element, the interpolation operator I h : Define Then operator I h : Denote The nonconforming CR element approximation of (2) is: where Korn's inequality for piecewise H 1 -vector fields (see [20]) plays an important role in the existence and uniqueness of the solution for the linear elasticity problem discreted by the discontinuous Galerkin method. For the nonconforming CR element, discrete Equation (7) has a unique solution because a h (·, ·) is positive definite (see page 325 in [4]). In fact, the nonconforming bilinear form a h (·, ·) is symmetric and positive definite on H(h). Because a h (v h , v h ) = 0 implies that v h is piecewise constant, and the zero boundary condition together with continuity at the midpoints imply v h ≡ 0.
Define the nonconforming energy norm · h and the norm | · | 1,h on H(h), respectively: Consider the associated boundary value problem of (2) and discrete form: where Let Ω ⊂ R 2 be a bounded convex polygonal domain, from (11.4.4) in [4] (or Lemma 2.2 in [5]) we have the elliptic regularity estimate: For any f ∈ L 2 (Ω), (9) exists a unique solution w ∈ H 2 (Ω) ∩ H 1 0 (Ω) satisfying where the positive constant C Ω is independent of Lamé parameters µ and λ. The above inequality (11) plays an essential role in showing the robustness of our numerical approximation to (9).
Brenner et al. in [4,5] studied and proved the following estimates:

Proposition 1.
Let Ω ⊂ R 2 be a bounded convex polygonal domain, w and w h be the kth eigenvalues of (9) and (10), respectively. There exists a positive constant C such that where C is independent of h and (µ, λ), which indicates that the nonconforming CR element method is locking-free.
Proof of Proposition 1. See the proof method of Theorem 3.1 on page 331 and Theorem 3.2 on page 332 in [5] or the proof method of Theorem 11.4.15 on page 325 in [4] to prove (12) and (13).
Define the operator T : and T h : Then both T and T h are self-adjoint and completely continuous operators. (2) and (7) have the following equivalent operator forms, respectively:

The Lower Bounds for the Eigenvalues of the Pure Displacement Problem
In [42] Liu established a framework that provides guaranteed lower eigenvalue bounds for the self-adjoint eigenvalue problems. We verify all conditions of the framework and apply it to obtain lower bounds for the eigenvalues of the planar linear elastic eigenvalue problem with the pure displacement boundary. A1 V is a Hilbert space of real function on Ω with the inner product M(·, ·) and the corresponding norm · M := M(·, ·). A2 N(·, ·) is another inner product of V. The corresponding norm · N := N(·, ·) is compact for V with respect to · M , i.e., every sequence in V which is bounded in · M has a subsequence, which is Cauchy in · N .
The assumption A4 means that M h (·, ·) and N h (·, ·) are inner products of V(h). In order to simplicity, the extended bilinear forms M h (·, ·) and N h (·, ·) are still denoted by M(·, ·) and N(·, ·), and the corresponding norms are denoted by · M and · N , respectively.
Under the above assumptions, in [42] Liu gave the following Lemma.
Lemma 1 (the Theorem 2.1 in [42]). Suppose that γ and γ h are the kth eigenvalues of (2) and (7), respectively, Suppose the following error estimation holds, ∀u ∈ V Then, there holds For the pure displacement problem, we take the following settings: It is easy to verify that the above settings satisfy the assumption A1 − A4. In the following discussion, for consistency of notations, we use Theorem 1 (guaranteed lower bounds for eigenvalues). Let γ k and γ h,k be the kth eigenvalue of (2) and (7), respectively, then the following guaranteed lower bounds holds where Proof of Theorem 1. For any u ∈ H(h), v h ∈ V h (Ω), using integration by parts and noticing that ∂v h ∂n is a constant function on edges (n is the unit outward normal vector), v h is a linear function on element κ and ∆ h v h = 0, which together with (6) we have thus the interpolation operator I h is the orthogonal projection. For κ ∈ π h , the edges of which are denoted by e 1 , e 2 , e 3 , let The following the Poincaré inequality with an explicit bound (see Theorem 3.2 in [42]) plays an crucial role in studying guaranteed lower eigenvalue bounds.

Remark 1.
Actually, when the angles of meshes meet contain condition (see Theorem 4.2 in [42] for details), such as uniform meshes used in our numerical experiments, the value of C h can be 0.1893h 1 √ µ .

The Pure Traction Problem
In this section, we present the nonconforming CR finite element for the planar linear elastic eigenvalue problem with the pure traction boundary condition. Unless explicitly noted in this section, we use the same notation as in Section 2.
We consider the planar linear elastic eigenvalue problem with the pure traction boundary: where n is the unit outward normal vector with respect to the domain Ω.
The weak formulation of (27) can be described as to where ω = γ + 1 and By the Korn's inequality (see [4], Theorem 11.2.16) ( we can know that the bilinear form a(·, ·) is H 1 (Ω)-coercive. We use a(·, ·) and · a = a(·, ·) as an inner product and norm on H 1 (Ω), respectively. Let V h (Ω) be the nonconforming CR element space: The nonconforming CR element approximation of (28) is: where ω k = γ k + 1 and It is easy to know that the nonconforming bilinear form a h (·, ·) is symmetric and positive definite on H(h).
Define the nonconforming energy norm · h and the norm | · | 1,h on H(h), respectively: Consider the associated boundary value problem of (28) and discrete form: Define the operator K : and Then both K and K h are self-adjoint and completely continuous operators.

The Lower Bounds for Eigenvalues of the Pure Traction Problem
For the elastic eigenvalue problem with the pure traction boundary condition, we also use , M(·, ·), N(·, ·), · M , · N , respectively, which also satisfy the assumption A1-A4. Let Using the argument in [44] we give the estimate between the interpolation and projection operators.
be the interpolation operator of the nonconforming CR element, for all u ∈ H(h), we have Proof of Lemma 2. For the convenience of readers, we refer to the Lemma 3.6 in [44] to give a detailed proof. From (34), for any Let v h := K h Φ h in the above formula, then From the definition of ω h,1 in (30) and the min-max principle, we get which together with (38) yields the estimate For any u ∈ H(h), using the same argument as (22) we get Let v h := I h u − P h u and Ψ h := K h v h ∈ V h (Ω). Combining (35), (37) and Let v h = Ψ h in (41), together with (40) and (42) From (43) we can immediately obtain (36).

Theorem 2 (guaranteed lower bounds for eigenvalues).
Let ω k and ω h,k be as defined in (28) and (30), respectively. Then, there holds where Proof of Theorem 2. For any u ∈ H(h), from (23) we get the following estimates: From the triangle inequality, (36) and (45) we obtain From Lemma 1 we can immediately obtain (44).

Remark 2.
Same as Remark 1, when the angles of meshes meet contain condition the value of C h can be ( 1

Numerical Experiments
In this section, we report some numerical experiments. In our computation, our program was completed under the package of iFEM [46], and the discrete eigenvalue problems were solved in MATLAB 2012a on a DELL inspiron 5480 PC with 8G memory, and the MATLAB codes are in Appendix A. The following notations are adopted in tables.
h: the meshes h = max γ h,k : the kth eigenvalue of (7) obtained by the nonconforming CR element. γ h,k : the approximation obtained by correcting γ h,k , i.e., the lower bounds of the kth eigenvalue of (1) or (27). γ c h,k : the kth eigenvalue obtained by the linear conforming finite element. γ glb : guaranteed lower bounds for the elastic eigenvalues. cond 2 (CR): the condition number of the stiffness matrix of (30) obtained by the nonconforming CR element. cond 2 (P1): the condition number of the stiffness matrix obtained by the linear conforming element.
Example 1. Consider the planar linear elastic eigenvalue problem with the pure displacement boundary condition (1), we take the mass density ρ = 1, Lamé parameter µ = 1 and take λ = 1, λ = 100, λ = 10 4 , λ = 10 8 , respectively. We use the nonconforming CR element to solve (1) on the unit square Ω S = [0, 1] 2 and L-shaped domain Ω L = [0, 1] 2 \ [ 1 2 , 1] 2 , respectively. On each domain, we select two eigenvalues to execute correction. One is the minimum eigenvalue, and the other is selected because it is an upper bound of the exact value on the coarsest mesh. For the uniform meshes, C h = 0.1893h 1 √ µ is used for the results in Tables 1-4. For the nonuniform meshes, the value of C h taken as 0.346h 1 √ µ . In addition, we show the meshes with geometrically triangular subdivision in Figure 1. In order to observe the influence of the Lamé parameter λ, we use the nonconforming CR element and the linear conforming finite element to solve (1) by taking λ = 1, λ = 10 3 , λ = 10 5 , λ = 10 8 , λ = 10 10 while µ = 1 and h = √ 2 256 , and depict the curves of the corrected eigenvalue γ h,1 and approximations γ h,1 and γ c h,1 of (1) in Figure 2. From Tables 1, 3 and 5, we see that γ h,k approximate the exact ones from below, and from Tables 2, 4 and 6, we see that on the coarsest triangulation of the square and the L-shaped domain, the discrete eigenvalues γ h,k computed by the nonconforming CR element cannot be a lower bound for the exact ones. But after correction, it must be the lower bound for eigenvalue. The corrected eigenvalues all converge to the exact ones from below. Combining with the analysis of the Morley finite element in [40], we know that the corrected eigenvalues are the guaranteed lower bounds, which coincide with our theoretical result of Theorem 1. Furthermore, from Figure 2, we can see that the approximations for the first eigenvalue of (1) obtained by the linear finite element become large as λ increases, but the approximations for the first eigenvalue of (1) obtained by the nonconforming CR element and the corrected eigenvalues become stable as λ increases, which indicates the nonconforming CR element method and the method of obtaining the guaranteed lower eigenvalue bounds are locking-free.      Tables 7-10. We also use the nonconforming CR element and the linear conforming finite element to solve (27) by taking λ = 1, λ = 10 3 , λ = 10 5 , λ = 10 8 , λ = 10 10 while µ = 1 and h = √ 2 256 , and depict the curves of the corrected eigenvalues and approximations for the first nonzero eigenvalue of (27) in Figure 3.     From Tables 7-10, we see that γ h,k approximate the exact ones from below. From Tables 8 and 10, we know that the eigenvalue obtained on the coarsest grid cannot be a lower bound for the exact ones, but after correction, it must be the lower bound for eigenvalue. The corrected eigenvalues converge to the exact ones from below, which coincide with our theoretical results. From Figure 3, we can see that curves of the corrected eigenvalues are parallel to that of the uncorrected eigenvalues, and the corrected eigenvalues are locking-free. As λ increases, the approximations for eigenvalue of (27) obtained by the nonconforming CR element and the linear conforming finite element are locking-free. Furthermore, from Figure 3, we find that the eigenvalues are abnormal as λ increases, so we compute the condition number of the stiffness matrix, and the results are shown in Tables 11 and 12. From Tables 11 and 12, we can see that the condition number increases as λ increases, we think it may be because the influence of the condition number and the rounding error.  cond 2 (CR) 1.004 × 10 7 3.408 × 10 9 3.395 × 10 11 3.416 × 10 14 3.179 × 10 16 cond 2 (P1) 1.681 × 10 6 5.708 × 10 8 5.636 × 10 10 5.618 × 10 13 5.191 × 10 15 Remark 3. In Tables 1-10, we can see that the eigenvalues selected are the guaranteed lower bounds after correction. We take the corrected eigenvalues on the smallest mesh size as the guaranteed lower bounds.

Conclusions
Generally, we know that the exact eigenvalues of the planar linear elastic eigenvalue problem are unknown, according to the min-max principle we can obtain upper bounds for eigenvalues by using conforming finite element. But it is generally more difficult to obtain lower bounds for the numerical eigenvalues. Noticing that the lower and upper bounds can produce intervals to which exact eigenvalue belongs. This is important for the design of the coefficient of safety in practical engineering. Therefore, we apply a locking-free nonconforming CR element to the planar linear elastic eigenvalue problem, and obtain the guaranteed lower bounds for eigenvalues.
In this paper, for the planar linear elastic eigenvalue problem with the pure displacement and the pure traction boundary conditions in two spacial dimensions, we prove that using the nonconforming CR element can obtain the guaranteed lower bounds for exact eigenvalues (see Theorems 1 and 2), and that is locking-free when the elliptic regularity estimate (11) holds. Besides, the discussion of the guaranteed lower eigenvalue bounds in this paper can be extended to mixed boundary-value problem and three spacial dimensions, and the schemes are locking-free when the elliptic regularity estimate holds, i.e., the constant C Ω is independent of µ and λ. Furthermore, it is a meaningful and difficult work to extend to higher order elements.