Abstract

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. The membrane elements with vertex rigid rotational freedom such as GQ12 and GQ12M 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.

1. 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 four-node 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–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 in-plane 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.

2. The Membrane Elements GQ12 and GQ12M

A quadrilateral membrane element with rotational degree of freedom (shown in Figure 1) was proposed in [25].

Figure 1 

The quadrilateral membrane element.

The nodal displacement vector over this element is given byThe freedoms at each node arein which and 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 is the conventional bilinear compatible displacement decided by the translations at the nodes and is the additional displacement determined by the rigid rotation at each node. The element GQ12 displacement can be described in the form ofwhere is the shape function matrix of element GQ12, which can be given in terms of in whichin whichwhere , and , 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, is the determinant of the Jacobian matrix, and is the element strain matrix which can be expressed as

In (9), 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 isin which is the bubble displacement. It can be assumed to bewhereIt 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. and in (14) are arbitrary independent parameters. Substituting (4) and (14) into (13), the element GQ12 displacement can be described in the form ofThe corresponding strain fields can be expressed asin which is the same as (10), and is the strain matrix of the bubble displacement , which can be expressed as

Figure 2 

The adding bubble node C.

According to (12) and (17), the strain energy of element GQ12M can be written aswherein which is the thickness of the element, is the determinant of the Jacobian matrix, and the matrices , , and are given by (10), (12), and (18), respectively.

From the stationary condition , the arbitrary parameters can be expressed in terms of the asSubstituting (21) into (16), the shape functions of the element GQ12M can be obtained as in which , , , and 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

3. 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 isin which and are the weight functions, and NG is the integration order. Generally, the element stiffness matrix of a 4-node quadrilateral element can be obtained using 2 × 2 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].

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 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 asand then

It is worth pointing out that in (25) and (26) is a substitution, which may be , , , and 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.

Figure 3 

The local Cartesian coordinate system in a traditional way.

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 isin which is the length of the vector.

A reference vector which can define the plane of the element with the vector is obtained by creating a vector passing through nodes and of the element as shown in Figure 3where , , 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 and , which isThe direction cosine for the local direction can be obtained by normalizing vector with respect to its length. The local axis is obtained by the vector cross product of the vectors and . The cross product of these two vectors will give the vector normal to the plane, which isThe direction cosine for the local direction is obtained by normalizing the vector with respect to its length, with reference to (28).

So the 3 × 3 transformation matrix for the transformation of coordinates from the local to the global axis can be expressed asin which , , and are the direction cosines for local , , and directions. The 12 × 12 transformation matrix for the 4-node quadrilateral element stiffness from the local to the global coordinate systems can be represented by , which isTo 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 asin 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 , and transformation matrix . 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].

4. A New 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 , , and are used to denote the symbols of the global Cartesian coordinate system, the local Cartesian coordinate system, and the intrinsic coordinates system, respectively.

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 and 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 is considered as the local direction. Then the vectors and define the plane which is tangential to the element surface. The normal vector to this plane can be obtained by the vector cross product of and , which is the local direction. At last, the local axis is obtained by the vector cross product of the vector in the local direction and vector in local direction. Now the local Cartesian coordinate system in which the axes and are tangential to the surface and is directed in the normal direction (shown in Figure 4) is established successfully.

Figure 4 

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

Figure 5 

The vectors and in the local tangent plane to the surface.

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–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.

Figure 6 

Local Cartesian coordinate systems with the origins at the Gauss points.

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.

Figure 7 

A curved element surface.

Figure 8 

The local Cartesian coordinate system established in the traditional method and the element local plane at the curved element surface.

Figure 9 

The local Cartesian coordinate systems established in the new method and the element local planes at the curved element surface.

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 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

4.2. 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 asin which and 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–34]; that is,where is the angle between axis and axis shown in Figure 4, andin which is the global coordinates of the element nodes, and , are the intrinsic coordinates.

In (39), is a substitution, which may be , , , and 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 and can be obtained by substituting the above equations into (10) and (18), respectively.

4.3. Determination of the Jacobian

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 is similar.

Referring to Figure 10, the Jacobian is equal to the magnitude of the vector cross product of the vectors and , 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 two-dimensional intrinsic coordinate system of the surface patch is introduced in the form of [26]wherein which , , and are the orthogonal unit basis vectors of the global coordinate axes and is normal to the surface. So the unit normal vector may be obtained from in which , , and are the magnitude of the components of the unit normal vector in the global coordinate axes.