Efficient Co-Rotational 3-Node Shell Element

Corresponding Author: Gil Rama Department of Structural Analysis, Berlin Institute of Technology, Berlin, Germany Email: gil.rama@tu-berlin.de Abstract: Geometric nonlinear FE-analysis of thin-walled shell structures, using a co-rotational formulation and the Cell Smoothed Discrete Shear Gap triangular shell element (CS-DSG3) is presented in this study. The CS-DSG3 element formulation uses the Mindlin-Reissner kinematic hypothesis to include transverse shear effects. In order to avoid locking effects Discrete Shear Gap (DSG) method is applied. In addition, cell based smoothing technique is adopted in order to improve accuracy and stability of the element. For the purpose of comparison, the Discrete Kirchhoff-Constant strain-Triangle (DKT-CST) is also implemented and studied in the linear static analysis. In the framework of the co-rotational FE-analysis rotations and displacements are adopted as finite, while strains are infinitesimal. Large rotation theory has been utilized to take into account the non-vectorial characteristic of rotations. Several static linear and nonlinear benchmark examples are presented and compared with commercial FE software Abaqus and analytical results. The presented approach, using CS-DSG3 element in co-rotational nonlinear analysis, illustrates very good results compared to reference solution and Abaqus results. The numerical effort can be reduced compared to Lagrange formulation with a similar accuracy for the studied cases. The formulation (including CS-DSG3 shell element) has been implemented into a test program.


Introduction
Thin-walled structures are widely used in engineering practice. This is the consequence of the optimization strategy to reduce the structural dead-load whereby the structural carrying capacity is kept at a very high level. The application of thin-walled structures becomes increasingly diverse ranging from tiny machine parts over aircrafts and ships (Rama, 2014) to bridges, buildings, storage vessels, etc. The continuous development of new complex structural designs (Samaei et al., 2011) demands reliable and efficient numerical tools for modeling and simulation of the elastic behavior of thin-walled structures.
The Finite Element Method (FEM) has established itself as the method of choice for problems in the field of structural analysis. Over the past 60 years, the FEM has experienced its evolution whereby its application in the field of structural analysis has been the major driving component of the development. Within the framework of the FEM, various shell elements represent the principle workhorse in modeling of thin-walled structures. The main requirements for shell elements are high efficiency, reliability and applicability over a wide range of thickness and curvature. Additionally, thinwalled structures are characterized by high susceptibility to geometrically nonlinear behavior, caused by large transverse deflections and therewith local rotations and the developed FEM formulations are also supposed to cover this aspect.
The requirement for high efficiency points to loworder elements. In its nature, this requirement is contrary to the requirements for high accuracy and reliability, so the short overview of available elements should focus only to high quality elements that offer a good compromise and therewith a conciliation of the requirements. The low-order shell elements that particularly stand out among numerous candidates are the "Discrete Kirchhoff Theory-Constant Strain Triangular element" (DKT-CST) and the Cell Smoothed Discrete Shear Gap triangular shell element (CS-DSG3) (Nguyen-Thoi et al., 2013). The DKT-CST shell element represents a superposition of the CST membrane element and the DKT plate element. The DKT approach, based on a single polynomial expression for the transverse deflection w, was first proposed by Dhatt (1969). The DKT plate has been extensively tested for its convergence properties and is still widely considered to be one of the best triangular elements for thin plates. The assumptions made in the DKT approach do not include transverse shear effects and therefore hold only for rather thin plates. In case of thicker plates, the appropriate firstorder approach would be the Reissner-Mindlin (RM) theory, in which the nodal deflection and rotations are independent from each other. If a rather thin structure is modeled by means of RM-based elements, the results are supposed to converge to those yielded by the Kirchhofbased approach. However, the MR-based elements often exhibit shear locking when slenderness ratio (thickness/with) is low. The CS-DSG3 is one of the latest Reissner-Mindlin triangular shell elements. It represents a further development of the DSG3 element originally developed by Bletzinger et al. (2000). The DSG3 eliminates parasitic locking effects by using the Discrete Shear Gap (DSG) method (Bletzinger et al., 2000).
Regarding the coverage of geometrical nonlinearities, the high efficiency request imposes the need to reconsider various formulations for description of nonlinear kinematics. Typical approach in commercially available software packages are the total and updated Lagrangian formulations, which offer rather high accuracy and reliability. Certain developments are essentially based on the updated Lagrangian formulation but use local reference frame attached to the structure for computation of the Cauchy stresses (Marinković et al., 2008). However, in a co-rotational approach (Rankin and Brogan, 1986;Izzuddin, 2005) the provision of such a coordinate system that is attached to the structure and performs the same rigid-body motion as the structure itself is essential. It provides means to separate the rigidbody motion from purely deformational motion. The corotational approach appears to be quite a promising solution to achieve high efficiency in cases involving large local rotations. The origins of the co-rotational description could be traced back to Wempner (1969;Belytschko and Hseih, 1973;Argyris et al., 1979). Recent developments and surveys of co-rotational shell elements are given in (Izzuddin, 2005) and (Wempner, 1969). Some of the advantages the CR-formulations may offer compared to the Lagrangian formulation (Izzuddin, 2005;Crisfield and Moita, 1996;Li and Vu-Quoc, 2007;Li et al., 2008; are: • A better convergence for large displacements and rotations, but small strain problems • The separation of geometrical and material nonlinearities, so that all formalisms for material nonlinearities used in combination with geometrically linear computations are readily applicable in CR-formulations as well • The numerical effort is less compared to Lagrangian formulations The CR formulation applied in this study is employing the element-base reference frame. This means that each single element is provided with an attached co-rotational frame and the resolution of accounting for rigid-body motion is element-wise (Marinković et al., 2012).

Element Formulations
The  (Fig. 1). The magnitude of the normal direction vector X h is t/2, with the shell thickness t: The vector U is given by the two parts U 0 and ξU h : where, u 0 , v 0 , w 0 are the displacements, r x , r y are the rotations of the shell mid-surface and ξ the transverse direction vector which is normal to the mid-surface. The linear strain field consisting of membrane strains ε m , flexural strains ξ κ b and shear strains γ can be expressed in the following form: , ,

CS-DSG3 Formulation
Firstly, an overview of the DSG3 formulation (Bletzinger et al., 2000) is given, followed by the further development (Nguyen-Thoi et al., 2013) which leads to the CS-DSG3 element formulation.

DSG Element
The shape functions related to the element coordinates x e and y e read (Bazeley et al. 1969): The approximation u e for the triangular element domain, associated with node i and shape function N i , can be expressed as follows: , , The strain displacement matrices B m and B b can be expressed as: where, N i,xe and N i,ye (i = 1,2,3) are the constant derivatives of the shape functions in the local element coordinate system (x e , y e , z e ; Fig. 2). The derivatives of the shape functions are obtained by the following geometric distances and the triangular element area A e : As previously mentioned, the DSG3 element eliminates shear locking effects by utilizing the 'shear gap' concept of displacements along the element edges. The shear strain field becomes to: The element stiffness matrix is then given in the global coordinate system by: where, T is the transformation matrix of coordinates from the global coordinate system (x,y,z) to the local element coordinate system (x e ,y e ,z e ), D m , D b , D s are the constitutive matrices, E is Young's modulus; t is the thickness of the shell and ν is the Poisson constant.
As previous mentioned has the element three rotational degrees of freedom at each node. The rotations r xi and r yi are directly associated with bending and twisting mode, in contrast the drilling degree of freedom r zi does not provide any resisting force or stiffness by itself. The resistance to the r zi rotation at each node comes from the coupling of the rotations of surrounding non planar elements (Erhart and Borrvall, 2013). Difficulty occurs if all the elements meet at a node are co-planar (e.g., flat segments). Then is the rotational degree of freedom totally free to spin and corresponding stiffness has a zero stiffness. In order to avoid this rotational singularity in the stiffness matrix the zero values corresponding to drilling degree of freedom r zi are replaced by an approximated value of 10 -3 times of maximal diagonal value in the element stiffness matrix (Nguyen-Thoi et al., 2013): 3 2 1 0 It was recommended (Bischoff et al., 2003) to add a stabilization term to the original DSG3 element to improve the accuracy of approximate solutions and to stabilize shear force oscillations. The basic idea is to modify the internal work parts which include parasitic locking effects (Llyly et al., 1993): where, Π y is the work done by shear deformation, Π y stab is the stabilization term and τ the stabilization parameter. The stabilization parameter can be expressed as: where, h e is the longest element edge, α is a weight factor and t is the shell thickness. There are two requirements which have to be met for τ. The first one, if the element edges become zero, τ has to be 1. The second one, the work done by shear deformation is zero for an infinitesimal shell thickness corresponding to the Kirchhoff-Love kinematics. Bischoff et al. (2003) transferred this idea directly to the DSG elements which leads to the resulting constitutive matrix D s which can be expressed as follows: In which k (k = 5/6) is the shear correction factor and α is a positive constant.

CS-DSG3 Element
In the CS-DSG3 formulation the element domain is divided into three DSG3 sub-triangle elements by connecting the central point O of the triangle with the three element nodes (Fig. 3). In each sub-triangle the DSG formulation is used to obtain the strain-displacement matrices of the sub-tringles. Afterwards, in order to smooth the strains in the sub-triangles, strain smoothing technique on the whole triangular element is applied. As a result the CS-DSG3 overcomes the drawback of the DSG3 element which depends on the sequence of node numbering and improves the accuracy and stability of the DSG3 element (Nguyen- Thoi et al., 2013).
The displacement vector of the central point is assumed as average of the displacement vectors of the element nodes: The displacement vectors corresponding to the subtriangles (Fig. 2) read: As a result of substituting d 0 and rearranging the strain displacement matrix, matrices for the third subtriangle are obtained as follows: 3  3  3  3  3  3  1  3  1  1  2   1  1  1  3 3 3  3  3  3  3  3  1  3  1  1  2   1  1  1  3 3 3 The calculation of the other sub-triangle strain displacement matrices proceeds in a quite analogous manner. The resulting smoothed strain gradient matrices are given by:

DKT-CST Formulation
As already mentioned the DKT-CST shell element is obtained by combining a membrane element (CST) with and a plate bending element (DKT) with their nodal degrees of freedom d mi and d pi :

DKT Plate Element
For the C 1 compatible (conforming) plate element a displacement interpolation function w consist of nine terms are used: These equations can be rewritten in the matrix form as: where, c is the vector of the coefficients ([c 1 …c 9 ] T ) The inverse of C matrix multiplied with the nodal degrees of freedom of the plate element d p leads to the unknown constants c (c 1 …c 9 ): As a result the bending strain displacement matrix B p can be expressed as: The plate stiffness matrix K p will be obtained by following integral, which can be solved by numerical integration:

CST Membrane Element
The strain displacement matrix of the CST element and the resulting stiffness matrix are already defined in Equation 10 and 18. For the in-plane degrees of freedom d m the matrix B m can be written as: The analytical integration over the element domain leads to the CST stiffness matrix K m : Afterwards the nodal matrices K mi,j and K pi,j (stiffness matrices of freedom i and j) are combined to the resulting DKT-CST stiffness matrix K DKT-CST,ij at Node i: As a result element stiffness matrix K DKT-CST is given by: In order to avoid rotational singularity in the stiffness matrix the values corresponding to drilling degree of freedom are replaced by an approximate value of 10 −3 times of maximal diagonal value in the element stiffness matrix (see Sec. 2.2).

Simplified Co-Rotational Geometrically Nonlinear FEM Formulation
A geometrically nonlinear co-rotational FEM formulation implies that the geometrical nonlinearities in structural behavior are accounted for by means of an auxiliary, local reference frame that is attached to the material and performs the same rigid-body motion as the structural material. Such an approach decouples rigidbody motion from deformational motion, thus allowing usage of engineering strain and stress measures in the formulation and, furthermore, decoupling geometrical from material nonlinearities. Theoretically, a corotational reference frame may be assigned to each material particle. In FEM formulations, this is typically done for each element integration point as those points are used in the evaluation of element tangential stiffness matrix. However, in the simplified corotational formulation, the essence of which is presented in this section, the rigid-body rotation is accounted for on a somewhat reduced level. Actually, an average rigid-body rotation is determined for each finite element and further used in the computation. In this manner, the elastic behavior of each element remains linear with respect to the local element frame attached to the element and following its rigid-body motion. Given the rotational matrix, R e , which describes the rigid-body rotation of an element, the element stiffness matrix is updated as follows: Hence, the deformation of an element with respect to the local, co-rotational frame is described by a linear model. The rotational matrix, current and initial configuration, x e and x 0e , are used to determine the rotation-free displacements with respect to the initial configuration as: Finally, the internal forces with respect to the initial configuration are computed using the linear stiffness matrix and rotation-free displacements and are furthermore rotated to obtain the internal forces with respect to the current configuration: Once the global tangential stiffness matrix and the vector of internal forces are updated, static or transient geometrically nonlinear FEM computation may proceed. As solid elements have only translational degrees of freedom and nodal forces as loads, Equation 50-52 give the essence of the proposed co-rotational FEM formulation for this type of finite elements. For other types of finite elements that additionally have rotations as degrees of freedom (beams, shells), the procedure needs to be adequately extended. The rotations and translations do not share the same properties and the update of rotational degrees of freedom is more demanding. The incremental rotations computed in each time-step are used to update the shell normals at each node. This is done by computing the incremental rotational matrix of a shell normal as (Argyris, 1982): and: With r g,x , r g,y and r g,z denoting the 3 incremental global rotations at the node. The rotational matrix Q is used to update the orientation of the shell node normal. Also, the initial normal at the node is rotated through the element rigid-body rotation. The so obtained direction is compared to the actual shell normal to get the deformational nodal rotations, free of rigid-body rotation. This enables the computation of internal nodal moments, which are needed to proceed with the geometrically nonlinear computation.
Thin-walled structures are known for their susceptibility to large rotations, whereby the strains remain small. This is why they represent a perfect candidate to check the accuracy of the proposed formulation. The formulation has been implemented with the CS-DSG3 shell.

Numerical Examples
Various examples have been studied to assess the performance of the CS-DSG3 in linear and co-rotational nonlinear static analysis. For comparison purpose the DKT-CST element and results of Abaqus using the S3 element are employed in the linear static cases. For the nonlinear static cases the results obtained by using the co-rotational approach and the CS-DSG3 element are compared with Abaqus solution (using S3 elements).

Linear Static Cases
Three linear static benchmark cases are studied using CS-DSG3, DKT-CST and Abaqus S3 elements. Figure 4 illustrates one quadrant of a pinched hemispherical shell with forces at point A and point B as shown in Fig. 4. The hemisphere has a radius of 10 m, thickness of 0.04 m, Young's modulus is 68.25 MPa, Poisson's ratio is 0.3 and an 18° hole at the top. The quadrant of the hemisphere is modeled utilizing symmetric boundary conditions. Five meshes (8, 32, 128, 512 and 2048 elements) are used. The utilized element orientation is presented in Fig. 4.

Pinched Hemispherical Shell with Hole
The absolute values of the resulting displacements in z-direction at point A (|Uz,a|), using CS-DSG3, DKT-CST and Abaqus S3 elements are presented in Fig. 5 and Table 1. For a better comparison displacements are normalized by the reference solution (Simo et al., 1989).

Skew Plate
The case geometry is shown Fig. 6. The analysis is performed for two different values of the skew angle, α: 30°, 90° (Fig. 6). Three meshes (32, 128 and 512 elements with element topology presented in Fig. 6) are used for each skew angle with CS-DSG3, DKT-CST and the Abaqus S3 elements. The plate is 0.01 m thick. All sides are 1.0 m long. The slenderness is 1/100 so that the plate is thin thus corresponding to the Kirchhoff-Love kinematics. Young's modulus is 30 MPa and Poisson's ratio is 0.3. The plate is loaded by a uniform pressure of 1.0E-6 MPa applied over the entire surface. The edges of the plate are all simply supported (u = v = w = 0).
The resulting vertical displacements at plate center are presented in Fig. 7 and Table 2. For a better comparison displacements are normalized by the reference solution (Morley, 1963).
Comparing the examples of the pinched hemispherical shell and the skew plate and considering each of the element formulations separately, one may note that the nature of convergence remains similar. Of course, with different element formulations different rate of convergence is obvious as a consequence of the fact that they account for different kinematics (transverse shear included or excluded) or resolve the locking effects in different manner.

Nonlinear Static Cases
Two geometric nonlinear static benchmark cases are solved using CS-DSG3 with a co-rotational formulation and Abaqus. The default automatic load incrementation scheme is adopted in Abaqus.

Bending of a Cantilever Plate
A flat cantilever plate (Fig. 8) is loaded with a bending moment at its tip. The solution to this problem, a beam 'rolled up' into a circular arc of radius r i , is given by the classical flexural formula (Shi and Voyiadjis 1991 where, L is the length of the cantilever and I the second moment of area. The analysis is performed with a mesh of 176 elements and with the following geometrical and physical properties: Young's modulus E = 2E6 MPa, Poisson's ratio of 0.0, L = 8.8 m, b = 0.4 m and h = 0.01 m. The load-displacement (|u|, |w|) curve for the middle tip node is presented in Fig. 9. The solution in Abaqus aborted (at load factor of = 0.64) before the load reaches its maximum (Fig. 9).

The Slit Annular Plate Loaded with Line Force
The slit annular plate is presented in Fig. 10. The problem has been studied in (Buechter and Ramm, 1992;Wriggers and Gruttmann, 1993;Kim et al., 2003) among others. The line force P is applied at one end of the slit while the other end of the slit is fixed.

Conclusion
The focus of the paper was highly efficient geometrically nonlinear analysis of shell structures. In order to reach that goal, two aspects were considered-the high quality low order elements and the co-rotational FEM formulation. The choice of the elements was explained and the basics of their formulation given. The element-base co-rotational FEM is chosen not only for the sake of numerical efficiency, but also due to high stability and therewith robustness it offers. With these properties, it may easily by the right choice for applications in the field of Multi-Body System (MBS) dynamics with elastic bodies included as well as in the field of real-time simulations, particularly when used in combination with the features of modern hardware (Nutti and Marinković, 2014). Popular linear and geometrically nonlinear benchmark cases were considered in order to demonstrate the applicability of the proposed solution.