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.
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 letter A, given as follows:
where is the change of configuration function, as depicted in Figure 1, and represents variation regarding initial position.
In Figure 1 and represent an infinitesimal fiber in the initial and current continuum configurations, respectively. Following Ogden [17], the Green strain tensor can be written as
in which index notation is adopted. The variables and are the right Cauchy-Green stretch tensor and the Kroenecker delta, respectively. The following quadratic strain energy per unit of initial volume is adopted:
resulting into a linear elastic constitutive law relating second Piola-Kirchhoff stress and Green strain, usually called Saint-Venant-Kirchhoff elastic law, that is,
The elastic tensor is given by
where is the shear modulus, given by
with 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. Figure 2 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.
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: , , ,
where are nondimensional coordinates. The classical continuum mapping is written as
where are the coordinates of any point of the mapped continuum, are the shape functions, and are the coordinates of nodes named position parameters.
A totally similar mapping is given by
where the mid line nodal coordinates and nodal generalized vectors are given by
Expression (3.3) shows that the equivalent vector mapping is done using nonunitary vectors as parameters. In order to generalize the formulation one writes the nodal vectors and as functions of their lengths and resulting
where brackets mean that summation is not implied. The first term of (3.5) describes the reference line approximation. The values 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 (), a usual assumption for bars and shells, (3.5) turns into the desired continuum mapping, that is,
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 Figure 2 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:
where is the ith coordinate of a generic point in the mid surface of the shell at initial configuration, is the ith coordinate of node , is the ith coordinate of a generic point in the mid surface of the shell at current configuration, and is the ith coordinate of node at current configuration.
One can see in Figure 3 that is the positional mapping from the auxiliary space to the initial mid-surface configuration, is the positional mapping from the auxiliary space to the current mid-surface configuration, is the positional mapping from the initial configuration to the current one (not to be written) and the values , , 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 vectors or , see Figure 4.
A general point of the shell can be defined by adding the position vectors to the corresponding mid-surface point, that is,
Following what has been described in Section 3, for constant strain regarding , one writes and as functions of nondimensional coordinates, as
where , , 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 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 . 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
The introduced quadratic term allows linear strain variation along the transverse direction of the shell. The final generalized nodal vectors are not constrained, so the final thickness of the shell is not the same as the initial one and can be recovered as
The final mapping is generalized and unconstrained, like (3.6), improved in the membrane and thickness sense and is written as
in which the linear rate of thickness variation (scalar) is parameterized by its nodal values as follows:
The current position has seven unknown parameters per each node , that is, three positions , three generalized nodal vectors , and the nodal rate of thickness variation .
Function is used to find while function is used to find (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,
It is worth to show the derivatives of regarding the nondimensional variables , constituting the gradient as follows:
5. Dynamic Nonlinear Equation
From this section on, a unified notation will be adopted for nodal parameters, they will be called simply for each element and for varying from 1 to 7. The correspondence is as follows, translations for varying from 1 to 3, rate of thickness variation for 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:
where, following Lánczos [20], can be stated as the quantity of energy withdrawn from the simple conservative idealized energy . As a consequence, is the remaining (current) mechanical energy of the system, and (5.1) can be rewritten as
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 for equilibrium analysis.
For a structural problem associated with a fixed reference system, Figure 5, the ideal potential energy function can be written as the composition of the strain energy (), the potential energy of applied conservative forces P, the kinetic energy () and dissipation (), as follows:
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 , (2.3), as
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
where represents forces applied in direction "i" and is the ith current position of the point where the load is applied, is the distributed force applied in direction "i" and is the current position of mid-surface points ( 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 represents the initial differential area of elements. The kinetic energy is written as
where are velocities and is the mass density, relative to the initial volume . The dissipative term, including normal distributed forces, is written in its differential form as
where is the specific dissipative functional, is a proportional damping constant, are velocities at any point, and are components of the normal distributed forces given by
where is the normal distributed force over the element and the generalized vector corresponding to the values for 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 , for usual applications of thin structures the value of .
Substituting (5.4), (5.5) and (5.6) in (5.3) results
This energy function can be written substituting the exact position field by its approximation described in Section 4, that is,
The minimum potential energy theorem is used on by differentiating (5.10) regarding a generic nodal position , resulting
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
or
The involved forces are, inertial force or , damping or and the external force, divided into conservative and following forces . 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
Consequently
where 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 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 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
The linear momentum expression for a total Lagrangian description is given as
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,
From the continuity of the body, Ogden [17], one concludes that, the following equation is valid
Using (6.5), written for two times and , into (6.2), it results
or by continuity
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
where 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
Form the continuity assumption the equality
must hold for any point of the continuum. It occurs in two situations, the trivial undesired one, that is, and the desired one,
where is the angular velocity of the body around the rotation axis and is, without loss of generality, the position vector of the point related to its projection over the rotation axis. Using (6.11) for time into the Newmark (6.1) and (6.2), one writes
Rearranging terms of (6.12), one has
Substituting (6.14) into (6.13) results
Post-vector-multiplying (6.14) and (6.15) by and , respectively, and subtracting results one achieves
Finally (6.16) simplifies to
By continuity one writes
or in a shorter notation
Therefore, the Newmark method is angular momentum conserving for , 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
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 as follows:
Substituting the Newmark approximations (6.1) and (6.2) into (7.2) results
where vectors QS and RS represent the dynamic contribution of the past and are given by
Equation (7.3) can be understood simply by and is clearly nonlinear regarding . A Taylor expansion to solve this nonlinear equation is necessary. The second derivative of the total energy potential is then given by
One builds the Taylor series of first order as:
and derives the Newton-Raphson procedure to solve the nonlinear (7.3), that is,
where is a trial position (usually ) for used in (7.3) to calculate . Solving one calculates a new trial for as
The acceleration must be corrected for each iteration by an expression obtained from (6.1), that is,
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,
It must be noted that, before the first time step, the initial acceleration must be calculated as follows:
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
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
and omitting, for simplicity, extra indices, one applies the proposed mapping, that is, , and writes
Remembering that is constant regarding the current position, the first derivative is performed as
and, from (4.11), one has the following (not null) values of the current positional mapping gradient:
The first index of is the element node, and the second is the degree of freedom. The second derivative of the Cauchy-Green stretch is given by
where are given by (8.5) and the second derivatives are straightforward. Recalling (2.2), one achieves directly
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 and The thickness of the shell is . The adopted discretization can also be seen in Figure 6.
(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 Figure 7 one can see an illustration of this situation for the first turn. The adopted rotation step is .
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 Figure 8.
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 is . 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.
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 Figure 10. 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.
At three quarters of the first turn the stress values are the ones depicted in Figure 11, exactly as expected.
Finally, after a complete revolution the stress values depicted in Figure 12 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.
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.
The adopted physical data are , , , and varying from 0 to 0.5. The applied load is . 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.
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 Figure 15. The adopted discretization ( 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: , , , . Two values are adopted for the Poisson ratio, and , respectively, in order to check the locking free behavior of the presented formulation.
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 Figure 16 the results for 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 even if the linear rate of thickness variation is totally constrained. A formulation suffering of shear locking is not able to run this example.
In Figure 17, the behavior of the proposed formulation (with the seventh parameter constrained or not) is depicted for . The reference value for this figure is the seven parameter result with . 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.
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 is subjected to an initial velocity in the vertical direction of value . 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, . 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 is mandatory the other one is adopted as .
As one can (see in Figure 18) 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 Figure 19), is subjected to a field of initial velocity generated by an angular velocity of . The dimensions of the cylinder are: , and . The physical parameters are: , and . The total energy introduced in the problem by this velocity field is (considering the thickness influence) . Eighteen finite elements with forty eight nodes and 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.
The exact solution is coincident with the numerical one for . 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 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 Figure 21. It intends to show the transfer behavior of kinetic and strain energies for a deformable body. The beam has young Modulus , mass density , length , thickness , width and Poisson’s ratio . The applied velocity is proportional to the distance from the clamped end and has a maximum value of , see Figure 20. The adopted time step is . The exact total energy of the system is
Forty finite elements and 217 nodes are employed in the discretization. In Figure 22, 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].
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 , Poisson’s ratio of , and mass density of . 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 and the diameter in its flat initial configuration is . Our initial discretization ( 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:
where t is time and r is the distance from the center of the airbag. The adopted time step is and the maximum tolerance for convergence is in positions. An initial random vertical defect with maximum amplitude of is assumed for the analysis.
Some selected deformed positions are depicted in Figures 25, 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.
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 Figure 29. 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 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 and . 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 . Using the proposed formulation we vary the time step from to and found not instability for the total Lagrangian formulation. The adopted physical properties are , , and thickness of .
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.