# Four-Node Generalized Conforming Membrane Elements with Drilling DOFs Using Quadrilateral Area Coordinate Methods

Authors: Xiao-Ming Chen; Song Cen; Jian-Yun  Sun; Yun-Gui Li
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Two 4-node generalized conforming quadrilateral membrane elements with drilling DOF, named QAC4θ and QAC4θM, were successfully developed. Two kinds of quadrilateral area coordinates are used together in the assumed displacement fields of the new elements, so that the related formulations are quite straightforward and will keep the order of the Cartesian coordinates unchangeable while the mesh is distorted. The drilling DOF is defined as the additional rigid rotation at the element nodes to avoid improper constraint. Both elements can pass the strict patch test and exhibit better performance than other similar models. In particular, they are both free of trapezoidal locking in MacNeal’s beam test and insensitive to various mesh distortions.


---

## Body

## 1. Introduction

It is well known that adding the drilling degree of freedom (DOF) at each node of a plane membrane element can enhance the element performance without increasing the number of the element nodes. Furthermore, a plane membrane element with drilling DOFs can be combined with a plate bending element to form a flat-shell element, which can avoid the problem of singular coefficients associated with the DOFs in the direction normal to the plane of the shell element.The research on drilling DOF can be traced to the 1960s. Olson and Bearden [1] proposed the first valuable triangular membrane element with drilling DOF. However, since the rotations of two adjacent edges are assumed to be equal, this element may not converge to the correct solution. Another model, a hybrid displacement triangular element with drilling DOF, was then proposed by Mohr [2], but its variational theory is not sufficient. The definition of the drilling DOFs proposed by Allman [3, 4] can be treated as a milestone for this topic, in which a quadratic displacement approximation was introduced to supplement the drilling DOFs at element nodes. Based on Allman’s work, numerous researches on plane elements with drilling degrees of freedom have been accomplished, such as the models proposed by Bergan and Felippa [5], Cook [6, 7], MacNeal and Harder [8], Yunus et al. [9], Hughes and Brezzi [10], Ibrahimbegovic et al. [11], Iura and Atluri [12], Cazzani and Atluri [13], Piltner and Taylor [14], Geyer and Groenwold [15], Pimpinelli [16], Groenwold et al. [17], Choi et al. [18], Choo et al. [19], Zhang and Kuang [20], Kugler et al. [21], and Cen et al. [22]. Long et al. [23–25] presented a new definition on the drilling DOFs. They treated these DOFs as the additional rigid rotations at the element nodes, so that the change of the angle between two adjacent edges along with the element deformation is allowed, and the rotation of the element edge has definite relation with the nodal drilling DOFs. Based on this assumption, the triangular element GT9 series and the quadrilateral element GQ12 series were developed. These elements all exhibit good performance [26], and among these models the quadrilateral element GQ12M8 is the best one.It is also well known that most quadrilateral elements use the isoparametric coordinates (ξ, η) to express their formulations. Lee and Bathe [27] have studied the influence of mesh distortions on the isoparametric membrane elements and showed that the serendipity family is quite sensitive to the mesh distortions. They concluded that the nonlinear transformation between the isoparametric (local) and the Cartesian (global) coordinates leads to such problem. Although the assumed displacement fields may contain high-order terms of ξ and η, their complete order in Cartesian coordinates x and y will degrade significantly once the meshes are distorted, which will lead to low accuracy. In order to make the isoparametric displacement fields satisfy second order completeness in Cartesian coordinates, even fourth order isoparametric terms should be introduced, such as the abovementioned element GQ12M8. This makes the element formulations quite complicated.For overcoming this inherent defect of the isoparametric coordinates, Long et al. [28, 29] developed the first kind of quadrilateral area coordinate method QACM-I. The QACM-I possesses an important character: the transformation between the area and the Cartesian coordinate systems is always linear. Thus, the disadvantage of the isoparametric coordinate system is avoided from the outset. Then, this new natural coordinate system was successfully applied to develop new finite element models. Soh et al. [30] constructed two 8-node plane quadrilateral generalized conforming elements which are insensitive to mesh distortion. Chen et al. [31] proposed two 4-node quadrilateral membrane elements AGQ6-I and AGQ6-II, which exhibit excellent performance in high-order benchmark examples; particularly, both can perfectly pass MacNeal’s thin beam test [32]. These two 4-node elements arouse the interests in further studies on the QACM-I. Cen et al. [33] derived out the analytical element stiffness matrix of AGQ6-I and developed a family of the quadrilateral plane membrane elements [34]. Du and Cen [35] extended the element AGQ6-I to geometrically nonlinear analysis. Cardoso et al. [36–38] introduced the element AGQ6-I to develop distortion-immune shell elements for linear, nonlinear, and dynamic fracture analyses. Wang and Sun [39] used the element AGQ6-II to formulate a new corotational nonlinear shell element. Chen et al. [40] modified the element AGQ6-I to make it pass the strict patch test. Li [41] improved the formulations and generalized them to simulate coupled solid-deformation/fluid-flow for porous geomaterials. Cardoso and Yoon [42], Prathap and Senthilkumar [43], and Flajs et al. [44] discussed the convergence for related AGQ6 models. Besides the above plane elements, the QACM-I has also been successfully employed to develop thin plate [45], Mindlin-Reissner plate [46], laminated composite plate [47], and shell models [48–51].Since the QACM-I contains four area coordinate components (L1, L2, L3, and L4), among which only two are independent, users may be confused on how to formulate a complete high-order polynomial. In view of this disadvantage, Chen et al. [52] proposed the second kind of quadrilateral coordinate method QACM-II. This QACM-II uses two midlines of opposite sides as the coordinate axes and defines only two independent coordinate components Z1 and Z2. The element formulations expressed by the QACM-II are quite simpler than those in terms of the QACM-I [52, 53]. In 2010, Long et al. [54] established the third kind of quadrilateral area coordinate method QACM-III. It takes two diagonals as the zero coordinate axes to define the two independent coordinate components T1 and T2. All the three kinds of area quadrilateral coordinate can be used simultaneously in one element, which will make the formulations quite simple and straightforward.In this paper, by combination with the definition of drilling DOFs proposed by Long et al. [23–25], a new plane membrane element with drilling DOFs, denoted by QAC4θ, was firstly developed by using the QACM-III. Then, by introducing a generalized bubble displacement field in terms of QACM-II into the element QAC4θ, a more accurate and robust element, denoted by QAC4θM, was constructed. Both elements can pass the strict patch test and exhibit better performance than other similar models. It is demonstrated again that the quadrilateral area coordinate methods are effective tools for developing high-performance quadrilateral finite element models.


## 2. Brief Reviews on the Quadrilateral Area Coordinate Methods

### 2.1. QACM-I [28, 29]

As shown in Figure1, the position of an arbitrary point P within a quadrilateral element 1234¯ is specified by the area coordinates L1, L2, L3, and L4, which are defined as(1)Li=AiA,i=1,2,3,4,where A is the area of the quadrilateral element; Ai (i=1,2,3,4) are the areas of the four triangles constructed by point P and four element sides 23¯, 34¯, 41¯, and 12¯, respectively. L1, L2, L3, and L4 can be expressed in terms of Cartesian coordinates (x,y) as follows:(2)Li=12Aai+bix+ciy,i=1,2,3,4with(3)ai=xjyk-xkyj,bi=yj-yk,ci=xk-xi,i=1,2,3,4;j=2,3,4,1;k=3,4,1,2.Figure 1
Definition of the quadrilateral area coordinatesLi of the QACM-I.Four dimensionless shape parametersg1, g2, g3, and g4 to each of the quadrangles, as shown in Figure 2, must be defined as(4)g1=A′A,g2=A′′A,g3=1-g1,g4=1-g2,0≤gi≤1,where A′ and A′′ are the areas of Δ124 and Δ123, respectively. Different values of these shape parameters mean different shapes of a quadrangle. And the area coordinates of four corner nodes can be obtained: (5)Node1g2,g4,0,0Node20,g3,g1,0Node30,0,g4,g2Node4g3,0,0,g1.Figure 2
Definitions ofg1, g2, g3, and g4.The relations between the area coordinatesLi and the isoparametric coordinates (ξ,η) are(6)L1=141-ξg21-η+g31+ηL2=141-ηg41-ξ+g31+ξL3=141+ξg11-η+g41+ηL4=141+ηg11-ξ+g21+ξ.


### 2.2. QACM-II [52]

As shown in Figure3, Mi (i=1,2,3,4) are the midside points of element sides 23¯, 34¯, 41¯, and 12¯, respectively. Thus, the position of an arbitrary point P within the quadrilateral element 1234¯ can be uniquely specified by the two-component area coordinates Z1 and Z2 (QACM-II), which are defined as(7)Z1=4Ω1A,Z2=4Ω2A,where Ω1 and Ω2 are the generalized areas of ΔPM2M4 and ΔPM3M1, respectively. It must be noted here that the values of generalized areas Ω1 and Ω2 can be both positive and negative: for ΔPM2M4 (or ΔPM3M1), if the permutation order of points P, M2, and M4 (or P, M3, and M1) is anticlockwise, a positive Ω1 (or Ω2) should be taken; otherwise, Ω1 (or Ω2) should be negative.Figure 3
Definition of the quadrilateral area coordinatesZi of QACM-II.Thus, the local coordinates of the corner nodes and midside points can be obtained:(8)Node1-g1-g2,-g1-g4Node2g1+g2,-g2-g3Node3g3+g4,g2+g3Node4-g3-g4,g1+g4M11,0M20,1M3-1,0M40,-1.The relations between the QACM-II and the QACM-I are(9)Z1=1Aa3-a1+b3-b1x+c3-c1y+g-1Z2=1Aa4-a2+b4-b2x+c4-c2y+g-2withg-1=g2-g1,g-2=g3-g2.And Z1 and Z2 can also be expressed in terms of ξ and η as follows:(10)Z1=ξ+g-2ξηZ2=η+g-1ξη.It can be seen that the new area coordinates Z1 and Z2 will degenerate to the isoparametric coordinates ξ and η for rectangular element cases.


### 2.3. QACM-III [54]

As shown in Figure4, 13¯ and 24¯ are the two diagonals of the quadrilateral 1234¯. Then, the position of an arbitrary point P within or outside the quadrilateral 1234¯ can be uniquely specified by the two-component area coordinates T1 and T2 (QACM-III), which are defined as(11)T1=S1A,T2=S2A,where S1 and S2 are the generalized areas of ΔP42 and ΔP13, respectively. The values of generalized areas S1 and S2 can be both positive and negative: for ΔP42 (or ΔP13), if the permutation order of the points P, 4, and 2 (or P, 1, and 3) is anticlockwise, a positive S1 (or S2) should be taken; otherwise, S1 (or S2) should be negative.Figure 4
Definition of the quadrilateral area coordinatesTi of the QACM-III.Then, the local coordinates of the corner nodes can be written as(12)node1-g1,0node20,-g2node3g3,0node40,g4.The relations between the QACM-III and the QACM-I are(13)T1=g3-L1-L2=L3+L4-g1T2=g4-L2-L3=L4+L1-g2.And T1 and T2 can also be expressed in terms of ξ and η as follows:(14)T1=14ξ+η+g3-g11+ξηT2=14-ξ+η+g4-g21-ξη.


## 2.1. QACM-I [28, 29]

As shown in Figure1, the position of an arbitrary point P within a quadrilateral element 1234¯ is specified by the area coordinates L1, L2, L3, and L4, which are defined as(1)Li=AiA,i=1,2,3,4,where A is the area of the quadrilateral element; Ai (i=1,2,3,4) are the areas of the four triangles constructed by point P and four element sides 23¯, 34¯, 41¯, and 12¯, respectively. L1, L2, L3, and L4 can be expressed in terms of Cartesian coordinates (x,y) as follows:(2)Li=12Aai+bix+ciy,i=1,2,3,4with(3)ai=xjyk-xkyj,bi=yj-yk,ci=xk-xi,i=1,2,3,4;j=2,3,4,1;k=3,4,1,2.Figure 1
Definition of the quadrilateral area coordinatesLi of the QACM-I.Four dimensionless shape parametersg1, g2, g3, and g4 to each of the quadrangles, as shown in Figure 2, must be defined as(4)g1=A′A,g2=A′′A,g3=1-g1,g4=1-g2,0≤gi≤1,where A′ and A′′ are the areas of Δ124 and Δ123, respectively. Different values of these shape parameters mean different shapes of a quadrangle. And the area coordinates of four corner nodes can be obtained: (5)Node1g2,g4,0,0Node20,g3,g1,0Node30,0,g4,g2Node4g3,0,0,g1.Figure 2
Definitions ofg1, g2, g3, and g4.The relations between the area coordinatesLi and the isoparametric coordinates (ξ,η) are(6)L1=141-ξg21-η+g31+ηL2=141-ηg41-ξ+g31+ξL3=141+ξg11-η+g41+ηL4=141+ηg11-ξ+g21+ξ.


## 2.2. QACM-II [52]

As shown in Figure3, Mi (i=1,2,3,4) are the midside points of element sides 23¯, 34¯, 41¯, and 12¯, respectively. Thus, the position of an arbitrary point P within the quadrilateral element 1234¯ can be uniquely specified by the two-component area coordinates Z1 and Z2 (QACM-II), which are defined as(7)Z1=4Ω1A,Z2=4Ω2A,where Ω1 and Ω2 are the generalized areas of ΔPM2M4 and ΔPM3M1, respectively. It must be noted here that the values of generalized areas Ω1 and Ω2 can be both positive and negative: for ΔPM2M4 (or ΔPM3M1), if the permutation order of points P, M2, and M4 (or P, M3, and M1) is anticlockwise, a positive Ω1 (or Ω2) should be taken; otherwise, Ω1 (or Ω2) should be negative.Figure 3
Definition of the quadrilateral area coordinatesZi of QACM-II.Thus, the local coordinates of the corner nodes and midside points can be obtained:(8)Node1-g1-g2,-g1-g4Node2g1+g2,-g2-g3Node3g3+g4,g2+g3Node4-g3-g4,g1+g4M11,0M20,1M3-1,0M40,-1.The relations between the QACM-II and the QACM-I are(9)Z1=1Aa3-a1+b3-b1x+c3-c1y+g-1Z2=1Aa4-a2+b4-b2x+c4-c2y+g-2withg-1=g2-g1,g-2=g3-g2.And Z1 and Z2 can also be expressed in terms of ξ and η as follows:(10)Z1=ξ+g-2ξηZ2=η+g-1ξη.It can be seen that the new area coordinates Z1 and Z2 will degenerate to the isoparametric coordinates ξ and η for rectangular element cases.


## 2.3. QACM-III [54]

As shown in Figure4, 13¯ and 24¯ are the two diagonals of the quadrilateral 1234¯. Then, the position of an arbitrary point P within or outside the quadrilateral 1234¯ can be uniquely specified by the two-component area coordinates T1 and T2 (QACM-III), which are defined as(11)T1=S1A,T2=S2A,where S1 and S2 are the generalized areas of ΔP42 and ΔP13, respectively. The values of generalized areas S1 and S2 can be both positive and negative: for ΔP42 (or ΔP13), if the permutation order of the points P, 4, and 2 (or P, 1, and 3) is anticlockwise, a positive S1 (or S2) should be taken; otherwise, S1 (or S2) should be negative.Figure 4
Definition of the quadrilateral area coordinatesTi of the QACM-III.Then, the local coordinates of the corner nodes can be written as(12)node1-g1,0node20,-g2node3g3,0node40,g4.The relations between the QACM-III and the QACM-I are(13)T1=g3-L1-L2=L3+L4-g1T2=g4-L2-L3=L4+L1-g2.And T1 and T2 can also be expressed in terms of ξ and η as follows:(14)T1=14ξ+η+g3-g11+ξηT2=14-ξ+η+g4-g21-ξη.


## 3. Definition of the Drilling DOFs

As shown in Figure5, Long et al. [23–25] defined the drilling DOFs θi as the additional rigid rotations at the element nodes.Figure 5
DOFs of a membrane element.
(a)
Translational freedom
(b)
Drilling freedomThe characteristics of this definition are as follows.(1)
The change of the angle between two adjacent sides along with the element deformation is allowed.(2)
The rotationθ of the element side has definite relation with the nodal drilling freedom θi.In this definition, the displacement fields within the domain of an element are assumed to include two parts:(15)u=u0+uθ=u0v0+uθvθ,where u0 are the displacement fields determined by the nodal translational displacements and uθ are the additional displacement fields only determined by the vertex rigid rotations.The element nodal displacement vectorqe is defined by(16)qe=u1v1θ1u2v2θ2u3v3θ3u4v4θ4T.


## 4. Formulations of the New Elements QAC4θ and QAC4θM

According to the definition of drilling DOFs, the element boundary displacement can be assumed as(17)u-ij=u-ij0+u-θijv-ij=v-ij0+v-θij,ij=12,23,34,41,where the translational displacements u-ij and v-ij can be interpolated by the nodal displacements(18)u-120v-120=-T1g1u1v1-T2g2u2v2,u-230v-230=-T2g2u2v2+T1g3u3v3,u-340v-340=T1g3u3v3+T2g4u4v4,u-410v-410=T2g4u4v4-T1g1u1v1.The element boundary displacements caused by the additional vertex rigid rotations can be assumed by using the QACM-III:(19)u-23θv-23θ=b1c1T1T2T2θ2g22g3+T1θ3g2g32,u-34θv-34θ=b2c2T1T2T1θ3g32g4-T2θ4g3g42,u-41θv-41θ=b3c3T1T2-T2θ4g42g1-T1θ1g4g12,u-12θv-12θ=b4c4T1T2-T1θ1g12g2+T2θ2g1g22.It can be seen that, at the element corners (nodes), these boundary additional displacements are always equal to zero, and their normal derivatives of each edge are given by(20)∂∂nuijθvijθi=1dmbmcmθi,∂∂nuijθvijθj=1dmbmcmθj,ijm=124,231,342,413.The element displacement fields can be assumed in QACM-III as follows:(21)u=α1+α2T1+α3T2+α4T1T2+α5T12+α6T22,v=β1+β2T1+β3T2+β4T1T2+β5T12+β6T22.In order to determine the constant αi, six generalized conforming conditions are introduced:(22)∑i=14u-u-i=0,∑i=14ξiηiu-u-i=0,∫liju-u-ds-=0ij=23,34,41,12.The conforming conditions for v are similar to those for u.Then, the constantsαi and βi in (21) can be solved:(23)α=L-1Ruquβ=L-1Rvqv,where(24)qu=u1u2u3u4θ1θ2θ3θ4Tqv=v1v2v3v4θ1θ2θ3θ4T,L=4g3-g1g4-g20g12+g32g22+g420g3-g1g2-g40g12+g32-g22-g421g32-g22-g2g36g323g2231g32g42g3g46g323g4231-g12g42-g4g16g123g4231-g12-g22g1g26g123g223,Ru=111100001-11-100000121200b112-b112000121200b212-b212120012-b31200b312121200b412-b41200,Rv=111100001-11-100000121200c112-c112000121200c212-c212120012-c31200c312121200c412-c41200.Let(25)Hu=L-1Ru,Hv=L-1Rv.Thus, the displacement fields can be written in the following:(26)u=∑i=14Ni0ui+∑i=14Nuθiθi,v=∑i=14Ni0vi+∑i=14Nvθiθi,where the shape functions of translational displacements are(27)Ni0=Hu,1i+Hu,2iT1+Hu,3iT2+Hu,4iT1T2+Hu,5iT12+Hu,6iT22=Hv,1i+Hv,2iT1+Hv,3iT2+Hv,4iT1T2+Hv,5iT12+Hv,6iT22,i=1,2,3,4.And the shape functions of additional displacement fields related to the vertex rigid rotations are(28)Nuθi=Hu,1j+Hu,2jT1+Hu,3jT2+Hu,4jT1T2+Hu,5jT12+Hu,6jT22,Nvθi=Hv,1j+Hv,2jT1+Hv,3jT2+Hv,4jT1T2+Hv,5jT12+Hv,6jT22,i=1,2,3,4,j=i+4,j>4.The element strain fields are given by(29)ε=Bqqe,where(30)Bq=B1B2B3B4,Bi=∂Ni0∂x0∂Nuθi∂x0∂Ni0∂y∂Nvθi∂y∂Ni0∂y∂Ni0∂x∂Nuθi∂y+∂Nvθi∂x,i=1,2,3,4,∂Ni0∂x∂Ni0∂y=Hu,2i2Ay4-y2x2-x4+Hu,3i2Ay1-y3x3-x1+Hu,4i2Ay4-y2T2+y1-y3T1x2-x4T2+x3-x1T1+Hu,5iAy4-y2T1x2-x4T1+Hu,6iAy1-y3T2x3-x1T2,∂Nuθi0∂x∂Nuθi0∂y=Hu,2j2Ay4-y2x2-x4+Hu,3j2Ay1-y3x3-x1+Hu,4j2Ay4-y2T2+y1-y3T1x2-x4T2+x3-x1T1+Hu,5jAy4-y2T1x2-x4T1+Hu,6jAy1-y3T2x3-x1T2,and ∂Nvθi0/∂x and ∂Nvθi0/∂y can be obtained following a similar procedure.Finally, the element stiffness matrix is given by(31)kqq=∬ABqTDBqtdA,where D is the elastic matrix. This element is named QAC4θ.In order to make further improvement on element QAC4θ, a generalized bubble displacement field uλ is introduced with the following generalized conforming conditions: (32)∫lijuλds=∫lijuλvλds=0,ij=12,23,34,41.u
λ is assumed to be expressed in terms of the QACM-II:(33)uλ=λ1Z12+λ2Z22+λ3Z1+λ4Z2+λ5vλ=λ1′Z12+λ2′Z22+λ3′Z1+λ4′Z2+λ5′.Substitution of (33) into (32) yields(34)λ2=λ1λ3=2g1-g23λ1λ4=2g2-g33λ1λ5=2g1g3+g2g4-53λ1.And λi′ (i=1~5) have similar relations. Thus, the shape functions of this additional displacement field can be obtained:(35)Nλ1=Z12+Z22+2g1-g23Z1+2g2-g33Z2+2g1g3+g2g4-53.Then, the corresponding strain vector can be written as(36)ελ=ελxελy=Bλλ1λ1′=Bλλ,where(37)Bλ=∂Nλ1∂x00∂Nλ1∂y∂Nλ1∂y∂Nλ1∂x,∂Nλ1∂x=2b3-b1AZ1+2b4-b2AZ2+2b3-b1g1-g23A+2b4-b2g2-g33A∂Nλ1∂y=2c3-c1AZ1+2c4-c2AZ2+2c3-c1g1-g23A+2c4-c2g2-g33A.The final element stiffness matrix of the element is(38)ke=kqq-kλqTkλλ-1kλq,where(39)kqq=∬BqTDBqtdAkλλ=∬BλTDBλtdAkλq=∬BλTDBqtdA.This element is named QAC4θM.


## 5. Numerical Examples

Seven benchmark problems, which are listed in Table1, have been used for evaluating the performance of the elements. The results solved by the other 14 element models listed in Table 2 are also given for comparison.Table 1
List of benchmark problems.
Number
Benchmark problems (figure number)
Results
1
Patch test
(Figure6)
2
Cook’s skew beam
(Figure7)
Table3
3
Beam divided by five quadrilateral elements
(Figure8)
Table4
4
Beam divided by four quadrilateral elements
(Figure9)
Table5
5
MacNeal’s thin beam
(Figure10)
Table6
6
Thin curving beam
(Figure11)
Table7
7
Beam divided by two elements with distortion parameter
(Figure12)
Table8Table 2
List of element models for comparison.
Number
Element model
Reference
1
Q4
4-node isoparametric element
2
Q6
4-node isoparametric element with internal parameters
Wilson et al. [55]
3
QM6
4-node isoparametric element with internal parameters
Taylor et al. [56]
4
P-S
Hybrid stress element
Pian and Sumihara [57]
5
QUAD4
4-node element in MSC/NASTRAN
MacNeal and Harder [32]
6
Q4S
Membrane element with drilling DOFs
MacNeal and Harder   [8]
7
GQ12
Membrane element with drilling DOFs
Long and Xu [23]
8
GQ12M8
Membrane element with drilling DOFs
Long and Xu [23]
10
D-type
Membrane element with drilling DOFs
Ibrahimbegovic et al. [11]
11
Groenwold1995
Membrane element with drilling DOFs
Groenwold and Stander [58]
12
AQR8
Membrane element with drilling DOFs
Aminpour [59]
13
RGD20
Refined hybrid element
Chen and Cheung [60]
14
Q8
8-node isoparametric elementTable 3
Results of Cook’s beam.
Element
V
C
σ
A
max
⁡
σ
B
min
⁡
2 × 2
4 × 4
8 × 8
2 × 2
4 × 4
8 × 8
2 × 2
4 × 4
8 × 8
Q4
11.80
18.29
22.08
0.1217
0.1873
0.2242
−0.0960
0.1524
−0.1869
Q6
21.61
23.04
23.69
0.1930
0.2237
0.2345
−0.1783
0.1867
−0.1992
D-type
20.68
22.98
23.63
—
—
—
—
—
—
GQ12
20.89
23.06
23.67
0.1802
0.2209
0.2315
−0.1784
−0.1950
−0.2007
GQ12M8
22.49
23.44
23.78
0.2083
0.2338
0.2361
−0.2216
−0.2045
−0.2028
QAC4θ
21.00
23.05
23.66
0.1917
0.2241
0.2318
−0.1877
−0.1938
−0.2009
QAC4θM
22.25
23.42
23.78
0.2147
0.2358
0.2364
−0.2092
−0.2033
−0.2027
Reference values
23.96
0.2362
−0.2023Table 4
Results of cantilever beam with five elements.
Elements
Load I
Load II
v
A
σ
x
B
v
A
σ
x
B
Q4
45.7
−1761
50.7
−2448
Q6
98.4
−2428
100.4
−3354
GQ12
95.5
−2989
96.0
−4096
GQ12M8
100.0
−3000
101.0
−4147
QAC4θ
100.0
−3000
98.6
−3931
QAC4θM
100.0
−3000
101.0
−3977
Exact
100.0
−3000
102.6
−4050Table 5
Results of cantilever with four elements.
Elements
Tip deflections
Normalized values
Point A
Point B
Average
Point A
Point B
Average
D-type
—
—
0.3065
—
—
0.861
Q4S
—
—
0.2978
—
—
0.837
Groenwold1995
—
—
0.3086
—
—
0.867
GQ12
0.3337
0.3324
0.3331
0.938
0.934
0.936
GQ12M
0.3420
0.3404
0.3412
0.961
0.957
0.959
QAC4θ
0.3523
0.3516
0.3520
0.990
0.988
0.989
QAC4θM
0.3523
0.3516
0.3520
0.990
0.988
0.989
Reference value
0.3558
1.000Example 1 (patch test).
The constant strain/stress patch test using irregular mesh is shown in Figure6. Let Young’s modulus E=1000, Poisson’s ratio μ=0.25, and thickness of the patch t=1. Both QAC4θ and QAC4θM can present exact solutions.Figure 6
Patch test of constant stress/strain state.Example 2 (Cook’s skew beam).
This example, in which a skew cantilever with shear distributed load at the free edge, as shown in Figure7, was proposed by Cook et al. [61]. The results of vertical deflection at point C, the maximum principal stress at point A, and the minimum principal stress at point B are listed in Table 3.Figure 7
Cook’s skew beam problem.Example 3 (cantilever beam divided by five quadrilateral elements).
The cantilever beam, as shown in Figure8, is divided by five irregular quadrilateral elements. And two loading cases are considered: (a) pure bending under moment M and (b) linear bending under transverse force P. Young’s modulus E=1500, and Poisson’s ratio ν=0.25. The results of the vertical deflection vA at point A and the stress σxB at point B are given in Table 4.Figure 8
Cantilever beam with five irregular elements.Example 4 (cantilever beam divided by five quadrilateral elements).
As shown in Figure9, the cantilever beam is divided by four irregular quadrilateral elements. The results of the deflections at the tip points A and B are shown in Table 5.Figure 9
Cantilever beam modeled with four irregular elements.Example 5 (MacNeal’s beam).
Consider the thin beams presented in Figure10. Three different mesh shapes, rectangular, parallelogram, and trapezoidal, are adopted. This example, proposed by MacNeal and Harder [32], is a famous benchmark for testing the sensitivity to mesh distortion of the 4-node quadrilateral membrane elements.Figure 10
MacNeal’s beam.
(a)
(b)
(c)Table 6
The normalized results of the tip deflection for MacNeal’s beam.
Elements
LoadP
LoadM
Mesh (a)
Mesh (b)
Mesh (c)
Mesh (a)
Mesh (b)
Mesh (c)
Q4
0.093
0.035
0.003
0.093
0.031
0.022
Q6
0.993
0.677
0.106
1.000
0.759
0.093
QM6
0.993
0.623
0.044
1.000
0.722
0.037
QUAD4
0.904
0.080
0.071
—
—
—
P-S
0.993
0.798
0.221
1.000
0.852
0.167
RGD20
0.981
0.625
0.047
—
—
—
AQR8
0.993
0.986
0.977
—
—
—
Q4S
0.993
0.986
0.988
—
—
—
QAC4θ
0.904
0.867
0.906
0.910
0.8804
0.930
QAC4θM
0.993
0.984
0.988
1.000
0.992
0.998
Exact
1.000 (−0.1081)
1.000 (−0.0054)Table 7
The tip deflection of a thin curving beam.
h/R
Q4
QM6
QUAD4
GQ12
GQ12M
QAC4θ
QAC4θM
Exact
0.03
0.024
0.339
0.615
0.670
0.897
0.712
1.000
1.000
0.006
0.001
0.022
0.163
0.612
0.896
0.645
1.008
1.000Table 8
Results of the tip deflection of a cantilever beam with distorted parametere.
e
0
0.5
1
2
3
4
4.9
Q4
28.0
21.0
14.1
9.7
8.3
7.2
6.2
Q8
100
100
99.3
89.3
59.7
31.6
19.0
QM6
100
80.9
62.7
54.4
53.6
51.2
46.8
P-S
100
81.0
62.9
55.0
54.7
53.1
49.8
GQ12
100
97.9
86.3
48.7
24.9
13.3
8.0
GQ12M
100
98.7
93.9
74.1
51.0
33.4
23.3
QAC4θ
100
99.9
98.9
99.8
102.0
102.2
100.3
QAC4θM
100
100
100
100
100
100
100
Exact
100There are two loading cases under consideration: pure bending and transverse linear bending. Young’s modulus of the beamE=107; Poisson’s ratio μ=0.3; the thickness of the beam t=0.1. The results of the tip deflection are shown in Table 6.Example 6 (thin curving beam).
As shown in Figure11, a cantilever thin curving beam is subjected to a transverse force at the tip. And it is also divided by five elements. Two thickness-radius ratios, (i) h/R=0.03 and (ii) h/R=0.006, are considered. The results of the tip displacement are listed in Table 7.Figure 11
Bending of a thin curved beam.Example 7 (cantilever beam divided by two elements containing a parameter of distortion).
The cantilever beam shown in Figure12 is divided by two elements. The shape of the two elements varies with the distorted parameter e. When e=0, both elements are rectangular. But with the increase of e, the mesh will be distorted more and more seriously. This is another famous benchmark for testing the sensitivity to the mesh distortion. For pure bending problem, the results of the tip deflection at point A are listed in Table 8.Figure 12
Cantilever beam divided by two elements.


## 6. Conclusions

In this paper, two membrane elements with drilling DOFs, named QAC4θ and QAC4θM, are developed by using the quadrilateral area coordinate methods QACM-II and QACM-III. In their formulations, the additional rigid rotations at the element nodes are considered as the drilling DOFs, so that these two elements can allow the change of the angle between two adjacent sides along with the element deformations. Furthermore, since the quadrilateral area coordinates can keep the order of the Cartesian coordinates unchangeable while the mesh is distorted, the new elements exhibit better performance than other similar models and insensitivity to mesh distortion. It is demonstrated again that the quadrilateral area coordinate methods are effective tools for developing high-performance quadrilateral finite element models.



---

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