Piecewise hierarchical p-version curved shell element for geometrically nonlinear behavior of laminated composite plates and shells

Pergamon

0045-7949(94)00506-o

Computers & Strucrures Vol. 55, No. I, pp. 4766. 1995 Copyright Q 1995 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0045-7949/95 99.50 + 0.00

PIECEWISE HIERARCHICAL p-VERSION CURVED SHELL ELEMENT FOR GEOMETRICALLY NONLINEAR BEHAVIOR OF LAMINATED COMPOSITE PLATES AND SHELLS J. H. Liu and K. S. Suranat The University

of Kansas,

Department

of Mechanical KS 66045-2234,

Engineering, U.S.A.

(Received 7 December

3013 Learned

Hall, Lawrence,

1993)

Abstract-This paper presents a nine-node three-dimensional laminated composite curved shell finite element formulation for geometrically nonlinear (GNL) analysis of laminated plates and shells where the displacement approximation for the laminate is piecewise hierarchical and is derived based on p-version. The displacement approximation for the element is developed first by establishing a hierarchical displacement approximation for each lamina of the laminate and then by imposing interlamina continuity conditions of displacements at the interfaces between the laminas. The nodal variables for the entire laminate are derived from the nodal variables of the laminas and the interlamina continuity conditions of displacements. The element formulation ensures Co continuity of displacements across the interelement as well as interlamina boundaries. The lamina stiffness matrices and the equivalent nodal force vectors are derived using the principle of virtual work and the hierarchical displacement approximation for the laminas. Interlamina continuity conditions are used to construct the transformation matrices for the individual laminas which permit transformation of the lamina degrees of freedom to the laminate degrees of freedom. The interlamina behavior incorporated in this formulation is in total agreement with the physics of laminate behavior for composite plates and shells. In formulating the properties of the element, complete three-dimensional stresses and strains are considered, hence the element is equally effective for very thin as well as extremely thick laminated shells and plates. Incremental equations of equilibrium are derived and solved using standard Newton’s method. The total load is divided in increments and for each increment of load equilibrium iterations are performed until each component of the residuals and the generalized nodal displacement vector are within preset tolerances. Numerical examples are presented to show the accuracy, p-convergence characteristics and overall advantages of the present formulation.

nent to this approach is discussed in detail. The displacement-based GNL shell formulations can be grouped loosely into three categories based on the manner in which the rotations of the normals are used in the kinematics of deformation. In the first approach, the element displacement field is a linear function of the nodal rotations. The element displacement approximation is based on the small rotation concept where the cosine of the angle is approximated by one and the sine of the angle is taken to be the angle itself [23]. It is rather obvious that with such linearization a true Lagrangian formulation is not possible and many load steps are required to simulate even moderate rotations. The solutions obtained from such formulations often drift from the true solution when large rotations are present due to the inaccurate representation of the kinematics of deformation. In the second category of GNL shell deformations, the element displacement field is expressed in terms of nonlinear functions of the nodal rotations [34]. In this approach the kinematics of deformation is exact and permits true Lagrangian GNL formulation. However, the noncumulative nature of nodal rotations presents some difficulties. This

INTRODUCTION

Many

curved

shell

finite

element

formulations

have

been developed and reported in the literature. A selected number of textbooks and papers are given in Refs [l-13]. Of special interest here are geometrically nonlinear (GNL) curved shell formulations. A majority of the basic GNL formulations were developed in the late 1970s and early 1980s. These formulations can be grouped into two broad categories: displacement based formulations 114-521 and hybrid or mixed formulations [53-681. Generally speaking displacement based formulations ensure only Co continuity. Due to their simplicity such formulations are very appealing and are commonly used. On the other hand, the hybrid or mixed formulations provide interelement continuity of stresses, but with a significant increase in the formulation complexity. An excellent discussion on both of these approaches can be found in Ref. [68]. The present formulation is based on p-version piecewise hierarchical displacement approximation, hence only the literature pertit Author to whom all correspondence should be addressed, 47

48

J. H. Liu and K. S. Surana

formulation removes the small load step limitation present in the first approach. Other approaches such as Euler angles and Rodrigues parameters [40,41] have also been reported to describe nonlinear nodal rotation functions. The third approach is similar to the three-dimensional solid element GNL formulations. The element geometry and the kinematics of deformation are described by either the nodes located on middle and top surfaces or middle and bottom surfaces of the element [22]. Though these formulations eliminate nonlinear nodal rotation functions, they introduce more complications in the use of the element, especially for thin shells. The GNL shell formulations based on the displacement and their first derivatives with respect to inplane natural coordinates [36-38,491 also do not have nonlinear nodal rotation functions and work well for thin shells. Since the inception of the p-version concept [69], numerous papers have been written on the subject. Surana and Sorem [79,81] presented p-version curved shell element formulations for linear and GNL behavior of shells. The hierarchical element nodal variables and the corresponding approximation functions were derived from the Lagrange family of interpolation functions of order pt. pq and p<. This was accomplished by first establishing onedimensional hierarchical approximation functions and the corresponding nodal variable operators in the 5, q and [ directions and then taking their products. The resulting approximation functions and the nodal variables were hierarchical. The element displacement approximation ensured Co continuity. The element geometry was described by the coordinates of the nodes located on the middle surface and the nodal vectors describing the bottom and top surfaces of the element. The element properties were established using the principle of virtual work and the hierarchical element approximation. In formulating the element properties complete three-dimensional stresses and strains were considered. A brief literature review of laminate plates and shell theories and finite element formulation for linear and GNL behavior of laminated plates and shells is presented in this section. The classical laminate plate theory is well established and is based on the Kirchhoff hypothesis [82, 831. Some analytical solutions to the laminate problem can be found in the work published by Pagan0 [84-861, and Srinivas and Rao [87]. As is well known, the shear deformation effects in laminated composites are much more pronounced than in monolithic materials. The theories incorporating transverse shear deformation effects are either based on the assumed stress field [88-911 or based on assumed displacement field [91-931. The extension of the displacement based theories to laminated composites has been accomplished by expressing displacements as higher order functions of the thickness coordinate [94-96,98-1031. A summary of various isotropic plate bending theories can be found in a recent paper by Szabo and Sahrmann [104].

The first successful p-approximation formulation was reported by Brombolich and Gould [105], however the authors of this paper did not recognize the hierarchical nature of the approximation. The idea of using the second and higher order derivatives of shape functions as coefficients of the polynomial shape functions has been reported by Basu et al. [106], and Basu and Lamprecht [107]. Since 1969, a large number of papers have been published on the p-version subject, and its merits over h-approximation are well documented in the literature [108]. In a series of papers by Surana and Kalim [ 109, I IO] and Surana and Phillips [l 11, 1121, shell and transition element formulations were presented for heat conduction where the temperature and temperature gradient in the thickness direction were retained as nodal variables, permitting a linear temperature distribution in the direction of the shell thickness. In these developments the temperature field in the element thickness direction was hierarchical, but the development stopped at a linear temperature approximation (in the thickness direction) which corresponds to a p-level of one. The method utilized by Surana and Kalim and Surana and Phillips, which used Lagrange interpolating polynomials in the thickness direction to derive the approximation functions and the hierarchical variables, offers a direct approach for constructing the approximation functions as well as nodal variables. Utilizing this approach, Surana and Sorem [ 113, 1141 developed a three-dimensional curved shell element with a higher-order hierarchical displacement approximation in the thickness direction of the shell. This work was further extended by Surana and Sorem [115, 1161 for laminated composite plates and curved shells. In these formulations the order of approximation in the plane of the element remained constant. Such formulations with hierarchical higherorder transverse effects, though very effective for many thick as well as thin (isotropic or laminated) plate and shell applications, fail to perform optimally in three-dimensional situations with concentrations, singularities, etc. In such applications, a variable p-level control in all three directions (5, q and [) is essential for optimal performance. In a recent paper Woo and Basu [ 1171 presented a cylindrical shell element formulation based on p-version in the plane of the element. Surana and Sorem [ 1181 also presented a completely hierarchical p-version curved shell finite element formulation for laminated composite plates and shells. However, the drawback of Surana and Sorem’s formulation [ 1181 for laminated composite plates and shells is that the displacement field for the entire laminate is described by a single higher order polynomial. This results in inter-lamina continuity of displacements as well as their derivative of first and higher orders, which is not in agreement with the physics of laminate behavior. All higher order shear deformation theory based formulations suffer from the same drawback.

Behavior of laminated composite plates and shells The motivation for the work presented in this paper is to remedy the shortcomings of the higher order shear deformation (including p-version) laminate theory-based finite element formulation and develop a p-version formulation for GNL behavior of laminated composite plates and shells, in which the interlamina behavior of displacements, stresses and strains is in agreement with the physics of laminate behavior. The result is the p-version nine node curved shell element based on piecewise hierarchical displacement approximation using total Lagrangian approach. The displacement approximation for the entire laminate is piecewise hierarchical and is derived based on p-version displacement approximation for individual laminas. The approximation functions and the nodal variables for each lamina are derived directly from the Lagrange family of interpolation functions of order pc, p,, and ‘p<. This is accomplished by first establishing one-dimensional hierarchical approximation functions and the corresponding nodal variables operators in the 5, 9 and kc directions for the three, three and node one equivalent configurations that correspond to pe:+ 1, p,, + 1 and “pi + 1 equally spaced nodes in the 5, 4 and kc directions and then taking their products. The displacement approximation for the entire element is developed by the established hierarchical displacement approximation for each lamina and then by imposing interlamina continuity of displacements across the interfaces between the laminas. The nodal variables of the entire laminate are derived from the nodal variables of the laminas and the interlamina continuity conditions of displacements. The element formulation is in complete agreement with the physics of deformation in laminated composite curved shells and ensures Co continuity or smoothness of displacements across inter-element as well as interlamina boundaries. The incremental equations of equilibrium are established for a single lamina using hierarchical displacement approximation and total Lagrangian approach. The incremental equations of equilibrium for the entire laminate are developed using an approach similar to the linear static case. Incremental equations of equilibrium are solved using Newton’s method. The total load is divided into increments. For each increment of load equilibrium iterations are performed until each component of the residual is within a preset tolerance. Numerical examples are presented for GNL behavior to demonstrate the effectiveness, modelling convenience, accuracy and overall superiority of the present formulation over existing element formulations for laminated composites. Convergence graphs showing the rates at which the errors are decreasing are presented. The element formulation produces very accurate results for extremely thin as well as very thick plates and shells.

49

ELEMENT GEOMETRY AND MAPPINGS

Consider a laminated curved shell element composted of an arbitrary number of generally orthotropic laminas with arbitrary lay-up pattern bonded together perfectly at the interlamina boundaries. Figure 1 shows a map of the nine-node curved shell element in the 5, q and c natural coordinate system. The element nodes are located at the mid-surface ([ = 0) and the bottom and the top surfaces are described by 5 = - 1 and 5 = + 1. The Cartesian coordinate of a point P can be defined by the Cartesian coordinates of the nodes and the nodal vectors.

where jVi(& q) are the serendipity parabolic interpolation functions and Vi are nodal vectors. In order to facilitate the development of the element properties, two different natural coordinate systems are used. The first coordinate system is &l[, in which the entire element is mapped into a two-unit cube with the origin of the coordinate system located at the center of the cube. The & coordinates are in the plane of the element, whereas the [ coordinate is in the transverse or the lay-up direction of the laminate. The second type of natural coordinate system is the lamina natural coordinate system 5, qk[ in which the kth lamina is mapped into a two-unit cube. By imposing laminate coordinate continuity across the interlaminas in the entire laminate, the relationship between [ and k[ can be described by:

(-kh(l-“0+2&h>.

[=-I+;

(2)

,=1

The inverse of the relation written as:

defined by eqn (2) can be

(3) where t is the total thickness

t = f/h;

of the laminate

m is the total number

of laminas.

,=I _*_--r----___

(4)

-_

~~~

..’

lamha

laYuP

direction

Fig. 1. A nine-node

laminated p-version element.

curved shell

J. H. Liu and K. S. Surana

50

The thickness in each lamina or the entire laminate may vary from node to node. The approximation for the lamina thickness can be described by the following:

where kh,; i = 1, 2, . , 8 are the nodal thicknesses of the k th lamina at the locations corresponding to the locations of the lamina nodes. Using eqns (2) and (5), interlamina values of the [ coordinate can be obtained. OUTLINE

OF GENERAL

PROCEDURE

There are several main steps involved in deriving the stiffness properties of the laminated curved shell. The first step is to develop a hierarchical displacement approximation based on the p-version to describe the deformation in each lamina of the entire laminate. The second step is to derive the stiffness matrix and equivalent load vector for each individual lamina of the laminate. In the last step, we impose the displacement continuity across inter-lamina boundaries, and thereby each lamina’s stiffness matrix and equivalent load vector can be transformed and summed to obtain the laminate stiffness matrix and equivalent load vector. The displacement approximation in each lamina can be of arbitrary polynomial order pc, p, and “pr, where 5 and q are the coordinates in the shell element’s plane and k< is the coordinate in the thickness direction of the lamina. The hierarchical approximation is constructed in the lamina coordinate system (5, r~,‘[). To ensure displacement continuity across the inter-lamina boundaries, the order of the approximation functions in the 5, r~ directions has to be the same for all laminas of the laminate, but the order of the approximation can be different and arbitrary in the transverse direction (k[) for each lamina. First, using the geometric description of the entire laminate given by eqn (1) and eqns (2)-(5) the 5, q, [ coordinate of the lamina nodes are located at the mid-surface of each lamina. Next we establish nodal vectors kV, at these nodes of the lamina k. The lamina k is mapped into its own natural coordinate spaces (5, r~,“0 and thus an equation similar to eqn (1) can be used to describe the geometry of the kth lamina. Next we construct a hierarchical displacement approximation for lamina k in the 5, r~and k[ coordinate system. The hierarchical approximation functions and corresponding variables are derived directly from the Lagrange family of interpolation functions. The lamina stiffness matrices and equivalent load vectors are derived using the principle of virtual work. Next we construct interlamina continuity conditions, i.e. the displacements km’6(t, ‘1,k-1<) and k~({, 1, kc) are equated at (pc + l)(p, + 1) equally spaced points at

the interface between the (k - 1)th and kth laminas. This process ensures that displacements are continuous everywhere in the interlamina plane. Using the inter-lamina continuity equations, (pt + 1) (p, + 1) degrees of freedom are eliminated for each component of displacement. It is worth noting that all the degrees of freedom of the first lamina are retained in the laminate degrees of freedom. For the second and subsequent laminas, (pc + l)(p, + 1) degrees of freedom are eliminated using interlamina continuity conditions and the remaining degrees of freedom are added to the laminate degrees of freedom. This is accomplished by using lamina transformation matrices which transform lamina degrees of freedom to the laminate degrees of freedom. For individual lamina tangent stiffness matrices, equivalent load vectors are transformed to entire laminates by using these transformation matrices. The transformed lamina matrices and vectors are then summed to give the tangent stiffness matrix and the equivalent load vector for the entire laminate. The numerical values of the lamina matrices and equivalent load vectors are calculated using Gaussian quadrature with the pt + 1, pq + 1 and “pi + I rule for 5, r) and k< directions. LAMINA DISPLACEMENT

APPROXIMATION

The geometry for the kth lamina can be described in the lamina natural coordinate system

The detail of deriving the hierarchical approximation can be found in the previous paper by Surana and Sorem [118]. A brief outline of how these approximation functions are used here to develop the lamina displacement approximation is presented in the following: let

{ks}T = {“u, ku,“w}.

(7)

The pt + 1 equally spaced nodes in the 5 direction between - 1 and + 1 with their corresponding Lagrange interpolation functions and nodal variables can be reduced to an equivalent three-node configuration for which we can write:

ikW)) = [%Hk~j5=

-1

(9)

51

Behavior of laminated composite plates and shells

55-r

yPe -_a - (P<)!

56-l

7’7’“”

1

...,rq)2,{k6};].

(10)

if pc is even, a = 1, if pc is odd, a = 5. The p,, + 1 equally spaced nodes in r] direction between - 1 and + 1 with their corresponding Lagrange interpolation functions and nodal variables can be reduced to an equivalent three node configuration, for which we can write:

(16)

The eliminated variables {kS,,} in the q direction for the kth lamina are defined as:

We define the retained as:

degrees for the k th laminas

. ...(w)].

$-‘I -z--

q6-1

’

6!

““’

!C.z! (p,)!

1

(18)

(13)

if pV is even, a = 1, if ps is odd, a = ye. The “pi + 1 equally spaced nodes in the k[ direction between - 1 and + 1 with their corresponding Lagrange interpolation functions and nodal variables can be reduced to an equivalent one-node configuration for which we can write the following:

{%Vki)} = ,~O~(+J(~)~,=O.

The p-version approximation functions for a lamina (in eliminated and retained forms) can be written as (for one component of displacement): ikNel = W&lP’;,I>

(14)

W:,IP;A

W:,lWb,L

W:,lP’;el>

D’:,lI%el~

P':el W:el,

[N~,l[J$,lt

Wt,lP$l~

W%ItN;eI)

= ([N:l,

P’,zl> W,31> P’:l,

W,51, W:l>

Let [“rn] = [[“rn, 1, [“&I,

W:l,W:l,

11

[Xl)

(19)

ikA731,rA741, f + . , rkflkpc

and IkNrl = Ck&lWtl,

(15) y is the coordinate along the thickness direction such that -khj/2 ,< y Ckhj/2. Equations (8)-(15) provide the required onedimensional approximation functions and the nodal variable operators along 5, q and k[ directions, whose products give the required approximation functions and the corresponding nodal variables for the kth lamina. In order to facilitate the derivation, the nodal degrees of freedom and the corresponding approximation functions are classified as retained (with subscript r) and eliminated (with subscript e). The eleminated variables [kGS,]Tin the 5 direction or the k th lamina are defined as:

[“w, 1[N:l,

[kf121tN:l>

fk&;

fkfi,lW:l,

[kfi,l[N:l,

..

>

[k~21[N~l~[k~21[N,Zl~

.

> [‘r3,lW:l,

lW,21, Lk&

1W,31, .

. . , ~k8~,;lIN~l,

. , Lkfllri 1[N:l). (20)

The corresponding nodal variables can be also classified as eliminated and retained. From eqns (16) and (17), the eliminated variables can be obtained:

{kSe)T=

[{x6}:, (q$q ,{“s}:, (qg),, L

\

-

I-

7

\

1

J. H. Liu and K. S. Surana

52

where

RI PI WI PI lkO,jT WI PI ikQ? )’ [Ol {kQr)T I”@,. )’ {“Q,1’ Ik%IT WI PI {k?vIT i”0.Y 1’

{“KY (21)

[k4,1 wherem=2,3,4 ,..., pc;n=2,3,4 ,..., p,,. The retained variables are shown in eqn (18). The element approximation can be written in a more compact form. For example, for displacement k~ for the k th lamina we can write:

r

(28)

00000000 00010000 00000001 (29)

00001010 Similar expressions can be written for displacement v and w. These can be combined to give the hierarchical displacement approximations for the k th lamina:

1

01000100 10100000 The incremental equations of equilibrium lamina can be written as:

for the kth

0

N,

0

0

kN,

0

0

N,

0

0

kN,

0

0

0

N,

0

0

kNr

= [kNl(kG 1,

dikBj = [kG] d(‘6)

(30)

d{ktO} = [H][kG] d[ks] = [kBo] d{k6) (23)

d{kt’} = [iAO][kG] d{k6} = [kBL] dik8

‘>

(31)

thus

where

d{kt} = [[HI + [kA,][kG]] d{k6} = [kB] d{kS}.

jk6)= = [(“S:}‘, {kS:}T, tk6: }T,

tks:)‘?{ks:}T2 fks:‘)‘l. (24) LAMINA TANGENT STIFFNESS MATRICES EQUIVALENT LOAD VECTORS

The equation written

be

AND

of equilibrium

[kB]T[kcr] s

(33) The residuals {“ul) can be visualized as the nodal forces required to bring the assumed displacement pattern into nodal equilibrium. Taking the variation of (33) yields

_

(25)

kL’,

d{ke} = [kD] d{ke}.

kW,]

(26)

Substituting

(34)

(27)

(35)

eqn (35) into eqn (34) produces:

dIkY} = d[kR] -

Then ik6} = {“e”} + {“E”} = [H]lkB} + f[kA,]{kB}

[kB]Td[ka] dR s rn+

and

of Refs [lo, 341, if {%,}T = [ku,

d[kB]T[k8] dR s me

{kg}T=[kax,kug,k~~,kZ1.,,kZZ,,kZr~]

the notations

dR = IkR] - [kP].

0

d{kY} = d[kR] -

Following

for the kth lamina can

as:

j”Y } = [kl?] -

Assume complete three-dimensional state of stress and strain for each lamina and further assume the material to be linear elastic. In the present formulation we consider Piola Kirchhoff stresses and Green’s strains in the total Lagrangian formulation. In the global xyz coordinate system the stresses and strains for the kth lamina are:

(32)

d[kB]T{ka} dR s rn<, - [kK”“] d{kS},

(36)

Behavior of laminated composite plates and shells where [kpL] =

[“B]T[kD] [kB] dR. s kc%

(37)

53

The p-level in the 5 and q directions are the same for all laminas, whereas thep-level in the ‘[ direction can be different for each lamina. First consider the interlamina continuity of incremental displacement du between laminas 1 and 2, i.e.

Since [kB] = [kP] + [ke”] i=l,2,...

9(Pt + l)(P, + 1) (43)

we can write d[kB] = dIkBo] + d[kBL]. But d[kBo] = [0], since [kBo] is not a function hence d[kB]T(ka} do =

of tk6},

[kSL]T{%) dR.

where <,, q, and ‘*& are the 5, q and < coordinates of the (pc + l)(p, + 1) equally spaced grid of points at the interlamina boundary between laminas 1 and 2. Equation (43) can be expanded using eqn (24) and then can be rearranged in the following form:

(38)

s /me It can be shown [I 181 that

s

[kG]T[kS][kG] dQ

d[‘BL]{%} dR = Kle

=

Now define

1

[kLld{k~l

d{kS} and

d{*x”} =

(39) Equation

(44) can be written

as:

where kflXKl [kS] =

“7&l

k7.J31

“g41

k7g:[I,1

Symmetric and [I,] is a 3 x 3 identity

(40)

k%Kl I (46)

matrix

d{kY} = d{kRj - [[“f@] + [kK,,,]]d{k6} = d[kR] - [kK] d@ ),

where (41)

where [kK] is the desired symmetric tangent stiffness matrix for the kth lamina. Comparing eqn (41) with the incremental equations of equilibrium derived for formulations based on nonlinear nodal rotation functions, we note that, in the present formulation matrices, [K,,] and [&I are null. This is due to the fact that [kBo] is not a function of {k6}, hence d[kB”] = (0). Numerical values of the [kK] and {“R) are calculated using Gaussian quadrature with the pc + 1, p, + 1 and “p; + 1 rule. INTERLAMINA CONTINUITY CONDITION, LAMINA TRANSFORMATION MATRICES AND LAMINATE PROPERTIES

For each lamina, using the hierarchical displacement approximation defined by (24) and taking variation of both sides we can write:

(42)

I’m = [2Pel-‘[‘P,l = VI [‘R”]= [‘p,]-‘[lp,] and

[‘A] = -[*pJ’[‘p,]

. (47)

Note that d(‘x”} are the total laminate degrees of freedom for u for lamina one and d{2x”) are the total degrees of freedom for u for the laminate consisting of laminas 1 and 2. Thus, the interlamina continuity relations (37) between laminas 1 and 2 allow us to eliminate d{*6:} degrees of freedom from lamina 2. Obviously all degrees of freedom for lamina one are retained in the laminate degrees of freedom. An expression similar to eqn (46) can be derived using the interlamina continuity conditions of displacement du and dw between laminas 1 and 2:

(48)

.

(49)

J. H. Liu and K. S. Surana

54

Equations (46), (48) and (49) can be combined the following:

to give

[‘RI [Ol WI [*@I PI PI PI [‘RI WI PI [‘RI WI PI KY [‘a PI PI t2R1

using eqn (50) the transformation 2 can be derived as:

matrix

(50)

for lamina

[‘RI PI

PI L2Rl PI [‘RI PI PI [‘iI WI [Ol M [‘RI PI PI Pii] WI PI to1 VI LOI PI M PI PI VI PI PI [Ol [Ol M PI

dr26} =

-

(51)

where

or d{26} = [*7-l d(*x}.

The matrix [2T] transforms the degrees of freedom for lamina two to the laminate degrees of freedom. Now consider the interlamina continuity conditions of du, dv and dw between laminas 2 and 3 (first, consider the displacement u), i.e.

dC3u(5. 9 231,)) = 4’45 (1

I>

I

I)

r2m= L3Pel -‘L2Pel = VI [‘WI= [3p,]m’[2p,]and

(52)

v 23i,)) I)

/

>

[‘I?] = -[jp,]

‘[3p,] t (56)

Substituting obtain:

for d{26:} from eqn (46) into eqn (55) we d{ ‘x”}

1i

d{36:} = [[2R][‘R], [[2iF][‘ri] + [‘I?]], [‘I?]] dj2x:‘} d36:}

(53)

(57) where &, ni and 23c, are the c, yeand [ coordinates the (pc + l)(p, + 1) equally spaced grid of points the interlamina boundary between lamina 2 and Equation (53) can be expanded using eqn (42) and can write:

of at 3. we

or (58) Similarly the interlamina continuity of da and dw between laminas 2 and 3 would yield:

Then (54) can be written

dOX1 = [[2Rl>Pffll

(59)

4’S: } = U2Rl,[‘ill

(60)

as:

d{36:} = [[‘RI, i2R], [3g]1

(55) Equations

d{*&} =

i2Rl

PI

101

12Rl

[Ol

PI

(58) (59) and (60) can be combined

t’fil 101 PI WI 13fil WI L2Rl PI WI LidI

to give

PI

101

(61)

Behavior of laminated composite plates and shells

55

where m is the total number of laminas in the laminate. For each lamina eqn (41) can be written as:

Now define

[‘K] d(‘6} - d(jR} = d(/Y}; Substituting multiplying

(62)

j = 1,2,

(68)

. . , m.

from eqn (62) into eqn (68) and then the entire equation by [jr]’ we obtain:

[‘T]T[‘K][‘T] d(‘x} - [‘Z-ITd{‘R} = [‘Z-ITd{“P}; Then using eqn (61) the transformation lamina 3 can be derived as:

d{36} =

-

matrix

for j=l,2

PI WI M PI

WI PI PI PI

12N PI [Ol PI

M VI PI PI

PI PI VI PI

,...,m,

(69)

[‘m ml PI [II -

(63)

where

or d(‘6) = [37-]d{3a}.

(64) \X

The matrix [37’J transforms degrees of freedom for lamina 3 to the degrees of freedom for the laminate consisting of laminas 1, 2, 3. Similar transformation matrices can be constructed for lamina 4, 5, . . , etc. In general, between the (k - 1)51and kth laminas the following can be written:

[“-‘RI

PI PI

d{k6,} =

-["-'RI

d{36) =

-

PI Dl PI M PI

represents the total degrees of freedom for the laminate. The matrices [‘TIT [‘K] [‘T] and the vector [IT] (‘R} provide the nonzero contribution for each lamina to the laminate tangent stiffness matrix and the equivalent nodal force vector. Thus the entire laminate can be symbolically written as

[“RI PI PI PI PI Ikm’Rl LOI PI [“Xl PI PI L”-‘RI WI WI [“fil PI PI lk-‘Rl PI Lkm ‘RI PI PI PI [Ol PI PI PI

and

,,

I

[“iI WI PI PI [“Al PI PI PI fk& VI [Ol PI PI VI PI PI PI VI (only the shown);

nonzero

(65)

contribution

of

each

lamina

[km’R]=[[k-‘~][k-2R],[[k-‘jTl[k-’i?]+[k--1R’]]], (66) thus

[‘~lTP~lWI] d{d) - ,$, WITd?R } m = ,;, [‘TIT d(“Y 1,

d(j8) = (‘T] d(,x}; CAS 5511-E

,$,

j = 2,3, .

,m

(67)

(70)

J.

56

H. Liu and K. S. Surana NUMERICAL

where d{6} = d{“m}. At this point several comments note that

[‘P,l= t*P,l=

(71) can be made. First

. = [mP,l.

(72)

Second, if each lamina has constant thickness (but the thickness may vary from lamina to lamina), that is if kh is not a function of r and r~, then it can be shown that the inverse of [‘p,] is not necessary (details omitted for the sake of brevity). This can be confirmed by closely examining the approximation function matrices defined earlier. Also note that all degrees of freedom for lamina one are retained in the laminate degrees of freedom, whereas for lamina 2, 3, . m only dp6,);j = 2, 3, . , m degrees of freedom contribute to the laminate degrees of freedom and d{‘S, j; j = 2, 3, . , m degrees of freedom are eliminated using interlamina continuity relations of displacements. Interlamina continuity equations are conveniently arranged in the form given by eqn (67) which permits transformation of lamina degrees of freedom to the lamina degrees of freedom easily. For each lamina the material properties can be defined in a material coordinate system which may be different for each lamina and may also vary in the 5 and q directions for each lamina. The details of transforming [kD] from material to the global coordinate system are omitted for the sake of brevity. SOLUTION

METHOD

Incremental equations of equilibrium resulting from the GNL formulation are solved using Newton’s method, the details of which are well documented in the literature. Many different convergence criteria have been proposed in the various papers. The criteria based on “norms” form a good test for overall global convergence of the solution. However, such criteria are less sensitive to local divergence or lack of convergence. In the present work the accuracy of the solution during equilibrium iterations is judged based on the following: I{Y,}l
1
(73)

where N* is the total number of degrees of freedom. Thus, based on eqn (73) the solution at each load step is considered converged when each component’s total residual is within a specified tolerance A. This criteria ensures the desired accuracy for the total residual corresponding to each nodal degree of freedom. The convergence of the incremental generalized nodal variables can also be used. However, for high p-levels, the higher order derivatives may fluctuate considerably without causing any significant change in the element displacement field.

EXAMPLES

Three numerical examples are presented in this section to demonstrate the accuracy and convergence characteristics of the present laminated composite shell element formulation for geometrically nonlinear behavior. For all three examples, linear analysis is performed first to determine the required mesh and the p-levels for the convergence of the quadratic functional. A complete GNL study is then performed using this mesh and these p-levels. Results are also computed for one order higher p-levels to ensure that the p-levels determined from the linear analysis indeed are adequate for GNL study. During the implementation or during the processing of data for the examples, conversion of units is not done. Hence any set of consistent units, either in FPS or metric systems, can be selected. For this reason units are not given for the examples presented in this paper. Computations were performed on an Apollo DN3500 computer using double precision arithmetic. Example

1: clamped

sandwich

square

plate

(thick-

ness = 1.03)

A clamped square laminated plate consisting of three orthotropic laminas subjected to a uniformly distributed conservative load p is considered. Figure 2 shows the material properties, plate dimensions and loading for this example. Because of the symmetry considerations, only a quarter of the entire plate needs to be modelled. Figure 3 shows a finite element model of the plate using a single p-version curved shell element. First, we conduct a linear static p-convergence study using p = 9.23076 to determine the necessary p-levels in the 5, n and ‘[ directions for convergence of the quadratic functional. Graph of the quadratic function vs degrees of freedom for various values of p:, pa and “pr are shown in Fig. 4. From Fig. 4 we observe that pc = 5, pV = 5 and kp; = 1 are the minimum p-levels needed for the convergence of the quadratic functional. The nonlinear load-deflection behavior of the plate is generated using these p-levels with force control. The total load of 92.3076 is applied in 10 equal increments of 9.23076. At each increment of load, equilibrium iterations are performed until each component of the residual is within a preset tolerance of 10m4. From Fig. 9.4, we observe that the results of one element model with p: = 5, p, = 5, “p; = 1 and with p: = 6, p, = 6, kpC= 2 are virtually identical. The results for a one element model with pc = 4, p,, = 4, “p: = 1 and pc = 3, pV = 3, “pC= 1 are also shown in Fig. 5. The results obtained from the present formulation are compared with those reported by Sorem 11191, Chang and Sawamiphakdi [120] and Schmidt and Monforton [121]. Further increase in the p-levels from pr = 5, p7 = 5, “p; = 1 and the further refinement in discretization had virtually no effect on the results.

Behavior of laminated composite plates and shells

57

,lamina

3

lemina 2

I edges

Thickness:

1 & 3 : 0.015

Lamina

Material

, Lamina 2 : 1.0

properties:

Lamina

1 & 3 :

Lamina

2

:

E = 10.5 x lo6

E, = E, = 0.0 G,,

v = 0.3 E, = 0.1 x lo6

= G,, = 0.05 x lo6

v = 0.3 G,,

= 0.0

Fig. 2. Three lamina clamped square plate (example 1). Example 2: clamped three Iamina shallow arch Here we consider a clamped shallow arch consisting of three laminas under a conservative distributed line load P (Fig. 6). Figure 6 shows the material properties, dimensions and loading of this example. Because of the symmetry considerations, only half of the arch needs to be modelled. Figure 7 shows a three element model of the arch. Y

I----I

line of symmetry-

_- . j---lx I i

Example 3: clamped sandwich square plate (thickness = 8.0)

line of symmetry Fig. 3. One element model of the quarter

First, we conduct a linear static p-convergence study to determine the necessary p-levels in the 5, q and k< directions. Graphs of the quadratic functional vs degrees of freedom for various values ps, p,, and “p; are shown in Fig. 8. From Fig. 8 we observe that pc = 5, p,, = 2, “pr = 1 are the minimum p-levels needed for the convergence of the quadratic functional. The nonlinear load-deflection behavior of the arch is generated using these p-levels with displacement control. The total central deflection of 2.0 is applied in 20 equal increments of 0.1. At each increment of central deflection equilibrium iterations are performed until each component of the residual is within a preset tolerance of 10e4. From Fig. 9 we observe that the results for p-levels of pc = 5, p, = 2, ‘pC= 1 and those for p-levels of pr: = 6, p,, = 3, “pi = 2 are virtually the same. The results for pc = 3, pa = 2, k pr = 1, and pc = 2, p,, = 2, “pi = 1 are also shown in Fig. 9. Further increase in the p-levels from pt = 5, pv = 2, “pr = 1 or refinement in the discretization beyond the one element model had virtually no effect on the results.

plate (example

1).

A thick laminated square plate consisting of two laminas under a uniformly distributed conservative

J. H. Liu and K. S. Surana

58

375 Degrees Fig. 4. p-Convergence

750 of Freedom

of the quadratic

load p is considered in this example. Figure 10 shows the material properties, plate dimensions and loading. Because of the symmetry considerations, only a quarter of the plate need to be modelled. Figure 11 shows a one element model of the quarter plate. The purpose of this example is to demonstrate that the present element formulation is equally effective for thick laminated composite plates and shells also. First, we conduct a linear static p-convergence study using p = 0.1 to determine the necessary plevels in the 5, q and kc directions for the convergence of the quadratic functional. The graphs of the

3.5

functional

1500

1125 (example

I).

quadratic function vs degrees of freedom for various values p:, pq and “pc are shown in Fig. 12. From Fig. 12 we observe that pt = 5, ps = 5, kpi = 3 are the minimum p-levels needed for the convergence of the quadratic functional for linear static behavior. Since the plate in this case is much thicker than in example 1, the necessary p-level in “pi direction for the convergence of the quadratic functional is higher as well. The non-linear load-deflection behavior of the plate is simulated using these p-levels with force control. The total load of p = 0.4 is applied in 10 equal increments of 0.04. At each increment of load equilibrium iterations are performed until each component of the residual is within a preset tolerance of 10-m4.Figure 13 presents the load-deflection plots for p< = 4, p, = 4, kpc= 3; pt = 4, ps = 4, “pi = 2; p< = 3, p,, = 3, “pi = 2; p: = 3, p,, = 3, “pi = 1; and pt = 2, pq = 2, kpr = 2 also. The converged non-linear solution is obtained for p-levels of pt. = 5, p,, = 5, kpr = 3 and further increase in the p-levels beyond these had virtually no effect on the results. SUMMARY AND CONCLUSIONS

A Schmit [1211

0.5

J

0.0

0

25

50

75

100

load P Fig.

5.

Load-deflection

curve for the laminated (example 1).

plate

A nine-node three-dimensional laminated composite curved shell finite element formulation has been presented for geometrically nonlinear analysis of laminated plates and shells where the displacement approximation for the laminate is piecewise hierarchical and is derived based on p-version. The displacement approximation for the element is developed first by establishing a hierarchical displacement approximation for each lamina of the laminate and then by imposing interlamina continuity conditions of displacements at the interface between the laminas. The approximation functions and the nodal variables for each lamina are derived directly from the Lagrange family of interpolation functions of order pt, p,, and “pi. This is accomplished by first establishing onedimensional hierarchical approximation functions

Lamina

(core ): E= l,OOO.OOO. v=O.2

v=O.2 thick=0.0625

thick=0.0625

6 = 7.33729’

Fig. 6. A three lamina clamped arch (example 2).

2

Lamina 1 & 3 (faces): E=10,000,000.

width =1 .O, total t= 0.1875,

Central radius R = 133.114 I lihe of symmetry

total load P

Fig. 7. A three element model of the arch (example 2).

Clomped

J. H. Liu and K. S. Surana

60

I

-------------

----

N-W-

II II II

II

0’ IF cr’ 2 * .* _=_* . Nmt4@4N

II

II

II

II

----____

Behavior of laminated composite plates and shells

61

lamina 2

P

,i: _I

Illlllllllllllll~~

I edges

Thickness: Material

I : 4.0

Lamina

Lamina

E, = E, = 1.1

vly =vrz =vyz = 0.25

Gxz = G,, = 0.66 Fig. 10. A two lamina

G,, - 0.55

thick square

and the corresponding nodal variable operators in the 5, rl and k[ directions for the three and one node equivalent configurations that correspond to pc + 1, p,, + 1 and “pi + 1 equally spaced nodes in the 5, r] and “[ directions and then taking their products. The nodal variables for the entire laminate are derived from the nodal variables of the laminas and the Y

line of symmetry- _- --x

line of sy-mmetry

11. A one element

4s”/-45O

Properties:

E,= 3.3

Fig.

2 : 4.0

model of (example 3).

the

quarter

plate

clamped

plate (example

3).

interlamina continuity conditions of displacements. The element formulation ensures Co continuity of displacements across the interelement as well as interlamina boundaries. The lamina stiffness matrices and the equivalent nodal force vectors are derived using the principle of virtual work and the hierarchical displacement approximation for the laminas. Interlamina continuity conditions are used to construct the transformation matrices for the individual laminas which permit transformation of the lamina degrees of freedom to the laminate degrees of freedom. There is no restriction on the number of laminas and their lay-up pattern. Each layer can be generally orthotropic. The material directions and the layer thickness may vary from point to point within each lamina. The geometry of the element is defined by the nodes located at the middle surface of the element and the lamina thicknesses. The formulation presented here removes virtually all of the drawbacks present in the existing GNL laminated curved shell formulations and has many additional benefits. The element displacement approximations in the existing formulations are either based on linearized (with respect to nodal rotation) displacement field in which case a true Lagrangian formulation is not possible and the load step size is

62

J. H. Liu and K. S. Surana

example 3 Q p=pc=pq=2-6, . p=p,=p,,=2-6, Q p=pr=ps=2-6, * p=pC=p1=2-6,

used

I

I

-1.8

300

0

Fig. 12. p-Convergence

I

5.4

3.6 c .P ‘;; t 8 E E s 1.8

Fig.

13. Load-deflection

3.0

2.0 Load

P

4.0

(x100)

curve for (example 3).

the

clamped

I 900

I

I 1200

I

I 1500

of Freedom

of the quadratic functional (example 3).

severely limited, or is based on the nonlinear nodal rotation function approach, in which case the kinematics of deformation is exact but additional complications arise due to the noncumulative nature of nonlinear nodal rotation functions. Such limitations and difficulties do not exist in the present formulation. The hierarchical displacement approximation used here does not involve traditional nodal rotations that have been used in the existing shell formulations, thus the difficulties associated with their use are not present in this formulation. p-Version formulation

1.o

for GNL

I

I 600 Degrees

‘p<=4 ‘p‘=3 ‘pC=2 ‘pc=l

plate

permits the selection of the desired interpolation order. In formulating the properties of the element complete three-dimensional stresses and strains are considered, hence the element is equally effective for very thin as well as extremely thick laminated shells and plates. Retaining true derivatives (i.e. derivative with respect to depth coordinate y as opposed to 5) in the 5 direction as nodal degrees of freedom enables modeling of junctions that exhibit step change in the shell thickness. Since the nodal vectors do not necessarily have to be normal to the middle surface of the beam, the element formulation permits nondifferentiable geometries (sharp junctions) in modeling. In the vicinity of such junctions the stress behavior is normally complicated and often stress risers also exist. Such situations present no difficulties in the present formulation. The interlamina behavior incorporated in this formulation is in complete agreement with the physics of laminate behavior for plates and shells which the higher order shear deformation theories are incapable of incorporating. Incremental equations of equilibrium are derived and solved using a standard Newton-Raphson method. The total load is divided in increments and, for each increment of load equilibrium, iterations are performed until each component of the residuals and the generalized nodal displacement vector are within preset tolerances. Numerical examples are presented to show the accuracy, efficiency and overall advantages of the present formulation. The formulation is insensitive of the load step size (within the limitations of Newton’s method). In all of the numerical examples presented here, first the linear static studies were conducted to establish the required p-levels in r, q and k[ directions for the convergence of the quadratic functional, then these p-levels were used (and kept fixed) to conduct GNL studies. Numerical studies show that any further increase in the p-levels

Behavior

of laminated

composite

(from these converged values for a linear static case) had virtually no effect on the results. All three examples employ very simple finite element models that yield converged solutions. Acknowledgements-The computing facilities provided by the Computational Mechanics Laboratory (CML) of the Department of Mechanical Engineering are gratefully acknowledged. The financial support provided to the first author by the Department of Mechanical Engineering through Graduate Teaching Assistantship is also acknowledged. This paper in whole or part has been accepted for presentation at the ETCE Symposium on Materials, Design and Analysis, New Orleans, LA, 1994. REFERENCES 1. R. H. Gallagher, Finite Element Analysis Fundamentals. Prentice-Hall, Englewood Cliffs, NJ (1973). 2. R. H. Gallagher, Shell Elements. In Proc. World Congr. an Finite Elements in Structural Mechanics, Vol. 1, Bournemouth (1975). 3. P. G. Bergan et al. (Eds), Finite Elements in Nonlinear Mechanics, Vol. 1. Tapir Publishers, Trondheim (1978). 4. T. J. R. Hughes, A. Pifko and A. Jay (Eds), Nonlinear Finite Element Analysis of Plates and Shells, AMDVol. 48. ASME, New York (1981). 5. R. D. Cook, Concepts and Applications of Finite Element Analvsis. 2nd edn. Wilev. New York (1981). 6. K. J. Bathe, Finite Element Procedures in Engineering Analysis, 2nd edn. Wiley, New York (1981). 7. T. J. R. Hughes and E. Hinton (Eds), Finite Element Methods for Plates and Shell Structures, Vol. 1, Element Technology and Vol. 2, Formulations and Algorithms. Pineridge Press, Swansea (1986). 8. T. Belytschko, A review of recent developments in plate and shell elements. In Computational Mechanics-Advances and Trends (Edited by A. K. Noor), AMD-Vol. 75.._. DD. 217-231. ASME. New York (1986). 9. T. J. R. Hughes, The Finite Element Method. PrenticeHall, Englewood Cliffs, NJ (1987). 10. 0. C. Zienkiewicz, The Finite Element Method, 2nd edn. McGraw-Hill, London (1977). 11. K. S. Surana, Theory of curved shell element, Technical Report. McDonnell Douglas Automation Company, Saint Louis, MI (1978). 12. S. Ahmad, B. M. Irons and 0. C. Zienkiewitcz, Analysis of thick and thin shell structures by curved elements. Int. J. numer. Meth. Engng 2, 4199451 (1970). 13. 0. C. Zienkiewicz, J. Too and R. L. Taylor, Reduced integration technique in general analysis of plates and shells. Int. J. numer. Meth. Engng 31 375-390 (1971). 14. A. B. Sabir and A. C. Lock, The apuhcation of finite __ elements to the large deflection geometrically nonlinear behavior of cylindrical shells. In Variational Methods in Engineering (Edited by C. A. Brebia and H. Tottenham), pp. 7/66-7/75. Southampton University Press Southhampton (1973). 15. J. H. Argyris and P. C. Dunn, Finite element applications to nonlinear and non-positive definite linear problems. U.S.-Germany Symp., M.I.T., Boston, MA (1976). 16. J. L. Batoz, A. Chattopadhvav and G. S. Dhatt. Finite element large deflection analysis of shallow shells. Int. J. numer. Meth. Engng 10, 39-58 (1976). 17. J. H. Argyris, P. C. Dunne, G. A. Malejannakis and E. Shelke, A simple triangular facet shell element with applications to linear and non-linear equilibrium and elastic stability problems. Comput. Meth. appl. mech. Engng 10, 371-403; 11, 97-131 (1977).

plates

and shells

63

18. E. Ramm, A plate/shell element for large deflections and rotations. In Formulations and Computational Algorithms in Finite Element Analysis (Edited by K. J. Bathe, J. T. Oden and W. Wunderlich) M.I.T. Press, Cambridge, MA (1977). and 19. J. H. Argyris, P. C. Dunne, G. A. Malejannakis Schelkle, A simple triangular facet shell element with applications to linear and non-linear equilibrium and stability problems. Comput. Meth. appl. Mech. Engng 10, 371-402; 11, 977131 (1977). 20. J. H. Argyris, P. C. Dunne and D. W. Scharpf, On large displacement-small strain analysis of structures with rotational degrees of freedom. Comput. Meth. appl. mech. Engng 14, 401-145; 15, 99-135 (1978). nonlinear analysis of shells. 21. H: Parisch, Geometrical Comma. Meth. awl. mech. Ennna 14. 159-178 (1978). R. L.- Taylor and T.‘ J. R. 22. W. Kanok-Nukuichai, Hughes, A large deformation formulation for shell analysis by the finite element method. Comput. Struct. 13, 19-27 (1981). A geometric and 23. K. J. Bathe and S. Bolourchi, material nonlinear plate and shell element. Comput. Struct. 11, 23-48 (1980). Large defor24. T. Y. Chang and K. Sawamiphakdi, mation analysis of laminated shells by finite element method. In Computational Methods in Nonlinear Strut tural and Solid Mechanics, Washington, DC (1980). behaviour of 25. A. Pica and R. D. Wood, Postbuckling plates and shells using a Mindlin shallow shell formulation. Comput. Sfruct. 12, 759-768 (1980). MAGNA: a finite element system 26. T. A. Bruockman, for three-dimensional nonlinear static and dynamics structural analysis. Comput. Struct. 13, 4155423 (1981). Nonlinear finite 27. J. H. Argyris and Sp. Symeonides, element analysis of elastic systems under nonconservative loading-natural formulation. Part I, quasistatic problems. Comput. Meth. appl. mech. Engng 26, 75-123; 373-383 (1981). A formulation of 28. K. J. Bathe and E. N. Dvorkin, general shell elements-the use of mixed interpolation of tensiorial components. Int. J. numer. meth. Engng 22, 697-722 (1986). finite 29. T. J. R. Hughes and W. K. Liu, Nonlinear element analysis of shells. Part I-three dimensional shells. Comput. Meth. appl. mech. Engng 26, 331-362 (1981). and 30. R. A. F. Martins and D. R. J. Owen, Elastoplastic geometrically nonlinear thin shell analysis by the semiloof element. Comput. Struct. 13, 5055513 (1981). of shells including 31. H. Parisch, Large displacements material nonlinearities. Comput. Meth. appl. mech. Engng 27, 183-214 (1981). finite 32. E. Ramm and H. Stegmuller, The displacement element method in nonlinear buckling analysis of shells. In Buckling of Shells (Edited by E. Ramm). Springer, Berlin (1982). 33. T. J. R. Hughes and E. Carnoy, Nonlinear finite element shell formulation accounting for large membrane strains. Comput. Meth. appt. mech. Engng 39, 69-82 (1983). nonlinear formulation for 34. K. S. Surana, Geometrically curve shell elements. Int. J. numer. Meth. Engng 19, 581-615 (1983). 35. J. Oliver and E. Onate, A total Lagrangian formulation for the geometrically nonlinear analysis of structures using finite elements. Part I, two-dimensional problems: shell and plate structures. Int. J. numer. Meth. Engng 20, 2253-2281 (1984). 36. T. Y. Yang and C. J. Moore, Geometrically nonlinear formulation of a 48 d.o.f. quadrilateral shell element with rational b-spline geometry. Int. J. numer. Meth. Engng 21, 317-328 (1985).

64

J. H. Liu and K. S. Surana

37. T. Y. Yang and S. Saigal, A curved quadrilateral element for static analysis of shells with geometric and material nonlinearities. Int. J. numer. Meth. Engng 21, 617-635 (1985). large 38. S. Saigal and T. Y. Yang, Elastic-viscoplastic, deformation formulations of a curved quadrilateral thin shell element. Inr. J. namer. Mefh. Engng 21, 1897-1909 (1985). 39. M. K. Nygard, J. Kjeken and P. G. Bergen, Nonlinear shell analysis based on free formulation elements. Proc. Int. Conf. Finite Elements in Computational Methods, p. 61. Pergamon Press, Oxford (1985). Large rotation in problems of 40. J. F. Besseling, structural mechanics in finite elements methods for nonlinear problems. In Proc. EuropeeU.S. Symp., The Norwegian Institute of Technology, Trondheim, pp. 25139. Springer, Berlin (1985). _41. M. Geradin. G. Robert and P. Buchet, Kinematic and dynamics analysis of mechanisms. A finite element approach based on Euler parameters, in finite element methods for nonlinear problems. In Proc. Europe-U.S. Symp., The Norwegian Institute of Technology, Trondheim, pp. 41-60. Springer, Berlin (1985). Nonlinear shell 42. P. G. Bergan and M. K. Nygard, analysis using free-formulation finite elements. In Finite Element Methods for Nonlinear Problems (Edited by P. G. Bergan et al.). Springer, Berlin (1985). and H. Stegmuller, Ulti43. E. Ramm, K. Schweizerhof mate load analysis of thin shells under pressure loads. In Finite Element Methods for Nonlinear Problems. Springer, Berlin (1986). ” 44. H. C. Huanr and E. Hinton, Elastic-elastic and geometrically-nonlinear analysis of plates and shells using new nine node elements. In Finite Element Methods for Nonlinear Problems. Springer, Berlin (1986). A formulation of 45. K. J. Bathe and K. N. Dvorkin, general shell elements-the use of mixed interpolation of tensorial components. Int. J. numer. Mefh. Engng 22, 679-722 (1986). Resultant 46. W. K. Liu, E. S. Lam and T. Belytschko, stress degenerated shell element. Comput. Meth. appl. mech. Engng 55, 259-300 (1986). Large deformation 47. E. Ramm and A. Matzenmiller, shell analysis based on the degeneration concept. In Finite Element Methods for Plates and Shell Structures, Vol. 2: Formulations and Algorithms (Edited by T. J. R. Hughes and E. Hinton), Chap. 15, pp. 3655393. Pineridge Press, Swansea (1986). 48. G. M. Stanley, K. C. Park and T. J. R. Hughes, Continuum-based resultant shell elements. In Finite Element Methods.for Plates and Shell Structures-Vol. 2: Formulations and Algorithms (Edited by T. J. R. Hughes and E. Hinton), Chap. 1, pp. 1-45. Pineridge Press, Swansea (1986). 49. H. T. Y. Yang and W. C. Wu, A geometrically non-linear tensorial formulation of a skewed quadrilateral thin shell element. Int. J. numer. Meth. Engng 28, 285552875 (1989). analy50. K. M. Hsiao and H. C. Hung, Large-deflection sis of shell structures by using corotational total Lanraneian formulation. Comuul. Meth. auul. ._ Mech. En&g 73, 2099225 (1989). . 51. D. W. White and J. F. Abel. Accurate and efficient non-linear formulation of a nine-node shell element with spurious mode control. Comput. Strucf. 35, 621-641 (1990). Total Lagrangian 52. B. Chang and A. A. Shabana, formulation for the large displacement analysis of rectangular plates. Ini. J. numer. Mefh. Engng 29, 73-103 (1990). 53. T. Matsui and 0. Matsuoka, New finite element

54.

55.

56.

57.

58.

59.

60.

61.

62.

63.

64.

65.

66.

67.

68.

69.

70.

71.

scheme for instability analysis of thin shells. Znt. J. numer. Merh. Engng 10, 145-170 (1976). A. K. Noor and S. J. Hartley, Nonlinear shell analysis via mixed isoparametric element. Comput. Struct. 7, 6155626 (1977). C. Tahiana and L. Lachance, Linear and nonlinear analysis of thin shallow shells by mixed finite elements. Comput. Struct. 5, 167-177 (1975). R. L. Boland and T. H. H. Pian, Large deflection analysis of thin elastic structures by the assumed stress hybrid finite element method. Comput. Sfruct. 1, l-12 (1977). G. Horrigmoe, Hybrid stress finite element model for nonlinear shell problems. Inf. J. numer. Meth. Engng 12, 1819-1839 (1978). G. Horrigmoe and P. G. Bergan, Nonlinear analysis of free form shells by flat finite elements. Comput. Meth. appl. mech. Engng 16, 11-35 (1978). N. Fukucki and S. N. Atluri, Finite deformation analysis of shells: a complementary energy hybrid approach. In Nonlinear Finite Element Analysis of Plates and Shells (Edited by T. J. R. Hughes et al.), AMD-Vol. 48, pp. 2233247. ASME, New York (1981). A. K. Noor and C. M. Anderson, Mixed models and reduced/selective integration displacement models for nonlinear shell analysis. Znt. J. Mefh. Engng 18, 142991454 (1982). S. N. Atluri, Alternate stress and conjugate strain measures, and mixed variational formulations involving rigid rotations, for computational analysis of finitely deformed solids with applications to plates and shells,. I: theory. Comput. St&t. 18, 93-116 (1984). L. Recke and W. Wunderlich. Rotations as orimarv unknowns in the nonlinear theory of sheils and corresponding finite element models. In Finite Rota/ions in Structural Mechanics (Edited by C. A. Brebbia and S. A. Orszag), pp. 239-258. Springer, Berlin (1985). A. F. Saleeb and T. Y. Chang, An efficient quadrilateral element for plate bending analysis. Znt. J. numer. Meth. Engng 24, 1123-1155 (1987). A. F. Saleeb, T. Y. Chang and W. Graf, A quadrilateral shell element using a mixed formulation. Comput. Struct. 26, 7877803 (1987). J. J. Thiu and S. W. Lee, A nine node finite element for analysis of geometrically nonlinear shells. ht. J. numer. ð. Engng 26, 194551962 (1988). J. C. Simo. D. D. Fox and M. S. Rafai, Formulation and computational aspects of a stress resultant geometrically exact shell model. In Computational Mechanics ‘88: Theory and Applications (Edited by S. N. Atluri and G. Yagawa), Vol. 1, pp. 26.11.1-26.ii.9. Springer, Berlin (1988). A. F. Saleeb, T. Y. Chang and S. Yingyeunyong, A mixed formulation of Co linear triangular plate/shell element-the role of edge shear constraints. Int. J. numer. Mefh. Engng 26, i101~1128 (1988). A. F. Saleeb. T. Y. Charm, W. Graf and S. Yinaveunyong, A hybrid/mixed model for nonlinear shell&alysis and its applications to large-rotation problems. Int. J. numer. Meth. Engng 29, 4077446 (1990). J. L. Brombolich and P. L. Gould, Finite element analysis of shells of revolution by minimization of the potential energy functional. Proc: Symp. Applied Finite Element Method Civil Engineering, Vanderbilt University, pp. 2799307 (1969). I. N. Katz. A. G. Peano and M. P. Morrow. Nodal variables for complete conforming finite elements of arbitrary polynomial order. Compuf. Math. Appl. 4, 85112 (1978). B. A. Szabo, P. K. Basu and M. P. Rossow, Adaptive finite element analysis based on p-convergence. NASA Publ. 2059, pp. 43-50 (1978).

Behavior

of laminated

composite

72. I. Babuska, I. N. Katz, B. A. Szabo and Greensfelder, Hierarchic families for the p-version of the finite element methods. Proc. Third IMACS Int. Symp, pp. 278-286 (1979). 73. I. Babuska, B. A. Szabo and I. N. Katz, The p-version of the finite element methods. SIAM J. numer. Anal. 18(3), 515-545 (1981). 74. D. W. Wang, I. N. Katz and B. A. Szabo, Implementation of a C’ triangular element based on the p-version of the finite element method. NASA Publ. 2245, pp. 153-170 (1982). 75. J. Jirousek, Hybrid-Trefftz plate bending elements with p-method capabilities. ht. J. numer. Meth. Engng 24, 1367-1393 (1987). 76. B. A. Szabo and G. J. Sahrmann, Hierarchic plate and shell models based on p-extensions. Int. J. numer. Meth. Engng 26, 1855-1881 (1988). 77. K. S. Woo and P. K. Basu, Analysis of singular cylindrical shells by p-version of FEM. Int. J. Solid S;truct 25(2), 151-165 (1989). 78. K. S. Surana and R. M. Sorem, Curved shell elements for elastostatics with p-version in the thickness direction. Comput. Struct. 36, 701-719 (1990). 79. K. S. Surana and R. M. Sorem, p-version hierarchical three dimensional curved shell element for elastostatits. Int. J. numer. Meth. Engng (in press). Large displacement 80. K. J. Bathe and S. Bolourchi, analysis of three-dimensional beam structures. Int. J. numer. Meth. Engng 14, 961-986 (1979). non81. K. S. Surana and R. M. Sorem, Geometrically linear formulation for three dimensional curved beam elements with large rotations. Int. J. numer. Meth. Engng 28, 43-73 (1989). 82. E. Reissner and Y. Stavsky, Bending and stretching of certain types of heterogeneous elastic plates. J. appl. Mech. 23, 402-408 (1961). 83. J. M. Whitney and A. W. Leissa, Analysis of heterogeneous anisotropic plates. J. appl. Mech. 28, 262-266 (1989). bi-direc84. N. J. Pagano, Exact solutions for rectangular tional composites and sandwich plates. J. Compos. Muter. 4, 20-34 (1970). 85. N. J. Pagano, Influence of shear coupling in cylindrical bending of anisotropic plates. J. Compos. Mater. 4, 330-343 (1970). 86. N. J. Pagan0 and S. J. Hatfield, Elastic behavior of multidirectional laminated composites. AIAA J. 10, 931-933 (1972). and 87. S. Srinivas and A. K. Rao, Bending vibration buckling of simply supported thick orthotropic rectangular plates and laminates. Int. J. Solids Struct. 6, 1463-1481 (1970). 88. E. Reissner, On the theory of bending of elastic plates. J. math. Phys. 23, 184-191 (1944). 89. E. Reissner, The effect of transverse shear deformation on the bending of elastic plates. J. appl. Mech. 12, A69-A77 (1945). 90. E. Reissner, On bending of elastic plates. Q. appl. Math. 5, 55-68 (1947). 91. R. D. Mindlin, Influence of rotary inertia and shear on flexural motions of isotropic elastic plates. J. appl. Mech. 18(l), 31-38 (1951). 92. F. B. Hildebrand, E. Reissner and G. B. Thomas Notes on the foundations of the theory of small displacements for orthotropic shells. NACA, Technical Note 1983 (1947). 93. A. B. Basset, On the extension and flexure of cylindrical and spherical thin elastic shells. Phi/. Trans. R. Sot. Lond. Ser. A 181, 433-480 (1980). 94. J. M. Whitney, The effects of transverse shear deformation on the bending of laminated plates. J. Compos. Mater. 3, 435-547 (1969). 95. J. M. Whitney and N. J. Pagano, Shear deformation

96.

97. 98.

99.

100.

101.

102.

103.

104.

105.

106.

107.

108.

109.

110.

111.

112.

113.

114.

115.

plates and shells

65

in heterogeneous anisotropic plates. J. appl. Mech. 37, 1031-1036 (1970). P. C. Yang, C. H. Norris and Y. Stavsky, Elastic wave propagation in heterogeneous plates. Int. J. Solids Struct. 2, 665-684 (1966). S. T. Mau, A refined laminated plate theory. J. aDa/. Mech. 40, 606-607 (1973). _ .. R. B. Nelson and D. R. Larch, A refined theorv of laminated orthotropic plates. J. appl. Mech.. 41, 177-183 (1974). K. H. Lo, R. M. Christensen and E. M. Wu, A higher-order theory of plate deformation. Part 1, homogeneous plates. J. appl. Mech. 44, 663-668 (1977). K. H. Lo, R. M. Christensen and E. M. Wu, A higher-order theory of plate deformation. Part 2, laminated plates. J. appt. Mech. 44, 669-676 (1977). M. Levinson, An accurate, sample theory of the statics and dynamics of elastic plates. Mech. Res. Commun. 7, 343-350 (1980). M. V. V. Murthy, An improved transverse shear deformation theor; for laminated anisotropic plates. NASA Technical Paper 1903 (1981). J. N. Reddy, A refined nonlinear theory of plates with transverse shear deformation. J. Solids Struct. 20, 881-896 (1984). B. A. Szabo and G. J. Sahrmann, Hierarchic plate and shell models based on p-extensions. Int. J. numer. Meth. Engng 26, 1855-1881 (1988). L. J. Brombolich and P. L. Gould, Finite element analysis of shells of revolution by minimization of the potential energy functional. Pro;. Symp. Appt. Finite Element Method Civil Enaineerinp. Vanderbilt University, pp. 279-307 (1960).P. K. Basu, B. A. Szabo and M. P. Rossow, Theoretical manual and user guide to COMET-X, an advanced computer program for stress analysis. Federal Railroad Administration. FRAiORD-77160 (1977). P. K. Basu and R. M. Lgmprecht, Some trends in computerized stress analysis. Proc. 7th Conf: Electronic. Computations, ASCE, St Louis, pp. 312-330 (1979). I. Babuska, 0. C. Zienkiewicz, I. Gago and E. R. de A. Oliveira, Accuracy Estimates and Adaptive Refinements in Finite Element Computations. Wiley, New York (1986). K. S. Surana and P. Kalim, Finite element formulation for axisymmetric shell heat conduction with temperature and gradients. Presented at the Winter Annual Meeting of the ASME, Miami Beach, FL. (1985). K. S. Surana and P. Kalim, Transition finite elements with temperature and temperature gradients as primary variables for axisymmetric heat conduction. Comput. Struct. 24, 197-211 (1986). K. S. Surana and R. K. Phillips, Three dimensional curved shell finite elements for heat conduction. Comput. Struct. 25, 775-785 (1987). K. S. Surana and R. K. Phillips, Three dimensional solid-shell transition finite element for heat conduction. Comput. Struct. 26, 941-950 (1987). K. S. Surana and R. M. Sorem, Higher-order hierarchical plates and curved shell elements based on p-approximation in the thickness direction. 7rh Int. Co& Mathematical and Computational Modelling, Chicago, IL (1989). K. S. Surana and R. M. Sorem, Curved shell element for elastostatics with p-version in the thickness direction. Cornput. Struct. 36, 701-719 (1990). K. S. Suiana and R. M. Sorem,‘p-Approximation curved shell elements for elastostatic analysis of laminated composite plates and shells. ASME Pressure Vessel and Piping Conference, Honolulu, HI (1989).

66

J. H. Liu and K. S. Surana

116. K. S. Surana and R. M. Sorem, Curved shell elements based on hierarchical p-approximation in the thickness direction for linear elastic analysis of laminated composites. Int. .I. numer. Meth. Engng 29, 139331420 (1990). 117. K. S. Woo and P. K. Basu, Analysis of cylindrical shells by p-version of finite element method. Int. J. Solids Struct. 25, 151-165 (1989). 118. K. S. Surana and R. M. Sorem, Completely hierarchical curved shell element for laminated composite plates and shells. Comput. Mech. 7, 237-251 (1991). 119. R. M. Sorem, Linear geometrically nonlinear p-

version three dimensional curved shell element for isotropic, general orthotropic and laminated composites materials. Ph.D. Dissertation, Department of Mechanical Engineering, University of Kansas, Lawrence (1991). 120. T. Y. Chang and K. Sawamiphakdi, Large deformation of laminated shells by finite element method. Proc. Comput. Meth. Nonlinear Struct. Solid Mech., Washington, DC (1980). 121. L. A. Schmit and G. R. Monforton, Finite element analysis of sandwich plate and cylindrical shells with laminate faces. AIAA J. 8, 1454-1461 (1970).