An approximate semi-analytical method for prediction of interlaminar shear stresses in an arbitrarily laminated thick plate

CkqmfrrsA Stnicfurts Voi. 23.No. 4. pp.627-636.1987 Printed in Gmi B&in.

AN APPROXIMATE SEMI-ANALYTICAL METHOD FOR PREDICTION OF TNTERLAMTNAR SHEAR STRESSES IN AN ARBITRARILY LAMINATED THICK PLATE &AZ A.

&AuDHlJRl

Department of Civil Engineering, Unjve~ity of Utah, Salt Lake City, UT 84112, U.S.A.

and PAUL SEIDE

Department of Civil Engineering, University of Southern California, Los Angeles, CA 90089-I 114, U.S.A. (Received 19 March 1986)

Abstract-An approximate semi-analytical method for determination of interlaminar shear stress distribution through the thickness of an arbitrarily laminated thick plate has been presented. The method is based on the assumptions of transverse inextensibility and layezwisecunstant shear angle theory (LCST) and utiiixes an assumed quadratic displacement potential energy based finite eiement method (FEM). Centroid of the triangular surface has been proved, from a rigorous methematical point of view ~Aubi~Ni~e theory), to be the point of exceptional accuracy for the ~~r~~r shear stresses. Numerical results indite close agreement with the avaiIable ~~irn~sio~al elasticity theory solutions. A comparison betvveen the present theory and that due to an assumed stress hybrid FEM suggests that the (normal) traction-free-e&e condition is not satisfied in the latter approach. Furthermore, the present paper is the first to present the results for interIaminar shear stresses in a two-layer thick square plate of balanced unsymmetric angle-ply construction. A comparison with the recently proposed Equilibrium Method (EM) indicates the superiority of the present method. because the latter assures faster convergence as well as simultaneous vanishing of the transverse shear stresses on both the exposed surfaces of the laminate. Superiority of the present method over the EM, in the case of a symmetric laminate, is limited to faster convergence alone. It has also been demonstrated that the combination of the present method and the reduced (quadratic order) numerical integration scheme yields convergence of the interlaminar shear stresses almost as rapidly as that of the nodal displacements, in the case of a thin plate.

a$), a;‘,

length of a square plate or width of an infinite strip strain+iisplacement relation matrix elastic (material) stiffness matrix of an anisotropic material nodal displacement vector for t.hejU triangular layer-element belonging to the i” layer mesh size for a regular or uniform mesh rlh component of unit normal vector P component of unit vector normal to the P side of the triangular element total number of layers number of subdivision applied surface load area of t&j* triangular layer eIement spa= of trial (~spla~cnt) functions total thickness of a laminated plate thickness of the P layer displaament vector displacement components in x, y, z directions, respectively displacement components in x and y directions respcaively at the i” interfaac of the triangular element vector of trial displacement functions Cartesian coordinatbz 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 the k* side of the j* element

(0 f(i) )I

7x??

7 !$

inplane normal and shear stresses in the i* layer transvcm shear stresses in the P kyer INTRODUCI'ION

Study of the distribution of interlaminar (transverse) shear stresses through the thickness of a laminated plate has assumed increasing importance because of the role these stresses play in causing delamination failures in composites. A brief look at the literature, however, suggests that with the exception of a few analytical solutions@-31 which are restricted to problems with simple shapes, loadings, boundary conditions, and cross-ply lamination, the problem of accurate prediction of in~ri~i~r shear stresses in an arbitrariiy Lninated anisotropic plate has received little attention. The finite element method (FEM) seems to be the only practical alternative, because of the ease with which irregular shapes, anisotropy, complex boundary conditions and general loadings can be handled by this method. A review of the FEM literature, however, suggests that short of conducting a highly refined threedimensional analysis (e.g. Barker et ~2.[4] wherein as many as three elements per layer have been used), few FEM-based methods are available which yield accu-

627

628

I&AZA. CHAUDHURI and PAUL SUDE

rate interlaminar shear stress variation through the thickness of an arbitrarily laminated plate. The quadrilateral element due to Mau er al. [15], based on assumed stress hybrid approach, seems to be an exception to this. This method considers the stresses as unknown nodal parameters and imposes constraints on the compatibility of these stresses at each interface and on vanishing of the interlaminar shear stresses on the top and bottom surfaces of the plate. This method, however, appears to be limited in its applicability, since the formulation of the element stiffness matrix involves too many matrix inversions at the element level and extension by Spilker et ai. [4] to triangular element shape has produced a very stiff element. More impo~antly, (normal) traction-freeedge formulation has not been sumssfu1[6], which leads to inaccuracy in the prediction of the transverse shear stresses near the edge [7]. Furthermore, a review of the literature suggests that this method has not been applied in the prediction of the interlaminar shear stresses in an asymmetrically laminated plate. A recent paper [7] has presented a method for prediction of interlaminar shear stress distribution through the thickness of a thick laminated plate. This method is based on finding the inplane stresses first and then util~g the first two equations of equiIibrium in terms of stresses-an idea originally due to Pryor and Barker [S], who introduced it for a quadrilateral symmetrically laminated plate element based on CST (Constant Shear Angle Theory). The triangular finite element formulation used by [7J utilizes an assumed quadratic displacement potential energy approach [9] and is based on the assumptions of transve.rse inextensibility and LCST (Layerwise Constant Shear Angie Theory, Fig. 1). The success of the Equilibrium Method (henceforth, referred to as EM), as presented in 171,hinges on either ~intwise fin the thr~-dimensional sense) of the second derivatives of convergence displacements or inplane force equilibrium {i.e. the through-thickness integrated form of the stressequilibrium equations at a point (x, JJ) on the reference surface). However, the fact remains that the first derivatives of displacements, as computed by FEM, converge in the mean-square sense [9, IO] while the

Fig- 2. LCST based composite plate element.

pointwise inplane force ~~lib~urn is generatly sat&% in case of homogeneous and symmet~cally laminated plates, which are characterized by the absence of mid-plane stretching. For asymmetrically laminated plates-even though it is possible that the EM may at times yield reasonably accurate resultsthe success of this method appears to be at best uncertain. The objective of the present paper is to alleviate this difficulty and develop a semi-analytical approximate method for prediction of the interlaminar shear stress distribution through the thickness of an arbitrarily laminated anisotropic plate.

Each triangular layer element belonging to the 2” layer is bounded by the i” and i f I* interfaces (Fig. 2). Each interface triangle is characterized by six nodes. The assumption of layerwise constant shear angle guarantees piecewise linear variation of inplane displacements through the plate thickness. By virtue of the assumption of transverse inextensibility, the transverse displacement does not vary through the thickness. The number of degrees of freedom per node is then 2N + 3 for an N-layer composite piate element. Details concerning formulation of element stiffness matrix and consistent load vector, effect of numerical integration on the convergence of displacements and stresses etc. are available in [9,11]. Once the displacements are determined, the element stresses can be obtained using the relation [l I]: {a”‘(z)}=[C~~~][~~‘~(z)][E:‘~](d~~),

(la)

where (C(i) (z))T= (a:‘(%), @F’(z), ,x 7~~~z),~~{z),i~~~z)~

Fig. 1. Possible severe cross-sectional warping in a thick laminated plate.

(lb)

and [PJ, [Lo], (df’l] and [Ej’)] are given in the Appendix [eqns (A3); (AGA7); and (A8) and (A9-A13), respectively]. It may be noted that while acazurate inplane stresses can be computed at the mid-points of the sides of the, triangular element interface, only through-layer-thickness average transverse shear stresses, r’zj and “Fj, are obtained using eqn (1). The transverse shear stresses, as computed by the

Prediction of interlaminar shear stresses EM of [7] will not, in general, simultaneously vanish at both the top and bottom surfaces of an asymmetrically laminated plate. The present paper, therefore, represents the first FEM-based attempt at computation of accurate transverse shear stress variation through the thickness of such a plate. Furthermore, as expected, the convergence of transverse shear stresses computed by the EM is much slower than that of the nodal displacements [9]. The second objective of the present study is then to devise a method that would force the transverse shear stresses, at special points, to converge almost as rapidly as the nodal displacements themselves. f$!j and f$J, as obtained from eqn (I), are known to be quadratic functions of x and y. The method to be proposed combines the desirable features of quadratic distribution of r’z and f$) over the element interface with the true variation through the thickness, as predicted by the equilibrium equations. Noting that the displacement functions are quadratic in x, y, the convergence of the transverse shear stresses can be expected to be almost as rapid as that of the nodal displacements. Besides, the proposed method, as demanded by the physics of the problem, forces the transverse shear stresses on the top and bottom surfaces of the transversely loaded plate to vanish-a goal that often seems to be so elusive to the EM, when applied to the asymmetrically laminated plates. METHOD OF ANALYSIS The assumption of layerwise constant shear angle implies parabolic variation of transverse shear strains and stresses through the thickness of a layer. This implies, for example, that 111;(z) is of the form

5:; (z) = N,
N,(z)=

3z 2zZ 1 -i_+2, I 1,

through-the-layer-thickness average of r$ (z) (N equations) and (iv) computing jump in r_ at each interface utilizing the first two equations of equilibrium in terms of the stresses (N - 1 equations). The conditions (i) and (ii) above imply lp’ =f$M -_ 0 fy)=f$‘-‘) The condition

(3a)

(4) for

i=2 ,.....,

N.

(5)

(iii) yields

6

1

to I r$j(z)dz

= [O 0 0 0 0 0 c$y C#] [B:“] {d:“}, for

(6)

i=l,...,N

which, with the help of qns

(2)-(S), become

G3s”+ f3f t I)+ t3yt I)) =[000000C~cg][g”]{d~‘}, for (W

Ua)

i=2,...,N-1

+ I,lS)) = [O 0 0 0 0 0 csf) CB] [B$ {dj”}

(7b)

C35”- ‘)+ ;3p> =[oooooocpcp][By]{d~~}.

(7c)

and

The left and right sides of the remaining N - 1 equations are given by the condition (iv) above and are obtained by using eqns (2) and (3) and the first equaton of equilibrium respectively. The jump at the (i + I)* interface, as computed from eqn (2), is

(2)

where x is the transverse coordinate, local to the i”’ layer, and is measured from its bottom surface.fy), fY’ and I$) d eno te r$ (z) at the bottom, middle and top surface of the i”’ layer respectively. N,(z), N2(z) and N,(z) are one-dimensional quadratic shape functions, defined by

629

Jut L

1)=

dr$t",(O)

ar:;(ti) _-

az

= ay3p

az

; afy3$i+ I,+ a;3 ,N+ “I,

[

-0

“y3p

_

; ~y’3p

; ayq

(8)

[

-h

which, with the help of eqn (3), becomes J’ft

1)

I

,(z)=;-7, I

+cfY”‘-&39”1. for r,, distribution through the thickness of a N-layer laminate, therefore, requires 3N unknown parameters, which, in turn, ask for 3N equations. The present method supplies these equations, by (i) forcing I, to vanish on the top and bottom surfaces of the laminate (2 equations), (ii) satisfying continuity of Tu at each interface (N - 1 equations), (iii) identifying <$!j, as computed by eqn (I), as the

i=2,...,N-2

(94 .

Jl” = ;3i1) - 3(i + ;)3$2) + ;3i2’ - ;3i2’ 1

(9b)

JiM = -&3y2)+-5-3y-I) -3

( &+l

1, fS”-” -$J. >

(SC)

630

&xz A. C~~AUDHURI and PAUL SEIDE

Fig. 3. j@ t~anguIar

plate element interface.

The first equation of equilibrium

art; tz)

-

az

=

whereJt+‘),i=ll,..,, N - 1, are given by eqn (9). Equations (4), (5) and (7) and a combination of eqns (9) and (13) determine the 3N unknown parameters required for describing the z,,(z) variation through the thickness of the N-layer arbitra~~y laminated a~isotropi~ plate. The procedure described above applies to both symmetricalty and asymmetrically laminated plates. Following an identical procedure and using the second equation of equilibrium, rYz(z) variation through the plate thickness can be determined.

is given by

-ug, (2) for

r = 1(=x),

2( =y),

(10)

where repeated index ‘r’ implies summation. Integrating over the surface area of the element and applying divergence theorem on the right side of the equation:

Since constant over the will not

5

U(‘)(Z) X, n , dr. (11) rj for quadratic shape functions o&(z) is with respect to x, y, integrating eqn (10) area of the element and dividing by the same alter anything. This will yield =-

Special cases For a homogeneous plate, this method reduces to finding fXXusing eqn (I). ‘t,(z) variation through the thickness being parabolic, (rXi)_ occurs at mid-surface and is equal to $iXz. For a two-layer uns~rnet~~ plate, the eqns (4), (51, (7) and (9) reduce to

1(14b) (144 (14a)

3~)= 3p = - llf2 4(1, f f2)

(1)= ; ?:;I_ 43p 32 3~)= f f!3’- f3S”

(144

where J$) is given by eqn (1) with i = 1. For a three-layer symmetric laminate,

3y =3$3)= ; 7’:; _ f3y’ 391= f fz _1;3p where nf’) for r = I (or x), 2(or y) and k = 1,. , . , 3 and r!‘) for k = 1, . . . ,3 (Fig. 3) for j* triangular elemeit are given by eqns (Al 5) and (Al 6). As the inplane stresses u,(‘)etc. are linear functions of x and y, one Gauss point is sufficient for exact integration on each side of the triangle. Besides, the inplane stresses at the midpoints of the sides are exceptionally accurate. Finally, the jump in

a%* (z) az at the (i + l)* interface may be computed as

3$”=Jp’=t0,

Wa>

(15b)

1

(15c) (15d)

Locarion of exceptional points fir the inieriaminar shear stresses The existence and location of special points, _ _. where the transverse shear stresses are exceptionally accurate can be explained from a physical as well as a rigorous mathematical point of view. The physical explanation for location of the exception points is as follows: the finite element method finds stresses by minimizing the error in the meansquare sense, which implies that equations of equilibrium are satisfied only in the weighted-average (Gale&in) sense. Only at the exceptional points, it is hypothesized, are the equations of equilibrium

Prediction of interhninar exactly satisfied. By virtue of the siphons of LCST and transverse inexten~bility, the variation of displacements through the thickness of each layer has been assumed and the finite element discretixation is restricted to the x-y plane. The consequence of this is that r$ now consists partly of displacement quantity [u(‘)(f,) - G(O)]/r, and, in part, of derivative of displacements w,* - both of which are now independent of z. Noting that quadratic shape functions have been used in the formulation and then writing the first equation of equilibrium as

it can be observed that the variation of the second term of the left side of eqn (I6) is constant with respect to x and y, but quadratic with respect to P while the variation of its first term is constant with respect to z, but quadratic with respect to x and y. Identifying the first term of the above equation to be the through-layer-thickness average, fg and the second term as the average of T:: over the area of the triangufar element interface will lead to the conclusion that the centroid of the element is the point of exceptional accuracy for the transverse shear stresses, when the quadratic shape functions are used for ~mpu~tion of displacements and inplane stresses. The same conclusion can be reached by using rigorous variational argument, known as Nitsche trick, due to Aubin and Nitsche and extended by &hultz, as discussed by Strang and Fix [IO]. They conclude that the average error is much smaller than the typical displacement error at a point. More precisely,

but

for displacement function a’e Sk of degree k’ - 1. This means that the error alternates rapidly in sign; the end practical implication is to even approximately discover the special points where these sign changes occur. Near such special points, the displacements II’ obtained by FJ2M will be of exceptional accuracy. In very simple cases, ‘super convergence* at the nodes has been observed [IO]. In actual ~rnpu~~on like the present work, the ~spla~ments at the nodes are computed and tested for convergences irrespective of their status as exceptional points. Our interest lies essentially in the following: once convergence in displacements has been achieved, find exceptional points which yield accurate transverse shear stresses with minimum further mesh refinement. In a secondorder elliptic problem, the inplane stresses a:), e:) and rfj, which consist of the first derivatives of displacements, are in error at a typical point by

shear stresses

631

0(/i"'I) and on the average

by O@“>. Therefore, these errors must also ahemate in sign and there must exist exceptional points for inplane stresses. Their presence in actual computation was first noticed by Barlow [12] and for quadratic shape functions, midpoints of the sides of the triangular element surfaces are the exceptional points for the inplane stresses. Extending the aforementioned variational argument further, it can be seen that the transverse shear stresses rz and r$ use equilibrium equations in computing jumps at the interfaces of layers and hence involve the second derivatives of displacements with respect to x and y and, therefore, are in error at a typical point by 0(/ry-2) and on the average by O(h”-‘). These errors must also alternate in sign and there must exist exceptional points for the interlaminar shear stresses, which are the loci of centroids of the triangular element interfaces. NUMERICAL EXAMPLES

In order to test the accuracy of the present method as well as the performance of the associated triangular element, the following example probiems will be investigated. Example 1. A square simply supported i.wtroJrrie thin

plate u&r

UniJormnormal pressure

It is essential to test the accuracy of the present method in the case of a homogeneous isotropic thin plate under uniform pressure p. , before applying it to more complex laminations. 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 2 cm. Young’s modulus, E, and Poisson’s ratio, v, of the plate material are 206 GPa and 0.3, respectively. hornets conditions about the plate centerlines and the diagonals permit the modeling to be limited to one-eighth of the plate surface. Antisymmetry about the mid-plane allows us to consider a plate, of half the thickness, under uniform pressure p,/2 so that the inplane displacements vanish at the middle surface (Fig. 4). The present method, when applied to a homogeneous isotropic thin plate, as has been discussed in the preceding section, requires the computation of only throughethickness average transverse shear stresses, 5= and ffl. Since the LCST assumption yields

it might, at first sight, seem reasonable to expect the exceptional points for ru (or z,.J to be the mid-points of the sides, This is because the displacement u(or V) and derivative W&are both exceptionally accurate at these mid-side nodes, when quadratic shape functions are employed [9, lo]. This conclusion has, however, been proved to be fallacious, because the exceptional point, as discussed earlier, requires equilibrium equations to be satisfied there. Furthermore, actual

632

zzurv=O&mLCond)

t12 1:+

Fig. 4. Finite element mesh and displacement boundary conditions of a square isotropic plate.

computations have revealed that the average of r, (or TJ at the three mid-side nodes yields exceptionally accurate results. This leads to the conclusion that the centroid of the triangular element surface is the exceptional point for the transverse shear stresses, which checks with the conclusion reached earlier from the theoretical standpoint. Convergence of (r,), (i.e. 2,(1/2) at x = a/2, y ii: 0), as computed by the EM and the present method, is displayed in Fig. 5. Most rapid converg ence is achieved by a combination of the reduced (i.e. quadratic order) integration and the present method, while the convergence is slowest for the combination of the full (i.e. quintic order) inte~ation and the EM. It is inte~ting to observe that the former combination yields (~xz~~xi~.~. tr,

hnax,exact

=

o 87,62 *

when W

“/F.E.=O.91384 W mu/cud

for a 2 x 2 mesh. w,,,,

and (T~,)~~,_, above are

Fig. 5. Convergence of maximum transverse shear stress for the isotropic plate.

as presented by Plantema [13]. If consideration is given to the fact that the exact theory solution for (r., ), is computed at the point x = a/2, y = 0, while its finite element counterpart has been computed at the centroid of the triangular element adjacent to this point, it can be concluded that T,~, as computed by the present method at the exceptional point, converges almost as’ rapidly as the nodal displacements themselves.

The problem under investigation, which involves severe cross-sectional warping, is a three-layer infinitely long strip of balanced cross-ply (9O”/O”/!W) construction loaded by a sinusoidal normal pressure q = qOsin (ny/a) [Fig. 6(a)]. The width, a, is 24 in,, while the thickness of each layer is 2 in. E,, and Ez2, the inplane Young’s moduli of the layer material in the directions 1 (i.e. fiber direction) and 2 (i.e. normal to the fiber dictions, are 25 x 106 lb/in’ and 1 x 106 lb/in” respectively. The inplane shear modulus, Gr2 and the transverse shear modulus, G,r, are assumed to be equal to 0.5 x IO6 lb@ while the remaining transverse shear modulus, Gr3, is equal to 0.2 x lo6 lb/in’. Major Poisson’s ratio vu = vlJ, defined to be the ratio of strains, perpendicular and parallel to the fiber direction, caused by the stress in the direction of the fibers, is taken equal to 0.25. Transverse Poisson’s ratio, vu, is also assumed to be the same. The problem is important because the accuracy of the present method needs to be established first for a symmetric laminate before applying it to an as~met~~lly laminated plate. Furthermore, it has exact elasticity and CLT (Classical Lamination Theory) solutions obtained by Pagan0 [If. An exact solution using the assumptions of LCST and transverse inextensibility has been given by Seide [4]. Another solution, based on the same theory and using the EM, has been obtained by Chaudhuri [7J. Solutions, utilizing assumed stress hybrid finite element formulation and quadrilateral element shape, have been obtained by Mau et 01. [4] and Spilker et al. [14], using LCST and CST respectively. The finite element analysis is performed utilizing sy~et~ conditions to permit the analysis of a strip of width a/t, which is su~ivided, as shown in Fig. 6(b). An element which fulfills the plane strain condition of independence of displacements and stresses of position along the length of the strip is obtained by combining two right-angled triangular elements as shown in Fig. 6(b) and imposing the conditions of equal displacement perpendicular to the plate (i.e. z-direction) and also in the y-direction at the three nodal points having the same y-coordinate, while all displacements in the x-direction vanish, Fig. 6(b). Antisymmetry about the mid-plane permits analysis of a laminate, of half the thickness, under sinusoidal normal pressure, q/2, so that inplane

Prediction of interlaminar shear stresses

633

9,

I

sol ci g

“-

I-

“;

”

0

aJts4

-

‘--‘---’

’t

c-

.

----

-

~-w--a*.

‘9:

I

I

I

I

00

4

8

12

16

--A Numbar of Subdivisions

Fig. 7. Convergence of interface transverse shear stress at the boundary to the corresponding elasticity solution.

L

I

Fig. 6. (a) Infinitely long three-layer cross-ply strip and (b) its finite element model.

displacement components vanish at the middle surface. First, the convergence of TV at the interface of layers is tested (Fig. 7), which demonstrates that by the present method, con(7yzyr)iumiace 3 computed verges faster than its EM-counterpart. Four divisions with 24 unconstrained degrees of freedom yield an error of approx. 1% while convergence to the elasticity solution is practically achieved (error 0.04%) when six divisions with 36 unconstrained degrees of freedom are used. Figure 8 presents the variation of r,/qO through half the laminate thickness and testifies to the accuracy of the present solution vis-d-vis the elasticity solution of Pagan0 [l]. rvz distribution through half the laminate thickness, as computed by the present method, is practically indistinguishable from their counterparts due to the EM [7] and the analytical solution of Seide [3]. That is why the last solution has not been presented in Fig. 8. As has been discussed also in [7j, the hybrid aproach due to Mau er al. [5], even though based on LCST, displays significant deviation from the elasticity solution of Pagan0 [ 11, because the normal-traction-free-edge condition is not satisfied by the hybrid formulation. Similar deviation also arises in the case of the CST solution due to Spilker et al. [14]. It is, moreover, noteworthy that the CST solution is generally closer to the CLT solution, while its LCST-based counterpart is closer to the elasticity solution.

angle-ply square laminate consisting of two layers of equal thickness, Fig. 9(a). The fiber orientations are - 45” (i.e. clockwise from x-axis) and 45” in the bottom and the top layer, respectively. The top surface of the plate is under uniform pressure, pO. The plate is simply supported on all four edges so that the transverse displacements as well as the mid-surface normal and parallel-edge displacements are not permitted. Normal and parallel-edge displacements are free to oazur at the top and bottom surfaces of the laminate boundaries. Symmetry about a diagonal allows the modeling to be limited to one-half the laminate surface, Fig.9(b). The elastic constants of the orthotropic layer material are the same as those of Example 2 except that E,, is 40 x lo6 lb/in*. Whitney (151 has presented an analytical solution for the problem using CLT. Solutions for displacements

-

Elasticity Theory

-A---

Resent

- -+__*_

Equilibrium Method

-+-

Hybrid Streaa FHte Element (CSn

-.-

uaaakalLamhatbnThecrY

Method

liYbridS~FH@EkllWtcLCS~

Zll

Example 3. Two-layer angle-ply square plate under uniform load A final problem, which establishes the accuracy of the present method as well as the EM of[7j in determination of interlaminar shear stress distribution through the thickness of an asymmetrically laminated plate, is that of a balanced unsymmetric

1121 0

0.4

I

I

0.8

1.2

1.6

I 2.0

TYtJq.

Fig. 8. Comparison of transverse shear stress variation throuab thickness at the boundarv.

REAZ A. CHAUDHURIand PAUL SUDE

634 Y

2

‘Iho v-5

‘Jr

--x--

present

---A-

Equilibrium

-

CLT

Method Method

X

Fig. 9. Geometry, finite element model and boundary conditions for a square two-layer balanced unsymmetric angle-ply plate.

and inplane

stresses

have been presented

T,,tha Fig. IO. Transverse shear stress variation through thickness of the unsymmetric angle-ply plate at x = 35a/36, y = 19a/36.

by Spilker

et nl. [6, 141using an assumed stress hybrid formulation, based on CST and quadrilateral element shape. These are in close agreement with the corrected versions of Whitney’s [ 151solutions. The interlaminar shear stress distribution through the thickness of this unsymmetric laminate has not, however, been presented by Spilker er al. [6, 141.This is a bit surprising, because the interlaminar shear stresses are the most important quantities to be investigated because of the critical role they play in causing delamination failure in layered composites. Figure 10 displays the results for interlaminar shear stress distribution through the laminate thickness at the point x = 35a/36, y = 19a/36. Twelve divisions on each side, with 288 layer elements, have been employed in obtaining these results. Moreover, reduced (quadratic order) numerical integration scheme has been used in the computation of the element stiffness matrices and consistent load vectors. The present solution is compared with that due to the EM [7j as well as the CLT solution. As expected, the EM solution for the transverse shear stress does not vanish on the top surface, while the present solution is forced to do so. The difference between the present solution and the CLT may be attributed largely to

the effect of transverse neglected

shear deformation which is in the CLT. The error in interlaminar shear

stress at the top surface (where it should vanish), computed by the EM at the location, x/a = 35136, u/a = 19/36, appears not too large and may even be considered acceptable. However, the accuracy of the EM, as applied to an unsymmetric angle-ply plate, caMOt he guaranteed as will be evident from Table 1, which exhibits three extremely poor situations wherein the EM error (i.e. the computed interlaminar shear stress at the top surface) is comparable or even larger than the maximum interlaminar shear stress at the location considered. CONCLUSIONS An approximate semi-analytical method for prediction of interlaminar shear stress distribution through the thickness of an arbitrarily laminated anisotropic plate has been presented. Even though the method is illustrated here for an assumed quadratic displacement triangular element, the principle behind the method is general enough for any element shape. Furthermore, this paper puts forward

Table 1. Comparison of normalized interlaminar shear stresses at certain locations computed by the EM and the nresent method Location y/a

Method

13/18

19136

25126

17136

29136

11/18

EM Present EM Present EM Present

x*=x/a

y*=

Normalized interlaminar shear stress 7: = 7-r /qa 7$‘(r3 7$(f2/2) 7$)(O) 7$(1,/2) 73r,) = 730) - 1.155 - 1.49 - 1.378 0 0.194 0 0.054 -0.519 -0.05 0 1.838 0 -0.699 0.402 1.607 0 -0.131 -0.080 -0.142 0 - 0.745 0 -0.051 -0.096 -0.700 0 -0.179 - 0.025 -0.279 0

Prediction of ~~r~~ar

the physical interpretation that the exceptional points for the interlaminar shear stresses are the points where the equations of equilibrium are nearly satisfied. The centroid of a triangular element surface has been proved, using the rigorous variational argument of Aubin and Nitsche, to be the point of exceptional accuracy for the interlaminar shear stresses. This paper essentially demonstrates that it is Rossible to obtain an accurate solution, based on LCST,

utilizing the ~nventional assumed displacement potential energy aproach (FEM) and that it is not necessary to resort to the more cumbersome assumed stress hybrid FEM. In fact, the superiority of the present method over the hybrid approach has been clearly demonstrated for both the symmetric and unsymmetric laminates. Deviation of the transverse shear stresses for a three-layer balanced symmetric cross-ply computed by the hybrid method, from that due to the exact elasticity solution, has been attributed to the fact that the normal-traction-freeedge condition along the simply supported edge has been, at best, app~~mately satisfled by the hybrid method. The close agreement between the present solution with that due to the exact elasticity theory suggests that this condition is automatically satisfied in the present case. Besides, the assumed stress hybrid method has not presented the interlaminar shear stress distribution through the thickness of the twolayer balanced unsymmetric angle-ply laminate, even though the displacements and the inplane stresses have been presented. The present method is the first

to present a reasonably accurate solution for the interlaminar shear stress distribution through the thickness of the afo~mention~ laminate. The present method has been demonstrated to be superior to the Equilibrium Method (EM) off;rl in the case of an unsymmetric laminate, because the former assures faster convergence as well as the vanishing of the transverse shear stresses on both the exposed surfaces of the laminate. The superiority of the present method over thd EM in the case of a

symmetric laminate is limited to faster convergence alone. It has also been demonstrated that, for a thin plate, the combination of the present method and the reduced integration scheme yields convergence of the transverse shear stresses almost as rapidly as that of the nodal displa~en~= It has further been concluded that, for a thick laminate, the LCST-based solution is closer to its counterpart due to the threedimensional elasticity theory, while the corresponding CST solution is closer to that due to the CLT.

shcu strews

635 -cl?s

1. N. J. Pagana,Exactsolutionfor composite laminates in cylindrical bending. /. camp. Morer. 3,39&-411 (1969). 2. N. J. Pagano, Exact solution for rectangular bidirectional composites and sandwich plates. 1. camp. Mater. 4, 2&25 (1970). 3. P. Seide,An improved approximate theory for the bending of laminated plates. In Mechanics Toalay (Edited by S. Nemat-Nasser), Vol. 5, pp. 451-465. Pergamon Press, New York (19gO). 4. R. M. Barker, J. R. Dana and C. W. Pryor, Stress ~n~~tions near holes in laminates. J. Engng Me&. Iliv., Am. SC. civ. Engrs. X00,(EM3), 477488 (1974). 5. S. T, Mau, P. Tong and T. H. H. Pian, Finite element solution for laminated thick plates. J. camp. Mater. 6, 304-311 (1972). 6. R. L. Spilker, 0. Orringer, E. A. Winner, S. Verbicse, S. E. French and A. Harris, Use of hybrid-stress finite element model for static and dynamic analysis of muhilayer composite plates and shells. Dept. AMMRC CTR 76:29, ASRL TR 181-2, M.I.T., Cambridge (1976). 7. R. A. Chaudhuri,An equilibriummethodfor prediction of transverseshear stressesin a thick laminated plate. Comrput.Struck. u, 139-146 (1986). 8. C. W. Pryor, Jr and R. M. Barker, A finite element analysis incIu&nS transverse shear etfects for applications to laminated plates. AMA Jtd 9,912-917 (1971). 9. R. A. Chaudhuri, Static anaiysis of fiber reinforced laminated p&s and shells with shear deformation using quadratic triangular elements. Ph.D. Dissertation, USC, Los Angeles (1983). 10. G. Strang and G. J. Fix, An Analysis of the Finite Element Method. Prentice-Hall, Englewood Cliffs, NJ (1973). 11. R. A. Chaudhuri and P. Seide, Triangular finite element for analysis of thick laminated plates. Int. 1. nwner, Meth. Engng (in press). 12. J. Barlow. Optimal stress isolations in finite element models. ht. 1. namer. Meth. Engng X$243-251 (1976). 13. F. J. Plan&ma, Sandwich Construction. John Wiley, New York (i966). t4, R. L. Spi]ker$S. C. Chou and 0. &ringer, Alternate hybrid-stress efrments for analysis of multi-layer composite plates, J. camp. Mater. 11, 51-70 (1977). 15. J. M. Whitney, Rending-extensional coupling in laminated plates under transverse loading. J. camp. Mater. 3, 20-28 (1969). APPENDIX This section presents some of the matrices and equations referred to in the section on background information and theoretical formulation of this paper. Stress-strain relations for an anisotropic Laminaare given by {@‘Of = ICC’91 (f(Q) (AlI where in the context of FEM formulation [ 12-131

while ~2 and 0 are inplane normal strains, y$i inphtne shear strains am&! and v’s are (through-thickness) average transverse shear strains in the i* layer, and

Acknowledgements-The research reported in this paper

was supported under a grant (NSG 1401)from NASA tangicy Research Center, Hampton, Virginia, U.S.A. it is a pleasure to thank Dr Manuel Stem, the NASA Technical Oftker monito~~g the grant. C.&s. 25/4-l.

[C”] =

IA31 ~.

Rwz

636

A.

CHAUDHURIand PAUL SEIDE

For a lamina with fiber orientation O,, [C”)] can be obtained from the orthotropic material matrix [C(‘)] by means of the relation pj [C”‘] = [T]T[C(“] [T]

[&I =

(A4)

while

(AlOb)

F-1= [VII.. . [Tt.ll

(Al la)

with cos2 Bi sin2 6,

fsin2 6,

0

sin'e, ax2 8, +inze,

- sin2 8, sin2 8,

VI= I

~0~2,~

0

-

0

0

0

0.

(Allb)

0

0

0

cos 8, sin 8,

0

0

0

-sin 8, cos ei

fork=1,...,6, [N’O] = [[NI”]]

. [Np]]

(AW

W)

with six non-vanishing elements of [Cc”] can be expressed in terms of the generalized Young’s moduli, Poisson’s ratios and shear moduli of the layer material [7]. [AQ] referred in eqn (I) is given by

(A12b)

and [kj”‘] = - [NW]

where

1-f [Ab”] =

0

0

’ l-f

0

0

for i=l,...,

0

0 I

(AW

1-f

(A13)

Nandk=l,...,6.

$k represents the shape function for the k' node of a triangular element interface. For the quadratic shape function chosen here, {4} is given by MI’=

ICl(2CI-

u.4Li2.

c2012-

1).

foe Y2L.

[AI’)]=

O;o

(AW

00;

WC,

-

11, U,

1

(At4)

where (,, c2, C, are the area coordinates of a point inside the triangular element interface. For q!“) referred to in eqn (5) and (Fig. 3): (AlSa)

and [A (‘J] = I

1 0 [ 0

I

I

(Al5b)

(A7c)

;

{d!)}, as referred to in eqn (I), is

(A154 (AlSd)

. . . P~+‘),d~+‘), w,, . . . , we} (A8) (

(AW

where the subscripts refer to the nodes of thejlh element. [By)], as referred to in eqn (l), is

(A15f) ]gj’)] = [ ];i]

,i]

i]

(A9)

ry =(x:

whose individual submatrices are given as

WI = H41.. . [&II with

where

(AlOa)

+y:p

w =[(x2-

x,)2+ (y2-YEW

ry=(~:+y:y'~.

(AlW (A16b) (A16c)