- Email: [email protected]
Implementation of assumed strain degenerated shell elements
~MFLEMENTATION
OF ASSUMED STRAIN SHELL ELEMENTS?
DEGENERATED
Hou CHENG HUANC$ Department of Mechanics, Tianjin University, China (Rereit$ed 23 firne 1986) Abstract-The implementation of the new degenerated shell elements with assumed transverse shear and membrane strains is described in this paper. An explanation of locking ~haviour and also the location of sampling points For the assumed strain fields is given in the present work. In the fo~u~ation of the new elements, assumed transverse shear strains in the naturat coordinate system are used to overcome the shear locking problem and assumed membrane strains in the orthogonal curvilinear coordinate system are applied to avoid membrane locking behaviour. Both tinear and nonii~e~r problems are considered~
1. INTRODUCTION originaf degenerated shell element performs reasonably well for moderately thick sheil situations [I]. However, for thin sheifs when full integration is used to evaluate the stiffness matrix, overstiff solutions are often produced due to shear or membrane locking. Attempts have been made to correct this behaviour by use of reduced or selective integration techniques[2-4]. Such schemes are not always successful in overcoming locking behaviour and the resulting solutions may still be overstiff for problems with highly constrained boundaries especially when coarse meshes are used. Furthermore, for problems with lightly constrained boundaries, element mechanisms or spurious zero energy modes may form. These mechanisms can spread from element to element causing either rank deficiency in the overall stiffness matrix and con~quently no solution, or the even more dangerous situation in which the solution obtained may be polluted by a near-m~hanism. Many research workers have attempted to improve the performance of Mindlin plate elements and degenerated shell elements [5]. From this work it appears that the ideal Mindlin-type element should: (1) converge, (2) not lock, (3) contain no mechanisms, (4) be capable of providing accurate displacements and stress resultants, (5) be insensitive to element distortions, (6) be invariant to the direction of coordinate system, (7) not be based on nume~cally adjusted factors, and (8) be easy to implement and The
Use.
In the formufation of the new degenerated she8 elements, assumed transverse shear strains in the natural coordinate system are used to overcome the shear locking problem and assumed membrane t’fhis is an extended version of the paper presented at the Conference of ICCMES-TOKYO, Japan, May 1986. $ Lecturer, Ph.D.. presently a research scholar at University College of Swansea, U.K. 147
strains in the orthogonal curvilinear coordinate system, are applied to avoid membrane locking behaviour. It has been shown that the new elements are free from any major defect [5].
According to ~indlin-ty~ theories, the rotations are independent of the displacements in both plates and shells. The transverse shear effects should gradually diminish as the plate/shell thickness becomes extremely thin. If we could solve the governing differential equations directly, shear Iocking behaviour should not appear. However, the shear strains obtained using a nurne~~l method, such as the finite element method, may be inaccurate locally, though average shear strains over a particular region are reasonable. Therefore, ii is possible to find some points at which values of the shear strains may represent the average distribution of the shear strain fields at the particular region. As we know, when the finite element method is used in conjunction with Mindlin-type formulations, the shear strain energy is always included and must be positive. Consequently, as the thickness of the plate and shell becomes extremely thin, the shear strain energy can be magnified unreasonably even though the average value of the shear strains over the area tends to zero 161. Consider for instance Fig. I, which shows a rectangular b&ear element subjected to pure bending moments at opposite edges. Apparently, there exists a ‘parasitic’ shear deformation, whereas in reality there should be no shear deformation. However, it is noticed that
Hou CHESGHCASG
148
Fig. 1. Rectangular bilinear element subjected to pure bending moments. whereas the element shear strain energy Gh
where yer, Y,,~are obtained from finite element soiutions and yf& Y; are the transverse shear strains obtained from exact anaiyticai soiutions. For thin plate and shell cases, we have
(2)
JA
y;: dx dy >>0.
It is also noted the values of the shear strain at the line x = 0 are correct. This explains the reason for the success of the use of reduced integration (1 x 1) for the shear strain energy [Z). According to Mindiin-type theory for general thin plate and shell problems, there is always some shear deformation at each point, even for the higher order Lagrangian elements. Therefore, when a full integration scheme is employed, the shear strain energy contribution will dominate the Total Potential Energy and the shear locking problem cannot be avoided. In the following, we take a g-node element as an example for the Lagrangian elements and an 8-node as an example for the Serendipity elements. 3. ASSUMED TRANSVERSE SHEAR STRAIN FIELD
An artificial method for the eitmination of shear locking is to interpolate new shear strain fields from the strain values at the sampling points which are appropriately located in individual etements. For the 9-node Lagrangian element the polynomial terms for the two contributions to each transverse shear strain evaluated from associated displacement fields do not match (51.The substitute shear strains yet and 7,’ should be at least taken as follows:
Yct+Y;‘O.
P,: -) r; + 0.
(5)
When a Mindiin-type hypothesis is applied, then equation (5) cannot be satisfied at ail points. However, equation (5) could be satisfied in an ‘average’ sense as follows:
J~,~?,,dtdl=J~~~~dtdn
-to.
(6)
In the finite element method, AA could be a particular area in an individual element. Now, we consider the special case of thin plates and shells. For the g-node Lagrangian element, the substitute shear strain field has been taken as a linear function of { for ycr and a linear function of 9 for Y,,~. It is assumed that
II J J -1
=0, J’‘r
ycr dc = 0,
0
-I
(7)
0
dn = 0.
y,,xdq==O,
(8)
0
Now, we try to find the value of ( at which the shear strain ycr is equal to zero and therefore we have to solve the following equation:
Y<,= b, t b,< + b,q + b‘& + b& + b&t$ F,, = cr + cr5*+ c,q + c&
+ c&2+
c&l.
(9) (3)
If the substitute strain fields are expressed as poiy nomials of a lower degree than those in equation (3), then zero energy modes will appear. This has happened in the element with seiective integration. It is noted that ‘it, is linear in 5 and quadratic in 9 and T,t is linear in 4 and quadratic in {.
which gives At2+Bt‘-t-C=O,
(10)
where A = H, (B. - 28, + @:,)I2f Hz@:4- 26, + B&2 f HA% - 25, f 4w B = H,[w, f w3 - 2w, f (4, - f&,)/2]
4. LOCATION OF THE SAMPLING POINTS FOR THE SHEAR STRAINS
In general, for plate and shell situations, following relations are expected: ‘J:, + $V
Y4t-+ r,“,
f H&q + wg - 2% -t (e,, - 4J/2] the
+ HJw, + w7 - 2w, + (k, - 4,),‘2] C=H,[(WJ-
(4)
-5 li,Kw,
d/2 -
+
w7)/2
421
f
+
461
ff*[(w,
-
d/2
f
e,1 (11)
of assumed strain degenerated shell elements
Implementation
a
A’I
a
(1.21 I 0 ;
_
b
0
/
-m
IO (2.1) I /
(3,21
0 (1,11
t ’ (**‘)
J
-
-
/ r---g (2.21
1,,i3.21
” !1,21
t 1
I
I
i
c”b
I
(1.1)
0 (3,ll --_-_+----_)c I*,*)
149
(3.11
;
J?r
‘655
Fig. 2. Interpolation points (i, j) for the 9-node Lagrangian element.
and H, = ~,(l + ~i)/2 + (1 - tl*)(l - tl!),
(12)
in which & and vi are the natural coordinates of the element nodes. However, from equation (7) it can be shown that B = 0, A = -3C and therefore equation (10) becomes t*= l/3. So, under the assumption of equation (7) it can be concluded that if 5 = fa (where a = 3-l’*, which are the two Gauss points in the { direction) then ycl = 0 for very thin plates and shells. Similarly, it can be shown that under the assumption of equation (8) if q = fa (which are the two Gauss points in the 7 direction) then yV = 0 for very thin plates and shells. As mentioned in equation (3), six sampling points are required for the representation of the substitute shear strains. For ytr the sampling points should be taken at the lines 5 = a and < = --(I on which three points must be taken. One of the three is q = 0. In order to guarantee the continuity of shear strains between two adjacent elements, the other two points should be 7 = 1 and q = - 1 at the lines 5 = a and 5 = --a. The six sampling points for ya are the three points < = 1, 5 = 0 and C = - 1 at the two lines n = a and n = -a (see Fig. 2). Therefore, we have
and the terms y$ and ;lyY are the transverse shear strains evaluated from the displacement field, a = 3-l’* and b = 1. Similarly, for the 8-node Serendipity element to avoid zero energy modes the substitute shear strains y. and yllr should be taken as polynomials of at the least the following degree [A: Fcr = b, + b25 + b,rj + b&j + b,q’ iisr= Cl + c25 + c,rJ + c&j + c55*.
(15)
Consequently, five sampling points are needed to construct the substitute shear strains. Apparently, along the edges of the element, deformations which are identical to those in the g-node element are expected, and therefore four of the sampling points can be chosen as two points q = 1 and n = - 1 at the two lines r = a and 5 = --a for y(,, and two points { = 1 and 5 = - 1 at the two lines q = a and q = --II. From the pure bending test mentioned in Sec. 2 it is found that the transverse shear strains are equal to zero at the centre of the element and therefore the centre point is introduced as the fifth point for both ycr and y,,*(see also [5]). In fact, for the 4-node element, if it is assumed that
I’-I
yet d5 = 0
(16)
~t,rdrl = 0,
(17)
and I I -1 where Pz(z)=l-
0
it can be shown that y. = 0 when C = 0 and yV = 0 when q = 0 for very thin plates. Therefore, the sampling points for the substitute shear strains can be chosen as those given in [8].
2
;
5. ASSUMED MEMBRANE STRAIN FIELD
Q,W=;(l+$
,,,)=;(I-~)
(14)
The membrane and shear stiffness terms are usually of the same order of magnitude. In shells, the stiffness associated with the membrane strain energy, though
150
HOU CHESG HUAHG Fig. 3). To define the local r and s directions at a typical midsurface point P. it is assumed that PS and PR are intersections of the plane tangential to the shell midsurface and the planes parallel to the r,, - to and s, - to planes drawn through point P respectively. Then the bisectors PB of angles R-P-S and r-P-s are made coincident in the tangent plane (see Fig. 4). The basic vectors gi of the orthogonal coordinate system are defined as follows:
Fig. 3. Symmetry between the natural coordinate system and the orthogonal curvilinear coordinate system.
important, should not dominate the totaf stiffness. It is known that the membrane strains can only be separated from the bending strains in the local Cartesian coordinate system (x’, J’, z’) [l, 91 in which the x’-y’ plane is tangent to the shell mid-surface. Therefore, it is expected that the use of an enhanced interpolation of the membrane strains in an orthogonal curvilinear coordinate system (r, s, t) [S] (in which the (r, s, 0) surface coincides with the mid-surface of the shell element) can help to eliminate membrane locking behaviour. It should be noted that symmetry between the natural coordinate system and the orthogonal curvilinear coordinate system must exist to ensure that this special coordinate system is uniquely defined. The origin of the orthogonal curvilinear coordinate system is set at the centre of each element. The axis t is normal to the shell mid-surface and the (r, s, 0) surface coincides with the ({, q, 0) surface at the shell midsurface. In order to guarantee the required symmetry, the bisectors of angles rO-o-s,, and &-o-~O are made coincident at the centre of each element (see
where g, are unit vectors [j]. From an assumption given in Sec. 4, it is deduced that the values of the membrane strains at lines q = a, -a (where a = 3-‘!?) give the q=-aandt=a,<= best representation of the membrane strain field in each element. Therefore, in the orthogonal curvilinear coordinate system, the enhanced interpolation for e,,,_,, and ie:L,, is the same as that for the y:, and the enhanced inte~oiation for e,_, and iegL,s is the same as that for y,,, but they are expressed in different coordinate systems. The details of the modification of the shell strains are given in [5]. However, the evaluation of the shell element stiffness is conducted in the local coordinate system (x’, y’, z’) in which the x’-J’ plane is tangential to the surface z = constant, as mentioned in [l]. For a shell element of uniform thickness the r, s, I and x’, y’, z’ coordinate systems may be taken as identical systems. All formulations are the same as those mentioned in [I] except for the definition of the local coordinate system. For a general shell element with the assumptions adopted in the present work, a transformation between the orthogonal curvilinear coordinate system and the local coordinate system should be
Fig. 4. ~e~nition of the orthogonal curvilinear coordinate system.
Implementation of assumed strain degenerated shell elements established. The tensor ~ansfo~ation foIIows: ax’ a.d ad
matrix #is as
z-z-z-
B=
dy'
ay'
ay'
drasdt azr at!
af
the stress symbols denote that the stresses at time t are referred to the configuration at time 0. Incremental second Piola-Kirchhoff stresses, ,,A$,, can be expressed as
(19)
(23)
$Si, = &k, OAEW
,,Aq,, are nonlinear incremental strains which are neglected in the *AE, expression and can be described as
arasat
6. GEOMETRICALLY
I51
NONLINEAl RI
(24)
In the geomet~~ally nonlinear analysis both Total Lagrangian (T.L.) and Updated Lagrangian (U.L.) formulations have been considered. In the T.L. formulation, all current static and kinematic variables are referred to the initial configuration before the deformation (say at time 0), whereas in the U.L. formulation, all the current variables are referred to the current configuration for each new load increment (say at time t). Both the T.L. and U.L. formulations include all kinematic nonlinear effects due to large displacements and rotations, but large strain behaviour is not taken into account in the present work f IO]. The local Cartesian coordinate system is used for both geometrically nonlinear formulations. The Green strain tensor can be expressed as follows
In equation (21) ‘+dW is the external incremental virtual work at the current load increment. According to the new element formulation, the shear strain should be interpolated from the values at certain sampling points as shown in equations (13) and (14). Therefore, it can be written that
(25)
p=lq=l
where oAEql = oAcF + oAeT oAEQ = ,A$ + oAe:r.
In the new element formulation the total shear strains have to be interpolated at every iteration in the Total Lagrangian formulation. It appears that the use of the Updated Lagrangian formulation is more appropriate with the interpolation concept used for the new element. 6.1 Total Lagra~gia~ ~or~ul~tio~ (T, L.) In the Total Lagrangian formulation, all variables are referred to the initial configuration at time 0 and the governing equation with incremental terms is
(26)
For the g-node element, p = 1, 2 and q = 1, 3. For the membrane strains, the interpolation in the orthogonal curvilinear coordinate system is only employed to the linear parts in equation (22) as shown in the linear case. ,AE$ and ,AE$ can be obtained from oAEij given by equation (22) by tensor transfo~ation. As the virtual displacements 6Au may be arbitrary, (21) can be rewritten in finite element form as
(I
&BTDoBodv+~oC?“]oCOdL.)Ad
=“SW-
-_ ‘+f -
s
;ISjjS oAEjj Ode, (2 1)
where &E,, is the incremental Green strain tensor which can be divided into two parts so that &E,
= ,AE{~+ ,Ae,
OAsCij = f (oA\u,,+ 0AUj.i)
,B”u dr,
where the incremental strain-displacement matrix ,,B = ,B + oBr,and the matrices oB and o& are associated with oAZand oA8 respectively. Therefore, the ,B matrix is different from that obtained in small deflection analyses. ‘+f are current applied loads and Ad is the incremental nodal displacement vector. Furthermore, it can be shown that KrAd = r,
deij = f(0~k.i o&j + 0Aut.i 0~)
(22) where KT is the tangential
and iSij are total second Piola-Kirchhoff stresses, the superscript 0 and subscript t given in front of
(27)
I
residual forces and K,=K,,-f-K,+K,
stiffness matrix,
(28) r is
152
Hoti CHENGHUANG
in which id=
which are referred to the configuration at time t. The Cauchy stresses are represented as rij. The incremental Green strains, ,A&, can be simply written as
oBTD,,g Odv
and *G is defined as follows oGijAcfi = oAuiQX,for
i = I, 3
oG,,Adi = &ui,,,
for
i = 1, 3
and
k = i+ 3
oGkiAdj = ,Aq,,
for
i = 1, 3
and
k = i + 6,
where ,Ar, and ,Aqi, are the linear and quadratic parts of ,AE,. As the quadratic strain expressions can be neglected in the first term of equation (34), the interpolations for the shear and membrane strains are the same as for the linear case. Equation (34) can be written in finite element form as
+
(301 where the maximum size of sum index j depends on the number of degrees of freedom in each element. For instance, for a 9-node shell eiement j is equal to 45. The stress vector u is
la1 = adds
;;;
:~I.
(32)
j,GT[rl,G'dc)Ad
,BTDIB’dv
= r+f-
,BTr ‘dv
where ,K is associated with ,Acij in which the shear and membrane strains have been substituted according to the formulation of the new shell elements. Therefore,
08)
K,Ad = I,
where S,, = 0 and
[ 1 1 0
I =o
0
10.
0
0
(37)
where KT = K0 + K, and (33)
1
,BT D J ‘do
K, = s
The second parts of shear strain in equation (26) can be neglected and the interpolation of shear strain is only required for the first parts of the shear strain as in the linear case.
,GT [r] ,G ‘do
s-
(39)
in which $G has the same form as in (30) and
6.2 Updated Lagrangian formulation (U.L.) Since dEcE,, and oAE,,, have to be interpolated at iteration, it appears that the use of the Updated Lagrangian formulation is more appropriate for use with the interpolation concept used for the new element. An incrementai form of the principle of virtual work based on the Updated Lagrangian formulation is as follows: every
rCi,k,,AEk, 6 ,AEij ‘do c I
where rj) = 0. 7.
‘~ij6 ,Aqij ‘dV s
7.1
= “6 W -
‘qj 6 ,AE, ‘du
(34)
I. where ‘+6 W is the external incremental virtual work at the current load increment and second PiolaKirchhofT stresses, &IS,, are AS, = rCijk/M’k/ 9
(35)
NUMERICAL
APPLICATIONS
Numerical example for linear onal_vsis
Numerical tests for element characteristics have been conducted in [S, 91. A more stringent example is the pinched hemisphe~cal she11subjected to concentrated forces shown in Fig. 5. TWOequal and opposite concentrated loads are applied in the x and y directions. By taking advantage of the double symmetry which exists in the problem, only a quadrant of the shell is considered with meshes of 1 x 1, 2 x 2, 4 x 4
Implementation
153
of assumed strain degenerated shell elements
2
QUAD 9" OUAO9S
1
--*-
Number of nodes per side
Fig. 6. Convergence of deflections at the loading point.
7;2 Numerical examples for geometrically nonlinear analysis Fig. 5. Spherical shell: radius = 10.0, thickness = 0.04, E = 6.825 x IO’, u = 0.3. and 8 x 8 elements (or n = 3, 5, 9, 17 nodes along each side). In this problem, QUAD9**, the new g-node Lagrangian element is used. A comparison is provided with the performance of QUAD9S, the g-node Lagrangian element with selective integration for the transverse shear stiffness. The locking of QUAD9S is evident in Fig. 6. For the 8 x 8 mesh, the results are less than 30% of the exact solution [I 11. The new element QUAD9** converges.
For the following numerical applications, a modified Newton-Raphson scheme is adopted and a displacement convergence criterion is employed. The preset tolerance is 0.001. 7.2. I Long cantilever. Five QUAD9** elements are used to idealize a cantilever subjected to an end moment as shown in Fig. 7. This problem is chosen to examine the behaviour of the QUAD9** when used for a geometrically nonlinear problem with large displacements and rotations. Here, an Updated Lagrangian formulation is adopted. In the analytical solution to this problem the slope is prescribed at
EI/L=l.O h/L =O.Ol v =o.o
a -1
-2 -3 -t -5 -6_ -1. -8 -2
I
-1
0
I
2
3
t
5
6
7
a
Fig. 7. Elastic cantilever beam, deformed shape.
9
18
Hoti CHENGHUASG
154
0.8 FOUR HETEROSIS .....O,,... FOUR WA0 9 ** -*_--_ WAY LINEAR SOLUTlON-
06
LO
30
20
1.0
00
CENTRE OlSPLACEflENTlmml
Fig. 8. Clamped square plate under uniform load.
the free end as 0 = ML/EI and the neutral axis is deformed into a circle with radius R = EI/M. L is the beam length and M is the end moment and El is the flexural rigidity. For a cantilever of length L = 10, the analytical solution indicates that the radius of the final circle is 3.183, whereas the radius obtained using the new element, QUAD9 **, is 3.185. In Fig. 7, it is shown that the error is reasonably small when the load approaches the value of 27rEZ/L which is divided into 100 load increments. 7.2.2 Clamped square plate. Four QUAD9** are used to idealize a symmetric quadrant of a clamped square plate subjected to a uniformly distributed load as shown in Fig. 8. The number of the increments
R=ZSLOmm L:ZSLmm ha12.7mm k310275 v.o.3
N/mm'
adopted is 15. Lateral displacements at the plate centre obtained using the new element are compared with the results obtained by the classical RayleighRitz solution presented by Way [12]. The results obtained by both T.L. and U.L. formulations are almost identical. Figure 8 shows that the central displacements obtained in the present study lie between those obtained by the classical solution and those obtained using heterosis elements [ 131. It is worth noting that the results obtained using the Lagrangian nine node element with selective integration, QUAD9S, are almost the same as those obtained using QUAD9**. 7.2.3 Hinged cylindrical shell. The hinged cylindrical shell. shown in Fig. 9, has been studied by Horrigmoe [14]. The straight edges are hinged and the curved edges are free. The shell is subjected to a concentrated central load on the convex side. A 2 x 2 mesh of QUAD9** elements is used to represent a symmetric quadrant of the shell. Here also, identical results are obtained using both T.L. and U.L. formulations with 20 load increments. In Fig. 9 the load-displacement response obtained using QUAD9** elements is compared with that given in Horrigmoe’s study in which a max. load of 2.22 kN was obtained. In the present study the max. load is 2.21 kN compared with a value of 2.20 kN obtained by Hughes [ 131 using four heterosis elements. REFERENCES
I. S. Ahmad, B. H. Irons and 0. C. Zienkiewicz, Analysis
of thick and thin shell structures by curved finite elements. Inr. J. numer. Merh. Engng 2, 419-451 (1970). 2. 0. C. Zienkiewicz, R. L. Taylor and J. M. Too, Reduced integration techniques in general analysis of plates and shells. Inr. J. numer. Meth. Engng 3, 275-290 (1971). 3. E. Onate, E. Hinton and N. Glover, Techniques for improving the performance of Ahmad shell elements. Proc. Int. Conf. on Applied Numerical Modelling, Madrid. Pentech Press (1978). 4. E. D. L. Pugh, E. Hinton and 0. C. Zienkiewicz, A study of quadrilateral plate bending elements with ‘reduced’ integration. Inr. J. numer. Meth. Engng 12, 1059-1079 (1978).
STRAIGHT MUOVABLE
EDGES ARE HINGE0 , CURVE0 EDGES
5. H. C. Huang. Defect-free shell elements. Ph.D. thesis, C/Ph/89/86. University College of Swansea. 6. H. C. Huang, Further discussion of the new Lagrangian and Serendipity degenerated shell elements. In Proc. Conf. ICCM86-Tokyo. 7. H. C. Huang and E. Hinton, A nine node Lagrangian Mindlin plate element with enhanced shear interpolation. Eng. Comput. 1, 369-379 (1984). a. K. J. Bathe and E. N. Dvorkin, A four-node plate bending element based on Mindlin/Reissner plate theory and a mixed interpolation. Inr. J. numer. Meth. Engng 21, 367-383 (1985). 9. H. C. Huang and E. Hinton, A new nine node degenerated shell element with enhanced membrane and shear intemolation. Inr. J. nutner. Merh. Engng 22, 73-92
AND
ARE
FREE 2.5
20
1.5
f g
1.0
0.5
(198i). 0.0
L/ 0.0
5.0
10.0
CENTRE OlSPLACEbIENT lmml
Fig. 9. Hinged cylindrical shell under concentrated
load.
IO. H. C. Huang and E. Hinton, Elasto-plastic and geometrically nonlinear anlaysis of plates and shells using a new nine node element. Proceedings of ‘Symposium of Finite Element Methods for Nonlinear Problems’, I. 3-1-3-15, Trondheim, Norway (1985).
Implementation
of assumed strain degenerated shell elements
11. L. S. D. Morley and A. J. Morris. Conflict between finite elements and shell theory. Royal Aircraft Establishment Report, London (1978). 12. S. Way, Uniformly loaded, clamped, rectangular plates with large deformation. Proceedings 5th Int. Congress of Applied Mechanics. Cambridge, Massachusetts (1938).
155
13. T. J. R. Hughes and W. K. Liu, Nonlinear finite element analysis of shells: part I. three-dimensional shells. Compur. Merhs appl. Mech. Engng 26, 331-362 (1981). 14. G. Honigmoe, Finite element instability analysis of free-form shells. Report 77-2, Dikision of Structure Mechanics. Norwegian Institute ofTechnology. University of Trondheim, Norway (May 1977).
OF ASSUMED STRAIN SHELL ELEMENTS?
DEGENERATED
Hou CHENG HUANC$ Department of Mechanics, Tianjin University, China (Rereit$ed 23 firne 1986) Abstract-The implementation of the new degenerated shell elements with assumed transverse shear and membrane strains is described in this paper. An explanation of locking ~haviour and also the location of sampling points For the assumed strain fields is given in the present work. In the fo~u~ation of the new elements, assumed transverse shear strains in the naturat coordinate system are used to overcome the shear locking problem and assumed membrane strains in the orthogonal curvilinear coordinate system are applied to avoid membrane locking behaviour. Both tinear and nonii~e~r problems are considered~
1. INTRODUCTION originaf degenerated shell element performs reasonably well for moderately thick sheil situations [I]. However, for thin sheifs when full integration is used to evaluate the stiffness matrix, overstiff solutions are often produced due to shear or membrane locking. Attempts have been made to correct this behaviour by use of reduced or selective integration techniques[2-4]. Such schemes are not always successful in overcoming locking behaviour and the resulting solutions may still be overstiff for problems with highly constrained boundaries especially when coarse meshes are used. Furthermore, for problems with lightly constrained boundaries, element mechanisms or spurious zero energy modes may form. These mechanisms can spread from element to element causing either rank deficiency in the overall stiffness matrix and con~quently no solution, or the even more dangerous situation in which the solution obtained may be polluted by a near-m~hanism. Many research workers have attempted to improve the performance of Mindlin plate elements and degenerated shell elements [5]. From this work it appears that the ideal Mindlin-type element should: (1) converge, (2) not lock, (3) contain no mechanisms, (4) be capable of providing accurate displacements and stress resultants, (5) be insensitive to element distortions, (6) be invariant to the direction of coordinate system, (7) not be based on nume~cally adjusted factors, and (8) be easy to implement and The
Use.
In the formufation of the new degenerated she8 elements, assumed transverse shear strains in the natural coordinate system are used to overcome the shear locking problem and assumed membrane t’fhis is an extended version of the paper presented at the Conference of ICCMES-TOKYO, Japan, May 1986. $ Lecturer, Ph.D.. presently a research scholar at University College of Swansea, U.K. 147
strains in the orthogonal curvilinear coordinate system, are applied to avoid membrane locking behaviour. It has been shown that the new elements are free from any major defect [5].
According to ~indlin-ty~ theories, the rotations are independent of the displacements in both plates and shells. The transverse shear effects should gradually diminish as the plate/shell thickness becomes extremely thin. If we could solve the governing differential equations directly, shear Iocking behaviour should not appear. However, the shear strains obtained using a nurne~~l method, such as the finite element method, may be inaccurate locally, though average shear strains over a particular region are reasonable. Therefore, ii is possible to find some points at which values of the shear strains may represent the average distribution of the shear strain fields at the particular region. As we know, when the finite element method is used in conjunction with Mindlin-type formulations, the shear strain energy is always included and must be positive. Consequently, as the thickness of the plate and shell becomes extremely thin, the shear strain energy can be magnified unreasonably even though the average value of the shear strains over the area tends to zero 161. Consider for instance Fig. I, which shows a rectangular b&ear element subjected to pure bending moments at opposite edges. Apparently, there exists a ‘parasitic’ shear deformation, whereas in reality there should be no shear deformation. However, it is noticed that
Hou CHESGHCASG
148
Fig. 1. Rectangular bilinear element subjected to pure bending moments. whereas the element shear strain energy Gh
where yer, Y,,~are obtained from finite element soiutions and yf& Y; are the transverse shear strains obtained from exact anaiyticai soiutions. For thin plate and shell cases, we have
(2)
JA
y;: dx dy >>0.
It is also noted the values of the shear strain at the line x = 0 are correct. This explains the reason for the success of the use of reduced integration (1 x 1) for the shear strain energy [Z). According to Mindiin-type theory for general thin plate and shell problems, there is always some shear deformation at each point, even for the higher order Lagrangian elements. Therefore, when a full integration scheme is employed, the shear strain energy contribution will dominate the Total Potential Energy and the shear locking problem cannot be avoided. In the following, we take a g-node element as an example for the Lagrangian elements and an 8-node as an example for the Serendipity elements. 3. ASSUMED TRANSVERSE SHEAR STRAIN FIELD
An artificial method for the eitmination of shear locking is to interpolate new shear strain fields from the strain values at the sampling points which are appropriately located in individual etements. For the 9-node Lagrangian element the polynomial terms for the two contributions to each transverse shear strain evaluated from associated displacement fields do not match (51.The substitute shear strains yet and 7,’ should be at least taken as follows:
Yct+Y;‘O.
P,: -) r; + 0.
(5)
When a Mindiin-type hypothesis is applied, then equation (5) cannot be satisfied at ail points. However, equation (5) could be satisfied in an ‘average’ sense as follows:
J~,~?,,dtdl=J~~~~dtdn
-to.
(6)
In the finite element method, AA could be a particular area in an individual element. Now, we consider the special case of thin plates and shells. For the g-node Lagrangian element, the substitute shear strain field has been taken as a linear function of { for ycr and a linear function of 9 for Y,,~. It is assumed that
II J J -1
=0, J’‘r
ycr dc = 0,
0
-I
(7)
0
dn = 0.
y,,xdq==O,
(8)
0
Now, we try to find the value of ( at which the shear strain ycr is equal to zero and therefore we have to solve the following equation:
Y<,= b, t b,< + b,q + b‘& + b& + b&t$ F,, = cr + cr5*+ c,q + c&
+ c&2+
c&l.
(9) (3)
If the substitute strain fields are expressed as poiy nomials of a lower degree than those in equation (3), then zero energy modes will appear. This has happened in the element with seiective integration. It is noted that ‘it, is linear in 5 and quadratic in 9 and T,t is linear in 4 and quadratic in {.
which gives At2+Bt‘-t-C=O,
(10)
where A = H, (B. - 28, + @:,)I2f Hz@:4- 26, + B&2 f HA% - 25, f 4w B = H,[w, f w3 - 2w, f (4, - f&,)/2]
4. LOCATION OF THE SAMPLING POINTS FOR THE SHEAR STRAINS
In general, for plate and shell situations, following relations are expected: ‘J:, + $V
Y4t-+ r,“,
f H&q + wg - 2% -t (e,, - 4J/2] the
+ HJw, + w7 - 2w, + (k, - 4,),‘2] C=H,[(WJ-
(4)
-5 li,Kw,
d/2 -
+
w7)/2
421
f
+
461
ff*[(w,
-
d/2
f
e,1 (11)
of assumed strain degenerated shell elements
Implementation
a
A’I
a
(1.21 I 0 ;
_
b
0
/
-m
IO (2.1) I /
(3,21
0 (1,11
t ’ (**‘)
J
-
-
/ r---g (2.21
1,,i3.21
” !1,21
t 1
I
I
i
c”b
I
(1.1)
0 (3,ll --_-_+----_)c I*,*)
149
(3.11
;
J?r
‘655
Fig. 2. Interpolation points (i, j) for the 9-node Lagrangian element.
and H, = ~,(l + ~i)/2 + (1 - tl*)(l - tl!),
(12)
in which & and vi are the natural coordinates of the element nodes. However, from equation (7) it can be shown that B = 0, A = -3C and therefore equation (10) becomes t*= l/3. So, under the assumption of equation (7) it can be concluded that if 5 = fa (where a = 3-l’*, which are the two Gauss points in the { direction) then ycl = 0 for very thin plates and shells. Similarly, it can be shown that under the assumption of equation (8) if q = fa (which are the two Gauss points in the 7 direction) then yV = 0 for very thin plates and shells. As mentioned in equation (3), six sampling points are required for the representation of the substitute shear strains. For ytr the sampling points should be taken at the lines 5 = a and < = --(I on which three points must be taken. One of the three is q = 0. In order to guarantee the continuity of shear strains between two adjacent elements, the other two points should be 7 = 1 and q = - 1 at the lines 5 = a and 5 = --a. The six sampling points for ya are the three points < = 1, 5 = 0 and C = - 1 at the two lines n = a and n = -a (see Fig. 2). Therefore, we have
and the terms y$ and ;lyY are the transverse shear strains evaluated from the displacement field, a = 3-l’* and b = 1. Similarly, for the 8-node Serendipity element to avoid zero energy modes the substitute shear strains y. and yllr should be taken as polynomials of at the least the following degree [A: Fcr = b, + b25 + b,rj + b&j + b,q’ iisr= Cl + c25 + c,rJ + c&j + c55*.
(15)
Consequently, five sampling points are needed to construct the substitute shear strains. Apparently, along the edges of the element, deformations which are identical to those in the g-node element are expected, and therefore four of the sampling points can be chosen as two points q = 1 and n = - 1 at the two lines r = a and 5 = --a for y(,, and two points { = 1 and 5 = - 1 at the two lines q = a and q = --II. From the pure bending test mentioned in Sec. 2 it is found that the transverse shear strains are equal to zero at the centre of the element and therefore the centre point is introduced as the fifth point for both ycr and y,,*(see also [5]). In fact, for the 4-node element, if it is assumed that
I’-I
yet d5 = 0
(16)
~t,rdrl = 0,
(17)
and I I -1 where Pz(z)=l-
0
it can be shown that y. = 0 when C = 0 and yV = 0 when q = 0 for very thin plates. Therefore, the sampling points for the substitute shear strains can be chosen as those given in [8].
2
;
5. ASSUMED MEMBRANE STRAIN FIELD
Q,W=;(l+$
,,,)=;(I-~)
(14)
The membrane and shear stiffness terms are usually of the same order of magnitude. In shells, the stiffness associated with the membrane strain energy, though
150
HOU CHESG HUAHG Fig. 3). To define the local r and s directions at a typical midsurface point P. it is assumed that PS and PR are intersections of the plane tangential to the shell midsurface and the planes parallel to the r,, - to and s, - to planes drawn through point P respectively. Then the bisectors PB of angles R-P-S and r-P-s are made coincident in the tangent plane (see Fig. 4). The basic vectors gi of the orthogonal coordinate system are defined as follows:
Fig. 3. Symmetry between the natural coordinate system and the orthogonal curvilinear coordinate system.
important, should not dominate the totaf stiffness. It is known that the membrane strains can only be separated from the bending strains in the local Cartesian coordinate system (x’, J’, z’) [l, 91 in which the x’-y’ plane is tangent to the shell mid-surface. Therefore, it is expected that the use of an enhanced interpolation of the membrane strains in an orthogonal curvilinear coordinate system (r, s, t) [S] (in which the (r, s, 0) surface coincides with the mid-surface of the shell element) can help to eliminate membrane locking behaviour. It should be noted that symmetry between the natural coordinate system and the orthogonal curvilinear coordinate system must exist to ensure that this special coordinate system is uniquely defined. The origin of the orthogonal curvilinear coordinate system is set at the centre of each element. The axis t is normal to the shell mid-surface and the (r, s, 0) surface coincides with the ({, q, 0) surface at the shell midsurface. In order to guarantee the required symmetry, the bisectors of angles rO-o-s,, and &-o-~O are made coincident at the centre of each element (see
where g, are unit vectors [j]. From an assumption given in Sec. 4, it is deduced that the values of the membrane strains at lines q = a, -a (where a = 3-‘!?) give the q=-aandt=a,<= best representation of the membrane strain field in each element. Therefore, in the orthogonal curvilinear coordinate system, the enhanced interpolation for e,,,_,, and ie:L,, is the same as that for the y:, and the enhanced inte~oiation for e,_, and iegL,s is the same as that for y,,, but they are expressed in different coordinate systems. The details of the modification of the shell strains are given in [5]. However, the evaluation of the shell element stiffness is conducted in the local coordinate system (x’, y’, z’) in which the x’-J’ plane is tangential to the surface z = constant, as mentioned in [l]. For a shell element of uniform thickness the r, s, I and x’, y’, z’ coordinate systems may be taken as identical systems. All formulations are the same as those mentioned in [I] except for the definition of the local coordinate system. For a general shell element with the assumptions adopted in the present work, a transformation between the orthogonal curvilinear coordinate system and the local coordinate system should be
Fig. 4. ~e~nition of the orthogonal curvilinear coordinate system.
Implementation of assumed strain degenerated shell elements established. The tensor ~ansfo~ation foIIows: ax’ a.d ad
matrix #is as
z-z-z-
B=
dy'
ay'
ay'
drasdt azr at!
af
the stress symbols denote that the stresses at time t are referred to the configuration at time 0. Incremental second Piola-Kirchhoff stresses, ,,A$,, can be expressed as
(19)
(23)
$Si, = &k, OAEW
,,Aq,, are nonlinear incremental strains which are neglected in the *AE, expression and can be described as
arasat
6. GEOMETRICALLY
I51
NONLINEAl RI
(24)
In the geomet~~ally nonlinear analysis both Total Lagrangian (T.L.) and Updated Lagrangian (U.L.) formulations have been considered. In the T.L. formulation, all current static and kinematic variables are referred to the initial configuration before the deformation (say at time 0), whereas in the U.L. formulation, all the current variables are referred to the current configuration for each new load increment (say at time t). Both the T.L. and U.L. formulations include all kinematic nonlinear effects due to large displacements and rotations, but large strain behaviour is not taken into account in the present work f IO]. The local Cartesian coordinate system is used for both geometrically nonlinear formulations. The Green strain tensor can be expressed as follows
In equation (21) ‘+dW is the external incremental virtual work at the current load increment. According to the new element formulation, the shear strain should be interpolated from the values at certain sampling points as shown in equations (13) and (14). Therefore, it can be written that
(25)
p=lq=l
where oAEql = oAcF + oAeT oAEQ = ,A$ + oAe:r.
In the new element formulation the total shear strains have to be interpolated at every iteration in the Total Lagrangian formulation. It appears that the use of the Updated Lagrangian formulation is more appropriate with the interpolation concept used for the new element. 6.1 Total Lagra~gia~ ~or~ul~tio~ (T, L.) In the Total Lagrangian formulation, all variables are referred to the initial configuration at time 0 and the governing equation with incremental terms is
(26)
For the g-node element, p = 1, 2 and q = 1, 3. For the membrane strains, the interpolation in the orthogonal curvilinear coordinate system is only employed to the linear parts in equation (22) as shown in the linear case. ,AE$ and ,AE$ can be obtained from oAEij given by equation (22) by tensor transfo~ation. As the virtual displacements 6Au may be arbitrary, (21) can be rewritten in finite element form as
(I
&BTDoBodv+~oC?“]oCOdL.)Ad
=“SW-
-_ ‘+f -
s
;ISjjS oAEjj Ode, (2 1)
where &E,, is the incremental Green strain tensor which can be divided into two parts so that &E,
= ,AE{~+ ,Ae,
OAsCij = f (oA\u,,+ 0AUj.i)
,B”u dr,
where the incremental strain-displacement matrix ,,B = ,B + oBr,and the matrices oB and o& are associated with oAZand oA8 respectively. Therefore, the ,B matrix is different from that obtained in small deflection analyses. ‘+f are current applied loads and Ad is the incremental nodal displacement vector. Furthermore, it can be shown that KrAd = r,
deij = f(0~k.i o&j + 0Aut.i 0~)
(22) where KT is the tangential
and iSij are total second Piola-Kirchhoff stresses, the superscript 0 and subscript t given in front of
(27)
I
residual forces and K,=K,,-f-K,+K,
stiffness matrix,
(28) r is
152
Hoti CHENGHUANG
in which id=
which are referred to the configuration at time t. The Cauchy stresses are represented as rij. The incremental Green strains, ,A&, can be simply written as
oBTD,,g Odv
and *G is defined as follows oGijAcfi = oAuiQX,for
i = I, 3
oG,,Adi = &ui,,,
for
i = 1, 3
and
k = i+ 3
oGkiAdj = ,Aq,,
for
i = 1, 3
and
k = i + 6,
where ,Ar, and ,Aqi, are the linear and quadratic parts of ,AE,. As the quadratic strain expressions can be neglected in the first term of equation (34), the interpolations for the shear and membrane strains are the same as for the linear case. Equation (34) can be written in finite element form as
+
(301 where the maximum size of sum index j depends on the number of degrees of freedom in each element. For instance, for a 9-node shell eiement j is equal to 45. The stress vector u is
la1 = adds
;;;
:~I.
(32)
j,GT[rl,G'dc)Ad
,BTDIB’dv
= r+f-
,BTr ‘dv
where ,K is associated with ,Acij in which the shear and membrane strains have been substituted according to the formulation of the new shell elements. Therefore,
08)
K,Ad = I,
where S,, = 0 and
[ 1 1 0
I =o
0
10.
0
0
(37)
where KT = K0 + K, and (33)
1
,BT D J ‘do
K, = s
The second parts of shear strain in equation (26) can be neglected and the interpolation of shear strain is only required for the first parts of the shear strain as in the linear case.
,GT [r] ,G ‘do
s-
(39)
in which $G has the same form as in (30) and
6.2 Updated Lagrangian formulation (U.L.) Since dEcE,, and oAE,,, have to be interpolated at iteration, it appears that the use of the Updated Lagrangian formulation is more appropriate for use with the interpolation concept used for the new element. An incrementai form of the principle of virtual work based on the Updated Lagrangian formulation is as follows: every
rCi,k,,AEk, 6 ,AEij ‘do c I
where rj) = 0. 7.
‘~ij6 ,Aqij ‘dV s
7.1
= “6 W -
‘qj 6 ,AE, ‘du
(34)
I. where ‘+6 W is the external incremental virtual work at the current load increment and second PiolaKirchhofT stresses, &IS,, are AS, = rCijk/M’k/ 9
(35)
NUMERICAL
APPLICATIONS
Numerical example for linear onal_vsis
Numerical tests for element characteristics have been conducted in [S, 91. A more stringent example is the pinched hemisphe~cal she11subjected to concentrated forces shown in Fig. 5. TWOequal and opposite concentrated loads are applied in the x and y directions. By taking advantage of the double symmetry which exists in the problem, only a quadrant of the shell is considered with meshes of 1 x 1, 2 x 2, 4 x 4
Implementation
153
of assumed strain degenerated shell elements
2
QUAD 9" OUAO9S
1
--*-
Number of nodes per side
Fig. 6. Convergence of deflections at the loading point.
7;2 Numerical examples for geometrically nonlinear analysis Fig. 5. Spherical shell: radius = 10.0, thickness = 0.04, E = 6.825 x IO’, u = 0.3. and 8 x 8 elements (or n = 3, 5, 9, 17 nodes along each side). In this problem, QUAD9**, the new g-node Lagrangian element is used. A comparison is provided with the performance of QUAD9S, the g-node Lagrangian element with selective integration for the transverse shear stiffness. The locking of QUAD9S is evident in Fig. 6. For the 8 x 8 mesh, the results are less than 30% of the exact solution [I 11. The new element QUAD9** converges.
For the following numerical applications, a modified Newton-Raphson scheme is adopted and a displacement convergence criterion is employed. The preset tolerance is 0.001. 7.2. I Long cantilever. Five QUAD9** elements are used to idealize a cantilever subjected to an end moment as shown in Fig. 7. This problem is chosen to examine the behaviour of the QUAD9** when used for a geometrically nonlinear problem with large displacements and rotations. Here, an Updated Lagrangian formulation is adopted. In the analytical solution to this problem the slope is prescribed at
EI/L=l.O h/L =O.Ol v =o.o
a -1
-2 -3 -t -5 -6_ -1. -8 -2
I
-1
0
I
2
3
t
5
6
7
a
Fig. 7. Elastic cantilever beam, deformed shape.
9
18
Hoti CHENGHUASG
154
0.8 FOUR HETEROSIS .....O,,... FOUR WA0 9 ** -*_--_ WAY LINEAR SOLUTlON-
06
LO
30
20
1.0
00
CENTRE OlSPLACEflENTlmml
Fig. 8. Clamped square plate under uniform load.
the free end as 0 = ML/EI and the neutral axis is deformed into a circle with radius R = EI/M. L is the beam length and M is the end moment and El is the flexural rigidity. For a cantilever of length L = 10, the analytical solution indicates that the radius of the final circle is 3.183, whereas the radius obtained using the new element, QUAD9 **, is 3.185. In Fig. 7, it is shown that the error is reasonably small when the load approaches the value of 27rEZ/L which is divided into 100 load increments. 7.2.2 Clamped square plate. Four QUAD9** are used to idealize a symmetric quadrant of a clamped square plate subjected to a uniformly distributed load as shown in Fig. 8. The number of the increments
R=ZSLOmm L:ZSLmm ha12.7mm k310275 v.o.3
N/mm'
adopted is 15. Lateral displacements at the plate centre obtained using the new element are compared with the results obtained by the classical RayleighRitz solution presented by Way [12]. The results obtained by both T.L. and U.L. formulations are almost identical. Figure 8 shows that the central displacements obtained in the present study lie between those obtained by the classical solution and those obtained using heterosis elements [ 131. It is worth noting that the results obtained using the Lagrangian nine node element with selective integration, QUAD9S, are almost the same as those obtained using QUAD9**. 7.2.3 Hinged cylindrical shell. The hinged cylindrical shell. shown in Fig. 9, has been studied by Horrigmoe [14]. The straight edges are hinged and the curved edges are free. The shell is subjected to a concentrated central load on the convex side. A 2 x 2 mesh of QUAD9** elements is used to represent a symmetric quadrant of the shell. Here also, identical results are obtained using both T.L. and U.L. formulations with 20 load increments. In Fig. 9 the load-displacement response obtained using QUAD9** elements is compared with that given in Horrigmoe’s study in which a max. load of 2.22 kN was obtained. In the present study the max. load is 2.21 kN compared with a value of 2.20 kN obtained by Hughes [ 131 using four heterosis elements. REFERENCES
I. S. Ahmad, B. H. Irons and 0. C. Zienkiewicz, Analysis
of thick and thin shell structures by curved finite elements. Inr. J. numer. Merh. Engng 2, 419-451 (1970). 2. 0. C. Zienkiewicz, R. L. Taylor and J. M. Too, Reduced integration techniques in general analysis of plates and shells. Inr. J. numer. Meth. Engng 3, 275-290 (1971). 3. E. Onate, E. Hinton and N. Glover, Techniques for improving the performance of Ahmad shell elements. Proc. Int. Conf. on Applied Numerical Modelling, Madrid. Pentech Press (1978). 4. E. D. L. Pugh, E. Hinton and 0. C. Zienkiewicz, A study of quadrilateral plate bending elements with ‘reduced’ integration. Inr. J. numer. Meth. Engng 12, 1059-1079 (1978).
STRAIGHT MUOVABLE
EDGES ARE HINGE0 , CURVE0 EDGES
5. H. C. Huang. Defect-free shell elements. Ph.D. thesis, C/Ph/89/86. University College of Swansea. 6. H. C. Huang, Further discussion of the new Lagrangian and Serendipity degenerated shell elements. In Proc. Conf. ICCM86-Tokyo. 7. H. C. Huang and E. Hinton, A nine node Lagrangian Mindlin plate element with enhanced shear interpolation. Eng. Comput. 1, 369-379 (1984). a. K. J. Bathe and E. N. Dvorkin, A four-node plate bending element based on Mindlin/Reissner plate theory and a mixed interpolation. Inr. J. numer. Meth. Engng 21, 367-383 (1985). 9. H. C. Huang and E. Hinton, A new nine node degenerated shell element with enhanced membrane and shear intemolation. Inr. J. nutner. Merh. Engng 22, 73-92
AND
ARE
FREE 2.5
20
1.5
f g
1.0
0.5
(198i). 0.0
L/ 0.0
5.0
10.0
CENTRE OlSPLACEbIENT lmml
Fig. 9. Hinged cylindrical shell under concentrated
load.
IO. H. C. Huang and E. Hinton, Elasto-plastic and geometrically nonlinear anlaysis of plates and shells using a new nine node element. Proceedings of ‘Symposium of Finite Element Methods for Nonlinear Problems’, I. 3-1-3-15, Trondheim, Norway (1985).
Implementation
of assumed strain degenerated shell elements
11. L. S. D. Morley and A. J. Morris. Conflict between finite elements and shell theory. Royal Aircraft Establishment Report, London (1978). 12. S. Way, Uniformly loaded, clamped, rectangular plates with large deformation. Proceedings 5th Int. Congress of Applied Mechanics. Cambridge, Massachusetts (1938).
155
13. T. J. R. Hughes and W. K. Liu, Nonlinear finite element analysis of shells: part I. three-dimensional shells. Compur. Merhs appl. Mech. Engng 26, 331-362 (1981). 14. G. Honigmoe, Finite element instability analysis of free-form shells. Report 77-2, Dikision of Structure Mechanics. Norwegian Institute ofTechnology. University of Trondheim, Norway (May 1977).