‘%“
[email protected]& vol. 3. No 5, pp. 451466. Printed in Great Britain.
09619526/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) AbstractLow 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 inplane 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 threelayered symmetric and a twolayered 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 fiberreinforcedplasticcompositelaminatesis 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.Suchgloballocal 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 threedimensionalfinite elementschemesto obtain interlaminar strains. The use of threedimensionalfinite element approach, is, perhaps, prohibitive for this purposeas it involves a large computational effort. A searchfor simpler approachesthan the threedimensionalfinite 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,severaltwodimensionalfinite 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 displacementbasedfinite 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 straindisplacement
=
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 bendingmembranecoupling 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(47)] 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 1cq 0 0 [Aj] = 0 (I 0 0 tf 0 0 1rq 0 (13) 0 0 1ta [ 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 inplane strains are computed using the straindisplacementrelations [eqns(47)]. The inplanestressesa,, 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(2325)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 (47), 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
ss SS 0
SS /B*
II "I=i
Fig. 3. The finite element mesh in a quarter of a crossply 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 inplanestrainsandcurvaturesarecalculatedat eachnode in the grid using the straindisplacementrelations. 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 inplanestressesareusedto 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 fourlayered symmetric crossply laminate [O/B], and a twolayeredasymmetric crossplylaminate [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 symmetriccrossplysquareplate are given in Table 1. The thicknesswisevariation of interlaminar stressesis shown in Figs 46. 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
O08
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 asymmetriccrossplysquareplate are given in Table 2. The thicknesswisevariation of interlaminar stressesare shown in Figs 78. It hasbeenfound that resultsby DDSF schemearevery erroneous.Becauseof
R. GOVMDARAJAN et al.
460
the asymmetry in the layup, the bendingstretchingcoupling terms exist and hence inplane displacementshave nonzero values. Thus both inplane displacementsand outofplanedisplacementscontributeto 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 crossply 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 crossply 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 o2 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 highprecisiontriangular element,basedon the classicallaminatedplatetheoryhas been adapted. Interlaminar stressesand strains are estimated using ply equilibrium equationsand stressstrainrelationsrespectively.The higherorderderivativesof inplane 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 asymmetriccrossplylaminatesobtainedby 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. AcknowledgementsThis 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 semianalytical method for prediction of interlaminar stresses in an arbitrarily laminated thick plate. Comput. Struct. 25, 627636. Dong, S. B., Pister, K. S. and Taylor, R. L. (1962). On the theory of laminated anisotropic shells and plates.
J. AerospaceSci.29, 969975. 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, 17591776. Harikumar, H. K. and Krishna Murty, A. V. (1991). On modelling of delaminations in symmetric composite laminates. Comp. Sci. Technol. 42, 393411. Hayashi, T. (1967). Analytical study of interlaminar shear stressesin a laminated composite plate. Trans.Japan
Sot. Aeronaut. SpaceSci. 10, 4348.
Jayachandrabose, C. and Kirkhope, J. (1985a). Explicit formulation of a high precision triangular laminated anisotropic thin plate finite element. Comput. Struct. 20, 9911007. Jayachandrabose, C. and Kirkhope, J. (1985b). A high precision triangular laminated anisotropic shallow thin shell finite element. Comput. Struct. 21, 701723. 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, 215246. Krishna Murty, A. V. and Vellaichamy, S. (1988). Higherorder theory of homogeneous plate flexure. AZAA J. 26, 719725. Lajczok. M. R. (1986). New approach in the determination of interlaminar shear stresses from the results of MSC/NASTRAN. Comput. Struct. 24, 651656. Lo, K. H., Chistenson, R. M. and Wu, E. M. (1977). A highorder theory of plate deformation, part I: Homogeneous plates, part II: Laminated plates. J. appl. Mech. 44, 663676. Mao, K. M. and Sun, C. T. (1992). Error estimators using globallocal methods. Znt. J. numer. Meth. Engng 35, 589599. 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, 171183. Owen, D. R. J. and Li, Z. H. (1987). A refined analysis of laminated plates by finite element displacement methodsI. Fundamentals and static analysis. Comput. Struct. 26, 907914. Pagano, N. J. (1970). Exact solutions for rectangular bidirectional composites and sandwich plates. J. Camp.
Mater. 4, 2034.
Pagano, N. J. and Hatfield, S. J. (1972). Elastic behavior of multilayered bidirectional composites. AZAA J. 10, 931933. Phan, N. D. and Reddy, J. N. (1985). Analysis of laminated composite plates using a higherorder shear deformation theory. Znt. J. numer. Meth. Enann 21. 22012219. Pryor Jr, C. W. and Barker, R. M. (1971). Afiniteelement analysis including transverse shear effects for applications to laminated plates. AZAA J. 9, 912917. Ransom, J. B. and Knight, N. F. Jr (1990). Global/local stress analysis of composite panels. Comput. Struct. 31, 375395. Reddy, J. N. (1984). A simple higherorder theory for laminated composite plates. J. appl. Mech. 51.745752. Reddy, J. N. (1989). On refined computational models of composite laminates. Znt. J. numer. Me&. Engng 27, 361382. 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, 22752292. Reissner, E. and Stavsky, Y. (1961). Bending and stretching of certain types of aeolotropic elastic plates. J. appl.
Mech. 28, 402408.
Shah, C. G. and Krishna Murty. A. V. (1991). Analysis of edge delaminations in laminates through combined use of quasithree dimensional, eightnoded, twonoded and transition elements. Comput. Struct. 39, 231242.
R. GO~INDARUAN et al.
462
Spilker, R. L. (1982). Hybridstress eight node elements for thin and thick multilayer laminated plates. ht. J. numer. Meth. Engng 18, 801828. 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, 14631481. 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, 1921 June. Viddussoni, A., Muheim Thompson, D. and Hayden GrifiIn. 0.. Jr (1990). Crossply laminates with holes in compression: straight freeedge stresses determined by two to threedimensional global/local ftite element analysis. J. Compos. Technol. Res. 20, 209216. Whitney, J. M. (1969). Rending extensional coupling in laminated plates under transverse loading. J. Camp. Mater. 3.2028. Whitney, J. M. and Leissa, A. W. (1969). Analysis of hetrogeneous anisotropic plates. J. appl. Mech. 36, 261266. Whitney, J. M. and Sun, C. T. (1973). A higher order theory for extensional motion of laminated composites. J. Sound Vibr. 30,8597. APPENDIX
I
RMatrix 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’: = h4
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=
CaE335f
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 adirection (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
’
;