- Email: [email protected]
Failure analysis of laminated composites by using iterative three-dimensional finite element method
FAILURE ANALYSIS OF LAMINATED COMPOSITES BY USING ITERATIVE THREE-DIMENSIONAL FINITE ELEMENT METHOD W. C. HWANG~ and C. T. SUNS tChung_Shan Institute of Science and Technology, Lung-Tan, Taiwan $Department of Aerospace Engineering, Mechanics and Engineering science, University of Florida, Gainesville, FL 3261 I, U.S.A. (Received 2 May 1988) Abstmct-In this paper, a failure analysis of laminated composites is accomplished by using an iterative three-dimensional finite element method. Based on Tsai-Wu failure theory, we first propose three different modes of failure, namely, fiber breakage, matrix cracking and delamination. The first-ply failure load is then evaluated. As the applied load exceeds the first-ply failure load, localized structural failure occurs and the global structural stiffness should change. We modify the global stiffness matrix by taking nonlinearity due to partial failures within a laminate into consideration. The first-ply failure load is analyzed by using an iterative mixed field in solving the linear part of the finite element equations. The progressive failure problem is solved numerically by using Ne~on-Raph~n iterative schemes for the solution of nonlinear finite element equations. Numerical examples include a three-layered cross-ply (O/90/0) E-glass-epoxy and an angle-ply four-layer Thomel 300 graphite-934 resin epoxy laminates under uniaxial tension in both cases. First-ply failure loads as well as the final loads are evaluated. Good correlation between analytical re.su1t.sand experimental data are observed. Numerical results also include the investigation of composite specimens with a centered hole, under uniaxial tension, First-ply failure loads and final failure loads are obtained by using the finite element program. Excellent correlation with the experimental data is observed.
failure problem. The increased use of high-strength, low-weight, fi~r-~infor~ composite materials, especially in aircraft and space vehicles, has heightened interest in understanding and characterizing their behavior. The trend toward higher design levels has led to the need to understand failure initiation and to predict strength. The overall effectiveness of the finite element analysis depends to a great extent on the numerical algorithms used in the solution of the algebraic system of equations resulting from the discretization process. As the finite element method is applied to larger and more complex problems, it becomes increasingly important that the solution process remains economical and accurate. We believe that finite element equations can be solved more cheaply than at present, but the failure analysis of multilayer composite laminates using thr~~mensional tinite element models often requires large computational resources [3]. The necessity to reduce computational cost or to run in minicomputers provided the incentive to formulate the effective algorithms presented in this research with the iterative mixed field method, which has already been described in [4]. Usually the sofution of these nonlinear equations is obtained by iteration [S]. In this paper, attention is focused on the failure analysis of composite laminates by using the iterative mixed field method in the linear part and a variable incremental load step scheme with the modified Newton-Raphson iteration in the nor&near part (progressive failure process). These procedures are carried
1. INTRODUCFION
The standard finite element discretization following system of algebraic equations:
yields the
in which (V) is the dispI~ment vector and considered as the basic unknown, (R} is the vector of applied “loads” and [K] is the assembled “stiffness” matrix. The interpretation of these quantities varies according to the physical problem under consideration. For structural appli~tions, the terms “load” and “stiffness” are directly applicable. In dynamic analysis, the vector {R) would include the inertia and damping forces. Hence, a dynamic analysis is basically a static analysis including inertia effects. We assumed that disp~a~ments of the finite element assemblage are small and that the material is linearly elastic. In addition, we also assumed that the nature of the boundary conditions remains unchanged during the application of the loads. Equation (1) corresponds to a linear analysis of a structural problem because {V) is a linear function of(R). For high-strength fiber-reinforced laminated composite, the fiber and matrix generally behave elastically until they have localized structural failures [I, 21. At this time the global structural stiffness should change; the response is only linear prior to localized failures, The nonlinearity comes from partial failures within a laminate. This situation arises in the analysis of a 41
W. C. HWANG and C. T. SUN
42
out on the numerical failure analysis of a T300/934 graphite-epoxy (GE) composite, using three-dimensional finite element models, and the results are compared with experimental results obtained at NASA Lewis Research Center.
2. MODIFIED NEWTON-RAPHSON METHOD fN COMPOSITE MATERIALS
The most popularly used iterative schemes for the solution of nonlinear finite element equations are forms of Newton-Raphson iteration [5]. This is often used in conjunction with an incremental loading procedure. The discretization of nonlinear problems can be written in the form (Y(U)} = {RI -WWl{U
-09
tb)
(2)
where we define the residual {r(u)} as the vector of the unbalanced nodal forces at any given deformed state (say IV)). The nonlinear soiution process consists of searching for a vector {CJ} that renders the residual as small as possible. The deformed state {V) is an exact solution of the discretized equations if, and only if, the unbalanced nodal forces vanish. Since the governing eqns (2) are nonlinear and generally of high order, they cannot be solved explicitly. An iterative algorithm, such as Newto~Raphson, can be used to solve the equation numerically. It is illustrated graphically in Fig. l(a) for problems with one variable u [6]. The instantaneous stiffness matrix [K,] is represented by the tangent to the curve at each step. For any known {U’), finite element modeling and analysis provides the numerical value of {y(V)} and [K,(V)] of the residual with respect to {vi) at the ith load step. Since evaluation of the stiffness matrix requires many calculations, it may be desirable to evaluate it less frequently. Figure l(b) shows the result of the Newton-Raphson algorithm with this modification. The choice of iteration steps, when the stiffness matrix should be updated, depends on the degree of nonlinearity in the system response, i.e. the more nonlinear the response, the more often the updating should be done. By using stress transformation equations we can determine the stress along the principal material axes in each lamina. Once the stress field is known, the application of a lamina failure criterion is appropriate at each “point” for initial failure analysis and at each “element” (averaged stress in each element if finite element method is used) to predict ultimate strength. In many cases, laminates have considembIe strength remaining after initial failure. If there are no failures in a load step, then the applied loading is increased until failures occur. Otherwise, the material is said to fail by the failure m~hanism. When a failure takes place, a change in the stiffness matrix, [AD], due to the localized structural failure is calculated. Using the change in
Pf PIR
fc) _~. 0
Ul
“2
U
U3U4U”
Fig. 1. Graphical iilustration of (a) Newton-Raphson, (b) modified Newton-Raphson, (c) iterative mixed field algorithms for a one variable u.
stiffness and the known strains, an imbalance vector {AQ} is found.
Ioad
(3) where n is the number of stiffness degradations during the iteration. The steps used in the modified Newton-Raphson iteration are for i = 1,2,3, . . . (for each iteration) [6]. (1) Solve the increased displacements
{AU}i
[K](ALr}j = (AQ}j. (2) Compute the new displacements {V)‘= {U)‘-‘+fAU)‘.
(4) (U)j (5)
(3) Execute convergence test. This test is the same as in the iterative mixed field method used in [4]. If convergence is not reached, this gives the required change in strains and stresses. (4) Extract nodal displacements {V}’ from {U}‘. (5) Determine the strains and the redistributed stresses. (c)‘= [e](uqi
[II*]’ = [b] - [AD]
(6) (7)
43
Failure analysis of laminated composites
(u)‘= p*l’((c)i-
(Lo])
{a}:=[rl{~}i,
(8) ments, since it permits identi~~tion (9)
where [b*]’ becomes a ~st-fails rotated elasti~ty matrix. The stresses, (a)‘, are changed by the degradation of matrix [D *I’. (6) Employ the failure criterion to check whether there are any failures or not. If there are no new failures, then add the incremental load. If the convergence rate becomes slower, a new global stiffness matrix, [K], is reassembled. For each iteration, the above-mentioned procedure is repeated until a converged state is reached or the laminate experiences gross failure. If we keep the finite element model and the structural stiffness unchanged, this iterative procedure is very similar to the iterative mixed field method, which is used in solving the linear part of the failure analysis. This similarity makes the programming implemented in the nonlinear analysis easier. The illustration of the iterative mixed field method with one variable u is shown in Fig. l(c).
3.
LAMINATE
FAILURE CRITERION
Failure criteria have been in use for structural materials for centuries. They can be divided into non-interactive criteria, such as the maximum stress or maximum strain; and interactive criteria, such as quadratic approximation. A common approach to formation of a failure criterion is to assume a quadratic polynomial in stress[7]. It is the simplest choice which can adequately describe the experimental data to consider the combined effects of stresses on failure. The quadratic failure criterion for composite materials has been proposed by many workers for at least 20 yr. A general quadratic, including linear terms, is the Tsai-Wu failure criterion [l). In this criterion, faihtre is assumed to occur when the following equations are satisfied: F,a, + Flja,al = 1 (i,j = 1,2, . . .6),
(10)
where F, and Fu are lamina strength tensors of the second and fourth order, respectively. The usual contracted stress notation is used except that Us = rz3, a, = tft, and a, = tn. The linear stress terms in (10) can be used to account for the difference between tensile and compressive strengths. Although it seems physically attractive [ 11,some components in the tensor polynomial require difficult and expensive biaxial testing for their determination. A major shortcoming of this quadratic failure criterion is that it predicts the imminence of failure but says nothing about the different failure modes and does not specify how the composite fails. This is of fundamen~l importance for design considerations of composite structural ele-
of the nature of the failure and finite element analysis of progressing structural failure. These problems can be avoided if the different failure modes of the fiber composite and the primary stresses contributing to them are identified and each mode is modeled separately by a quadratic. The principal failure modes are (1) fiber mode: tensile fiber mode is described by fiber breakage and compressive fiber mode is described by fiber buckling; (2) matrix mode: tensile matrix mode is described by plane failure surface parallel to fibers with cr + g3 > 0 and compressive matrix mode with a, + or < 0; and (3) delamination mode: produced by the interlaminar stresses. Experimental evidence for some of these modes obtained with off-axis specimens of various composites has been described in [S] and [9]. Each of these may be expected to have a different failure criterion; and further, these criteria may be of different form for tensile and compressive stresses. Failure criteria based on this approach have been presented by Hashin [IO], in which the general three-dimensional stress is considered. Each of the failure modes described can be modeIed separately by a quadratic polynomial containing the appropriate stress components. The advantages of this approach for the biaxial test results are not needed and the failure mode is identified. The following criteria make use of the following material strength parameters, which are fiber direction (longitudinal) tensile strength fiber direction compressive strength transverse tensile strength transverse compressive strength shear strength in the 1-2 plane shear strength in the 1-3 plane shear strength in the 2-3 plane interlaminar (through-the-thickness) tensile strength interlaminar shear strength. 2, and S, are assumed as the tensile strength and pure shear strength of the matrix, respectively. The faiiure criteria are separated into three different modes. For fiber breakage, the criteria are
at failure for 6, > 0, and
0
2
01
x,
(12) =I
W. C. HWANGand C. T. SUN
44
at failure for u, c 0. For matrix cracking, the criteria are
+(zJ+($r=
1 (13)
at failure for cr + c3 > 0, and
$[($-
Il(o,+o,)+&(a,+~,)’
at failure for g2 + 6) < 0. Delaminations are predicted utilizing a simple quadratic interactive formula [3]. The delamination occurs if either a,zZ,
or
(tf,+r
)-2s :3:_
I’
Yet another form used to identify delamination follows [ 111:
0
7:3+7:3>
1
s2+7/. I
(1%
is as
the response of composite laminated structural elements up to collapse. The reason is that failure criteria are not valid after FPF [7], strictly speaking. Until this impasse is overcome, the design of laminated structural elements with no high stress concentrator is likely to be based on the onset of initial failure in an individual lamina of the laminate, even though such a criterion may be very conservative from the strength standpoint [l]. The case of biaxial loading is a, in the x direction and ka, with a fixed load factor k in the y direction. Based on the formulae of the failure criteria, the failure mode, the prediction of the externally applied loads which will initiate failure at an opening, and the location of the failure can be obtained. The progressive failure occurs when the applied loading is increased beyond the initial failure load. The treatment of the case of bending problems is similar. The mode of failure determines how the stiffness matrix will be modified to reflect the effects of the predicted damage. The way in which the lamina properties are changed depends upon the type of ply damage. A simple, post-failure reduced stiffness approach is used here. From laminate constitutive equations, if the damage includes fiber failure in an element, D,,, Di2, D13,D,, and D, are set equal to zero for that element in this study. A change in the stiffness matrix of that element occurs.
(16)
I
In the prediction of failure, (15) and ( 16) are evaluated at the critical locations in the ply interfaces of a laminate only. It is very useful for progressive damage analysis to use a finite element method to analyze a fiber composite structure [lo]. Such a computational procedure starts with a finite element stress analysis. Certain elements will fail when their stresses satisfy the failure criterion. But in order to continue the analysis, it is essential to know how the elements have failed by fiber rupture, matrix crack, or delamination from an adjacent lamina. With this information, it is possible to set up a new finite element configuration with appropriate stiffness changes or new free surfaces, and the analysis can proceed to identify the new failures.
4,
42
D,3
0
0
0
0
0
0
0
0
0
0
0
0
4s
0
[AD] = Symmetry
D,
1
1 .
(17)
If matrix failure is the only failure mode, it is assumed that the fiber can still carry load; thus, only the corresponding transverse and shear properties of the stiffness matrix are removed. The stiffness D,* = D,, = D2, = Du = D, = 0, or
[AD] =
4. INITIAL FAILURE
ANALYSIS AND THE PREDICTION OF STRENGTH
For a laminate which is loaded by in-plane forces and/or bending moments, the stresses in the layers can be obtained. If there is no bending and in-plane loads are uniform distributions on the edges, then this is a special case where the stresses in each layer are constant and planar. The layer for which the criteria are satisfied by the lowest external load will define the layer that fails initially and the nature of the failure (failure mode). This is called first-ply failure (FPF). Despite continuing research effort, satisfactory analytical techniques are not yet available for predicting
If damage includes fiber and matrix modes, then the entire stiffness matrix of that element is removed, or [AD] = [D], The components of [AD] are calculated by use of where A6, terms are substituted for D4 terms. Also, damage in one element will not affect the stiffness of other elements. If delamination occurs, we remodel a new finite element configuration with a new free surface at the ply interface. The area around the delamination fronts can be modeled by singularity crack tip elements [ 121to reduce the need for extremely fine meshes. Since the delamination mode is checked on each element basis, the average stress as calculated
Failure analysis of laminated composites
with the singular element differs little from that calculated with a conventional element. However, this mode of delamination can be easily eliminated by using the techniques which are discussed in [13]. The laminates analyzed were 50.8mm wide by 457.2 mm long specimens, which were tabbed at each end with 50.8 mm wide beveled aluminum tabs and loaded in uniaxial tension. The following types were studied: three-layered cross-ply (O~~/O) E-gIassepoxy solid specimens, hi = h3= 0.127 mm, h, = 1.27 mm; four-ply panels with orientation ( f 0), Fiberite 1034E prepreq (934 resin matrix with thornel 300 graphite fibers) solid specimensand specimensnotched with centered 6.35 mm diameter through holes. The nominal thickness of one ply is 0.127 mm. In all cases, uniaxial tension was applied at one end while the other end of the specimen was fixed. The material properties of the first simple E-glassepoxy lamina are [l] E, = 53.8 GPa (7.8 x 10”psi), E2 = Ex = 18 GPa (2.6 x 10’ psi), Gn = Gu = 8.6GPa (1.25 x 106psi), G,, = 7.2 GPa (1.04 x lo6 PSI), vu = v,j = v, = 0.25, a, = 3.5 x lO+/“F, a2= a3 = 11.4 x 10a6pF, X,=X,= 1034MPa (150 ksi), Yr = 27.6 MPa (4 ksi), Y, = 138MPa (20 ksi), S,, = S,, = S,, = S, = 41.4 MPa (6 ksi), 2, = 48.3 MPa (7 ksi), AT = -2OO“F, AM = 0.
45
Fig. 2. A comparison of the prcdii and measured angle-ply laminates ( f f?), failure loads in graphite-epoxy composites.
The thy-dimensional finite element model used for solid specimens is a 20 x 4 x 1 uniform mesh of eight-node elements in each layer. Figure 2 shows the comparison between the numerical results and experimental results [15]. Because the specimens are very thin, they can be considered as plane stress state problems. The strengths in Fig. 2 are essentially those of FPF predictions. The laminates fail to carry more load after the FPF. The dominant modes of initial failure predicted by laminate failure criteria are fiber breakage in solid laminates with ply orientation (O), to ( f 30), and matrix cracking in solid specimens with (19) ply orientations > 30”. The finite element mesh used to model the compoThe solutions are obtained by using a quadrant of site specimens with centered holes is as shown in the laminates because of symmetry. The finite element Fig. 3. Results of the numerical analyses, which are model is a 10 x 2 x 1 uniform mesh of eight-node shown in Fig. 4, agree well with experimental results elements in each layer. Results of the initial failure 1151.The predicted strengths are higher than the initial stress, a,,, and the ultimate failure stress, efi are failure stresseswhich are induced by the stress concene,, = 23.2 MPa and cr/= 195MPa. Those results are trations around holes. The laminates continue to carry reasonably close to the experimental values where a0 load after the initial failures occur. It also indicates is about 23.4 MPa and ur is about 200 MPa [14]. The that the results are conservative in their prediction of initial failure mode, which is the tensile matrix mode, strengths. Table 1 shows that the initial delamination is predicted to occur in the 90” layer. failure mode will occur when the number of plies per The assumed material properties of the second exlayer, n, increases from two to five. Delamination, ample, T300/934, are [ 15; Irvine, personal communicaas a form of internal damage, is det~men~l to the tion, 1985; Malvem, personal ~mmunication, 19851 laminate strength. This mode of failure can be easily eliminated when it is predicted numerically. The E, = 138GPa (20 x lo6 psi), predicted initial failure stresses, a,, are listed in Fig. 5. Et = Es = 10.3 GPa (1.5 x lO”psi), The values of a, will decrease as n increases from two G,* = G,, = 6.6 GPa (0.95 x 10'psi), to six. Figure 6 indicates the maximum &,, which is G2,= 2.6 GPa (0.37 x 106psi), seen in the (-+45), laminate. This is the reason why v,z = v,j = 0.21, v*3= 0.47, Xr = 1393MPa (202 ksi), Xc = 1447.9MPa (210 ksi), Yr = 44.8 MPa (6.5 ksi), Y, = 172.4MPa (25 ksi), S ,2= S13= S, = 62.1 MPa (9 ksi), S, = 34.5 MPa (5 ksi), Fig. 3. The top view of the finite element mesh of the Z, = 49 MPa (7.1 ksi). (20) angle-ply laminate specimen ( f O), with a through hole.
W. C.
46
HWANG
and C. T. SUN
Table 1. The predictions of initial delamination failure mode of (t?,,/-6,)s G/E laminates with holes Initial failure n 0 2
3
4
5
(deiree) 1.5 30 45 15 30 45 15 30 45 15 30 45 60
(gi) 89.8 35.2 31.6 69.6 27.1 24.3 57.3 22.6 20.2 50.2 19.8 17.7 22.4
Approx.? location
Approx.7 locations
(dek+e)
Maximum SCF
(de:ee)
-75 -75 -75 -75 -75 -75 -15 -15 -75 -75 -75 -75 -60
5.2 5.3 5.7 5.3 5.4 5.7 5.3 5.4 5.8 5.4 5.4 5.8 5.1
-105 -75 -60 -105 -75 -60 -105 -75 -60 -105 -75 -60 -60
Approx.? locations Maximum rd
(de&e)
-0.34 -0.95 -1.10 -0.45 - 1.23 -1.42 -0.55 -1.48 -1.71 -0.63 -1.70 -2.00 -1.40
-75 -15 -15 -75 -75 -75 -75 -75 -75 -75 -15 -75 -60
t Because of symmetry, there are always two locations on the circumference; one is located at B = tl (or &2} and another is at @= OL+ 180 (or #3,y).
the initial failure stresses are lower around the (f45), laminate. The numerical results presented herein indicate that a three-dimensional finite element progressive failure procedure can be used to predict the initiation and propagation of damage in solid and notched composite laminates. More accurate numerical results might be obtained if more refined models are used and the effects of the residual thermal stresses [16] and the material nonlinear behavior [I 71in which the stiffness matrix, [D], can be a function of the state variables (i.e. stress, strain, etc.) are considered. 5. CONCLUDING
120
” , ee 1 g so 3 so i
40 20 to
20
20
40
66
60
70
66
60
PLY ANGLE (degree)
Two
different approaches are used for the failure of the laminated composites. The iterative mixed field method is used in solving the linear part and the modified Newton-Raphson method is used in
.:NUMERICAL, 0; EXPERIYENTAL 1 ,I ,I,
-140
e p
0
REMARKS
analysis
I,,,,,,,‘,
Es[IlCHE5 ANGELPIED, IM
160
Fig. 5. The predicted initial failure stresses of angle-ply laminate specimens (@I,/- (3,)s with a through hole.
NOTCHED ANGELPIED, itr(l I
.~~$
Fig 4. A comparison of the predicted and measured angleply laminates (k@), failure loads in notched (with hole) graphit~~xy composites.
Fig. 6. The maximum ?*, distributions of angle-ply laminate specimen (6 / - 6)s with a through hole.
Failure analysis of laminated composites
solving the nonlinear part for modeling the initiation and progression of damage. All models are seen to agree reasonably well with experimental results for solid and notched tensile specimens. Many laminated composites will continue to carry load long after the matrix or fiber has experienced many local failures. The failure analysis used herein underestimates the load carrying capability in angle-plied laminates with holes. More refined models around stress concentrators might be needed. The effect of layer thickness on the initial failure stresses and modes could also be identified.
REFERENCES
R. M. Jones. Mechanics of Composite Materials. Scripta Book Co., Washington (1975). B. D. Agarwal and L. J. Broutman. Analysis and Performance of Fiber Comttosite. John Wilev. _I New York (1980). ” . J. D. Lee. Three-dimensional finite element analysis of damage accumulation in composite laminate. Comput. Struct. 15, 335-350 (1982). W. C. Hwang and C. T. Sun. A finite element iterative approach for analysis of laminated composite structural elements. Conference on Finite Methods in Engineering, Melbourne, Australia, 19-21 August 1987, pp. 1623.
CA.9 31,1--D
47
K. J. Bathe. Finite Element Procedures in Engineering Analysis. Prentice-Hall, Englewood Cliffs, NJ (1982). G. Dhatt and G. Touzot. The Finite Element Dtkplayed. John Wiley, New York (1984). S. W. Tsai. A survey of macroscopic failure criteria for composite materials. J. Reinforced Plast. Compos. 3, 40-62 (1984).
8. R. B. Pipes and B. W. Cole. On the off-axis strength test for anisotropic materials. J. Compos. Mater. 7, 245-256 (1974). 9. R. Y. Kim. On the off-axis and angle-ply strength of composites. Test methods and design allowables for fibrouscomposites. ASTM STP 734, pp. 281-337 (1981). 10. Z. Hashin. Failure criteria for unidirectional fiber composites. J. appl. Mech. 47, 329-334 (1980). 11. F. K. Chang and G. S. Springer. The strengths of fiber reinforced composite bends. J. Compos. Mater. 20, 3@45 (1986). 12. 0. C. Zienkiewicz. The Finite Element Method, 3rd Edn. McGraw-Hill, New York (1977). 13. W. C. Hwang. Three-dimensional finite element based, iterative approach for analysis of laminated composites. Ph.D. Dissertation, Department of Engineering sciences, University of Florida. (1986). 14. S. W. Tsai. Strength characteristics of composite materials. NASA CR-224 (1965). 15. T. B. Irvine and C. A. Ginty. Progressive fracture of fiber composites. NASA TM-83701 (1984). 16. H. T. Hahn and N. J. Pagano. Curing stresses in composite laminates. J. Compos. Mater. 9, 91-106. (1975). 17. A. Rotem and Z. Hashin. Failure modes of angle ply laminates. J. Compos. Mater. 9, 191-206 (1975).
failure problem. The increased use of high-strength, low-weight, fi~r-~infor~ composite materials, especially in aircraft and space vehicles, has heightened interest in understanding and characterizing their behavior. The trend toward higher design levels has led to the need to understand failure initiation and to predict strength. The overall effectiveness of the finite element analysis depends to a great extent on the numerical algorithms used in the solution of the algebraic system of equations resulting from the discretization process. As the finite element method is applied to larger and more complex problems, it becomes increasingly important that the solution process remains economical and accurate. We believe that finite element equations can be solved more cheaply than at present, but the failure analysis of multilayer composite laminates using thr~~mensional tinite element models often requires large computational resources [3]. The necessity to reduce computational cost or to run in minicomputers provided the incentive to formulate the effective algorithms presented in this research with the iterative mixed field method, which has already been described in [4]. Usually the sofution of these nonlinear equations is obtained by iteration [S]. In this paper, attention is focused on the failure analysis of composite laminates by using the iterative mixed field method in the linear part and a variable incremental load step scheme with the modified Newton-Raphson iteration in the nor&near part (progressive failure process). These procedures are carried
1. INTRODUCFION
The standard finite element discretization following system of algebraic equations:
yields the
in which (V) is the dispI~ment vector and considered as the basic unknown, (R} is the vector of applied “loads” and [K] is the assembled “stiffness” matrix. The interpretation of these quantities varies according to the physical problem under consideration. For structural appli~tions, the terms “load” and “stiffness” are directly applicable. In dynamic analysis, the vector {R) would include the inertia and damping forces. Hence, a dynamic analysis is basically a static analysis including inertia effects. We assumed that disp~a~ments of the finite element assemblage are small and that the material is linearly elastic. In addition, we also assumed that the nature of the boundary conditions remains unchanged during the application of the loads. Equation (1) corresponds to a linear analysis of a structural problem because {V) is a linear function of(R). For high-strength fiber-reinforced laminated composite, the fiber and matrix generally behave elastically until they have localized structural failures [I, 21. At this time the global structural stiffness should change; the response is only linear prior to localized failures, The nonlinearity comes from partial failures within a laminate. This situation arises in the analysis of a 41
W. C. HWANG and C. T. SUN
42
out on the numerical failure analysis of a T300/934 graphite-epoxy (GE) composite, using three-dimensional finite element models, and the results are compared with experimental results obtained at NASA Lewis Research Center.
2. MODIFIED NEWTON-RAPHSON METHOD fN COMPOSITE MATERIALS
The most popularly used iterative schemes for the solution of nonlinear finite element equations are forms of Newton-Raphson iteration [5]. This is often used in conjunction with an incremental loading procedure. The discretization of nonlinear problems can be written in the form (Y(U)} = {RI -WWl{U
-09
tb)
(2)
where we define the residual {r(u)} as the vector of the unbalanced nodal forces at any given deformed state (say IV)). The nonlinear soiution process consists of searching for a vector {CJ} that renders the residual as small as possible. The deformed state {V) is an exact solution of the discretized equations if, and only if, the unbalanced nodal forces vanish. Since the governing eqns (2) are nonlinear and generally of high order, they cannot be solved explicitly. An iterative algorithm, such as Newto~Raphson, can be used to solve the equation numerically. It is illustrated graphically in Fig. l(a) for problems with one variable u [6]. The instantaneous stiffness matrix [K,] is represented by the tangent to the curve at each step. For any known {U’), finite element modeling and analysis provides the numerical value of {y(V)} and [K,(V)] of the residual with respect to {vi) at the ith load step. Since evaluation of the stiffness matrix requires many calculations, it may be desirable to evaluate it less frequently. Figure l(b) shows the result of the Newton-Raphson algorithm with this modification. The choice of iteration steps, when the stiffness matrix should be updated, depends on the degree of nonlinearity in the system response, i.e. the more nonlinear the response, the more often the updating should be done. By using stress transformation equations we can determine the stress along the principal material axes in each lamina. Once the stress field is known, the application of a lamina failure criterion is appropriate at each “point” for initial failure analysis and at each “element” (averaged stress in each element if finite element method is used) to predict ultimate strength. In many cases, laminates have considembIe strength remaining after initial failure. If there are no failures in a load step, then the applied loading is increased until failures occur. Otherwise, the material is said to fail by the failure m~hanism. When a failure takes place, a change in the stiffness matrix, [AD], due to the localized structural failure is calculated. Using the change in
Pf PIR
fc) _~. 0
Ul
“2
U
U3U4U”
Fig. 1. Graphical iilustration of (a) Newton-Raphson, (b) modified Newton-Raphson, (c) iterative mixed field algorithms for a one variable u.
stiffness and the known strains, an imbalance vector {AQ} is found.
Ioad
(3) where n is the number of stiffness degradations during the iteration. The steps used in the modified Newton-Raphson iteration are for i = 1,2,3, . . . (for each iteration) [6]. (1) Solve the increased displacements
{AU}i
[K](ALr}j = (AQ}j. (2) Compute the new displacements {V)‘= {U)‘-‘+fAU)‘.
(4) (U)j (5)
(3) Execute convergence test. This test is the same as in the iterative mixed field method used in [4]. If convergence is not reached, this gives the required change in strains and stresses. (4) Extract nodal displacements {V}’ from {U}‘. (5) Determine the strains and the redistributed stresses. (c)‘= [e](uqi
[II*]’ = [b] - [AD]
(6) (7)
43
Failure analysis of laminated composites
(u)‘= p*l’((c)i-
(Lo])
{a}:=[rl{~}i,
(8) ments, since it permits identi~~tion (9)
where [b*]’ becomes a ~st-fails rotated elasti~ty matrix. The stresses, (a)‘, are changed by the degradation of matrix [D *I’. (6) Employ the failure criterion to check whether there are any failures or not. If there are no new failures, then add the incremental load. If the convergence rate becomes slower, a new global stiffness matrix, [K], is reassembled. For each iteration, the above-mentioned procedure is repeated until a converged state is reached or the laminate experiences gross failure. If we keep the finite element model and the structural stiffness unchanged, this iterative procedure is very similar to the iterative mixed field method, which is used in solving the linear part of the failure analysis. This similarity makes the programming implemented in the nonlinear analysis easier. The illustration of the iterative mixed field method with one variable u is shown in Fig. l(c).
3.
LAMINATE
FAILURE CRITERION
Failure criteria have been in use for structural materials for centuries. They can be divided into non-interactive criteria, such as the maximum stress or maximum strain; and interactive criteria, such as quadratic approximation. A common approach to formation of a failure criterion is to assume a quadratic polynomial in stress[7]. It is the simplest choice which can adequately describe the experimental data to consider the combined effects of stresses on failure. The quadratic failure criterion for composite materials has been proposed by many workers for at least 20 yr. A general quadratic, including linear terms, is the Tsai-Wu failure criterion [l). In this criterion, faihtre is assumed to occur when the following equations are satisfied: F,a, + Flja,al = 1 (i,j = 1,2, . . .6),
(10)
where F, and Fu are lamina strength tensors of the second and fourth order, respectively. The usual contracted stress notation is used except that Us = rz3, a, = tft, and a, = tn. The linear stress terms in (10) can be used to account for the difference between tensile and compressive strengths. Although it seems physically attractive [ 11,some components in the tensor polynomial require difficult and expensive biaxial testing for their determination. A major shortcoming of this quadratic failure criterion is that it predicts the imminence of failure but says nothing about the different failure modes and does not specify how the composite fails. This is of fundamen~l importance for design considerations of composite structural ele-
of the nature of the failure and finite element analysis of progressing structural failure. These problems can be avoided if the different failure modes of the fiber composite and the primary stresses contributing to them are identified and each mode is modeled separately by a quadratic. The principal failure modes are (1) fiber mode: tensile fiber mode is described by fiber breakage and compressive fiber mode is described by fiber buckling; (2) matrix mode: tensile matrix mode is described by plane failure surface parallel to fibers with cr + g3 > 0 and compressive matrix mode with a, + or < 0; and (3) delamination mode: produced by the interlaminar stresses. Experimental evidence for some of these modes obtained with off-axis specimens of various composites has been described in [S] and [9]. Each of these may be expected to have a different failure criterion; and further, these criteria may be of different form for tensile and compressive stresses. Failure criteria based on this approach have been presented by Hashin [IO], in which the general three-dimensional stress is considered. Each of the failure modes described can be modeIed separately by a quadratic polynomial containing the appropriate stress components. The advantages of this approach for the biaxial test results are not needed and the failure mode is identified. The following criteria make use of the following material strength parameters, which are fiber direction (longitudinal) tensile strength fiber direction compressive strength transverse tensile strength transverse compressive strength shear strength in the 1-2 plane shear strength in the 1-3 plane shear strength in the 2-3 plane interlaminar (through-the-thickness) tensile strength interlaminar shear strength. 2, and S, are assumed as the tensile strength and pure shear strength of the matrix, respectively. The faiiure criteria are separated into three different modes. For fiber breakage, the criteria are
at failure for 6, > 0, and
0
2
01
x,
(12) =I
W. C. HWANGand C. T. SUN
44
at failure for u, c 0. For matrix cracking, the criteria are
+(zJ+($r=
1 (13)
at failure for cr + c3 > 0, and
$[($-
Il(o,+o,)+&(a,+~,)’
at failure for g2 + 6) < 0. Delaminations are predicted utilizing a simple quadratic interactive formula [3]. The delamination occurs if either a,zZ,
or
(tf,+r
)-2s :3:_
I’
Yet another form used to identify delamination follows [ 111:
0
7:3+7:3>
1
s2+7/. I
(1%
is as
the response of composite laminated structural elements up to collapse. The reason is that failure criteria are not valid after FPF [7], strictly speaking. Until this impasse is overcome, the design of laminated structural elements with no high stress concentrator is likely to be based on the onset of initial failure in an individual lamina of the laminate, even though such a criterion may be very conservative from the strength standpoint [l]. The case of biaxial loading is a, in the x direction and ka, with a fixed load factor k in the y direction. Based on the formulae of the failure criteria, the failure mode, the prediction of the externally applied loads which will initiate failure at an opening, and the location of the failure can be obtained. The progressive failure occurs when the applied loading is increased beyond the initial failure load. The treatment of the case of bending problems is similar. The mode of failure determines how the stiffness matrix will be modified to reflect the effects of the predicted damage. The way in which the lamina properties are changed depends upon the type of ply damage. A simple, post-failure reduced stiffness approach is used here. From laminate constitutive equations, if the damage includes fiber failure in an element, D,,, Di2, D13,D,, and D, are set equal to zero for that element in this study. A change in the stiffness matrix of that element occurs.
(16)
I
In the prediction of failure, (15) and ( 16) are evaluated at the critical locations in the ply interfaces of a laminate only. It is very useful for progressive damage analysis to use a finite element method to analyze a fiber composite structure [lo]. Such a computational procedure starts with a finite element stress analysis. Certain elements will fail when their stresses satisfy the failure criterion. But in order to continue the analysis, it is essential to know how the elements have failed by fiber rupture, matrix crack, or delamination from an adjacent lamina. With this information, it is possible to set up a new finite element configuration with appropriate stiffness changes or new free surfaces, and the analysis can proceed to identify the new failures.
4,
42
D,3
0
0
0
0
0
0
0
0
0
0
0
0
4s
0
[AD] = Symmetry
D,
1
1 .
(17)
If matrix failure is the only failure mode, it is assumed that the fiber can still carry load; thus, only the corresponding transverse and shear properties of the stiffness matrix are removed. The stiffness D,* = D,, = D2, = Du = D, = 0, or
[AD] =
4. INITIAL FAILURE
ANALYSIS AND THE PREDICTION OF STRENGTH
For a laminate which is loaded by in-plane forces and/or bending moments, the stresses in the layers can be obtained. If there is no bending and in-plane loads are uniform distributions on the edges, then this is a special case where the stresses in each layer are constant and planar. The layer for which the criteria are satisfied by the lowest external load will define the layer that fails initially and the nature of the failure (failure mode). This is called first-ply failure (FPF). Despite continuing research effort, satisfactory analytical techniques are not yet available for predicting
If damage includes fiber and matrix modes, then the entire stiffness matrix of that element is removed, or [AD] = [D], The components of [AD] are calculated by use of where A6, terms are substituted for D4 terms. Also, damage in one element will not affect the stiffness of other elements. If delamination occurs, we remodel a new finite element configuration with a new free surface at the ply interface. The area around the delamination fronts can be modeled by singularity crack tip elements [ 121to reduce the need for extremely fine meshes. Since the delamination mode is checked on each element basis, the average stress as calculated
Failure analysis of laminated composites
with the singular element differs little from that calculated with a conventional element. However, this mode of delamination can be easily eliminated by using the techniques which are discussed in [13]. The laminates analyzed were 50.8mm wide by 457.2 mm long specimens, which were tabbed at each end with 50.8 mm wide beveled aluminum tabs and loaded in uniaxial tension. The following types were studied: three-layered cross-ply (O~~/O) E-gIassepoxy solid specimens, hi = h3= 0.127 mm, h, = 1.27 mm; four-ply panels with orientation ( f 0), Fiberite 1034E prepreq (934 resin matrix with thornel 300 graphite fibers) solid specimensand specimensnotched with centered 6.35 mm diameter through holes. The nominal thickness of one ply is 0.127 mm. In all cases, uniaxial tension was applied at one end while the other end of the specimen was fixed. The material properties of the first simple E-glassepoxy lamina are [l] E, = 53.8 GPa (7.8 x 10”psi), E2 = Ex = 18 GPa (2.6 x 10’ psi), Gn = Gu = 8.6GPa (1.25 x 106psi), G,, = 7.2 GPa (1.04 x lo6 PSI), vu = v,j = v, = 0.25, a, = 3.5 x lO+/“F, a2= a3 = 11.4 x 10a6pF, X,=X,= 1034MPa (150 ksi), Yr = 27.6 MPa (4 ksi), Y, = 138MPa (20 ksi), S,, = S,, = S,, = S, = 41.4 MPa (6 ksi), 2, = 48.3 MPa (7 ksi), AT = -2OO“F, AM = 0.
45
Fig. 2. A comparison of the prcdii and measured angle-ply laminates ( f f?), failure loads in graphite-epoxy composites.
The thy-dimensional finite element model used for solid specimens is a 20 x 4 x 1 uniform mesh of eight-node elements in each layer. Figure 2 shows the comparison between the numerical results and experimental results [15]. Because the specimens are very thin, they can be considered as plane stress state problems. The strengths in Fig. 2 are essentially those of FPF predictions. The laminates fail to carry more load after the FPF. The dominant modes of initial failure predicted by laminate failure criteria are fiber breakage in solid laminates with ply orientation (O), to ( f 30), and matrix cracking in solid specimens with (19) ply orientations > 30”. The finite element mesh used to model the compoThe solutions are obtained by using a quadrant of site specimens with centered holes is as shown in the laminates because of symmetry. The finite element Fig. 3. Results of the numerical analyses, which are model is a 10 x 2 x 1 uniform mesh of eight-node shown in Fig. 4, agree well with experimental results elements in each layer. Results of the initial failure 1151.The predicted strengths are higher than the initial stress, a,,, and the ultimate failure stress, efi are failure stresseswhich are induced by the stress concene,, = 23.2 MPa and cr/= 195MPa. Those results are trations around holes. The laminates continue to carry reasonably close to the experimental values where a0 load after the initial failures occur. It also indicates is about 23.4 MPa and ur is about 200 MPa [14]. The that the results are conservative in their prediction of initial failure mode, which is the tensile matrix mode, strengths. Table 1 shows that the initial delamination is predicted to occur in the 90” layer. failure mode will occur when the number of plies per The assumed material properties of the second exlayer, n, increases from two to five. Delamination, ample, T300/934, are [ 15; Irvine, personal communicaas a form of internal damage, is det~men~l to the tion, 1985; Malvem, personal ~mmunication, 19851 laminate strength. This mode of failure can be easily eliminated when it is predicted numerically. The E, = 138GPa (20 x lo6 psi), predicted initial failure stresses, a,, are listed in Fig. 5. Et = Es = 10.3 GPa (1.5 x lO”psi), The values of a, will decrease as n increases from two G,* = G,, = 6.6 GPa (0.95 x 10'psi), to six. Figure 6 indicates the maximum &,, which is G2,= 2.6 GPa (0.37 x 106psi), seen in the (-+45), laminate. This is the reason why v,z = v,j = 0.21, v*3= 0.47, Xr = 1393MPa (202 ksi), Xc = 1447.9MPa (210 ksi), Yr = 44.8 MPa (6.5 ksi), Y, = 172.4MPa (25 ksi), S ,2= S13= S, = 62.1 MPa (9 ksi), S, = 34.5 MPa (5 ksi), Fig. 3. The top view of the finite element mesh of the Z, = 49 MPa (7.1 ksi). (20) angle-ply laminate specimen ( f O), with a through hole.
W. C.
46
HWANG
and C. T. SUN
Table 1. The predictions of initial delamination failure mode of (t?,,/-6,)s G/E laminates with holes Initial failure n 0 2
3
4
5
(deiree) 1.5 30 45 15 30 45 15 30 45 15 30 45 60
(gi) 89.8 35.2 31.6 69.6 27.1 24.3 57.3 22.6 20.2 50.2 19.8 17.7 22.4
Approx.? location
Approx.7 locations
(dek+e)
Maximum SCF
(de:ee)
-75 -75 -75 -75 -75 -75 -15 -15 -75 -75 -75 -75 -60
5.2 5.3 5.7 5.3 5.4 5.7 5.3 5.4 5.8 5.4 5.4 5.8 5.1
-105 -75 -60 -105 -75 -60 -105 -75 -60 -105 -75 -60 -60
Approx.? locations Maximum rd
(de&e)
-0.34 -0.95 -1.10 -0.45 - 1.23 -1.42 -0.55 -1.48 -1.71 -0.63 -1.70 -2.00 -1.40
-75 -15 -15 -75 -75 -75 -75 -75 -75 -75 -15 -75 -60
t Because of symmetry, there are always two locations on the circumference; one is located at B = tl (or &2} and another is at @= OL+ 180 (or #3,y).
the initial failure stresses are lower around the (f45), laminate. The numerical results presented herein indicate that a three-dimensional finite element progressive failure procedure can be used to predict the initiation and propagation of damage in solid and notched composite laminates. More accurate numerical results might be obtained if more refined models are used and the effects of the residual thermal stresses [16] and the material nonlinear behavior [I 71in which the stiffness matrix, [D], can be a function of the state variables (i.e. stress, strain, etc.) are considered. 5. CONCLUDING
120
” , ee 1 g so 3 so i
40 20 to
20
20
40
66
60
70
66
60
PLY ANGLE (degree)
Two
different approaches are used for the failure of the laminated composites. The iterative mixed field method is used in solving the linear part and the modified Newton-Raphson method is used in
.:NUMERICAL, 0; EXPERIYENTAL 1 ,I ,I,
-140
e p
0
REMARKS
analysis
I,,,,,,,‘,
Es[IlCHE5 ANGELPIED, IM
160
Fig. 5. The predicted initial failure stresses of angle-ply laminate specimens (@I,/- (3,)s with a through hole.
NOTCHED ANGELPIED, itr(l I
.~~$
Fig 4. A comparison of the predicted and measured angleply laminates (k@), failure loads in notched (with hole) graphit~~xy composites.
Fig. 6. The maximum ?*, distributions of angle-ply laminate specimen (6 / - 6)s with a through hole.
Failure analysis of laminated composites
solving the nonlinear part for modeling the initiation and progression of damage. All models are seen to agree reasonably well with experimental results for solid and notched tensile specimens. Many laminated composites will continue to carry load long after the matrix or fiber has experienced many local failures. The failure analysis used herein underestimates the load carrying capability in angle-plied laminates with holes. More refined models around stress concentrators might be needed. The effect of layer thickness on the initial failure stresses and modes could also be identified.
REFERENCES
R. M. Jones. Mechanics of Composite Materials. Scripta Book Co., Washington (1975). B. D. Agarwal and L. J. Broutman. Analysis and Performance of Fiber Comttosite. John Wilev. _I New York (1980). ” . J. D. Lee. Three-dimensional finite element analysis of damage accumulation in composite laminate. Comput. Struct. 15, 335-350 (1982). W. C. Hwang and C. T. Sun. A finite element iterative approach for analysis of laminated composite structural elements. Conference on Finite Methods in Engineering, Melbourne, Australia, 19-21 August 1987, pp. 1623.
CA.9 31,1--D
47
K. J. Bathe. Finite Element Procedures in Engineering Analysis. Prentice-Hall, Englewood Cliffs, NJ (1982). G. Dhatt and G. Touzot. The Finite Element Dtkplayed. John Wiley, New York (1984). S. W. Tsai. A survey of macroscopic failure criteria for composite materials. J. Reinforced Plast. Compos. 3, 40-62 (1984).
8. R. B. Pipes and B. W. Cole. On the off-axis strength test for anisotropic materials. J. Compos. Mater. 7, 245-256 (1974). 9. R. Y. Kim. On the off-axis and angle-ply strength of composites. Test methods and design allowables for fibrouscomposites. ASTM STP 734, pp. 281-337 (1981). 10. Z. Hashin. Failure criteria for unidirectional fiber composites. J. appl. Mech. 47, 329-334 (1980). 11. F. K. Chang and G. S. Springer. The strengths of fiber reinforced composite bends. J. Compos. Mater. 20, 3@45 (1986). 12. 0. C. Zienkiewicz. The Finite Element Method, 3rd Edn. McGraw-Hill, New York (1977). 13. W. C. Hwang. Three-dimensional finite element based, iterative approach for analysis of laminated composites. Ph.D. Dissertation, Department of Engineering sciences, University of Florida. (1986). 14. S. W. Tsai. Strength characteristics of composite materials. NASA CR-224 (1965). 15. T. B. Irvine and C. A. Ginty. Progressive fracture of fiber composites. NASA TM-83701 (1984). 16. H. T. Hahn and N. J. Pagano. Curing stresses in composite laminates. J. Compos. Mater. 9, 91-106. (1975). 17. A. Rotem and Z. Hashin. Failure modes of angle ply laminates. J. Compos. Mater. 9, 191-206 (1975).