# A Conforming Triangular Plane Element with Rotational Degrees of Freedom

Authors: Xiang-Rong Fu; Ming-Wu Yuan; Chen Pu
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

This paper presents a novel way to formulate the triangular plane element with rotational degrees of freedom (RDOF). The linear distribution of rotational displacement is assumed. The conforming displacement along the sides based on the rotational displacement assumption is derived, and the triangular plane element TR3 for isotropic material is formulated. By using the explicit integral formulae of the triangular element, the matrices used in the proposed plane element TR3 are calculated efficiently. The benchmark examples showed thier high accuracy and high efficiency.


---

## Body

## 1. Introduction

The triangular plane element is widely used [1–3]. In order to improve the performance of the triangular plane element, rotational degrees of freedom (RDOF) has been introduced. The other advantage of plane elements with RDOF is that the singularity of global stiffness can be avoided in the analysis of shell structures. Mervyn and Terrence introduced the RDOF to describe the displacement of the flat shell element [4]. Mohr constructed the hybrid membrane element with RDOF [5]. Allman et al. presented some rational displacement-type plane elements with RDOF [6–10]. Many researchers improve accuracy of the membrane element by this way [11–15], and the plane elements with RDOF have been widely used in analysis of shell structure.Long et al. studied the generalized conforming plane elements [16–20] with RDOF. Taking advantage of the conforming boundary displacement, the nonconforming displacement model can be involved freely by the generalized conforming equations [17–20]. Chen et al. proved that the hybrid element using the balanced stress fields can play the same role in the element formulation as the generalized conforming element introducing the associated generalized conforming equations [20]. On the other hand, the hybrid element employed analytical trial functions of stress needs conforming boundary displacement [21–26].This paper studies the triangular plane element based on the rational assumption of the rotation displacement field. There are mainly four steps as follows in Section2.(a)
The first step assumes the rational distribution of the rotational displacement (in Section2.1).(b)
The second step gives the displacement field based on the assumed rotational displacement (in Section2.2).(c)
The third step formulates the stiffness matrix of the element TR3 based on the associated strain fields (in Section2.3).(d)
The fourth step calculates the stiffness matrix using the explicit integral formulae of the triangular element (in Section2.4).In Section3, some numerical examples are shown as benchmark to study the accuracy and the efficiency of the proposed element model TR3.In Section4, some conclusions are given.


## 2. The Formulation of the Triangular Element TR3 with RDOF

### 2.1. The Rational Distribution of the Rotational Displacement

The rotational displacement field in the triangular element with three nodes can be assumed as follows:(1)θz=L1θZ1+L2θZ2+L3θZ3,where(2)L1=a1+b1x+c1y2A,L2=a2+b2x+c2y2A,L3=a3+b3x+c3y2A,where(3)a1=x2y3-x3y2,a2=x3y1-x1y3,a3=x1y2-x2y1,b1=y2-y3,b2=y3-y1,b3=y1-y2,c1=x3-x2,c2=x1-x3,c3=x2-x1,2A=x1y2-x2y1-x1y3+x3y1+x2y3-x3y2,in which x1, x2, and x3 and y1, y2, and y3 are coordinates of three nodes in the triangular element. θZ1, θZ2, and θZ3 are rotational displacement of three nodes.


### 2.2. The Distribution of the Displacement

According to the definition of the rotation displacement in the elasticity, the relationship between displacement symbolsu, v and the rotational displacement θZ can be denoted as(4)θz=12∂u∂y-∂v∂x.From (1),(5)θz=a1θz1+a2θz2+a3θz3+b1θz1+b2θz2+b3θz3x+c1θz1+c2θz2+c3θz3y=Az+Bzx+Czy.From (4) and (5), we can assume(6)ux,y=Az2+Bzxy2+Cz4y2+fux,y,vx,y=-Az2-Bz4x2-Czxy2+fvx,y,where we can give two assumptions,(7)fux,y=u0+uxx+uyy,fvx,y=v0+vxx+vyy.On the other hand, we can define the displacement fields in the triangular element as(8)u=guu1,u2,u3,θZ1,θZ2,θZ3,x,y,v=gvv1,v2,v3,θZ1,θZ2,θZ3,x,y.According to (6)–(8), the nodal conforming requirement can be given as(9)u1=ux1,y1,u2=ux2,y2,u3=ux3,y3,v1=vx1,y1,v2=vx2,y2,v3=vx3,y3.Then, we can find the relationship between nodal displacement symbolsu1, u2, u3, v1, v2, and v3 and inner parameters u0, ux, uy, v0, vx, and vy; it is a linear translation matrix. By solving the associated equations (9), we can get the distribution of the displacement as follows:(10)u=u1L1+u2L2+u3L3+θz1a12y+c14y2+b12xy+θz2a22y+c24y2+b22xy+θz3a32y+c34y2+b32xy+θz142Ay1-y22xy2-y32x1y1y3-y211111111111111111111.1+x2y1+y2y1-y3111111111111111111111+x3y1+y3y2-y1+θz142Ah+yx2-x3+θz142Ah·x1-x2y1-y32+x3-x1y1-y22+θz142Ah-x22x1x2y32y1-y2y2-y3+θz142Ah+x1x3y2y3-2y1y2-y3+θz142Ah+x22y1y3y3-y1+x32y1y2y2-y1+θz142Ah+x2x3y1y1+y2+y3y2+y3+θz142Ah+x2xhh3y1-6y2y3x22y1-y22+θz242Ay1-y22y2-y12xy3-y12x2y2y1-y32x2y2y1-ywwee.eW+x3y2+y3y2-y12x2y2y1-ywweeeW+x1y2+y1y3-y2+θz242Ah+yx3-x1+θz242Ah·x2-x3y2-y12+x1-x2y2-y32+θz242Ah-x32x2x3y12y2-y3y3-y1+θz242Ah+x1x2y3y1-2y2y3-y1+θz242Ah+x32y2y1y1-y2+x12y2y3y3-y2+θz242Ah+x3x1y2y1+y2+y3y1+y3hhhhhhhhhhhhhhhhhhhh-6y1y3x32y2-y12+θz342Ay3-y22xy1-y22x3y3y2-y12x3y3y2-ywwee.eW+x1y3+y1y3-y22x3y3y2-ywweeeW+x2y3+y2y1-y3+θz342Ah+yx1-x2+θz342Ah·[x3-x1y3-y22+x2-x3y3-y12]+θz342Ah-x12x3x1y22y3-y1y1-y2+θz342Ah+x3x2y1y2-2y3y1-y2+θz342Ah+x12y3y2y2-y3+x22y3y1y1-y3+θz342Ah+x2x1y3y1+y2+y3y1+y2x2x1y3y1+y2+y301-6y1y2x12y3-y22,v=v1L1+v2L2+v3L3+θz1a12x+b14x2+c12xy-θz2a22x+b24x2+c22xy-θz3a32x+b34x2+c32xy+θz142Ax1-x32-yx2-x32x1y1x3-x22x1y1x3-xwweee.eWy2+y2x1+x2x1-x32x1y1x3-xwweeeeWy2+y3x1+x3x2-x1+θz142Ah-xy2-y3+θz142Ah·y1-y2x1-x32+y3-y1x1-x22+θz142Ah+y22y1y2x32x1-x2x2-x3+θz142Ah+y1y3x2x3-2x1x2-x3+θz142Ah+y22x1x3x3-x1+y32x1x2x2-x1+θz142Ah+y2y3x1x1+x2+x3x2+x3-6x2x3y22x1-x32+θz242Ax2-x12-yx3-x12y2x2x1-x32y2x2x1-xwweeee.Wy3+y3x2+x3x2-x12y2x2x1-xwweeeeWy3+y1x2+x1x3-x2+θz242Ah-xy3-y1+θz242Ah·[y2-y3x2-x12+y1-y2x2-x32]+θz242Ah+y32y2y3x12x2-x3x3-x1+θz242Ah+y1y2x3x1-2x2x3-x1+θz242Ah+y32x2x1x1-x2+y12x2x3x3-x2+θz242Ah+y3y1x2x1+x2+x3x1+x3-6x1x3y32x2-x12+θz342Ax3-x12-yx1-x22y3x3x2-x12y3x3x2-xwweee.eWy2+y1x3+x1x3-x22y3x3x2-xwweeeeWy2+y2x3+x2x1-x3+θz342Ah-xy1-y2+θz342Ah·[y3-y1x3-x22+y2-y3x3-x12]+θz342Ah+y12y3y1x22x3-x1x1-x2+θz342Ah+y3y2x1x2-2x3x1-x2+θz342Ah+y12x3x2x2-x3+y22x3x1x1-x3+θz342Ah+y2y1x3x1+x2+x3x1+x2-6x1x2y12x3-x12.It can be proved that the assumed displacement is conforming.


### 2.3. The Stiffness Matrix of the Plane Element TR3 with RDOF

From (10), taking into account the definition of the strain field, we can get the strain matrix about the nodal displacement. It can be denoted as(11)B=B1B2B3,B1=b12A0b12y-y1-bx40c12A-c12x-x1-by4c12Ab12Ax1b1+x2b3+x3b24,B2=b22A0b22y-y2-bx40c22A-c22x-x2-by4c22Ab22Ax1b3+x2b2+x3b24,B3=b32A0b32y-y3-bx40c32A-c32x-x3-by4c32Ab32Ax1b2+x2b1+x3b34,in which(12)bx=x1y1b1+x2y2b2+x3y3b32A,by=x1y1c1+x2y2c2+x3y3c32A.Then, the stiffness matrix can be obtained:(13)K=∬ABTDBdA,where(14)D=E1-μ21μ0μ10001-μ2.


### 2.4. The Stiffness Matrix Using the Explicit Integral Formulae

The strain matrix of the triangular plane element TR3 with RDOF can be presented in the natural coordinates of the triangular element. The stress trial functions contain the variablesx and y, which can be written as(15)x=x1L1+x2L2+x3L3,y=y1L1+y2L2+y3L3.Then, the matrixK can be expressed in the natural coordinates L1, L2, and L3. To the integral formula of the natural coordinates in the triangular element, we have(16)∬AL1aL2bL3cdxdy=a!b!c!a+b+c+2!2A,∫lL1aL2bds=a!b!a+b+1!l12,where A is the area of the triangular element. Utilizing the explicit integral formulae of (16) in the triangular element, the explicit formula of the proposed elements in this paper can be obtained.


## 2.1. The Rational Distribution of the Rotational Displacement

The rotational displacement field in the triangular element with three nodes can be assumed as follows:(1)θz=L1θZ1+L2θZ2+L3θZ3,where(2)L1=a1+b1x+c1y2A,L2=a2+b2x+c2y2A,L3=a3+b3x+c3y2A,where(3)a1=x2y3-x3y2,a2=x3y1-x1y3,a3=x1y2-x2y1,b1=y2-y3,b2=y3-y1,b3=y1-y2,c1=x3-x2,c2=x1-x3,c3=x2-x1,2A=x1y2-x2y1-x1y3+x3y1+x2y3-x3y2,in which x1, x2, and x3 and y1, y2, and y3 are coordinates of three nodes in the triangular element. θZ1, θZ2, and θZ3 are rotational displacement of three nodes.


## 2.2. The Distribution of the Displacement

According to the definition of the rotation displacement in the elasticity, the relationship between displacement symbolsu, v and the rotational displacement θZ can be denoted as(4)θz=12∂u∂y-∂v∂x.From (1),(5)θz=a1θz1+a2θz2+a3θz3+b1θz1+b2θz2+b3θz3x+c1θz1+c2θz2+c3θz3y=Az+Bzx+Czy.From (4) and (5), we can assume(6)ux,y=Az2+Bzxy2+Cz4y2+fux,y,vx,y=-Az2-Bz4x2-Czxy2+fvx,y,where we can give two assumptions,(7)fux,y=u0+uxx+uyy,fvx,y=v0+vxx+vyy.On the other hand, we can define the displacement fields in the triangular element as(8)u=guu1,u2,u3,θZ1,θZ2,θZ3,x,y,v=gvv1,v2,v3,θZ1,θZ2,θZ3,x,y.According to (6)–(8), the nodal conforming requirement can be given as(9)u1=ux1,y1,u2=ux2,y2,u3=ux3,y3,v1=vx1,y1,v2=vx2,y2,v3=vx3,y3.Then, we can find the relationship between nodal displacement symbolsu1, u2, u3, v1, v2, and v3 and inner parameters u0, ux, uy, v0, vx, and vy; it is a linear translation matrix. By solving the associated equations (9), we can get the distribution of the displacement as follows:(10)u=u1L1+u2L2+u3L3+θz1a12y+c14y2+b12xy+θz2a22y+c24y2+b22xy+θz3a32y+c34y2+b32xy+θz142Ay1-y22xy2-y32x1y1y3-y211111111111111111111.1+x2y1+y2y1-y3111111111111111111111+x3y1+y3y2-y1+θz142Ah+yx2-x3+θz142Ah·x1-x2y1-y32+x3-x1y1-y22+θz142Ah-x22x1x2y32y1-y2y2-y3+θz142Ah+x1x3y2y3-2y1y2-y3+θz142Ah+x22y1y3y3-y1+x32y1y2y2-y1+θz142Ah+x2x3y1y1+y2+y3y2+y3+θz142Ah+x2xhh3y1-6y2y3x22y1-y22+θz242Ay1-y22y2-y12xy3-y12x2y2y1-y32x2y2y1-ywwee.eW+x3y2+y3y2-y12x2y2y1-ywweeeW+x1y2+y1y3-y2+θz242Ah+yx3-x1+θz242Ah·x2-x3y2-y12+x1-x2y2-y32+θz242Ah-x32x2x3y12y2-y3y3-y1+θz242Ah+x1x2y3y1-2y2y3-y1+θz242Ah+x32y2y1y1-y2+x12y2y3y3-y2+θz242Ah+x3x1y2y1+y2+y3y1+y3hhhhhhhhhhhhhhhhhhhh-6y1y3x32y2-y12+θz342Ay3-y22xy1-y22x3y3y2-y12x3y3y2-ywwee.eW+x1y3+y1y3-y22x3y3y2-ywweeeW+x2y3+y2y1-y3+θz342Ah+yx1-x2+θz342Ah·[x3-x1y3-y22+x2-x3y3-y12]+θz342Ah-x12x3x1y22y3-y1y1-y2+θz342Ah+x3x2y1y2-2y3y1-y2+θz342Ah+x12y3y2y2-y3+x22y3y1y1-y3+θz342Ah+x2x1y3y1+y2+y3y1+y2x2x1y3y1+y2+y301-6y1y2x12y3-y22,v=v1L1+v2L2+v3L3+θz1a12x+b14x2+c12xy-θz2a22x+b24x2+c22xy-θz3a32x+b34x2+c32xy+θz142Ax1-x32-yx2-x32x1y1x3-x22x1y1x3-xwweee.eWy2+y2x1+x2x1-x32x1y1x3-xwweeeeWy2+y3x1+x3x2-x1+θz142Ah-xy2-y3+θz142Ah·y1-y2x1-x32+y3-y1x1-x22+θz142Ah+y22y1y2x32x1-x2x2-x3+θz142Ah+y1y3x2x3-2x1x2-x3+θz142Ah+y22x1x3x3-x1+y32x1x2x2-x1+θz142Ah+y2y3x1x1+x2+x3x2+x3-6x2x3y22x1-x32+θz242Ax2-x12-yx3-x12y2x2x1-x32y2x2x1-xwweeee.Wy3+y3x2+x3x2-x12y2x2x1-xwweeeeWy3+y1x2+x1x3-x2+θz242Ah-xy3-y1+θz242Ah·[y2-y3x2-x12+y1-y2x2-x32]+θz242Ah+y32y2y3x12x2-x3x3-x1+θz242Ah+y1y2x3x1-2x2x3-x1+θz242Ah+y32x2x1x1-x2+y12x2x3x3-x2+θz242Ah+y3y1x2x1+x2+x3x1+x3-6x1x3y32x2-x12+θz342Ax3-x12-yx1-x22y3x3x2-x12y3x3x2-xwweee.eWy2+y1x3+x1x3-x22y3x3x2-xwweeeeWy2+y2x3+x2x1-x3+θz342Ah-xy1-y2+θz342Ah·[y3-y1x3-x22+y2-y3x3-x12]+θz342Ah+y12y3y1x22x3-x1x1-x2+θz342Ah+y3y2x1x2-2x3x1-x2+θz342Ah+y12x3x2x2-x3+y22x3x1x1-x3+θz342Ah+y2y1x3x1+x2+x3x1+x2-6x1x2y12x3-x12.It can be proved that the assumed displacement is conforming.


## 2.3. The Stiffness Matrix of the Plane Element TR3 with RDOF

From (10), taking into account the definition of the strain field, we can get the strain matrix about the nodal displacement. It can be denoted as(11)B=B1B2B3,B1=b12A0b12y-y1-bx40c12A-c12x-x1-by4c12Ab12Ax1b1+x2b3+x3b24,B2=b22A0b22y-y2-bx40c22A-c22x-x2-by4c22Ab22Ax1b3+x2b2+x3b24,B3=b32A0b32y-y3-bx40c32A-c32x-x3-by4c32Ab32Ax1b2+x2b1+x3b34,in which(12)bx=x1y1b1+x2y2b2+x3y3b32A,by=x1y1c1+x2y2c2+x3y3c32A.Then, the stiffness matrix can be obtained:(13)K=∬ABTDBdA,where(14)D=E1-μ21μ0μ10001-μ2.


## 2.4. The Stiffness Matrix Using the Explicit Integral Formulae

The strain matrix of the triangular plane element TR3 with RDOF can be presented in the natural coordinates of the triangular element. The stress trial functions contain the variablesx and y, which can be written as(15)x=x1L1+x2L2+x3L3,y=y1L1+y2L2+y3L3.Then, the matrixK can be expressed in the natural coordinates L1, L2, and L3. To the integral formula of the natural coordinates in the triangular element, we have(16)∬AL1aL2bL3cdxdy=a!b!c!a+b+c+2!2A,∫lL1aL2bds=a!b!a+b+1!l12,where A is the area of the triangular element. Utilizing the explicit integral formulae of (16) in the triangular element, the explicit formula of the proposed elements in this paper can be obtained.


## 3. Benchmark of the Proposed Element TR3

### 3.1. Cantilever Beam

As shown in Figure1, a slender cantilever beam has the length b of 32 m, the height h of 2 m, and the thickness t of 1 m. It is made of the isotropy material, which has the material parameters E=768Pa, µ=0.25. Figure 2 shows four schemes of the mesh to the proposed cantilever beam.Figure 1
Cantilever beam subjected to force/moment.Figure 2
Schemes of the mesh.The load is an in-plane momentMx = 1 N·m on the end side C. Table 1 gives the results of the displacement VCz obtained by the proposed element TR3. Comparing with other elements, TR3 can give higher precise results.Table 1
DeflectionVCz/m of the cantilever beam subjected to Mx.
Element
Mesh
2 × 2
2 × 4
2 × 8
2 × 16
CST
0.0128
0.0489
0.1654
0.4058
GT9
0.1045
0.3156
0.6371
0.8510
GT9M
0.1049
0.3191
0.6504
0.8692
TR3
0.9906
0.9903
0.9911
0.9961


### 3.2. Cook’s Skew Beam

As showed in Figure3, the Cook’s skew beam [10] is studied. The thickness t of the shell is 1 m. It is also made of the isotropy, whose material parameters E = 1 Pa, µ = 0.3333, and is subjected to the uniformly distributed load which has the sum of P = 1 N. Table 2 gives the deflection of the central point C in different mesh scheme.Table 2
Deflections (m) of the central pointC in different mesh scheme.
Element
Mesh
2 × 2
4 × 4
8 × 8
16 × 16
32 × 32
64 × 64
CST
11.99
18.28
22.02
23.41
/
/
GT9
20.08
22.71
23.61
23.87
23.94
23.96
GT9M
20.36
22.80
23.63
23.88
23.94
23.96
TR3
19.66
22.46
23.47
23.81
23.92
23.95Figure 3
Cook’s skew shell.


## 3.1. Cantilever Beam

As shown in Figure1, a slender cantilever beam has the length b of 32 m, the height h of 2 m, and the thickness t of 1 m. It is made of the isotropy material, which has the material parameters E=768Pa, µ=0.25. Figure 2 shows four schemes of the mesh to the proposed cantilever beam.Figure 1
Cantilever beam subjected to force/moment.Figure 2
Schemes of the mesh.The load is an in-plane momentMx = 1 N·m on the end side C. Table 1 gives the results of the displacement VCz obtained by the proposed element TR3. Comparing with other elements, TR3 can give higher precise results.Table 1
DeflectionVCz/m of the cantilever beam subjected to Mx.
Element
Mesh
2 × 2
2 × 4
2 × 8
2 × 16
CST
0.0128
0.0489
0.1654
0.4058
GT9
0.1045
0.3156
0.6371
0.8510
GT9M
0.1049
0.3191
0.6504
0.8692
TR3
0.9906
0.9903
0.9911
0.9961


## 3.2. Cook’s Skew Beam

As showed in Figure3, the Cook’s skew beam [10] is studied. The thickness t of the shell is 1 m. It is also made of the isotropy, whose material parameters E = 1 Pa, µ = 0.3333, and is subjected to the uniformly distributed load which has the sum of P = 1 N. Table 2 gives the deflection of the central point C in different mesh scheme.Table 2
Deflections (m) of the central pointC in different mesh scheme.
Element
Mesh
2 × 2
4 × 4
8 × 8
16 × 16
32 × 32
64 × 64
CST
11.99
18.28
22.02
23.41
/
/
GT9
20.08
22.71
23.61
23.87
23.94
23.96
GT9M
20.36
22.80
23.63
23.88
23.94
23.96
TR3
19.66
22.46
23.47
23.81
23.92
23.95Figure 3
Cook’s skew shell.


## 4. Conclusion

This paper present a new type of plane element with RDOF: the assumptions of the rotational displacement and the displacement are proved to be useful. The accuracy and the efficiency of the proposed element model TR3 are shown in the benchmark problems.



---

*Source: 274709-2015-05-20.xml*