Operational Research in Engineering Sciences: Theory and Applications Vol. 3, Issue 2, 2020, pp. 74-86 ISSN: 2620-1607 eISSN: 2620-1747 DOI: https://doi.org/10.31181/oresta2003074m * Corresponding author. dragan.marinkovic@tu-berlin.de (D. Marinković), manfred.zehn@tu-berlin.de (M. Zehn), ana.pavlovic@unibo.it (A. Pavlović) HIGHLY EFFICIENT FE SIMULATIONS BY MEANS OF SIMPLIFIED COROTATIONAL FORMULATION Dragan Marinković 1*, Manfred Zehn 1, Ana Pavlović 2 1 Berlin Institute of Technology, Department of Structural Analysis, Germany 2 Department Industrial Engineering, Alma Mater Studiorum University of Bologna, Italy Received: 14 June 2020 Accepted: 15 July 2020 First online: 15 July 2020 Original scientific paper Abstract: Finite Element Method (FEM) has deservedly gained the reputation of the most powerful numerical method in the field of structural analysis. It offers tools to perform various kinds of simulations in this field, ranging from static linear to nonlinear dynamic analyses. In recent years, a particular challenge is development of FE formulations that enable highly efficient simulations, aiming at real-time dynamic simulations as a final objective while keeping high simulation fidelity such as nonlinear effects. The authors of this paper propose a simplified corotational FE formulation as a possible solution to this challenge. The basic idea is to keep the linear behavior of each element in the FE assemblage, but to extract the rigid-body motion on the element level and include it in the formulation to cover geometric nonlinearities. This paper elaborates the idea and demonstrates it on static cases with three different finite element types. The objective is to check the achievable accuracy based on such a simplified geometrically nonlinear FE formulation. In the considered examples, the difference between the results obtained with the present formulation and those by rigorous formulations is less than 3% although fairly large deformations are induced. Key words: Structural analysis, Co-rotational FEM, Geometric nonlinearity, Solid, Shell 1. Introduction Structural analysis is an important engineering discipline encountered in various fields of mechanical and civil engineering. Reliable, accurate and efficient predictions of structural behavior in general, and deformations in particular, are of crucial importance for successful design and optimization of structures, testing their functionality, prediction of their load-carrying capacity and life-time, etc. Recently, this aspect started gaining in importance in some modern fields as well, such as Highly efficient FE simulations by means of simplified corotational formulation 75 interactive simulations, where physics-based simulation are supposed to increase the realism of various applications. Until several decades ago, computations in the field of structural analysis were performed mainly analytically by implementing significant simplifications. Those simplifications have made the models mathematically tractable, but seriously affected the accuracy of the obtained results and therewith their reliability. However, the development of modern hardware tools has set the ambient for development of modern numerical methods for this purpose and the development of computer aided engineering (CAE) software packages skyrocketed. Those programs offered a great assistance to engineers in the previous decades. Among different methods, the Finite Element Method (FEM) has established itself as the method of choice for all problems characterized by complex domains, arbitrary boundary conditions and described by partial differential equations. Problems in the field of structural analysis fit perfectly well into this description and this is why structural analysists have initiated the development of this powerful method and made the very first steps (Turner et al., 1956). Its general applicability to many other engineering fields was later recognized and richly used. Initial developments of FEM were done for the least complicated but quite often encountered problems of structural analysis, namely the linear static problems. Those problems are characterized by slowly increasing loads of constant direction, constant geometric boundary conditions and quite small deformations, so that the initial and deformed structural configurations are almost identical. Those assumptions imply that the balance can be considered over the initial configuration (Bathe, 1996). It is a straightforward task to extend the FE formulations from linear static to linear dynamic cases and it basically comes down to extending the equations by including inertial and damping effects. However, high level of structural optimization implies exploitation of structures to the levels quite close to their limits. In such cases, structural deformations are not small any more and more sophisticated FE formulations were needed to meet the objectives. Total Lagrange and updated Lagrange formulations have set the standards in commercially available FE codes. The essential difference between the two lies in the choice of the reference configuration. Principally it could be any configuration between the initial one and the last determined one, but the common sense choice would be to use either the initial one (total Lagrange formulation) or the last determined one (updated Lagrange formulation). Different strain and stress measures are used in those two formulations and building the tangential stiffness matrix also reflects those differences, but numerics of the two formulations is essentially the same and the choice between the two is basically a matter of taste. Another interesting formulation, namely the corotational FE formulation, appeared several decades ago. Related to FEM the term ‘corotational’ was used for the first time in a paper by Belytschko and Hseih (1979). The idea to cover geometric nonlinearities by attaching a corotational frame to single elements was introduced by Horrigmoe and Bergan (1978). The work in this direction continued under the supervision of Bergan and the developments were summarized in a survey article by Nygard and Bergan (1989). Crisfield (1990, 1997) and Crisfield & Moita (1996) introduced “consistent CR formulation” by developing the stiffness matrix as the Marinković et al./Oper. Res. Eng. Sci. Theor. Appl. 3 (2) (2020) 74-86 76 actual variation of the internal force. Rankin and Brogan (1986) proposed the concept of element independent corotational formulation. A high-quality survey of these developments and a detailed analysis of their properties was provided by Felippa and Haugen (2005). This paper suggests a corotational (CR) FE approach that offers a trade-off between numerical efficiency and achievable accuracy by simplifying the rigorous corotational FE formulation. The idea is to offer a FE formulation that would be of high interest for specific applications such as multi-body system dynamics, where parts exhibit large rigid-body deformations but small strains, or applications involving real-time simulations. In this paper, the achievable accuracy will be tested with a few basic solid and shell elements in cases involving large local rigid-body rotations. 2. Simplified CR formulation and implemented elements 2.1. Basic principles of the simplified CR formulation While the linear FE formulations offer very efficient and stable computations, the nonlinear formulations are very time consuming and prone to computational stability issues, as they might not necessarily produce converged solutions. On the other hand, linear formulations are accurate only for small deformations, but geometrically nonlinear formulations offer engineering accuracy for large deformations involving arbitrarily large rigid-body rotations. While one would wish to have advantages of both formulations in one formulation, it is certainly not possible to have all the advantages to the full extent. But a formulation may offer a kind of trade-off or a compromise between those. The formulation that will be explained here follows the idea of element-based CR formulation. Hence, the basic concept is to attach a coordinate system to an element and keep the linear FE formulation of the element with respect to this coordinate system. The attached coordinate system follows the element in its rigid-body motion. As the elastic behavior of the element remains linear with respect to the attached coordinate system, this implies that the element matrices are computed only once for this coordinate system. As deformation proceeds, it is necessary to determine the motion of the attached coordinate system, or, in other words, to determine the element rigid-body motion. Once this is described by the element rotational matrix, Re, the related element matrices and vectors can be rotated to the current configuration and the assemblage of the global matrices and vectors can be done for further computation. Hence, the element elastic behavior is described as linear with respect to the attached corotational frame and the element stiffness matrix with respect to this frame is not updated. In this manner, the local element deformation and the stress state is neglected from the consideration of geometrically nonlinear effects, thus simplifying the formulation significantly compared to the rigorous nonlinear formulations. The formulation keeps the very important aspect of geometric nonlinearity, namely the rigid-body rotation that is accounted for on the element level. In continuum every point exhibits its own rigid-body rotation, generally speaking. Obviously, this aspect is described here in a coarser way, as it is always the Highly efficient FE simulations by means of simplified corotational formulation 77 case with methods that apply discretization. But one may arbitrarily adjust the ‘resolution’ of accounting for rigid-body rotation by performing finer or coarser FE meshing. Hence, assuming the element rigid-body rotation between the initial and the current configuration is known and given as the rotation matrix Re, the element stiffness matrix, Ke, is update with respect to the global coordinate system in a straightforward manner: Tt0tt eeee RKRK = , (1) where the left superscript denotes the moment in time at which the term is given. This simple way of updating the element stiffness matrix is where the efficiency of the method resides. Not only is the tangential stiffness matrix efficiently updated, but also its condition number does not change dramatically in this manner, so that the stability of computation is kept to a large extent. This is not always the case with rigorous nonlinear FE formulations in which single elements may suffer significant deformations, and, as a consequence, the solution may not converge. In order to perform nonlinear computations, one needs the tangential stiffness matrix, the update of which was elaborated above, and the internal forces. In order to determine the internal forces, deformational displacements and rotations are required. Those are obtained when the rigid-body rotation is removed from the overall element displacements. This procedure is best explained using Figure 1. In this figure, a tetrahedron element is shown in its original and an arbitrarily deformed configuration. Again, it is assumed that the rigid-body rotation of the element is known. It is sufficient to rotate the deformed element back to the original element configuration. By comparing so obtained element configuration with the initial one, one obtains deformational displacements. Figure 1. Extraction of deformational displacements for a tetrahedron element Hence, the expression for the deformational displacements reads: eee xxRu 0tTt R t 0 −= , (2) where xe denotes the element configuration, i.e. those are the nodal coordinates of all element nodes. With the known deformational displacements, one may simply multiply those with the element stiffness matrix for the initial configuration to obtain the internal forces with respect to the initial configuration. The internal forces are Marinković et al./Oper. Res. Eng. Sci. Theor. Appl. 3 (2) (2020) 74-86 78 the rotated to the current configuration. Those steps are summarized in the following expression: eeeeeeeeee xKRxRKRFRF 00ttTt0tt 0 tt t −== , (3) If finite elements contain rotational degrees of freedom, the procedure is essentially the same for rotations, as illustrated in Figure 2. Figure 2. Extraction of deformational rotations for a shell element With the incremental nodal rotations, t-ti1, t-ti2 and t-ti3, the element nodal normals are updated by means of the rotation matrix Qi (Argyris, 1982): ( ) ( ) 2 i tΔt 2 i tΔt i tΔt t i tΔt i tΔt i tΔt 2/γ 2/γsin 2 1 γ γsin SSIQ − − − − − −         ++= , (4) where 2 3i tΔt2 2i tΔt2 1i tΔt i tΔt θΔθΔθΔγ −−−− ++= , (5)           − − − = −− −− −− − 0θΔθΔ θΔ0θΔ θΔθΔ0 1i tΔt 2i tΔt 1i tΔt 3i tΔt 2i tΔt 3i tΔt i tΔt S , (6) so that i tΔttΔt i t nQn −− = . (7) After rotating the element from the current to initial orientation, it is a straightforward task to compute the internal moments and rotate them again to the current configuration in an analogous manner as done above with the forces. 2.2. Finite elements implemented into the formulation So far three finite elements have been implemented into the proposed simplified CR formulation – two solid elements and one shell element. Highly efficient FE simulations by means of simplified corotational formulation 79 The solid elements are the linear tetrahedron element and the quadratic hexahedron element. The linear tetrahedron element is notorious for its too stiff behavior and is, therefore, often avoided in modeling. However, it has two nice properties. It is numerically very efficient and can discretize any geometry. Actually, the second advantage makes it often inevitable in FE models in order to model some Figure 3. Linear tetrahedron element (left) and quadratic hexahedron element (right) areas of the model that are otherwise too difficult to discretize. In addition, this element is characterized by unambiguity of the rotation matrix. There is a single, unique rotation matrix describing the rigid-body rotation of this element, which is not the case with most of the finite elements. For any two given configurations of the element there is a unique matrix that transforms the element from one configuration to the other one. This is due to the fact that the element employs the linear shape functions, so that the deformation gradient has a constant value over the whole element domain. Polar decomposition of this transformation matrix yields the rotation matrix. In order to obtain reasonable results with this element, a quite fine discretization is required. But this also increases the “resolution” of accounting for the rigid-body rotation, which is a positive aspect regarding the corotational FE approach. Nguyen et al. (2016) have used this element in combination with a corotational FE approach that implements the projector matrix for the sake of better result convergence. The quadratic hexahedron element is in most aspects the opposite of the linear tetrahedron element. It offers the best accuracy among solid elements (apart from those that use special techniques), but is numerically very demanding and requires partitioning of complex geometries for successful meshing, whereby the ‘corners of the geometry’ will still require tetrahedron elements. The rotational matrix is not unique for the element, i.e. it differs for different points within the element domain. Hence, it is ambiguous and one has to decide what strategy to use in order to determine it. It may be determined by local coordinate systems defined in a special ways by using the current nodal positions. A better option would be to use the deformational gradient at some point of the element, whereby the element centroid appears to be a natural choice – exactly this option was applied in this work. The best, but also the most demanding option would be to obtain some kind of average rotational matrix of the element. It is so far, however, an open question with respect to what criteria the averaging is to be performed. The implemented shell element is a linear triangular shell element (Figure 4) recently developed (Rama et al., 2018, 2018a, 2018b, Marinkovic et al., 2019). Essentially, the element is a combination of a plate element and a membrane element. It implements the Mindlin-Reissner kinematics and uses the Discrete Shear Marinković et al./Oper. Res. Eng. Sci. Theor. Appl. 3 (2) (2020) 74-86 80 Gap (DSG) method (Li et al., 2019) in combination with the strain smoothing technique to alleviate the notorious shear locking. The membrane part of the element is actually the ANDES membrane formulation developed by Felippa and Militello (1992). Similarly to the linear tetrahedron element, this one also has a unique rotational matrix, can discretize any surface geometry and is numerically highly efficient. Due to the flat shape of the element, the discretized geometry is actually faceted, which affects the achievable accuracy. Figure 4. Linear triangular shell element 3. Numerical examples In what follows, three examples of large deformations, each for one type of implemented elements, will be considered in order to investigate the achievable accuracy by means of the proposed corotational FE formulation. The major purpose is the comparison of computed displacements obtained by rigorous geometrically nonlinear FE formulation (computed in Abaqus) and those obtained by the presented development. In accordance with this objective, all quantities will be given as dimensionless. The selected examples are of academic nature involving structures of rather simple geometry. Sufficiently large loading will be chosen to produce geometrically nonlinear deformations, i.e. those that significantly differ from deformations computed by the linear formulation. 2.1. Solid elements The same structure, which may be referred to as a block, with dimensions 10101.5 and clamped over one surface with dimensions 101.5 will be discretized with both tetrahedron and hexahedron element. The geometry with kinematic boundary conditions is depicted in Figure 5, left. The material is linear elastic with the following properties: Young’s modulus E=21011 and Poisson ration =0.3. The load cases are chosen to be different for the discretization with the tetrahedron element and for the discretization with the hexahedron element. In both cases the force is set to be F=1010 in order to produce sufficiently large, geometrically nonlinear deformations. As shown in Figure 5, middle, in the case of discretization with the tetrahedron element, the force acts only at one corner of the structure so as to bend and twist it at the same time. Figure 5, right, shows discretization with the hexahedron element and the load case with a pair of oppositely oriented forces that cause twisting of the considered structure. Highly efficient FE simulations by means of simplified corotational formulation 81 Figure 5. Block structure: geometry and kinematic boundary conditions (left), load cases and discretization with tetrahedrons (middle) and hexahedrons (right) In order to visualize the deformations and get the feeling for the magnitude of deformation, the unscaled deformations (i.e. scale factor set to 1) are depicted in Figures 6 (the FE model with tetrahedron elements and one force) and 7 (the FE model with hexahedron elements and two forces) together with the undeformed structure. The structure is shown from different perspectives. Obviously, the magnitude of deformation is well beyond the realm of linearity. Figure 6. Deformed and undeformed block structure under single force load, discretization with tetrahedron elements, three different perspectives Figure 7. Deformed and undeformed block structure under force couple load, discretization with hexahedron elements, three different perspectives As a representative point to follow its displacements with the gradually increasing loading, the point at which the force acts in Figure 5, middle, is selected. Its displacements in all three global directions are considered in both cases and compared with the linear and geometrically nonlinear results from Abaqus. The Marinković et al./Oper. Res. Eng. Sci. Theor. Appl. 3 (2) (2020) 74-86 82 results are given in Figure 8, for the case shown in Figure 5, middle, and Figure 9, for the case shown in Figure 5, right. One may notice that the results for the same global displacements show a similar trend in both considered cases. Comparing the displacements in the same directions in those two cases, one would notice that the major distinction is in the relative difference between linear and geometrically nonlinear results. There is actually a significant difference between the linear and geometrically nonlinear results, and it goes even up to 50%. This was expected and, Figure 8. Model with tetrahedron elements – displacements in three global directions Figure 9. Model with hexahedron elements – displacements in three global directions in fact, the loading was chosen with this objective. On the other hand, there is also a good agreement between the nonlinear results by Abaqus and present formulation. The difference is observable in the last 30-40% of the loading but stays in the limits of up to 2%, which is an acceptable result for many different applications. In addition, the highest difference is at the full loading, where local element deformations start to kick in and this is an effect not accounted for by the present formulation. As long as this effect is not present, the difference in the results is practically negligible. 2.2. Shell element The example considered for the shell element is a typical benchmark case used in development of shell elements for geometrically nonlinear analysis. It a straight beam, with one end clamped, while the free end is exposed to a bending moment of such a magnitude that the beam bends into a circle. The moment required to produce such a deformation can be computed analytically assuming beam kinematics and is given as M = Ebh3/6l (Bathe and Bolourchis, 1979), where E is the Young’s modulus, while b, l and h are the width, length and thickness of the beam, respectively. In this Highly efficient FE simulations by means of simplified corotational formulation 83 case, the Young’s modulus is the same as in the previous cases with solid elements, while Poisson ratio is set to zero, so that the shell element reproduces the beam behavior (Poisson effect is neglected over the beam’s cross-section). Dimensions of the beam can be seen in Figure 10, left, while Figure 10, right, shows the FE mesh applied. Figure 10. Beam model and discretization with triangular shell elements Interestingly, Abaqus encounters a problem to complete the computation with its 3-node shell element. The computation runs until approximately 95% of the load, and when this level is reached, a converged solution is not found any more (automatic stepping was used to facilitate the computation). Figure 11 shows the initial and deformed configuration as computed by Abaqus. The relatively coarse mesh is the reason for the faceted deformed geometry and could be one of the reasons for the computational issues encountered by Abaqus. Figure 11. Beam model – deformed (Abaqus) and undeformed configuration The results for the displacements along the x- and y-axes are given in Figure 12. The diagrams include only nonlinear results as the large difference between those and linear results would make the inclusion of linear results unreasonable. The computation with the proposed corotational formulation proceeds till 100% of the load. One may notice a good agreement of the results up to the load level computed by Abaqus. Marinković et al./Oper. Res. Eng. Sci. Theor. Appl. 3 (2) (2020) 74-86 84 Figure 12. Beam model – free end displacements in the x- and y-directions 4. Conclusions The paper elaborates a simplified corotational FE formulation in which the element local behavior is described as linear, but the element rigid-body rotation is accounted for. Hence, it neglects the effect of the element pure deformation and element stress state on the element tangential stiffness matrix. This is where the numerical savings are made, thus rendering the formulation very efficient in terms of computational effort and also numerically stable. At the same time, this means that it delivers results that are an approximation compared to the results delivered by the rigorous geometrically nonlinear FE formulations, such as total and updated Lagrange formulation, or the rigorous corotational FE formulation. The examples were focused on accuracy of predicting the structural displacements, which is equivalent to the accuracy of predicting the deformed structural configuration. It was shown in the considered examples that the discrepancy between the rigorous results and those obtained by the proposed formulation is only a few percent for fairly large deformations. Of course, the achievable accuracy certainly depends on the nature of the deformation and is expected to be better in cases where local rigid-body rotation dominates. Furthermore, this means that the formulation can successfully be used for certain engineering simulations where this level of accuracy is acceptable. For instance, it can be a very attractive alternative for consideration of elastic bodies in Multi-Body System (MBS) simulations, which is currently mainly done based on the modal- superposition technique thus covering only linear deformations with respect to the local frame of the whole structure. The proposed formulation would offer better accuracy and fidelity of the full-scale FE model, while keeping the numerical effort in acceptable limits. Another interesting field of application would be Virtual Reality (VR) where physics-based real-time simulations have always played a challenging task (Marinkovic et al., 2018, Marinkovic & Zehn, 2019, Zehn & Marinkovic , 2019). In this field, the presented formulation can be successfully used for various types of simulators such as surgery (Marinkovic & Zehn, 2018), assembly planning and practicing assembling of various complex products, thus improving the productivity, etc. Highly efficient FE simulations by means of simplified corotational formulation 85 References Argyris, J. (1982). An excursion into large rotations. Computer Methods in Applied Mechanics and Engineering, 32(1-3), 85-155. Bathe, K. J. (1996). Finite element procedures. New York : Prentice Hall. Bathe, K. J., & Bolourchis, S. (1979). Large displacement analysis of three- dimensional beam structures. International Journal of Numerical Methods in Engineering 14(7), 961-986. Belytschko, T., & Hsieh, B. J. (1979). Application of higher order corotational stretch theories to nonlinear finite element analysis. Computers & Structures, 11, 175–182. Crisfield, M. A. (1990). A consistent corotational formulation for nonlinear three- dimensional beam element. Comp. Meths. Appl. Mech. Engrg., 81, 131–150. Crisfield, M. A. (1997). Nonlinear finite element analysis of solids and structures. Vol. 2: Advanced Topics. Chichester: Wiley. Crisfield, M. A., & Moita, G. F. (1996). A unified co-rotational for solids, shells and beams. Int. J. Solids Struc., 33, 2969–2992. Felippa, C., & Haugen, B. (2005). A unified formulation of small-strain corotational finite elements: I. Theory. Computer Methods in Applied Mechanics and Engineering, 194(21-24), 2285-2335. Felippa, C., & Militello, C. (1992). Membrane triangles with corner drilling freedoms – II The ANDES element. Finite Elem. Anal. Des., 12(3–4), 189–201. Horrigmoe, G., & Bergan, P. G. (1978). Instability analysis of free-form shells by flat finite elements. Comp. Meths. Appl. Mech. Engrg., 16, 11–35. Li, S., Zhang, J. & Cui, X. (2019). Nonlinear dynamic analysis of shell structures by the formulation based on a discrete shear gap. Acta Mechanica (230), 3571–3591. Marinkovic, D., Rama, G., & Zehn, M. (2019). Abaqus implementation of a corotational piezoelectric 3-node shell element with drilling degree of freedom. Facta Universitatis, series: Mechanical Engineering, 17(2), 269-283. Marinkovic, D., & Zehn, M. (2018). Corotational finite element formulation for virtual-reality based surgery simulators. Physical Mesomechanics, 21(1), 15-23. Marinkovic, D., & Zehn, M. (2019). Survey of finite element method-based real-time simulations. Applied Sciences, 9(14), art. no. 2775. Marinkovic, D., Zehn, M., & Rama, G. (2018). Towards real-time simulation of deformable structures by means of co-rotational finite element formulation. Meccanica, 53(11-12), 3123-3136. Nguyen, V. A., Zehn, M., & Marinkovic, D. (2016). An efficient co-rotational FEM formulation using a projector matrix. Facta Universitatis, series: Mechanical Engineering, 14(2), 227-240. Nygard, M. K., & Bergan, P. G. (1989). Advances in treating large rotations for nonlinear problems. In A. K. Noor, & J. T. Oden (Eds.), State-of the Art Surveys on Computational Mechanics (pp. 305–332). New York: ASME. Rama, G., Marinkovic, D., & Zehn, M. (2018). Efficient three-node finite shell element for linear and geometrically nonlinear analyses of piezoelectric laminated structures. Journal of Intelligent Material Systems and Structures, 29(3), 345-357. Marinković et al./Oper. Res. Eng. Sci. Theor. Appl. 3 (2) (2020) 74-86 86 Rama, G., Marinkovic, D., & Zehn, M. (2018a). A three-node shell element based on the discrete shear gap and assumed natural deviatoric strain approaches. Journal of the Brazilian Society of Mechanical Sciences and Engineering, 40(7), art. no. 356. Rama, G., Marinkovic, D., & Zehn, M. (2018b). High performance 3-node shell element for linear and geometrically nonlinear analysis of composite laminates. Composites Part B: Engineering, 151, 118-126. Rankin, C. C., & Brogan, F. A. (1986). An element-independent corotational procedure for the treatment of large rotations, ASME J. Pressure Vessel Technology, 108, 165– 174. Turner, M. J., Clough, R. W., Martin, H. C., & Topp, L. J. (1956). Stiffness and deflection analysis of complex structures, J. Aeronaut. Sci., 23, 805-824. Zehn, M.W., Marinkovic, D. (2019). Chances of real-time simulation in FE analyses with conventional hardware. Advances in Engineering Materials, Structures and Systems: Innovations, Mechanics and Applications - Proceedings of the 7th International Conference on Structural Engineering, Mechanics and Computation, 2019, Cape Town, South Africa, pp. 531-536. © 2020 by the authors. Submitted for possible open access publication under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). highly efficient FE simulations by means of simplified corotational formulation Dragan Marinković 1*, Manfred Zehn 1, Ana Pavlović 2 1. Introduction 2. Simplified CR formulation and implemented elements 2.1. Basic principles of the simplified CR formulation 2.2. Finite elements implemented into the formulation 3. Numerical examples 2.1. Solid elements 2.2. Shell element 4. Conclusions References