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

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.


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 Mathematical Problems in Engineering Pian et al. [12][13][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 quasiconforming 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][31][32][33][34][35][36][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][40][41][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][45][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.

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   ( = 1 ∼ 18) and resulting stresses and displacements for isotropic case are given in Table 1.

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

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][40][41][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); C is the elasticity matrix of compliances, it possesses different forms for isotropic and anisotropic cases, or for plane strain and plane stress states; t is the ) and  → /(1 − ) for plane strain problem.E and  are Young's modulus and Poisson's ratio, respectively.thickness of the element; , the element stress vector; T, the traction force vector along the element boundaries; u, the displacement vector along element boundaries, which can be interpolated by the element nodal displacement vector q  : where matrix N| Γ 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);   ( = 4 ∼ +3) are  fundamental analytical solutions (in Cartesian coordinates) of the Airy stress function , which must be selected from complete second-order terms ( 2 ,  2 , ) to complete higherorder terms (see Table 1);   ( = 1 ∼ ) are  unknown constants.Obviously, such trial functions will directly lead to more reasonable stress fields satisfying both equilibrium and compatibility conditions.
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 q  : Substitution of ( 16) into the Π *  in (12) yields From the viewpoint of the element definition given in [12], matrix K * 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 q  is solved, the element stresses can be given by  = SM −1 Hq  .
(18) [40].Up to present, several plane HSF element models for both isotropic and anisotropic materials, including 8-node and 12-node quadrilateral elements [39][40][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.

A Typical HSF Element HSF-Q8-15𝛽
Differing from the usual models, the element shapes can be either convex or concave.The element nodal displacement vector q  is given by Assume that the displacements along each element edge are interpolated by the nodal displacements of each edge.Therefore, , V and corresponding matrix N (see (6) and ( 12)) of each element edge can be given as follows: with in which  0  (, ) ( = 1 ∼ 8) are the shape functions of the plane 8-node isoparametric serendipity elements where (  ,   ) are the isoparametric coordinates of node .
Fifteen analytical solutions   ( = 4 ∼ 18) 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 S in ( 14) is a 3×15 matrix.It can be seen that the stress fields possess thirdorder completeness in both  and .The resulting element model is denoted as HSF-Q8-15.
In order to evaluate the matrices M and H by Gaussian numerical integration procedure, the Cartesian coordinates  and  within an element should be expressed in terms of local coordinates (isoparametric coordinates).Let where (  ,   ) ( = 1 ∼ 8) 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.[42].Consider a 4-node quadrilateral element with drilling degrees of freedom shown in Figure 3. Differing from the usual displacement-based Table 2: Results at selected locations for the pure bending problem (Figure 2).models, the element shape is allowed to be either convex or concave.The element nodal displacement vector q  is defined as

A 4-Node HSF Quadrilateral Plane Membrane Element with Drilling Degrees of Freedom
where   ( = 1 ∼ 4) 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  = 7, that is, the first seven analytical solutions   ( = 1 ∼ 7) 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.
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   (, ) ( = 1 ∼ 8) are not complete.[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

The Shape-Free 8-Node Unsymmetric Element US-ATFQ8
in which Table 3: HSF-Q4-7 results of radial stress at selected points for a wedge subjected to a uniformly distributed load (Figure 4).where   and   ( = 1 ∼ 18) are the displacement solutions derived from the fundamental analytical solutions of the Airy stress function, and have been given in Table 1;   ( = 1 ∼ 14) are 14 unknown parameters related to   and V  ;   ( = 1 ∼ 4) are 4 unknown parameters related to   and V  ; are assumed displacements which possess third-order completeness in  and y; and are assumed fourth-order displacement fields.
For the fourth-order displacement fields   and V  , two additional relaxed point bubble conditions are introduced So,  2 and  4 can be expressed in terms of  1 and  3 as follows in which  1 = 0.5 Let Then, the displacement fields in (31) can be rewritten as Substitution of coordinates of eight nodes into (39) yields  where   ,   and   ( = 4 ∼ 18) 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.

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][31][32][33][34][35][36][37][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.

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

Figure 2 :
Figure 2: Pure bending problem for a cantilever beam.

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

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

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