An equilibrium method for prediction of transverse shear stresses in a thick laminated plate

Compu:rrr & Sirucrvr~s Prnrcd I” Great Bnran.

Vul

13. No

1. pp. 139-146.

1936 t

Oo-li--94956 53.00 - M) 1986 Pergamon Press Lrd

AN EQUILIBRIUM METHOD FOR PREDICTION OF TRANSVERSE SHEAR STRESSES IN A THICK LAMINATED PLATE Department

of Civil Engineering.

RE:\~ A. Universit)

CH-\C’I)HL1KI

of Utah. Salt Lake Cify.

Utah WI 12. U.S.;\.

Abstract-First two equations of equilibrium are utilized to compute the transverse shear wes\ variation through thickness of a thick laminated plate after in-plane stresses have been computed using an assumed quadratic displacement triangular element based on transverse inextensibility and la)erwise coWant shear angle theory I LCST). Centroid of the triangle is the point ofexceptional accuracy for transverse shear stresses. Numerical results indicate close agreement with elabticit) theory. An interesting comparison between the present theory and that based on ashumed glress hybrid finite element approach suggests that the latter does not satisfy the condition of free normal traction at the edge. Comparison with numerical results obtained by using constant shear angle theory suggests that LCST is close to the elasticity solution while the CST is clober to classical tCLT) solution. Ir i5 ai40 demonstrated that the reduced integration gives faster convergence when the present theory i> applied to a thin plate.

likely to fail due to delamination caused by the scissoring effect of these stresses. Analytical solutions are few and are usually restricted to problems with simple loading and boundary conditions[ I, 21. Finite element method (FEM) is. therefore. a practical alternative for solution to more complex problems. Most of the finite element analyses are either based on classical lamination theory (CLT) which ignores shear deformation effect altogether. or constant shear angle theory (CST) which assumes constant shear angle through the entire thickness[3-81. The finite element analyses that are based on layerwise constant shear angle theory (LCST) and transverse inextensibility are due to Mau er trI.[9]. Mawenya et N/.[ IO], and Seide and Chang[ I I] (Fig. I). However. only Pryor et rr/.[3]. Spilker ct (11.14. 51, and Mau ef 01.191have investigated transverse shear stress variation through thickness. Assumed linear displacement triangular element due to Seide and Chang[ I I] can only compute average transverse shear stress through thickness of a layer. The quadrilateral element due to Mau (gt n/.[9] based on assumed stress hybrid approach has the capability of obtaining the variation of transverse shear stresses through thickness of a thick multilayer plate by considering stresses as unknown nodal parameters and imposing constraints on the compatibility of these stresses at each interface. However, this method appears to be limited in its applicability. as the formulation of the stiffness matrix involves too many matrix inversions at the element level. More importantly. extension by Spilker et a/.[41 to triangular element shape has produced a stiff element, and traction free edge condition is not satisfied. While quadrilateral element shape is adequate for plates with regular boundaries. use of triangular

NOhlEYCLATURE

Length of a square plate or width of an infinite strip Elastic (material) stiffness matrix of the material Nodal displacement vector for the ,jth triangular layer element belonging to the ith layer Young’s modulus for isotropic material Elastic moduli of anisotropic material (stretching) Shear moduli of anisotropic material Strain-displacement relation matrix Mesh size for a regular or uniform mesh rth component of unit normal vector rth component of unit vector normal to the kth side of the triangular element Total number of layers Number of subdivisions Applied surface load Area of thejth triangular layer element Total thickness of a laminated plate. Thickness of layer i Displacement components in .s. x, ; direclions. respectively Cartesian coordinate-: is measured from the bottom surface of a layer: Z is measured from the bottom (reference) surface of the laminate Boundary of the triangular element Length of&h side of the jth element Poisson’s ratio of an isotropic material Major and minor Poisson’s ratios in the plane of the ftbers 0:” . u:’ ’ . 2:‘: In-plane normid and shear slresses at interface i 7:‘;. r:!.! Transverse shear stresses in layer i

INTRODUCTION

of transverse shear stress variation through thickness of a laminated plate has assumed increasing importance because such a plate is more Determination

I ?Y

I

z.2 Deformed

Normal

Fig. 2. LCST

Fig.

I.

Possible

severe cross-wctional laminated plate.

warping

in a thick

element is more convenient to represent such irregular boundaries as holes and cutouts. Seide and Chaudhuri[ 12. 131have developed a highly efficient and accurate triangular plate/shell element based on LCST and assumed quadratic displacement field. The present paper will discuss the use of this element to obtain accurate transverse shear stress variation through thickness.

R.ACKCROtiND

based cumpo\ite

plate clement.

stress variation through thickness of a thick laminated plate using the stresses computed by assumed quadratic displacement potential energy approach and first two equations of equilibrium. The concept of use of equilibrium equations to compute transverse shear stresses is due to Pryor et 01.131 who utilized it for a quadrilateral element based on CST.

METHOD

OF ASALI’SIS

The first t\\‘o equations of equilibrium inside layer i can be written as

at a point

IXFORXIATION

Each triangular element belonging to the layer i is bounded by the ith and i + Ith interfaces (Fig. 2). Each interface triangle is characterized by six nodes, each node being associated with three nodal displacement parameters. By virtue of the assumption of transverse inextensibility. ~1’does not vary with thickness. The assumption of layerwise constant shear angle guarantees piecewise linear variation of inplane displacements through the plate thickness. The number of degrees of freedom per node is then 7N + 3 for an N-layer composite plate element. Details concerning formulation of element stiffness matrix and consistent load vector. effect of numerical integration on convergence of displacements and stresses. solution of global equations and element stresses, etc. have been presented by Chaudhuri[ 121. Once displacements are determined. the element stresses can be obtained using the relation[ 12. 131 {u’} = [C’][A’][G:][d;J.

(1)

for r. 1 = Itor x). 2tor ~5). Repeated indices here indicate summation. Performing integration with respect to : will yield

The assumption leads to

u);‘(c) =

of layerwise

(,- fI

&’

constant

+

;

a):-

(3)

shear angle

“,

(3)

For assumed quadratic displacement field. u)::,.(,-I is constant with respect to .V and _v. Use of this knowledge and eqn (4) coupled with applying divergence theorem on the right-hand side of eqn (3) will give

T)!‘( . -) c

where [C’], [A’J. and [G;] are presented in the Appendix [eqns (A3): tA7) and tA8): and (A9) and (AlO), respectivelyj. It is noted that while accurate in-plane stresses can be computed at midpoints of the sides of the triangle, only average (through layer-thickness) transverse shear stresses. T:ii and 7:‘: are determined by using eqn (I). The following section will present a method of determination of transverse shear

u):!,-(3d,.

$(,7) - T;;)(O) = -

-

$‘(O)

A square. simply supported isotropic plate subjected to uniform normal pressure. [J,) is studied first. The edges of the plate are assumed to be attached to rigid end plates which are free to rotate. Each side of the plate is 100 cm while the plate thickness is 1 cm. E and u of the plate material are 2.1 x IO6 kg/cm’ and 0.3 respectively. The exact solution is availabie in Ref. [13]. Symmetry conditions about the plate centerlines and diagonals allow the modeling to be limited to one-eighth of the plate surface. Antisymmetry about the middle surface permits us to consider a plate of half the thickness subjected to uniform pressure pJ:!. so that at the middle surface, the inplane displacements vanish (Fig. 4) Y,._is given by

Fig. 1. .$h triangular plate element interface.

T,;

In the above equation. quantities are referred to XThside of thejth triangular layer element belonging to layer i. n: and fyi are given by (refer to Fig. 3) eqn (AIZ) and (Al3). As element stresses al:‘. al,I’+ ” are linear in .Vand .v, one Gauss point is sufficient for numerical integration on each side of the triangle! 131. It is noteworthy that stresses o$’ etc. at midpoints are exceptionally accurate. Equation (5) on numerical integration finally yields

+ c$l node 4 -

CT:‘/ node

(Yz - Yl) -

??;I

Y 2 + 7:; i

X21

6

node 4

node 6

=

(7)

-

The expression for r,, consists of second derivatives of displacements. which are constants with respect to .Y and y. The transverse shear stresses are then automatically computed at the centroid of the triangular element. Centroid is then the point of exceptional accuracy for T,_ and ?,.._. Convergence of (T.,:)~;~~. i.e. T,, at .Y - d?. y - 0 is presented in Fig. 5 for both full integration (quintic order, error 0(1r6)f14]) and reduced integration (quadratic order, error 0(h3)[ 141). As expected, for a thin plate like the one under consideration (a/r = SO), the reduced integration yields faster convergence.

u=v=o

(x2 - Xl)

(Symm. Cond.)

t/2 .-L ‘-?-

(6)

x, node

Y

x (II:!-XI)

- UY)lnodeYZ +

T~;‘)lnodeX2)1

6

Similar expression

can be obtained

6

for T$.

NUMERICAL EXAMPLES

In order to test the accuracy of the present theory as welt as the performance of the associated triangular element. the following examples will be considered.

Fig. 4. Finite element mesh and displ~cemenl boundary conditions of a square isotropic plate.

remaining equal

transverse

to 0.1

~1,;. defined parallel

in the direction

of’the

tibcr\. ratio.

The problem

for the present

theory

gano[ II.

An

exact

of the present Solutions.

formulation

have

been obtained

rrl.151 using (LCST)

volves layer

problem severe

under

cross-sectional

balanced

finitely loaded The

long

s),mmetric fiber

is 2 in. E,, IOh and

I x

shear

modulus

G,:

are assumed

of width

IO’ psi. and

equal

the thickness

transverse

dx,y)

x

shear

Q

are

In-plane modulus

IO” psi while

=

6(a)].

of each

material

respectively.

to 0.5

(I and

[Fig.

and E1: of the layer

25 x GII

is a three-

pressure

while

in-

ply (YO”iO”/YO”) in-

strip

normal

LI is 24 in..

which

warping.

cross

reinforced

by 21 sinusoidal

width

layer

consideration.

sin(ny/a)

the

the

and

strip

length

iriangular

imposing points ments

in

about

the

with

presents

the

and testifies method.

present

tical to its analytical by

Seide[l61.

present

theory

compared both

It

yields

of the present is iden-

(not shown). to

note

given

that

the

result

to that due to hlau

CI rrl.[XI.

even though

fact that

normal

satisfied

when

edge condition

rtrcss

hybrid

is not

approach

of

Classical Lamination Theory Hvbrid Stress Finite Element (CST) H&rid Stress Finite Element (LCST) Present Finite Element (LCST) Elasticity Theory

-8r----

q(x.y)/2 2

for yi;

This is per-hap< due to the

traction-free assumed

0

0.8

8

thickness

accurate

A

tz

solu-

Figure

:I more

-.-

t/6 gz.!--r-------=i 113 L.L

success

interesting

are based on LCST.

1.4

is used.

FEM based solution counterpart

is

7). Ten

of’ freedom

to the elasticity scheme

to the complete

The

the con-

(Fig.

of T,,ip,,through

variation

at the

First.

degrees

convergence

f’ull integration

and per-

;\nti\ymmetrq

is tested

60 unconstrained

for

when

right-nn6(b))

\\ hile displace-

vanish.

divisions

the

displacement

is assumed.

surface

of T,, at interface

are needed

trio

and in the J’ direction

.V direction

the middle

of’

along

[Fig.

the same 1’ coordinate

vergence

tion

independence

of position

of’ equal

to the plate

having

in

\I hich fulfills

as shown

the condition

pendicular

of’

of’ :I

~15shotvn

by combining

elements

utilizing

the analysis

element

condition

is obtained

(CST).

is perf’ormed

to permit

(‘I

theor)

shear angle theor!,

analysis

displacements

gled

ele-

shape.

angle

(r/2. bvhich is subdivided

plane-strain

tinite

element shear

A one-dimensional

stresses

h) Scicle[ 161.

stress hybrid

constant

conditions

Fig. 6(b).

Pa-

by Xlau c’f rr/.[81 and Spiller

element

of width

Lamib),

thr: :t\\umption\i

and quadrilateral

and also constant

symmetry

element.

obtained

using

to

ground

and cI;i\sical

has been giLen

layer-wise

The finite

to 0.25.

ir\sunwJ

i\ an ideal te>tlng

assumed

ment

i\ tahen equal

solutions

theory

;md

b) the \trez\

~a::. i\ ;IIW

solution

utilizing

~8,: =

perpcndicul;~r cLtu4

elasticity

tCLT)

Theory

to be

ratio

and the a~\v~iatecl

it has exact

nation

The

directic)n.

Poisson’s

bt: the same.

because

c‘:: i\ ;l\\urneJ Poi\son‘\

to be ratio af>tr;~in\.

to the tiber

Transverse

strip

modulu\

X IOh phi. .1I;ijor

4

6

8

10

12 14

16

18 20

-.-

Classical

__+_-.

Hybrid

Stress

Finite

Element

(CST)

__G__.

Hybrid

Stress

Finite

Element

(LCST)

Present

Finite

Element

___x___

Lamination

Elasticity

Theory

(LCST)

Theory

u=v=o (Symm.

112’ o

Mau

0.4

CI tr/.[8]

Similar

0.8

is used

deviation

of Spilker

is closer

ticity

solution.

ratio

(a/r)

shear

stress

A final

the

which

solution

can be obtained

is that of

angle-ply

laminate

uniformly

same as those

of example

Symmetry the

subjected

of

of freedom.

2

layer

that the in-plane

and reduced in Fig.

results

lO.Oi

!

‘=. I-”

7.5 /

0

art: the

half

reduce

Variation

k--~ -+--

Full

.._-___ Integration

(Quintic)

!

using

to converge.

Classical

vanish

I

-e---

on

Reduced

Integration!

computed

the Both

this tipre full

all

using

schemes

=

is 112

that

-.-_--..___

._.

integration

the results

ob-

N

Lamination Theory

Present Finite Element

20

30

Elasticitv

Theory

40

60

50

/

(Quadratic)

the number

of T,,(Z) at .Y = rri

thickness. from

.

Cl_T

the middle

displacements

I I. It is evident

_____- -

us to

I -.-.-

I

load p,,12. This and the

integration

obtained

10

IO. The

numerical

are yet

6

in Fig.

____-

I

/j,,.

(W-ply) permits

and

considerbly

17tri36 through

scheme

sim-

pressure

lamination

to uniform

of degrees

the

has edges.

that E, , is 10 x

7 except

bottom

surface

shown

is shown

of an orthotropic

the middle

full

and fat

square

plate

distributed

configuration

properties

36. .Y =

the effect plate

in the same sense as that of example

Geometrical

condition

of aspect

transverse

in a thin

The square

is under

layer.

effect

demonstrates

symmetric

ply supported

only

than to rhe elas-

the

deformation

( -45”/45”/-45”).

IOh psi.

solution

9 shows

nondimensionalized

problem

an unbalanced

model

solution that CST

at the interface.

shear

elastic

of T,;I:).

It is also noteworthy

to CLT

analytical

I and

2.0

in the computation

Figure

on

of layer

1.6

arises also in case 0fCST

er a1.[5].

solution

which

1.2

Cod.)

70

60 a/t

90

100

2T,*

fqo

50

rained using the present theor! are compared with CLT solution. obtained using Green’s [ 17-701 dwDie Fourier series 122 x 2 terms) approach. It is interesting to observe that the ;,._ at the middle surface obtained using reduced integration scheme converges to the corresponding CLT result. Ho\\eier. the present result utilizing LCST differs considerably from the CLT solution in the sense that the latter monotonically increases with ; and reaches its maximum at the middle surfxe. LCST solution. on the contrary. reaches its maximum at the interface of layers and demonstrates a marked change of slope of r,; vs ; curve there. Because of simplifying assumptions, CLT fails to detect this behavior and will be in error. in predicting (:,:I interface, by more than 30%. This example demonstrates the importance of LCST in predicting interlaminar shear stresses even in a thin (l/r = 33.33) angle-ply laminate like the one under investigation. It is also worthwhile to note that unlike the threetayer symmetric cross-ply strip, (T,:),,,;,, occurs at the interface. where interlaminar shear strength assumes the lowest value. This creates ;I more dungerous potential for delamination failure in the angle-ply laminate.

CONCI.C’SIOSS An accurate method for prediction of transverse [interiamin~lr} shear stress variation through thickness ofa laminated plate is presented. Even though the method is demonstrated here for an assumed quadratic displacement triangular element. the principle behind the method is general enough to be utilized for any element shape. This paper also demonstrates that it is possible to obtain an xcurate solution. based on LCST. with a conventional assumed potential energy approach. It is also shown that centroid of an interface triangle is the point of exceptional accuracy for transverse shear stresses. Deviation ONTO.computed by MXI cl trl.fYj from the corresponding elasticity solution (example 2) is perhaps due to the fact that the free normal traction condition along the edge is only approximately satisfied by the assumed stress hybrid method. The close correlation of the present solution suggests that the traction-free edge condition is automatically satisfied in the present case. The present study also demonstrates that transverse shear deformation effect is significant in an angleply laminate. even if such a laminate can be otherIvise categorized as a thin plate on the basis ofwidth to thickness ratio. CLT solution may be in error by more than 30% in prediction ot’ transverse shear stresses at the interface of layers for such a laminate. It should be noted that the success of the method. presented herein, depends on either pointwise convergence of derivatives of displacements or in-plane force equilibriLim at a point is. x). HOWever. derivatives of displacements in FEltI con-

vtfrge onl) in the mean >quare $ense and in-plane force equilibrium may. in general. be satislied for homogeneous and ~!mmetrically laminated plates alone. An alternative method has been developed by Chaudhuri[ 12) for accurate computation of transverse shear s[resc variation through thickness of an arbitrarily laminated plats. r\c.l(/~rr~~,j~,t/i’/~~~,~rr~-~to~tof the revolts reported here were obtained during investigations. supported by NASA Langley. under supervision of Profezwr Paul Seitle of USC. The author also acknowledges the assistance of Sfr. Kamal .-\bu-.-\rja of the L:niverGty of Utah in carrying the CLT computations of Example 3.

KEFERESCES 1. N. J. Pagano. Exact wfution fur composite laminates in cylindrical bending. J. Compos. Mater. 3, 398-41 I (1969). 2. N. J. Pagano. Exact solution for rectangular bidirectional composites nnd sandnich platef. J. Cortrp~~v. IW~//PI..-1, X-35 f 1970). 3. C. W. Pryor Jr. and R. iLl. Barker. .-\ finite element analysis inclttding transverse shear effects for applications to larnjn~~ted plate5. AllrlJ 3.912-913 i t971). 4. R. L, Spilker. 0. Orringer. E. A. Witmcr. S. Verbiese. S. E. French and ,A. Harris. Use ofh\,brid stress finite element model for the static and dynamic analysis of multilayer composite plates and shells. Dept. AMMRC CTR 76-29. ASRL TR 1X1-2. ASHL. MIT. Cambridge. MA. Sept. (lY76l. 3 . R. I_. Sprlker. S. C. Chou ;md 0. Orrin_eer. Altern~lte hybrid-itress elements for analysis of m;itilayer comoositr nlates. J. Cf>rrrno.\. .Iltrlcr. I I, 51-70 11977). 6. ‘s. C. 61nda and R. X\;at;lrajan. Finite elemenl analysis of laminaled composite plate>. /jr/. J. ,Vrwrcr Mcth. Erf‘&vl.c i-6, 69-79 ( IY7Y 1. 7. J. N. Reddy. A penalty plate bending eiemcnt for the ;inalysis of laminalrd anisotropic composite plates,

Itrf. J. iVit777~.itfcrir. I-ltlgrry 15, i IX?- 1206(19X01. 8, G. R. Bhashyam and R. H. Gallagher. A triangular

9.

IO.

I I.

12.

i3.

shear-tlexible finite element for moder;~tely thick laminaled plates. Co/trprrt. .I/ct/r. Appl. .\lrt 11. Li:,,p/r,g JO, 309-326 f 19X3). S. T. Mau. P. Tone and T. H. H. Pian. Finite element solutions for Iam-inated thick plates. J. Co!trpw. ivtrrc*r. 6, 304-3 I I ( 1972). A. S. hlawenva and J. D. lXtvi\. Finite clement bending analysis o+mullil;~yer plates. Irli. i. .V/rrrrc,r. ,Llrrl~. &,Is//L’ 8, 2 15-226 I 19741. P. Se’ide and P. H. H. Chang. Finite element analysis of laminated platss and jhells. NASA-CR-157106. 157107 (197X1. R. A. Chaudhuri. Static anafysix of fiber reinforced laminated plates and \hellz \t ith shear deform~~tion using quadratic trisneul;rr elements. Ph.D. dis!erta{ion. Dept. Civil Engng. USC. Lo< Angeles. CA

Il9X31. P. Sride and R. A. Chaudhuri. Triangulw’ finite element for analysis of thick laminated shells. Int. J. fVtVccmer. Meth. Engng (to appear).

~~~I~.s~~~~~,~;~~/~ t Tk ~~,~~~/if~~ 1-I. F. J Plantema. S~~/?~/~i.i(./f rrtrtl ~ti~.~ii~l~ i?~‘S~j~r(~~~,~i,~t 3etrirr.s. 1%7lc.\trrrcfSlWli\l. John Wiley. New York ( IYhhl. T/r,, Fi/ri/c E/r,/rrt*/rt ,\lct/wtl. 15. 0. C. Zienkiewicz. McGraw-Hill. London (1977). 16. P. Scide. An improved approximate theory t’or the bending oflaminatedplates. .Ilrc.lrtr!ric-.y 7;dcr\ I Edited by S. Ncmat-Na\sert. Vol. .i. pp. J?l-165. Perg;tmon. New York i l%(t). 17 ,A. E. Green. Double Fourier sericq and boundar>

IX.

.\/tr:‘tr:ir,c,.Ser.

7. Vol.

36. pp. h!Y-6x8

Strucrufibibehavior

19. R. .4. Chaudhuri.

I - VI:“:,

I IY-lil.

uf FRP reclangular plates and cylindrical \hells. X1.S. thesis. Dept. .-\ero. Engng. I.I.T.. Madra\. India llY7J1. 20. K. .A. Chaudhuri. V. N. Kunuhkazseril and K. Balwarnan. Free vibration\ of rectangular multila~er anibotropic plate?. Presented a~ First Conf. Reint. PIa+ tic3 ,~nd Their Aerospace Applic.. ;II Vikram Sarabhai Space Center. ISRO. Trivandrum. India. .4ug. f 1972).

CLZ =

c,, [A’]

V12EZI

kIE,l

I - v,:vz, =

= G:,:

referred

I - “I:“:, :

T,,

to in eqn

I -

= G,,:

Y,:“:, c,,

( I I ih given

(A6)

:

= G:J. b)

t.471 ,\PPESDIS This section presents some of the matrice and equalions referred 10 in the background information and theoretical formulation of this paper. Stress-strain relations for a lamina are given by

where

(If/ = Ic’l(c’t. where 131

in the context

of FEM

formulation

I AXa)

0 of Refs.

1I2

and

(Mb)

f‘4XCl

L

C.VJ

[G:]

as referred

I [Cl =ITl’lc]lTl

For a lamina with fiber orientation H. &I can be obtained from the orthotropic ma(erial matrix ICI by means of [he relalionl I I

whose

individual 6 and i =

to in eqn

( I)

submalrices

is

are defined

for

I

:

1.2.

I. 2. . . IY b!

(A-l)

(*4 I oa 1

wbile cos’ sin’ (I’]

=

-sin’ 0 0

H

sin’

H

H

cosz H

H

sin’

l41Obl

H

0 0 IAS) IAIOC)

Six nonvanishiny elements of ICI can terms ofthe generalized Youn_e’s moduli.

be eupre3red in Poisson’s ratios.

0 L.

(bi 1,

& 2

J

.

i:--\,

l.-ll2CI

“i = I. ,,; =

‘: - ‘1 1-5

l.-\IYl