Finite element estimation of elastic interlaminar stresses in laminates

Finite element estimation of elastic interlaminar stresses in laminates

‘%“[email protected]& vol. 3. No 5, pp. 451466. Printed in Great Britain. 0961-9526/93 s6.00+ .oo 0 1993 Pergamon Pra8 Ltd 1993. FINITE ELEMENT ESTIMAT...

896KB Sizes 3 Downloads 94 Views

‘%“[email protected]& vol. 3. No 5, pp. 451466. Printed in Great Britain.

0961-9526/93 s6.00+ .oo 0 1993 Pergamon Pra8 Ltd

1993.

FINITE ELEMENT ESTIMATION OF ELASTIC INTERLAMINAR STRESSESIN LAMINATES R. GO~INDARAJAN, A. V. KRISHNA MUIZTY, K. P. v. RAGHuaAM

VIJAYAKUMAR

and

Department of Aerospace Engineering, Indian Institute of Science, Bangalore 560012, India (Received 19 November 1991; revised version accepted 4 September 1992) Abstract-Low interlaminar strength and the consequent possibility of interlaminar failures in composite laminates demand an examination of interlaminar stressesand/or strains to ensure their satisfactory performance. As a first approximation, these stressescan be obtained from thicknesswise integration of ply equilibrium equations using in-plane stresses from the classical laminated plate theory. Implementation of this approach in the finite element form requires evaluation of third and fourth order derivatives of the displacement functions in an element. Hence, a high precision element developed by Jayachandrabose and Kirkhope (1985) is used here and the required derivatives are obtained in two ways. (i) from direct differentiation of element shape functions; and (ii) by adapting a finite difference technique applied to the nodal strains and curvatures obtained from the finite element analysis. Numerical results obtained for a three-layered symmetric and a two-layered asymmetric laminate show that the second scheme is quite effective compared to the first scheme particularly for the case of asymmetric laminates.

INTRODUCTION

Interlaminar strengthof fiber-reinforcedplasticcompositelaminatesis known to bevery much lower thanthe ply strength.Consequently,onehasto designlaminatessuchthat the possibility of interlaminar failure is avoided.With regardto relatively weakinterlaminar zonesin a laminate,they areproneto developdamagessuchassmall delaminationseither dueto largeinterlaminar stressesinducednearfree edges,boundariesof rivet holes,ply drops etc. or during serviceand/or manufacture. Mechanismsof damageonset are complex,involving simultaneousoccurrenceof matrix cracks,fiber breaks,delaminations etc. and so prediction of damageonset is, indeed, highly complicated as it requires considerationof many details. However,it is desirableto retain a proceduresufficiently simple to direct the initial designprocess. A simple approachmay be initially usedto identify the interfaceof failure and the minimum load at which the failure occurs at discretepoints, followed by a detailed analysisof damagesiteto gainan understandingof the effect of damageon the behaviour of the laminate. The first stepusedis calledthe global approachand the second,the local approach.Suchglobal-local conceptsaregainingincreasedacceptanceasusefulanalytical tools (Stroudet al., 1985;Vidussoniet al., 1990;Mao and Sun, 1991;RansomandKnight Jr, 1990;Shah and Krishna Murty, 1991). Apart from the availability of a reliable interlaminar failure criteria, a method of estimation of interlaminar stressesand strains is a necessaryprerequisitein order to identify critical sitesby global approach.Interlaminar failure criteria could be basedon stressesor strains or both. Therefore, it is necessaryto ensurethat both stressesand strainsareavailablefrom theanalysis.It is known that the interlaminartransversestrains, unlike thecorrespondingstresses,arediscontinuousacrossthe interface.This featuremay createspecialdifficulties evenin three-dimensionalfinite elementschemesto obtain interlaminar strains. The use of three-dimensionalfinite element approach, is, perhaps, prohibitive for this purposeas it involves a large computational effort. A searchfor simpler approachesthan the three-dimensionalfinite elementanalysisis worthwhile. In the early 196Os,the classical plate theories (Reissnerand Stavsky, 1961; Dong et al., 1962;Whitney and Leissa, 1969;Whitney, 1969)have been usedfor the analysisof laminated composites.[A few elasticity solutions (Srinivas and Rao, 1970; 451

452

R.

GOVINDARNAN

et al.

Pagano,1970;Pagan0and Hatfield, 1972)werealsoobtainedfor particularplatebending problems.] Hayashi (1967)presentedthe first analyticalmodel for treating interlaminar shearstresses.In general,higher order theories(Whitney and Sun, 1973;Nelson and Larch, 1974;Lo et al., 1977;Murthy, 1981;Reddy,1984;Phan andReddy, 1985;Krishna Murty and Vellaichamy, 1988;Harikumar and Krishna Murty, 1991)are expectedto do betterthan the classicallaminatedplatetheoryin estimatinginterlaminar stresses.Higher order theoriesaregenerallymuch more complicatedand someof them (Reddy,1989)are attractivealternatesfor local analysis.However, for the purposeof global analysis,it is preferableto keepthe analysisas simple as possible,sincea more detailedlocal analysis of the damagedsite is anyway necessaryat a later stage.With this view, the classical laminatedplatetheory is consideredto beadequatefor the global analysis.This approach will also allow us to take advantageof finite elementscurrently availablein the literature. In recentyears,severaltwo-dimensionalfinite elementmodelshavebeendeveloped to estimate interlaminar stresses.Broadly, such elementscan be classified into two categoriesbasedon the variableschosenfor the formulation. The first, perhapsthe most popular one, has displacementsas the primary variables(Pryor Jr and Barker, 1971; Engblom and Ochoa, 1985;Owen and Li, 1987;Kant and Pandya, 1988;Reddyet al., 1989) and the second type may be generically called hybrid models as they allow independentassumptionsof stressand displacementswithin eachlayer, Spilker (1982). In displacement-basedfinite elements,interlaminar stressesare estimated either directly by using the constitutiverelationsbetweenstressesand strainsor by integration of the ply equilibrium equations.The latter is generallymore accurateand preservesthe continuity of interlaminar stressesacrossthe layer interfaces. Jayachandrabose and Kirkhope (1985a,b) presentedan explicit formulation of the 36 d.f. triangular plate elementbased on classicallaminated plate theory, originally developedby Cowperet al. In the presentstudy, this elementhasbeenadaptedto study its capabilityto estimateinterlaminar stressesin laminatesundertransverseloading.With this element, it is possible to obtain interlaminar stressesand strains using the two approachesnamely: (i) from direct differentiation of elementshapefunctions; and (ii) by adapting a finite difference technique applied to the nodal strains and curvaturesobtained from the finite element analysis. Numerical results are presentedfor both symmetric and asymmetriclaminates. FORMULATION

Figure 1 showsthe geometryof the elementand the coordinatesystemused. The elementhasthreenodesat its verticesand 12 degreesof freedomper node. The d.f. are the inplane displacementsu and v and their first derivativesu,, , Us, v,~, Us and the normal displacementw and its first and secondderivatives,namely w,* wy , w,= w,~, We, The geometryof the laminate is shownin Fig. 2. The displacementfield of the elementis written in terms of its area coordinates< and 1 (definedin Fig. 1) as:

u = ~~11vw.l (1) A1 and AZ the polynomial terms describingthe displacementfield, givenby:

R and S are the transformation matrices(seeAppendix I). 6,, 8, and 8, are the nodal

Finite elementestimation of elasticinterlaminar stresses

eJl4

453

are Area Coordinates

5=

Area Area

203 123

?=

Area 103 Ana !23

c=

Area 102 Area 123

Fig. 1. Elementgeometry.

k+l k

Fig. 2. Geometryof the laminate.

displacement vectors containing the degrees of freedom, given by: (6JT

=

[Ul

U,l

uyl

u2

(s,)=

=

[Ill

II,,

uyl

v2

fW

=

tw1

Wxl

WY1

***

243

*-a

u3

w**,

WAyl

*** *a-

UC]

(3)

v,]

&yl

w2

**-

w3

“’

yvy31.

u, and u, are the centroidal degrees of freedom and

au

4i

The strain-displacement

=

0 G

i;

wyvi

=

a2w &If i; ( >

i = 1,2,3;

etc.

relations may be written as: (4 = bO1 - ZIX)

(4) (5)

454

R. Gommm

et al.

and (6)

(7) The strain energyof the laminate may be written as:

N and M are membraneand bendingstressresultants,respectively,given by: (9) where [A], [II], [B], the effectivemembrane,bendingand bending-membranecoupling stiffnessmatricesof the laminate are given by: zk+l

k=l s i&

@,I, U,z,2) dz,

(10)

n is the total numberof layers(seeFig. 2). Qu is the transformedreducedstiffnessmatrix. Using relations[eqns(4-7)] and the displacementfield [eqn(l)] in the strain energy expression[eqn (8)] and simplifying, onegets:

(11) [K] is the elementalstiffnessmatrix. Explicit expressionsfor the elementsof the stiffness matrix [K] are availablein Jayachandrabose and Kirkhope (1985a,b). Consistentload vector P,, Py and Pz representthe appliedload over the elementin X, Y and 2 directions, respectively.To determinea consistentload vector, the appliedload is assumedto vary linearly over the element.The appliedload distribution may be written as:

PI = Px py e1= = kul~,l

(12)

1

{ 0 0 tj 0 0 1-c-q 0 0 [Aj] = 0 (I 0 0 tf 0 0 1-r-q 0 (13) 0 0 1-t-a [ 00~00~ (14) IspI= = w..l py1 PZl px2 --* PLPI. Pxl representsthe magnitudeof the load appliedat node1, in X direction. The subscripts of the otherterms shouldalsobe interpretedthe sameway. The virtual work W, doneby the applied loadsis givenby: w=

ss

(P)=(s)ds =

ss

vN=m ds

(15)

Finite element estimation of elastic interlaminar

stresses

where[a]’ = [u v w}. Using eqn (12) for (P) and eqn (1) for (S], Wean be written as

From eqn (16), the load vector ? can be identified as:

ROOT ? = 2A 0 R 0 fiti,] 0 0 s

[ 1

(17)

where l&r,]is the load intensity vector, given by: and

(19) and A is the areaof the triangle. The elementsof the matrix P aregivenin Appendix II. The consistentload vectorderivedabovecorrespondsto degreesof freedomarranged as ~lr~xl,I(yl,~2,...,~,,~l,~..,tr,,~l,~~~,~yu3.

The terms correspondingto tl, and v, in the stiffnessmatrix and load vector are later condensedout, to obtain a 36 x 36 elementstiffnessmatrix and a 36 x 1 load vector.\ STRESS ESTIMATION

Utilizing the element stiffness and load vectors derived earlier, the method of evaluatingthe nodal displacementsis standard.Having obtainedthe elementaldegreesof freedom, the in-plane strains are computed using the strain-displacementrelations [eqns(4-7)]. The in-planestressesa,, a,,, oWarethencalculatedfrom thesestrains,using constitutive relations.

wherethe superscriptk indicatesa typical layerin a n layeredlaminate(seeFig. 2). oti is the transformedreducedstiffnessmatrix. The local equilibrium equationsof a ply are given by:

=ux,x

+

QY,Y

a,,,

+

ows

+ +

~YZ,Z

=

0

qz

=

0

(21)

and it follows from the aboveeqns:

(22)

R. GOVLUDARAIM et al.

456

The boundaryconditionson a,,, oyr, a, are: uxz= t7xz =

=YZ

=oz=o

at 2 = z1

and 0, = Q at 2 = z~+~. =0 In evaluatinga,, , cryyr , a,, the condition of a,, = a,, = a, = 0 at z = z1is imposed initially. Using eqns(4 and 20) in eqn (22), we obtain: GYZ

%lk,k+l

- (z- Zk)b~,, +&I I +4.

(25)

After lengthy algebra,it can be shown that the secondboundary condition gets satisfied automatically, when the governingequation of equilibrium of the laminated platetheoryis exactlysatisfied.However,in the finite elementform, the secondboundary condition cannotbe satisfiedexactlydue to discretizationerrors. The derivativeof strainswhich occur in eqns(23-25)are obtainedin the following ways: (i) direct differentiation of the shapefunctions describingthe displacementfield; (ii) a finite difference scheme based on first and second order Taylor series expansion. 1. Direct differentiation

of the shape functions

(DDSF)

Using eqns (4-7), the strains and curvaturescan be expressedin terms of displacementsandtheir derivatives.Thesedisplacementscanbe written in termsof the shape functions describing the displacementfield in an element [eqn (l)]. Higher order derivativesof strainsandcurvaturesareobtainedby differentiatingtheseshapefunctions

Finite element estimation of elastic interlaminar

stresses

451

s-s S-S 0

S-S /B*

II "I-=i

Fig. 3. The finite element mesh in a quarter of a cross-ply square plate with simply supported edges.

and calculatingtheir valuesat a given location inside an element,specifiedby its area coordinates< and q within an element.Basedon Chaudhuriand Seide’s(1987)work, in which it was indicated that the centroid of the triangular element is the point of exceptional accuracy, the stressesare calculatedat the centroid of eachelement. 2. Finite difference scheme (FD)

A finite differencegrid is generatedtreatingthe nodesof the finite elementmeshas the grid points (seeFig. 3). The in-planestrainsandcurvaturesarecalculatedat eachnode in the grid using the strain-displacementrelations. The higher order derivativesof the strainsand curvaturesarecalculatedusingfirst and secondorderTaylor seriesexpansion. Dependingon the location of the grid point, three difference schemesare used viz. forward difference, central difference and backward difference. Forward differenceis usedin the left and bottom boundaryof the domain, backwarddifferenceis usedin the top and right boundaryof the domain and centraldifferenceis usedin the interior of the domain, Lajzcok (1986).It has been found that secondorder terms are necessaryto obtain accurateresultsat the boundaryand first order terms areadequatein the interior. The relevant expressionsof the three difference schemesare given in Appendix III. Thoughthe numericalschememight presentdifficulties in obtaining higherderivativesin complicated geometries,it is believedthat thesedifficulties could be overcomeby a combination of establishinga convenientcoordinatesystemand a transformation to a Cartesiancoordinatesystem(cf. computationalfluid mechanicsproblems). STRAIN ESTIMATION

Classicallaminated plate theory assumesthat the transversestrains(E,,, e,,,and e& are equalto zero. While this assumptionsimplifies the theory and makesit attractivefor a preliminary global analysis,an estimationof interlaminar strains,as indicatedearlier, is necessary. Thesestrainscan be computedusing the constitutiverelationsbetweenstrainsand stresses.Estimated values of interlaminar stresses(by integrating ply equilibrium equations)and in-planestressesareusedto estimateinterlaminarstrains.The compliance

R. GOVINDARAJANet al.

458

matrix, which relatesthe strainsto stressesfor an orthotropic material, is given by:

Sl, s,, s,, 0 0 0

s,2 s,, $2 0 0 0

$3 0 g3 0 s33 0 0 s// 0 0 0 0

0 0 0 0 0 0 0 0 $3, 0 0 s,,

sti is the transformedcompliancematrix eqn (24). NUMERICAL

I(

(26)

RESULTS AND DISCUSSIONS

Numerical results have been obtained using the presentapproach for a simply supportedsquareplate of sidelength~1,subjectedto transversesinusoidalloading: q

= q,sinysinT.

(27)

A four-layered symmetric cross-ply laminate [O/B], and a two-layeredasymmetric cross-plylaminate [O/90]are consideredwith properties: E, = 25 x lO,psi(174.6GPa);

G,2 = G13= 0.5 x lo6psi (3.5GPa); Vl2

=

V23

=

013

Ez = Es = lo6psi (7 GPa)

G,, = 0.2

x

lo6 psi (1.4GPa)

(28)

- 0.25.

The interlaminarstressesarecomparedwith the valuesof CLPT analyticalsolutions. The resultsare normalizedas follows: axz=

0 qoS cos(m/Z)

sin(ny/b) ’

cr,

a, = q0S2

sin(nx/a) sin(ny/b) ’

& “’ = S sin(nx/~~cos(rry/b)’

0 oyz = q. S sin(rrx/z cos(ny/b) & eXz= S cos(nx/0~sin(ny/b) &Z

” = S2sin(nx/a) sin(lry/b)

(29) (30)

Case 1

The valuesof interlaminar stressesat midplane for symmetriccross-plysquareplate are given in Table 1. The thicknesswisevariation of interlaminar stressesis shown in Figs 4-6. It can be noted from theseresultsthat resultsof DDSF are closeenoughto theoreticalvaluesonly in the caseof interlaminarshearstresses,whereasresultsof FD are very accurate.For examplein the value aXz(for z = 0) the percentageerror of DDSF is 4.13 whereasFD resultsare accurateup to threesignificant digits. Similarly in the value BYE(for z = 0) the percentageerror is 2.9 in caseof DDSF and is 1.45in caseof FD. The satisfactionof the zeroshearboundarycondition at the top surfaceof the laminate (at z = z”+~, oxz= vyz = 0) is automaticallyachieved.However,DDSF schemedoesnot satisfy the top surfaceboundary condition for a, (at z = z,,+~,o, = q). The error is believedto be due to discretization.While the local equilibrium equationsare satisfied through integration, the laminate governing differential equations is only satisfied approximatelyin finite elementmethod. It is believedthat this error can greatly be reducedby refining the meshfurther. However,this erroneousstressdistribution can be

Finite elementestimationof elasticinterlaminar stresses 0.61 ---

CLPT Analytical FEM /DDSF

459

8 FEM /FD

( W~l~)Laminate

I

Or

-0.6

i 1 I 0.4 Ii =

I 0.2

I

I 0.6

1

I

Fig. 4. Variation of a;* through the thickness.

0

O-08

0.16

0.24

%

Fig. 5. Variation of a& through the thickness.

0.1

0.2 Fz

0.3

x lo+

Fig. 6. Variation of u; through the thickness.

corrected with referenceto known surface values (thus ensuring the satisfaction of laminate governingdifferential equations)along with an interpolation at ply level. This modification producesaccuratevaluesof oz. The resultsby FD schemefor a, distribution is quite accurate without any modification as mentioned previously. But the results becomeevenbetterwith modification (seeTable 1). Case2 The valuesof interlaminar stressesat midplanefor asymmetriccross-plysquareplate are given in Table 2. The thicknesswisevariation of interlaminar stressesare shown in Figs 7-8. It hasbeenfound that resultsby DDSF schemearevery erroneous.Becauseof

R. GOVMDARAJAN et al.

460

the asymmetry in the lay-up, the bending-stretchingcoupling terms exist and hence in-plane displacementshave nonzero values. Thus both in-plane displacementsand out-of-planedisplacementscontributeto the error. However,the resultsby FD scheme areaccuratefor interlaminar shearstressesand for normal stress(seeFig. 8) the error in the unmodified variation is less,but with modification the resultsbecomevery accurate (seeTable 2). Table 1. a,, and ayz at Z = 0.0, for symmetric cross-ply square plate

0x2 QYZ

uz

CLPT Analytical

FEM FD

0.122(0.1221)* 0.122(0.1221)* (OS)*

0.118 0.118 0.500(0.486)*+

*Exact elasticity solution for S = 50. **Unmodified values. Table 2. a,, and a,, at Z = 0.0, for asymmetric cross-ply square plate FEM

uxz QYZ

%

CLPT Analytical

DDSF

FD

0.339(0.328)* 0.138(0.156)* -

0.353 0.134 0.125(6.213)**

0.339 0.136 0.125(0.117)**

* Exact elasticity solution for S = 50. ** Unmodified values. ----

CLPT FEM

Analytical /FD

-0.4 -0.2 5 0 o-2 0.4

Fig. 7. Variation of uiz and uiz through the thickness.

( [O/901

) Lominot

e

Fig. 8. Variation of u; through the thickness.

Finite element estimation of elastic interlaminar

stresses

461

CONCLUSIONS

A high-precisiontriangular element,basedon the classicallaminatedplatetheoryhas been adapted. Interlaminar stressesand strains are estimated using ply equilibrium equationsand stress-strainrelationsrespectively.The higherorderderivativesof in-plane strainsandcurvaturesrequiredfor the calculationof interlaminar stressesareobtainedby two schemes,namely, direct differentiation of the shapefunctions (DDSF) and a finite difference scheme (FD). By comparing the interlaminar results for symmetric and asymmetriccross-plylaminatesobtainedby DDSF and FD, it is found that the latter is more accurate.The first approachmay alsobe expectedto giveaccurateresults,provided a very refined meshis usedfor finite elementanalysis. The approachdiscussedand presentedin this paper is perhapsthe simplest of all currently availablemethodsand henceis attractiveas a global methodin the global local analysisof compositelaminates. Acknowledgements-This work was supported by Aeronautics Research and Development Board, Ministry of Defence, Government of India, underproject no. 604. Also, it is a pleasure to acknowledge the discussions with Mr H. K. Harikumar. REFERENCES Chaudhuri, R. A. and Seide, P. (1987). An approximate semi-analytical method for prediction of interlaminar stresses in an arbitrarily laminated thick plate. Comput. Struct. 25, 627-636. Dong, S. B., Pister, K. S. and Taylor, R. L. (1962). On the theory of laminated anisotropic shells and plates.

J. AerospaceSci.29, 969-975. Engblom, J. J. and Ochoa, 0. 0. (1985). Through the thickness stress prediction for laminated plates of advanced composite materials. Znt. J. numer. Meth. Engng 21, 1759-1776. Harikumar, H. K. and Krishna Murty, A. V. (1991). On modelling of delaminations in symmetric composite laminates. Comp. Sci. Technol. 42, 393-411. Hayashi, T. (1967). Analytical study of interlaminar shear stressesin a laminated composite plate. Trans.Japan

Sot. Aeronaut. SpaceSci. 10, 43-48.

Jayachandrabose, C. and Kirkhope, J. (1985a). Explicit formulation of a high precision triangular laminated anisotropic thin plate finite element. Comput. Struct. 20, 991-1007. Jayachandrabose, C. and Kirkhope, J. (1985b). A high precision triangular laminated anisotropic shallow thin shell finite element. Comput. Struct. 21, 701-723. Kant, T. and Pandya, B. N. (1988). A simple finite element formulation of a higher order theory for unsymmetrically laminated composite plate. Compos. Struct. 9, 215-246. Krishna Murty, A. V. and Vellaichamy, S. (1988). Higher-order theory of homogeneous plate flexure. AZAA J. 26, 719-725. Lajczok. M. R. (1986). New approach in the determination of interlaminar shear stresses from the results of MSC/NASTRAN. Comput. Struct. 24, 651-656. Lo, K. H., Chistenson, R. M. and Wu, E. M. (1977). A high-order theory of plate deformation, part I: Homogeneous plates, part II: Laminated plates. J. appl. Mech. 44, 663-676. Mao, K. M. and Sun, C. T. (1992). Error estimators using global-local methods. Znt. J. numer. Meth. Engng 35, 589-599. Murthy, M. V. V. (1981). An Zmproved TransverseShearDeformation Theoryfor Laminated Anisotropic Plates.NASA technical paper, 1903. Nelson, R. B. and Larch, D. R. (1974). A refined theory of laminated orthotropic plates. J. appl. Mech. 41, 171-183. Owen, D. R. J. and Li, Z. H. (1987). A refined analysis of laminated plates by finite element displacement methods-I. Fundamentals and static analysis. Comput. Struct. 26, 907-914. Pagano, N. J. (1970). Exact solutions for rectangular bidirectional composites and sandwich plates. J. Camp.

Mater. 4, 20-34.

Pagano, N. J. and Hatfield, S. J. (1972). Elastic behavior of multilayered bidirectional composites. AZAA J. 10, 931-933. Phan, N. D. and Reddy, J. N. (1985). Analysis of laminated composite plates using a higher-order shear deformation theory. Znt. J. numer. Meth. Enann 21. 2201-2219. Pryor Jr, C. W. and Barker, R. M. (1971). A-finiteelement analysis including transverse shear effects for applications to laminated plates. AZAA J. 9, 912-917. Ransom, J. B. and Knight, N. F. Jr (1990). Global/local stress analysis of composite panels. Comput. Struct. 31, 375-395. Reddy, J. N. (1984). A simple higher-order theory for laminated composite plates. J. appl. Mech. 51.745-752. Reddy, J. N. (1989). On refined computational models of composite laminates. Znt. J. numer. Me&. Engng 27, 361-382. Reddy, J. N., Barbero, E. J. and Teply, J. L. (1989). A plate bending element based on a generalized laminate plate theory. Znt. J. numer. Meth. Engng 28, 2275-2292. Reissner, E. and Stavsky, Y. (1961). Bending and stretching of certain types of aeolotropic elastic plates. J. appl.

Mech. 28, 402-408.

Shah, C. G. and Krishna Murty. A. V. (1991). Analysis of edge delaminations in laminates through combined use of quasi-three dimensional, eight-noded, two-noded and transition elements. Comput. Struct. 39, 231-242.

R. GO~INDARUAN et al.

462

Spilker, R. L. (1982). Hybrid-stress eight node elements for thin and thick multilayer laminated plates. ht. J. numer. Meth. Engng 18, 801-828. Srinivas, S. and Rae, A. K. (1970). Rending, vibration and buckling of simply supported thick orthotropic rectangular plates and laminates. Int. J. Soli& Structs 6, 1463-1481. Stroud, W. .I., Housner, J. M., Tanner, J. A. and Hayduk, R. .I. (Editors), (1985). Local/global nonlinear stress analysis. Proc. Session I Workshop on Computational Methods for Structural Mechanics and Dynamics, NASA, Langley Research Center, Hampton, VA, 19-21 June. Viddussoni, A., Muheim Thompson, D. and Hayden GrifiIn. 0.. Jr (1990). Cross-ply laminates with holes in compression: straight free-edge stresses determined by two- to three-dimensional global/local ftite element analysis. J. Compos. Technol. Res. 20, 209-216. Whitney, J. M. (1969). Rending extensional coupling in laminated plates under transverse loading. J. Camp. Mater. 3.20-28. Whitney, J. M. and Leissa, A. W. (1969). Analysis of hetrogeneous anisotropic plates. J. appl. Mech. 36, 261-266. Whitney, J. M. and Sun, C. T. (1973). A higher order theory for extensional motion of laminated composites. J. Sound Vibr. 30,85-97. APPENDIX

I

R-Matrix 00000010

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

b,

0

0

0

3

-c,

0

-c,

b,

0

-3

-2q

24

0

f;

f;

-13

f;

f9’

27

cl

-bl

-3

2c,

-2b1

0

c2 -b,

0

0

0

2

cz

-4

0

7

f,”

f!

7

fs”

f:

13

fe8

f;

-27

7

f2’ 0

f,’ 0

7

f,’ -cl

f,’ b,

13

fl

f;

-27

2

-q

0 -2

-0

fi 0

-7

-4

3

-7

f;

c,

0

0

-2

b,

b, = Yj - Yk c, = x, - xj, wherei=l,2,3;j=2,3,l;k=3,1,2.

f2' = f; = f2” = f; =

c, + 2c2

f; = b, - 2bz

f; = -c, - 2c,

bz + 26,

fe! = 3(c, - c2) f,” = 2(b, + b3 f; = 3c, - 2c,

f; = -c, - 2cz

f,’ = b, + 2bz

19’ = f,” = f; = f,’ =

fl = -2(b, + 9)

f; = 2~ - 3c,

f9’ = -2b, + 3b,

-2(c, + CJ -b, - 26,

Xbz - 4) c, + 2c, -39

+ 2b,

2(c, + CJ

0

where

i =

1, 2. 3;j

0

- 3c,

g;’

A”

6

-3od,

d”

g:'

0

0

0

0

= 2, 3, 1; k = 3. 1. 2.

g,

0

0

0

0

0

0

0

0

0

d'

$"

gi'

e3

0

0

$'

id2

-YS

0

0

g6

8

iY3

0

0

0

0

0

0

bt + ct

;

-

6

- 30d2 4%

&

&

4:

xj,

Yk

0

0

dd

g::

0

0

+xX,

do

0

0

0

0

0

0

0

0

0

0

*x2

‘4:

g::

s::

-x2

g::

g::

0

0

+x2

ST*

0

0

0

0

0

0

0

0

bf+cf

b, x b2 + c, x c,

-

Yj x,

=

c,

-3b,

do

&’

d”

-x,

bi =

d, =

3c,

ido

d’

d”

18 g7 d’

0

0 0

0

0

0

g? 76,

gY

d’

0

0

- 4b,

d

0

0

0

0

0

0

0

0

7c,

-15

m2

g7

A’

0

0 13

0

k,

z

0

10

30d2

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

b, x b, + c, x c,

0

0

g:'

g:*

d'

iv2

0

0

g:'

d"

0

0

gs

iY2 8

-y2

d, = -

0

0

g:'

g:*

g:"

d’

id'

G*

0

0

id'

d'

-Y,

0

0

g4

iY,

g:'

34

0

0

g:)

g:"

0

19

0

d’

0

13

0

g,

d'

6od,

0

0

0

-74

0

0

7c*

g3

Lz;

0

-15

-3w3

10 8

0

0

0

4th 8

0

0

0

0

-4c*

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

;

+x3

d3 =

3c,

-3b,

$4

bf+cf

-

-+ ‘x L

&

g:: gl619

g::

s::

-iY,

h,

&

-r ‘x

& 2

19 g17

ST:",

4:

-iY2

+x2

g :*,

g::

12 g17 g::

SY2 gl612

- *x2

id7

d7

-$Y2

ix2

-P2

iY2

0

0

0

G,

- ;x,

d6

d6

-G,

h,

-P,

iY,

0

0

0

g::

g::

3bz

b, x b2 + c, x c,

-6

&

g:‘: 3od, &

g14 19 g14

id: 18

- 3c2

84

s::

- 8c,

4:

g::

- 8b2

-6b,

I&

d5

14 g14

g::

3od,

-6

15

0

0

0

b,

662

-b,

0

13 g14

12 g14

8~2

%

d4

,9 g12

0

0

-6od,

g::

- 60d3

15

8

0

0

0

g14

-SC2

-Cl

c2

0

g::

-x3

d’2

4:

0

0

ix3

3% - 10

3od,

9 g12

- 10

0

0

0

0

0

0

0

0

-r ‘X 3

&

g::

4:

4:

- L3

*x3

&

A’8

12 gl8

G3

- %x3

d8

g:*

-G,

ix3

-P3

L3

0

0

0

R. GOVINDARAJANet al.

464 x, = c:

Yl = c:

z, =c:

x2 = -2b,c,

y2 = -2bZc2

z2 = -2b3c,

~2 = -(b,cz

x3 = bf

Y, = b:

z3 =

~3 = b,b2

g: = 15c,d3 - 3c,

g; = 15b2d3

gi = ~2 - jy,d,

$

g$

= -15bld3

da

= bA

do

= P,

- 3b,

g;

+ 2~3

g;,

= &,

g;’

= -2b,

g::

= -P,

= -2b,

=

$’

= 5yy,d,

- PP

g;:

= -30c2d3

g::

= -p2

- 2q

-

g:;

30b,d2

82

= ~3 - jx3

&

= b,

jx,4

g:”

= 2c, -

30c2d3

-

2c,

- 5y2d3 + 2b2

d’

= 5r2d3

- ~2

g;:

=30b2d3

+ 2b,

g:;

= -P,

-

5&

gt’o = 5x,d2

- PI

- ~3

g;:

= 2c2 + 30cld2

- 5x,d2

g;:

= -p2

- 5x2d2

g;’

= 15c2d3

g;’

= -15b2d3

= hh

4’

= j~3d3

= -15b2d3

d:

= fy,d,

s:: = by24

g::

= hh

do = -15c,d,

g”

= 15b,d2

do, = -jx,d,

d; = -h%

82

=

-+x+4

di = -15c,d,

g$

=

&

=

h4

d’: = h-4

d;

= hh

- lSb,d,

g;’

= -jz,d,

+ 5y,d,

-

3p,

- x2

g;’

= -15b3d, = f(x2

- 3Ob,d, - z&

- 6b2 -

+ jv2

lJb,

+ 2~2 + 5x24

g;; = -30[1 + 2(d2 + d3)] = -6(b,

- b2) - 30(b,d,

-Y&

+

;” = -3O(d,

54)

-

- b,d,)

2p,

-

x2(+

+

+ 15b3d,

8:”

= 3P2 + jx2

18 87

= 3O(d,

89‘*

= l5(b,d,

sit

= jz2dl

-

- hd, d2 -

s::

= -94

+ 66,

s::

= ~2(j

+ 54)

-3W

gy

E

W,d,

gj'

=

&2

g;’

= 3O(d,

+

= - jz,d,

+ 5y3d3

- 3~3 - x3

+ 9b2

= j(x,

- z,d,)

+ +Y, + 2p,

=

- z,d,)

+ ju3

g;:

= 6(cl

g::

= -Y&

+ Sd,)

g::

= -u,(j

+

g;”

= 3tk2d3

g;: - 3~2

+ 15b,d2

+ $dJ

+ 6c2 + 15c,

g::

&x3

-

jzzJ d,

-

+ 54

Zd,) - 2~33 - s(j

+ 54)

-

- tp,

c2dJ

- x,6

15c3d,

+ jx,

= jz,d,

+ 5x14

+ 2~3 + 5x,4

c2) + 30(c1d2

-

5yA

- x,

g::

= 3p,

30b2d3

+ x2(1

- !k,

- jz, 4

- 5~1d3

- jzt,d,

-

+ c,d2

5y,d,

+ c,)

- 9c2

- jx,(d,

+ 1) - ju,

-

-

+

- 3P3

jx,(d,

d”2

=

g;:

= 9c2 - 6c, + 30c2d3

1) -

b3

15c,d2

-

s:: = ~18 + 54) + x,(1 + fds, g::

= u3(j

+ 54)

+ x3(1

+ fdz,

g;9 = 15(c,d, - c3d1) - 6c,

4)

~2 -

g;’

g:’ = -15(c,d,

- b2d3) +

- 3p,

= 3~3 + jx3

+ 1) - jvz

-

+ 5y,d3

g:’

+ b,)

- $x,(4

= -jz,d,

g:”

1)

+ b,d2

g;’

+ 9b,

s:: = 30(1 + dd + 6Od,

19 81 =

54)

+ 2d3)

= -30b2d3

- 30c2d3

g:’ = 15c,d, + 30cld2

g;’ = 30(1 - d, + 2d2)

s::

15bld2

g;’ = 9c, + 15qd, + 3Ob,d,

4

+ h%

g;;

= -9b,

2~32

15c,d2

d’

g;’

g;*

+

+ 15b1d2

= $Y,&

d’

8

$y,d,

g; = -34

- P,

= 5x34

g;“, = -p,

=

+ 3c,

g;: = 15c2d3

+ 2d3)

817 =

2~33 +

gi4 = -3Ob,d,

- 5x3d2

= 3O(d,

13

+

da

- ~2

g;‘8 = -p3

g;;

9:’

lSc,d,

+ 3Ob,d,

- 5y,d,

g;;

+ 2p,

= 3b,

gi4 = 30c,d2 = 5x24

= h4 = 3c2 -

= ~2 - jx2d2

= 5~34

= P, - jy,d,

gf4 = 15c2d3

gt

g:’

d

+ 4c,)

= ~3 - h4

8,

15c,d,

= c,c2

+ 3b,

-

+ +xd,

g::

b;

- jx,d,

gt4 = -3c,

d’

do

PI

+ 6b1

fd,z2

- 24f2 -

-

I)

j4yz

gi9 = Zp,

+ x,

-

jd,z,

d’

=

+ x3 -

jd,z3

g:’

= -15(c,d,

?p3

+ cl) -

- jd,y, -

j4n

30cld2

- 6~2

3P1

Finite element estimationof elasticinterlaminar stresses

465

g19 9 = lS(b,d, + b,) + 3Ob,d, + 6b 2

s:‘a = &z,d, - x,) - 5x14 - ~1 - 2pl

g:; = #(z&l -

g:; = 30(1 + d,) + 60d2

g:: = $63 4 - x3) - 5x,4 - ~3 g;: = 6c2 - 9c, + 1Sc2d3- 30c,d2

g;; = -6b2 + 96, - 15b2d3 + 30b,d2

g:: = Y# + $4) + x,(t + 54)

d’, = M

x2)

- 5x24 -

+ td& +

x2($

~2

+ 54)

-

2~2

g:; = ~30 + td,) + APPENDIX II

p=

CaE335-f

x3(+

+ 54)

2p3

466

R. GOWNDAMJAN et 01. APPENDIX a.-zta.-1ra”.a.+1~a.+2

III

are the coordinates of successive nodes in the a-direction (see Fig. 3). a denotes

a typical direction x, y, r or 0. ha = am+2 - a.+, = a,+, - a,. Let 6 be a typical quantity such as eX. The value of derivative of + with respect to a at node n, can be obtained by using any of the following finite difference schemes: (i) forward difference (from second order Taylor series expansion)

#,.,W =

-9(a,+3

+ Wa,,,) 2Aa

- %(a,) ;

(ii) backward difference (from second order Taylor series expansion) 4

( ) = 4ta,d rO an

- Ma,-,)

+ 34taJ

2Aa

(iii) central difference (from first order Taylor series expansion)

&,,(a,) =

Ha,+, 1 - 4ta.J 2Aa



;