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

1China State Construction Technical Center, Beijing 101300, China 2Department of Engineering Mechanics, School of Aerospace Engineering, Tsinghua University, Beijing 100084, China 3High Performance Computing Center, School of Aerospace Engineering, Tsinghua University, Beijing 100084, China 4Key Laboratory of Applied Mechanics, School of Aerospace Engineering, Tsinghua University, Beijing 100084, China


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][24][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 formu-lations. 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 highorder terms of and , their complete order in Cartesian coordinates and 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 4node 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][37][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][49][50][51].
Since the QACM-I contains four area coordinate components ( 1 , 2 , 3 , and 4 ), 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 1 and 2 . 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 1 and 2 . 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][24][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. [28,29]. As shown in Figure 1, the position of an arbitrary point within a quadrilateral element 1234 is specified by the area coordinates 1 , 2 , 3 , and 4 , which are defined as

QACM-I
where is the area of the quadrilateral element; ( = 1, 2, 3, 4) are the areas of the four triangles constructed by point and four element sides 23, 34, 41, and 12, respectively. 1 , 2 , 3 , and 4 can be expressed in terms of Cartesian coordinates ( , ) as follows: (3) Four dimensionless shape parameters 1 , 2 , 3 , and 4 to each of the quadrangles, as shown in Figure 2, must be defined as where and 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: The relations between the area coordinates and the isoparametric coordinates ( , ) are 2.2. QACM-II [52]. As shown in Figure 3, ( = 1, 2, 3, 4) are the midside points of element sides 23, 34, 41, and 12, respectively. Thus, the position of an arbitrary point within the quadrilateral element 1234 can be uniquely specified by the two-component area coordinates 1 and 2 (QACM-II), which are defined as where Ω 1 and Ω 2 are the generalized areas of Δ 2 4 and Δ 3 1 , respectively. It must be noted here that the values of generalized areas Ω 1 and Ω 2 can be both positive and negative: for Δ 2 4 (or Δ 3 1 ), if the permutation order of points , 2 , and 4 (or , 3 , and 1 ) is anticlockwise, a positive Ω 1 (or Ω 2 ) should be taken; otherwise, Ω 1 (or Ω 2 ) should be negative.  Thus, the local coordinates of the corner nodes and midside points can be obtained: The relations between the QACM-II and the QACM-I are And 1 and 2 can also be expressed in terms of and as follows: It can be seen that the new area coordinates 1 and 2 will degenerate to the isoparametric coordinates and for rectangular element cases. 2.3. QACM-III [54]. As shown in Figure 4, 13 and 24 are the two diagonals of the quadrilateral 1234. Then, the position of an arbitrary point within or outside the quadrilateral 1234 can be uniquely specified by the two-component area coordinates 1 and 2 (QACM-III), which are defined as where 1 and 2 are the generalized areas of Δ 42 and Δ 13, respectively. The values of generalized areas 1 and 2 can be both positive and negative: for Δ 42 (or Δ 13), if the permutation order of the points , 4, and 2 (or , 1, and 3) is anticlockwise, a positive 1 (or 2 ) should be taken; otherwise, 1 (or 2 ) should be negative. Then, the local coordinates of the corner nodes can be written as The relations between the QACM-III and the QACM-I are And 1 and 2 can also be expressed in terms of and as follows:

Definition of the Drilling DOFs
As shown in Figure 5, Long et al. [23][24][25] defined the drilling DOFs as the additional rigid rotations at the element nodes.
The 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 .
In this definition, the displacement fields within the domain of an element are assumed to include two parts: where u 0 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 vector q is defined by

Formulations of the New Elements QAC4 and QAC4 M
According to the definition of drilling DOFs, the element boundary displacement can be assumed as where the translational displacements and V can be interpolated by the nodal displacements The element boundary displacements caused by the additional vertex rigid rotations can be assumed by using the QACM-III: 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 The element displacement fields can be assumed in QACM-III as follows: In order to determine the constant , six generalized conforming conditions are introduced: The conforming conditions for V are similar to those for .
Then, the constants and in (21) can be solved:     , Let Thus, the displacement fields can be written in the following: The element strain fields are given by where ] , ( = 1, 2, 3, 4) , and 0 V / and 0 V / can be obtained following a similar procedure.
Finally, the element stiffness matrix is given by 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: ( = 12, 23, 34, 41) . (32) u is assumed to be expressed in terms of the QACM-II: Substitution of (33) into (32) yields And ( = 1 ∼ 5) have similar relations. Thus, the shape functions of this additional displacement field can be obtained: Then, the corresponding strain vector can be written as where The final element stiffness matrix of the element is where This element is named QAC4 M.

Numerical Examples
Seven benchmark problems, which are listed in Table 1, 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.
Example 1 (patch test). The constant strain/stress patch test using irregular mesh is shown in Figure 6. Let Young's modulus = 1000, Poisson's ratio = 0.25, and thickness of the patch = 1. Both QAC4 and QAC4 M can present exact solutions.

Number
Benchmark problems (figure number) Results 1 Patch test ( Figure 6) 2 Cook's skew beam ( Figure 7) T a b l e 3 3 Beam divided by five quadrilateral elements ( Figure 8) T a b l e 4 4 Beam divided by four quadrilateral elements ( Figure 9) T a b l e 5 5 MacNeal's thin beam ( Figure 10) T a b l e 6 6 Thin curving beam ( Figure 11) T a b l e 7 7 Beam divided by two elements with distortion parameter ( Figure 12) T a b l e 8 Example 2 (Cook's skew beam). This example, in which a skew cantilever with shear distributed load at the free edge, as shown in Figure 7, 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.
Example 3 (cantilever beam divided by five quadrilateral elements). The cantilever beam, as shown in Figure 8, is  Table 4.
Example 4 (cantilever beam divided by five quadrilateral elements). As shown in Figure 9, the cantilever beam is    Table 5.
Example 5 (MacNeal's beam). Consider the thin beams presented in Figure 10. 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.
There are two loading cases under consideration: pure bending and transverse linear bending. Young's modulus of the beam = 10 7 ; Poisson's ratio = 0.3; the thickness of the beam = 0.1. The results of the tip deflection are shown in Table 6.
Example 6 (thin curving beam). As shown in Figure 11, 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) ℎ/ = 0.03 and (ii) ℎ/ = 0.006, are considered. The results of the tip displacement are listed in Table 7.
Example 7 (cantilever beam divided by two elements containing a parameter of distortion). The cantilever beam shown in Figure 12 is divided by two elements. The shape of the two elements varies with the distorted parameter . When = 0, both elements are rectangular. But with the increase of , 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.

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.