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

Authors: Song Cen; Ming-Jue Zhou; Yan Shang
Journal: Mathematical Problems in Engineering
 (2013)

---

## 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 calledshape-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.


---

## Body

## 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 thehybrid 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 calledshape-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 theshape-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 theimproved unsymmetric element methods. For plane problem without body forces (or with only constant body forces), the analytical solutions ϕi 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
(1)∇4ϕi=∂4ϕi∂x4+2∂4ϕi∂x2∂y2+∂4ϕi∂y4=0,
and for the anisotropic case, the equation becomes
(2)C^11∂4ϕi∂y4+C^22∂4ϕi∂x4+(2C^12+C^66)∂4ϕi∂x2∂y2-2C^16∂4ϕi∂x∂y3-2C^26∂4ϕi∂x3∂y=0,
where C^ij=C^ji are the reduced elastic compliances and have been defined in [41].The fundamental analytical solutionsϕi 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 ϕi  (i=1~18) 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
ϕ
i
σ
x
i
σ
y
i
τ
x
y
i
U
-
i
V
-
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(μx2 + y2)/E
6xy/E
8
x
2
y
0
2y
−2x
−2μxy/E
[y2 − (2 + μ)x2]/E
9
xy
2
2x
0
−2y
[
x
2
-
(
2
+
μ
)
y
2
]
/
E
−2μxy/E
10
y
3
6y
0
0
6xy/E
−3(μy2 + x2)/E
11
x
3
y
0
6xy
−3x2
−(3μx2y + y3)/E
[
3
x
y
2
-
(
2
+
μ
)
x
3
]
/
E
12
xy
3
6xy
0
−3y2
[3x2y − (2 + μ)y3]/E
−(3μx2y + y3)/E
13
x
4
-
y
4
−12y2
12x2
0
−4(3xy2+ μx3)/E
4(3x2y + μy3)/E
14
6
x
2
y
2
-
x
4
-
y
4
12(x2-y2)
−12(x2-y2)
−24xy
4(1 +μ)(x3-3xy2)/E
4(1 +μ)(y3-3x2y)/E
15
x
3
y
2
-
x
y
4
2x(x2-6y2)
6xy2
−2y(3x2-2y2)
[x4-6(2 + μ)x2y2 + (3 + 2 μ)y4]/2E
[2(1 + 2μ)xy3-2 μx3y]/E
16
5
x
3
y
2
-
x
5
10x3
−10x(2x2-3y2)
−30x2y
2.5[(1 + 2μ)x4-6 μx2y2 − y4]/E
10xy[y2-(2+μ)x2]/E
17
x
2
y
3
-
x
4
y
6x2y
−2y(6x2-y2)
2x(2x2-3y2)
[2(1 + 2μ)x3y-2 μxy3]/E
[y4-6(2 + μ)x2y2 + (3+2 μ)x4]/2E
18
5
x
2
y
3
-
y
5
10y(3x2-2y2)
10y3
−30xy2
10xy[x2 − (2 + μ)y2]/E
2.5[(1 + 2μ)y4-6 μx2y2 − x4]/E
*Note: E→E/(1-μ2) and μ→μ/(1-μ) 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]:
(3)ΠC=ΠC*+VC*=12∬AeσTCσtdA-∫ΓeTTu-tds,(4)ΠC*=12∬AeσTCσtdA,VC*=-∫ΓeTTu-tds,(5)σ={σxσyτxy},T={TxTy},u-={u-v-},
where ΠC* is the complementary energy within the element; VC* 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 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 qe:
(6)u-=N-|Γqe,
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,
(7)σ={σxσyτxy}={∂2ϕ∂y2∂2ϕ∂x2-∂2ϕ∂x∂y}=R~(ϕ).
And the traction force vector T in (5) can be written as
(8)T={TxTy}=[l0m0ml]{σxσyτxy}=LR~(ϕ)withL=[l0m0ml],
where l and m are the direction cosines of the outer normal n of the element boundaries.Substituting (7) and (8) into (3) yields
(9)ΠC=ΠC*+VC*=12∬AeR~(ϕ)TCR~(ϕ)tdA-∫Γ[LR~(ϕ)]Tu-tds.
Let
(10)ϕ=∑i=1Nϕi+3βi=φβ,
with
(11)φ=[ϕ4ϕ5ϕ6⋯⋯ϕN+3],β=[β1β2β3⋯⋯βN]T,
where N is the number of the fundamental analytical solutions used for stress function ϕ in (10); ϕi(i=4~N+3) are N fundamental analytical solutions (in Cartesian coordinates) of the Airy stress function ϕ, which must be selected from complete second-order terms (x2, y2, xy) to complete higher-order terms (see Table 1); βi(i=1~N) are N 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
(12)ΠC=ΠC*+VC*=12βTMβ-βTHqe,(13)M=∬AeSTCStdA,H=∫ΓSTLTN-tds
with
(14)S=[∂2ϕ4∂y2∂2ϕ5∂y2⋯∂2ϕN+3∂y2∂2ϕ4∂x2∂2ϕ5∂x2⋯∂2ϕN+3∂x2-∂2ϕ4∂x∂y-∂2ϕ5∂x∂y⋯-∂2ϕN+3∂x∂y]=[σx4σx5σx6⋯σx(N+3)σy4σy5σy6⋯σy(N+3)τxy4τxy5τxy6⋯τxy(N+3)],
where σxi, σyi and τxyi (i=4~N+3) are the corresponding stresses of ϕi(i=4~N+3) (see Table 1), and all are expressed in Cartesian coordinates.According to the principle of minimum complementary energy, we require(15)∂ΠC∂β=0.
Then, the unknown constant vector β can be expressed in terms of the nodal displacement vector qe:
(16)β=M-1Hqe.Substitution of (16) into the ΠC* in (12) yields
(17)ΠC*=12qeTK*qe,K*=(M-1H)TH=HTM-1H.
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 vectorqe is solved, the element stresses can be given by
(18)σ=SM-1Hqe.


### 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 Figure1, 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 qe is given by(19)qe=[u1v1u2v2u3v3u4v4u5v5u6v6u7v7u8v8]T.8-node quadrilateral plane element with arbitrary shapes.
(a)
(b)
(c)
(d)
(e)Assume that the displacements along each element edge are interpolated by the nodal displacements of each edge. Therefore,u-, v- and corresponding matrix N- (see (6) and (12)) of each element edge can be given as follows:(20)u-12¯={u-v-}12¯=N-|η=-1qe,u-23¯={u-v-}23¯=N-|ξ=1qe,u-34¯={u-v-}34¯=N-|η=1qe,u-41¯={u-v-}41¯=N-|ξ=-1qe,
with
(21)N-=[N100N200⋯⋯N8000N100N20⋯⋯0N80],
in which Ni0(ξ,η) (i=1~8) are the shape functions of the plane 8-node isoparametric serendipity elements
(22)Ni0={-14(1+ξiξ)(1+ηiη)(1-ξiξ-ηiη)(i=1,2,3,4)12(1-ξ2)(1+ηiη)(i=5,7)12(1-η2)(1+ξiξ)(i=6,8),
where (ξi, ηi) are the isoparametric coordinates of node i.Fifteen analytical solutionsϕi (i=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 βi are introduced. Then, the matrix S in (14) is a 3×15 matrix. It can be seen that the stress fields possess third-order completeness in both x and y. The resulting element model is denoted as HSF-Q8-15β.In order to evaluate the matricesM and H by Gaussian numerical integration procedure, the Cartesian coordinates x and y within an element should be expressed in terms of local coordinates (isoparametric coordinates). Let
(23)x=∑i=18Ni0(ξ,η)xi,y=∑i=18Ni0(ξ,η)yi,
where (xi,yi) (i=1~8) are the Cartesian coordinates of the node i.Example 1 . (pure bending for a cantilever beam (Figure2)).
As shown in Figure2, a cantilever beam under plane stress condition is subjected to a constant bending moment M. 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]:
(24)σx=240cy-120,σy=0,τxy=0,u=-120Ex+240cExy,v=36Ey-120cEx2-36cEy2.
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 (Figure2).
Q8
Q9L
HSF-Q8-15β
Exact
Mesh 1
σ
x(0,10)
120.000
120.000
120.000
120.0
σ
x(0,0)
−120.000
−120.000
−120.000
−120.0
v(100,0) × 103
−12.000
−12.000
−12.000
−12.0
Mesh 2
σ
x(0,10)
56.447
120.000
120.000
120.0
σ
x(0,0)
−74.863
−120.000
−120.000
−120.0
v(100,0) × 103
−2.328
−12.000
−12.000
−12.0
Mesh 3
σ
x(0+,10)
13.665
120.000
120.000
120.0
σ
x(0,10-)
5.262
120.000
120.000
120.0
σ
x(0,0+)
−5.665
−120.000
−120.000
−120.0
σ
x(0+,0)
−14.299
−120.000
−120.000
−120.0
v(100,0) × 103
−0.477
−12.000
−12.000
−12.0
Mesh 4
σ
x(0,10)
—
—
120.000
120.0
σ
x(0,0)
—
—
−120.000
−120.0
v(100,0) × 103
—
—
−12.000
−12.0
Mesh 5
σ
x(0,10)
—
—
120.000
120.0
σ
x(0,0)
—
—
−120.000
−120.0
v(100,0) × 103
—
—
−12.000
−12.0
Mesh 6
σ
x(0+,10)
—
—
120.000
120.0
σ
x(0,10-)
—
—
120.000
120.0
σ
x(0,0+)
—
—
−120.000
−120.0
σ
x(0+,0)
—
—
−120.000
−120.0
v(100,0) × 103
—
—
−12.000
−12.0Figure 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 Figure3. Differing from the usual displacement-based models, the element shape is allowed to be either convex or concave. The element nodal displacement vector qe is defined as (25)qe=[u1v1θz1u2v2θz2u3v3θz3u4v4θz4]T,Figure 3
A 4-node HSF quadrilateral plane membrane element with drilling degrees of freedom.whereθzi(i=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
(26){u-v-}ij¯=[N-i0-(yj-yi)2N-iN-jN-j0(yj-yi)2N-iN-j0N-i(xj-xi)2N-iN-j0N-j-(xj-xi)2N-iN-j]×{uiviθziujvjθzj},(i=1,2,3,4j=2,3,4,1),N-1=12(1-s),N-2=12(1+s),-1≤s≤1.For matrixS defined in (14), let N=7, that is, the first seven analytical solutions ϕi(i=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 x and y. The resulting element model is denoted as HSF-Q4θ-7β.Example 2 . (a wedge subjected to a uniformly distributed load (Figure4)).
As shown in Figure4, a cantilever wedge is subjected to a uniformly distributed load q. 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]
σ
r at point A (0, 5)
7.6835(1.38%)
7.5894 (0.13%)
7.5806 (0.02%)
7.5792
σ
r at point B (1, 5)
−7.7920 (1.47%)
−7.7088 (0.39%)
−7.6832 (0.05%)
−7.6792Figure 4
A wedge subjected to a uniformly distributed load.


## 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]:
(3)ΠC=ΠC*+VC*=12∬AeσTCσtdA-∫ΓeTTu-tds,(4)ΠC*=12∬AeσTCσtdA,VC*=-∫ΓeTTu-tds,(5)σ={σxσyτxy},T={TxTy},u-={u-v-},
where ΠC* is the complementary energy within the element; VC* 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 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 qe:
(6)u-=N-|Γqe,
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,
(7)σ={σxσyτxy}={∂2ϕ∂y2∂2ϕ∂x2-∂2ϕ∂x∂y}=R~(ϕ).
And the traction force vector T in (5) can be written as
(8)T={TxTy}=[l0m0ml]{σxσyτxy}=LR~(ϕ)withL=[l0m0ml],
where l and m are the direction cosines of the outer normal n of the element boundaries.Substituting (7) and (8) into (3) yields
(9)ΠC=ΠC*+VC*=12∬AeR~(ϕ)TCR~(ϕ)tdA-∫Γ[LR~(ϕ)]Tu-tds.
Let
(10)ϕ=∑i=1Nϕi+3βi=φβ,
with
(11)φ=[ϕ4ϕ5ϕ6⋯⋯ϕN+3],β=[β1β2β3⋯⋯βN]T,
where N is the number of the fundamental analytical solutions used for stress function ϕ in (10); ϕi(i=4~N+3) are N fundamental analytical solutions (in Cartesian coordinates) of the Airy stress function ϕ, which must be selected from complete second-order terms (x2, y2, xy) to complete higher-order terms (see Table 1); βi(i=1~N) are N 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
(12)ΠC=ΠC*+VC*=12βTMβ-βTHqe,(13)M=∬AeSTCStdA,H=∫ΓSTLTN-tds
with
(14)S=[∂2ϕ4∂y2∂2ϕ5∂y2⋯∂2ϕN+3∂y2∂2ϕ4∂x2∂2ϕ5∂x2⋯∂2ϕN+3∂x2-∂2ϕ4∂x∂y-∂2ϕ5∂x∂y⋯-∂2ϕN+3∂x∂y]=[σx4σx5σx6⋯σx(N+3)σy4σy5σy6⋯σy(N+3)τxy4τxy5τxy6⋯τxy(N+3)],
where σxi, σyi and τxyi (i=4~N+3) are the corresponding stresses of ϕi(i=4~N+3) (see Table 1), and all are expressed in Cartesian coordinates.According to the principle of minimum complementary energy, we require(15)∂ΠC∂β=0.
Then, the unknown constant vector β can be expressed in terms of the nodal displacement vector qe:
(16)β=M-1Hqe.Substitution of (16) into the ΠC* in (12) yields
(17)ΠC*=12qeTK*qe,K*=(M-1H)TH=HTM-1H.
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 vectorqe is solved, the element stresses can be given by
(18)σ=SM-1Hqe.


## 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 Figure1, 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 qe is given by(19)qe=[u1v1u2v2u3v3u4v4u5v5u6v6u7v7u8v8]T.8-node quadrilateral plane element with arbitrary shapes.
(a)
(b)
(c)
(d)
(e)Assume that the displacements along each element edge are interpolated by the nodal displacements of each edge. Therefore,u-, v- and corresponding matrix N- (see (6) and (12)) of each element edge can be given as follows:(20)u-12¯={u-v-}12¯=N-|η=-1qe,u-23¯={u-v-}23¯=N-|ξ=1qe,u-34¯={u-v-}34¯=N-|η=1qe,u-41¯={u-v-}41¯=N-|ξ=-1qe,
with
(21)N-=[N100N200⋯⋯N8000N100N20⋯⋯0N80],
in which Ni0(ξ,η) (i=1~8) are the shape functions of the plane 8-node isoparametric serendipity elements
(22)Ni0={-14(1+ξiξ)(1+ηiη)(1-ξiξ-ηiη)(i=1,2,3,4)12(1-ξ2)(1+ηiη)(i=5,7)12(1-η2)(1+ξiξ)(i=6,8),
where (ξi, ηi) are the isoparametric coordinates of node i.Fifteen analytical solutionsϕi (i=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 βi are introduced. Then, the matrix S in (14) is a 3×15 matrix. It can be seen that the stress fields possess third-order completeness in both x and y. The resulting element model is denoted as HSF-Q8-15β.In order to evaluate the matricesM and H by Gaussian numerical integration procedure, the Cartesian coordinates x and y within an element should be expressed in terms of local coordinates (isoparametric coordinates). Let
(23)x=∑i=18Ni0(ξ,η)xi,y=∑i=18Ni0(ξ,η)yi,
where (xi,yi) (i=1~8) are the Cartesian coordinates of the node i.Example 1 . (pure bending for a cantilever beam (Figure2)).
As shown in Figure2, a cantilever beam under plane stress condition is subjected to a constant bending moment M. 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]:
(24)σx=240cy-120,σy=0,τxy=0,u=-120Ex+240cExy,v=36Ey-120cEx2-36cEy2.
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 (Figure2).
Q8
Q9L
HSF-Q8-15β
Exact
Mesh 1
σ
x(0,10)
120.000
120.000
120.000
120.0
σ
x(0,0)
−120.000
−120.000
−120.000
−120.0
v(100,0) × 103
−12.000
−12.000
−12.000
−12.0
Mesh 2
σ
x(0,10)
56.447
120.000
120.000
120.0
σ
x(0,0)
−74.863
−120.000
−120.000
−120.0
v(100,0) × 103
−2.328
−12.000
−12.000
−12.0
Mesh 3
σ
x(0+,10)
13.665
120.000
120.000
120.0
σ
x(0,10-)
5.262
120.000
120.000
120.0
σ
x(0,0+)
−5.665
−120.000
−120.000
−120.0
σ
x(0+,0)
−14.299
−120.000
−120.000
−120.0
v(100,0) × 103
−0.477
−12.000
−12.000
−12.0
Mesh 4
σ
x(0,10)
—
—
120.000
120.0
σ
x(0,0)
—
—
−120.000
−120.0
v(100,0) × 103
—
—
−12.000
−12.0
Mesh 5
σ
x(0,10)
—
—
120.000
120.0
σ
x(0,0)
—
—
−120.000
−120.0
v(100,0) × 103
—
—
−12.000
−12.0
Mesh 6
σ
x(0+,10)
—
—
120.000
120.0
σ
x(0,10-)
—
—
120.000
120.0
σ
x(0,0+)
—
—
−120.000
−120.0
σ
x(0+,0)
—
—
−120.000
−120.0
v(100,0) × 103
—
—
−12.000
−12.0Figure 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 Figure3. Differing from the usual displacement-based models, the element shape is allowed to be either convex or concave. The element nodal displacement vector qe is defined as (25)qe=[u1v1θz1u2v2θz2u3v3θz3u4v4θz4]T,Figure 3
A 4-node HSF quadrilateral plane membrane element with drilling degrees of freedom.whereθzi(i=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
(26){u-v-}ij¯=[N-i0-(yj-yi)2N-iN-jN-j0(yj-yi)2N-iN-j0N-i(xj-xi)2N-iN-j0N-j-(xj-xi)2N-iN-j]×{uiviθziujvjθzj},(i=1,2,3,4j=2,3,4,1),N-1=12(1-s),N-2=12(1+s),-1≤s≤1.For matrixS defined in (14), let N=7, that is, the first seven analytical solutions ϕi(i=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 x and y. The resulting element model is denoted as HSF-Q4θ-7β.Example 2 . (a wedge subjected to a uniformly distributed load (Figure4)).
As shown in Figure4, a cantilever wedge is subjected to a uniformly distributed load q. 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]
σ
r at point A (0, 5)
7.6835(1.38%)
7.5894 (0.13%)
7.5806 (0.02%)
7.5792
σ
r at point B (1, 5)
−7.7920 (1.47%)
−7.7088 (0.39%)
−7.6832 (0.05%)
−7.6792Figure 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]
(27)Ke=∬AeB-TDB^tdA=∬-11B-*T|J|DB^|J|tdξdη=∬-11B-*TDB^tdξdη,
where |J| is the Jacobian determinant; D is the elasticity matrix; B- is the element strain matrix for conventional 8-node plane isoparametric element Q8:
(28)B-=1|J|B-*=1|J|×[N1,x0N2,x0N3,x0⋯N8,x00N1,y0N2,y0N3,y⋯0N8,yN1,yN1,xN2,yN2,xN3,yN3,x⋯N8,yN8,x],
in which Ni(ξ,η)(i=1~8) have been given by (22); B^ is given by
(29)B^=[M1,x0M2,x0M3,x0⋯M8,x00M1,y0M2,y0M3,y⋯0M8,yM1,yM1,xM2,yM2,xM3,yM3,x⋯M8,yM8,x],
where Mi(x,y)(i=1~8) are Cartesian shape functions, given by
(30)[11111111x1x2x3x4x5x6x7x8y1y2y3y4y5y6y7y8x12x22x32x42x52x62x72x82x1y1x2y2x3y3x4y4x5y5x6y6x7y7x8y8y12y22y32y42y52y62y72y82x12y1x22y2x32y3x42y4x52y5x62y6x72y7x82y8x1y12x2y22x3y32x4y42x5y52x6y62x7y72x8y82]×{M1M2M3M4M5M6M7M8}={1xyx2xyy2x2yxy2},
in which (xi, yi) are the Cartesian coordinates of node i.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 Mi(x,y)(i=1~8) 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
(31)ue={uv}=[U-U^V-V^]{α^β^}={uα+uβvα+vβ},
in which(32)U-=[U-1U-2U-3U-4U-5U-6U-7U-8U-9U-10U-11U-12U-13U-14]V-=[V-1V-2V-3V-4V-5V-6V-7V-8V-9V-10V-11V-12V-13V-14],U^=[U-15U-16U-17U-18],V^=[V-15V-16V-17V-18],α^=[α1α2α3α4α5α6α7α8α9α10α11α12α13α14]T,β^=[β1β2β3β4]T,whereU-i and V-i(i=1~18) are the displacement solutions derived from the fundamental analytical solutions of the Airy stress function, and have been given in Table 1; αi(i=1~14) are 14 unknown parameters related to uα and vα; βi  (i=1~4) are 4 unknown parameters related to uβ and vβ;
(33){uαvα}={U-α^V-α^},
are assumed displacements which possess third-order completeness in x and y; and
(34){uβvβ}={U^β^V^β^},
are assumed fourth-order displacement fields.For the fourth-order displacement fieldsuβ and vβ, two additional relaxed point bubble conditions are introduced
(35)∑i=58uβ(xi,yi)=∑i=58[U-15(xi,yi)β1+U-16(xi,yi)β2U-17+U-17(xi,yi)β3+U-18(xi,yi)β4]=0,∑i=58vβ(xi,yi)=∑i=58[V-15(xi,yi)β1+V-16(xi,yi)β2V-17+V-17(xi,yi)β3+V-18(xi,yi)β4]=0.
So, β2 and β4 can be expressed in terms of β1 and β3 as follows
(36){β2β4}=-[a2a4b2b4]-1[a1a3b1b3]{β1β3}=[r11r12r21r22]{β1β3},
in which
(37)a1=0.5∑i=58xi4-(6+3μ)∑i=58xi2yi2+(1.5+μ)∑i=58yi4,a2=(2.5+5μ)∑i=58xi4-15μ∑i=58xi2yi2-2.5∑i=58yi4,a3=(2+4μ)∑i=58xi3yi-2μ∑i=58xiyi3,a4=10∑i=58xi3yi-(20+10μ)∑i=58xiyi3,b1=-2μ∑i=58xi3yi+(2+4μ)∑i=58xiyi3,b2=-(20+10μ)∑i=58xi3yi+10∑i=58xiyi3,b3=(1.5+μ)∑i=58xi4-(6+3μ)∑i=58xi2yi2+0.5∑i=58yi4,b4=-2.5∑i=58xi4-15μ∑i=58xi2yi2+(2.5+5μ)∑i=58yi4.
Let
(38)β1=α15β3=α16.
Then, the displacement fields in (31) can be rewritten as
(39)ue={uv}=[U-U̿V-V̿]{α^{α15α16}}=Pα,
with
(40)U̿=[(U-15+r11U-16+r21U-18)(U-17+r12U-16+r22U-18)]V̿=[(V-15+r11V-16+r21V-18)(V-17+r12V-16+r22V-18)],P=[U-U̿V-V̿],α=[α1α2α3⋯α14α15α16]T.Substitution of coordinates of eight nodes into (39) yields
(41)d-α=qe,
in which
(42)d-=[P(x1,y1)P(x2,y2)⋮P(x8,y8)]16×16,whereqe is the element nodal displacement vector and has been given by (19). Then, αi(i=1~16) can be solved by
(43)α=d--1qe.
It should be noted that so long as the element nodes are not in coincidence with each other, the matrix d- 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
(44)Ke=∬AeB-TS^d--1tdA=∬-11B-*T|J|S^d--1|J|tdξdη=∬-11B-*TS^tdξdη·d--1,
where B- and B-* are given by (28), and (45)S^=[000σx4σx5⋯σx14σx15+r11σx16+r21σx18σx17+r12σx16+r22σx18000σy4σy5⋯σy14σy15+r11σy16+r21σy18σy17+r12σy16+r22σy18000τxy4τxy5⋯τxy14τxy15+r11τxy16+r21τxy18τxy17+r12τxy16+r22τxy18]3×16,whereσxi, σyi and τxyi(i=4~18) have been given in Table 1. And the element stresses can be evaluated by
(46)σ={σxσyτxy}=S^d--1qe.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 (Figure5)).
As shown in Figure5, a cantilever beam under plane stress condition is subjected to a linear bending moment caused by a shear force P 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 (Figure5).
Q8  (3 × 3)
HSF-Q8-15β  (5 × 5)
US-QUAD8  (3 × 3)
US-ATFQ8(5 × 5 or 3 × 3)
Exact
σ
x (0,10)
−66.871
−120.00
−66.871
−120.000
−120.0
Mesh 1
σ
x (0,0)
66.871
120.00
66.871
120.000
120.0
v(100,0) × 103
6.282
8.558
6.282
8.0460
8.046
σ
x (0,10)
−35.248
−113.63
−74.317
−120.000
−120.0
Mesh 2
σ
x (0,0)
67.040
119.63
71.781
120.000
120.0
v(100,0) × 103
1.611
7.104
6.415
8.0460
8.046
σ
x
(0+,10)
−11.307
−68.645
−66.349
−120.000
−120.0
σ
x (0,10-)
−55.699
−66.980
−67.974
−120.000
−120.0
Mesh 3
σ
x (0,0+)
55.305
66.146
67.995
120.000
120.0
σ
x (0+,0)
11.073
69.328
66.328
120.000
120.0
v(100,0) × 103
4.520
6.350
6.292
8.0460
8.046
σ
x(0,10)
−72.471
−120.00
−120.0
Mesh 4
σ
x(0,0)
Failed
76.265
Failed
120.00
120.0
v(100,0) × 103
6.663
8.0460
8.046
σ
x(0,10)
−77.515
−58.158
−120.000
−120.0
Mesh 5
σ
x(0,0)
Failed
66.306
59.158
120.000
120.0
v(100,0) × 103
6.729
5.471
8.0460
8.046
σ
x(0+,10)
−70.346
−42.671
−120.000
−120.0
σ
x(0,10-)
−66.835
−44.636
−120.000
−120.0
Mesh 6
σ
x(0,0+)
Failed
66.877
40.974
120.000
120.0
σ
x(0+,0)
68.512
41.048
120.000
120.0
v(100,0)  × 103
6.373
3.752
8.0460
8.046
σ
x(0,10)
−105.69
−120.00
−105.69
−120.000
−120.0
Mesh 7
σ
x(0,0)
105.69
120.00
105.69
120.000
120.0
v(20,0) × 104
3.485
3.975
3.485
3.6600
3.66
σ
x(0,10)
−115.25
−119.02
−118.99
−120.000
−120.0
Mesh 8
σ
x(0,0)
115.25
119.02
118.99
120.000
120.0
v(20,0) × 104
3.3423
3.796
3.5570
3.6600
3.66
σ
x (0,20)
−118.68
−120.00
−118.68
−120.000
−120.0
Mesh 9
σ
x(0,0)
118.68
120.00
118.68
120.000
120.0
v(10,0) × 105
11.1963
15.750
11.1963
13.200
13.2
σ
x(0,20)
−97.066
−121.61
−118.94
−120.00
b/−119.65a
−120.0
Mesh 10
σ
x(0,0)
138.21
124.61
117.74
120.00
b/119.36a
120.0
v(10,0) × 105
12.8162
13.738
12.9157
13.200
b/13.190a
13.2
σ
x(0,10)
−24.732
−36.104
−120.000
−120.0
Mesh 11
σ
x(0,0)
Failed
20.048
36.230
120.000
120.0
v(100,0) × 103
1.517
3.1470
8.0460
8.046
a3 × 3 integration scheme; b5 × 5 integration scheme.Figure 5
Linear bending problem for a cantilever beam.


## 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]
(27)Ke=∬AeB-TDB^tdA=∬-11B-*T|J|DB^|J|tdξdη=∬-11B-*TDB^tdξdη,
where |J| is the Jacobian determinant; D is the elasticity matrix; B- is the element strain matrix for conventional 8-node plane isoparametric element Q8:
(28)B-=1|J|B-*=1|J|×[N1,x0N2,x0N3,x0⋯N8,x00N1,y0N2,y0N3,y⋯0N8,yN1,yN1,xN2,yN2,xN3,yN3,x⋯N8,yN8,x],
in which Ni(ξ,η)(i=1~8) have been given by (22); B^ is given by
(29)B^=[M1,x0M2,x0M3,x0⋯M8,x00M1,y0M2,y0M3,y⋯0M8,yM1,yM1,xM2,yM2,xM3,yM3,x⋯M8,yM8,x],
where Mi(x,y)(i=1~8) are Cartesian shape functions, given by
(30)[11111111x1x2x3x4x5x6x7x8y1y2y3y4y5y6y7y8x12x22x32x42x52x62x72x82x1y1x2y2x3y3x4y4x5y5x6y6x7y7x8y8y12y22y32y42y52y62y72y82x12y1x22y2x32y3x42y4x52y5x62y6x72y7x82y8x1y12x2y22x3y32x4y42x5y52x6y62x7y72x8y82]×{M1M2M3M4M5M6M7M8}={1xyx2xyy2x2yxy2},
in which (xi, yi) are the Cartesian coordinates of node i.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 Mi(x,y)(i=1~8) 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
(31)ue={uv}=[U-U^V-V^]{α^β^}={uα+uβvα+vβ},
in which(32)U-=[U-1U-2U-3U-4U-5U-6U-7U-8U-9U-10U-11U-12U-13U-14]V-=[V-1V-2V-3V-4V-5V-6V-7V-8V-9V-10V-11V-12V-13V-14],U^=[U-15U-16U-17U-18],V^=[V-15V-16V-17V-18],α^=[α1α2α3α4α5α6α7α8α9α10α11α12α13α14]T,β^=[β1β2β3β4]T,whereU-i and V-i(i=1~18) are the displacement solutions derived from the fundamental analytical solutions of the Airy stress function, and have been given in Table 1; αi(i=1~14) are 14 unknown parameters related to uα and vα; βi  (i=1~4) are 4 unknown parameters related to uβ and vβ;
(33){uαvα}={U-α^V-α^},
are assumed displacements which possess third-order completeness in x and y; and
(34){uβvβ}={U^β^V^β^},
are assumed fourth-order displacement fields.For the fourth-order displacement fieldsuβ and vβ, two additional relaxed point bubble conditions are introduced
(35)∑i=58uβ(xi,yi)=∑i=58[U-15(xi,yi)β1+U-16(xi,yi)β2U-17+U-17(xi,yi)β3+U-18(xi,yi)β4]=0,∑i=58vβ(xi,yi)=∑i=58[V-15(xi,yi)β1+V-16(xi,yi)β2V-17+V-17(xi,yi)β3+V-18(xi,yi)β4]=0.
So, β2 and β4 can be expressed in terms of β1 and β3 as follows
(36){β2β4}=-[a2a4b2b4]-1[a1a3b1b3]{β1β3}=[r11r12r21r22]{β1β3},
in which
(37)a1=0.5∑i=58xi4-(6+3μ)∑i=58xi2yi2+(1.5+μ)∑i=58yi4,a2=(2.5+5μ)∑i=58xi4-15μ∑i=58xi2yi2-2.5∑i=58yi4,a3=(2+4μ)∑i=58xi3yi-2μ∑i=58xiyi3,a4=10∑i=58xi3yi-(20+10μ)∑i=58xiyi3,b1=-2μ∑i=58xi3yi+(2+4μ)∑i=58xiyi3,b2=-(20+10μ)∑i=58xi3yi+10∑i=58xiyi3,b3=(1.5+μ)∑i=58xi4-(6+3μ)∑i=58xi2yi2+0.5∑i=58yi4,b4=-2.5∑i=58xi4-15μ∑i=58xi2yi2+(2.5+5μ)∑i=58yi4.
Let
(38)β1=α15β3=α16.
Then, the displacement fields in (31) can be rewritten as
(39)ue={uv}=[U-U̿V-V̿]{α^{α15α16}}=Pα,
with
(40)U̿=[(U-15+r11U-16+r21U-18)(U-17+r12U-16+r22U-18)]V̿=[(V-15+r11V-16+r21V-18)(V-17+r12V-16+r22V-18)],P=[U-U̿V-V̿],α=[α1α2α3⋯α14α15α16]T.Substitution of coordinates of eight nodes into (39) yields
(41)d-α=qe,
in which
(42)d-=[P(x1,y1)P(x2,y2)⋮P(x8,y8)]16×16,whereqe is the element nodal displacement vector and has been given by (19). Then, αi(i=1~16) can be solved by
(43)α=d--1qe.
It should be noted that so long as the element nodes are not in coincidence with each other, the matrix d- 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
(44)Ke=∬AeB-TS^d--1tdA=∬-11B-*T|J|S^d--1|J|tdξdη=∬-11B-*TS^tdξdη·d--1,
where B- and B-* are given by (28), and (45)S^=[000σx4σx5⋯σx14σx15+r11σx16+r21σx18σx17+r12σx16+r22σx18000σy4σy5⋯σy14σy15+r11σy16+r21σy18σy17+r12σy16+r22σy18000τxy4τxy5⋯τxy14τxy15+r11τxy16+r21τxy18τxy17+r12τxy16+r22τxy18]3×16,whereσxi, σyi and τxyi(i=4~18) have been given in Table 1. And the element stresses can be evaluated by
(46)σ={σxσyτxy}=S^d--1qe.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 (Figure5)).
As shown in Figure5, a cantilever beam under plane stress condition is subjected to a linear bending moment caused by a shear force P 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 (Figure5).
Q8  (3 × 3)
HSF-Q8-15β  (5 × 5)
US-QUAD8  (3 × 3)
US-ATFQ8(5 × 5 or 3 × 3)
Exact
σ
x (0,10)
−66.871
−120.00
−66.871
−120.000
−120.0
Mesh 1
σ
x (0,0)
66.871
120.00
66.871
120.000
120.0
v(100,0) × 103
6.282
8.558
6.282
8.0460
8.046
σ
x (0,10)
−35.248
−113.63
−74.317
−120.000
−120.0
Mesh 2
σ
x (0,0)
67.040
119.63
71.781
120.000
120.0
v(100,0) × 103
1.611
7.104
6.415
8.0460
8.046
σ
x
(0+,10)
−11.307
−68.645
−66.349
−120.000
−120.0
σ
x (0,10-)
−55.699
−66.980
−67.974
−120.000
−120.0
Mesh 3
σ
x (0,0+)
55.305
66.146
67.995
120.000
120.0
σ
x (0+,0)
11.073
69.328
66.328
120.000
120.0
v(100,0) × 103
4.520
6.350
6.292
8.0460
8.046
σ
x(0,10)
−72.471
−120.00
−120.0
Mesh 4
σ
x(0,0)
Failed
76.265
Failed
120.00
120.0
v(100,0) × 103
6.663
8.0460
8.046
σ
x(0,10)
−77.515
−58.158
−120.000
−120.0
Mesh 5
σ
x(0,0)
Failed
66.306
59.158
120.000
120.0
v(100,0) × 103
6.729
5.471
8.0460
8.046
σ
x(0+,10)
−70.346
−42.671
−120.000
−120.0
σ
x(0,10-)
−66.835
−44.636
−120.000
−120.0
Mesh 6
σ
x(0,0+)
Failed
66.877
40.974
120.000
120.0
σ
x(0+,0)
68.512
41.048
120.000
120.0
v(100,0)  × 103
6.373
3.752
8.0460
8.046
σ
x(0,10)
−105.69
−120.00
−105.69
−120.000
−120.0
Mesh 7
σ
x(0,0)
105.69
120.00
105.69
120.000
120.0
v(20,0) × 104
3.485
3.975
3.485
3.6600
3.66
σ
x(0,10)
−115.25
−119.02
−118.99
−120.000
−120.0
Mesh 8
σ
x(0,0)
115.25
119.02
118.99
120.000
120.0
v(20,0) × 104
3.3423
3.796
3.5570
3.6600
3.66
σ
x (0,20)
−118.68
−120.00
−118.68
−120.000
−120.0
Mesh 9
σ
x(0,0)
118.68
120.00
118.68
120.000
120.0
v(10,0) × 105
11.1963
15.750
11.1963
13.200
13.2
σ
x(0,20)
−97.066
−121.61
−118.94
−120.00
b/−119.65a
−120.0
Mesh 10
σ
x(0,0)
138.21
124.61
117.74
120.00
b/119.36a
120.0
v(10,0) × 105
12.8162
13.738
12.9157
13.200
b/13.190a
13.2
σ
x(0,10)
−24.732
−36.104
−120.000
−120.0
Mesh 11
σ
x(0,0)
Failed
20.048
36.230
120.000
120.0
v(100,0) × 103
1.517
3.1470
8.0460
8.046
a3 × 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.



---

*Source: 491626-2013-09-05.xml*