F i n i t e E l e m e n t s i n A n a l y s i s a n d D e s i g n 1 ( 1 9 8 5 ) 3 - 20 North-Holland
A PRO POS ED
STANDARD
S E T O F P R O B L E M S T O T E S T F IN IN I T E E L E M E N T
ACCURACY Richard H. MACNEAI. and Robert L. HARDER The
ac Ne al-S chw end ler Corpor Corporati ation. on. 815 Colorado Colorado Boule Boulevard. vard. Los Angeles. CA 90041, U.S.A.
Received February 1984 A b s t r a c t . A p r o p o s e d s t a n d a r d s e t o f t e s t p r o b l e m s i s d e s c r i b e d a n d a p p l i e d t o r e p r e s e n t a t i v e q u a d r i l at e r a l p l a t e a n d s o l i d b r i c k f i n i t e e l e m e n t s . T h e p r o b l e m s e t c o n t a i n s p a t c h t e s t s a n d b e a m , p l a t e , a n d s h e l l p r o b le m s , s o m e of which have become de facto standards for comparing the accuracy of finite elements. Although few in number, the tests are able to display most of the parameters which affect finite elemen t accuracy.
Introduction
The intended purpose of the proposed problem set is to help users and developers of finite element programs to a~certain the accuracy of particular finite elemen1~ in various applications. It is not intended that the problems be used as benchmarks for cost co mparisons since the p r o b l e m s a r e , in i n g e n e r a l , t o o s m a l l t o b e m e a n i n g f u l f o r t h is is p u r p o s e . Nothing is as important to the success of a finite element analysis as the a ccuracy of the elements. Indeed, in a linear static analysis, the finite elements embody all of the discretizing assumptions: the rest of the calculations are exact except for their lack of precisio n. Thus the accuracy of the finite elements should be a matter of primary concern to thos e who perform f i n itit e e le le m e n t a n a l y s e s a n d t o t h o s e w h o a r e r e s p o n s i b l e f o r c o n c lu l u s i o n s d e r iv iv e d t h e r e f r o m . E v e r y n e w f i n i te t e e l e m e n t i s te t e s t ed e d w i t h o n e o r m o r e s m a l l p r o b l e m s a n d t h e r e s u l ts ts o f t h e s e tests are generally published in the open literature or in FEM program docum entation. Such test results are all that the user has to help him evaluate the elements prior to actual use. They have invariably proven to be woefully inadequate for this purpose: (a) becau se they test an insufficient number of conditions. (b) because few, if any, bad results are reported (not by design but because the developer fixed only the bugs he discovered), and (c) because they usually cannot be compared with results for other elements, particularly with those in other p r o g r a m s . T h e s e s a m e d e f e c ts ts w o u l d n o t b e p r e s e n t i n a ca c a r e f u l ly ly d e s i g n e d s et e t o f s t a n d a r d t es es t p r o b l e m s a p p l i e d t o m a n y d i f f e re re n t f i n i te t e e l e m e n t s a n d w i d e l y c i rc r c u l a te te d . The need for verifying finite element accuracy by independent testing and f or compariag finite element results is becoming more widely recognized. A recent effort to evaluate the plate e l e m e n t s i n c o m m e r c i a l p r o g r a m s [ 1 ,2 ,2 ] r e v e a l e d r e s u lt lt s r a n g i n g f r o m e x c e llll e n t to to e x t r e m e l y p o o r a n d m i s l e a d in in g . G o v e r n m e n t a l c o n c e r n f o r t h e a c c u r a cy c y o f f i n i te te e l e m e n t a n a l y s is i s is i s e v id id e n c e d b y t h e N u c l e a r R e g u l a t o r y C o m m i s s i o n ' s re re q u , e m e n t f o r s t r u c t u r a l a n a l y s is is c o m p u t e r p r o g r a validation, and abroad by the recent formation in the United Kingdom of a Na tional Agency for Finite Element Methods and Standards (NAFEMS). While the former agency relies on verification problems supplied by the developer, the latter shows promise of doin g their own testing. The authors can confirm from personal experience that the design of an adequate set of problems for finite element testing requires careful planning. Too often initial te sting is done 0 1 6 8 - 8 7 4 X / 8 5 / $ 3 . 3 0 © 1 98 98 5 , E ls ls e v i e r S c i e n c e P u b l i s h e r s B .V .V . ( N o r t h - H o l l a n d
4
R.H. MacNeal,
. L H a r d e r / P r o p o s e d s t a n d a r d s e t o f p r o b l e m s t o t e st F E a c c ur a c
w i t h o n e o r t w o r e l a t i v e ly e a s y p r o b l e m s w i t h k n o w n a n s w e r s . F o r e x a m p l e , t h e o n l y p r o b l e m u s ed i n d e v e l op m e n t t e s ti n g o f t h e N A S T R A N ® p l at e b e n d i n g e l em e n t s T R P L T a n d Q D P L T w a s t h e l a te r a l l o a d i n g o f a r e c t a n g u l a r p l a t e [ 3 ]. T h e s u b s e q u e n t u s e o f Q D P L T a n d i t membrane counterpart QDMEM in the solution of verification problems did not reveal any weaknesses in these elements because the verification problems were chosen to gi ve excellent c o m p a r i s o n s w i t h t h e o r y r a t h e r t h a n t o a c t a s c r it i c al e l e m e n t t e s t s . W e w i ll s h o w i n t h i s p a p e r t h a t t h e se tw o e l e m e n t s ( w h ic h w h e n c o m b i n e d f o rm t h e Q U D 2 e l e m e n t) in f a c t h a v e m a n y w e a k n es s es . T h i s e x a m p l e is n o t a n i s o la t e d on e . E v er y e l e m e n t i n M S C / N S T N has bee revised in response to difficulties encountered in the field. In that process, we have gradually built up a library of element test problems which clearly demonstrate frequen tly encountered element failure modes. The problems to be described in this paper represent a sub stantial part o f t h a t l i b r a ry . The most important symptoms of accuracy failure in modern finite elements are spurious mechanisms, also known as rank deficiencies, and a phenomenon known as lock ing in which e x c es s iv e s t i ff n e s s is e x h i b it e d f o r p a r t i c u l a r l o a d i n g s a n / o r i r r e g u l a r s h a p e s . M o s t e l e m e n t s display one or the other of these symptoms, but not usually both. An important stat e-of-the-art problem is the design of elements which are free from spurious mechanisms and lock ing in situati ons. Element ary ,de fects of element design, such as viol ation of rigi d body p roperty and noninvariance to node r~umbering, are less frequently encountered nowadays, but are de vastating when they occur. The design of a comprehensive set of element test problems should, of c ourse, take into account the parameters which affect accuracy. These parameters can be clas sified under the headings of loading, element geometry, problem geometry, and material properties. With regard to loading, the problem set should, as a whole, provide significant loading for
Aspect Ratio
t
b
w
Skew
Taper (2 Directions)
Warp
Fig. 1. Ty pes of geometricdistortion from a square plate
a/b
R.H. MacN eai,
. L H a r d e r / P r o p o s e d s t a n d a rd s e t o f p r o b l e m s t o t e st F E a c cu r ac y
each of the types of deformation which the elements can exhibit. For exampl e, a three-noded shell element should be subjected to extension, in-plaae shear, and out-of-plane be nding. For a four-noded shell element, add in-plane bending and twist. For an eight-noded elem ent, add the o t i o n o f t h e e d g e n o d e s r e l a t iv e to t h e c o r n e r s . T h e l a t t e r a r e le s s i m p o r t a n t , b u t s h o u l d n o t b e neglecte d entir ely. E a c h e le m e n t h a s a s t a n d a r d s h a p e w h i c h m a y b e t h e o n l y s h a p e t h a t th e d e v e l o p e r h a s tested. In the case of a quadrilateral, the standard shape is a square; in the case of a hexahedron, the standard shape is a cube; and in the case of a triangle, th e standard shape is u s u a l l y a n i s o c e le s r i g h t t r i a n g l e. C a r e s h o u l d b e t a k e n t o t e s t n o n s t a n d a r d s h a p e s . F i g . 1 s h o w s the four basic modes of distortion of a square, each of which should be exhibited in the test problems and tested with several kinds of loading. Geometric parameters which are not isolated to single elements can also affect dement a c c u r a c y . C u r v a t u r e i s t h e m o s t i m p o r t a n t s u c h p a r a m e t e r . I t i s n o t s u f f i c ie n t t o t e s t o n l y s i n g l e c u r v a t u r e s i n c e s o m e e l e m e n t s w h i c h b e h a v e w e ll f o r s i ng l e c u r v a t u r e b e h a v e p o o r l y f o r d o a b l e c u r v a t u r e . T h e s l e n d e r n e s s r a t io a n d t h e m a n n e r o f s u p p o r t o f a s t r u c t u r e a f f ec t t h e c o n d i t io n ing of the stiffness matrix and therefore can be used to check element failures related to precision. Poisson's ratio has a strong effect on element accuracy as its value app roaches 0.5. Such values should be included in the problem set if the use of nearly incompr essible materials is contemplated. Plasticity affects element accuracy in much the same wa y as incompressible material. Plasticity and all other nonlinear effects are outside the scope o f the paper. Anisotropic material properties also have a significant effect on element accuracy whic h will not be examined here T h e t e s t p r o b le m s
The names of the proposed test problems are listed in Table 1, which also indicates the suitability of each problem for testing various types of elements . The geometry, material p r o p e rt i e s, b o u n d a r y c o n d i t io n s , l o a d in g , a n d e l e m e n t m e s h i n g f o r e a c h p r o b l e m a r e d e s c r ib e d in Figs. 2 through 10 in sufficient detail to permit construction o f a finite element model
l
4
F i g . 2 . P a t c h t e s t f o r p l a t e s , a = 0 . 1 2 ; b = 0 . 2 4; t - - 0 . 0 0 1 ; E = 1 . 0 × 10 6 ; p = 0 . 2 5 . B o u n d a r y c o n d i t i o n s : s e e T a b l e 2 . Location of inner nodes: x 1
0.04
0.02
2 3
0.18 0.16
0.03 0.08
4
0.08
0.08
6
R.H . Mac Neal, R.I., Hard er / Proposed standard set of problems to test FE accuracy
Table 1 S u m m a r y o f p r o p o s e d t e s t p r o b le m s Test problem
S u i t a b i li t y o f p r o b l e m f o r e l e m e n t t y p e Beam
Putch tests Straight cantilever beam Curved beam Twisted beam Rectan gular plate S c o r d e l .; s - L o r o o f Spherical shell Thick-walled cylind er
× x x
Me mbra ne plate
Bending plate
× × x
x x x
S h el l
Solid
× x × × × b
a A s h e ll e l e m e n t is d e f i n e d h e r e as a n e l e m e n t t h a t c o m b i n e s m e m b r a n e a n d b e n d i n g p r o p e r t ie s . b Using plane strain option.
F i g . 3 . P a t c h t e s t f o r s o ! i d s. O u t e r d i m e n s i o n s : u n i t c u b e ; E = 1 .0 x 1 0 6 ; p = 0 . 2 5 . B o u n d a r y c o n d i t i o n s : s e e T a b l e 2 . Location of inner nodes:
1 2 3 4 5 6 7 8
x
y
0.249 0.826 0.850 0.273 0.320 0.677 0.788 0.165
0.342 0.288 0.649 0,750 0.186 0.305 0.693 0.745
0.192 0,288 0,263 0,230 0,643 0.683 0.644 0.702
R . H . M a c N e a l , R . L . H a rd e r
P r o p o s e d s t a n d a r d s e t o f p r o b l e m s t o t e s t F E a cc ur ac y
I
t
I
I
",,,
/
",,,
/
/
/
I
/
, ~ 5 /°
I
/•45
",,,
/
,
Fig. 4. Straight cantilever beam. (a) Regular shape elements; (b) Trapezoidal shape elements; (c) Parallelogram sh ape elements. Length -- 6.0; ~.,idth -~ 0.2; depth -- 0.1; E -- 1.0 × 107; p -- 0.30: mesh -- 6 x 1. L oading: unit forces a t free end. ( Note: All elements have equal volume.) c o n s i s t i n g o f b e a m , q u a d r i l a t e r a l p l a t e , s h el l, o r b r i c k e l e m e n t s . A n a p p r o p r i a t e m e s h i n g
for
t r ia n g l e s a n d w e d g e e l e m e n t s c a n b e o b t a i n e d b y s u b d i v i d in g t h e q u a d s a n d b r i ck s . T h e o r e t i c a r e s u l ts f o r th e p r o b l e m s a r e g iv e n in T a b l e s 2 t h r o u g h 5 .
FIXED
Fig. 5. Curved beam. Inner radius-- 4.12: outer radius = 4.32: a rc -~ 90° : thickness = 0 .] ; E = I . 0 x l 0 7 ; p = 0 . 2 5 : m e s h = 6 Loading: unit forces at tip. ×
8
R.H. MacNeal, R.L. Harder
Proposed standard set o f problems to test FE accuracy
FIXED EN
F i g . 6. T w i s t e d b e a m . L e n g t h = 1 2 .0 ; w i d t h
--
1.1; dep th = 0.32; twi st = 90 ° (roo t to tip); E = 29.0 × 106;
, -- 0.22:
mesh -1 2 x 2. Load ing: unit forces at tip.
o c o m p r e h e n s i v e s e t o f f in i t e e le m e n t t e s t p r o b l e m s w o u l d b e c o m p l e t e if i t d i d n o t i n c l u d e patch tests for plate and solid problems. The patch test that we propose for pla tes, shown in F i g . 2 , h a s b e e n u s e d b y R o b i n s o n [ 1 , 2 ] t o t e s t c o m m e r c i a l f i n it e e le m e n t s . N o t e t h a t t h e arbitrarily distorted element shapes are an essential part of the test. On th e other hand, the rectangular exterior shape of the plate makes it easy to provide boundary conditions corresponding to constant membrane strains or constant bending curvatures, independent of ele ment s h a p e . W e h a v e e l e c te d t o u s e d i s p l a c e m e n t b o u n d a r y c o n d i t i o n s ( s e e T a b l e 2 ) b e c a u s e t h e y a r e easier to specify for a variety of elements than the force and raoment b oundary conditions
--T
I• _! sym
.I Fig. 7. Re cta ng ula r plat e, a -- 2.0; b -- 2.0 or 10.0; thickn ess = 0.0001 (p lates ); th icknes s = 0.01 (solids) ; E = 1.7472 x 10T; /4 o f pl a t e) . L o a d i n g : u n i f o r m p r e s s u r e = 0 . 3; b o u n d a r i e s = s i m p l y s u p p o r t e d o r c l a m p e d ; m e s h - N x N ( o n q = 10 -4, or centra l load P = 4 .0 x 10 -4.
YX
F i g. 8 . S c o r d e l i s - L o r o o f , R a d i u s = 2 5. 0 ; len gth --- 50.0; th ick ne ss --- 0.25 ; x 10s; area
in
E .-, 4.32
p = 0.0; loa din g-- 90.0 per - Z
dire ctio n;
curved edges; mesh: area.
unit
ux ~- u= - 0 on
N x N on sh aded
R.H. MacNeai, R. L Harder / Proposed standard set of problems to test FE accuracy
quadrant)
Loading: concentrated forces as shown.
employed by Robinson. If the latter are used, the load distribution rules s~own in Fig. 11 are a p p r o p r i a t e f o r i s o p a r a m e t r i c e l em e n t s . T h e p r i n c i p a l v i r t u e o f a p a t c h t e s t i s t h a t , i f a n e l e m e n t p r o d u c e s c o rr ex :t r e s u lt s f o r t h e t e s t , t h e r e s u l t s f o r a n y p r o b l e m s o l v ed w i t h t h e e l e m e n t w i ll c o n v e r g e t o w a r d t h e c o r r e c t s o l u t i o n a s
Radius
Fig. 10. Thick-wall ed cylin der. I nner radiu s = 3.0; out er r adius = 9.0: thickn ess = 1.0; E -- 1000; , -- 0.49, 0.499. 0.4999: plane strain condition; mesh: 5 × 1 (as shown above). Loading: unit pressure at inner radius.
10
R.H. MacNeal, R.I., Harder / Proposed standard set of problems to test FE accuracy
Table 2 Boundary conditions ~nd theoretical solutions for patch tests (a) Membrane plate patch test Boundary conditions: u = 1 0 - 3 ( x v = 10-3(y Theoretical solution: x= Cy =¥ =1 0- 3;
y/2)
x/2)
x =Oy=1333. ; % y= 400
(b) Bending plate patch test Boundary conditions: w = 10 -3( x 2 xy y2)/2 x = att'/~y = 10-3(y + x / 2 ) o,, =
-aw/ax = 10-3(-
x
-
y/ 2)
Theoretical solution: Bending moments per unit length: =1.111 )<10-7; mxy =10 -7 mx= Surface stresses: o~ = oy = +0.667; ~'xy= +0.200 (c) Solid patch test Boundary conditions: u = 10-3(2x + y + z)/2 v = 10-3(x +2 y z)/2 w =10-3(x + y + 2z)/2 Theoretical solution: x = ¢y = ~: = Yxy= Yyz Vzx = 1 0 - 3 x = Oy =.o: = 2000; ~-,,y= ~-y:= ~':x= 400
Table 3 Theoretical solutions for beam problems Tip load direction
Extension In-plane shear Out-of-plane shear Twist
Displacement in direction of load Straight beam
Curved beam
Twisted beam
3.0 × 100.108 0.432 0.03208
0.08734 0.5022
0.005424 0.001754
Table 4 Theoretical solutions for rectangular plate Boundary supports
Simple Simple Cla mped Cla mped
Aspect ratio
Displacement at center of plate b Uniform pressure
Concent rated force
1.0 5.0 1.0 5.0
4.062 12.97 1.26 2.56
11.60 16.96 5.60 7.23
See [8] for the method. b Values shown are for shell elements; multipl y by 10 -6 for solid element s because the thickne ss used for solid elements is 100 times greater. Erratic results were obtained with the solid elements at the lesser thickness.
R . H . M a c N e a l , R . L . H a r d e r / P r o p o s e d s t a n d ar d s e t o f p r o b l e m s t o te s t F E a c c u ra c y
11
Table 5 Theoretical solutions for shell and cylinder problems (a S c o r d e l i s - L o r o o f The value for the midside vertical displacement quoted in [5] is 0.3086. an y finite elements converse to a slightly smaller value. We hav e used the value 0.3024 for normalization of our results. (b Spherical shell Reference [7] computes a theoretical lower bound o f 0.0924 for the displacement under lo ad in the case w here the hole at the center is not present. W e hav e used the value of 0.0940 for normalization of ou r results. (c Th ick. walled cylinder Form ula for radial displacement: u =
0 +
R2 _ R!
[R
/r
where p = pressure; R I -- inner radius; R 2 = outer radius. poisson's ratio
Radial displacer, ent at r = R
0.49 0.499 0.4999
5.0399 × 10 5.0602 x 10 5.0623 x 10 -
t h e e l e m e n t s a r e s u b d i v i d e d . T h e r e a s o n , o f c o u rs e , i s t h a t t h e s t r e s s w i t h in e a c h e l e m e n t t e n d s to a uniform value in the limit. Many authorities, including Irons [4], feel that an element that d o e s n o t p a s s t h e p a t c h t e s t s h o u l d n o t b e t r u s te d . O n t h e o t h e r h a n d , p a s s i n g t h e p a t c h t e s d o e s n o t g u a r a n t e e s a t i s fa c t i o n s i nc e t h e r a t e o f c o n v e r g e n c e m a y b e t o o s l ow f o r p ra c t i c a l u s e . T h e p r o p o s e d p a t c h t e s t f o r so l i d e l e m e n t s ( s e e F i g . 3) i s s e e n t o b e a n e x t e n s i o n o f Robinson's patch test to three dimensions. Displacement boundar y conditions for the test are also given in Table 2. Th e straight cantilever beam (see Fig. 4) is a frequen tly used test proble m w hich can be a p p l i e d t o b e a m , p l a t e , a n d s o l id e l e m e n t s , t I t s v i r t u e s a r e i t s s i m p l i c i t y a n d t h e f a c t t h a t a l l o f t h e p r i n c i p a l e le m e n t d e f o r m a t i o n m o d e s d e s c r i b e d e a r l i er ( c o n s t a n t a n d l in e a r ly v a r y i n g s t r a in s a n d c u r v a t u r e s ) c a n b e e v o k e d b y l o a d s a p p l i e d t o t h e f re e e n d . W e h a v e ~ a d e d i r re g u l a r e l e m e n t s h a p e s p r i n c i p a l l y t o t e s t t h e c o m b i n a t i o n o f s u c h s h a p e s w i t h l i n e a rl y v a r y in g s t ra i n s . I n t h e c u r v e d c a n t i l e v e r b e a m ( s ee F i g . 5 ), c o m b i n a t i o n s o f t h e p r i n c ip a l d e f o r m a t i o n m o d e s a r e e v o k e d b y a s i n g le i n - p l a n e o r o u t - o f - p l a n e s h e a r l o a d a t t h e t i p . N o t e a l s o t h a t t h e e l e m e n t s h a p e i s n o t q u i t e r e c t a n g u l a r , w h i c h w i l l t e s t t h e e f fe c t o f s l ig h t i r r e g u l a r it y . T h e o r e t i c a l t i p displacement results for all of the beam problems are shown in Ta ble 3. T h e t w i s te d b e a m e l e m e n t is t h e o n l y o n e i n o u r p r o p o s e d s e t t h a t t e s t s th e e ff e c t o f w a r p o n p l a t e e l e m e n t s . T h e w a r p o f e a c h e l e m e n t i s o n l y 7 . 5 ° , b u t e v e n t h i s s m a l l a m o u n t w i ll p r o d u c a su rprising result. A s n o t e d e a r l i e r , t h e l a t e r a l l y l o a d e d r e c t a n g u l a r p l a t e ( s e e F i g . 7) is a s u r v i v o r f r o m t h e o r i g in a l t e s ti n g o f N A S T N ® p l a t e e le m e n t s . I t h a s b e c o m e a d e f a c t o s t a n d a r d t e s t a n d h a s b e e n s e e n f r e q u e n t l y in t h e t e c h n i c a l l it e r a t u r e . T h e o r e t i c a l r e s u l t s fo r l a te r a l d i s p l a c e m e n t a t the center are provided in Table 4 for all of the combinations of boundary supports, aspect r a t i o , a n d l o a d i n g c o n d i t i o n s l i s t e d in F i g . 7 . T h i s i s th e f i r st p r o b l e m i n w h i c h c o n v e r g e n c e w i t h decreasing mesh spacing will he studied. Note that the proper wa y to apply free end load depends on element type (see Fig. 11).
12
R.H. MacNeal, R .L Harder / Pr~posedstandard set of problems to test FE accuracy
The Scordelis-Lo roof [5] (see Fig. 8) also has achieved the status of a de facto standard te st p r o b l e m a n d h a s f o u n d i ts w a y i n t o t e x tb o o k s [ 6 ] . T h e t e s t r e su l t m o s t f r e q u e n t l y d i s p l a y e d i s the vertical displacement at the midpoint of the free edge. The theoretical value for this result is 0 .3 0 8 6 , b u t m o s t e l e m e n t s c o n v e r g e t o a s l i gh t ly lo w e r v a l u e . M e m b r a n e a n d b e n d i n g d e f o r m a tions both contribute significantly to it. The Scordelis-Lo roof is the only singly-c urved shell problem in the proposed standard problem set. T h e s p h e r i c a l sh e l l s h o w n i n F i g . 9 is o u r p r o p o s e d d o u b l y - c u r v e d s h e l l p r o b l e m . N o t e t h a t the equator is a free edge so that the problem represents a hemisphere with fo ur point loads a l t e r n a t in g i n s i g n a t 9 0 ° i n t e rv a l s o n t h e e q u a t o r . T h e h o l e a t t h e t o p h a s b e e n i n t r o d u c e d t o avoid the us e of triangles ne:~r the axis of revolution. Conv ergen ce can be studied by vary ing m e s h s i ze . B o t h m e m b r a n e a n d b e n d i n g s t r a i n s c o n t r i b u t e s i g n if i ca n t ly t o t h e r a d i a l d i s p l a ce ment at the load point. A theoretical value of the displacement under load has been computed for a slightly different co nfigu ration 7] in wh ich the h ole at the axis is close(;.
1
,
12
1
4
I
4
1-'2
HEXA (8
1"2
HEX20 (R) .
z
C
4
12
12
QUAD8 C
QUAD2 r QUAD4
O
h~
0 I
0 I
~
~ I
6-~
I
-2--~
2-
HEX20 6 -E "
HE XA (8
HEX20 (R)
QUAD8 C
1
0
QUAD2 r QUAD 0
1
0
I
Fig. 11. Consistent load distribution on beam end. (a) Tension or shear oad. (b) Mo me nt load.
13
R . H . M a c N e a l , R . L H a r d e r / P r o p o s e d s t a n d a rd s et o f p r o b l e m s t o t e s t F E a c c u r ac y
Table 6 Summ ary of element properties Shape
Name
Connections Num ber of points
2
~
OUAP4
4
QUADS
~
)
0
HEX20(R) ~
a Note:
Description
Reference
7", R
Com posed of two pairs of overlapping triangles. Each triangle has constant mem brane strain and lateral deflection described by a n incomplete (nine-term) cubic polynomial
[9
r,
Isoparam etric with selective reduced order integration. Transve rse shear uses string-net approximation and augmented shear flex ibilit
ll0.111
7", R
lsoparam etric with selective reduced order integration
D21
Components
~
8
lsoparam etric with selective reduced order integration and internal strain functions equivalent to b ubb le functions
D3d41
~
0
~?tandard isoparam etric
[15]
Standard isoparametric with reduced order integration.
1161
20
T = all three components o f translation; R -- all three components of rotation.
The section of a thick-walled cylinder shown in Fig. 10 has been chosen to tes t the effect of nearly incompressible material. Note that a plane strain condition is assumed which, along with the radial symmetry, confines the material in all but the radial direct ion and intensifies the numerical d~fficulty caused by near incompressibility.
The elements
T h e a b i l it y o f t h e p r o p o s e d t e s t p r o b l e m s t o d i s c r i m i n a t e g o o d a n d p o o r e l e m e n t a c c u ra c y w i ll b e ex a m i n e d b y t e st in g a s e le c ti o n o f N A S T R A N ® a n d S C / N S T R A N e l em e n t s . T h e s e le c t io n i n c lu d e s t h r e e q u a d r i l a te r a l s h e ll e l e m e n t s ( Q U A D 2 , Q U A L M , a n d Q U A D S ) a n d t h r e e s o l id b r ic k e l e m e n t s ( H E X A ( 8 ) , H E X 2 0 , a n d H E X 2 0 ( R ) ) . T h e p r o p e r t ie s o f e a c h o f t h e si x elements 2 are su mm arized in Table 6. The QUAD!2 elem ent is the original NA STR AN ® quadrilateral plate element. It was conceived in 1967 an d it is still heavily used in N ST RA N ®. It has be en replaced in MSC/NASTRAN by the QUAD4, which is probably the most frequently used S C / N A S T R A N e le m e n t, a n d b y t he e ig h t n o de d Q U A D 8 . 2 One of the authoi!s of Ithis paper, R.H. MacNeal, is also the author of four of the elements (QUA D2, QUA D4, QU AD S, an d HF~I~,A(8))and therefore need not apologize to anyone, except perhaps to the users of these elements. for an y unk ind wol~'ds hat m ;,.y be said he re abo ut them.
14
R . H . M a c N e a l , R . L . H a r d e r / P r o p o s e d s t a n d a r d s e t o f p r o b l e m s t o t e st F E a c c u r a c y
T h e M S C / N S T R A N b ri ck e le m e n t, H E X A , h a s b o t h a n e ig h t- n od e d an d a t w e n ty - n od e d f o rm . O n l y t h e m o r e f r e q u e n t l y u s e d e i g h t - n o d e d f o r m , H E X A ( 8 ) , w i l l b e t e s te d . T h e o t h e r t w o b r ic k e l e m e n t s t o b e te s te d a r e t h e i n fr e q u e n tl y u s e d S C / N S T N HEX 20 and HEX 20(R), w h i c h a r e i n c l u d e d h e r e b e c a u s e " ~h ey a r e r e s p e ct iv e l y i d e n t ic a l t o t h e w e l l k n o w n a n d w i d e l y used standard isoparametric brick element with full and reduced order integration. The QUAD4, QUADS, and HEXA(8) elements are modified occasionally, and som e of the m o d i f i c a t i o n s m a y a f f e c t a cc u r a c y . T h e r e s u l t s to b e p r e s e n t e d w e r e o h : a i n e d w i t h V e r s io n 6 3 o f MSC/NASTRAN and may not agree in all details with results from earlier or later versions.
Test results
T h e t e s t r e s u l t s a r e d e t a il e d i n T a b l e s 7 t h r o u g h 1 5 a n d s u m m a r i z e d i n T a b l e s 1 6 a n d 1 7. T h e test data have been restricted to one data point per case in order to save s pace. The quantity recorded in the patch test results (Table 7) is the maximum stress error. In cases with point loads, the recorded quantity is the displacement under load, which is a direc t measure of the t o t a l s t r a in e n e r g y , i n c a s e s w i t h d i s t r i b u t e d l o a ds , t h e l a r g e s t d i s p l a c e m e n t c o m p o n e n t h a s b e e n recorded. M o s t o f t h e d a t a i n T a b l e s 7 t h r o u g h 1 5 a r e p r e s e n t e d i n n o r m a l i z e d f o r m , i .e ., t h e r e c o rd e d v a lu e is t he r a t io o f t h e v a l u e o b t a i n e d i n a n M S C / N S T N f in i te el e m e n t a n a ly s i s t o t h e t h e o r e t ic a l v al u e . T h e o n e e x c e p t i o n o c c u r s i n t h e p a t c h t e s t w h e r e t h e p e r c e n t e r r o r i s r ec o r d e d . I n T a b l e s 1 6 a n d 1 7, t h e t e s t d a t a a r e r e d u c e d t o a l e t t e r g r a d e b y t h e f o l l o w i n g r u le s : Graae
Rule
A B C D F
27o >i er ro 107o >I er ror > 2% 20¢g 1er ro r > 107O 50~ >I error > err or > 507o
In situations where more than one case contributes to a letter grade, the abso lute errors have been avera ged before assigni ng a letter grade. Tables 16 and 17 reveal a great deal about the elements and about the te sts. It should be noted, first of all, that t he propose d test problems are deceptively difficult and that a fai ling grade in one or more problems should not disqualify the element for appli cations where the combination of parameters that caused the failure do not occur. Note, for ex ample, that both t h e Q U A D 2 a n d t h e Q U A D 8 f a il t h e p a t c h t e st fo r o u t - o f - p la n e l o ad i ng . T h e Q U A D 2 t h e n p r o c e e d s t o f a i l fo u r o f t h e e i g h t o t h e r t e s ts t h a t i n v o l v e o u t - o f - p l a n e l o a d i n g w h i l e t h e Q U A D f a il s n o ne . T h e t e s t sc o r e s o f a l l t h e e l e m e n t s c o u l d , o f c o u r s e , b e i m p r o v e d b y m e s h r e f i n e m e n t , particularly in the beam problems where the mesh is only a single element deep in the slender d i r e ct i o n s. O n t h e o t h e r h a n d , m o s t u s e r s e x p e c t t o g e t a c c u r a t e r e s u l ts w i t h c o a r s e m e s h e s , a n d t h e a u t h o r s h a v e c o n s t r u c t e d t h e t e s t p r o b l e m s i n t h i s s p ir i t .
Table 7 Patch test results Maximum error in stres QUAD2 Constan!ostress loading Constant-curvature loading
QUAD4
QUAD8
HEXA(8)
HEX20
HEX20(R)
0.0
0.0
18.0~
0.0
0.0
0.0
30,7~
0.0
51.6~
N/ A
N/ A
N/
R.H. MacNeal, R .L Harder / Proposed standard set of problems to test FE accuracy
15
The most distu rbing failure of the QU 2 e l e m e n t is i ts in a b i l i t y to g et a p a s s i n g g r a d e f o r i n - p l a n e b e n d i n g o f t h e s t r a i g h t b e a m w i t h r e g u l a r e l e m e n t s h a p e . T h i s f a i lu r e i s e c h o e d i n t h e curved beam and the Scordelis-Lo roof problem s and was responsible for the developm ent of t he Q D M E M I a n d Q D E M 2 m e m b r a n e e le m e n t s n N A S T R A N ~ , n e it he r o f w h i c h p r ov e d t o b e a s i g n i f ic a n t i m p r o v e m e n t . T h e s p e c t a c u l a r f a il u re o f t h e Q U 2 i n th e t w i s te d b e a problem s due to the misalignm ent of the mom ents at the boundaries between adjacent e l e m e n t s w h e n t h e e l e m e n t s a re s ke w e d . T h e m i s a l i g n m e n t p r o d u c e s a r e s u l t an t m o m e n t a b o u t h e n o r m a l t o t h e s u r f a c e w h i c h is n o t r e s is t e d . In s p i t e o f t h e s e w e a k n e s s e s , th e Q r e c e i v e d a n e x c e l l e n t g r a d e o n t h e s p h e r i c a l sh e l l p r o b l e m i n c o n t r a s t t o t h e h i g h e r o r d e r elements which received C's or worse. T h e o n l y f a il u r e s o f t h e Q U
4 e l e m e n t a r e l o ck i n g f o r t he s t r ai g h t b e a m w i t h i n - p l a n e
Table 8 Results for straight cantilever beam Tip loading direction
Norm alized tip displacement in direction of load Q UA D2
Q UA D4
QU A'D 8
HE XA (8)
HEX20
HEX20( R )
Extension In-plane shear Ou t-of-p lane shea r Tw ist
(a Rectangular elements 0.992 0.995 0.032 0.904 0.971 0.986 0.566 0.941
0.999 0.987 0.991 0.950
0.988 0.981 0.981 0.910
0.994 0.970 0.961 0,904
0.999 0.984 0.972 0,911
Extension In-p lane shear Out-of-plan e shear Tw ist
(b Trapezoidal elements 0.992 0.996 0.016 0.071 0.963 0.968 0.616 0.951
0.999 0.946 0.998 0.943
0.989 0.069 0.051 0.906
0.994 0.886 0.920 0.904
0.999 0.964 0.964 0.918
Extension In-plane shear Ou t-of-p lane shea r Tw ist
(c Parallelogramelements 0.992 0.996 0.014 0.080 0.961 0.977 0.615 0.945
0.999 0.995 0.985 0.965
0.989 0.080 0.055 0.910
0.994 0.967 0.941 0.904
0.999 0.994 0.961 0.913
a The g ood to excellent results show n here for H EX 20(R ) w ere obtained in spite of singularity in the be am 's stiffness matrix. Table 9 Results for curved beam Tip loading
Normalized tip displacement in direction of load
di ec on
QUAD2
QUAD4
QUAD8
In-p lane (vertical)
0.025
0.833
1.007
0.880
0.875
1.00
Ou t-of-p lane
0.594
0.951
"0.971
0.849
0.946
0.959
(8 )
20
HE X20 (
Table 10 Results for twisted beam Tip loading
Norm alized tip displacement in direction of load
direction
QUA D2
QUA Da
QUAD 8
In-p lane
100.4
0.993
0.998
Out-of-plane
228.9
0.985
0.998
(8 )
HEX20
H EX20( R
0.983
0.991
0.993
0.977
0.995
0.999
R . H . M a c N e a l, R . L . H a r d e r / P r o p o s e d s t a n d a r d s e t o f p r o b le m s to t es t F E a c c u ra c y
Table I1 Results for rectang~L~-~.=p l a ~ e s i m p l e s u p p o r t s : u n i f o r m l o a d (a A s p e c t r a t i o = 1 . Nu mb er of node spaces a per edge of model
Aspect ratio =
N o r m a l i z e d l a te r a l d e f l ec t i o n a t c e n t e r QUA D2
QUA D4
QUA D8
H E X A ( 8)
HEX20
H E X 2 0( R
0.971 0.995 0.998 0.999
0.981 1.004 1.003 1.002
0.927 0.996 0.999 1.000
0.98 9 0.998 0.999 1.000
0.023 0.738 0.967 0.991
1.073 0.993 1.011 1.008
.0
N u m b e r o f n o d e s p ac e s a per edge of mod el
N o r m a l i z e d l a t e r a l d e f le c t i o n a t c e n t e r
2 4 6 8
QUAD 2
QUA D4
QUAD8
H E X A ( 8)
HEX20
HEX20(R)
0.773 0.968 0.993 0.998
1.052 0.991 0.997 0.998
1.223 1.003 1.000 1.000
0.955 0.978 0.990 0.995
0.028 0.693 1.066 1.026
1.139 0.995 1.024 1.006
a F o r e l e m e n t s w i t h m i d s i d e n o d e s, t h e n u m b e r o f e l e m e n t s p e r e d g e o f m o d e l i s e q u a l t o o n e - h a l f t h e n u m b e r o f n o d e spaces.
Table 12 R e s u lt s o f r e c t a n g u l a r p l a t e c l a m p e d s u p p o r t s : c o n c e n t r a t e d l o a d (a A s p e c t r a t i o = 1 . 0 N u m b e r o f n o d e s p a ce s per edge of model 2 4 6 8
N o r m a l i z e d l a t e ra l d e f l e c t i o n a t c e n t e r QUAD 2
QUA D4
QUAD 8
H E X A ( 8)
HEX20
H E X 20 ( R
0.979 1.008 1.006 1.005
0.934 1.010 1.012 1.010
1.076 0.969 0.992 0.997
0.885 0.972 0.988 0.994
0.002 0.072 0.552 0.821
0.983 0.433 0.813 0.942
( b ) A s p e c t r a t i o -- 5.0
Nu mb er of node spaces per edge of m odel
N o r m a l i z e d l a te r a l d e f l ec t i on a t c e n t e r QUAD 2
QUA D4
QUAD 8
H E X A ( 8)
HEX20
HEX20(R
2 4 6 8
0.333 0.512 0.638 0.723
0.519 0,~63 t£~ ,~ 0.972
0.542 0.754 0.932 0.975
0.321 0.850 0.927 0.957
0.001 0.041 0.220 0.374
0.363 0.447 0.721 0.867
Table 13 R e s u l ts f o r S c o r d e l i s - L o r o o f N u m b e r o f n o d e s p a ce s per edge of mod el
2 4 6 8 10
N o r m a l i z e d v e r t ic a l d e f l e c t io n at m idpoint of free edge QUA D2
QUA D4
QUAD S
HEXA(8)
HEX20
HEX20(R
0.784 0.665 0.781 0.854 0.897
1.376 1.050 1.018 1.008 1.004
1.021 0.984 1.002 0.997 0.996
1.320 1.028 1.012 1.005 -
0.092 0.258 0.589 0.812 -
1.046 0.967 1.003 0.999
17
R . H . M a c N e a l , R . L . H a r d e r / P r o p o s e d s t a n d a r d s e t o / p r o b l e m s t o t e st F E a c cu r a c
bending and irregular elements and locking for the thick-walled cylinder with nearly incomp r e s s ib l e m a t e r i a l. T h e s e a r e s t a t e - o f -t h e - a r t p r o b le m s f o r f o u r - n o d e d e l em e n t s . T h e o n l y o t h e r test where QUAD4 shows weakness is the curved beam with in-plane loadin g which is a m a n i f e s t a t io n o f t h e f i rs t m e n t i o n e d l o c k i n g p r o b l e m . s u r p r i s in g i n v ie w o f t h e i m p o r t a n c e u s u a l l y a t t a c h e d t o t h e p a t c h t e st . O n t h e o t h e r h a n d , t h e C
Table 14 Results for spherical shell Nu ber of node spaces per edge of model
N o r m a l i z e d r a d ia l d e f l e c t io n a t l o a d p o i n t QUA D2
QUA D4
QUAD 8
HEXA(8)
HEX20
2
0.928
6 8 10 12
0.990 0.986 0.98 4 0.982
0.972 1.024 1.013 1.005 1.001 0.998
0.025 0,121 0.494 0.823 0.955 0.992
0.039 0.730 0.955
0.001 0.021 0,097
0.990
HE×20( R
00162 0.776 0.972
Table 15 Results for thick-wailed cylinder P o i s s o n ' s r a t io 0.49 0.499 0.4999
N o r m a l i z e d r a di a l d i s p l a c e m e n t a t i n n e r b o u n d a r y QU AD 2
QUA D4
QUA D8
HEXA(8)
HEX20
H E X 2 0( R
~.6-¢3 0.156 0,'318
0.846 0.359 0.053
!.000 0.997 0.967
0.986 0.980 0.986
0 .999 0 .986 0.879
1.000 1.000 1.000
Table 16 S u m m a r y o f t e s t r e s u l ts f o r s h e ll e l e m e n t s Test
(1) (2) (3) (4) (5) (6)
Patch test Patch test Straight beam , extens ion S t r a ig h t b e a m , b e n d i n g S t r a ig h t b e a m , b e n d i n g S t r a ig h t b e a m , b e n d i n g (7 S t r a ig h t b e a m , b e n d i n g (8) Straight beam , twist (9) Cur ved be am (10) Curved beam (11) Twisted beam ( 1 2) R e c t a n g u l a r p l a t e ( N = 4 ) (13) Scordelis-Lo roof (N -- 4) ( 1 4 ) S p h e r i c al s h e ll ( N = 8 ) (15) Thick-walled cylin der (r -- 0.4999)
E l e m e n t l o a d in g
Element
In-plane
shape
Out-of-plane
× × × x × ×
x
×
x x × x
x
x
x
x
Irregular Irregular All Regular I r r e g u l ar Regular Irregular All R_,~ular Regular Regular Regular Regular Regula
Regu lar
N u m b e r o f f a i le d t e s ts ( D ' s a n d F ' s ) a R e g u l a r m e a n s t h a t e l e m e n t s h a p e h a s n o t b e e n i n t e n ti o n a l ly d is t o r te d .
QU AD
QUAD4
A D A F F B B D F D F C D A F
A A A B F A B R C B A B B A F
9
2
QUAD8
18
R.H. MacNeal,
. L H a r d e r / P r o p o s e d s t a n d a r d s e t o f p r o b l e m s t o t e s t F E a c c u ra c y
n o t e d f o r t h e s p h e r i c a l s h e ll i s d i s t u r b i n g i n v i e w o f t h e p r a c t i c a l i m p o r t a n c e o f t h i s p r o b l e m T h e H E X A ( 8 ) e l e m e n t p a s se s t h e p a t c h t e s t a n d p e r f o r m s w e l l w h e n t h e e l e m e n t s a r e r eg u l a r but behaves poorly for bending loads on irregular elements. This is probably why it also converges slowly in the spherical shell problem and gets C's in the curved beam test where the elements deviate slightly from rectangular shape. It is noteworthy that the HEXA(8) gets an excellent mark for the nearly incompressible thick-walled cylinder a s the result of a design feature recently built into the element. The standard isoparametric HEX20 does fairly well on all the beam pr oblems but fails the i m p o r t a n t p l a t e a n d s h el l p r o b l e m s . T h e s e f a i lu r e s a r e s o m e w h a t a l l ev i a te d b y t h e i n t ro d u c t i o n o f re d u c e d o r d e r i n t e g r a ti o n a s e v i d e n c ed i n t h e re s u l ts f o r H E X 2 0 ( R ) . O n t h e o th e r h a n d , reduced order integration is responsible for the singularity of the stiffn ess matrix noted in the straight beam tests. The spuri'ous modes associated with this singularity inv olve a second order distortion of the cross-section which has no component in the direction of th e applied load. The rather erratic results for HEX20(R) shown in Table 12(a) are probably related to the reduced order of integration. Thus, while reduced order integration appears to im prove performance for many applications, it can lead to sudden unexpected failure in others. Turning our attention away from the elements to the test problems, it is seen that the p r o b l e m s a r e c a p a b l e o f ev o k i n g b o t h g o o d a n d b a d r e s u lt s . I n o n l y o n e p r o b l e m , e x t en s i o n o f a straight beam, did all of the elements get A's. (This is a throwaway t est in the sense that any element that cannot pass it should be thrown away.) All but three of the remaining fourteen tests produced failing grades for one or more elements and all but three of the remaining tests p r o d u c e d a g r a d e o f A f o r o n e o r m o r e e l em e n t s . Two important questions are whether any of the tests are d uplicates of other tests and whether the tests are sufficient. The first question can be answered by noting that the letter grade test results are not duplicated across the six elements by any pair of tests. The second question is more difficult to answer because the list of error causing pa rameters is not closed. W e c a n , h o w e v e r , e x a m i n e w h e t h e r t h e p a r a m e t e r s d e s c r i b ed i n t h e i n t r o d u c t i o n a s t h o se w h i c h affect accuracy are present in one or more problems. The result of the examination is
Table 17 Summ ary of test resu lts for solid elements Test
Element shape
HEXA(8)
(1,2) Patch test (3) Straigh t beam , extension (4,6) Straight beam, bending (5) Straigh t beam , bending (7) Straigh t beam , bending (8) Straigh t beam , twist (9) Cu rved beam , in-plane loading (10) Curved beam, out-of-plane loading (11) Twistedbeam 02) Rectangularplate (N -- 4) (13) Scordelis-Lo roof(N = 4) (14) Sphericalshell (N = 8) (15) Thick-walledcylinder (v = 0.4999)
Irregular All Regular Irregular a Irregular b All Regular Regular Regular Regular Regular Regular Regular
A A
Num ber of failed tests (D's and F's) a Bending in p lane of irregularity. b Bending out of plane of irregularity. c In spite of singu lar stiffness matrix.
A
0 A A B
F F B
B B B
C
C
B
F
C B
B B
B D A
F F C
3
3
20 (R
A B
A B B
R.H. MacNeai. R.L. Harder
Proposed standard s et of problems to test FE accurac.v
19
T a b l e 18 Accuracy affecting parameters vs. test problems P a r a m e t e r s which affect accuracy
Test problems where illustrated
Response excited bv loading Extension In-plane shear and bending Out-of-plane shear and bending Twist Higher order stress gradients
I. 3,8.9. 10. 13. 14. 15 1.4.5.9. 11.13. 14. 15 2. 6. 7. 10. 11.12. 13, 14 2. 8. 10. 12. 13. 14 12. 13. 14. 15
Element shape Aspect ratio Skew Taper Warp Displaced edge nodes
ll 1.2.3,5.7.8 i.2,3,5.7. 8.9. 10. 13, 14. 15 11 10. 13. 14. 15
Problem geometry Single curvature Doub le curvatur Free boundary Simply-supported boundary Clamped boundary Symmetric boundary Enforced motion and boundary Thinness Slenderness
All except 15 3 . 4 . 5 . 6 . 7 . 8 . 9 , 1 0.
Material Nearly incompressible Anisotropic
15 None
13 11.14 3 . 4 , 5 , 6 . 7 , 8 . 9 , 1 0. 1 1, 1 3. 1 4. 1 5 12. 13 3. 4. 5, 6. 7. 8. 9. 10. 11, 12 12. 13, 14. 15
summarized in Table 18, where it is seen that all of the parameters affecti ng accuracy except anisotropic material are present in one or more of the proposed test problems and more of ten in several of them. Concluding remarks The present paper has described a proposed standard set of finite element test problems in sufficient detail to allow their reproduction in any general purpose fini te element code. The problems, while few in number and relatively simple, have been designed to include, collectively, nearly all of the parameters which have important effects on element accu racy. The proposed test problems have been exercised on a representative set of quadrilateral shell and brick elements with different orders of complexity, and it has been shown that al most every problem is capable of evoking results ranging from excellent to poor. As developers of finite elements, this is exactly the sort of test problem library we would wa nt to help us find and correct element errors and weaknesses. It is also hoped that the user community will be able to profit from them in selecting elements and mesh sizes. Acknowledgment The work reported here has been prepared in response to a request from the AIAA Standards Committee on Structural Analysis, chaired by Kevin J. Forsberg.
20
R.H. MacNeal, R.L. Harder / Proposedstandard set of problems to test FE accuracy
References [1] ROBINSON,J., and S. B ACKHAM. An e valuation of lower order m embran es as contained in the M SC /N AS TR AN ASAS, and PAFEC FEM systems", Robinson and Associates, Dorset, England, 1979. [2] ROBINSON,J., and S. BLACKHAM," A n e v a l u a t i o n o f p l a t e b e n d i n g e l e m e n t s : M S C / N A S T RA N , A S A S . P A F E C , ANSYS, and SAP4", Robinson and Associates, Dorset, England, 1981. [3] The NA ST RA N Theoretical Manual level 16.0, NASA SP-221(03), Section 15.2, 1976. [4] IRONS,B. an d S. AHMAD, Techniques of Finite Elements, Wiley, New York, p. 155, 1980. [5] S¢ORDVUS,A.C. a nd K.S. Lo, " Co pu ter analysis of cylindrical shells", J. Amer. Concr. Inst. 61, pp. 539-561, 1969. [6] ZIENKIEW_'CZ,O.C., The Finite Element Method, McGraw-Hill, Londo n, p. 417, 1977. [7] MORLEY,L.S.D. and A.J. MORRIS, Conflict Between Finite Elements and Shell Themy, Royal Aircraft Establishment Report, London, 1978. [8] TIMOSH~NKO,S., and S. WOINOWSKY-KRIEGER, Theory of Plates and Shells 2nd ed., McG raw-H ill, New York, pp. 120, 143, 202, 206, 1959. [9] The NAS TR AN Theoretical Manual Level 16.0, NAS A S P-221(03), Section 5.8.3, 1976. [10] MACNEAL,R.H., "A simple quadrilateral shell element", Computers and Structures 8, pp. 175-183, 1978. [11] M^cNvA L, R.H., "D eriva tion of element stiffness matrices by assum ed strain distributions", Nuclear Engineering and Design 70, pp. 3-12, 1982. [12] MAcNEAL, R.H., "Sp ecifications for the QU AD 8 qu adrilateral curved shell elemen t," MacN eaI-Schw endler Corp. Memo. RHM-46B, 1980. [13] MACNEAL,R.H., "S pecifications for the HE XA element", M acN eai-Sch wen dler Corp. Memo. RHM -38C, 1976. [14] HARDER,R.L., "M odifie d isoparametric elemen ts and patch tests", Ma cNeaI-S chwen dler Corp. M emo. RLH -35, 1982. [15] IRONS,B.M., "E ngin eerin g application of numerical integration in stiffness method" J. AIAA !4, pp. 2035-2037, 1966. [16] ZI~NKIEWlCZ,O.C., R.L. TAYLOR~nd J.M. To o, " Red uced integration techniques in general analysis of plates and shells", Internat. J. Numer. Method. Engrg. 3, pp. 275-290, 1971.