# Several Treatments on Nonconforming Element Failed in the Strict Patch Test

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

---

## Abstract

For nonconforming finite elements, it has been proved that the models whose convergence is controlled only by the weak form of patch tests will exhibit much better performance in complicated stress states than those which can pass the strict patch tests. However, just because the former cannot provide the exact solutions for the patch tests of constant stress states with a very coarse mesh (strict patch test), their usability is doubted by many researchers. In this paper, the non-conforming plane 4-node membrane element AGQ6-I, which was formulated by the quadrilateral area coordinate method and cannot pass the strict patch tests, was modified by three different techniques, including the special numerical integration scheme, the constant stress multiplier method, and the orthogonal condition of energy. Three resulting new elements, denoted by AGQ6M-I, AGQ6M-II, and AGQ6M, can pass the strict patch test. And among them, element AGQ6M is the best one. The original model AGQ6-I and the new model AGQ6M can be treated as the replacements of the well-known models Q6 and QM6, respectively.


---

## Body

## 1. Introduction

In order to overcome the over-stiffness problem existing in conforming finite elements and improve their performances in regular and distorted meshes, the appearance of the nonconforming elements seems to be inevitable, such as the famous 4-node quadrilateral plane element Q6 proposed by Wilson et al. [1], the rectangular thin plate bending element ACM proposed by Melosh [2], and the triangular thin plate bending element BCIZ proposed by Bazeley et al. [3]. Since the conforming requirements along element edges are relaxed, the stiffness matrices of these elements become much softer. So, they can exhibit better performance than usual conforming element models in most cases. Furthermore, the construction procedures of these formulations are also easier. But, from the viewpoint of the variational principle, these elements violate the principle of minimum potential energy, which leads to that most of them cannot pass the strict patch tests, that is, their convergence cannot be guaranteed.The patch test proposed by Irons et al. [3, 4] has been taken as a standard for examining the convergence of the nonconforming elements for a long time. In such test, the nonconforming elements are required to produce exact solutions for a constant strain/stress patch with a very coarse mesh. Although it has been adopted very broadly, there are still some further discussions [5–8] on it. The patch test usually includes two types: the “strict” form and the “weak” form [9]. The strict patch test requests that the element can yield exact solutions under fixed coarse meshes and has been accepted by most researchers, while the weak form only requires the element to provide convergent solution for constant stress/strain problem as the mesh is repeatedly subdivided. Actually, the thought of the weak form of the patch test is more consistent with the concept of convergence.Many efforts have been made for developing 4-node nonconforming plane elements without any problem in convergence, such as the element QM6 proposed by Taylor et al. [10], QP6 by Wachspress [11], NQ6 by Pian and Wu [12], the generalized conforming element GC-Q6 by Long and Huang [13], the quasiconforming element QC6 by Chen and Tang [14], the hybrid-stress element P-S by Pian and Sumihara [15], and the Hamilton hybrid-stress element HH4 by Cen et al. [16]. All these elements can pass the strict form of patch test and possess much better performance than usual 4-node conforming isoparametric element Q4. In 2004, Chen et al. [17] formulated a 4-node nonconforming quadrilateral plane element AGQ6-I by using the quadrilateral area coordinates [18, 19]. This model exhibits excellent performance and is quite insensitive to mesh distortions in higher order problems. For example, it can keep the same high precision for famous MacNeal’s beam bending problem [20] using various distorted meshes, while other models cannot achieve this. This distinguished property can also be shown in nonlinear problems [21]. Furthermore, the explicit expressions of the element stiffness matrix can be obtained according to the integral formulae of the quadrilateral area coordinates [22], which can bring better computation efficiency. So, the element AGQ6-I attracts other researchers to make further developments and applications [23–25]. On the other hand, this element failed in strict patch test and can only pass the weak patch test. So, its convergence also arose some discussions [26–28].In this paper, three treatments were tried to improve the convergence of the element AGQ6-I, including the special numerical integration scheme, the constant stress multiplier method, and the orthogonal condition of energy. Three resulting new elements, denoted by AGQ6M-I, AGQ6M-II, and AGQ6M, respectively, can pass the strict patch test. Numerical examples show that, among the three new models, the element AGQ6M is the best one. The original model AGQ6-I and the new model AGQ6M can be treated as the replacements of the well-known models Q6 and QM6, respectively.


## 2. Brief Reviews on Element AGQ6-I

The 4-node quadrilateral plane membrane element AGQ6-I was formulated by the quadrilateral area coordinate methods [17]. Its displacement fields contain two parts: a generalized conforming low-order displacement part which is interpolated by the nodal displacements, and an additional second-order nonconforming part.


### 2.1. Low-Order Displacement Interpolation Formulae

For a 4-node quadrilateral plane element, the degrees of freedom can be expressed by the nodal displacement vector:(1)q=[u1v1u2v2u3v3u4v4]T.Then, the low-order displacement fields can be assumed as(2)u0=α1+α2(L3-L1)+α3(L4-L2)+α4(L3-L1)(L4-L2),v0=β1+β2(L3-L1)+β3(L4-L2)+β4(L3-L1)(L4-L2)
in which Li (i=1,2,3,4) are quadrilateral area coordinates and have been defined in [18]. In order to solve the undetermined parameters αi and βi (i=1,2,3,4), two nodal conforming conditions and two generalized conforming conditions along the perimeter ∂Ae of the element are introduced for both u0 and v0. Taking u0 as the example, these conditions are as follows:
(3)∑i=14(u0-u~)i=0,∑i=14(u0-u~)iξiηi=0,∮∂Ael(u0-u~)ds=0,∮∂Aem(u0-u~)ds=0.
Thus, the shape functions of the low-order displacement fields can be obtained:
(4)Ni0=-1-gi2+Li+Li+1+ξiηi(1-gi)P-,(i=1,2,3,4)
with
(5)P-=1g2g4+g1g3[(L3-L1)(L4-L2)-12(g2g4-g1g3)]
in which gi (i=1,2,3,4) are the shape parameters for a quadrangle and have been defined in [18].Then, the corresponding strain matrix can be expressed by(6)Bq=[Bq1Bq2Bq3Bq4],
where
(7)Bqi=[∂Ni∂x00∂Ni∂y∂Ni∂y∂Ni∂x],(i=1,2,3,4),(8a)∂Ni0∂x=bi2A+bj2A+ξiηigk2A(1+g1g3+g2g4)Sx,∂Ni0∂y=ci2A+cj2A+ξiηigk2A(1+g1g3+g2g4)Sy,(i=1,2,3,4;j=2,3,4,1;k=3,4,1,2),(8b)Sx=∑i′=14bi′ξi′ηi′[3(Lj′-Lm′)+(gj′-gk′)],Sy=∑i′=14ci′ξi′ηi′[3(Lj′-Lm′)+(gj′-gk′)],(i′=1,2,3,4⃡j′=2,3,4,1⃡k′=3,4,1,2⃡m′=4,1,2,3⃡),(9)bi=yj-yk,ci=xk-xj,(i,j,k⃡=1,2,3,4⃡)
in which (xi,yi) (i=1,2,3,4) are the Cartesian coordinates of the four corner nodes.


### 2.2. The Additional Second-Order Displacement Interpolation Formulae

The additional displacements are assumed to be(10)uλ=λ1L1L3+λ2L2L4,vλ=λ1′L1L3+λ2′L2L4.
Thus, the c matrix Bλ is given by
(11)Bλ=[b1L3+b3L12A0b2L4+b4L22A00c1L3+c3L12A0c2L4+c4L22Ac1L3+c3L12AcbL3+b3L12Ac2L4+c4L22Ab2L4+b4L22A].


### 2.3. The Element Stiffness Matrix

Finally, the element stiffness matrix of the element can be expressed by(12)Ke=Kqq-KλqTKλλ-1Kλq
with
(13)Kqq=∬ABqTDBqtdA,Kλλ=∬ABλTDBλtdA,Kλq=∬ABλTDBqtdA,
where t is the thickness of element; D is the elasticity matrix
(14)D=E1-μ2[1μ0μ10001-μ2],
where E and μ are Young’s modulus and Poisson’s ratio, respectively. For plane strain problems, the E and μ in (14) should be replaced by E/(1-μ2) and μ/(1-μ), respectively.The resulting element model is the element AGQ6-I. It can be easily shown that the values of the additional displacement fields given by (10) are zero at four element nodes. However, without any other restrictions along element edges, element AGQ6-I is a nonconforming element. On the other hand, since the relationship between area and Cartesian coordinates is linear, the displacement fields of AGQ6-I are actually a second-order complete polynomial in Cartesian coordinates. This makes the element AGQ6-I insensitive to mesh distortion in higher-order problems. So, it cannot only present the exact solution for pure bending problem but also be free of MacNeal’s trapezoid locking problem.Unfortunately, as a nonconforming element like Wilson’s Q6, element AGQ6-I cannot give the exact solutions for strict patch test either. Although the results can converge to the exact solutions by subdividing the mesh (weak patch test), further improvements on its convergence are valuable.


## 2.1. Low-Order Displacement Interpolation Formulae

For a 4-node quadrilateral plane element, the degrees of freedom can be expressed by the nodal displacement vector:(1)q=[u1v1u2v2u3v3u4v4]T.Then, the low-order displacement fields can be assumed as(2)u0=α1+α2(L3-L1)+α3(L4-L2)+α4(L3-L1)(L4-L2),v0=β1+β2(L3-L1)+β3(L4-L2)+β4(L3-L1)(L4-L2)
in which Li (i=1,2,3,4) are quadrilateral area coordinates and have been defined in [18]. In order to solve the undetermined parameters αi and βi (i=1,2,3,4), two nodal conforming conditions and two generalized conforming conditions along the perimeter ∂Ae of the element are introduced for both u0 and v0. Taking u0 as the example, these conditions are as follows:
(3)∑i=14(u0-u~)i=0,∑i=14(u0-u~)iξiηi=0,∮∂Ael(u0-u~)ds=0,∮∂Aem(u0-u~)ds=0.
Thus, the shape functions of the low-order displacement fields can be obtained:
(4)Ni0=-1-gi2+Li+Li+1+ξiηi(1-gi)P-,(i=1,2,3,4)
with
(5)P-=1g2g4+g1g3[(L3-L1)(L4-L2)-12(g2g4-g1g3)]
in which gi (i=1,2,3,4) are the shape parameters for a quadrangle and have been defined in [18].Then, the corresponding strain matrix can be expressed by(6)Bq=[Bq1Bq2Bq3Bq4],
where
(7)Bqi=[∂Ni∂x00∂Ni∂y∂Ni∂y∂Ni∂x],(i=1,2,3,4),(8a)∂Ni0∂x=bi2A+bj2A+ξiηigk2A(1+g1g3+g2g4)Sx,∂Ni0∂y=ci2A+cj2A+ξiηigk2A(1+g1g3+g2g4)Sy,(i=1,2,3,4;j=2,3,4,1;k=3,4,1,2),(8b)Sx=∑i′=14bi′ξi′ηi′[3(Lj′-Lm′)+(gj′-gk′)],Sy=∑i′=14ci′ξi′ηi′[3(Lj′-Lm′)+(gj′-gk′)],(i′=1,2,3,4⃡j′=2,3,4,1⃡k′=3,4,1,2⃡m′=4,1,2,3⃡),(9)bi=yj-yk,ci=xk-xj,(i,j,k⃡=1,2,3,4⃡)
in which (xi,yi) (i=1,2,3,4) are the Cartesian coordinates of the four corner nodes.


## 2.2. The Additional Second-Order Displacement Interpolation Formulae

The additional displacements are assumed to be(10)uλ=λ1L1L3+λ2L2L4,vλ=λ1′L1L3+λ2′L2L4.
Thus, the c matrix Bλ is given by
(11)Bλ=[b1L3+b3L12A0b2L4+b4L22A00c1L3+c3L12A0c2L4+c4L22Ac1L3+c3L12AcbL3+b3L12Ac2L4+c4L22Ab2L4+b4L22A].


## 2.3. The Element Stiffness Matrix

Finally, the element stiffness matrix of the element can be expressed by(12)Ke=Kqq-KλqTKλλ-1Kλq
with
(13)Kqq=∬ABqTDBqtdA,Kλλ=∬ABλTDBλtdA,Kλq=∬ABλTDBqtdA,
where t is the thickness of element; D is the elasticity matrix
(14)D=E1-μ2[1μ0μ10001-μ2],
where E and μ are Young’s modulus and Poisson’s ratio, respectively. For plane strain problems, the E and μ in (14) should be replaced by E/(1-μ2) and μ/(1-μ), respectively.The resulting element model is the element AGQ6-I. It can be easily shown that the values of the additional displacement fields given by (10) are zero at four element nodes. However, without any other restrictions along element edges, element AGQ6-I is a nonconforming element. On the other hand, since the relationship between area and Cartesian coordinates is linear, the displacement fields of AGQ6-I are actually a second-order complete polynomial in Cartesian coordinates. This makes the element AGQ6-I insensitive to mesh distortion in higher-order problems. So, it cannot only present the exact solution for pure bending problem but also be free of MacNeal’s trapezoid locking problem.Unfortunately, as a nonconforming element like Wilson’s Q6, element AGQ6-I cannot give the exact solutions for strict patch test either. Although the results can converge to the exact solutions by subdividing the mesh (weak patch test), further improvements on its convergence are valuable.


## 3. Treatments on the Convergence of Element AGQ6-I

### 3.1. Special Numerical Integration Scheme from QM6 [10]

In order to make the element Q6 [1] present exact solution for constant stress/strain problem with fixed coarse mesh, the following restriction condition was imposed by Taylor et al. [10] on its additional strain matrix:
(15)∫VeBλTdV≡∫-11∫-11BλT|J|tdξdη=0.
For satisfying the previous equation, the value of determinant |J| of Jacobi matrix must be a constant. So, when evaluating the differentials of x and y to ξ and η in matrix Kλq (see (13)) by numerical integration, the isoparametric coordinates of element centroid, ξ=η=0, were used to replace the coordinates of Gauss’s integral points. Then, the value of |J| would be always a constant, and the resulting element, denoted by QM6, can pass the strict patch test.In order to employ the previous technique to modify element AGQ6-I, the quadrilateral area coordinates can be rewritten in terms of the isoparametric coordinates:(16)L1=14(1-ξ)[g2(1-η)+g3(1+η)],L2=14(1-η)[g4(1-ξ)+g3(1+ξ)],L3=14(1+ξ)[g1(1-η)+g4(1+η)],L4=14(1+η)[g1(1-ξ)+g2(1+ξ)].
Then, after substituting (16) into (11), the additional strain matrix Bλ of element AGQ6-I in (11) can be rewritten in terms of the isoparametric coordinates.It is obvious that the additional strain matrixBλ of element AGQ6-I is different with that of element Q6. For element Q6, there are only terms of ξ, η, and ξη available in Bλ, and the integral interval of (15) is [-1,1]. So long as the |J| is a constant, (15) will be satisfied, and the modified version of Q6, element QM6, can pass the strict patch test. Compared with Q6, the additional strain matrix Bλ of AGQ6-I contains not only the terms of ξ, η, and ξη but also some constant terms which Q6 does not have. This leads to that the result of (15) is not a zero matrix. Therefore, the special numerical integration scheme mentioned previously cannot be directly used by element AGQ6-I. The influence of these constants must be considered.The constants in the additional strain matrixBλ are as follows:
(17)C1=b1(g1+g4)+b3(g2+g3),C2=b2(g1+g2)+b4(g3+g4),C3=c1(g1+g4)+c3(g2+g3),C4=c2(g1+g2)+c4(g3+g4).
In order to eliminate the influence of these constants, the modified additional strain matrix B-λ of element AGQ6-I can be rewritten as
(18)B-λ=Bλ-[C10C200C30C4C3C1C4C2].
Using this modified additional strain matrix and the numerical integral method of QM6, the resulting element will pass the strict patch test.The new model is named as AGQ6 M-I.


### 3.2. The Constant Stress Multiplier Method

The second approach is the constant stress multiplier method proposed by Pian and Wu [12], and can also make a nonconforming element pass the strict patch test. After introducing constant stress multiplier to eliminate the nonconforming energy on element sides, the final energy functional of the element can be written as
(19)∏mpe=∏pe-∮∂VeσCTnTuλds.
Finally, the new form of the element stiffness matrix can be obtained as
(20)KPe=[KqqKqλKqλTKλλ]-[0GGT0],
where
(21)G=∫VeBCTDBλdV
in which BC is the constant strain part of the matrix Bq given in (6).The final element stiffness matrix can be derived by usual static condense procedure. The resulting new element is denoted by AGQ6 M-II and can also pass the strict form of patch test.


### 3.3. The Orthogonal Condition of Energy

In order to ensure the element to pass the strict patch test, the strain energy under constant stress state caused by the nonconforming strains is required to be zero:(22)∫VeσcTελdV=0,
where ελ=nTuλ is the nonconforming strains of element. Equation (22) was called the orthogonal condition of energy in reference [29].By eliminating the constant factor in (22), the patch test of constant stress/strain for any element can be rewritten as
(23)∮∂VenTuλds=∫VeελdV=0.
Equation (23) can be further decoupled into three cases corresponding to the states subjected to one stress component of σx, σy, and τxy.Thus, four constantsCi (i=1,2,3,4) can be added to modify the additional strain field of AGQ6-I as follows:
(24)εx=(b12AL3+b32AL1+C1)λ1+(b22AL4+b42AL2+C2)λ2,εy=(c12AL3+c32AL1+C3)λ1′+(c22AL4+c42AL2+C4)λ2′,εxy=(c12AL3+c32AL1+C1)λ1+(c22AL4+c42AL2+C2)λ2+(b12AL3+b32AL1+C3)λ1′+(b22AL4+b42AL2+C4)λ2′.
From ∫AeεxdA=0, we obtain
(25)∫Ae[(b12AL3+b32AL1+C1)λ1+(b22AL4+b42AL2+C2)λ2]dA=0.
Its stronger forms can be written as
(26)∫Ae(b12AL3+b32AL1+C1)dA=0,∫Ae(b22AL4+b42AL2+C2)dA=0.
Then, C1 and C2 can be solved:
(27)C1=-16A[b1(1-g2g3)+b3(1-g4g1)],C2=-16A[b2(1-g3g4)+b4(1-g1g2)].
Similarly, from ∫AeεydA=0, C3 and C4 can be obtained:
(28)C3=-16A[c1(1-g2g3)+c3(1-g4g1)],C4=-16A[c2(1-g3g4)+c4(1-g1g2)].
And it can be proved that all previous Ci (i=1,2,3,4) satisfy ∫AeεxydA=0 automatically.Substitution of (27) and (28) into (24) yields the additional strain matrix of a new element by modifying element AGQ6-I. This new model is denoted by AGQ6 M.


## 3.1. Special Numerical Integration Scheme from QM6 [10]

In order to make the element Q6 [1] present exact solution for constant stress/strain problem with fixed coarse mesh, the following restriction condition was imposed by Taylor et al. [10] on its additional strain matrix:
(15)∫VeBλTdV≡∫-11∫-11BλT|J|tdξdη=0.
For satisfying the previous equation, the value of determinant |J| of Jacobi matrix must be a constant. So, when evaluating the differentials of x and y to ξ and η in matrix Kλq (see (13)) by numerical integration, the isoparametric coordinates of element centroid, ξ=η=0, were used to replace the coordinates of Gauss’s integral points. Then, the value of |J| would be always a constant, and the resulting element, denoted by QM6, can pass the strict patch test.In order to employ the previous technique to modify element AGQ6-I, the quadrilateral area coordinates can be rewritten in terms of the isoparametric coordinates:(16)L1=14(1-ξ)[g2(1-η)+g3(1+η)],L2=14(1-η)[g4(1-ξ)+g3(1+ξ)],L3=14(1+ξ)[g1(1-η)+g4(1+η)],L4=14(1+η)[g1(1-ξ)+g2(1+ξ)].
Then, after substituting (16) into (11), the additional strain matrix Bλ of element AGQ6-I in (11) can be rewritten in terms of the isoparametric coordinates.It is obvious that the additional strain matrixBλ of element AGQ6-I is different with that of element Q6. For element Q6, there are only terms of ξ, η, and ξη available in Bλ, and the integral interval of (15) is [-1,1]. So long as the |J| is a constant, (15) will be satisfied, and the modified version of Q6, element QM6, can pass the strict patch test. Compared with Q6, the additional strain matrix Bλ of AGQ6-I contains not only the terms of ξ, η, and ξη but also some constant terms which Q6 does not have. This leads to that the result of (15) is not a zero matrix. Therefore, the special numerical integration scheme mentioned previously cannot be directly used by element AGQ6-I. The influence of these constants must be considered.The constants in the additional strain matrixBλ are as follows:
(17)C1=b1(g1+g4)+b3(g2+g3),C2=b2(g1+g2)+b4(g3+g4),C3=c1(g1+g4)+c3(g2+g3),C4=c2(g1+g2)+c4(g3+g4).
In order to eliminate the influence of these constants, the modified additional strain matrix B-λ of element AGQ6-I can be rewritten as
(18)B-λ=Bλ-[C10C200C30C4C3C1C4C2].
Using this modified additional strain matrix and the numerical integral method of QM6, the resulting element will pass the strict patch test.The new model is named as AGQ6 M-I.


## 3.2. The Constant Stress Multiplier Method

The second approach is the constant stress multiplier method proposed by Pian and Wu [12], and can also make a nonconforming element pass the strict patch test. After introducing constant stress multiplier to eliminate the nonconforming energy on element sides, the final energy functional of the element can be written as
(19)∏mpe=∏pe-∮∂VeσCTnTuλds.
Finally, the new form of the element stiffness matrix can be obtained as
(20)KPe=[KqqKqλKqλTKλλ]-[0GGT0],
where
(21)G=∫VeBCTDBλdV
in which BC is the constant strain part of the matrix Bq given in (6).The final element stiffness matrix can be derived by usual static condense procedure. The resulting new element is denoted by AGQ6 M-II and can also pass the strict form of patch test.


## 3.3. The Orthogonal Condition of Energy

In order to ensure the element to pass the strict patch test, the strain energy under constant stress state caused by the nonconforming strains is required to be zero:(22)∫VeσcTελdV=0,
where ελ=nTuλ is the nonconforming strains of element. Equation (22) was called the orthogonal condition of energy in reference [29].By eliminating the constant factor in (22), the patch test of constant stress/strain for any element can be rewritten as
(23)∮∂VenTuλds=∫VeελdV=0.
Equation (23) can be further decoupled into three cases corresponding to the states subjected to one stress component of σx, σy, and τxy.Thus, four constantsCi (i=1,2,3,4) can be added to modify the additional strain field of AGQ6-I as follows:
(24)εx=(b12AL3+b32AL1+C1)λ1+(b22AL4+b42AL2+C2)λ2,εy=(c12AL3+c32AL1+C3)λ1′+(c22AL4+c42AL2+C4)λ2′,εxy=(c12AL3+c32AL1+C1)λ1+(c22AL4+c42AL2+C2)λ2+(b12AL3+b32AL1+C3)λ1′+(b22AL4+b42AL2+C4)λ2′.
From ∫AeεxdA=0, we obtain
(25)∫Ae[(b12AL3+b32AL1+C1)λ1+(b22AL4+b42AL2+C2)λ2]dA=0.
Its stronger forms can be written as
(26)∫Ae(b12AL3+b32AL1+C1)dA=0,∫Ae(b22AL4+b42AL2+C2)dA=0.
Then, C1 and C2 can be solved:
(27)C1=-16A[b1(1-g2g3)+b3(1-g4g1)],C2=-16A[b2(1-g3g4)+b4(1-g1g2)].
Similarly, from ∫AeεydA=0, C3 and C4 can be obtained:
(28)C3=-16A[c1(1-g2g3)+c3(1-g4g1)],C4=-16A[c2(1-g3g4)+c4(1-g1g2)].
And it can be proved that all previous Ci (i=1,2,3,4) satisfy ∫AeεxydA=0 automatically.Substitution of (27) and (28) into (24) yields the additional strain matrix of a new element by modifying element AGQ6-I. This new model is denoted by AGQ6 M.


## 4. Numerical Examples

### 4.1. Strict Patch Test

The constant strain/stress patch test using irregular mesh is shown in Figure1. Let Young’s modulus E=1000, Poisson’s ratio μ=0.25, and thickness of the patch t=1. Element AGQ6-I cannot give the exact solutions under such coarse meshes, while all three new elements with full integration scheme (2 × 2) can produce exact solutions without any problem. It shows that the proposed three methods are all workable for modifying nonconforming elements.Figure 1
Patch test.


### 4.2. Cook’s Skew Beam

This example was proposed by Cook et al. [9]. As shown in Figure 2, a skew cantilever beam is subjected to distributed shear load along its free edge. The results of vertical deflection at point C are listed in Table 1.Table 1
Displacement at point C of Cook’s beam.
Element
V
C
2 × 2
4 × 4
8 × 8
Q4
11.80
18.29
22.08
Q6
22.94
23.48
23.80
QM6
21.05
23.02
—
AGQ6-I
23.06
23.68
23.87
AGQ6M-I
20.86
23.00
23.69
AGQ6M-II
11.75
18.28
22.08
AGQ6M
20.74
22.99
23.69
Reference value
23.96Figure 2
Cook’s skew beam.Compared with the other elements, the accuracy of AQ6 M-II is even poor as Q4 element, while AGQ6 M-I and AGQ6 M are as accurate as QM6.


### 4.3. Cantilever Beam Divided by Five Quadrilateral Elements

The cantilever beam, as shown in Figure3, is divided by five irregular quadrilateral elements. Two loading cases are considered: (a) pure bending under moment M; (b) linear bending under transverse force P. The Young’s modulus E=1500, and Poisson’s ratio μ=0.25. Numerical results of the vertical deflection vA at point A and the stress σxB at point B are given in Table 2.Table 2
The deflections at point A of a cantilever beam.
Element type
Load (a)
Load (b)
v
A
σ
x
B
v
A
σ
x
B
Q4
45.7
−1761
50.7
−2448
Q6
98.4
−2428
100.4
−3354
QM6
96.1
−2497
98.0
−3235
QC6
96.1
−2439
98.1
−3339
NQ6
96.1
−2439
98.0
−3294
AGQ6-I
100.0
−3000
102.0
−4151
AGQ6M-I
91.9
—
94.5
—
AGQ6M-II
44.4
—
49.4
—
AGQ6M
96.0
−3015
97.9
−4135
Exact
100
−3000
102.6
−4050Figure 3
Cantilever beam with five irregular elements.From the results listed in Table2, AGQ6 M-II still exhibits poor performance as Q4. Although AGQ6 M-I has passed the patch test, it suffers from the mesh distortion now; only AGQ6 M can keep the same accuracy as other elements.


### 4.4. MacNeal’s Thin Cantilever Beam with Distorted Meshes

Consider the thin beams presented in Figure4. Three different mesh shapes, rectangular, parallelogram and trapezoidal, are adopted. This example, proposed by MacNeal and Harder [20], is a famous benchmark for testing the sensitivity to mesh distortion of quadrilateral membrane elements. In the three kinds of mesh, the trapezoidal shapes will cause the element locking easily.MacNeal’s thin beam.
(a)
(b)
(c)There are two loading cases under consideration: pure bending and transverse linear bending. The Young’s modulus of the beamE=107; the Poisson’s ratio μ=0.3; and the thickness of the beam t=0.1.The results of the tip deflection are shown in Table3. From the results listed in Table 3,Table 3
The normalized deflections at the free end of MacNeal’s beam.
Element type
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
AGQ6-I
0.993
0.994
0.994
1.000
1.000
1.000
AGQ6M-I
0.993
0.631
0.050
1.000
0.722
0.044
AGQ6M-II
0.093
0.034
0.027
0.093
0.031
0.022
AGQ6M
0.993
0.632
0.051
1.000
0.726
0.046
Exact
−0.1081
−0.0054(1) both AGQ6 M-II and Q4 suffer from locking problems for the three kinds of mesh distortions ((a) length-width ratio distortion, (b) parallelogram distortion, and (c) trapezoidal distortion) of three different meshes.(2) As a nonconforming element, Q6 still cannot overcome the locking problem for the trapezoidal mesh. Benefited from the area coordinate methods, AGQ6-I is locking free from each kind of mesh.(3) All those elements which can pass the strict patch test are locked in the trapezoidal mesh distortion. This verified the theory proposed in [30] again; that is, the trapezoidal locking is unavoidable for an element passed the strict patch test.


### 4.5. Mesh Distortion

The cantilever beam shown in Figure5 is divided by two elements. The shape of the two elements varies with the variety of 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. Pure bending M=2000 is considered. The results of the tip deflection at point A are listed in Table 4.Table 4
Results of the tip deflection of a cantilever beam subjected to a pure bendingM.
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
Q6
100
78.0
56.1
42.5
41.5
44.2
47.4
QM6
100
80.9
62.7
54.4
53.6
51.2
46.8
AGQ6-I
100
100
100
100
100
100
100
AGQ6M-I
100
80.8
55.4
24.1
12.9
9.4
7.6
AGQ6M-II
28.0
21.2
14.1
9.2
7.3
5.9
4.9
AGQ6M
100
83.8
66.5
60.1
61.4
60.3
56.0
Exact
100
100
100
100
100
100
100Figure 5
Cantilever beam divided by two elements with distorted parametere.The accuracy of element Q4 and AGQ6 M-II is the poorest. Their relative error reaches 72% whene=0. Except Q4, AGQ6 M-I is much more sensitive to mesh distortion than any other element although it can present exact solution while e=0. Compared with other elements which can pass the strict patch test, AGQ6 M is the most robust element against mesh distortion.


## 4.1. Strict Patch Test

The constant strain/stress patch test using irregular mesh is shown in Figure1. Let Young’s modulus E=1000, Poisson’s ratio μ=0.25, and thickness of the patch t=1. Element AGQ6-I cannot give the exact solutions under such coarse meshes, while all three new elements with full integration scheme (2 × 2) can produce exact solutions without any problem. It shows that the proposed three methods are all workable for modifying nonconforming elements.Figure 1
Patch test.


## 4.2. Cook’s Skew Beam

This example was proposed by Cook et al. [9]. As shown in Figure 2, a skew cantilever beam is subjected to distributed shear load along its free edge. The results of vertical deflection at point C are listed in Table 1.Table 1
Displacement at point C of Cook’s beam.
Element
V
C
2 × 2
4 × 4
8 × 8
Q4
11.80
18.29
22.08
Q6
22.94
23.48
23.80
QM6
21.05
23.02
—
AGQ6-I
23.06
23.68
23.87
AGQ6M-I
20.86
23.00
23.69
AGQ6M-II
11.75
18.28
22.08
AGQ6M
20.74
22.99
23.69
Reference value
23.96Figure 2
Cook’s skew beam.Compared with the other elements, the accuracy of AQ6 M-II is even poor as Q4 element, while AGQ6 M-I and AGQ6 M are as accurate as QM6.


## 4.3. Cantilever Beam Divided by Five Quadrilateral Elements

The cantilever beam, as shown in Figure3, is divided by five irregular quadrilateral elements. Two loading cases are considered: (a) pure bending under moment M; (b) linear bending under transverse force P. The Young’s modulus E=1500, and Poisson’s ratio μ=0.25. Numerical results of the vertical deflection vA at point A and the stress σxB at point B are given in Table 2.Table 2
The deflections at point A of a cantilever beam.
Element type
Load (a)
Load (b)
v
A
σ
x
B
v
A
σ
x
B
Q4
45.7
−1761
50.7
−2448
Q6
98.4
−2428
100.4
−3354
QM6
96.1
−2497
98.0
−3235
QC6
96.1
−2439
98.1
−3339
NQ6
96.1
−2439
98.0
−3294
AGQ6-I
100.0
−3000
102.0
−4151
AGQ6M-I
91.9
—
94.5
—
AGQ6M-II
44.4
—
49.4
—
AGQ6M
96.0
−3015
97.9
−4135
Exact
100
−3000
102.6
−4050Figure 3
Cantilever beam with five irregular elements.From the results listed in Table2, AGQ6 M-II still exhibits poor performance as Q4. Although AGQ6 M-I has passed the patch test, it suffers from the mesh distortion now; only AGQ6 M can keep the same accuracy as other elements.


## 4.4. MacNeal’s Thin Cantilever Beam with Distorted Meshes

Consider the thin beams presented in Figure4. Three different mesh shapes, rectangular, parallelogram and trapezoidal, are adopted. This example, proposed by MacNeal and Harder [20], is a famous benchmark for testing the sensitivity to mesh distortion of quadrilateral membrane elements. In the three kinds of mesh, the trapezoidal shapes will cause the element locking easily.MacNeal’s thin beam.
(a)
(b)
(c)There are two loading cases under consideration: pure bending and transverse linear bending. The Young’s modulus of the beamE=107; the Poisson’s ratio μ=0.3; and the thickness of the beam t=0.1.The results of the tip deflection are shown in Table3. From the results listed in Table 3,Table 3
The normalized deflections at the free end of MacNeal’s beam.
Element type
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
AGQ6-I
0.993
0.994
0.994
1.000
1.000
1.000
AGQ6M-I
0.993
0.631
0.050
1.000
0.722
0.044
AGQ6M-II
0.093
0.034
0.027
0.093
0.031
0.022
AGQ6M
0.993
0.632
0.051
1.000
0.726
0.046
Exact
−0.1081
−0.0054(1) both AGQ6 M-II and Q4 suffer from locking problems for the three kinds of mesh distortions ((a) length-width ratio distortion, (b) parallelogram distortion, and (c) trapezoidal distortion) of three different meshes.(2) As a nonconforming element, Q6 still cannot overcome the locking problem for the trapezoidal mesh. Benefited from the area coordinate methods, AGQ6-I is locking free from each kind of mesh.(3) All those elements which can pass the strict patch test are locked in the trapezoidal mesh distortion. This verified the theory proposed in [30] again; that is, the trapezoidal locking is unavoidable for an element passed the strict patch test.


## 4.5. Mesh Distortion

The cantilever beam shown in Figure5 is divided by two elements. The shape of the two elements varies with the variety of 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. Pure bending M=2000 is considered. The results of the tip deflection at point A are listed in Table 4.Table 4
Results of the tip deflection of a cantilever beam subjected to a pure bendingM.
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
Q6
100
78.0
56.1
42.5
41.5
44.2
47.4
QM6
100
80.9
62.7
54.4
53.6
51.2
46.8
AGQ6-I
100
100
100
100
100
100
100
AGQ6M-I
100
80.8
55.4
24.1
12.9
9.4
7.6
AGQ6M-II
28.0
21.2
14.1
9.2
7.3
5.9
4.9
AGQ6M
100
83.8
66.5
60.1
61.4
60.3
56.0
Exact
100
100
100
100
100
100
100Figure 5
Cantilever beam divided by two elements with distorted parametere.The accuracy of element Q4 and AGQ6 M-II is the poorest. Their relative error reaches 72% whene=0. Except Q4, AGQ6 M-I is much more sensitive to mesh distortion than any other element although it can present exact solution while e=0. Compared with other elements which can pass the strict patch test, AGQ6 M is the most robust element against mesh distortion.


## 5. Conclusions

In order to make the nonconforming element AGQ6-I pass the strict patch test, three treatments are used and tested. These methods are special numerical integration scheme used in QM6, constant stress multiplier method, and the orthogonal condition of energy. Numerical results of numerical examples show that each method can achieve the goal, but the performance for complicated stress problem of new elements will degenerate inevitably. Among the three treatments, the constant stress multiplier method made the element so rigid that the accuracy will deteriorate even in rectangular meshes. The special numerical integration scheme can make element QM6 exhibit good performance, but it is not always suitable for all other elements. The orthogonal condition of energy can effectively improve the compatibility of element AGQ6-I without losing too much accuracy.



---

*Source: 901495-2013-08-26.xml*