The use of 2D enriched elements with bubble functions for finite element analysis

Computers and Structures 84 (2006) 2081–2091 www.elsevier.com/locate/compstruc

The use of 2D enriched elements with bubble functions for finite element analysis Shi-Pin Ho *, Yen-Liang Yeh Department of Mechanical Engineering, National Cheng Kung University, Tainan, Taiwan, ROC Received 11 January 2005; accepted 11 April 2006 Available online 2 October 2006

Abstract In this paper, the concept that adds the interior nodes of the Lagrange elements to the serendipity elements is described and a family of enriched elements is presented to improve the accuracy of finite element analysis. By the use of the static condensation technique at the element level, the extra computation time in using these elements can be ignored. Plane stress problems are used as examples in this paper. The numerical results show that these enriched elements are more accurate than the traditional serendipity elements. The convergence rate of the proposed elements is the same as the traditional serendipity elements. The error norm of the second and third order proposed elements can be reduced from 40% to 60% when compared with the use of the traditional serendipity elements. In the numerical examples, the use of the second and third order proposed elements not only give an improvement in element accuracy but also save computation time, when the precondition conjugate gradient method is used to solve the system of equations. The saving of computation time is due to the decrease of iteration number in iteration method.  2006 Elsevier Ltd. All rights reserved. Keywords: Finite element; Enriched element; Bubble function; Serendipity element; Lagrange element; Static condensation

1. Introduction The selection and use of elements is very important in finite element analysis, affecting both accuracy and computation time. Finite element methodology has been employing Lagrange elements and serendipity elements for three decades. However, research for improvements in computational speed and accuracy has continued. In the 1970’s, Herrmann [1] presented a two-dimensional incompatible element and evaluated its efficiency. In the same decade Wilson et al. [2] presented the incompatible models method, aiming to improve the behavior of elements in bending situations. Wilson’s method was later improved by Taylor et al. [3]. Chen and Cheung [4] applied Taylor’s method to solve the axisymmetric solid problem. Taiebat *

Corresponding author. Tel.: +886 6 2757575x62164; fax: +886 6 2352973. E-mail address: [email protected] (S.-P. Ho). 0045-7949/$ - see front matter  2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2006.04.008

and Carter [5] extended the idea to produce three dimensional non-conform elements and Alves de Sousa et al. [6] recommended low order elements for three-dimensional analysis. The use and improvement of higher order elements is another way to improve the accuracy of finite element analysis. Rathod and Sridevi [7,8] derived the interpolation polynomials for the general serendipity elements which allow arbitrarily placed nodes along edges. Celia and Gray [9] showed that side nodes positioned at the same relative distance from corner nodes in both local and global space and interior nodes positioned such that they do not affect the Jacobian of the coordinate transformation provide improved accuracy over nodes positioned locally without regard for the global configuration. Universal serendipity elements [10,11] defined as isoparametric elements at element edges in an arbitrary manner is another application. Kikuchi et al. [12] developed a general principle to derive eight-node elements by eliminating the interior node of

2082

S.-P. Ho, Y.-L. Yeh / Computers and Structures 84 (2006) 2081–2091

a nine-node element with all the Cartesian quadratic polynomials preserved. Adding bubble functions can also greatly improve element performance. A bubble function is defined as a function that disappears along the element boundaries. Accordingly, the inter-element compatibility requirements for displacement are maintained for the element, and the bubble function parameters can be eliminated at the element level by static condensation [13] and can be recovered for subsequent calculations. Bubble functions have been introduced to construct plate element models [14– 17]. For example, Pinsky and Jasti [14] developed fournode and nine-node plate bending elements which incorporated bubble function displacements. Kemp et al. [16] used the bubble function displacements to develop a four-node solid shell element based on the assumed strain formulation and verified that the bubble function displacement alleviates sensitivity to mesh distortion as well as element locking. On the other hand, bubble functions are employed to solve advection–diffusion problems by Brezzi, Franca and Farhat [18–21] and are shown to help in stabilizing the advective operator without recourse to upwinding or any other numerical strategy. Furthermore, the limitation of bubble functions is discussed by Franca and Farhat [22] and error analysis of residual-free bubbles is discussed by Brezzi et al. [23] and Sanglli [24]. Brezzi and Marini [25] introduced a possible way of mimicking the effect of residual-free bubbles just by constructing a suitable subgrid. They proposed sufficient conditions on the subgrid discretization for convectiondomain problem in order to obtain error estimates of the same quality as one could obtain by solving the fine-level equations in an exact way. However, the introduction of the higher order bubble functions may result in element models that require a higher-order numerical integration rule for the element stiffness matrix. Therefore accuracy and computation time must be considered when we evaluate elements. Presently, the second order serendipity element is the most popular choice in industrial applications. The goal of this paper is to derive a family of enriched elements that improve the accuracy of serendipity elements and do not increase computation time. To do so, the interior nodes of Lagrange elements are added to serendipity elements and then the polynomial order of the shape functions is increased to improve accuracy. In the finite element calculation process, static condensation is used to reduce the degree of freedom at interior nodes. Since static condensation is performed at the element level, only a small increase in computation time is needed. The structure of the paper is as follows. The shape functions of the enriched quadrilateral elements are derived in Section 2. Static condensation of the enriched elements is discussed in Section 3. In Section 4, the performance of the proposed elements is compared on the solid problems. Conclusions are provided in Section 5.

2. Formulation of the enriched elements In this paper, a nodal based basis is constructed by adding the interior nodes of Lagrange basis to a serendipity basis. The serendipity basis is then corrected such that the Kronecker delta property is satisfied at the interior nodes, and a family of enriched elements is presented. ðLÞ We use the notation N i to indicate the shape function which corresponds to the ith node of the Lagrange element. After the numbers of all edge nodes are assigned, the interior node is assigned in turn. On the other hand, serendipity elements have no interior nodes. In this methodology ðSÞ the notation N i is used to indicate the shape function which corresponds to ith node of the serendipity element. For the case of an nth order enriched element, the interior nodes of the mth order Lagrange element are added to the nth order serendipity elements. This method can provide higher order shape functions for the interior of lower order serendipity elements. We need to keep the relation Ni(nj, gj) = dij for every shape function within the element, so the serendipity shape functions corresponding to the edge nodes must be modified. The enriched element, designated nSmL, is derived by combining the nth order serendipity element and the interior nodes of mth order Lagrange element (m P 2). The nth order serendipity elements have 4n nodes and the shape functions are repreðSÞ sented as N i ; i ¼ 1 to 4n. The mth order Lagrange elements have (m  1)2 interior nodes. Since the 4m nodes on the edge are numbered first, the (m  1)2 interior nodes of the mth order Lagrange element are numbered from 4m + 1 to (m + 1)2. The shape functions of the Lagrange ðLÞ element interior nodes, N i ; i ¼ 4m þ 1 to ðm þ 1Þ2 , are used as the interior shape functions of the nSmL element. Since the edges of the nSmL element have 4n nodes, these shape functions are numbered from 4n + 1 to 4n + (m  1)2. Therefore, the shape functions of the nSmL element can be represented as

N i ðn; gÞ ¼

8 ðLÞ > > < N ð4m4nþiÞ ðSÞ > > : Ni 

4nþðm1Þ P

ðSÞ

2

for i ¼ 4n þ 1 to 4n þ ðm  1Þ ; 2

aij N j

for i ¼ 1 to 4n

j¼4nþ1

ðLÞ

where aij ¼ N j ðni ; gi Þ. Because N j ; j ¼ 4m þ 1 to 2 ðm þ 1Þ is zero at all edge nodes, the modification of aijNj makes Ni satisfy the Ni(nj,gj) = dij condition. For example, the shape functions of 1S2L element are 1 1 N 1 ðn; gÞ ¼ ð1  nÞð1  gÞ  ð1  n2 Þð1  g2 Þ 4 4 1 1 N 2 ðn; gÞ ¼ ð1 þ nÞð1  gÞ  ð1  n2 Þð1  g2 Þ 4 4 1 1 N 3 ðn; gÞ ¼ ð1 þ nÞð1 þ gÞ  ð1  n2 Þð1  g2 Þ 4 4 1 1 N 4 ðn; gÞ ¼ ð1  nÞð1 þ gÞ  ð1  n2 Þð1  g2 Þ 4 4 N 5 ðn; gÞ ¼ ð1  n2 Þð1  g2 Þ

S.-P. Ho, Y.-L. Yeh / Computers and Structures 84 (2006) 2081–2091

"

#

(



It will be seen that the use of N5(n, g) = (1  n2)(1  g2) increases the accuracy of the element. Because the shape function of the 1S2L element interior node vanishes on all edge sides, the 1S2L element is conforming.

where

3. Static condensation

K ee ¼ K ee  K ei K 1 ii K ie

0

K ie

K ii

ue ui

¼

Fe

)

Fi

F e ¼ F e  K ei K 1 ii F i

The degrees of freedom of an nSmL element can be classified into two groups. Those that are not connected to the freedoms of another element are named ‘‘interior freedoms’’ and these that are connected to at least one other element are named ‘‘edge freedoms’’. When the 1S2L element is used in the finite element model, the element coefficient matrix Ke and the element loading vector Fe are 38 e 9 8 e 9 2 e F > u > > K 11 K e12 K e13 K e14 K e15 > > > > > > 1e > > 1e > > > > 6 K e K e K e K e K e 7> > > > > 6 21 22 23 11 25 7> = > = < u2 > 7 6 e e e e e e e 6 K 31 K 32 K 33 K 34 K 35 7 u3 ¼ F 3 7> > > > 6 > > > > e> 7> e > 6 e u4 > F4> > > 4 K 41 K e42 K e43 K e44 K e45 5> > > > > > > > > ; ; : : K e51 K e52 K e53 K e54 K e55 ue5 F e5 Because the interior freedoms are not coupled with the other elements, the entries of the row and column that correspond to the interior freedoms are not used in the algorithm for the following assembly procedure. When the K e55 is used as the pivot for the Gaussian elimination method to solve the system of equations, the modified coefficient matrix and the loading vector are 38 e 9 8 e 9 2 e K 11 K e12 K e13 K e14 0 > >F1> > > > u1 > > > > > > > > 7> 6 e e> e e e > > > u > > 0 7> F e2 > 6 K 21 K 22 K 23 K 11 2> = = < < 7 6 e e e e 7 ue3 ¼ F e 6K K K K 0 32 33 34 7> > > 3 > 6 31 > > 7> e > 6 e > u4 > > > > > 4 K 41 K e42 K e43 K e44 0 5> F e4 > > > > > > > > ; : e; : e> e e e e e u5 K 51 K 52 K 53 K 54 K 55 F5 Ke

K ee

2083

Ke

and use them to assemble the global coefficient matrix and force. After ue is solved, we get ui as 1 ui ¼ K 1 ii F i  K ii K ie ue

Therefore we can see that Gaussian elimination is used to produce the K ee and F e , then store the parts which are 1 needed to solve ui, K 1 ii F i and K ii K ie . In the coefficient matrix K ee ¼ K ee  K ei K 1 ii K ie it can be seen that Kee applies partial factorization and this procedure can be regarded as a preconditioner [26]. Notably, this preconditioner is applied directly to the entities of the coefficient matrix, so the procedure for solving the system of equations can be applied to another preconditioner if the precondition conjugate gradient method is used. Therefore, in the following numerical experiment, the enriched elements, static condensation and precondition conjugate gradient method cooperate in the finite element analysis procedure. 4. Numerical examples Plane stress problems are considered in this section to verify the efficiency and accuracy of the enriched elements with bubble functions. All the computational work is carried out on a PC running an AMD Duron 800 MHz with 384MB RAM. These problems are solved using first order, second order and third order quadrilateral elements. Isotropic material is used in the numerical example and the stress-strain relationship is

where K eij ¼ K eij  K eii K e5j and F ei ¼ F ei  K eii F e5 . 55 55 Therefore we use the interior freedoms to do the Gaussian elimination procedure called static condensation at the element level. After the edge freedoms are assembled and solved, we get the interior freedom by

r ¼ Ce 2 C 11 6 C ¼ 4 C 12 0

ue5 ¼ ðF e5  K e51 ue1  K e52 ue2  K e53 ue3  K e54 ue4 Þ=K e55

C 11 ¼ C 22

In general, all the interior freedoms can be made invisible to the user by applying a small dose of static condensation. The element coefficient matrix can be reduced at the element level and divided into several partitions      K ee K ei ue Fe ¼ K ie K ii ui Fi where ue is the undetermined variables vector of the edge nodes and ui is the undetermined variables vector of the interior nodes. Kee, Kei, Kie and Kii are submatrices. Because the interior nodes have no relation with other elements, we can let

C 12

0

3

7 C 22 0 5 0 C 66 E mE ¼ ; C 12 ¼ ; 1  m2 1  m2

C 66 ¼ G

where E is the elastic modulus of the material, m is the Poisson ratio of the material and the G is the shear modulus of the material. Stress vector r includes the rx, ry and sxy components. Strain vector e includes the ex, ey and cxy components. We will next compare the performance of both the traditional elements used in finite element analysis (herein designated 1L, 2S, 2L, 3S and 3L) and the proposed elements (1S2L, 1S3L, 2S3L, 2S4L, 3S2L, and 3S4L). Diagrams of the proposed elements are shown in Table 1. The precondition conjugate gradient method (PCG) is used to solve the system of equations. After comparing

2084

S.-P. Ho, Y.-L. Yeh / Computers and Structures 84 (2006) 2081–2091

Table 1 The elements used in the numerical experiment

Serendipity

First order element

Second order element

Third order element

1S

2S

3S

2L

3L

1S2L

2S3L

3S2L

1S3L

2S4L

3S4L

Lagrange

Enriched

different preconditioner, the incomplete Cholesky factorization preconditioner [27] is used in the following numerical examples.

formula and the 1S3L element needs at least the 9-point formula. The 2S element only needs the 4-point formula but The 2S4L element needs the 16-point formula. In the third order elements, the 3L and the 3S4L elements need at least the 16-point formula. However, the deflection does not vary significantly with a different Gauss–Legendre quadrature formula if the Gauss–Legendre quadrature formula can be successfully used for the targeted element. Therefore, we should select the lower limit formula in order to save computation time. The formulae that are used in the later numerical experiment are marked by the symbol * in the Table 2. Table 3 shows the deflection w of point A and the PCG iteration number NPCG when this example is solved for regular mesh and different element sizes h (0.2, 0.1, 0.05 and 0.025). The reference deflection 9.8547E09 is given by the model of the third order Lagrange element and element size 0.01. The numerical results show that the deflections of the enriched elements are more accurate than the comparable serendipity elements and the PCG iteration numbers are reduced by the static condensation. 4.2. Plane stress problem which devises an exact solution

4.1. Clamped plane with applied linear pressure The first example is a clamped plane with applied linear pressure and is shown in Fig. 1. The deflection of point A is compared in the following numerical experimentation. We consider which formula of the Gauss–Legendre quadrature is applicable when using the proposed element. We use this numerical example in which element size h is 0.2 to check the deflection when adopting different Gauss–Legendre quadrature formulae (herein included 1-point, 4-point, 9-point, 16-point and 25-point formulae). The numerical results are presented in Table 2 and the symbol X means that the Gauss–Legendre quadrature formula cannot be used for the targeted element. From the Table 2, the 1S and 1S2L elements need at least the 4-point

A plane stress problem which devises exact solution is taken as an example. A small deformation field is assumed for the region bounded by 0 6 x 6 1 and 0 6 y 6 1. This example is solved for regular mesh and different element sizes h (0.2, 0.1, 0.05 and 0.025). The exact solutions of the displacement field are sine functions of the Cartesian coordinates (x, y), shown as u ¼ 0:001  sinðpxÞ  sinðpyÞ v ¼ 0:001  sinðpxÞ  sinðpyÞ where u and v are displacement functions in the x and y directions. These displacement fields satisfy the condition that the Jacobian is positive in the field.

Fig. 1. The clamped plane with applied linear pressure.

S.-P. Ho, Y.-L. Yeh / Computers and Structures 84 (2006) 2081–2091

2085

Table 2 The deflection of point A vs. Gauss–Legendre quadrature formulae 1-Point

4-Point

9-Point

16-Point

25-Point

1S 1S2L 1S3L

X X X

8.627E09* 9.259E09* X

8.627E09 9.146E09 9.169E09*

8.627E09 9.146E09 9.162E09

8.627E09 9.146E09 9.162E09

2S 2L 2S3L 2S4L

X X X X

9.766E09* X X X

9.750E09 9.794E09* 9.829E09* X

9.750E09 9.794E09 9.824E09 9.829E09*

9.750E09 9.794E09 9.824E09 9.828E09

3S 3L 3S2L 3S4L

X X X X

X X X X

9.779E09* X 9.814E09* X

9.777E09 9.837E09* 9.813E09 9.846E09*

9.777E09 9.837E09 9.813E09 9.844E09

(Unit: meter.)

Table 3 The deflection and PCG iteration number vs. various element size h

0.2

0.1

0.05

0.025

w

NPCG

w

NPCG

w

1S 1S2L 1S3L

8.627E09 9.259E09 9.169E09

10 11 11

9.476E09 9.664E09 9.639E09

16 16 16

9.745E09 9.796E09 9.789E09

27 28 28

9.822E09 9.836E09 9.834E09

51 52 52

2S 2L 2S3L 2S4L

9.766E09 9.802E09 9.829E09 9.829E09

25 24 19 19

9.827E09 9.840E09 9.845E09 9.846E09

52 40 32 32

9.845E09 9.850E09 9.852E09 9.852E09

114 75 60 60

9.851E09 9.853E09 9.854E09 9.854E09

245 147 121 121

3S 3L 3S2L 3S4L

9.779E09 9.830E09 9.814E09 9.846E09

36 29 29 20

9.828E09 9.846E09 9.840E09 9.852E09

67 51 49 35

9.845E09 9.852E09 9.849E09 9.854E09

138 96 93 68

9.851E09 9.854E09 9.853E09 9.855E09

284 194 183 134

NPCG

w

NPCG

(Unit: meter.)

When the deformation functions are substituted into the equations of stress-strain and strain-deformation, the stress components are

The results show solution error as error L2 norm kekL2 and error energy norm kekeng with the following definitions [28]:

rx ¼ 0:001  p  ½C 12  cosðpyÞ  sinðpxÞ þ C 11  cosðpxÞ  sinðpyÞ

Error energy norm: kekeng ¼

sxy ¼ 0:001  p  ½C 66  sinðpðx þ yÞÞ When the stress components are substituted into the equilibrium equations, we get the body force functions as fx ¼ 0:001  p  ½C 11  sinðpxÞ  sinðpyÞ  C 12  cosðpxÞ  cosðpyÞ  C 66  cosðpðx þ yÞÞ fy ¼ 0:001  p  ½C 22  sinðpxÞ  sinðpyÞ  C 12  cosðpxÞ  cosðpyÞ  C 66  cosðpðx þ yÞÞ

Thus, this example of body force applied to a fixed boundary has been solved by the proposed nSmL finite element method.

ðe  eh ÞT ðr  rh Þdxdy

12

X

ry ¼ 0:001  p  ½C 22  cosðpyÞ  sinðpxÞ þ C 21  cosðpxÞ  sinðpyÞ

Z Z

Error L2 norm: kekL2 ¼

Z Z

ju  uh j2 dxdy

12

X

where u is the exact solution vector which includes ux, uy and uz components. r is the exact stress vector which includes rx, ry and sxy components. e is the exact strain vector which includes ex, ey and cxy components. uh is the finite element approximate solution vector, rh is the finite element approximate stress vector and eh is the finite element approximate strain vector. The error norms of the proposed elements are compared with the traditional serendipity and Lagrange elements. The ratio of the reduced error norm is defined such that the difference in the error norm of the elements is divided by the error norm of the serendipity elements or

2086

S.-P. Ho, Y.-L. Yeh / Computers and Structures 84 (2006) 2081–2091

the Lagrange elements. In the finite element analysis procedure, we must deal with producing the element coefficient matrix, assembling the global coefficient matrix and solving the system of equations. Therefore, computation times are separated into three groups in the following tables; matrix formation time, equation solution time and total computation time. The iteration numbers of the precondition conjugate gradient method are listed in the parenthesis. Computation time for static condensation is included in matrix formation time. Because the computation time for calculating the interior freedoms was much less than the others, computation time for the calculation of interior freedoms is also included in matrix formation time. 4.2.1. Error and computation time of first order quadrilateral elements The error energy norm of the first order quadrilateral elements are shown in Table 4. From Table 4, the 1S2L element, with one extra node in the center, reduces the error energy norm about 54%, while the 1S3L element, with four extra nodes, reduces the energy norm about 61% when compared with the 1S element. The error L2 norm in Table 5 shows that 1S2L and 1S3L elements reduce the error L2 norm about 78%. The three computation times for the first order quadrilateral elements are shown in Table 6. From the table we can see the equation solution time for the 1S2L and 1S3L elements is quite similar to that of the 1S element. The total computation time of the 1S2L element increases only 9% when compared with the 1S element. However, the matrix formation time of the 1S3L element, with four additional interior nodes, is much more than the 1S element because each additional interior freedom must undergo static condensation and the 1S3L element needs 9-point Gauss– Legendre quadrature formulation.

Table 6 The computation time of the first order quadrilateral elements h (m)

1S

1S2L

1S3L

Matrix formation time 0.2 0.01 0.1 0.01 0.05 0.04 0.025 0.08

<0.01 0.01 0.03 0.12

<0.01 0.03 0.04 0.42

Equation solution time 0.2 <0.01 (6) 0.1 <0.01 (10) 0.05 0.04 (18) 0.025 0.39 (35)

<0.01 <0.01 0.03 0.38

Total computation time 0.2 0.01 0.1 0.01 0.05 0.06 0.025 0.47

<0.01 0.01 0.07 0.51

(6) (10) (18) (34)

<0.01 (6) <0.01 (10) 0.04 (18) 0.371 (34)

<0.01 0.03 0.15 0.801

(Unit: second.)

4.2.2. Error and computation time of second order quadrilateral elements The error energy norm of the second order quadrilateral elements is shown in Table 7 and the error L2 norm of the second order quadrilateral elements is shown in Table 8. The error norm of the 2L element is very close to the error norm of the 2S element. From Table 7, the 2S3L element, with four extra nodes, reduces the error energy norm about 46%, while the 2S4L element, with nine extra nodes, reduces the energy norm about 51% when compared with the 2S element. Table 8 shows that the 2S3L element reduces the error L2 norm about 54% and 2S4L element reduces the error L2 norm about 57%. Table 9 shows the three computation times for the second order quadrilateral elements. The equation solution time for the 2S3L and 2S4L elements is smaller than that for the 2S element. The equation solution time is smaller

Table 4 Error energy norm of the first order quadrilateral elements

Table 7 Error energy norm of the second order quadrilateral elements

h (m)

1S

1S2L

1S3L

h (m)

0.2 0.1 0.05 0.025

2.187E+02 1.096E+02 5.486E+01 2.743E+01

1.030E+02 4.968E+01 2.461E+01 1.228E+01

8.656E+01 4.249E+01 2.115E+01 1.056E+01

0.2 0.1 0.05 0.025

(Unit:

pffiffiffiffiffiffiffiffiffiffi N m.)

Table 5 Error L2 norm of the first order quadrilateral elements

2S

1.818E+01 4.472E+00 1.114E+00 2.783E01 pffiffiffiffiffiffiffiffiffiffi (Unit: N m.)

2L

2S3L

2S4L

1.785E+01 4.453E+00 1.113E+00 2.782E01

9.369E+00 2.397E+00 6.029E01 1.510E01

8.364E+00 2.165E+00 5.462E01 1.369E01

Table 8 Error L2 norm of the second order quadrilateral elements

h (m)

1S

1S2L

1S3L

h (m)

2S

2L

2S3L

2S4L

0.2 0.1 0.05 0.025

2.800E05 7.076E06 1.774E06 4.439E07

5.345E06 1.206E06 2.934E07 7.286E08

6.542E06 1.568E06 3.880E07 9.675E08

0.2 0.1 0.05 0.025

1.392E06 1.772E07 2.225E08 2.785E09

1.382E06 1.769E07 2.225E08 2.785E09

5.940E07 8.057E08 1.029E08 1.292E09

5.651E07 7.499E08 9.526E09 1.196E09

(Unit: meter.)

(Unit: meter.)

S.-P. Ho, Y.-L. Yeh / Computers and Structures 84 (2006) 2081–2091

2087

Table 9 The computation time of the second order quadrilateral elements

Table 12 The computation time of the third order quadrilateral elements

h (m)

h (m)

2S

2L

Matrix formation time 0.2 0.01 0.1 0.02 0.05 0.06 0.025 0.25 Equation solution time 0.2 <0.01 (21) 0.1 0.08 (40) 0.05 0.951 (82) 0.025 8.202 (168)

0.01 0.02 0.08 0.3

<0.01 (12) 0.05 (20) 0.45 (38) 3.825 (78)

Total computation time 0.2 0.01 0.1 0.1 0.05 1.011 0.025 8.452

0.01 0.07 0.53 4.125

2S3L

2S4L

0.02 0.06 0.21 0.841

0.05 0.16 0.631 2.514

<0.01 (11) 0.05 (21) 0.461 (41) 4.166 (81)

<0.01 (11) 0.04 (21) 0.43 (40) 3.966 (80)

0.02 0.11 0.681 5.027

0.05 0.21 1.081 6.53

3S

3L

Matrix formation time 0.2 0.01 0.1 0.06 0.05 0.22 0.025 0.881

0.02 0.09 0.36 1.442

Equation solution time 0.2 0.02 (25) 0.1 0.28 (48) 0.05 2.233 (96) 0.025 19.278 (191)

0.02 (13) 0.13 (23) 1.062 (45) 9.003 (88)

Total computation time 0.2 0.03 0.1 0.34 0.05 2.453 0.025 20.159

0.04 0.23 1.432 10.475

3S2L 0.02 0.07 0.25 1.012

0.02 (19) 0.19 (33) 1.482 (63) 12.588 (124)

0.04 0.26 1.732 13.61

3S4L 0.06 0.24 0.941 3.766

0.01 (12) 0.13 (23) 1.042 (44) 8.923 (88)

0.07 0.38 2.003 12.759

(Unit: second.)

(Unit: second.)

for the second order enriched elements relative to the 2S element since the iteration number is reduced. However, the matrix formation time of the 2S3L and 2S4L element is slightly greater than the 2S element, so the 2S3L and 2S4L total computation times are slightly less than the serendipity elements.

node in the center, is as same as the traditional Lagrange 3L element. The three computation times of the third order quadrilateral elements are shown in Table 12. It is seen that the equation solution times of the enriched elements are less than the serendipity elements because the iteration number of the precondition conjugate gradient method is greatly reduced. The equation solution time for the 3S4L element is about 53% smaller than the 3S element. Therefore the total computation time of the third quadrilateral elements decreases significantly as the number of interior nodes increases. Use of the enriched third order quadrilateral elements not only improves accuracy but also saves total computation time.

4.2.3. Error and computation time of third order quadrilateral elements The error energy norm and error L2 norm of the third order quadrilateral elements are shown in Tables 10 and 11. From Tables 10 and 11, the 3S2L element reduces the error energy norm and error L2 norm about 65% when compared with the 3S element. The 3S4L element reduces the error energy norm and error L2 norm significantly when compared with the 3S element and the 3L element. The error norm of the 3S2L element, with only one extra

Table 10 Error energy norm of the third order quadrilateral elements h (m)

3S

0.2 0.1 0.05 0.025

3.533E+00 4.265E01 5.208E02 7.706E03 pffiffiffiffiffiffiffiffiffiffi (Unit: N m.)

3L

3S2L

3S4L

9.416E01 1.207E01 1.220E02 4.422E03

1.021E+00 1.233E01 1.230E02 4.423E03

1.101E01 7.111E03 6.992E04 5.847E04

Table 11 Error L2 norm of the third order quadrilateral elements h(m)

3S

3L

3S2L

3S4L

0.2 0.1 0.05 0.025

1.743E07 1.016E08 6.164E10 4.460E11

4.061E08 2.636E09 1.337E10 2.425E11

4.594E08 2.718E09 1.352E10 2.426E11

2.561E09 7.242E11 9.091E12 4.371E12

(Unit: meter.)

4.2.4. Summary From above numerical testing, we find that the 1S2L element is better than the 1S element and the 2S3L element is better than the 2L element, showing an obvious improvement of approximately 50% when compared with the traditional Lagrange element. Although the 1S3L and 2S4L elements improved about 60%, they still require more interior nodes and cost more computation time than the 1S2L and 2S3L elements. It is also seen that the proposed 2S3L and 3S4L elements provide significantly more accurate when compared with the traditional comparable Lagrange elements. From the error norm of the first order elements, the 1S element, the 1S2L element and the 1S3L element show the same convergence rate. The 2S, 2L, 2S3L and 2S4L elements have similar convergence rates. Similarly, the 3S, 3L, 3S2L and 3S4L elements have equivalent convergence rates. As above, we found that the iteration number of the conjugate gradient method is decreased by static condensation. This relation can be expected to extend to higher orders, a useful relation in terms of saving equation

2088

S.-P. Ho, Y.-L. Yeh / Computers and Structures 84 (2006) 2081–2091 2.50

log(||e||)

2.00 1S

1.50

1S2L

1.00

1S3L

0.50 0.00

0

0.2

0.4

0.6

0.8

Fig. 2(a) shows that the 1S2L element and 1S3L element own similar trend. Fig. 2(b) and (c) shows that the 2S3L element is the best one in the cases of the second order elements and the 3S4L is the best one in the cases of the third order elements, respectively. The dependence on error L2 norm vs. computation time shows the same trend as the dependence on error energy norm vs. computation time.

1

time

4.3. The plane stress analysis of the plate with a hole

1.50 1.00

log(||e||)

2S

0.50

2L 2S3L

0.00

0

2

4

6

8

10

15

20

25

2S4L

-0.50 -1.00

time 1.00 0.00

log(||e||)

0

5

10

-1.00

3S 3L 3S2L

-2.00

3S4L

-3.00 -4.00

time

Fig. 2. The dependence ‘‘error energy norm vs. computation time’’ for various elements. (a) First order elements, (b) second order elements and (c) third order elements.

solution time. As the degrees of freedom increase, solving the system of equations uses progressively more of the computation time. Therefore increase in static condensation computation time can be ignored when the equation solution time is much greater than matrix formation time. Fig. 2 shows the dependence on error energy norm vs. computation time for the various elements. The error energy norms are took the logarithmic operation in order to show the trend clearly. In Fig. 2, all enriched elements shown a better trend between accuracy and computation time than the one of traditional serendipity and Lagrange elements.

A simple plane stress example, which is a plate to be applied loading Fx = 1000Nt, is considered in this section to verify the performance of the proposed elements. Since the problem is symmetric, the geometry is simplified to quarter and the symmetric boundary conditions are applied. The finite element mesh and the boundary condition are shown in Fig. 3. The total number of elements is 234. The material properties are set as E = 2 · 1011 and m = 0.3. The computation time and the x direction deformation of point A are shown in Table 13. A reference deformation 3.1439E09 is pointed out by ANSYS, using element type is plane82 and number of elements is 14,232 (total 43,293 nodes). The DOF means the total degree of freedom, including both interior freedoms and edge freedoms. From the numerical results, the solution procedure dominates the computation time on a complex geometry, distortion elements and large number of nodes. The proposed elements show the advantage on both the accuracy and efficient. Especially, the 2S3L and 3S4L elements can provide a more accurate solution and save computation time when comparing with the comparable serendipity elements. 4.4. Extend the concept to enrich the hexahedral elements The concept that adds the interior nodes of the Lagrange elements to the serendipity elements can also be extended to enrich the hexahedral elements. The nodes of hexahedral serendipity elements only exist at the edge and do not have any nodes on the surface and in the interior of element. In order to distinguish the quadrilateral enriched elements and the hexahedral enriched elements, the hexahedral enriched element which the interior nodes of mth order hexahedral Lagrange elements are added to

Fig. 3. The mesh and boundary conditions of the plate with a hole.

S.-P. Ho, Y.-L. Yeh / Computers and Structures 84 (2006) 2081–2091

2089

Table 13 The numerical results of the plate with a hole analysis DOF

Formation time

Solution time

Total time

Ux (m)

1S 1S2L 1S3L

552 1020 2424

0.01 0.02 0.06

0.07 (39) 0.07 (38) 0.07 (38)

0.08 0.09 0.13

2.8364E09 2.9113E09 2.9005E09

2S 2L 2S3L 2S4L

1570 2038 3442 5782

0.02 0.07 0.11 0.33

0.76 0.42 0.38 0.42

(126) (69) (64) (64)

0.78 0.49 0.49 0.75

3.1041E09 3.1270E09 3.1442E09 3.1478E09

3S 3L 3S2L 3S4L

2588 4460 3056 6800

0.11 0.31 0.13 0.41

1.99 0.94 1.32 0.90

(145) (74) (97) (71)

2.10 1.25 1.45 1.31

3.1142E09 3.1426E09 3.1351E09 3.1438E09

The reference deformation is 3.1439E09. (Time unit : second.)

the nth order hexahedral serendipity elements is designated Hex-nSmL. Table 14 shows the diagrams of the first and second order hexahedral enriched elements which are tested in following numerical experimentation. We will examine the performance of the hexahedral serendipity elements (herein designated Hex-1S and Hex-2S), hexahedral Lagrange elements (herein designated Hex-2L) and the hexahedral enriched elements (herein designated Hex1S2L, Hex-1S3L, Hex-2S2L and Hex-2S3L) used in finite element analysis. 4.4.1. Curved beam with a shear loading A curved cantilever beam loaded with an in-plane shear load is considered (Fig. 4) and is solved using the proposed elements. We also check which formulae of Gauss–Legendre quadrature should be used when using the proposed elements. The Gauss–Legendre quadrature formulae are frequently used in the finite element method because they require fewer base points to achieve the same accuracy.

Table 14 The first and second order hexahedral enriched elements

Serendipity

First order

Second order

Hex-1S

Hex-2S

Lagrange

Enriched

Hex-2L

Hex-1S2L

Hex-2S2L

Hex-1S3L

Hex-2S3L

Fig. 4. Curved beam: inner radius = 4.12; outer radius = 4.32; arc = 90; thickness = 0.1; E = 1.0 · 107; m = 0.25; mesh = 6 · 1 · 1. A unit-shear force F is applied at the tip.

Different Gauss–Legendre quadrature formulae (herein including 1-point, 8-point, 27-point and 64-point formulae) are adopted in the following numerical tests. The elements are tested for its principal deformation modes under curved geometry. Normalized tip deflections with respect to theoretical values are shown in Table 15. The symbol X in the following table means that the Gauss–Legendre quadrature formula cannot be used for the targeted element because a serious error occurred. The resultant values which are listed in the following table are the computational results to be divided by theoretical solution 0.08734 [29] and the number of iterations for precondition conjugate gradient method is listed in the parenthesis. In the numerical results, the Hex-2S3L shows a more accurate solution than the other second order element in the curved beam test when using the 27-point Gauss– Legendre quadrature formula. However, we find that the

Table 15 Normalized tip deflection ratio of curved beam

Hex-1S Hex-1S2L Hex-1S3L Hex-2S Hex-2L Hex-2S2L Hex-2S3L

1-Point

8-Point

27-Point

64-Point

X X X X X X X

0.07324 (27) 0.06389 (26) X X X X X

0.07324 0.06886 0.07373 0.87455 0.87782 0.87486 0.88207

0.07324 0.06683 0.07374 0.87468 0.87796 0.87499 0.87967

(27) (26) (27) (99) (113) (100) (98)

(27) (26) (27) (98) (112) (100) (97)

The above results are the computational tip deflections to be divided by the theoretical solution and the theoretical solution is 0.08734 [29].

2090

S.-P. Ho, Y.-L. Yeh / Computers and Structures 84 (2006) 2081–2091

for saving the equation solution time since the iteration number of the conjugate gradient method is decreased. 5. Conclusion

Fig. 5. Twisted beam: length = 12.0; width = 1.1; depth = 0.32; twist = 90 (root to tip); E = 29.0 · 106; m = 0.22; mesh = 12 · 2 · 2. A unit-shear force F is applied at the tip and the beam is fixed at the root as shown.

Hex-1S element shows its inability to represent the bending state and the Hex-1S2L and Hex-1S3L elements also present the same deficiency. When both accuracy and efficiency is considered, one should select the lowest and most usable formula in order to save computation time. The Hex-1S and Hex-1S2L elements should adopt the 8-point Gauss– Legendre quadrature formula, but the Hex-1S3L element should adopt the 27-point Gauss–Legendre quadrature formula. All the second order elements should adopt the 27-point Gauss–Legendre quadrature formula. 4.4.2. Twisted cantilever beam with a shear load A twisted cantilever beam loaded with a unit shear load is considered. Fig. 5 shows the mesh and boundary condition of the twisted cantilever beam. The enriched hexahedral elements are tested for its principle deformation modes under angular distortion. Table 16 shows that the computed tip deflections are normalized with respect to the theoretical solution. The resultant values which are listed in Table 16 are the computational results to be divided by theoretical solution 0.005424 [29] and the number of iterations for precondition conjugate gradient method is listed in the parenthesis. The second order enriched hexahedral elements also show the advantage Table 16 The numerical results of the twisted cantilever beam DOF

Formation time

Hex-1S Hex-1S2L Hex-1S3L

351 495 1503

0.02 0.02 0.09

Hex-2S Hex-2L Hex-2S2L Hex-2S3L

1143 1875 1287 2295

0.20 0.30 0.20 0.29

Solution time 0.17 (90) 0.17 (90) 0.17 (90) 7.52 21.16 7.17 6.75

(579) (872) (552) (505)

Total time

Deflection ratio

0.19 0.19 0.26

0.20257 0.20461 0.20463

7.71 21.46 7.37 7.04

0.89672 0.89720 0.89688 0.89700

The above results are the computational tip deflections to be divided by the theoretical solution and the theoretical solution is 0.005424 [29].

We have presented a family of enriched elements for finite element analysis. These elements are a combination of several of the characteristics of traditional Lagrange and serendipity elements. The proposed elements have demonstrated improved accuracy and computation time relative to their Lagrange and serendipity counterparts for at least certain cases. The same concept can also be extended to enrich the hexahedral elements. Specifically, the 2S3L and 3S4L elements are more accurate than their serendipity and Lagrange counterparts, while the 1S2L and 3S2L elements are more accurate than their serendipity counterparts. Numerical experiments showed that error norms are reduced from 40% to over 60% when compared with the serendipity elements. In addition, the convergence rate of the nSmL element is the same as the traditional serendipity element. Since the enriched elements, static condensation and precondition conjugate gradient method cooperate, the accuracy is improved and the computation time is saved in finite element analysis. By the use of the static condensation technique at the element level, extra computation time using these enriched elements can be ignored and the procedure can be regarded as a preconditioner. When the precondition conjugate gradient method is used to solve the system of equations, the iteration number can be reduced. Hence the computation time can be cut even though the dimension of coefficient matrix is as same as the serendipity element. Notably, use of the second and third order proposed elements not only improves accuracy but also reduces total computation time relative to third order serendipity elements. Therefore, the 2S3L element is recommended for the second quadrilateral element and the 3S4L element is recommended for the third quadrilateral element according to their high performance. References [1] Herrmann LR. Efficiency evaluation of a two-dimensional incompatible finite element. Comput Struct 1973;3:1377–95. [2] Wilson EL, Taylor RL, Doherty WP, Ghaboussi J. Incompatible displacement models. In: Fenves SJ et al., editors. Numerical and computer methods in structural mechanics. New York: Academic Press; 1973. [3] Taylor RL, Beresford PJ, Wilson EL. A non-conforming element for stress analysis. Int J Numer Methods Eng 1976;10:1211–9. [4] Chen W, Cheung YK. The nonconforming element method and refined hybrid element method for axisymmetric solid. Int J Numer Methods Eng 1995;39:2509–29. [5] Taiebat HH, Carter JP. Three-dimension non-conforming elements, Research Report R808, Department of Civil Engineering, The University of Sydney, 2001. [6] Alves de Sousa RJ, Natal Jorge RM, Areias PMA, Fontes Valente RA, Ce´sar de Sa´ JMA. Low order elements for 3D analysis. In: Mang

S.-P. Ho, Y.-L. Yeh / Computers and Structures 84 (2006) 2081–2091

[7]

[8]

[9]

[10] [11] [12] [13]

[14]

[15]

[16]

[17]

HA, Rammerstorfer FG, Eberhardsteiner J, editors. Proceedings of the fifth world congress on computational mechanics, Vienna, Austria; 2002. Rathod HT, Sridevi K. General complete Lagrange family for the cube in finite element interpolations. Comput Methods Appl Mech Eng 2000;181:295–344. Rathod HT, Sridevi K. General complete Lagrange interpolations with applications to three-dimensional finite element analysis. Comput Methods Appl Mech Eng 2001;190:3325–68. Celia MA, Gray WG. An improved isoparametric transformation for finite element analysis. Int J Numer Methods Eng 1984;20: 1443–59. Citipitioglu E. Universal serendipity elements. Int J Numer Methods Eng 1983;19:803–10. Utku M. An improved transformation for universal serendipity elements. Comput Struct 1999;73:199–206. Kikuchi F, Okabe M, Fujio H. Modification of the 8-node serendipity element. Comput Methods Appl Mech Eng 1999;179:91–109. Cook RD, Malkus DS, Plesha ME, Witt RJ. Concepts and applications of finite element analysis. New York: John Wiley and Sons; 2002. Pinsky PM, Jasti RV. A mixed finite element formulation for Reissner–Mindlin plates based on the use of bubble functions. Int J Numer Methods Eng 1989;28:1677–702. Hong WI, Kim YH, Lee SW. An assumed strain triangular solid shell element with bubble function displacements for analysis of plates and shells. Int J Numer Methods Eng 2001;52:455–69. Kemp BL, Cho C, Lee SW. A four-node solid shell element formulation with assumed strain. Int J Numer Methods Eng 1998; 43:909–24. Auricchio F, Taylor RL. A triangular thick plate finite element with an exact thin limit. Finite Elem Anal Design 1995;19:57–68.

2091

[18] Brezzi F, Bristeau MO, Franca LP, Mallet M, Roge G. A relationship between stabilized finite element methods and the Galerkin method with bubble functions. Comput Methods Appl Mech Eng 1992;96: 117–29. [19] Baiocchi C, Brezzi F, Franca LP. Virtual bubbles and Galerkin-leastsquares type methods (Ga.L.S). Comput Methods Appl Mech Eng 1993;105:125–41. [20] Franca LP, Farhat C. Bubble functions prompt unusual stabilized finite element methods. Comput Methods Appl Mech Eng 1995;123: 299–308. [21] Brezzi F, Russo A. Choosing bubbles for advection–diffusion problems. Math Models Methods Appl Sci 1994;4(4):571–87. [22] Franca L, Farhat C. On the limitation of bubble functions. Comput Methods Appl Mech Eng 1994;117:225–30. [23] Brezzi F, Hughes TJR, Marini LD, Russo A, Suli E. A priori error analysis of residual-free bubbles for advection–diffusion problems. SIAM J Numer Anal 1999;36(6):1933–48. [24] Sangalli G. Global and local error analysis for the residual-free bubbles method applied to advection-dominated problems. SIAM J Numer Anal 2000;38(5):1496–522. [25] Brezzi F, Marini LD. Augmented spaces, two-level methods, and stabilizing subgrids. Int J Numer Methods Eng 2002;40:31–46. [26] Farhat C, Sobh N. A coarse/fine preconditioner for very ill-conditioned finite element problems. Int J Numer Methods Eng 1989;28: 1715–23. [27] Golub GH, Van Loan CF. Matrix computations. London: Johns Hopkins; 1996. [28] Zienkiewicz OC, Taylor RL. The finite element method—volume 1. Basic formulation and linear problems. London: McGraw-Hill; 1989. [29] Macneal RH, Harder RL. A proposed standard set of problems to test finite element accuracy. Finite Elem Anal Design 1985;1:3–20.