mixed formulation

Finite Elements in Analysis and Design 2 (1986) 357-375 North-Holland

357

ON THE NUMERICAL PERFORMANCE OF THREE-DIMENSIONAL THICK SHELL ELEMENTS USING A HYBRID/MIXED FORMULATION * W. GRAF, T.Y. CHANG and A.F. SALEEB Department of Civil Engineering, Universityof Akron, Akron, OH 44325, U.S.A. Received 30 September 1985 Revised 2 December 1985 Abstract. Three-dimensional thick shell elements with 8, 16, and 18 nodes are formulated by using the hybrid/mixed method. In bending applications, these elements are free from locking effect and give improved stress predictions. Finite element equations are derived from the Hellinger-Reissner variational principle in which both the displacement and stress fields are approximated by independent interpolation functions. For the assumption of stress parameters, three guidelines are followed: (i) suppression of kinematic deformation modes, (ii) invariant element property, and (iii) the constraint index exhibited by the element, when applied to constrained-media problems, must be greater than or equal to one. Numerical results are presented to show the element's behavior characteristics regarding sensitivity to locking, distortion effect (patch tests), mesh convergence and the accuracy of stress evaluation.

Introduction Within the framework of finite element analysis of shells or shell-like structures, their geometric configurations are often represented by mid-surface thin shell elements. A great variety of such elements has been developed and published in the literature for different application purposes (see, e.g., [1,2,10,32]). However, in many cases two-surfaced thick shell (or solid) elements offer several advantages over mid-surface elements. For example, in transition regions, where the shell is part of a solid structure or at the intersection between two cylinders, use of mid-surface shell elements will present special boundary conditio~, problems, which do not exist when two-surfaced elements are utilized. Additionally, two-surfaced elements can easily capture thickness deformation effects whereas mid-surface elements cannot. Moreover, the effect of shearing deformation is generally included in the thick shell formulations. In current practice, many thick shell structures are simply modeled by three-dimensional solid elements with 8, 16, or 20 nodes. However, when the thickness of the structure becc,.qles small, numerical difficulties may be encountered in two different forms: (i) solid elements may exhibit shear locking phenomenon, and (ii) stiffness coefficients associated with the thickness direction become much larger than other coefficients and this may induce ill-conditioning of the system equations to be solved [10]. Several solutions to the locking problem have been proposed by various researchers. One alternative is to use a reduced integration technique for the evaluation of element stiffness matrices [12,13,14,26,29,33]. However, solid elements will still eventually lock when the thick* This work was supported by the NASA Lewis Research Center under Grant No. NAG3-307. 0168-874X/86/$3.50 © 1986, Elsevier Science Publishers B.V. (North-Holland)

358

W. Graf et al. / Three-dimensional thick shell elemolts

ness is reduced further. Mofe,:~'~'~:~,~l~e use of reduced integration techniques can cause rank deficiency of the element stiffr,o ~ matrix (also known as spurious singular or kinematic deformation modes), which ~t t~t~ may yield questionable analysis results under certain boundary or loading condRio~s ~ t ~ ~6]. ~'~oxe recently, Belyigchko et al. [4,5] have introduced a projection (stabilization) m¢~.~_~ ~o be used in conjunction with reduced integration to alleviate the locking phenomenon and also eliminate the formation of kinematic deformation modes. Another noteworthy approach was that taken by Wilson et al. [31] who proposed the introduction of incompatible displacement modes at the element level to improve the numerical performance of solid elements for thick shell analysis. While Wilson's incompatible element approach has indeed resolved the locking problem, the associated elements have nevertheless been found to be sensitive to geometric distortions. In fact, a consistently formulated incompatible element of this type will not pass the patch test [16] unless some adjustments are made [10,16]. More specifically, it involves the use of a 'constant' Jacobian when determining the contribution of incompatible modes to element stiffness matrices, and disregarding the loads associated with these modes in stress calculations. A fourth alternative to improve element performance is to adopt a hybrid or hybrid/mixed formulation, the initial development of which is attributed to Pian [18]. Two different versions of the hybrid formulation are generally available [19,20]. The first involves assumed equilibrating stresses only within each element and compatible displacements along the interelement boundary. A new formulation presented by Pian et al. [22,23,24] termed the hybrid/mixed method relaxes the stress equilibrium condition initially, but then enforces equilibrium through constraint equations by using internal displacement variables as Lagrange multipliers in the Hellinger-Reissner principle. Indeed, there is a close resemblance between the hybrid method and other techniques mentioned above. For instance, some discussions on the equivalence and differences bet~°een the hybrid method and the incompatible displacement models are given in [11,15,25]. Different forms of the hybrid/mixed variational principles were also referenced to justify both the utilization of a reduced integration technique [14,17] and the projection stabilization procedures [4]. In the hybrid/mixed method, one has the freedom to independently interpolate both the stress (or strain) and displacement fields. By carefully selecting appropriate parameters for the assumed fields, it is possible to arrive at computationally efficient and numerically reliable elements for application purposes. However, it appears that the potential and versatility of this method have not yet been fully realized. Hence, further research in this direction is still needed. In the present paper, the hybrid/mixed method is used to formulate three-dimensional thick shell elements with 8, 16, and 18 nodes for the bending analysis of shell-like structures. To determine the numerical performance of these elements, four parameters have been investigated: (i) element sensitivity to locking (shear and membrane), (ii) effect of element distortion (patch test), (iii) mesh convergence characteristics, and (iv) accuracy of stress evaluation. Numerical results obtained for the hybrid/mixed elements are then compared with those of equivalent displacement-based elements.

Finite element formulation

In this study, the formulation is restricted to linear elastic materials with small deformations. All kinematic and kinetic quantities (e.g., displacements, strains, and stresses) are referred to a rectangular Cartesian system, x~ (i ffi 1, 2, 3), unless defined otherwise. Beginning with the Hellinger-Reissner variational principle, which may be derived from

W. Grafet al. / Three-dimensionalthickshellelements

359

either the potential or complementary energy principle via the Lagrange multiplier technique [30], the following functional ~rR for an element should be maintained stationary:

ftR-" fv[--½oTso+oT(Du)] d V - / s o T T U d S - f s T T ( u - u )

dS,

(1)

where = element stresses, -elastic compliance matrix, element displacements, boundary tractioas, T () prescribed quantity, Du defines the strain-displacement relation, V -- element volume, &,S ffi portions of the total element boundary surface area S over which tractions and displacements, respectively, are prescribed.

0

S

=

In the above, the effect of body forces is ignored. Invoking the stationarity of ¢tR with respect to variations in o and .u, the following Euler equations are recovered:

Strain-displacement relations: • -- S o ffi D u in V.

(2)

Equilibrium equations: Do---0 in V. Boundary conditions:

(3)

lr=T u-~

(4) (5)

OnSo. on Su"

In finite element approximations, the displacements u and stresses o within an element are interpolated in terms of nodal displacements, q, and stress parameters, fl, as

.=Nq,

o=P~,

(6), (7)

where N and P are the matrices of interpolation functions for element displacements and stresses, respectively; they are polynomial functions of local coordinates n~ ( i - 1, 2, 3). Some comments on the rationale for the particular choices of the interpolation functions for stresses are given in the following section. Here, we note that the stress parameters may be chosen in such a way that the stress equilibrium conditions of equation (3) are satisfied pointwis,.~ v,:tZ~2a ~ ; ~:'-ment. The resulting formulation is then referred to as an 'equilibrium' hybrid element [19,20]. Tiffs has been the most popular approach in the early developments of hybrid e!e:,.,ea~s. However, recently in a new formulation by Pian et al. [23,24], the stress equilibrium conditions are not enforced initially, but are brought in, in an 'approximate' (integral or variational) sense, through the use of incompatible (or internal) displacement parameters. In our study, we chose to use only compatible displacement functions within an element. For the approximation of stresses, two alternatives are considered: (i) stress functions satisfy the equilibrium conditions identically in the element (local) coordinates, and (ii) the equilibrium requirement is completely relaxed. By substituting equations (6) and (7) into (1) and grouping appropriate terms, the functional ¢rx can be written as ~ g ffi -- ½ f l V H f l ÷ f l r G q - ~ T q ,

(8)

W. Graf et al. / Three-dimensional thick shell elements

360

where H, the flexibility matrix, and G are defined by

H= fvPTSPdV, G---f P (ON) dV,

(9) (10)

and Q is an equivalent nodal force vector,

Q= fsNXTdS.

(11)

Invoking the stationarity of (8) with respect to the stress parameters yields fl in terms of q,

flffi H-1Gq.

(12)

Substitution of (12) into (8) yields "R =

1

TK q -

QTq,

(13)

from which the stiffness matrix for the hybrid/mixed element is given by

K=GXH-~G.

(14)

In the above, K is symmetric due to the symmetry of H.

Approximations for displacements and stresses The general formulations for hybrid/mixed elements were outlined in the previous section. To derive specific finite element equations, three element types are considered: an 18-node element and two sub-elements with 8 and 16 nodes, respectively. The 18-node (HM18) element is defined by 9 nodes each on the upper and lower surfaces as shown in Fig. 1, with each node having three translational degrees of freedom. The 16-node (HM16) and 8-node (HM8) elements are obtained by removing the upper and lower central nodes, and all mid-side nodes, respectively. As mentioned in the preceding section, only compatible displacement functions are used, these being the familiar interpolation functions for three-dimensional isoparametric solid elements [2,32]. The assumption on the stress functions, i.e., entries of the P matrix, however, is not as straightforward. In this section, a discussion on the basic strategy used in selecting an appropriate set of stress parameters is outlined. Many of the useful concepts, procedures, and criteria proposed in the literature regarding this point are utilized here as guidelines. Three important considerations should be carefully examined ,vhen selecting the stress parameters: (i) all kinematic deformation modes must be avoided, (ii) whether or not the stress components should satisfy the equilibrium conditions, and (iii) ability of the resulting element to handle constrained-media problems (specifically, this implies herein the avoidance of locking phenomenon for 'thin' elements) based on a convenient 'constraint counting' measure. Each of these items is discussed separately below. Suppression of kinematic deformation modes is perhaps the most important requirement in the hybrid/mixed formulation. Early in the development of the assumed stress hybrid elements [18], it was recognized that a necessary condition for the stiffness matrix to be of sufficient rank is that the number of stress parameters should be greater than or equal to n - r, where n is the total number of degrees of freedom in the element, and r the number of rigid body degrees of freedom. Noting that examples exist where kinematic deformation modes develop even when the stress parameters satisfy this minimum requirement, Plan and Chert [21] offer additional

W. Graf et al. / Three-dimensional thick shell elements

361

!

f

t?@ / %

X2,U 2

/

XI,U I

X~l ,U s Fig. 1. A thick shell element.

insight as to the criterion for choosing the proper set of stress terms. Based on deformation energy considerations, a scheme was proposed whereby the total number of stress parameters are minimized while simultaneously suppressing all kinematic deformation modes. The basic procedure is to choose the stress parameters in such a manner so as to have a /]-term corresponding to each of the terms in the strain expressions obtained from the strain-displacement relations. If needed, additional stress terms may then be introduced in order to satisfy equilibrium conditions. A second consideration involves the question of whether or not it is necessary to completely satisfy the pointwise equilibrium requirement in the initial assumption for the stress polynomials. Since the global equilibrium is, after all, enforced through the variation of the functional ~rR, it is our view that the local (or pointwise) equilibrium requirement can be relaxed either partially or fully. Much can be gained if indeed the local equilibrium requirement is relaxed. First, the tedious task of selecting the stress polynomials in trying to satisfy the equilibrium equations is simplified. We note here that in some cases it is extremely difficult, or even impossible, to satisfy the 'exact' equilibrium conditions imposed on the physical stress components; for example, in situations involving elements with distorted geometries or in nonlinear applications. Further, it has been found that introduction of additional stress terms to satisfy equilibrium may again induce kinematic deformation modes, which in turn requires that more stress terms be added [22]. In order to shed some light on the question of equilibrium requirement, a numerical study is carried out by both including and deleting those stress terms involving equilibrium. The analysis results, to be described in the next section, actually indicate the favorable effect of relaxing the equilibrium constraints. Adhering to the guidelines set forth by Plan and Chen [21], namely, minimizing the total number of stress parameters while simultaneously suppressing all kinematic deformation

362

W. Graf et aL / Three.dimensional thick shell elements

Table 1 Stress parameters and corresponding displacement modes to be suppressed for the 18-node hybrid/mixed element

1 2 S

Or,

°ss

°u

Ors

#,

#2

#3

#a

u----r

v=S

w=t

u=s;

#7 u=r 2

#16

#8

#17

u

t rs

St

=

rst

0=5

2

w=

#32 rt

u=t;

w=r

r2

w=

r2

ll = s t

W=S

--(#35+#40)

#34

--(#7+#33 )

--(#17 + #32)

#+l

#19

#27

#37

o = rs 2

w ---- r s t

w = r2s

#11

#20

-

#21

U = r2t

0 = rst

#12

1

#22

#28

-- i#19

0 = r2s

w = r2t

#29

u = rs 2

w = S2t

=

rs 2

#20

-

#13

U

-

w=

2

s2t

#12

#14

#40

#36

ZI "~ S 2

u = r2s

=

rs

#33

#10

0

w=

w = st

v = st

rst

W=$

#26

#is

=

O-~l;

#39

#35

o=

u = rt

r2

s2

#25 rs

#6

#5 v=r

#9

U

rt

rs

o=

°"t

1

-- ~#41

#42 U ---- r 2 t

-

-

½#10

#38

1

--

2#37

u = s2t

#23 r2st

r2s

r 2t

u

~-

rs2t

#44

#30

#45

O m r25 2

w m r2$t

w = r25 2

#24

- #45

- -~#23

- ~

V = r2st

.2

s2t

&~ U m r2s 2

#31 rs2t

W--

#is

- ~#,4

- #43

u " rs2t

rs2t r 2st

#46 U m r252t

#4~ 0 = r25 2t

r 2s 2

#4~ W

s2t 2 r':

=

r 2 s 2t

- I#46 - ½#4,

modes, an expression for the stress polynomials of an HM18-element is obtained in terms of the natural coordinates (r, s, t), where r = rt, s - r2, t = r3. The form of the P matrix as defined in equation (7) is given in Table 1, which contains the stress parameters and corresponding displacement modes to be suppressed for the HM18-element (with regular

W. Graf et al. / Three-dimensional thick shell elements

363

geometry). Here, the equilibrium terms are included as indicated by those blocks not containing displacement terms. For instance, six terms must be added to %s in order to satisfy the equilibrium conditions: "rrs . . . . .

~~12st-//2ort-

½//19r 2 - ½//10$2- l//23r2t- 1//14S2t.

(15)

With the stress parameters listed in Table 1, all columns of the G matrix are independent, and hence no kinematic deformation modes are present in the HM18-element. By deleting those terms corresponding to the basic deformation modes given by u, v, and w (where u - u l, v - - u 2, and w - u3) equal to r2s 2 and r2s2t, i.e., //43 through//48 along with their associated equilibrating terms, the stress pattern for the HM16-element is obtained. In a similar fashion, the stress parameters for the HM8-element can also be obtained directly from Table 1. The third consideration in our choice of the stress parameters is to check the element behavior so that any locking problem is precluded for its application to 'thin' structures. A convenient way to judge such behavior is to use the method of constraint counting, often termed the 'constraint index' as suggested by Hughes et al. [14], Malkus and Hughes [17] and Pugh et al. [261. As mentioned in the above, in order to suppress all kinematic deformation modes it is necessary to assign at least one //-term in the stress polynomials to each of the basic deformation modes involved in the corresponding displacement field. Whether the aforementioned //-term is assigned to any one of the normal or shear stress components is, however, immaterial. It is precisely this freedom in assigning the//-terms to various stress components which enables us to control the locking phenomenon induced by the 'troublesome' stresses. By exercising this freedom, the number of element constraints is minimized and the 'constraint index' is thus improved [14,17]. Based on the above considerations, stress parameters for elements HM18, HM16 and HM8 were selected and are listed in Table 1. For illustrative purposes, we consider the HM18-element shown in Fig. 1 (with regular geometry), to model the behavior of a thin shell. With t being the through-thickness direction, in the limit as the thickness approach zero, the stiffness contribution associated with the normal stress component ¢, and out-of-plane shears ¢rt and ¢, should vanish. This requirement will thus impose (approximately) 20 constraints, in reference to the stress parameters given in Table 1. However, when adding an HM18-element to an existing finite element mesh, the resulting constraint index is 4 (i.e., 2 4 - 2 0 ) . Since the constraint index of the element is greater than or equal to one, no locking is expected i~, this case. Since the P matrix is expressed in terms of the natural coordinates instead of the Cartesian system, the resulting element stiffness is frame-indifferent. However, the stress components so defined in the natural coordinates are not the physical quantities. In order to obtain the physical stresses, a transformation is necessary: ax~ ~xj

(16)

The above transformation is evaluated at a constant location, namely the centroid of an element. If evaluated otherwise, the order of the stress polynomial would be altered which may either trigger the kinematic deformation modes or result in an increase in the number of constraints so that locking may once again become a problem.

Case studies The aforementioned thick shell elements have been implemented into a general purpose finite element program NFAP [7] for computational purposes. In order to determine the

364

W. Graf et al. / Three-dimensional thick shell elements

numerical performance of these elements, the following items are investigated: (i) shear locking vs. large aspect ratio, (ii) membrane locking, (iii) patch test, (iv) mesh convergence, and (v) stress evaluation. Results are compared with those obtained from the equivalent displacementbased elements. In addition, solution accuracy in connection with relaxation of equilibrium conditions in the hybrid/mixed elements is also noted. All computations are performed in double-precision on the IBM 3033 Computer at the University of Akron.

Shear locking As mentioned before, when conventional, displacement-based, lower-order isoparametric elements are used to analyze thin plates or shells, there is a tendency for the element to lock if its thickness is relatively small. This degradation in element performance is attributed to its inability to accurately represent the transverse shear behavior in the element stiffness formulation. With this background in mind, it is our intention to determine whether or not the hybrid/mixed elements also exhibit the shear locking phenomenon. For this purpose we consider a cantilevel beam subjected to an end shear force. The beam is modeled by one HM16-element or two HMS-elements with the parameter varied being the beam's aspect ratio, L/h, where L is the span of the beam and h the height. Analysis results for the HM16-element are shown in Fig. 2 along with those of the 16-node standard brick element. The hybrid/mixed element yields tip deflections corresponding to elementary beam theory for aspect ratios as high as 1000. For L/h greater than 2000, element deterioration is observed, primarily due to the large stiffness coefficients in the thickness direction. Results obtained using a 16-node standard brick element exhibit locking phenomenon for the entire range of aspect ratios. 0 la iu4

go • 4.

U

L[G[NO II;-liOOE OZBIq.llgiB~N'f |I.F.NIKN'! I ~ l t l l - l L I I l l H l r 1|8U|LJlillIUH M ' l r l l l g l [ O |

=.'.. laJ

2¢ 00

¢.,t bJ _J I.i..

,-,~. n_

4" ,

lid

4

4,

4,

4,

•

•

•

B

t,,q

I-. Om

Wm r,,Ic:;-

•

4.'

.J (¢

• O

z~..

.

,o.

.....

......

ELEHENT RSPECT-RRTIO

(--L/H)

.......

io,

Fig, 2, Effects of element aspect ratio on tip deflection of a cantilever beam subjected to end shear.

365

W. Graf et aL / Three.dimensional thick shell elements

Membrane locking .lust as element reliability is suspected in thin plate applications, when dealing with curved structures, membrane locking may result due to the increase in bending stiffness coupled with membrane terms. Stolarski and Belytschko [29] have indicated that when the predominant deformation in a curved structure is bending, a membrane-dominated response results if lower-order displacement elements are used. In this context, a fixed-fixed curved beam subjected to a concentrated force at the center, as shown in Fig. 3, is analyzed using the HM16-element. For this problem, again there is no need to use HM18-elements since the bending action is only in one direction. One-half of the curved beam is modeled, where now two cases of stress parameters are considered: (a) local equilibrium conditions are satisfied (in the natural coordinate system only), and (b) local stress equilibrium is completely ignored. In Fig. 4, the vertical deflection beneath the applied force is plotted as a function of the number of elements used. The HM16-element shows better mesh convergence rate than the 'equivalent' displacement-based element for both of the stress assumptions considered. Moreover, the improved element

P

E = 22000. Y= 0.00 R= 4.22 t -0.20 b= 0.10

~,. , I " = %

I #w

Fig. 3. Axial force and bending moment at a typical section of a curved beam subjected to concentrated load.

366

W. Graf et aL / Three.dime,sional thick shell elements I.SO ..

t.EO

a • 4-)

LEGENO

,=-=o=

=t=,,",--

NNtE-ELEHENIr IEOUILJIIIIUII SIIJSF|EOl HIItO-ELEHEli'f IEIIU|LllllllJII RELIIXEBI

t.tO t.O0 0.g0

Z 0

0.110

I-,. U

0.70 "

,.J b. W Q

O. 60 -

D bJ It4

a

Q

+ •

4.

0.60 o.qo

.3

E E 0 Z

O.SO 0.20 O.tO O

0.00

I

I

I

I

I

I

t

~

S

II

g

li

NUHBER OF ELEHENT$ Fig, 4, Deflection beneath point of load application for a fixed-fixed curved beam.

performance due to relaxation of the equilibrium condition is quite apparent. In either case, the hybrid/mixed element does not appear to show any sign of membrane locking. Patch test

An important convergence requirement for an element is to pass the patch test for arbitrary geometry. In bending analysis, the element assemblage must be able to assume a state of constant curvature [10,16]. For this purpose, we consider a cantilevel beam modeled by two

E = ?-£000. V= 0 , 0 0

O LVI

--7~

o o

m~. 'p

,

L

I

IO

r

I -1

Fig, 5. Element geometry for the patch test.

367

W. Graf et al. / Three-dimensional thick shefl elements I .~10 - - - -

01~ .b.

1.80-

LEG[H0 O-NOOE OlSPt~CEflEifl ELEIII[N1 HIM~.ELENEN1

I.tO X

- -I

1.00

4-

O.O0 -

Z 0

i i i

I-C.J L~ ..J ta. ILl Q

0.00 0.70 O. 60

4. 4.

O.

0.60 o.qo

UJ

0.30 2: n" O Z

0.20

' &

0 •

•

O.LO

•

0.00

•

•

i

O;SO

0.00

•

I

"i 1.10

I./60

i

I!. 011

O|$T0fl'[JON

It~60

i

• ..,

-,

$~0|

.

1~60

PAflAH['fEfl

. =..

,~..

U.O0 u

..............

|

I

tl.fm

s.N

(al

Fig. 6. Effects of element distortion on tip deflection of a cantilever beam subjected to pure bending (8-node). t.SO

L[G[HD :: ~ e - m t olm..r.,em t L u ~ Hlltll-lLillff liWll, lilllUII IAI|IPIEOI • • Iltltll-|l, l l q m I [ W l I . I I I I I U I I MI.IIIXI~Ill

t,t0 4.)

.

.

'

•

1.to X

...

1,00

A

•L

A

• Z

II

I

I

I

.

...,

--

4'

6 4.

O.liU

Q

÷

4.

0

0,I0-

I." C,J U,I 0 , 7 0 ° ,J ia. ta.I C:)

O.IO-

0.. .-,

0.60-

I--

o

UJ J~4 ..J

nZ nO Z

O.qO0.$0" 0.20O, tO 0.00 0.00

~ ' , OS0 I.O0

, I.S0

, = ~. 11.oe l s o

O]$T0flT|0H

,

~

~

,

s.ee

ago

u.oo

i.so

PgflgHE'fEfl

~ s.oo

!al

Fig. 7. Effects of element distortion on lip deflection of a cantilever beam subjected to pure bending (16.node).

W. Graf et al. / Three-dimensional thick shell elements

368

Rigid Diaphragm S u p p o ~ E v t R L

= = = = =

22000 0.30 1.0 100 200

V

J

Rigid Diaphragm Support Fig. 8. A pinched cylindrical shell.

distorted HM8- and HM16-elements, respectively, as shown in Fig. 5. The normalized tip deflection is plotted as a function of the distortion parameter, a, for the HM8-element in Fig. 6, and for the HM16-element in Fig. 7. The HM8-element shows a maximum of 62% stiffening when it approaches to a tetrahedron. On the other hand, the HM16-element is rather insensitive to geometric distortion. When the equilibrium condition is imposed, 17% error is observed if fully distorted. When relaxed, however, the error is now only 5%. By comparison, numerical results obtained using displacement-based elements deteriorate much more rapidly as evidenced by the 90% stiffening observed in the fully distorted 16-node element.

Mesh convergence for a pinched cylindrical shell To test the convergence characteristics of the elements developed, a cylindrical shell subjected to two concentrated forces (see Fig. 8) is analyzed using HM18-elements with mesh sizes ranging from 2 x 2 to 10 x 10. This problem has been one of the most frequently used for convergence studies [1,6,10]. Taking advantage of the triple symmetry condition, only one-eight of the cylinder is modeled. The load is calculated so as to produce a unit deflection beneath the point of application. The normalized vertical deflections vs. mesh sizes are plotted in Fig. 9. As evident from the results in Fig. 9, unlike the displacement-based element, the HM18-element does not exhibit any sign of locking. Although the convergence rate is slightly improved when the equilibrium conditions are enforced, the solution corresponding to the 10 x 10 mesh appears to slightly overshoot theory (U/Ue,c, = 1.04). In this example, it shows that the analysis results are relatively insensitive to the equilibrium requirement on stresses.

W. Graf et al. / Three-dimensional thick shell elements i

I.SO 1 I,EO

-,

~rs >t

369

j

.. +.LEGEflO I 1 - . o o ( ozm.,co~m i L r . . e m HNLli'.ELEHEWI IEBUILJIAIUII SAT|SF|EOI HHLE'ELItHENT IEOU|L]IIIIIUll HELAXE01

6 •

4)

.,

I,LO 4. ..

1,00

,.

A

, ,

1

4. 0.II0 -

2~ o.oo ' I-C.3 ..J LI. UJ t"J

0,70

I

0.60 0,60

LI3 1'14

O, YO

..I

3~

O.'JO

tD Z

O,~O

n~

O,LO .~ 0 0.00

|

0

~t

II

O

l

HESPI

SIZE

tN X

I0

N)

Fig. 9. Displacement convergence for a pinched cylindrical shell. II

0

t,OO

z

l.YO

n__

1.80

%

1.00

~" IIJ

O.IO

L~

.....

ee 4,+.

"

LEGEND

11-140|| OJmq.ACEHEWI |LFJIEWI HH'I,I-ILIIIMI IIWIL, JIBIUII M"IJSFJIIOI

• • HHLO-ItL|HIEIll I|BUIL.IIIIIUH RELAXEOI -- AHAI.'III CIL $|LUT|OH O

0.80

E O b.

O.qO

..J O:

0.20

•

•

iml

X O:

0.00

I

II

o "-O.l~O I1.1 j,t,i -O.YO .J C Z -O.80 O m 05 -O,Oil

l.t,,l -;.120 E

•

_

_

~

_

Ol

-+

-.-

_ ;++

i

_

++

+

""

.... +

_-.

=,;.

el

iml

c~ - I ,a0 . O :am 10 Z -!.~O

• O '"'

|

|O

i

|O ANGLE

' '

|

'

'

30 HERSURED

i

IIO FI:iOH

|

Ill

"

I

II0

VEflT!Cl:iL

Fig. 10. Axial force in a fixed-fixed curved beam.

I

70 ( O !

,=L

|

80

'

IlO

W. Graf et al. / Three.dimensional thick shell elements

370

Stress evaluation

As is well known, there is an overall reduction in accuracy when the stresses are evaluated from the gradient of assumed displacement functions in the displacement-based finite element formulation. In the hybrid/mixed formulation, however, stresses are determined directly from the assumed stress functions; hence, improved accuracy in the evaluation of stresses is anticipated. The curved beam problem of the section on "Membrane Locking" (see p. 365) is considered here again. Both fixed-fixed (180 ° subtended angle) and fixed-free (90 ° subtended angle) boundary conditions are analyzed. The beam is modeled by six HM16-elements with two variations as before: (a) equilibrium conditions are imposed, and (b) equilibrium conditions are relaxed. Depicted in Figs. 10-13 are normalized stress distributions (axial forces and bending moments) along the centerline of the beam. As was the case in predicting deflections, the stress patterns calculated from the displacement-based elements are greatly in error, oscillating about the analytical solution for both the fixed-fixed (Figs. 10 and 11) and fixed-free (Fig. 12) boundary conditions. On the other hand, results obtained using the HM16-elements with relaxed equilibrium condition match theory for both axial forces and bending moments in the fixed-fixed case. Less accurate, but still quite good results are obtained when equilibrium constraints are imposed.

O.tO CC G.

O. 18 . t8 -

00. IC I-

el

L(~I[ND 4,

a • "--

IS.HOOf DIIPl.figlCHENI ILl[HEM

HHtE-ILIHENI I | O U | L I I A I U N M T J I F | | O I H H t S ' | L | H | M t|OU|I.JIJtJUH RELAX|01 HHIN,YIJClY, $OLUT|OH

O, t q

kJJ O. t 2 E X:

O. tO . 4 ,

1:0.00 k,,e

0

isJ m

O. O0 "

4. Q

O.Oq •. I: W r,,4 O.Ot m t,,,,q

..J gc

0.00

O -O.U u'; z 16j -O, Oq 1C Omq

I:~ - 0 . O0

&

z

I0 Z - o . 08 0

, tO

, |O

"

, ......... tO

i qO

6 '

i SO

-

, . . . . . 80

RNGL( HEI:ISUflEO FflOl'l VERTICRL

, 70

(01

Fig, I I , Bending m o m e n t in a fixed-fixed curved beam.

"''

i

80

|

|0

W. Graf et al. / Three-dimensional thick shell elements 0

tO.O0

ac

8.oo =

"

I • • ,, ----

| . 00

371

I.[GENO 1E-NODE O|SPLIICF.HENI ELEIIEWT HI4LG-ELEHEWI I E O U | L I I I I i l M SRTlSF|EOl AIWL~rll CPL SOLUTION

7.00

_°-" ,.oo w

•

•

•

•

•

•

6.00

n-

O

.J "".

G::

11.003.00~.00

"

t o00 -

C:) 0,,1212

±

±

"

-"

•

•

LL --

"

!

A

AL

I-=,I

_j

-1,,00

Z

-re. O0 =

0

-

e,,,,o

Ir,lr'l -qJ,, O0 I,a,JI

:E: "q. O0 -

t,,,q

z o Z

-S.O0

- •

- I I , O0

, tl

l

00 , 10

O0 i ..... tO

i

10

QO

O0 k' El

"

, ' iO

ANGLE HEASUPIEO FROM V E R T I C A L

|

|

70

(0

II0

90

I

Fig. 12. A x i a l f o r c e i n a f i x e d - f r e e c u r v e d beam.

t.tS

Q E

t,t0

I m • ,, -.--.

,"

LEGEND t ~ . N O O [ 0JSPLIICITJ~ITT [ L [ H | W T HHLG-ELEHENI IEOU|L||RJUH SATJSFJE01 AN~Y5|¢~ 8OLUTXON

I=,, m ILl E ID Z

t,OS-

z

I.oo

:,

A

A

•

-

A

,,L

,-

,,

,,',

& &

&

O 2g

O " .J G: Z O "

0.80

•

•

O0

•

0(O

•

O0

Q

0

•

•

BO

•

•

0.U

I,IJ 2::

o

o, eo

Z

Z 0,75 0

!

!

|

t0

|0

20

RNGL[

HERSUREO

|

110

|

lle

i

|

|

li0

70

FROH V E R T I C R L

tO)

Fig. 13. B e n d i n g m o m e n t i n a f i x ~ l - f r e e c u r v e d b e a m .

,i

!

ll0

80

372

w. Graf et al. / Three.dimensionalthick shell elements

For the fixed-free boundary conditions, again accurate stress predictions are obtained as seen in Figs. 12 and 13. Comparison with the incompatible displacement model

It was mentioned earher, when the incompatible displacement model is employed, certain modifications have to be made in the evaluation of stiffness matrices in order for the element to be reliable when distorted. This includes the use of a 'constant' Jacobian for stiffness evaluation and discarding the load contribution corresponding to the imcompatible modes [10,32]. By doing so, it can generally be shown that the incompatible displacement model is equivalent to a certain class of hybrid/mixed elements [11,15,25].

E = 1200

B

M= 0.10

L=4.0 h=4.0

Pressure I o o d = 4 . 0

.k ' t X

j.J

A

A L

v I

r a

b Fig, 14. A 16-node element subjected to pressure load, (a) Regular geometry. (b) Distorted geometry.

W. G r a f et al. / Three-dimensional thick shell elements

Table 2 Displacements and stresses for the pressure loaded element of Fig. Case

uA ( - 0.oo6oo)

uB (-

0.01200)

373

14

Oxx a ( - 4.ooo)

O;T a (o.ooo)

~: a ( - 4.ooo)

1

- 0.00603

- 0.01189

- 3.973

0.OOl

- 4.077

2

- 0.00667

- 0.01340

- 4.423

- 0.025

- 4.031

3

-0.00600

-0.01200

-4.000

0.OO0

-4.000

4

- 0.00600

- 0.01200

- 4.000

0.OO0

- 4.000

5

- 0.00637

- 0.01274

- 4.229

- 0.026

- 4.060

6

- 0.00600

- 0.01200

- 4.000

0.OO0

- 4.000

-

0.00600

- 0.01200

- 4.000

0.OO0

- 4.000

- 0.00600

- 0.01200

- 4.000

0.OO0

- 4.000

HM16-element (undistorted) HM16-element (distorted)

a Computed at the element centroid. Table

2a

Case

Incompatible mode load included

Distorted

1

Yes

No

2

Yes

Yes No

3

No

4

No

No

5

No

6

No

Yes Yes

Jacobian Variable Variable Variable Constant Variable Constant

In order to investigate the effect of the modifications mentioned above on element performance, a single 16-node displacement-based solid element with the incompatible mode (1 - t 2) is considered in Fig. 14. The element is subjected to a pressure loading (P ffi 4.00 units) on its + r, + t, and - t faces as indicated in the figure. Distorted as well as regular geometries are considered in the analysis. Considering a 'consistent' formulation, in which contributions of the incompatible modes to the nodal load vector are included, it is observed that the results are in error for both undistorted and distorted geometry (Cases 1 and 2 of Table 2). Discarding this load contribution, an element with regular geometry gives exact results for constant as well as variable Jacobian (Cases 3 and 4). However, even when the internal degrees of freedom are not loaded, a distorted element with the incompatible mode's stiffness contribution evaluated using variable Jacobian (Case 5) still gives inaccurate results. For a valid solution in this case, a constant Jacobian (at element centroid) must therefore be used, as evident from the results for Case 6. For comparison, the solution using an HM16-element with relaxed equilibrium is given in Table 2a, where good results are obtained for both regular and distorted models. It is noted that the modification introduced in the formulation of incompatible mode displacement elements may not add significantly to the computational effort in linear analysis. However, this will have a significant effect on computations when the model is extended for solving nonlinear problems, and thus makes the model unattractive. Conclusion

A hybrid/mixed formulation of three.dimensional thick shell elements with 8, 16, and 18 nodes based on the Hellinger-Reissner variational principle is presented. It is important to

374

W. Graf et al. / Three-dimensional thick shell elements

note that the key element for developing a successful hybrid/mixed lies in the careful selection of stress interpolation functions. In this connection, three guidelines should be considered: (i) suppression of kinematic deformation modes, (ii) element property must be invariant with respect to the reference coordinates, and (iii) the constraint index of the element must be equal to or greater than one for constrained-media applications. The thick shell elements reported herein have been developed by following these guidelines. Several test problems were analyzed to determine the numerical performance of the aforementioned thick shell elements. Based on the results obtained, the following conclusions are reached: (1) The elements are free from either shear or membrane locking phenomenon even when the aspect ratio becomes fairly large (on the order of 2000). (2) The elements are relatively insensitive to geometric distortion as compared to their equivalent displacement-based elements. (3) Excellent stress predictions are obtainable. (4) The elements show good mesh convergence characteristics. (5) For the most part, the element behavior appears to be rather insensitive to the equilibrium constraint. However, in some cases the elements with relaxed equilibrium outperform those with the equilibrium constraint imposed in the natural coordinates. (6) Although the incompatible displacement model may give qualitatively equivalent behavior, the modifications introduced will make its nonlinear application unattractive due to added computational efforts. In conclusion, it is anticipated that the hybrid/mixed formulation approach has great potential for the development of various efficient finite elements, including plates and shells~ and their applications to nonlinear problems. Some of our research results concerning several of these developments will be reported in the related write-ups [8,9,27,28].

References [|] ASHWELL,D.G. and R.H. GALLAGHElt,editors, Finite Elementsfor Thin Shells and Curved Members, Wiley, New York, 1976. [2] B^'rHl~,K.J., Finite Element Procedures in Engineering Anal.vsis, Prentice.Hall, Englewood Cliffs, N J, 1982. [3] BELYTSCHKO,T., C.S. TSAYand W.K. LIu, "A stabilization matrix for the bilinear Mindlin rate element", Comput. Meths. AppL Mech. Engrg. 29, pp. 313-327, 1981. [4] BELYTSCHKO,T., W.K. LIu, J. I~NHI~DYand J.S. ONO, "Hourglass control for linear and nonlinear problems", Comput. Meths. AppL Mech. Engr$. 43, pp. 251-276, 1984. [5] BELYTSCHKO,T., W. LIU, J,S. ONO and D. LAM, "Implementation and application and a 9-node La&range shell element with spurious mode control", Computers & Structures, to appear. [6] BOLOURCm,S., "On finite element nonlinear analysis of general shell structures", Ph.D. Dissertation, Massachusetts Institute of Technology, 1979. [7] CHANG, T.Y. and S. PSACHUKTAM,"NFAP--A nonlinear finite element analysis program", Rept. No. SE76-3, The Univ. of Akron, Akron, Ohio, October 1976. [8] CHANG,T.Y. and A.F. SALE~n,A hybrid/mixed quadrilateral element for plane stress and plane strain analysis", In preparation. [9] CHANG, T.Y., A.F. SALEEnand Y.C. ZHANG, "On the hybrid/mixed formulation of isoparametric elements for nonlinear bending analysis", Proc. Invitational China-American Workshop on Finite Element Methods, Chengde, China, June 2-6, 1986. [10] COOK, R.D., Concepts and Applications of Finite Element Analysis, 2nd ed., Wiley, New York, 1981. [11] FxOIER, M., L. NXLSSONand A, SAMU~LSSOH,"The rectangular plane stress element by Turner, Plan and Wilson", lnternat. J. Numer. Meths. Engrg. 8, pp. 0,33-437, 1974. [12] HIHTGN, E., E.M. SALON~Nand N. BlCANXC,"A study in locking phenomena in isoparametric elements", Third MAFELAP Cot~, Brunel Univ., UxbHdge, 1978. [13] HUGHES,T.J.R. and M. COHEN,"The 'Heterosis' finite element for plate bending", Computers & Structures 9, pp. 445-450, 1978.

W. Graf el al. / Three-dimensional thick shell elements

375

[14] Huom's, T.J.R., M. COHENand M. HAROUN,"Reduced and selective integration techniques in the finite element analysis of plates", Nuclear Engrg. & Design 46, pp. 203-222, 1978. [151 IRONS, B.M., "An assumed stress version of the Wilson 8-node isoparametric brick", Rept. No. C N M E / C R / 5 6 , Center for Numerical Methods in Engineering, Dept. of Civil Engineering, Univ. of Wales, Swansen, U.K., May 1972. [16] IRONS, B.M. and S. AHMAD, Techniques of Finite Elements, Ellis Horwood, Ltd., Chichester, 1980. [17] MALKUS,D.S. and T.J.R. HUGHES,"Mixed finite element methods--reduced and selective integration techniques: A unification of concepts", Comput. Meths. Appi. Mech. Engrg. 15, pp. 63-81, 1978. [18] PLAN, T.H.H., "Derivation of element stiffness matrices by assumed stress distributions", AIAA Journal 2, pp. 1333-1336, 1964. [19] PLAN, T.H.H., "Recent advances in hybrid/mixed finite elements", Proc. Internat. Conf. of Finite Element Methods, Shanghai, China, August 2-6, 1982. [20] PLAN,T.H.H. and D.P. CHEN, "Alternative ways for formulation of hybrid stress elements", Internat. 3. Numer. Meths. Engrg. 18, pp. 1679-1684, 1982. [211 PLAN,T.H.H. and D.P. CHEN, "On the suppression of zero energy deformation modes", Internat. 3. Numer. Meths. Engrg. 19, pp. 1741-1752, 1983. [22] PLAN, T.H.H., D.P. CHEN and D. KANO, "A new formulation of hybrid/mixed finite element", Computers & Structures 16, pp. 81-87, 1983. [23] PLAN, T.H.H. and K. SUMmARA, "Rational approach for assumed stress finite elements", Internat. 3. Numer. Meths. Engrg. 20, pp. 1685-1695, 1984. [24] PLAN, T.H.H., K. SUMIHARAand D. KANO, "New variational formulations of hybrid stress elements", Proc. Nonlinear Structural Analysis Workshop, Cleveland, Ohio, April 19-20, 1983. [251 PLAN, T.H.H. and P. TONG, "Relations between incompatible displacement model and hybrid stress model", Internat. J. Numer. Meths. Engrg., to appear. [26] PUGH, E.D.L., E. HINTONand O.C. ZIENKtEWlCZ,"A study of quadrilateral plate bending elements with reduced integration", lnteruat. J. Numer. Meths. Engrg. 12, pp. 1059-1079, 1978. [271 SArEe, A.F. and T.Y. CHANG, "On the mixed formulation of curved beam elements", Comput. Meths. Appl. Mech. Engrg. 60, 1987, to appear. [2s] SALEEe, A.F., T.Y. CHANG and C. JIANO,"A hybrid/mixed quadrilateral element for plate bending analysis", In preparation. [29] STOLARSKk H. and T. BELYTSCHKO,"Membrane locking and redued integration for curved elements", J. Appi. Mech. 49, pp. 172-176, 1982. [3o1 WASHlZU, K., Variational Methods in Elasticity and Plasticity, 3rd ed., Pergamon Press, Oxford, 1982. 1311 WILSON,E.L., R.L. TAYLOR,W. DOHERTYand J. GHABOUSSk"Incompatible displacement models", in: Numerical and Computer Methods in Structural Mechanics, Academic Press, New York, 1973. I321 ZIENKIEWtCZ,O.C., The Finite Element Method, 3rd ed., McGraw-Hill, London, 1977. [331 ZIENKIEWiCZ,O,C., R,L. TA~LORand J.M, Too, "Reduced integration technique in general analysis of plates and shells", lnternat. J. Numer. Meths. Engrg. 3, pp. 275-290, 1971.