# A New Method Applied to the Quadrilateral Membrane Element with Vertex Rigid Rotational Freedom

Authors: Xiaowei Gao; Yunfei Liu; Jun Lv
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

In order to improve the performance of the membrane element with vertex rigid rotational freedom, a new method to establish the local Cartesian coordinate system and calculate the derivatives of the shape functions with respect to the local coordinates is introduced in this paper. The membrane elements with vertex rigid rotational freedom such as GQ12 and GQ12M based on this new method can achieve higher precision results than traditional methods. The numerical results demonstrate that the elements GQ12 and GQ12M with this new method can provide better membrane elements for flat shell elements. Furthermore, this new method presented in this paper offers a new approach for other membrane elements used in flat shell element to improve the computing accuracy.


---

## Body

## 1. Introduction

The finite element method has been used for solving problems in different fields of engineering. For now, a large number of different finite elements have been developed and the shell element is one of these elements to solve the multiscale problems [1, 2]. The shell element can be categorized into two types: curved shell element and flat shell element [3]. The curved shell element can be used to make a good description of the geometrical shape of the shell structure, so it is with a good calculation precision. But the formula expressions of the curved shell element are very complicated so it is limited in practical application. In this case, the flat shell element is becoming more attractive in deriving efficient numerical accuracy and its concise theory due to its superposition theory of the membrane element and plate bending element [4].The membrane elements are among the simplest elements to develop, which are used for analyzing structures subjected to in-plane forces. The membrane elements are usually used to model the behavior of shear wall, stiffened sheet construction, and membrane action in shells. Some plane elements can be considered as membrane elements, such as the CST (constant strain triangle) element and the four-node isoparametric quadrilateral plane element (QUAD4) [5]. In finite element methods, many plate bending elements also have been developed. Bazeley et al. [6] developed the confirming and nonconfirming plate bending elements in 1966. Batoz et al. [7] developed three types of plate bending elements, that is, the DKT element, the HSM element, and the SRI element for the analysis of plates and shells in 1980. In 1982, the quadrilateral plate bending element (DKQ) was formulated by Batoz and Tahar [8] based on the DKT element. Other plate bending elements have been developed in the following years [9–12].The flat shell elements are developed by combining membrane elements containing two in-plane translational degrees of freedom and plate bending elements containing two rotational degrees of freedom and one out of plane translational degree of freedom. Since the in-plane rotational degrees of freedom are not included, the null values for the in-plane rotational degrees of freedom will lead to singularity in structure stiffness matrix if all the elements are coplanar. The simplest method adopted to remove the singularity is to add the in-plane rotational degrees of freedom to the membrane elements, which can also improve their performance [5]. After several triangular elements [13, 14] with vertex rigid rotational freedom are proposed, the quadrilateral membrane elements including corner rotations have been developed by Allman [15] and MacNeal and Harder [16] as well. And then, numerous researches focused on membrane elements with drilling degrees of freedom have been accomplished, such as the models proposed by Iura and Atluri [17], Cazzani and Atluri [18], Piltner and Taylor [19], Geyer and Groenwold [20], Groenwold et al. [21], Choi et al. [22], Kugler et al. [23], and Cen et al. [24]. In 1994, two new arbitrary quadrilateral membrane elements called GQ12 and GQ12M with vertex rotation were proposed by Long and Xu [25], resulting in more reasonable compatible conditions between adjoining elements and a more simple formulation.In order to further improve the computing accuracy of the elements GQ12 and GQ12M, a new method for establishing the local Cartesian coordinate system and calculating the derivatives of the shape function with respect to the local coordinates is presented in this paper. This new method can be applied to the membrane elements which may provide a constituent part for flat shell elements. The elements GQ12 and GQ12M have more accurate numerical results compared to the bilinear quadrilateral element Q4 [3, 25]. In this paper, the applications of the new method to the elements GQ12 and GQ12M were proposed to examine the performance of this new method described in this paper.


## 2. The Membrane Elements GQ12 and GQ12M

A quadrilateral membrane element with rotational degree of freedom (shown in Figure1) was proposed in [25].Figure 1
The quadrilateral membrane element.The nodal displacement vector over this element is given by(1)
q
e
=
q
1
T
q
2
T
q
3
T
q
4
T
T
.The freedoms at each node are(2)
q
i
=
u
i
v
i
θ
i
T
,
i
=
1,2
,
3,4
,in which u
i and v
i are the translations and θ
i is the additional rotation at each corner. A membrane element is usually located in a three-dimensional space, so its orientation may be in any directions. For convenience, the membrane element is commonly studied in the local Cartesian coordinate system represented by x
′
,
y
′
,
z
′ and located on the membrane element. Once a quantity is formed in the local system, it can be transformed to the global coordinate system.In the local system, the element GQ12 is formulated by the displacement field which includes two parts as follows [25]:(3)
u
=
u
0
+
u
θ
,in which u
0 is the conventional bilinear compatible displacement decided by the translations u
i
v
i
T at the nodes and u
θ is the additional displacement determined by the rigid rotation θ
i
i
=
1,2
,
3,4 at each node. The element GQ12 displacement can be described in the form of(4)
u
=
u
0
+
u
θ
=
N
q
e
=
∑
i
=
1
4
N
i
q
i
,where N is the shape function matrix of element GQ12, which can be given in terms of (5)
N
i
=
N
i
0
0
N
u
θ
i
0
N
i
0
N
v
θ
i
,
i
=
1,2
,
3,4
,in which(6)
N
i
0
=
1
4
1
+
ξ
i
ξ
1
+
η
i
η
,
i
=
1,2
,
3,4
,
(7)
N
u
θ
i
=
1
8
ξ
i
1
-
ξ
2
b
1
+
b
3
η
i
1
+
η
i
η
+
η
i
1
-
η
2
b
2
+
b
3
ξ
i
1
+
ξ
i
ξ
,
N
v
θ
i
=
-
1
8
ξ
i
1
-
ξ
2
a
1
+
a
3
η
i
1
+
η
i
η
+
η
i
1
-
η
2
a
2
+
a
3
ξ
i
1
+
ξ
i
ξ
,
i
=
1,2
,
3,4
,in which(8)
a
1
=
1
4
∑
i
=
1
4
ξ
i
x
i
′
,
a
2
=
1
4
∑
i
=
1
4
η
i
x
i
′
,
a
3
=
1
4
∑
i
=
1
4
ξ
i
η
i
x
i
′
,
b
1
=
1
4
∑
i
=
1
4
ξ
i
y
i
′
,
b
2
=
1
4
∑
i
=
1
4
η
i
y
i
′
,
b
3
=
1
4
∑
i
=
1
4
ξ
i
η
i
y
′
,where ξ
i, η
i and x
i
′, y
i
′
i
=
1,2
,
3,4 denote the intrinsic coordinates and the local Cartesian coordinates of the element nodes, respectively.In the local system, the element stiffness matrix of GQ12 can be written as [25](9)
K
G
Q
12
′
e
=
∬
A
e
B
T
D
B
t
d
A
=
∬
-
1
1
B
T
D
B
J
t
d
ξ
d
η
,where t is the thickness of the element, J is the determinant of the Jacobian matrix, and B is the element strain matrix which can be expressed as(10)
B
=
B
1
B
2
B
3
B
4
,
(11)
B
i
=
∂
∂
x
′
0
0
∂
∂
y
′
∂
∂
y
′
∂
∂
x
′
N
i
0
0
N
u
θ
i
0
N
i
0
N
v
θ
i
=
∂
N
i
0
∂
x
′
0
∂
N
u
θ
i
∂
x
′
0
∂
N
i
0
∂
y
′
∂
N
v
θ
i
∂
y
′
∂
N
i
0
∂
y
′
∂
N
i
0
∂
x
′
∂
N
u
θ
i
∂
y
′
+
∂
N
v
θ
i
∂
x
′
,
i
=
1,2
,
3,4
.In (9), D is the elasticity matrix and it can be expressed as [4](12)
D
=
E
1
-
μ
2
1
μ
0
μ
1
0
0
0
1
-
μ
2
,in which E is the elastic modulus and μ is Poisson’s ratio.In order to further improve the precision of the calculation, the element GQ12M based on the element GQ12 was presented in [25] by Long and Xu. The element GQ12M is formulated with an assumed displacement given by adding a bubble displacement u
^ to GQ12, which is(13)
u
=
u
0
+
u
θ
+
u
^in which u
^ is the bubble displacement. It can be assumed to be(14)
u
^
=
u
^
v
^
=
N
^
ρ
ρ
=
N
^
0
0
N
^
ρ
1
ρ
2
,where(15)
N
^
=
1
-
ξ
2
1
-
η
2
.It can be considered as the shape function at the center node C of the element, as shown in Figure 2, which will be eliminated in the future. ρ
1 and ρ
2 in (14) are arbitrary independent parameters. Substituting (4) and (14) into (13), the element GQ12 displacement can be described in the form of(16)
u
=
u
0
+
u
θ
+
u
^
=
N
q
e
+
N
^
ρ
ρ
.The corresponding strain fields can be expressed as(17)
ε
=
B
q
e
+
B
^
ρin which B is the same as (10), and B
^ is the strain matrix of the bubble displacement u
^, which can be expressed as(18)
B
^
=
∂
∂
x
′
0
0
∂
∂
y
′
∂
∂
y
′
∂
∂
x
′
N
^
0
0
N
^
=
∂
N
^
∂
x
′
0
0
∂
N
^
∂
y
′
∂
N
^
∂
y
′
∂
N
^
∂
x
′
.Figure 2
The adding bubble node C.According to (12) and (17), the strain energy of element GQ12M can be written as(19)
U
=
t
2
∬
ε
T
D
ε
d
A
=
1
2
q
e
T
K
q
q
′
e
q
e
+
ρ
T
K
ρ
q
′
e
q
e
+
1
2
ρ
T
K
ρ
ρ
′
e
ρ
,where(20)
K
q
q
′
e
=
∬
-
1
1
B
T
D
B
J
t
d
ξ
d
η
,
K
ρ
q
′
e
=
∬
-
1
1
B
^
T
D
B
J
t
d
ξ
d
η
,
K
ρ
ρ
′
e
=
∬
-
1
1
B
^
T
D
B
^
J
t
d
ξ
d
ηin which t is the thickness of the element, J is the determinant of the Jacobian matrix, and the matrices B, D, and B
^ are given by (10), (12), and (18), respectively.From the stationary condition∂
U
/
∂
ρ
=
0, the arbitrary parameters ρ can be expressed in terms of the q
e as(21)
ρ
=
-
K
ρ
ρ
′
e
-
1
K
ρ
q
′
e
q
e
.Substituting (21) into (16), the shape functions of the element GQ12M can be obtained as (22)
N
∗
=
N
-
N
^
ρ
K
ρ
ρ
′
e
-
1
K
ρ
q
′
ein which N, N
^
ρ, K
ρ
ρ
′
e
-
1, and K
ρ
q
′
e are given by (5), (14), and (20), respectively.According to the shape functions of GQ12M expressed in (22), the element stiffness matrix of element GQ12M can be written as(23)
K
G
Q
12
M
′
e
=
K
q
q
′
e
-
K
ρ
q
′
e
T
K
ρ
ρ
′
e
-
1
K
ρ
q
′
e
.


## 3. Calculating the Stiffness Matrix of the Membrane Element under the Global Coordinate System in the Traditional Manner

The element stiffness matrices of GQ12 and GQ12M under the local coordinate system are described in the above section, and these matrices can be numerically calculated by the Gaussian integral method. It can be seen from (9) that the GQ12 element stiffness matrix is(24)
K
G
Q
12
′
e
=
∬
-
1
1
B
T
D
B
J
t
d
ξ
d
η
=
∑
i
=
1
N
G
∑
j
=
1
N
G
W
i
W
j
B
T
D
B
J
t
,in which W
i and W
j are the weight functions, and NG is the integration order. Generally, the element stiffness matrix of a 4-node quadrilateral element can be obtained using 2 × 2 Gauss quadrature which can achieve the exact solution of the integral equation. The roots and weight functions for 2 × 2 Gauss quadrature are given in Table 1 [3].Table 1
Roots and weight functions for 2 × 2 Gauss quadrature.
Roots
Weight functionsW
i (W
j)
±0.577350269189626
1.0It can be seen from (11) and (24) that the shape functions with respect to the local coordinates ∂
N
i
/
∂
x
′ and ∂
N
i
/
∂
y
′ and the determinant of the Jacobian matrix J need to be firstly determined to calculate the GQ12 element stiffness matrices. The conventional method [3, 4] to compute ∂
N
i
/
∂
x
′ and ∂
N
i
/
∂
y
′ is to calculate the product of the inverse Jacobi matrix and the derivatives matrix of the shape functions with respect to the intrinsic coordinates, which can be expressed as(25)
∂
N
i
∂
ξ
∂
N
i
∂
η
=
∂
x
′
∂
ξ
∂
y
′
∂
ξ
∂
x
′
∂
η
∂
y
′
∂
η
∂
N
i
∂
x
′
∂
N
i
∂
y
′
=
J
∂
N
i
∂
x
′
∂
N
i
∂
y
′and then(26)
∂
N
i
∂
x
′
∂
N
i
∂
y
′
=
J
-
1
∂
N
i
∂
ξ
∂
N
i
∂
η
,
J
=
∂
x
′
∂
ξ
∂
y
′
∂
ξ
∂
x
′
∂
η
∂
y
′
∂
η
.It is worth pointing out thatN
i in (25) and (26) is a substitution, which may be N
i
0, N
u
θ
i, N
v
θ
i, and N
^ in (6), (7), and (15), respectively.Once the element local stiffness matrix for membrane element is obtained, it is necessary to transform it from the local coordinate system to the global coordinate system. So the element local coordinate system needs to be established first. In the current application of the membrane elements, a traditional manner has been adopted to establish the local Cartesian coordinate system over an element [5].Figure3 shows the global and local coordinate axes for the quadrilateral membrane element. The midpoints of sides 1-2, 2-3, 3-4, and 4-1 are represented by i, j, k, and l, respectively, which can be determined by the shape functions of the four-node isoparametric element. The element local plane is defined by creating two vectors intersecting each other and passing through the midpoints of the sides 2-3 and 3-4 of the quadrilateral as shown in Figure 3.Figure 3
The local Cartesian coordinate system in a traditional way.Assuming that the localx
′ axis of the quadrilateral is parallel to the vector passing through nodes l and j, the vector passing through these nodes can be given by [5](27)
V
x
′
=
V
l
j
=
x
j
-
x
l
y
j
-
y
l
z
j
-
z
l
=
x
j
l
y
j
l
z
j
l
,where x
j, y
j, and z
j and x
l, y
l, and z
l represent the global coordinates of nodes j and l, respectively. The direction cosine λ
x
′ for the local x
′ direction can be obtained by normalizing the vector with respect to its length, which is(28)
λ
x
′
=
1
l
j
l
x
j
l
y
j
l
z
j
lin which l
j
l
=
x
j
l
2
+
y
j
l
2
+
z
j
l
2 is the length of the vector.A reference vectorV
i
k which can define the plane of the element with the vector V
x
′ is obtained by creating a vector passing through nodes i and k of the element as shown in Figure 3
(29)
V
i
k
=
x
k
-
x
i
y
k
-
y
i
z
k
-
z
i
=
x
k
i
y
k
i
z
k
i
,where x
k, y
k, and z
k and x
i, y
i, and z
i represent the global coordinates of nodes k and i, respectively. The normal vector to the plane can be obtained by the vector cross product of V
x
′ and V
i
k, which is(30)
V
z
′
=
V
x
′
×
V
i
k
.The direction cosine λ
z
′ for the local z
′ direction can be obtained by normalizing vector V
z
′ with respect to its length. The local y
′ axis is obtained by the vector cross product of the vectors V
z
′ and V
x
′. The cross product of these two vectors will give the vector V
y
′ normal to the x
′
o
′
z
′ plane, which is(31)
V
y
′
=
V
z
′
×
V
x
′
.The direction cosine λ
y
′ for the local y
′ direction is obtained by normalizing the vector V
y
′ with respect to its length, with reference to (28).So the 3 × 3 transformation matrixL for the transformation of coordinates from the local to the global axis can be expressed as(32)
L
=
λ
x
′
3
×
1
λ
y
′
3
×
1
λ
z
′
3
×
1
Tin which λ
x
′, λ
y
′, and λ
z
′ are the direction cosines for local x
′, y
′, and z
′ directions. The 12 × 12 transformation matrix R for the 4-node quadrilateral element stiffness from the local to the global coordinate systems can be represented by L, which is(33)
R
=
L
L
L
L
.To calculate the structure stiffness matrix, the element stiffness matrix must be transformed to the global coordinate system. The transformation of the GQ12 element stiffness matrix from the local to the global coordinate system is given by [3, 4](34)
K
G
Q
12
e
=
R
T
K
G
Q
12
′
e
R
=
R
T
∑
i
=
1
N
G
∑
j
=
1
N
G
W
i
W
j
B
T
D
B
J
t
R
.Similarly, from (20) and (23), the GQ12M element stiffness matrix under the global coordinate system can be described as(35)
K
G
Q
12
M
e
=
R
T
K
G
Q
12
M
′
e
R
=
R
T
K
q
q
′
e
R
-
K
ρ
q
′
e
R
T
K
ρ
ρ
′
e
-
1
K
ρ
q
′
e
Rin which(36)
R
T
K
q
q
′
e
R
=
R
T
∑
i
=
1
N
G
∑
j
=
1
N
G
W
i
W
j
B
T
D
B
J
t
R
,
(37)
K
ρ
q
′
e
R
=
∑
i
=
1
N
G
∑
j
=
1
N
G
W
i
W
j
B
^
T
D
B
J
t
R
.It can be seen from the above equations that the key task of computing the membrane stiffness matrix is to accurately evaluate the derivatives of the shape functions with respect to the local coordinates∂
N
i
/
∂
x
′ and ∂
N
i
/
∂
y
′, the determinant of the Jacobian matrix J, and transformation matrix R. The membrane element in the three-dimensional space may be not regular and coplanar. Thus, the calculation accuracy of the traditional way to establish the element local Cartesian coordinate system is heavily dependent on the element regularity.In order to improve the performance of the membrane elements GQ12 and GQ12M, a new method to establish the element local Cartesian coordinate system is proposed. And the corresponding techniques to calculate the derivatives of the shape functions with respect to the local coordinates, the transformation matrix from the local to the global coordinate systems, and the determinant of the Jacobian matrix are also provided in this paper [26, 27].


## 4. A New Method to Calculate the Stiffness Matrix of the Membrane Element in the Global Coordinate System

### 4.1. Establishing the Local Coordinate System

In an effort to avoid confusion and awkward phrasing, the signsx
1
,
x
2
,
x
3, x
1
′
,
x
2
′
,
x
3
′, and ξ
,
η are used to denote the symbols of the global Cartesian coordinate system, the local Cartesian coordinate system, and the intrinsic coordinates system, respectively.In order to improve the calculation accuracy, a new way [26] to establish the local Cartesian coordinate system is suggested in this paper. As shown in Figure 4, the origin of the local Cartesian coordinate system is set to the origin of the intrinsic coordinate system. Firstly, the vectors r
ξ and r
η from the origin need to be built along the tangent directions of the axes ξ and η, respectively, as illustrated in Figure 5. The direction of the vector r
ξ is considered as the local x
1
′ direction. Then the vectors r
ξ and r
η define the plane which is tangential to the element surface. The normal vector n to this plane can be obtained by the vector cross product of r
ξ and r
η, which is the local x
3
′ direction. At last, the local x
2
′ axis is obtained by the vector cross product of the vector in the local x
1
′ direction and vector in local x
3
′ direction. Now the local Cartesian coordinate system in which the axes x
1
′ and x
2
′ are tangential to the surface and x
3
′ is directed in the normal direction (shown in Figure 4) is established successfully.Figure 4
Local orthogonal set of axes over a surface element in the new way.Figure 5
The vectorsr
ξ and r
η in the local tangent plane to the surface.However, when the curvature of element surface is large, the numerical results are still not very accurate in this local Cartesian coordinate system. Therefore, in order to further enhance the precision of the calculation, the origin of the local Cartesian coordinate system can be set at the Gauss points. This idea came from the traction recovery method [28–31] in Boundary Element Method (BEM) which establishes the local Cartesian coordinate systems at each node of the element in order to improve the calculation accuracy. In finite element method (FEM), the numerical integration is performed on Gauss points, so the local Cartesian coordinate systems can be established at Gauss points. Take the 2 × 2 Gauss point integration as an example, as indicated in Figure 6; the local Cartesian coordinate systems with the origins at the four Gauss points can be established in a similar way as shown in Figure 4.Figure 6
Local Cartesian coordinate systems with the origins at the Gauss points.It is reasonable to establish the local Cartesian coordinate systems at Gauss points since the numerical integration is performed on them. In this case, the calculation accuracy can be improved, because the local Cartesian coordinate systems which are established on the tangent plane to the element surface at Gauss points can adapt to the curved element surface better. In order to better describe this problem, an example shown in Figure7 which is a curved element surface is introduced to illustrate the performance of this new local Cartesian coordinate systems established at Gauss points. Figure 8 describes the local Cartesian coordinate system established in the traditional method at the curved element surface. It can be observed that the element local plane has a large difference with the curved element surface. By comparison, the element local planes defined by the local Cartesian coordinate systems established at 2 × 2 Gauss points are more accordant with the curved element surface, as illustrated in Figure 9. In essence, Gaussian integration is the summation of the numerical results at Gauss points, so it is reasonable to believe that the calculation accuracy of the numerical integration can be improved by establishing the local Cartesian coordinate system at each Gauss point.Figure 7
A curved element surface.Figure 8
The local Cartesian coordinate system established in the traditional method and the element local plane at the curved element surface.Figure 9
The local Cartesian coordinate systems established in the new method and the element local planes at the curved element surface.As shown in Figure9, the local Cartesian coordinate systems are different since they are established in different tangent planes to the curved element surface. So the direction cosines λ
x
′, λ
y
′, and λ
z
′ of the local coordinate system established at different Gauss point with respect to the global coordinate system are different with each other. The detailed formulas and computation steps of λ
x
′, λ
y
′, and λ
z
′ will be given in the following sections. According to (31) and (32), it can be seen that the transformation matrix R for the transformation of element stiffness matrix from the local to the global coordinate systems may be different. So the element stiffness matrices of GQ12 and GQ12M under the global coordinate system expressed as (33), (35), and (36) can be modified to be(38)
K
G
Q
12
e
=
∑
i
=
1
N
G
∑
j
=
1
N
G
W
i
W
j
B
R
T
D
B
R
J
t
,
R
T
K
q
q
′
e
R
=
∑
i
=
1
N
G
∑
j
=
1
N
G
W
i
W
j
B
R
T
D
B
R
J
t
,
K
ρ
q
′
e
R
=
∑
i
=
1
N
G
∑
j
=
1
N
G
W
i
W
j
B
^
T
D
B
R
J
t
.


### 4.2. Derivatives of Shape Functions with respect to Local Coordinates

According to the derivation rule for compound function, the derivatives of the shape functions with respect to the local coordinates can be written as(39)
∂
N
i
∂
x
1
′
∂
N
i
∂
x
2
′
=
∂
ξ
∂
x
1
′
∂
η
∂
x
1
′
∂
ξ
∂
x
2
′
∂
η
∂
x
2
′
∂
N
i
∂
ξ
∂
N
i
∂
η
,
i
=
1,2
,
3,4in which ∂
ξ
/
∂
x
j
′ and ∂
η
/
∂
x
j
′
j
=
1,2 are the derivatives of the intrinsic coordinates with respect to the local coordinates, which can be evaluated using the method proposed by Lachat [26, 32–34]; that is,(40)
∂
ξ
∂
x
1
′
=
1
m
1
,
∂
ξ
∂
x
2
′
=
-
cos
⁡
θ
m
1
sin
⁡
θ
,
∂
η
∂
x
1
′
=
0
,
∂
η
∂
x
2
′
=
1
m
2
sin
⁡
θ
,where θ is the angle between x
1
′ axis and η axis shown in Figure 4, and(41)
m
1
=
∂
x
1
∂
ξ
2
+
∂
x
2
∂
ξ
2
+
∂
x
3
∂
ξ
2
,
m
2
=
∂
x
1
∂
η
2
+
∂
x
2
∂
η
2
+
∂
x
3
∂
η
2
,
(42)
cos
⁡
θ
=
1
m
1
m
2
∂
x
i
∂
ξ
∂
x
i
∂
η
,
i
=
1,2
,
3in which x
i
i
=
1,2
,
3 is the global coordinates of the element nodes, and ξ, η are the intrinsic coordinates.In (39), N
i is a substitution, which may be N
i
0, N
u
θ
i, N
v
θ
i, and N
^ in (6), (7), and (15), respectively. So separately substituting (6), (7), and (15) into (39) gives(43)
∂
N
i
0
∂
x
′
∂
N
i
0
∂
y
′
=
∂
N
i
0
∂
x
1
′
∂
N
i
0
∂
x
2
′
=
∂
ξ
∂
x
1
′
∂
η
∂
x
1
′
∂
ξ
∂
x
2
′
∂
η
∂
x
2
′
∂
N
i
0
∂
ξ
∂
N
i
0
∂
η
,
i
=
1,2
,
3,4
,
∂
N
u
θ
i
∂
x
′
∂
N
u
θ
i
∂
y
′
=
∂
N
u
θ
i
∂
x
1
′
∂
N
u
θ
i
∂
x
2
′
=
∂
ξ
∂
x
1
′
∂
η
∂
x
1
′
∂
ξ
∂
x
2
′
∂
η
∂
x
2
′
∂
N
u
θ
i
∂
ξ
∂
N
u
θ
i
∂
η
,
i
=
1,2
,
3,4
,
∂
N
v
θ
i
∂
x
′
∂
N
v
θ
i
∂
y
′
=
∂
N
v
θ
i
∂
x
1
′
∂
N
v
θ
i
∂
x
2
′
=
∂
ξ
∂
x
1
′
∂
η
∂
x
1
′
∂
ξ
∂
x
2
′
∂
η
∂
x
2
′
∂
N
v
θ
i
∂
ξ
∂
N
v
θ
i
∂
η
,
i
=
1,2
,
3,4
,
∂
N
^
∂
x
′
∂
N
^
∂
y
′
=
∂
N
^
∂
x
1
′
∂
N
^
∂
x
2
′
=
∂
ξ
∂
x
1
′
∂
η
∂
x
1
′
∂
ξ
∂
x
2
′
∂
η
∂
x
2
′
∂
N
^
∂
ξ
∂
N
^
∂
η
.The derivatives of the shape functions with respect to the local coordinates can be obtained by substituting (40) into the above equations. Then the strain matrices B and B
^ can be obtained by substituting the above equations into (10) and (18), respectively.


### 4.3. Determination of the JacobianJ

For the convenience of description, the local Cartesian coordinate system involved in this section is based on the way as shown in Figure4. For the local Cartesian coordinate systems with the origins at the Gauss points, the approach to determine J is similar.Referring to Figure10, the Jacobian is equal to the magnitude of the vector cross product of the vectors r
ξ and r
η, which are located in the local tangent plane to the surface as can also be seen in Figure 5. The Jacobian of the transformation from the global three-dimensional coordinate system to the two-dimensional intrinsic coordinate system of the surface patch is introduced in the form of [26](44)
J
′
=
n
∗
=
r
ξ
×
r
η
,where(45)
r
ξ
=
∂
x
1
ξ
,
η
∂
ξ
i
+
∂
x
2
ξ
,
η
∂
ξ
j
+
∂
x
3
ξ
,
η
∂
ξ
k
,
r
η
=
∂
x
1
ξ
,
η
∂
η
i
+
∂
x
2
ξ
,
η
∂
η
j
+
∂
x
3
ξ
,
η
∂
η
k
,
(46)
n
∗
=
r
ξ
×
r
η
=
∂
x
2
ξ
,
η
∂
ξ
·
∂
x
3
ξ
,
η
∂
η
-
∂
x
3
ξ
,
η
∂
ξ
·
∂
x
2
ξ
,
η
∂
η
i
+
∂
x
3
ξ
,
η
∂
ξ
·
∂
x
1
ξ
,
η
∂
η
-
∂
x
1
ξ
,
η
∂
ξ
·
∂
x
3
ξ
,
η
∂
η
j
+
∂
x
1
ξ
,
η
∂
ξ
·
∂
x
2
ξ
,
η
∂
η
-
∂
x
2
ξ
,
η
∂
ξ
·
∂
x
1
ξ
,
η
∂
η
k
,in which i, j, and k are the orthogonal unit basis vectors of the global coordinate axes and n
∗ is normal to the surface. So the unit normal vector may be obtained from (47)
n
=
n
∗
n
∗
=
n
1
i
+
n
2
j
+
n
3
k
,in which n
1, n
2, and n
3 are the magnitude of the components of the unit normal vector in the global coordinate axes.Figure 10
Normal vector to a surface element.The Jacobian of the transformation from the local orthogonal coordinate system to the intrinsic coordinate system is the same as that from the global orthogonal to the intrinsic coordinate systems. So the value ofJ in (38) can be computed using (44).


## 4.1. Establishing the Local Coordinate System

In an effort to avoid confusion and awkward phrasing, the signsx
1
,
x
2
,
x
3, x
1
′
,
x
2
′
,
x
3
′, and ξ
,
η are used to denote the symbols of the global Cartesian coordinate system, the local Cartesian coordinate system, and the intrinsic coordinates system, respectively.In order to improve the calculation accuracy, a new way [26] to establish the local Cartesian coordinate system is suggested in this paper. As shown in Figure 4, the origin of the local Cartesian coordinate system is set to the origin of the intrinsic coordinate system. Firstly, the vectors r
ξ and r
η from the origin need to be built along the tangent directions of the axes ξ and η, respectively, as illustrated in Figure 5. The direction of the vector r
ξ is considered as the local x
1
′ direction. Then the vectors r
ξ and r
η define the plane which is tangential to the element surface. The normal vector n to this plane can be obtained by the vector cross product of r
ξ and r
η, which is the local x
3
′ direction. At last, the local x
2
′ axis is obtained by the vector cross product of the vector in the local x
1
′ direction and vector in local x
3
′ direction. Now the local Cartesian coordinate system in which the axes x
1
′ and x
2
′ are tangential to the surface and x
3
′ is directed in the normal direction (shown in Figure 4) is established successfully.Figure 4
Local orthogonal set of axes over a surface element in the new way.Figure 5
The vectorsr
ξ and r
η in the local tangent plane to the surface.However, when the curvature of element surface is large, the numerical results are still not very accurate in this local Cartesian coordinate system. Therefore, in order to further enhance the precision of the calculation, the origin of the local Cartesian coordinate system can be set at the Gauss points. This idea came from the traction recovery method [28–31] in Boundary Element Method (BEM) which establishes the local Cartesian coordinate systems at each node of the element in order to improve the calculation accuracy. In finite element method (FEM), the numerical integration is performed on Gauss points, so the local Cartesian coordinate systems can be established at Gauss points. Take the 2 × 2 Gauss point integration as an example, as indicated in Figure 6; the local Cartesian coordinate systems with the origins at the four Gauss points can be established in a similar way as shown in Figure 4.Figure 6
Local Cartesian coordinate systems with the origins at the Gauss points.It is reasonable to establish the local Cartesian coordinate systems at Gauss points since the numerical integration is performed on them. In this case, the calculation accuracy can be improved, because the local Cartesian coordinate systems which are established on the tangent plane to the element surface at Gauss points can adapt to the curved element surface better. In order to better describe this problem, an example shown in Figure7 which is a curved element surface is introduced to illustrate the performance of this new local Cartesian coordinate systems established at Gauss points. Figure 8 describes the local Cartesian coordinate system established in the traditional method at the curved element surface. It can be observed that the element local plane has a large difference with the curved element surface. By comparison, the element local planes defined by the local Cartesian coordinate systems established at 2 × 2 Gauss points are more accordant with the curved element surface, as illustrated in Figure 9. In essence, Gaussian integration is the summation of the numerical results at Gauss points, so it is reasonable to believe that the calculation accuracy of the numerical integration can be improved by establishing the local Cartesian coordinate system at each Gauss point.Figure 7
A curved element surface.Figure 8
The local Cartesian coordinate system established in the traditional method and the element local plane at the curved element surface.Figure 9
The local Cartesian coordinate systems established in the new method and the element local planes at the curved element surface.As shown in Figure9, the local Cartesian coordinate systems are different since they are established in different tangent planes to the curved element surface. So the direction cosines λ
x
′, λ
y
′, and λ
z
′ of the local coordinate system established at different Gauss point with respect to the global coordinate system are different with each other. The detailed formulas and computation steps of λ
x
′, λ
y
′, and λ
z
′ will be given in the following sections. According to (31) and (32), it can be seen that the transformation matrix R for the transformation of element stiffness matrix from the local to the global coordinate systems may be different. So the element stiffness matrices of GQ12 and GQ12M under the global coordinate system expressed as (33), (35), and (36) can be modified to be(38)
K
G
Q
12
e
=
∑
i
=
1
N
G
∑
j
=
1
N
G
W
i
W
j
B
R
T
D
B
R
J
t
,
R
T
K
q
q
′
e
R
=
∑
i
=
1
N
G
∑
j
=
1
N
G
W
i
W
j
B
R
T
D
B
R
J
t
,
K
ρ
q
′
e
R
=
∑
i
=
1
N
G
∑
j
=
1
N
G
W
i
W
j
B
^
T
D
B
R
J
t
.


## 4.2. Derivatives of Shape Functions with respect to Local Coordinates

According to the derivation rule for compound function, the derivatives of the shape functions with respect to the local coordinates can be written as(39)
∂
N
i
∂
x
1
′
∂
N
i
∂
x
2
′
=
∂
ξ
∂
x
1
′
∂
η
∂
x
1
′
∂
ξ
∂
x
2
′
∂
η
∂
x
2
′
∂
N
i
∂
ξ
∂
N
i
∂
η
,
i
=
1,2
,
3,4in which ∂
ξ
/
∂
x
j
′ and ∂
η
/
∂
x
j
′
j
=
1,2 are the derivatives of the intrinsic coordinates with respect to the local coordinates, which can be evaluated using the method proposed by Lachat [26, 32–34]; that is,(40)
∂
ξ
∂
x
1
′
=
1
m
1
,
∂
ξ
∂
x
2
′
=
-
cos
⁡
θ
m
1
sin
⁡
θ
,
∂
η
∂
x
1
′
=
0
,
∂
η
∂
x
2
′
=
1
m
2
sin
⁡
θ
,where θ is the angle between x
1
′ axis and η axis shown in Figure 4, and(41)
m
1
=
∂
x
1
∂
ξ
2
+
∂
x
2
∂
ξ
2
+
∂
x
3
∂
ξ
2
,
m
2
=
∂
x
1
∂
η
2
+
∂
x
2
∂
η
2
+
∂
x
3
∂
η
2
,
(42)
cos
⁡
θ
=
1
m
1
m
2
∂
x
i
∂
ξ
∂
x
i
∂
η
,
i
=
1,2
,
3in which x
i
i
=
1,2
,
3 is the global coordinates of the element nodes, and ξ, η are the intrinsic coordinates.In (39), N
i is a substitution, which may be N
i
0, N
u
θ
i, N
v
θ
i, and N
^ in (6), (7), and (15), respectively. So separately substituting (6), (7), and (15) into (39) gives(43)
∂
N
i
0
∂
x
′
∂
N
i
0
∂
y
′
=
∂
N
i
0
∂
x
1
′
∂
N
i
0
∂
x
2
′
=
∂
ξ
∂
x
1
′
∂
η
∂
x
1
′
∂
ξ
∂
x
2
′
∂
η
∂
x
2
′
∂
N
i
0
∂
ξ
∂
N
i
0
∂
η
,
i
=
1,2
,
3,4
,
∂
N
u
θ
i
∂
x
′
∂
N
u
θ
i
∂
y
′
=
∂
N
u
θ
i
∂
x
1
′
∂
N
u
θ
i
∂
x
2
′
=
∂
ξ
∂
x
1
′
∂
η
∂
x
1
′
∂
ξ
∂
x
2
′
∂
η
∂
x
2
′
∂
N
u
θ
i
∂
ξ
∂
N
u
θ
i
∂
η
,
i
=
1,2
,
3,4
,
∂
N
v
θ
i
∂
x
′
∂
N
v
θ
i
∂
y
′
=
∂
N
v
θ
i
∂
x
1
′
∂
N
v
θ
i
∂
x
2
′
=
∂
ξ
∂
x
1
′
∂
η
∂
x
1
′
∂
ξ
∂
x
2
′
∂
η
∂
x
2
′
∂
N
v
θ
i
∂
ξ
∂
N
v
θ
i
∂
η
,
i
=
1,2
,
3,4
,
∂
N
^
∂
x
′
∂
N
^
∂
y
′
=
∂
N
^
∂
x
1
′
∂
N
^
∂
x
2
′
=
∂
ξ
∂
x
1
′
∂
η
∂
x
1
′
∂
ξ
∂
x
2
′
∂
η
∂
x
2
′
∂
N
^
∂
ξ
∂
N
^
∂
η
.The derivatives of the shape functions with respect to the local coordinates can be obtained by substituting (40) into the above equations. Then the strain matrices B and B
^ can be obtained by substituting the above equations into (10) and (18), respectively.


## 4.3. Determination of the JacobianJ

For the convenience of description, the local Cartesian coordinate system involved in this section is based on the way as shown in Figure4. For the local Cartesian coordinate systems with the origins at the Gauss points, the approach to determine J is similar.Referring to Figure10, the Jacobian is equal to the magnitude of the vector cross product of the vectors r
ξ and r
η, which are located in the local tangent plane to the surface as can also be seen in Figure 5. The Jacobian of the transformation from the global three-dimensional coordinate system to the two-dimensional intrinsic coordinate system of the surface patch is introduced in the form of [26](44)
J
′
=
n
∗
=
r
ξ
×
r
η
,where(45)
r
ξ
=
∂
x
1
ξ
,
η
∂
ξ
i
+
∂
x
2
ξ
,
η
∂
ξ
j
+
∂
x
3
ξ
,
η
∂
ξ
k
,
r
η
=
∂
x
1
ξ
,
η
∂
η
i
+
∂
x
2
ξ
,
η
∂
η
j
+
∂
x
3
ξ
,
η
∂
η
k
,
(46)
n
∗
=
r
ξ
×
r
η
=
∂
x
2
ξ
,
η
∂
ξ
·
∂
x
3
ξ
,
η
∂
η
-
∂
x
3
ξ
,
η
∂
ξ
·
∂
x
2
ξ
,
η
∂
η
i
+
∂
x
3
ξ
,
η
∂
ξ
·
∂
x
1
ξ
,
η
∂
η
-
∂
x
1
ξ
,
η
∂
ξ
·
∂
x
3
ξ
,
η
∂
η
j
+
∂
x
1
ξ
,
η
∂
ξ
·
∂
x
2
ξ
,
η
∂
η
-
∂
x
2
ξ
,
η
∂
ξ
·
∂
x
1
ξ
,
η
∂
η
k
,in which i, j, and k are the orthogonal unit basis vectors of the global coordinate axes and n
∗ is normal to the surface. So the unit normal vector may be obtained from (47)
n
=
n
∗
n
∗
=
n
1
i
+
n
2
j
+
n
3
k
,in which n
1, n
2, and n
3 are the magnitude of the components of the unit normal vector in the global coordinate axes.Figure 10
Normal vector to a surface element.The Jacobian of the transformation from the local orthogonal coordinate system to the intrinsic coordinate system is the same as that from the global orthogonal to the intrinsic coordinate systems. So the value ofJ in (38) can be computed using (44).


## 5. The Matrix Algorithm for the Membrane Element Stiffness Transformation from the Local to the Global Coordinate Systems

As described in Figure4, the direction of local x
1
′ axis is the tangent direction of the axis ξ at the origin, which is the direction of the vector r
ξ. So the direction cosine λ
x
1
′ for the local x
1
′ direction is obtained by normalizing the vector r
ξ with respect to its length. According to (41) and (45), the direction cosine λ
x
1
′ can be written as [26](48)
λ
x
′
=
λ
x
1
′
=
1
m
1
∂
x
1
∂
ξ
∂
x
2
∂
ξ
∂
x
3
∂
ξ
,in which x
i
i
=
1,2
,
3 is the global coordinates of the element nodes and m
1 is the mold of the vector r
ξ expressed in (41). The local x
3
′ direction is the direction of the normal vector n to the element surface which can be obtained from (47). So the direction cosine λ
x
3
′ can be written as(49)
λ
z
′
=
λ
x
3
′
=
n
1
n
2
n
3in which n
1, n
2, and n
3 are given by (47). Because the local x
2
′ axis can be obtained by the vector cross product of the vector in the local x
1
′ direction and vector in local x
3
′ direction, the direction cosine λ
x
2
′ can be written as(50)
λ
y
′
=
λ
x
2
′
=
λ
z
′
×
λ
x
′
.So the 3 × 3 transformation matrixL for the transformation of coordinates from the local to the global axis can be obtained from (32), which can be expressed in terms of L
i
j
i
,
j
=
1,2
,
3 as [26](51)
L
1
j
=
1
m
1
∂
x
j
∂
ξ
,
j
=
1,2
,
3
,
L
21
=
n
2
L
13
-
n
3
L
12
,
L
22
=
n
3
L
11
-
n
1
L
13
,
L
23
=
n
1
L
12
-
n
2
L
11
,
L
3
j
=
n
j
,
j
=
1,2
,
3
.The 12 × 12 transformation matrixR for the 4-node quadrilateral element stiffness from the local to the global coordinate systems can be obtained from (33). It is worth mentioning that the matrix R is just used for the membrane element stiffness transformation from the local to the global coordinate systems. For the 4-node quadrilateral shell element which has six degrees of freedom in each node, the 24 × 24 transformation matrix T for the transformation of shell element stiffness from the local to the global coordinate systems is written as [4](52)
T
=
R
R
.


## 6. Numerical Examples

In order to examine the performance of the new method for GQ12 and GQ12M described in this paper, four numerical examples are proposed in this section. And the computational results are compared with those calculated by the traditional method based on the elements GQ12 and GQ12M.


### 6.1. Example 1: Cook’s Problem

The problem defined in Figure11 was proposed by Cook [36] as a test case of plane stress elements, which is a standard example for accuracy testing of plane stress problem. Young’s modulus of the element is E
=
1Pa, and Poisson’s ratio is v
=
1
/
3. The structure is meshed with 2 × 2, 4 × 4, and 8 × 8 elements. Here the new methods based on the elements GQ12 and GQ12M proposed in this paper are utilized to solve this problem. The best known answers taken from [35] are used for comparison, since there is no analytical solution available for this problem. The results of the vertical deflection at C, the maximum stress at A, and the minimum stress at B on different meshes of the structure are given in Tables 2 and 3. The results demonstrate that the new method for the elements GQ12 and GQ12M have the desirable numerical accuracy, both for the displacement and for the stress.Table 2
The vertical deflection results of Cook’s problem.
Element
Vertical deflection at C (m)
(2
×
2) mesh
(4
×
4) mesh
(8
×
8) mesh
Deflection
Deflection
Deflection
Q4 [4]
11.80
18.29
22.08
Allman [15]
20.27
22.78
23.56
GQ12 [25]
20.89
23.06
23.67
GQ12 (present method)
21.27
23.07
23.67
GQ12M [25]
21.69
23.30
23.74
GQ12M (present method)
22.47
23.43
23.78
Comparison solution [35]
23.96Table 3
The stress results of Cook’s problem.
Element
Maximum stress at A (Pa)
Minimum stress at B (Pa)
(2
×
2) mesh
(4
×
4) mesh
(8
×
8) mesh
σ
A
σ
B
σ
A
σ
B
σ
A
σ
B
Q4 [4]
0.1217
−0.0957
0.1873
−0.1524
0.2242
−0.1869
Allman [15]
0.1825
−0.1716
0.2261
−0.1921
0.2340
−0.2004
GQ12 [25]
0.1802
−0.1784
0.2209
−0.1950
0.2315
−0.2007
GQ12 (present method)
0.1841
−0.1810
0.2214
−0.1948
0.2315
−0.2008
GQ12M [25]
0.1861
−0.2088
0.2232
−0.2017
0.2305
−0.2017
GQ12M (present method)
0.2070
−0.2142
0.2333
−0.2037
0.2361
−0.2027
Comparison solution [35]
σ
A
=
0.2362
σ
B
=
-
0.2023Figure 11
Cook’s problem.


### 6.2. Example 2: MacNeal Beam

In order to examine the membrane locking phenomena, a typical test example [37, 38] defined in Figure 12 was proposed by MacNeal. The thickness of the membrane elements is t
=
0.1m, Young’s modulus is E
=
1
0
7 Pa, and Poisson’s ratio is v
=
0.3. The beam end displacements obtained with several elements under the unit shear load at the free end are compared with the theoretical values taken from [38]. The results under different meshes in Table 4 demonstrate that GQ12 and GQ12M with the new method described in this paper can delete membrane locking and have the same accuracy as the traditional method.Table 4
The results of the MacNeal beam under different meshes.
Element
The beam end deflection (m)
Mesh (a)
Mesh (b)
Mesh (c)
Q4 [4]
−0.0904
−0.0710
−0.0800
GQ12 [25]
−0.0904
−0.0767
−0.0860
GQ12 (present method)
−0.0977
−0.0871
−0.0944
GQ12M [25]
−0.0977
−0.0886
−0.0923
GQ12M (present method)
−0.1073
−0.0980
−0.1050
Theoretical value [38]
−0.1081Figure 12
MacNeal beam.
(a)
Rectangular mesh
(b)
Trapezoid mesh
(c)
Rhomboid mesh


### 6.3. Example 3: Uniform Stretching (Patch Test) and Bending of a Plate

This test example [25] is a rectangular plate under two loading cases, which are the uniform stretching under Load 1 and the pure bending under Load 2, respectively. The uniform stretching case is the patch test problem under constant strain. Because of the symmetry of the model and loads, only a quarter of the plate with the irregular mesh shown in Figure 13 is considered. Young’s modulus of the element is E
=
1Pa, and Poisson’s ratio is v
=
1
/
3. Table 5 shows the results of displacement at corner A which are obtained with different elements under these two loading cases. It demonstrates that the elements GQ12 and GQ12M with present method in this paper pass the patch test for a general quadrilateral mesh and achieve more accurate results.Table 5
The displacement at corner A under two loading cases.
Element
Load 1
Load 2
U
A (m)
V
A (m)
Q4 [4]
6.00
−17.00
GQ12 [25]
6.00
−17.72
GQ12 (present method)
6.00
−17.75
GQ12M [25]
6.00
−17.77
GQ12M (present method)
6.00
−17.82
Exact [39]
6.00
−18.00Figure 13
A rectangular plate under two loading cases.
(a)
Uniform stretching problem (patch test problem)
(b)
Pure bending problem


### 6.4. Example 4: Linear Analysis of a Hemispherical Shell

The aim of this test example [40] is to examine the performance of GQ12 and GQ12M elements with the present method in this paper as a constituent part of the flat shell element to compute a curved shell structure. As known, the flat shell elements are formed by superimposing the stiffness of membrane and plate bending elements. In this section, the plate bending element is based on the Mindlin theory which includes transverse shear deformations, and the membrane element may be selected from GQ12 and GQ12M with present method in this paper, respectively.Figure14 shows a hemispherical shell structure with an 18-degree hole at the top [40]. The radius of the shell is 10 m, thickness is 0.04 m, Young’s modulus is 6.825
×
1
0
7 Pa, and Poisson’s ratio is 0.3. The top and bottom circumferential edges of the hemisphere are free and the shell is subjected to two radial unit point loads. Only a quarter of the hemispherical shell with the meshes shown in Figure 15 is separated out for research due to the geometric symmetry. The radial displacement at point A from different meshes is compared with the theoretical solution u
A
=
0.094 [38] in Table 6. It can be seen that the shell elements based on the membrane elements with present method in this paper can achieve more accurate results than those based on a traditional manner.Table 6
The radial displacement at point A under different meshes.
Flat shell element
⁢The radial displacement at A under different meshes (m)
4
×
4
8
×
8
12
×
12
16
×
16
Mindlin + GQ12
0.02234
0.05708
0.07222
0.08019
Mindlin + GQ12 (present method)
0.03321
0.08975
0.09275
0.09329
Mindlin + GQ12M
0.02365
0.06000
0.07470
0.08197
Mindlin + GQ12M (present method)
0.03642
0.09027
0.09287
0.09335
Comparison solution [38]
0.094Figure 14
A hemispherical shell structure with an 18-degree hole on the top.Figure 15
A quarter of the hemispherical shell structure with the meshes.


## 6.1. Example 1: Cook’s Problem

The problem defined in Figure11 was proposed by Cook [36] as a test case of plane stress elements, which is a standard example for accuracy testing of plane stress problem. Young’s modulus of the element is E
=
1Pa, and Poisson’s ratio is v
=
1
/
3. The structure is meshed with 2 × 2, 4 × 4, and 8 × 8 elements. Here the new methods based on the elements GQ12 and GQ12M proposed in this paper are utilized to solve this problem. The best known answers taken from [35] are used for comparison, since there is no analytical solution available for this problem. The results of the vertical deflection at C, the maximum stress at A, and the minimum stress at B on different meshes of the structure are given in Tables 2 and 3. The results demonstrate that the new method for the elements GQ12 and GQ12M have the desirable numerical accuracy, both for the displacement and for the stress.Table 2
The vertical deflection results of Cook’s problem.
Element
Vertical deflection at C (m)
(2
×
2) mesh
(4
×
4) mesh
(8
×
8) mesh
Deflection
Deflection
Deflection
Q4 [4]
11.80
18.29
22.08
Allman [15]
20.27
22.78
23.56
GQ12 [25]
20.89
23.06
23.67
GQ12 (present method)
21.27
23.07
23.67
GQ12M [25]
21.69
23.30
23.74
GQ12M (present method)
22.47
23.43
23.78
Comparison solution [35]
23.96Table 3
The stress results of Cook’s problem.
Element
Maximum stress at A (Pa)
Minimum stress at B (Pa)
(2
×
2) mesh
(4
×
4) mesh
(8
×
8) mesh
σ
A
σ
B
σ
A
σ
B
σ
A
σ
B
Q4 [4]
0.1217
−0.0957
0.1873
−0.1524
0.2242
−0.1869
Allman [15]
0.1825
−0.1716
0.2261
−0.1921
0.2340
−0.2004
GQ12 [25]
0.1802
−0.1784
0.2209
−0.1950
0.2315
−0.2007
GQ12 (present method)
0.1841
−0.1810
0.2214
−0.1948
0.2315
−0.2008
GQ12M [25]
0.1861
−0.2088
0.2232
−0.2017
0.2305
−0.2017
GQ12M (present method)
0.2070
−0.2142
0.2333
−0.2037
0.2361
−0.2027
Comparison solution [35]
σ
A
=
0.2362
σ
B
=
-
0.2023Figure 11
Cook’s problem.


## 6.2. Example 2: MacNeal Beam

In order to examine the membrane locking phenomena, a typical test example [37, 38] defined in Figure 12 was proposed by MacNeal. The thickness of the membrane elements is t
=
0.1m, Young’s modulus is E
=
1
0
7 Pa, and Poisson’s ratio is v
=
0.3. The beam end displacements obtained with several elements under the unit shear load at the free end are compared with the theoretical values taken from [38]. The results under different meshes in Table 4 demonstrate that GQ12 and GQ12M with the new method described in this paper can delete membrane locking and have the same accuracy as the traditional method.Table 4
The results of the MacNeal beam under different meshes.
Element
The beam end deflection (m)
Mesh (a)
Mesh (b)
Mesh (c)
Q4 [4]
−0.0904
−0.0710
−0.0800
GQ12 [25]
−0.0904
−0.0767
−0.0860
GQ12 (present method)
−0.0977
−0.0871
−0.0944
GQ12M [25]
−0.0977
−0.0886
−0.0923
GQ12M (present method)
−0.1073
−0.0980
−0.1050
Theoretical value [38]
−0.1081Figure 12
MacNeal beam.
(a)
Rectangular mesh
(b)
Trapezoid mesh
(c)
Rhomboid mesh


## 6.3. Example 3: Uniform Stretching (Patch Test) and Bending of a Plate

This test example [25] is a rectangular plate under two loading cases, which are the uniform stretching under Load 1 and the pure bending under Load 2, respectively. The uniform stretching case is the patch test problem under constant strain. Because of the symmetry of the model and loads, only a quarter of the plate with the irregular mesh shown in Figure 13 is considered. Young’s modulus of the element is E
=
1Pa, and Poisson’s ratio is v
=
1
/
3. Table 5 shows the results of displacement at corner A which are obtained with different elements under these two loading cases. It demonstrates that the elements GQ12 and GQ12M with present method in this paper pass the patch test for a general quadrilateral mesh and achieve more accurate results.Table 5
The displacement at corner A under two loading cases.
Element
Load 1
Load 2
U
A (m)
V
A (m)
Q4 [4]
6.00
−17.00
GQ12 [25]
6.00
−17.72
GQ12 (present method)
6.00
−17.75
GQ12M [25]
6.00
−17.77
GQ12M (present method)
6.00
−17.82
Exact [39]
6.00
−18.00Figure 13
A rectangular plate under two loading cases.
(a)
Uniform stretching problem (patch test problem)
(b)
Pure bending problem


## 6.4. Example 4: Linear Analysis of a Hemispherical Shell

The aim of this test example [40] is to examine the performance of GQ12 and GQ12M elements with the present method in this paper as a constituent part of the flat shell element to compute a curved shell structure. As known, the flat shell elements are formed by superimposing the stiffness of membrane and plate bending elements. In this section, the plate bending element is based on the Mindlin theory which includes transverse shear deformations, and the membrane element may be selected from GQ12 and GQ12M with present method in this paper, respectively.Figure14 shows a hemispherical shell structure with an 18-degree hole at the top [40]. The radius of the shell is 10 m, thickness is 0.04 m, Young’s modulus is 6.825
×
1
0
7 Pa, and Poisson’s ratio is 0.3. The top and bottom circumferential edges of the hemisphere are free and the shell is subjected to two radial unit point loads. Only a quarter of the hemispherical shell with the meshes shown in Figure 15 is separated out for research due to the geometric symmetry. The radial displacement at point A from different meshes is compared with the theoretical solution u
A
=
0.094 [38] in Table 6. It can be seen that the shell elements based on the membrane elements with present method in this paper can achieve more accurate results than those based on a traditional manner.Table 6
The radial displacement at point A under different meshes.
Flat shell element
⁢The radial displacement at A under different meshes (m)
4
×
4
8
×
8
12
×
12
16
×
16
Mindlin + GQ12
0.02234
0.05708
0.07222
0.08019
Mindlin + GQ12 (present method)
0.03321
0.08975
0.09275
0.09329
Mindlin + GQ12M
0.02365
0.06000
0.07470
0.08197
Mindlin + GQ12M (present method)
0.03642
0.09027
0.09287
0.09335
Comparison solution [38]
0.094Figure 14
A hemispherical shell structure with an 18-degree hole on the top.Figure 15
A quarter of the hemispherical shell structure with the meshes.


## 7. Concluding Discussion

In this paper a new method for establishing the local Cartesian coordinate system and calculating the derivatives of the shape function with respect to local coordinates is applied to the membrane elements GQ12 and GQ12M. When the membrane elements are introduced to be the component of the flat shell element, this new method can establish the local Cartesian coordinate system on each Gauss point over the tangent plane to the element surface, and the precision is improved because the membrane elements are able to adapt the curved surface better. The numerical results of the test problems show that the elements GQ12 and GQ12M with the new method can obtain comparatively high accuracies. This new method can bring some new ideas and approaches to improve the computing accuracies of other membrane elements and flat shell elements, so it has a good application prospect in the future.



---

*Source: 1045438-2016-02-17.xml*