- Email: [email protected]
Nonlinear analysis of shell structures by degenerated isoparametric shell element
#4$.W9/89
s3.00 + WI0
Q 1989 Pcrgamon Press pit
Km-MO
HSAO
(Received
and
YEH-REN CHEN
15 January 1988)
Abstract-Two rotation strategies termed the finite rotation method atrd the mixed rotation method are proposed to describe the rotation of the shell normal and four rotation strategies in the literature are reviewed. The rotadan variables of the f&&erotation method are chosen to be the incremental rotations with respect to the ICY and X, axes of 8 moving coordinate system rigid& tied to the she& Both the rotation increments be?ween two sacce&ve increments and the rotation corrections be&veer3 tw62 successive iterations are used as the incrementfll rotations. The previous convergent stress is employed to update the geometric stiffness matrix and its pt?rformance is compared with that of the standard geometric stiffness matrix update method. The six types of rotation variables are embedded into a shell element whose formulation is based on that of the degenerated isoparametric shell element proposed by Ramm. In order to test and compare the performance of the various rotation strategies, ~urnero~s numerical examples are studied. The n&d rttgorithm ufed here is an in~~~a~~~~t~ve method based on the ~~to~R~~sa~ method combined with the constant arc length method.
INTRODUCTlON
Since
the
~~~~~r~t~
sfieff
eI=ent
was
fhi
intro-
duced for Enear analysis hy Ahmad et af. fl], the popularity of the element has grown enarmousIy. The element and numeraus versions of it are now applied not anly in the linear analysis of shell structures, but also in the geometric and material nonlinear analysis of $h& strXz@ures[2-l ‘13% The displacement components, referred to the global coordinate system, of any generic point of the degenerated shell elements can be separated into the displacement of the she11 middle surface and the relative displacement related to the rcrlations of the shell normal. The ~presentat~on of the relative dispja~ent is based on the kineznatic a~~rnpt~on that the shell normal of the middle surface remains straight and inextensional during deformation. Due to this kinematic assumption only two components of the norma vector are independent. TINS, only two i~de~ndent variables, referrred to as rotation variables, are re@ed to uterine the ~~enta~on ofthe normal vector during deformation. ‘it is a weit knowa fact that finite rotation components are not vectorial and that the sum of finite rotations depends on the order of rotations. fn geometric no&near analysis this has led to the adoption of many different types of ro~~on variables and different rotation strategies of the shed normal to accommodate large rotation capability during the large deformation process, Ramm [2] chose tha rotation of the shell normal around the global X, axis and the angle of the shell normal with respect to the global X, axis as the rotation variables. Oliver and Onate 191detined a set
of fixed local axe$ xi (i = t, 2, 3)+ associated with a set of orthonormai vectors e, fi = 1,2, 3)>where e3 is normal to the u~defo~~ shelt surface, e, is perpen&cular to e, and the glohal X, axis, and e, is the cross product of e3 and e,. The rotation variables a and fl are taken as the rotations about the positive local xI axis and the negative local x2 axis, respectively. The rotation of the shell normal corresponding to the rotation variables is defhzed as the finite rotation ff2] resulting from the appEcation of the rotation vector 112, 131B = ae, - fie2 to the initial shell normal W, = e3. Parisch [3] and Surana [S] also define a set of fixed orthogonal local coordinate system and rotation variables in a similar way as that given by Oliver and Onate [9], but they chose different rotation axes and ~te~~t~ the rotations of the shed normal corresponding to the two rotation variables in a different way. Tho orientation of the shell normal is determined as follows: the initial shell normal OV,= e:, and the rotation vectar ae, are rotated to their new positions V; and e;, respectively by the appi~~t~on of the rotation vector &zr; then the vector V; is transparted to V, by the application of the rotation vector are{. Surana [8] aSsa tried different sequences of rotations to determine the orientation of the shell normal and compared their numerical efficiency. Bathe and Bolourchi fq, Hughes and Liu [6], and Milford and ~h~~b~~h 1%i] adopted the ~n~r~en~l rotations with respect to the x1 and x, axes of a moving coordinate system as the rotation variables. Tha x,-direction of”the moving coordinate system is chosen to coincide with the shell normal during deformation, the x1 axis is perpendicular to the x3 axis and the global X2 axis, the x2 axis defines a Cartesian
428
Kuo-MO HSIAOand YEH-RENCHEN
right-handed system. In order to avoid the use of the trigonometric representation for the rotation, the motion of the shell normal corresponding to the rotation variables is determined by the ‘vectorial rotation increments’ used in the linear case. The newly obtained normal vector is then projected radially to maintain the constant length assumption of the shell normal. It should be mentioned that the above procedure does not define rotations in the strict sense, but Cartesian displa~ments. In IS], it was mentioned that the rotation of the shell normal may be evaluated accurately by dividing the incremental rotation into subincrements. However, a large number of subincrements are required for large increment rotations. Instead of using the subincrement method, in this paper a strategy termed the finite rotation method is proposed to describe the motion of the shell normal. In this method the rotation variables tl and 0 are chosen to be the incremental rotations with respect to the x, and x, axes of a moving coordinate system which is rigidly tied to the middle surface of the shell. The initial xj direction of the moving coordinate system is taken as the normal of the undeformed shell surface, the initial x, axis is perpendicular to the x3 axis and tangential to a proper element side, and the xz axis is established by taking the cross product of the unit vectors associated with the xj and xi axes. The rotation of the moving coordinate system corresponding to the incremental rotation is determined by the application of the rotation vector B = ae, + /?e_,to the triad of the unit vectors associated with the moving coordinate system. Here, both the rotation increment between two successive increments and the rotation correction between two successive iterations are used as the incremental rotation. Since the choice of rotation axes may affect the numerical efficiency of nonlinear analysis, a mixed rotation method, which is a combination of the proposed finite rotation method and Bathe’s method [S], is also proposed to compare the ~rfo~an~e of the same rotation strategy with different rotation axes. Theoretically, it seems that all six types of rotation strategies mentioned above can handle arbitrarily large rotations and permit large rotations between successive increments. However, because they are employed in different shell elements, their convergence characteristics have not been compared in the literature. In order to test and compare the performance of these rotation strategies, they are embedded into a shell element whose formulation is based on that given by Ramm [2] using a total Lagrangian approach, and numerous numerical examples are studied. In [13] the previous convergent stress method was employed to update the geometric stiffness matrix for nonlinear analysis using a flat triangular shell element. It was reported that this method can improve the convergence characteristics for nonlinear analysis. Here this method is referred to as the PS
Fig. I. Nine-node degenerated shell element.
method and its performance is compared with that of the standard geometric stiffness update method, referred to as the CS method here. The numerical algorithm used here is an incremental-iterative method based on the Newton-Raphson method combined with a constant arc length of incremental displacement vector (7, 141. FINITE ELEMENT FORMULATION
In the following only the basic relations of the shell element are derived. A more detailed description may be obtained from [2]. Geometry and kinematics of the element
The geometry of the shell element is described in Fig. 1 by the coordinates and normal vectors at the nodes on the middle surface. The geometry and displacement representation are based on the kinematic assumption that the normal to the middle surface remains straight and inextensional during deformation. On the basis of the first assumption;the current configutation of the global Cartesian coordinate vector X = I&‘,, X,, X,) of a generic point of the shell element is expressed by X= 5
N’(r,s)X’+
f c
k=l
‘k=l
h&N&(r,s)V:,
(1)
where M is the total number of nodes of the shell, Nk is the matrix of shape functions associated with node k, Xk is the Cartesian coordinate vector of node k, r and s are the local curvilinear coordinates of the middle surface, t is the coordinate in the direction of the shell normal, hk is the thickness of the shell at node k, and Vi = fYfs, Vt3, Vi,) is the unit normal vector at node k. The displacement vector U = {iJ1, U,, U3f between the initial configuration 0 and the current configuration are specified to be consistent with eqn (I) U = 5 Nk(r, s)Uk + G 5 hkNk(r, s)(V; - “V$), (2) k=l
k-l
where Uk = {CJ:, Vi, .Y:} is the current displacement vector of node k, and “v: is the unit normal vector at node k for the undeformed configuration. It should be noted that because of the inextensional
Nanbar
anaiysis of shell structures
assumptions of the shell normal only two components of the nodal normal vector are independent. Thus, only two independent variables, referred to as nodal rotation variables, are required to determine the orientation of the nodal normal vector during deformation. So, i$ (I’= 1,2, 3) can he expressed as a function of two rotation variables c1and /I, and thus each node of the shell element has only five degrees of freedom: three translations in the direction of global axes and two rotation variabies ta be speci%d later. When the rotation variables and the rotation axes are chosen, the variations of displacements, which are required for the derivation af the element internal nodal force vector and element stiffness matrix, can be expressed as a function of the variations of nodal variables. From eqn (2), the variation of the disptacement vector can be written as SU = g N’(x, s)6Uk
kzl
where vectors f and g are obtained by differentiation of Vi with respect to the nodal rotation variables a“ and flk, respectively, and are given by f = 8V$/c3aR g =Ea~~]~~~.
14)
Assume that the incremental iterative method is used for the solution of nonlinear equilibrium equa. tions and that configuration I is known. Here ~~gu~t~on I may denote the equ~~ib~urn ~~g~t~oo of the previous increment or the iteration. configuration of the previous Let A& = {AU:, AU& AU:} and A+k = {A.ak, A/?“} be the given incremental nodal displacements and incremental nodal rotations, respectively, at nodes k (k = 1-M). In this paper when configuration I denotes the equilibrium con~~ra~o~ of the previous increment, AU’ and A#k denote the incremental nodal displacement and rotation increments between two successive increments. When the configuration I denotes the configuration of the previous iteration, Auk and Adik denote the nodal displacement and rotation correction between two successive iterations. Let Ixk, %I*, ‘v! and ‘+” = {‘a”, ‘/I”] denote the nodal coordinate vectors, nodal displacement vectors, unit shell normal vectors and nodal rotation vectors, respectively, at nodes k (k = 1-M) of an element at configuration I; the current values of nodal coordinates, nodal dispiacements and nodal rotation variables can be obtained by adding the corresponding incrments, so Xk=‘Xk+AUk Uk=%Jk+AUk Qlk= ‘rlik+ A#.
(5)
a!? x , ‘1
Dyformed normal
Fig. 2. Description of rotation for Ramm’s method.
At this point, an interesting question arises. Given the incremental nodat rotations, how are the current orientations of the shell normal Vi dete~in~~ It is a weil known fact that finite rotation components are aat vectorial and that the sum of finite rotations depends on the order of rotations. In geometric nonlinear analysis this has led to the adoption of many different types of rotation variables and different rotation strategies of the sheil normal to accommodate large rotation capability during the large deformation process. In the following, four different types of nodal rotation variables and rotation strategies proposed in the literature are reviewed and two sets of nodal rotation variables and rotation strategies are introduced. Note Chat ek in eqn (5) is only used for the A, 3 and C types of rotation strategies to determine the orientation of the shell normal; q: and A4k are used for the D, E and F types of rotation strategies. For convenience, the subscript k used to indicate node k is dropped throughout this paper to avoid confusion. (A) Ramm’s method. fn [Z] 01(Fig. 2), the angle between the global X,-X, plane and the plane formed by the shell normal and the global X, axis, and /3, the angle between the shell normal and the giobal X, axis are selected as the nodal rotational variables. The unit normal vector at node k in eqns (1) and (2) may now be written as V, = (cos b, sin /? cos a, sin @sin u }.
(6)
From eqns (4) and (fr), vectors f and g may be expressed as f=(O,
-siaasinjY,cosasinfij
g={-sinfl,cosacosB,sinacas~).
(7)
(B) Oliver’s m&ad. Oliver and Onate [PI defined a set of fixed local axes x, (i = f , 2, 3), associated with B set of orthonormal vectors ei (i = I, 2, 3), where e3 is normal to the undeformed shell surface, and e, is perpendicular to e3 and the global X2 axis. Note that in this study the vector eI is replaced by the unit vector which is perpendicdar to e5 and tangential to a proper element side. Vector e, is the cross product
430
KIJO-MO HSIAO and
t 0
normal and tangential to a proper element side. The rotation variables LY and /I are chosen as the rotations about the positive local xi and x, axes, respectively. The orientation of the shell normal is determined as follows: the initial shell normal qj = e3 and the rotation vector e, are rotated to their new positions Vi and ae;, respectively, as shown in Fig. 4, by the application of the rotation vector Be,; then the vector V; is transported to its final position V, by the application of the rotation vector ae;. The unit normal vector at node k may now be expressed as
x3
$ v3
Ov3
- P2 /
e
/
-
/
YEH-REN CHEN
x2
gel
V, = cos CIsin Be, - sin ae2 + cos a cos /3e3.
(10)
Xi
From eqns (3) and (lo), vectors f and g may be expressed as
Fig. 3. Description of rotation for Oliver’s method.
f = -sin CIsin De, - cos aeZ - sin a cos /3e3 of e, and e, . The rotation variables a and B are taken as the rotations about the positive local x, axis and the negative local x2 axis, respectively. The rotation of the shell normal corresponding to the rotation variables is defined as the finite rotation [12, 131 resulting from the application of the rotation vector 6 = ae, - Be2 to the initial shell normal q, = e3 as shown in Fig. 3. The unit normal vector at node k may now be written as V,=cosBe,+(l
-cos@(e;e,)e,
g = cos a cos Be, - cos a sin /?e3.
(11)
(D) &the’s method. Bathe and Bolourchi [S], Hughes and Liu [6] and Milford and Schnobrich [1 1] adopted the incremental rotations a and j with respect to the x, and x, axes of a moving coordinate system as the rotation variables. Let vectors e, (i = 1, 2, 3) be unit vectors associated with the moving coordinate axes xi (i = 1, 2, 3). Vector e, is chosen to coincide with the shell normal V, during deformation; vectors e, and e, are defined by
+ sin 8 (e, x e3) B cos Be, - i cos 8e, + cos 8e3, 8
=--
e, = I, x e3/ III, x e3 II (8)
and e2=e3 xe,,
where 6 = (a2 + fi2)“2 and e, = 0/t?. By substituting eqn (3) into (8), vectors f and g may be given as f= - $(cosO
g=-
e,-
-$(cost?-y)e,-
a
8sm0e,
e,
$costI+$y (
where i, is the unit vector in the global X2 direction. Let ‘V3= ‘es be the normal vector of configuration I at element node k, and a and B be the given incremental rotations. In order to avoid the use of the trigonometric representation for the rotation, the current shell normal V, is determined by the gener-
-y)e,
p2 sine a2 ~cosB+8’o
-
(12)
)
isin*e,.
(9)
(C) Parisch’s method. Parisch [3] and Surana [8] define a set of fixed orthogonal local coordinate systems and rotation variables in a similar way to that given by Oliver and Onate [9]. The local xi axis is chosen to be perpendicular to the initial shell
Xl
Fig. 4. Description of rotation for Parisch’s method.
Nonlinear analysis of shell
431
structW%i
incremental rotations at node k. A rotation vector is defined by 6 = a’e, + fl’e, = Bn,
(13
where 0 = 11 B 11is the angle of rotation and n e e/0 is the unit vector along the axis of the rotation, The current base vectors are determined by the application of the rotation vector 8 to base vectors ‘ei and expressed as (see Fig. 6) [l 11
e,scos@‘e,+ (1 -cosB)(r~-~e&r + sin @I x ‘e,) Fig. 5. Description of rotation for Bathe% method.
alization of the ‘vectorial rotation increments’ used in the linear case and given as follows (see Fig. 5): ‘V, -
” =e3= I/‘v3 -
a’e,4” B’q are2 -t_jPe,/I*
f= -e,
g=ei.
i = 1, 2, 3.
As mentioned earlier the current shell normal vector Vj - e,. From the definition of the rotations of the normal vector about the rotation axes corresponding to the variation of the rotation variables, vectors f and g in eqn (3) may be given by
(13)
It should be mentioned that the above procedure does not de&e rotations in the strict sense, but Cartesian displacements. After the vector VS is determined, the current rotation axes e, and eZ can be determined using eqn (12). From the definition of the rotations of the normal vector about the rotation axes corresponding to the variation of the rotation variables, vectors f and g in eqn (3) may be given by
06)
f= -e,
g=ei, where 9 and ei are given in eqn (16). (-F) Mixed rotationmethod. The choice of rotation axes about which the rotation variables are defined may affect the convergence rate of equilibrium iterations. In order to investigate this effect, a mixed rotation method is proposed here to determine base vectors. The shell normal vector V, = e3 is determined by eqn (16), but e, and e, are determined using eqn (12). Vectors f and g in eqn (3) are the same as those given in eqn (14).
(14)
(E) Finite rotation ~e~~u~. I&e a finite rotation method is proposed to describe the motion of the shell normal. In this method the rotation variabies are chosen to be the incremental rotations, a and jl with respect to the x1 and x2 axes of a moving coordinate system which is rigidly tied to the middle surface of the shell. The initial xg direction of the moving coordinate system is taken as the normal of the undeformed shell surface, the initial xi axis is perpendicular to the x, axis and tangential to a proper element side, and the xi axis is established by taking the cross product of the unit vectors associated with the X, and xi axes. Let vectors e, fi = 1,2,3) be base vectors, unit vectors associated with the moving coordinate axes xj (i = 1, 2, 3). Due to the definition of the moving coordinate system, vector e3 coincides with the shell normal Vj during deformation. The rotation of the moving coordinate system corresponding to the incremental rotations is determined by a finite rotation scheme similar to that proposed by Hsiao [19. Let ‘ei denote the base vectors at node k of configuration I, and 01 and /I be the given
The nonlinear expressed by
equilibrium
$ =F-IP=O,
equations
may
be
(W
where # is the unbalanced force between the internal nodal force vector F and the external nodal force dP;
Fig, 6, Deacriptian of rotation for finite rotation method.
1 is a loading parameter and P is a no~aiiz~ loading vector. The internal nodal force vector is obtained by summing up the element nodal force VeCtOorS. In this paper, a weighted Euchdean norm of the unbalanced Force [ f 5j is employed as the error measure of the equiiibr~um state during the e~uil~b~um iterations, and the convergence criterion is given by
where N is the number of degrees of fr~dom for the d~s~~et~~d structure and ptOSis a prescribed value of error tolerance. Unless it is stated otherwise, the error tolerance is set to 10 -4 in this paper,
An incremental-ite~t~ve method based on the ~ewto~~aphs~n method is adopted here. In order to deai with the hmit points and snap through, the arc length of the incrementat d~sp~a~m~nt vector is kept constant during the ~u~~~br~~rn iteration using Crisfield’s method f7, 141. The PS and CS methods are employed to usage the geometric st~~ness matrix and their ~rfo~an~es are compared. if the ~ui~ibrium con~gurat~on of the Ith increment is assumed to be known, the system infect stiffness matrix K+ then can be calculated at this configuration and an initial displacement increment Aq for the next increment may be obtained by using the Euter predictor as Aq = Ariq,,
QQ)
where Al is the initial in~rementa1 loading param&r and qr = K; ‘P is the tangential displa~ment of unit goading P. For all ~ncreme~fs other than the first, A;t is obtained in much the same way as that mentioned in f143, and is given by A.;1= f Aa(qq;q,)‘;“,
(zx)
while the sign is chosen following an ap~roa&h due to Bergan and Soreide [Id] in which the sign fatlows that of the previous ~nc~ment unbss the d~tE~j~~~~ of the tangent sti&ess matrix has changed the s&n, in which case a sign reversal is applied. rt\ais the incremental arc length used for next inurement, and is dete~ined by
number of itara~ons~ the cutout para~te~ C, and C, are chosen tu be 03 and 1.5, resp&iveIy, to prevent the yielding of an incremental displacement which is too large or too smalt. Using the djs~~~rn~nt &rement obtained in eqn (20) and ane oF the six rotation strategies of the shell normat deemed in the pre~nt paper, the current eonfi~rat~o~ of the shell structure can be determined. Then, the internal noda force vector F in eqn (tgj donated with the current e~~~gu~at~on can be caleutated, The foadirtg parameter d corresponding to the current ~on~gura~on is given by ,I-- li+ 81, where t’ is the convergent loading parameter at the fth ~n~~ment and Atl the loading parameter increment. Then the unbalanced force I& can be obtained from eqn (IQ If the convergence criterion feqn (I!?)] is not satisfied, a ~sp~~rn~~t correction r and B loading parameter correction SI 17, 141 are added to the previous Aq and A\A respectively to obtajn a new in~r~e~~ displacement and ~~~rnenta~ foading parameter for the next ~te~t~o~. The values of r and lil may be deterred by
and Aa 2 = (Aq + @If&q+ r)y
(zs)
where K, may be the tangent stifiness matrix at same known configuration. Two different st~ate~es are used to determine the Cflrrent Confi~ffrat~on. fn the first Stt’stegy, referred to as the dj~~acem~nt increment (LX) method, the current d~s~~~~m~nt i~~r~~~~~69 and the previous eq~~~b~~rn cor&guration are used Fcr eonstro.ct tb&! current configuration, fn the second strategy, referred to as the displa~m~nt correction (IX) method, the current con~guratio~ is determined using the displacement correction r and the configuration of the previous iteration. When these two methods are in~~~orat~ with any one of the first three rotation strategies of the shell normal mentioned in the present paper, ide~tjcal configuration wilt be obtained, because the rotation variables of these three methods are de&ed in fixed coordinate systems. But when these two methods are ~ornbjn~ with rotation strategy with rotation variabks de&d in axes of rotating coordinate system, diiYerent ~o~guratio~s may be obtained. This may a&et the convergence ~haracte~stics of ~u~~jb~urn iteration and wilt be investigated in the next section. Then, the internal nodal force corresponding to the new ~o~~~~tio~ obtained can be calcufated to initiate the next iteration. This procedure is repeated until the convergerme criterion is satisfied.
433
Nonlinear analysis af sheli structures
In the following, for convenience, XYYZZ is used to denote that X (A-F) type of rotation variables combined with W @I or DC) and ZZ (I’S or CS) methods is employed. E = 1.2X108kN/m2 II=0 L = 10.0 m b = 1.0 m h = 0.1 m l
r
Analyticul Present
Example 1. Cantilever beam with an end moment
solution
{ ) : NO. of ilaration
Displacement
(m)
Fig* 7. Cantilever beam with an end momeat.
geometric stiffness matrix for iteration strategies based on the Newton-Raphson method, In this paper, this method is referred to as the current stress (CS) method. As mentioned in [13], during the first few equilibrium iterations, the values of the stresses obtained from the current deformation may be several orders larger than their convergent values for certain problems with very large rotation increments. This may cause difficulty in convergence or even divtkgence for a large rotation increment using the CS method to update geometric stiffness matrix. in order to overcome this difficulty, a previous convergent stress (PS) method was proposed in [13] to update the geometric stiffness matrix. In this method the convergent stress of the previous increment is adopted to &c&ate the geometric sti&ess matrix. The convergence characteristics of the CS and PS methods will be investigated and compared using numerical examples in the next section.
This example is a cantilever beam subjected to a concentrated moment at the free end as shown in Fig. 7. The beam is idealized using five elements. For this two dimensional problem, the orientation of the shell normal determined by all types of the rotation variables and rotation strategies are identical except that determined by the DDI or the DDC method. It is found that if the CS method is used to update the geometric stiffness matrix, the convergence of the equilibrium iterations cannot be achieved even for small in~~ment steps. Thus, only the results obtained using the FS method are given for this example The results shown in Fig. 7 are obtained by using the APS method. It is seen that five increments are used. The number of iterations for each increment is also given in parenthesis beside each point on the graph. Nearly identicaI results are obtained for alI methods except that obtained by the DDIPS method. For the DDIPS method, 14 increments are used to reach M = 6.2 kNm and the average number of iterations per increment is 1I _
Exczmpfe2. ~~arnped~~~~edcirculaf a& The flamed-hinted deep arch depicted in Fig. 8 is subjected to a paint load at the crown. A 16 x 1 mesh is used to represent one half of the arch. The
E R @ t b
= = = = =
6Xi06
tb/‘in?
100 in. 215Q 1.0 in. 2.0 in,
j/j” section at A-A
041
9
u
8
In order to compare the performance of the six types of rotation variables and rotation strategies combined with the DI or DC method and PS or CS method, many numerical examples of plate and shell structures are studied here. The nine-node Lagrange element is employed for nonlinear analysis. The element matrices are calculated using the 2 x 2 x 2 Gauss integration. Crisfield’s constant arc length method [7, 141 combined with the Newton-Raphson method is used for equilibrium iteration exclusively. Unless otherwise stated, the error tolerance pIol is set to 10-4,
7 * (26)
cF
3
Present
2
( ;
: No. of iteration
1
0 Fig. 8.
0.2
0.4
0.6
0.8 U/R,V/R
1.0
1.2
1.4
Clamp&-hinged circular arch.
Kuo-MO HSIAOand
434 Table
1. Comparison of solution methods clamncd-hinaed circular arch Max.
ACS DDCPS DDIPS FDIPS
Error tolerance
No. of increments
Total No. of iterations
1.198 1.201 1.184 1.203
3 x IO-4 10-d lO-4 IO-4
30 6 25 6
223 55 283 56
in Fig. 8 are obtained by the FDIPS method using six increments. For this two dimensional case, except for the D type of rotation variables all types of rotation variables are equivalent. When the CS methods are used, the convergence of equilibrium iterations becomes very difficult and thus a looser error tolerance ptO,= 3 x 10m4is used. A comparison of performance among different methods is given in Table 1. The PS methods give significantly better convergence characteristics than the CS methods for this problem. As in the previous example, the DDCPS method shows significantly better performance than the DDIPS method. results
Table 2. Comparison of CS and CS matrix update method for column under a vertical load (B = 1)
for
U ii
Solution methods
YEH-REN CHEN
cr
The column shown in Fig. 9 is subjected to a vertical end load which is higher than the elastic buckling load. Two cases, B = 1 and B = 10 are considered. For the case when B = 1. five elements
PI L = 100 Lo= 20 h=l B = 1, 10 E = 12 Thickness = 1 p, = o.247x1o-3
B =
B=l ----
1.9580
.
3.3970 3.7554
0.0173 0.0864 0.2550 0.3115 0.5340 0.853 1 0.9067 1.2609 1.3058
ACS I
APS
; 8 6 8 5 10 9 5
are used and for the case when B = 10, six elements are used. To avoid numerical difficulty, the out of plane degrees of freedom are restrained for the case B = 10. Here the shell elements are used in the ‘plane stress mode’, therefore no rotation of the shell normal occurs during deformation. Hence, the choice of the type of the rotation variables does not affect the results. Only the performance of the CS and PS stiffness matrix update methods are compared for this example. Details of their convergence characteristics are given in Tables 2 and 3. As with the previous examples, the PS method gives much better convergence characteristics than the CS method. The results obtained using the APS method are shown in Fig. 9 together with the results given in [8] in which 20 load steps are used. Example 4. Cantilever plate with end load
t
The example considered here is a cantilever plate subjected to a concentrated end force as shown in Fig. 10. A four by three mesh is used for discretization. The plate is first analyzed using the APS and ACS methods to compare the convergence characteristics of the PS and CS matrix update methods. The results are given in Table 4. Clearly the PS method gives significantly better convergence characteristics than the CS method as with the previous
iTI s
r
10
-
Surana[B] Present
P
0.9317 1.0294 1.1417 1.1834 1.3829 1.8445
shown
Example 3. Column under a vertical end load
No. of iterations
P
Surana [ 81 Present
.
Table 3. Comparison of CS and CS matrix update method for column under a vertical load (B = 10) P
0
I
0.2
I
!I!
0.4
0.6
I
l
0.8
II
1.0
I!
1.2
U/L, V/L
Fig. 9. Column under a vertical load.
1
I
1.4
II
PC,
r,
1.4676 1.5628 1.7071 1.7245 2.0375 2.3708 2.6335 3.8676 3.9304
0.0177 0.0790 0.1941 0.2073 0.4061 0.5575 0.6479 0.8983 0.9066
No. of iterations ACS
APS
7 6 6 7 8 5 8 7 5
Nonlinear anafysis of sheif stndctures Table 5. Iterative performance for cantilever plate with end load (the table gives the number of iterations for each
30 m 40 m 0.4 m 1.2x1 d
r V.XIO~
40
u
7’ W . u = 0.3
iilCrement)
kN/mz G
APS BPS CPS DDIPS DDCPS EDIPS EDCPS FDIPS FDCPS
‘;24 c a x 4 16
4 4 4 8 4 7 4 7 4
4 4 4 6 4 4 % 4
4 4 4 6 4 5 4 5 4
Example 5. Simply supported plate subjected to lateral pressure
Deflection
(m>
Fig. 10. cantilever plate with end
load.
examples. Three i~~ements are used for the APS method and the load versus end displacements curves are shown in Fig. 10, together with those reported by Hsiao [ 131. The plate is later reanalyzed using different types of rotation variables in conjunction with the PS stiffness matrix update method, In order to compare their performance, the same fixed incrementaI arc length is used for each increment. The details are given in Table 5. Clearly the three DC methods give better convergence than the three cor~sponding DI methods as with the previous examples. The same ~rfo~an~ is obtained for the XPS (X = A, B, C) methods and the XDCPS (X = D, E, F) methods, From the results obtained using the EYYPS and FWPS methods (Table 5), it seems that the choice of the rotation axes only has a small effect on the convergence characteristics.
The pIate is subjected to a uniform pressure loading. Figure 11 shows the geometry and physical data of the plate. A quarter of the plate is analyzed using a two by two element mesh. It is found that the performance of all methods is similar for this example. However, the buoyance of the DC method is slightly better than that of the DI method, and the performance of the I’S method is mar~na~iy worse than that of the CS method. Note that the order of performance between the PS and the CS methods for this example is contrary to that for the previous examples in which the deflections are up
q = uniform pressure
/---a---l’ 3CD *s : ::
Levy 1171 Present
( )
l
:
No. of iteration l
(31
I
f
Table 4, Comparison of CS and Ps matrix update methods for cantilever nlate with end load No. of iterations
6.492 12.650 16.307 19.559 27.662 36.336 37.357
w&N
A(Js
7.897 13.903 16.719 18.810 22.730 25.548 25.818
6 5
APS 4 4
: 4 5
Centrat Deflection
Ratio W,/t
Fig. Il. Simply-supPorted Plate with a uniform lead.
436
E R L h Y
60
= = = = =
68.95 N/mm2 2540 mm 1569.8 mm 99.44 mm 0.3
r
40 s? t 30 a. 10 9w 20
10
0
Fig. 12.
-----
Surano [8] Horrigmoe and Bergan [Is]
_
Present ( 1
50
No. of iteration
100 $50 ZOO 250 Central Deflection ~c(mm~
300
Hinged spherical shell (vertical cutting).
Also shown in Fig. 12 are the solutions given in [&l&19], and in Fig. 13 are the solutions given in (81.The solutions given in [8] are obtained using 12 increments with about 50 iterations. Example 7. Hinged cylindrical shell (h = 12.7 mm) Figure 14 shows a circular cylindrical shell subjected to a concentrated central load on the convex side. The longitudinal boundaries are hinged and immovable, whereas the curved edges are complete& free, Here the symmet~~i response is studied by using a two by two mesh for one-quarter of the shell. This example is analyzed using all the six types of rotation variables. It is found that the performance of all methods is similar for this example. However, the performance of the DC method is slightly better than that of the DI method as with the previous examples, and the performance of the CS method is marginally better than that of the PS method as with examples 5 and 6. The results obtained using the EDICS method are given in Fig, 14 with the solutions reported by Surana [81 using the C type of rotation variables and Oliver and Onate [9] using the B type of rotation variables. The results of the present study are obtained by using three increments and 15 iterations. The solutions given in [8] are obtained using 12 increments and 54 iterations.
to the order of ma~itude of the major st~ctural Example 8. Hinged ~yl~nd~~~al shell (h = 6.35 mm) dimensions. In this example the shell thickness of the previous The results obtained using the ACS method are given in Fig. 11 with Levy’s classical solution [17]. example is reduced by one-half, i.e. k = 6.35 mm. The The results of the present study are obtained by using three increments and 10 iterations,
I-he we consider a spherical shell subjected to a concentrated load at the crown. All four edges are hinged and immovable. Two cases are considered: (1) the spherical shell is formed by four vertical cutting planes (Fig. 12) and (2) the spherical shell is formed by four cutting planes passing through the center of the sphere (Fig. 13). A two by two element mesh is used for one-quarter of the shell. Only the 3, C and E types of rotation variables can describe the boundary conditions of this example properly, therefore only these three types of rotation variables are used. Nearly the same performance is obtained for these three types of rotation variables, The results obtained using the EDICS method for both cases are shown in Figs 12 and 13, respectively. It is seen that only three increments with 13 iterations are used for the first case and three increments with 14 iterations are used for the second case. The solid curves shown in Figs 12 and 13 are also obtained by the present study using 17 increments to obtain complete loadingdeflection curves. As with the previous example, the CS method gives a marginally better performance than the PS method.
& R L ;
= = = 7 -
68.95 N/mm’ 2540 mm 1569.8 mm ;9;45 mm
.
I’
P Xx
:: 20 f
____
I’
_
~~
y
.
( , ,;;;,;
0
f3.
Present
;lym;,
3~*
50 Central
Fig.
Surano [ 8}
Deflection
W,(mm>
Hinged spherical shell Qsyramid cutting).
Nonlinear anaIysis of
R = 2540 mm L = 508 mm tl = 12.7 mm
R
437
shell structures ----^-
tlorrigmoe
---.
Ramm [20]
and @argon [18]
.
6 = 0.1 rodions E = &I0275 l&/mm2
V = 0.3 4
-
Oliver [ 9 ]
-----
Sunara [S]
(5) ( >
No. of itemtion
3
5
10
15
20
25
30
Deflection (mm)
Fig. 15. Hinged cylindrical shell (h = 6.35 mm).
Fig. 14. Hinged cylindrical shell (h = 12.7 mm).
remaining geometric data, physical data and bout& at-y conditions are unchanged. A quarter of the shell is analyzed using a two by two element mesh. The shell is first analyzed using the APS and ACS methods to compare the convergence characteristics of the PS and CS matrix update methods. It is found that the CS method gives s~gni~~antly better convergence ~hara~~e~stj~ than the PS method far this highly nonlinear problem. Three increments with 21 iterations are used for the ACS method and 13 increments with 73 iterations are used for the APS method. The load versus end displacements curves obtained using the ACS method are shown in Fig. 15, together with those reported in f18, ZO].The sohd fine shown in Fig. 15 is also obtained by the present study using a 3 x 3 element mesh and 27 increments to obtain complete loading-deflection curves. This example was studied by Surana [8] using the C type of rotation variables The rest&s reported in fg] (not shown in Fig. 15) are obtained by using 2 1 increments and 98 iterations. The shell is later reanalyzed using different types of rotation variables in conjunction with the CS stiffness matrix update method. In order to compare their performance. The same fixed incremental arc length is used for each increment. The details are given in Table 6. Clearly the three DC methods give better convergence than the three corresponding Df methods as with the previous examples. Nearly the same performances are obtained for the XCS (X = A, B, C) methods and the XDCCS (X = D, E, F) methods. From the results obtained using the EYYPS
(YV = DI, DC) and FYYPS methods (Table 5), it seems that the choice of the rotation axes has no effect on the convergence characteristics For this example. CONCLUSIONS
Six types of rotation variables and rotation strategies are employed to describe the rotation of the shell normal. In particular, a finite rotation method is proposed and tested. In this method the rotation variables are chosen to be the incremental rotations with respect to the x1 and x, axes of a moving coordinate system which is rigidly tied to the middle surface of the shell. Both the rotation increments between two successive increments and the rotation
Table 6. Iterative performance for hinged cyfindricai bad (the table gives the number of iterations for each imxement)
G APS BPS CPS DDIPS
7 6 6 8
7 7 7 8
7 7 7
EDIPS DDCPS EDCPS FDIPS FDCPS
7 7 8 7
:
; 7 g 7
; 7
8
438
Kuo-MO Hsuo and YEH-RENCHEN
corrections between two successive iterations are used as the incremental rotations (rotation variables) to update the orientation of the shell normal. The previous convergent stress (PS) method and the current stress (CS) method are employed to update the geometric stiffness matrix and their performances are compared. Prom the numerical experiments that were conducted the following conclusions can be drawn. (1) The six types of rotation variables give a similar performance. (2) The displacement correction (DC) method always gives a better performance than the displacement increment (DI) method. (3) The previous convergent stress (PS) method gives significantly better convergence characteristics than the current stress (CS) method for the problems with very large displacements and rotations up to the order of magnitude of the major structural dimensions, but for the problems in which the deflections are the order of several thickness, the CS method is more efficient and is generally a little more efficient than the PS method. (4) The choice of the rotation axes has little effect on the convergence characteristics. (5) Only the B, C and E types of rotation variables can describe the boundary conditions of the spherical shells given in example 6. (6) Although no example is given here, it seems that only the E type of rotation variables proposed in the present paper can handle moving boundary conditions.
REFERENCES 1. S. Ahmad, B. M. Iron and 0. C. Zienkiewicz, Analysis
of thick and thin shell structures by curved finite element. Inr. J. Numer. Meth. Engng 2, 419-451 (1970). 2. R. Ramm, A plate/shell element for large deflections and rotations. In Formulations and Computational Algorithms in Finite Element Analysis, U.S.-Germany Symp. (Edited by K. J. Bathe, J. T. Cklen and W. Wunderlich), pp. 264293. MIT Press, Cambridge, MA (1977). 3. II. Par&h, Geometrical nonlinear analysis of shells. Comput. Meth. appl. Mech. Engng 14, 159-178 (1978). 4. H. Parisch, Large displacements of shells including
5. 6.
7. 8.
9.
material nonlinearities. Comput. Merh. appl. Mech. Engng 27, 183-214 (1981). K. J. Bathe and S. Bolourchi, A geometric and material nonlinear plate and shell element. Compur. Strucf. 11, 2348 (1980). T. J. R. Hughes and W. K. Liu, Nonlinear finite element analysis of shells: Part I. Three-dimensional shells; Part II. Two-dimensional shells. Comnut. Mesh. anol. Mech. Engag 26, 331-362 and 27, 1671181 (1981)” T. Y. Chang and K. Sawamiphakdi, Large defIection and post-buckling analysis of shell structure. Compur. Merh. appl. Mech. Engng 32, 311-326 (1982). K. S. Surana, Geometrically nonlinear fo~ulation for the curved shell elements. Int. J. Numer. Merh. Engng 19, 581-615 (1983). J. Oliver and E. Onate, A total Lagrangian formulation for the geometrically nonlinear analysis of structures using finite elements. Part I. Two dimensional problems: shell and plate structures. Inr. J. Numer. Meth. Engng 20, 2253-228 1 (1984).
10. K. J. Bathe and E. N. Dvorkin, A formulation of general shell ‘elements-the use of mixed interpolation of tensorial components. Int. J. Numer. Meth. Engng 22, 697-722 (1986).
11. R. V. Milford and W. C. Schnobrich, Degenerated isoparametric finite elements using explicit integration. Inf. J. Numer. Meth. Engng 23, 133-154 (1986). 12. H. Goldstein, Ciassicai Mechanics, 2nd Bdn. Addision-Wesley, London (1980). 13. K. M. Hsiao, Nonlinear analysis of general shell structures by flat triangular shell element. Compuf. Srrucr. 25, 665675 (1987). 14. M. A. Crisfield, A fast incremental/iterative solution procedure that handles ‘snap-through’. Comput. Strucr. 13, 5562 (1981).
15. A. K. Noor and J. M. Peters, Tracing post-limit point with reduced basis technique. Comput. Meth. appl. Mech. Engng 28, 217-240 (1981). 16. P. G. Bergan and T. Soreide, Solution of large displacement and instability problems using the current stiffness parameter. In Finite Elements in Non-linear Me&attics (Edited by P. G. Bergan et al.). Tapir Press, Trondheim, Norway (1978). 17. S. Levy, Rending of rectangular plates with large deflection, NACA Tech. Note 846 (1942). 18. G. Horrigmoe and P. G. Bergan, Nonlinear analysis of free-form shells by flat finite elements. Comput. Merh. appf. Mech. Engng 16, 11-35 (1978). 19. A. Pica and R. D. Wood, Postbuckling behaviour of plates and shells using a Mindlin shallow shell formulation. Comput. Struct. 12, 759-768 (1980). 20. E. Ramm, Strategies for tracing the nonlinear response near limit points. In Nonlinear Finite EIement Analysis in Structural Mechanics (Edited by W. Wunderlich, E. Stein and K. J. Bathe), pp. 63-89. Springer, Berlin (1981).
s3.00 + WI0
Q 1989 Pcrgamon Press pit
Km-MO
HSAO
(Received
and
YEH-REN CHEN
15 January 1988)
Abstract-Two rotation strategies termed the finite rotation method atrd the mixed rotation method are proposed to describe the rotation of the shell normal and four rotation strategies in the literature are reviewed. The rotadan variables of the f&&erotation method are chosen to be the incremental rotations with respect to the ICY and X, axes of 8 moving coordinate system rigid& tied to the she& Both the rotation increments be?ween two sacce&ve increments and the rotation corrections be&veer3 tw62 successive iterations are used as the incrementfll rotations. The previous convergent stress is employed to update the geometric stiffness matrix and its pt?rformance is compared with that of the standard geometric stiffness matrix update method. The six types of rotation variables are embedded into a shell element whose formulation is based on that of the degenerated isoparametric shell element proposed by Ramm. In order to test and compare the performance of the various rotation strategies, ~urnero~s numerical examples are studied. The n&d rttgorithm ufed here is an in~~~a~~~~t~ve method based on the ~~to~R~~sa~ method combined with the constant arc length method.
INTRODUCTlON
Since
the
~~~~~r~t~
sfieff
eI=ent
was
fhi
intro-
duced for Enear analysis hy Ahmad et af. fl], the popularity of the element has grown enarmousIy. The element and numeraus versions of it are now applied not anly in the linear analysis of shell structures, but also in the geometric and material nonlinear analysis of $h& strXz@ures[2-l ‘13% The displacement components, referred to the global coordinate system, of any generic point of the degenerated shell elements can be separated into the displacement of the she11 middle surface and the relative displacement related to the rcrlations of the shell normal. The ~presentat~on of the relative dispja~ent is based on the kineznatic a~~rnpt~on that the shell normal of the middle surface remains straight and inextensional during deformation. Due to this kinematic assumption only two components of the norma vector are independent. TINS, only two i~de~ndent variables, referrred to as rotation variables, are re@ed to uterine the ~~enta~on ofthe normal vector during deformation. ‘it is a weit knowa fact that finite rotation components are not vectorial and that the sum of finite rotations depends on the order of rotations. fn geometric no&near analysis this has led to the adoption of many different types of ro~~on variables and different rotation strategies of the shed normal to accommodate large rotation capability during the large deformation process, Ramm [2] chose tha rotation of the shell normal around the global X, axis and the angle of the shell normal with respect to the global X, axis as the rotation variables. Oliver and Onate 191detined a set
of fixed local axe$ xi (i = t, 2, 3)+ associated with a set of orthonormai vectors e, fi = 1,2, 3)>where e3 is normal to the u~defo~~ shelt surface, e, is perpen&cular to e, and the glohal X, axis, and e, is the cross product of e3 and e,. The rotation variables a and fl are taken as the rotations about the positive local xI axis and the negative local x2 axis, respectively. The rotation of the shell normal corresponding to the rotation variables is defhzed as the finite rotation ff2] resulting from the appEcation of the rotation vector 112, 131B = ae, - fie2 to the initial shell normal W, = e3. Parisch [3] and Surana [S] also define a set of fixed orthogonal local coordinate system and rotation variables in a similar way as that given by Oliver and Onate [9], but they chose different rotation axes and ~te~~t~ the rotations of the shed normal corresponding to the two rotation variables in a different way. Tho orientation of the shell normal is determined as follows: the initial shell normal OV,= e:, and the rotation vectar ae, are rotated to their new positions V; and e;, respectively by the appi~~t~on of the rotation vector &zr; then the vector V; is transparted to V, by the application of the rotation vector are{. Surana [8] aSsa tried different sequences of rotations to determine the orientation of the shell normal and compared their numerical efficiency. Bathe and Bolourchi fq, Hughes and Liu [6], and Milford and ~h~~b~~h 1%i] adopted the ~n~r~en~l rotations with respect to the x1 and x, axes of a moving coordinate system as the rotation variables. Tha x,-direction of”the moving coordinate system is chosen to coincide with the shell normal during deformation, the x1 axis is perpendicular to the x3 axis and the global X2 axis, the x2 axis defines a Cartesian
428
Kuo-MO HSIAOand YEH-RENCHEN
right-handed system. In order to avoid the use of the trigonometric representation for the rotation, the motion of the shell normal corresponding to the rotation variables is determined by the ‘vectorial rotation increments’ used in the linear case. The newly obtained normal vector is then projected radially to maintain the constant length assumption of the shell normal. It should be mentioned that the above procedure does not define rotations in the strict sense, but Cartesian displa~ments. In IS], it was mentioned that the rotation of the shell normal may be evaluated accurately by dividing the incremental rotation into subincrements. However, a large number of subincrements are required for large increment rotations. Instead of using the subincrement method, in this paper a strategy termed the finite rotation method is proposed to describe the motion of the shell normal. In this method the rotation variables tl and 0 are chosen to be the incremental rotations with respect to the x, and x, axes of a moving coordinate system which is rigidly tied to the middle surface of the shell. The initial xj direction of the moving coordinate system is taken as the normal of the undeformed shell surface, the initial x, axis is perpendicular to the x3 axis and tangential to a proper element side, and the xz axis is established by taking the cross product of the unit vectors associated with the xj and xi axes. The rotation of the moving coordinate system corresponding to the incremental rotation is determined by the application of the rotation vector B = ae, + /?e_,to the triad of the unit vectors associated with the moving coordinate system. Here, both the rotation increment between two successive increments and the rotation correction between two successive iterations are used as the incremental rotation. Since the choice of rotation axes may affect the numerical efficiency of nonlinear analysis, a mixed rotation method, which is a combination of the proposed finite rotation method and Bathe’s method [S], is also proposed to compare the ~rfo~an~e of the same rotation strategy with different rotation axes. Theoretically, it seems that all six types of rotation strategies mentioned above can handle arbitrarily large rotations and permit large rotations between successive increments. However, because they are employed in different shell elements, their convergence characteristics have not been compared in the literature. In order to test and compare the performance of these rotation strategies, they are embedded into a shell element whose formulation is based on that given by Ramm [2] using a total Lagrangian approach, and numerous numerical examples are studied. In [13] the previous convergent stress method was employed to update the geometric stiffness matrix for nonlinear analysis using a flat triangular shell element. It was reported that this method can improve the convergence characteristics for nonlinear analysis. Here this method is referred to as the PS
Fig. I. Nine-node degenerated shell element.
method and its performance is compared with that of the standard geometric stiffness update method, referred to as the CS method here. The numerical algorithm used here is an incremental-iterative method based on the Newton-Raphson method combined with a constant arc length of incremental displacement vector (7, 141. FINITE ELEMENT FORMULATION
In the following only the basic relations of the shell element are derived. A more detailed description may be obtained from [2]. Geometry and kinematics of the element
The geometry of the shell element is described in Fig. 1 by the coordinates and normal vectors at the nodes on the middle surface. The geometry and displacement representation are based on the kinematic assumption that the normal to the middle surface remains straight and inextensional during deformation. On the basis of the first assumption;the current configutation of the global Cartesian coordinate vector X = I&‘,, X,, X,) of a generic point of the shell element is expressed by X= 5
N’(r,s)X’+
f c
k=l
‘k=l
h&N&(r,s)V:,
(1)
where M is the total number of nodes of the shell, Nk is the matrix of shape functions associated with node k, Xk is the Cartesian coordinate vector of node k, r and s are the local curvilinear coordinates of the middle surface, t is the coordinate in the direction of the shell normal, hk is the thickness of the shell at node k, and Vi = fYfs, Vt3, Vi,) is the unit normal vector at node k. The displacement vector U = {iJ1, U,, U3f between the initial configuration 0 and the current configuration are specified to be consistent with eqn (I) U = 5 Nk(r, s)Uk + G 5 hkNk(r, s)(V; - “V$), (2) k=l
k-l
where Uk = {CJ:, Vi, .Y:} is the current displacement vector of node k, and “v: is the unit normal vector at node k for the undeformed configuration. It should be noted that because of the inextensional
Nanbar
anaiysis of shell structures
assumptions of the shell normal only two components of the nodal normal vector are independent. Thus, only two independent variables, referred to as nodal rotation variables, are required to determine the orientation of the nodal normal vector during deformation. So, i$ (I’= 1,2, 3) can he expressed as a function of two rotation variables c1and /I, and thus each node of the shell element has only five degrees of freedom: three translations in the direction of global axes and two rotation variabies ta be speci%d later. When the rotation variables and the rotation axes are chosen, the variations of displacements, which are required for the derivation af the element internal nodal force vector and element stiffness matrix, can be expressed as a function of the variations of nodal variables. From eqn (2), the variation of the disptacement vector can be written as SU = g N’(x, s)6Uk
kzl
where vectors f and g are obtained by differentiation of Vi with respect to the nodal rotation variables a“ and flk, respectively, and are given by f = 8V$/c3aR g =Ea~~]~~~.
14)
Assume that the incremental iterative method is used for the solution of nonlinear equilibrium equa. tions and that configuration I is known. Here ~~gu~t~on I may denote the equ~~ib~urn ~~g~t~oo of the previous increment or the iteration. configuration of the previous Let A& = {AU:, AU& AU:} and A+k = {A.ak, A/?“} be the given incremental nodal displacements and incremental nodal rotations, respectively, at nodes k (k = 1-M). In this paper when configuration I denotes the equilibrium con~~ra~o~ of the previous increment, AU’ and A#k denote the incremental nodal displacement and rotation increments between two successive increments. When the configuration I denotes the configuration of the previous iteration, Auk and Adik denote the nodal displacement and rotation correction between two successive iterations. Let Ixk, %I*, ‘v! and ‘+” = {‘a”, ‘/I”] denote the nodal coordinate vectors, nodal displacement vectors, unit shell normal vectors and nodal rotation vectors, respectively, at nodes k (k = 1-M) of an element at configuration I; the current values of nodal coordinates, nodal dispiacements and nodal rotation variables can be obtained by adding the corresponding incrments, so Xk=‘Xk+AUk Uk=%Jk+AUk Qlk= ‘rlik+ A#.
(5)
a!? x , ‘1
Dyformed normal
Fig. 2. Description of rotation for Ramm’s method.
At this point, an interesting question arises. Given the incremental nodat rotations, how are the current orientations of the shell normal Vi dete~in~~ It is a weil known fact that finite rotation components are aat vectorial and that the sum of finite rotations depends on the order of rotations. In geometric nonlinear analysis this has led to the adoption of many different types of rotation variables and different rotation strategies of the sheil normal to accommodate large rotation capability during the large deformation process. In the following, four different types of nodal rotation variables and rotation strategies proposed in the literature are reviewed and two sets of nodal rotation variables and rotation strategies are introduced. Note Chat ek in eqn (5) is only used for the A, 3 and C types of rotation strategies to determine the orientation of the shell normal; q: and A4k are used for the D, E and F types of rotation strategies. For convenience, the subscript k used to indicate node k is dropped throughout this paper to avoid confusion. (A) Ramm’s method. fn [Z] 01(Fig. 2), the angle between the global X,-X, plane and the plane formed by the shell normal and the global X, axis, and /3, the angle between the shell normal and the giobal X, axis are selected as the nodal rotational variables. The unit normal vector at node k in eqns (1) and (2) may now be written as V, = (cos b, sin /? cos a, sin @sin u }.
(6)
From eqns (4) and (fr), vectors f and g may be expressed as f=(O,
-siaasinjY,cosasinfij
g={-sinfl,cosacosB,sinacas~).
(7)
(B) Oliver’s m&ad. Oliver and Onate [PI defined a set of fixed local axes x, (i = f , 2, 3), associated with B set of orthonormal vectors ei (i = I, 2, 3), where e3 is normal to the undeformed shell surface, and e, is perpendicular to e3 and the global X2 axis. Note that in this study the vector eI is replaced by the unit vector which is perpendicdar to e5 and tangential to a proper element side. Vector e, is the cross product
430
KIJO-MO HSIAO and
t 0
normal and tangential to a proper element side. The rotation variables LY and /I are chosen as the rotations about the positive local xi and x, axes, respectively. The orientation of the shell normal is determined as follows: the initial shell normal qj = e3 and the rotation vector e, are rotated to their new positions Vi and ae;, respectively, as shown in Fig. 4, by the application of the rotation vector Be,; then the vector V; is transported to its final position V, by the application of the rotation vector ae;. The unit normal vector at node k may now be expressed as
x3
$ v3
Ov3
- P2 /
e
/
-
/
YEH-REN CHEN
x2
gel
V, = cos CIsin Be, - sin ae2 + cos a cos /3e3.
(10)
Xi
From eqns (3) and (lo), vectors f and g may be expressed as
Fig. 3. Description of rotation for Oliver’s method.
f = -sin CIsin De, - cos aeZ - sin a cos /3e3 of e, and e, . The rotation variables a and B are taken as the rotations about the positive local x, axis and the negative local x2 axis, respectively. The rotation of the shell normal corresponding to the rotation variables is defined as the finite rotation [12, 131 resulting from the application of the rotation vector 6 = ae, - Be2 to the initial shell normal q, = e3 as shown in Fig. 3. The unit normal vector at node k may now be written as V,=cosBe,+(l
-cos@(e;e,)e,
g = cos a cos Be, - cos a sin /?e3.
(11)
(D) &the’s method. Bathe and Bolourchi [S], Hughes and Liu [6] and Milford and Schnobrich [1 1] adopted the incremental rotations a and j with respect to the x, and x, axes of a moving coordinate system as the rotation variables. Let vectors e, (i = 1, 2, 3) be unit vectors associated with the moving coordinate axes xi (i = 1, 2, 3). Vector e, is chosen to coincide with the shell normal V, during deformation; vectors e, and e, are defined by
+ sin 8 (e, x e3) B cos Be, - i cos 8e, + cos 8e3, 8
=--
e, = I, x e3/ III, x e3 II (8)
and e2=e3 xe,,
where 6 = (a2 + fi2)“2 and e, = 0/t?. By substituting eqn (3) into (8), vectors f and g may be given as f= - $(cosO
g=-
e,-
-$(cost?-y)e,-
a
8sm0e,
e,
$costI+$y (
where i, is the unit vector in the global X2 direction. Let ‘V3= ‘es be the normal vector of configuration I at element node k, and a and B be the given incremental rotations. In order to avoid the use of the trigonometric representation for the rotation, the current shell normal V, is determined by the gener-
-y)e,
p2 sine a2 ~cosB+8’o
-
(12)
)
isin*e,.
(9)
(C) Parisch’s method. Parisch [3] and Surana [8] define a set of fixed orthogonal local coordinate systems and rotation variables in a similar way to that given by Oliver and Onate [9]. The local xi axis is chosen to be perpendicular to the initial shell
Xl
Fig. 4. Description of rotation for Parisch’s method.
Nonlinear analysis of shell
431
structW%i
incremental rotations at node k. A rotation vector is defined by 6 = a’e, + fl’e, = Bn,
(13
where 0 = 11 B 11is the angle of rotation and n e e/0 is the unit vector along the axis of the rotation, The current base vectors are determined by the application of the rotation vector 8 to base vectors ‘ei and expressed as (see Fig. 6) [l 11
e,scos@‘e,+ (1 -cosB)(r~-~e&r + sin @I x ‘e,) Fig. 5. Description of rotation for Bathe% method.
alization of the ‘vectorial rotation increments’ used in the linear case and given as follows (see Fig. 5): ‘V, -
” =e3= I/‘v3 -
a’e,4” B’q are2 -t_jPe,/I*
f= -e,
g=ei.
i = 1, 2, 3.
As mentioned earlier the current shell normal vector Vj - e,. From the definition of the rotations of the normal vector about the rotation axes corresponding to the variation of the rotation variables, vectors f and g in eqn (3) may be given by
(13)
It should be mentioned that the above procedure does not de&e rotations in the strict sense, but Cartesian displacements. After the vector VS is determined, the current rotation axes e, and eZ can be determined using eqn (12). From the definition of the rotations of the normal vector about the rotation axes corresponding to the variation of the rotation variables, vectors f and g in eqn (3) may be given by
06)
f= -e,
g=ei, where 9 and ei are given in eqn (16). (-F) Mixed rotationmethod. The choice of rotation axes about which the rotation variables are defined may affect the convergence rate of equilibrium iterations. In order to investigate this effect, a mixed rotation method is proposed here to determine base vectors. The shell normal vector V, = e3 is determined by eqn (16), but e, and e, are determined using eqn (12). Vectors f and g in eqn (3) are the same as those given in eqn (14).
(14)
(E) Finite rotation ~e~~u~. I&e a finite rotation method is proposed to describe the motion of the shell normal. In this method the rotation variabies are chosen to be the incremental rotations, a and jl with respect to the x1 and x2 axes of a moving coordinate system which is rigidly tied to the middle surface of the shell. The initial xg direction of the moving coordinate system is taken as the normal of the undeformed shell surface, the initial xi axis is perpendicular to the x, axis and tangential to a proper element side, and the xi axis is established by taking the cross product of the unit vectors associated with the X, and xi axes. Let vectors e, fi = 1,2,3) be base vectors, unit vectors associated with the moving coordinate axes xj (i = 1, 2, 3). Due to the definition of the moving coordinate system, vector e3 coincides with the shell normal Vj during deformation. The rotation of the moving coordinate system corresponding to the incremental rotations is determined by a finite rotation scheme similar to that proposed by Hsiao [19. Let ‘ei denote the base vectors at node k of configuration I, and 01 and /I be the given
The nonlinear expressed by
equilibrium
$ =F-IP=O,
equations
may
be
(W
where # is the unbalanced force between the internal nodal force vector F and the external nodal force dP;
Fig, 6, Deacriptian of rotation for finite rotation method.
1 is a loading parameter and P is a no~aiiz~ loading vector. The internal nodal force vector is obtained by summing up the element nodal force VeCtOorS. In this paper, a weighted Euchdean norm of the unbalanced Force [ f 5j is employed as the error measure of the equiiibr~um state during the e~uil~b~um iterations, and the convergence criterion is given by
where N is the number of degrees of fr~dom for the d~s~~et~~d structure and ptOSis a prescribed value of error tolerance. Unless it is stated otherwise, the error tolerance is set to 10 -4 in this paper,
An incremental-ite~t~ve method based on the ~ewto~~aphs~n method is adopted here. In order to deai with the hmit points and snap through, the arc length of the incrementat d~sp~a~m~nt vector is kept constant during the ~u~~~br~~rn iteration using Crisfield’s method f7, 141. The PS and CS methods are employed to usage the geometric st~~ness matrix and their ~rfo~an~es are compared. if the ~ui~ibrium con~gurat~on of the Ith increment is assumed to be known, the system infect stiffness matrix K+ then can be calculated at this configuration and an initial displacement increment Aq for the next increment may be obtained by using the Euter predictor as Aq = Ariq,,
QQ)
where Al is the initial in~rementa1 loading param&r and qr = K; ‘P is the tangential displa~ment of unit goading P. For all ~ncreme~fs other than the first, A;t is obtained in much the same way as that mentioned in f143, and is given by A.;1= f Aa(qq;q,)‘;“,
(zx)
while the sign is chosen following an ap~roa&h due to Bergan and Soreide [Id] in which the sign fatlows that of the previous ~nc~ment unbss the d~tE~j~~~~ of the tangent sti&ess matrix has changed the s&n, in which case a sign reversal is applied. rt\ais the incremental arc length used for next inurement, and is dete~ined by
number of itara~ons~ the cutout para~te~ C, and C, are chosen tu be 03 and 1.5, resp&iveIy, to prevent the yielding of an incremental displacement which is too large or too smalt. Using the djs~~~rn~nt &rement obtained in eqn (20) and ane oF the six rotation strategies of the shell normat deemed in the pre~nt paper, the current eonfi~rat~o~ of the shell structure can be determined. Then, the internal noda force vector F in eqn (tgj donated with the current e~~~gu~at~on can be caleutated, The foadirtg parameter d corresponding to the current ~on~gura~on is given by ,I-- li+ 81, where t’ is the convergent loading parameter at the fth ~n~~ment and Atl the loading parameter increment. Then the unbalanced force I& can be obtained from eqn (IQ If the convergence criterion feqn (I!?)] is not satisfied, a ~sp~~rn~~t correction r and B loading parameter correction SI 17, 141 are added to the previous Aq and A\A respectively to obtajn a new in~r~e~~ displacement and ~~~rnenta~ foading parameter for the next ~te~t~o~. The values of r and lil may be deterred by
and Aa 2 = (Aq + @If&q+ r)y
(zs)
where K, may be the tangent stifiness matrix at same known configuration. Two different st~ate~es are used to determine the Cflrrent Confi~ffrat~on. fn the first Stt’stegy, referred to as the dj~~acem~nt increment (LX) method, the current d~s~~~~m~nt i~~r~~~~~69 and the previous eq~~~b~~rn cor&guration are used Fcr eonstro.ct tb&! current configuration, fn the second strategy, referred to as the displa~m~nt correction (IX) method, the current con~guratio~ is determined using the displacement correction r and the configuration of the previous iteration. When these two methods are in~~~orat~ with any one of the first three rotation strategies of the shell normal mentioned in the present paper, ide~tjcal configuration wilt be obtained, because the rotation variables of these three methods are de&ed in fixed coordinate systems. But when these two methods are ~ornbjn~ with rotation strategy with rotation variabks de&d in axes of rotating coordinate system, diiYerent ~o~guratio~s may be obtained. This may a&et the convergence ~haracte~stics of ~u~~jb~urn iteration and wilt be investigated in the next section. Then, the internal nodal force corresponding to the new ~o~~~~tio~ obtained can be calcufated to initiate the next iteration. This procedure is repeated until the convergerme criterion is satisfied.
433
Nonlinear analysis af sheli structures
In the following, for convenience, XYYZZ is used to denote that X (A-F) type of rotation variables combined with W @I or DC) and ZZ (I’S or CS) methods is employed. E = 1.2X108kN/m2 II=0 L = 10.0 m b = 1.0 m h = 0.1 m l
r
Analyticul Present
Example 1. Cantilever beam with an end moment
solution
{ ) : NO. of ilaration
Displacement
(m)
Fig* 7. Cantilever beam with an end momeat.
geometric stiffness matrix for iteration strategies based on the Newton-Raphson method, In this paper, this method is referred to as the current stress (CS) method. As mentioned in [13], during the first few equilibrium iterations, the values of the stresses obtained from the current deformation may be several orders larger than their convergent values for certain problems with very large rotation increments. This may cause difficulty in convergence or even divtkgence for a large rotation increment using the CS method to update geometric stiffness matrix. in order to overcome this difficulty, a previous convergent stress (PS) method was proposed in [13] to update the geometric stiffness matrix. In this method the convergent stress of the previous increment is adopted to &c&ate the geometric sti&ess matrix. The convergence characteristics of the CS and PS methods will be investigated and compared using numerical examples in the next section.
This example is a cantilever beam subjected to a concentrated moment at the free end as shown in Fig. 7. The beam is idealized using five elements. For this two dimensional problem, the orientation of the shell normal determined by all types of the rotation variables and rotation strategies are identical except that determined by the DDI or the DDC method. It is found that if the CS method is used to update the geometric stiffness matrix, the convergence of the equilibrium iterations cannot be achieved even for small in~~ment steps. Thus, only the results obtained using the FS method are given for this example The results shown in Fig. 7 are obtained by using the APS method. It is seen that five increments are used. The number of iterations for each increment is also given in parenthesis beside each point on the graph. Nearly identicaI results are obtained for alI methods except that obtained by the DDIPS method. For the DDIPS method, 14 increments are used to reach M = 6.2 kNm and the average number of iterations per increment is 1I _
Exczmpfe2. ~~arnped~~~~edcirculaf a& The flamed-hinted deep arch depicted in Fig. 8 is subjected to a paint load at the crown. A 16 x 1 mesh is used to represent one half of the arch. The
E R @ t b
= = = = =
6Xi06
tb/‘in?
100 in. 215Q 1.0 in. 2.0 in,
j/j” section at A-A
041
9
u
8
In order to compare the performance of the six types of rotation variables and rotation strategies combined with the DI or DC method and PS or CS method, many numerical examples of plate and shell structures are studied here. The nine-node Lagrange element is employed for nonlinear analysis. The element matrices are calculated using the 2 x 2 x 2 Gauss integration. Crisfield’s constant arc length method [7, 141 combined with the Newton-Raphson method is used for equilibrium iteration exclusively. Unless otherwise stated, the error tolerance pIol is set to 10-4,
7 * (26)
cF
3
Present
2
( ;
: No. of iteration
1
0 Fig. 8.
0.2
0.4
0.6
0.8 U/R,V/R
1.0
1.2
1.4
Clamp&-hinged circular arch.
Kuo-MO HSIAOand
434 Table
1. Comparison of solution methods clamncd-hinaed circular arch Max.
ACS DDCPS DDIPS FDIPS
Error tolerance
No. of increments
Total No. of iterations
1.198 1.201 1.184 1.203
3 x IO-4 10-d lO-4 IO-4
30 6 25 6
223 55 283 56
in Fig. 8 are obtained by the FDIPS method using six increments. For this two dimensional case, except for the D type of rotation variables all types of rotation variables are equivalent. When the CS methods are used, the convergence of equilibrium iterations becomes very difficult and thus a looser error tolerance ptO,= 3 x 10m4is used. A comparison of performance among different methods is given in Table 1. The PS methods give significantly better convergence characteristics than the CS methods for this problem. As in the previous example, the DDCPS method shows significantly better performance than the DDIPS method. results
Table 2. Comparison of CS and CS matrix update method for column under a vertical load (B = 1)
for
U ii
Solution methods
YEH-REN CHEN
cr
The column shown in Fig. 9 is subjected to a vertical end load which is higher than the elastic buckling load. Two cases, B = 1 and B = 10 are considered. For the case when B = 1. five elements
PI L = 100 Lo= 20 h=l B = 1, 10 E = 12 Thickness = 1 p, = o.247x1o-3
B =
B=l ----
1.9580
.
3.3970 3.7554
0.0173 0.0864 0.2550 0.3115 0.5340 0.853 1 0.9067 1.2609 1.3058
ACS I
APS
; 8 6 8 5 10 9 5
are used and for the case when B = 10, six elements are used. To avoid numerical difficulty, the out of plane degrees of freedom are restrained for the case B = 10. Here the shell elements are used in the ‘plane stress mode’, therefore no rotation of the shell normal occurs during deformation. Hence, the choice of the type of the rotation variables does not affect the results. Only the performance of the CS and PS stiffness matrix update methods are compared for this example. Details of their convergence characteristics are given in Tables 2 and 3. As with the previous examples, the PS method gives much better convergence characteristics than the CS method. The results obtained using the APS method are shown in Fig. 9 together with the results given in [8] in which 20 load steps are used. Example 4. Cantilever plate with end load
t
The example considered here is a cantilever plate subjected to a concentrated end force as shown in Fig. 10. A four by three mesh is used for discretization. The plate is first analyzed using the APS and ACS methods to compare the convergence characteristics of the PS and CS matrix update methods. The results are given in Table 4. Clearly the PS method gives significantly better convergence characteristics than the CS method as with the previous
iTI s
r
10
-
Surana[B] Present
P
0.9317 1.0294 1.1417 1.1834 1.3829 1.8445
shown
Example 3. Column under a vertical end load
No. of iterations
P
Surana [ 81 Present
.
Table 3. Comparison of CS and CS matrix update method for column under a vertical load (B = 10) P
0
I
0.2
I
!I!
0.4
0.6
I
l
0.8
II
1.0
I!
1.2
U/L, V/L
Fig. 9. Column under a vertical load.
1
I
1.4
II
PC,
r,
1.4676 1.5628 1.7071 1.7245 2.0375 2.3708 2.6335 3.8676 3.9304
0.0177 0.0790 0.1941 0.2073 0.4061 0.5575 0.6479 0.8983 0.9066
No. of iterations ACS
APS
7 6 6 7 8 5 8 7 5
Nonlinear anafysis of sheif stndctures Table 5. Iterative performance for cantilever plate with end load (the table gives the number of iterations for each
30 m 40 m 0.4 m 1.2x1 d
r V.XIO~
40
u
7’ W . u = 0.3
iilCrement)
kN/mz G
APS BPS CPS DDIPS DDCPS EDIPS EDCPS FDIPS FDCPS
‘;24 c a x 4 16
4 4 4 8 4 7 4 7 4
4 4 4 6 4 4 % 4
4 4 4 6 4 5 4 5 4
Example 5. Simply supported plate subjected to lateral pressure
Deflection
(m>
Fig. 10. cantilever plate with end
load.
examples. Three i~~ements are used for the APS method and the load versus end displacements curves are shown in Fig. 10, together with those reported by Hsiao [ 131. The plate is later reanalyzed using different types of rotation variables in conjunction with the PS stiffness matrix update method, In order to compare their performance, the same fixed incrementaI arc length is used for each increment. The details are given in Table 5. Clearly the three DC methods give better convergence than the three cor~sponding DI methods as with the previous examples. The same ~rfo~an~ is obtained for the XPS (X = A, B, C) methods and the XDCPS (X = D, E, F) methods, From the results obtained using the EYYPS and FWPS methods (Table 5), it seems that the choice of the rotation axes only has a small effect on the convergence characteristics.
The pIate is subjected to a uniform pressure loading. Figure 11 shows the geometry and physical data of the plate. A quarter of the plate is analyzed using a two by two element mesh. It is found that the performance of all methods is similar for this example. However, the buoyance of the DC method is slightly better than that of the DI method, and the performance of the I’S method is mar~na~iy worse than that of the CS method. Note that the order of performance between the PS and the CS methods for this example is contrary to that for the previous examples in which the deflections are up
q = uniform pressure
/---a---l’ 3CD *s : ::
Levy 1171 Present
( )
l
:
No. of iteration l
(31
I
f
Table 4, Comparison of CS and Ps matrix update methods for cantilever nlate with end load No. of iterations
6.492 12.650 16.307 19.559 27.662 36.336 37.357
w&N
A(Js
7.897 13.903 16.719 18.810 22.730 25.548 25.818
6 5
APS 4 4
: 4 5
Centrat Deflection
Ratio W,/t
Fig. Il. Simply-supPorted Plate with a uniform lead.
436
E R L h Y
60
= = = = =
68.95 N/mm2 2540 mm 1569.8 mm 99.44 mm 0.3
r
40 s? t 30 a. 10 9w 20
10
0
Fig. 12.
-----
Surano [8] Horrigmoe and Bergan [Is]
_
Present ( 1
50
No. of iteration
100 $50 ZOO 250 Central Deflection ~c(mm~
300
Hinged spherical shell (vertical cutting).
Also shown in Fig. 12 are the solutions given in [&l&19], and in Fig. 13 are the solutions given in (81.The solutions given in [8] are obtained using 12 increments with about 50 iterations. Example 7. Hinged cylindrical shell (h = 12.7 mm) Figure 14 shows a circular cylindrical shell subjected to a concentrated central load on the convex side. The longitudinal boundaries are hinged and immovable, whereas the curved edges are complete& free, Here the symmet~~i response is studied by using a two by two mesh for one-quarter of the shell. This example is analyzed using all the six types of rotation variables. It is found that the performance of all methods is similar for this example. However, the performance of the DC method is slightly better than that of the DI method as with the previous examples, and the performance of the CS method is marginally better than that of the PS method as with examples 5 and 6. The results obtained using the EDICS method are given in Fig, 14 with the solutions reported by Surana [81 using the C type of rotation variables and Oliver and Onate [9] using the B type of rotation variables. The results of the present study are obtained by using three increments and 15 iterations. The solutions given in [8] are obtained using 12 increments and 54 iterations.
to the order of ma~itude of the major st~ctural Example 8. Hinged ~yl~nd~~~al shell (h = 6.35 mm) dimensions. In this example the shell thickness of the previous The results obtained using the ACS method are given in Fig. 11 with Levy’s classical solution [17]. example is reduced by one-half, i.e. k = 6.35 mm. The The results of the present study are obtained by using three increments and 10 iterations,
I-he we consider a spherical shell subjected to a concentrated load at the crown. All four edges are hinged and immovable. Two cases are considered: (1) the spherical shell is formed by four vertical cutting planes (Fig. 12) and (2) the spherical shell is formed by four cutting planes passing through the center of the sphere (Fig. 13). A two by two element mesh is used for one-quarter of the shell. Only the 3, C and E types of rotation variables can describe the boundary conditions of this example properly, therefore only these three types of rotation variables are used. Nearly the same performance is obtained for these three types of rotation variables, The results obtained using the EDICS method for both cases are shown in Figs 12 and 13, respectively. It is seen that only three increments with 13 iterations are used for the first case and three increments with 14 iterations are used for the second case. The solid curves shown in Figs 12 and 13 are also obtained by the present study using 17 increments to obtain complete loadingdeflection curves. As with the previous example, the CS method gives a marginally better performance than the PS method.
& R L ;
= = = 7 -
68.95 N/mm’ 2540 mm 1569.8 mm ;9;45 mm
.
I’
P Xx
:: 20 f
____
I’
_
~~
y
.
( , ,;;;,;
0
f3.
Present
;lym;,
3~*
50 Central
Fig.
Surano [ 8}
Deflection
W,(mm>
Hinged spherical shell Qsyramid cutting).
Nonlinear anaIysis of
R = 2540 mm L = 508 mm tl = 12.7 mm
R
437
shell structures ----^-
tlorrigmoe
---.
Ramm [20]
and @argon [18]
.
6 = 0.1 rodions E = &I0275 l&/mm2
V = 0.3 4
-
Oliver [ 9 ]
-----
Sunara [S]
(5) ( >
No. of itemtion
3
5
10
15
20
25
30
Deflection (mm)
Fig. 15. Hinged cylindrical shell (h = 6.35 mm).
Fig. 14. Hinged cylindrical shell (h = 12.7 mm).
remaining geometric data, physical data and bout& at-y conditions are unchanged. A quarter of the shell is analyzed using a two by two element mesh. The shell is first analyzed using the APS and ACS methods to compare the convergence characteristics of the PS and CS matrix update methods. It is found that the CS method gives s~gni~~antly better convergence ~hara~~e~stj~ than the PS method far this highly nonlinear problem. Three increments with 21 iterations are used for the ACS method and 13 increments with 73 iterations are used for the APS method. The load versus end displacements curves obtained using the ACS method are shown in Fig. 15, together with those reported in f18, ZO].The sohd fine shown in Fig. 15 is also obtained by the present study using a 3 x 3 element mesh and 27 increments to obtain complete loading-deflection curves. This example was studied by Surana [8] using the C type of rotation variables The rest&s reported in fg] (not shown in Fig. 15) are obtained by using 2 1 increments and 98 iterations. The shell is later reanalyzed using different types of rotation variables in conjunction with the CS stiffness matrix update method. In order to compare their performance. The same fixed incremental arc length is used for each increment. The details are given in Table 6. Clearly the three DC methods give better convergence than the three corresponding Df methods as with the previous examples. Nearly the same performances are obtained for the XCS (X = A, B, C) methods and the XDCCS (X = D, E, F) methods. From the results obtained using the EYYPS
(YV = DI, DC) and FYYPS methods (Table 5), it seems that the choice of the rotation axes has no effect on the convergence characteristics For this example. CONCLUSIONS
Six types of rotation variables and rotation strategies are employed to describe the rotation of the shell normal. In particular, a finite rotation method is proposed and tested. In this method the rotation variables are chosen to be the incremental rotations with respect to the x1 and x, axes of a moving coordinate system which is rigidly tied to the middle surface of the shell. Both the rotation increments between two successive increments and the rotation
Table 6. Iterative performance for hinged cyfindricai bad (the table gives the number of iterations for each imxement)
G APS BPS CPS DDIPS
7 6 6 8
7 7 7 8
7 7 7
EDIPS DDCPS EDCPS FDIPS FDCPS
7 7 8 7
:
; 7 g 7
; 7
8
438
Kuo-MO Hsuo and YEH-RENCHEN
corrections between two successive iterations are used as the incremental rotations (rotation variables) to update the orientation of the shell normal. The previous convergent stress (PS) method and the current stress (CS) method are employed to update the geometric stiffness matrix and their performances are compared. Prom the numerical experiments that were conducted the following conclusions can be drawn. (1) The six types of rotation variables give a similar performance. (2) The displacement correction (DC) method always gives a better performance than the displacement increment (DI) method. (3) The previous convergent stress (PS) method gives significantly better convergence characteristics than the current stress (CS) method for the problems with very large displacements and rotations up to the order of magnitude of the major structural dimensions, but for the problems in which the deflections are the order of several thickness, the CS method is more efficient and is generally a little more efficient than the PS method. (4) The choice of the rotation axes has little effect on the convergence characteristics. (5) Only the B, C and E types of rotation variables can describe the boundary conditions of the spherical shells given in example 6. (6) Although no example is given here, it seems that only the E type of rotation variables proposed in the present paper can handle moving boundary conditions.
REFERENCES 1. S. Ahmad, B. M. Iron and 0. C. Zienkiewicz, Analysis
of thick and thin shell structures by curved finite element. Inr. J. Numer. Meth. Engng 2, 419-451 (1970). 2. R. Ramm, A plate/shell element for large deflections and rotations. In Formulations and Computational Algorithms in Finite Element Analysis, U.S.-Germany Symp. (Edited by K. J. Bathe, J. T. Cklen and W. Wunderlich), pp. 264293. MIT Press, Cambridge, MA (1977). 3. II. Par&h, Geometrical nonlinear analysis of shells. Comput. Meth. appl. Mech. Engng 14, 159-178 (1978). 4. H. Parisch, Large displacements of shells including
5. 6.
7. 8.
9.
material nonlinearities. Comput. Merh. appl. Mech. Engng 27, 183-214 (1981). K. J. Bathe and S. Bolourchi, A geometric and material nonlinear plate and shell element. Compur. Strucf. 11, 2348 (1980). T. J. R. Hughes and W. K. Liu, Nonlinear finite element analysis of shells: Part I. Three-dimensional shells; Part II. Two-dimensional shells. Comnut. Mesh. anol. Mech. Engag 26, 331-362 and 27, 1671181 (1981)” T. Y. Chang and K. Sawamiphakdi, Large defIection and post-buckling analysis of shell structure. Compur. Merh. appl. Mech. Engng 32, 311-326 (1982). K. S. Surana, Geometrically nonlinear fo~ulation for the curved shell elements. Int. J. Numer. Merh. Engng 19, 581-615 (1983). J. Oliver and E. Onate, A total Lagrangian formulation for the geometrically nonlinear analysis of structures using finite elements. Part I. Two dimensional problems: shell and plate structures. Inr. J. Numer. Meth. Engng 20, 2253-228 1 (1984).
10. K. J. Bathe and E. N. Dvorkin, A formulation of general shell ‘elements-the use of mixed interpolation of tensorial components. Int. J. Numer. Meth. Engng 22, 697-722 (1986).
11. R. V. Milford and W. C. Schnobrich, Degenerated isoparametric finite elements using explicit integration. Inf. J. Numer. Meth. Engng 23, 133-154 (1986). 12. H. Goldstein, Ciassicai Mechanics, 2nd Bdn. Addision-Wesley, London (1980). 13. K. M. Hsiao, Nonlinear analysis of general shell structures by flat triangular shell element. Compuf. Srrucr. 25, 665675 (1987). 14. M. A. Crisfield, A fast incremental/iterative solution procedure that handles ‘snap-through’. Comput. Strucr. 13, 5562 (1981).
15. A. K. Noor and J. M. Peters, Tracing post-limit point with reduced basis technique. Comput. Meth. appl. Mech. Engng 28, 217-240 (1981). 16. P. G. Bergan and T. Soreide, Solution of large displacement and instability problems using the current stiffness parameter. In Finite Elements in Non-linear Me&attics (Edited by P. G. Bergan et al.). Tapir Press, Trondheim, Norway (1978). 17. S. Levy, Rending of rectangular plates with large deflection, NACA Tech. Note 846 (1942). 18. G. Horrigmoe and P. G. Bergan, Nonlinear analysis of free-form shells by flat finite elements. Comput. Merh. appf. Mech. Engng 16, 11-35 (1978). 19. A. Pica and R. D. Wood, Postbuckling behaviour of plates and shells using a Mindlin shallow shell formulation. Comput. Struct. 12, 759-768 (1980). 20. E. Ramm, Strategies for tracing the nonlinear response near limit points. In Nonlinear Finite EIement Analysis in Structural Mechanics (Edited by W. Wunderlich, E. Stein and K. J. Bathe), pp. 63-89. Springer, Berlin (1981).