# A fully incremental simple triangular multilayer Kirchhoff-Love shell element

Gustavo Canário Gomes, Paulo de Mattos Pimenta, Matheus Lucci Sanchez, Adnan Ibrahimbegovic

► To cite this version:

Gustavo Canário Gomes, Paulo de Mattos Pimenta, Matheus Lucci Sanchez, Adnan Ibrahimbegovic. A fully incremental simple triangular multilayer Kirchhoff-Love shell element. Latin American Journal of Solids and Structures, 2024, 21 (9), 10.1590/1679-78258159. hal-05384306

HAL Id: hal-05384306

https://hal.science/hal-05384306v1

Submitted on 26 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A simple triangular multilayer Kirchhoff-Love shell element with shear stress analysis

Gustavo Canário Gomesª* https://orcid.org/0000-0002-0859-3440, Paulo de Mattos Pimentaᵇ
https://orcid.org/0000-0002-7808-8319, Adnan Ibrahimbegovic, https://orcid.org/0000-0002-6502-0198ᶜ

## Abstract

This paper presents a new triangular multi-layer nonlinear shell finite element, suitable for large displacements and rotations, and a formulation to obtain through the thickness shear stress considering both geometrical and material nonlinearity. This is a nonconforming element with 10 nodes, cubic displacement interpolation and enforcement of the rotation field based on Rodrigues rotation parameters and lagrangean parameters in the 6 side nodes, with a total of 42 degrees of freedom. Associated with the new element, another novelty of this work is the development of a multilayer kinematical model with properties from Kirchhoff-Love theory, approximating the shell director across layers as constant and the formulation to obtain through the thickness shear stress considering both geometrical and material nonlinearity. The model is numerically implemented, and results are compared to different references in multiple examples, showing the capabilities of the formulation. Although lagrangean parameters are used to enforce $C^1$ continuity, it is believed that due to it being possibly the simplest multilayer extension, combined with fully incremental DoFs, simple kinematic, a relatively small number of DoFs, possibility to use various 3D material models, easily connected with multiple branched shells and beams, and geometric exact theory, this is a simple yet powerful shell element. The shear stress formulation, capable of being used with nonlinear materials, presents itself as a novelty, deviating from the standart use of linear material behaviors with co-rotational formulations, with good results.

## Keywords

Triangular Shell Element, Multilayer shell, Shear stress, Nonlinear Shell Formulation, Kirchhoff-Love shell, Stress analysis, Large Strain.

## 1 INTRODUCTION

The finite element method is, until this day, an important numerical method for solving partial differential equations with many special applications in solid mechanics and engineering. Among the possible elements, shells still are widely used in structural engineering due to their relatively simple models (in terms of degrees of freedom (DoF) quantities) capable of simulating complex strucutres with thin aspect ratios and very few elements.

In the present work, a rigid shear model based on the Kirchhoff-Love shell theory is considered. In a succinct manner, the theory does not exhibit shear strain due to three fundamental kinematic hypotheses (see (REDDY, 2006)): (1) Straight lines perpendicular to the mid-surface (i.e., transverse normal) before deformation remain straight after deformation; (2) Inextensible transverse normal; (3) The transverse normals rotate such that they remain perpendicular to the middle surface after deformation. As presented by (SANCHEZ, PIMENTA and IBRAHIMBEGOVIC, 2023) and (Gomes et al., 2024), the advantages of using Kirchhoff-Love shells when compared with other possibilities are the better numerical convergence in very thin shells and the fewer DoF. The theory requires a $C^1$-continuity in the interface between

<footnote>ª Polytechnic School at University of São Paulo, Av. Prof. Almeida Prado, trav. 2 n° 83 Edifício Paula Souza, 05424-970, São Paulo, Brazil.
gustavocanario@usp.br</footnote>

<footnote>ᵇ Polytechnic School at University of São Paulo, Av. Prof. Almeida Prado, trav. 2 n° 83 Edifício Paula Souza, 05424-970, São Paulo, Brazil.
ppimenta@usp.br</footnote>

<footnote>ᶜ Laboratoire Roberval de Mécanique, Université de Technologie de Compiègne/Alliance Sorbonne Université, Rue du docteur Schweitzer
CS 60319,60203 Compiègne, France. adnan.ibrahimbegovic@utc.fr</footnote>

<footnote>* Corresponding author</footnote>
adjacent elements, which may be a challenge. Among the possibilities to deal with such ordeal, there is the well stab-
lished TUBA family elements. initially developed for plates by (ARGYRIS, 1968), expanded for shells by (ARGYRIS, 1972)
and recently revisited by (IVANNIKOV, 2015). Yet another possibilities are the one presented by (CIRAK F, 2001), which
is the use of subdivision surfaces, with progressive division of the mesh until the C¹ continuity is achieved and another
being the use of a rotation scalar, located at the mid-points, associated with a rotation in the direction of the sides of
the shell, as presented in (SANCHEZ, PIMENTA and IBRAHIMBEGOVIC, 2023). A different approach, proposed by
(Neunteufel and Schörbel, 2019) with the development of a mixed shell finite element with the Hellan–Herrmann–
Johnson method and latter improved by (Neunteufel and Schörbel, 2021) with Regge interpolation and tying points is
also possible. For the meshless case, as seen in (IVANNIKOV, 2014), the requirement can be achieved with the multiple
fixed least squares (MFLS) approximation functions and their generalized version, considering the use of proper
weights and base functions. The presented options are just a few of the possibilities, as the research field for Kirchhoff-
Love plates and shells has been active in the past two decades, with some additional works being by (PIMENTA, NETO
and CAMPELLO, 2010), (RABCZUK T, 2007), (KIENDL, 2009), (Bohinc, 2014), (SILVA et al., 2019) and (Wu, 2024).

In this work we present an extension of the triangular Kirchhoff-Love shell element developed in (VIEBAHN, PIMENTA and SCHRÖDER, 2017) and (Sanchez, Costa e Silva and Pimenta, 2020) to the multi-layer consideration, creating a new element, capable of representing composite material behavior. In the reference work, in order to guarantee C1 continuity, a rotation continuity is imposed with the use of penalty parameters at the midside node and with mesh refinement. Here, we use lagrangean multipliers in order to obtain the same effect.

In the current work, a consideration that greatly simplifies the multi-layer kinematics and eliminates the need of
multiple independent variables along the thickness is used, this being based on the thin nature of Kirchhoff-Love shells.
It is assumed that rotation along the layers remains constant, regardless of possible change of material and thickness.
The theory defines the energy of the system by conjugating the first Piola-Kirchhoff stress tensor and the deformation
gradient, allowing for contribution of each layer in the internal energy of the shell.

We assume a plane reference configuration for the mid-surface, and the rotation is considered with the Rodrigues formula, in a pure lagrangian way, considering Rodrigues parameter and with the use of the shell displacement. In the present work, we have a nonconforming element with 10 nodes, cubic displacement and enforcement of the rotation field based on Rodrigues rotation parameters and lagrangean parameter in the 6 points along the sides, giving in total 42 degrees of freedom. The numerical implementation is done using AceGen and AceFEM, fully utilizing the automation of the packages as a mean to obtain physical quantities such as the residual force vector and the stiffness matrix.

Associated with the element, we present an equilibrium based formulation to obtain the through the thickness
shear stresses with the use of a simplified Crank-Nicholson scheme, once again taking advantage of the thin nature of
the Kirchhoff-love shell. This formulation fully considers nonlinearities, allowing geometrical nonlinearity and the use of
nonlinear material formulations, such as the neo-hookean hyperelastic material model. The idealization for the current
element came from the necessity of a third degree polynomial interpolation of the displacement field to properly rep-
resent the shear stress field with the equilibrium based stress formulation here presented. The need of such interpola-
tion degree was identified when the stress formulation was used in the element presented in (Gomes et al., 2024),
which uses a quadratic displacement field with a linear rotation field.

In regard to the notation considered throughout the text, italic Latin or Greek lowercase letters (a, b, ..., α, β) de-
note scalars quantities, bold italic Latin or Greek lowercase letters (a, b, ..., α, β) denote vectors, bold italic Latin or
Greek capital letters (A, B, ...) denote matrix and second order tensors. Summation convention over repeated indices,
with values {1, 2, 3} for latin letters and {1, 2} for greek letters

2 SHELL FORMULATION

2.1 Multi-Layer consideration

The main hypothesis that guides the current multi-layer theory is that, due to the thin nature of the Kirchhoff-Love
shell, the shell director does not change direction along each layer, see Figure 1. This greatly simplifies the multi-layer
kinematics, as there is no need of specific inter-layer displacement equations to deal with the shell director variation
along the thickness, as usually is the case for Reissner-Mindlin models.

In the model, we assume the middle plane of the shell as the reference to obtain the kinematics equations and to
apply the FEM mesh and external loads.
<figure><img src="image_2.png" /><figcaption>Figure 1</figcaption></figure>

Figure 1 - a) Shell director along layers; b) Shell description and kinematics

In order to obtain the position along the shell thickness, we initially establish a general mid-surface for all the layers, by taking a direct average of the sum of the layer's height

$$ \zeta_{Med} = \frac{\sum_{i=1}^{i=n_{layers}} \zeta_i}{n_{layers}}, \qquad (1) $$

For the current position along the thickness of the shell, due to the simplification of a constant director, we use a simple procedure, described by

$$ \zeta_{posiGen} = \zeta_{LayerAcum} - \zeta_{Med} + \zeta_{PosiLayerCurr}, \qquad (2) $$

In the equation above, we assume a starting point from the bottom to the top of distribution of the plies. As such, $\zeta_{LayerAcum}$ is the accumulated lengths of the previous layers, $\zeta_{Med}$ is the reference mid-plane for all the plies and $\zeta_{layerCurr}$ is the current layer position. When working with numerical integration, in which the domain is discretized from $-1 \le \zeta_{GP} \le 1$, we have that

$$ \zeta_{PosiLayerCurr} = \frac{(\zeta_{GP} + 1) \zeta_{LayerCurr}}{2}, \qquad (3) $$

with $h_{LayerCurr}$ being the current layer length.

The director vectors $\mathbf{a}$ and $\mathbf{a}^r$ at current and reference configurations, for any point in the length of the shell are then given by

$$ \mathbf{a} = \zeta_{posiGen} e_3 \quad \text{and} \quad \mathbf{a}^r = \zeta_{posiGen} e_3^r. \qquad (4) $$

## 2.2 Kinematics

It is assumed that the middle plane is in the reference configuration of the shell and it is flat, see Figure 1, with an orthogonal coordinate system {$e_1^r, e_2^r, e_3^r$} where $e_1^r$ and $e_2^r$ are in the middle plane of the shell and $e_3^r$ is normal to the plane. In the current configuration, the orthogonal coordinate system is {$e_1, e_2, e_3$}, with the relation $e_i = Qe_i^r$ valid, where $Q$ is the rotation tensor.

The definition of any point in the reference configuration will then be given by

$$ \xi = \zeta + \mathbf{a}^r, \zeta = \xi_{\alpha} \mathbf{e}_{\alpha}^r \in \Omega^r \quad \text{and} \quad \mathbf{a}^r = \zeta \mathbf{e}_3^r \in H^r, \qquad (5) $$

where $\zeta$ are the coordinates of a material point in the shell middle plane and $\mathbf{a}^r$ is the fiber direction of the shell at this point. The thickness $h$ of the shell in the reference configuration is considered through the parameter $\zeta \in H = [-h^b, h^t]$, where $h = h^b + h^t$.

In a similar manner to the definition of a material point in the reference configuration, the definition of any point in the current is given by
$$x = z + a, z = \zeta + u \quad \text{and} \quad a = Qa^r, \qquad (6)$$

$z$ represents the current position of a point in the middle surface and $a$ the position along the thickness of the shell, in $e_3$ direction.

Adopting the notation $\partial(\cdot)/\partial(\xi_\alpha) = (\cdot)_{,\alpha}$, the deformation gradient can be obtained from:

$$F = x_{,\alpha} \otimes e_{\alpha}^r + \frac{\partial x}{\partial \zeta} \otimes e_3^r = f_a \otimes e_{\alpha}^r + f_3 \otimes e_3^r = (z_{,\alpha} + \kappa_{\alpha} \times a) \otimes e_{\alpha}^r + Q e_3^r \otimes e_3^r. \qquad (7)$$

where

$$K_{\alpha} = Q_{,\alpha}Q^{T} \text{ and } \kappa_{\alpha} = \text{axial}(K_{\alpha}). \qquad (8)$$

The back-rotated deformation gradient and the strains are defined as

$$F^r = Q^T F = I + \gamma^r \otimes e^r_{\alpha}, \qquad (9)$$

with

$$\gamma_{\alpha}^{r} = \eta_{\alpha}^{r} + \kappa_{\alpha}^{r} \times a^{r}; \quad \eta_{\alpha}^{r} = Q^{T}z_{,\alpha} - e_{\alpha}^{r}; \quad \kappa_{\alpha}^{r} = Q^{T}\Gamma_{\beta}u_{\beta\alpha}. \qquad (10)$$

$$\Gamma_1 = (e_1 \cdot z_{,1})^{-1} [\text{Skew}(e_1) - (e_1 \cdot z_{,2})(e_2 \cdot z_{,2})^{-1}(e_1 \otimes e_3)]. \qquad (11)$$

$$\Gamma_2 = (e_2 \cdot z_{,2})^{-1}(e_1 \otimes e_3). \qquad (12)$$

In the previous equations, the vectors $\eta_{\alpha}^r$ and $\kappa_{\alpha}^r$ are the cross-sectional generalized strains.

It is important to remark that, due to the Kirchhoff-Love assumption, $\gamma_{\alpha}^r \cdot e_3^r = \eta_{\alpha}^r \cdot e_3^r = 0$.

As we are working with a Kirchhoff-Love type shell, the director vector $a$ remains orthogonal to the middle surface of the shell. As such, the rotation tensor $Q$ can be obtained with the use of the orthogonal coordinate systems as

$$Q = e_i \otimes e_i^r, \qquad (13)$$

with the current orthogonal basis obtained in relation to the geometry of the deformed configuration, as

$$e_1 = \|z_{,1}\|^{-1} z_{,1}, \quad e_3 = \|z_{,1} \times z_{,2}\|^{-1} (z_{,1} \times z_{,2}), \quad e_2 = e_3 \times e_1 \qquad (14)$$

## 2.3 Variational formulation and material equations

The stiffness matrix of the current work is obtained through the variational formulation and with the help of AceGen. The usage of this software allows for multiple material models readily implemented with automatic differentiation, as long as the material model is given with a strain energy function. In the present work, in order to assure the convergence of the simulations, we work here with two polyconvex material models suitable for large displacements and rotation, this being in alignment with (Brank, Korelc and Ibrahimbegović, 2002).

As the materials considered are hyperelastic, there exists a total strain energy function $\psi$ that describes the elastic energy stored in the body. Here, we define the functional

$$\prod_{B} (u) = \int_{B} (\psi(C(u)) - \rho_{0}\bar{b} \cdot u) dV - \int_{\partial B} (\bar{t} \cdot u) dA. \qquad (15)$$

The directional derivative of the functional is obtained as

$$\delta \prod_{B} (u) = D \prod_{B} (u) \cdot v = \frac{\partial}{\partial \alpha} \prod_{B} (u + \alpha v)_{\alpha=0} = 0. \qquad (16)$$

One can define the internal total energy for the shell as

$$\prod_{\Omega} (u) = U^{\text{int}} = \int_{\Omega} \int_{H^r} \psi(\xi_1^r, \xi_2^r) dH d\Omega = \int_{\Omega} \tilde{\psi}(\xi_1^r, \xi_2^r) d\Omega. \qquad (17)$$

The strain energy is, obviously, material dependent, with three options presented in the current work. The internal stress and material tangent moduli are obtained by
$$ P = \frac{\partial \psi}{\partial F} \quad \text{and} \quad D_{\alpha\beta}^{r} = \frac{\partial \sigma_{\alpha}^{r}}{\partial \varepsilon_{\alpha}^{r}} = \frac{\partial^2 \bar{\psi}}{\partial \varepsilon_{\alpha}^{r} \partial \varepsilon_{\beta}^{r}}. \qquad (18) $$

The first material described is the Ciarlet-Simo Neo-Hookean isotropic material, defined by the polyconvex strain energy function

$$ \psi = \frac{1}{2} \lambda \left( \frac{1}{2}(J^2 - 1) - \ln(J) \right) + \frac{1}{2} \mu (I_1 - 3 - 2 \ln(J)), \qquad (19) $$

with the invariants, Cauchy-Green stress tensor and Lamé coefficients given as

$$ I_1 = \operatorname{tr} C = f_i \cdot f_i, \quad I_2 = \operatorname{tr}(C \circ f C) = g_i \cdot g_i \quad \text{and} \quad I_3 = \det C = J^2 = (f_1 \cdot (f_2 \times f_3))^2, \qquad (20) $$

$$ C = F^T F, \qquad (21) $$

$$ \lambda = \frac{Ev}{(1 + \nu)(1 - 2\nu)} \quad \text{and} \quad \mu = \frac{E}{2(1 + \nu)}. \qquad (22) $$

The second material is anisotropic, polyconvex in both the isotropic and anisotropic contributions, as presented by (VIEBAHN, PIMENTA and SCHRÖDER, 2017), (SCHRÖDER and NEFF, 2010) and (SCHRÖDER et al., 2016) and is given as

$$ \psi(C) = \psi^{i-p}(C) + \psi^{a-p}(C) = \psi^{i-p}(C) + \sum_{i}^{2} \alpha_{1}^{(i)} \langle I_{4}^{(i)} - 1 \rangle_{2}^{\alpha_{2}^{(i)}}, \qquad (23) $$

where $\alpha_1^{(i)}$ and $\alpha_2^{(i)}$ are variables in relation to material properties in the $m^{(i)}$ direction, "⟨⟩" are Macaulay brackets and $I_4^{(i)}$ is an invariant defined as

$$ I_4^{(i)} = \operatorname{tr}[A_{\alpha\beta}^{(i)} f_{\alpha} \cdot f_{\beta}], \qquad (24) $$

with

$$ A_{\alpha\beta}^{(i)} = (e_{\alpha} \cdot m^{(i)})(e_{\beta} \cdot m^{(i)}) = e_{\alpha} \cdot M^{(i)} e_{\beta} \quad \text{and} \quad M^{(i)} = m^{(i)} \otimes m^{(i)}. \qquad (25) $$

In relation to the isotropic part, we use the neo-Hookean material model previously presented.

The material model presented can also be simplified in order to obtain a transversely isotropic material, this being described by

$$ \psi(C) = \psi^{i-p}(C) + \psi^{a-p}(C) = \psi^{i-p}(C) + \alpha_{1}^{(1)} \langle I_{4}^{(1)} - 1 \rangle_{2}^{\alpha_{2}^{(1)}}. \qquad (26) $$

## 2.4 Plane stress

The plane stress condition is a necessity when the shell theory does not consider thickness variation. As such, the present work enforces this condition analytically for the presented materials by adding the term $\gamma_{33}^r = \gamma_{33} e_3^r$ to eq. (5), rewriting it as

$$ F^r = I + \gamma_{\alpha}^r \otimes e_{\alpha}^r + \gamma_{33}^r \otimes e_3^r. \qquad (27) $$

The plane stress is then enforced by $\tau_{33} = \tau_3^r \cdot e_3^r = 0$.

With eq. (11), eq. (12) and eq. (19), one arrives at

$$ \tau_{33} = 2(1 + \gamma_{33}) \frac{\partial \psi}{\partial I_1} + \frac{\partial \psi}{\partial J} \bar{J}, \qquad (28) $$

where $\bar{J} = e_3^r \cdot (f_1^r \times f_2^r)$. After some manipulations, one arrives at

$$ \gamma_{33} = \sqrt{\frac{\lambda + 2\mu}{\lambda \bar{J}^2 + 2\mu}} - 1. \qquad (29) $$
In the case that the process of determining the value of $\gamma_{33}$ is not as direct or simple, one can use the newton method algorithm described in (CAMPELLO, PIMENTA and WRIGGERS, 2003) in order to ensure the plane stress condition.

## 3 THE FINITE ELEMENT

The element created is an extrapolation of the element developed by (Sanchez, Costa e Silva and Pimenta, 2020). The new element, entitled T10LKLM is a 16-noded displacement based triangular element (see Figure 2a). The discretization of the element coordinates and displacement field is done with area coordinates, cubic interpolation and, additional rotation parameters ($\phi_{\Delta}$) and lagrangean parameters, located at the sides of the element as presented in Figure 2b, used in order to guarantee the C1 continuity with mesh refinement. As such, this element presents 42 DoF.

<figure><img src="image_3.png" /><figcaption>Figure 2 - T10-6iKLM</figcaption></figure>

In the present element, we use the simplification of a constant shell director through the layers of the element due to the small thickness.

The element uses 7 integration points at the shell plane (located at the edge nodes, bubble node and at the mid-sides) and, when integrating along the thickness, each layer used 3 Lobatto integration points. Lobatto points are used in order to have integration points located at each layer limits. Obviously, by integrating each layer separately, one allows for use of different materials and thickness at each layer. The use of 3 integration points along the layer is also compeled by the format of the finite difference method used for the shear stresses obtained along the thickness direction. For the method described in a posterior section, one uses points at the top and bottom layer boundary and 1 point at the middle of the layer.

We here give a brief presentation of implementation details. The AceGen and AceFem computer codes from the software Wolfram Mathematica were used as the implementation environment. AceGen allows for the implementation of the element itself, with a symbolic type language and a vast quantity of pre-programmed capabilities, such as multiple Gauss points distributions, elements geometry and multiple commands already present in Mathematica readily available. Another remarkable quality of AceGen, fully used throughout the implementation of the element here presented, is the automated differentiation. The main use of such feature in the current research was to obtain the element residual vector and its stiffness matrix directly from the differentiation of the energy polyconvex strain energy function in relation to the DoFs. Regarding AceFem, this is another package that also possesses multiple pre-programed routines which greatly facilitate the geometry construction, meshing, boundary conditions imposition, analysis and post-processing. As some of its remarkable characteristics are fully programed Newton and arc-length methods, associated with multiple
possible linear system solution methods and, high quality images and animations obtained at post-processing (see section 6 for examples, as the meshes, deformed configurations and graphs were generated within AceFem).

## 4 KINKING ANGLE ENFORCEMENT

The kinking angle enforcement here presented follows the methodology of (Sanchez, Costa e Silva and Pimenta, 2020) and (VIEBAHN, PIMENTA and SCHRÖDER, 2017), with a different approach at the end when obtaining the energy contribution as the original work uses a penalty parameter and here, we use a lagrangean multiplier.

The shell kinematics here adopted is based on a Kirchhoff-Love model, with the necessity of continuity of the displacement gradient. As such, the finite element must guarantee $C^1$ continuity in the whole domain, internal domain and edges of the element. Different from the reference work, here, the kinking angle between adjacent elements is constrained with the help of lagrangean multipliers, based on the Rodrigues rotation parameters developed by (SILVA, 2020).

The enforcement of $C^1$ continuity is asymptotically satisfied as $h \to 0$ (mesh refinement) as long as the kinking angle (at the edges) between adjacent elements does not change during deformation. The kinking angle may be understood as the angle between mid-surface normal vectors of adjacent elements.

The rotation tensor $Q$ represents the rotation of the normal vector of the shell from reference to current coordinate system. As such, if one imposes that the rotation tensor (or its incremental part) is the same for both adjacent elements, there is continuity of the kinking angle. In order to achieve that, the authors follow the steps presented by (Sanchez, Costa e Silva and Pimenta, 2020).

We begin using Rodrigues rotation parameters and with the description of the rotation tensor as

$$ \hat{Q}(\alpha) = \left(I - \frac{1}{2}A\right)^{-1} \left(I + \frac{1}{2}A\right) = I + \frac{4}{4 + \alpha^2} \left(A + \frac{1}{2}A^2\right), \qquad (30) $$

where

$$ \alpha = 2 \tan(\theta/2), \quad \alpha = \alpha_e, \quad A = \text{Skew}(\alpha). \qquad (31) $$

With such description, the incremental description of rotation is obtained as

$$ Q_{i+1} = Q_{\Delta}Q_{i}, \quad Q_{i+1} = \hat{Q}(\alpha_{i+1}), \quad Q_{\Delta} = \hat{Q}(\alpha_{\Delta}), \quad Q_{i} = \hat{Q}(\alpha_{i}) \qquad (32) $$

where

$$ \alpha_{i+1} = \frac{4}{4 + \alpha_i \cdot \alpha_{\Delta}} \left( \alpha_i + \alpha_{\Delta} - \frac{1}{2} \alpha_i \times \alpha_{\Delta} \right). \qquad (33) $$

In the previous equations, $\alpha_i$, $\alpha_{\Delta}$ and $\alpha_{i+1}$ are the rotation from the previous, increment and next step, respectively. The same logic follow for the rotation tensor $Q_i$, $Q_{\Delta}$, $Q_{i+1}$.

With such description, one can write the coordinate system in sequential steps $e_1^i, e_2^i, e_3^i$ in step $i$ and $e_1^{i+1}, e_2^{i+1}, e_3^{i+1}$ in step $i+1$ as

$$ e_3^{i+1} = Q_{\Delta}e_3^i, \quad e_3^{i+1} - e_3^i = \alpha_{\Delta} \times e_3^m \quad \text{and} \quad e_3^m = \frac{1}{2}(e_3^{i+1} + e_3^i) \qquad (34) $$

Based on this description, two equations were presented originally in (SILVA, 2020) and used in (Sanchez, Costa e Silva and Pimenta, 2020), those being

$$ \alpha_{\Delta} = \|e_3^m\|^{-2}(e_3^i \times e_3^{i+1}) + \varphi_{\Delta} \|e_3^m\|^{-1} \qquad (35) $$

$$ \varphi_{\Delta} = \|e_3^m\|^{-1} \alpha_{\Delta} \cdot e_3^m \qquad (36) $$

The equations above show that it is possible to determine the normal vector $e_3$ at the edge of the element at step $i+1$ if the parameter $\varphi_{\Delta}$ is known. With this, we have two different ways to determine $e_3$, one given by the geometry (see (14)) and one given by rotation parameters, presented in this section. These approach, however, have no guarantee to be equal. In order to do so, a lagrangean multiplier approach is used to ensure that the displacement field is compatible with $\varphi_{\Delta}$ independent parameters, as a mechanism to ensure $C^1$ continuity with sufficient mesh refinement. In order to do so, one uses (14) and the expressions for the derivatives of $u$ in relation to directions 1 and 2 (obtained automatically for the cubic interpolation with the use of AceGEN).
With this, one may write the incremental rotation vector based on the displacement field as

$$ \alpha_{\Delta} = \frac{2(e_j^i \otimes e_j^{i+1})}{1 + e_k^i \cdot e_k^{i+1}} \qquad (37) $$

and, consequently, have a rotation parameter

$$ \phi_{\Delta} = \alpha_{\Delta} \cdot \frac{e_1^m}{\|e_1^m\|}, \quad \text{with} \quad e_1^m = \frac{1}{2}(e_1^i + e_1^{i+1}) \qquad (38) $$

With this, one can add the contribution of the terms needed to make the rotation parameters $\phi_{\Delta}$ and $\phi_{\Delta}$ equal to the Helmholtz free energy, obtaining

$$ \Pi_{Lag} = \psi + \lambda_4(\phi_{\Delta}^{(4)} - \phi_{\Delta}^{(4)}) + \lambda_5(\phi_{\Delta}^{(5)} - \phi_{\Delta}^{(5)}) + \lambda_6(\phi_{\Delta}^{(6)} - \phi_{\Delta}^{(6)}) + \lambda_7(\phi_{\Delta}^{(7)} - \phi_{\Delta}^{(7)}) + \lambda_8(\phi_{\Delta}^{(8)} - \phi_{\Delta}^{(8)}) \\ + \lambda_9(\phi_{\Delta}^{(9)} - \phi_{\Delta}^{(9)}) \qquad (39) $$

## 5 SHEAR STRESS ALONG LAYERS

We begin the shear stress analysis with the local force balance given by

$$ \mathrm{div}P = -\bar{b}. \qquad (40) $$

As $P$ is the first Piola-kirchhoff stress tensor, one can write

$$ \tau_{i,i} = (Q\tau_i^r),_i = Q_i\tau_i^r + Q\tau_{i,i}^r = -\bar{b}, \qquad (41) $$

with the relation

$$ Q_i Q^T Q \tau_i^r = Q \kappa_\alpha^r \times \tau_\alpha^r, \qquad (42) $$

one can write

$$ \tau_{i,i} = Q(\tau_{i,i}^r + \kappa_\alpha^r \times \tau_\alpha^r) = -\bar{b} \qquad (43) $$

We can further expand the expression (43) as

$$ \tau_{i,i}^r + \kappa_\alpha^r \times \tau_\alpha^r = \tau_{\alpha,\alpha}^r + \kappa_\alpha^r \times \tau_\alpha^r + \tau_{3,3}^r = -Q^T \bar{b} = -\bar{b}^r. \qquad (44) $$

Once we apply the plane stress condition, it renders $\tau_3^r = \tau_{3,\alpha}e_\alpha^r = \tau_{31}e_1^r + \tau_{32}e_2^r$. On the other hand, we have

$$ \tau_1^r = \tau_{11}e_1^r + \tau_{12}e_2^r + \tau_{13}e_3^r, \qquad (45) $$

$$ \tau_2^r = \tau_{21}e_1^r + \tau_{22}e_2^r + \tau_{23}e_3^r, \qquad (46) $$

and

$$ \kappa_1^r = \kappa_{11}e_1^r + \kappa_{12}e_2^r + \kappa_{13}e_3^r, \qquad (47) $$

$$ \kappa_2^r = \kappa_{21}e_1^r + \kappa_{22}e_2^r + \kappa_{23}e_3^r. \qquad (48) $$

Hence, considering the plane stress condition and substituting (45) to (48) in (44) renders

$$ (\tau_{11}e_1^r + \tau_{12}e_2^r + \tau_{13}e_3^r),_1 + (\tau_{21}e_1^r + \tau_{22}e_2^r + \tau_{23}e_3^r),_2 + (\kappa_{11}e_1^r + \kappa_{12}e_2^r + \kappa_{13}e_3^r) \times (\tau_{11}e_1^r + \tau_{12}e_2^r + \tau_{13}e_3^r) \\ + (\kappa_{21}e_1^r + \kappa_{22}e_2^r + \kappa_{23}e_3^r) \times (\tau_{21}e_1^r + \tau_{22}e_2^r + \tau_{23}e_3^r) + \tau_{31,3}e_1^r + \tau_{32,3}e_2^r = -\bar{b}^r. \qquad (49) $$

Applying the derivatives, and after some manipulations, one finds the equations

$$ \tau_{11,1} + \tau_{21,2} + (\kappa_{12}\tau_{13} - \kappa_{13}\tau_{12}) + (\kappa_{22}\tau_{23} - \kappa_{23}\tau_{22}) + \tau_{31,3} = -\bar{b}_1, \qquad (50) $$

$$ \tau_{12,1} + \tau_{22,2} + (\kappa_{13}\tau_{11} - \kappa_{11}\tau_{13}) + (\kappa_{23}\tau_{21} - \kappa_{21}\tau_{23}) + \tau_{32,3} = -\bar{b}_2. \qquad (51) $$

Let us briefly focus on the local balance with equation
$$ \text{Skew}(PF^T) = \frac{1}{2} (PF^T - (PF^T)^T) = 0. \qquad (52) $$

As the skew must be zero, we have

$$ (PF^T) = (PF^T)^T. \qquad (53) $$

By focusing on the left term,

$$ (PF^T) = Q(\tau_1^r \otimes e_1^r + \tau_2^r \otimes e_2^r + \tau_3^r \otimes e_3^r)F^{rT}Q^T. \qquad (54) $$

With the rotation tensor property, we have

$$ (PF^T) = (\tau_1^r \otimes e_1^r + \tau_2^r \otimes e_2^r + \tau_3^r \otimes e_3^r)F^{rT}. \qquad (55) $$

Substituting the back-rotated transformation tensor and transposing, we have

$$ (PF^T) = (\tau_1^r \otimes e_1^r + \tau_2^r \otimes e_2^r + \tau_3^r \otimes e_3^r)[e_1^r \otimes (e_1^r + \gamma_1^r) + e_2^r \otimes (\gamma_2^r + e_2^r) + e_3^r \otimes e_3^r]. \qquad (56) $$

After some manipulations and considering that

$$ \gamma_1^r = \eta_1^r + \zeta \kappa_1^r \times e_3^r, \qquad (57) $$

$$ \gamma_2^r = \eta_2^r + \zeta \kappa_2^r \times e_3^r, \qquad (58) $$

one finds the components

$$ (PF^T)_{11} = [\tau_{11}(1 + \eta_{11} + \zeta\kappa_{12}) + \tau_{21}(\eta_{21} + \zeta\kappa_{22})]e_1^r \otimes e_1^r, \qquad (59) $$

$$ (PF^T)_{12} = [\tau_{11}(\eta_{12} - \zeta\kappa_{11}) + \tau_{21}(1 + \eta_{22} - \zeta\kappa_{21})]e_1^r \otimes e_2^r, \qquad (60) $$

$$ (PF^T)_{13} = \tau_{31}e_1^r \otimes e_3^r, \qquad (61) $$

$$ (PF^T)_{21} = [\tau_{12}(1 + \eta_{11} + \zeta\kappa_{12}) + \tau_{22}(\eta_{21} + \zeta\kappa_{22})]e_2^r \otimes e_1^r, \qquad (62) $$

$$ (PF^T)_{22} = [\tau_{12}(\eta_{12} - \zeta\kappa_{11}) + \tau_{22}(1 + \eta_{22} - \zeta\kappa_{21})]e_2^r \otimes e_2^r, \qquad (63) $$

$$ (PF^T)_{23} = \tau_{32}e_2^r \otimes e_3^r, \qquad (64) $$

$$ (PF^T)_{31} = [\tau_{13}(1 + \eta_{11} + \zeta\kappa_{12}) + \tau_{23}(\eta_{21} + \zeta\kappa_{22})]e_3^r \otimes e_1^r, \qquad (65) $$

$$ (PF^T)_{32} = [\tau_{13}(\eta_{12} - \zeta\kappa_{11}) + \tau_{23}(1 + \eta_{22} - \zeta\kappa_{21})]e_3^r \otimes e_2^r, \qquad (66) $$

Focusing in (65) and (66), and due to the symmetry described in (53),

$$ \tau_{13}(1 + \eta_{11} + \zeta\kappa_{12}) + \tau_{23}(\eta_{21} + \zeta\kappa_{22}) = \tau_{31}, \qquad (67) $$

$$ \tau_{13}(\eta_{12} - \zeta\kappa_{11}) + \tau_{23}(1 + \eta_{22} - \zeta\kappa_{21}) = \tau_{32}. \qquad (68) $$

Adapting to matrix form,

$$ \begin{bmatrix} 1 + \eta_{11} + \zeta\kappa_{12} & \eta_{21} + \zeta\kappa_{22} \\ \eta_{12} - \zeta\kappa_{11} & 1 + \eta_{22} - \zeta\kappa_{21} \end{bmatrix} \begin{bmatrix} \tau_{13} \\ \tau_{23} \end{bmatrix} = \begin{bmatrix} \tau_{31} \\ \tau_{32} \end{bmatrix} \qquad (69) $$

or,

$$ \begin{bmatrix} \tau_{13} \\ \tau_{23} \end{bmatrix} = \begin{bmatrix} 1 + \eta_{11} + \zeta\kappa_{12} & \eta_{21} + \zeta\kappa_{22} \\ \eta_{12} - \zeta\kappa_{11} & 1 + \eta_{22} - \zeta\kappa_{21} \end{bmatrix}^{-1} \begin{bmatrix} \tau_{31} \\ \tau_{32} \end{bmatrix}. \qquad (70) $$

One can define

$$ \beta = \begin{bmatrix} \beta_{11} & \beta_{12} \\ \beta_{21} & \beta_{22} \end{bmatrix} = \begin{bmatrix} 1 + \eta_{11} + \zeta\kappa_{12} & \eta_{21} + \zeta\kappa_{22} \\ \eta_{12} - \zeta\kappa_{11} & 1 + \eta_{22} - \zeta\kappa_{21} \end{bmatrix}^{-1} \qquad (71) $$

and, from (70) and (71),
$$ \begin{bmatrix} \tau_{13} \\ \tau_{23} \end{bmatrix} = \begin{bmatrix} \beta_{11} & \beta_{12} \\ \beta_{21} & \beta_{22} \end{bmatrix} \begin{bmatrix} \tau_{31} \\ \tau_{32} \end{bmatrix}. \qquad (72) $$

or

$$ \tau_{13} = \beta_{11}\tau_{31} + \beta_{12}\tau_{32}, \qquad (73) $$

$$ \tau_{23} = \beta_{21}\tau_{31} + \beta_{22}\tau_{32}, \qquad (74) $$

Introducing (73) and (74) in (50) and (51),

$$ \tau_{11,1} + \tau_{21,2} + (\kappa_{12}(\beta_{11}\tau_{31} + \beta_{12}\tau_{32}) - \kappa_{13}\tau_{12}) + (\kappa_{22}(\beta_{21}\tau_{31} + \beta_{22}\tau_{32}) - \kappa_{23}\tau_{22}) + \tau_{31,3} = -\bar{b}_1, \qquad (75) $$

$$ \tau_{12,1} + \tau_{22,2} + (\kappa_{13}\tau_{11} - \kappa_{11}(\beta_{11}\tau_{31} + \beta_{12}\tau_{32})) + (\kappa_{23}\tau_{21} - \kappa_{21}(\beta_{21}\tau_{31} + \beta_{22}\tau_{32})) + \tau_{32,3} = -\bar{b}_2. \qquad (76) $$

Grouping in relation to $\tau_{31}$ and $\tau_{32}$

$$ (\beta_{11}\kappa_{12} + \beta_{21}\kappa_{22})\tau_{31} + (\beta_{12}\kappa_{12} + \beta_{22}\kappa_{22})\tau_{32} + \tau_{31,3} = -(\bar{b}_1 + \tau_{11,1} + \tau_{21,2} - \tau_{12}\kappa_{13} - \tau_{22}\kappa_{23}), \qquad (77) $$

$$ (-\beta_{11}\kappa_{11} - \beta_{21}\kappa_{21})\tau_{31} + (-\beta_{12}\kappa_{11} - \beta_{22}\kappa_{21})\tau_{32} + \tau_{32,3} = -(\bar{b}_2 + \tau_{12,1} + \tau_{22,2} + \tau_{11}\kappa_{13} + \tau_{21}\kappa_{23}). \qquad (78) $$

By fitting the expression

$$ \frac{\partial t}{\partial \zeta} - Bt = b, \qquad (79) $$

one finds

$$ B_{11} = -\beta_{11}\kappa_{12} - \beta_{21}\kappa_{22}, \qquad (80) $$

$$ B_{12} = -\beta_{12}\kappa_{12} - \beta_{22}\kappa_{22}, \qquad (81) $$

$$ B_{21} = \beta_{11}\kappa_{11} + \beta_{21}\kappa_{21}, \qquad (82) $$

$$ B_{22} = \beta_{12}\kappa_{11} + \beta_{22}\kappa_{21}, \qquad (83) $$

$$ b_1 = -(\bar{b}_1 + \tau_{11,1} + \tau_{21,2} - \tau_{12}\kappa_{13} - \tau_{22}\kappa_{23}), \qquad (84) $$

$$ b_2 = -(\bar{b}_2 + \tau_{12,1} + \tau_{22,2} + \tau_{11}\kappa_{13} + \tau_{21}\kappa_{23}), \qquad (85) $$

and thus assembling the matrixes

$$ t = \begin{bmatrix} \tau_{31} \\ \tau_{32} \end{bmatrix}, \qquad (86) $$

$$ B = \begin{bmatrix} B_{11} & B_{12} \\ B_{21} & B_{22} \end{bmatrix}. \qquad (87) $$

$$ b = \begin{bmatrix} b_1 \\ b_2 \end{bmatrix}. \qquad (88) $$

we may write the system of two ordinary differential equations presented in (79).

Applying Crank-Nicholson method for each layer ($k = 1,2,...,n_{Layer}$), with $\zeta^0 = -h^{bottom}$, $\zeta^k = \zeta^{k-1} + h^k$, where $h^k$ is the thickness of the k-th layer, we get from (79),

$$ \frac{1}{h^k}(t^{k+1} - t^{k-1}) - B^k \frac{1}{2}(t^{k+1} + t^{k-1}) = b^k, \qquad (89) $$

where $B^k$ and $b^k$ are computed at the middle of the layers, i.e. $B^k = B(\zeta^{k-1} + \frac{1}{2}h^k)$ and $b^k = b(\zeta^{k-1} + \frac{1}{2}h^k)$. As for the stress, $t^{k+1} = t(\zeta^k)$ with the initial condition given as $t^0 = \bar{t}^{bottom}$. Isolating the current stress vector, one finds
$$ t^{k+1} = \left(I - \frac{1}{2}h^k B^k\right)^{-1} \left[ \left(I + \frac{1}{2}h^k B^k\right) t^{k-1} + h^k b^k \right], \quad k = 1, 2, \dots, n_{\text{layers}}. \qquad (90) $$

In here, we consider that the above equation can be used starting from the bottom or from the top, as such we have the same equation starting from the top as

$$ t^{j+1} = \left(I + \frac{1}{2}h^j B^j\right)^{-1} \left[ \left(I - \frac{1}{2}h^j B^j\right) t^{j-1} - h^j b^j \right], \quad j = 1, 2, \dots, n_{\text{layers}}. \qquad (91) $$

with the initial condition $t^0 = \bar{t}^{\text{top}}$.

In order to consider the two boundary conditions, here we consider the approach presented by (Soltani and Kordkheili, 2021), in which we apply an interpolation from the results obtained by finite difference starting from the bottom and from the top as

$$ t = \frac{(1 - \zeta_s)}{2} t_{\uparrow} + \frac{(1 + \zeta_s)}{2} t_{\downarrow}. \qquad (92) $$

The variable $\zeta_s$ is referenced to the middle layer of the shell, varying from $-1 \le \zeta_s \le 1$ at bottom and top of the shell thickness, respectively.

# 6 RESULTS

## 6.1 Cantilever beam (locking check)

The current example is based on the example proposed by Simo, Fox and Rifai (1990), revisited by Wriggers and Gruttmann (1993), Campello, Pimenta and Wriggers (2003) and Sanchez, Pimenta and Ibrahimbegovic (2023). In the first situation, however, a linear elastic behavior is used to verify if the element presents a locking behavior. As for the geometry, shown in Figure 3, a cross section of $b = 0.1$ is used, with length $L = 1.0$, and with varying values of thickness, ranging from $10^{-5} \le h \le 1$. As for the mesh, we have as an example the one presented in Figure 3, equivalent to 40 elements. With this, we check from a range of thick up to considerably thin shell, with increasing amount of elements. Regarding the boundary conditions, one end of the beam is clamped, while the other is subjected to a concentrated force $P \sim h^3$. This load-thickness relation is used in order to maintain the example allways at the linear elastic range, as in (Campello, Pimenta and Wriggers, 2003). The beam material has an Young's modulus $E = 10^7$ and Poisson's ratio $\nu = 0.3$ for all the layers, with Neo-Hookean material model. In the current situation, 1 and 9 layers where considered, with results for each case presented in Figure 3. The two cases where used to verify the stability of the element with different layer distributions. From the results, one can see the convergence behavior as the mesh is refined for all the cases, with no locking identified.
<figure><img src="image_4.png" /><figcaption>Figure 3 - Mesh refinement 1 layer and 9 layers – Linear case</figcaption></figure>

The example was also tested in a nonlinear range, with the load of $P = 1000$ and thickness $h = 0.1$ used. For such case, we find the results in Figure 10, compared with the T6-3i and T63iKL elements. One can see that the element still does not present any locking behavior.

<figure><img src="image_5.png" /><figcaption>Figure 4 - Mesh refinement 1 layer and 9 layers – Nonlinear case</figcaption></figure>
## 6.2 Cantilever Beam

This example uses the same mesh and general geometry from the previous example, however, with nonlinear behavior present. In the current case, the mesh is the same as in Figure 3, with fixed cross-section dimensions of $b = 0.1 \times h = 0.1$. One end of the beam is clamped, while the other is subjected to a concentrated force $P = 1000$. The beam material has an Young's modulus $E = 10^7$ and Poisson's ratio $\nu = 0.3$ for all the layers, with Neo-Hookean material model.

The values for comparison are the results of displacement of the midpoint of the free end of the beam, seen in Figure 5. In the current work, the results for 1 layer and 50 layers are presented, in order to exemplify the stability of the element with different layer distributions. From the results, one can see that solution path closely resembles the one obtained for the T6-3iKL from (SANCHEZ, PIMENTA and IBRAHIMBEGOVIC, 2023) and, in the same way as in this reference, deviates slightly from the results of the T6-3i element from (CAMPELLO, PIMENTA and WRIGGERS, 2003). This is expected, as the T6-3iKL is a Kirchhoff-Love type element, while the T6-3i is a Reissner-Mindlin type element.

<figure><img src="image_6.png" /><figcaption>Figure 5 - In and out of plane cantilever beam</figcaption></figure>

## 6.3 Raasch Challenge

The Raasch challenge is an example that, due to the coupling of four different modes of deformation (Shear, Twist, stretch and bending), is considered a challenge. The problem is composed of a curved strip hook, with first radius $R_1 = 14$", second radius $R_2 = 46$" and thickness $t = 2$". The left side of the hook is clamped, and the right is subjected to an upward in-plane total shear force of 1.0 lb. Regarding the mesh, we assume a variable $N$, which is the number of lateral divisions and, along the hook curves, there are a total of 6N divisions, equally divided between each curve (see Figure 6 for geometry, example mesh and deformed configurations). As presented by (SANCHEZ, PIMENTA and IBRAHIMBEGOVIC, 2023), the vertical expected displacement is 4.9352 in, obtained with the finite element method and with a different kinematical model, which explain the slight difference in results from the ones presented by (SANCHEZ, PIMENTA and IBRAHIMBEGOVIC, 2023) and the new proposed element. The results for multiple refinements with different layers layout (1 and 5 layers), proved to be acceptable as they have a small divergence in relation to the expected value (see Figure 7 for the results and comparison with multiple other authors and elements).
$w$ (width) = 20 in = 50,8 cm
$R_1$ = 46 in = 116,84 cm; $E$ = 3300 psi = 22,753 MPa

$R_2$ = 14 in = 35,56 cm; $v$ = 0.35;

$t$ = 2 in = 5,08 cm; $q$ = $0.05 \frac{\text{lbf}}{\text{in}} = 8.756 \frac{\text{N}}{\text{m}}$

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

Mesh

Deformed configuration
and mesh

<figure><img src="image_9.png" /></figure>

Shear force

<figure><img src="image_10.png" /><figcaption>Figure 7 - Raasch challenge vertical displacements</figcaption></figure>

## 6.4 Lateral Buckling cantilever beam

The buckling cantilever beam example analyzes the ability of the element to describe buckling behaviour. This was proposed by (SILVA, 2020), adapted from the work of (PIMENTA and Yojo, 1993) As for the geometry and the mesh, shown in Figure 8, a cross section of $b$ = 11.64 mm x $h$ = 100 mm is used, with length $L$ = 2400 mm. One end of the beam is clamped, while the other, parallel to $h$, is subjected to a concentrated force $P$ = 4500 N and a lateral load, $P_L$ = 1 N, is applied just to impose the initial instability into the system. The beam material has an Young's modulus $E$ = 210 GPa and Poisson's ratio $v$ = 0.3125 for all the 3 layers adopted, with Neo-Hookean material model. The geometry with mesh, deformed configuration and the lateral displacements of the tip of the beam are presented in Figure 8. As seen in the results, the solution path closely resembles the other elements, predicting essentially the same buckling load.
<figure><img src="image_11.png" /><figcaption>Figure 8 - Lateral and vertical displacements / Undeformed and deformed configurations of buckling cantilever beam</figcaption></figure>

## 6.5 Ring plate loaded at free edge

This example is a classic in the literature. A circular ring plate with a radial rip is loaded at the free edge by $q = 0.1\text{fkN/m}$, with the other edge being fully clamped. The external and internal plate radius are $R = 10.0\text{ m}$ and $r = 6.0\text{ m}$, respectively, with thickness $h = 0.03\text{ m}$. The material considered is Neo-Hookean, with $E = 2.1 \cdot 10^8\text{kN/m}^2$ and $v = 0$. The mesh along with the deformed configuration and the vertical displacement in two points (A, B) are presented in Figure 9. As expected, the results obtained for the 1 layer and 5 layers cases are closer to the Kirchhoff-Love type elements (T6-3iKL and T6-3iKLPen) than from the Reissner-Mindlin type element (T6-3i).

<figure><img src="image_12.png" /></figure>
Figure 9 - Vertical Displacements / Undeformed and deformed configurations of ring plate loaded at free edge

## 6.6 Pull-out of an open cylinder

This example consists of applying two diametrically opposite point forces in an open cylinder (see Figure 10). The two layers, modeled with Neo-Hookean material, have $E = 10.5 \cdot 10^6$ and $\nu = 0.3125$. In Figure 10 the deflection at point A versus pulling force graph is presented along with the final deformed configuration. From the results, one can see that the element is capable of representing slight snap-through behaviour of the solution when $P \approx 2.0 \cdot 10^4$, in the same fashion as the solid shell, T6-3iKLPen and the T6-3i elements, showing improved behavior in relation to the results from the T6-3iKL. In the current example, the new element presents a higher deformability when compared with the other kirchoff- love models. This could be due to the higher interpolation polynomial from the T10 when compared with the T6 types, allowing for more flexibility and deformation. We remark that the solution path still closely follows the same pattern seen in the other elements.

<figure><img src="image_13.png" /><figcaption>Figure 10 - Vertical Displacement / Undeformed and deformed configurations of pull out of open cylinder</figcaption></figure>

## 6.7 Pinched cylinder

This is another classic in the literature. The cylinder has a radius of $R = 100$, length of $L = 200$ and thickness $h = 1$. Additionally, the cylinder has a rigid diaphragm on both sides, in a similar fashion to a metal can. The material used is Neo-hookean for all the 3 layers, with Young's modulus $E = 3 \cdot 10^4$ and Poisson's coefficient $\nu = 0.3$. The pinch force applied symmetrically at top and bottom is $P = 20000$. The mesh and results for vertical and lateral displacement of points A and B, are presented in Figure 11. From the results, the current element presents itself as more flexible than the other elements, possibly due to the higher order polynomial interpolation. We remark that the solution path still closely follows the same pattern seen in the other elements.
<figure><img src="image_14.png" /><figcaption>Figure 11 - Pinched cylinder</figcaption></figure>

## 6.8 Anisotropic cantilever beam

The current example was developed to demonstrate the multilayer capabilities of the multilayer theory and firstly presented in (Gomes et al., 2024). The geometry, load and base material properties for the isotropic part of the material are the same as the configuration in the cantilever beam example, presented in section 6.2., with 5 equally sized layers composing the shell thickness. The difference in relation to the previous example lies in the additional anisotropic material contribution. Regarding the anisotropic parameters throughout the section, we always assume that $\alpha_1 = 5 \cdot 10^6$ and $\alpha_2 = 2.3$ are associated with orthotropic direction vector $m_1$ and that $\alpha_1 = 1 \cdot 10^5$ and $\alpha_2 = 2.3$ are associated with orthotropic direction vector $m_2$.

We assume initially that there are two possible configurations for the orthotropic directions, the first being $m_1 = (1,0,0)$ and $m_2 = (0,1,0)$ and the second being $m_1 = (0,1,0)$ and $m_2 = (1,0,0)$. With these directions, four possible plies configurations, presented in Figure 12 and named from ANISO1 up to ANISO4, are tested. The first and fourth configurations are chosen as they represent the extremes of anisotropic behavior. The intermediate configurations are adopted as to verify the increasingly lower value of stiffness due to the assumed directions of anisotropy.

<figure><img src="image_15.png" /><figcaption>Figure 12 - Plies configurations along the thickness – ANISO1 up to ANISO4</figcaption></figure>

The results for the relation between vertical deflection and the current beam X axis position for all plies configurations and for the purely Neo-Hookean material are presented in Figure 13. As expected, all the anisotropic cases rendered more stiff final configurations when compared with the Neo-Hookean case. As also expected, configurations ANISO1 and ANISO4 respectively resulted in the most and least stiff behaviors when compared with the other anisotropic cases.

It is worth mentioning that, although the vertical deflection for ANISO1 and ANISO2 were close, they were not identical. The number of layer configuration 1 along the thickness of ANISO2 greatly influenced the general behavior of the displacement when compared with ANISO3. As seen in Figure 13, The behavior of the cases are almost identical as it was observed in (Gomes et al., 2024), with the T6-3iKLMI element.
<figure><img src="image_16.png" /><figcaption>Figure 13 - Vertical deflection x Current beam axis position – ANISO1 up to ANISO4</figcaption></figure>

Additionally to verifying the influence of different configurations of plies when applying the anisotropy in the direction of the axis, we also verified the influence of rotating the vectors $m_1$ and $m_2$. As such, two new configurations are used for the orthotropic directions, the first being $m_1 = (\frac{\sqrt{2}}{2}, \frac{\sqrt{2}}{2}, 0)$ and $m_2 = (\frac{\sqrt{2}}{2}, -\frac{\sqrt{2}}{2}, 0)$ and the second being $m_1 = (\frac{\sqrt{2}}{2}, -\frac{\sqrt{2}}{2}, 0)$ and $m_2 = (\frac{\sqrt{2}}{2}, \frac{\sqrt{2}}{2}, 0)$. These new directions are numbered as 3 and 4 for compatibility with the previous model. With these directions, three possible plies configurations, presented in Figure 14 and named ANISO5 up to ANISO7, are tested.

<figure><img src="image_17.png" /><figcaption>Figure 14 - Plies configurations along the thickness – ANISO5 up to ANISO7</figcaption></figure>

The first and third configurations in Figure 14 are chosen as they represent the extremes of anisotropic directions for such case. The intermediate configuration is adopted as to verify the effect in the stiffness due to the assumed direction of anisotropy.

The results for the relation between vertical deflection and the current beam axis position for plies configurations ANISO5 up to ANISO7 and for the purely neo-Hookean material are presented in Figure 15.
<figure><img src="image_18.png" /><figcaption>Figure 15 - Vertical deflection x Current beam axis position – ANISO5 up to ANISO7</figcaption></figure>

As expected, due to boundary conditions directions and the assumed angle of anisotropy, the results for ANISO5 and ANISO7 are identical. Also, due to the same factors, one can see that for ANISO6 there is an increase in stiffness in relation to ANISO5 and ANISO7. This is due to the preferential stiffness direction associated with direction $m_1$. As ANISO5 and ANISO7 assume a specific 45° angle, they generate a smaller stiffness in relation to the loading and displacement directions when compared with the case in which we have a variation of such 45° angle, such as ANISO6. This behavior is also identified in the reference element (T6-3iKLMI) as all the cases results are nearly identical to the new element.

Another interesting behavior to observe is the lateral deflection difference, due to the anisotropic behavior. Obviously, when observing the Neo-hookean material reference, there is no noticeable lateral deflection for any side, as the material is isotropic. When observing Figure 16 for the extreme cases (ANISO5 and ANISO7), however, one can notice that the preferential stiffness resultant from the anisotropy generates a lateral deflection (in the case mirrored due to the proposed configuration of the vectors $m_1$ and $m_2$). As expected for the intermediate case, ANISO6, one can see an intermediate lateral deflection, closer to ANISO5 as there are more plies 3 than 4 in such case. The same behavior is seen when an inverted mesh is used.
<figure><img src="image_19.png" /><figcaption>Figure 16 - Lateral deflections</figcaption></figure>

## 6.9 Wrinkling of membrane

This is another example used to verify the anisotropic capability of the shell multi-layered element. This example is divided in two situations. Both of these follow and expand on the results presented in (Gomes et al., 2024) and (SANCHEZ, PIMENTA and IBRAHIMBEGOVIC, 2023).

As general information, the square membrane has side equal to 1 and beveled corners to reduce stress concentration. The vertical corners are clamped, and the horizontal corners are stretched by a distributed force $q = 0.1$. The isotropic parameters for the Neo-Hookean part are $E = 200$ and $v = 0.3$. For a visual reference of geometry and general informations, see Figure 17.

<figure><img src="image_20.png" /><figcaption>Figure 17 - Wrinkling of a membrane: Mesh and boundary conditions</figcaption></figure>
To ensure the wrinkling behaviour, only at the first step of the analysis, alternating positive and negative displacements of 0.00003 are applied normal to the shell plane, in some points along the middle line of the mesh (see Figure 17). Regarding the anisotropic parameters, we assume that $\alpha_1 = 40$ and $\alpha_2 = 2.3$ are associated with orthotropic direction vector $m_1$ and that $\alpha_1 = 1$ and $\alpha_2 = 2.3$ are associated with orthotropic direction vector $m_2$. We assume initially that there are two possible configurations for the orthotropic directions, the first being $m_1 = (0,1,0)$ and $m_2 = (1,0,0)$ and the second being $m_1 = (1,0,0)$ and $m_2 = (0,1,0)$. With these directions, three possible plies configurations, presented in Figure 18 and named from ANISO1 up to ANISO3, are tested with five layers. The first and third configurations are chosen as they represent the extremes of anisotropic behavior. The intermediate configuration is adopted as to verify the increasingly lower value of stiffness due to the assumed directions of anisotropy.

<figure><img src="image_21.png" /><figcaption>Figure 18 - Wrinkling of a membrane: Plies configurations along the thickness – ANISO1 up to ANISO3</figcaption></figure>

The results of displacement for points along the middle line of the mesh are presented in Figure 19. As it can be seen, the general behavior of the deformation of the fabric is close in relation to the results from (Gomes et al., 2024), and in the plies configurations equivalent to the ones presented by (SANCHEZ, PIMENTA and IBRAHIMBEGOVIC, 2023), these being ANISO1 and ANISO3. As also expected, configurations ANISO1 and ANISO3 respectively resulted in the least and most stiff behaviors when compared with the intermediate value ANISO2. The location of ANISO2, between both and closer to ANISO1 is logical, as there are more layers from the first than the second direction (see Figure 18).

<figure><img src="image_22.png" /><figcaption>Figure 19 - Vertical displacement of middle line in wrinkling membrane</figcaption></figure>
## 6.10 Stress analysis of simply supported plate with Neo-hookean material

This example is dedicated to analyse the performance of the shear stress formulation here presented. It consists of a simply supported square plate, with length $L = 2000$ mm and thickness $h = 0.2$ mm, equally divided along 50 layers, subjected to a distributed load $q = h^3 \cdot 10^{-14}$. The material of the plate is determined with Young's modulus $E = 10^3 N/mm^2$ and Poisson's coefficient $\nu = 0.3$. A Neo-hookean material law is used. For a reference of the geometry, boundary conditions and mesh, see Figure 20 a) and b).

<figure><img src="image_23.png" /><figcaption>Figure 20 – Simply supported plate: a)General description; b) Mesh and stress reference point c) vertical deflection results</figcaption></figure>

In order to verify the displacements, the equation

$$ w_{maxLin} = 0.00406 \frac{qL^4}{(Eh^3/(12(1-\nu^2)))'} \quad (93) $$

Obtained from (TIMOSHENKO and WOINOWSKY-KRIEGER, 1959) for a linear elastic simply supported plate submitted to a distributed load is used as a reference. When normalized with such expression, we obtain

$$ \frac{w_{max}}{w_{maxLin}} = 1.000766, \quad (94) $$

Indicating excellent agreement from the new element and linear elastic theory for the current example.

With the displacement verified, we now focus on the stress. In order to obtain reference value, one uses the equation

$$ w = \frac{qL^4}{(Eh^3/(12(1-\nu^2)))} \sum_{m=1}^{21} \left( \frac{4}{\pi^5 m^5} + A_m \cosh \frac{m\pi y}{L} + B_m \frac{m\pi y}{L} \sinh \frac{m\pi y}{L} \right) \sin \frac{m\pi x}{L}, \quad m = 1, 3, 5, \dots, 21 \quad (95) $$

$$ A_m = - \frac{2(\alpha_m \tanh \alpha_m + 2)}{\pi^5 m^5 \cosh \alpha_m} \quad (96) $$

$$ B_m = \frac{2}{\pi^5 m^5 \cosh \alpha_m} \quad (97) $$

$$ \alpha_m = \frac{m\pi L}{2L} \quad (98) $$

once again from (TIMOSHENKO and WOINOWSKY-KRIEGER, 1959), considering that the plate is square. From such equation, we use the relation
$$V_x = \frac{Eh^3}{12(1 - v^2)} (w_{xx} + w_{yy})_{,x}, \qquad (99)$$

$$V_y = \frac{Eh^3}{12(1 - v^2)} (w_{xx} + w_{yy})_{,y}, \qquad (100)$$

with the derivatives determined using the software Wolfram Mathematica, and determine the stresses

$$\tau_{31} = -\frac{3V_x}{2h} \left[ 1 - \left( \frac{z}{h/2} \right)^2 \right], \qquad (101)$$

$$\tau_{32} = -\frac{3V_y}{2h} \left[ 1 - \left( \frac{z}{h/2} \right)^2 \right], \qquad (102)$$

With z being the point analysed along the thickness of the shell. The shear stresses are determined at four points in the plate surface $P_1 = \{\frac{L}{4}, \frac{L}{2}, 0\}$, $P_2 = \{\frac{3L}{4}, \frac{L}{2}, 0\}$, $P_3 = \{\frac{L}{2}, \frac{L}{4}, 0\}$ and $P_4 = \{\frac{L}{2}, \frac{3L}{4}, 0\}$, indicated at Figure 20 b). As such, we have the shear stresses results in Figure 21 and Figure 22.

<figure><img src="image_24.png" /><figcaption>Figure 21 – Shear stress $\tau_{31}$ at points</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Figure 22 – Shear stress $\tau_{32}$ at points</figcaption></figure>

One can see in the results that, for $\tau_{31}$ in points 1 and 2 and $\tau_{32}$ in points 3 and 4, the diagrams from both the current element and the reference linear theory are in agreement, with a slight difference. In relation to the other situations, one can see that the new element presents stress values, these being around $1/10^{th}$ of the stress in the other direction. When analysing the linear theory, these other points give exactly zero, due to the nature of the assumed displacement and consequent stresses obtained. This discrepancy does not eliminate the validity of the new shear stress formulation, as the displacement field assumed in the current element is more precise than the reference theory, and the new formulation still presents a considerable difference between the stresses $\tau_{31}$ and $\tau_{32}$ at each point, as expected.

## 6.11 Stress analysis of cantilever Beam with Neo-hookean material (Same Material)

The cantilever beam presented in 6.1 is used as a reference for the analysis of the shear stress. For this, we use 5 points along the beam (see Figure 23). These points are checked for the stress, considering a set of multiple possible loads applied at the edge of the cantilever beam, there are $P = \{0.0001,0.001,0.01,0.1,1,10,100,200,300,400,500\}$. The increase in load is considered in order to verify the influence of change in stress distribution as the example loads start in a linear domain and increasingly move towards the nonlinear domain. For this example, 50 plies are considered along the thickness.

<figure><img src="image_26.png" /><figcaption>Figure 23 - Analysed points in beam for shear stress distribution</figcaption></figure>
As comparison, we use the linear elastic beam shear formula $\tau = \frac{VQ}{IT}$ in the beam axis direction for each load. Due to the exact nature of the shell formulation, it is immediately visible that there is a variation of the stress behavior along the beam axis in the results. The results closer to the linear elastic reference is obtained at point 5. Obviously, the smaller the load, the closer the result obtained from the current element and the linear elastic reference.

<figure><img src="image_27.png" /><figcaption>Figure 24 – Shear stress $\tau_{31}$: Loads from $P = 0.0001$ to $P = 10.0$</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Figure 25 – Shear stress $\tau_{31}$: Loads from $P = 100.0$ to $P = 500.0$</figcaption></figure>

It is worth to remark that, as the load increases, we see how the nonlinear material behavior influences the shear stress. In the smaller loads, the shape of the shear stress diagram is that of a parabole, with the maximum value at the middle of the thickness. When the loads begin to increase, it is perceptible a change in the point of maximum shear stress, this being distorted to a point below the middle of the thickness, indicating that the nonlinear material and geometric aspect of the theory becomes increasingly more pronounced.

In order to exemplify the behavior of the normal stress, we take $\sigma_{11}$ obtained directly from the first Piola-Kirchhoff stress tensor at point 1 with the new element and compare it with the linear elastic beam normal stress formula $\sigma_{11} = \frac{Mc}{l}$. As seen in Figure 26, for the loads between 0.0001 and 10, we see a diagram close related to the standard linear
elastic theory, with the change from compression to traction happening at the middle of the thickness of the shell and with general results close to the reference.

<figure><img src="image_29.png" /><figcaption>Figure 26 – Normal stress $\sigma_1$ at point 1: Loads from $P = 0.0001$ to $P = 10.0$</figcaption></figure>

With the load $P = 100$, one can identify the translation from the origin for the change from traction to compression (see Figure 27), indicating the influence of not only bending, but also normal forces in the normal stress. This behavior coincides with the change of the position of shear stress maximum from the middle thickness to a position below the middle thickness (see Figure 25).

<figure><img src="image_30.png" /></figure>
<caption>Figure 27 – Normal stress $σ_1$ at point 1: Loads from $P = 0.0001$ to $P = 10.0$</caption>

As for the other loads results, the change continuously increases, as seen in Figure 28.

<figure><img src="image_31.png" /><figcaption>Figure 28 – Normal stress $σ_1$ at point 1: Loads from $P = 100.0$ to $P = 500.0$</figcaption></figure>

## 6.12 Stress analysis of cantilever Beam with Neo-hookean material (Variable Material Properties)

The cantilever beam presented in 6.1 is used as a reference for the analysis of the shear stress, however, in the current case, we analyze a case in which different material properties compose the thickness. For this example, 50 plies are considered along the thickness, with intercalated material properties of Young's modulus $E = 10^7$ and $E = 10^6$ at groups of 5 plies, with all being composed of Poisson's ratio $v = 0.3$. In this case, we use 4 points along the beam (see Figure 23). These points are checked for the stress, considering a set of multiple possible loads applied at the edge of the cantilever beam, there are $P = \{0.0001,0.001,0.01,0.1,1,10,100,200,300,400,500\}$. The increase in load is to check the influence of change in stress distribution as the example loads start in a linear domain and increasingly move towards the nonlinear domain.

<figure><img src="image_32.png" /><figcaption>Figure 29 - Analysed points in beam for variable shear stress distribution</figcaption></figure>

Due to the exact nature of the shell formulation, it is immediately visible that there is a variation of the stress behavior along the beam axis in the results. As observed in Figure 30 and Figure 31, it is possible to clearly identify the variation in the shear stress due to the different material properties at each group of 5 plies.
<figure><img src="image_33.png" /><figcaption>Figure 30 – Variable material - Shear stress $\tau_{31}$: Loads from $P = 0.0001$ to $P = 10.0$</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Figure 31 – Variable material - Shear stress $\tau_{31}$: Loads from $P = 100.0$ to $P = 500.0$</figcaption></figure>

It is worth to remark that, as the load increases, we see how the nonlinear material behavior influences the shear stress. In the smaller loads, the shape of the shear stress diagram is more stable, with better defined behavior along the shell thickness. When the loads begin to increase, it is perceptible a change in general behavior of the diagrams, this being distorted from the smaller loads more stable behavior, indicating that the nonlinear material and geometric aspect of the theory becomes increasingly more pronounced.

In order to exemplify the behavior of the normal stress, we take $\sigma_{11}$ obtained directly from the first Piola-Kirchhoff stress tensor at point 1. For all the load cases, one can see the clear steps formed due to the different material properties assumed at each group of 5 plies.
<figure><img src="image_35.png" /><figcaption>Figure 32 – Normal stress $\sigma_1$ at point 1: Loads from P = 0.0001 to P = 500.0</figcaption></figure>

## 6.13 Stress analysis of cantilever Beam with transversely isotropic material

In this example we check the behavior of the shear stress formulation when used with transversely isotropic materials in three different cases of structural vectors (fibers) directions. The geometry, mesh, boundary conditions, isotropic material properties and evaluated points in this example are the same as in the previous example for all the three different cases. In the current case, a transversely isotropic material energy contribution is added in the way described by equation (26) with anisotropic parameters $\alpha_1 = 5 \cdot 10^6$ and $\alpha_2 = 2.0$ applied to all layers in direction $m_1$.

The first case considered is constructed with all layers assuming $m_1 = \{1,0,0\}$. The results for the shear stress $\tau_{31}$ for such configuration are presented in Figure 33 and Figure 34. The results in point 1 of $\sigma_{11}$, obtained directly from the first Piola-Kirchhoff stress tensor, is presented in Figure 35.
<figure><img src="image_36.png" /><figcaption>Figure 33 – Shear stress $\tau_{31}$ (Case 1): Loads from $P = 0.0001$ to $P = 10.0$</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>Figure 34 – Shear stress $\tau_{31}$ (Case 1): Loads from $P = 100.0$ to $P = 500.0$</figcaption></figure>
<figure><img src="image_38.png" /><figcaption>Figure 35 – Normal stress $\sigma_1$ at point 1 (Case 1): Loads from $P = 0.0001$ to $P = 500.0$</figcaption></figure>

In the results one can see the influence of the fiber addition to the material. The material energy assumed for the transversely isotropic material assumes contribution only when there is traction at the point. If the point is under compression, no material contribution associated with the fiber is considered. This is due to the Macaulay brackets present in equation (26). One can clearly see the influence of such behavior in Figure 35, in the sharp change in the normal stress along the thickness of the shell. This is also clearly observed in the sudden change of shear stress behavior in Figure 33 and Figure 34, at the same point along the thickness of the normal stress.

The second case considered is constructed with alternating layers assuming $m_1 = \{1,0,0\}$ and $m_1 = \{0,1,0\}$, respectively, starting from the bottom of the thickness and anisotropic parameters $\alpha_1 = 5 \cdot 10^6$ and $\alpha_2 = 2.0$. The results for the shear stress $\tau_{31}$ for such configuration are presented in Figure 36 and Figure 37. The results in point 1 of $\sigma_{11}$, obtained directly from the first Piola-Kirchhoff stress tensor, is presented in Figure 38.
<figure><img src="image_39.png" /><figcaption>Figure 36 – Shear stress $\tau_{31}$ (Case 2): Loads from $P = 0.0001$ to $P = 10.0$</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>Figure 37 – Shear stress $\tau_{31}$ (Case 2): Loads from $P = 100.0$ to $P = 500.0$</figcaption></figure>
<figure><img src="image_41.png" /><figcaption>Figure 38 – Normal stress $\sigma_1$ at point 1 (Case 2): Loads from $P = 0.0001$ to $P = 500.0$</figcaption></figure>

In this case, the results also show the influence of the fiber addition to the material. Once again we see the influence of the material fibers being activated only in traction. One can clearly see the influence of such behavior in Figure 38, in the sharp change in the normal stress along the thickness of the shell, with clear steps seen in the traction portion of the diagram. This is also clearly observed in the sudden changes of shear stress behavior in Figure 36 and Figure 37, at the same point along the thickness of the normal stress. Due to the extremely complex interactions of the material both in linear and non-linear domains, the stresses pattern is extremely variable along the shell thickness.

The third case considered is constructed with 10 layers assuming $\alpha_1 = 5 \cdot 10^7$ and $\alpha_2 = 2.0$ at the top and bottom and $\alpha_1 = 5 \cdot 10^6$ and $\alpha_2 = 2.0$ at the rest of the layers. The results for the shear stress $\tau_{31}$ for such configuration are presented in Figure 39 and Figure 40. The results in point 1 of $\sigma_{11}$, obtained directly from the first Piola-Kirchhoff stress tensor, is presented in Figure 41.
<figure><img src="image_42.png" /><figcaption>Figure 39 – Shear stress $\tau_{31}$ (Case 3): Loads from $P = 0.0001$ to $P = 10.0$</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>Figure 40 – Shear stress $\tau_{31}$ (Case 3): Loads from $P = 100.0$ to $P = 500.0$</figcaption></figure>
<figure><img src="image_44.png" /><figcaption>Figure 41 – Normal stress $\sigma_1$ at point 1 (Case 3): Loads from P = 0.0001 to P = 500.0</figcaption></figure>

In this case, the results also show the influence of the fiber addition to the material. Once again we see the influence of the material fibers being activated only in traction, with no influence of the fibers in the compression portion of the diagram. One can clearly see the influence of such behavior in Figure 41, in the sharp change in the normal stress along the thickness of the shell. This is also clearly observed in the sudden changes of shear stress behavior in Figure 39 and Figure 40, at the same point along the thickness of the normal stress.

# 7 CONCLUSION

As seen in the current work, the context of a thin shell with Kirchhoff-Love theory and consequent assumption of a constant direction of the shell director connected with a multi-layer approach and incremental DoFs, create an extremely simple multi-layer finite element. Although lagrangean parameters are used in order to obtain $C^1$ continuity, the results from the numerical implementation indicate that, although simple, the element performed exceedingly well with a singular or multiple layers in various non-linear situations including large bending in and out-of-plane movement, buckling phenomena, complex load-displacement curves (with snaps involved) and anisotropic membrane. It is also worth to remark that, in the context of composite finite elements, the amount of degrees of freedom used is still quite small, as they do not increase with the number of layers present. Regarding the shear stress formulation, although presented here in a limited capability due to the considerable length of the work, it proved itself capable of being used in a fully nonlinear context, with both geometric and material nonlinearities present. In a future dedicated work, a more thorough analysis of the stress formulation is to be made. The authors highlight that the current work is a step toward a dynamic formulation capable of simulation of composites and their stress behavior.

**Author's Contributions:** Conceptualization, Gustavo Canário Gomes, Paulo de Mattos Pimenta and Adnan Ibrahimbegovic; Methodology, Gustavo Canário Gomes, Paulo de Mattos Pimenta and Adnan Ibrahimbegovic;
Investigation, Gustavo Canário Gomes, Paulo de Mattos Pimenta and Adnan Ibrahimbegovic; Writing - original draft, Gustavo Canário Gomes and Paulo de Mattos Pimenta; Writing - review & editing, Gustavo Canário Gomes, Paulo de Mattos Pimenta and Adnan Ibrahimbegovic; Funding acquisition, Paulo de Mattos Pimenta; Supervision, Paulo de Mattos Pimenta and Adnan Ibrahimbegovic.

**Editor:** (Filled by: Journal Staff)

## REFERENCES

Argyris, J. H., Fried, I., Scharpf, D. W. (1968). The TUBA family of plate elements for the matrix displacement method. The Aeronautical Journal of the Royal Aeronautical Society, v. 72, p. 701-709. https://doi.org/10.1017/S000192400008489X

Argyris, J. H., Lochner, N. (1972). On the application of the SHEBA shell element. Computer Methods in Applied Mechanics and Engineering, v. 1, p. 317-347. https://doi.org/10.1016/0045-7825(72)90012-6.

Brank, B., Korelc, J., Ibrahimbegović, A. (2002) Nonlinear shell problem formulation accounting for through-the-thickness stretching and its finite element implementation. Comput Struct 80(9). P. 699-717. https://doi.org/10.1016/S0045-7949(02)00042-1.

Campello, E. M. B.; Pimenta, P. M.; Wriggers, P. (2003). A triangular finite shell element based on a fully nonlinear shell formulation. Computational Mechanics, v. 31, p. 505-518. https://doi.org/10.1007/s00466-003-0458-8.

Campello, E. D. M. B.; Pimenta, P. M., Wriggers, P. (2011). An exact conserving algorithm for nonlinear dynamics with rotational dofs and general hyperelasticity. part 2: shells. Computational Mechanics, p. 195-211. https://doi.org/10.1007/s00466-011-0584-7.

Cirak F, Ortiz M. (2001). Fully C1-conforming subdivision elements for finite deformation thin-shell analysis. International Journal for Numerical Methods in Engineering, v. 51, p. 813-833. https://doi.org/10.1002/nme.182.

Ibrahimbegovic A., B. Brank, P. Courtois, (2001) Stress resultant geometrically exact form of classical shell model and vector-like parameterization of constrained finite rotations, International Journal for Numerical Methods in Engineering, 52 p.1235-1252. https://doi.org/10.1002/nme.247.

Ivannikov, V., Tiago, C., Pimenta, P.M. (2014). Meshless implementation of the geometrically exact Kirchhoff-Love shell theory. International Journal for Numerical Methods in Engineering, v. 100, p. 1-39. https://doi.org/10.1002/nme.4687.

Ivannikov, V., Tiago, C., Pimenta, P.M. (2015). Generalization of the C1 TUBA plate finite elements to the geometrically exact Kirchhoff-Love shell model. Computer Methods in Applied Mechanics and Engineering, v. 294, p. 210-244. https://doi.org/10.1016/j.cma.2015.05.018.

Kiendl, J, Bletzinger, K-U, Linhard, J, Wüchner, R. (2009). Isogeometric shell analysis with Kirchhoff-Love elements. Computer Methods in Applied Mechanics and Engineering, v. 198, p. 3902-3914. https://doi.org/10.1016/j.cma.2009.08.013.

Pimenta, P. M.; Neto, E. S. A.; Campello, E. M. B. (2010). A Fully Nonlinear Thin Shell Model of Kirchhoff-Love type. New Trends in Thin Structures: Formulation, Optimization and Coupled Problems, p. 29-58. https://doi.org/10.1007/978-3-7091-0231-2_2.

Pimenta, P. M. ; Yojo, T. (1993). Geometrically exact analysis of spatial frames. Applied Mechanics Review, v. 46, n. 11, p. 118-128. http://dx.doi.org/10.1115/1.3122626.

Rabczuk T, Areias P.M.A., Belytschko T. (2007). A meshfree thin shell method for non-linear dynamic fracture. International Journal for Numerical Methods in Engineering, v. 72, p. 524-548. https://doi.org/10.1002/nme.2013.

Sanchez, M. L.; Pimenta, P. M., Silva, C. C. (2020). A simple fully nonlinear Kirchhoff-Love shell finite element. Latin American Journal of Solids and Structures, v. 17, n. 8. https://doi.org/10.1590/1679-78256120.

Sanchez, M. L.; Pimenta, P. M.; Ibrahimbegovic, A (2023). A simple geometrically exact finite element for thin shells — Part 1: statics. Computational Mechanics, 72 p. 1119-1139. https://doi.org/10.1007/s00466-023-02339-2.

Schröder, J. ; Neff, P. (2010). Poly-, quasi-and rank-one convexity in applied mechanics. Springer Science & Business. https://doi.org/10.1007/978-3-7091-0174-2.
Schröder, J., Viebahn, N., Balzani D., Wriggers P. (2016). A novel mixed finite element for finite anisotropic elasticity; the SKA-element Simplified Kinematics for Anisotropy. Computer Methods in Applied Mechanics and Engineering, p. 475-494. https://doi.org/10.1016/j.cma.2016.06.029.

Silva, C. C. (2020). Modelo geometricamente exato rígido ao cisalhamento de cascas e barras. Ph.D. Thesis, University of São Paulo, São Paulo. https://doi.org/10.11606/T.3.2020.tde-12042021-143412.

da Costa e Silva, C., Maassen, S.F., Pimenta, P.M. and J. Schröder (2019). A simple finite element for the geometrically exact analysis of Bernoulli–Euler rods. Computational Mechanics, p. 1-19. https://doi.org/10.1007/s00466-019-01800-5.

Silva, C.C., Maassen, S. F., Pimenta, P. M., Schröder, J. (2021). On the simultaneous use of simple geometrically exact shear-rigid rod and shell finite elements, Computational Mechanics, vol. 67, no. 3, pp. 867-881. https://doi.org/10.1007/s00466-020-01967-2.

Simo, J. C.; Fox, D. D. ; Rifai, M. S. (1990). On a stress resultant geometrically exact shell model. Part III: computational aspects of the nonlinear theory. Computer Methods in Applied Mechanics and Engineering., p. 21-70. https://doi.org/10.1016/0045-7825(90)90094-3.

Reddy, J. N. (2006). Theory and analysis of elastic plates and shells. Boca Raton: CRC Press. 568 p.
https://doi.org/10.1201/9780849384165.

Bohinc, U. . Brank, B.. Ibrahimbegovic, A. (2014). Discretization error for the Discrete Kirchhoff plate finite element approximation.
Computer Methods in Applied Mechanics and Engineering, v. 269, p. 415-436. https://doi.org/10.1016/j.cma.2013.11.011.

Viebahn, N.; Pimenta, P. ; Schröder, J. (2017). A simple triangular finite element for nonlinear thin shells: statics, dynamics and anisotropy. Computational Mechanics, p. 281-297. https://doi.org/10.1007/s00466-016-1343-6.

Wu, T.P., Pimenta, P.M. & Wriggers, P. (2024). On triangular virtual elements for Kirchhoff–Love shells. Arch Appl Mech. https://doi.org/10.1007/s00419-024-02591-9.

Wriggers, P. ; Gruttmann, F. (1993). Thin shells with finite rotations formulated in Biot stresses: theory and finite element formulation. Int. J. Numer. Meth. Engrg., p. 2049-2071. https://doi.org/10.1002/nme.1620361207.