Plane Thermoelastic Waves in Infinite Half-Space Caused FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 14, N o 2, 2016, pp. 227 - 240 Original scientific paper 1AN EFFICIENT CO-ROTATIONAL FEM FORMULATION USING A PROJECTOR MATRIX UDC 531:519.6 Viet Anh Nguyen, Manfred Zehn, Dragan Marinković Structural Mechanics Department, Berlin Institute of Technology, Germany Abstract. Co-rotational finite element (FE) formulations can be seen as a very efficient approach to resolving geometrically nonlinear problems in the field of structural mechanics. A number of co-rotational FE formulations have been well documented for shell and beam structures in the available literature. The purpose of this paper is to present a co-rotational FEM formulation for fast and highly efficient computation of large three-dimensional elastic deformations. On the one hand, the approach aims at a simple way of separating the element rigid-body rotation and the elastic deformational part by means of the polar decomposition of deformation gradient. On the other hand, a consistent linearization is introduced to derive the internal force vector and the tangent stiffness matrix based on the total Lagrangian formulation. It results in a non- linear projector matrix. In this way, it ensures the force equilibrium of each element and enables a relatively straightforward upgrade of the finite elements for linear analysis to the finite elements for geometrically non-linear analysis. In this work, a simple 4-node tetrahedral element is used. To demonstrate the efficiency and accuracy of the proposed formulation, nonlinear results from ABAQUS are used as a reference. Key Words: Co-rotational FEM, Tetrahedral Element, Projector Matrix, Polar Decomposition, Newton-Raphson-Method 1. INTRODUCTION Beside the updated (UL) and total (TL) Lagrangian formulation, the co-rotational (CR) finite element (FE) formulations represent an efficient approach to handle large non-linear deformations of flexible structures. As it allows a relatively straightforward upgrade of the finite elements for linear analysis to the finite elements for geometrically Received February 15, 2016 / Accepted July 16, 2016 Corresponding author: Viet Anh Nguyen TUB, Structural Mechanics Department, Strasse des 17. Juni 135, 10623 Berlin E-mail: viet.a.nguyen@tu-berlin.de 228 V. A NGUYEN, M. ZEHN, D. MARINKOVIĆ non-linear analysis, the co-rotational approach has attracted significant interest over the past years [1]. The main idea of the approach is to decompose the total motion of each finite element into a rigid-body rotation and a moderate deformational part, which can be measured within a local element coordinate frame. The local frame can be defined at an arbitrary point of the element. It is attached to the element and performs the same rigid- body rotation and translation as the element itself. Material deformational response is determined with respect to this local frame, while geometric non-linearities are accounted for by means of the large rigid-body rotation. For sufficiently small local strains, the constant strain-displacement matrix can be applied [2]. In the past decade, many co-rotational frameworks for the geometrically non-linear analysis of flexible structures were developed. Even some classic, Lagrangian formulation based developments used certain aspects of the co-rotational formulation such as the element attached co-rotational frame to determine the true strains and stresses [3]. A rather thorough and systematic description of the small-strain co-rotational finite shell and beam elements is available in the literature [1, 4, 5]. In those developments, a local frame was defined at an arbitrary node of the element. The orientation change of the local frame was used to represent the element rigid-body rotation. On this basis, the rotation- free nodal displacements are determined by filtering out the rigid body rotation from the global displacements. Furthermore, the internal element forces were derived using a consistent linearization of the internal element energy with respect to the global element displacements. Marinkovic and Zehn [2] presented a simplified co-rotational FE formulation using the linear tetrahedral element for highly efficient computation of geometrically nonlinear deformations in the field of multibody systems (MBS) and real- time simulations. This development was based on a single co-rotational frame per finite element and a linear element stiffness matrix with respect to the co-rotational frame. In the work of Crisfield and Moita [6, 7], a co-rotational hexahedral element has been developed that enabled computation of deformations involving large strains. For this purpose, incompatible displacement modes were used to enhance the accuracy of computing the local element displacements. In addition, the element force vector and stiffness matrix were estimated using the assumption that the spin at the element centroid was zero. Espath et al. [8] used the co-rotational approach in combination with the NURBS-based isogeometric FEM to solve geometrically nonlinear dynamic problems. In this work, a co-rotational formulation for the non-linear static FE analysis is given using the tetrahedral element with the volume coordinates as linear shape functions. The rigid-body rotation is estimated by the polar decomposition of the deformation gradient, which is determined based on the initial and deformed element configurations. In addition, the internal element forces and element stiffness matrix are computed using directly the consistent variations of the internal element energy with respect to the global displacements. A resulting projector matrix is used to improve the results. The paper is organized as follows. After the introduction, in the second section, the co-rotational kinematics and the method to handle the rigid-body rotation for the tetrahedral element are presented. The derivation of the internal force vector and the tangent stiffness matrix are given in the third section. The fourth section gives the numerical procedure of the co-rotational FEM as a flow chart followed by a couple of numerical examples to verify the accuracy and to assess the method efficiency. Finally, based on the obtained results, conclusions are drawn. An Efficient Co-Rotational FEM formulation Using a Projector Matrix 229 2. NON-LINEAR CO-ROTATIONAL ELEMENT KINEMATICS In this section, important kinematical relations for the calculation of the purely deformational part are given. In addition, a method to calculate the rigid-body rotation by means of the Newton-Raphson iteration is described in detail. 2.1. Rotation-free element displacements Two coordinate systems are applied to fully describe the motion of an element. A fixed global coordinate system (Xg, Yg, Zg) measures the global nodal positions, whereas a local co-rotational system (Xl, Yl, Zl) is defined at a node and rigidly translated and rotated with the element. While it moves with the element, it remains an orthogonal frame. The motion of an element from the initial (undeformed) to the current (deformed) configuration is depicted in Fig. 1. Origin 0 of the local element system is here defined at one of the nodes. In the initial configuration, the element has an orientation matrix T0 that remains constant during the analysis. The current deformed configuration can be transformed back to the initial one by matrix (T0 t+Δt R (i) ) T , where t+Δt R (i) is the rigid-body rotation matrix determined in each iteration i within the increment at time t+Δt. Fig. 1 The co-rotational description of the tetrahedral element kinematics By comparing the back rotated with the initial element configuration, the pure deformational displacements of node j are computed as shown in [9]: t Δt (i) t Δt (i) T t Δt (i) t Δt (i) T j 0 j 0 0 j 0 ( ) ( ) ( )        u T R x x T X X , (1) where t+Δt xj (i) and t+Δt x0 (j) are the actual global nodal position of node j and the origin of the co-rotational system (CRS), respectively. Similarly, Xj and X0 represent their counterparts in the initial element configuration. 2.2. Rigid-body rotation In general, rigid-body rotation R at a point of a deformable 3D-continuum can be calculated by the polar decomposition of deformation gradient, F, which is given as: RUF  . (2) 230 V. A NGUYEN, M. ZEHN, D. MARINKOVIĆ where U is the right stretch matrix that describes the non-linear deformation of an infinitesimal small volume around the considered point. For 3D solid finite elements, the deformation gradient is typically determined at the integration points of the element. However, since a linear 3D tetrahedral element is used in this paper, whose shape functions derivatives are constant over the volume of the element, there is a unique rotational matrix that describes the rigid-body rotation between two element configurations. There are various approaches to compute the deformation gradient [1, 6]. In this work, the Newton-Raphson method, as described by Rankin [5], is applied. In the first step, the deformation gradient is estimated at the center of each element, which generally requires a choice of alternative shape functions at the center of the element. The deformation gradient matrix of an element in iteration i within the increment at time t+∆t is given as:     k 1j T(i) j tt j 0T(i)tt xvF , (3) where the summation runs over all element nodes, while constant interpolation function derivative 0 vj of a node j is calculated as: T jjj10 j 0 ggg               Jv , (4) and it depends only on shape functions gj and the initial element geometry [10] and finally, (i) j tt x  is the local nodal coordinate vector that can be given as: t Δt (i) T t Δt (i) j 0 j 0 j ( )     x T X X u (5) In Eq. (4) 0 J is the constant Jacobian matrix of the dimension 3×3, which transforms an infinitesimal small volume element from the natural (ξ, η, ) to co-rotational Cartesian frame (Xl, Yl, Zl) in the undeformed configuration [11]. A 3D finite rotation matrix about a rotation axis r in space can be expressed as the exponential of a 3×3 skew-symmetric matrix S [12]: 2 2 2 sin sin 1 2 exp 2 2                      R S I S S (6) where skew-symmetric matrix S is defined by:               0χ 0ψ χψ0 θspin S (7) and where so-called rotational pseudo-vector rθ  , which represents a finite rotation θ along axis r, is defined by rotations ,  and  around the axes of the Cartesian system, ex, ey and ez, respectively, so that: An Efficient Co-Rotational FEM formulation Using a Projector Matrix 231 zyx eeeθ  (8) Variation of Eq. (6) yields: spin( )  R θ R (9) Using the polar decomposition, Eq. (2), and the orthogonality property of the rotation matrix, the variation of the right stretch matrix is obtained in the same way: spin( )  U θ U (10) and the following relationship can be derived [5]: 2 (axial ) ( trace )    U I U U θ (11) This equation is a basis to develop an iterative process in order to perform the polar decomposition. Generally, the axial operator in Eq. (11) is the inverse of the spin operator and it yields an axial vector a from the skew-symmetric part of a quadratic matrix A by:         T 2 1 axial)spin (axial AAaa (12) Eq. (11) describes a way to calculate the change of pseudo vector δθ with a given right stretch matrix U. Hence, the change of rotation matrix R and the actual right stretch matrix in the next iteration step can be obtained using Eqs. (6) and (10). This iterative process is terminated when U becomes symmetric: 0U  axial (13) The method is performed through the following steps. In the first step, deformation gradient F is computed using Eq. (3) and the right stretch matrix is initialized by: FU  (0) (14) The iterative process is performed as follows: (n) (n) (n) (n 1) 2(axial ) ( trace )     U I U U θ (n 1) (n 1) exp(spin( ))     R θ (15) (n)1)(n1)(n URU   The iteration is terminated, if the magnitude of vector axialU (n) is smaller than a predefined tolerance t: taxial (n) U (16) 232 V. A NGUYEN, M. ZEHN, D. MARINKOVIĆ 3. INTERNAL ELEMENT FORCES AND TANGENT STIFFNESS MATRIX To guarantee the quadratic convergence of the Newton-Raphson method, the internal force vector of each element can be formulated using variation of the internal element energy. On this basis, the tangent element stiffness matrix is estimated by a consistent linearization of the element force with respect to the global nodal displacements. 3.1. Internal element force vector The incremental/iterative Newton-Raphson solution (NRS) for non-linear static structural mechanical problems requires a known solution at time t and estimates the next configuration at time t+Δt using an iterative process. The equations for an assembled finite element structure in the iteration step i can be written as: (i) Gint, Δtt ext tt1)(i G Δtt(i) G Δtt FFuK   1)(i G Δtt(i) G Δtt1)(i G Δtt   uuu , (17) where t+Δt KG (i) is the assembled global stiffness matrix, t+Δt uG (i) and t+Δt ∆uG (i+1) represent the global and incremental nodal displacements while t+Δt Fext and (i) Gint, Δtt F  are the global external and internal force vectors, respectively. In general, the internal force vector can be computed by the variation of the internal element energy with respect to global element displacement vector t+Δt ue (i) as: (i) e tt (i) inte, tt (i) Ge, Δtt W u F       , (18) where t+Δt ue (i) is represented as the global element displacement vector in the form: T T Tt t (i) t t (i) t t (i) t t (i) T e 1 2 k [ ... ]     u u u u . (19) It should be noted here that internal element energy (i) inte, tt W  is unaffected by the rigid-body motion. In the TL formulation, the 2 nd Piola-Kirchhoff stress and the Green- Lagrange strain are used to compute the internal element energy. Because they are invariant with respect to the rigid-body rotation, it can be written [10]: 0 e T t t (i) t t (i) t t (i) 0 e,int e V W ( ) ( ) d V       E S , (20) where (i)tt S  is the 2 nd Piola-Kirchhoff stress matrix expressed in the actual CR frame and d 0 Ve is an infinitesimal small volume in the initial element configuration. In addition, (i)tt E  denotes the Green-Lagrange strain with respect to the actual CR frame as: t t (i) t t (i) t t (i) t t (i) L NL e { }        E B B u , (21) where (i) NL tt B  is the non-linear strain-displacement matrix and is used to describe large strain problems [11]. Matrix (i) L tt B  can be decomposed in two parts: An Efficient Co-Rotational FEM formulation Using a Projector Matrix 233 (i) L1 tt L0 (i) L tt BBB   , (22) where L0 B only depends on the derivatives of the interpolation functions with respect to the natural coordinates in the initial CR system and, hence, is constant. On the other hand, (i) L1 tt B  is estimated using the local nodal displacements from the previous iteration step [11]. In this paper, only small elastic strains are considered and hence, (i) L1 tt B  and (i) NL tt B  can be neglected. In case of small strains, the 2 nd Piola-Kirchhoff stress components are equal to the Cauchy stress components expressed in the CR element frame: (i)tt(i)tt σS   [11]. Hence, the global internal element forces, Eq. (18), can be rewritten as: e 0(i)tt V T L0T(i) e tt T(i) e tt (i) inte, Δtt Vd e 0 σB u u F         . (23) The derivative in front of the integral in Eq. (23) describes the change of the approximate rotation-free element displacement vector with respect to the global nodal displacements and is introduced as a non-linear projector matrix into the global system. Detailed discussions about the roles and derivation of the projector matrix can be found in [1, 4, 13]. Finally, the projector matrix, expressed in the co-rotational frame, is given as: T(i)tt(i) S tt(i)tt    PPIP , (24) where I is the 3k×3k identity matrix and (i) S tt P  is the lever-arm matrix of the dimension 3k×3: t t (i) t Δt (i) t Δt (i) t Δt (i) T S 1 2 k [spin( ) spin( )... spin( )]     P x x x (25) and (i)tt   P represents the non-linear rotation-projector matrix of the dimension 3×3k, which allows to find the best-fit co-rotational element configuration as [1]:                        (i) k tt (i)tt (i) 2 tt (i)tt (i) 1 tt (i)tt T(i)tt ... uuu P . (26) It can be rewritten as: S 0-1(i)ttT(i)tt GUP      (27) with the 3×3 matrix: k t t (i) 0 T t Δt (i) t Δt (i) 0 T j j j j j 1 { ( ) ( )}       U I v x x v (28) and the constant 3×3k matrix: 0 0 0 0 S 1 2 k [spin( ) spin( )... spin( )]G v v v . (29) 234 V. A NGUYEN, M. ZEHN, D. MARINKOVIĆ In the case of linear elastic materials, the co-rotational Cauchy stress is computed using the material constitutive law: (i)tt(i)tt εHσ   , (30) with constant Hooke’s matrix H and linear local strain matrix (i)tt ε  , which is obtained using the rotation-free displacements [7]: (i) e tt L0 (i)tt uBε   . (31) Finally, the internal element forces for a linear tetrahedral element can be written in the closed form: (i) e tt e T(i)tt(i)Δtt 0 (i) inte, Δtt uKPRTF   . (32) The operator A is a 3k×3k diagonal matrix of submatrices A, and eK denoting the element stiffness matrix with: eL0 T L0e VBHBK  (33) where Ve is the volume of a tetrahedral element. Matrix eK is constant and, therefore, calculated once and saved. In this way, the computational effort is reduced. 3.2. Tangent element stiffness matrix The tangent element stiffness matrix, expressed in the global coordinate system, can be calculated by the variation of the global internal element forces with respect to the global element displacements [1]: (i) e tt (i) inte, Δtt (i) Te, tt u F K       . (34) Applying the product rule to Eq. (32), one obtains: (i) e tt (i) e tt e T(i)tt(i)Δtt 0 (i) e tt e(i) e tt T(i)tt (i)Δtt 0 (i) e tt e T(i)tt (i) e tt (i)Δtt 0(i) Te, tt u u KPRT uK u P RTuKP u RT K                     (35) More details on Eq. (35) are available in [1, 4, 5]. The variations in Eq. (35) are: (i) e tt T (i)Δtt 0 (i)tt(i) e tt uRTPu   , (36) t Δt (i) t Δt (i) t t (i) spin( )       R R , (37) T(i)tt(i) S tt(i)tt(i)tt    PPPP . (38) An Efficient Co-Rotational FEM formulation Using a Projector Matrix 235 Substituting Eqs. (36), (37) and (38) into Eq. (35), the resulting global tangent stiffness matrix is generally asymmetric and given as: T (i)Δtt 0 (i) Te, tt(i)Δtt 0 (i) Te, tt RTKRTK   , (39) where (i) Te, tt K  is the local tangent element stiffness that can be expressed as: T Tt t (i) t t (i) t t (i) t Δt (i) t t (i) t t (i) t Δt (i) T t t (i) e,T e e,int e,int [spin ] [spin ]             K P K P F P P F P (40) where t Δt (i) e,int [spin ]  F represents a 3k×3 matrix that contains k spin-matrices of the local nodal force vectors: t Δt (i) t Δt (i) T t Δt (i) T t Δt (i) T T e,int 1,int 2,int k,int [spin ] [spin( ) spin( ) ... spin( ) ]     F F F F , (41) with local element forces (i) inte, Δtt F  computed as follows: (i) e tt e T(i)tt(i) inte, Δtt uKPF   . (42) Simo [14] has shown that the symmetric part of the global tangent element stiffness matrix is sufficient for a quadratic rate of convergence of the Newton-Raphson method, hence: T t t (i) t Δt (i) t t (i) t Δt (i) e,T sym 0 e,T sym 0 ( ) ( )     K T R K T R , (43) where the symmetric part of the local tangent element stiffness is given as: Tt t (i) t t (i) t t (i) t t (i) t Δt (i) T e,T sym e e,int T Tt t (i) t Δt (i) t t (i) e,int 1 ( ) [spin ] 2 1 [spin ] 2               K P K P P F P F P (44) A detailed discussion about the tangent stiffness matrix in Eq. (44) can be found in [1]. 4. CO-ROTATIONAL COMPUTATIONAL PROCEDURE AND NUMERICAL EXAMPLES The proposed co-rotational framework for the calculation of small elastic strains but large rotation structural problems is summarized in Fig. 2 and has been implemented in the commercial FE software package ABAQUS using the user-defined-element subroutine UEL [15]. All constant variables are computed once and saved in the initialization step. The equilibrium of the structure in each increment is computed by an iterative process. The non-linear projector matrix and the element rigid-body rotation are estimated using the global nodal positions in each iteration step. Furthermore, the global tangent element stiffness matrix and the global internal element forces are computed in order to assemble the global system of equations for the whole structure. The solution in the current load increment is obtained when the convergence condition is satisfied. 236 V. A NGUYEN, M. ZEHN, D. MARINKOVIĆ Fig. 2 The presented co-rotational FEM with the projector matrix This paper is focused on the geometrically non-linear static analysis with linear elastic material. All the numerical examples are computed in ABAQUS. The results obtained by the presented CR formulation are compared to the ABAQUS results to verify the accuracy and assess the efficiency of the presented approach. An Efficient Co-Rotational FEM formulation Using a Projector Matrix 237 4.1. A clamped solid block In the first example, the static geometrically non-linear analysis of a clamped solid block exposed to a single force is represented. The solid block, Fig. 3, is made of steel (E=210000 N/mm 2 and ν = 0.3) with the dimensions 400×200×50 mm. The single force F acts at point P in the global z-direction and has the magnitude of 1.5×10 7 N. The FE mesh contains 1420 linear tetrahedral elements C3D4. The resulting deformation is bending-dominated but also accompanied by torsion. Fig. 3 Geometry and boundary conditions of the solid block (left), FE mesh with 1420 C3D4elements (right) In Fig. 4, the development of global displacements of point P with the increasing force is presented. A very good agreement of the obtained results by the presented CR formulation and those from ABAQUS is observable. In addition, the displacements in x-, y- and z-direction of point P are given in Tab. 1. The column entitled “Error” gives the relative difference with respect to the ABAQUS solution. It should be emphasized that this example involves the maximum principal logarithmic strains of 20%. Fig. 4 Displacements in three directions of the solid block for the load case 1.5×10 7 N (left), deformation state with deformation scale factor 1 (right) 238 V. A NGUYEN, M. ZEHN, D. MARINKOVIĆ To assess the efficiency of this formulation, the number of iteration steps used by ABAQUS as well as by the CR formulation is also listed in Tab. 1. The geometrically non- linear analysis was performed in ABAQUS once setting the automatic load increment size (initial increment size: 0.2, minimum size: 10 -5 , maximum size 1) and then again with the fixed load increment size of 0.4. ABAQUS needed 38 iterations with the automatic increment size but did not converge with the fixed increment size used. In contrast, the CR formulation has converged after 24 iterations in the case of automatic increment size and after only 12 iterations with the fixed increment size. The result is reduction of computational time of about 68%. It can be noted that the presented approach also allows for relatively large increments. Table 1 Results for the clamped solid block Cases Abaqus CR-FEM Error X direction -30.78 -30.70 0,3% Y direction -132.12 -132.81 0.5% Z direction 262.28 261.34 0.4% Automatic 38 Iter. 24 Iter. Direct Fails 12 Iter. Max. Log. Strains 20% 20% 4.2. A three dimensional curved beam In the second example, a three dimensional curved beam with a single force F of 5.5×10 5 N at point P is considered. The geometry and the boundary conditions are given in Fig. 5. The curved beam is made of steel (E = 21000 N/mm 2 and ν = 0.3) with the arc length of 0.5×π×450 mm. The cross-section has dimensions of 40 mm × 40 mm. The curved beam is clamped on a one end and discretized by 3652 linear tetrahedral elements C3D4. Again, the deformational behavior can be seen as a combination of bending and torsion. To evaluate the results, the displacements of point P from the geometrically non-linear calculations in ABAQUS and from this CR formulation are given in Fig. 6. A very good agreement is obvious in this case as well. The detailed values are given in Tab. 2. Fig. 5 Geometry and boundary conditions of the curved beam for the load case 5.5×10 5 N (left), mesh with 3652 C3D4 elements (right) An Efficient Co-Rotational FEM formulation Using a Projector Matrix 239 Fig. 6 Displacements [mm] in three directions of the curved beam with 3652 C3D4 elements (left), deformation state with deformation scale factor 1 (right) The maximal error in the displacement components is in the y-direction and reads 1.9%. As for the efficiency, performing the computation with an automatic (initial increment size: 0.2, minimum size: 10 -5 , maximum size 1) and fixed load increment size of 0.3, ABAQUS needed 31 iterations for the both cases. In contrast, the CR-formulation required 24 iterations for the automatic increment size and only 13 iterations for the fixed one. In addition, the computed maximum principal logarithmic strain was in this case 7%. The computational time was reduced by 58%. Table 2 Results for the curved clamped beam Cases Abaqus CR-FEM Error X direction 198.62 201.85 1.6% Y direction -142.98 -140.27 1.9% Z direction 470.03 470.75 0.2% Automatic 31 Iter. 24 Iter. Direct 31 Iter. 13 Iter. Max. Log. Strains 7% 7% 5. CONCLUSIONS The paper presented a co-rotational FE formulation to simulate geometrically non-linear elastic behavior of 3D structures. A systematic derivation of the numerical procedure is given as well as a couple of examples. The basis of the proposed calculation framework is a co-rotational FE-formulation for 3D geometrically non-linear deformations, where the large element rigid body rotation and small elastic displacements can be separately treated in an elegant manner. The element rigid-body rotation is efficiently computed using the polar decomposition of the deformation gradient, where the Newton-Raphson iterative method is applied. The rotation-free nodal displacements, given with respect to the local CR element coordinate system, are assumed to be small. The internal element force vector and the tangent element stiffness matrix are 240 V. A NGUYEN, M. ZEHN, D. MARINKOVIĆ estimated using a consistent linearization of the internal element energy with respect to the global element displacements. This results in a non-linear projector matrix, which generally improves the obtained results. Finally, the presented approach allows the linear tetrahedral element to be easily upgraded into the element for the geometrically non-linear analysis, requiring rather small modifications. A very good agreement between the results from ABAQUS and the proposed CR formulation validates the high accuracy of the presented approach. It was also demonstrated that the computational effort could be significantly reduced by using large load increments. This fact particularly gains importance if the method is used in combination with the advantages of modern hardware tools [16]. The problems in the field of structural mechanics with large rigid-body rotation and with relatively large strains can be efficiently solved using the projector matrix. With those properties, the presented formulation could be embedded into the software packages for multi-body system simulation to provide a greater versatility compared to existing solutions in describing flexible bodies undergoing large rigid-body rotation. REFERENCES 1. Felippa, C. A., Haugen, B., 2005, Unifield formulation of small-strain corotational finite elements: I. Theory, Center for aerospace structures, College of engineering, University of Colorado. 2. Marinković, D., Zehn, M., 2012, Finite element formulations for effective computations of geometrically nonlinear deformations, Advances in Engineering Software, 50, pp. 3-11. 3. Marinković D., Köppe, H., Gabbert, U., 2008, Degenerated shell element for geometrically nonlinear analysis of thin-walled piezoelectric active structures, Smart Materials and Structures, 17(1), pp. 1-10. 4. Nour-Omid, B., Rankin, C., 1991, Finite rotation analysis and consistent linearization using projectors, Computer Methods in Applied Mechanics and Engineering, 93, pp. 353-384. 5. Rankin, C., 2006, Application of linear finite elements to finite strain using corotation, Rhombus Consultants Group, Inc., Palo Alto. 6. Crisfield, M., Moita, F., 1996, A unified co-rotational framework for solids, shells and beams, International Journal of Solids and Structures, 33(20-22), pp. 2969-2992. 7. Moita, G.F., Crisfield M.A., 1996, A finite element formulation using the corotational technique, International Journal for Numerical Methods in Engineering, 39(22), pp. 3775-3792. 8. Espath, L.F.R, Braun, A.L., Awruch, A.L., Dalcin, L.D., 2015, A NURBS-based finite element model applied to geometrically nonlinear elastodynamics using a corotational approach, International Journal for Numerical Methods in Engineering, 102(13), pp. 1839-1868. 9. Nguyen, V.A., Zehn, M., Marinković, D., 2014, Effiziente Berechnung von geometrischen und materiellen Nichtlinearitäten mit einer co-rotationalen Finite-Elemente-Formulierung, Deutschsprachige NAFEMS Konferenz, Bamberg, Germany. 10. Crisfield, M.A., 2000, Non-linear finite element analysis of solids and structures, Volume I, essentials, Johns Wiley & sons LTD, Baffins Lane, Chichester. 11. Bathe, K.J., 1996, Finite element procedures, Prentice-Hall, Inc., Englewood Cliffs, New Jersey. 12. Argyris, J.H., 1982, An excursion into large rotations, Computer Methods in Applied Mechanics and Engineering, 32(1-3), pp, 85-155. 13. Mostafa, M., Sivaselvan, M.V., Felippa, C.A., 2013, Reusing linear finite elements in material and geometrically nonlinear analysis-Applications to the plane stress problems, Finite Elements in Analysis and Design, 69, pp. 62-72. 14. Simo, J.C., 1985, A finite strain formulation. The three-dimensional problem. Part 1, Computer Methods in Applied Mechanics and Engineering, 49, pp. 55-70. 15. Park, K., Paulino, G.H., 2012, Computational implementation of PPR potential-based cohesive model in ABAQUS: Educational perspective, Engineering Fracture Mechanics, 93, pp. 239-262. 16. Nutti, B., Marinković, D., 2014, An approach to efficient FEM simulations on graphics processing units using CUDA, Facta Universitatis series: Mechanical Engineering, 12(1), pp. 15-25.