Mathematical Problems in Engineering

Volume 2013 (2013), Article ID 491626, 14 pages

http://dx.doi.org/10.1155/2013/491626

Research Article

Shape-Free Finite Element Method: Another Way between Mesh and Mesh-Free Methods

Song Cen,1,2,3 Ming-Jue Zhou,¹ and Yan Shang¹

¹Department of Engineering Mechanics, School of Aerospace, Tsinghua University, Beijing 100084, China

²High Performance Computing Center, School of Aerospace, Tsinghua University, Beijing 100084, China

³Key Laboratory of Applied Mechanics, School of Aerospace, Tsinghua University, Beijing 100084, China

Received 25 June 2013; Accepted 18 July 2013

Academic Editor: Zhiqiang Hu

Copyright © 2013 Song Cen et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Performances of the conventional finite elements are closely related to the mesh quality. Once distorted elements are used, the accuracy of the numerical results may be very poor, or even the calculations have to stop due to various numerical problems. Recently, the author and his colleagues developed two kinds of finite element methods, named hybrid stress-function (HSF) and improved unsymmetric methods, respectively. The resulting plane element models possess excellent precision in both regular and severely distorted meshes and even perform very well under the situations in which other elements cannot work. So, they are called shape-free finite elements since their performances are independent to element shapes. These methods may open new ways for developing novel high-performance finite elements. Here, the thoughts, theories, and formulae of above shape-free finite element methods were introduced, and the possibilities and difficulties for further developments were also discussed.

1. Introduction

As the cornerstone of computational mechanics and mathematics, the finite element method (FEM) has been considered as one of the greatest academic achievements in the last century [1, 2]. Over the past 60 years, with the progress in techniques of computers, the FEM has obtained remarkable developments in theories and applications, and became one of the main tools for computations and numerical simulations of science and engineering. However, it should be pointed out that, the FEM is a kind of numerical piecewise interpolation methods, so that its performance often depends on the mesh quality. For example, the accuracy of most conventional finite elements may drop dramatically once the mesh is distorted [3, 4], and this phenomenon will lead to incorrect results or interruption of computations due to numerical difficulties. Since the mesh distortions are unavoidable for complex geometry and large deformation problems, the sensitivity problem to mesh distortion has been regarded as one of the severest defects existing in the FEM.

In order to circumvent above trouble caused by mesh, numerous researchers began to develop so-called meshless or mesh-free approaches to replace the FEM, such as the element-free Galerkin method [5], the reproducing kernel particle method [6], the meshless local Petrov-Galerkin (MLPG) and the local boundary integral equation (LBIE) methods [7], the boundary node method [8], the hybrid boundary node method [9], the least-squares collocation meshless method [10], the PU-based meshless Shepard interpolation method [11]. These methods can produce excellent results without meshing (only distributed points are needed for most cases), so that they are free of many numerical problems aroused by the FEM mesh. However, since higher order interpolations and more complicated techniques are often used, these methods usually need much more computation costs, which is not acceptable for applications of large-scale engineering problems, especially for nonlinear problems. Therefore, at present, besides several special problems, the meshless or mesh-free methods have not completely taken the place of the FEM yet.

On the other hand, during the history of the FEM itself, numerous efforts have been also made for improving performance and robustness of the traditional finite element models, such as the hybrid stress method proposed by Pian et al. [12–14], the incompatible displacement modes proposed by Wilson et al. [15] and Taylor et al. [16], the enhanced strain method proposed by Simo and Rifai [17], the stabilization method proposed by Belytschko and Bacharch [18], the selectively reduced integration scheme proposed by Hughes [19], the assumed strain formulations proposed by MacNeal [20] and Piltner and Taylor [21], the quasi-conforming element method proposed by Tang et al. [22], the generalized conforming method proposed by Long and Huang [23], the Alpha finite element method (αFEM) [24] and the smoothed finite element method (S-FEM) [25, 26] proposed by Liu et al., the new spline finite element method [27] proposed by Chen et al., the FE-mesh-free element method proposed by Rajendran et al. [28, 29], the new natural coordinate methods proposed by Long et al. [30–37], and the Hamilton hybrid tress element method proposed by Cen et al. [38]. These works made great contributions on the finite element method. However, the sensitivity problem to mesh distortion has never been overcome from the outset. For example, once a quadrilateral element almost degenerates into a triangle or becomes a concave quadrangle, the accuracy will be lost greatly, or even the computation has to be stopped due to the numerical difficulties.

Recently, two mesh distortion immune techniques were proposed by the author and his coworkers. The first one is the hybrid stress-function (HSF) element method [39–42], which can be treated as the improved version for the initial hybrid stress element method proposed by Pian [12]. It starts from the principle of minimum complementary energy and employs the fundamental analytical solutions of the Airy stress function as the trial functions (analytical trial function method). The second one is the improved unsymmetric element method [43], which can solve the interpolation failure and rotational frame dependence problems existing in the original unsymmetric element method proposed by Rajendran et al. [44–46] and can produce more accurate results. It is based on the virtual work principle and used two different strain matrices in the final formulae: one is from the conventional isoparametric element, and the other is derived from the fundamental analytical solutions of elasticity.

The element models constructed by above methods possess excellent precision in both regular and severely distorted meshes, and even perform very well under the situations in which other elements cannot work (e.g., a quadrilateral element degenerates into triangular or concave quadrangular shapes). So, they are called shape-free finite elements since their performances are independent to the element shapes. It is very interesting that effective ways for completely avoiding sensitivity problem to mesh distortion may be found.

In the following sections, the thoughts, theories, and formulae of the shape-free finite element methods will be introduced, and the possibilities and difficulties for further developments were also discussed.

2. Fundamental Analytical Solutions for Trial Functions

The fundamental analytical solutions of the mechanics play important roles in both the hybrid stress-function (HSF) and the improved unsymmetric element methods. For plane problem without body forces (or with only constant body forces), the analytical solutions of the Airy stress function should satisfy the Beltrami-Michell equation (the compatibility equation expressed in terms of stress function). For the isotropic case, the equation is written as

and for the anisotropic case, the equation becomes

where are the reduced elastic compliances and have been defined in [41].

The fundamental analytical solutions can be used as the trial function in finite element method or other numerical methods. In order to choose appropriate number of the solutions for the construction of new element model, two principles must be followed: (i) the fundamental analytical solutions should be selected in turn from the lowest-order to higher-order; and (ii) the resulting stress or displacement fields should possess completeness in Cartesian coordinates. In plane elasticity, there are 3 solutions corresponding to the rigid-body displacement state, 3 for the constant stress state (linear displacement state), and 4 for each other higher order stress or displacement state. As an example, the first 18 solutions    and resulting stresses and displacements for isotropic case are given in Table 1.

Table 1: Eighteen fundamental analytical solutions of the Airy stress function and resulting stress and displacement solutions for plane stress problems*.

i

1
1
0
0
0
1/E
0

2
x
0
0
0
0
1/E

3
y
0
0
0
y/E
−x/E

4
x 2
0
2
0
−2 x/E
2y/E

5
xy
0
0
−1
− (1 + )y/E
−(1 + )x/E

6
y 2
2
0
0
2x/E
−2 y/E

7
x 3
0
6x
0
−3(x ² + y ²)/E
6xy/E

8
x 2y
0
2y
−2x
−2 xy/E
[y ² − (2 + )x ²]/E

9
xy 2
2x
0
−2y

−2 xy/E

10
y 3
6y
0
0
6xy/E
−3(y ² + x ²)/E

11
x 3y
0
6xy
−3x ²
−(3 x ²y + y ³)/E

12
xy 3
6xy
0
−3y ²
[3x ²y − (2 + )y ³]/E
−(3 x ²y + y ³)/E

13

−12y ²
12x ²
0
−4(3xy ²+ x ³)/E
4(3x ²y + y ³)/E

14

12()
−12()
−24xy
4(1 + )()/E
4(1 + )()/E

15

2x()
6xy ²
−2y()
[(2 + )x ²y ² + (3 + 2 )y ⁴]/2E
[2(1 + 2 ) x ³y]/E

16

10x ³
−10x()
−30x ²y
2.5[(1 + 2 ) x ²y ² − y ⁴]/E
10xy[]/E

17

6x ²y
−2y()
2x()
[2(1 + 2 )x ³ xy ³]/E
[(2 + )x ²y ² + ( )x ⁴]/2E

18

10y()
10y ³
−30xy ²
10xy[x ² − (2 + )y ²]/E
2.5[(1 + 2 )  x ²y ² − x ⁴]/E

Note: and for plane strain problem. E and are Young’s modulus and Poisson’s ratio, respectively.

3. Shape-Free FEM I: Plane Hybrid Stress-Function (HSF) Element

3.1. The Construction Procedure for 2D Hybrid Stress-Function Elements

For a plane finite element model, its complementary energy functional can be written in the following matrix form [12, 39–42]:

where is the complementary energy within the element; is the complementary energy along the kinematic boundaries (here, all element boundaries are treated as the kinematic boundaries because the boundary displacements will be prescribed in (6); is the elasticity matrix of compliances, it possesses different forms for isotropic and anisotropic cases, or for plane strain and plane stress states; is the thickness of the element; , the element stress vector; , the traction force vector along the element boundaries; , the displacement vector along element boundaries, which can be interpolated by the element nodal displacement vector :

where matrix is the interpolation function matrix for element boundary displacements.

The stress vector in (5) can be derivable from the Airy stress function , that is,

And the traction force vector T in (5) can be written as

where and are the direction cosines of the outer normal of the element boundaries.

Substituting (7) and (8) into (3) yields

Let

with

where is the number of the fundamental analytical solutions used for stress function in (10); are fundamental analytical solutions (in Cartesian coordinates) of the Airy stress function , which must be selected from complete second-order terms (, , ) to complete higher-order terms (see Table 1); are unknown constants. Obviously, such trial functions will directly lead to more reasonable stress fields satisfying both equilibrium and compatibility conditions.

Substitution of (6) and (10) into (9) yields

with

where , and () are the corresponding stresses of (see Table 1), and all are expressed in Cartesian coordinates.

According to the principle of minimum complementary energy, we require

Then, the unknown constant vector can be expressed in terms of the nodal displacement vector :

Substitution of (16) into the in (12) yields

From the viewpoint of the element definition given in [12], matrix mentioned earlier can be considered as the stiffness matrix of the hybrid stress-function element, and therefore, it can readily be incorporated into the standard finite element program framework.

Once the element nodal displacement vector is solved, the element stresses can be given by

3.2. A Typical HSF Element HSF-Q8-15β [40]

Up to present, several plane HSF element models for both isotropic and anisotropic materials, including 8-node and 12-node quadrilateral elements [39–41], have been successfully developed. As an example, the 8-node HSF element using 15β, denoted by HSF-Q8-15β [40], is introduced here.

Consider an 8-node quadrilateral element shown in Figure 1, and any edge of the element can be either straight or curved, and (, ) are the usual isoparametric coordinates. Differing from the usual models, the element shapes can be either convex or concave. The element nodal displacement vector is given by

(a)

(b)

(c)

(d)

(e)

Figure 1: 8-node quadrilateral plane element with arbitrary shapes.

Assume that the displacements along each element edge are interpolated by the nodal displacements of each edge. Therefore, , and corresponding matrix (see (6) and (12)) of each element edge can be given as follows:

with

in which () are the shape functions of the plane 8-node isoparametric serendipity elements

where (, ) are the isoparametric coordinates of node .

Fifteen analytical solutions () for the stress function (see (10)), which have been listed in Table 1, are taken as the trial functions. That is to say, fifteen unknown constants are introduced. Then, the matrix in (14) is a matrix. It can be seen that the stress fields possess third-order completeness in both and . The resulting element model is denoted as HSF-Q8-15β.

In order to evaluate the matrices and by Gaussian numerical integration procedure, the Cartesian coordinates and within an element should be expressed in terms of local coordinates (isoparametric coordinates). Let

where () () are the Cartesian coordinates of the node .

Example 1 . (pure bending for a cantilever beam (Figure 2)). As shown in Figure 2, a cantilever beam under plane stress condition is subjected to a constant bending moment . This problem was earlier used by Lee and Bathe [4] to assess the distortion sensitivity of the quadratic plane elements. The theoretical solutions for this problem are given by [47]:

Six regular and distorted mesh divisions are employed (see Figure 2), in which meshes 4, 5, and 6 are so severely distorted that some quadrilateral elements degenerate into triangles or concave quadrangles. Numerical results, obtained by present element HSF-Q8-15β, 8-node isoparametric serendipity element Q8, 9-node isoparametric Lagrangian element Q9L, are listed in Table 2. It can be seen that exact solutions can always be obtained by HS-F elements, no matter the meshes are distorted or not, and no matter the element shapes are convex or concave.

Table 2: Results at selected locations for the pure bending problem (Figure 2).

Q8
Q9L
HSF-Q8-15
Exact

Mesh 1
()
120.000
120.000
120.000
120.0

()
−120.000
−120.000
−120.000
−120.0

v() × 10³
−12.000
−12.000
−12.000
−12.0

Mesh 2
()
56.447
120.000
120.000
120.0

()
−74.863
−120.000
−120.000
−120.0

v() × 10³
−2.328
−12.000
−12.000
−12.0

Mesh 3
()
13.665
120.000
120.000
120.0

(
5.262
120.000
120.000
120.0

(
−5.665
−120.000
−120.000
−120.0

()
−14.299
−120.000
−120.000
−120.0

v() × 10³
−0.477
−12.000
−12.000
−12.0

Mesh 4
()
—
—
120.000
120.0

()
—
—
−120.000
−120.0

v() × 10³
—
—
−12.000
−12.0

Mesh 5
()
—
—
120.000
120.0

()
—
—
−120.000
−120.0

v() × 10³
—
—
−12.000
−12.0

Mesh 6
()
—
—
120.000
120.0

()
—
—
120.000
120.0

()
—
—
−120.000
−120.0

()
—
—
−120.000
−120.0

v() × 10³
—
—
−12.000
−12.0

Figure 2: Pure bending problem for a cantilever beam.

Actually, element HSF-Q8-15β can produce the exact solutions for constant stress/strain problems, no matter the element edges are straight or curved, and no matter the shapes of the elements are convex or concave. And so long as all element edges keep straight, element HSF-Q8-15β can also produce the exact solutions for linear stress/strain problems, no matter the shapes of the elements are convex or concave. For higher-order problems, HSF-Q8-15β elements possess much better convergence than the corresponding isoparametric elements.

3.3. A 4-Node HSF Quadrilateral Plane Membrane Element with Drilling Degrees of Freedom [42]

Consider a 4-node quadrilateral element with drilling degrees of freedom shown in Figure 3. Differing from the usual displacement-based models, the element shape is allowed to be either convex or concave. The element nodal displacement vector is defined as

Figure 3: A 4-node HSF quadrilateral plane membrane element with drilling degrees of freedom.

where are not the physical rotations of the element nodes. Instead, the definitions of the drilling degrees of freedom given by Allman [48] are employed. So, the element boundary displacements can be written as

For matrix S defined in (14), let , that is, the first seven analytical solutions for the stress function , which have been listed in Table 1, are taken as the trial functions. Then, the matrix S is a 3 × 7 matrix. It can be seen that the stress fields possess linear completeness in both and . The resulting element model is denoted as HSF-Q4θ-7β.

Example 2 . (a wedge subjected to a uniformly distributed load (Figure 4)). As shown in Figure 4, a cantilever wedge is subjected to a uniformly distributed load . Because of its triangular shape, the wedge can not be modeled without the use of triangular and/or distorted quadrilateral elements. Since the present quadrilateral element HSF-Q4θ-7β can still perform when its shape degenerates into triangle, it can therefore be readily used to model such wedge problem. Numerical results and the percentage errors of the radial stresses at selected points are listed in Table 3. The present element HSF-Q4θ-7β performs very well for such high-order bending problem.

Other examples can be found in [42]. Actually, element HSF-Q4θ-7β may be the best model among all the 4-node quadrilateral plane membrane elements with drilling degree of freedom.

Table 3: HSF-Q4-7 results of radial stress at selected points for a wedge subjected to a uniformly distributed load (Figure 4).

Mesh
1 × 6
2 × 12
4 × 24
Exact [48]

at point A (0, 5)
7.6835

(1.38%)
7.5894 (0.13%)
7.5806 (0.02%)
7.5792

at point B (1, 5)
−7.7920 (1.47%)
−7.7088 (0.39%)
−7.6832 (0.05%)
−7.6792

Figure 4: A wedge subjected to a uniformly distributed load.

4. Shape-Free FEM II: Improved Unsymmetric Element Method [43]

4.1. Brief Reviews on Original 8-Node Unsymmetric Element US-QUAD8 [44]

Based on the virtual work principle [44], the final element stiffness matrix (unsymmetric) can be written as [44–46]

where is the Jacobian determinant; D is the elasticity matrix; is the element strain matrix for conventional 8-node plane isoparametric element Q8:

in which have been given by (22); is given by

where are Cartesian shape functions, given by

in which (, ) are the Cartesian coordinates of node .

The resulting element model is denoted as US-QUAD8. It can keep producing exact solutions for constant and linear stress/strain problems by using both regular and severely distorted meshes. However, once the element is distorted from quadrilateral to a certain shape, such as a triangle, the first matrix on the left hand of (30) will be singular, that is, interpolation failure will occur. Furthermore, this element exhibits rotational frame dependence in higher-order problems since the displacement fields constructed by are not complete.

4.2. The Shape-Free 8-Node Unsymmetric Element US-ATFQ8 [43]

In [43], the above unsymmetric element method was improved by analytical trial function method. The element displacement fields in terms of Cartesian coordinates are assumed as

in which

where and are the displacement solutions derived from the fundamental analytical solutions of the Airy stress function, and have been given in Table 1; are 14 unknown parameters related to and ;    are 4 unknown parameters related to and ;

are assumed displacements which possess third-order completeness in and ; and

are assumed fourth-order displacement fields.

For the fourth-order displacement fields and , two additional relaxed point bubble conditions are introduced

So, and can be expressed in terms of and as follows

in which

Let

Then, the displacement fields in (31) can be rewritten as

with

Substitution of coordinates of eight nodes into (39) yields

in which

where is the element nodal displacement vector and has been given by (19). Then, can be solved by

It should be noted that so long as the element nodes are not in coincidence with each other, the matrix is hardly singular for various element shapes.

The displacement fields in terms of Cartesian coordinates of the element US-QUAD8 can be replaced by (39). Thus, the element stiffness matrix in (27) can be rewritten as

where and are given by (28), and

where , and have been given in Table 1. And the element stresses can be evaluated by

The resulting element model is denoted as US-ATFQ8. It can avoid the interpolation failure and rotational frame dependence problems existing in the original unsymmetric element US-QUAD8, and possesses much better accuracy for higher-order problems.

Example 3 . (linear bending for a cantilever beam (Figure 5)). As shown in Figure 5, a cantilever beam under plane stress condition is subjected to a linear bending moment caused by a shear force at the free end. This problem is often used to assess the distortion sensitivity of the cubic plane elements [4]. The theoretical solutions for this problem are given by [47]. The meshes used for this example are also given in Figure 5. Table 4 lists results using full or suggested integration schemes. It is astonishing that the present element US-ATFQ8 can almost produce the exact solutions under all meshes. This phenomenon has not been found before for any other 8-node plane elements.

The element US-ATFQ8 can nearly produce the exact solutions for constant, pure bending and linear bending problems, no matter the element edges are straight or curved, and no matter the shapes of the elements are convex or concave. This important advantage cannot be achieved by other 8-node, 16-DOF elements. Compared with the original unsymmetric 8-node plane element US-QUAD8, the present element possesses obvious advantages: (i) element US-ATFQ8 can still work well when interpolation failure modes for US-QUAD8 occur; (ii) element US-ATFQ8 naturally eliminates the rotation dependency which exists in element US-QUAD8 for higher-order problems; and (iii) element US-ATFQ8 can provide more accurate results than those obtained by element US-QUAD8 in most complicated problems. Therefore, element US-ATFQ8 is a truly shape-free finite element model.

Table 4: Results at selected locations for the linear bending problem (Figure 5).

Q8  (3 × 3)
HSF-Q8-15  (5 × 5)
US-QUAD8  (3 × 3)
US-ATFQ8  

(5 × 5 or 3 × 3)
Exact

()
−66.871
−120.00
−66.871
−120.000
−120.0

Mesh 1
()
66.871
120.00
66.871
120.000
120.0

v() × 10³
6.282
8.558
6.282
8.0460
8.046

()
−35.248
−113.63
−74.317
−120.000
−120.0

Mesh 2
()
67.040
119.63
71.781
120.000
120.0

v() × 10³
1.611
7.104
6.415
8.0460
8.046

()
−11.307
−68.645
−66.349
−120.000
−120.0

()
−55.699
−66.980
−67.974
−120.000
−120.0

Mesh 3
()
55.305
66.146
67.995
120.000
120.0

()
11.073
69.328
66.328
120.000
120.0

v() × 10³
4.520
6.350
6.292
8.0460
8.046

()

−72.471

−120.00
−120.0

Mesh 4
()
Failed
76.265
Failed
120.00
120.0

v() × 10³

6.663

8.0460
8.046

()

−77.515
−58.158
−120.000
−120.0

Mesh 5
()
Failed
66.306
59.158
120.000
120.0

v() × 10³

6.729
5.471
8.0460
8.046

()

−70.346
−42.671
−120.000
−120.0

()

−66.835
−44.636
−120.000
−120.0

Mesh 6
()
Failed
66.877
40.974
120.000
120.0

()

68.512
41.048
120.000
120.0

v(  × 10³

6.373
3.752
8.0460
8.046

()
−105.69
−120.00
−105.69
−120.000
−120.0

Mesh 7
()
105.69
120.00
105.69
120.000
120.0

v() × 10⁴
3.485
3.975
3.485
3.6600
3.66

()
−115.25
−119.02
−118.99
−120.000
−120.0

Mesh 8
()
115.25
119.02
118.99
120.000
120.0

v() × 10⁴
3.3423
3.796
3.5570
3.6600
3.66

()
−118.68
−120.00
−118.68
−120.000
−120.0

Mesh 9
()
118.68
120.00
118.68
120.000
120.0

v() × 10⁵
11.1963
15.750
11.1963
13.200
13.2

()
−97.066
−121.61
−118.94
−120.00 b/−119.65a
−120.0

Mesh 10
()
138.21
124.61
117.74
120.00 b/119.36a
120.0

v() × 10⁵
12.8162
13.738
12.9157
13.200 b/13.190a
13.2

()

−24.732
−36.104
−120.000
−120.0

Mesh 11
()
Failed
20.048
36.230
120.000
120.0

v() × 10³

1.517
3.1470
8.0460
8.046

3 × 3 integration scheme; b5 × 5 integration scheme.

Figure 5: Linear bending problem for a cantilever beam.

5. Discussions and Concluding Remarks

Since there is no inverse of the Jacobian matrix existing in the expressions of the element stiffness matrix, most sensitivity problems to mesh distortion can be automatically avoided. And the use of the fundamental analytical solutions of elasticity significantly improved the element accuracy. Have we already found effective ways for developing new finite elements whose performances are independent to element shapes? If it was true, the users of FEM would not worry about the mesh quality in practical analyses, and it was not necessary for researchers to make more efforts on mesh generation techniques.

However, the answer is not optimistic. We have to encounter many additional difficulties when generalize above methods for developing other models.

When one develops 3D 20-node hexahedral element, in order to construct a complete fourth-order displacement fields (similar with (31) for 2D case), 75 fundamental analytical solutions must be considered, which is too complicated and not economical for practical applications. On the other hand, when one wants to develop lower-order element models (such as plane 4-node quadrilateral element with 8 DOFs, or 3D 8-node hexahedral element with 24 DOFs) directly using the above method, the number of the effective fundamental analytical solutions will be limited by the number of the element DOFs, so that there is no noteworthy improvements on accuracy.

In HSF element method, in order to guarantee its robust performance, the assumed exact displacements of the element boundaries may be required. However, it is difficult to obtain such displacements in 3D problems.

Furthermore, how to select the appropriate forms of the fundamental analytical solutions [49], how to construct shape-free plate and shell elements, how to improve the computational efficiency of the elements with unsymmetric stiffness matrices, how to deal with the problems in which there are no fundamental analytical solutions (such as nonlinear problems), and so forth, are all important issues that we have to face to.

Obviously, new techniques must be developed for solving above problems. But at the same time, do not forget that some existing achievements, such as the new natural coordinate methods [30–38], generalized conforming method [2, 23], assumed strain method [17], could play important roles in the constructions of the new shape-free finite elements, especially for low-order high-performance models.

So, for establish systematic shape-free finite element method, there is still a long way to go.

Acknowledgments

The authors would like to acknowledge the financial support of the National Natural Science Foundation of China (11272181), the Specialized Research Fund for the Doctoral Program of Higher Education of China (20120002110080), and the National Basic Research Program of China (Project no. 2010CB832701).

References

O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method for Solid and Structural Mechanics, Elsevier, Oxford, UK, 6th edition, 2005.

Y.-Q. Long, S. Cen, and Z.-F. Long, Advanced Finite Element Method in Structural Engineering, Tsinghua University Press, Beijing, China; Springer, Berlin, Germany, 2009.

J. A. Stricklin, W. S. Ho, E. Q. Richardson, and W. E. Haister, “On isoparametric vs linear strain triangular elements,” International Journal for Numerical Methods in Engineering, vol. 11, pp. 1041–1043, 1977.

N.-S. Lee and K.-J. Bathe, “Effects of element distortions on the performance of isoparametric elements,” International Journal for Numerical Methods in Engineering, vol. 36, no. 20, pp. 3553–3576, 1993.

T. Belytschko, Y. Y. Lu, and L. Gu, “Element-free Galerkin methods,” International Journal for Numerical Methods in Engineering, vol. 37, no. 2, pp. 229–256, 1994.

W. K. Liu, S. Jun, and Y. F. Zhang, “Reproducing kernel particle methods,” International Journal for Numerical Methods in Fluids, vol. 20, no. 8-9, pp. 1081–1106, 1995.

S. N. Atluri, H.-G. Kim, and J. Y. Cho, “A critical assessment of the truly Meshless Local Petrov-Galerkin (MLPG), and Local Boundary Integral Equation (LBIE) methods,” Computational Mechanics, vol. 24, no. 5, pp. 348–372, 1999.

Y. X. Mukherjee and S. Mukherjee, “The boundary node method for potential problems,” International Journal for Numerical Methods in Engineering, vol. 40, no. 5, pp. 797–815, 1997.

J.-M. Zhang, Z.-H. Yao, and H. Li, “A hybrid boundary node method,” International Journal for Numerical Methods in Engineering, vol. 53, no. 4, pp. 751–763, 2002.

X. Zhang, X.-H. Liu, K.-Z. Song, and M.-W. Lu, “Least-squares collocation meshless method,” International Journal for Numerical Methods in Engineering, vol. 51, no. 9, pp. 1089–1100, 2001.

Y. C. Cai and H. H. Zhu, “A PU-based meshless Shepard interpolation method satisfying delta property,” Engineering Analysis with Boundary Elements, vol. 34, no. 1, pp. 9–16, 2010.

T. H. H. Pian, “Derivation of element stiffness matrices by assumed stress distributions,” AIAA Journal, vol. 2, no. 7, pp. 1333–1336, 1964.

T. H. H. Pian and D. P. Chen, “Alternative ways for formulation of hybrid stress elements,” International Journal for Numerical Methods in Engineering, vol. V 18, no. 11, pp. 1679–1684, 1982.

T. H. H. Pian and K. Sumihara, “Rational approach for assumed stress finite elements,” International Journal for Numerical Methods in Engineering, vol. 20, no. 9, pp. 1685–1695, 1984.

E. L. Wilson, R. L. Taylor, W. P. Doherty, and T. Ghabussi, “Incompatible displacement models,” in Numerical Computational Methods in Structural Mechanics, S. T. Fenven, Ed., pp. 43–57, Academic Press, New York, NY, USA, 1973.

R. L. Taylor, P. J. Beresford, and E. L. Wilson, “A non-conforming element for stress analysis,” International Journal for Numerical Methods in Engineering, vol. 10, pp. 1211–1219, 1976.

J. C. Simo and M. S. Rifai, “A class of mixed assumed strain methods and the method of incompatible modes,” International Journal for Numerical Methods in Engineering, vol. 29, no. 8, pp. 1595–1638, 1990.

T. Belytschko and W. E. Bachrach, “Efficient implementation of quadrilaterals with high coarse-mesh accuracy,” Computer Methods in Applied Mechanics and Engineering, vol. 54, no. 3, pp. 279–301, 1986.

T. J. R. Hughes, “Generalization of selective integration procedures to anisotropic and nonlinear media,” International Journal for Numerical Methods in Engineering, vol. 15, no. 9, pp. 1413–1418, 1980.

R. H. MacNeal, “Derivation of element stiffness matrices by assumed strain distributions,” Nuclear Engineering and Design, vol. 70, no. 1, pp. 3–12, 1982.

R. Piltner and R. L. Taylor, “A systematic construction of B-bar functions for linear and non-linear mixed-enhanced finite elements for plane elasticity problems,” International Journal for Numerical Methods in Engineering, vol. 44, no. 5, pp. 615–639, 1999.

L.-M. Tang, W.-J. Chen, and Y.-X. Liu, “Formulation of Quasi-Conforming element and Hu-Washizu principle,” Computers & Structures, vol. 19, no. 1-2, pp. 247–250, 1984.

Y.-Q. Long and M.-F. Huang, “A generalized conforming isoparametric element,” Applied Mathematics and Mechanics, vol. 9, no. 10, pp. 929–936, 1988.

G. R. Liu, T. Nguyen-Thoi, and K. Y. Lam, “A novel FEM by scaling the gradient of strains with factor α (αFEM),” Computational Mechanics, vol. 43, no. 3, pp. 369–391, 2009.

G. R. Liu, K. Y. Dai, and T. T. Nguyen, “A smoothed finite element method for mechanics problems,” Computational Mechanics, vol. 39, no. 6, pp. 859–877, 2007.

G. R. Liu, H. Nguyen-Xuan, and T. Nguyen-Thoi, “A theoretical study on the smoothed FEM (S-FEM) models: properties, accuracy and convergence rates,” International Journal for Numerical Methods in Engineering, vol. 84, no. 10, pp. 1222–1256, 2010.

J. Chen, C.-J. Li, and W.-J. Chen, “A family of spline finite elements,” Computers & Structures, vol. 88, no. 11-12, pp. 718–727, 2010.

S. Rajendran and B. R. Zhang, “A “FE-meshfree” QUAD4 element based on partition of unity,” Computer Methods in Applied Mechanics and Engineering, vol. 197, no. 1–4, pp. 128–147, 2007.

J. P. Xu and S. Rajendran, “A partition-of-unity based “FE-Meshfree” QUAD4 element with radial-polynomial basis functions for static analyses,” Computer Methods in Applied Mechanics and Engineering, vol. 200, no. 47-48, pp. 3309–3323, 2011.

Y.-Q. Long, J.-X. Li, Z.-F. Long, and S. Cen, “Area co-ordinates used in quadrilateral elements,” Communications in Numerical Methods in Engineering, vol. 15, no. 8, pp. 533–545, 1999.

Z.-F. Long, J.-X. Li, S. Cen, and Y.-Q. Long, “Some basic formulae for area co-ordinates in quadrilateral elements,” Communications in Numerical Methods in Engineering, vol. 15, no. 12, pp. 841–852, 1999.

X.-M. Chen, S. Cen, X.-R. Fu, and Y.-Q. Long, “A new quadrilateral area coodinate method (QACM-II) for developing quadrilateral finite element models,” International Journal for Numerical Methods in Engineering, vol. 73, no. 13, pp. 1911–1941, 2008.

Z.-F. Long, S. Cen, L. Wang, X.-R. Fu, and Y.-Q. Long, “The third form of the quadrilateral area coordinate method (QACM-III): theory, application, and scheme of composite coordinate interpolation,” Finite Elements in Analysis and Design, vol. 46, no. 10, pp. 805–818, 2010.

H.-G. Li, S. Cen, and Z.-Z. Cen, “Hexahedral volume coordinate method (HVCM) and improvements on 3D Wilson hexahedral element,” Computer Methods in Applied Mechanics and Engineering, vol. 197, no. 51-52, pp. 4531–4548, 2008.

X.-M. Chen, S. Cen, Y.-Q. Long, and Z.-H. Yao, “Membrane elements insensitive to distortion using the quadrilateral area coordinate method,” Computers & Structures, vol. 82, no. 1, pp. 35–54, 2004.

S. Cen, X.-M. Chen, and X.-R. Fu, “Quadrilateral membrane element family formulated by the quadrilateral area coordinate method,” Computer Methods in Applied Mechanics and Engineering, vol. 196, no. 41–44, pp. 4337–4353, 2007.

A.-K. Soh, Y.-Q. Long, and S. Cen, “Development of eight-node quadrilateral membrane elements using the area coordinates method,” Computational Mechanics, vol. 25, no. 4, pp. 376–384, 2000.

S. Cen, T. Zhang, C.-F. Li, X.-R. Fu, and Y.-Q. Long, “A hybrid-stress element based on Hamilton principle,” Acta Mechanica Sinica, vol. 26, no. 4, pp. 625–634, 2010.

X.-R. Fu, S. Cen, C. F. Li, and X.-M. Chen, “Analytical trial function method for development of new 8-node plane element based on the variational principle containing airy stress function,” Engineering Computations, vol. 27, no. 4, pp. 442–463, 2010.

S. Cen, X.-R. Fu, and M.-J. Zhou, “8- and 12-node plane hybrid stress-function elements immune to severely distorted mesh containing elements with concave shapes,” Computer Methods in Applied Mechanics and Engineering, vol. 200, no. 29–32, pp. 2321–2336, 2011.

S. Cen, X.-R. Fu, G.-H. Zhou, M.-J. Zhou, and C. F. Li, “Shape-free finite element method: the plane hybrid stress-function (HS-F) element method for anisotropic materials,” Science China Physics, Mechanics & Astronomy, vol. 54, no. 4, pp. 653–665, 2011.

S. Cen, M.-J. Zhou, and X.-R. Fu, “A 4-node hybrid stress-function (HS-F) plane element with drilling degrees of freedom less sensitive to severe mesh distortions,” Computers & Structures, vol. 89, no. 5-6, pp. 517–528, 2011.

S. Cen, G.-H. Zhou, and X.-R. Fu, “A shape-free 8-node plane element unsymmetric analytical trial function method,” International Journal for Numerical Methods in Engineering, vol. 91, no. 2, pp. 158–185, 2012.

S. Rajendran and K. M. Liew, “A novel unsymmetric 8-node plane element immune to mesh distortion under a quadratic displacement field,” International Journal for Numerical Methods in Engineering, vol. 58, no. 11, pp. 1713–1748, 2003.

S. Rajendran, “A technique to develop mesh-distortion immune finite elements,” Computer Methods in Applied Mechanics and Engineering, vol. 199, no. 17-20, pp. 1044–1063, 2010.

E. T. Ooi, S. Rajendran, and J. H. Yeo, “Remedies to rotational frame dependence and interpolation failure of US-QUAD8 element,” Communications in Numerical Methods in Engineering, vol. 24, no. 11, pp. 1203–1217, 2008.

S. P. Timoshenko and J. N. Goodier, Theory of Elasticity, McGraw-Hill, New York, NY, USA, 3rd edition, 1934.

D. J. Allman, “A compatible triangular element including vertex rotations for plane elasticity analysis,” Computers & Structures, vol. 19, no. 1-2, pp. 1–8, 1984.

X.-R. Fu, M.-W. Yuan, S. Cen, and G. Tian, “Characteristic equation solution strategy for deriving fundamental analytical solutions of 3D isotropic elasticity,” Applied Mathematics and Mechanics, vol. 33, no. 10, pp. 1253–1264, 2012.