# Study on the Explicit Formula of the Triangular Flat Shell Element Based on the Analytical Trial Functions for Anisotropy Material

Authors: Xiang-Rong Fu; Li-Na Ge; Ge Tian; Ming-Wu Yuan
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

This paper presents a novel way to formulate the triangular flat shell element. The basic analytical solutions of membrane and bending plate problem for anisotropy material are studied separately. Combining with the conforming displacement along the sides and hybrid element strategy, the triangular flat shell elements based on the analytical trial functions (ATF) for anisotropy material are formulated. By using the explicit integral formulae of the triangular element, the matrices used in proposed shell element are calculated efficiently. The benchmark examples showed the high accuracy and high efficiency.


---

## Body

## 1. Introduction

The triangular flat shell element is included widely in finite element software packages. In order to improve the performances of the triangular flat shell element, most of studies focused on two challenges [1–15].One challenge is the drilling degree of freedom in the membrane element. Introducing the in-plane drilling degree of freedom, the singularity of global stiffness can be avoided in the analysis of shell structures when the neighboring elements around a common node are close to coplanar. Olson introduced the drilling degree of freedom to describe the displacement of the flat shell element [4]. Mohr constructed the hybrid membrane element with the drilling degree of freedom [5]. Allman presented a rational interpolation function to construct the displacement fields associated with the drilling degree of freedom in the triangular membrane elements [6]. Many researchers made different possibilities to improve accuracy of the membrane element by introducing the drilling degree of freedom to the corner nodes of triangular or quadrilateral elements [7–9].Another challenge is the shear locking of some bending plate elements based on the Mindlin/Reissner thick plate theory [10–15]. The pioneering work of Adini, Clough, and Melosh studied the famous thin plate element ACM, which showed the difficulties in constructing the C1 conforming displacement along the boundaries of the thin plate elements [10]. The thick plate element based on the Mindlin/Reissner theory can avoid the strict C1 conforming condition [11–15], but the shear locking problem in these elements became the main challenge, which raised the interest of many well-known researchers, including Zienkiewicz et al. [11], Hughes et al. [12], Hinton and Huang [13], Belytschko et al. [14], Bathe and Dvorkin [15].Long gave a general strategy to deal with the nonconforming displacement along the boundary between the neighboring elements, which was named as the generalized conforming element [16]. Taking advantage of the conforming boundary displacement, the nonconforming displacement model can be used freely to describe the inner field of the element. The generalized conforming equations which are utilized to determine the parameters of the proposed displacement model were studied systematically in the past decades [17–20]. Cen and Long gave a rational way to define the displacement field in the thick plate elements [19]. Fu et al. proved that the hybrid element using the balanced stress fields and the generalized conforming element introducing the associated generalized conforming equations can derive the same transform matrix between the inner parameters and the nodal displacements in the element [20].Fu et al. and Cen studied the hybrid element based on the variational principle containing the stress function [21]. A series of works showed that elements based on the analytical trial functions (ATF) would be insensitive to the distortion of the element shape [22–25]. The hybrid element strategy was improved through turning to the stress trial functions which are derived from the stress functions (e.g., the Airy stress functions for plane problem). The idea of ATF can be traced back to the first element proposed by Clough which is named as CST (constant stress triangular) element [26]. Some works showed that the element CST can be derived from three rigid body displacements and three displacements of constant stress [1]. On the other hand, the first hybrid element also employed five analytical trial functions of stress [27].This paper studies the triangular flat shell element based on the analytical trial functions of anisotropy. There are mainly four parts as follows.(a)
The first part extracts the general analytical solutions of plane and bending plate problem in anisotropy elasticity (in Section2).(b)
The second part gives the hybrid element strategy based on the analytical trial functions of stress derived in the first part (in Section3).(c)
The third part studies the explicit formulae of the hybrid element based on explicit integral formulae of the triangular element (in Section4).(d)
The fourth part shows some numerical examples as benchmark to study the accuracy and the efficiency of the proposed element model (in Section5).The proposed shell element of anisotropy has high accuracy and computational efficiency.


## 2. The Analytical Trial Functions (ATF)

### 2.1. The Analytical Trial Functions of the Plane Problem of Anisotropy

The general analytical solutions play very important roles in many methods of computational mechanics [21–27]. Article [28] provided an efficient approach named characteristic differential equation method (CDEM) to derive the general analytical solutions from the governing differential equations directly.Following the strategy of CDEM to study the plane problem of anisotropy, stress-displacement relationship can be described as(1){σxσyτxy}=[C11C12C16C12C22C26C16C26C66]{∂u∂x∂v∂y∂v∂x+∂u∂y}.
And the equilibrium equations in terms of displacement can be written as
(2)C11∂2u∂x2+C66∂2u∂y2+2C16∂2u∂x∂y+C16∂2v∂x2+C26∂2v∂y2+(C66+C12)∂2v∂x∂y=0,(3)C16∂2u∂x2+C26∂2u∂y2+(C66+C12)∂2u∂x∂y+C66∂2v∂x2+C22∂2v∂y2+2C26∂2v∂x∂y=0.
Utilizing the determinant expansion equation of the differential operator matrix in (2), the characteristic differential equation can be given as
(4)(C11C66-C612)∂4φ∂x4+(C22C66-C622)∂4φ∂y4+(2C11C62-2C61C21)∂4φ∂x3∂y+(2C61C22-2C62C21)∂4φ∂x∂y3+(C11C22-2C66C21-C122+2C61C62)×∂4φ∂x2∂y2=0.Then the general analytical solutions of the displacement can be derived from the characteristic solutionφ of (4):
(5){uv}={C66∂2∂x2+C22∂2∂y2+2C26∂2∂x∂y-(C16∂2∂x2+C26∂2∂y2+(C66+C12)∂2∂x∂y)}φ.
Substituting (5) into (1), the general stress solutions can be derived from the characteristic solutions φ. Table 1 gives seven fundamental analytical stress solutions of the anisotropy plane problem, which can be selected as trial functions of the following numerical method.Table 1
The fundamental stress solutions of the anisotropic plane problem.
Item
1
2
3
4
5
6
7
σ
x
1
0
0
(
C
11
C
66
-
C
61
2
)
(
(
C
22
C
66
-
C
62
2)x+(C21C62-C61C22)y)
(
C
11
C
66
-
C
61
2)((C21C66-C62C61)x +2(C62C11-C21C61)y)
(
C
11
C
66
-
C
61
2)((C61C22-C62C21)x +(C21C66-C22C11+C212-C61C62)y)
(
C
11
C
66
-
C
61
2
)
2
y
σ
y
0
1
0
(
C
21
C
66
-
C
61
C
62)(C66C22-C622)x
(
2
C
66
2
C
21
2
+
C
66
C
21
3
-
4
C
21
C
66
C
61
C
62
-
C
11
C
66
2
C
22
-
C
22
C
21
C
11
C
66
+
C
61
2
C
66
C
22
+
C
11
C
66
C
62
2
+
C
61
2
C
62
2
+
C
62
C
61
C
11
C
22
-
C
62
C
61
C
21
2
)
x
+
(
C
21
C
62
C
11
C
66
-
C
21
C
62
C
61
2
-
C
61
C
22
C
11
C
66
+
C
22
C
61
3
)
y
2
(
C
61
C
62
2
C
21
-
C
66
C
62
C
21
2
+
C
21
C
61
C
66
C
22
-
C
22
C
62
C
61
2)x  + (C22C662C11-C22C66C612-C622C11C66+C622C612)y
(
2
C
61
C
21
2
C
66
-
C
11
C
66
C
62
C
21
-
3
C
21
C
62
C
61
2
-
C
11
C
61
C
66
C
22
+
C
22
C
61
3
+
2
C
61
C
62
2
C
11
)
x +(C62C613-C61C62C11C66+C21C662C11-C21C66C612)y
τ
x
y
0
0
1
-
(
C
11
C
66
-
C
61
2)(C66C22-C622)y
(
C
11
C
66
-
C
61
2
) ((C61C22-C62C21)x +(C61C62-C21C66)y)
(
C
61
2
-
C
11
C
66)((C22C66-C622)x+(C61C22-C21C62)y)
(
C
11
C
66
-
C
61
2
) (C61C62-C21C66)x


### 2.2. The Analytical Trial Functions of the Bending Plate Problem of Anisotropy

The characteristic differential equation of the displacement solutions can be written as [29]
(6)-t312[C11∂4w∂x4+(2C12+4C66)∂4w∂x2∂y2,mm,,+4C61∂4w∂x3∂y+4C62∂4w∂x∂y3+C22∂4w∂y4]+q(x,y)=0.Deriving from the general analytical displacement solutions of (6), Table 2 presents seven general analytical stress solutions of the anisotropy plate problem, which will be selected as trial functions of the following hybrid element in Section 3.Table 2
The fundamental stress solutions of anisotropic material used in bending plate problem.
Item
1
2
3
4
5
6
7
M
x
1
0
0
C
11
x
C
12
y
C
11
y
+
2
C
61
x
C
12
x
+
2
C
61
y
M
y
0
1
0
C
12
x
C
22
y
C
12
y
+
2
C
62
x
C
22
x
+
2
C
62
y
M
x
y
0
0
1
C
61
x
C
62
y
C
61
y
+
2
C
66
x
C
62
x
+
2
C
66
y
Q
x
0
0
0
C
11
C
26
3
C
16
C
12
+
2
C
66
Q
y
0
0
0
C
16
C
22
C
12
+
2
C
66
3
C
26


## 2.1. The Analytical Trial Functions of the Plane Problem of Anisotropy

The general analytical solutions play very important roles in many methods of computational mechanics [21–27]. Article [28] provided an efficient approach named characteristic differential equation method (CDEM) to derive the general analytical solutions from the governing differential equations directly.Following the strategy of CDEM to study the plane problem of anisotropy, stress-displacement relationship can be described as(1){σxσyτxy}=[C11C12C16C12C22C26C16C26C66]{∂u∂x∂v∂y∂v∂x+∂u∂y}.
And the equilibrium equations in terms of displacement can be written as
(2)C11∂2u∂x2+C66∂2u∂y2+2C16∂2u∂x∂y+C16∂2v∂x2+C26∂2v∂y2+(C66+C12)∂2v∂x∂y=0,(3)C16∂2u∂x2+C26∂2u∂y2+(C66+C12)∂2u∂x∂y+C66∂2v∂x2+C22∂2v∂y2+2C26∂2v∂x∂y=0.
Utilizing the determinant expansion equation of the differential operator matrix in (2), the characteristic differential equation can be given as
(4)(C11C66-C612)∂4φ∂x4+(C22C66-C622)∂4φ∂y4+(2C11C62-2C61C21)∂4φ∂x3∂y+(2C61C22-2C62C21)∂4φ∂x∂y3+(C11C22-2C66C21-C122+2C61C62)×∂4φ∂x2∂y2=0.Then the general analytical solutions of the displacement can be derived from the characteristic solutionφ of (4):
(5){uv}={C66∂2∂x2+C22∂2∂y2+2C26∂2∂x∂y-(C16∂2∂x2+C26∂2∂y2+(C66+C12)∂2∂x∂y)}φ.
Substituting (5) into (1), the general stress solutions can be derived from the characteristic solutions φ. Table 1 gives seven fundamental analytical stress solutions of the anisotropy plane problem, which can be selected as trial functions of the following numerical method.Table 1
The fundamental stress solutions of the anisotropic plane problem.
Item
1
2
3
4
5
6
7
σ
x
1
0
0
(
C
11
C
66
-
C
61
2
)
(
(
C
22
C
66
-
C
62
2)x+(C21C62-C61C22)y)
(
C
11
C
66
-
C
61
2)((C21C66-C62C61)x +2(C62C11-C21C61)y)
(
C
11
C
66
-
C
61
2)((C61C22-C62C21)x +(C21C66-C22C11+C212-C61C62)y)
(
C
11
C
66
-
C
61
2
)
2
y
σ
y
0
1
0
(
C
21
C
66
-
C
61
C
62)(C66C22-C622)x
(
2
C
66
2
C
21
2
+
C
66
C
21
3
-
4
C
21
C
66
C
61
C
62
-
C
11
C
66
2
C
22
-
C
22
C
21
C
11
C
66
+
C
61
2
C
66
C
22
+
C
11
C
66
C
62
2
+
C
61
2
C
62
2
+
C
62
C
61
C
11
C
22
-
C
62
C
61
C
21
2
)
x
+
(
C
21
C
62
C
11
C
66
-
C
21
C
62
C
61
2
-
C
61
C
22
C
11
C
66
+
C
22
C
61
3
)
y
2
(
C
61
C
62
2
C
21
-
C
66
C
62
C
21
2
+
C
21
C
61
C
66
C
22
-
C
22
C
62
C
61
2)x  + (C22C662C11-C22C66C612-C622C11C66+C622C612)y
(
2
C
61
C
21
2
C
66
-
C
11
C
66
C
62
C
21
-
3
C
21
C
62
C
61
2
-
C
11
C
61
C
66
C
22
+
C
22
C
61
3
+
2
C
61
C
62
2
C
11
)
x +(C62C613-C61C62C11C66+C21C662C11-C21C66C612)y
τ
x
y
0
0
1
-
(
C
11
C
66
-
C
61
2)(C66C22-C622)y
(
C
11
C
66
-
C
61
2
) ((C61C22-C62C21)x +(C61C62-C21C66)y)
(
C
61
2
-
C
11
C
66)((C22C66-C622)x+(C61C22-C21C62)y)
(
C
11
C
66
-
C
61
2
) (C61C62-C21C66)x


## 2.2. The Analytical Trial Functions of the Bending Plate Problem of Anisotropy

The characteristic differential equation of the displacement solutions can be written as [29]
(6)-t312[C11∂4w∂x4+(2C12+4C66)∂4w∂x2∂y2,mm,,+4C61∂4w∂x3∂y+4C62∂4w∂x∂y3+C22∂4w∂y4]+q(x,y)=0.Deriving from the general analytical displacement solutions of (6), Table 2 presents seven general analytical stress solutions of the anisotropy plate problem, which will be selected as trial functions of the following hybrid element in Section 3.Table 2
The fundamental stress solutions of anisotropic material used in bending plate problem.
Item
1
2
3
4
5
6
7
M
x
1
0
0
C
11
x
C
12
y
C
11
y
+
2
C
61
x
C
12
x
+
2
C
61
y
M
y
0
1
0
C
12
x
C
22
y
C
12
y
+
2
C
62
x
C
22
x
+
2
C
62
y
M
x
y
0
0
1
C
61
x
C
62
y
C
61
y
+
2
C
66
x
C
62
x
+
2
C
66
y
Q
x
0
0
0
C
11
C
26
3
C
16
C
12
+
2
C
66
Q
y
0
0
0
C
16
C
22
C
12
+
2
C
66
3
C
26


## 3. The Hybrid Element Strategy Based on the Analytical Trial Functions

### 3.1. The Hybrid Membrane Element Based on the Analytical Trial Functions of Anisotropy

In the sense of novel hybrid element strategy [21], the basic stress solutions of the anisotropic plane problem in Table 1 can be employed as the trial function of the inner stress field in the element, which is denoted as
(7){σ~}={σxσyτxy}=[σx1σx2⋯σxNσy1σy2⋯σyNτxy1τxy2⋯τxyN]{β~1β~2⋮β~N}=[S~]{β~},
where [S~] is matrix of the analytical trial functions, whose dimension is 3×N. N will be determined by the number of trial functions involved in the element. {β~} is the vector of the unknowns of the inner parameters. The modified complementary energy can be presented as [21]
(8)Π~=Π~C-H~PC,
where the complementary energy Π~C is given by
(9)Π~C=12∬Ω{σ~}T[D~]-1{σ~}tdA.
The constitutive matrix [D~] is
(10)[D~]=[C11C12C16C12C22C26C16C26C66].
Substituting (7) into (9) yields
(11)Π~C=12∬Ω{β~}T[S~]T[D~]-1[S~]{β~}tdA=12{β~}T[V~]{β~},
where
(12)[V~]=∬Ω[S~]T[D~]-1[S~]tdA.
The additional complementary energy H~PC is given by
(13)H~PC=∫Γ(Txu~+Tyv~)tds,
in which Tx and Ty  are the boundary forces along the side Γ around the element, determined by the inner stress field in (7):
(14){TxTy}=[l0m0ml]{σxσyτxy}=[L~]{σ~}=[L~][S~]{β~},
where
(15)[L~]=[l0m0ml],
in which
(16)l=sin(θ)=-yi-yjlij,m=-cos(θ)=xi-xjlij.As showed in Figure1, θ is the angle between the side li-j and the axis x and lij is the length of the side li-j.Figure 1
Normal direction and tangential direction.The conforming displacementsu~ and v~ along the side li-j can be obtained from the displacements of u~n and u~s in Figure 1, which can be assumed as
(17)u~n(s)=H1(s)u~ni+H2(s)θ~i+H3(s)u~nj+H4(s)θ~j,u~s(s)=G1(s)u~si+G2(s)u~sj,
where s is the natural coordinate Lj along the side li-j in the triangular element; its value is defined as si=0 and sj=1.H
i
(
s
) are the Hermite functions, which have the expressions as follows:
(18)H1(s)=(1-s)2(2s+1),H2(s)=lij(1-s)2s,H3(s)=s2(3-2s),H4(s)=lijs2(s-1).Gi(s) is the linear shape function
(19)G1(s)=1-s,G2(s)=s.
Involving the relationship between displacement u~n/u~s and displacement u~/v~, the conforming displacement along the side li-j can be defined as
(20)u~(s)=u~i(l2H1(s)+m2G1(s))-v~ilm(G1(s)-H1(s))+θ~ilH2(s)+u~j(l2H3(s)+m2G2(s))-v~jlm(G2(s)-H3(s))+θ~jlH4(s),(21)v~(s)=-u~ilm(G1(s)-H1(s))+v~i(l2G1(s)+m2H1(s))+θ~imH2(s)-u~jlm(G2(s)-H3(s))+v~j(l2G2(s)+m2H3(s))+θ~jmH4(s).
It can also be denoted as
(22){u~v~}=[N~]{q~ij},
where [N~] is the matrix of the trial functions of the displacements defined in side li-j as (23)[N~]=[l2H1+m2G1-lm(G1-H1)lH2l2H3+m2G2-lm(G2-H3)lH4-lm(G1-H1)l2G1+m2H1mH2-lm(G2-H3)l2G2+m2H3mH4].And{q~ij} is part of the unknown nodal displacement vector,
(24){q~ij}=[u~iv~iθ~iu~jv~jθ~j]T.
Substituting (14) and (22) into (13), we have
(25)H~PC=∫T{β}T[S~]T[L~]T[N~]{q~}tds={β}T[H~]{q~},
where
(26)[H~]=∫T[S~]T[L~]T[N~]tds.Thus the modified complementary energyΠ~ can be presented as
(27)Π~=Π~C-H~PC=12{β~}T[V~]{β~}-{β~}T[H~]{q~}.With the principle of the modified complementary energy,δΠ~=0,
(28)∂Π~∂{β~}={0},∂Π~∂{q~}={0}.From the first part of (28), the relationship between the inner parameters {β~} and the nodal displacement {q~} can be determined as
(29){β~}=[V~]-1[H~]{q~}.Substituting (29) into (27),
(30)Π~=-12{q~}T[H~]T[V~]-1[H~]{q~}.The stiffness matrix can be calculated by(31)[K~]=[H~]T[V~]-1[H~].Two triangular membrane elements based on the analytical trial functions of anisotropy are studied in this paper. The one containing five items of basic stress solutions in Table1 is named as ATF-TR5; the other containing seven items is named as ATF-TR7.


### 3.2. The Hybrid Plate Element Based on the Analytical Trial Functions of Anisotropy

The basic stress solutions of the anisotropic plate problem in Table2 can be employed as the trial function of the inner stress field in the plate element, which is denoted as
(32){F}={MxMyMxyQxQy}=[Mx1Mx2⋯MxNMy1My2⋯MyNMxy1Mxy2⋯MxyNQx1Qx2⋯QxNQy1Qy2⋯QyN]{β-1β-2⋮β-N}=[S-]{β-},
where [S-] is matrix of the analytical trial functions, whose dimension is 5×N. N  will be determined by the number of trial functions involved in the element. {β-} is the vector of the unknowns of the inner parameters.The modified complementary energy of the plate problem can be presented as(33)Π-=Π-C-H-PC=∬Ae12{F}T[A-]-1{F}dxdy-∫Se(ψ-nM-n+ψ-sM-ns-w-Q-n)ds=∬Ae12{F}T[A-]-1{F}dxdy-∫Se{F-}T{d-}ds,
where [A-] is the constitutive matrix of the anisotropy plate problem and {F-} is the vector of the boundary forces along the side of the element; it can be denoted as
(34){F-}={M-nM-ns-Q-n}=[L-]{F}=[l2m22lm00-lmlml2-m200000-l-m]×{MxMyMxyQxQy}=[L-][S-]{β-}.{d-} is the conforming displacement; it can be denoted as
(35){d-}={ψ-nψ-sw-}.As showed in Figure2, the conforming displacement along the side li-j can be defined as
(36)w-ij=H1(s)wi+H3(s)wj-H2(s)(xijψxi-yijψyi)lij-H4(s)(xijψxj-yijψyj)lij,ψ-sij=Z2(s)(wj-wi)lij-Z3(s)(xijψxi+yijψyi)lij-Z4(s)(xijψxj+yijψyj)lij,ψ-nij=G1(s)(-yijψxi+xijψyi)lij+G2(s)(-yijψxj+xijψyj)lij,
where
(37)Z2=6s(1-s),Z3=1-4s+3s2,Z4=3s2-2s.
And it is defined as
(38)yij=yi-yj,xij=xi-xj.
Equation (36) can also be denoted as
(39){d-ij}={ψ-nijψ-sijw-ij}=[N-]{q-ij},
where [N-] is the matrix of the trial functions of the displacements defined in side li-j as
(40)[N-]=[0-yijG1lijxijG1lij0-yijG2lijxijG2lij-Z2lij-xijZ3lij-yijZ3lijZ2lij-xijZ4lij-yijZ4lijH1-xijH2lijyijH2lijH3-xijH4lijyijH4lij].
And {q-ij} is part of the unknown nodal displacement vector:
(41){q-ij}=[wiψxiψyiwjψxjψyj]T.Figure 2
The conforming displacement along the side of the element.Following the hybrid element strategy proposed in Section3.1, the modified complementary energy Π- can be presented as
(42)Π-=12{β-}T[V-]{β-}-{β-}T[H-]{q-}.Similarly, the stiffness matrix can be calculated by(43)[K-]=[H-]T[V-]-1[H-].Two triangular plate elements based on the analytical trial functions of anisotropy are studied in this paper. The one containing five items of basic stress solutions in Table2 is named as ATF-TP5; the other containing seven items is named as ATF-TP7.


## 3.1. The Hybrid Membrane Element Based on the Analytical Trial Functions of Anisotropy

In the sense of novel hybrid element strategy [21], the basic stress solutions of the anisotropic plane problem in Table 1 can be employed as the trial function of the inner stress field in the element, which is denoted as
(7){σ~}={σxσyτxy}=[σx1σx2⋯σxNσy1σy2⋯σyNτxy1τxy2⋯τxyN]{β~1β~2⋮β~N}=[S~]{β~},
where [S~] is matrix of the analytical trial functions, whose dimension is 3×N. N will be determined by the number of trial functions involved in the element. {β~} is the vector of the unknowns of the inner parameters. The modified complementary energy can be presented as [21]
(8)Π~=Π~C-H~PC,
where the complementary energy Π~C is given by
(9)Π~C=12∬Ω{σ~}T[D~]-1{σ~}tdA.
The constitutive matrix [D~] is
(10)[D~]=[C11C12C16C12C22C26C16C26C66].
Substituting (7) into (9) yields
(11)Π~C=12∬Ω{β~}T[S~]T[D~]-1[S~]{β~}tdA=12{β~}T[V~]{β~},
where
(12)[V~]=∬Ω[S~]T[D~]-1[S~]tdA.
The additional complementary energy H~PC is given by
(13)H~PC=∫Γ(Txu~+Tyv~)tds,
in which Tx and Ty  are the boundary forces along the side Γ around the element, determined by the inner stress field in (7):
(14){TxTy}=[l0m0ml]{σxσyτxy}=[L~]{σ~}=[L~][S~]{β~},
where
(15)[L~]=[l0m0ml],
in which
(16)l=sin(θ)=-yi-yjlij,m=-cos(θ)=xi-xjlij.As showed in Figure1, θ is the angle between the side li-j and the axis x and lij is the length of the side li-j.Figure 1
Normal direction and tangential direction.The conforming displacementsu~ and v~ along the side li-j can be obtained from the displacements of u~n and u~s in Figure 1, which can be assumed as
(17)u~n(s)=H1(s)u~ni+H2(s)θ~i+H3(s)u~nj+H4(s)θ~j,u~s(s)=G1(s)u~si+G2(s)u~sj,
where s is the natural coordinate Lj along the side li-j in the triangular element; its value is defined as si=0 and sj=1.H
i
(
s
) are the Hermite functions, which have the expressions as follows:
(18)H1(s)=(1-s)2(2s+1),H2(s)=lij(1-s)2s,H3(s)=s2(3-2s),H4(s)=lijs2(s-1).Gi(s) is the linear shape function
(19)G1(s)=1-s,G2(s)=s.
Involving the relationship between displacement u~n/u~s and displacement u~/v~, the conforming displacement along the side li-j can be defined as
(20)u~(s)=u~i(l2H1(s)+m2G1(s))-v~ilm(G1(s)-H1(s))+θ~ilH2(s)+u~j(l2H3(s)+m2G2(s))-v~jlm(G2(s)-H3(s))+θ~jlH4(s),(21)v~(s)=-u~ilm(G1(s)-H1(s))+v~i(l2G1(s)+m2H1(s))+θ~imH2(s)-u~jlm(G2(s)-H3(s))+v~j(l2G2(s)+m2H3(s))+θ~jmH4(s).
It can also be denoted as
(22){u~v~}=[N~]{q~ij},
where [N~] is the matrix of the trial functions of the displacements defined in side li-j as (23)[N~]=[l2H1+m2G1-lm(G1-H1)lH2l2H3+m2G2-lm(G2-H3)lH4-lm(G1-H1)l2G1+m2H1mH2-lm(G2-H3)l2G2+m2H3mH4].And{q~ij} is part of the unknown nodal displacement vector,
(24){q~ij}=[u~iv~iθ~iu~jv~jθ~j]T.
Substituting (14) and (22) into (13), we have
(25)H~PC=∫T{β}T[S~]T[L~]T[N~]{q~}tds={β}T[H~]{q~},
where
(26)[H~]=∫T[S~]T[L~]T[N~]tds.Thus the modified complementary energyΠ~ can be presented as
(27)Π~=Π~C-H~PC=12{β~}T[V~]{β~}-{β~}T[H~]{q~}.With the principle of the modified complementary energy,δΠ~=0,
(28)∂Π~∂{β~}={0},∂Π~∂{q~}={0}.From the first part of (28), the relationship between the inner parameters {β~} and the nodal displacement {q~} can be determined as
(29){β~}=[V~]-1[H~]{q~}.Substituting (29) into (27),
(30)Π~=-12{q~}T[H~]T[V~]-1[H~]{q~}.The stiffness matrix can be calculated by(31)[K~]=[H~]T[V~]-1[H~].Two triangular membrane elements based on the analytical trial functions of anisotropy are studied in this paper. The one containing five items of basic stress solutions in Table1 is named as ATF-TR5; the other containing seven items is named as ATF-TR7.


## 3.2. The Hybrid Plate Element Based on the Analytical Trial Functions of Anisotropy

The basic stress solutions of the anisotropic plate problem in Table2 can be employed as the trial function of the inner stress field in the plate element, which is denoted as
(32){F}={MxMyMxyQxQy}=[Mx1Mx2⋯MxNMy1My2⋯MyNMxy1Mxy2⋯MxyNQx1Qx2⋯QxNQy1Qy2⋯QyN]{β-1β-2⋮β-N}=[S-]{β-},
where [S-] is matrix of the analytical trial functions, whose dimension is 5×N. N  will be determined by the number of trial functions involved in the element. {β-} is the vector of the unknowns of the inner parameters.The modified complementary energy of the plate problem can be presented as(33)Π-=Π-C-H-PC=∬Ae12{F}T[A-]-1{F}dxdy-∫Se(ψ-nM-n+ψ-sM-ns-w-Q-n)ds=∬Ae12{F}T[A-]-1{F}dxdy-∫Se{F-}T{d-}ds,
where [A-] is the constitutive matrix of the anisotropy plate problem and {F-} is the vector of the boundary forces along the side of the element; it can be denoted as
(34){F-}={M-nM-ns-Q-n}=[L-]{F}=[l2m22lm00-lmlml2-m200000-l-m]×{MxMyMxyQxQy}=[L-][S-]{β-}.{d-} is the conforming displacement; it can be denoted as
(35){d-}={ψ-nψ-sw-}.As showed in Figure2, the conforming displacement along the side li-j can be defined as
(36)w-ij=H1(s)wi+H3(s)wj-H2(s)(xijψxi-yijψyi)lij-H4(s)(xijψxj-yijψyj)lij,ψ-sij=Z2(s)(wj-wi)lij-Z3(s)(xijψxi+yijψyi)lij-Z4(s)(xijψxj+yijψyj)lij,ψ-nij=G1(s)(-yijψxi+xijψyi)lij+G2(s)(-yijψxj+xijψyj)lij,
where
(37)Z2=6s(1-s),Z3=1-4s+3s2,Z4=3s2-2s.
And it is defined as
(38)yij=yi-yj,xij=xi-xj.
Equation (36) can also be denoted as
(39){d-ij}={ψ-nijψ-sijw-ij}=[N-]{q-ij},
where [N-] is the matrix of the trial functions of the displacements defined in side li-j as
(40)[N-]=[0-yijG1lijxijG1lij0-yijG2lijxijG2lij-Z2lij-xijZ3lij-yijZ3lijZ2lij-xijZ4lij-yijZ4lijH1-xijH2lijyijH2lijH3-xijH4lijyijH4lij].
And {q-ij} is part of the unknown nodal displacement vector:
(41){q-ij}=[wiψxiψyiwjψxjψyj]T.Figure 2
The conforming displacement along the side of the element.Following the hybrid element strategy proposed in Section3.1, the modified complementary energy Π- can be presented as
(42)Π-=12{β-}T[V-]{β-}-{β-}T[H-]{q-}.Similarly, the stiffness matrix can be calculated by(43)[K-]=[H-]T[V-]-1[H-].Two triangular plate elements based on the analytical trial functions of anisotropy are studied in this paper. The one containing five items of basic stress solutions in Table2 is named as ATF-TP5; the other containing seven items is named as ATF-TP7.


## 4. The Explicit Formulae of the Hybrid Element Based on the Analytical Trial Functions

The triangular flat shell elements based on the analytical trial functions of anisotropy can be presented in the natural coordinates of the triangular element. For example, the matrix[S~] and, [S-] of the stress trial functions contain the variables x  and  y, which can be written as
(44)x=xiLi+xjLj+xmLm,y=yiLi+yjLj+ymLm,
while the matrix [N~] and matrix [N-] of the conforming displacements contain the variables s, which can be denoted as the natural coordinate too, that is in the side li-j(45)s=Lj.According to (16), the matrix [L~] and [L-] can be expressed in the coordinates of the nodes in the element. For example, in the side li-j,
(46)[L~]=[-yijlij0xijlij0xijlij-yijlij],(47)[L]=[(yijlij)2(xijlij)2-2xijyijlij200xijyijlij2-xijyijlij2(yijlij)2-(xijlij)200000yijlij-xijlij].In utilized equations (44)~(47), all variables in the matrices [H~], [V~], [K~], [H-], [V-],  and [K-] can be expressed in the natural coordinates Li, Lj, and Lk. To the integral formula of the natural coordinates in the triangular element, we have
(48)∬ALiaLjbLmcdxdy=a!b!c!(a+b+c+2)!2A,∫lLiaLjbds=a!b!(a+b+1)!lij,
where A is the area of the triangular element. Utilizing the explicit integral formulae of (48) in the triangular element, the explicit formula of the proposed elements in this paper can be obtained.For example, in the plane elements ATF-TR5 and ATF-TR7, we have(49)V~(1,1)=At(C22C66-C262)ΔD,V~(1,2)=-At(C12C66-C16C26)ΔD,V~(2,2)=At(C11C66-C162)ΔD,V~(1,3)=At(C12C26-C16C22)ΔD,V~(2,3)=-At(C11C26-C12C16)ΔD,V~(3,3)=At(C11C22-C122)ΔD,V~(1,4)=t(xj+xm)(C22C66-C262)C663,V~(2,4)=At(yj+ym)(C11C66-C162)C263,ΔD=|D~|=C11C22C66-C11C262-C122C66+2C12C16C26-C162C22.In the plate elements ATF-TP5 and ATF-TP7, we have(50)V-(1,1)=12AtD(C22C66-C262),V-(1,2)=-12AtD(C12C66-C16C26),V-(2,2)=12AtD(C11C66-C162),V-(1,3)=12AtD(C12C26-C16C22),V-(2,3)=-12AtD(C11C26-C12C16),V-(3,3)=12AtD(C11C22-C122),V-(1,4)=4A(xj+xm)t3,V-(2,4)=0,AtD=A/t3ΔD.


## 5. Numerical Examples

### 5.1. Cantilever Shell

As shown in Figure3, a slender cantilever shell has the length b of 32 m, the height h of 2 m, and the thickness t of 0.01 m. It is made of the anisotropy T700, which has the material parameters C11=137.54 GPa, C12=4.14 GPa, C22=9.79 GPa, and C66=4.47 GPa, C44=4.47 GPa, C55=4.47 GPa. There is an angle θ between the principle axis 1 and the axis y, designated as 0°, 45°, and 90°. Figure 4 shows four schemes of the mesh to the proposed cantilever shell.Figure 3
Cantilever beam subjected to force/moment.Schemes of the mesh.
(a)
(b)
(c)
(d)Two load conditions are studied in the example. One loading is the unit in-plane momentMx=1 N·m on the end side C. Table 3 gives the results of the displacement VCz obtained from the proposed elements ATF-TR5 and ATF-TR7, respectively. Comparing with the element S3, which is employed in the software ABAQUS, ATF-TR5 and ATF-TR7 are more accurate.Table 3
DeflectionVCz/10−6 m of the cantilever shell subjected to Mx.
Element
Mesh
θ
2
×
2
2
×
4
2
×
8
2
×
16
S3
0
0.0726
0.1855
0.3035
0.3604
45
0.0126
0.0483
0.2076
1.0908
90
0.0883
0.3400
1.1821
3.1053
ATF-TR5
0
0.4838
0.4814
0.5408
0.6396
45
0.9997
1.7910
3.4030
5.0761
90
5.9177
4.2597
4.3973
5.7013
ATF-TR7
0
0.2864
0.42617
0.5069
0.5627
45
0.2191
0.6552
1.8858
3.7485
90
0.5630
1.6336
3.5687
5.3100The other loading is the unit out-plane loadPx=1 N  on the end point C; Table 4 gives the results of the displacement VCx obtained from the proposed element, ATF-TP5 and ATF-TP7. Comparing with the element S3, ATF-TP5 and ATF-TP7 are more accurate.Table 4
DeflectionVCx/m of the cantilever shell subjected to Px.
Element
Mesh
θ
2
×
2
2
×
4
2
×
8
2
×
16
S3
0
0.4516
0.4748
0.4811
0.4867
45
2.2091
2.8138
3.6630
4.5418
90
7.9508
8.3589
6.7483
8.4965
ATF-TP5
0
0.4793
0.4806
0.4815
0.4818
45
2.3397
3.4545
4.3604
4.8783
90
6.4248
6.6769
6.7493
6.7700
ATF-TP7
0
0.4785
0.4804
0.4814
0.4818
45
2.3199
3.4114
4.3236
4.8473
90
6.3856
6.6638
6.7446
6.7681


### 5.2. Clamped Square Plate

As shown in Figure5, a clamped square plate has the length of 1 m in each side and the thicknesses t of 0.001 m, 0.01 m, and 0.1 m. It is also made of the anisotropy T700 and is subjected to the uniformly distributed loading q=1 Pa. Table 5 gives the deflection wC of the central point C in different mesh schemes.Table 5
DeflectionswC of the central point C.
Thicknesst
Element
Mesh
θ
4
×
4
8
×
8
16
×
16
32
×
32
0.001 mwC/10−3 m
S3
0
0.0010
0.0344
0.1976
0.2370
45
0.0010
0.0292
0.2023
0.3555
ATF-TP5
0
0.2404
0.2381
0.2367
0.2364
45
0.2630
0.3392
0.3519
0.3551
ATF-TP7
0
0.2344
0.2367
0.2363
0.2363
45
0.2526
0.3313
0.3496
0.3546
0.01 mwC/10−6 m
S3
0
0.0423
0.2053
0.2372
0.2413
45
0.0375
0.2231
0.3572
0.3792
ATF-TP5
0
0.2435
0.2412
0.2396
0.2389
45
0.2652
0.3414
0.3542
0.3572
ATF-TP7
0
0.2374
0.2398
0.2392
0.2388
45
0.2549
0.3339
0.3523
0.3568
0.1 mwC/10−9 m
S3
0
0.4583
0.5469
0.5679
0.5732
45
0.4249
0.6214
0.6799
0.6961
ATF-TP5
0
0.4636
0.4548
0.4699
0.4933
45
0.4231
0.5101
0.5727
0.6037
ATF-TP7
0
0.4431
0.4450
0.4612
0.4882
45
0.4136
0.5071
0.5719
0.6035Figure 5
Schemes of the mesh.


### 5.3. Cook’s Skew Beam

As showed in Figure6, Cook’s skew flat shell [30] is studied. The thickness t of the shell is 1 m. It is also made of the anisotropy T700 and is subjected to the uniformly distributed load which has the summation of P=1 N. Table 6 gives the deflection of the central point C in different mesh schemes.Table 6
Deflections (10−9 m) of the central point C in different mesh scheme.
Element
Mesh
θ
2
×
2
4
×
4
8
×
8
16
×
16
32
×
32
S3
0
0.6225
0.8210
1.0105
0.9906
1.1225
45
0.3286
0.4335
0.4932
0.5017
0.5184
90
1.0912
1.7895
2.2273
2.1889
2.3824
ATF-TR5
0
0.8643
1.0041
1.0830
1.1153
1.1282
45
0.5402
0.4987
0.5054
0.5086
0.5103
90
1.4926
1.6447
1.7202
1.7496
1.7643
ATF-TR7
0
0.7102
0.8929
1.0015
1.0483
1.0660
45
0.4576
0.4801
0.4934
0.4974
0.4984
90
1.4055
1.6099
1.7015
1.7410
1.7591Figure 6
Cook’s skew shell.


## 5.1. Cantilever Shell

As shown in Figure3, a slender cantilever shell has the length b of 32 m, the height h of 2 m, and the thickness t of 0.01 m. It is made of the anisotropy T700, which has the material parameters C11=137.54 GPa, C12=4.14 GPa, C22=9.79 GPa, and C66=4.47 GPa, C44=4.47 GPa, C55=4.47 GPa. There is an angle θ between the principle axis 1 and the axis y, designated as 0°, 45°, and 90°. Figure 4 shows four schemes of the mesh to the proposed cantilever shell.Figure 3
Cantilever beam subjected to force/moment.Schemes of the mesh.
(a)
(b)
(c)
(d)Two load conditions are studied in the example. One loading is the unit in-plane momentMx=1 N·m on the end side C. Table 3 gives the results of the displacement VCz obtained from the proposed elements ATF-TR5 and ATF-TR7, respectively. Comparing with the element S3, which is employed in the software ABAQUS, ATF-TR5 and ATF-TR7 are more accurate.Table 3
DeflectionVCz/10−6 m of the cantilever shell subjected to Mx.
Element
Mesh
θ
2
×
2
2
×
4
2
×
8
2
×
16
S3
0
0.0726
0.1855
0.3035
0.3604
45
0.0126
0.0483
0.2076
1.0908
90
0.0883
0.3400
1.1821
3.1053
ATF-TR5
0
0.4838
0.4814
0.5408
0.6396
45
0.9997
1.7910
3.4030
5.0761
90
5.9177
4.2597
4.3973
5.7013
ATF-TR7
0
0.2864
0.42617
0.5069
0.5627
45
0.2191
0.6552
1.8858
3.7485
90
0.5630
1.6336
3.5687
5.3100The other loading is the unit out-plane loadPx=1 N  on the end point C; Table 4 gives the results of the displacement VCx obtained from the proposed element, ATF-TP5 and ATF-TP7. Comparing with the element S3, ATF-TP5 and ATF-TP7 are more accurate.Table 4
DeflectionVCx/m of the cantilever shell subjected to Px.
Element
Mesh
θ
2
×
2
2
×
4
2
×
8
2
×
16
S3
0
0.4516
0.4748
0.4811
0.4867
45
2.2091
2.8138
3.6630
4.5418
90
7.9508
8.3589
6.7483
8.4965
ATF-TP5
0
0.4793
0.4806
0.4815
0.4818
45
2.3397
3.4545
4.3604
4.8783
90
6.4248
6.6769
6.7493
6.7700
ATF-TP7
0
0.4785
0.4804
0.4814
0.4818
45
2.3199
3.4114
4.3236
4.8473
90
6.3856
6.6638
6.7446
6.7681


## 5.2. Clamped Square Plate

As shown in Figure5, a clamped square plate has the length of 1 m in each side and the thicknesses t of 0.001 m, 0.01 m, and 0.1 m. It is also made of the anisotropy T700 and is subjected to the uniformly distributed loading q=1 Pa. Table 5 gives the deflection wC of the central point C in different mesh schemes.Table 5
DeflectionswC of the central point C.
Thicknesst
Element
Mesh
θ
4
×
4
8
×
8
16
×
16
32
×
32
0.001 mwC/10−3 m
S3
0
0.0010
0.0344
0.1976
0.2370
45
0.0010
0.0292
0.2023
0.3555
ATF-TP5
0
0.2404
0.2381
0.2367
0.2364
45
0.2630
0.3392
0.3519
0.3551
ATF-TP7
0
0.2344
0.2367
0.2363
0.2363
45
0.2526
0.3313
0.3496
0.3546
0.01 mwC/10−6 m
S3
0
0.0423
0.2053
0.2372
0.2413
45
0.0375
0.2231
0.3572
0.3792
ATF-TP5
0
0.2435
0.2412
0.2396
0.2389
45
0.2652
0.3414
0.3542
0.3572
ATF-TP7
0
0.2374
0.2398
0.2392
0.2388
45
0.2549
0.3339
0.3523
0.3568
0.1 mwC/10−9 m
S3
0
0.4583
0.5469
0.5679
0.5732
45
0.4249
0.6214
0.6799
0.6961
ATF-TP5
0
0.4636
0.4548
0.4699
0.4933
45
0.4231
0.5101
0.5727
0.6037
ATF-TP7
0
0.4431
0.4450
0.4612
0.4882
45
0.4136
0.5071
0.5719
0.6035Figure 5
Schemes of the mesh.


## 5.3. Cook’s Skew Beam

As showed in Figure6, Cook’s skew flat shell [30] is studied. The thickness t of the shell is 1 m. It is also made of the anisotropy T700 and is subjected to the uniformly distributed load which has the summation of P=1 N. Table 6 gives the deflection of the central point C in different mesh schemes.Table 6
Deflections (10−9 m) of the central point C in different mesh scheme.
Element
Mesh
θ
2
×
2
4
×
4
8
×
8
16
×
16
32
×
32
S3
0
0.6225
0.8210
1.0105
0.9906
1.1225
45
0.3286
0.4335
0.4932
0.5017
0.5184
90
1.0912
1.7895
2.2273
2.1889
2.3824
ATF-TR5
0
0.8643
1.0041
1.0830
1.1153
1.1282
45
0.5402
0.4987
0.5054
0.5086
0.5103
90
1.4926
1.6447
1.7202
1.7496
1.7643
ATF-TR7
0
0.7102
0.8929
1.0015
1.0483
1.0660
45
0.4576
0.4801
0.4934
0.4974
0.4984
90
1.4055
1.6099
1.7015
1.7410
1.7591Figure 6
Cook’s skew shell.


## 6. Conclusions

This paper presents a novel way to formulate the flat shell element of anisotropy. The basic analytical solutions and the proposed conforming displacement are employed to formulate the triangular elements ATF-TR5, ATF-TR7, ATF-TP5, and ATF-TP7. The high accuracy and computational efficiency of these elements are proved by the benchmark.



---

*Source: 486453-2013-09-16.xml*