- Email: [email protected]
Nonlinear static and dynamic analysis of shell structures with finite rotation
Computer methods in applied mechanics and engineering Comput. Methods Appl. Mech. Engrg. 162 (1998) 301-315
ELSEYIER
Nonlinear
static and dynamic analysis of shell structures with finite rotation J.L. Meek*, Yicai Wang Civil Engineering
Department.
The University
Received
of Queensland,
10 October
Brisbane.
Australia
1997
Abstract A simple, effective finite element incremental formulation and procedure for geometrically nonlinear static and dynamic analysis of a shell structure with finite rotations are presented, based on both the total Lagrangian and updated Lagrangian description of motion. The co-rotation formulation is used for finite rotation. The element employed here to implement the present method is a faceted shell element with Loof nodes (DKL + LST). An incremental iterative method based on the constant arc length method in conjunction with Newton-Raphson method is implemented in the present study for static and Newmark’s integration scheme for dynamic analysis. To demonstrate the accuracy and efficiency of the formulations and to compare the difference between the total and updated Lagrangian formulation, numerical studies are presented. 0 1998 Elsevier Science S.A. All rights reserved.
1. Introduction the total Lagrangian and updated Lagrangian For structural problems with geometric nonlinearities, formulations come as the two natural choices [l]. The total Lagrangian formulation has found wide application in problems involving geometrical nonlinearity and elastic stability. The updated Lagrangian formulation may be particularly useful for slender structures such as beams, plates and shells. The dominant factors in the geometrical nonlinearities of shells are attributed to large rigid body translational and rotational displacements coupled with elastic deformation while the strains remain small. If the rigid body motion part is eliminated from the total displacements the deformational part of the motion is always a small quantity relative to the local element axes, so that the faceted shell elements developed for the small displacement analysis may be applied to the nonlinear analysis of shell problems. Applying a co-rotation formulation it is possible to separate the rigid body displacements from the total deformation. The co-rotational formulation was presented by Argyris and Wempner [2,3].
2. Formulation 2. I. Nonlinear
description
The equilibrium displacements [I]:
* Corresponding
oj’ motion
of the body
at time
t + At is readily
author
0045-7825 /98/$19.00 0 1998 Else&r PII: SOO45-7825(97)00349-6
Science S.A. All rights reserved.
expressed
by means
the principle
of virtual
302
J.L. Meek, Y. Wang I Comput. Methods Appl. Mech. Engrg. 162 (1998) 301-315
I
(1)
“,+~, (r;,:A’)(S~;,:A’) dV= R’+”
where rb+Af are the components of Cauchy stress tensor at time t + At, and &;,+A’ are the corresponding virtual variations in strains. RrCAr is the external virtual work. The direct solution of Eq. (1) is not possible since the configuration at time t + At is not known. To obtain an approximate solution, Eq. (1) must be rewritten by refering all variables to a previously known calculated equilibrium configuration. In practice, the choice lies essentially between two formulations: (1) Total Lagrangian formulation where all the static and kinematic variables are referred to the initial configuration; (2) updated Lagrangian formulation where all variables refer to the updated previously calculated equilibrium position at time t. The governing nonlinear equilibrium equation for the total Lagrangian formulation is given:
I
v. (;+A’Si,)(8;+A’e~j) dV” = R’+Ar
(2)
where aA‘Sij and yA’eij are the Second Piola-Kirchhoff stresses and Green-Lagrange strains at time t + At, respectively. The governing nonlinear equilibrium equation for the updated Lagrangian formulation is given:
I
“~ (:+A’Sij)(S:+A’~;j) dV’ = Rt+”
(3)
where :+A’Sjj and :+A’e,j are the Second Piola-Kirchhoff stresses and Green-Lagrange strains, respectively, from the configuration t to the configuration t + At referred to the configuration at time t. 2.2. Finite element solution for large displacement For large displacements,
the strain-displacement
with moderate
rotation
matrix has the form:
1 4
=
y
('i.1
+
uj,i +
'k,i
The strain tensor components
lij = e:,
(4)
. uk,~)
are thus given in terms of linear and nonlinear
+ qij
terms: (5)
In Eq. (5), e:, represents the linear part and vii the nonlinear part of the strain tensor. Meek and Tan [4] developed an efficient and reliable triangular plate bending element (DKL+LST) used herein for large displacement theory, which has six nodes with three degrees of freedom: ui, u,, wi as displacements and six Loof nodes, each with one degree of rotation freedom, pi, along the edge in which node i is embedded. This flat shell element is illustrated in Fig. 1. Since small strain but finite rotation conditions are assumed, the membrane and bending contributions can be evaluated separately in the body-attached coordinate system, much in the same way as in linear analysis. The matrices are then transformed to the global coordinate system. The extension of linear theory ignores the shear strains and follows the Von Karman assumptions, so that, the nonlinear strain-displacement matrix for the faceted shell element has the form
(6)
where {e/1= represents
W”1 + Z[~bl>w~ = ~B,lW~ the linear part, and
(7)
J.L. Meek, Y. Wang 1 Comput. Methods Appl. Mech. Engrg. 162 (1998) 301-315
303
X
Fig. 1. Flat faceted shell element DKL+LST.
(8)
the nonlinear part of the strain tensor. From Eqs. (6), (7) and (8)
2.2.1. Total Lqrangian formulation In the total Lagrangian formulation, time t is written.
the 2nd Piola-Kirchhoff
stress tensor referred to the configuration
@I= [ml4
at the
(10)
From Eq. (2) the principle
of virtual displacements
is written
(11)
{F’}={Gu’} = /-“,)[S]={Se} dV Using Eqs. (9) and ( 10) finally
(12) {AF’} = Iv,,(([B,l= + 2[&JTXA9 + 2W”lTW> dV for total Lagrangian
formulation,
the incremental
stress-strain
(13) relationship
is:
WI = PI-&~ = rDw,l+
W,l)W'I
and the second term in Eq. (13) can be rewritten
i “0
(14) as
2[AB”]={S} dV= [K’,]{Au’}
Substituting
Eqs. (14) and (15) into Eq. (13) results in
(1%
J.L. Meek, Y. Wang I Comput. Methods Appl. Mech. Engrg. 162 (1998) 301-315
304
w;l + [K:l + W’,lHW
(16)
= if’> - 10
where [5]
[KFI = j-” ( t[B”]T[D][Bm] + & ,Bb]T,Dl,Bbl > dA r
[K;] = j, 2t([B”lT[Dl[Bnl + [B”lT~~l[~“l + 2[~“lTPlUf’l)
1 VJ =I,o K’,l =
A”
PI PI
[Ol PI
PI to1
PI WI
[o]
[o]
La]
lo]
PI
WI
WI [Ol
4
dA
dA
+ z[Bb]= + 2[BJT>{S’} dV
w7=
[KY], [Kz] and [K’,] are, respectively, the element linear, nonlinear and geOmetriC Stiffness v’} is the external force and the vector {re} is the element nodal force vector. The vector {S’} is the 2nd Piola-Kirchhoff stress at time t, and
The matrices ma&es.
The
vector
[aI = q~~,,1kh,.J=
+ ~“[~J[~W,,lT
+ ~x,wJqvl
The stress resultants [N] = [N,, I$, N,,] are obtained by integrating the stresses through the thickness of the element and they are in general related to the 2nd Piola-Kirchhoff stresses at time t. Use of the total Lagrangian formulation implies that the reference configuration for the field variables is stationary, and as such the local-global transformation matrix, [T] does not change during the deformation process. 2.2.2. Updated Lqrangian formulation In the updated Lagrangian formulation, the 2nd Piola-Kirchhoff the time t + At is expressed in terms of the two components: $A
stress tensor referred to the configuration
= r:, + ASIj
(17)
where [r’] is the Cauchy stress at time t. Notice, for the updated Lagrangian
method, {AE} = {e} and {Au} = {u},
(18)
{AS] = [WWJ,l + PW~W Substituting
at
Eqs. (9), (17) and (18) into Eq. (3)
(19)
(W;l + KZI + K’,l)@u’) = if’) - be> where
K;l = I(A’
t[B”]=[D][B”]
+
g
dA
,BblT,D][Bb] >
r
[K;] = j,, t([B”]=[D][B”]
K',l
=
A’
WI PI
Ku PI
[o]
[a]
f”]
[Ol
PI
[o]
{re} = Iv, ([fIrnIT +
The mat&es
[K;],
2[B”l=[~1[@7 + 2[B”1T[~1[B”1>dA
WI PI WI [Ol LOI to1
4
+
[Kz]
1 dA
z[BblT>[~‘l dV and [K’,]
are, respectively,
the element
linear,
nonlinear
and geOIIU?tIiC Stiffness
J.L. Meek, Y.
Wang I Compur.Methods Appl. Mech. Engrg. 162 (1998) 301-315
305
matrices. The vector {re} is the element nodal force vector, and the stress resultants [N] are related to the Cauchy stresses. The use of the updated Lagrangian formulation implies that the reference coordinate system is defined by the element coordinate system in the position at time t. Therefore, the local stiffness matrices as well as force and displacement vectors referred to this coordinate system must be transformed to the global coordinate system prior to assemblage. Thus, the local-global transformation matrix [T] needs to be changed after every step of the displacement. Also, the transformation between Piola-Kirchhoff and Cauchy stresses must be performed in every step. 2.3, Finite element solution for large displacement
with finite rotation
The use of the Von Karman assumptions to formulate the nonlinear characteristics of the faceted shell element means that the in-plane displacements u and u are infinitesimal so that their derivatives are small and hence their quadratic variations can be neglected. Thus only the second degree terms in the deflection w are retained. Unfortunately, in geometric nonlinear analysis, the dominant effect of the nonlinearities can be attributed to the finite rotations with small strains. This implies that the motion of an individual element will consist of rigid body motion. When the total or updated Lagrangian formulation in conjunction with Von Karman assumptions is employed, the solution of the governing equations obtained may result in large displacements and rotations which in turn give strains incorrectly because the rigid body motion does not affect stresses/strains in the structure. Fig. 2 shows a vector r which as a result of application of a rotation vector n is transported to a new position r’. The relation between r and r’ may be expressed as [6]: r’=cos@+(l-cos~)(n*r)n+sin#QzXr)
(20)
If it is assumed that the change is infinitesimal between the current and the previous position, the relation between the vectors r’ and r can be written as {r’} = [l + E](T)
(21)
where {E} is the infinitesimal operator which is defined by the antisymmetric matrix:
Thus, the change of the components of a vector can be represented as {dr} = {r’} - {r} = E(T)
(22)
That is {d4 = [dfil x k>
(23)
Fig. 2. Rotation
formula.
J.L. Meek, Y. Wang I Compui. Methods Appl. Mech. Engrg. 162 (1998) 301-315
306
Also, Eq. (20) becomes k’) = k> + in> X {r) d@
(24)
&I = -h>X {r) d4
(25)
Comparing Eqs. (23) and (25), it is seen that
[dfll = b> d4
(26)
2.3.1. Total Lagrangian formulation In the case of the total Lagrangian formulation, where the reference configuration is always the initial state, the use of the Von I&man assumptions for handling the large rotations is incorrect since the displacement of the element must be so restricted that the element only experiences moderate rotations. If the rigid body motion part can be eliminated from the total displacements, the deformational part of the element is always small relative to the local element axes. Thus, in the context of an incremental iterative procedure, the element internal forces are calculated using these incremental deformational displacements and incremental mid-side rotations. The incremental deformation of the element is calculated from the difference in the local element coordinate system at the current and the initial configuration illustrated in Fig. 3. This assumption imposes a restriction on the maximum element size that can be employed. Incremental
membrane
deformations
From Fig. 3 it is seen that {Xi’“‘} = {X’,‘“‘}+ [T]{x;+*‘}
(27)
where {x:+*~}1s the coordinate vector of node i in the local coordinate system at time t + At. That is
{x:+*‘}= [T]T({X:+A’}
-{X’,+“‘})
(28)
and in which i = 2, . . . ,6 is the node number. The previously defineh element is then rigidly moved to this updated position and the local incremental membrane deformations are calculated as
(29)
Fig. 3. Incremental
membrane
deformations
for T.L.
J.L. Meek, Y. Wang J Comput. Methods Appl. Mech. Engrg.
162 (1998) 301-315
307
(30)
These local incremental
in-plane
deformations
of the element
Incremental bending deformations The incremental bending deformations include and the incremental rotational deformations. (A) Incremental transverse Similar to the membrane
deformations deformations,
are shown in Fig. 3.
two parts, namely,
the vector of transverse
the incremental
deformations
transverse
deformations
is given by
(31)
(B) Incremental rotational deformations If the element incremental rotation matrix is defined as
e,&, -exP(il
0
1
(32)
where {e} = [e, eY eZIT is the direction cosine of the side associated with node i, and &, is the element incremental rotation associated with node i obtained from the solution of the incremental equilibrium equations Eq. (16) a joint orientation matrix [(w] is defined for each node with a rotation d.o.f. It is initially assumed to be parallel to the global axes and hence
[LUO]=
1 0 [0
0 1 0
0 0 1
1
(33)
[da’] = [dL?][cu’]
(34)
]ff ‘+“] = [(Y’] + [da’]
(35)
Therefore,
the rotational
deformation
of the node i is
sin p, = -{e}‘{p”‘} where the matrix rotations.
[p] describing
[p] = [a’+“][r”]
(36) the rotational
deformation
of the element
at node i contains
the rigid body
(37)
and [r’] is the matrix of direction cosines of the triad of orthogonal unit vectors [4]. As mentioned above, the local updated x - y axes pass through the three corner nodes, so only third column of [p] is used. For moderate rotations it is possible to use the approximation p, = -{e}‘{p’“}
(38)
308
J.L. Meek. Y. Wang I Comput. Methods Appl. Mech. Engrg. 162 (1998) 301-315
2.3.2. Updated Lagrangian formulation Application of the co-rotational updated Lagrangian formulation requires updating of the nodal coordinates for each new position of the element. The incremental deformation of the element is calculated from the difference in the local element coordinate system at the current and the previous configuration shown in Fig. 4. As with the total Lagrangian formulation, it is assumed that the membrane and bending deformations are uncoupled. Then, the calculation of the incremental nodal deformations can be divided into those for membrane deformations and for bending deformations. The bending deformations are discussed in terms of transverse nodal deformations and deformational nodal rotations. Incremental membrane deformations From Fig. 4 it is seen that {Xi’“‘} = {x;‘“‘} + [Tr+Ar]{x;+A’} where {.x~+~‘}IS ’ th e coordinate
(39)
vector of node i in the local coordinate
system at time t + At. That is
{x;+AI} = [Tf+A’]T({X;+A’} - {Xf’“‘})
(40)
and in which i = 2, . . . ,6 is the node number. The previously defined element is then rigidly updated position and the local incremental membrane deformations are calculated as ‘u,’
0 r+Ar
r
X2 -x2 r+Ar f x3 -x3 r+Ar I x4 --x4 r+Ar I x5 --x5 r+A.r f .x6 -x6.
u2 u3
{u”} = U4 u5 '6.
moved to this
(41)
0 0
=
r+Ar I Y3 -Y3 t+Ar f Y4 -Y4 t+Ar f Y5 -Y5 r+Ar .Y, -Y:,
(42)
These local incremental in-plane deformations of the element are given in Fig. 4. By computing the local membrane deformations in this manner, the rigid body motion of the element is automatically removed from the displacement vectors un and II”.
Fig. 4. Incremental
membrane
deformations
for U.L
J.L. Meek, Y. Wang / Comput. Methods Appl. Mech. Engrg. 162 (1998) 301-315
Incremental bending deformations The incremental bending deformations include and the incremental rotational deformations. (A) Incremental transverse Similar to the membrane
deformations deformations,
two parts, namely,
the vector of transverse
the incremental
deformations
transverse
309
deformations
is given by
(43)
Due to the local updated x - y coordinate axes being defined as passing through the three comer nodes, the local transverse displacements at these nodes are zero. (B) Incremental rotational deformations As with the total Lagrangian formulation,
[da,
=[;;?
;
if the element
incremental
rotation
matrix is defined as
+X]
(44)
where cu,, aY and LY,are the element incremental rigid body rotations results in a change of the element vector r, according to Eq. (24):
about x, y and z axes, respectively.
Firstly, the t-element
conditions
is rotated about z axis by setting the following
This
in Eq. (45)
a; = 0 a,, = 0 $2 = 0 and this results in the calculation a:=--
of the element
rotation
about z axis as
Yl2 X12
Similarly
it is possible
to calculate
the rotations
about y and x axes as (47)
(48) Since the final nodal parameters of the element include the rotations along the sides, it is therefore necessary to project the vector of rigid rotation onto the sides of the triangle and from geometrical considerations this gives the incremental rigid body rotation along side k as (49) where the indexes
i, j, k = 1,2,3
correspond
to the three sides of the triangle
and
J.L. Meek, Y. Wang I Comput. Methods Appl. Mech. Engrg. 162 (1998) 301-315
Then the incremental rotational deformation at Loof node i is obtained by subtracting this part of the rigid body rotation ‘Ye,from the calculated rotational displacement at Loof node i from equilibrium equations, PC,, Pi=&,-~kr
i=l,2,3
2.4. Step by step integration
,....,
6
k=l,2,3.
(50)
method
In the current study, the Newmark integration scheme is used to analyse the nonlinear The flow chart for the above integration procedure is presented in Fig. 5.
dynamic behaviour
givenattimet=O: displacements r0 velocities +O I accelerations: f. = M-l F(O) - Cko - R(rO)]
initialisation of accelerations: fk = EL-1 I displacements and velocities (predictor): rk = r&l + At&l + Ata(4- CY) iike1 + At%& ik = ik-1 + At (I- 6) fk-1 + Lit&k
-
I
t
dynamic out-of-balance forces: R&k) = F - Mfk - Cir, - R(Q)
diiplacement correction: [-g + &M + AC + K&h
I
=
%(rk) 1
I
t
adjusted displacements, velocities and accelerations (corrector): rk = rk + Ark +k = fk =
ik
+
fk
+
’ Ark ;;pi -Ark
Fig. 5. Flow chart of the Newmark
procedures
applied to nonlinear
dynamic
analyses.
[7].
J.L. Meek, Y. Wang I Compui. Methods Appl. Mech. Engrg. 162 (1998) 301-315
311
3. Examples The numerical examples that are presented in this section demonstrate the capability of the formulation and solution methods used in this paper. Also, the difference between the total Lagrangian and the updated Lagrangian formulation is shown to be as follows. Using the updated Lagrangian formulation can take larger load step than that of the total Lagrangian formulation, especially for slender structures. The structures were discretized by the finite element method with the DKL+LST element developed by Meek and Tan [4], using a basic rectangular element composed of four sub-triangles.
3.1. Examples for static analysis 3.1.1. Hinged cylindrical panel under a concentrated load The first example considers the large displacement behaviour of a cylindrical panel under a central vertical load. The cylindrical panel depicted in Fig. 6 has been analysed by various researchers for different loading and boundary conditions, well into the advanced nonlinear range. As shown in Fig. 6, the radius is R = 2540.0 mm, thickness h = 6.35 mm, length L = 2 X 2540.0 mm. The material properties are E = 3.10275 kN/mm*, v = 0.3. Taking into account the symmetry of geometry, boundary conditions and loading, one quarter of the panel was modelled by finite elements. The mesh used of 65 nodes and 16 elements. Fig. 7 shows the load-deflection behaviour of the panel at load point for both the total and updated Lagrangian formulations. Fig. 8 shows the the shapes of the cylindrical panel at different loads. Both total and updated Lagrangian formulation used a load increment of 0.01 kN.
‘” I o.I1
R
Fig. 6. Cylindrical
(a) Load=O.kN
shell.
Fig. 7. Load-displacement
(b) Load=-0.5687kN
(c) Load=O.O29kN
Fig. 8. The shapes of panel at different
loads.
curve of the panel shell at central point.
(d) Load=O.3616kN
J.L. Meek, Y. Wang I Comput. Methods Appl. Mech. Engrg. 162 (1998) 301-315
312
h
M
Js: b
L
U lW
Fig. 9. A cantilever
beam subjected
Fig. 10. Load-displacement
to an end moment.
curve for the cantilever
beam.
3.1.2. Cantilever beam under an end moment
This classical elastic problem has been widely used by many investigators as a benchmark test for large deformation analysis. As shown in Fig. 9, the geometrical properties of the beam are length L = 10.0 m, width b = 1.0 m and thickness h = 0.1 m. The Young’s modulus of the beam is E = 1.2 X lo5 N/m2 with Poisson’s ratio Y = 0.0. The load increments AP = 0.01 kN. Analytical solutions show that under end moment M, the beam rolls up into a circular arc of radius p: p=~=
L
1.59155m
(a) Lload=o.kN
(b) Load=3.47kN Fig. Il.
(c) Load=6.666kN
The shapes of the cantilever
beam at different
loads.
Fig. 12. Load-displacement
curve for T.L. with 20 elements.
Fig. 13. Load-displacement
curve for T.L. with 40 elements.
(d) Load=6.538?‘kN
J.L. Meek, Y. Wang I Comput. Methods Appl. Mech. Engrg. 162 (1998) 301-315
Fig. 14. Step load.
313
Fig. 15. Time history of the panel shell at central point.
In the finite element discretization, a 103 node, 20 element mesh was used. For the updated Lagrangian formulation, the load increments were AP = 0.01 kN. Figs. 10 and 11 include results obtained from the present investigation. These indicate that the beam rolls up to a complete circle when M = 6.5 kN. In this configuration, a maximum lateral displacement of 7.247 m is observed at the mid point of the beam. For total Lagrangian formulation, when the element number is taken 20 and 40 the results are not convergent at load level of 3.78 kN and 5.62 kN, respectively. They are shown in Figs. 12 and 13. 3.2. Examples of dynamic analysis 3.2.1. Hinged cylindrical panel under a concentrated load This cylindrical panel has been studied in 3.1.1 and in here E = 0.310275 kN/mm’, mass density p = 3210.05 kg/m3. The apex load is applied as a step load shown in Fig. 14. P, = 400 N and t, = 0.01 s. The time step size At = 0.001 s. Fig. 15 shows the dynamic response curve of the panel at the load point and Fig. 16 shows the shapes of the cylindrical panel at different times. It can be seen that, at time = 0.048 s, the displacement first reaches the maximum value which is -37.4 mm, the shape of the panel becomes a concave shell. The displacements obviously include finite rigid rotations.
(a) Time=O.sec.
(b) Tie=O.O4sec.
(c) Time=O.OSsec.
(d) Time=O.l2sec.
(e) Time=O.%ec.
(f) Time=O.2aec.
Fig. 16. The shapes of panel at different
times.
J.L. Meek, Y. Wang I Comput. Methods Appl. Mech. Engrg. 162 (1998) 301-315
314
Fig. 17. The rectangular
plate.
Fig. 18. Time history
of centre displacement
of the rectangular
plate.
(b) Time=O.O8sec.
(a) Time=O.aee.
(c) Tiie=O.l6sec.
Fig. 19. The shapes of plate at different
3.2.2.
Simply supported
(d) Tiie=0.32sec.
times.
plate
This example considers the geometric nonlinear dynamic analysis of a rectangular plate with simply supported edges. The plate is shown in Fig. 17. Its dimensions are a = 1016 mm, b = 1524 mm and thickness h = 25.4 mm. The material properties are: Young’s modulus E = 2.0955 X lo8 N/m*, Poisson’s ratio v = 0.25, and mass density per unit volume p = 3210.05 kg/m3. Owing to symmetry of geometry, boundary condition and loading, one quarter of the plate was discretized. For direct integration the time step size was 0.004 s, twice that suggested in [8]. The step function shown in Fig. 14 was considered, with PO = 44.54 N and t, = 0.006 s. A linear and nonlinear dynamic analysis was performed. The results are plotted in Figs. 18 and 19. Comparing the results with those of Wilson et al. [8], very good agreement is observed. The nonlinear results had smaller amplitude and shorter period compared with the linear solution, which indicates stiffness hardening due to the inplane forces and deformed shape.
4. Conclusions
The objective in this study was to derive and evaluate finite element formulations for general nonlinear static and dynamic analysis which have been implemented in the search for the most effective procedure. A simple, practical formulation and program to remove the restriction of small rotations for geometrically nonlinear dynamic analysis of shells using flat triangular shell elements are presented. The conceptual difference between the formulations is the reference configuration that is used for the linearization of the incremental equations of motion. In the total Lagrangian formulation the initial configuration is used as reference, whereas in the updated Lagrangian formulations, the reference configuration corresponds to the last calculated configuration. An advantage of the total Lagrangian formulation is that the derivatives of the interpolation functions are with respect to the initial configuration, and therefore need be formed only once. On the other hand, for updated Lagrangian formulation, the convergence rates are better than that of total Lagrangian formulation, especially for slender structures. It is very clearly shown in Example 3.1.2.
J.L. Meek, Y. Wang I Compuf. Methods Appl. Mech. Engrg.
162 (1998) 301-31.5
315
References [l] K.J. Bathe, E. Rramm and E.L. Wilson, Finite element formulations for large deformation dynamic analysis, Int. J. Numer. Methods Engrg. 9 (1975) 353-386. [2] J.H. Argyris and P.C. Dunne, On the application of the natural node technique to small strain large displacement problems, World Congress on Finite Element Methods in Structural Mechanics, Boumemouth, 1975. [3] G. Wempner, Finite elements finite rotations and small strains of flexible shells, Int. J. Solids Struct. 5 (1969) 117-153. [4] J.L. Meek and H.S. Tan, A faceted shell element with Loof nodes, hit. J. Numer. Methods Engrg. 23 (1986) 49-67. [5] J.L. Meek and H.S. Tan, A discrete Kirchgoff Plate bending element with Loof nodes, Comput. Stmct. 21 (1985) 1197-1212. [6] H. Goldstein, Classical Mechanics (Addison-Wesley, 1980). [7] J. Argyris and HP Mlejnek. Dynamics of Structures (North-Holland, Amsterdam, New York, Oxford, Tokyo, 1991). [8] E.L. Wilson, K.J. Bathe and R.H. Iding, A structural analysis program for static and dynamic response of nonlinear systems, SESM 74-3, Structural Engineering Laboratory, University of California, Berkeley, CA, 1974.
ELSEYIER
Nonlinear
static and dynamic analysis of shell structures with finite rotation J.L. Meek*, Yicai Wang Civil Engineering
Department.
The University
Received
of Queensland,
10 October
Brisbane.
Australia
1997
Abstract A simple, effective finite element incremental formulation and procedure for geometrically nonlinear static and dynamic analysis of a shell structure with finite rotations are presented, based on both the total Lagrangian and updated Lagrangian description of motion. The co-rotation formulation is used for finite rotation. The element employed here to implement the present method is a faceted shell element with Loof nodes (DKL + LST). An incremental iterative method based on the constant arc length method in conjunction with Newton-Raphson method is implemented in the present study for static and Newmark’s integration scheme for dynamic analysis. To demonstrate the accuracy and efficiency of the formulations and to compare the difference between the total and updated Lagrangian formulation, numerical studies are presented. 0 1998 Elsevier Science S.A. All rights reserved.
1. Introduction the total Lagrangian and updated Lagrangian For structural problems with geometric nonlinearities, formulations come as the two natural choices [l]. The total Lagrangian formulation has found wide application in problems involving geometrical nonlinearity and elastic stability. The updated Lagrangian formulation may be particularly useful for slender structures such as beams, plates and shells. The dominant factors in the geometrical nonlinearities of shells are attributed to large rigid body translational and rotational displacements coupled with elastic deformation while the strains remain small. If the rigid body motion part is eliminated from the total displacements the deformational part of the motion is always a small quantity relative to the local element axes, so that the faceted shell elements developed for the small displacement analysis may be applied to the nonlinear analysis of shell problems. Applying a co-rotation formulation it is possible to separate the rigid body displacements from the total deformation. The co-rotational formulation was presented by Argyris and Wempner [2,3].
2. Formulation 2. I. Nonlinear
description
The equilibrium displacements [I]:
* Corresponding
oj’ motion
of the body
at time
t + At is readily
author
0045-7825 /98/$19.00 0 1998 Else&r PII: SOO45-7825(97)00349-6
Science S.A. All rights reserved.
expressed
by means
the principle
of virtual
302
J.L. Meek, Y. Wang I Comput. Methods Appl. Mech. Engrg. 162 (1998) 301-315
I
(1)
“,+~, (r;,:A’)(S~;,:A’) dV= R’+”
where rb+Af are the components of Cauchy stress tensor at time t + At, and &;,+A’ are the corresponding virtual variations in strains. RrCAr is the external virtual work. The direct solution of Eq. (1) is not possible since the configuration at time t + At is not known. To obtain an approximate solution, Eq. (1) must be rewritten by refering all variables to a previously known calculated equilibrium configuration. In practice, the choice lies essentially between two formulations: (1) Total Lagrangian formulation where all the static and kinematic variables are referred to the initial configuration; (2) updated Lagrangian formulation where all variables refer to the updated previously calculated equilibrium position at time t. The governing nonlinear equilibrium equation for the total Lagrangian formulation is given:
I
v. (;+A’Si,)(8;+A’e~j) dV” = R’+Ar
(2)
where aA‘Sij and yA’eij are the Second Piola-Kirchhoff stresses and Green-Lagrange strains at time t + At, respectively. The governing nonlinear equilibrium equation for the updated Lagrangian formulation is given:
I
“~ (:+A’Sij)(S:+A’~;j) dV’ = Rt+”
(3)
where :+A’Sjj and :+A’e,j are the Second Piola-Kirchhoff stresses and Green-Lagrange strains, respectively, from the configuration t to the configuration t + At referred to the configuration at time t. 2.2. Finite element solution for large displacement For large displacements,
the strain-displacement
with moderate
rotation
matrix has the form:
1 4
=
y
('i.1
+
uj,i +
'k,i
The strain tensor components
lij = e:,
(4)
. uk,~)
are thus given in terms of linear and nonlinear
+ qij
terms: (5)
In Eq. (5), e:, represents the linear part and vii the nonlinear part of the strain tensor. Meek and Tan [4] developed an efficient and reliable triangular plate bending element (DKL+LST) used herein for large displacement theory, which has six nodes with three degrees of freedom: ui, u,, wi as displacements and six Loof nodes, each with one degree of rotation freedom, pi, along the edge in which node i is embedded. This flat shell element is illustrated in Fig. 1. Since small strain but finite rotation conditions are assumed, the membrane and bending contributions can be evaluated separately in the body-attached coordinate system, much in the same way as in linear analysis. The matrices are then transformed to the global coordinate system. The extension of linear theory ignores the shear strains and follows the Von Karman assumptions, so that, the nonlinear strain-displacement matrix for the faceted shell element has the form
(6)
where {e/1= represents
W”1 + Z[~bl>w~ = ~B,lW~ the linear part, and
(7)
J.L. Meek, Y. Wang 1 Comput. Methods Appl. Mech. Engrg. 162 (1998) 301-315
303
X
Fig. 1. Flat faceted shell element DKL+LST.
(8)
the nonlinear part of the strain tensor. From Eqs. (6), (7) and (8)
2.2.1. Total Lqrangian formulation In the total Lagrangian formulation, time t is written.
the 2nd Piola-Kirchhoff
stress tensor referred to the configuration
@I= [ml4
at the
(10)
From Eq. (2) the principle
of virtual displacements
is written
(11)
{F’}={Gu’} = /-“,)[S]={Se} dV Using Eqs. (9) and ( 10) finally
(12) {AF’} = Iv,,(([B,l= + 2[&JTXA9 + 2W”lTW> dV for total Lagrangian
formulation,
the incremental
stress-strain
(13) relationship
is:
WI = PI-&~ = rDw,l+
W,l)W'I
and the second term in Eq. (13) can be rewritten
i “0
(14) as
2[AB”]={S} dV= [K’,]{Au’}
Substituting
Eqs. (14) and (15) into Eq. (13) results in
(1%
J.L. Meek, Y. Wang I Comput. Methods Appl. Mech. Engrg. 162 (1998) 301-315
304
w;l + [K:l + W’,lHW
(16)
= if’> - 10
where [5]
[KFI = j-” ( t[B”]T[D][Bm] + & ,Bb]T,Dl,Bbl > dA r
[K;] = j, 2t([B”lT[Dl[Bnl + [B”lT~~l[~“l + 2[~“lTPlUf’l)
1 VJ =I,o K’,l =
A”
PI PI
[Ol PI
PI to1
PI WI
[o]
[o]
La]
lo]
PI
WI
WI [Ol
4
dA
dA
+ z[Bb]= + 2[BJT>{S’} dV
w7=
[KY], [Kz] and [K’,] are, respectively, the element linear, nonlinear and geOmetriC Stiffness v’} is the external force and the vector {re} is the element nodal force vector. The vector {S’} is the 2nd Piola-Kirchhoff stress at time t, and
The matrices ma&es.
The
vector
[aI = q~~,,1kh,.J=
+ ~“[~J[~W,,lT
+ ~x,wJqvl
The stress resultants [N] = [N,, I$, N,,] are obtained by integrating the stresses through the thickness of the element and they are in general related to the 2nd Piola-Kirchhoff stresses at time t. Use of the total Lagrangian formulation implies that the reference configuration for the field variables is stationary, and as such the local-global transformation matrix, [T] does not change during the deformation process. 2.2.2. Updated Lqrangian formulation In the updated Lagrangian formulation, the 2nd Piola-Kirchhoff the time t + At is expressed in terms of the two components: $A
stress tensor referred to the configuration
= r:, + ASIj
(17)
where [r’] is the Cauchy stress at time t. Notice, for the updated Lagrangian
method, {AE} = {e} and {Au} = {u},
(18)
{AS] = [WWJ,l + PW~W Substituting
at
Eqs. (9), (17) and (18) into Eq. (3)
(19)
(W;l + KZI + K’,l)@u’) = if’) - be> where
K;l = I(A’
t[B”]=[D][B”]
+
g
dA
,BblT,D][Bb] >
r
[K;] = j,, t([B”]=[D][B”]
K',l
=
A’
WI PI
Ku PI
[o]
[a]
f”]
[Ol
PI
[o]
{re} = Iv, ([fIrnIT +
The mat&es
[K;],
2[B”l=[~1[@7 + 2[B”1T[~1[B”1>dA
WI PI WI [Ol LOI to1
4
+
[Kz]
1 dA
z[BblT>[~‘l dV and [K’,]
are, respectively,
the element
linear,
nonlinear
and geOIIU?tIiC Stiffness
J.L. Meek, Y.
Wang I Compur.Methods Appl. Mech. Engrg. 162 (1998) 301-315
305
matrices. The vector {re} is the element nodal force vector, and the stress resultants [N] are related to the Cauchy stresses. The use of the updated Lagrangian formulation implies that the reference coordinate system is defined by the element coordinate system in the position at time t. Therefore, the local stiffness matrices as well as force and displacement vectors referred to this coordinate system must be transformed to the global coordinate system prior to assemblage. Thus, the local-global transformation matrix [T] needs to be changed after every step of the displacement. Also, the transformation between Piola-Kirchhoff and Cauchy stresses must be performed in every step. 2.3, Finite element solution for large displacement
with finite rotation
The use of the Von Karman assumptions to formulate the nonlinear characteristics of the faceted shell element means that the in-plane displacements u and u are infinitesimal so that their derivatives are small and hence their quadratic variations can be neglected. Thus only the second degree terms in the deflection w are retained. Unfortunately, in geometric nonlinear analysis, the dominant effect of the nonlinearities can be attributed to the finite rotations with small strains. This implies that the motion of an individual element will consist of rigid body motion. When the total or updated Lagrangian formulation in conjunction with Von Karman assumptions is employed, the solution of the governing equations obtained may result in large displacements and rotations which in turn give strains incorrectly because the rigid body motion does not affect stresses/strains in the structure. Fig. 2 shows a vector r which as a result of application of a rotation vector n is transported to a new position r’. The relation between r and r’ may be expressed as [6]: r’=cos@+(l-cos~)(n*r)n+sin#QzXr)
(20)
If it is assumed that the change is infinitesimal between the current and the previous position, the relation between the vectors r’ and r can be written as {r’} = [l + E](T)
(21)
where {E} is the infinitesimal operator which is defined by the antisymmetric matrix:
Thus, the change of the components of a vector can be represented as {dr} = {r’} - {r} = E(T)
(22)
That is {d4 = [dfil x k>
(23)
Fig. 2. Rotation
formula.
J.L. Meek, Y. Wang I Compui. Methods Appl. Mech. Engrg. 162 (1998) 301-315
306
Also, Eq. (20) becomes k’) = k> + in> X {r) d@
(24)
&I = -h>X {r) d4
(25)
Comparing Eqs. (23) and (25), it is seen that
[dfll = b> d4
(26)
2.3.1. Total Lagrangian formulation In the case of the total Lagrangian formulation, where the reference configuration is always the initial state, the use of the Von I&man assumptions for handling the large rotations is incorrect since the displacement of the element must be so restricted that the element only experiences moderate rotations. If the rigid body motion part can be eliminated from the total displacements, the deformational part of the element is always small relative to the local element axes. Thus, in the context of an incremental iterative procedure, the element internal forces are calculated using these incremental deformational displacements and incremental mid-side rotations. The incremental deformation of the element is calculated from the difference in the local element coordinate system at the current and the initial configuration illustrated in Fig. 3. This assumption imposes a restriction on the maximum element size that can be employed. Incremental
membrane
deformations
From Fig. 3 it is seen that {Xi’“‘} = {X’,‘“‘}+ [T]{x;+*‘}
(27)
where {x:+*~}1s the coordinate vector of node i in the local coordinate system at time t + At. That is
{x:+*‘}= [T]T({X:+A’}
-{X’,+“‘})
(28)
and in which i = 2, . . . ,6 is the node number. The previously defineh element is then rigidly moved to this updated position and the local incremental membrane deformations are calculated as
(29)
Fig. 3. Incremental
membrane
deformations
for T.L.
J.L. Meek, Y. Wang J Comput. Methods Appl. Mech. Engrg.
162 (1998) 301-315
307
(30)
These local incremental
in-plane
deformations
of the element
Incremental bending deformations The incremental bending deformations include and the incremental rotational deformations. (A) Incremental transverse Similar to the membrane
deformations deformations,
are shown in Fig. 3.
two parts, namely,
the vector of transverse
the incremental
deformations
transverse
deformations
is given by
(31)
(B) Incremental rotational deformations If the element incremental rotation matrix is defined as
e,&, -exP(il
0
1
(32)
where {e} = [e, eY eZIT is the direction cosine of the side associated with node i, and &, is the element incremental rotation associated with node i obtained from the solution of the incremental equilibrium equations Eq. (16) a joint orientation matrix [(w] is defined for each node with a rotation d.o.f. It is initially assumed to be parallel to the global axes and hence
[LUO]=
1 0 [0
0 1 0
0 0 1
1
(33)
[da’] = [dL?][cu’]
(34)
]ff ‘+“] = [(Y’] + [da’]
(35)
Therefore,
the rotational
deformation
of the node i is
sin p, = -{e}‘{p”‘} where the matrix rotations.
[p] describing
[p] = [a’+“][r”]
(36) the rotational
deformation
of the element
at node i contains
the rigid body
(37)
and [r’] is the matrix of direction cosines of the triad of orthogonal unit vectors [4]. As mentioned above, the local updated x - y axes pass through the three corner nodes, so only third column of [p] is used. For moderate rotations it is possible to use the approximation p, = -{e}‘{p’“}
(38)
308
J.L. Meek. Y. Wang I Comput. Methods Appl. Mech. Engrg. 162 (1998) 301-315
2.3.2. Updated Lagrangian formulation Application of the co-rotational updated Lagrangian formulation requires updating of the nodal coordinates for each new position of the element. The incremental deformation of the element is calculated from the difference in the local element coordinate system at the current and the previous configuration shown in Fig. 4. As with the total Lagrangian formulation, it is assumed that the membrane and bending deformations are uncoupled. Then, the calculation of the incremental nodal deformations can be divided into those for membrane deformations and for bending deformations. The bending deformations are discussed in terms of transverse nodal deformations and deformational nodal rotations. Incremental membrane deformations From Fig. 4 it is seen that {Xi’“‘} = {x;‘“‘} + [Tr+Ar]{x;+A’} where {.x~+~‘}IS ’ th e coordinate
(39)
vector of node i in the local coordinate
system at time t + At. That is
{x;+AI} = [Tf+A’]T({X;+A’} - {Xf’“‘})
(40)
and in which i = 2, . . . ,6 is the node number. The previously defined element is then rigidly updated position and the local incremental membrane deformations are calculated as ‘u,’
0 r+Ar
r
X2 -x2 r+Ar f x3 -x3 r+Ar I x4 --x4 r+Ar I x5 --x5 r+A.r f .x6 -x6.
u2 u3
{u”} = U4 u5 '6.
moved to this
(41)
0 0
=
r+Ar I Y3 -Y3 t+Ar f Y4 -Y4 t+Ar f Y5 -Y5 r+Ar .Y, -Y:,
(42)
These local incremental in-plane deformations of the element are given in Fig. 4. By computing the local membrane deformations in this manner, the rigid body motion of the element is automatically removed from the displacement vectors un and II”.
Fig. 4. Incremental
membrane
deformations
for U.L
J.L. Meek, Y. Wang / Comput. Methods Appl. Mech. Engrg. 162 (1998) 301-315
Incremental bending deformations The incremental bending deformations include and the incremental rotational deformations. (A) Incremental transverse Similar to the membrane
deformations deformations,
two parts, namely,
the vector of transverse
the incremental
deformations
transverse
309
deformations
is given by
(43)
Due to the local updated x - y coordinate axes being defined as passing through the three comer nodes, the local transverse displacements at these nodes are zero. (B) Incremental rotational deformations As with the total Lagrangian formulation,
[da,
=[;;?
;
if the element
incremental
rotation
matrix is defined as
+X]
(44)
where cu,, aY and LY,are the element incremental rigid body rotations results in a change of the element vector r, according to Eq. (24):
about x, y and z axes, respectively.
Firstly, the t-element
conditions
is rotated about z axis by setting the following
This
in Eq. (45)
a; = 0 a,, = 0 $2 = 0 and this results in the calculation a:=--
of the element
rotation
about z axis as
Yl2 X12
Similarly
it is possible
to calculate
the rotations
about y and x axes as (47)
(48) Since the final nodal parameters of the element include the rotations along the sides, it is therefore necessary to project the vector of rigid rotation onto the sides of the triangle and from geometrical considerations this gives the incremental rigid body rotation along side k as (49) where the indexes
i, j, k = 1,2,3
correspond
to the three sides of the triangle
and
J.L. Meek, Y. Wang I Comput. Methods Appl. Mech. Engrg. 162 (1998) 301-315
Then the incremental rotational deformation at Loof node i is obtained by subtracting this part of the rigid body rotation ‘Ye,from the calculated rotational displacement at Loof node i from equilibrium equations, PC,, Pi=&,-~kr
i=l,2,3
2.4. Step by step integration
,....,
6
k=l,2,3.
(50)
method
In the current study, the Newmark integration scheme is used to analyse the nonlinear The flow chart for the above integration procedure is presented in Fig. 5.
dynamic behaviour
givenattimet=O: displacements r0 velocities +O I accelerations: f. = M-l F(O) - Cko - R(rO)]
initialisation of accelerations: fk = EL-1 I displacements and velocities (predictor): rk = r&l + At&l + Ata(4- CY) iike1 + At%& ik = ik-1 + At (I- 6) fk-1 + Lit&k
-
I
t
dynamic out-of-balance forces: R&k) = F - Mfk - Cir, - R(Q)
diiplacement correction: [-g + &M + AC + K&h
I
=
%(rk) 1
I
t
adjusted displacements, velocities and accelerations (corrector): rk = rk + Ark +k = fk =
ik
+
fk
+
’ Ark ;;pi -Ark
Fig. 5. Flow chart of the Newmark
procedures
applied to nonlinear
dynamic
analyses.
[7].
J.L. Meek, Y. Wang I Compui. Methods Appl. Mech. Engrg. 162 (1998) 301-315
311
3. Examples The numerical examples that are presented in this section demonstrate the capability of the formulation and solution methods used in this paper. Also, the difference between the total Lagrangian and the updated Lagrangian formulation is shown to be as follows. Using the updated Lagrangian formulation can take larger load step than that of the total Lagrangian formulation, especially for slender structures. The structures were discretized by the finite element method with the DKL+LST element developed by Meek and Tan [4], using a basic rectangular element composed of four sub-triangles.
3.1. Examples for static analysis 3.1.1. Hinged cylindrical panel under a concentrated load The first example considers the large displacement behaviour of a cylindrical panel under a central vertical load. The cylindrical panel depicted in Fig. 6 has been analysed by various researchers for different loading and boundary conditions, well into the advanced nonlinear range. As shown in Fig. 6, the radius is R = 2540.0 mm, thickness h = 6.35 mm, length L = 2 X 2540.0 mm. The material properties are E = 3.10275 kN/mm*, v = 0.3. Taking into account the symmetry of geometry, boundary conditions and loading, one quarter of the panel was modelled by finite elements. The mesh used of 65 nodes and 16 elements. Fig. 7 shows the load-deflection behaviour of the panel at load point for both the total and updated Lagrangian formulations. Fig. 8 shows the the shapes of the cylindrical panel at different loads. Both total and updated Lagrangian formulation used a load increment of 0.01 kN.
‘” I o.I1
R
Fig. 6. Cylindrical
(a) Load=O.kN
shell.
Fig. 7. Load-displacement
(b) Load=-0.5687kN
(c) Load=O.O29kN
Fig. 8. The shapes of panel at different
loads.
curve of the panel shell at central point.
(d) Load=O.3616kN
J.L. Meek, Y. Wang I Comput. Methods Appl. Mech. Engrg. 162 (1998) 301-315
312
h
M
Js: b
L
U lW
Fig. 9. A cantilever
beam subjected
Fig. 10. Load-displacement
to an end moment.
curve for the cantilever
beam.
3.1.2. Cantilever beam under an end moment
This classical elastic problem has been widely used by many investigators as a benchmark test for large deformation analysis. As shown in Fig. 9, the geometrical properties of the beam are length L = 10.0 m, width b = 1.0 m and thickness h = 0.1 m. The Young’s modulus of the beam is E = 1.2 X lo5 N/m2 with Poisson’s ratio Y = 0.0. The load increments AP = 0.01 kN. Analytical solutions show that under end moment M, the beam rolls up into a circular arc of radius p: p=~=
L
1.59155m
(a) Lload=o.kN
(b) Load=3.47kN Fig. Il.
(c) Load=6.666kN
The shapes of the cantilever
beam at different
loads.
Fig. 12. Load-displacement
curve for T.L. with 20 elements.
Fig. 13. Load-displacement
curve for T.L. with 40 elements.
(d) Load=6.538?‘kN
J.L. Meek, Y. Wang I Comput. Methods Appl. Mech. Engrg. 162 (1998) 301-315
Fig. 14. Step load.
313
Fig. 15. Time history of the panel shell at central point.
In the finite element discretization, a 103 node, 20 element mesh was used. For the updated Lagrangian formulation, the load increments were AP = 0.01 kN. Figs. 10 and 11 include results obtained from the present investigation. These indicate that the beam rolls up to a complete circle when M = 6.5 kN. In this configuration, a maximum lateral displacement of 7.247 m is observed at the mid point of the beam. For total Lagrangian formulation, when the element number is taken 20 and 40 the results are not convergent at load level of 3.78 kN and 5.62 kN, respectively. They are shown in Figs. 12 and 13. 3.2. Examples of dynamic analysis 3.2.1. Hinged cylindrical panel under a concentrated load This cylindrical panel has been studied in 3.1.1 and in here E = 0.310275 kN/mm’, mass density p = 3210.05 kg/m3. The apex load is applied as a step load shown in Fig. 14. P, = 400 N and t, = 0.01 s. The time step size At = 0.001 s. Fig. 15 shows the dynamic response curve of the panel at the load point and Fig. 16 shows the shapes of the cylindrical panel at different times. It can be seen that, at time = 0.048 s, the displacement first reaches the maximum value which is -37.4 mm, the shape of the panel becomes a concave shell. The displacements obviously include finite rigid rotations.
(a) Time=O.sec.
(b) Tie=O.O4sec.
(c) Time=O.OSsec.
(d) Time=O.l2sec.
(e) Time=O.%ec.
(f) Time=O.2aec.
Fig. 16. The shapes of panel at different
times.
J.L. Meek, Y. Wang I Comput. Methods Appl. Mech. Engrg. 162 (1998) 301-315
314
Fig. 17. The rectangular
plate.
Fig. 18. Time history
of centre displacement
of the rectangular
plate.
(b) Time=O.O8sec.
(a) Time=O.aee.
(c) Tiie=O.l6sec.
Fig. 19. The shapes of plate at different
3.2.2.
Simply supported
(d) Tiie=0.32sec.
times.
plate
This example considers the geometric nonlinear dynamic analysis of a rectangular plate with simply supported edges. The plate is shown in Fig. 17. Its dimensions are a = 1016 mm, b = 1524 mm and thickness h = 25.4 mm. The material properties are: Young’s modulus E = 2.0955 X lo8 N/m*, Poisson’s ratio v = 0.25, and mass density per unit volume p = 3210.05 kg/m3. Owing to symmetry of geometry, boundary condition and loading, one quarter of the plate was discretized. For direct integration the time step size was 0.004 s, twice that suggested in [8]. The step function shown in Fig. 14 was considered, with PO = 44.54 N and t, = 0.006 s. A linear and nonlinear dynamic analysis was performed. The results are plotted in Figs. 18 and 19. Comparing the results with those of Wilson et al. [8], very good agreement is observed. The nonlinear results had smaller amplitude and shorter period compared with the linear solution, which indicates stiffness hardening due to the inplane forces and deformed shape.
4. Conclusions
The objective in this study was to derive and evaluate finite element formulations for general nonlinear static and dynamic analysis which have been implemented in the search for the most effective procedure. A simple, practical formulation and program to remove the restriction of small rotations for geometrically nonlinear dynamic analysis of shells using flat triangular shell elements are presented. The conceptual difference between the formulations is the reference configuration that is used for the linearization of the incremental equations of motion. In the total Lagrangian formulation the initial configuration is used as reference, whereas in the updated Lagrangian formulations, the reference configuration corresponds to the last calculated configuration. An advantage of the total Lagrangian formulation is that the derivatives of the interpolation functions are with respect to the initial configuration, and therefore need be formed only once. On the other hand, for updated Lagrangian formulation, the convergence rates are better than that of total Lagrangian formulation, especially for slender structures. It is very clearly shown in Example 3.1.2.
J.L. Meek, Y. Wang I Compuf. Methods Appl. Mech. Engrg.
162 (1998) 301-31.5
315
References [l] K.J. Bathe, E. Rramm and E.L. Wilson, Finite element formulations for large deformation dynamic analysis, Int. J. Numer. Methods Engrg. 9 (1975) 353-386. [2] J.H. Argyris and P.C. Dunne, On the application of the natural node technique to small strain large displacement problems, World Congress on Finite Element Methods in Structural Mechanics, Boumemouth, 1975. [3] G. Wempner, Finite elements finite rotations and small strains of flexible shells, Int. J. Solids Struct. 5 (1969) 117-153. [4] J.L. Meek and H.S. Tan, A faceted shell element with Loof nodes, hit. J. Numer. Methods Engrg. 23 (1986) 49-67. [5] J.L. Meek and H.S. Tan, A discrete Kirchgoff Plate bending element with Loof nodes, Comput. Stmct. 21 (1985) 1197-1212. [6] H. Goldstein, Classical Mechanics (Addison-Wesley, 1980). [7] J. Argyris and HP Mlejnek. Dynamics of Structures (North-Holland, Amsterdam, New York, Oxford, Tokyo, 1991). [8] E.L. Wilson, K.J. Bathe and R.H. Iding, A structural analysis program for static and dynamic response of nonlinear systems, SESM 74-3, Structural Engineering Laboratory, University of California, Berkeley, CA, 1974.