# A Novel Shape-Free Plane Quadratic Polygonal Hybrid Stress-Function Element

Authors: Pei-Lei Zhou; Song Cen
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

A novel plane quadratic shape-free hybrid stress-function (HS-F) polygonal element is developed by employing the principle of minimum complementary energy and the fundamental analytical solutions of the Airy stress function. Without construction of displacement interpolation function, the formulations of the new model are much simpler than those of the displacement-based polygonal elements and can be degenerated into triangular or quadrilateral elements directly. In particular, it is quite insensitive to various mesh distortions and even can keep precision when element shape is concave. Furthermore, the element does not show any spurious zero energy modes. Numerical examples show the excellent performance of the new element, denoted by HSF-AP-19β, in both displacement and stress solutions.


---

## Body

## 1. Introduction

The conventional triangular and quadrilateral elements are widely used in 2D finite element analyses [1, 2]. However, when simulating microstructure problems, inconveniences may still exist. For instance, the shape of one crystal of plane polycrystalline materials is usually regarded as a polygon. If these microstructures are meshed by those conventional triangular and quadrilateral elements, the density of the elements will be extremely high. But if each crystal is modeled only by one polygonal element, the related computation costs can be dramatically reduced [3]. Actually, the polygonal elements are more flexible for meshing structures with complex shapes [4–6] or connecting different kinds of meshes [7, 8], so that they possess obvious advantages for simulating heterogeneous materials [9–13].Compared with the conventional finite elements, it is more difficult to formulate the interpolation functions for polygonal elements because of the arbitrariness of the element shape. Many researchers have made great efforts on developing effective techniques to solve this problem. In 1975, Wachspress [14] firstly suggested to take the rational functions as the element shape functions. These functions satisfy Kronecker-delta property and continuity requirement along the polygon and can keep linear on the adjacent edges and zero on the opposite edge. However, their construction procedures were relatively complicated. Later, a series of improvements on the Wachspress interpolation were successfully proposed [2, 15–23]. Beside the Wachspress method, the Laplace functions were also utilized to build the interpolation functions for polygonal elements, such as the works proposed by Wang and Li [24] and Sukumar and Tabarraei [25]. Another mentionable approach is the barycentric coordinate method [26–28] which can be applied to develop both conventional and polygonal elements. During recent years, new polygonal elements are still emerging in many literatures. Dohrmann et al. [7] developed a transition element for connecting dissimilar finite element meshes; Peters and Heymsfield [29] developed constant strain elements consisting of an arbitrary number of nodes; Liu et al. [30, 31] constructed several linear polygonal elements by the smoothed finite element method; Chen et al. [32] proposed a quadratic spline element based on area coordinates and B-net method; Sukumar [33] proposed an arbitrary quadratic polygonal element by applying quadratic maximum-entropy serendipity shape functions. All the elements mentioned above belong to displacement-based elements, in which the displacement interpolation functions satisfying compatibility requirements must be considered.In order to avoid the difficulties in constructions of displacement-based polygonal elements, some researchers also utilized the hybrid finite element method to develop polygonal elements. Based on the principle of minimum complementary energy, Ghosh and Mallett [34] proposed the Voronoi cell finite element method (VCFEM), in which no displacement interpolation function is needed. By applying the Hellinger-Reissner variational principle, Zhang and Katsube [35, 36] developed the hybrid polygonal element (HPE) to simulate materials with inclusions and holes. Peng et al. [37, 38] proposed a novel base force element method by using the principle of minimum complementary energy and developed quadratic polygonal elements to solve problems with concave polygonal meshes.Recently, Cen et al. [39–44] proposed a hybrid stress-function (HS-F) element method based on the principle of minimum complementary energy and the fundamental analytical solutions of Airy stress function. They constructed several elements with excellent performance, such as 8-node and 12-node plane quadrilateral elements and 4-node plane quadrilateral element with drilling degrees of freedom. These quadrilateral elements are not sensitive to severely distorted meshes and can work well even when the element shapes degenerate into concave quadrilaterals or triangles. So, they are called shape-free finite elements since their performances are independent to the element shapes. Furthermore, this HS-F method can be directly and easily extended to develop arbitrary polygonal element due to its flexible theoretical frame. In the paper, a quadratic polygonal element will be developed by the above HS-F approach.The arrangements of this paper are as follows. In Section2, the 2D HS-F element method is briefly reviewed. Then, a new quadratic polygonal HS-F element, denoted by HSF-AP-19β, is formulated in Section 3. Several standard numerical examples are performed in Section 4 to validate the high performance of the new element. Finally, some conclusion remarks are given in Section 5.


## 2. Brief Reviews on the Plane HS-F Element Method

For a 2D finite element model, the complementary energy functional can be written as [39–44](1)
Π
C
=
1
2
∬
A
e
σ
T
C
σ
t
d
A
-
∫
Γ
e
T
T
u
-
t
d
s
,with(2)
σ
=
σ
x
σ
y
τ
x
y
,
T
=
T
x
T
y
=
l
0
m
0
m
l
σ
x
σ
y
τ
x
y
=
L
σ
,
u
-
=
u
-
v
-
,in which t is the thickness of the element; A
e is the element area; Γ
e is the element boundary; σ is the element stress vector; T is the traction force vector along the element boundaries; l and m are the direction cosines of the outer normal n
→ of the element boundaries; C is the elasticity matrix of compliances, and, for isotropic cases, it can be expressed by(3)
C
=
1
E
′
1
-
μ
′
0
-
μ
′
1
0
0
0
2
(
1
+
μ
′
)
,
plane
stress:
E
′
=
E
,
μ
′
=
μ
,
plane
strain:
E
′
=
E
1
-
μ
2
,
μ
′
=
μ
1
-
μ
,
E
,
Young’s modulus;
μ
,
Poisson’s ratio
.
u
- is the assumed displacement vector along element boundaries and can be interpolated by the element nodal displacement vector q
e:(4)
u
-
=
u
-
v
-
=
N
-
Γ
q
e
,where N
-
Γ is the interpolation function matrix for element boundary displacements.The element stress fields are assumed as follows:(5)
σ
=
σ
0
+
σ
*
=
σ
x
0
σ
y
0
τ
x
y
0
+
σ
x
*
σ
y
*
τ
x
y
*
=
S
β
+
σ
*
,
β
=
β
1
β
2
β
3
⋯
⋯
β
N
T
,where β
i (i
=
1
~
N) are N unknown stress parameters; S is the stress solution matrix:(6)
S
=
σ
x
1
σ
x
2
σ
x
3
⋯
σ
x
k
σ
y
1
σ
y
2
σ
y
3
⋯
σ
y
k
τ
x
y
1
τ
x
y
2
τ
x
y
3
⋯
τ
x
y
k
3
×
k
.The components (stress interpolation functions) in S are all derived from k fundamental solutions of stress function ϕ. The first nineteen solutions are given in Table 1.Table 1
Nineteen fundamental analytical solutions of the Airy stress function and resulting stress solutions for plane problem.
(a)
i
1
2
3
4
5
6
7
8
9
10
ϕ
i
x
2
x
y
y
2
x
3
x
2
y
x
y
2
y
3
x
3
y
x
y
3
x
4
-
y
4
σ
x
i
0
0
2
0
0
2
x
6
y
0
6
x
y
-
12
y
2
σ
y
i
2
0
0
6
x
2
y
0
0
6
x
y
0
12
x
2
τ
x
y
i
0
-
1
0
0
-
2
x
-
2
y
0
-
3
x
2
-
3
y
2
0
(b)
i
11
12
13
14
15
16
ϕ
i
6
x
2
y
2
-
x
4
-
y
4
x
3
y
2
-
x
y
4
5
x
3
y
2
-
x
5
x
2
y
3
-
x
4
y
5
x
2
y
3
-
y
5
x
5
y
-
x
y
5
σ
x
i
12
(
x
2
-
y
2
)
2
x
(
x
2
-
6
y
2
)
10
x
3
6
x
2
y
10
y
(
3
x
2
-
2
y
2
)
-
20
x
y
3
σ
y
i
-
12
(
x
2
-
y
2
)
6
x
y
2
-
10
x
(
2
x
2
-
3
y
2
)
-
2
y
(
6
x
2
-
y
2
)
10
y
3
20
x
3
y
τ
x
y
i
-
24
x
y
-
2
y
(
3
x
2
-
2
y
2
)
-
30
x
2
y
2
x
(
2
x
2
-
3
y
2
)
-
30
x
y
2
-
5
x
4
+
5
y
4
(c)
i
17
18
19
ϕ
i
10
x
3
y
3
-
3
x
5
y
-
3
x
y
5
x
6
-
10
x
4
y
2
+
5
x
2
y
4
y
6
-
10
x
2
y
4
+
5
x
4
y
2
σ
x
i
60
x
3
y
-
60
x
y
3
-
20
x
4
+
60
x
2
y
2
30
y
4
-
120
x
2
y
2
+
10
x
4
σ
y
i
60
x
y
3
-
60
x
3
y
30
x
4
-
120
x
2
y
2
+
10
y
4
-
20
y
4
+
60
x
2
y
2
τ
x
y
i
-
90
x
2
y
2
+
15
x
4
+
15
y
4
80
x
3
y
-
40
x
y
3
80
x
y
3
-
40
x
3
yσ
* is a particular solution corresponding to body forces, and, for constant body forces b
x and b
y, it can be given by(7)
σ
*
=
σ
x
*
σ
y
*
τ
x
y
*
=
0
0
-
b
y
x
-
b
x
y
.Therefore, these components are also fundamental stress solutions which satisfy all governing equations. This is the key point for developing the hybrid stress-function (HS-F) element models. According to the principle of minimum complementary energy [39–44], we have(8)
β
=
M
-
1
H
q
e
-
M
*
with
M
=
∬
A
e
S
T
C
S
t
d
A
,
M
*
=
∬
A
e
S
T
C
σ
*
t
d
A
.By applying the principle of the minimum complementary energy again, the final finite element equation of the hybrid stress-function (HS-F) element method can be written as(9)
∑
e
K
e
q
e
=
∑
e
p
e
,where K
e is the element stiffness matrix of 2D hybrid stress-function elements:(10)
K
e
=
H
T
M
-
1
H
,with(11)
H
=
∫
Γ
e
S
T
L
T
N
-
Γ
t
d
s
;
p
e is the element equivalent load vector. Its part caused by concentrated and distributed line forces can be determined by the standard procedure for the conventional finite elements. And the part caused by body forces is given by(12)
p
e
=
H
T
M
-
1
M
*
-
V
T
with
V
=
∫
Γ
e
σ
*
T
L
T
N
-
Γ
t
d
s
.Once the element nodal displacement vectorq
e is solved, the stresses at any point within the element can be given by(13)
σ
=
σ
x
σ
y
τ
x
y
=
S
M
-
1
(
-
M
*
+
H
q
e
)
+
σ
*
.


## 3. Formulations of a New Plane Quadratic Polygonal HS-F Element

As shown in Figure1, consider an arbitrary n-sided quadratic polygonal HS-F element; the element nodal displacement vector q
e is defined as(14)
q
e
=
u
1
v
1
u
2
v
2
⋯
⋯
u
2
n
-
1
v
2
n
-
1
u
2
n
v
2
n
T
,in which u
i and v
i (i
=
1
-
2
n) are the nodal displacements in x- and y-directions, respectively.Figure 1
An arbitraryn-sided polygonal HS-F element.


### 3.1. Displacement Interpolation along the Element Boundary

The element boundary displacement vector can be described by the element nodal displacement vector (see (4)). For element edge with three nodes i, j
=
i
+
1, and k
=
j
+
1 (if k
>
2
n, let k
=
1), let(15)
x
=
N
-
1
ξ
x
i
+
N
-
2
ξ
x
j
+
N
-
3
ξ
x
k
,
y
=
N
-
1
ξ
y
i
+
N
-
2
ξ
y
j
+
N
-
3
ξ
y
k
,with(16)
N
-
1
=
-
0.5
ξ
1
-
ξ
,
N
-
2
=
1
-
ξ
2
,
N
-
3
=
0.5
ξ
1
+
ξ
,where ξ is the local coordinate along each element boundary and -
1
≤
ξ
≤
1.Then, forith edge Γ
i, the matrix N
-
Γ in (4) can be expressed by(17)
2
i
-
1
2
i
2
j
-
1
2
j
2
k
-
1
2
k
∣
∣
∣
∣
∣
∣
N
-
Γ
i
=
0
0
⋯
⋯
N
-
1
0
0
N
-
1
N
-
2
0
0
N
-
2
N
-
3
0
0
N
-
3
⋯
⋯
0
0
2
×
4
n
if
k
>
2
n
,
let
k
=
1
.


### 3.2. Evaluation Procedure for the MatrixH in (11)

For an arbitraryn-sided polygonal element, the evaluation of matrix H in (11) should be performed along the whole element boundaries. Thus, (11) can be rewritten as (18)
H
=
∫
Γ
S
T
L
T
N
-
Γ
d
s
=
∑
i
=
1
n
∫
Γ
i
S
x
,
y
T
L
T
N
-
Γ
i
d
s
.And the direction cosines of the outer normal of each element edge, l and m in (2), are given by(19)
l
=
d
y
d
s
,
m
=
-
d
x
d
s
.Substitution of (15) and (16) into (18) yields(20)
H
=
∑
i
=
1
n
∫
-
1
1
S
ξ
T
L
~
T
N
-
Γ
i
d
ξ
,with(21)
L
~
=
d
y
d
ξ
0
-
d
x
d
ξ
0
-
d
x
d
ξ
d
y
d
ξ
,
d
x
=
d
N
-
1
ξ
x
i
+
d
N
-
2
ξ
x
j
+
d
N
-
3
ξ
x
k
,
d
y
=
d
N
-
1
ξ
y
i
+
d
N
-
2
ξ
y
j
+
d
N
-
3
ξ
y
k
.Five Gauss integration points are used to evaluate (18).


### 3.3. Evaluation Procedure for the MatrixM in (8)

As shown in Figure2, any polygonal element can be divided into some subtriangles. The matrix M in (8) will be firstly evaluated in each subtriangle by standard numerical integration technique. And the sum of these values is the final matrix M. The whole procedure can be performed automatically by computer code.Figure 2
Domain divisions of a dodecagonal element.The resulting new polygonal element is denoted by HSF-AP-19β.


## 3.1. Displacement Interpolation along the Element Boundary

The element boundary displacement vector can be described by the element nodal displacement vector (see (4)). For element edge with three nodes i, j
=
i
+
1, and k
=
j
+
1 (if k
>
2
n, let k
=
1), let(15)
x
=
N
-
1
ξ
x
i
+
N
-
2
ξ
x
j
+
N
-
3
ξ
x
k
,
y
=
N
-
1
ξ
y
i
+
N
-
2
ξ
y
j
+
N
-
3
ξ
y
k
,with(16)
N
-
1
=
-
0.5
ξ
1
-
ξ
,
N
-
2
=
1
-
ξ
2
,
N
-
3
=
0.5
ξ
1
+
ξ
,where ξ is the local coordinate along each element boundary and -
1
≤
ξ
≤
1.Then, forith edge Γ
i, the matrix N
-
Γ in (4) can be expressed by(17)
2
i
-
1
2
i
2
j
-
1
2
j
2
k
-
1
2
k
∣
∣
∣
∣
∣
∣
N
-
Γ
i
=
0
0
⋯
⋯
N
-
1
0
0
N
-
1
N
-
2
0
0
N
-
2
N
-
3
0
0
N
-
3
⋯
⋯
0
0
2
×
4
n
if
k
>
2
n
,
let
k
=
1
.


## 3.2. Evaluation Procedure for the MatrixH in (11)

For an arbitraryn-sided polygonal element, the evaluation of matrix H in (11) should be performed along the whole element boundaries. Thus, (11) can be rewritten as (18)
H
=
∫
Γ
S
T
L
T
N
-
Γ
d
s
=
∑
i
=
1
n
∫
Γ
i
S
x
,
y
T
L
T
N
-
Γ
i
d
s
.And the direction cosines of the outer normal of each element edge, l and m in (2), are given by(19)
l
=
d
y
d
s
,
m
=
-
d
x
d
s
.Substitution of (15) and (16) into (18) yields(20)
H
=
∑
i
=
1
n
∫
-
1
1
S
ξ
T
L
~
T
N
-
Γ
i
d
ξ
,with(21)
L
~
=
d
y
d
ξ
0
-
d
x
d
ξ
0
-
d
x
d
ξ
d
y
d
ξ
,
d
x
=
d
N
-
1
ξ
x
i
+
d
N
-
2
ξ
x
j
+
d
N
-
3
ξ
x
k
,
d
y
=
d
N
-
1
ξ
y
i
+
d
N
-
2
ξ
y
j
+
d
N
-
3
ξ
y
k
.Five Gauss integration points are used to evaluate (18).


## 3.3. Evaluation Procedure for the MatrixM in (8)

As shown in Figure2, any polygonal element can be divided into some subtriangles. The matrix M in (8) will be firstly evaluated in each subtriangle by standard numerical integration technique. And the sum of these values is the final matrix M. The whole procedure can be performed automatically by computer code.Figure 2
Domain divisions of a dodecagonal element.The resulting new polygonal element is denoted by HSF-AP-19β.


## 4. Numerical Examples

In this section, several problems are solved to test the performance of the new quadratic polygonal HSF-AP-19β element. The results are compared with those obtained by the following two elements:PS2: the quadratic polygonal spline element based on area coordinates and B-net method, Chen et al. [32];BFEM: the base force element method for quadratic polygonal mesh problems, Peng et al. [37, 38].


### 4.1. Patch Test

A small patch is divided into some arbitrary polygonal elements in Figure3. The displacement fields and the exact stress solution corresponding to the constant strain are(22)
u
=
1
0
-
3
x
+
y
2
,
v
=
1
0
-
3
y
+
x
2
,
σ
x
=
σ
y
=
1333.3333
,
τ
x
y
=
400.0
.Figure 3
Constant stress/strain patch test.The displacements of boundary nodes are taken as the displacement boundary conditions. No matter whether the inner element edges are straight or curved and no matter whether the shapes of the elements are convex or concave, the exact results of the displacements at each node and stresses at any point can be obtained by the present HSF-AP-15β element. This demonstrates that the new element passes the patch test and thus can ensure solution convergence.


### 4.2. Pure Shear Test of a Rectangular Plate

As shown in Figure4, several polygonal meshes are applied to analyze a cantilever rectangular plate subjected to uniformly distributing shearing forces. Consider the plane stress condition; the parameters of the model are taken as a
=
8, b
=
5, E
=
1, μ
=
0.3, and τ
0
=
1.0. It can be seen that no matter whether the element shapes in the meshes are convex or concave polygonal, the following exact solutions [2],(23)
τ
x
y
=
1
,
v
x
=
a
=
7.8
,can be obtained.Figure 4
Rectangular plate subjected to shearing forces.


### 4.3. Pure Bending Problem of a Cantilever Beam

As shown in Figure5, a cantilever beam under plane stress condition is subjected to a constant bending moment M. The theoretical solutions for this problem are given by [48] (24)
σ
x
=
12
c
y
-
6
,
σ
y
=
0
,
τ
x
y
=
0
,
u
=
-
6
E
x
+
12
c
E
x
y
,
v
=
6
μ
E
y
-
6
c
E
x
2
-
6
μ
c
E
y
2
.The dimensions of the beam are L
=
10 and c
=
1. Other parameters are given in Figure 5. Eighteen mesh divisions containing polygonal elements with concave and convex shapes are employed for the calculations. The displacement and stress results of select points are listed in Table 2.Table 2
Results at selected locations for the pure bending problem (Figure5).
Models
Meshes
σ
x
0,0
σ
x
0,1
v
10,0
×
1
0
-
4
HSF-AP-19β
Mesh 1
−6.000
6.000
−6.0000
Mesh 2
−6.000
6.000
−6.0000
Mesh 3
−6.000
6.000
−6.0000
Mesh 4
−6.000
6.000
−6.0000
Mesh 5
−6.000
6.000
−6.0000
Mesh 6
−6.000
6.000
−6.0000
Mesh 7
−6.000
6.000
−6.0000
Mesh 8
−6.000
6.000
−6.0000
Mesh 9
−6.000
6.000
−6.0000
Mesh 10
−6.000
6.000
−6.0000
Mesh 11
−6.000
6.000
−6.0000
Mesh 12
−6.000
6.000
−6.0000
Mesh 13
−6.000
6.000
−6.0000
Mesh 14
−6.000
6.000
−6.0000
Mesh 15
−6.000
6.000
−6.0000
Mesh 16
−6.000
6.000
−6.0000
Mesh 17
−6.000
6.000
−6.0000
Mesh 18
−6.048
6.089
−5.9972
Mesh 19
—
—
−6.0339
BFEM [37]
Mesh 20
—
—
−6.0359
Mesh 21
—
—
−6.0379
Exact
−6.000
6.000
−6.0000Figure 5
Pure bending problem for a cantilever beam.It can be seen that, so long as all element edges keep straight, exact displacements and stress solutions can always be obtained by the element HSF-AP-19β. However, even if a large number of meshes are employed to solve this problem, the exact solutions cannot be obtained by the BFEM [37].


### 4.4. Linear Bending for a Cantilever Beam

As shown in Figure6, a cantilever under plane stress condition beam is subjected to a linear bending moment caused by a shear force P at the free end. The theoretical solutions for this problem are given by [48] (25)
σ
x
=
240
c
L
x
y
-
120
L
x
-
240
c
y
+
120
,
σ
y
=
0
,
τ
x
y
=
-
120
c
y
2
+
120
L
y
,
u
=
120
c
L
x
2
y
-
92
c
L
y
3
-
60
L
x
2
-
240
c
x
y
m
m
c
+
138
L
y
2
+
120
x
-
46
c
L
y
E
-
1
,
v
=
-
40
c
L
x
3
-
36
c
L
x
y
2
+
120
c
x
2
+
36
L
x
y
m
m
c
+
36
c
y
2
+
46
c
L
x
-
36
y
E
-
1
.The mesh divisions employed are given in Figure 6. All numerical results are listed in Table 3. For this linear bending problem, the exact solutions cannot be directly obtained by the present HSF-AP-19β element. However, once the upper and lower boundaries of the beam are divided into more segments, the results will rapidly converge to the exact solutions.Table 3
Results obtained by the element HSF-AP-19β at selected locations for the linear bending problem (Figure 6).
v
(
100,0
)
×
1
0
-
3
σ
(
0,0
)
σ
(
0,10
)
Numerical results
Normalized results
Numerical results
Normalized results
Numerical results
Normalized results
Mesh 1
6.340
0.788
67.57
0.563
−67.34
0.561
Mesh 2
7.834
0.974
102.16
0.851
−101.94
0.850
Mesh 3
7.915
0.984
112.57
0.938
−113.12
0.943
Mesh 4
7.911
0.983
112.88
0.941
−112.56
0.938
Mesh 5
8.026
0.998
114.18
0.952
−113.66
0.947
Mesh 6
8.035
0.999
115.74
0.965
−116.11
0.968
Mesh 7
8.033
0.998
113.44
0.945
−113.10
0.943
Mesh 8
8.062
1.002
117.59
0.980
−117.53
0.979
Mesh 9
8.055
1.001
115.67
0.964
−118.93
0.991
Mesh 10
8.066
1.002
117.82
0.982
−117.84
0.982
Exact
8.046
120.00
−120.00Figure 6
Linear bending problem for a cantilever beam.


### 4.5. MacNeal’s Thin Cantilever Beam Problem [45]

As shown in Figure7, six new mesh divisions are designed to solve thin beam problem proposed by MacNeal and Harder [45]. Two loading cases are considered: pure bending and transverse linear bending. Young’s modulus E
=
10
7, Poisson’s ratio μ
=
0.3, and the thickness of the beam t
=
0.1. The results of the tip deflection are given in Table 4. It can be seen that the present HSF-AP-19β element exhibits good performance: it can not only provide the exact solutions for the pure bending case but also produce the high precision results for the linear bending case.Table 4
Results at selected locations for MacNeal’s thin beam problem (Figure7).
LoadP
LoadM (×10−02)
Numerical results
Normalized results
Numerical results
Normalized results
Mesh 1
−0.1060
0.9806
−0.5400
1.0000
Mesh 2
−0.1060
0.9806
−0.5400
1.0000
Mesh 3
−0.1060
0.9806
−0.5400
1.0000
Mesh 4
−0.1060
0.9806
−0.5400
1.0000
Mesh 5
−0.1060
0.9806
−0.5400
1.0000
Mesh 6
−0.1064
0.9843
−0.5400
1.0000
[45]
−0.1081
−0.5400Figure 7
MacNeal’s thin beam problem.


### 4.6. Cook’s Skew Beam Problem [49]

As shown in Figure8, a skew cantilever under plane stress condition is subjected to a shear distributed load at the free edge. The typical meshes are given in Figure 8. The results of vertical deflection at point C (48, 52), the maximum principal stress at point A (24, 22), and the minimum principal stress at point B (24, 52) are all listed in Table 5, in which the stresses are the average values at the corresponding nodes. It can be seen that the new HSF-AP-19β element has a good convergence.Table 5
Results at selected locations for Cook’s skew beam problem (Figure8).
Mesh
v
C
σ
A
max
⁡
σ
B
min
⁡
2 × 2
23.79
0.2473
−0.1761
4 × 4
23.94
0.2406
−0.1946
6 × 6
23.95
0.2386
−0.2018
8 × 8
23.96
0.2378
−0.2033
Reference solutiona
23.96
0.2362
−0.2023
aResults of the element GT9M8 [46] using 64 × 64 mesh.Figure 8
Cook’s skew beam problem.


### 4.7. Curved Beam

As shown in Figure9, a curved cantilever beam under plane stress condition is subjected to a transverse force P at the upper end in the radial direction. The dimensions of the model are as follows: the inner radius a
=
0.7, the outer radius b
=
1.0, and other constants are given in Figure 9. Two Poisson’s ratio cases are considered: μ
=
0.3 and μ
=
0.4999. Furthermore, two mesh division types are adopted: elements with curved edges and straight edges. The theoretical radial displacement of point A at the upper end is given by [48] (26)
v
A
=
-
P
π
(
a
2
+
b
2
)
E
[
(
a
2
-
b
2
)
+
(
a
2
+
b
2
)
log
⁡
(
b
/
a
)
]
.Figure 9
Curved beam problem.The results of the tip displacementv
A are given in Table 6. It can be observed that the new HSF-AP-19β element exhibits good convergence. Moreover, the errors produced by the HSF-AP-19β element are quit small, no matter whether the meshes have curved or straight edges.Table 6
The tip displacementsv
A of the curved beam (Figure 9).
Curved edges
Straight edges
μ
=
0.3
μ
=
0.4999
μ
=
0.3
μ
=
0.4999
1 × 2
−0.0211
−0.0211
−0.0223
−0.0222
1 × 4
−0.0217
−0.0216
−0.0218
−0.0217
1 × 6
−0.0217
−0.0216
−0.0217
−0.0216
1 × 8
−0.0217
−0.0216
−0.0217
−0.0216
1 × 10
−0.0217
−0.0217
−0.0217
−0.0217
Analytical solution
−0.0218


### 4.8. Infinite Plate with a Circular Hole [3]

As shown in Figure10, a plate with a central circular hole is under plane stress state and is subjected to a uniform tensile load of 1.0 in x-direction. In order to simulate the properties of the infinite plate, the plate-to-hole aspect ratio of 50 is applied, and the radius of hole is taken as 1.0, while the half-width of plate is 50. The related parameters are as follows: E
=
1.0
×
10
3, μ
=
0.3, and thickness t
=
1.0. The meshes employed are given in Figure 10.Figure 10
Infinite plate with a circular hole.The results at pointsA and B are listed in Table 7. Besides the HSF-AP-19β element, the results obtained by the PS2 [32] element are also given for comparison. It can be seen that the accuracy of the results obtained by the HSF-AP-19β element is higher than those obtained by the PS2 element, especially for stress solutions.Table 7
The results at pointsA and B of the infinite plate with a circular hole (Figure 10).
Mesh
u
A
/
1
0
-
3
v
B
/
1
0
-
3
σ
y
A
σ
x
B
Number of nodes
HSF-AP-19β
4 × 4
2.992
−0.999
−1.058
3.045
65
8 × 8
3.002
−1.002
−1.007
3.008
225
12 × 12
3.003
−1.003
−1.004
3.006
481
16 × 16
3.003
−1.002
−1.003
3.003
883
20 × 20
3.003
−1.002
−1.002
3.002
1281
4 × 4
2.917
−0.950
−1.025
3.095
65
PS2 [32]
8 × 8
2.985
−0.993
−1.026
3.044
225
12 × 12
2.994
−0.998
−1.013
3.021
481
Exact
3.000
−1.000
−1.000
3.000


### 4.9. A Cantilever Beam under a Concentrated Load

As shown in Figure11, a cantilever beam with large aspect ratio is subjected to a concentrated load P
=
1 at free end. The geometric parameters of the cantilever beam are as follows: L
=
5, t
=
1, and h
=
0.1. Young’s modulus is E
=
10
9, and Poisson’s ratio is zero. Six mesh divisions are shown in Figure 11, in which Meshes a, b, and c are employed for the HSF-AP-19β element and Meshes d, e, and f for the BFEM [37].Figure 11
A cantilever beam under a concentrated load.The vertical displacements of pointA are presented in Table 8. It can be seen that the HSF-AP-19β element can provide better results by using much less elements.Table 8
The vertical displacements at pointA of the cantilever beam under a concentrated load (Figure 11).
Mesh
HSF-AP-19β
BFEM
Analytical solution
Numerical results
Relative errors
Numerical results
Relative errors
a
−4.7915
4.1696%
b
−4.9561
0.8772%
—
c
−4.9789
0.4230%
−5.0000
d
−4.9787
0.4260%
e
—
−5.0159
0.3180%
f
−4.9734
0.5320%


### 4.10. Cook’s Skew Beam Problem with Heterogeneous Material

As shown in Figure12, a Cook’s skew beam with heterogeneous material (under plane stress condition) is subjected to a shear distributed load P at the free edge. The beam is heterogeneous and made of two materials: (i) E
=
1.0, μ
=
0.25 and (ii) E
=
1.5, μ
=
0.3. Two heterogeneous models are considered here, and the corresponding mesh divisions are given in Figures 13(a) and 13(b), in which the gray domains are made of the first material, while other domains are made of the second material. The results of vertical deflection at point C (48, 52) are listed in Table 9.Table 9
The vertical displacements at pointC of Cook’s skew beam problem with heterogeneous materials (Figures 12 and 13).
Heterogeneous model
Mesh
Number of elements/nodes
HSF-AP-19β
CPS8 [47]
CPS8R [47]
The first case (a)
Mesh 1
9/48
18.7076
—
—
Mesh 2
15/64
—
18.6115
18.5847
Mesh 3
50/177
18.6890
18.7002
Mesh 4
77/264
18.6957
18.7000
Mesh 5
112/383
18.7227
18.7200
Mesh 6
424/1365
18.7408
18.7422
The second case (b)
Mesh 1
27/128
19.4072
—
—
Mesh 2
112/382
—
19.3741
19.3846
Mesh 3
245/802
19.3992
19.4037
Mesh 4
424/1365
19.4062
19.4080
Mesh 5
764/2413
19.4103
19.4114Figure 12
Cook’s skew beam problem with heterogeneous material.Figure 13
Meshes divisions for the heterogeneous models.
(a)
Meshes for the first heterogeneous model
(b)
Meshes for the second heterogeneous modelWith the mesh refinement, the results obtained by Abaqus 8-node elements CPS8 and CPS8R [47] will tend to the exact solutions. As to a specified precision, the total element number needed by Abaqus is much more than that required by the HSF-AP-19β element. This example is a very simple problem of heterogeneous materials. However, when solids with more complicated materials are modeled, such as the crystal structures, the computational cost of the HSF-AP-19β element will be much less than that of CPS8 and CPS8R in Abaqus. Therefore, it can be said that the new polygonal element exhibits obvious advantages in the analysis of complex heterogeneous material problems.


## 4.1. Patch Test

A small patch is divided into some arbitrary polygonal elements in Figure3. The displacement fields and the exact stress solution corresponding to the constant strain are(22)
u
=
1
0
-
3
x
+
y
2
,
v
=
1
0
-
3
y
+
x
2
,
σ
x
=
σ
y
=
1333.3333
,
τ
x
y
=
400.0
.Figure 3
Constant stress/strain patch test.The displacements of boundary nodes are taken as the displacement boundary conditions. No matter whether the inner element edges are straight or curved and no matter whether the shapes of the elements are convex or concave, the exact results of the displacements at each node and stresses at any point can be obtained by the present HSF-AP-15β element. This demonstrates that the new element passes the patch test and thus can ensure solution convergence.


## 4.2. Pure Shear Test of a Rectangular Plate

As shown in Figure4, several polygonal meshes are applied to analyze a cantilever rectangular plate subjected to uniformly distributing shearing forces. Consider the plane stress condition; the parameters of the model are taken as a
=
8, b
=
5, E
=
1, μ
=
0.3, and τ
0
=
1.0. It can be seen that no matter whether the element shapes in the meshes are convex or concave polygonal, the following exact solutions [2],(23)
τ
x
y
=
1
,
v
x
=
a
=
7.8
,can be obtained.Figure 4
Rectangular plate subjected to shearing forces.


## 4.3. Pure Bending Problem of a Cantilever Beam

As shown in Figure5, a cantilever beam under plane stress condition is subjected to a constant bending moment M. The theoretical solutions for this problem are given by [48] (24)
σ
x
=
12
c
y
-
6
,
σ
y
=
0
,
τ
x
y
=
0
,
u
=
-
6
E
x
+
12
c
E
x
y
,
v
=
6
μ
E
y
-
6
c
E
x
2
-
6
μ
c
E
y
2
.The dimensions of the beam are L
=
10 and c
=
1. Other parameters are given in Figure 5. Eighteen mesh divisions containing polygonal elements with concave and convex shapes are employed for the calculations. The displacement and stress results of select points are listed in Table 2.Table 2
Results at selected locations for the pure bending problem (Figure5).
Models
Meshes
σ
x
0,0
σ
x
0,1
v
10,0
×
1
0
-
4
HSF-AP-19β
Mesh 1
−6.000
6.000
−6.0000
Mesh 2
−6.000
6.000
−6.0000
Mesh 3
−6.000
6.000
−6.0000
Mesh 4
−6.000
6.000
−6.0000
Mesh 5
−6.000
6.000
−6.0000
Mesh 6
−6.000
6.000
−6.0000
Mesh 7
−6.000
6.000
−6.0000
Mesh 8
−6.000
6.000
−6.0000
Mesh 9
−6.000
6.000
−6.0000
Mesh 10
−6.000
6.000
−6.0000
Mesh 11
−6.000
6.000
−6.0000
Mesh 12
−6.000
6.000
−6.0000
Mesh 13
−6.000
6.000
−6.0000
Mesh 14
−6.000
6.000
−6.0000
Mesh 15
−6.000
6.000
−6.0000
Mesh 16
−6.000
6.000
−6.0000
Mesh 17
−6.000
6.000
−6.0000
Mesh 18
−6.048
6.089
−5.9972
Mesh 19
—
—
−6.0339
BFEM [37]
Mesh 20
—
—
−6.0359
Mesh 21
—
—
−6.0379
Exact
−6.000
6.000
−6.0000Figure 5
Pure bending problem for a cantilever beam.It can be seen that, so long as all element edges keep straight, exact displacements and stress solutions can always be obtained by the element HSF-AP-19β. However, even if a large number of meshes are employed to solve this problem, the exact solutions cannot be obtained by the BFEM [37].


## 4.4. Linear Bending for a Cantilever Beam

As shown in Figure6, a cantilever under plane stress condition beam is subjected to a linear bending moment caused by a shear force P at the free end. The theoretical solutions for this problem are given by [48] (25)
σ
x
=
240
c
L
x
y
-
120
L
x
-
240
c
y
+
120
,
σ
y
=
0
,
τ
x
y
=
-
120
c
y
2
+
120
L
y
,
u
=
120
c
L
x
2
y
-
92
c
L
y
3
-
60
L
x
2
-
240
c
x
y
m
m
c
+
138
L
y
2
+
120
x
-
46
c
L
y
E
-
1
,
v
=
-
40
c
L
x
3
-
36
c
L
x
y
2
+
120
c
x
2
+
36
L
x
y
m
m
c
+
36
c
y
2
+
46
c
L
x
-
36
y
E
-
1
.The mesh divisions employed are given in Figure 6. All numerical results are listed in Table 3. For this linear bending problem, the exact solutions cannot be directly obtained by the present HSF-AP-19β element. However, once the upper and lower boundaries of the beam are divided into more segments, the results will rapidly converge to the exact solutions.Table 3
Results obtained by the element HSF-AP-19β at selected locations for the linear bending problem (Figure 6).
v
(
100,0
)
×
1
0
-
3
σ
(
0,0
)
σ
(
0,10
)
Numerical results
Normalized results
Numerical results
Normalized results
Numerical results
Normalized results
Mesh 1
6.340
0.788
67.57
0.563
−67.34
0.561
Mesh 2
7.834
0.974
102.16
0.851
−101.94
0.850
Mesh 3
7.915
0.984
112.57
0.938
−113.12
0.943
Mesh 4
7.911
0.983
112.88
0.941
−112.56
0.938
Mesh 5
8.026
0.998
114.18
0.952
−113.66
0.947
Mesh 6
8.035
0.999
115.74
0.965
−116.11
0.968
Mesh 7
8.033
0.998
113.44
0.945
−113.10
0.943
Mesh 8
8.062
1.002
117.59
0.980
−117.53
0.979
Mesh 9
8.055
1.001
115.67
0.964
−118.93
0.991
Mesh 10
8.066
1.002
117.82
0.982
−117.84
0.982
Exact
8.046
120.00
−120.00Figure 6
Linear bending problem for a cantilever beam.


## 4.5. MacNeal’s Thin Cantilever Beam Problem [45]

As shown in Figure7, six new mesh divisions are designed to solve thin beam problem proposed by MacNeal and Harder [45]. Two loading cases are considered: pure bending and transverse linear bending. Young’s modulus E
=
10
7, Poisson’s ratio μ
=
0.3, and the thickness of the beam t
=
0.1. The results of the tip deflection are given in Table 4. It can be seen that the present HSF-AP-19β element exhibits good performance: it can not only provide the exact solutions for the pure bending case but also produce the high precision results for the linear bending case.Table 4
Results at selected locations for MacNeal’s thin beam problem (Figure7).
LoadP
LoadM (×10−02)
Numerical results
Normalized results
Numerical results
Normalized results
Mesh 1
−0.1060
0.9806
−0.5400
1.0000
Mesh 2
−0.1060
0.9806
−0.5400
1.0000
Mesh 3
−0.1060
0.9806
−0.5400
1.0000
Mesh 4
−0.1060
0.9806
−0.5400
1.0000
Mesh 5
−0.1060
0.9806
−0.5400
1.0000
Mesh 6
−0.1064
0.9843
−0.5400
1.0000
[45]
−0.1081
−0.5400Figure 7
MacNeal’s thin beam problem.


## 4.6. Cook’s Skew Beam Problem [49]

As shown in Figure8, a skew cantilever under plane stress condition is subjected to a shear distributed load at the free edge. The typical meshes are given in Figure 8. The results of vertical deflection at point C (48, 52), the maximum principal stress at point A (24, 22), and the minimum principal stress at point B (24, 52) are all listed in Table 5, in which the stresses are the average values at the corresponding nodes. It can be seen that the new HSF-AP-19β element has a good convergence.Table 5
Results at selected locations for Cook’s skew beam problem (Figure8).
Mesh
v
C
σ
A
max
⁡
σ
B
min
⁡
2 × 2
23.79
0.2473
−0.1761
4 × 4
23.94
0.2406
−0.1946
6 × 6
23.95
0.2386
−0.2018
8 × 8
23.96
0.2378
−0.2033
Reference solutiona
23.96
0.2362
−0.2023
aResults of the element GT9M8 [46] using 64 × 64 mesh.Figure 8
Cook’s skew beam problem.


## 4.7. Curved Beam

As shown in Figure9, a curved cantilever beam under plane stress condition is subjected to a transverse force P at the upper end in the radial direction. The dimensions of the model are as follows: the inner radius a
=
0.7, the outer radius b
=
1.0, and other constants are given in Figure 9. Two Poisson’s ratio cases are considered: μ
=
0.3 and μ
=
0.4999. Furthermore, two mesh division types are adopted: elements with curved edges and straight edges. The theoretical radial displacement of point A at the upper end is given by [48] (26)
v
A
=
-
P
π
(
a
2
+
b
2
)
E
[
(
a
2
-
b
2
)
+
(
a
2
+
b
2
)
log
⁡
(
b
/
a
)
]
.Figure 9
Curved beam problem.The results of the tip displacementv
A are given in Table 6. It can be observed that the new HSF-AP-19β element exhibits good convergence. Moreover, the errors produced by the HSF-AP-19β element are quit small, no matter whether the meshes have curved or straight edges.Table 6
The tip displacementsv
A of the curved beam (Figure 9).
Curved edges
Straight edges
μ
=
0.3
μ
=
0.4999
μ
=
0.3
μ
=
0.4999
1 × 2
−0.0211
−0.0211
−0.0223
−0.0222
1 × 4
−0.0217
−0.0216
−0.0218
−0.0217
1 × 6
−0.0217
−0.0216
−0.0217
−0.0216
1 × 8
−0.0217
−0.0216
−0.0217
−0.0216
1 × 10
−0.0217
−0.0217
−0.0217
−0.0217
Analytical solution
−0.0218


## 4.8. Infinite Plate with a Circular Hole [3]

As shown in Figure10, a plate with a central circular hole is under plane stress state and is subjected to a uniform tensile load of 1.0 in x-direction. In order to simulate the properties of the infinite plate, the plate-to-hole aspect ratio of 50 is applied, and the radius of hole is taken as 1.0, while the half-width of plate is 50. The related parameters are as follows: E
=
1.0
×
10
3, μ
=
0.3, and thickness t
=
1.0. The meshes employed are given in Figure 10.Figure 10
Infinite plate with a circular hole.The results at pointsA and B are listed in Table 7. Besides the HSF-AP-19β element, the results obtained by the PS2 [32] element are also given for comparison. It can be seen that the accuracy of the results obtained by the HSF-AP-19β element is higher than those obtained by the PS2 element, especially for stress solutions.Table 7
The results at pointsA and B of the infinite plate with a circular hole (Figure 10).
Mesh
u
A
/
1
0
-
3
v
B
/
1
0
-
3
σ
y
A
σ
x
B
Number of nodes
HSF-AP-19β
4 × 4
2.992
−0.999
−1.058
3.045
65
8 × 8
3.002
−1.002
−1.007
3.008
225
12 × 12
3.003
−1.003
−1.004
3.006
481
16 × 16
3.003
−1.002
−1.003
3.003
883
20 × 20
3.003
−1.002
−1.002
3.002
1281
4 × 4
2.917
−0.950
−1.025
3.095
65
PS2 [32]
8 × 8
2.985
−0.993
−1.026
3.044
225
12 × 12
2.994
−0.998
−1.013
3.021
481
Exact
3.000
−1.000
−1.000
3.000


## 4.9. A Cantilever Beam under a Concentrated Load

As shown in Figure11, a cantilever beam with large aspect ratio is subjected to a concentrated load P
=
1 at free end. The geometric parameters of the cantilever beam are as follows: L
=
5, t
=
1, and h
=
0.1. Young’s modulus is E
=
10
9, and Poisson’s ratio is zero. Six mesh divisions are shown in Figure 11, in which Meshes a, b, and c are employed for the HSF-AP-19β element and Meshes d, e, and f for the BFEM [37].Figure 11
A cantilever beam under a concentrated load.The vertical displacements of pointA are presented in Table 8. It can be seen that the HSF-AP-19β element can provide better results by using much less elements.Table 8
The vertical displacements at pointA of the cantilever beam under a concentrated load (Figure 11).
Mesh
HSF-AP-19β
BFEM
Analytical solution
Numerical results
Relative errors
Numerical results
Relative errors
a
−4.7915
4.1696%
b
−4.9561
0.8772%
—
c
−4.9789
0.4230%
−5.0000
d
−4.9787
0.4260%
e
—
−5.0159
0.3180%
f
−4.9734
0.5320%


## 4.10. Cook’s Skew Beam Problem with Heterogeneous Material

As shown in Figure12, a Cook’s skew beam with heterogeneous material (under plane stress condition) is subjected to a shear distributed load P at the free edge. The beam is heterogeneous and made of two materials: (i) E
=
1.0, μ
=
0.25 and (ii) E
=
1.5, μ
=
0.3. Two heterogeneous models are considered here, and the corresponding mesh divisions are given in Figures 13(a) and 13(b), in which the gray domains are made of the first material, while other domains are made of the second material. The results of vertical deflection at point C (48, 52) are listed in Table 9.Table 9
The vertical displacements at pointC of Cook’s skew beam problem with heterogeneous materials (Figures 12 and 13).
Heterogeneous model
Mesh
Number of elements/nodes
HSF-AP-19β
CPS8 [47]
CPS8R [47]
The first case (a)
Mesh 1
9/48
18.7076
—
—
Mesh 2
15/64
—
18.6115
18.5847
Mesh 3
50/177
18.6890
18.7002
Mesh 4
77/264
18.6957
18.7000
Mesh 5
112/383
18.7227
18.7200
Mesh 6
424/1365
18.7408
18.7422
The second case (b)
Mesh 1
27/128
19.4072
—
—
Mesh 2
112/382
—
19.3741
19.3846
Mesh 3
245/802
19.3992
19.4037
Mesh 4
424/1365
19.4062
19.4080
Mesh 5
764/2413
19.4103
19.4114Figure 12
Cook’s skew beam problem with heterogeneous material.Figure 13
Meshes divisions for the heterogeneous models.
(a)
Meshes for the first heterogeneous model
(b)
Meshes for the second heterogeneous modelWith the mesh refinement, the results obtained by Abaqus 8-node elements CPS8 and CPS8R [47] will tend to the exact solutions. As to a specified precision, the total element number needed by Abaqus is much more than that required by the HSF-AP-19β element. This example is a very simple problem of heterogeneous materials. However, when solids with more complicated materials are modeled, such as the crystal structures, the computational cost of the HSF-AP-19β element will be much less than that of CPS8 and CPS8R in Abaqus. Therefore, it can be said that the new polygonal element exhibits obvious advantages in the analysis of complex heterogeneous material problems.


## 5. Concluding Remarks

In the paper, a novel quadratic hybrid stress-function polygonal element model, HSF-AP-19β, is developed based on the principle of minimum complementary energy and the fundamental analytical solutions of the Airy stress function. The present element has the following characters.(1)
It is simple to build this new element without construction of displacement interpolation functions, which are very complex for polygonal elements compared with quadrilateral elements. This polygonal HSF-AP-19β element can be degenerated into triangular or quadrilateral elements directly.
(2)
Spurious zero energy modes can be avoided by selecting suitable number of fundamental analytical solutions for stresses (from constant terms to higher-order terms). For quadratic HS-F model, the first nineteen stress solutions must be employed.
(3)
The present HSF-AP-19β element exhibits better accuracy and convergence in both displacement and stress solutions. Even though meshes containing convex or concave elements are employed, the precisions of the results will not be affected. So, it is a shape-free element model. For constant stress/strain problem, the exact solutions can be obtained for all cases. So long as element edges keep straight, all elements can produce the exact solutions for pure bending problem. Moreover, the new element can produce better results using relatively coarse meshes for some problems, such as high-order problems, problems with circular holes, and heterogeneous material problems.



---

*Source: 491325-2015-05-19.xml*