Locking-free nonconforming finite elements for three-dimensional elasticity problem

Applied Mathematics and Computation 217 (2011) 5790–5797

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Locking-free nonconforming finite elements for three-dimensional elasticity problem Liu-Chao Xiao a,⇑, Yong-Qin Yang b, Shao-Chun Chen b a b

College of Science, Henan University of Technology, Zhengzhou 450001, China Department of Mathematics, Zhengzhou University, Zhengzhou 450001, China

a r t i c l e

i n f o

Keywords: Three-dimensional elasticity Locking-free Nonconforming finite element Error estimate Uniform convergence

a b s t r a c t Two locking-free nonconforming finite elements are presented for three-dimensional elasticity problem with pure displacement boundary condition. Convergence rate of the elements are uniformly optimal with respect to k. The energy norm and L2 norm errors are O(h2) and O(h3), respectively. Lastly, a numerical experiment is carried out, which coincides with the theoretical analysis. Ó 2010 Elsevier Inc. All rights reserved.

1. Introduction For the homogeneous isotropic linear elasticity, it is well known that the convergence rate for the standard displacement method using normal lower order conforming finite elements deteriorates as the Lamé constant k ? 1, i.e., as the elastic material is incompressible. This is known as the phenomenon of locking [1–5], which is caused by that the kernel of the discrete div space is too small, or the factor in the error estimate tends to 1 as k ? 1. In order to avoid the phenomenon of locking, we should construct some special elements such that the finite element solution uniformly converge to the true solution with respect to k. Various approaches have been proposed which work uniformly well for all k. Many works have been done using mixed finite element methods [3,4,6,7], which have one drawback involving large number of variables. Another method resorts to nonconforming finite element methods [5,8–11]. The resulted matrix is positive definite and easier to be solved. However, there are a few papers dealing with the locking phenomenon of three-dimensional elasticity by finite element methods [10,12]. In this paper, we construct two locking-free nonconforming finite elements, and prove that the convergence rate of the finite elements is uniformly optimal with respect to k for three-dimensional elasticity with pure displacement boundary condition. The energy norm and L2 norm errors are O(h2) and O(h3), respectively. Lastly, we present a numerical experiment to verify our theoretical results. 2. Nonconforming locking-free elements The pure displacement boundary value problem for a three-dimensional homogeneous isotropic material is given by

(

div rð~ uÞ ¼ ~ f in X; ~ u ¼~ 0 on @ X;

⇑ Corresponding author. E-mail addresses: [email protected] (L.-C. Xiao), [email protected] (Y.-Q. Yang), [email protected] (S.-C. Chen). 0096-3003/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2010.12.061

ð1Þ

L.-C. Xiao et al. / Applied Mathematics and Computation 217 (2011) 5790–5797

5791

where X  R3 is a bounded convex polyhedron with the boundary oX, ~ u ¼ ðu1 ; u2 ; u3 Þt is the displacement, ~ f 2 ðL2 ðXÞÞ3 is the body force, k 2 (0, 1), l 2 [l1, l2] are Lamé constants, 0 < l1 < l2,

rð~ uÞ ¼ 2lð~ uÞ þ kdiv~ uI; div~ u¼

@u1 @u2 @u3 þ þ ; @x @y @z

ð~uÞ ¼

1 ðgrad~ u þ ðgrad~ uÞt Þ; 2

where I is the unit matrix. (1) is equivalent to the following partial differential equation

(

lM~ u  ðl þ kÞgradðdiv~ uÞ ¼ ~ f in X; ~ u ¼~ 0 on @ X;

ð2Þ

whose weak form is

(

Find ~ u 2 V; such that ~ ~ aðu; v Þ ¼ ð~ f;~ v Þ 8~ v 2 V;

ð3Þ

 3 R where V ¼ H10 ðXÞ ; ð~ f;~ v Þ ¼ X ~f  ~ v dxdydz,

að~ u; ~ vÞ ¼

Z

flgrad~ u : grad~ v þ ðl þ kÞðdiv~uÞðdiv ~ v Þgdxdydz ) Z ( X 3 ¼ l gradui  gradv i þ ðl þ kÞðdiv~ v Þ dxdydz: uÞðdiv ~ X

X

ð4Þ

i¼1

By Poincáre’s inequality, the bilinear form in (4) is symmetric, coercive and continuous on V. So there exists a unique solution to problem (3) by Lax–Milgram Theorem [13,14]. For the sake of simplicity, we do not discriminate the positive constant C, which may denote different constant in the different position. We denote by Pk the space of polynomials with degree 6k, and by Qk the space of polynomials that are of degree 6 k with respect to each variable. 2.1. Tetrahedron element Let T h be a regular tetrahedron subdivision for the polyhedral domain X. For any T 2 T h ; Ai ðxi ; yi ; zi Þ and Fi (i = 1, 2, 3, 4) are the vertices and faces, respectively. (k1, k2, k3, k4) is the capacity coordinate. We denote hT ¼ max16i64 diamF i ; h ¼ b be the tetrahedron in the (n, g, f) space. The four vertices and faces are maxT2T h hT . Let the reference element T b 1 ð0; 0; 0Þ; A b 2 ð1; 0; 0Þ; A b 3 ð0; 1; 0Þ; A b 4 ð0; 0; 1Þ; b A F i ði ¼ 1; 2; 3; 4Þ, respectively. It is obvious that

^k1 ¼ 1  n  g  f;

^k2 ¼ n;

^k3 ¼ g;

^k4 ¼ f:

b Þ as following: the freedoms are b ; P; b R We define finite element ð T

X d

8 R 1 ~ > F i Þ i ¼ 1; 2; 3; 4; < b bF v^ pd^s 8p 2 P1 ð b jF ij i ¼ R > : 1 b div ~ v^ qdndgdf 8q 2 P1 ð Tb Þ jb Tj T

ð5Þ

and the shape function space is

b Tb Þ ¼ f~ Pð v^ ; ~ v^ 2 ðWð Tb ÞÞ3 ; div ~ v^ 2 P1 ð Tb Þg;

ð6Þ

n b Þ ¼ P2 ð T bÞ  ^ k2  ^ k22  ^ k2 þ ^ k22 ; ^ k3  ^ k23  ^ k2 þ ^ k22 ; ^ k4  ^ k24  ^ k2 þ ^ k22 ; ^ k1  ^ k21  ^ k2 þ where Wð T k21 ^ k1 ^ k24 ^ k4 ^ k22 ^ k2 ^ k24 ^ k4 ^ k23 ^ k3 ^ k24 ^ k4 ^ k24 ^ k4 ^ k24 ^ ^ k23 ^ k3 ^ k24 ^ k4 ^ k22 ; ^ k1  ^ k21  ^ k2 þ ^ k22 g. k4 ^ The freedom and the shape space are 39 dimensions, and the element is well defined by simple calculation. 2.2. Cuboid element Assume X to be a cuboid domain and T h to be one of its regular cuboid subdivision. For any given T 2 T h , we denotes A0(x0, y0, z0), Fi (i = 1, 2, . . . , 6), 2h1, 2h2, 2h3 are the central point, faces of T, the length of sides corresponding to x, y and z direction, respectively. Let hT ¼ maxðh1 ; h2 ; h3 Þ; h ¼ maxT2T h hT . We denote ~ m ¼ ðm1 ; m2 ; m3 Þt as the outer normal. Suppose reference b ¼ ½1; 13 ; b b defined as element T F i ð1 6 i 6 6Þ are faces of T

5792

L.-C. Xiao et al. / Applied Mathematics and Computation 217 (2011) 5790–5797

b F 1 ¼ fðn; g; fÞ 2 b F 3 ¼ fðn; g; fÞ 2 b F 5 ¼ fðn; g; fÞ 2

Tb j n ¼ 1g; b F 2 ¼ fðn; g; fÞ 2 Tb j g ¼ 1g; Tb j f ¼ 1g; b F 4 ¼ fðn; g; fÞ 2 Tb j n ¼ 1g; Tb j g ¼ 1g; b F 6 ¼ fðn; g; fÞ 2 Tb j f ¼ 1g:

b as following: We define spaces of polynomials on T

W 1 ð Tb Þ ¼ ðP3 n ff3 gÞ  fn2 g2 ; n2 f2 ; n3 g; n3 f; g3 f; gf3 ; ng2 f; ngf2 ; n4 g; n4 f; ng3 f; ngf3 g; W 2 ð Tb Þ ¼ ðP3 n fn3 gÞ  fg2 f2 ; g2 n2 ; g3 f; g3 n; f3 n; n3 f; ngf2 ; n2 gf; g4 f; g4 n; ngf3 ; n3 gfg; W 3 ð Tb Þ ¼ ðP3 n fg3 gÞ  ff2 n2 ; f2 g2 ; f3 n; f3 g; n3 g; ng3 ; n2 gf; ng2 f; f4 n; f4 g; n3 gf; ng3 fg: The shape function space is taken as

b Tb Þ ¼ f~ v^ ; ~ v^ 2 W 1 ð Tb Þ  W 2 ð Tb Þ  W 3 ð Tb Þ; div ~ v^ 2 P1 ð Tb Þg: Pð

ð7Þ

b T b Þ ¼ 75. Let ~ It is easy to see that dim Pð v^ ¼ ðv^ 1 ; v^ 2 ; v^ 3 Þt , the degrees of freedom are

X d

8 R 1 ~ ^ > < b bF v^ pd^s 8p 2 Q 1 ðF i Þ; 1 6 i 6 6; jF i j i ¼ R > : 1 b div ~ v^ qdndgdf 8q 2 P1 ð Tb Þ: b jT j T

ð8Þ

The above 75-freedom finite element is also well defined. b ! T be: Let affine transformation F T : T

ðx; y; zÞt ¼ Bðn; g; fÞt þ ðx1 ; y1 ; z1 Þt ; 0

x2  x1 where B ¼ @ y2  y1 z2  z1

ð9Þ

1 x4  x1 y4  y1 A for the tetrahedron element, or z4  z1

x3  x1 y3  y1 z3  z1

ðx; y; zÞt ¼ Bðn; g; fÞt þ ðx0 ; y0 ; z0 Þt ; 0

h1 where B ¼ @ 0 0

0 h2 0

ð10Þ

1

0 0 A or the cuboid element. h3

For any ~ v^ on Tb , Piola’s transformation reads

~ v ¼ B~ v^ ;

ð11Þ

~ b b ^  F 1 ^ b b then the shape function space on T is PðTÞ ¼ f~ p; ~ p ¼ B~ p T ; p 2 Pð T Þg, here Pð T Þ is defined by (6) or (7). b T b Þ and P(T), respectively. Obviously, the b Suppose P and PT are the finite element interpolation operators introduced by Pð b~ interpolation operators are affine equivalent, i.e. Pd v ¼ B1 PT ~ v¼P v^ , and the following properties hold for PT: T~

1 jF i j 1 jTj

Z Z

Fi

T

~ v pds ¼

1 jF i j

Z

ðPT ~ v Þpds 8p 2 PðF i Þ;

Fi

div ~ v  qdxdydz ¼

1 jTj

Z

ð12Þ

div PT ~ v  qdxdydz 8q 2 P1 ðTÞ;

ð13Þ

T

where P(Fi) = P1(Fi) (1 6 i 6 4) for the tetrahedron element, P(Fi) = Q1(Fi) (1 6 i 6 6) for the cuboid element. For any ~ v 2 ðH1 ðXÞÞ3 , the global finite element interpolation operator Ph is defined as

Ph~ v jT ¼ PT ~ v 8T 2 T h :

ð14Þ

The nonconforming finite element space corresponding to Ph is

  Z Vh ¼ ~ v h; ~ v h jT 2 PðTÞ; ~ v h  ~qds ¼ 0; 8F  @ X; 8~q 2 ðPðFÞÞ3 ; 8T 2 T h :

ð15Þ

F

It follows easily from (12) that

Z

½~ v h  ~qds ¼ 0 8~ vh 2 V h;

~ q 2 ðPðFÞÞ3 ;

ð16Þ

F

where we denote [v] as the jump of the function

v on the face F.

L.-C. Xiao et al. / Applied Mathematics and Computation 217 (2011) 5790–5797

5793

We define the discrete bilinear form and energy norm as

ah ð~ uh ; ~ v hÞ ¼

XZ T

flgrad~ uh : grad~ v h þ ðl þ kÞdiv~uh  div ~ v h gdxdydz;

ð17Þ

T

1

k~ v h kh ¼ ah ð~ vh; ~ v h Þ2 . Then the approximation of (3) is

(

Find ~ uh 2 V h ; such that ~ v h Þ ¼ ð~f ; ~ v h Þ 8~ v h 2 V h: ah ðuh ; ~

ð18Þ

Clearly, k~ v h kh is a norm on Vh, so there exists a unique solution to problem (18).

3. Error estimates

v ¼ div ~ v^ ; 8~ v 2 ðH1 ðXÞÞ3 . Lemma 3.1. div ~ ^ then r ^ ¼ grad, ^ ¼ Bt r, Proof. Let r ¼ grad; r

^ ~ div ~ v^ ¼ r v^ ¼ Bt r  ~ v^ ¼ r  B~ v^ ¼ r  ~ v ¼ div ~ v: By well known technique [14] we can obtain the following estimates.

h

v 2 ðW m;p ðTÞÞ2 ; then ~ v^ 2 ðW m;p ð Tb ÞÞ2 , and Lemma 3.2. If ~ j~ v^ j

1

m;p;b T

6 CkBkm kB1 kjdetBjp j~ v jm;p;T ;

ð19Þ

where kBk is the Euclid norm of B. Similarly, 1 j~ v jm;p;T 6 CkB1 km kBkjdetBjp j~ v^ j

m;p;b T

ð20Þ

:

Define the projection operator cT : L2(T) ? P1(T) as following: For given w 2 L2(T), cTw 2 P1(T) and

ðw  cT w; pÞ ¼ 0;

8p 2 P1 ðTÞ:

ð21Þ

Using Bramble–Hilbert lemma [13] and the affine equivalent finite element methods [14], one can prove the following result. Lemma 3.3. Suppose cT be the operator of (21), then 2

kw  cT wkL2 ðTÞ 6 ChT jwj2;T

8w 2 H2 ðTÞ;

T 2 T h;

ð22Þ

where C = Const > 0 independent of hT. Lemma 3.4. For any T 2 T h , we have

div PT ~ v ¼ cT div ~ v 8~ v 2 ðH1 ðTÞÞ3 :

ð23Þ

Proof. From the definition of PT and Lemma 3.1, we have div PT ~ v 2 P1 ðTÞ. By (12) and (13), using Green’s formula, we have

Z T

div PT ~ v  pdxdydz ¼

Z

div ~ v  pdxdydz 8p 2 P1 ðTÞ:

T

By the definition of cT, (23) can be easily gotten. h Lemma 3.5. There exists a positive C, independent of hT and T, such that

k~ v  PT ~ v k0;T þ hT j~ v  PT ~ v j1;T 6 Ch3T j~ v j3;T

8~ v 2 ðH3 ðTÞÞ3 :

Proof. By Lemma 3.2 and interpolation theorem [14], we have

ð24Þ

5794

L.-C. Xiao et al. / Applied Mathematics and Computation 217 (2011) 5790–5797 3

3

b~ b~ k~ v  PT ~ v k0;T þ hT j~ v  PT ~ v j1;T 6 Ch2T kBkðk~ v^  P v^ k0;bT þ j~ v^  P v^ j1;bT Þ 6 Ch2T kBkj~ v^ j3;bT 6 Ch3T j~ v j3;T : We assume that

k~ uk2;X þ kkdiv~ f k0;X ; uk1;X 6 Ck~

ð25Þ

k~ uk3;X þ kkdiv~ uk2;X 6 C;

ð26Þ

where ~ u is the solution to problem (2), C = Const > 0 independent of k. h Remark 3.6. (25) and (26) are true for the planar elasticity when the boundary oX is smooth sufficiently. As for the threedimensional case (see [12]), we have not found the results. We use them as assumptions. Furthermore, a numerical example is provided to verify our theoretical analysis. uh be the solution to problem (3) and (18) respectively, ~ f 2 ðL2 ðXÞÞ3 . Then we have Theorem 3.7. Let ~ u 2 ðH3 ðXÞ \ H10 ðXÞÞ3 and ~ the following error estimate 2 k~ u ~ uh kh 6 Ch ;

ð27Þ

where C = Const > 0 independent of h and k. Proof. By the second Strang lemma [7],

(

) ~ h Þ  ð~ ~ h Þj jah ð~ u; w f;w : k~ u ~ uh kh 6 C inf k~ u~ v h kh þ sup ~ h kh ~ v h 2V h kw ~ h 2V h 0–w

ð28Þ

Using Lemma (3.3)–(3.5) and the regularity assumption (26), we obtain

inf k~ u~ v h k2h 6 k~u  Ph~uk2h ¼ l

X

~ v h 2V h

6 Ch

4

X

j~ u  PT ~ uj21;T þ ðl þ kÞ

T

j~ uj23;T þ ðl þ kÞ

X

T

X

kdiv~ u  div PT ~ uk20;T

T 4 4 kdiv~ uj23;X þ kjdiv~ u  cT div~ uk20;T 6 Ch ðj~ uj22;X Þ 6 Ch :

ð29Þ

T

The following is the nonconforming error estimate. By Green’s formula,

~ h Þ ¼ ah ð~ ~ h Þ  ð~ ~h Þ ¼ l Eh ð~ u; w u; w f;w

XZ T

@T

~ h ds þ ðl þ kÞ uw @ m~

XZ T

~h  ~ uw div~ mds:

ð30Þ

@T

R 1 Following the standard nonconforming error estimate method [15], let P F0 ðwÞ ¼ jFj wds. For the subdivision T h on X, let F Ih(w) be the piecewise linear (for tetrahedron element) or cubic-linear (for cuboid element) interpolation operator defined by the value of w 2 H2(X) at the element vertices, then Ih(w) 2 C0(X). From (15) and (16), by trace theorem and interpolation theorem [14] we have

    X Z  X X X  3 Z     F ~ ~ ~ ~ ~ ~ @ m u  wh ds ¼  ½@ i u  Ih ð@ i uÞ  ½wh  P0 ðwh Þ  mi ds   T @T   T F@T i¼1 F  6

3 X XX T

~ h  PF0 ðw ~ h Þk0;@T 6 C u  Ih ð@ i~ uÞk0;@T kw k@ i~

i¼1 F@T

3 XX T

3

1

~ h j1;T uj2;T h2T jw h2T j@ i~

i¼1

2

~ h kh : 6 Ch j~ uj3;X kw

ð31Þ

Similarly,

    X Z  X X Z     F ~h  ~ ~ h  P0 ðw ~ h Þ  ~ ~ h kh : uw u  Ih ðdiv~ uÞ  ½w uj2;X kw div~ mds ¼  ½div~ mds 6 Ch2 jdiv~   T @T   T F@T F 

ð32Þ

Substituting (31), (32) into (30) we have

~ h Þj 6 Ch2 ðj~ ~ h kh : jEh ð~ u; w uj3;X þ kjdiv~ uj2;X Þkw By the regularity assumption (26), then

ð33Þ

L.-C. Xiao et al. / Applied Mathematics and Computation 217 (2011) 5790–5797

sup ~ h 2V h 0–w

~hÞ Eh ð~ u; w 2 6 Ch : ~ h kh kw

5795

ð34Þ

The theorem now follows by combining (28), (29) and (34). h Next, we present the L2 norm error estimate by the dual argument. Lemma 3.8. Let ~ u 2 ðH3 ðTÞÞ2 ; ~ v 2 ðH1 ðTÞÞ2 ; T 2 T h ; e  @T, then there exists a constant C independent of T such that

 Z    6 Ch3 j~  ~ ~ ~ ~ v  ð u  P u Þds T T v j1;T juj3;T ;   e Z    3  ~ u  PT ~ v j2;T j~uj2;T : uÞds 6 ChT j~  v  ð~

ð35Þ ð36Þ

e

b ; e ¼ F T ð^ b to T. Define Proof. Let ^ e  @T eÞ; B is the transformation matrix from T

bð~ v^ ; ~u^Þ ¼

Z

b~ ^ Þd^s: u B~ v^  Bð~u^  P

ð37Þ

^e

By trace theorem and interpolation theorem [14], we have

b~ b~ b~ ^ Þk0;^e 6 CkBk2 k~ ^ k0;^e 6 CkBk2 k~ ^ k 6 CkBk2 k~ u u u jbð~ v^ ; ~u^Þj 6 kB~ v^ k0;^e kBð~u^  P v^ k0;^e k~u^  P v^ k1;bT k~u^  P v^ k1;bT k~u^k3;bT ; 1;b T

ð38Þ

b ÞÞ2  ðH3 ð T b ÞÞ2 , and kbk 6 CkBk2. It follows from the definition of interpolation i.e., b(, ) is a bounded bilinear form on ðH1 ð T b that operator P

^ 2 ðP1 ð Tb ÞÞ2 ; p 8~

^ 2 ðH3 ð Tb ÞÞ2 ; u 8~

^; ~ ^ Þ ¼ 0: bð~ p u

ð39Þ

8~ v^ 2 ðH1 ð Tb ÞÞ2 ;

^ 2 ðP2 ð~ q 8~ Tb ÞÞ2 ;

bð~ v^ ; ~q^Þ ¼ 0:

ð40Þ

jbð~ v^ ; ~u^Þj 6 Ckbkj~ v^ j1;bT j~u^j3;bT 6 CkBk2 j~ v^ j1;bT j~u^j3;bT ; jbð~ v^ ; ~u^Þj 6 Ckbkj~ v^ j2;bT j~u^j2;bT 6 CkBk2 j~ v^ j2;bT j~u^j2;bT :

ð41Þ

By the bilinear lemma [14],

ð42Þ

Using lemma 3.2 and regularity condition, we have

Z  Z   jej ~ ~   jej  ~ ~ ~ b   ¼ jbðv^ ; u  ~  ^ ^ ~ ^ ^ ^ Þj 6 Ch2T kBk2 j~ ~ ¼ u Þd s v  ð u  P B v  Bð u  P v^ j1;bT j~u^j3;bT 6 Ch3T j~ v j1;T j~uj3;T : u Þds T  j^ej   j^ej  ^ e e Similarly, we have j

R

e

~ v  ð~u  PT~uÞdsj 6 Ch2T kBk2 j~ v^ j2;bT j~u^j2;bT 6 Ch3T j~ v j2;T j~uj2;T . The lemma is proved. h

By the regularity assumption (25), using the general nonconforming error estimate methods [13,14], we also have the following estimate. u 2 ðH2 ðXÞ \ H10 ðXÞÞ3 is the solution to problem (3), ~ uh is the solution to approximation Lemma 3.9. Suppose ~ f 2 ðL2 ðXÞÞ3 ; ~ problem (18). Then

k~ u ~ uh kh 6 Chk~ f k0;X ;

k~ u  Ph~ f k0;X ; ukh 6 Chk~

ð43Þ

where C = Const > 0 independent of h and k. u 2 ðH3 ðXÞ \ H10 ðXÞÞ3 be the solution to problem (3), ~ uh be the solution to approximation Theorem 3.10. Suppose ~ f 2 ðL2 ðXÞÞ3 ; ~ problem (18). Then 3 k~ u ~ uh k0;X 6 Ch ;

ð44Þ

where C = Const > 0 independent of h and k. Proof. From the definition of norm, we have

k~ u ~ uh k0;X ¼

sup ~ 2ðL2 ðXÞÞ3 0–w

j

R

~~ ~ dxj uh Þ  w : ~ k0;X kw

X ðu

ð45Þ

5796

L.-C. Xiao et al. / Applied Mathematics and Computation 217 (2011) 5790–5797

~ 2 ðL2 ðXÞÞ3 , let ~ For given w f 2 ðH2 ðXÞ

(

T

H10 ðXÞÞ3 be the solution of the following equation

~ in X; fÞ ¼ w lM~ f  ðl þ kÞgradðdiv~ ~ f ¼~ 0 on @ X:

ð46Þ

Then ~ f satisfies the variational equation

að~ f; ~ vÞ ¼

Z

~ ~ v dx 8~ v 2 V: w

ð47Þ

X

Suppose ~ fh 2 V h be the solution to the following discrete equation

ah ð~ fh ; ~ vÞ ¼

Z

~ ~ w v dx 8~ v 2 V h:

ð48Þ

X

From Lemma 3.9 and (25), there exists a positive constant C independent of h and k, such that

~ k0;X fk1;X 6 Ckw k~ fk2;X þ kkdiv~

ð49Þ

~ k0;X ; k~ f ~ fh kh 6 Chkw

ð50Þ

~ k0;X : fkh 6 Chkw k~ f  Ph~

By (47) and (48) we have

Z     ð~ ~ dx ¼ jah ð~ f; ~ u ~ uh Þ þ ah ðPh~ f; ~ u ~ uh Þj uh Þ  w f ~ fh ; ~ u  Ph~ f ~ fh ; Ph~ fh  Ph~ uÞ þ ah ð~ uÞ þ ah ð~  u ~ X

fkh k~ f; ~ u ~ uh Þj: 6 k~ f ~ fh kh k~ u  Ph~ fh  Ph~ u ~ uh kh þ jah ð~ f ~ fh ; Ph~ ukh þ k~ uÞj þ jah ðPh~

ð51Þ

Using (29) and (50), we obtain 3 ~ k0;X ; k~ f ~ fh kh k~ u  Ph~ ukh 6 Ch kw

ð52Þ

3 ~ k0;X : fkh k~ fkh Þk~ u ~ uh kh 6 ðk~ fh  ~ fkh þ k~ f  Ph~ u ~ uh kh 6 Ch kw k~ fh  Ph~

ð53Þ

Using Green’s formula, Lemma 3.8, (25) and (26) for the last two items of (51), we have

 Z      ~ ~  ~  Ph~ uÞ ¼ ah ð~ uÞ  w udx f; Ph~ ah ðf  fh ; Ph~ X       X Z X Z     6 l fÞ~ ðgrad~ fÞ~ m  ðPh~ ðdiv~ m  ðPh~ u ~ uÞds þ ðl þ kÞ u ~ uÞds    T @T  T @T 3 3 3 3 ~ k0;X : 6 C lh k~ fk1;X j~ fk1;X Þ 6 Ch kw fk2;X j~ uj3;X þ Cðl þ kÞh kdiv~ uj3;X 6 Ch j~ uj3;X ðk~ fk2;X þ kkdiv~

ð54Þ

Similarly, 3 3 3 3 ~ k0;X : f; ~ u ~ uh Þj 6 C lh kgrad~ jah ðPh~ uk2;X j~ fj2;X þ Cðl þ kÞh kdiv~ fj2;X 6 Ch j~ fj2;X ðk~ uk3;X þ kkdiv~ uk2;X j~ uk2;X Þ 6 Ch kw

ð55Þ Substituting (51)–(55) into (45), (44) can be easily gotten. h

4. Numerical experiment In this section, we only list the numerical information of the cuboid nonconforming finite element. The tetrahedron element is also numerically convergent uniformly with respect to k 2 (0, +1). Let X ¼ ½0; 13 ; l ¼ 1; ~ f ¼ ðf1 ; f2 ; f3 Þt , where

f1 ¼ 400lð2y  1Þð2z  1Þ½3ðx2  xÞ2 ðy2  y þ z2  zÞ þ ð1  6x þ 6x2 Þðy2  yÞðz2  zÞ 3p2 þ sinðpxÞ sinðpyÞ sinðpzÞ þ p2 sinðpxÞ sinðpyÞ sinðpzÞ 1þk  p2 cosðpxÞ cosðpyÞ sinðpzÞ  p2 cosðpxÞ sinðpyÞ cosðpzÞ; f2 ¼ 200lð2x  1Þð2z  1Þ½3ðy2  yÞ2 ðx2  x þ z2  zÞ þ ð1  6y þ 6y2 Þðx2  xÞðz2  zÞ 3p2 þ sinðpxÞ sinðpyÞ sinðpzÞ  p2 cosðpxÞ cosðpyÞ sinðpzÞ 1þk þ p2 sinðpxÞ sinðpyÞ sinðpzÞ  p2 sinðpxÞ cosðpyÞ cosðpzÞ; f3 ¼ 200lð2x  1Þð2y  1Þ½3ðz2  zÞ2 ðx2  x þ y2  yÞ þ ð1  6z þ 6z2 Þðx2  xÞðy2  yÞ 3p2 þ sinðpxÞ sinðpyÞ sinðpzÞ  p2 cosðpxÞ sinðpyÞ cosðpzÞ 1þk  p2 sinðpxÞ cosðpyÞ cosðpzÞ þ p2 sinðpxÞ sinðpyÞ sinðpzÞ:

5797

L.-C. Xiao et al. / Applied Mathematics and Computation 217 (2011) 5790–5797 Table 1 Numerical results of 75-freedom FEM for k = 1. h

k~ u ~ uh k0;X

a

k~ u ~ uh kh

b

0.5 0.25 0.125 0.1

0.03854021 0.00393362 0.00040486 0.00020012

n 3.2924 3.2804 3.1577

0.50864544 0.12266192 0.02721298 0.01691425

n 2.0520 2.1723 2.1311

Table 2 Numerical results of 75-freedom FEM for k = 108. h

k~ u ~ uh k0;X

a

k~ u ~ u h kh

b

0.5 0.25 0.125 0.1

0.01287767 0.00156254 0.00017444 0.00008853

n 3.0429 3.1631 3.0395

0.23559395 0.05248456 0.01187056 0.00749438

n 2.1663 2.1445 2.0610

Then the exact solution ~ u ¼ ðu1 ; u2 ; u3 Þ of (2) is:

1 sinðpxÞ sinðpyÞ sinðpzÞ; 1þk 1 sinðpxÞ sinðpyÞ sinðpzÞ; u2 ¼ 100ðy  y2 Þ2 ð2x3  3x2 þ xÞð2z3  3z2 þ zÞ þ 1þk 1 sinðpxÞ sinðpyÞ sinðpzÞ: u3 ¼ 100ðz  z2 Þ2 ð2y3  3y2 þ yÞð2x3  3x2 þ xÞ þ 1þk

u1 ¼ 200ðx  x2 Þ2 ð2y3  3y2 þ yÞð2z3  3z2 þ zÞ þ

The domain X is divided into cubes uniformly. Let h be the edge length of each element. We assume a and b be the convergence order of k~ u ~ uh k0;X and k~ u ~ uh kh , respectively. The numerical results are presented in Tables 1 and 2. Acknowledgements This work was supported by the NSF of China (No. 11071226), the NSF of the Education Department of Henan Province (No. 2010A110005), and the Doctoral Foundation of Henan University of Technology (No. 2008BS013). References [1] M. Vogelius, An analysis of the p-version of the finite element method for nearly incompressible materials. Uniformly valid, optimal error estimates, Numer. Math. 41 (1983) 39–53. [2] I. Babuska, M. Suri, Locking effects in the finite element approximation of elasticity problems, Numer. Math. 62 (1992) 439–463. [3] D.N. Arnold, J. Douglas, C.P. Gupta, A family of higher order mixed finite element methods for plane elasticity, Numer. Math. 45 (1984) 1–22. [4] D.N. Arnold, R. Winther, Mixed finite elements for elasticity, Numer. Math. 92 (2002) 401–419. [5] S.C. Brenner, L.Y. Sung, Linear Finite element methods for planar linear elasticity, Math. Comput. 220 (1992) 321–330. [6] R. Stenberg, A family of mixed finite elements for the elasticity problem, Numer. Math. 53 (1988) 513–538. [7] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Spring-Verlag, 1991. [8] R.S. Falk, Nonconforming finite element methods for the equations of linear elasticity, Math. Comput. 57 (1991) 529–550. [9] L.H. Wang, H. Qi, On locking-free finite element schemes for the pure displacement boundary value problem in the planar elasticity, Mathematica Numerica Sinica 24 (2002) 243–256. [10] H. Qi, L.H. Wang, W.Y. Zheng, On locking-free finite element schemes for three-dimension elasticity, J. Comput. Math. 23 (2005) 101–112. [11] P. Ming, Nonconforming finite element method VS locking problems, Doctor thesis, The Academy of Mathematics and System Sciences, CAS, 1991. [12] P.G. Ciarlet, Mathematical Elasticity, Vol. I, Three-dimensional elasticity, North-Holland, Amsterdam, 1988. [13] S.C. Brenner, L.R. Scott, The Mathematical Theory of Finite Element Methods, Springer-Verlag, New York, 1994. [14] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, New York, 1978. [15] Z.C. Shi, The F-E-M-test for convergence of nonconforming finite elements, Math. Comput. 49 (1991) 391–405.