A New Method Applied to the Quadrilateral Membrane Element with Vertex Rigid Rotational Freedom

In order to improve the performance of the membrane element with vertex rigid rotational freedom, a new method to establish the local Cartesian coordinate system and calculate the derivatives of the shape functions with respect to the local coordinates is introduced in this paper.Themembrane elements with vertex rigid rotational freedom such as GQ12 andGQ12M based on this new method can achieve higher precision results than traditional methods. The numerical results demonstrate that the elements GQ12 and GQ12M with this new method can provide better membrane elements for flat shell elements. Furthermore, this new method presented in this paper offers a new approach for other membrane elements used in flat shell element to improve the computing accuracy.


Introduction
The finite element method has been used for solving problems in different fields of engineering.For now, a large number of different finite elements have been developed and the shell element is one of these elements to solve the multiscale problems [1,2].The shell element can be categorized into two types: curved shell element and flat shell element [3].The curved shell element can be used to make a good description of the geometrical shape of the shell structure, so it is with a good calculation precision.But the formula expressions of the curved shell element are very complicated so it is limited in practical application.In this case, the flat shell element is becoming more attractive in deriving efficient numerical accuracy and its concise theory due to its superposition theory of the membrane element and plate bending element [4].
The membrane elements are among the simplest elements to develop, which are used for analyzing structures subjected to in-plane forces.The membrane elements are usually used to model the behavior of shear wall, stiffened sheet construction, and membrane action in shells.Some plane elements can be considered as membrane elements, such as the CST (constant strain triangle) element and the fournode isoparametric quadrilateral plane element (QUAD4) [5].In finite element methods, many plate bending elements also have been developed.Bazeley et al. [6] developed the confirming and nonconfirming plate bending elements in 1966.Batoz et al. [7] developed three types of plate bending elements, that is, the DKT element, the HSM element, and the SRI element for the analysis of plates and shells in 1980.In 1982, the quadrilateral plate bending element (DKQ) was formulated by Batoz and Tahar [8] based on the DKT element.Other plate bending elements have been developed in the following years [9][10][11][12].
The flat shell elements are developed by combining membrane elements containing two in-plane translational degrees of freedom and plate bending elements containing two rotational degrees of freedom and one out of plane translational degree of freedom.Since the in-plane rotational degrees of freedom are not included, the null values for the inplane rotational degrees of freedom will lead to singularity in structure stiffness matrix if all the elements are coplanar.The simplest method adopted to remove the singularity is to add the in-plane rotational degrees of freedom to the membrane elements, which can also improve their performance [5].After several triangular elements [13,14] with vertex rigid rotational freedom are proposed, the quadrilateral membrane elements including corner rotations have been developed by Allman [15] and MacNeal and Harder [16] as well.And then, numerous researches focused on membrane elements with drilling degrees of freedom have been accomplished, such as the models proposed by Iura and Atluri [17], Cazzani and Atluri [18], Piltner and Taylor [19], Geyer and Groenwold [20], Groenwold et al. [21], Choi et al. [22], Kugler et al. [23], and Cen et al. [24].In 1994, two new arbitrary quadrilateral membrane elements called GQ12 and GQ12M with vertex rotation were proposed by Long and Xu [25], resulting in more reasonable compatible conditions between adjoining elements and a more simple formulation.
In order to further improve the computing accuracy of the elements GQ12 and GQ12M, a new method for establishing the local Cartesian coordinate system and calculating the derivatives of the shape function with respect to the local coordinates is presented in this paper.This new method can be applied to the membrane elements which may provide a constituent part for flat shell elements.The elements GQ12 and GQ12M have more accurate numerical results compared to the bilinear quadrilateral element Q4 [3,25].In this paper, the applications of the new method to the elements GQ12 and GQ12M were proposed to examine the performance of this new method described in this paper.

The Membrane Elements GQ12 and GQ12M
A quadrilateral membrane element with rotational degree of freedom (shown in Figure 1) was proposed in [25].
The nodal displacement vector over this element is given by The freedoms at each node are in which   and V  are the translations and   is the additional rotation at each corner.A membrane element is usually located in a three-dimensional space, so its orientation may be in any directions.For convenience, the membrane element is commonly studied in the local Cartesian coordinate system represented by (  ,   ,   ) and located on the membrane element.Once a quantity is formed in the local system, it can be transformed to the global coordinate system.In the local system, the element GQ12 is formulated by the displacement field which includes two parts as follows [25]: in which u 0 is the conventional bilinear compatible displacement decided by the translations [  V  ] T at the nodes and u  is the additional displacement determined by the rigid rotation   ( = 1, 2, 3, 4) at each node.The element GQ12 displacement can be described in the form of where N is the shape function matrix of element GQ12, which can be given in terms of in which in which where   ,   and    ,    ( = 1, 2, 3, 4) denote the intrinsic coordinates and the local Cartesian coordinates of the element nodes, respectively.
In the local system, the element stiffness matrix of GQ12 can be written as [25] where  is the thickness of the element, |J| is the determinant of the Jacobian matrix, and B is the element strain matrix which can be expressed as In (9), D is the elasticity matrix and it can be expressed as [4] in which  is the elastic modulus and  is Poisson's ratio.In order to further improve the precision of the calculation, the element GQ12M based on the element GQ12 was presented in [25] by Long and Xu.The element GQ12M is formulated with an assumed displacement given by adding a bubble displacement û to GQ12, which is in which û is the bubble displacement.It can be assumed to be where It can be considered as the shape function at the center node C of the element, as shown in Figure 2, which will be eliminated in the future. 1 and  2 in (14) are arbitrary independent parameters.Substituting (4) and ( 14) into (13), the element GQ12 displacement can be described in the form of The corresponding strain fields can be expressed as in which B is the same as (10), and B is the strain matrix of the bubble displacement û, which can be expressed as According to (12) and (17), the strain energy of element GQ12M can be written as where BT DB |J|   , in which  is the thickness of the element, |J| is the determinant of the Jacobian matrix, and the matrices B, D, and B are given by ( 10), (12), and (18), respectively.From the stationary condition / = 0, the arbitrary parameters  can be expressed in terms of the q  as Substituting ( 21) into ( 16), the shape functions of the element GQ12M can be obtained as in which N, N , K   −1 , and K   are given by ( 5), (14), and (20), respectively.
According to the shape functions of GQ12M expressed in (22), the element stiffness matrix of element GQ12M can be written as

Calculating the Stiffness Matrix of the Membrane Element under the Global Coordinate System in the Traditional Manner
The element stiffness matrices of GQ12 and GQ12M under the local coordinate system are described in the above section, and these matrices can be numerically calculated by the Gaussian integral method.It can be seen from ( 9) that the GQ12 element stiffness matrix is in which   and   are the weight functions, and NG is the integration order.Generally, the element stiffness matrix of It can be seen from ( 11) and ( 24) that the shape functions with respect to the local coordinates   /  and   /  and the determinant of the Jacobian matrix |J| need to be firstly determined to calculate the GQ12 element stiffness matrices.The conventional method [3,4] to compute   /  and   /  is to calculate the product of the inverse Jacobi matrix and the derivatives matrix of the shape functions with respect to the intrinsic coordinates, which can be expressed as and then It is worth pointing out that   in ( 25) and ( 26) is a substitution, which may be  0  ,   ,  V , and N in ( 6), (7), and (15), respectively.
Once the element local stiffness matrix for membrane element is obtained, it is necessary to transform it from the local coordinate system to the global coordinate system.So the element local coordinate system needs to be established first.In the current application of the membrane elements, a traditional manner has been adopted to establish the local Cartesian coordinate system over an element [5].
Figure 3 shows the global and local coordinate axes for the quadrilateral membrane element.The midpoints of sides 1-2, 2-3, 3-4, and 4-1 are represented by , , , and , respectively, which can be determined by the shape functions of the four-node isoparametric element.The element local plane is defined by creating two vectors intersecting each other and passing through the midpoints of the sides 2-3 and 3-4 of the quadrilateral as shown in Figure 3.  Assuming that the local   axis of the quadrilateral is parallel to the vector passing through nodes  and , the vector passing through these nodes can be given by [5] where   ,   , and   and   ,   , and   represent the global coordinates of nodes  and , respectively.The direction cosine    for the local   direction can be obtained by normalizing the vector with respect to its length, which is 2 is the length of the vector.
A reference vector V  which can define the plane of the element with the vector V   is obtained by creating a vector passing through nodes  and  of the element as shown in Figure 3 V where   ,   , and   and   ,   , and   represent the global coordinates of nodes  and , respectively.The normal vector to the plane can be obtained by the vector cross product of V   and V  , which is The direction cosine    for the local   direction can be obtained by normalizing vector V   with respect to its length.The local   axis is obtained by the vector cross product of the vectors V   and V   .The cross product of these two vectors will give the vector V   normal to the       plane, which is The direction cosine    for the local   direction is obtained by normalizing the vector V   with respect to its length, with reference to (28).So the 3 × 3 transformation matrix L for the transformation of coordinates from the local to the global axis can be expressed as in which    ,    , and    are the direction cosines for local   ,   , and   directions.The 12 × 12 transformation matrix R for the 4-node quadrilateral element stiffness from the local to the global coordinate systems can be represented by L, which is To calculate the structure stiffness matrix, the element stiffness matrix must be transformed to the global coordinate system.The transformation of the GQ12 element stiffness matrix from the local to the global coordinate system is given by [3,4] Similarly, from ( 20) and ( 23), the GQ12M element stiffness matrix under the global coordinate system can be described as in which It can be seen from the above equations that the key task of computing the membrane stiffness matrix is to accurately evaluate the derivatives of the shape functions with respect to the local coordinates   /  and   /  , the determinant of the Jacobian matrix |J|, and transformation matrix R. The membrane element in the three-dimensional space may be not regular and coplanar.Thus, the calculation accuracy of the traditional way to establish the element local Cartesian coordinate system is heavily dependent on the element regularity.
In order to improve the performance of the membrane elements GQ12 and GQ12M, a new method to establish the element local Cartesian coordinate system is proposed.And the corresponding techniques to calculate the derivatives of the shape functions with respect to the local coordinates, the transformation matrix from the local to the global coordinate systems, and the determinant of the Jacobian matrix are also provided in this paper [26,27].), and (, ) are used to denote the symbols of the global Cartesian coordinate system, the local Cartesian coordinate system, and the intrinsic coordinates system, respectively.

A New
In order to improve the calculation accuracy, a new way [26] to establish the local Cartesian coordinate system is suggested in this paper.As shown in Figure 4, the origin of the local Cartesian coordinate system is set to the origin of the intrinsic coordinate system.Firstly, the vectors r  and r  from the origin need to be built along the tangent directions of the axes  and , respectively, as illustrated in Figure 5.The direction of the vector r  is considered as the local   1 direction.Then the vectors r  and r  define the plane which is tangential to the element surface.The normal vector n to this plane can be obtained by the vector cross product of r  and r  , which is the local   3 direction.At last, the local   2 axis is obtained by the vector cross product of the vector in the local   1 direction and vector in local   3 direction.Now the local Cartesian coordinate system in which the axes   1 and   2 are tangential to the surface and   3 is directed in the normal direction (shown in Figure 4) is established successfully.
However, when the curvature of element surface is large, the numerical results are still not very accurate in this local Cartesian coordinate system.Therefore, in order to further enhance the precision of the calculation, the origin of the local Cartesian coordinate system can be set at the Gauss points.This idea came from the traction recovery method [28][29][30][31] in Boundary Element Method (BEM) which establishes the local Cartesian coordinate systems at each node of the element in order to improve the calculation accuracy.In finite element method (FEM), the numerical integration is performed on Gauss points, so the local Cartesian coordinate systems can be established at Gauss points.Take the 2 × 2 Gauss point integration as an example, as indicated in Figure 6; the local Cartesian coordinate systems with the  origins at the four Gauss points can be established in a similar way as shown in Figure 4.
It is reasonable to establish the local Cartesian coordinate systems at Gauss points since the numerical integration is performed on them.In this case, the calculation accuracy can be improved, because the local Cartesian coordinate systems which are established on the tangent plane to the element surface at Gauss points can adapt to the curved element surface better.In order to better describe this problem, an example shown in Figure 7 which is a curved element surface is introduced to illustrate the performance of this new local Cartesian coordinate systems established at Gauss points.Figure 8 describes the local Cartesian coordinate system established in the traditional method at the curved element surface.It can be observed that the element local plane has a large difference with the curved element surface.By comparison, the element local planes defined by the local Cartesian coordinate systems established at 2 × 2 Gauss points are more accordant with the curved element surface,  as illustrated in Figure 9.In essence, Gaussian integration is the summation of the numerical results at Gauss points, so it is reasonable to believe that the calculation accuracy of the numerical integration can be improved by establishing the local Cartesian coordinate system at each Gauss point.
As shown in Figure 9, the local Cartesian coordinate systems are different since they are established in different tangent planes to the curved element surface.So the direction cosines    ,    , and    of the local coordinate system established at different Gauss point with respect to the global coordinate system are different with each other.The detailed formulas and computation steps of    ,    , and    will be given in the following sections.According to (31) and (32), it can be seen that the transformation matrix R for the transformation of element stiffness matrix from the local to the global coordinate systems may be different.So the element stiffness matrices of GQ12 and GQ12M under the global coordinate system expressed as (33), (35), and ( 36) can be modified to be

Derivatives of Shape Functions with respect to Local
Coordinates.According to the derivation rule for compound function, the derivatives of the shape functions with respect to the local coordinates can be written as in which /   and /   ( = 1, 2) are the derivatives of the intrinsic coordinates with respect to the local coordinates, which can be evaluated using the method proposed by Lachat [26,[32][33][34]; that is, where  is the angle between   1 axis and  axis shown in Figure 4, and In (39),   is a substitution, which may be  0  ,   ,  V , and N in ( 6), (7), and (15), respectively.So separately substituting ( 6), (7), and ( 15) into (39) gives The derivatives of the shape functions with respect to the local coordinates can be obtained by substituting (40) into the above equations.Then the strain matrices B and B can be obtained by substituting the above equations into (10) and (18), respectively.

Determination of the Jacobian |J|.
For the convenience of description, the local Cartesian coordinate system involved in this section is based on the way as shown in Figure 4.For the local Cartesian coordinate systems with the origins at the Gauss points, the approach to determine |J| is similar.
Referring to Figure 10, the Jacobian is equal to the magnitude of the vector cross product of the vectors r  and r  , which are located in the local tangent plane to the surface as can also be seen in Figure 5.The Jacobian of the transformation from the global three-dimensional coordinate system to the twodimensional intrinsic coordinate system of the surface patch is introduced in the form of [26] where in which i, j, and k are the orthogonal unit basis vectors of the global coordinate axes and n * is normal to the surface.So the unit normal vector may be obtained from in which  1 ,  2 , and  3 are the magnitude of the components of the unit normal vector in the global coordinate axes.The Jacobian of the transformation from the local orthogonal coordinate system to the intrinsic coordinate system is the same as that from the global orthogonal to the intrinsic coordinate systems.So the value of |J| in (38) can be computed using (44).

The Matrix Algorithm for the Membrane Element Stiffness Transformation from the Local to the Global Coordinate Systems
As described in Figure 4, the direction of local   1 axis is the tangent direction of the axis  at the origin, which is the direction of the vector r  .So the direction cosine    1 for the local   1 direction is obtained by normalizing the vector r  with respect to its length.According to (41) and (45), the direction cosine    1 can be written as [26] in which   ( = 1, 2, 3) is the global coordinates of the element nodes and | 1 | is the mold of the vector r  expressed in (41).The local   3 direction is the direction of the normal vector n to the element surface which can be obtained from (47).So the direction cosine    3 can be written as in which  1 ,  2 , and  3 are given by (47).Because the local   2 axis can be obtained by the vector cross product of the vector in the local   1 direction and vector in local   3 direction, the direction cosine    2 can be written as So the 3 × 3 transformation matrix L for the transformation of coordinates from the local to the global axis can be obtained from (32), which can be expressed in terms of   (,  = 1, 2, 3) as [26]  (51) The 12 × 12 transformation matrix R for the 4-node quadrilateral element stiffness from the local to the global coordinate systems can be obtained from (33).It is worth mentioning that the matrix R is just used for the membrane element stiffness transformation from the local to the global coordinate systems.For the 4-node quadrilateral shell

Numerical Examples
In order to examine the performance of the new method for GQ12 and GQ12M described in this paper, four numerical examples are proposed in this section.And the computational results are compared with those calculated by the traditional method based on the elements GQ12 and GQ12M.

Example 1:
Cook's Problem.The problem defined in Figure 11 was proposed by Cook [36] as a test case of plane stress elements, which is a standard example for accuracy testing of plane stress problem.Young's modulus of the element is  = 1 Pa, and Poisson's ratio is V = 1/3.The structure is meshed with 2 × 2, 4 × 4, and 8 × 8 elements.Here the new methods based on the elements GQ12 and GQ12M proposed in this paper are utilized to solve this problem.The best known answers taken from [35] are used for comparison, since there is no analytical solution available for this problem.
The results of the vertical deflection at C, the maximum stress at A, and the minimum stress at B on different meshes of the structure are given in Tables 2 and 3.The results demonstrate that the new method for the elements GQ12 and GQ12M have the desirable numerical accuracy, both for the displacement and for the stress.

Example 2: MacNeal Beam.
In order to examine the membrane locking phenomena, a typical test example [37,38] defined in Figure 12 was proposed by MacNeal.The thickness of the membrane elements is  = 0.1 m, Young's modulus is  = 10 7 Pa, and Poisson's ratio is V = 0.3.The beam end displacements obtained with several elements under the unit shear load at the free end are compared with the theoretical values taken from [38].The results under different meshes in Table 4 demonstrate that GQ12 and GQ12M with the new method described in this paper can delete membrane locking and have the same accuracy as the traditional method.shows the results of displacement at corner A which are obtained with different elements under these two loading cases.It demonstrates that the elements GQ12 and GQ12M with present method in this paper pass the patch test for a general quadrilateral mesh and achieve more accurate results.

Example 4:
Linear Analysis of a Hemispherical Shell.The aim of this test example [40] is to examine the performance of GQ12 and GQ12M elements with the present method in this paper as a constituent part of the flat shell element to compute a curved shell structure.As known, the flat shell elements are formed by superimposing the stiffness of membrane and plate bending elements.In this section, the plate bending element is based on the Mindlin theory which includes transverse shear deformations, and the membrane element may be selected from GQ12 and GQ12M with present method in this paper, respectively.Figure 14 shows a hemispherical shell structure with an 18-degree hole at the top [40].The radius of the shell is 10 m, thickness is 0.04 m, Young's modulus is 6.825 × 10 7 Pa, and Poisson's ratio is 0.3.The top and bottom circumferential edges of the hemisphere are free and the shell is subjected to two radial unit point loads.Only a quarter of the hemispherical shell with the meshes shown in Figure 15 is separated out for research due to the geometric symmetry.The radial displacement at point A from different meshes is compared with the theoretical solution  A = 0.094 [38] in Table 6.It can be seen that the shell elements based on the membrane elements with present method in this paper can achieve more accurate results than those based on a traditional manner.

Concluding Discussion
In this paper a new method for establishing the local Cartesian coordinate system and calculating the derivatives of the shape function with respect to local coordinates is applied to the membrane elements GQ12 and GQ12M.When the membrane elements are introduced to be the component of the flat shell element, this new method can establish the local Cartesian coordinate system on each Gauss point over the tangent plane to the element surface, and the precision is improved because the membrane elements are able to adapt the curved surface better.The numerical results of the test

Figure 2 :
Figure 2: The adding bubble node C.

Figure 3 :
Figure 3: The local Cartesian coordinate system in a traditional way.

Figure 4 :
Figure 4: Local orthogonal set of axes over a surface element in the new way.

Figure 5 :
Figure 5: The vectors r  and r  in the local tangent plane to the surface.

Figure 6 :
Figure 6: Local Cartesian coordinate systems with the origins at the Gauss points.

1 Figure 8 :Figure 9 :
Figure 8: The local Cartesian coordinate system established in the traditional method and the element local plane at the curved element surface.

Figure 10 :
Figure 10: Normal vector to a surface element.

Figure 13 :
Figure 13: A rectangular plate under two loading cases.

Figure 14 :
Figure 14: A hemispherical shell structure with an 18-degree hole on the top.

Figure 15 :
Figure 15: A quarter of the hemispherical shell structure with the meshes.

Table 1 :
[3]ts and weight functions for 2 × 2 Gauss quadrature.Gauss quadrature which can achieve the exact solution of the integral equation.The roots and weight functions for 2 × 2 Gauss quadrature are given in Table 1[3].

Method to Calculate the Stiffness Matrix of the Membrane Element in the Global Coordinate System 4
.1.Establishing the Local Coordinate System.In an effort to avoid confusion and awkward phrasing, the signs ( 1 ,  2 ,  3 ), (  1 ,   2 ,   3 ( = 1, 2, 3) is the global coordinates of the element nodes, and ,  are the intrinsic coordinates.

Table 2 :
The vertical deflection results of Cook's problem.

Table 3 :
The stress results of Cook's problem.

Table 4 :
The results of the MacNeal beam under different meshes.
The uniform stretching case is the patch test problem under constant strain.Because of the symmetry of the model and loads, only a quarter of the plate with the irregular mesh shown in Figure13is considered.Young's modulus of the element is  = 1 Pa, and Poisson's ratio is V = 1/3.Table5

Table 5 :
The displacement at corner A under two loading cases.