Abstract
Two 4-node generalized conforming quadrilateral membrane elements with drilling DOF, named QAC4 and QAC4M, 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.
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 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 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 (, , , and ), 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 and . 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 and . 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 Figure 1, the position of an arbitrary point within a quadrilateral element is specified by the area coordinates , , , and , which are defined aswhere is the area of the quadrilateral element; () are the areas of the four triangles constructed by point and four element sides , , , and , respectively. , , , and can be expressed in terms of Cartesian coordinates as follows:with
Four dimensionless shape parameters , , , and to each of the quadrangles, as shown in Figure 2, must be defined aswhere and are the areas of and , 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, () are the midside points of element sides , , , and , respectively. Thus, the position of an arbitrary point within the quadrilateral element can be uniquely specified by the two-component area coordinates and (QACM-II), which are defined aswhere and are the generalized areas of Δ and Δ, respectively. It must be noted here that the values of generalized areas and can be both positive and negative: for Δ (or Δ), if the permutation order of points , , and (or , , and ) is anticlockwise, a positive (or ) should be taken; otherwise, (or ) 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 areAnd and can also be expressed in terms of and as follows:It can be seen that the new area coordinates and will degenerate to the isoparametric coordinates and for rectangular element cases.
2.3. QACM-III [54]
As shown in Figure 4, and are the two diagonals of the quadrilateral . Then, the position of an arbitrary point within or outside the quadrilateral can be uniquely specified by the two-component area coordinates and (QACM-III), which are defined aswhere and are the generalized areas of Δ and Δ, respectively. The values of generalized areas and can be both positive and negative: for Δ (or Δ), if the permutation order of the points , 4, and 2 (or , 1, and 3) is anticlockwise, a positive (or ) should be taken; otherwise, (or ) 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 areAnd and can also be expressed in terms of and as follows:
3. Definition of the Drilling DOFs
As shown in Figure 5, Long et al. [23–25] defined the drilling DOFs as the additional rigid rotations at the element nodes.
(a) Translational freedom
(b) Drilling freedom
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 are the displacement fields determined by the nodal translational displacements and are the additional displacement fields only determined by the vertex rigid rotations.
The element nodal displacement vector is defined by
4. Formulations of the New Elements QAC4θ and QAC4θM
According to the definition of drilling DOFs, the element boundary displacement can be assumed aswhere the translational displacements and 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 are similar to those for .
Then, the constants and in (21) can be solved:where
Let Thus, the displacement fields can be written in the following:where the shape functions of translational displacements areAnd the shape functions of additional displacement fields related to the vertex rigid rotations are
The element strain fields are given bywhereand and can be obtained following a similar procedure.
Finally, the element stiffness matrix is given bywhere is the elastic matrix. This element is named QAC4θ.
In order to make further improvement on element QAC4θ, a generalized bubble displacement field is introduced with the following generalized conforming conditions:
is assumed to be expressed in terms of the QACM-II:Substitution of (33) into (32) yieldsAnd () have similar relations. Thus, the shape functions of this additional displacement field can be obtained:
Then, the corresponding strain vector can be written aswhere
The final element stiffness matrix of the element iswhere
This element is named QAC4θM.
5. 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 , Poisson’s ratio , and thickness of the patch . Both QAC4θ and QAC4θM can present exact solutions.
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 divided by five irregular quadrilateral elements. And two loading cases are considered: (a) pure bending under moment and (b) linear bending under transverse force . Young’s modulus , and Poisson’s ratio . The results of the vertical deflection at point A and the stress at point B are given in Table 4.
Example 4 (cantilever beam divided by five quadrilateral elements). As shown in Figure 9, 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.
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.
(a)
(b)
(c)
There are two loading cases under consideration: pure bending and transverse linear bending. Young’s modulus of the beam ; Poisson’s ratio ; the thickness of the beam . 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) and (ii) , 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 , 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.
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.
Conflict of Interests
The authors declare no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors would like to acknowledge the financial support of the National Natural Science Foundation of China (11272181), the Specialized Research Fund for the Doctoral Program of Higher Education of China (20120002110080), and the National Basic Research Program of China (Project no. 2010CB832701).