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. [1214], 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. [3037], 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 [3942], 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. [4446] 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.

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, 3942]: 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 [3941], 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

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.

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

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.

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 [4446] 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.

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 [3038], 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).