A simple fully nonlinear Kirchhoff-Love shell finite element

*Corresponding author https://doi.org/10.1590/1679-78256120 Abstract The current paper implementates a simple fully non-linear Kirchhoff-lovel shell penalty based finite element. The 6 nodes and 21 DoF triangular element developed in this work has a quadratic displacement field associated to it and the C 1 continuity required by Kirchhoff-Love Hyphotesis is approximated by an internal penalty. The biggest novelty in this article is the simultaneous use of penalty and a Rodrigues incremental Rotation parameter (scalar DOF) between neighboring elements further explained in the text. The nonlinear finite element model developed in this article is compared to analytical results, commercial finite element code and another FEM model developed in bibliography. Simulations have demonstrated consistency when comparing results to other models and it is deemed that reliable mesh generation together with a powerfull triangular finite element is a good option for trustworthy Love shell,


INTRODUCTION
A numerical approach solution usually is an effective method of obtaining the structural response for situations where the load distribution, material non-linearity or geometrical properties are complex, thus turning the attainment of an analytical solution into a real challenge. Shells are a very important structural model, which can be regarded as a MOR (Model Order Reduction) of Solid Mechanics (Nigro et al., 2019) and are consistently applicable to engineering structural problems.
Theories for Plates, Shells and Membranes have been developed in the past decades by many researchers due to its importance in mechanics. Kirchhoff-Love and Reissner-Mindlin theories are example of shell/plate models used in engineering. In the current paper, the authors address a non-linear Kirchhoff-Love element implemented in the AceGEN / AceFEM code environment, which is a Mathematica based code suitable for developing optimized code for computational mechanics (Korelc, J.A., Wriggers, P., 2016).
In Finite Element Methods, this structural problem may be implemented and studied by three-dimensional elements. As stated in Korelc and Wriggers (2016) this strategy is, however, inadequate for certain occasions. It is shown that for especially thin shells, the locking-phenomena may occur, causing the simulated structure to appear stiffer than what it really is. Hence, special attention must be taken regarding the interpolation of 3-D elements for Shell modeling.
Another approach that have been used by some researches is the use of kinematical simplification for shells. With these models, the system can be modeled as a surface (2D) with a specified thickness, consequently with less degrees of freedom, but can yet deal with nonlinear phenomena and with finite deflections and Rotations. Some pioneer work has been developed by Simo and Fox (1989) and Pimenta (1993) developing new shell models before 2000's and departing from the traditional so-called degenerated solid approach. Campello et al. (2003) and Pimenta et al. (2004) developed shear flexible models (Reissner-Mindling) based on kinematics proposed by Pimenta (1993). These works had successfully implemented a shear flexible nonlinear six-parameter triangular shell finite element. In the same work, finite rotations were treated by using Euler-Rodrigues formula in a very conviniet way (due to its incremental formulation). Viebahn et al. (2017), Costa e Silva et al. (2018), Costa e Silva (2020) and Costa e Silva et al. (2020) built finite element models based on Kirchoff-Love Hypothesis and are the basis for the current article. The main advantage for using Kirchoff-Love shells is that the element will carry less degrees of freedom per element, simplifying the numerical problem. However, as shown in the next sections, one must ensure C1-continuity in domain due to kinematical formulation. As the displacement field is quadratic, this continuity is guaranteed inside the element. The present work and the work of Viebahn et al. (2017) and Costa e Silva (2020) differ in the methodology to enforce C1 continuity in the element interfaces. The procedure adopted in this article can be regarded as a discrete Galerking approach (see Brezzi and Fortin, 1991).

OBJECTIVES
In the current paper, a simple non-linear finite element for Kirchhoff-Love shells undergoing finite deflections and rotations is developed and it is implemented with the use of AceGen and AceFem. In the first part of this paper the kinematical description of the model is made and in the second part three numerical simulations are done and results are comparared to outputs from commercial software (Nastran) and from the model developed in Costa e Silva (2020). The simulations are performed with different dimensions, mesh discretization and loads in order to evaluate its reliability and response for different system conditions.

Shell Kinemátics
Plates and shell theories are very diverse, and many models may be found in the bibliography considering different kinematical hypothesis. The current work is based on the Kirchhoff Plate theory which has three kinematical constrains, also known as Kirchhoff hypothesis (see Reddy, 2006): 1. Straight perpendicular lines (transverse normals) to the middle surface stay straight after any deformation.
3. These transverse normals remain perpendicular to the middle surface during any deformation imposed to the shell.
There is no shear strain in any shell cross section.
The Kirchhoff-Love shell model, also known as shear rigid because of the absence of shear strain in consequence of its kinematical assumptions, developed in this paper may be understood as a continuation of the work developed by Costa e Silva et al. (2018), Costa e Silva et al. (2020). and Viebahn et al. (2017). It includes a new methodology to approcimate C1 continuity between adjacent elements with penalty and is also is valid for finite deflections, rotations and strains.
A plane reference configuration is assumed, however initial curved shells may be represented by a a stress-free deformation (Pimenta and Campello, 2009). The model is illustrated by Figure 1: The main kinematic properties of the shell developed in this article are based on Pimenta, P.M., et al. (2010) and are defined by:

•
An arbitraty point is parametrized in the reference configuration by the vector ∈ ℝ 3 and in the current configuration by the vector ∈ ℝ 3 .

•
The projection of the arbitrary point in Middle surface in the reference configuration is parametrized by the vector ∈ Ω ⊂ ℝ 2 and in the current configuration by the vector ∈ Ω ⊂ ℝ 3 .

•
The vector that connects the arbitrary point to its projection in the middle surface is and in the current configuration by . It is important to notice that and are perpendicular to Ù and Ù respectively.

•
The shell volume and thickness of the reference shell is defined such as and = [−ℎ , ℎ ]. The Shell total thickness is ℎ = ℎ + ℎ An orthonormal right-handed coordinate system is defined in the reference configuration such as any point in the shell may be defined by equation (1) (Pimenta, P.M., et al., 2010) = + , = , ∈ Ω r and = 3 , 3 ∈ . (1) Equivalently, an arbitrary point in the current configuration can also be defined in an associated orthonormal coordinate system by position vector , rotation tensor and displacement with expressions = + , = and = + .
It is important to mention that tensor rotates vector to , which are always perpendicular to shell middle surface, consequently is important to describe shell curvature and related kinematical quantities. The approach of rotational quantities is better explained in further sections of this paper.
As one of the kinematical assumptions of Kirchhoff-Love theory is that the director vector remains orthogonal to the middle surface of the shell, then the rotation tensor may be defined by The local orthogonal coordinate system in the current configuration is expressed by and its equation is defined in further section(See equation (30)).
The nonlinear deformation map that takes a vector in the reference configuration and takes it to the current configuration at any positive time is defined by : ξ → x and the deformation gradient = t ( ) is defined by where the vectors are defined by = , + , and = = .
The Curvature vectors are defined by where can yet be rewritten as Consequently, one may write the components of deformation gradient as and the Jacobian and cofactor of may be written as Back-rotated counterparts of the deformation gradient and of strains may also be defined by Latin American Journal of Solids and Structures, 2020, 17(8), e325 5/16 The plane stress condition is forced considering vanishing stresses in the normal direction of the mid-plane by the expression where, τ is the decomposition of the first Piola-Kirchhoff stress tensor on the cartesian axes. It is defined by In the previous equation, = ( ) is the so-called Helmholtz free energy per unit of reference volume.

Weak form of equilibrium
In solid mechanics, the balance of linear momentum in statics (LAI et al., 2009) may be defined by the differential equation (14) as However, classically, the mathematical model implemented on FEM code is based on the weak form of the equilibrium equation. The weak form (in this case, variational formulation) can be obtained applying the principle of virtual work (Wriggers, 2008) (Viebahn et al., 2017) and consequentily we have where: • δw is virtual work of external forces, • δw is the virtual work of internal forces, • is the first Piola-Kirchhoff stress tensor • δ is the virtual deformation gradient (funcition of virtual displacements) • is the displacement; δ is the virtual displacement.
• is the boundary forces at Domain's boundary.
• is the volumetric external force. The internal virtual work at Equation (15) may be further developed (Pimenta et al., 2010) by Latin American Journal of Solids and Structures, 2020, 17(8), e325 6/16 and the external virtual work at Equation (15) may also be rewritten as where has been considered

Constitutive equations (Elasticity)
The model implemented in this article uses an elastic material model (neo-Hookean material) suitable for large displacements and deformations (Pimenta, P.M, 2006). The strain energy is given by equation (21) = where are the principal invariants of the deformation Cauchy-Green right tensor "C", defined by with: • and μ -Lamé coefficients (Material properties) • = ( ) = Helmholtz free energy • = = deformation Cauchy-Green (right) tensor.

Spatial Discretization
The finite element developed is triangular. Because of the Kirchhoff-Love shell theory, it would be required 1 continuous approximations in order to preserve continuity of the displacement gradient among all domain. This model, however, imposes 1 displacement field inside the element and enforces rotation between adjacent elements. Emforcement of kinking angle between adjacent elements is better explained in further section.
In the current model, shape functions are based on Area coordinates (Papadopoulos P., 2015). The interpolations for position vector (reference configuration), displacement vector and derivatives are defined by equations (23) and (24), which are similarly represented in Viebahn et al. (2017).
Expressions (24) represents the discretization of the variation and linearization of displacements and spatial derivatives It is important to remark that the current Paper does not focus time-dependent problems. Further development in this topic will be approached in future paper under development. Figure 2 illustrates the finite element implemented in this article with its degrees of freedom in its 6 nodes. It is assumed a quadratic displacement field (with nodal values of (1) . . (6) ) and 3 independent rotation parameters ( Δ (4) , Δ (5) , Δ (6) ).

Kinking angle enforcement
The kinking angle enforcement in the present article follows a methodology similar but different from Costa e Silva et al. (2018) and Viebahn et al. (2017). The shell kinematics based on Kirchhoff-Love model imposes the necessity of continuity of the displacement gradient (See equations (8), (16) and (18)). Consequently, the finite element construction should guarantee 1 continuity in the whole domain, include inside and the edges of the elements. The finite element structure implemented in this article penalizes the kinking angle(See Figure 3) between adjacent elements based on one of the so-called Rodrigues rotation parameters developed in Costa e Silva (2020).
As explained in Viebahn et al. (2017), the enforcement of C1 continuity at element edges is expected to be satisfied asymptotically as ℎ → 0 (mesh refinement) if one guarantees that the kinking angle(at the midpoint) between adjacent elements does not change during the deformation. Geometrically, the kinking may be understood as the angle between mid-surface normal vectors of adjacent elements.
As seen in equation (2), the tensor rotates the normal vector of the shell from reference to current configuration. As will be shown in further equations, if one imposes that the rotation tensor (or its incremental part) is the same for both adjacents elements, one may guarantee the continuity of the kinking angle. In this sense, the author uses one of the characteristics of the Rodrigues rotation parameters, in which the tensor may be written as function of : When working with Rodrigues parameter in rotations, rotation increment is easily treated by equation where: With the previous expressions, one may write the coordinate system in sequential time steps ({ 1 , 2 , 3 } at and { 1 +1 , 2 +1 , 3 +1 } at +1 ) by 3 +1 = ∆ 3 , 3 +1 − 3 = ∆ × 3 and 3 = 1 2 � 3 +1 + 3 � .
In (Costa e Silva, 2020), two important expressions are achieved regarding rotation parameters. It is shown that The previous equations shows that the normal vector at the edge of the element at time +1 may be find if the parameter Ä is known, but at the same time one may determine with the displacement field and so far no condition has been imposed to guarantee that both normal vectors are the same. In this sense, the penalty approach is used to guarantee that the displacement field is compatible with the Ä independent parameters, witch have been introduced, before, as a mechanism to guarantee 1 continuity with sufficient mesh refinement. The following equations are used , , With the previous equations, one may write based on the displacement field only and consequently Δ may be defined by Here, the penalty is summed to Helmholtz free energy to create the penalty function: As a note, it is important to remark that it is used similar letters in this paper ( and ) because the represent the same parameter, however, as explainded before, they are calculated in different ways and are forced to be the same by equation (34).
In the previous equation, it is used the only artificial numerical factor of the formulation ("k").
In the examples, it is assumed as a multiple of the bending stiffness ( ℎ 3 /�12(1 − 2 )�). Figure 4 illustrates adjacent elements with its degrees of freedom (DoF).

RESULTS
The FEM model have been implemented in AceGen / AceFEM environment and tested for 3 simple cases: Simply supported square plate, Cantiliever Beam and Pinchided Cylinder. For convenience, the finite element model developed in this article has been nominated as "KL T63i Pen" by similarity to the model developed in Campello et al. (2003). It is used the letters Pen as "Penalty", in contrast to Lag in "Lagrangean" -Model developed in Costa e Silva (2020).

Square Plate
The square plate is simply supported among all 4 sides and an external normal force is applied in the center of the plate. This is a classical problem in theory of elasticity and analytical results may be found in bibliography for small displacements. The formulas used in the simulations have been refered to Timoshenko and Krieger (1959) and Young et al. (2002). The simulation has been executed for different mesh discretization and for different Thickness and its output is compared to analytical results and to commercial software Nastran Triangular Shell Element -Ctria6 (See Autodesk Inc, 2019).
In the simulations and calculations, mechanical properties, dimensions and loads are defined by • Center down Force = 10 N; • a(width) =1 m; • b(length) = 1 m; • h(thickness) = 1mm and 0,2mm;    (Figure 6) plots the total displacement of the center of the plate for different load factors. Each line represents different models used to simulate. The second (Figure 7) represents tha same simulation but with different thickiness, the last graph (Figure 8) in this subsection compares maximum displacement (at center) for different mesh discretization.   In these simulations it has been observed that for K-L Lagrange model and Nastran Ctria6 the simulations took longer time to converge with finer meshes. In addition, it has been noticed that the K-L Penalty based element has performed the faster simulations (compared to the other models with same parameters). Another interesting fact is that the author could not simulate successfully the problem using 3D solid tetrahedral Nastran elements for given thickness. In this case, the software was unable to solve the numerical problem and presented several warnings and errors. This was expected due to the locking phenomena discussed in section 1 of this paper. Figure 6 demonstrate great results agreement between different models. For 1mm thickness the model developed in this article had practically the same results as the model developed by Costa e Silva (2020) and to the commercial software Nastran, with differences smaller than 5%. The analytical solution could only approximate the results at the beginning of the curve due to its linearized form. For the thinner(0,2mm) plate (Figure 7) it may be seen that the models present some slight differences. The (vertical) displacements at plate center are very similar, however Nastran FEM have shown as the stiffer result. For a load factor of 1, the displacements are less than 10% different from one another. Here again the analytical formula is suitable only for very low displacements. Figure 8 shows that all the models tend to converge to a maximum value when using a bigger number of elements, as it would be expected. However, in the graphic the model developed in this article (KL T63i Pen) have presented the fastest convergence in this simulation. Converging to the more accurate result with fewer elements is a desirable important feature for a finite element model.

Cantilever Beam
A Cantilever Beam is subject to a transverse force at the free end. The simulation has been executed for different mesh discretization and for different thicknesses in order to analyze the FEM model built in this article. The results are compared to the Lagrange shell FEM model from Costa e Silva (2020) and to commercial software (Nastran Triangular Shell Element -CTRIA6).
Mechanical properties and dimensions of the beam that have been used are displayed below.
• E = 210 GPa; • = 0,3; • L (beam length) =2400 mm; • h(height) = 100 mm; • b(width of cross section) = 11,64 mm; The beam is subject to vertical down force up to = 10000 N and an horizontal perturbation force is also introduced ( = ⋅ 10 −4 ) in order to induce the buckling phenomena.  The analytical results (Linear elastic) for this problem may be found in Solid mechanics books and is illustrated by the formula where • m is the maximum vertical displacement at the center of the free side of the beam.
• is the concentrated force at the center of the free side of the beam.
• is the length of the beam • E is the Young modulus of the material.
• is the moment of inertia of the rectangular cross section • b is the thickness of the cross section.
• h is the height of the cross section.   The results shown above demonstrates that both Lagrange and penalty methods led to considerable similar (difference smaller than 5% from one another) and consistent results, once they could capture the flexural buckling of the beam. Both methods have close results to the analytical model only for low load factors due to the Linear elastic limitation of the analytical model. For high loads, the beam buckles and the analytical results loses its physical meaning. For this simulation, the commercial software (Nastran Triangular Shell Element -CTRIA6) was unable to converge for load factor bigger than 0,24, presenting numerical errors and warnings. This load factor is approximately the step that the beam starts to buckle. The simulations were accomplished with approximately 160 elements and in Nastran, it has been used the non-linear material setup with the same mechanical properties as described above. The capacity to successfully capture the buckling phenomena is an important and desirable feature that was observed in this simulation.

Pinched Cylinder
This example is reproduced in other articles (Ivannikov et al., 2015) (Campello et al, 2003) (Sansour and Kollmann, 2000) (Costa e Silva, 2020) and is a good benchmark to evaluate the shell model for big displacements and curvatures. It consists of a cylinder with rigid ends and subject to a radial vertical force in its center. The Figure 13   In the simulations, it has been used the same parameters of the bibliography, then one can compare results. The mechanical properties, dimensions and loads, that are used in the simulation, are • E = 3 ⋅ 10 10 Pa • R(Cylinder Radius) = 100 mm; • b(thickness) =1 mm; Figure 15 illustrates the vertical displacement of point A and Figure 16 plots the horizontal displacement of point B for different load factors. Each line in the graphics represents a different Finite element model. Results from some sources in bibliography have been plotted in the same Graph for comparison.  Figure 15 and Figure 16 illustrates that all models have similar results for low loading, however for bigger forces, and big displacements, Nastran CTRIA6 has different and unstable results. Yet in the same figures, the model developed in this article has demonstrated very close results to the lagrange based model (Costa e Silva, 2020) with results differing less than 10% from one another even for higher loads. Nastran CTRIA6 presented non-smooth results regarding the displacements and presented some warnings in the numerical simulation regarding numerical convergence quality. In the simulations, it has been used about 5000 shell elements and in Nastran it has been used non-linear material setup with the same mechanical properties as specified above.

CONCLUSIONS
The present article develops analytically a very simple finite element and presents results in 3 different simulations comparing it to other FEM models and to linear analytical models. All the kinematics is based on Kirshoff-Love hyphothesis and the C1 continuity expected between adjacent elements because of the K-L hypothesis is expect to be fulfilled by applying penalty method internally to the element enforcing agreement between a rotation degree of freedom (that is shared between adjacent elements) and the displacement field (defined by others DoFs). One of the advantages of the current model compared to others K-L based elements is that a second mesh for enforcing penalty or lagrange multiplyers is not necessary. All the approximation is accomplished internally to the element.
The promising results obtained so far are coherent and satisfactory, but more tests are going to be made in further studies to have a stronger idea of the positives/negatives aspects of the model. An important conclusion of these first simulations is that all the simulations presented in this article shows compatible and reliable results regarding the model developed. Some differences, however, are observed. The results obtained with commercial software (Nastran CTRIA6) were observed as more stiffer than the K-L models (Penalty and Lagrange based models) in the first and third example (Square plate and Pinched cylinder) and it is also unstable for large displacements (Second and third Examples). It is also possible to check by Figure 8 that the K-L Penalty based model has shown better convergence with fewer elements. It is important to remark as well that the analitical models in all examples could not depict the non-linear phenomena that is captured with the non-linear Finite element and consequently presents similar results only for small displacements due to its linear characteristic.
For the next studies, the authors plan also to develop a shell element with out the use of penalty or Lagrange multipliers as a mechanism to guarantee C1 continuity. The theory behind this new element is still under developement and more information is going to be presented in forthcoming work. It is believed that a very reliable and consistent thiangular shell finite element as presented in this article combined with powerfull mesh generators is a great option for simulations of thin complex structures.