A simple fully nonlinear Kirchhoff-Love shell finite element

Abstract

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

Keywords:
Finite Element Method; Kirchoff Love shell; Non-linear

Graphical Abstract

1 INTRODUCTION

A numerical approach solution usually is an effective method of obtaining the structural response for situations where the load distribution, material non-linearity or geometrical properties are complex, thus turning the attainment of an analytical solution into a real challenge. Shells are a very important structural model, which can be regarded as a MOR (Model Order Reduction) of Solid Mechanics (Nigro et al., 2019) and are consistently applicable to engineering structural problems.

Theories for Plates, Shells and Membranes have been developed in the past decades by many researchers due to its importance in mechanics. Kirchhoff-Love and Reissner-Mindlin theories are example of shell/plate models used in engineering. In the current paper, the authors address a non-linear Kirchhoff-Love element implemented in the AceGEN / AceFEM code environment, which is a Mathematica based code suitable for developing optimized code for computational mechanics (Korelc, J.A., Wriggers, P., 2016).

In Finite Element Methods, this structural problem may be implemented and studied by three-dimensional elements. As stated in Korelc and Wriggers (2016) this strategy is, however, inadequate for certain occasions. It is shown that for especially thin shells, the locking-phenomena may occur, causing the simulated structure to appear stiffer than what it really is. Hence, special attention must be taken regarding the interpolation of 3-D elements for Shell modeling.

Another approach that have been used by some researches is the use of kinematical simplification for shells. With these models, the system can be modeled as a surface (2D) with a specified thickness, consequently with less degrees of freedom, but can yet deal with nonlinear phenomena and with finite deflections and Rotations. Some pioneer work has been developed by Simo and Fox (1989) and Pimenta (1993) developing new shell models before 2000’s and departing from the traditional so-called degenerated solid approach. Campello et al. (2003) and Pimenta et al. (2004) developed shear flexible models (Reissner-Mindling) based on kinematics proposed by Pimenta (1993). These works had successfully implemented a shear flexible nonlinear six-parameter triangular shell finite element. In the same work, finite rotations were treated by using Euler-Rodrigues formula in a very conviniet way (due to its incremental formulation). Viebahn et al. (2017), Costa e Silva et al. (2018), Costa e Silva (2020) and Costa e Silva et al. (2020) built finite element models based on Kirchoff-Love Hypothesis and are the basis for the current article. The main advantage for using Kirchoff-Love shells is that the element will carry less degrees of freedom per element, simplifying the numerical problem. However, as shown in the next sections, one must ensure C1-continuity in domain due to kinematical formulation. As the displacement field is quadratic, this continuity is guaranteed inside the element. The present work and the work of Viebahn et al. (2017) and Costa e Silva (2020) differ in the methodology to enforce C1 continuity in the element interfaces. The procedure adopted in this article can be regarded as a discrete Galerking approach (see Brezzi and Fortin, 1991).

2 OBJECTIVES

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

3 METHODOLOGY

3.1 Shell Kinemátics

Plates and shell theories are very diverse, and many models may be found in the bibliography considering different kinematical hypothesis. The current work is based on the Kirchhoff Plate theory which has three kinematical constrains, also known as Kirchhoff hypothesis (see Reddy, 2006):

1. Straight perpendicular lines (transverse normals) to the middle surface stay straight after any deformation.

2. These transverse normals are inextensible. Their length (thickness) are constant.

3. These transverse normals remain perpendicular to the middle surface during any deformation imposed to the shell. There is no shear strain in any shell cross section.

The Kirchhoff-Love shell model, also known as shear rigid because of the absence of shear strain in consequence of its kinematical assumptions, developed in this paper may be understood as a continuation of the work developed by Costa e Silva et al. (2018), Costa e Silva et al. (2020). and Viebahn et al. (2017). It includes a new methodology to approcimate C1 continuity between adjacent elements with penalty and is also is valid for finite deflections, rotations and strains.

A plane reference configuration is assumed, however initial curved shells may be represented by a a stress-free deformation (Pimenta and Campello, 2009). The model is illustrated by Figure 1:

Figure 1
Non-linear finite element for Kirchhoff-Love shells undergoing finite deflections and rotations. Shell description and basic kinematical quantities.

The main kinematic properties of the shell developed in this article are based on Pimenta, P.M., et al. (2010) and are defined by:

• An arbitraty point is parametrized in the reference configuration by the vector $\mathit{\xi }\in {\mathbb{R}}^3$ and in the current configuration by the vector $\mathit{x}\in {\mathbb{R}}^3$.

• The projection of the arbitrary point in Middle surface in the reference configuration is parametrized by the vector $\mathit{\zeta }\in {\mathrm{\Omega }}^r\subset {\mathbb{R}}^2$ and in the current configuration by the vector $\mathit{z}\in \mathrm{\Omega }\subset {\mathbb{R}}^3$.

• The vector that connects the arbitrary point to its projection in the middle surface is ${\mathit{a}}^{\mathit{r}}$ and in the current configuration by $\mathit{a}$. It is important to notice that ${\mathit{a}}^{\mathit{r}}$ and $\mathit{a}$ are perpendicular to ${\mathrm{Ù}}^r$ and $\mathrm{Ù}$ respectively.

• The shell volume and thickness of the reference shell is defined such as $V^r$ and $H^r=\left[-h_b^r,h_t^r\right]$. The Shell total thickness is $h^r=h_b^r+h_t^r$

An orthonormal right-handed coordinate system ${\mathit{e}}_{\mathit{i}}^{\mathit{r}}$ is defined in the reference configuration such as any point in the shell may be defined by equation (1) (Pimenta, P.M., et al., 2010)

$\xi =\zeta +a^r,\zeta ={\xi }_{\alpha }{e}_{\alpha }^r,{\xi }_{\alpha }\in {\Omega }^{r}and{a}^r={\xi }_3{e}_3^r,{\xi }_3\in H^r.$ (1)

Equivalently, an arbitrary point in the current configuration can also be defined in an associated orthonormal coordinate system ${\mathit{e}}_{\mathit{i}}$ by position vector $\mathit{x}$, rotation tensor $\mathit{Q}$ and displacement $\mathit{u}$ with expressions

$z=\zeta +u,a=Q{a}^{r}andx=z+a.$ (2)

It is important to mention that tensor $\mathit{Q}$ rotates vector ${\mathit{a}}^{\mathit{r}}$ to $\mathit{a}$, which are always perpendicular to shell middle surface, consequently $\mathit{Q}$ is important to describe shell curvature and related kinematical quantities. The approach of rotational quantities is better explained in further sections of this paper.

First and second spatial derivatives of $\mathit{z}$ are defined by

$z_{,\alpha }=\frac{\partial \left({\xi }_{\alpha }{e}_{\alpha }^r+u\right)}{\partial {\xi }_{\alpha }}=e_{\alpha }^r+u_{,\alpha },z{,}_{\alpha \beta }=u{,}_{\alpha \beta }and{(\cdot )}_{\alpha }=\frac{\partial (\cdot )}{\partial {\xi }_{\alpha }}.$ (3)

As one of the kinematical assumptions of Kirchhoff-Love theory is that the director vector $\mathit{a}$ remains orthogonal to the middle surface of the shell, then the rotation tensor may be defined by

$\mathit{Q}={\mathit{e}}_{\mathit{i}}⨂{\mathit{e}}_{\mathit{i}}^{\mathit{r}}$ (4)

The local orthogonal coordinate system in the current configuration is expressed by ${\mathit{e}}_{\mathit{i}}$ and its equation is defined in further section(See equation (30)).

The nonlinear deformation map that takes a vector in the reference configuration and takes it to the current configuration at any positive time $t$ is defined by ${\phi }_t:\mathrm{\xi }\to \mathrm{x}$and the deformation gradient $\mathit{F}=Grad{\mathbf{\phi }}_{\mathrm{t}}\left(\mathbf{\xi }\right)$ is defined by

$F=\partial x/\partial \xi ,F=\frac{\partial \left(z+Q{a}^r\right)}{\partial {\xi }_{\alpha }}\otimes e_{\alpha }^r+\frac{\partial \left(z+Q{a}^r\right)}{\partial {\xi }_3}\otimes e_3^{r}andF=f_{\alpha }\otimes e_{\alpha }^r+f_3\otimes e_3^r,$ (5)

where the vectors ${\mathit{f}}_i$ are defined by

$f_{\alpha }=z_{,\alpha }+Q_{,\alpha }{a}^{r}and{f}_3=Q{e}_3^r=e_3.$ (6)

The Curvature vectors are defined by

$K_{\alpha }=Q_{,\alpha }{Q}^{T}and{\kappa }_{\alpha }=axial\left(K_{\alpha }\right),$ (7)

where ${\mathit{\kappa }}_{\mathit{\alpha }}$ can yet be rewritten as

$\begin{array}{l}{\kappa }_{\alpha }={\Gamma }_{\beta }{u}_{,\beta \alpha },\\ {\Gamma }_1={\left(e_1\cdot z_{\mathrm{,1}}\right)}^{-1}\left[Skew\left(e_1\right)-\left(e_1\cdot z_{\mathrm{,2}}\right){\left(e_2\cdot z_{\mathrm{,2}}\right)}^{-1}\left(e_1\otimes e_3\right)\right]and\\ {\Gamma }_2={\left(e_2\cdot z_{\mathrm{,2}}\right)}^{-1}\left(e_1\otimes e_3\right).\end{array}$ (8)

Consequently, one may write the components of deformation gradient as

${\mathit{f}}_{\mathit{\alpha }}={\mathit{z}}_{,\mathit{\alpha }}+{\mathit{\kappa }}_{\mathit{\alpha }}\times \mathit{a}$ (9)

and the Jacobian and cofactor of $\mathit{F}$ may be written as

$J=\det F=f_1\cdot \left(f_2\times f_3\right),CofF=J{F}^{-T}=g_i\otimes e_{i}and\begin{array}{c}{g}_1=f_2\times f_3\\ g_2=f_3\times f_1\\ g_3=f_1\times f_2\end{array}$ (10)

Back-rotated counterparts of the deformation gradient and of strains may also be defined by

$F^r=Q^{T}F=I+{\gamma }_{\alpha }^r\otimes e_{\alpha }^r+{\gamma }_{33}^r\otimes e_3^r,{\gamma }_{\alpha }^r={\eta }_{\alpha }^r+k_{\alpha }^r\times a^{r}and\begin{array}{c}{\eta }_{\alpha }^r=Q^T{z}_{,\alpha }-e_{\alpha }^r\\ {\kappa }_{\alpha }^r=axial\left(Q^T{Q}_{,\alpha }\right)\end{array}$ (11)

The plane stress condition is forced considering vanishing stresses in the normal direction of the mid-plane by the expression

$\left(\mathbf{\tau }{\mathit{e}}_3\right)\cdot {\mathit{e}}_3=0,$ (12)

where, $\mathrm{\tau }$ is the decomposition of the first Piola-Kirchhoff stress tensor on the cartesian axes. It is defined by

$P=\frac{\partial \psi }{\partial F}=:{\partial }_F\psi andP={\tau }_i\otimes e_i^r,\tau =\frac{\partial \psi }{\partial f_i}.$ (13)

In the previous equation, $\psi =\psi \left(\mathit{F}\right)$ is the so-called Helmholtz free energy per unit of reference volume.

3.2 Weak form of equilibrium

In solid mechanics, the balance of linear momentum in statics (LAI et al., 2009) may be defined by the differential equation (14) as

$Div\mathit{P}+{\rho }_0\left(\stackrel{-}{\mathit{b}}\right)=0$ (14)

However, classically, the mathematical model implemented on FEM code is based on the weak form of the equilibrium equation. The weak form (in this case, variational formulation) can be obtained applying the principle of virtual work (Wriggers, 2008) (Viebahn et al., 2017) and consequentily we have

$\begin{array}{l}\delta W=\delta W_{\mathrm{int}}-\delta W_{ext}=0,\forall \delta uand\delta W_{\mathrm{int}}=\underset{{\rm B}}{{\displaystyle \int }}P:\delta FdV\\ \delta W_{ext}=\underset{\partial {\rm B}}{{\displaystyle \int }}\overline{t}\cdot \delta xdA+\underset{{\rm B}}{{\displaystyle \int }}\overline{f}\cdot \delta xdV,\end{array}$ (15)

where:

• ${\mathrm{\delta }\mathrm{w}}_{ext}$ is virtual work of external forces,

• ${\mathrm{\delta }\mathrm{w}}_{int}$ is the virtual work of internal forces,

• $\mathit{P}$ is the first Piola-Kirchhoff stress tensor

• $\mathrm{\delta }\mathit{F}$ is the virtual deformation gradient (funcition of virtual displacements)

• $\mathit{u}$ is the displacement; $\mathrm{\delta }\mathit{u}$ is the virtual displacement.

• $\overline{\mathit{t}}$ is the boundary forces at Domain’s boundary.

• $\overline{\mathit{f}}$ is the volumetric external force.

The internal virtual work at Equation (15) may be further developed (Pimenta et al., 2010) by

$\delta W_{int}={\int }_{{\mathrm{\Omega }}^r}{(\mathit{\sigma }}_{\mathit{\alpha }}^{\mathit{r}}\cdot {\delta \mathit{\epsilon }}_{\mathit{\alpha }}^{\mathit{r}})d{\mathrm{\Omega }}^r,$ (16)

where

$n_{\alpha }^r=\underset{H}{{\displaystyle \int }}{\tau }_{\alpha }^{r}dHand{m}_{\alpha }^r=\underset{H}{{\displaystyle \int }}(a^r\times {\tau }_{\alpha }^r)dH,$ (17)

with

• $H=\left[-\frac{h}{2},\frac{h}{2}\right]$

• ${\mathit{n}}_{\alpha }^r=$ Back-rotated Forces per unit length

• ${\mathit{m}}_{\alpha }^r=$ Back-rotated moments per unit length

• $\stackrel{-}{M}$ = Inertia property of cross section

• ${\mathit{\kappa }}_{\alpha }^r$= Back-rotated curvature vector

• ${\mathit{\eta }}_{\alpha }^r$ = Back-rotated membrane strains

For convenience, ${\mathit{\sigma }}_{\mathit{\alpha }}^{\mathit{r}}$ and ${\mathit{\epsilon }}_{\mathit{\alpha }}^{\mathit{r}}$ are defined by

${\sigma }_{\alpha }^r={\left[\begin{array}{cc}{n}_{\alpha }^r& m_{\alpha }^r\end{array}\right]}^T,{\epsilon }_{\alpha }^r={\left[\begin{array}{cc}{\eta }_{\alpha }^r& {\kappa }_{\alpha }^r\end{array}\right]}^T$ (18)

and the external virtual work at Equation (15) may also be rewritten as

$\delta W_{ext}={\int }_{{\mathrm{\Omega }}^r}\left({\stackrel{-}{\mathit{t}}}^{\mathit{t}}+{\stackrel{-}{\mathit{t}}}^{\mathit{b}}+{\int }_{H^r}\stackrel{-}{f}d{H}^r\right)\cdot \delta \mathit{x}d{\mathrm{\Omega }}^r+{\int }_{{\partial \mathrm{\Omega }}^r}{\int }_{H^r}{\stackrel{-}{\mathit{t}}}^{\mathit{l}}d{H}^{r}d{\partial \mathrm{\Omega }}^r,$ (19)

where has been considered

${\int }_{\partial B}\stackrel{-}{\mathit{t}}dA={\int }_{{\mathrm{\Omega }}^r}{(\stackrel{-}{\mathit{t}}}^{\mathit{t}}+{\stackrel{-}{\mathit{t}}}^{\mathit{b}})d{\mathrm{\Omega }}^r+{\int }_{{\partial \mathrm{\Omega }}^r}{\int }_{H^r}{\stackrel{-}{\mathit{t}}}^{\mathit{l}}d{H}^{r}d{\partial \mathrm{\Omega }}^r,$ (20)

with

• $\stackrel{-}{\mathit{t}}$ = Surface traction

• ${\stackrel{-}{\mathit{t}}}^{\mathit{t}}$ = Top Surface traction

• ${\stackrel{-}{\mathit{t}}}^{\mathit{b}}$ = Bottom Surface traction

• ${\stackrel{-}{\mathit{t}}}^{\mathit{l}}$ = Lateral Surface traction

• $\overline{\mathit{f}}=$Vector of ext. volumetric body forces.

3.3 Constitutive equations (Elasticity)

The model implemented in this article uses an elastic material model (neo-Hookean material) suitable for large displacements and deformations (Pimenta, P.M, 2006). The strain energy is given by equation (21)

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

where $I_i$ are the principal invariants of the deformation Cauchy-Green right tensor “C”, defined by

$I_1=trC=f_i\cdot f_i,I_2=tr\left[CofC\right]=g_i\cdot g_i,I_3=\det C=J^2={\left(f_1\cdot \left(f_2\times f_3\right)\right)}^2,$ (22)

with:

• $\lambda$ and $\mathrm{\mu }$ - Lamé coefficients (Material properties)

• $\psi =\psi \left(\mathit{F}\right)$ = Helmholtz free energy

• $\mathit{C}={\mathit{F}}^T\mathit{F}=$deformation Cauchy-Green (right) tensor.

3.4 Spatial Discretization

The finite element developed is triangular. Because of the Kirchhoff-Love shell theory, it would be required $C^1$ continuous approximations in order to preserve continuity of the displacement gradient among all domain. This model, however, imposes $C^1$ displacement field inside the element and enforces rotation between adjacent elements. Emforcement of kinking angle between adjacent elements is better explained in further section.

In the current model, shape functions are based on Area coordinates (Papadopoulos P., 2015). The interpolations for position vector (reference configuration), displacement vector and derivatives are defined by equations (23) and (24), which are similarly represented in Viebahn et al. (2017).

$\mathit{\zeta }\simeq {\mathit{\zeta }}^h={\sum }_I^{nen}{N}_I{\mathit{\zeta }}_I,$ $\mathit{u}\simeq {\mathit{u}}^h={\sum }_I^{nen}{N}_I{\mathit{d}}_{I}and{{\mathit{u}}_{,\alpha }\simeq \mathit{u}}_{,\alpha }^h={\sum }_I^{nen}{N}_{I,\alpha }{d}_I,$ ${\mathit{u}}_{,\alpha \beta }\simeq {\mathit{u}}_{,\alpha \beta }^h={\sum }_I^{nen}{N}_{I,\alpha \beta }{d}_I.$ (23)

Expressions (24) represents the discretization of the variation and linearization of displacements and spatial derivatives

${\mathrm{\Delta }\mathit{u}}^h={\sum }_I^{nen}{N}_I{\mathrm{\Delta }\mathit{d}}_I,$ ${\delta \mathit{u}}^h={\sum }_I^{nen}{N}_I{\delta \mathit{d}}_I,{\mathrm{\Delta }\mathit{u}}_{,\alpha }^h={\sum }_I^{nen}{N}_{I,\alpha }{\mathrm{\Delta }\mathit{d}}_I,$ ${\delta \mathit{u}}_{,\alpha }^h={\sum }_I^{nen}{N}_{I,\alpha }{\delta \mathit{d}}_{I}and{\mathrm{\Delta }\mathit{u}}_{,\alpha \beta }^h={\sum }_I^{nen}{N}_{I,\alpha \beta }{\mathrm{\Delta }\mathit{d}}_I,$ ${\delta \mathit{u}}_{,\alpha \beta }^h={\sum }_I^{nen}{N}_{I,\alpha \beta }{\delta \mathit{d}}_I.$ (24)

It is important to remark that the current Paper does not focus time-dependent problems. Further development in this topic will be approached in future paper under development.

3.5 Finite element definition

Figure 2 illustrates the finite element implemented in this article with its degrees of freedom in its 6 nodes. It is assumed a quadratic displacement field (with nodal values of $u^{\left(1\right)}..u^{\left(6\right)}$) and 3 independent rotation parameters (${\phi }_{\mathrm{\Delta }}^{\left(4\right)},{\phi }_{\mathrm{\Delta }}^{\left(5\right)},{\phi }_{\mathrm{\Delta }}^{\left(6\right)}$).

Figure 2
K-L shell finite element.

3.6 Kinking angle enforcement

The kinking angle enforcement in the present article follows a methodology similar but different from Costa e Silva et al. (2018) and Viebahn et al. (2017).

Figure 3
Enforcement of C1-continuity; local coordinate system in reference (a) and current (b) configuration. Source: Viebahn et al. (2017)

The shell kinematics based on Kirchhoff-Love model imposes the necessity of continuity of the displacement gradient (See equations (8),(16) and (18)). Consequently, the finite element construction should guarantee $C^1$ continuity in the whole domain, include inside and the edges of the elements. The finite element structure implemented in this article penalizes the kinking angle(See Figure 3) between adjacent elements based on one of the so-called Rodrigues rotation parameters developed in Costa e Silva (2020).

As explained in Viebahn et al. (2017), the enforcement of C1 continuity at element edges is expected to be satisfied asymptotically as $h\to 0$ (mesh refinement) if one guarantees that the kinking angle(at the midpoint) between adjacent elements does not change during the deformation. Geometrically, the kinking may be understood as the angle between mid-surface normal vectors of adjacent elements.

As seen in equation (2), the tensor $\mathit{Q}$ rotates the normal vector of the shell from reference to current configuration. As will be shown in further equations, if one imposes that the rotation tensor $\mathit{Q}$ (or its incremental part) is the same for both adjacents elements, one may guarantee the continuity of the kinking angle. In this sense, the author uses one of the characteristics of the Rodrigues rotation parameters, in which the tensor $\mathit{Q}$ may be written as function of $\mathbf{\alpha }$:

$\mathrm{\alpha }=2\mathrm{t}\mathrm{a}\mathrm{n}\left(\mathrm{\theta }/2\right),\mathbf{\alpha }=\mathrm{\alpha }\mathbf{e}\mathbf{},\mathbf{A}=\mathrm{s}\mathrm{k}\mathrm{e}\mathrm{w}\left(\mathbf{\alpha }\right),$ (25)

$\widehat{\mathbf{Q}}\left(\mathbf{\alpha }\right)={\left(\mathbf{I}-\frac{1}{2}\mathbf{A}\right)}^{-1}\left(\mathbf{I}+\frac{1}{2}\mathbf{A}\right)=\mathbf{I}+\frac{4}{4+{\mathrm{\alpha }}^2}\left(\mathbf{A}+\frac{1}{2}{\mathbf{A}}^2\right).$

When working with Rodrigues parameter in rotations, rotation increment is easily treated by equation

$Q_{i+1}=Q_{\Delta }{Q}_i,Q_{i+1}=\widehat{Q}\left({\alpha }_{i+1}\right),Q_{\Delta }=\widehat{Q}\left({\alpha }_{\Delta }\right),Q_i=\widehat{Q}\left({\alpha }_i\right),$ (26)

where:

${\mathbf{\alpha }}_{\mathbf{i}+1}=\frac{4}{4-{\mathbf{\alpha }}_{\mathbf{i}}\cdot {\mathbf{\alpha }}_{∆}}\left({\mathbf{\alpha }}_{\mathbf{i}}+{\mathbf{\alpha }}_{∆}-\frac{1}{2}{\mathbf{\alpha }}_{\mathbf{i}}\times {\mathbf{\alpha }}_{∆}\right).$ (27)

With the previous expressions, one may write the coordinate system in sequential time steps ($\{{\mathit{e}}_1^i,{\mathit{e}}_2^i,{\mathit{e}}_3^i\}$ at $t_i$ and $\{{\mathit{e}}_1^{i+1},{\mathit{e}}_2^{i+1},{\mathit{e}}_3^{i+1}\}$ at $t_{i+1}$) by

$e_3^{i+1}=Q_{\Delta }{e}_3^i,e_3^{i+1}-e_3^i={\alpha }_{\Delta }\times e_3^{m}and{e}_3^m=\frac{1}{2}\left(e_3^{i+1}+e_3^i\right).$ (28)

In (Costa e Silva, 2020), two important expressions are achieved regarding rotation parameters. It is shown that

${\alpha }_{\Delta }={\Vert e_3^m\Vert }^{-2}\left(e_3^i\times e_3^{i+1}\right)+{\phi }_{\Delta }{\Vert e_3^m\Vert }^{-1}{e}_3^{m}and{\phi }_{\Delta }={\Vert e_3^m\Vert }^{-1}{\alpha }_{\Delta }\cdot e_3^m.$ (29)

The previous equations shows that the normal vector ${\mathit{e}}_3$ at the edge of the element at time $t_{i+1}$ may be find if the parameter ${\phi }_{\mathrm{Ä}}$ is known, but at the same time one may determine ${\mathit{e}}_3$ with the displacement field and so far no condition has been imposed to guarantee that both normal vectors are the same. In this sense, the penalty approach is used to guarantee that the displacement field is compatible with the ${\phi }_{\mathrm{Ä}}$ independent parameters, witch have been introduced, before, as a mechanism to guarantee $C^1$ continuity with sufficient mesh refinement. The following equations are used

$e_1={\Vert z_{\mathrm{,1}}\Vert }^{-1}{z}_{\mathrm{,1}},e_3={\Vert z_{\mathrm{,1}}\times z_{\mathrm{,2}}\Vert }^{-1}\left(z_{\mathrm{,1}}\times z_{\mathrm{,2}}\right),e_2=e_3\times e_1,$ (30)

${\mathit{u}}_{,1}^{\left(4\right)}=\frac{1}{l_3}\left({\mathit{u}}^{\left(2\right)}-{\mathit{u}}^{\left(1\right)}\right),{\mathit{u}}_{,1}^{\left(5\right)}=\frac{1}{l_1}\left({\mathit{u}}^{\left(3\right)}-{\mathit{u}}^{\left(2\right)}\right),{\mathit{u}}_{,1}^{\left(6\right)}=\frac{1}{l_2}\left({\mathit{u}}^{\left(1\right)}-{\mathit{u}}^{\left(3\right)}\right),$ (31)

${\mathit{u}}_{,2}^{\left(4\right)}=\left(\frac{C_1}{2A}\right){\mathit{u}}^{\left(1\right)}+\left(\frac{C_2}{2A}\right){\mathit{u}}^{\left(2\right)}+\left(\frac{C_3}{2A}\right){\mathit{u}}^{\left(3\right)}+\left(\frac{-C_3}{A}\right){\mathit{u}}^{\left(4\right)}+\left(\frac{C_3}{A}\right){\mathit{u}}^{\left(5\right)}+\left(\frac{C_3}{A}\right){\mathit{u}}^{\left(6\right)},$ $where{C}_1=x_3-x_2,C_2=x_1-x_3,C_3=x_2-x_1.$ (32)

With the previous equations, one may write ${\mathit{\alpha }}_{\mathit{\Delta }}$ based on the displacement field only and consequently ${\varphi }_{\mathrm{\Delta }}$ may be defined by

${\alpha }_{\Delta }=\frac{2\left(e_j^i\times e_j^{i+1}\right)}{1+e_k^i\cdot e_k^{i+1}},{\varphi }_{\text{Δ}}={\alpha }_{\Delta }\cdot \frac{e_1^m}{\Vert e_1^m\Vert }and{e}_1^m=\frac{1}{2}\left(e_1^i+e_1^{i+1}\right).$ (33)

Here, the penalty is summed to Helmholtz free energy to create the penalty function:

${\Pi }^{pen}=\psi +\frac{1}{2}k{\left({\phi }_{\Delta }^{\left(4\right)}-{\varphi }_{\Delta }^{\left(4\right)}\right)}^2+\frac{1}{2}k{\left({\phi }_{\Delta }^{\left(5\right)}-{\varphi }_{\Delta }^{\left(5\right)}\right)}^2+\frac{1}{2}k{\left({\phi }_{\Delta }^{\left(6\right)}-{\varphi }_{\Delta }^{\left(6\right)}\right)}^2.$ (34)

As a note, it is important to remark that it is used similar letters in this paper (${\phi }_{\Delta }^{}$ and ${\varphi }_{\Delta }^{}$) because the represent the same parameter, however, as explainded before, they are calculated in different ways and are forced to be the same by equation (34).

In the previous equation, it is used the only artificial numerical factor of the formulation (“k”). In the examples, it is assumed as a multiple of the bending stiffness ($E{h}^3/\left(12\left(1-{\upsilon }^2\right)\right)$). Figure 4 illustrates adjacent elements with its degrees of freedom (DoF).

Figure 4
Adjacent elements.

4 RESULTS

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

4.1 Square Plate

The square plate is simply supported among all 4 sides and an external normal force is applied in the center of the plate. This is a classical problem in theory of elasticity and analytical results may be found in bibliography for small displacements. The formulas used in the simulations have been refered to Timoshenko and Krieger (1959) and Young et al. (2002). The simulation has been executed for different mesh discretization and for different Thickness and its output is compared to analytical results and to commercial software Nastran Triangular Shell Element - Ctria6 (See Autodesk Inc, 2019).

In the simulations and calculations, mechanical properties, dimensions and loads are defined by

• $\mathrm{E}=200$ GPa;

• $\nu =\mathrm{0,3}$;

• Center down Force = 10 N;

• a(width) =1 m;

• b(length) = 1 m;

• h(thickness) = 1mm and 0,2mm;

Figure 5 illustrates the simulated model.

Figure 5
FEM model. rectangular plate with concentrated normal center force.

Figure 6, Figure 7 and Figure 8 illustrate the main results of this first simulated model. The first (Figure 6) plots the total displacement of the center of the plate for different load factors. Each line represents different models used to simulate. The second (Figure 7) represents tha same simulation but with different thickiness, the last graph (Figure 8) in this subsection compares maximum displacement (at center) for different mesh discretization.

Figure 6
Square Plate - Thk=1mm.

Figure 7
Square Plate - Thk=0,2mm.

Figure 8
Square Plate, different mesh refinements - Thk=0,2mm.

In these simulations it has been observed that for K-L Lagrange model and Nastran Ctria6 the simulations took longer time to converge with finer meshes. In addition, it has been noticed that the K-L Penalty based element has performed the faster simulations (compared to the other models with same parameters). Another interesting fact is that the author could not simulate successfully the problem using 3D solid tetrahedral Nastran elements for given thickness. In this case, the software was unable to solve the numerical problem and presented several warnings and errors. This was expected due to the locking phenomena discussed in section 1 of this paper.

Figure 6 demonstrate great results agreement between different models. For 1mm thickness the model developed in this article had practically the same results as the model developed by Costa e Silva (2020) and to the commercial software Nastran, with differences smaller than 5%. The analytical solution could only approximate the results at the beginning of the curve due to its linearized form. For the thinner(0,2mm) plate (Figure 7) it may be seen that the models present some slight differences. The (vertical) displacements at plate center are very similar, however Nastran FEM have shown as the stiffer result. For a load factor of 1, the displacements are less than 10% different from one another. Here again the analytical formula is suitable only for very low displacements. Figure 8 shows that all the models tend to converge to a maximum value when using a bigger number of elements, as it would be expected. However, in the graphic the model developed in this article (KL T63i Pen) have presented the fastest convergence in this simulation. Converging to the more accurate result with fewer elements is a desirable important feature for a finite element model.

4.2 Cantilever Beam

A Cantilever Beam is subject to a transverse force at the free end. The simulation has been executed for different mesh discretization and for different thicknesses in order to analyze the FEM model built in this article. The results are compared to the Lagrange shell FEM model from Costa e Silva (2020) and to commercial software (Nastran Triangular Shell Element - CTRIA6).

Mechanical properties and dimensions of the beam that have been used are displayed below.

• $\mathrm{E}=210$ GPa;

• $\nu =\mathrm{0,3}$;

• L (beam length) =2400 mm;

• h(height) = 100 mm;

• b(width of cross section) = 11,64 mm;

The beam is subject to vertical down force up to $F_y=10000\mathrm{N}$ and an horizontal perturbation force is also introduced ($F_z=F_y\cdot {10}^{-4}$) in order to induce the buckling phenomena.

Figure 9 illustrates the simulated model and Figure 10 shows FEM model.

Figure 9
Cantilever Beam. Image Source: Costa e Silva (2020)

Figure 10
Lateral view, Cantilever Beam.

The analytical results (Linear elastic) for this problem may be found in Solid mechanics books and is illustrated by the formula

$y_{\mathrm{m}ax}=-F_y{L}^3/3EIandI=b{h}^3/12,$ (35)

where

• $y_{\mathrm{m}ax}$ is the maximum vertical displacement at the center of the free side of the beam.

• $F_y$ is the concentrated force at the center of the free side of the beam.

• $L$ is the length of the beam

• E is the Young modulus of the material.

• $I$ is the moment of inertia of the rectangular cross section

• b is the thickness of the cross section.

• h is the height of the cross section.

Figure 11 and Figure 12 illustrate the vertical and lateral displacement of the free side of the beam for different forces. Each line represents different models. For certain high transversal forces, it is observed flexural buckling of the structure.

Figure 11
Vertical maximum displacement vs load factor.

Figure 12
Lateral maximum displacement vs load factor.

The results shown above demonstrates that both Lagrange and penalty methods led to considerable similar (difference smaller than 5% from one another) and consistent results, once they could capture the flexural buckling of the beam. Both methods have close results to the analytical model only for low load factors due to the Linear elastic limitation of the analytical model. For high loads, the beam buckles and the analytical results loses its physical meaning. For this simulation, the commercial software (Nastran Triangular Shell Element - CTRIA6) was unable to converge for load factor bigger than 0,24, presenting numerical errors and warnings. This load factor is approximately the step that the beam starts to buckle. The simulations were accomplished with approximately 160 elements and in Nastran, it has been used the non-linear material setup with the same mechanical properties as described above. The capacity to successfully capture the buckling phenomena is an important and desirable feature that was observed in this simulation.

4.2 Pinched Cylinder

This example is reproduced in other articles (Ivannikov et al., 2015) (Campello et al, 2003) (Sansour and Kollmann, 2000) (Costa e Silva, 2020) and is a good benchmark to evaluate the shell model for big displacements and curvatures. It consists of a cylinder with rigid ends and subject to a radial vertical force in its center. The Figure 13 illustrates the simulated model and Figure 14 FEM simulation.

Figure 13
Pinched cylinder - Scheme of boundary value problem. Image Source: Viebahn et al (2017)

Figure 14
Pinched Cylinder - FEM simulation (current model).

In the simulations, it has been used the same parameters of the bibliography, then one can compare results. The mechanical properties, dimensions and loads, that are used in the simulation, are

• $\mathrm{E}=$$3\cdot {10}^{10}$ Pa

• $\nu =\mathrm{0,3}$

• $\mathit{F}\left(baseforce\right)$ = 1000 N

• L (Cylinder length) =200 mm;

• R(Cylinder Radius) = 100 mm;

• b(thickness) =1 mm;

Figure 15 illustrates the vertical displacement of point A and Figure 16 plots the horizontal displacement of point B for different load factors. Each line in the graphics represents a different Finite element model. Results from some sources in bibliography have been plotted in the same Graph for comparison.

Figure 15
Maximum vertical displacement at A

Figure 16
Maximum Horizontal displacement at B

Figure 15 and Figure 16 illustrates that all models have similar results for low loading, however for bigger forces, and big displacements, Nastran CTRIA6 has different and unstable results. Yet in the same figures, the model developed in this article has demonstrated very close results to the lagrange based model (Costa e Silva, 2020) with results differing less than 10% from one another even for higher loads. Nastran CTRIA6 presented non-smooth results regarding the displacements and presented some warnings in the numerical simulation regarding numerical convergence quality. In the simulations, it has been used about 5000 shell elements and in Nastran it has been used non-linear material setup with the same mechanical properties as specified above.

5 CONCLUSIONS

The present article develops analytically a very simple finite element and presents results in 3 different simulations comparing it to other FEM models and to linear analytical models. All the kinematics is based on Kirshoff-Love hyphothesis and the C1 continuity expected between adjacent elements because of the K-L hypothesis is expect to be fulfilled by applying penalty method internally to the element enforcing agreement between a rotation degree of freedom (that is shared between adjacent elements) and the displacement field (defined by others DoFs). One of the advantages of the current model compared to others K-L based elements is that a second mesh for enforcing penalty or lagrange multiplyers is not necessary. All the approximation is accomplished internally to the element.

The promising results obtained so far are coherent and satisfactory, but more tests are going to be made in further studies to have a stronger idea of the positives/negatives aspects of the model. An important conclusion of these first simulations is that all the simulations presented in this article shows compatible and reliable results regarding the model developed. Some differences, however, are observed. The results obtained with commercial software (Nastran CTRIA6) were observed as more stiffer than the K-L models (Penalty and Lagrange based models) in the first and third example (Square plate and Pinched cylinder) and it is also unstable for large displacements (Second and third Examples). It is also possible to check by Figure 8 that the K-L Penalty based model has shown better convergence with fewer elements. It is important to remark as well that the analitical models in all examples could not depict the non-linear phenomena that is captured with the non-linear Finite element and consequently presents similar results only for small displacements due to its linear characteristic.

For the next studies, the authors plan also to develop a shell element with out the use of penalty or Lagrange multipliers as a mechanism to guarantee C1 continuity. The theory behind this new element is still under developement and more information is going to be presented in forthcoming work. It is believed that a very reliable and consistent thiangular shell finite element as presented in this article combined with powerfull mesh generators is a great option for simulations of thin complex structures.

6 ACKNOWLEDGEMENTS

Matheus L. Sanchez thankfully acknowledges Clark Solutions, where he works as a project engineer providing financial support, knowledge, applicability and industrially experience. C. Costa e Silva gratefully acknowledges the Federal Institute of Science and Technology Education of São Paulo for financial support. P. M. Pimenta acknowledges the support by CNPq under the grant 308142/2018-7 and the Alexander von Humboldt Foundation for the Georg Forster Award that made possible his stays at the Universities of Duisburg-Essen and Hannover in Germany as well as the French and Brazilian Governments for the Chair CAPES-Sorbonne that made possible his stay at Sorbonne Universités on a leave from the University of São Paulo.

References

• Autodesk, Inc (2019, March 11) CTRIA6 Autodesk Nastran. https://knowledge.autodesk.com/support/nastran/learn-explore/caas/CloudHelp/cloudhelp/2020/ENU/NSTRN-Reference/files/GUID-956704D9-26C5-4FF8-85F2-955F383C5ECD-htm.html
  » https://knowledge.autodesk.com/support/nastran/learn-explore/caas/CloudHelp/cloudhelp/2020/ENU/NSTRN-Reference/files/GUID-956704D9-26C5-4FF8-85F2-955F383C5ECD-htm.html
• Brezzi, F.; Fortin, M. (1991). Mixed and hybrid finite element methods. New York: Springer-Verlag.
• Campello, E.M.B.; Pimenta, P.M.; Wriggers, P. (2003). A triangular finite shell element based on a fully nonlinear shell formulation. Computational Mechanics. Comput. Methods Appl. Mech. Engrg
• Costa e Silva, C., Maassen, S., Viebahn, N., Pimenta, P.M., Schröder, J. (2018). Connection Between Bernoulli-Euler Rods And Kirchhoff-Love Shells. Cilamce Congress.
• Costa e Silva, C., Maassen, S., Viebahn, N., Pimenta, P.M., Schröder, J (2020). On simultaneous use of geometrically exact shear-rigid rod and shell finite elements. Submitted to Computational Mechanics.
• Costa e Silva, C. (2020) Geometrically exact shear-rigid shell and rod models. Ph.D. Thesis, University of Sao Paulo, Brazil
• Ivannikov, V., Tiago, C., Pimenta, P.M. (2015). Generalization of the C₁ TUBA plate finite elements to the geometrically exact Kirchhoff-Love shell model. Comput Methods Appl Mech Eng
• Korelc, J.A., Wriggers, P. (2016). Automation of Finite Element Methods. [S.l.]: Springer.
• LAI, W. Michael et al. (2009) Introduction to continuum mechanics. Butterworth-Heinemann.
• Nigro, P.S.B., Simões, E.T., Pimenta, P.M. ; Schröder, J. (2019) Model Order Reduction with Galerkin Projection Applied to Nonlinear Optimization with Infeasible Primal-Dual Interior Point Method. Int J Numer Methods Eng, v. 1, nme.6181
• Papadopoulos P. (2015) Introduction to the Finite Element Method. University of California, Berkeley
• Pimenta, P.M, Campello, E.M.B., Wriggers, P. (2004). A fully nonlinear multi-parameter shell model with thickness variation and a triangular shell finite element.
• Pimenta, P.M (2006) Fundamentos da Mecânica dos Sólidos e das Estruturas. São Paulo.
• Pimenta, P.M, Campello, E.M.B. (2009) Shell curvatures as an initial deformation: a geometrically exact finite element approach. Int J Numer Methods Eng 78:1094-1112.
• Pimenta, P.M., Neto, E.S.A, Campello, E.M.B. (2010). Fully Nonlinear Thin Shell Model of Kirchhoff-Love type.
• Pimenta, P.M. (1993). On the geometrically-exact finite-strain shell. Proceeding of the third Pan-American congress on.
• Reddy, J. N. (2006) Theory and analysis of elastic plates and shells. CRC press.
• Sansour, C., Kollmann, F. (2000) Families of 4-node and 9-node finite elements for a finite deformation shell theory. An assessment of hybrid stress, hybrid strain and enhanced strain elements. Comput Mech 24:435-447
• Simo, J.C., Fox, D.D. (1989). On a stress resultant geometrically exact. Comput Methods Appl Mech.
• Timoshenko, S. P., Krieger, S.W. (1959). Theory of plates and shells. [S.l.]: McGraw-hill, 1959.
• Viebahn, N.; Pimenta, P.M.; Schröder, J.A. (2017). A simple triangular finite element for nonlinear thin shells: statics, dynamics and anisotropy. Computational Mechanics.
• Wriggers, P. (2008). Nonlinear finite element methods. Springer Science & Business Media.
• Young, W.C., Budynas, R.G., Sadegh, A.M. (2002). Roark's formulas for stress and strain. McGraw-Hill New York.

• SYMBOLS AND ABBREVIATIONS

  As a remark note, throughout the text italic Greek or Latin lowercase letters ($a,b,\dots ,\alpha ,\beta ,\dots$) represents scalar quantities, bold Greek or Latin lowercase letters ($\mathit{a},\mathit{b},\dots ,\mathbf{\alpha },\mathbf{\beta },...$) indicate vectors and bold Greek or Latin capital letters ( A,B, . . .) signify second-order tensors. In summations, Greek indices indicates range from 1 to 2, Latin indices indicate from 1 to 3. ${\mathbb{R}}^3$ denotes 3-dimension Euclidian space.

Edited by

Editor:

Marco L. Bittencourt.

Publication Dates

• Publication in this collection
  11 Dec 2020
• Date of issue
  2020

History

• Received
  17 May 2020
• Reviewed
  20 Sept 2020
• Accepted
  30 Sept 2020
• Published
  06 Oct 2020