A fully incremental simple triangular multilayer Kirchhoff-Love shell element

Abstract

This paper presents a new triangular multi-layer nonlinear shell finite element with incremental degrees of freedom, suitable for large displacements and rotations. This is a nonconforming element with 6 nodes, quadratic displacement and linear rotation field based on Rodrigues incremental rotation parameters, with a total of 21 DoFs. The novelty of this element is the extension to a multilayer, fully incremental situation of the T6-3iKL element, a kinematical model with properties from Kirchhoff-Love theory, approximating the shell director across layers as constant. The model is numerically implemented, and results are compared to different references in multiple examples, showing the capabilities of the formulation. It is believed that the possibly simplest multilayer extension, combined with fully incremental DoFs, simple kinematic, no necessity of artificial parameters such as penalties, 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 create a simple yet powerful shell element.

KeywordsTriangular Shell Element; Multilayer shell; Nonlinear Shell Formulation; Kirchhoff-Love shell; Large Strain; Incremental degrees of freedom

Graphical Abstract

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 structures 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), 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 ${\mathrm{C}}^1$-continuity in the interface between adjacent elements, which may be a challenge. Among the possibilities to deal with such ordeal, there is the well stablished TUBA family elements, initially developed for plates by Argyris, Fried and Scharpf (1968), expanded for shells by Argyris and Lochner (1972) and recently revisited by Ivannikov, Tiago and Pimenta (2015). Other possible approach is to impose rotation continuity with the use of penalty parameters at the midside node and with mesh refinement, such as done by Viebahn, Pimenta and Schröder (2017) and Sanchez, Pimenta and Silva (2020). Yet another possibility, presented by Cirak and Ortiz (2001), is the use of subdivision surfaces, with progressive division of the mesh until the ${\mathrm{C}}^1$ continuity is achieved. For the meshless case, as seen in Ivannikov, Tiago and Pimenta (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, Campello (2010), Rabczuk, Areias, and Belytschko (2007), Kiendl, Bletzinger, Linhard and Wüchner (2009), Bohinc, Brank and Ibrahimbegovic (2014), Silva, Maassen, Pimenta and Schröder (2021) and Wu, Pimenta and Wriggers (2024).

In this work we present an extension of the triangular Kirchhoff-Love shell element T6-i3KL developed by Sanchez, Pimenta and Ibrahimbegovic (2023) to the multi-layer consideration with incremental degrees of freedom, creating a new element.

The main multi-layer consideration used to build the element is 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. This consideration greatly simplifies the multi-layer kinematics and eliminate the need of multiple independent variables along the 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 incremental rotation. In the present work, we have a nonconforming element with 6 nodes, a quadratic displacement and a linear rotation field based on Rodrigues incremental rotation parameters, giving in total 21 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. The degrees of freedom of the element are all incremental. The adaptation from a full displacement formulation, presented by Sanchez, Pimenta and Ibrahimbegovic (2023), to an incremental displacement is a step towards a formulation in dynamics (following and improving upon the previous work in Campello, Pimenta and Wriggers (2011)). We note that the main novelty of this work is the multilayer formulation. The fully incremental formulation is an additional uniqueness when compared with the initial semi-incremental formulation present at the reference work.

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

2 SHELL FORMULATION

2.1 Multi-Layer consideration

The main consideration 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.

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

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.

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 $\left\{{\mathit{e}}_1^r{,\mathrm{}\mathit{e}}_2^r,{\mathit{e}}_3^r\right\}$ where ${\mathit{e}}_{1\mathrm{}}^r$ and ${\mathit{e}}_2^r$ are in the middle plane of the shell and ${\mathit{e}}_3^r$ is normal to the plane. In the current configuration, the orthogonal coordinate system is $\left\{{\mathit{e}}_1,{\mathit{e}}_2,{\mathit{e}}_3\right\}$, with the relation ${\mathit{e}}_i=\mathit{Q}{\mathit{e}}_i^r$ valid, where $\mathit{Q}$ is the rotation tensor.

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

$\mathbf{\xi }=\mathbf{\zeta }+{\mathit{a}}^r,\mathrm{}\mathbf{\zeta }={\mathrm{\xi }}_{\alpha }{\mathit{e}}_{\alpha }^r\mathrm{}\in {\Omega }^r\mathrm{}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{}{\mathit{a}}^r=\mathrm{\zeta }{\mathit{e}}_3^r\mathrm{}\in H^r,$

(1)

where $\mathbf{\zeta }$ are the coordinates of a material point in the shell middle plane and ${\mathit{a}}^r\mathrm{}$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$\mathrm{}\mathrm{\zeta }\in \mathrm{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

$\mathit{x}=\mathit{z}+\mathit{a},\mathrm{}\mathit{z}=\mathit{\zeta }+\mathit{u}\mathrm{}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{}\mathit{a}=\mathit{Q}{\mathit{a}}^{\mathit{r}},\mathrm{}$

(2)

$\mathit{z}$ represents the current position of a point in the middle surface and $\mathit{a}$ the position along the thickness of the shell, in ${\mathit{e}}_3$ direction.

Adopting the notation $\partial \left(\bullet \right)/\partial \left({\xi }_{\alpha }\right)={\left(\bullet \right)}_{,\alpha }$, the deformation gradient can be obtained from:

$\mathbf{F}={\mathit{x}}_{,\alpha }\otimes {\mathit{e}}_{\alpha }^r+\frac{\partial \mathit{x}}{\partial \zeta }\otimes {\mathit{e}}_3^r={\mathit{f}}_a\otimes {\mathit{e}}_{\alpha }^r+{\mathit{f}}_3\otimes {\mathit{e}}_3^r=\left({\mathit{z}}_{,\alpha }+{\mathit{\kappa }}_{\alpha }\times \mathit{a}\right)\otimes {\mathit{e}}_{\alpha }^r+\mathit{Q}{\mathit{e}}_3^r\otimes {\mathit{e}}_3^r.$

(3)

where

${\mathit{K}}_{\alpha }={\mathit{Q}}_{,\alpha }{\mathit{Q}}^{\mathrm{T}}\mathrm{}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{}{\mathit{\kappa }}_{\alpha }=\mathrm{a}\mathrm{x}\mathrm{i}\mathrm{a}\mathrm{l}\left({\mathit{K}}_{\alpha }\right).$

(4)

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

${\mathit{F}}^r={\mathit{Q}}^{\mathrm{T}}\mathit{F}=\mathit{I}+{\mathit{\gamma }}_{\alpha }^r\otimes {\mathit{e}}_{\alpha }^r,\mathrm{}$

(5)

with

${\mathit{\gamma }}_{\alpha }^r={\mathbf{\eta }}_{\alpha }^{\mathit{r}}+{\mathit{\kappa }}_{\alpha }^r\times {\mathit{a}}^r\mathrm{};\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{\mathit{\eta }}_{\alpha }^r={\mathit{Q}}^T{\mathit{z}}_{,\alpha }-{\mathit{e}}_{\alpha }^r;{\mathit{\kappa }}_{\alpha }^r=\mathrm{a}\mathrm{x}\mathrm{i}\mathrm{a}\mathrm{l}\left({\mathit{Q}}^{\mathrm{T}}{\mathit{Q}}_{,\alpha }\right).$

(6)

In the previous equations, the vectors ${\mathit{\eta }}_{\alpha }^r$ and ${\mathit{\kappa }}_{\alpha }^r$ are the cross-sectional generalized strains.

It is important to remark that, due to the Kirchhoff-Love assumption, ${\mathit{\gamma }}_{\alpha }^r\bullet {\mathit{e}}_3^r={\mathit{\eta }}_{\alpha }^r\bullet {\mathit{e}}_3^r=0$.

2.3 Rodrigues rotation vector

The rotation vector can be represented as $\mathit{\theta }=\theta \mathit{e}$, with $\theta$ being the rotation around the axis $\mathit{e}$. As presented in Campello, Pimenta and Wriggers (2003), the rotation tensor can be a function of the Rodrigues rotation vector as $\mathit{\alpha }=\frac{\mathit{t}\mathit{a}\mathit{n}\left(\theta /2\right)}{\theta /2}\mathit{\theta }$, restricting the rotation at $-\pi <\theta <+\pi$. As such, the rotation tensor can be defined as

$\mathit{Q}=\mathit{I}+h\left(\alpha \right)\left(\mathit{A}+\frac{1}{2}{\mathit{A}}^2\right),\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}h\left(\alpha \right)=\frac{4}{4+{\alpha }^2},$

(7)

$\alpha =‖\mathit{\alpha }‖$ and $\mathit{A}=\mathrm{S}\mathrm{k}\mathrm{e}\mathrm{w}\left(\mathit{\alpha }\right)$.

For other possible descriptions of the Rodrigues rotation dealing with large displacements and rotations, we refer to the work of Ibrahimbegovic, Brank and Courtois (2001).

2.4 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 \left(\mathit{u}\right)={\int }_B^{}\left(\psi \left(\mathit{C}\left(\mathit{u}\right)\right)-{\rho }_0\stackrel{-}{\mathit{b}}\bullet \mathit{u}\right)dV-{\int }_{\partial B}^{}\left(\stackrel{-}{\mathit{t}}\bullet \mathit{u}\right)dA.$

(8)

The directional derivative of the functional is obtained as

$\delta \prod \left(\mathit{u}\right)=D\prod \left(\mathit{u}\right)\bullet \mathit{v}=\frac{\partial }{\partial \alpha }\prod {\left(\mathit{u}+\alpha \mathit{v}\right)}_{\alpha =0}=0.$

(9)

One can define the internal total energy for the shell as

$\prod \left(\mathit{u}\right)=U^{int}={\int }_{\Omega }^{}{\int }_{H^r}^{}\psi \left({\xi }_1^r,{\xi }_2^r\right)dHd\Omega ={\int }_{\Omega }^{}\widehat{\psi }\left({\xi }_1^r,{\xi }_2^r\right)d\Omega .$

(10)

The strain energy is, obviously, material dependent, with two options presented in the current work. The internal stress and material tangent moduli are obtained by

$\mathit{P}=\frac{\partial \psi }{\partial \mathit{F}}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{\mathit{D}}_{\alpha \beta }^r=\frac{\partial {\mathit{\sigma }}_{\alpha }^r}{\partial {\mathit{\epsilon }}_{\alpha }^r}=\frac{{\partial }^2\stackrel{-}{\psi }}{\partial {\mathit{\epsilon }}_{\alpha }^r\partial {\mathit{\epsilon }}_{\beta }^r}.$

(11)

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}\left(J^2-1\right)-\ln \left(\mathrm{}J\mathrm{}\right)\right)+\frac{1}{2}\mu \left(I_1-3-2\ln \left(\mathrm{}J\right)\right),$

(12)

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

$I_1=tr\mathrm{}\mathit{C}={\mathit{f}}_i\bullet {\mathit{f}}_i,\mathrm{}\mathrm{}\mathrm{}\mathrm{}{I}_2=\mathit{t}\mathit{r}\left(\mathrm{C}\mathrm{o}\mathrm{f}\mathrm{}\mathit{C}\right)={\mathit{g}}_i\bullet {\mathit{g}}_i\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{}\mathrm{}\mathrm{}{I}_3=\det \mathit{C}=J^2={\left({\mathit{f}}_1\bullet \left({\mathit{f}}_2\times {\mathit{f}}_3\right)\right)}^2,$	(13)
$\mathit{C}={\mathit{F}}^{\mathrm{T}}\mathit{F},$	(14)
$\lambda =\frac{E\nu }{\left(1+\nu \right)\left(1-2\nu \right)}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mu =\frac{E}{2\left(1+\nu \right)}.$	(15)

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 \left(\mathit{C}\right)={\psi }^{i\_p}\left(\mathit{C}\right)+{\psi }^{a\_p}\left(\mathit{C}\right)={\psi }^{i\_p}\left(\mathit{C}\right)+\sum _i^2{\alpha }_1^{\left(i\right)}{〈I_4^{\left(i\right)}-1〉}^{{\alpha }_2^{\left(i\right)}},$

(16)

where ${\alpha }_1^{\left(i\right)}$ and ${\alpha }_2^{\left(i\right)}$ are variables in relation to material properties in the ${\mathit{m}}^{\left(i\right)}$ direction, $"〈〉"$ are Macaulay brackets and $I_4^{\left(i\right)}$ is an invariant defined as

$I_4^{\left(i\right)}=tr\left[A_{\alpha \beta }^{\left(i\right)}{\mathit{f}}_{\alpha }\bullet {\mathit{f}}_{\beta }\right],$

(17)

with

$A_{\alpha \beta }^{\left(i\right)}=\left({\mathit{e}}_{\alpha }\bullet {\mathit{m}}^{\left(i\right)}\right)\left({\mathit{e}}_{\beta }\bullet {\mathit{m}}^{\left(i\right)}\right)={\mathit{e}}_{\alpha }\bullet {\mathit{M}}^{\left(i\right)}{\mathit{e}}_{\beta }\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{\mathit{M}}^{\left(i\right)}={\mathit{m}}^{\left(i\right)}\otimes {\mathit{m}}^{\left(i\right)}.$

(18)

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

2.5 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 ${\mathit{\gamma }}_{33}^r={\gamma }_{33}{\mathit{e}}_3^r$ to eq. (5), rewriting it as

${\mathit{F}}^r=\mathit{I}+{\mathit{\gamma }}_{\alpha }^r\otimes {\mathit{e}}_{\alpha }^r+{\mathit{\gamma }}_{33}^r\otimes {\mathit{e}}_3^r.\mathrm{}$

(19)

The plane stress is then enforced by ${\tau }_{33}={\mathit{\tau }}_3^r\bullet {\mathit{e}}_3^r=0$.

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

${\tau }_{33}=2\left(1+{\gamma }_{33}\right)\frac{\partial \psi }{\partial I_1}+\frac{\partial \psi }{\partial J}\stackrel{-}{J},$

(20)

where $\stackrel{-}{J}={\mathit{e}}_3^r\bullet \left(f_1^r\times f_2^r\right)$. After some manipulations, one arrives at

${\gamma }_{33}=\sqrt{\frac{\lambda +2\mu }{\lambda {\stackrel{-}{J}}^2+2\mu }}-1.$

(21)

3 INCREMENTAL DESCRIPTION OF DEGREES OF FREEDOM

3.1 General incremental relations

The incremental nature of the formulation takes into account some general relations. We consider a mean value correlation

${(\bullet )}_{1/2}=\frac{1}{2}\left({(\bullet )}_{i+1}+{(\bullet )}_i\right).\mathrm{}$

(22)

and

${(\bullet )}_{1/2}+\frac{1}{2}∆\left(\bullet \right)={(\bullet )}_{i+1},\mathrm{}$	(23)
${(\bullet )}_{1/2}-\frac{1}{2}∆\left(\bullet \right)={(\bullet )}_i.\mathrm{}$ (24)

We also define a general rule for a generic product of generic quantities as

$∆\left(\mathrm{A}\mathrm{*}\mathrm{B}\right)={\mathrm{A}}_{1/2}\mathrm{*}∆\mathrm{B}+∆\mathrm{A}\mathrm{*}{\mathrm{B}}_{1/2}.\mathrm{}$ (25)

These relations are used throughout the section. It is important to remark that the incremental aspects of the DoF in the theory follow the same logic of the presented by Campello, Pimenta and Wriggers (2011), with alterations when required.

3.2 Incremental description of rotation

The incremental rotational variables are presented as

${\mathit{Q}}_{\mathrm{\Delta }}=\mathit{I}+h\left({\alpha }_{\mathrm{\Delta }}\right)\left({\mathit{A}}_{\mathrm{\Delta }}+\frac{1}{2}{\mathit{A}}_{\mathrm{\Delta }}^2\right),\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}h\left({\alpha }_{\mathrm{\Delta }}\right)=\frac{4}{4+{\alpha }_{\mathrm{\Delta }}^2},$

(26)

where ${\alpha }_{\mathrm{\Delta }}=‖{\mathit{\alpha }}_{\mathrm{\Delta }}‖$ and ${\mathit{A}}_{\mathrm{\Delta }}=\mathrm{S}\mathrm{k}\mathrm{e}\mathrm{w}\left({\mathit{\alpha }}_{\mathrm{\Delta }}\right).$

The update of the rotation vector is done by

${\alpha }_{\mathrm{i}+1}=\frac{4}{4-{\mathit{\alpha }}_{\mathrm{i}}\bullet {\mathit{\alpha }}_{\mathrm{\Delta }}}\left({\mathit{\alpha }}_{\mathrm{i}}+{\mathit{\alpha }}_{\mathrm{\Delta }}-\frac{1}{2}{\mathit{\alpha }}_{\mathrm{i}}\times {\mathit{\alpha }}_{\mathrm{\Delta }}\right),$

(27)

where ${\mathit{\alpha }}_{i+1}$ is the next step of ${\mathit{\alpha }}_{\mathrm{i}}$. We also update the rotation tensor following the relations

${\mathit{Q}}_{\mathrm{i}+1}={\mathit{Q}}_{\mathrm{\Delta }}{\mathit{Q}}_{\mathrm{i}},\mathrm{}\mathrm{}\mathrm{}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{}\mathrm{}\mathrm{}{\mathit{Q}}_{\mathrm{i}+1}=\widehat{Q}\left({\mathit{\alpha }}_{\mathrm{i}+1}\right),\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{\mathit{Q}}_{\mathrm{\Delta }}=\widehat{Q}\left({\mathit{\alpha }}_{\mathrm{\Delta }}\right)\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{}\mathrm{}\mathrm{}{\mathit{Q}}_i=\widehat{Q}\left({\mathit{\alpha }}_{\mathrm{i}}\right).$

(28)

As $\mathit{Q}$ applies to the shell director, one may also define

${\mathit{e}}_3^{i+1}={\mathit{Q}}_{\mathrm{\Delta }}{\mathit{e}}_3^i.$

(29)

Here we bring the description of the rotation vector in relation to a scalar parameter ${\phi }_{\mathrm{\Delta }}$, first presented by da Costa e Silva et al (2019) in relation to vector ${\mathit{e}}_3$ and, here, described in relation to ${\mathit{e}}_1$ to obtain the rotation vector in agreement with Sanchez, Pimenta and Ibrahimbegovic (2023). This is done as ${\mathit{e}}_1$ can be calculated only by displacement. Thus, we have

${\mathit{\alpha }}_{\mathrm{\Delta }}=\frac{{\mathit{e}}_1^i\times {\mathit{e}}_1^{i+1}}{{‖{\mathit{e}}_1^m‖}^2}+{\phi }_{\mathrm{\Delta }}\frac{{\mathit{e}}_1^m}{‖{\mathit{e}}_1^m‖},\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{\mathit{e}}_1^m=\frac{1}{2}\left({\mathit{e}}_1^{i+1}+{\mathit{e}}_1^i\right).$

(30)

3.3 Rotation tensor correlations

In this section we focus in some useful rotation tensor correlations. These are initially presented for $\mathit{Q}$, but also hold for its incremental equivalents. Afterwards, we deal with specific definitions for the incremental quantities.

We begin with the Cayley transform of (7), resulting in

$\mathit{Q}=\left(\mathit{I}+\frac{1}{2}\mathit{A}\right){\left(\mathit{I}-\frac{1}{2}\mathit{A}\right)}^{-1}={\left(\mathit{I}-\frac{1}{2}\mathit{A}\right)}^{-1}\left(\mathit{I}+\frac{1}{2}\mathit{A}\right).$

(31)

Based on the previous equation and with some basic manipulations, one can obtain

$\frac{1}{2}\left(\mathit{I}+\mathit{Q}\right)={\left(\mathit{I}-\frac{1}{2}\mathit{A}\right)}^{-1}$

(32)

and

$\mathit{Q}-\mathit{I}=\mathit{A}{\left(\mathit{I}-\frac{1}{2}\mathit{A}\right)}^{-1}={\left(\mathit{I}-\frac{1}{2}\mathit{A}\right)}^{-1}\mathit{A}.$

(33)

Let us now focus on specific equations. We define the back-rotated incremental rotation tensor and the relation with the skew tensor of the incremental rotation, respectively

${\mathit{Q}}_{∆}^r={\mathit{Q}}_{i+1}^{\mathrm{T}}{\mathit{Q}}_{∆}{\mathit{Q}}_{i+1}={\mathit{Q}}_i^{\mathrm{T}}{\mathit{Q}}_{∆}{\mathit{Q}}_{\mathit{i}},$	(34)
${\mathit{Q}}_{∆}^r-\mathit{I}={\mathit{A}}_{∆}^r{\left(\mathit{I}-\frac{1}{2}{\mathit{A}}_{∆}\right)}^{-1}.$	(35)

We also need to obtain the increment of the rotation tensor

$∆\mathit{Q}={\mathit{Q}}_{i+1}-{\mathit{Q}}_i={\mathit{Q}}_i\left({\mathit{Q}}_{∆}^r-\mathit{I}\right).$

(36)

Considering (35), the equation above becomes

$∆\mathit{Q}={\mathit{Q}}_i{\mathit{A}}_{∆}^r{\left(\mathit{I}-\frac{1}{2}{\mathit{A}}_{∆}^r\right)}^{-1}={\mathit{Q}}_i{\mathit{A}}_{∆}^r{\mathit{Q}}_i^{\mathrm{T}}{\mathit{Q}}_i{\left(\mathit{I}-\frac{1}{2}{\mathit{A}}_{∆}^r\right)}^{-1}={\mathit{A}}_{∆}{\mathit{Q}}_{1/2}.$

(37)

With $∆\mathit{Q}$ defined, we now present ${\mathit{Q}}_{1/2}$. With the help of (22), (28) and (34), one can define

${\mathit{Q}}_{1/2}=\frac{1}{2}\left({\mathit{Q}}_{i+1}+{\mathit{Q}}_i\right)=\frac{1}{2}\left({\mathit{Q}}_{∆}{\mathit{Q}}_i+{\mathit{Q}}_i\right)=\frac{1}{2}\left({\mathit{Q}}_i{\mathit{Q}}_{∆}^r+{\mathit{Q}}_{\mathrm{i}}\right)=\frac{1}{2}{\mathit{Q}}_i\left({\mathit{Q}}_{∆}^r+\mathit{I}\right).\mathrm{}$

(38)

With equation (34), (38) becomes

${\mathit{Q}}_{1/2}=\frac{1}{2}{\mathit{Q}}_i\left({\mathit{Q}}_i^{\mathbf{T}}{\mathit{Q}}_{∆}{\mathit{Q}}_i+\mathit{I}\right)=\mathrm{}\frac{1}{2}\left({\mathit{Q}}_{∆}+\mathit{I}\right){\mathit{Q}}_i$

(39)

and with (32),

${\mathit{Q}}_{1/2}={\left(\mathit{I}-\frac{1}{2}{\mathit{A}}_{∆}\right)}^{-1}\mathrm{}{\mathit{Q}}_i.$

(40)

Using (31), the previous equations becomes

${\mathit{Q}}_{1/2}={\left(\mathit{I}+\frac{1}{2}{\mathit{A}}_{∆}\right)}^{-1}\mathrm{}{\mathit{Q}}_{i+1}\mathrm{}$

(41)

additionally,

$\mathrm{d}\mathrm{e}\mathrm{t}\left({\mathit{Q}}_{1/2}\right)=h\left({\alpha }_{∆}\right)=\frac{4}{4+{\alpha }^2}.$

(42)

3.4 Incremental strains

In the present work, the increment of displacements and rotations are considered as DOFs. This consideration facilitates the development and convergence of dynamic formulations, a possible next step of the current formulation.

The curvature vectors presented in eq. (4) and eq. (6) can be rewritten in relation to Rodrigues rotation parameters as

${\mathit{\kappa }}_{\alpha }^r={\mathbf{\Xi }}^{\mathrm{T}}{\mathit{\alpha }}_{,\alpha }\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{\mathit{\kappa }}_{\alpha }=\mathbf{\Xi }{\mathit{\alpha }}_{,\alpha }\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathbf{\Xi }=\frac{4}{4+{\alpha }^2}\left(\mathit{I}+\frac{1}{2}\mathit{A}\right).$

(43)

In the same fashion as Campello, Pimenta and Wriggers (2011), we use the incremental formulation of the curvature vectors as

${\mathit{\kappa }}_{i+1}^r={\mathit{\kappa }}_{\alpha i}^r+\mathrm{\Delta }{\mathit{\kappa }}_{\alpha }^r,\mathrm{}\mathrm{}\mathrm{}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{}\mathrm{}\mathrm{\Delta }{\mathit{\kappa }}_{\alpha }^r=\mathrm{d}\mathrm{e}\mathrm{t}\left({\mathit{Q}}_{1/2}\right){\mathit{Q}}_{1/2}^{-1}{\mathit{\alpha }}_{\mathrm{\Delta },\alpha }.$

(44)

Regarding the rotation, it is implemented incrementally and as so, at every step the rotation vector $\mathit{\alpha }$ is updated as

${\mathit{\alpha }}_{i+1}=\frac{4}{4-{\mathit{\alpha }}_i\bullet {\mathit{\alpha }}_{\mathrm{\Delta }}}\left({\mathit{\alpha }}_i+{\mathit{\alpha }}_{\mathrm{\Delta }}-\frac{1}{2}{\mathit{\alpha }}_i\times {\mathit{\alpha }}_{\mathrm{\Delta }}\right),$

(45)

with ${\mathit{\alpha }}_{\mathrm{\Delta }}$ updated as in eq. (30). Regarding the curvature vector, it is incremented as in eq. (44).

At every step, the increment parameters are reset as ${\mathit{\alpha }}_{\mathrm{\Delta }}=0,\mathrm{}\mathrm{}{\phi }_{\mathrm{\Delta }}=0,\mathrm{}\mathrm{}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{}\mathrm{\Delta }{\mathit{\kappa }}_{\alpha }^r=0$.

The displacements are updated according to the following expression

${\mathit{u}}_{i+1}={\mathit{u}}_{\mathrm{i}}+{\mathit{u}}_{∆},$

(46)

in which ${\mathit{u}}_{i+1}$ is the displacement of the next step, ${\mathit{u}}_{\mathrm{i}}$ the current and ${\mathit{u}}_{∆}$ the increment. One also need to define the vectors

${\mathit{z}}_{1/2,\alpha }=\frac{1}{2}\left({\mathit{z}}_{i+1,\alpha }+{\mathit{z}}_{i,\alpha }\right).\mathrm{}$	(47)
${\mathit{z}}_{i+1,\alpha }={\mathit{z}}_{1/2,\alpha }+\frac{1}{2}∆{\mathit{z}}_{,\alpha }.\mathrm{}$ (48)

With this, one can write the incremental strain as

$∆{\mathit{\eta }}_{\alpha }^{\mathrm{r}}=∆\left({\mathit{Q}}^{\mathrm{T}}{\mathit{z}}_{,\alpha }-{\mathit{e}}_{\alpha }^r\right)=∆{\mathit{Q}}^{\mathrm{T}}{\mathit{z}}_{1/2,\alpha }+{\mathit{Q}}_{1/2}^{\mathrm{T}}∆{\mathit{z}}_{,\alpha }.\mathrm{}$

(49)

considering (47), (48) and (37), and that ${\mathit{A}}^{\mathrm{T}}=-\mathit{A}$ the equation above can be rewritten as

$∆{\mathit{\eta }}_{\alpha }^{\mathrm{r}}=∆\left({\mathit{Q}}^{\mathrm{T}}{\mathit{z}}_{,\alpha }-{\mathit{e}}_{\alpha }^r\right)=-{\mathit{Q}}_{1/2}^{\mathrm{T}}{\mathit{A}}_{∆}{\mathit{z}}_{1/2,\alpha }+{\mathit{Q}}_{1/2}^{\mathrm{T}}{\mathit{u}}_{∆,\alpha }.\mathrm{}$

(50)

or, considering the cross product property and rearranging

$∆{\mathit{\eta }}_{\alpha }^{\mathrm{r}}={\mathit{Q}}_{1/2}^{\mathrm{T}}\left({\mathit{u}}_{{∆}_{,\alpha }}+{\mathit{z}}_{1/2,\alpha }\times {\mathit{\alpha }}_{∆}\right).\mathrm{}$

(51)

At every step, the increment parameters are reset as ${\mathit{u}}_{\mathrm{\Delta }}=0\mathrm{}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{}\mathrm{\Delta }{\mathit{\eta }}_{\alpha }^r=0$.

4 THE FINITE ELEMENT

The element created is an extrapolation of the T6-3iKL developed by Sanchez, Pimenta and Ibrahimbegovic (2023). The new element, entitled T6-3iKLMI is a six-noded displacement based triangular element (see Figure 2). The discretization of the element coordinates and displacement field is done with area coordinates and quadratic interpolation and, an additional scalar rotation parameter $\left({\phi }_{\mathrm{\Delta }}\right)$. The rotation vector ${\mathit{\alpha }}_{∆}$ is a linear non-conform field over the shell element obtained indirectly by the scalar DoF ${\phi }_{\mathrm{\Delta }}$ shared between adjacent elements. This model renders a quite efficient element, with a reduced number of DoF when compared with usual models that consider the rotation using three-dimensional vectors as DoF. This characteristic associated with the simplification of the continuous shell director along the multi-layers, constitutes an incredibly simple and efficient multi-layer element with incremental DoF.

Figure 2
T6-3iKLMI

This recent description of the rotation, proposed by Sanchez, Pimenta and Ibrahimbegovic (2023) is a quite interesting adoption, reducing the number of DoF involved in the analysis while not reducing the quality of results. Such advantage is fully used in the present work with the adoption of a constant rotation through the layers of the element due to the small thickness.

The element uses 3 integration points at the mid-side nodes of the shell plane and, when integrating along the thickness, each layer has 3 integration points, using Simpson’s 1/3 rule. Obviously, by integrating each layer separately, one allows for use of different materials and thickness at each 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 5 for examples, as all the mesh, deformed configurations and graphs were generated within AceFem).

5 RESULTS

5.1 Cantilever Beam

The cantilever beam was an 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) to analyze the behavior of the element when used under conditions of high rotations and displacement in and out of plane, in this work, with multiple layers possibility. As for the geometry and the mesh, shown in Figure 3, a cross section of $b=0.1\times h=0.1$ is used, with length $L=1.0$. One end of the beam is clamped, while the other is subjected to different values of a concentrated force$\mathrm{}P$, variable from 0 to 1000. The beam material has an Young’s modulus $E={10}^7$ and Poisson’s ratio $v=0.3$ for all the layers, with neo-Hookean material model.

Figure 3
– In and out of plane cantilever beam

The values for comparison are the results of displacement of the midpoint of the free end of the beam, which are compared with that obtained in the work of Wriggers and Gruttmann (1993), where the maximum deflection are $u_{in-plane}=0.8425$ and $u_{out-of-plane}=0.8406$.

5.2 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 geometry of the problem is presented in Figure 4. The problem is compose of a curved strip hook, with first radius ${\mathrm{R}}_1=14"$, second radius ${\mathrm{R}}_2=46"$ and thickness $\mathrm{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\mathrm{}\mathrm{l}\mathrm{b}$. 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. As presented by Sanchez, Pimenta and Ibrahimbegovic (2023), the vertical expected displacement is $4.3966\mathrm{}\mathrm{i}\mathrm{n}$, 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). The results for multiple refinements with multiple layers (from 1 up to 9 layers), seem to be acceptable as they have a small divergence in relation to the expected value (see Figure 4 for the normalized result in relation to the expected value for all values of N and number of layers).

Figure 4
– Raasch challenge geometry and normalized vertical displacements

5.3 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 5, a cross section of $b=11.64\mathrm{}\mathrm{m}\mathrm{m}\times h=100\mathrm{}\mathrm{m}\mathrm{m}$ is used, with length $L=2400\mathrm{}\mathrm{m}\mathrm{m}$. One end of the beam is clamped, while the other, parallel to $h$, is subjected to a concentrated force$\mathrm{}P=10000\mathrm{}\mathrm{N}$ and a lateral load, $P_L=1\mathrm{}\mathrm{N}$, is applied just to impose the initial instability into the system. The beam material has an Young’s modulus $E=210\mathrm{}\mathrm{G}\mathrm{P}\mathrm{a}$ and Poisson’s ratio $v=0.3125$ for all the 3 layers adopted, with neo-Hookean material model. The deformed configuration and the, respectively from left to right, lateral and vertical displacements of the tip of the beam are presented in Figure 5.

Figure 5
– Lateral and vertical displacements / Undeformed and deformed configurations of buckling cantilever beam

5.4 Ring plate loaded at free edge

This example is a classic in the literature, referring to Campello, Pimenta and Wriggers (2003) for results comparison. A circular ring plate with a radial rip is loaded at the free edge by $p=0.1\mathrm{f}\mathrm{k}\mathrm{N}/\mathrm{m}$, with the other edge being fully clamped. The external and internal plate radius are $R=10.0\mathrm{}\mathrm{m}$ and $r=6.0\mathrm{}\mathrm{m}$, respectively, with thickness $h=0.03\mathrm{}\mathrm{m}$. The material considered is neo-Hookean, with $E=2.1\bullet {10}^8\mathrm{k}\mathrm{N}/{\mathrm{m}}^2$ and $v=0$. The mesh along with the deformed configuration and the vertical displacement in three points (A, B, C) are presented in Figure 6.

Figure 6
– Vertical Displacements / Undeformed and deformed configurations of ring plate loaded at free edge

5.5 Pull-out of an open cylinder

This example consists of applying two diametrically opposite point forces in an open cylinder (see Figure 7). The three layers, modeled with neo-Hookean material, have $E=10.5\bullet {10}^6$ and $v=0.3125$. In Figure 7 the deflection at point A versus pulling force graph is presented along with the final deformed configuration. It is worth to remark the slight snap-through behaviour of the solution when $P\approx 2.0\bullet {10}^4$. We refer to Campello, Pimenta and Wriggers (2003) for results comparison.

Figure 7
– Vertical Displacement / Undeformed and deformed configurations of pull out of open cylinder

5.6 Pinched cylinder

This is another classic in the literature, referring to Sanchez, Pimenta and Ibrahimbegovic (2023) for results comparison. 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\bullet {10}^4$ and Poisson’s coefficient $v=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, respectively, are presented in Figure 8.

Figure 8
– Pinched cylinder

5.7 Anisotropic cantilever beam

The current example was developed to demonstrate the multilayer capabilities of the multilayer theory. The geometry, load and base material properties for the isotropic part of the material are the same as the out-of-plane configuration in the first example, presented in section 5.1, 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\bullet {10}^6$ and ${\alpha }_2=2.3$ are associated with orthotropic direction vector $m_1$ and that ${\alpha }_1=1\bullet {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=\left(\mathrm{1,0},0\right)$ and $m_2=\left(\mathrm{0,1},0\right)$ and the second being $m_1=\left(\mathrm{0,1},0\right)$ and $m_2=\left(\mathrm{1,0},0\right)$. With these directions, four possible plies configurations, presented in Figure 9 and named from ANISO1 up to ANISO4, are tested. The first and forth 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 9
– Plies configurations along the thickness – ANISO1 up to ANISO4

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 10. 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.

Figure 10
– Vertical deflection x Current beam axis position – ANISO1 up to ANISO4

It is worth mentioning that, although the vertical deflection for ANISO1 and ANISO2 were close, with tip vertical displacement respectively of 0.743629 and 0.746318, 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.

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=\left(\frac{\sqrt{2}}{2},\frac{\sqrt{2}}{2},0\right)$ and $m_2=\left(\frac{\sqrt{2}}{2},-\frac{\sqrt{2}}{2},0\right)$ and the second being $m_1=\left(\frac{\sqrt{2}}{2},-\frac{\sqrt{2}}{2},0\right)$ and $m_2=\left(\frac{\sqrt{2}}{2},\frac{\sqrt{2}}{2},0\right)$. The previous directions are numbered as 3 and 4 for compatibility with the previous model. With these directions, three possible plies configurations, presented in Figure 11 and named ANISO5 up to ANISO7, are tested.

Figure 11
– Plies configurations along the thickness – ANISO5 up to ANISO7

The first and third configurations in Figure 11 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 ANISO 7 and for the purely neo-Hookean material are presented in Figure 12.

Figure 12
– Vertical deflection x Current beam axis position – ANISO5 up to ANISO7

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.

5.8 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 cases. The first analysed case, presented in Viebahn, Pimenta and Schröder (2017), focus on verifying the capabilities of the element to represent the stretching phenomena with multiple layers, with the same material. This is used to verify the stability of the theory when dividing the thickness in multiple smaller domains. For the second case, we use the proposal presented by Sanchez, Pimenta and Ibrahimbegovic (2023) as two of the plies configurations and expand the analyzes to other possibilities that arise due to the multilayer approach. It is also worth to mention that different mesh configurations and methods to generate the instability needed to initiate the wrinkling process are used for each case.

As general information for both cases, the square membrane has side equal to 1 and beveled corners to avoid 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 $\nu =0.3.$ For a visual reference of geometry and general informations, see Figure 13.

Figure 13
– Wrinkling of a membrane

With the general information given, let us focus on the first case. To ensure the wrinkling behaviour, some of the points in the mesh are moved, in order to remove symmetry and allow the beginning of wrinkling. The material used is the anisotropic model presented in a previous section, with orthotropic direction $m_1={\left[\mathrm{0,1},0\right]}^{\mathrm{T}}$ and $m_2={\left[\mathrm{1,0},0\right]}^{\mathrm{T}}$. The anisotropic parameters are ${\alpha }_1=4$ and ${\alpha }_2=2.3$ for the warp direction and ${\alpha }_1=1$ and ${\alpha }_2=2.3$ for the weft direction. For comparison, we refer to Viebahn, Pimenta and Schröder (2017). The example was modeled with one to nine layers and the results for deformation of one and nine layers are presented in Figure 13, along with the mesh adopted. It is perceptible that the displacements are identical for both configurations.

With the first case complete, we shift to the second. 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 14). 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=\left(\mathrm{0,1},0\right)$ and $m_2=\left(\mathrm{1,0},0\right)$ and the second being $m_1=\left(\mathrm{1,0},0\right)$ and $m_2=\left(\mathrm{0,1},0\right)$. With these directions, three possible plies configurations, presented in Figure 14 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 14
– Wrinkling of a membrane: a) Mesh and boundary conditions b) Plies configurations along the thickness – ANISO1 up to ANISO3

The results of displacement for points along the middle line of the mesh are presented in Figure 15. As expected, the plies configurations equivalent to the ones presented by Sanchez, Pimenta and Ibrahimbegovic (2023), these being ANISO1 and ANISO3, rendered the same result as theirs. As also expected, configurations ANISO1 and ANISO3 respectively resulted in the most and least 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 14).

Figure 15
– Vertical displacement of middle line in wrinkling membrane

6 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. The fact that there is no need of penalty parameters and any numeric tool in the element formulation is also remarkable and desirable. 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. The authors highlight that the current work is a step toward the most appropriate formulation in dynamics, with the incremental aspect of the formulation being extremely advantageous for such endeavor.

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
  » 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
  » 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
  » 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
  » 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
  » 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
  » 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
  » 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
  » 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
  » 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
  » 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
  » 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
  » 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
  » 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
  » 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
  » 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
  » 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
  » 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
  » 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
  » 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
  » 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
  » 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
  » 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
  » 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
  » 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
  » 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
  » https://doi.org/10.1002/nme.1620361207

Edited by

• Editor:
  Marco L. Bittencourt

Publication Dates

• Publication in this collection
  07 Oct 2024
• Date of issue
  2024

History

• Received
  16 Apr 2024
• Reviewed
  03 July 2024
• Accepted
  13 Aug 2024
• Published
  19 Aug 2024