Finite element analysis of progressive failure in laminated composite plates

Compurers & S~luclures Vol. 38, No. 3, pp. 361-376, Printed in Great Britain.

1991 0

o!Ms-7949/91 s3.00 + 0.00 1991 Perpnon Press pk

FINITE ELEMENT ANALYSIS OF PROGRESSIVE FAILURE IN LAMINATED COMPOSITE PLATES S. TOLS~N and N. ZASARAS Department of Mechanical Engineering, University of Minnesota, 111Church Street SE., Minneapolis, MN 55455, U.S.A. (Received 8 February 1990)

Abstract-Acceptance and utilization of composite materials require confidence in their load carrying capacity. Therefore, it is desirable to develop a computational model capable of determining the ultimate strength of laminated composite plates under conditions of complex loading. A new seven degree of freedom finite element model for laminated composite plates is developed. The model utilizes three ~sp~~ments, two rotations of normals about the plate midplane, and two warps of the normals, to accurately and efficiently determine the laminate- stresses, Based on these stresses a failure model for determining first ply failure (FPF) and last ply failure (LPF) by a progressive stiffness reduction technique has been developed. The progressive failure model produces results in good agreement with experimental data. The calculated FPF and LPF results form lower and upper bounds within which the true load carrying capacity lies.

INTRODUCTION

The objective of this paper is to construct a methodology for determining when failure of a composite plate occurs and to develop the mathematical and computations tools to carry out the analysis. Research into the analysis of plates has been extensive. Kirchhoff was one of the first to develop a comprehensive theory for plates. His classical plate theory yields inaccurate results for thick plates (S = a/h < 20) and for thin plates with holes. Reissner was one of the first to develop a plate theory that considered transverse shear deformations in static analysis [ 1,2]. Mindlin later expanded the shear deformation theory to accommodate rotary inertia terms [3]. These first order shear deformation theories fail to produce accurate results for very thick plates (S c 10) and they do not satisfy the condition of vanishing transverse shear stress on the top and bottom plate surfaces. They also require the use of a shear correction factor to obtain accurate stress results. Higher order plate theories seem to begin with Levinson [4]. The displacement field for his theory includes up to cubic variation in the throughthickness direction. The theory predicts zero transverse shear stresses on the top and bottom plate surfaces, does not require a shear correction factor, and is accurate for both thin and thick piates. The preceding theories were all developed for isotropic materials, but each one has been generalized for anisotropic materials. Reissner and Stavsky developed an orthotropic laminated plate theory based on classical plate assumptions [S]. A first order shear deformation theory for laminated anisotropic plates was presented by Yang ef al. [6]. Finally, higher order

shear deformation theories for laminated composite plates were developed concurrently by Reddy [A and Murthy [8]. Reddy’s penalty plate bending fmite element technique is one of the first finite element models for first order shear deformation theory [9]. The theory uses eight-noded quadrilateral elements with five degrees of freedom at each node. The degrees of freedom include three displacements and two rotations. Higher order theories contain second order derivatives of the transverse displacement, Therefore, a finite element model requires C’ continuity. Construction of such an element requires many degrees of freedom at each node and, therefore, much computer processing time. In an attempt to efficiently solve this problem, Reddy developed a mixed formulation for his simple higher order theory [IO]. The resulting finite element model consists of 11 degrees of freedom per node, three displacements, two rotations, and six moment resultants. While this may have been an improvement, it remains computationally intensive and a more efficient alternative is required. Failure of a composite material can be described in a variety of ways. First ply failure (FPF) occurs when initial failure of a single layer in a laminate fails in either the fiber direction or in the direction perpendicular to the fibers. Last ply failure (LPF) occurs after the structure has degraded to the point where it is no longer capable of carrying additional load. Most authors beheve failure is caused by (1) lon~tudinal tensile loads in the fiber direction, (2) longitudinal compressive loads in the fiber direction, (3) tensile loads transverse to the fibers, (4) compressive loads transverse to the fibers, or (5) shear loads [l 1,121. Utilizing these five loading cases, the failure modes of 361

362

S. TOLSON and N. ZABARAS

a composite material can be described as (1) breaking of fibers, (2) cracking of the matrix, (3) separation of the fiber and the matrix (debonding), and (4) separation of one iamina from another (delamination) (131. The consequence of individual failure modes is not of great interest for FPF, but will be of great importance in determining LPF. Standard laminate strength analysis is the most common and oldest of the composite analysis methods [II, 121. The method neglects local effects such as fiber misalignment, material discontinuities, and free edge effects, and assumes that the stiffness of the laminate receives no contribution from failed layers. The first finite element based failure analysis of composite materials was performed by Lee [14]. Lee used his own direct mode dete~ning failure criterion and standard laminate strength analysis methods to determine the ultimate strength of plates with circular holes. The major drawback of a threedimensional failure analysis such as Lee’s is the tremendous amount of memory space and calculation time required. Hwang and Sun attempted to improve the computational aspects of the three-dimensional formulation by incorporating a Newton-Raphson type iteration process [ 151. The search for more efficient finite element analysis of composite plates, therefore, leads to twodimensional plate fo~ulations. Reddy and Pandey developed a first ply failure analysis of composite laminates based on first order shear deformation plate theory [16]. The limiting factor of this analysis is the inadequacy of first order shear deformation theory for thick composite plates. Engblom and Ochoa developed a two-dimensional plate analysis similar to the above, but with increased interpolation in the through thickness direction [17, 181. Their analysis is carried out to LPF. The stiffness reduction and progressive damage accumulation are treated in a manner similar to standard laminate analysis. The generated stresses are less accurate than those obtained from higher order shear deformation plate theory. Also, since the interpolation is of a different order through the thickness than in the plane, the analysis cannot be easily adapted to elements other than the eightnoded variety discussed in their paper. The objective of this paper is to develop a twodimensional finite element failure analysis for composite plates that is more accurate and more flexible than previously developed plate analyses, while at the same time more efficient than current higher order shear deformation or three-dimensional formulations. The plan of this paper is as follows. First, a new seven degree of freedom finite element composite plate formulation for laminate stress calculation will be presented. Then a brief review of some failure theories will be given. An algorithm for analyzing progressive failure in laminated plates will be given and tested on several problems with various stacking

sequences and load conditions. experiments will be reported.

Comparison

with

A SEVEN DEGREE OF FIRM FINITE ELBOW MODEL FOR LAMINATED COMPOSITF. PLATES

Kinematic assumptions The following displacement field is assumed u(x,Y,zf=u~x,Y)-zY~Y,(x,Y)+2*5,(X,Yt + Z3t;,kY) V(x,Y,z)=u(x,Y)-zY:(x,Y)+z2r,(X,Y) + Z35&*Y) W&

Y) = we,

YX

(1)

where u, v, and w denote the displacements at the midplane (z = 0), and Y, r, and [ are appropriately selected functions of x and y. The above displacement field can be simplified by utilizing the condition that the transverse shear stresses a,, and cry2vanish on the top and bottom surfaces of the plate. The final expressions for the displacement field in an orthotropic plate are as follows: u=U(X,Y)-z~~Y,(x,Y)+Z3t;,tX,Yf ~=~(x,Y)-zy:(x,Y)+~3ry(x,Y) w=

WCG Y),

(2)

where YXand Y,, are rotations of the normals to the midplane about the x and y axes respectively, and 6, and {, describe the warping of the normal in the x and y directions, respectively. It can be seen that for this displacement field, normals to the midplane of the plate before deformation do not necessarily remain normal or straight after deformation. The strains can now be derived using the final form of the ~splacement field, in eqn (2).

au

ay,

au

ay if&

C,=z$-Z-g+Z

4’&i-z

ay

3

ai,

z

ay

-G=o

y,,= -YX+;

+3.r21,

_a,+“_z(!T$+!s)+zf!!+~). yxy ay ax (3)

363

Progressive failure in laminated composite plates

The reduced stiffness matrix, Q’, from the material coordinates to the global coordinates, can be calculated as

Note that the assumption 6, = 0 is retained from first order shear deformation theory. Constitutive model

(6)

If 6, = 0 and a transversely isotropic material is assumed, the stress-strain relationship can be stated as[ll]

[Q’l= P’IIQW-‘I.

where TV,= o,, 02 = (rm, a, = au, tr5 = Cm, be = 6L* (L and T are the fiber and normal to fiber inplane directions, respectively), and the reduced stiffness components, Q,, are as follows:

Therefore, the stress-strain relationship in the global, xyz coordinate system can be stated as the following:

(;j=(;;

i;

Z)($

Q,,=E’ 1 -

Q,2-

h2V21

v12E2 1 -

Q22=E’ 1Q.M =

_

Vl2V21

vzlE1 1 -

h2V2I

h2V2I

G23

Qs = G2

=

Equations (7a) and (7b) are considered as the fundamental constitutive equations in terms of the global coordinate system.

Cab)

Gs.

The relationships used to transform stress and strain from one coordinate system to another are as follows:

P-l{~j

(54

16’1= [Tl+L

WI

(6’) =

Finite element model

Isoparametric finite elements are introduced, the transformation of space is as follows:

and

and

y= i

where 0’ and 6’ are the stress and strain components in the xyz coordinate system and

i=l

r

PI =

7

2sinO

cor.0

cos20

sin’ 0

0

0

sin2 0

co? 0

0

0

0

0

cos 0

-sin

0

0

sin 0

cos 0

0

0

0

cos2 0 - sin2 8

L -cosOsinO

cos@sinQ

-2sinBcosO 8

0

,

y;w(r,rl),

@b)

S. TOLK~N and N. ZABAIW

364

where

where XT and Y; are the global coordinates of the element nodes, NT are the shape functions and n is the number of nodes in element e. The displacement field for each element can be written in terms of the shape functions as follows: u(x,Y) = t U,N:,

6: is E~, 6: is f, and ej’ is yXy c4 is the extensional strain component cfi is the Hexural strain component cWiis the warping strain component.

(9)

i=l

and C: = es, + e,,iz*

where

(i = 4, 5),

(11)

where 6: 1s yyr, and 6: is yXz cs, is the first order shear strain component &swIis the shear-warp strain component

is the displacement vector and

The strain-displacement relations for higher order theory can be determined for each of the strain components in eqns (10) and (11). Based on the assumption of small displacements, one can write for each element e the following:

are the nodal displacements. The eight-noded serendipity and nine-noded heterosis elements are used in this analysis. The strains for higher order plate theory can now be derived in terms of the element nodal displacements. First, it is helpful to subdivide the total strains from (3) into components associated with extension, flexure, warp, shear, and shear-warp, as follows: &: = Cri+ Cf,Z+ t,,z’

(i=1,2,3),

(10)

where

aN, 0

0

-a&

0

aNI --

aN, ax

ax [BcI@)=

0

[&I”’ =

0

0

0

0

0

0

0000

0

0

0

,

0

0

0

0

0

0

0

..*

0

0

0

0

0

0

o...

ooo-2

0

-$

0

0

‘..

000

-z

-2

0

0

... 0

0

0

2

0

0

0

0

‘..

-...o aN,

ay

aN, aN,

--~..oooo ay

ax

0

0 0 0

0

0

0000

[B,,.]“’ =

0

-$

000

0

0

fs ay ooooo...~~ooooo ,

ay

ay

..$

. o

0

0

0

0

0

0

0

0

0

0

0

_alv,

0

0

-$

0

0

ay

-2

0

aN

2

0

aN, 0

00

oaN”

ay dN,

Fax

J

I

00%-N, 00% ay

[B J I@)=

3

003-N ay

-N,OO...OO 0 OO...

0

ax [B SW l(e)=

365

Progressive failure in laminated composite plates

00000 000003N,

0

3N,...OOO 0 . ..OOO

Higher than linear order derivatives of the displacements do not exist in the strain field. This means that only C,, element continuity is required. The principle of virtual work can now be used to derive the stiffness and load matrices. Assuming that the only external force acting on the system is the distributed load, q, the following equation can be obtained:

-N, 0

0

00 0 0

n 3N,

0

0

0

0

0

3N,

0 >’

(Sl,,

s2,,

S3,) =

s hi*

L,(l, z*, z4)dz

(i,j = 1,2)

- hi2

(16b)

and NT the set of shape functions for the element being used. All terms involved in eqn (14) have been calculated using a reduced 2 x 2 integration scheme. The frontal solution method was used to solve for the unknown nodal displacements without explicitly assembling the stiffness and load matrices. Interlaminar stress calculation

+ z46~~[‘l]c,,)=

SuTq dx dy,

(13)

sr

where V is the material volume occupied by the plate and r is the top surface. Substituting in the straindisplacement relations [eqn (3)] and integrating with respect to the z coordinate the above expression finally becomes Ku = F,

The inplane stresses are calculated from the constitutive equations, eqns (7a), but the transverse stresses are evaluated from the three-dimensional equilibrium equations. This method gives results very similar to the three-dimensional elasticity solution. The transverse shear stresses are very accurate, and they are continuous through the thickness of the laminate. Also, this method provides a way to calculate the transverse normal stress. Substitution of the constitutive equations into the stress equilibrium equations in the x and y axes, and integration in the z axis, finally gives the following expressions: Aa,z= -

where,

a% Q,,g+Q,,ax ay

+~B,IT[BI[Be1+ P$IT[~IP/l+ Pfl~FIP,.l +PL1T[~l[Be1+PwIVlP,l+

PwlT[~l~B,,J

+~~~1~~~1~~~1+~~~1~~~1~~,,1+~~,,.1~~~1~~,1 + Psw17s 31[BwI) d-x dy

(14)

and F=

NTqdxdy,

(15)

s Iwith (A,, Bij, D,, Eij, Fi,, Cii) hi* = Qij(l,z,z2,z3,z4,z6)dz s -h/2

(i,j=l,2,3) WM

(174

366

S. TOLS~N and N. ZABARAS

Aq,z= -

ah Q,$+Q,2axay

where

.,,,($+2&)+Q**$

+P33(&+~)]z

P is the value of the pressure load on the top surface, and aZj is the transverse normal stress acting at the interface of the (j - 1)th and jth layers. Thus, a continuous transverse normal stress distribution can be found at each Gauss point in the laminated plate, by solving the above equations set using a least squares procedure.

+;[ld$

a*y -+a,3+%(~+3) ax ay

+QQ,*

Numerical examples

The current finite element formulation will be compared with classical plate theory, three-dimensional elasticity solutions, and other finite element formulations. The orthotropic material properties used for the comparison are those of a graphite/epoxy compound and will be specified as $+2g31z4c"', X

UW 2,

where zi and zi+ , are the z coordinates at the bottom and top of the ith lamina, respectively. Since the transverse shear stresses are zero on both the top and bottom surfaces of the laminate, it is a simple matter to sum the shear stress contributions from each lamina to obtain the transverse shear stress distribution through the thickness of the plate [19]. The calculations indicated in eqns (17) are performed at the integration points of each element and then transferred to the four corner nodal points using a local least squares technique [19]. The values of the derivatives of the transverse shear stresses can now be determined at the Gauss points. It is then possible to write the through-thickness variation of the transverse normal stress in numerical form at a single Gauss point as

E, /E2 = 25.0 G,*/G2, = 2.5 v,* = v2, = 0.25.

All plates analyzed will be square with planar dimensions, a x a, and total thickness, h. The nine-noded heterosis element is used in all the calculations reported in this section. The stresses from the current theory are evaluated at the Gauss points, then extrapolated to the nodes for direct comparison with the elasticity solutions. The stresses reported in the following table and figures will be stated in their normalized form. This is done to remove the effects of varying applied loads and changing aspect ratios, S = a/h. The normalizing equations used are as follows:

o* =-

x

where k denotes the layer or portion of a layer under consideration. For a laminate of n layers, n similar equations can be written in terms of both the unknown normal stresses at the layer interfaces and the normal stresses at the laminate surfaces. Assuming the bottom surface of the laminate to be stress free and the top surface to have a normal stress equal to the pressure load, n equations in n - 1 unknowns can be written. The equation set is over determined and has the form specified below:

-1

1 -1

1

-1

= I

)+I*=---

Qx

PO&s

*=Y

a

ay PoS

wuh’ Poa4 ’

where

n4 ‘tG,* + [E, + (1 + ~v,~)E~]

a=E

(1 - v12v21)

)

and P,, is the maximum value of the pressure load applied on the top surface.

367

Progressive failure in laminated composite plates Table 1. Normalized stresses and displacements for a simply supported O/90/90/0 square plate Approach

S

w!!O 2’2’

aah @.t 2’2’* --

a a h --=’ 2’2’4

o,, O,$O

a 00 ovz j’ 9

0.665 0.720

0.701 0.678 0.615 0.632 0.666

0.0495 0.0472 0.0341 0.0440 0.0467

0.219 0.237 0.300 0.206 0.270

0.280 0.291 0.336 0.239 0.292

1.693 1.671 1.790 1.534 1.643 1.709

0.605 0.575 0.490 0.484 0.546 0.559

0.429 0.401 0.382 0.350 0.389 0.403

0.0299 0.0280 0.0257 0.0234 0.0268 0.0276

0.291 0.299 0.325 0.264 0.301

0.184 0.191 0.209 0.153 0.196

20

1.188 1.177 1.216 1.136 1.163 1.189

0.586 0.557 0.528 0.511 0.539 0.543

0.324 0.313 0.299 0.287 0.304 0.309

0.0249 0.0235 0.0228 0.0214 0.0228 0.0230

0.313 0.326 0.339 0.283 0.328

0.148 0.153 0.167 0.123 0.156

Present (2 x 2) Present (4 x 4) Ref. f17] (6 x 6) Ref. [9] (2 x 2) Ref. [7] Elasticity [2 l]

50

1.036 1.026 I.045 1.019 1.031

0.581 0.543 0.542 0.520 0539

0.296 0.282 0.264 0.265 0.276

0.0235 0.0220 0.0217 0.0207 0.0216

0.318 0.334 0.345 0.337

0.135 0.139 0.152 0.141

Present (2 x 2) Present (4 x 4) Ref. [17] (6 x 6) Ref. [9] (2 x 2) Ref. [7] Elasticity [21]

LOO

1.013 1.006 1.021 1.005 I.001 1.008

0.580 0.542 0.544 0.523 0.539 0.539

0.286 0.277 0.259 0.263 0.271 0.271

0.0233 0.0218 0.0216 0.0207 0.0213 0.0214

0.322 0.335 0.346 0.290 0.339

0.133 0.137 0.146 0.112 0.139

1BOO

0.539

0.269

0.0213

0.339

0.138

Present (2 x 2) Present (4 x 4) Ref. [17] (6 x 6) Ref. [9] (2 x 2) Ref. [7] Elasticity [21]

4

Present (2 x 2) Present (4 x 4) Ref. [17] (6 x 6) Ref. [9] (2 x 2) Ref. [7] Elasticity [21]

4.498 4.393 5.245 4.355 4.491

0.763 0.723 0.387

10

Present (2 x 2) Present (4 x 4) Ref. [t7] (6 x 6) Ref. [9] (2 x 2) Ref. [7] Elasticity [2 l]

Classical plate

The example considered is a simply supported (0~90~90/0) symmetric cross-ply laminated plate subjected to a transverse sinusoidally distributed pressure of the following form: P = POsin m/a sin ny/a.

(20)

The origin of the plate is located at the lower left comer on the midplane. The boundary conditions used are as follows [20]: 0(X, a/2) = Y,(X, a/2) =&(x, a/2) = 0 uWZY)

= ~.&/ZY)

V(O,Y) = MJ,Y)=

=L@/~,Y>

=O

yly(O,Y) = 5,(&Y) =o

u(x, 0) = w (x, 0) = YX(X, 0) = rX(x, 0) = 0.

The normalized stresses and deflections for the simply supported O/90/90/0 plate are compared to elasticity and other finite element solutions for aspect ratios ranging from four to 100 in Table 1. Convergence is demonstrated by providing two levels of mesh refinement. It is seen that a 4 x 4 mesh of elements provides better results than a 2 x 2 mesh when compared to the exact elasticity solution of Pagan0 [al]. The re-

suits of the present theory are at least as accurate, in comparison with the elasticity solution, as any of the other theories for which results could be obtained. Results from all analyses are equally accurate at high aspect ratios, that is for thin plates. The results from Reddy’s penalty plate bending analysis [9] at low aspect ratios are the worst of those presented in Table 1 as compared to the elasticity soiution. The inaccuracy is to be expected since the penalty plate bending analysis is a first order shear deformation analysis and does not account for warping of the plate. The analysis of Engblom and Ochoa provides somewhat better results, but even with a fine mesh of elements, (6 x 6), the stresses are still not tremendously accurate at low aspect ratios 1171.Their analysis is essentially a first order shear deformation analysis with increased interpolation in the throughthickness dimension. The analytical results from Reddy’s simple higher order theory [7] compare quite closely with the present analysis and the elasticity solutionf21j. The present theory is especially accurate in calculating the transverse shear stresses due to its direct determination of the warp terms. Finally, Fig. 1 shows the calculated transverse normal stress. No comparison is available.

S. TOWN

368

and N.

ZABARAS

Fig. 1. Transverse normal stress distribution for a O/90/90/0simply supported plate. In conclusion, it can safely be said that the present theory provides the most accurate results of those compared for a O/90/90/0 laminated plate. More examples and comparisons are given in [19]. PROGRESSIVE

FAILURE

ANALYSIS

Review of failure criteria

Failure criteria fall into four basic categories: (1) limit theories, (2) polynomial theories, (3) strain energy theories, and (4) direct mode determining theories. Several papers can be found which list the most commonly used composite failure theories [22-241. The limit theories compare the value of each stress or strain component to a corresponding ultimate value. The polynomial theories use a polynomial in stress to describe the failure surface. The strain energy theories attempt to use a nonlinear energy based criteria to define failure. Finally, the direct mode determining theories are usually polynomials in stress and use separate equations to describe each mode of failure. The maximum stress failure theory is the dominant member of the limit failure theory category. The maximum stress failure criterion for a transversely isotropic material can be expressed as follows:

;=I

04 -=

ST

or

$=I

s*

F,ai+Fijaiaj=

1.

(22)

Hoffman proposed a fracture criterion for brittle orthotropic materials based on the work of von Mises and Hill [26]. For a transversely isotropic material the Hoffman criterion has the following form [27]:

+-

ala2

xx-

+

0 2

‘+wzl,

s,

(23)

A

The above equation includes both the tensile and compressive strengths of the material. Finally, let us examine the Tsai-Wu failure theory [25]. Because of its general nature, this theory contains almost all other polynomial theories as special cases. For the case of a transversely isotropic material the criterion is as follows:

I

05 -=1 S” a6 -=

where the stresses a,-u6 have been defined earlier. In this and all following expressions, X and X-, are the tensile and compressive strength in the fiber direction, Y and Y-, are the tensile and compressive strength in the direction transverse to the fibers, S, is the axial shear strength associated with the l-2 and l-3 planes, and ST is the transverse shear strength associated with the 2-3 plane. The vast majority of the polynomial theories are quadratic and have the form [25]

F,a, +F2a2+F,,a:+2F,,a,a2

1,

(21)

+F22a:+F,a:+F,(a:+a:)=l,

(24)

369

Progressive failure in laminated composite plates Tensile fiber mode:

where

F,=_:+& F,,

=-

F2=_:+$ F,=-

-1

or

0f 2@:+u:)=

XX-

61~1 x

-1

Compressive fiber mode:

1.

(26a)

A

YY01 --=X-

(26b)

1.

Tensile matrix mode: The value of the interaction coefficient, Flz, is determined from biaxial test data. Wu determined that FL2 can be approximated in the worst case as &,/(F,,F,,). The strain energy failure criteria are primarily used with nonlinear theory and are of no concern in this paper. The direct mode determining failure criteria are possibly the most useful of the four categories. Two theories proposed by Hashin and Lee provide separate failure equations for each mode of failure [12, 141. The mode determining ability of the limit theories is of great importance in the cumulative failure stage of composite failure because it allows for easy stiffness reduction. Hashin developed his failure criterion by considering the fundamental stress invariants for a transversely isotropic material. In a three-dimensional stress space he determined that four distinct quadratic polynomials describe the four modes of failure. They are as follows: Tensile fiber mode:

0 %

2

or

++(o:+o:)=l

A

0,=X.

(25a)

Compressive fiber mode: 6, = -x-.

(25b)

Tensile matrix mode:

0

_: 2(u2+uJ+&+~2c))+~(cr:+c+1. T

A

cw

Compressive matrix mode:

a2

-=l Y

(25d)

Lee also proposed a direct mode determining failure criterion [14]. Lee’s criterion, unlike Hashin’s, is based entirely on empirical reasoning. Lee’s failure criterion in terms of the four failure modes is as follows:

2(u:+u:)=

1.

(26~)

Compressive matrix mode: --=

02 y-

1.

(26d)

Algorithm for incremental analysis of failure

As stated in the introduction, the procedure for determining the strength of a laminate involves an incremental load analysis. For a given load, the stresses in each lamina can be calculated with respect to the material coordinates. These stresses are inserted into the appropriate failure criterion to determine if failure has occurred within a lamina. When failure occurs, the stiffness is modified and, upon satisfaction of equilibrium, the load is increased until final failure is reached. After failure at one point, the load continues to be carried by the remaining fibers and matrix of the lamina. Failure of a portion of one lamina is compensated for by an increase in the load carried by adjacent laminae. The analysis progresses in the following manner. The appropriate failure criterion is chosen from those described previously. The stresses are transformed into the material coordinate system using eqn (7), and substituted into the failure criterion to determine if failure occurs in any lamina of any element at the initial load. If no failure occurs at the initial load, then the FPF load is calculated. Since the analysis is elastic until failure, it is possible to determine the failure load by simply scaling up the stresses until the failure ratio (the value found by evaluating the failure criterion) is equal to 1. The initial load is multiplied by this factor to give the first ply failure load. If failure occurs at the initial load, the analysis can be restarted at a lower initial load. After first ply failure has occurred, it is necessary to determine the mode of failure so the stiffness can be reduced in the correct manner. The maximum stress, Lee, and Hashin failure criteria automatically determine the mode of failure. However, the Hoffman and Tsai-Wu failure criteria do not. For the latter, it is necessary to determine the relative sizes of the contributions from shear terms, transverse (to the fiber) direction terms, and fiber direction terms.

&((3- +~2+.3)+&u2+u3)2 qo,)+~(u:+4)=1. +&T A

or

S. TOWN and N.

370

ZABARAS

1500

z

1000-

4

6

MAXIMUM STRESS FPF

l

MAXIMUM STRESS LPF

l

HOFFMANFPF

X

HOFFMANLPF

0

EXPERIMENTAL (231

si k t;; 500-

0

r 40

I 20

0

I 60

I 80

100

ANGLE (DEGREES)

Fig. 2. First and last ply failure curves for (0,/O,/-O,), plate using the maximum stress and Hoffman failure criteria.

The mode of failure is said to be matrix mode if either the transverse contributions or the shear contributions are the largest. If the fiber direction contributions are the largest, the failure is said to be fiber mode. Once the mode of failure is known, the stiffness may be reduced in the appropriate way. Since all elements used have four Gauss points the element stiffness properties associated with a given lamina and direction are reduced by one-fourth of their original value for each location within the element that fails. A fiber mode failure at all four Gauss points of an element would reduce the stiffness matrix of the failed lamina within that element. Consequently, eqn (4a) would become 0

0

0

00

0

Q22

0

00

0

0

Q44

0

0

0

0

00

0

0

0

00

reduced

to

zero

fbr

the

=

failed

0

0

0

00

0

0

0

00

0

0

0

0 0 Qss 0

0

00

0

(28)

0

because a2, ad, and c6 (a,, bTZ, and cLT) would be reduced to zero for this failed lamina. Again, intermediate stiffness values would be prescribed if failure were to occur at 1, 2, or 3 Gauss points. Delamination is the tinal mode of failure. It is characterized by the interlaminar stresses acting between adjacent layers. An interface of two adjacent layers is identified as a delamination failure if either [ 141

0

$=l

(4 + w2

or

= 1

’

sz

because 6,) c5, and rr6_ (aL, _aLr,_._ a ncI bLT) would

be

Q,, 0 0

lamina.

If

failure was detected at less than four Gauss points, the appropriate stiffness components would be given intermediate values between those of eqns (4a) and (27). The intermediate values would be found by proportionally reducing the stiffness components by the fraction of failed area in the element. A matrix failure at all four Gauss points of an element would reduce the stiffness for that element in a different manner. Equation (4a) would then become

(29)

where S, is the through-thickness shear strength. A delamination failure at all four Gauss points of an element would have yet another effect on the stiffness matrix of eqn (4a). For both of the lamina adjacent to the delamination, the new stiffness would be as follows: QH

Q,2

0

0

0

Q12

Q22

0

0

0

0

0

00

0

=5

0

0

00

0

u6

0

0

0

02 *4

=

0

Qss

371

Progressive failure in laminated composite plates

q

HASHIN FPF

l

HASHINLPF

l

LEEFPF

z

X

LEELPF

%

0

EXPERIMENTAL 1231

8 III E 500 -

0

20

40

60

60

100

ANGLE (DEGREES)

Fig. 3. First and last ply failure curves for (0,/O,/-0,),

In this case the stresses u), cd, and o5 (a,, on, and ati) are reduced to zero in both laminae adjacent to the failed interface. Obviously, when all the members of the element stiffness matrix have been reduced to zero the element makes absolute no further contribution to the structure and is considered to have undergone total failure. The composite analysis program written to carry out the required calculations stores the lamina stiffness properties in a three-dimensional array. The three indices of the array are property number, layer number, and element number. These stiffness properties can, therefore, be reduced one at a time as failure of a specific location is determined. The reduced material property matrices in global coordinates, Qj must be recalculated every time new stresses are calculated to insure that the stresses take the failed lamina into effect. After first ply failure the stresses are recalculated at the FPF load using the newly calculated element stiffness matrices. These new stresses are once again inserted into the failure criteria to determine if further failure will occur. Assuming this happens, the stiffness is reduced appropriately and the stresses are recalculated once again at the FPF load. These ‘equilibrium iterations’ continue until no further failures occur at the FPF load. The load is incremented for the second time. Again, the stresses are calculated, and the laminate is checked for failure. Stiffness reduction and load increment continue until the stiffness matrix has been reduced to zero for all laminae at a single (x, y) location. This is considered to constitute LPF or the ultimate strength. THEORY-EXPERIMENT CORRELATION Vniaxial tension

The numerically calculated FPF and LPF results will be compared with experimental results obtained

plate using the Hashin and Lee failure criteria.

by Soni for uniaxial tension of two composite laminates [23]. The analysis was performed using a one eight-noded element plate model subjected to uniaxial tension. The laminate studied consisted of 24 layers of T300/5208 graphite-epoxy having an aspect ratio (S = a/h) of 150. The material and strength properties for T300/5208 are as follows: vu = 0.3

Y- = 210 MPa

vu = 0.3

S, = 93 MPa

v*, = 0.3

ST = 93 MPa

Gu = 6.46 GPa

S, = 93 MPa

Gu = 6.46 GPa G,, = 6.46 GPa. The layup studied is the (0,/O,,/-f3,),. The laminate was studied using all five of the failure criteria discussed earlier. Figure 2 shows the FPF and LPF curves using the maximum stress and Hoffman criteria. The experimental results fall near the failure curves with the exception of the 15” experimental point which is considerably lower than the predicted values. Soni states that the failure mechanism of the 15” specimen was extensive delamination at the free edges, followed by tensile fiber failure. Note that the current analysis neglects free-edge effects. It is interesting to note that for both failure criteria the FPF and LPF values coincide for angles less than 30”. This is due to the fact that the fibers run primarily in the direction of the applied load for these angles. The first failure in each case was fiber mode failure of the 0” and -6” laminae. Since the FPF values in each case are greater than the fiber direction strength of the two remaining 0” laminae, complete

312

S. TOLS~Nand N. ZABARAS

0

EXPERIMENTAL (241 40 X STRESS

60

80

(MPA)

Fig. 4. First and last ply failure curves for a (O/90),plate using the Lee failure criterion.

failure follows immediately after first ply failure. The LPF curve determined using maximum stress criterion levels off for angles of 0 greater than 45”; however, the LPF curve determined using the Hoffman criterion continuously decreases with increasing angle. This is because the Hoffman criterion is a combined stress criterion and includes transverse stress contributions not present in the maximum stress criterion. The results obtained from the Tsai-Wu criterion are very close to those obtained using the Hoffman criterion. Therefore, they are not shown. Figure 3 shows the first and last ply failure curves determined using the Lee and Hashin failure criteria. Again the FPF and LPF curves coincide

0 0 0

for angles less than 30” and the LPF curves are horizontal lines from 60 to 90”. The predicted strengths are again less than the experimental values in the 60-90” range by a small percentage. The results using both criteria are again close to the experimental data with the exception of the 15” data point. In all, Hashin’s criterion seems to give results which are slightly less accurate than the Lee and maximum stress criteria, but more accurate than Hoffman and Tsai-Wu criteria for this particular laminate type. Interestingly, the FPF curves for all failure criteria are quite similar, and differ by a maximum of 12% (at 15’). The LPF curves are less similar and differ by as much as 30% (at 90’).

EXPERIMENTAL [24] I 20

I 40 X STRESS

1 60

80

(MPA)

Fig. 5. First and last ply failure curves for a (O/k45/90), plate using the Lee failure criterion.

313

Progressive failure in laminated composite plates

•I FPFCURE

0

l

LWCURVE

0

EXPERIMENTAL [24]

X STRESS

80

60

40

20

(MPA)

Fig. 6. First and last ply failure curves for a (O/*45), plate using the Lee failure criterion.

Biaxial tension The laminates are all of boron-epoxy following properties:

having the

E, = 204 GPa

X+ = 1260 MPa

ES= 19GPa

X- = 2500 MPa

E, = 19 GPa

Y+ = 61 MPa

vu = 0.25

Y- = 202 MPa

v,, = 0.25

S, = 76 MPa

vZ3= 0.25

S, = 76 MPa

G,* = 5.6 GPa

S, = 76 MPa

G,, = 5.6 GPa GZ3= 5.6 GPa. One eight-noded element plate model was used to perform the analysis. The aspect ratio of all laminates analyzed was S = a/h = 100. The theoretical strength predictions were compared to experimental biaxial strength results obtained from [24]. The Lee criterion was chosen for this study. The first layup studied is the (O/90),. Figure 4 shows the first and last ply failure curves determined at various biaxial load ratios, PX : P,. As can be seen, significant scatter is present in the experimental data. For x-stresses of 0, 18, and 32 MPa two or more y-stresses are shown. This type of material behavior seems quite unusual and must be due to processing defects and testing inconsistencies. Figure 5 shows the first and last ply failure curves for a (O/+45/90), laminate. The theoretical FPF and LPF curves are much closer to bounding the experimental results for this case. One experimental data point is still outside the curves, but only by a small margin. Thus, the theory-experiment correlation for this laminate is very good.

Figure 6 shows the first and last ply failure curves for a (O/&45), laminate. The theoretical FPF and LPF curves for this case also compare quite well with the experimental data. An interesting feature of this graph is the fact that much of the FPF and LPF curves are coincident. This is due to the fact that the +45” layers make up a large volumne percentage of the laminate. If they fail, the remaining layers are no longer capable of supporting the load. Also, the failure of a f45” ply in either the matrix of the fiber mode implies that the layer has failed in both the x and y directions. Thus, biaxial load ratios greater than 1 : 1 lead to matrix failure in the f45” layers which weakens the structure to the point that failure of the remaining layers follows immediately. Plates under transverse loading The composite analysis computer program will be used to investigate the first and last ply failure loads of a O/90/90/0 laminate subjected to a sinusoidally distributed transverse pressure as specified in eqn (20). The plate is simply supported and the boundary conditions are as given earlier. The material and strength properties for the carbon-epoxy composite used are as follows: E, = 180GPa

X+ = 1500 MPa

Ez = 10.6 GPa

X- = 1500 MPa

E, = 10.6 GPa

Y+ = 40 MPa

vn = 0.28

Y- = 250 MPa

v,, = 0.28

S, = 68 MPa

v2, = 0.28

S, = 68 MPa

G,, = 7.56 GPa Gu = 7.56 GPa Gzs = 7.56 GPa.

Sz = 68 MPa

S. TOLS~N and N. ZABARAS

314

q

6

0

I

0

I

20

I

I

40 ASPECT

7DoFLPF

0

SDOFFPF

X

5DOFLF’F

I

60

80 RATIO,

7DOFFPF

l

I

100

1 3

S

Fig. 7. First and last ply failure curves for a O/90/90/0 simply supported plate with transverse pressure load.

Only a quarter of the plate needs to be modeled due to the geometric symmetry of the problem. Four eight-noded elements were used to insure reasonable accuracy in the stress calculations. The Lee failure criterion was used in the analysis. Figure 7 shows the first and last ply failure curves for the O/90/90/0 laminated plate as a function of aspect ratio. The failure loads plotted in the graph have been normalized as FPF* =T

and

LPF*=ik!?& Since no experimental data could be found to compare with the failure predictions of the composite analysis program, a new failure analysis was carried out using a first order shear deformation plate theory to determine the stresses within the laminate. The set of curves designated ‘5 DOF’ represent the failure results obtained using this particular plate theory. The five degree of freedom shear deformation plate theory used in the analysis is essentially equivalent to Reddy’s penalty plate bending theory and provides equal results [9]. The set of curves designated ‘7 DOF represent the failure results obtained using the higher order plate theory developed in this paper. A comparison of the stresses obtained for the two plate models was presented previously in Table 1. The stress values calculated using the first order deformation theory were found to be less accurate than the stress values calculated using the higher order theory developed in this paper when compared to the three-dimensional elasticity solution. There-

fore, the results of the failure analysis determined using the five degree of freedom first order shear deformation model are expected to be less accurate than the results determined using the seven degree of freedom higher order shear deformation model. Further, the results found using the two theories are expected to converge as aspect ratios increase, because the stress values were found to converge as aspect ratios increased. The results shown in Fig. 7 converge as expected for the FPF curves, but the LPF curves maintain about an equal percent difference throughout. The FPF curves have a maximum difference of 23% at an aspect ratio of 4, and converge to within 1% at an aspect ratio of 100. The LPF curves do not seem to converge or diverge. The maximum difference of 19% occurs at an aspect ratio of 4 and the minimum difference of 15% occurs at an aspect ratio of 100. Since the stresses calculated by the two models at an aspect ratio of 100 are nearly equal prior to failure, the above differences must be associated with the sensitivity of the two models to stiffness changes that occur after FPF. The progression of failure within the O/90/90/0 laminated plate is shown in Figs 8a through 8d for an aspect ratio of S = 10. Figure 8a shows that FPF occurs in the area closest to the center of the plate at the bottom 0” lamina. The initial failure is matrix mode caused by tensile stresses. Figure 8b shows the additional failure sights that occur during the ‘equilibrium iterations’ at the FPF load. The small area of initial failure causes a sufficient reduction in the strength of the laminate to trigger failure of a much larger region. The progression of failure stabilizes after three iterations, allowing the load to be incremented. Figure 8c shows the additional failure that occurs immediately after the increased load is applied. Failure has now spread to the interior 90”

37.5

Progressive failure in laminated composite plates

q

Matdx

q M&lx&de Falkire q Matiix Mode Failun

Mode Fallwe (T%=I@W)

Layer 3

Matrix Made Failure (Tensile)

q Matrix Mode Failurn

(Compressive)

q Fltw Mcda Failure (Tensile) q Flbsr Mode Failure (Comprasslve)

Layer 4

Fig. 8a, Failure progression of a O/90/90/0 plate under transverse sinusoidal pressure (FPF-no~aiized Ioad 1513 MPa). &j

(Team)

(Compressive)

q Fiber Mode Failure (Tensile) •m

Layer 3

Layer 4

Fig. 8b. Failure progression of a O/~~~~O plate under transverse sinusoidal pressure (after eq~librium iterationsnormalized load 1513 MPa).

Fiber h+ode Failure (Compulsive)

Layer1 Layer 2

Layer 3

Layer 4

Fig. 8c. Failure progression of a O/90/90/0 plate under transverse sinusoidal pressure (after first load incrementnormatized load 1590 MPal.

layers and a large percentage of the laminate has failed. Due to the extensive damage present, the laminate fails compietely during the equilibrium iterations. The failure pattern at the point of last ply failure is shown in Fig. 8d. The region of zero stiffness through the thickness that characterizes LPF occurs near the center of the plate in the same location as FPF occurred. Thus it has been shown that composite analysis based on cumulative stiffness reduction techniques can be used to estimate the first and last ply failure loads for an arbitrary Iaminated plate. The analysis works equally well regardless of the laminate layup or

Layer 3

Layer 4

Fig. 8d. Failure progression of a O/90190/0 plate under transverse sinusoidal pressure (LPF-no~alized load 1590 MPa).

thickness. The progression of damage between FPF and LPF can also be estimated based on the same stiffness reduction principles. While the theoretical progression of damage is not expected to compare exactly with the progression of damage in a real material, it is expected to show the major trends.

SUMMARY This paper has presented a computational model for dete~ining the ultimate strength of an arbitrary laminated composite plate. A new higher order shear

S. TOLKINand N. ZABARAS

376

deformation plate theory was developed. The theory utilizes seven degrees of freedom at each node. An improvement in the accuracy of the transverse shear stresses was obtained by calculating these stresses using three-dimensional elasticity equilibrium equations. The method proved so successful that only the inplane stresses were calculated from constitutive equations in the final form of the computational model. The composite failure analysis computer program used to determine first and last ply failure of a laminated composite plate was developed based on the seven degree of freedom higher order shear deformation plate theory. The program utilizes a system of failure mode determined stiffness reduction to decrease the strength of a failed laminate. The FPF and LPF strengths determined by the analysis are expected to bound the actual failure strengths. Five failure criteria from various categories were used in the analysis. For the cases studied, the Lee criterion gave the best results. The model is quite effective in providing these limits, and would be a valuable aid in the design of laminated composite plates. Acknow~ledgemenr-The computational part of this project has been supported with a grant from the University of Minnesota Supercomputer Institute. REFERENCES I. E. Reissner, The effect of transverse shear deformation on the bending of elastic plates. ASME J. appl. Mech.

67, A69-A77 (1945). 2. E. Reissner, On bending of elastic plates. Q. appl. Math. 5, 5567

(1947).

3. R. D. Mindlin, Influence of rotatory inertia and shear on flexural motions of isotrouic. elastic mates. ASME J. appl. Mech. 18, 31-39 (195;). _ 4. M. Levinson, An accurate simple theory of the statics and dynamics of elastic plates. Mech. Rex Commun. 7, 343-350 (1980).

5. E. Reissner and Y. Stavsky, Bending and stretching on certain types of heterogeneous aelotropic elastic plates. ASME J. appl. Mech. 28, 402408 (1961). 6. P. C. Yang, C. H. Norris and Y. Stavsky, Elastic wave propagation in heterogeneous plates. Inr. J. Numer. Merh. Engng 2, 665484 (1966). I. J. N. Reddy, A simple higher order theory for laminated composite plates. ASME J. appl. Mech. 51, 745-152

(1984). 8. M. Murthy, An improved transverse shear deformation theory for laminated anisotropic plates. NASA Technical Paper 1903 (1981).

9. J. N. Reddy, A penalty plate bending element for the analysis of laminated anisotropic composite plates. Int. J. Numer. Merh. Engng 15(8), 1187-1206 (1980). 10. N. S. Putcha and J. N. Reddy, A refined mixed shear flexible finite element for the nonlinear analysis of laminated plates. Compur. Srrucr. 22, 529-538 (1986).

11. R. M. Jones, Mechanics of Composite Materials. Scripta, Washington (1975). 12. Z. Hashin, Failure criteria for unidirectional fiber composites. ASME J. appl. Mech. 47, 329-334 (1980). 13. B. D. Agarwal and L. J. Broutman, Analysis and Performance of Fiber Composites. John Wiley, New York (1980). 14. J. D. Lee, Three dimensional finite element analysis of damage accumulation in composite laminate. Compur. Srrucr. 15, 335-350 (1982).

15. W. C. Hwang and C. T. Sun, Failure analysis of laminated composites by using iterative three-dimensional finite element method. Compur. Srrucr. 33,4147 (1989).

16. J. N. Reddy and A. K. Pandey, A first-ply failure analysis of composite laminates. Compur. Srrucr. 25, 371-393 (1987).

17. J. J. Engblom and 0. 0. Ochoa, Finite element formulation including interlaminar stress calculations. Comput. Srrucr. 23, 241-249 (1986).

18. J. J. Engblom and 0. 0. Ochoa, Analysis of progressive failure in composites. Comp. Sci. Technol. Zs, 87-102 (1987).

19. S. A. Tolson, Finite element analysis of progressive failure in laminated composite pltes. M.S. thesis, University of Minnesota, MN (1990). 20. J. N. Reddy, A note on symmetry considerations in the transient response of unsymmetrically laminated anisotropic plates. Inr. J. Numer. Merh. Engng 20, 175-180 (1984). 21. N. J. Pagano,

Exact solutions for rectangular bidirectional composites and sandwich plates. J. Camp.

Mater. 4, 20-34 (1970). 22. M. N. Nahas, Survey of failure and post-failure theories of laminated fiber-reinforced composites. J. Camp. Technol. Res. 8(4), 138-153 (1986). 23. S. R. Soni, A new look at commonly used failure

laminates. 24th AIAA/ theories in composite ASME/ASCE/AHS Structures, Structural Dynamics and Materials Conference, Procedures, pp. 171.-179 (1983). 24. G. C. Sih and A. M. Skudra, Failure Mechanics of Composites: Handbooks of Composites, Vol. 3. Elsevier Science, Amsterdam (1985). 25. S. W. Tsai and E. M. Wu, A general theory of strength for anisotropic materials. J. Camp. Mater. 5, 58-80 (1971).

26. R. Hill, The Marhemarical Theory of Plasriciry. Oxford University Press, London (1950). 21. 0. Hoffman, The brittle strength of orthotropic materials. J. Comp. Marer. 1, 200-206 (1967).