# Unconstrained Finite Element for Geometrical Nonlinear Dynamics of Shells

Authors: Humberto Breves Coda; Rodrigo Ribeiro Paccola
Journal: Mathematical Problems in Engineering
 (2009)

---

## Abstract

This paper presents a positional FEM formulation to deal with geometrical nonlinear dynamics of shells. The main objective is to develop a new FEM methodology based on the minimum potential energy theorem written regarding nodal positions and generalized unconstrained vectors not displacements and rotations. These characteristics are the novelty of the present work and avoid the use of large rotation approximations. A nondimensional auxiliary coordinate system is created, and the change of configuration function is written following two independent mappings from which the strain energy function is derived. This methodology is called positional and, as far as the authors' knowledge goes, is a new procedure to approximated geometrical nonlinear structures. In this paper a proof for the linear and angular momentum conservation property of the Newmarkβ algorithm is provided for total Lagrangian description. The proposed shell element is locking free for elastic stress-strain relations due to the presence of linear strain variation along the shell thickness. The curved, high-order element together with an implicit procedure to solve nonlinear equations guarantees precision in calculations. The momentum conserving, the locking free behavior, and the frame invariance of the adopted mapping are numerically confirmed by examples.


---

## Body

## 1. Introduction

An accurate analysis of structures that exhibit large deflections is of great importance for structural design. The increasing search for economy and optimal material application leads to the conception of very flexible structures. As a consequence, the equilibrium analysis in the  nondeformed position is no more acceptable for most applications. This affirmation is confirmed by the large amount of research regarding this subject. One can see pioneering studies related to nonlinear analysis of structures as the works of Bisshopp and Drucker [1], Mattiasson [2], Goto et al. [3], Jenkins et al. [4], Kerr [5], Mondkar and Powell [6], Belytschko et al. [7], among others.Moreover, some structures are naturally geometrically nonlinear as balloons, airbags, cables, membranes and so on. The design of this kind of structures requires more sophisticated theories than the linear ones. One can see, for instance, the works of Stein and Hedgepath [8], Baginski et al. [9], Pipkin [10], Bonet et al. [11] among others.This study is concerned with the development of a new Finite Element methodology to solve geometrical nonlinear dynamics of shells. In order to achieve a robust formulation, the resulting element should be free of shear and volumetric locking. This problem is solved here by the natural presence of the transverse shear strain in the proposed kinematics. The novelty of the proposed formulation is the use of positions and generalized unconstrained vector mapping, resulting in a naturally objective continuum representation of the shell, free of large rotation descriptions and locking. There exists another kind of rotation-free elements as proposed in the works of Oñate and Zárate [12] and Flores and Oñate [13], developed specifically for thin shells and based on curvature considerations.There are different approaches regarding time integration for transient nonlinear FEM dynamics, one can mention, as an example, three different approaches. The first is the explicit time integration for fast solutions, analyzed in details by Argyris et al. [14] and works cited therein, where small time steps are adopted and an indirect control of errors is made. The second is the so-called variational energy conserving algorithms, necessary for corotational like formulations; see for instance Simo et al. [15]. Finally, the implicit time integration for total Lagrangian formulations as described by Lopez [16].The formulation proposed here is total Lagrangian and, due to its unconstrained vector mapping, it presents constant mass matrix. It is possible to apply the Newmakβ integrator as a momentum conserving algorithm. A simple proof of the momentum conserving property of the Newmark β method is given in this paper. The proof is restricted to total Lagrangian formulation (not extended to corotational formulations) and trivially fulfils the energy conserving property for rigid bodies. All required features of the formulation as: locking free, frame invariance, and momentum conserving (linear and angular) are checked in the numerical examples section.


## 2. Strain Measure and Specific Strain Energy Potential

This section summarizes simple concepts used to derive the proposed formulation. The Green strain tensor is derived directly from the gradient of the change of configuration function, represented by letterA, given as follows:(2.1)Aij=∂fi∂xj,
where fi is the change of configuration function, as depicted in Figure 1, and xj represents variation regarding initial position.Figure 1
Change of configuration.In Figure1dxi and dyi represent an infinitesimal fiber in the initial and current continuum configurations, respectively. Following Ogden [17], the Green strain tensor can be written as(2.2)Eij=12[AkiAkj-δij]=12[Cij-δij],
in which index notation is adopted. The variables Cij and δij are the right Cauchy-Green stretch tensor and the Kroenecker delta, respectively. The following quadratic strain energy per unit of initial volume is adopted:(2.3)ue=12EijCijklEkl,
resulting into a linear elastic constitutive law relating second Piola-Kirchhoff stress and Green strain, usually called Saint-Venant-Kirchhoff elastic law, that is,(2.4)Sij=∂ue∂Eij=CijklEkl.
The elastic tensor is given by(2.5)Cijkl=2Gν1-2νδijδkl+G(δikδjl+δilδjk),
where G is the shear modulus, given by(2.6)G=E2(1+ν)
with E being the well-known Young modulus and ν  Poisson’s ratio.The true stress (Cauchy stress) is achieved directly from the Second Piolla-Kirchhoff stress following simple expressions given by Ogden [17] or Ciarlet [18], for instance. For the sake of completeness it is necessary to recall that the right Cauchy-Green stretch tensor is positive definite, symmetric, and has six independent values.


## 3. Equivalence between Classical and Generalized Unconstrained Mapping of Solids

This section presents the equivalence of a classical finite element solid mapping and its counterpart, the generalized vector mapping. Figure2 shows a solid (following plane stress or plane strain condition) element of quadrangular shape classically mapped from a usual nondimensional coordinate system. Figure 2 also shows the same solid mapped from the same nondimensional coordinate system, this time using generalized unconstrained vector parameters.Figure 2
Classical and generalized unconstrained vector discretizations.Without loss of generality the equivalence is shown for a low-order mapping and in the next section extended to consider high-order interpolations. In order to be complete, let us take the following shape functions:ϕ1=(1/4)(1-ξ1)(1-ξ2), ϕ2=(1/4)(1+ξ1)(1-ξ2), ϕ3=(1/4)(1+ξ1)(1+ξ2), ϕ4=(1/4)(1-ξ1)(1+ξ2),(3.1)φ1=12(1-ξ1),φ2=12(1+ξ1),
where ξi are nondimensional coordinates. The classical continuum mapping is written as(3.2)xi=ϕℓ(ξ1,ξ2)Xℓifori=1,2,ℓ=1,2,3,4,
where xi are the coordinates of any point of the mapped continuum, ϕℓ are the shape functions, and Xℓi are the coordinates of nodes Pℓ named position parameters.A totally similar mapping is given by(3.3)xi=φℓXℓim+Hℓiφℓ2ξ2fori=1,2,ℓ=1,2,
where the mid line nodal coordinates Xℓim and nodal generalized vectors Hℓi are given by(3.4)X1im=X1i+X4i2,X2im=X2i+X3i2,H1i=X4i-X1i,H2i=X3i-X2i.
Expression (3.3) shows that the equivalent vector mapping is done using nonunitary vectors Hℓi as parameters. In order to generalize the formulation one writes the nodal vectors H1i and H2i as functions of their lengths H1 and H2 resulting(3.5)xi=φℓXℓim+H(ℓ)2ξ2Viℓφℓfori=1,2,ℓ=1,2,
where brackets mean that summation is not implied. The first term of (3.5) describes the reference line approximation. The values Viℓ are the generalized nodal vectors. These vectors are not orthogonal to the reference line and could not be unitary, if desired, see (3.3). However, for the initial configuration they are made unitary, but in the current configuration these vectors assume nonunitary values. This feature is the original of the generalized unconstrained vector mapping of solids. Finally, if one adopts constant height (h0), a usual assumption for bars and shells, (3.5) turns into the desired continuum mapping, that is,(3.6)xi=φℓXℓim+h02ξ2Viℓφℓfori=1,2,ℓ=1,2.
This mapping is totally similar to the classical one, given in (3.2) and, consequently, has the same objectiveness, Crisfield and Jelenić [19].


## 4. Improved Finite Element Kinematics

Improving the solid description of Figure2 by a three-dimensional representation of a shell, one can approximate the mid-surface positions of a shell, (see Figure 3), by the following mapping:Figure 3
Mid-surface mapping.(4.1)fim0=xim=ϕℓXℓi,(4.2)fim1=yim=ϕℓYℓi,
where xim is the ith coordinate of a generic point in the mid surface of the shell at initial configuration, Xℓi is the ith coordinate of node ℓ, yim is the ith coordinate of a generic point in the mid surface of the shell at current configuration, and Yℓi is the ith coordinate of node ℓ at current configuration.One can see in Figure3 that fm0 is the positional mapping from the auxiliary space to the initial mid-surface configuration, fm1 is the positional mapping from the auxiliary space to the current mid-surface configuration, fm is the positional mapping from the initial configuration to the current one (not to be written) and the values Am0, Am1, Am are their respective gradients. Expression (4.1) is totally similar to the reference line description of the first term of (3.6). This time, high-order approximations for a surface are used.To complete the shell kinematic description for both initial and current configurations, one realizes that the difference between a point out of the mid-surface, and its corresponding belonging to the mid-surface generates position vectorsv⃗0 or v⃗1, see Figure 4.Figure 4
Position vectors.A general point of the shell can be defined by adding the position vectors to the corresponding mid-surface point, that is,(4.3)xi=xim+vi0,yi=yim+vi1.Following what has been described in Section3, for constant strain regarding ξ3, one writes vi0 and vi1 as functions of nondimensional coordinates, as(4.4)vi0=h02ϕℓVℓi0ξ3,(4.5)vi1=h02ϕℓVℓi1ξ3,
where h0, vi0, vi1 are, respectively, the initial constant thickness of the shell, the normal unit vector to the initial mid-surface and the current generalized vector (not necessarily normal to the mid-surface).To consider linear strain variation regardingξ3 an improved kinematics is required. This is done adopting a new scalar variable called here the linear rate of thickness variation and denoted by letter a. It is not necessary to introduce this new variable for initial configuration, so expression (4.4) does not change, however, expression (4.5) turns into(4.6)vi1=h02ϕℓVℓi1[ξ3+aξ32].
The introduced quadratic term allows linear strain variation along the transverse direction of the shell. The final generalized nodal vectors Vℓi1 are not constrained, so the final thickness of the shell is not the same as the initial one and can be recovered as(4.7)vi1=ϕℓVℓi1,h=h0vi1vi1.
The final mapping is generalized and unconstrained, like (3.6), improved in the membrane and thickness sense and is written as(4.8)fi0=xi=ϕℓXℓi+h02ξ3ϕℓViℓ0,fi1=yi=ϕℓYℓi+h02[ξ3+ϕℓAℓξ32]ϕkVik1,
in which the linear rate of thickness variation (scalar) is parameterized by its nodal values Aℓ as follows:(4.9)a=ϕℓAℓ.
The current position has seven unknown parameters per each node ℓ, that is, three positions Yℓi, three generalized nodal vectors Viℓ, and the nodal rate of thickness variation Aℓ.Functionfi0 is used to find A0 while function fi1 is used to find A1 (trial). The composition of these two values for each integration point gives the numerical value of the gradient of the change of configuration for any initial geometry (curved), that is,(4.10)A=A1(A0)-1.
It is worth to show the derivatives of fi1 regarding the nondimensional variables ξj, constituting the gradient Aij1 as follows:(4.11)Ai11=ϕℓ,1Yℓi+h02{[ξ3+ϕℓAℓξ32]ϕk,1Vik+[ϕℓ,1Aℓξ32]ϕkVik},Ai21=ϕℓ,2Yℓi+h02{[ξ3+ϕℓAℓξ32]ϕk,2Vik+[ϕℓ,2Aℓξ32]ϕkVik}Ai31=h02[1+2ϕℓAℓξ3]ϕkVik.,


## 5. Dynamic Nonlinear Equation

From this section on, a unified notation will be adopted for nodal parameters, they will be called simplyYℓi for each element ℓ and  for  i  varying from 1 to 7. The correspondence is as follows, translations for i varying from 1 to 3, rate of thickness variation for i=4 and generalized vectors for i varying from 5 to 7.The conservation of energy in a mechanical system is guaranteed if the input and output of energy are at balance. If there is some kind of dissipation, the total energy  of the system changes along time. It can be understood  by writing the total potential energy of a system as follows:(5.1)Π=Π0-Q(t,x),
where, following Lánczos [20], Q(t,x) can be stated as the quantity of energy withdrawn from the simple conservative idealized energy Π0. As a consequence, Π is the remaining (current) mechanical energy of the system, and (5.1) can be rewritten as(5.2)Π0=Π+Q(t,x).This equation can be understood as recovering the possibility of using stationary properties for the mechanical system analysis, that is, one can use the minimum potential energy theorem on the ideal energy functionΠ0 for equilibrium analysis.For a structural problem associated with a fixed reference system, Figure5, the ideal potential energy function can be written as the composition of the strain energy (Ue), the potential energy of applied conservative forces P, the kinetic energy (K) and dissipation (Q), as follows:Figure 5
Deformation of a beam under the action of a force.(5.3)Π0=Ue-P+K+Q.
In this work the nonconservative forces will be considered as part of the dissipative potential. The strain energy function of the body, shell for instance, is considered to be stored in the initial volume of the body V0 and is written as an integral of the specific strain energy value ue, (2.3), as(5.4)Ue=∫V0uedV0.The strain energy is assumed to be zero at the initial position, called nondeformed. The potential energy of the applied conservative forces is written as(5.5)P=FiYi+∫s0tiyids0,
where Fi represents forces applied in  direction "i" and Yi is the ith current position of the point where the load is applied, ti is the distributed force applied in direction "i" and yi is the current position of mid-surface points (i=1,2,3 only). Gravitational force has not been mentioned; however, as it is a conservative force, t can be introduced directly into the integral of (5.5). The letter ds0 represents the initial differential area of elements. The kinetic energy is written as(5.6)K=12∫V0ρ0ẏiẏidv0,
where ẏi are velocities and ρ0 is the mass density, relative to the initial volume V0. The dissipative term, including normal distributed forces, is written in its differential form as(5.7)∂∂yiQ(t,y)=∫v0∂∂yiq̅(y,t)dv0=∫v0λmρ0ẏidv0-∫s0qids0,
where q̅ is the specific dissipative functional, λm is a proportional damping constant, ẏi are velocities at any point, and qi are components of the normal distributed forces given by(5.8)qi=qϕℓVℓi1,
where q is the normal distributed force over the element and Vℓi1the generalized vector corresponding to the values Yℓi for i varying from 5 to 7. The integral of these forces respects the direction of the current position but its integral is performed over the initial surface as the load magnitude is written regarding initial position. The current load is easily known by multiplying these magnitudes by ds0/ds, for usual applications of thin structures the value of ds0/ds≅1.Substituting  (5.4), (5.5) and (5.6) in (5.3) results(5.9)Π0=∫V0uedV0-FiYi-∫s0tiyids0+12∫V0ρ0ẏiẏidV0+Q.
This energy function can be written substituting the exact position field by its approximation described in Section 4, that is,(5.10)Π0=∫V0ue(ξ,Yi)dV0-FiYi-∫s0tiyi(ξ,Yi)ds0+12∫V0ρ0ẏi(ξ,Yi)ẏi(ξ,Yi)dV0+Q(ξ,Yi).
The minimum potential energy theorem is used on Π0 by differentiating (5.10) regarding a generic nodal position Yℓj, resulting(5.11)∂Π0∂Yℓj=∫V0∂ue(ξ,Yi)∂YℓjdV0+∫V0ρ0ϕℓj(ξ)ϕki(ξ)dV0Ÿki+∫V0λmρ0ϕℓj(ξ)ϕki(ξ)dV0Ẏki-Fℓj-∫s0ϕℓj(ξ)ϕki(ξ)ds0tki-∫s0qjϕℓj(ξ)ds0=0.
It is worth noting that in this equation the dissipative potential is differentiated regarding nodal positions, differently from (5.7), so (5.8) should be introduced in the last integral of (5.11) to perform the numerical integration. Moreover, as the vibration frequency of the thickness variation is too high, when compared to the other movements, the mass matrix is generated neglecting this term.One can rewrite (5.11) in a simple vector form as(5.12)gℓj=∂Ue∂Yℓj+Fℓjinert.+Fℓjdamp.-Fℓjc-Fℓjnc=0,
or(5.13)gℓj=∂Ue∂Yℓj+MŸℓj+CẎℓj-Fℓjc-Fℓjnc=0.
The involved forces are, inertial force Fℓjinert.or MŸℓj, damping Fℓjdamp.or CẎℓj and the external force, divided into conservative Fℓjc and following forces Fℓjnc. It is important to note that the second representation of the inertial and damping forces is possible because the simple vector mapping described in Sections 3 and 4 generates constant mass matrix. Splitting the derivative of the specific strain energy, one writes(5.14)12∂∂Yℓj(EklCklimEim)=12∂∂Eαβ(EklCklimEim)∂Eαβ∂Yℓj=CαβimEim∂Eαβ∂Yℓj=Sαβ∂Eαβ∂Yℓj.
Consequently(5.15)Fℓjint.=∫V0CαβimEim∂Eαβ∂YℓjdV0,
where Fℓjint is the first gradient vector of the strain energy potential regarding positions, understood as internal force. Equation (5.12) represents the dynamic equilibrium of the shell in the D’Alambert sense. If not, vector gℓj can be understood as the unbalanced force of the mechanical system.The current position is the unknown of the problem, so it is necessary to solve the nonlinear (5.12) regarding Yℓj and time. The first solution step is to integrate (5.12) regarding time. For an implicit approach this step is of great importance regarding the momentum conserving properties of the adopted time integrator. In this work a proof (alternative to the one given by Kane et al. [21]) that the Newmark β method conserves linear momentum and angular momentum for any adopted time step is given. This proof is restrict to total Lagrangian formulation (not corotational) and is trivially extended to energy conserving property for rigid bodies.


## 6. Proof of the Linear and Angular Momentum Conserving Properties of the Newmarkβ Method

In this section no indexes are used, so the variables, as they appear, are vectors. It is important to mention that the proofs given here are restricted to total Lagrangian formulations. It is not extended to nonlagrangian or corotational formulations.The Newmarkβ approximations, following the notation given by Argyris and Mlejnek [22], for position description are(6.1)YS+1=YS+ΔtẎS+Δt2[(12-β)ŸS+βŸS+1],(6.2)ẎS+1=ẎS+Δt(1-γ)ŸS+γΔtŸS+1.
The linear momentum expression for a total Lagrangian description is given as(6.3)Q=∫V0ρ0ẎdV0.
If the body does not develop any deformation and the external forces are zero, then the linear momentum does not vary along time, that is,(6.4)∂Q∂t=∫V0ρ0ŸdV0.
From the continuity  of the body, Ogden [17], one concludes that, the following equation is valid(6.5)Ÿ=0.
Using (6.5), written for two times ts and ts+1, into (6.2), it results(6.6)ẎS+1=ẎS,
or by continuity(6.7)∫V0ρ0Ẏs+1dV0=∫V0ρ0ẎsdV0,
that is, the linear momentum is conserved for any adopted time step and Newmark constants.To prove the conservation of angular momentum more steps are required. The Lagrangian expression of angular momentum is(6.8)J=∫V0ρ0ẎxYdV0,
where x is the vector product.The angular momentum is constant when there is a fixed axis around which a rigid body turns with constant angular velocity. As the body is considered rigid, no transfer of energy from kinetics to strain energy occurs; as a consequence the conservation of momentum means the conservation of energy for an isothermal situation. Assuming this hypothesis one writes(6.9)∂J∂t=∫V0ρ0(ŸxY+ẎxẎ)dV0=∫V0ρ0(ŸxY)dV0=0.
Form the continuity assumption the equality(6.10)ŸxY=0
must hold for any point of the continuum. It occurs in two situations, the trivial undesired one, that is, Ÿ=0 and the desired one,(6.11)Ÿ=-ω2Y,
where ω is the angular velocity of the body around the rotation axis and Y is, without loss of generality, the position vector of the point related to its projection over the rotation axis. Using (6.11) for time ts+1 into the Newmark β (6.1) and (6.2), one writes(6.12)YS+1=YS+ΔtẎS-Δt2(12-β)ω2YS-Δt2βω2YS+1,(6.13)ΔtẎS+1=ΔtẎS-Δt2(1-γ)ω2YS-Δt2γω2YS+1.
Rearranging terms of (6.12), one has(6.14)ΔtẎS=(Δt2βω2+1)YS+1-(1-Δt2(12-β)ω2)YS.
Substituting (6.14) into  (6.13) results(6.15)ΔtẎS+1=(Δt2βω2+1-Δt2γω2)YS+1-(1-Δt2(12-β)ω2+Δt2(1-γ)ω2)YS.
Post-vector-multiplying (6.14) and (6.15) by Ys and Ys+1, respectively, and subtracting results one achieves(6.16)Δt(Ẏs+1xYs+1-ẎsxYs)=(1-Δt2(12-β)ω2+Δt2(1-γ)ω2-Δt2βω2-1)(Ys+1xYs).
Finally (6.16) simplifies to(6.17)(Ẏs+1xYs+1-ẎsxYs)=Δtω2(12-γ)(Ys+1xYs).
By continuity one writes(6.18)∫V0ρ0Ẏs+1xYs+1dV0=∫V0ρ0ẎsxYsdV0+∫V0ρ0Δtω2(12-γ)(Ys+1xYs)dV0,
or in a shorter notation(6.19)Js+1=Js+∫V0ρ0Δtω2(12-γ)(Ys+1xYs)dV0.
Therefore, the Newmark β method is angular momentum conserving for γ=1/2, despite the adopted time step, angular velocity or β parameter.


## 7. Time Marching Process and Newton Rapson Procedure

From the previous developments (5.12) can be written in a simpler form as(7.1)g=∂Ue∂Y-Fc-Fnc+MŸ+CẎ=0.
Expression (7.1) represents the dynamic equilibrium equation at any time and has to be solved. In order to do so the first step is to write this equilibrium for a specific instant  ts+1 as follows:(7.2)∂Π∂Y|S+1=∂Ue∂Y|S+1-FS+1+MŸS+1+CẎS+1=0.
Substituting the Newmark approximations (6.1) and (6.2) into (7.2) results(7.3)g(YS+1)=∂Π∂Y|S+1=∂Ue∂Y|S+1-FS+1+MβΔt2YS+1-MQS+CRS+γCβΔtYS+1-γΔtCQS=0,
where vectors QS and RS represent the dynamic contribution of the past and are given by(7.4)QS=YSβΔt2+ẎSβΔt+(12β-1)ŸS,RS=ẎS+Δt(1-γ)ŸS.
Equation (7.3) can be understood simply by g(YS+1)=0 and is clearly nonlinear regarding (YS+1). A Taylor expansion to solve this nonlinear equation is necessary. The second derivative of the total energy potential is then given by(7.5)∂2Π∂Y2|S+1=∇g(Ys+1)=∂2Ue∂Y2|S+1+MβΔt2+γCβΔt.
One builds the Taylor series of first order as:(7.6)0=g(Y)≅g(Y0)+∇g(Y0)ΔY
and derives the Newton-Raphson procedure to solve the nonlinear  (7.3), that is,(7.7)∇g(Y0)ΔY=-g(Y0),
where Y0 is a trial position (usually Ys) for Ys+1 used in (7.3) to calculate g(Y0). Solving ΔY one calculates a new trial for Ys+1 as(7.8)YS+1=Y0+ΔY.
The acceleration must be corrected for each iteration by an expression obtained from (6.1), that is,(7.9)ŸS+1=YS+1βΔt2-QS.
This equation is used in (6.2) to correct velocity. The stop criterion is given by (7.10), when a chosen tolerance (TOL) is satisfied, that is,(7.10)∥g(Y0)∥≤TOLor∥ΔY∥≤TOL.
It must be noted that, before the first time step, the initial acceleration must be calculated as follows:(7.11)Ÿ0=M-1[F0-∂Ue∂Y|0-CẎ0].


## 8. The Derivatives of the Specific Strain Energy

In order to conclude the description of the formulation the second derivatives of the strain energy regarding nodal positions should be given as it has been done for the first derivative in (5.15). From (5.14) and (5.15) one writes(8.1)∂2Ue∂Yk∂Yj=∫V0∂∂Yk(CαβimEim∂Eαβ∂Yj)dV0=∫V0(∂Eim∂YkCαβim∂Eαβ∂Yj+EimCαβim∂2Eαβ∂Yj∂Yk)dV0.
Finally, the first and second derivatives of the Green strain regarding current nodal positions should be done. Firstly the necessary derivatives of the Cauchy-Green stretch tensor are presented. Next the derivatives of strains are straightforwardly achieved. Recalling that the Cauchy-Green stretch tensor is given by(8.2)C=AtA
and omitting, for simplicity, extra indices, one applies the proposed mapping, that is,  A=A1(A0)-1, and writes(8.3)C=[(A0)t]-1(A1)t(Yi)A1(Yi)(A0)-1.Remembering thatA0 is constant regarding the current position, the first derivative is performed as(8.4)∂C∂Yj=[(A0)t]-1∂(A1)t(Yi)∂YjA1(Yi)(A0)-1+[(A0)t]-1(A1)t(Yi)∂A1(Yi)∂Yj(A0)-1
and, from (4.11), one has the following (not null) values of the current positional mapping gradient:(8.5)∂A111∂Yℓ1=ϕℓ,1(ξ1,ξ2),∂A121∂Yℓ1=ϕℓ,2(ξ1,ξ2),∂A211∂Yℓ2=ϕℓ,1(ξ1,ξ2),∂A221∂Yl2=ϕl,2(ξ1,ξ2),∂A311∂Yℓ3=ϕℓ,1(ξ1,ξ2),∂A321∂Yℓ3=ϕℓ,2(ξ1,ξ2),∂Ai11∂Yℓ4=h02ξ32{ϕℓ,1(ξ1,ξ2)Vkϕk(ξ1,ξ2)+ϕℓ(ξ1,ξ2)Vkϕk,1(ξ1,ξ2)},∂Ai21∂Yℓ4=h02ξ32{ϕℓ,2(ξ1,ξ2)Vkϕk(ξ1,ξ2)+ϕℓ(ξ1,ξ2)GVkϕk,2(ξ1,ξ2)},∂Ai31∂Yℓ7=h0ξ3ϕℓ(ξ1,ξ2)Vkϕk(ξ1,ξ2),∂A111∂Yℓ5=h02{ϕℓ,1(ξ1,ξ2)ξ3+[ϕk,1(ξ1,ξ2)Akϕℓ(ξ1,ξ2)+ϕk(ξ1,ξ2)Akϕℓ,1(ξ1,ξ2)]ξ32},∂A121∂Yℓ6=h02{ϕℓ,2(ξ1,ξ2)ξ3+[ϕk,2(ξ1,ξ2)Akϕℓ(ξ1,ξ2)+ϕk(ξ1,ξ2)Akϕℓ,2(ξ1,ξ2)]ξ32},∂A131∂Yℓ7=h02{ϕℓ(ξ1,ξ2)+2ξ3ϕk(ξ1,ξ2)Akϕℓ(ξ1,ξ2)},∂A211∂Yℓ5=h02{ϕℓ,1(ξ1,ξ2)ξ3+[ϕk,1(ξ1,ξ2)Akϕℓ(ξ1,ξ2)+ϕk(ξ1,ξ2)Akϕℓ,1(ξ1,ξ2)]ξ32},∂A221∂Yℓ6=h02{ϕℓ,2(ξ1,ξ2)ξ3+[ϕk,2(ξ1,ξ2)Akϕℓ(ξ1,ξ2)+ϕk(ξ1,ξ2)Akϕℓ,2(ξ1,ξ2)]ξ32},∂A231∂Yℓ7=h02{ϕℓ(ξ1,ξ2)+2ξ3ϕk(ξ1,ξ2)Akϕℓ(ξ1,ξ2)},∂A311∂Yℓ5=h02{ϕℓ,1(ξ1,ξ2)ξ3+[ϕk,1(ξ1,ξ2)Akϕℓ(ξ1,ξ2)+ϕk(ξ1,ξ2)Akϕℓ,1(ξ1,ξ2)]ξ32},∂A321∂Yℓ6=h02{ϕℓ,2(ξ1,ξ2)ξ3+[ϕk,2(ξ1,ξ2)Akϕℓ(ξ1,ξ2)+ϕk(ξ1,ξ2)Akϕℓ,2(ξ1,ξ2)]ξ32},∂A331∂Yℓ7=h02{ϕℓ(ξ1,ξ2)+2ξ3ϕk(ξ1,ξ2)Akϕℓ(ξ1,ξ2)}.The first index ofYℓm is the element node, and the second is the degree of freedom. The second derivative of the Cauchy-Green stretch is given by(8.6)∂2C∂Yj∂Yk=[(A0)t]-1∂(A1)t(Yi)∂Yj∂A1(Yi)∂Yk(A0)-1+[(A0)t]-1∂(A1)t(Yi)∂Yk∂A1(Yi)∂Yj(A0)-1+[(A0)t]-1∂2(A1)t(Yi)∂Yj∂YkA1(Yi)(A0)-1+[(A0)t]-1(A1)t(Yi)∂A1(Yi)∂Yj∂Yk(A0)-1,
where ∂A1/∂Yk are given by  (8.5) and the second derivatives ∂2A1/(∂Yk∂Yj) are straightforward. Recalling  (2.2), one achieves directly(8.7)∂E∂Yj=12∂C∂Yj,∂2E∂Yj∂Yk=12∂2C∂Yj∂Yk.It is important to mention that the present technique can be applied to any strain measure based on the Cauchy-Green stretch. Equation (7.5) indicates that the proposed procedure can be operated by creating the Hessian matrix and internal forces for finite elements and building the global matrix and internal force vector by summation of coincident degrees of freedom, as it is done for usual FEM procedures. One should remember that all nodal parameters follow the global system of reference, avoiding the use of rotation schemes.


## 9. Numerical Examples

This section provides eight examples covering selected tests to confirm the generality and accuracy of the proposed formulation for static and dynamic situations. No mention to units is made to keep a coherence with references. The important features are objectivity, locking free behavior, Linear momentum conserving, angular momentum conserving, total energy conserving, and generality in applications. More examples regarding static analysis can be seen in Coda and Paccola [23, 24].


### 9.1. Objectivity of the Formulation Regarding Rotations

As mentioned in the introduction, this formulation is tested regarding mapping objectivity. The employed way to test this property follows well-known methodologies; see for instance Criesfield and Jelenić [19] and Ibrahimbergovic and Taylor [25]. A clamped vertical plate (shell in deformed configuration) is subject to a transverse load at its free end as depicted in Figure 6. Static conditions are considered, so inertial forces are neglected, gravitational forces are also not considered. The physical properties of the structure are E=100000 and ν=0. The thickness of the shell is h=0.1. The adopted discretization can also be seen in Figure 6.Geometrical characteristics of the problem and discretization.
(a)(b)Two situations are created. The first consists into applying a rotation over the shell regarding the clamping axis without applying any load. The objective is to show that no stress will be generated at any stage of rotation. One hundred turns are applied and no stress appears; moreover, the positions are exactly the same after each turn. In Figure7 one can see an illustration of this situation for the first turn. The adopted rotation step is 0.1π.Figure 7
Stress values for the first turn—unloaded.In the second situation the process is divided into two phases. First, the load is increased to its final value in ten equal steps. The resulting stress, following the longitudinal direction of the shell, at the superior face of the shell, is depicted in Figure8.Figure 8
Stress values for the first deformed configuration—no rotation.Then the load is kept constant and acts, in the same sense and direction, and the rotation, similar to the one used in the first situation, is applied. The adopted rotation step is0.01π. At the beginning of the rotation process the action of rotation is against the action of the loading. At a quarter of the first turn the loading is compressing the shell and the stress values, following the longitudinal direction of the shell, are depicted in Figure 9.Figure 9
Stress values for a quarter of the first turn (π/2).At the half of the first turn the shell is in opposite position to the beginning of the rotation process, and the initially superior face of the shell is now the inferior one. The stress values at this face are negative and their values are depicted in Figure10. The difference in stress magnitude from the first deformed configuration and this one is due to the normal traction force that increases the values in Figure 8 and decreases the absolute values in Figure 10.Figure 10
Stress values for the half turn (π).At three quarters of the first turn the stress values are the ones depicted in Figure11, exactly as expected.Figure 11
Stress values for three quarters of a turn (3π/2).Finally, after a complete revolution the stress values depicted in Figure12 are exactly the same ones present in Figure 8. Ninety nine more turns were performed, and the results are repeated for each turn, revealing the total objectiveness of the generalized vector mapping.Figure 12
Stress values for one complete turn (2π).


### 9.2. Shear and Volumetric Locking Analysis

This example is extracted from Bucalem et al. [26] and is a benchmark to check FEM formulations regarding shear and volumetric locking. It is an important example as the solution for plates is very sensitive to the Poisson Ratio, and some shell theories fail to reproduce the analytical solution by the presence of shear locking. It is the analysis of a simple supported square plate subjected to a static transverse concentrated load at its center. The numerical results are compared with the analytical solution obtained using Navier’s series for Kirchhoff plate theory. The thin plate geometry is depicted in Figure 13.Figure 13
Analyzed plate and adopted discretization.The adopted physical data areL=2, E=2.1×109, h=0.002, and υ varying from 0 to 0.5. The applied load is P=0.4×10-2. In Figure 14 the results obtained using the presented improved formulation are compared to the analytical solution and the solution for constrained linear rate of thickness variation. This last situation is called simply six-parameter shell formulation.Figure 14
Displacement at the center of the plate versus Poisson Ratio.As one can observe the presented formulation is free from shear and volumetric locking and reproduces perfectly the analytical solution. For more examples regarding the shear locking free behavior of the proposed kinematics one is referred to Coda and Paccola [23, 24], where an extensive static analysis for a six parameter shell element is presented.


### 9.3. Pinched Cylinder with Rigid Diaphragms (Static)

This benchmark consists in a cylinder with rigid diaphragms pinched by concentrated loads at two opposite points at its top and bottom, see Figure15. The adopted discretization (2×18×6 mesh) is also depicted in Figure 15 comprising 1045 nodes. This example is also used to test the formulation regarding shear and Poisson locking. It should be noted that if a formulation suffers from any locking the correct results cannot be achieved. The use of symmetry can hide some buckling modes of the problem; however, to compare with the reference papers the symmetry is assumed. The number of modes one can get with accuracy (excluding nonsymmetric ones) is about 30 for a total of 220 node for the total circumference and cubic approximations. The results are compared to the ones achieved by Sansour and Kollmann [27] with a mesh of 1681 nodes. Taking advantage of symmetry only one octant of the cylinder is discretized. The adopted physical parameters are: R=100, h=1, E=3×104, L=200. Two values are adopted for the Poisson ratio, υ=0.3 and υ=0.49, respectively, in order to check the locking free behavior of the presented formulation.Figure 15
Pinched Cylinder geometry, loading and the adopted discretization.The problem is solved with the proposed formulation following two strategies. The first, called 7 parameters, does not constrain the linear rate of thickness. The second, called 6 parameters, totally constrains the linear rate of thickness variation. In Figure16 the results for υ=0.3 are compared  with Sansour and Kollmann [27] that employed an enhanced strain quadrilateral element. As one can see the formulation proposed here can capture the flexibility of the pinched shell for ν=0.3 even if the linear rate of thickness variation is totally constrained. A formulation suffering of shear locking is not able to run this example.Figure 16
Displacements for pointsA and B, υ=0.3.In Figure17, the behavior of the proposed formulation (with the seventh parameter constrained or not) is depicted for υ=0.49. The reference value for this figure is the seven parameter result with υ=0.3. As expected, the proposed formulation does not lock for large Poisson ratio. However, the results for the totally constrained rate of thickness variation lock completely. Shell formulations based on six parameters and full constitutive relations can be free from shear locking, but not from Poisson locking.Figure 17
Displacements for pointsA and B, υ=0.49.From this result it  is obvious that the positional improved formulation based on generalized vectors, together with high-order curved elements, is able to solve geometrically nonlinear shell problems with precision and reduced mesh. Additional and practical information is that the pinched cylinder benchmark problem is not very sensitive to the Poisson ratio intensity.


### 9.4. Linear Momentum Conservation

This example is used to confirm the proof given for the linear momentum conserving property of the Newmarkβ method when used with the generalized vector mapping to develop the positional formulation. A plate with the same dimensions of Section 9.1 with no displacement restrictions and mass density ρ0=1 is subjected to an initial velocity in the vertical direction of value V2=1. As the gravitational force is neglected and no displacement restrictions are imposed the plate does not deform and moves with constant velocity. The total applied kinetic energy is, of course, k=0.2. Only two finite elements were used to run this problem. In Figure 18 the numerical kinetic energy calculated for 1000 time steps using different time steps is depicted. Remembering that the Newmark parameter γ=0.5 is mandatory the other one is adopted as β=1/4.Figure 18
Energy conserving for linear momentum.As one can (see in Figure18)  the Newmark β method preserves the total energy of the system when a linear momentum is applied. No strains were developed in this example.


### 9.5. Angular Momentum Conservation

In this example a vertical circular thick cylinder, (see Figure19), is subjected to a field of initial velocity generated by an angular velocity of w=1. The dimensions of the cylinder are: h=0.4, l=1.0 and r=1.0. The physical parameters are: E=100, ν=0 and ρ0=0.1. The total energy introduced in the problem by this velocity field is (considering the thickness influence)  k=1.3069. Eighteen finite elements with forty eight nodes and β=0.25 are used to run this problem, as shown in Figure 19. One should note that the elements are curved despite the straight lines that appear in Figure 19. The achieved total energy is depicted along 1000 time steps in Figure 20 for two different time steps. It is interesting noting that a small strain energy was computed for this analysis.Figure 19
Geometry and discretization.Figure 20
Energy conserving for angular momentum.The exact solution is coincident with the numerical one forΔt=0.01. The numerical result for the large time step is also energy conserving; the difference in results is floating and less than 1%. This difference is due to the better accuracy achieved when using a small time step. In Figure 20 the number of time steps is the same, however for the larger one 15.9 turns are depicted and only 1.59 turns are depicted for the smaller one. The same total energy is found after 1000000 time steps for the smaller time step, that is, performing 1590 turns. Even better results are achieved for thinner shells. As a consequence, the proof given for the momentum conserving property of the Newmark β method is confirmed also for angular momentum.


### 9.6. Transverse Dynamic Load over a Clamped Beam—Energy Conservation Check

This example adopts the proposed nonlinear shell element to run the nonlinear transversal dynamic vibration of a clamped beam subjected to an initial field of velocity, see Figure21. It intends to show the transfer behavior of kinetic and strain energies for a deformable body. The beam has young Modulus E=0.2×109, mass density ρ0=500, length L=1, thickness h=0.01, width b=0.20 and Poisson’s ratio ν=0. The applied velocity is proportional to the distance from the clamped end and has a maximum value of v0=1, see Figure 20. The adopted time step is Δt=0.001. The exact total energy of the system is k0=0.16667.Figure 21
Geometry, initial velocity and discretization.Forty finite elements and 217 nodes are employed in the discretization. In Figure22, the kinetic, strain and total energies are depicted. The transversal displacement of the beam at the free end is also depicted in Figure 22. As one can see the energy is completely conserved for deformable situations. It is important to mention that in nonlinear applications the coupling of vertical and horizontal movements is present, and this is the reason why the peaks of the analysis  do not repeat with the same shape, also there is not necessarily an instant for which the kinetic energy becomes zero, as it is expected in some simple linear analysis. Although it is not mentioned, the results for bar examples were extensively tested against the ones presented by  Greco and Coda [28].Figure 22
Energy and free end displacement along time.


### 9.7. Simple Airbag Simulation

This simple example is inspired in the airbag simulation presented by Cirak and Radovitzky [29] that includes fluid structure interaction. In the present work only the structure (airbag) is modeled subjected to a deterministic applied load in order to demonstrate the possibilities of the proposed formulation regarding the analysis of very thin membranes and general problems. A formulation suffering of shear locking is not able to run this example. The load is orthogonal to the airbag surface and, as no comparisons can be made, the initial structure discretization is simplified. The simulation corresponds to an initially-flat airbag made of an elastic fabric with  Young’s modulus of E=6.0×109, Poisson’s ratio of ν=0.3, and mass density of ρ0=1000. The initial position of the real airbag can be seen in Figure 23, extracted from the work of Cirak and Radovitzky [29]. The thickness of the airbag is h=7.3×10-4 and the diameter in its flat initial configuration is D=0.74. Our initial discretization (1/8th of the airbag due to the assumed double symmetry) consists of 800 cubic elements of ten nodes resulting in 3721 nodes, Figure 24. The adopted boundary condition at the curved board is simply supported. The pressure of the fluid is simulated by the following function:Figure 23
Airbag geometry, Cirak and Radovitzky [29].Figure 24
Adopted discretization, initially flat.(9.1)100te(0.999t-0.01)r25066.25for0<t<0.01,5066.25fort>0.01,
where t is time and r is the distance from the center of the airbag. The adopted time step is Δt=0.0002 and the maximum tolerance for convergence is tol=0.0005 in positions. An initial random vertical defect with maximum amplitude of δ=0.001 is assumed for the analysis.Some selected deformed positions are depicted in Figures25, 26, and 27. In Figure 28 the time story of the displacement of the top of the airbag is depicted. As one can see in Figure 28 the top of the airbag stabilizes, to the final value, by geometrical accommodation.Figure 25
Airbag at 0.0042s.Figure 26
Airbag at 0.0084s.Figure 27
Airbag at  0.0314s.Figure 28
Top displacement along time.


### 9.8. Cylindrical Shell with Dynamic Snap Through

Following Argyris et al. [14], the second example is a cylindrical shell exhibiting dynamic snap through, a severe nonlinearity. Important theoretical studies, related to severe nonlinearities are presented by Breslavsky et al. [30] and others cited therein. This is a typical benchmark example that has been used extensively as a test for all nonlinear shell dynamics formulations presented so far. Snap through problems in shells produces higher dynamical modes and this is the reason why it is believed that the standard integration schemes such as the Newmark method are not adequate to produce a stable and accurate solution and that only algorithms with numerical dissipation and energy decaying schemes can be applied with an acceptable time step. Despite this widespread belief, it will be shown that using the proposed shell element it is possible to obtain stable and accurate solutions with the Newmark integration for reasonable time step. Results will be compared to the TRIC element and an ABACUS solution, both presented by Argyris et al. [14].The geometry of the cylindrical shell is shown in Figure29. The two straight edges of the shell are simply supported, while the two curved edges are free. A concentrated load is applied at the central node of the shell. The value of this load increases linearly from 0 to 50×106 in a time of 0.2, after that it is held constant. To avoid any distortion in results the structure is totally discretized. The mesh used in the analysis is shown in Figure 30 and consists of 32 curved finite elements with cubic approximation resulting in 169 nodes and 1183 degrees of freedom. Argyris et al. [14] used (for a symmetric quarter of the structure) 128 TRIC elements with 407 degrees of freedom to run this example, see Figure 30. The adopted Newmark Parameters are β=0.25 and γ=0.5. The TRIC element is designed to run this example using large time step within the context of corotational formulation. by the coupling of rigid body movements, for the overall element, and strain modes. Argyris et al. [14] found a maximum stable time step for their co rotational formulation of Δt=1×10-3. Using the proposed formulation we vary the time step from Δt=0.0625×10-3 to Δt=4×10-3 and found not instability for the total Lagrangian formulation. The adopted physical properties are E=200×109, v=0.25, ρ=10.000 and thickness of h=0.1.Figure 29
Geometry.Figure 30
Discretizations.Figure 31
Present solution for the apex displacement.Figure 32
Results presented by Argyris et al. [14].As one can see, the results converge to the smaller time step, while for large time steps the accuracy is lost, but the stability is kept.


## 9.1. Objectivity of the Formulation Regarding Rotations

As mentioned in the introduction, this formulation is tested regarding mapping objectivity. The employed way to test this property follows well-known methodologies; see for instance Criesfield and Jelenić [19] and Ibrahimbergovic and Taylor [25]. A clamped vertical plate (shell in deformed configuration) is subject to a transverse load at its free end as depicted in Figure 6. Static conditions are considered, so inertial forces are neglected, gravitational forces are also not considered. The physical properties of the structure are E=100000 and ν=0. The thickness of the shell is h=0.1. The adopted discretization can also be seen in Figure 6.Geometrical characteristics of the problem and discretization.
(a)(b)Two situations are created. The first consists into applying a rotation over the shell regarding the clamping axis without applying any load. The objective is to show that no stress will be generated at any stage of rotation. One hundred turns are applied and no stress appears; moreover, the positions are exactly the same after each turn. In Figure7 one can see an illustration of this situation for the first turn. The adopted rotation step is 0.1π.Figure 7
Stress values for the first turn—unloaded.In the second situation the process is divided into two phases. First, the load is increased to its final value in ten equal steps. The resulting stress, following the longitudinal direction of the shell, at the superior face of the shell, is depicted in Figure8.Figure 8
Stress values for the first deformed configuration—no rotation.Then the load is kept constant and acts, in the same sense and direction, and the rotation, similar to the one used in the first situation, is applied. The adopted rotation step is0.01π. At the beginning of the rotation process the action of rotation is against the action of the loading. At a quarter of the first turn the loading is compressing the shell and the stress values, following the longitudinal direction of the shell, are depicted in Figure 9.Figure 9
Stress values for a quarter of the first turn (π/2).At the half of the first turn the shell is in opposite position to the beginning of the rotation process, and the initially superior face of the shell is now the inferior one. The stress values at this face are negative and their values are depicted in Figure10. The difference in stress magnitude from the first deformed configuration and this one is due to the normal traction force that increases the values in Figure 8 and decreases the absolute values in Figure 10.Figure 10
Stress values for the half turn (π).At three quarters of the first turn the stress values are the ones depicted in Figure11, exactly as expected.Figure 11
Stress values for three quarters of a turn (3π/2).Finally, after a complete revolution the stress values depicted in Figure12 are exactly the same ones present in Figure 8. Ninety nine more turns were performed, and the results are repeated for each turn, revealing the total objectiveness of the generalized vector mapping.Figure 12
Stress values for one complete turn (2π).


## 9.2. Shear and Volumetric Locking Analysis

This example is extracted from Bucalem et al. [26] and is a benchmark to check FEM formulations regarding shear and volumetric locking. It is an important example as the solution for plates is very sensitive to the Poisson Ratio, and some shell theories fail to reproduce the analytical solution by the presence of shear locking. It is the analysis of a simple supported square plate subjected to a static transverse concentrated load at its center. The numerical results are compared with the analytical solution obtained using Navier’s series for Kirchhoff plate theory. The thin plate geometry is depicted in Figure 13.Figure 13
Analyzed plate and adopted discretization.The adopted physical data areL=2, E=2.1×109, h=0.002, and υ varying from 0 to 0.5. The applied load is P=0.4×10-2. In Figure 14 the results obtained using the presented improved formulation are compared to the analytical solution and the solution for constrained linear rate of thickness variation. This last situation is called simply six-parameter shell formulation.Figure 14
Displacement at the center of the plate versus Poisson Ratio.As one can observe the presented formulation is free from shear and volumetric locking and reproduces perfectly the analytical solution. For more examples regarding the shear locking free behavior of the proposed kinematics one is referred to Coda and Paccola [23, 24], where an extensive static analysis for a six parameter shell element is presented.


## 9.3. Pinched Cylinder with Rigid Diaphragms (Static)

This benchmark consists in a cylinder with rigid diaphragms pinched by concentrated loads at two opposite points at its top and bottom, see Figure15. The adopted discretization (2×18×6 mesh) is also depicted in Figure 15 comprising 1045 nodes. This example is also used to test the formulation regarding shear and Poisson locking. It should be noted that if a formulation suffers from any locking the correct results cannot be achieved. The use of symmetry can hide some buckling modes of the problem; however, to compare with the reference papers the symmetry is assumed. The number of modes one can get with accuracy (excluding nonsymmetric ones) is about 30 for a total of 220 node for the total circumference and cubic approximations. The results are compared to the ones achieved by Sansour and Kollmann [27] with a mesh of 1681 nodes. Taking advantage of symmetry only one octant of the cylinder is discretized. The adopted physical parameters are: R=100, h=1, E=3×104, L=200. Two values are adopted for the Poisson ratio, υ=0.3 and υ=0.49, respectively, in order to check the locking free behavior of the presented formulation.Figure 15
Pinched Cylinder geometry, loading and the adopted discretization.The problem is solved with the proposed formulation following two strategies. The first, called 7 parameters, does not constrain the linear rate of thickness. The second, called 6 parameters, totally constrains the linear rate of thickness variation. In Figure16 the results for υ=0.3 are compared  with Sansour and Kollmann [27] that employed an enhanced strain quadrilateral element. As one can see the formulation proposed here can capture the flexibility of the pinched shell for ν=0.3 even if the linear rate of thickness variation is totally constrained. A formulation suffering of shear locking is not able to run this example.Figure 16
Displacements for pointsA and B, υ=0.3.In Figure17, the behavior of the proposed formulation (with the seventh parameter constrained or not) is depicted for υ=0.49. The reference value for this figure is the seven parameter result with υ=0.3. As expected, the proposed formulation does not lock for large Poisson ratio. However, the results for the totally constrained rate of thickness variation lock completely. Shell formulations based on six parameters and full constitutive relations can be free from shear locking, but not from Poisson locking.Figure 17
Displacements for pointsA and B, υ=0.49.From this result it  is obvious that the positional improved formulation based on generalized vectors, together with high-order curved elements, is able to solve geometrically nonlinear shell problems with precision and reduced mesh. Additional and practical information is that the pinched cylinder benchmark problem is not very sensitive to the Poisson ratio intensity.


## 9.4. Linear Momentum Conservation

This example is used to confirm the proof given for the linear momentum conserving property of the Newmarkβ method when used with the generalized vector mapping to develop the positional formulation. A plate with the same dimensions of Section 9.1 with no displacement restrictions and mass density ρ0=1 is subjected to an initial velocity in the vertical direction of value V2=1. As the gravitational force is neglected and no displacement restrictions are imposed the plate does not deform and moves with constant velocity. The total applied kinetic energy is, of course, k=0.2. Only two finite elements were used to run this problem. In Figure 18 the numerical kinetic energy calculated for 1000 time steps using different time steps is depicted. Remembering that the Newmark parameter γ=0.5 is mandatory the other one is adopted as β=1/4.Figure 18
Energy conserving for linear momentum.As one can (see in Figure18)  the Newmark β method preserves the total energy of the system when a linear momentum is applied. No strains were developed in this example.


## 9.5. Angular Momentum Conservation

In this example a vertical circular thick cylinder, (see Figure19), is subjected to a field of initial velocity generated by an angular velocity of w=1. The dimensions of the cylinder are: h=0.4, l=1.0 and r=1.0. The physical parameters are: E=100, ν=0 and ρ0=0.1. The total energy introduced in the problem by this velocity field is (considering the thickness influence)  k=1.3069. Eighteen finite elements with forty eight nodes and β=0.25 are used to run this problem, as shown in Figure 19. One should note that the elements are curved despite the straight lines that appear in Figure 19. The achieved total energy is depicted along 1000 time steps in Figure 20 for two different time steps. It is interesting noting that a small strain energy was computed for this analysis.Figure 19
Geometry and discretization.Figure 20
Energy conserving for angular momentum.The exact solution is coincident with the numerical one forΔt=0.01. The numerical result for the large time step is also energy conserving; the difference in results is floating and less than 1%. This difference is due to the better accuracy achieved when using a small time step. In Figure 20 the number of time steps is the same, however for the larger one 15.9 turns are depicted and only 1.59 turns are depicted for the smaller one. The same total energy is found after 1000000 time steps for the smaller time step, that is, performing 1590 turns. Even better results are achieved for thinner shells. As a consequence, the proof given for the momentum conserving property of the Newmark β method is confirmed also for angular momentum.


## 9.6. Transverse Dynamic Load over a Clamped Beam—Energy Conservation Check

This example adopts the proposed nonlinear shell element to run the nonlinear transversal dynamic vibration of a clamped beam subjected to an initial field of velocity, see Figure21. It intends to show the transfer behavior of kinetic and strain energies for a deformable body. The beam has young Modulus E=0.2×109, mass density ρ0=500, length L=1, thickness h=0.01, width b=0.20 and Poisson’s ratio ν=0. The applied velocity is proportional to the distance from the clamped end and has a maximum value of v0=1, see Figure 20. The adopted time step is Δt=0.001. The exact total energy of the system is k0=0.16667.Figure 21
Geometry, initial velocity and discretization.Forty finite elements and 217 nodes are employed in the discretization. In Figure22, the kinetic, strain and total energies are depicted. The transversal displacement of the beam at the free end is also depicted in Figure 22. As one can see the energy is completely conserved for deformable situations. It is important to mention that in nonlinear applications the coupling of vertical and horizontal movements is present, and this is the reason why the peaks of the analysis  do not repeat with the same shape, also there is not necessarily an instant for which the kinetic energy becomes zero, as it is expected in some simple linear analysis. Although it is not mentioned, the results for bar examples were extensively tested against the ones presented by  Greco and Coda [28].Figure 22
Energy and free end displacement along time.


## 9.7. Simple Airbag Simulation

This simple example is inspired in the airbag simulation presented by Cirak and Radovitzky [29] that includes fluid structure interaction. In the present work only the structure (airbag) is modeled subjected to a deterministic applied load in order to demonstrate the possibilities of the proposed formulation regarding the analysis of very thin membranes and general problems. A formulation suffering of shear locking is not able to run this example. The load is orthogonal to the airbag surface and, as no comparisons can be made, the initial structure discretization is simplified. The simulation corresponds to an initially-flat airbag made of an elastic fabric with  Young’s modulus of E=6.0×109, Poisson’s ratio of ν=0.3, and mass density of ρ0=1000. The initial position of the real airbag can be seen in Figure 23, extracted from the work of Cirak and Radovitzky [29]. The thickness of the airbag is h=7.3×10-4 and the diameter in its flat initial configuration is D=0.74. Our initial discretization (1/8th of the airbag due to the assumed double symmetry) consists of 800 cubic elements of ten nodes resulting in 3721 nodes, Figure 24. The adopted boundary condition at the curved board is simply supported. The pressure of the fluid is simulated by the following function:Figure 23
Airbag geometry, Cirak and Radovitzky [29].Figure 24
Adopted discretization, initially flat.(9.1)100te(0.999t-0.01)r25066.25for0<t<0.01,5066.25fort>0.01,
where t is time and r is the distance from the center of the airbag. The adopted time step is Δt=0.0002 and the maximum tolerance for convergence is tol=0.0005 in positions. An initial random vertical defect with maximum amplitude of δ=0.001 is assumed for the analysis.Some selected deformed positions are depicted in Figures25, 26, and 27. In Figure 28 the time story of the displacement of the top of the airbag is depicted. As one can see in Figure 28 the top of the airbag stabilizes, to the final value, by geometrical accommodation.Figure 25
Airbag at 0.0042s.Figure 26
Airbag at 0.0084s.Figure 27
Airbag at  0.0314s.Figure 28
Top displacement along time.


## 9.8. Cylindrical Shell with Dynamic Snap Through

Following Argyris et al. [14], the second example is a cylindrical shell exhibiting dynamic snap through, a severe nonlinearity. Important theoretical studies, related to severe nonlinearities are presented by Breslavsky et al. [30] and others cited therein. This is a typical benchmark example that has been used extensively as a test for all nonlinear shell dynamics formulations presented so far. Snap through problems in shells produces higher dynamical modes and this is the reason why it is believed that the standard integration schemes such as the Newmark method are not adequate to produce a stable and accurate solution and that only algorithms with numerical dissipation and energy decaying schemes can be applied with an acceptable time step. Despite this widespread belief, it will be shown that using the proposed shell element it is possible to obtain stable and accurate solutions with the Newmark integration for reasonable time step. Results will be compared to the TRIC element and an ABACUS solution, both presented by Argyris et al. [14].The geometry of the cylindrical shell is shown in Figure29. The two straight edges of the shell are simply supported, while the two curved edges are free. A concentrated load is applied at the central node of the shell. The value of this load increases linearly from 0 to 50×106 in a time of 0.2, after that it is held constant. To avoid any distortion in results the structure is totally discretized. The mesh used in the analysis is shown in Figure 30 and consists of 32 curved finite elements with cubic approximation resulting in 169 nodes and 1183 degrees of freedom. Argyris et al. [14] used (for a symmetric quarter of the structure) 128 TRIC elements with 407 degrees of freedom to run this example, see Figure 30. The adopted Newmark Parameters are β=0.25 and γ=0.5. The TRIC element is designed to run this example using large time step within the context of corotational formulation. by the coupling of rigid body movements, for the overall element, and strain modes. Argyris et al. [14] found a maximum stable time step for their co rotational formulation of Δt=1×10-3. Using the proposed formulation we vary the time step from Δt=0.0625×10-3 to Δt=4×10-3 and found not instability for the total Lagrangian formulation. The adopted physical properties are E=200×109, v=0.25, ρ=10.000 and thickness of h=0.1.Figure 29
Geometry.Figure 30
Discretizations.Figure 31
Present solution for the apex displacement.Figure 32
Results presented by Argyris et al. [14].As one can see, the results converge to the smaller time step, while for large time steps the accuracy is lost, but the stability is kept.


## 10. Conclusions

A new, simple and robust formulation to solve dynamic geometrical nonlinear problems with large deflections applied to shells is proposed and implemented. The formulation is based on unconstrained vector mapping of the continuum, called here position description, simplifying the understanding and the implementation of total Lagrangian geometrical nonlinear analysis when compared to typical FEM shell formulations. The Newmarkβ method has been proved to be linear and angular momentum conserving, for the proposed total Lagrangian formulation. The high-order curved triangular element with improved transverse position field is free from locking and does not need reduced integration or relaxed constitutive relation to reproduce, with accuracy and small degrees of freedom, the static benchmarks of plate and shell analysis. Moreover, the formulation proved to be objective regarding rotation for unloaded and loaded structures. The general dynamic analysis of thin shells (airbag) and  the snap through benchmark indicate that the formulation is promising and should be extended to include physical nonlinearities (hyperelasticity and plasticity), fluid-structure iteration and impact.



---

*Source: 575131-2009-06-04.xml*