t... o; .1 A. .72.. . : c: . 1.11.: 93" a r . f I In“. \ptulr 4 . F . 1.1.5.... . zn 1.9.3.1 1.2.1.. 1.; 3.3. . on!) .uu.1n1n. .r....... 1...“! .. DRUNHUNHVO. ootihru. no!!! Iv .. (filo. .) fi...’:tlclolhl.ll¢£u . . . y .. . .r . . . . . . .. . 3K . . . . "flaw“; . . , . r O , u - . v . i. F ‘ A. .r 4 . VIII-IN .r “710: - . . u v En .4. . _ . . ”£24.. .fifia- .33.». 5W)»; 1 m... .. Mm... flan... . «an Art? 2 fill. . 153.1. .. . .... g... . .. . . -. .. ya: havvufiafi . .h u... . lad-.0"; . «mi-gm... .w .. .. . . . . ,1 bk”... :1“... Ivan. . . w- a» , . .l‘m . a . .353 v.1. . $.W.ft~.nhemww% .991. . ._ vfllt... “Dumdu... .3?! a. . a N u .. .9 5.39.1. - «5:... m ,. m . 1 - wax. 2%....3... ME... 15% u .. fiufifimsmwfi... b 3%.. 2w. . .4.” .w....w.......m...wiq...... i... 531?: , 34? a ..n-fi..«..u..~ .53....4. mmvammvfizfimwmm . “gamma”? -2»..va 5 . pfiuflflwfi! .. . 5M . . . u. «a. .u... ,N..o-.v . v; .. 1. . .mmzz . -.a€:.......... fl... 2.... 5:..wawuwwwfi»... flammam.” .. lwfiumhflwa, \. .911 ...J.....znm..m.n,.fi.i.. i . nfitdrvuhtfl‘u... Lynda... vuihuuu . u :9: can“... . ‘1‘. ulhuucfnumwh an... 5.... Prep .. .U-I-‘nv. .. \J% . .. [III . .r"~-f.vto!i\ul$0llvhplll)ht 1-11. '5'. C... i . .OIIIIII’ID . . . .uv‘\vcv\»1'95|tllul.flvhd§_3 . an a r . . . flfl‘ciJfiv-I-‘fi . v I n .3 1| rhna 1.3.3.1. : Il.ivxr.i .x.‘u..9 I I. . W .(afirufi 1.1. . ....r. 13.. 34.1,...) \ v 0 I'D. {‘s I! . 1m want. . 8531.. I «31".-- :52? vé . . . . . o 1.!" WWW .\. n3 :1: I mg" dbglme, li‘ll‘tv69< I. Ingl‘nli 1. tllovlb‘. .A.IV¢.11IIJ‘.¥‘I§J u THES!S 31293 01570 0879 This is to certify that the thesis entitled Investigation of Interlaminar Stresses in Laminated Composites using the Resin-Layer Concept and Cubic Zig-Zag Plate Theory_ presented by Venugopal Amineni has been accepted towards fulfillment of the requirements for Master's degree in Mechanics mm Major professor Date 7/l47/97 RONALD AVERILL 0-7539 MS U is an Affirmative Action/Equal Opportunity Institution LIBRARY Mlchlgan State University PLACE IN RETURN 30X to remove this checkout from your record. TO AVOID FINES return on or bdoro duo duo. DATE DUE DATE DUE DATE DUE I. MAGIC 2 1 {‘9‘ ‘ #4.— v-# , , .~: "‘_; . r , I _. MSU Is An Man-two mm Opportunity Inflation gamma INVESTIGATION OF INTERLAMINAR STRESSES IN LAMINATED COMPOSITES USING THE RESIN-LAYER CONCEPT AND CUBIC ZIG-ZAG PLATE THEORY By Venugopal Amineni A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Department of Materials Science and Mechanics 1997 ABSTRACT INVESTIGATION OF INTERLAMINAR STRESSES IN LAMINATED COMPOSITES USING THE RESIN-LAYER CONCEPT AND CUBIC ZIG-ZAG PLATE THEORY By Venugopal Amineni The first section of this thesis describes the utilization of a finite element model based on cubic zig-zag sublarninate theory to analyze free edge effect in both angle ply and cross ply laminates subjected to uniform extension. Results are found to be in good agreement ' with other models in the literature. The model proved computationally more efficient com- pared to the other three-dimensional continuum based models. Resin-rich zones which exist at interlaminar boundaries of plies provide paths for crack to initiate and to propagate because of their low strength. So, it is necessary to predict the stresses and strains in resin regions accurately and efficiently. The present model is able to predict the response of simply supported composite beams with thin resin and delaminated layers subjected to sinusoidal loading very accurately compared to equivalent single layer theories. A new failure criterion based on stress state in the resin-rich zone is developed to pre- dict the failure loads in quasi—isotropic laminates under off-axis loading. Predictions of the cubic zig-zag model are found to be in good agreement with the experimental results. It is also observed that a small shift in loading direction relative to fiber direction cause a dis- tinct change in the strength of the laminates. to my beloved, Parents and Brother iii ACKNOWLEDGEMENTS Firstly, I would like to express my deep gratitude and appreciation to my advisor, Dr. Ronald C. Averill, for his constant support and guidance during this research. Without his encouragement and help this work would not have been possible. The guidance and comments of the member of my thesis advisory committee, Dr. Dashin Liu and Dr. Patrick Kwon are greatly appreciated. I must also thank the Computational Mechanics Research Group for their support and helpful discussions during the course of this work. Special thanks to my colleagues Yip, Aitharaju, Cho and Kudapa for their friendship and helpful discussions. I would like to thank my parents and brother, without their support and patience, nothing would have become reality. iv Table of Contents List of Figures ................................................................................................................... vii List of Tables ..................................................................................................................... xi CHAPTER 1 Introduction ......................................................................................................................... 1 1.1 Introduction ................................................................................................. l 1.2 Literature Review on Laminated Plate Theories ......................................... 3 1.3 Interlaminar Stresses ................................................................................... 6 1.4 Failure in Laminated Composites ............................................................. 10 1.5 Delamination and Growth in Composite Laminates ................................. 11 1.5.1 Fracture Mechanics Concept ......................................................... 13 1.5.2 Strain Energy Release Rate ........................................................... 15 1.5.4 Crack Closure Technique for Finite Element ................................ 17 1.6 Present Work ............................................................................................. 20 CHAPTER 2 Review of Cubic Zig-Zag Plate Theory and Finite Element Model ................................. 21 2. 1 Introduction ............................................................................................... 21 2.2 Theory and Formulation ........................................................................... 22 2.3 Finite Element Formulation ...................................................................... 29 CHAPTER 3 Application of Cubic Zig-Zag Plate Theory to Predict Free Edge Stresses in Composite Laminates ........................................................................................... 32 3. 1 Introduction ............................................................................................... 32 3.2 Interlaminar Stresses ................................................................................. 33 3.3 Problem Definition and Boundary Conditions .......................................... 40 3.3.1 Angle-Ply Laminate ...................................................................... 43 3.3.2 Cross-Ply Laminate ....................................................................... 44 3.3.3 Quasi-Isotropic Laminate .............................................................. 44 3.4 Results and Discussions ............................................................................ 45 CHAPTER 4 v fir Resin Rich Region Concept .............................................................................................. 58 4. 1 Introduction ............................................................................................... 58 4.2 Concept of Resin Rich Region .................................................................. 59 4.3 Results and Discussions ............................................................................ 63 4.3.1 Case I ............................................................................................ 64 4.3.2 Case H ........................................................................................... 67 4.3.3 Case III .......................................................................................... 69 4.3.4 Case IV .................. 71 4.4 Conclusions ............................................................................................... 72 CHAPTER 5 Failure Analysis of Laminate Composite with Free Edge Effects .................................... 92 5. 1 Introduction ............................................................................................... 92 5.2 Description of the Problem ....................................................................... 93 5.2.1 Finite Element Mesh .................................................................... 96 5.3 Interlaminar Stresses ................................................................................. 98 5.4 Failure Criterion for Laminates ................................................................ 99 5.5 Result and Discussions ........................................................................... 101 5.6 Conclusion .............................................................................................. 104 CHAPTER 6 Summary and Conclusions ............................................................................................. 111 6.1 Recommendation for Future Work ......................................................... 1 13 List of References ........................................................................................................... 115 vi Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Figure 8 Figure 9 Figure 10 Figure 1 1 Figure 12 Figure 13 List of Figures Energy Release Rate versus Crack Length Curve for Unstable Crack Growth. ...................................................................................................... 16 Energy Release Rate versus Crack Length for Stable Crack Growth... ..... 17 Schematic of Crack Closure Technique. ....................................... 18 Finite Element Representation for the Crack Closure Technique (a) before extension and (b) after extension ............................................... l9 Schematic Diagram of Division of Sublaminates. ..................................... 23 Topology of Element .................................................................................. 30 Geometry of the Problem. .......................................................................... 36 Schematic Diagram of the Interlaminar Stresses in An Angle Ply Laminate. ................................................................................................... 38 Geometry and Boundary Conditions. ........................................................ 4O Schematic Representation of Finite Element Mesh for (45/-45)s Laminate. ................................................................................................... 43 Schematic Diagram of Finite Element Mesh for (O/90)s Laminate. .......... 44 Schematic Representation of Finite Element Mesh for (90/O/-45/45)s Laminate. ................................................................................................... 44 Axial Displacement Across Top Surface (z=2h) of (45/—45)s Laminate. ................................................................................................... 47 vii Figure 14 Figure 15 Figure 16 Figure 17 Figure 18 Figure 19 Figure 20 Figure 21 Figure 22 Figure 23 Figure 24 Figure 25 Figure 26 Figure 27 Distribution of 0 Along Center of Top Layer (z=3/2h) of (45l—45)s Laminate. ................................................................................................... 48 Distribution of 1x), Along the Center of Top Layer (z=3/2h) of (45l-45)s Laminate. .................................................................................... 49 Distribution of sz Along (45l-45) Interface (z=h) of (45/-45)s Laminate. ................................................................................................... 50 Distribution of tyz at [4S/-45] Interface (z=h) of (45/-45)s Laminate. ................................................................................................... 51 Interlaminar Shear Stress Distribution Through the Laminate Thickness of (45-45)s Laminate. ............................................................... 52 Axial Displacement Distribution Through the Laminate Thickness of (45/-45)s Laminate. ............................................................................... 53 Distribution of tyz Along [0/90] Interface (z=h) of (O/90)s Laminate. ................................................................................................... 54 Transverse Displacement Across Top Surface (z=2h) of (0/90)s Laminate. ................................................................................................... 55 Shear Stresses Distributions at Different Interfaces for (90/0/—45/45)s Laminate ............................................................................. 56 Through-Thickness Distribution of 1: x2 for (90/Ol-45/45)s Laminate. ................................................................................................... 57 Micrograph of AS4/3501-6 Carbon Fiber Epoxy Laminate Showing Resin-Rich Zone between the Plies [55] .................................................... 59 Schematic Diagram of Geometry with Loading and Boundary Conditions. ................................................................................................. 65 Effect of Resin Layer Thickness on Transverse Deflection (w) at Middle of Simply Supported Beam. .......................................................... 74 Variation of Shear Strain sz With Resin Thickness at Simply Supported End of Beam. ............................................................................................. 75 viii Figure 28 Figure 29 Figure 30 Figure 31 Figure 32 Figure 33 Figure 34 Figure 35 Figure 36 Figure 37 Figure 38 Figure 39 Figure 40 Through the Thickness Distribution of In-plane Displacement vs. Normalized Thickness, Layupl. ................................................................ 76 Through the Thickness Distribution of 0' H vs. Normalized Thickness, Layupl. ...................................................................................................... 77 Through the Thickness Distribution of Shear Strain 7“ vs. Normalized Thickness, Layupl. ................................................................ 78 Through the Thickness Distribution of Normal Strain ezz vs. Normalized Thickness, Layupl. ................................................................ 79 Through the Thickness Distribution of In-plane Displacement vs. Normalized Thickness, Layup2. ................................................................ 80 Through the Thickness Distribution of on vs. Normalized Thickness, Layup2. ...................................................................................................... 81 Through the Thickness Distribution of Shear Strain sz vs. Normalized Thickness, Layup2. ................................................................ 82 Through the Thickness Distribution of Normal Strain 822 vs. Normalized Thickness, Layup2. ................................................................ 83 Through the Thickness Distribution of In-plane Displacement vs. Normalized Thickness, Layup3. ................................................................ 84 Through the Thickness Distribution of on vs. Normalized Thickness, Layup3. ...................................................................................................... 85 Through the Thickness Distribution of Shear Strain sz vs. Normalized Thickness, Layup3. ................................................................ 86 8 Through the Thickness Distribution of Normal Strain 5522 vs. Normalized Thickness, Layup3. ................................................................ 87 Through the Thickness Distribution of In-plane Displacement vs. Normalized Thickness, Layup4. ................................................................ 88 ix Figure 41 Figure 42 Figure 43 Figure 44 Figure 45 Figure 46 Figure 47 Figure 48 Figure 49 Figure 50 Figure 51 Through Thickness Distribution of on vs. Normalized Thickness, Layup4. ...................................................................................................... 89 Through the Thickness Distribution of Shear Strain 'sz vs. Normalized Thickness, Layup4. ................................................................ 90 Through the Thickness Distribution of Normal Strain 822 vs. Normalized Thickness, Layup4. ................................................................ 91 Geometry and Other Details of the Problem for Tensile Loading. ............ 94 Finite Element Mesh. ................................................................................. 97 Convergence of U-Displacement at Free Edge as a Function of Number of Elements Along y-axis. ..................................................... 105 Convergence of U-Displacement at Free Edge as a Function of Number of Elements Along x-axis .......................................................... 106 Interlaminar stress I n Distribution Along the Width of (0/90/45/-45)s Laminate at (45/-45) Interface of 7.50 Off-Axis Loading Case. .......................................................................................... 107 Interlaminar stress I H Distribution Along the Width of (0/60/-60)s Laminate at (60/—60) Interface of 100 Off-Axis loading Case. ................ 108 Failure Loads of (0/90/45/-45)s Quasi-Isotropic Laminate vs. Loading Angle. ........................................................................................ 109 Failure Loads of (0/60/-60)s Quasi-Isotropic Laminate vs. Loading Angle. ........................................................................................ 110 Table 1. Table 2. Table 3. Table 4. Table 5. Table 6. Table 7. Table 8. Table 9. Table 10. Table l 1. Table 12. Table 13. Table 14. List of Tables Material Properties of and Dimensions Graphite-Epoxy Plies .................. 41 Material Properties of the Beam ................................................................ 65 Layupl .................................................................................... - ................... 67 Layup2 ....................................................................................................... 69 Layup3 ....................................................................................................... 70 Layup4 ................................................................................. . ..................... 72 Material Properties of AS4/3501-6 Graphite/Epoxy (in psi) ..................... 95 Material Properties of Epoxy Resin ........................................................... 95 Laminate Stacking Sequence and Loading Angle for (0/90/-45/45)s Laminate .................................................................................................... 96 Laminate Stacking Sequence and Loading Angle for (0l60/-60)s Laminate .................................................................................................... 96 Predicted Failure Loads (lb/in) for (0l90/45/-45)s Laminate at Different Critical Lengths ........................................................................ 101 Predicted Failure Loads (lb/in) for (0/60/-60)s Laminate at Different Critical Lengths ........................................................................................ 102 Comparison Failure Loads (lb/in) for (0/90/45/-45)s Laminate ............. 102 Comparison of Failure Loads (lb/in) for (0/60/-60)s Laminate ............... 103 xi CHAPTER 1 Introduction 1.1 Introduction The use of advanced composites is increasing in aerospace and automobile structures where light weight and high strength are essential for achieving high performance designs. Composites offer also the flexibility of integrating manufacturing with design. For the first time, designers can use elastic tailoring to achieve specific design objectives that conven- tional structural materials cannot. However, because of the heterogeneity of the composite materials through the thickness and the anisotropy in the individual layers, the design and analytical techniques developed for conventional materials and structures cannot be used for composite materials. For example, stress concentrations around a cutout in a laminated composite must account for the boundary-layer effects. The low ratio of transverse shear modulus to in-plane tensile modulus renders composite laminates more vulnerable to transverse shear deformation and the coupling effects among the in-plane loading, in- plane shear deformation and out of plane deformation make the prediction of composite behavior more complicated and induces unique failure modes. A thorough understanding of the behavior and efficient analysis tools are imperative for the use of composites to their full potential. In this regard accurate structural theories and efficient computational mod- els have been developed on the basis of various kinematic assumptions. Anisotropy and inhomogenity associated with laminated composites greatly increase the difficulties in finding analytical solutions. The finite element method provides a power- ful tool for laminates having complex geometries, arbitrary loadings and boundary condi- tions. Conventional three-dimensional (3-D) displacement based finite elements are computationally very expensive. However, the cubic zig-zag theory and its associated finite element model [20] have several important features that are ideally suited for this purpose. Firstly, it posssess only a fixed number of degrees of freedom, regardless of the number of layers in the laminate, which makes it computationally efficient compared to conventional continuum 3-D elements. Secondly, it is accurate and satisfies conditions such as transverse shear stress continuity through the thickness and homogeneous shear traction conditions. In composite laminate analysis where local effects are important such as free edge effects or delamination, it is necessary to use theories which are computationally viable and based on the 3-D kinematics to capture complex three dimensional states of stress. Interlaminar stress is the driving force behind edge delamination. Therefore, composite laminates under uniform extension have been the subject of many investigations because they are known to develop these interlaminar stresses at free edges. Approximate analyti- cal solutions have indicated that complex stress states with rapid change of gradients occur along the edges of composites. And experimental results have shown that when the failure occurs due to interlaminar stress, the failure loads are low compared to other modes of failure. 1.2 Literature Review on Laminated Plate Theories Many plate theories have been developed for accurate analysis of laminate composites based on various kinematic assumptions. These theories can be classified based on the assumption of their displacement fields into three categories: equivalent single layer theo- ries, discrete layer-wise theories and zig-zag theories. The first theory used in the analysis of laminated composites is the classical laminate theory (CLPT) based on Kirchhoff-Love classical plate theory (CPD [1]. In this theory, the in-plane displacements are assumed to vary linearly through the thickness, while the transverse displacement is constant through the thickness. When the aspect ratio of the structure becomes smaller, the CLPT produces erroneous results [9]. This shortcoming mainly comes from the Kirchhoff ’s hypothesis that lines originally straight and normal to the reference surface remain straight and normal to the reference surface during the defor- mation, which does not allow transverse shear deformation. In recent years, many investi- gations have been focused on the development of new and refined laminate theories to overcome the shortcoming of CPLT and to improve the prediction of the behaviors of lam- inated composites with various types of geometry and loading conditions. To improve the deficiency of CLPT, the first-order shear deformation theory (FSDT) [4] was developed which is also called Mindlin plate theory. The first order shear deforma- tion theory yields a constant value of transverse shear strain through the thickness of the plate, and thus requires a shear correction factor. This theory yields good response of lam- inate behavior, provided that the material properties of adjacent layers do not vary drasti- cally. The first order shear deformation theory was improved upon by so called higher order shear deformation theories (HSDT) [5-8]. In these theories, displacements are 4 higher order polynomial form and Cl continuous through the thickness. All the theories (CPLT, FSDT and HSDT) are termed as equivalent single layer theories, differ only in their choice of displacement fields. Generally, these equivalent single layer approaches give good results in global responses, such as deflection, vibration frequency, critical buckling load, etc. However these theories have a common drawback. They are unable to account for the severe discontinuity in the transverse shear strains that occur at the inter- faces between the two adjacent layers with drastically vary material properties. Therefore, it eliminates the satisfaction of continuity of transverse shear stresses at the interfaces of adjacent layers. In view of the drawbacks of discontinuous transverse shear stresses at the lamina inter- faces, discrete layerwise theories have been proposed by Reddy [10], Toledano and Murakami [11] and Lu and Liu [12-13]. These theories assumes a unique displacement field in each layer, and enforce the interlaminar continuity of displacements and possibil- ity of achieving continuous transverse shear stresses at all layer interfaces. Reddy’s theory treats each layer within the multi-layered laminate as separate individual layers or com- bines several layers into sublaminates. He assumes that the displacement field in each layer is expanded through the thickness direction using Lagrangian shape function employed in the conventional finite element method. Several layers can combined into a sublaminate for computational efficiency. Toledano and Murakami [11] assume indepen- dent displacement and stress fields in each layer, and used Reissner’s mixed variational principle [21] in the formulation of the model. Lu and Liu extended Reddy’s theory by adding nodal rotations in addition to nodal displacement as independent variables. The theory allows the prediction of transverse shear stresses from the constitutive equation rather than resorting to stress recovering technique which is usually achieved by using equilibrium equations. During the stress recovery process, the numerical differentiation can worsen the accuracy of the results. Most significant feature of these discrete layerwise theories compared to conventional 3-D displacement finite element model is saving com- putational time while giving the same results for comparable meshes. However, drawback of these theories is large number of degrees of freedom that is proportional to the number of layers in the laminate. Recently, new theories were developed called zig-zag theories. These zig-zag theories were extremely useful because they account for layer wise distribution of in-plane dis- placement and at the same they overcome the problem of large number of degrees of free- dom associated with discrete layer-wise laminate theories. DiSciuva [14] developed first order zig-zag theory (FZZT) where in-plane displace- ments in a laminate are assumed to be piecewise linear and continuous through the thick- ness, yet the total number of degrees of freedom is only five. Although the theory enforces interlaminar stress continuity at each interface analytically, the transverse shear stress are constant through the thickness because of the low order assumed displacement. In addi- tion, it. does not satisfy the traction conditions at the top and bottom surfaces. This draw- back was overcome by extending the theory to include higher-order terms for the in-plane displacement field [17,18-20]. An improvement to first order zig-zag theory was achieved by superimposing a piecewise linear variation of in-plane displacements on a continuous cubic function through the thickness coordinate. A new theory was developed by Averill and Yip [16,17], which combines the discrete layerwise and zig-zag theories. This theory employs the sub-laminate approach in which each sublaminate may consist of several physical layers. An accurate displacement field is assumed that accounts for discrete layer effects without increasing the number of degrees of freedom as the number of layers increases. This is accomplished by satisfying analyti- cally the continuity of transverse shear stresses at layer interfaces as well as shear traction conditions at top and bottom surfaces of each sublaminate. Since the in-plane displace- ment takes a piecewise linear variation, the theory is called a zig-zag theory. Based on this theory an eight node brick element was developed by Averill and Yip which permits verti- cal stacking through the thickness of a laminate with all the proper continuity conditions. The operative degrees of freedom are located at top and bottom surfaces of each sublami- nate to facilitate the satisfaction of continuity conditions between sublaminate. 1.3 Interlaminar Stresses The response of a multi-layer fiber-reinforced composite laminate near its geometric boundaries has been a subject of intensive investigation during the last decade. Both experimental studies and approximate analytical solutions have indicated that complex stress states with a rapid change of gradients occur along the edges of composite lami- nates. This phenomenon is considered to result from the presence and interactions of geo- metric discontinues of the composite material through the laminate thickness. The irregularity has been found to occur only within a very local region near the geometric boundaries of a composite laminate, and is therefore referred as boundary layer effect or free edge effect. The high interlaminar stresses coupled with low interlaminar strength have been found to be of critical significance in the failure of laminated composite struc- tures. It has been further observed that the boundary layer effect is three dimensional in nature and not predictable by classical lamination theory. The boundary layer effect is apparently one of the most fundamental and important problems in the mechanics and mechanical behavior of composite laminates. Boundary layer stresses have been observed to be responsible for the initiation and growth of local heterogeneous damage in the forms of interlaminar (delamination) and interlaminar fracture in composite laminates under static loading. There are two stress components of interlaminar stress, i.e., normal stress and shear stress. Many researchers have concentrated their efforts on opening mode delarnination that is dominated by interlaminar normal stress. Experimental results indicate that if free edge delarnination is caused by interlaminar normal stress, the laminate can withstand loads beyond the onset of delamination. On the other hand, the free edge failure caused by the interlaminar shear stresses often results in immediate laminate failure. The first approximation solution of interlaminar stresses was proposed by Puppo and Evensen [25]. They modeled a laminate as a set of anisotropic layers under the assumption of generalized plane stress, i.e. the normal stress is zero, and the in-plane stresses are the thickness averages of the actual values. Pipes and Pagano [26,37] are the first to address all the three components of interlaminar stresses. They derived a reduced form of the elas- ticity field equations for symmetric laminates for uniform axial strain and employed the finite difference solution technique. Few years later Pagano [37] proposed a new theory and a global local model to define the complete stress field. The theory is based on Reiss- ner’s variational principle [21]. In the global and local model they divided the laminate into two parts: i) the local region that is of interest ii) the global region, remaining part. The field equations in this model are based upon an assumed thickness distribution of dis- placement component in the global region. An apparent loss of accuracy occurs in the cal- culation of transverse shear stress components in the global model, while increasing the number of sub-laminates will obviously enhance both the accuracy of transverse stress prediction as well as the local predictions. Wang and Choi [30-31] developed analytical solution for composite laminates under uniform axial strain based on Lekhnitskii’s stress potentials and the theory of anisotropic elasticity. The formulation of the problem lead to a pair of coupled governing partial dif- ferential equations. An ei gen function expansion method is developed to obtain the homo- geneous solution for the governing partial differential equations. They also showed that boundary layer or free edge stresses in composites are generally singular in nature due to the geometric and material discontinuities. Yin [46-47] used a variational method based on Lekhnitskii’s stress function to determine the free edge interlaminar stresses in multi-layer laminates subject to combined axial extension, bending and twisting loads. Wang and Crossman [39] are one of the first researchers to employ the finite element method to study free edge effects of laminates under uniform axial strain. The analysis is formulated on the basis of constant strain, triangular finite elements. They calculate inter- laminar stresses by smearing two or more laminae as a quasi-homogenous laminate. This approach can reduce the number of elements required but one should not smear two adja- cent layers with drastic change in material properties and where stresses are of particular importance. And they concluded that 622 and In are singular at the free edge of a lami- nated composite. Whitcomb and Raju [39] studied stresses in the thick laminates. They used quasi-three dimensional finite element analysis to calculate stresses in (45/0/-45/90) ply groups. They found out the interlaminar normal stress distribution in the outer group plies were insensitive to the total thickness of laminate. But the free-edge boundary layer width is related to the thickness of the total laminate. Jones et al. [40] reported three- dimensional super-element schemes which replace several elements by a super element. This scheme can reduce substantially the computational cost for the analyses of multi- layer larninates, but a region of high stress concentration needs to be modeled by ordinary elements. A sublaminate approach was utilized to construct thick laminates by Kim et al. [41]. Moment equilibrium approach is adopted with the coupling properties of unsymmet- ric laminates. The objective of the work is to rank thick laminates rather than accurate cal- culation of interlaminar stresses. Recently many 3-D finite element models are developed to analyze free edge stresses with less number of elements to reduce the computational cost. In this category, Yang and He developed a 24-node, 64-degree of freedom, solid hexahedronal finite element to ana- lyze the interlaminar free edge stresses and delarnination of composite laminate plate. The difficulty due to the requirement of large storage space is resolved by using the precondi- tion conjugate gradient method. One of the advantages of this method is that there is no need to store the fill-in elements in the stiffness matrix decomposition process, which leads to a significant reduction in storage space required. One drawback of this model is still the computational time because of large number of degree of freedom involved with each node. Lessard et al. developed a three-dimensional 20-node quadratic brick element to simulate cross-ply composite laminate. The advantage of so called slice model is less computational time and memory. 10 1.4 Failure in Laminated Composites Interlaminar stresses at free edge play important role in the failure of composite lami- nates. In classical lamination theory, a state of plane stress is assumed and the interlaminar stress are not included. It has been shown by many researchers that a three-dimensional state of stress exists near the free edge in laminated composites. These stresses give rise to edge delamination and result in failure of the laminate. Classical laminate theory does not include interlaminar stresses. Thus, the associated failure analysis predicts, failure due only to in-plane stresses but not due to free edge interlaminar-stresses which are very high in magnitude compared to in-plane stresses. If the failure of a laminate is controlled by free edge stresses, then the prediction of laminate strength using classical laminate theory is often substantially higher than the actual strength of the laminate. Thus, it is absolutely necessary to carry out failure analysis based on the theory which has required kinematics to evaluate interlaminar stresses. Soni and Kim [48] are one of first to address the failure due to predominate interlami- nar stress In. They used the global-local laminate variational model developed by Pagano [36] to analyze a large number of laminates made of T300/1034-C graphite epoxy. Predic- tion of the onset of failure was on the basis of the average stress failure criterion of In. Experimental results conducted by them have also shown significantly lower strength than predicted by in-plane first ply failure strength. Joo et al. [49] have shown that the average interlaminar shear stress In is linearly related to the mismatch of the extension-shear coupling of the top and bottom sublami- nates. A simple failure criterion was developed based on the mismatch of extension shear coupling to predict failure stresses and strains. Herakovich et al. [60,64] pointed out the ll interlaminar shear stress In dominates the initiation of failure in angle laminates for fiber angle smaller that 370. They also demonstrated the influence of thickness on the strength of angle-ply laminate. Experimental results by Zhou [62] indicated that if free edge delamination is caused by interlaminar normal stress, the laminate can withstand loads beyond the onset of delamination. But on the other hand, the free edge failure caused by interlaminar shear stresses often results in immediate failure. Sun et al. [65] introduced a failure criterion based on all three interlaminar stresses components at the free edge. They showed that the- oretical results obtained using this approach were in good agreement with experimental results for quasi-isotropic laminates under various loading directions. Experiments con- ducted by Sun revealed that interlaminar strength is strongly dependent on the loading angle. Further, it was observed that the classical laminated plate and Tsai-Hill failure crite- rion [66] predict much higher loads for off-axis laminates compared to the experimental results. Sun performed a quasi-three dimensional stress analysis on a single cross-section of the laminates. Symmetry was assumed about xy and 1:2 planes and the stress analysis was performed only on one quarter cross-section of the plate. 1.5 Delamination and Growth in Composite Laminates Laminated composites gained importance in structural applications due to their very high strength-to-weight ratios. These applications are, however, still very limited due to the lack of a thorough understanding of the failure mechanisms involved in laminated struc- tures. Particularly the phenomenon of progressive failure in laminated structures is yet to be understood, and as a result reliable strategies for designing optimal laminated compos- 12 ite structures for desired life and strength are not yet available. It is known that laminated structures, in spite of their excellent mechanical properties are particularly prone to delarnination type of failures. Such delarnination are know to occur due to weakness in planes that are parallel to the lamination plane. This weakness usually manifests itself in the form of delarnination (or delayering) failure, a failure mode that has been and continues to be investigated thoroughly. Mode 1, mode II and mixed mode fracture tests have been performed by many researchers to observe and understand delamination failure [68-70]. Presence of these failures leads to flaws and cracks and their growth cause great concern because they may, in some cases, lead to catastrophic laminate failure. In most other case they induce long term laminate property degradation. Both are detrimental to the structural reliability and durability of the laminates. One of the modes of matrix dominated failure that involves ply structural interactions and has not been investigated as intensively as the other failure modes. Cause of delamina- tion for this type of failure have attributed generally to the existence of interlaminar stresses usually found near the free edge of laminates [71-72]. Ply delarnination is found to depend on the thickness of the ply which is stressed inter- laminarly, The characteristic of the thickness effect is that the occurrence of ply delamina- tion may be Suppressed by decreasing the ply thickness as shown by Rodini, et al. [73] on edge delamination of graphite-epoxy laminates. Rodini used a probabilistic argument that the laminate with thick plies contains statistically more defects than laminates with thin plies. Consequently, a thicker ply is likely to fail at lower stress levels. This approach requires the knowledge of the interlaminar stress distributions which may become singular and the strength distribution of material in a unit volume. Wu [74] presents a more general 13 approach based on the Weibull statistical strength theory which takes into account the effect of stress gradient near the site of stress concentration. These statistical approaches are certainly plausible, which do provide an explanation of the thickness effect. However, the observed thickness effect on the threshold stress for ply delarnination may be explained readily from an energy point of view. In this method the actual amount of strain energy trapped in the plies of the laminate and the manner of its release during a crack process play an important role in crack initiation and in its growth behavior. So, the thickness of the plies that contain the crack constitutes an important parameter in the crack propagation. Furthermore, the strain energy release mechanism is controlled at least par- tially by the structural interaction between plies during loading the laminate. Since this interaction can be altered by the kinematics of the crack, the energy argument provides not only a criterion for crack growth but also for the kinematic effects such as growth stability. 1.5.1 Fracture Mechanics Concept Delamination is a crack separating adjacent laminae, and the plane of the crack lies in the plane of the interface between laminae. In order to understand the initiation of delarni- nation, it is necessary to study the interlaminar stresses which are part of a complex three- dimensional state of stress. Such a complex state of stress at the crack tip inhibits the effective use of the stress intensity factor approach, and makes the problem ideally suited for the strain energy release rate approach. Rybicki, et al. [75], described the interlaminar delarnination problem by a crack growth analysis based on the strain energy release rate. The premise for analysis is the interlaminar flaws or defects existing near the ply interface region cause initiation of ply delarnination whenever a certain condition is reached. The 14 growth of the delarnination is stable initially, but the applied load must be increased in order to extend the crack. Critical stain energy Gc is the measure of material resistance against ply delarnination or sometimes called fracture toughness of the material which can be obtained experimentally. During the stable growth, strain energy is released when new crack surface is created under the applied load. The rate of the energy release per unit crack surface, G, can be calculated by stress analysis. The driving force can extend the crack further when G=Gc. If GGc. The strain energy release rate G, is generally a complicated function of crack location, crack geometry, ply properties and the applied loads. Rybicki employed the finite element method in conjunction with a crack closure procedure to calculate G. The main assump- tions in this method is the delarnination involves matrix dominated fracture and the crack path is parallel to ply interfaces. The numerical method involves the introduction of a crack with a known length and computing the work done to close it. This concept is an extension of Irwin’s [76] crack closure integral using the finite element method. Since the work done in closing the crack is computed directly from the finite element nodal force and nodal displacement thereby eliminating the analysis of stress field is eliminated. Hence, the advantage of the energy approach compared to stress analysis is that, for the same level of accurancy, a lower order of precision suffices for the energy calculation which is based on the product of force and displacement, as opposed to the stress calcula- tion which is based on the derivatives of displacement. Furthermore, there is no need to 15 use any singular stress element formulation in order to obtain a solution for the crack tip nodal force and nodal displacement. 1.5.2 Strain Energy Release Rate In classical fracture mechanics, it is assumed that strain energy is released when a crack surface is created in a stressed body. The rate of strain energy release during the crack extension depends on the material involved. However, the basic character of the energy release rate coupled with the fact that it can be evaluated with a simple computa- tional scheme make it a potentially useful tool for studying damage growth in composite materials. The energy release rate, denoted by G is defined for a virtual crack extension as G = — -— (1.5.3) Where W is the work of the external tractions per unit thickness, U is the strain energy of the body per unit thickness and a is the crack length. The physical significance of the energy release rate corresponds to the rate of change of energy of the system per unit area of crack length. The utility of the energy release rate as a criterion depends on the energy dissipation in the fracture process being a material property. When taken as criterion for crack growth, the energy release rate also provides a means of identifying stable crack growth which is the ability of a material configuration to sustain higher loads after some crack growth. Stable growth can be considered in terms of the shape of the energy release rate versus crack length curve. A steady increasing curve as shown in Figure 1 indicates unstable crack growth. 16 P T P = Constant curves O _ Increasing values 2. a0 of P a 0.1 U 2 i 1% Critical Value of G a; — Ta" CiaEk Extension “2 I m I | | I I I l I 1 5 > a0 a0 + Aa Crack Length, a Figure 1. Energy Release Rate versus Crack Length Curve for Unstable Crack Growth. Consider a flaw of size (10 as illustrated in Figure 1. The critical value of G for crack exten- sion is reached when the load is increase to P0. As the crack extends to a0 + Aa , the curve shows that the amount of energy available to drive the crack exceeds the critical value. Thus, the crack will continue to grow unstably. An illustration of stable crack growth is shown in Figure 2. 17 A 1* P P = Constant curves Increasing values G - r P " O 9:; an a: § ; P P“ / 'b’ Critical Value of G . :3 _______ -Fo_r' Crack Extension :3 C: LU I . . l Stable ' Unstable Regime l Regime l l n > 30 a0 + Aa Crack Length, a Figure 2. Energy Release Rate versus Crack Length for Stable Crack Growth. Again, the critical value of P for crack extension is denoted by P0. However, as the crack grows to a length of a0 + Aa, the amount of energy available to drive the crack decreases to a value below the critical value of G. Thus, the crack will not propagate until the load is increased above P0 and the crack is stable until the curve starts an upward trend. 1.5.4 Crack Closure Technique for Finite Element Ribicki and Kanninen [75] presented a very simple method of evaluating the Irwin’s [76] crack closure integral 18 An G = lim Lje-Aada (1.5.5) 0 M -> 02M where AI) is a displacement vector between the crack faces, and 8 is a vector defining stress components per unit area which are required to close the cracks as shown in Figure 3. E4— Figure 3. Schematic of Crack Closure Technique. This technique simplifies the energy release rate computations because the knowledge of the singular stress field near the crack tip is not required. Refering to Figure 4(a), the nodal forces are evaluated at nodes b and c for a crack size a. Then nodes b and c are released and crack is allowed to extend an amount Aa as shown in Figure 4(b). After the crack extends, the relative displacements between nodes b' and e' , and nodes c' and f' are measured. The energy release rate can now be expressed as the work required to close the crack by an amount of An by: 19 1.5.6 G, = 5:3[sz(wb.—we.)+Fn(wC.—wf)] ( ) 1 G" = i—A—a-[Fyb(vb' — ve') + ch(vc' — vf’)] 1 Gill = mlFxbO‘b' - ue') + Fxc(uc' '- “f” where the components G , , G ,1 and G H, are referred to as mode 1, mode 11 and mode III strain energy release rates, respectively. In the above equation, F zb is the nodal force in the z direction at node b, and wb’ is the displacement in the z-direction at the node b’. Then the total strain energy release rate can be calculated as the sum of all the possible components, i.e., f, . . +Aa Lu . (a) ‘__a (b) «A3 a Figure 4. Finite Element Representation for the Crack Closure Technique (a) before extension and (b) after extension. 20 1.6 Present Work In the present work, a finite element model is developed to analyze composite lami- nates under uniform axial strain. The model is based on the cubic zig-zag theory and the sublaminate approach developed by Averill and Yip [19]. And the same model is used to predict free edge interlaminar stresses for composite laminates subjected to uniform extension. The second part of this study involves investigation of resin-rich regions in lam- inated composites. It is observed that there often exists a finite region of resin with few if any fibers at the interface of plies of laminated composites. Because of inherent weakness of the strength of resin layers, they are prone to shear failure. So, it is important to com- pute the stresses and strains in resin layer accurately and efficiently. A new failure crite- rion is developed based on the maximum interlaminar shear stress in the resin layers to predict the failure loads in quasi-isotropic laminates under off-axis loading. CHAPTER 2 Review of Cubic Zig-Zag Plate Theory and Finite Element Model 2.1 Introduction Many theories have been developed for accurate analysis of composite structures. Most of these theories were obtained from three-dimensional elasticity theory by making assumptions regrading the variation displacement field through the thickness direction of the laminate. The first advancement in this regard was first order shear deformation theory (FSDT) [2-3], which accounted for transverse shear effects ignored by classical plate the- ory. First order shear deformation theory requires a shear correction factor when applied to cases of drastic change in the material properties in the composite laminates. Then, the higher order theories were developed [5-8] which involve higher order expansions of the displacement field which overcame the problem of shear correction factor of FSDT. The prominent common feature of the theories discussed above is an assumed dis- placement field that is continuous across the entire laminate thickness. The theories differ only in the their specific choice of the assumed displacements. This displacement assump- tion, however, usually guarantees discontinuous interlaminar stresses or tractions at inter- faces between layers of different material properties. A different class of approximate laminate theories were developed called discrete 21 22 layer theories, in which a unique displacement field is assumed within each layer. These theories allow the possibility of achieving continuous transverse shear stresses at layer interfaces. But, the major drawback is the number of degrees of freedom is proportional to the number of layer in the laminates which makes it computationally very expensive. Finally laminate theories called zig-zag theories were developed which account for all the features of discrete layer theories except they are layer independent in terms of the number of degrees of freedom. This reduction in the number of degrees of freedom over- comes many of the computational difficulties associated with stress analysis of discrete- layer theories. In this chapter, cubic zig-zag theory developed by Yip and Averill [20] will be described. This theory satisfies transverse shear stress continuity between the layers of the composite laminate and traction conditions at top and bottom surfaces of the laminate. So, it is ideally suited for the investigation of free edge stresses and free edge delamination due to interlaminar shear stresses. 2.2 Theory and Formulation Consider a laminate consisting of N perfectly bonded layers, with each layer being of independent thickness and having independent stiffness properties stacked together in the thickness direction. The laminate is modeled as M sublaminates with each sublaminate consisting of N”, layers, where m is the sublaminate number as shown in Figure 5. Mathematical representation is: where N is the total number of layers in the laminate. (2.2.1) / zN(2) /// / i / . Z- /////:/; :20) ' Sb ' / / / / 1(2) ._,._._._JilammateM___-L ///;/:// 20(2) Z2 b ___________ ,//// / / . ///// /// -—--Sublammate-2-—-——-4r/ / / / F ___________ ]’/’ //// 21+ ——————————— 3’ //// / / p--- _- ——__lr / .___S_I.lbIaImnateI__._./ Zo x Figure 5. Schematic Diagram of Division of Sublaminates. In each sublaminate, the transverse displacement field is assumed to vary quadratically through the thickness. The in-plane displacement field is assumed to be a cubic polyno- mial expansion in the local thickness coordinate and a linear piecewise function is super imposed. The displacement field for the nth layer of mth sublaminate is given by [20]: 24 3 n-l k. “ uin)(x. y. z) = 2 z uk+ 2 (z-Z,)§i k=0 i=1 3 n-l n kA A u§ ’(x.y. z) = 2 z vk+ X (z—z,)n, (2.2.2) k=0 i=1 aim = will -%l+w:(%)+rs(%XI-i) E; and fl,- are the zigzag functions and u x, uy and uz are displacements in global x, y, and 2 directions. h is the total thickness of sublaminate m, subscripts t and b refer to top and bottom surfaces of the sublaminate and n is the layer numbers. Poisson’s ratio stiffening effect is eliminated by introducing the non-conforming bub- ble function associated with displacement field variable uz which is observed in contin- uum based elements when only one layer is used through the thickness. This stiffening effect arises from the bending energy terms and is due to the inability of approximate functions to interpolate linearly varying transverse normal strain through the thickness of the plate. In the computation of transverse shear strain, W3 terms are ignored because of the fact that it is used only to eliminate Poisson’s stiffening effect. This eventually leads to W3 being a nodeless variable and which is condensed out at the element level. The constitutive relations for a three dimensional stress state for the kth layer of the mth sublaminate with respect to the laminate coordinate axes is given by: 25 - w) Cfik’c‘k’c‘“ o o C‘k’r w.) 6"" C(k) C5“ 0 0 C00 8” a” (k) (k) a” 6 C33 0 0 C36 8 a = 22 (2.2.3) 1:yz C(k) CY) O yyz T xz Sym C(Sk) 0 M fry. C221 .70.. The thirteen anisotropic material constants in Eq. (2.2.3) are related to the nine elastic material constants with respect to material symmetry axes through the well known stiff- ness transformation law [44]. The strain-displacement relations are as follows: E = xx ZZ mlcu a = to II 6: ll <2le N S _a a a a _a 3 (2°14) w v _ _w if v u 'sz =8y +37 V“ - 3x +32: 7” =8—x+ 8y Both Ei and 1”],- degrees of freedoms can be eliminated by satisfying the transverse shear stress continuity at each layer interface. The conditions for transverse shear stress continuity at the nth interface is; (n) - rm”) I“) - 10”” (2.2.5) sz _ xz yz — yz Using the strain-displacement and constitutive relations along with Eq. (2.2.2), Eq. (2.2.5) can be solved for E,- and 1"], in terms of the remaining degrees of freedom. 26 .. .. Bwb .. . 8w dwb Tli = ali(vl +3—y )+02iv2+a3iv3+a4i(ay_ 3y ) 3w Bwb) .. Bwb +b1i(ul +— )+ b21u2+b3lu3+b4,( - 3y ax ax (2.2.6) 3 (. awb) .. .. (3w, awb) I = ' +_ + . + o + . — -— I Ci. Vi 8y C21V2 C3iV3 ‘34; 3y 8y .. awb . 8w, awb +d1i(“1+3—yb)+d21“2+d31“3+d4z(a”‘3; ) where pk k-l ak[f(pk )q+ZIapq]+bkL cpq k- l bpk = bk[f(k)+: 2 dpq]+ak[vz hm) (2.27) pk Ek[f(k )+ :2 llapq]+dk ] k- l ci[f""+qk- 2 dpq]+c‘ X b m] 32]} Zk q=l (k) fik)=1 f3 = (k) (k) Zk k = 1,2,...,Nm; p =1,2,...,4 and 27 (k+ 1)C(k) a, = Ak+1(cg§+1)c<,:>-c,,,)-1 5k _ Ak1+l(Css+”C(k)-C§s+l)c(k)) c“, = Ak——1—+1(C(k+1)CC(§)—Cg§+an4k)) (2.2.8) 4:. = A]:1(C(k+l)C(k)—Cf,§+l)C(k))-l kl kl kl Ak+l_ C( + )C( + )_(C(5 + ))2 The above conditions reduce the number of degrees of freedom by 2( N m — l) for each sublaminate, leaves ten degrees of freedom per sublaminate. Further simplification of Eq. (2.2.2) is possible by satisfying the applied transverse shear traction boundary conditions at the top and bottom surfaces of the sublaminate as follows: (1) _ (N...) sz z=zo - Txb sz z=h " Txt (2.2.9) T(l)| - I IW’”) - T 2:20 yb yz z=h yt where Nm is the number of layers in the sublaminate, z=0 is bottom surface and z= h is the top surface of the sublaminate, and In, , be , I and In are applied tractions (or inter- xt ’ laminar shear stresses when more than one sublaminate is used) at the bottom and top sur- faces of the sublaminate. The final form of the displacement field is obtained by introducing the surface variables: “(1) 2 c- I 2:0 0 (N...) (2.2.10) u h ’ y 28 and xb-g; “—3; (2211) awb aw! . . 9),!) = 'a—y' 6),, = 3; The rotational degrees of freedom are introduced in the above form to allow the assumed displacement form to remain CO continuous. The constraints in Eq. (2.2.11) are enforced by a penalty formulation in the finite element model. The final form of displacement field for each sublaminate is expressed in terms of operative degrees-of-freedom {2}” by solv- ing Eq. (2.2.9) and Eq. (2.2.10) either analytically or numerically [see 20], where: {fi}’" = ub, vb, wb, 6n), Oyb, In), In), u,, v,, w,, 9n, 9 I I (2.2.12) Y" xt’ W The number of degrees of freedom is not coupled with the number of layers in each sublaminate. The displacement field in Eq. (2.2.2) for the mth sublaminate can be repre- sented in terms of through-thickness shape functions as; u‘“ - a"’(x.y)¢f;"’(z) x (I. r= 1,... 7 ’ (2.2.13) a = 1, 2 u”) = a$’(x.y)~vf;"’(z) u = Wa(X, Y)Ma(z) + W3M3(Z) In the above Eq. (2.2.13) indicial notation is used. Index 0t is used to represent the top ((1:1) and bottom ((1:2) surfaces. (1):"), ‘1’?) and M a are shape functions of the thick- ness coordinate 2 only and are known functions obtained either analytically or numerically 29 by the process mentioned earlier. This zig-zag theory has several important features that are ideally suited for composite laminate analysis. Firstly, it. possesses only a fixed number of degrees of “freedom, regardless of the number of layers in the laminate. Secondly, it is accurate and satisfies conditions such as: 1. The continuity of transverse shear stress through the thickness of the laminate. 2. The homogeneous shear traction boundary conditions. 2.3 Finite Element Formulation According to the principle of minimum potential energy, the potential energy is sta- tionary at equilibrium, i.e: 511 = 5U+5V == 0 (2.3.1) where U is the strain energy is stored in the system and V is the work done by the external forces. The strain energy stored in the system can be written as: M M 5U = Z 8W0 = 2 [ cijae,jdv""’ (i, j): 1,2,3 (2.3.2) =1 m =17”) where V'" is the volume of the mth sublaminate. The potential energy of the external forces: M M iv = 2 SW") = 2 —[ tisuidr‘m) (i, j) = 1,2,: (2.3.3) =1 m=l r(”" 30 where I, = oil-n]. and I‘M) is the boundary where surface tractions are specified. In order to reduce the continuity requirements on wb , w, , the rotational degrees-of- freedom 0n, , 9n, Byb and By, were introduced, thus eliminating all second-order deriva- tives of the transverse deflection degrees-of—freedom in the strain energy function to make the associated finite element C0 continuous. Substituting the Eq. (2.2.13) into the expression of the first variation of the total poten- tial energy and collecting the coefficients of i7; , the stiffness matrix and force vector can ° be obtained. For details see [20]. The transverse deflection wb , w, is approximated using quadratic interpolation functions, while all other degrees-of freedom are approximated by using the linear Lagrange interpolation functions. This unequal interpolation scheme is a consequence of trying to satisfy the Kirchhoff ’3 constraints in a consistent manner. a) Eight node Brick element (b) Element with degrees of freedom 7 (“t’ vt’ wt’ 9xr’ eyt’ TJuct’ Tyr) ‘ 3 2 (“1y Vb, Wb’ exb’ eyb’ Txb’ Tyb) Figure 6. Topology of Element. The finite element model uses the interdependent interpolation concept of Tessler and Hughes [42-43]. This unequal order of interpolations typically give rises to elements with nodes having different degrees of freedom at comer and mid side nodes. Tessler and 31 Hughes overcame this problem by condensing out the undesirable middle nodes through the use of four differential edge constraints, i.e the vanishing of the linear transverse shear strain. The transverse deflection degree of freedom at mid-side nodes are condensed out to achieve an eight node brick element as shown in Figure 6. This topology of the elements allows the concept of the sublaminate approach to be employed via for mesh refinement in the thickness direction. Furthermore, this eight node brick element allows coupling with conventional three dimensional elements. The concept of interdependent interpolation alleviates the shear locking problem‘so that an exact order of numerical integration may be used. The model described here has been shown to be very well-suited for predicting struc- tural response and bending stresses in very complicated laminates panels [20]. In the sequel, its ability to predict interlaminar stresses at free edges will be investigated. CHAPTER 3 Application of Cubic Zig-Zag Plate Theory to Predict Free Edge Stresses in Composite Laminates 3.1 Introduction Composite materials are being used extensively in many applications due to their excellent properties of high strength-to-weight and high stiffness-to—weight ratios. Both thin and thick composite laminates are prevalent in many aerospace and automobile appli- cations which are of considerable interest to structural analyst. So, it is very important for scientists to analyze local stresses in these structures accurately and efficiently. The response of a fiber-reinforced laminated composite near the free edges has been a subject of many intensive investigations. Both experimental studies and approximate ana- lytical solutions have indicated that complex stress states with rapid change in gradients occur especially in the interlaminar shear stresses along the edges of composite laminates. This phenomenon results from the presence and interactions of geometric discontinuities and mismatch of the composite material pr0perties through the laminate thickness. The effect is found to occur only within a limited region near the geometric boundaries of a composite laminate, and is therefore referred as to a boundary layer or free edge effect. So, these high interlaminar stresses at the edges of composite laminates coupled with low interlaminar strength have been found to be a critical factor in the free edge delamination. 32 33 3.2 Interlaminar Stresses The classical lamination plate theory (CLPT) based on Kirchhoff-Love classical plate theory (CPT) [1] provides a simple and direct way to calculate stresses and strains in lam- inated composites. However, it is mathematically not rigorous because of Kirchhoff ’s hypothesis that lines originally straight and normal to the reference surface remain straight and normal to the reference surface during deformation. This can be illustrated by a simple example as shown in [23]. Let us consider a four layer laminated composite plate sub- jected to a uniaxial tension in the x-direction (Nx=No, Ny=ny=0) where Nx, N), and N x), represents force per unit length in the respective directions. The relationship between mid- plane strains and curvatures are given by: —l 0 8 = A B N (3.2. 1) K B D M where A, B and D are extensional, coupling and bending stiffness matrices. Since the lam- inate is subjected to uniaxial tension only, the Eq. (3.2.1) reduces to: [2"] = [Al—'[N] (3.2.2) where 4 k AU = 2 [Q] [hk—hk_1] (3.2.3) k =1 a is the transformed stiffness matrix, subscript k refers to the kth layer of the laminate and h is the thickness of the kth ply. 34 For angle ply symmetric laminate (+0, -0) as show in Figure 7, the extensional stiff- ness matrix is given by: (3.2.4) [A] = 2hii§i°+i§r°i where the superscript 6 and —0 represent orientation of plies with respect to global coor- dinates as shown in Figure 7. an 512 {—216 Q11 E212 "Q16 [A] = 2" Q22 Q26 + Q22 'Q26 (3'25) _sym 66.; gym @611 (211312 0 [A] = 4h Q22 0 (3.2.6) L Q66. The inverse of [A] can be written as: 22.2 .213 0 S S [A] ‘ — —1— 'sz 0 (3.2.7) .h 3" MI I'—‘ 35 — — —2 Wheres = Q11Q22‘Q12 Substituting Eq. (3.2.7) into Eq. (3.2.2), we get: ex 1 _ Nx o = _ Q22 3y 4h T 0 0 (3.2.8) Yo I 0 - ”— sym -_— _ Qé. or: '- 0‘ QZZNx 8x 4hS o = — 3y -Q12Nx (3.2.9) 0 S Lyxxi 0 So, from the constitutive relation the following stress in +0 layers can be obtained: P n I"_ _ _ " QZZNx 5; Qii Q12 Q16 4hS 6y = Q22 Q26 -Q12NJr (32-10) fix); _sym Q66 5 .- 0 _ (‘2 a N 6’ N 12 22 12 a, = 4125 ‘— 4115’ (3.2.11) ‘ ‘ N ’ ‘ N _ Q16Q22 J.,__Q26Q12 x (32.12) I” ‘ 4hS 4hS 36 It can be seen from Eq. (3.2.11) and Eq. (3.2.12), lamina stresses 0’), and In, are nonzero even on a free edge boundary parallel to the x-axis shown in Figure 7, which violates the boundary conditions at edges of laminates. Similarly it can also be shown that classical lamination plate theory violates boundary conditions for cross-ply laminates. Free Edge Fiber orientation of an angle-ply symmetric laminate (+0l-9)s Figure 7. Geometry of the Problem. Within a laminate, adjacent layers tend to slide over each other because of the differences in their material properties. Since the layers are elastically connected through their inter- faces, transverse shear stresses are developed at the interface of each layers. The shear stresses thus produced are negligible in the region away from the laminate boundaries or edges. Therefore, the laminate analysis by classical lamination plate theory which ignores shear deformation of layers gives adequate results in the domain away from the edges of laminates. A state of plane stress actually does not exist at the free edges. Rather, a three-dimen- sional stress state exists. The behavior of interlaminar stresses near a free edge in a lami- nate can be demonstrated by using equilibrium equations. ac, 8er BIZJr _ 0 3 2 13 a. +5; +37 " ‘°° ’ BI Bo BI _xy _y _yz _ ax +ay +32 0 (3.2.14) BIn 81:2), 30': ax +5; +3.2. 0 (3.2.15) Let us consider again the uniaxially loaded composite laminate as shown in Figure 7. 0' Assuming that the stresses do not vary along the loading direction, e.g. —x = 0, the 3x interlaminar shear stress, In from Eq. (3.2.13) can be written as: Z a: 1.4:) = -ja—y‘ydz (3.2.16) t where t is the total thickness of the laminate and 2 represents distance in the thickness direction. A constant in-plane stress, I is obtained by the classical lamination plate theory in xy the interior regions of the laminate as shown in Eq. (3.2.12). Now, as we move along the y direction toward a free edge, I xy must vanish rapidly to satisfy stress-free boundary con- ditions at y = ib. Thus as y —> ib , 2%” must increase. It follows from Eq. (3.2.16) that In must increase from zero in the interior region to a very large value at the free edge as shown in Figure 8. The region where these rapid changes take place is referred to as boundary layer region as shown in Figure 8. From Eq. (3.2. 14) and Eq. (3.2.15) respectively, the other interlaminar stresses are: 1,42) = -j 3’22 (3.2.17) t '2 and Z BI o,(z) = -1 5;?de (3.2.18) t (boundary Layer region) Plane Stress 4 bm > Try 1X2 ‘ b * Figure 8. . Schematic Diagram of the Interlaminar Stresses in An Angle Ply Laminate. Free edge delamination is mainly caused by the existence of these interlaminar shear stresses at free edge. Experimental results show that the boundary layer region depends on 39 the laminate thickness [24]. A great amount of analytic works [26-31] have been reported on interlaminar stress distribution in laminates. Many approximate analytical method [30] or numerical approaches such as finite difference [26] and finite element [34-37] methods have been used for the analysis of free edge stress problems. Numerical methods provide the option of having a very fine mesh to capture the stresses near regions of stress concen- tration. However, one major difficulty in conventional finite element analysis with brick elements is that the number of elements in the mesh can be excessive for accurate analysis of multi-layer composite laminates consisting of a large number of layers with varying material properties. But the present theory and associated finite element model [20] can account for many different layers in an element and can predict stresses with a less number of elements through the thickness. Wang and Crossman [35] used a simplified method to calculate interlaminar stresses by smearing two or more laminae as a quasi-homogenous laminate. This approach can reduce the number of elements required but one should not smear two adjacent layers with drastic change in material properties and where stresses are of particular importance. Whitcomb and Raju [39] studied the stresses in thick laminates. They used quasi-three dimensional finite element analysis to calculate the stresses in (45/0/~45/90) ply groups. They found out the interlaminar normal stress distribution in the outer group plies were insensitive to the total thickness of the laminate. But the free-edge boundary layer width is related to the thickness of the total laminate. Jones et al. [40] reported on a three-dimensional super-element scheme which replaces several elements by a super element. This scheme was used to reduce substantially the computational cost for the analysis of multi-layer laminates, but a region of high stress 40 concentration needs to be modeled explicitly by ordinary elements. The objective of this chapter is to show that the cubic zig-zag plate theory [20] and associated finite element model can predict free edges stresses more accurate with fewer elements than conventional finite elements in laminated composites.The theory considers all the six stress components as non-zero in general and satisfies interlaminar shear stress and displacement continuity conditions at interfaces between adjacent layers. 3.3 Problem Definition and Boundary Conditions A graphite-epoxy laminated rectangular plate subjected to uniform axial strain 60 in the x- direction is considered. The length and width of the laminate are 21 and 2b, respectively with 21=16h and 2b=8h, and ply thickness h is equal to 0.25mm. Figure 9 illustrates the I dimension and the stress boundary conditions associated with the model: '__I_y T1221: Figure 9. Geometry and Boundary Conditions. (1) On the top and bottom surfaces, 2 = it 41 In=In=oz=O (2) On the two edge-lateral surfaces, y=0 and y=2b Iyz = I xy = 0y = 0 (3) On the two end-lateral surfaces, x=0 and x=21. I =I =0, =O,u| =8021 qu=0 x=21 To demonstrate the solution scheme and the fundamental nature of boundary layer stresses, angle-ply, cross-ply and qausi-isotropic laminates of graphite-epoxy composites under uniform axial extension ex = 80 are studied. Ply elastic properties of each graphite- epoxy composite and geometric variables of the composite laminates used in present study is shown in Table 1. Table 1. Material Properties of and Dimensions Graphite-Epoxy Plies Material properties Values 13,r 20 x 106 psi Ey=Ez 2.1 x 106 psi ny=Gyz=Gn 0.85 x 106 psi 150:1),sz 0.21 h1=h2 0.25 in b=8h 4.0 in 42 These particular material and geometric constants are selected because they have been used extensively in previous studies of free edge problems [30]. Thus, solutions obtained in this study can readily be compared with existing numerical solutions available in the lit- erature. Let Dx, D), Dz represent the number of mesh division in the x, y, and 2 directions of laminate, respectively. 3.3.1 Angle-Ply Laminate A (45/-45)s graphite-epoxy laminated rectangular plate subjected to a uniform axial strain 80 in the x-direction is modeled. This laminate has symmetry with respect to the x-y plane at z=0. Therefore, only half of the laminate is modeled in this analysis. Material properties from the Table 1. are used. The uniform axial strain is applied through displacement boundary conditions. In such a case, the displacement boundary conditions are u(21, y, z) = 212°, u(0, y, z) = o, v(l, b, 0) = o and w(x, y, 0) = o. The laminate was modeled with two different meshes to study the convergence of the results. The first mesh contains Dx=10, where DJr is the number of gradually refined spaces towards the center of the laminate from edge as shown in Figure 10. Similarly Dy=15, are gradually refined elements in y-axis with the smallest element being at free edge. Only two elements in z-axis were used in both meshes because the present element has quadratic variation in shear through the thickness and has interlaminar shear stresses as degrees of freedom. The second mesh consists of Dx=10, Dy=20 and DZ=2 elements with a total 400 elements and 693 nodes. A convergence study of two the meshes has shown that the second mesh is more accurate compare to the first mesh, so the second mesh is used in all the analyses. 43 a) The front view b) Top view of actual mesh Figure 10. Schematic Representation of Finite Element Mesh for (45/-45)s Laminate. 3.3.2 Cross-Ply Laminate A [0/90]S cross-ply graphite-epoxy laminated rectangular plate subjected to an uniform axial strain 30 is considered. The material properties are the same as shown in Table 1. The displacement w at the mid-plane is assumed to zero. A finite element n 10x20x2 as shown in Figure 11 is used, where the mesh along x-axis was gradually refined toward the center of the quarter of the plate into 10 divisions, the length along y-axis was divided into 20 unequal divisions with finer mesh near the free edges and only two ele- ments in the thickness direction. Owing to the symmetry about the x, y, z axes, the analysis has been performed for the quarter of the laminate only as shown by the shade region of Figure 11. and the finite element mesh is shown in figure. The displacement boundary conditions are: u(l, y, z) = [80, u(0, y, z) = 0, v(x,0,z) = 0 and w(x, y, 0) = 0. L 2” J a) The front view b) Top view of actual mesh Figure 11. Schematic Diagram of Finite Element Mesh for (0l90)s Laminate. 3.3.3 Quasi-Isotropic Laminate A [90/0/-45/45]s graphite-epoxy laminated rectangular plate subjected to a uniform axial strain 80 as shown Figure 12.The material properties are same as shown in Table 1. The laminate consists of eight layers each of thickness h=0.25mm. Owing to symmetry about the x-y plane only half of the laminate is modeled. The laminate is modeled with a mesh of 10x20x4 elements in x, y, and 2 directions. Element in x and y direction are refined towards center and the free edge surfaces respectively. L 2” J a) The front view b) Top view of actual mesh Figure 12. Schematic Representation of Finite Element Mesh for (90/0/-45/45)s Laminate. 45 Four elements through the thickness are used to capture interlaminar stresses at each inter- face between adjacent layers as shown in Figure 12 3.4 Results and Discussions Boundary layer stresses in a composite laminate determined by the finite element model based on present cubic zig-zag plate theory are compared with existing numerical solu- tions available in the literature. I Figure 13 shows the axial displacement distribution across the width of the top surface of (45/-45)s laminate under uniform axial strain with two different meshes. Figure l4-23, distribution of in-plane and interlaminar shear stresses along the ply interface for graphite- epoxy laminate under uniform axial strain 80 and finite element solution obtained by Wang, et al., [34] are also shown for the purpose of comparison. It can be seen that the present solutions are in good agreement with Wang’s results. Figure 14 and 15 show the distribution of 0' x and I respectively along the width of xy ’ laminate at the center of the top layer. The solution given by the present model is in good agreement with the results of Wang et al., [34]. The distribution of the most dominant interlaminar stress In at the [45/-45] interface is shown in Figure 16. It can be seen that the present solution converges very closely to Wang’s result. The state of a stress singular- ity of In can be captured by using finer and finer meshes at the free edge. The present results demonstrate that the current model is able to capture the presence stress singulari- ties at the free edge. The steep increase in the gradient at the free edge as shown in Figure 16, demonstrates that In reaches a very high value and is the major contributor to mode II delarnination in laminated composites. The distribution of the interlaminar shear stress 46 through the laminate thickness at different locations from the free edge is shown in Figure 18. The shear stress vanishes at the free surface, z/ h = 2.0 and at the laminate mid plane, I/ h = 0 , while attaining maximum magnitude at the interface, z/ h = 1.0. The varia- tion of axial displacement distribution through the thickness of laminate at different loca- tions at from the free edge is shown in Figure 19. Figure 19 shows that g—Z- is very large at the interface which also suggests that In has singular behavior at the free edge of the [45/ -45] interface of the laminate. Figure 17 shows the variation of other component of inter- laminar stress In at the [45/-45] interface. Wang [34] showed that In is also singular at the free edge of the laminate. But, from the mechanics. point of view In should go to zero at the free edge because of a traction free boundary condition. The present results agrees well with Wang’s results and also satisfy the traction free boundary conditions at the free edge. Important free edge effects are also found in cross-ply laminates. Figure 20 shows the distribution of In at the [0/90] interface of a [0/90]s laminate. It is seen that In rises toward the free edge, but decrease rather suddenly to zero at y=b. This latter is, of course, required by the stress-free edge boundary conditions. Results obtained are in good agree- ment with Wang’s [34] solution and show the effectiveness of the present model.Figure 21 shows the variation of transverse displacement across the top surface of the [0/90]s lami- nate. Results of a quasi-isotropic laminate [90/0/-45/45]S under uniform tensile stretch are shown in Figure 22 and Figure 23. Figure 22 shows the in-plane shear stress In, at z = h , and the interlaminar stresses In at z = h, In at z = 3h. The distribution of In, at z=h indicates that the edge effects extend far inside the laminate than thinner four layer lami- 47 nates as shown by Whitcomb [39]. In at z=h displays a singularity at the free edge. On the other hand In at z=3h is relatively smaller in magnitude. Through the thickness distri- bution of In at the free edge is shown in Figure 23. A Significant stress singularity in In exists again at the {45/45} interface. 1.0 . 1 """ Mesh] — Mesh2 JED 0 5 w I- B 0‘00.0 L05 1.0 y/b Figure 13. Axial Displacement Across Top Surface (z=2h) of (45/-45)S Laminate. 48 3.5 7 1 T I ' f V I '—-' A.S.D Wang — Present 3.0 . a .5; 3 3f 3 3r 3 3 3 3 3 3 3 3 3 3 3 3 3 - Q. a? \O b 2 \ x e 2.5 - - b , 2 o A 1 i l a l n 1 a ' 0.0 0.2 0.4 0.6 0.8 1.0 y/b Figure 14. Distribution of 0’ Along Center of Top Layer (z=3/2h) of (457-45), Laminate. Txy/ 106 £0, psi 49 1.5 ' I ' 1 ' T ‘ r 1.0 ’ 0.5 ’ 0.0 0.0 Figure 15. Distribution of I Along the Center of Top Layer (z=3/2h) of (25/45), Laminate. 50 3.0 f I ' I ‘ I ' I , '——‘ A.S.D Wang — Present 2.0 " '6 o. 55 I co 0 v- \ 3 1.0 ~ I 0.0 ‘ ‘ l ‘ ‘ ‘ 0.0 0.2 0.4 Figure 16. Distribution of In Along (45/-45) Interface (z=h) of (45l-45)S Laminate. 51 0.40 0.20 " 0.00 Iyz/ 10680, psi I -0.20 ~0.40 °—-* A.S.D Wang -— Present 0.0 Figure 17 0.2 0.4 0.6 0.8 y/b . Distribution of In at {45/45} Interface (z=h) of (45/-45)s Laminate. 1.0 52 2.0 y=1.00b y=0.89b 1.0 0.0 0.0 1.0 2.0 1,, / 106130, psi Figure 18. Interlaminar Shear Stress Distribution Through the Laminate Thickness of (45-45)s Laminate. 53 2.0 r 1 1 y=0.0 Um”: U @ y=b, z=2h y= . 6b 1.5 r .1 y=0. 9b y=1.0b 1.0 ” ‘ 0.5 I ‘ 0.0 1 l A a 1 4 -1.0 -0.5 0.0 0.5 1.0 U/Umax Figure 19. Axial Displacement Distribution Through the Laminate Thickness of (45/-45)s Laminate. Iyz/ 1062:0’ psi 54 0.40 r 1 T 1 F °-—° A.S.D Wang —— Present 0.30 r '1 0.20 r _ 1 0.10 _ 0 00 l 1 L l L ' 0.0 0.2 0.4 0.6 0.8 Figure 20. Distribution of In Along [0/90] Interface (z=h) of (0/90)S Laminate. 1.0 55 0.8 fl -— present 1 —— A.S.D Wang 0.6 ' I I 3 0.4 I r 0.2 ” ‘ 0.0 ‘ 0.0 0.5 y/b Figure 21. Transverse Displacement Across Top Surface (z=2h) of (0/90)S Laminate. 1.0 56 2.0 0.5 " 0.0 -0.5 " 0.2 0.4 0.6 0.8 y/b Figure 22. Shear Stresses Distributions at Different Interfaces for (90/0/-45/45)s Laminate. 1.0 z\h 57 1 .00 I I . I 90 0.75 0 0 0.50 .45 , -45 0.25 45 0.00 -1.5 0.5 1,,2/106 80, psi Figure 23. Through-Thickness Distribution of In for (90/0/-45/45)S Laminate. CHAPTER 4 Resin Rich Region Concept 4.1 Introduction Laminated composites possess a number of different failure modes. These include fiber breakage, fiber buckling, matrix cracking, fiber-matrix interface debonding, interlam- inar tension failure, and interlaminar shear failure. Many researchers [50-53] concentrated their effort on open mode delarnination that is dominated by interlaminar normal stress, but very few papers dealt with the mode 11 delarnination caused by interlaminar shear stresses. Several common loadings give rise to mode H delarnination in composite struc- tures. Examples where this is the case include transverse impact of panels. Mode 11 delam- ination is characterized by relative sliding of adjacent plies after their interlaminar bond has been broken. So, this failure mode is governed by the interlaminar shear stresses in the interfacial region between the plies. Unlike metals, laminated composites with traction free edges are known to develop very high interlaminar stresses near the free edge region as illustrated in Chapter 3. Such stresses play a vital role in the failure of composite lami- nates. Moreover, in multi—directional fibrous composites, the failure process reveals inter- laminar shear stresses cause catastrophic failure [62]. To successfully predict failure in composite laminates is difficult. Since Tsai [54] first 58 59 suggested the ply-by-ply failure method to determine laminate strength, this procedure has been widely practiced. The main advantage of this procedure is its simplicity. While quite adequate for strength prediction for unidirectional composites, classical failure criteria fail to achieve the desired accuracy in most laminate strength prediction. One reason for this lack of success is the fact that lamina strength properties are not the same as the in-situ lamina properties in the laminate. Furthermore, classical lamination theory assumes that no interlaminar stresses exist between plies, which is not true. 4.2 Concept of Resin” Rich Region Resin-rich areas are often observed in composite materials even though the part has been processed according to the manufacturer’s specifications. A typical example of a resin-rich zone in composite laminate is shown in the micrograph presented in Figure]. Mllrix cracking . 5,21.- 3,51: -, 1' 1", I u .9 is if ya!“ Resin-rich zone Figure 24. Micrograph of AS4/3501-6 Carbon Fiber Epoxy Laminate Showing Resin-Rich Zone between the Plies [55]. 60 An examination of the micrograph reveals a resin-rich zone adjacent to the fibers. Resin-rich layers transfer the load between adjacent plies when load is applied to the com- posite laminate. But when the shear stresses in these regions exceed the material strength of the resin layer, mode 11 delarnination occurs and is responsible for reducing the mechanical performance of the composite laminate. Interlaminar shear failure manifests itself as a crack in the resin between the plies. This is because delarnination cracks prefer to grow along a resin path to avoid the necessity of breaking fibers and such a path is provided at interlaminar boundaries. In general, lami- nated composite structures are usually modeled as perfectly bonded layers stacked together in the laminates. The presence of a resin layer between the plies is neglected from the analysis point of view because of the fact that previous analytical models couldn’t account for the thin resin rich regions. In order to predict interlaminar failure modes, it is necessary to predict the interface shear stresses and compare these stresses with the inter- face strength property of the resin. Accurate prediction of failure requires a robust theory and associated finite element which can take into account a very thin compliant resin rich layers. In the conventional continuum based element approach, it is computationally intractable due to element aspect ratio difficulties associated with resin layers. If the resin layer thickness is increased to alleviate the element aspect ratio problems, then bending and shear prediction become sensitive to the thickness of resin layers. But cubic zig-zag theory has necessary kinematics not only to account for resin rich layers but also satisfies interply shear continuity. It has been shown that the constant normal strain assumption for thick beams wherein the layer properties vary drastically or with compliant layers doesn’t predict satisfactory 61 transverse normal stress [56]. Therefore, to improve the distribution of normal stress and its related energy in the laminate, the present theory assumes a constant normal stress through the thickness of the laminate. Normal stress, (In is calculated by using Reissner’s mixed variational principle [57], which results in following condition: N Zr 2 j (efz—ef§)dz = 0 (4.2.1) k=l ZK_l where 8:: is the strain derived from the displacement field and is given by: k _(W1_Wb) zz " h (4.2.2) 8 8:: is the transverse strain in the k’h layer evaluated from the constitutive equation assum- ing a constant variation of transverse normal stress. w, , wb are transverse displacement at top and bottom of the laminate and h is the thickness of laminate. Please refer to [56] for more detail on the evaluation of normal strain. The resulting final form of normal strain is no longer constant through the thickness but depends on in-plane and transverse stiffness and the displacement field. This procedure improves the transverse normal strain energy and gives good normal strains and average normal stress. It also eliminates the Poisson’s ratio locking in the associated finite element. Since cracks prefer to initiate and grow in thin resin rich regions of a composite lami- nate because of inherent weakness of the resin material properties, one of the following failure criteria can be used to predict the failure of laminates based on the stress state in the resin-rich regions. 1. Maximum Principal stress criterion: In maximum principal stress criterion, failure of a given material is assumed to occur if one of the principal stresses satisfies the following condition: 6’ = 6° (4.2.3) 0'2 = -00 Where do is the yield stress of the material in simple tension. 0', and 0'2 are maximum principal stress and minimum principal stress. 2. Maximum Principal Strain (or Saint-Venant) criterion: In maximum strain criterion, failure of the material is assumed to occur if the maxi- mum value of the principal strain equals the value of the yield strain in simple tension (or compression), so = co/ E . Thus if 81 is assumed to be the largest strain in absolute value, failure will occur when: 13581 = (51 —u(oz+o3) = ioo (4.2.4) Where 0’1 , (5'2 and 0'3 are the principal stresses. Go and so are yield stress and strain of material. 3. Maximum shear stress (Tresca) criterion: In the maximum shear stress criterion, failure of the material is assumed to occur if any the maximum shear stress reaches the value of the maximum shear stress. Failure will occur if one of the following condition is satisfied: where 01 , 0'2 and 0'3 are the principal stresses. The objective of the remainder of this chapter is to illustrate the effectiveness of cubic zig-zag plate theory in accommodating the resin rich layer and predicting stresses and strains in resin and delaminated resin layers of composite laminated beam under sinusoi- dal loading. The goal of this study is to determine if such an approach can eventually be used to predict mode 11 failure and model of progressive failure in laminated composites. 4.3 Results and Discussions A simply supported beam of (0/90/0)s laminate is modeled with a sinusoidally varying load of unit amplitude as shown in Figure 25. Material properties used for the laminate are given in Table 2. The aspect ratio of the beam is Llh=20. Figure 26 shows the variation of transverse deflection of the beam normalized with respect to the deflection when no resin present, and Figure 27 shows the variation of maximum shear strain in the middle layer plotted as a function of the thickness of resin layer. It can be seen from the figures that when the resin layer thickness exceeds approximately 0.01, the bending deflections and resin layer shear strains become very dependent on the thickness of the resin layer. For thicknesses below 0.0001, these results are insensitive to the resin layer thickness used. ’ Resin layer £9)— ,1 l‘ 20” Figure 25. Schematic Diagram of Geometry with Loading and Boundary Conditions. Thus, in the present study a resin layer thickness of 0.0001 (0.6% of the ply thickness) is used for each resin layer between each composite ply. Four different cases are simulated to illustrated the efficiency of cubic zig-zag theory in predicting the response of the com- posite beam under sinusoidal loading. 4.3.1 Case I In this case a simply supported (0/90/0)S beam of with resin layers at all the five ply interfaces along the beam length is modeled as shown in Figure 25. Material properties used are given in Table 2. The lamination scheme for this case is called Layupl, given in Table 3. 65 Table 2. Material Properties of the Beam Materiall Material 2 Material3 (AS4/3501-6) (in psi) (in psi) Property (in psi) (resin) (delaminated) E1 1 20.16e06 0.62e06 0.62e04 E22 1 .43e06 0.62e06 0.62e04 E33 1 .43e06 0.62e06 0.62e04 v12=v23=v31 0.3 0.34 0.34 612 = G13 0.75e06 0.23e06 0.23e04 623 0.47e06 0.23e06 0.23e04 66 Through the thickness distribution of in-plane displacement at the simply supported end of the beam is shown in Figure 28, where results presented are based on the elasticity solution of Pagano [9], the first order shear deformation theory (FSDT) [2-3], and the cubic zig-zag theory. Figures 2829 show the variation of in-plane displacement on through the thickness of the beam, wherein stresses predicted by zig-zag theory is indi- cated by circles shows good agreement with elasticity solution indicated by solid line and FSDT is denoted by long dashed lines. Through the thickness distribution of transverse shear strains are shown in Figure 30. FSDT'predicts a constant variation of shear strain and thus does not account for the presence of the resin layers.The shear strain distribution predicted by the elasticity solution has a parabolic variation in the composite plies with large discontinuity in the shear strains at the interface of plies and resin layers. The shear strain in the resin layers is represent by the horizontal lines emanating from the parabolic curve as predict by Pagano’s solution. The magnitude of the shear strains in the resin lay- ers predicted by the zigzag model is indicated by circles. Shear strain predicted in the middle layer at zlh=0.5 is 0.3944E-04 by Pagano and 0.3823E-04 by zig-zag model. It can be seen that the cubic zig-zag model represents the variation of shear strain very well in comparison with the elasticity solution. The maximum error in the shear strain at other interfaces is approximately 8%. Figure 31 shows through the thickness variation of trans- verse normal strain. At the middle of the beam the zig-zag theory matches with the elastic- ity solution. As we move towards the top and bottom of the beam, the cubic zig-zag model under-predicts and over-predicts respectively with maximum error of approximately 13%. This is due to the fact that the theory assumes a constant normal stress. 67 Table 3. Layupl Layer No. Material Thickness(in) orientation 1 1 1 t 01667 0 2 2 0.0001 - 3 1 0.1667 90 4 2 0.0001 - 5 1 0.1667 0 6 2 0.0001 - 7 1 0.1667 0 8 2 0.0001 - 9 1 0.1667 90 10 2 0.0001 - 1 l 1 0.1667 0 4.3.2 Case 11 When interlaminar shear failure is predicted to occur, delarnination failure can be rep- resented in the cubic zig-zag model by significantly reducing the stiffness of the corre- Sponding resin layers. In this case all five ply interfaces are modeled as delaminated resin layfirs along the length of the (0/90/0)S simply supported beam. Material properties and 68 lamination scheme called Layup2 are given in Table 2 and Table 4. respectively. Figure 32 shows through the thickness distribution of in-plane displacement at the simply supported end of the (0/90/0)s beam. It can be seen that the FSDT solution could not account for the very thin delaminated layers in the beam. Whereas, the cubic zig—zag models is able to predict small kinks in the delaminated resin layers and shows good agreement with the elasticity solution of Pagano. Through the thickness distribution of on is shown in Figure 33. Again cubic zig-zag model shows excellent agreement with the elasticity solution. Figure 34 shows the numerical predictions of through the thickness dis- tribution of transverse shear strains. It can be seen that the shear strain in the reduced stiff- ness resin layers has increased by two orders of magnitude over the undamaged case. Shear strain predicted by the cubic zig-zag model at z/h=0.5 is 0.3734E-02 and 0.3926E- 02 by the elasticity solution of Pagano. There is excellent agreement between the cubic Zig-zag model and the elasticity solution, whereas the FSDT doesn’t recognize the pres- ence of thin compliant delaminated resin layers. Figure 35 shows through the thickness distribution of normal strain at the middle of the composite beam. Because of the constant normal stress assumption, the zigzag model undel’Pl'CClicts and over-predicts at the top and bottom of the beam respectively. But at the Center of the middle of the beam it shows good agreement with the elasticity solution. The maximum error in the top most delaminated layer is approximately 32%. 69 Table 4. Layup2 Layer No. Material Thickness(in) orientation 1 1 0.1667 0 2 3 0.0001 - 3 1 0.1667 90 4 3 0.0001 - 5 1 0.1667 0 6 3 0.0001 - 7 1 0.1667 0 8 3 0.0001 - 9 1 0.1667 90 10 3 0.0001 - 1 1 1 0.1667 0 4.3.3 Case In In this case the (0/90/0)s simply supported beam is modeled with only the middle interface delaminated along the length of the beam. This is accomplished by reducing the Stiffness in just the middle resin layers. The lamination scheme, called Layup3 is given in Table 5, 70 Table 5. Layup3 Layer No. Material Thickness(in) orientation 1 — l = 0.1667 0 2 2 0.0001 - 3 1 0.1667 90 4 2 0.0001 - 5 1 ' 0.1667 0 6 3 0.0001 - 7 1 0.1667 0 8 2 0.0001 - 9 1 0.1667 90 10 2 0.0001 - 1 l 1 0.1667 0 Through the thickness distribution of in-plane displacement is shown in Figure 36. It can be seen that FSDT predicts linear variation of displacement without accounting for the Presence of the thin compliant delaminated resin layer at the middle of the beam. The cubic zig~zag and the elasticity solutions agree well. Figure 37 shows good agreement of on between the elasticity and zig-zag solution through the thickness at the mid-span of the beam. Figm‘e 38 shows the through thickness distribution of shear strain in the delaminated 71 resin layer. Shear strain predicted in the middle delaminated resin layer is 0.3897E-02 (Pagano) and 0.3805E-02 (Cubic zig-zag) which is approximately 2% in difference pre- dict by the models. Whereas, FSDT predicted constant shear strain. Normal strain distri- bution through the thickness of the beam is show in Figure 39. There is excellent agreement between the elasticity solution and zig-zag model in the delaminated resin layen 4.3.4 Case IV An unsymmetrical delamination failure is simulated in the (0/90/0)s simple supported beam. The lamination scheme called Layup4 is given in Table 6. Again there is excellent agreement of in-plane displacement and on through the thickness of the beam as shown in Figure 40 and Figure 41 respectively. Figure 42 shows the shear strain distribution of through the thickness in the unsym- metrically delaminated resin layers. Both cubic zig-zag model and the elasticity solutions show good agreement, whereas FSDT predicts a constant value. The maximum error in the cubic zig-zag results was approximately 18%. Normal strain also shows good agree- ment at the middle of the beam but has error of approximately 15% in the top delaminated layer compared to the elasticity solution as illustrated in Figure 43. 72 Table 6. Layup4 Layer No. Material Thickness(in) orientation = J ____=__==; l 1 0.1667 0 2 3 0.0001 - 3 1 0.1667 90 4 3 0.0001 - 5 1 0.1667 0 6 2 0.0001 - 7 1 0.1667 0 8 2 0.0001 - 9 1 0.1667 90 10 2 0.0001 - l l 1 0.1667 0 4.4 Conclusions Results for four different layups have shown that cubic zig-zag theory is accurate com- Pared to FSDT in predicting stresses and strains. Cubic zig-zag theory is capable of accommodating an extremely thin resin layer between plies without increasing the number 0f degrees of freedom and without numerical problems. This layer can be included in a 73 laminate model at minimal computational cost, providing results with no sensitivity to the resin layer thickness. Whereas, FSDT didn’t account for the presence of thin resin layers at the ply interfaces. Zig-Zag models are able to capture the discontinuous transverse shear strain and continuous shear stress at the ply interfaces of the laminate Since the cubic zig-zag model has necessary kinematics to include resin layers, it will be used to perform failure analysis due to free edge stresses in composite laminates in the next chapter. wresin\wwithout resin 74 18 v r vaval f v vvvvvvl fi 1 vavvvl V v YV‘V’fYTI 1 I 1.6 r 1.4 - 1.2 r 1.0 0.8 r 0.6 - 0.4 ~ 0.2 ~ A ALAIJ g 0.0 ‘ J .J 4 0.000001 A 0.000010 0.000100 0.001000 0.010000 lOgaresin) 0.100000 Figure 26. Effect of Resin Layer Thickness on Transverse Deflection (w) at Middle of Simply Supported Beam. 1112 75 0.400 I 0.390 0.380 I 0.370 I 4 A L L L r rrvvvv A J 0.360 Figure 27. Variation of Shear Strain 7 Supported En 0.0000 0.0000 0.0001 log(tresin) 0.0010 0.0100 0.1000 With Resin Thickness at Simply of Beam. 76 1.0 - 1 —— Pagano's 0 3 L °——-° Cubic 219-239 A ' ——- FSDT 0.6 e a 0.4 r a 0.2 — a 0.0 : 1 ‘ -0.00010 0.00000 0.00010 11 displacement Figure 28. Through the Thickness Distribution of In-plane Displacement vs. Normalized Thickness, Layup 1. 77 1.0 l I I A 0.8 ~ V - —— Pagano o——o Cubic Zig-Zag — — - FSDT 0.6 ~ - 0.4 ~ ~ 0.2 e a 0.0 v A g A l A 1 -400.0 -200.0 0.0 200.0 400.0 GXX Figure 29. Through the Thickness Distribution of on vs. Normalized Thickness, Layupl. 78 1.0 A T I r I : +—+- Pagano’s , °-—o Cubic Zig-Zag i — - — FSDT 0.8 L - Y I ' ‘ . l I l .. I - . - l T ' .. 0.6 ~ I d 0.4 - 7g, 0 g I 1 ¢ 4 C l T I l 0.2 b I -i /'/‘ i ' r = i t I ‘ i l 0.0 v 1 l l 1 0.000000 0.000010 0.000020 0.000030 0.000040 Figure 30. Through the Thickness Distribution of Shear Strain 'sz vs. Normalized Thickness, Layup 1. 79 1.0 1 T l T l’ —— Pagano’s 0.3 - °--° Cubic Zig-Zag _ 0.6 ~ , 0.4 F _ r 0.2 ~ . 0.0 4 l . 1 4 1 ~0.000010 0000006 -0.000002 0.000002 0.000006 82.2 0.000010 Figure 31. Through the Thickness Distribution of Normal Strain 822 vs. Normalized Thickness, Layupl. 80 1.0 v x —— Pagano’s °——o Cubic Zig-Zag 0.8 ~ —-- FSDT. - . 0.6 L 0.4 ~ 0.2 L 0.0 * ‘ - -0.00010 0.00000 0.00010 u displacement Figure 32. Through the Thickness Distribution of In-plane Displacement vs. Normalized Thickness, Layup2. 81 400.0 1.0 - . f T . A 3- v/ 0.8 ~ - — Pagano’s r 0—0 Cubic Zig-Zag 4 — — - FSDT 0.6 - 1 0.4 ~ ~ 0.2 ~ _ 0.0 L A 1 . 1 L 400.0 -200.0 0.0 200.0 GXX Figure 33. Through the Thickness Distribution of (1'M vs. Normalized Thickness, Layup2. 82 1.0 ,. , , - , a d —— Paganos :3 0—0 Cubic Zig-Zag ’ C, -——- FSDT 0.8 ~ 0 .. i ‘1 C1 Ci J c .. 0.6 ~ cl _ <5 <1 q W “l 0.4 - c - 21 ° cl <1 0.2 ~ ci a :, c ‘1 < < 00 fi 1 A l . 1 . -'0.0010 0.0000 0.0010 0.0020 0.0030 0.0040 sz Figure 34. Through the Thickness Distribution of Shear Strain 'sz vs. Normalized Thickness, Layup2. 83 1.0 . ”X: a T . fl b 0.8 ~ o———o Cubic Zig-zag A v 0.5 p g — Pagano’s l 0.4 ~ . 0.2 » X: _ 0.0 . l . L . 1 -0.00005 0.00000 0.00005 0.00010 832 Figure 35. Through the Thickness Distribution of Normal Strain ezZ vs. Normalized Thickness, Layup2. 84 1.0 . , .. —— Pagano’s 0.8 ~ ~ e——+> Cubic Zig-Zag . ’ — - - FSDT 0.4 ~ 0.2 ~ 0.0 . . '0-00010 0.00000 0.00010 u displacement Figure 36. Through the Thickness Distribution of In-plane Displacement vs. Normalized Thickness, Layup3. 85 1.0 r 1 I f 7 0.8 ~ ' J —— Pagano’s °——° Cubic Zig-Zag 0.6 — ——- FSDT - / 0.4 ~ — 0.2 ~ « L 00 V . m . l . 4 . -400.0 -200.0 0.0 200.0 400.0 Figure 37. Through the Thickness Distribution of on vs. Normalized Thickness, Layup3. 86 1.0 f 1 I I I <1 ’_2) + _f + .5 = 1 (5.4.1) sz Syz S SC 100 . . t c . . where 0' 6y: are average 1nterlam1nar shear stresses and oz , oz are average 1nterlam1- XZ ’ nar tensile and compressive normal stresses, respectively. sz, Syz are interlaminar shear strengths and Sz', SZ" are interlaminar tensile and compressive strengths respectively. The average of each components of stresses can be defined as: 1 Aavg 6,]. = 2. I 0‘.de (5.4.2) avg 0 where A. is the distance from the free edge or critical length, Gij are stress components and overbar denotes its average value. For angle ply laminates the dominate stress is interlaminar shear stress Tu and the second, third and fourth terms in Eq. (5.4.1) are negligible for mode II failure, so the qua- dratic failure criterion can be reduced to a simple form: 0’ (4‘) = 1 (5.4.3) where 6'” is the average interlaminar shear stresses as described by Eq. (5.4.2) and sz is the in-situ shear strength of the resin layer present between the plies, not the bulk strength assumed by Sun and Zhou [65]. The value of resin shear strength can be obtained by car- rying out experimental testing for one laminate and then the same value can be used to cal- culate the failure loads for the remaining laminates. 101 5.5 Result and Discussions The interlaminar failure loads are calculated by using the failure criterion as described in previous section and average interlaminar shear stress 1: xz. In both laminates the in-situ strength of the resin layers present between the plies is backed out from experimental results of Sun and Zhou [65]. The in-situ resin strength value obtained for one loading direction is used to predict the failure loads for remaining off-axis loadings. It is also observed that when the in-situ strength of resin is used, the failure loads are not sensitive to the critical length chosen as shown by the-Sun and Zhou [65]. The loads calculated for different critical lengths for (0/90/45/1-45)s and (0/60/-60)s laminates are shown in Table 11 and Table 12. Table 11. Predicted Failure Loads (lb/in) for (0l90/45l-45),3 Laminate at Different Critical Lengths Critical Length Critical Length Critical Length Loading Angle 1.0t 1.5t 2.0t 7.5 3211 3138 3158 15 2911 2842 2848 22.5 2692 2675 2684 30.0 2722 2693 2706 37.5 2996 2996 2996 102 Table 12. Predicted Failure Loads (lb/in) for (0l60/-60)s Laminate at Different Critical Lengths Critical Length Critical Length Critical Length Loading Angle 1.0t 1.5t 2.0t 10 2842 2833 2819 20 1866 1870 1874 30 1822 1822 1822 For off-axis (0/90/45l-45)s laminates, the third interface (45l-45) dictates the laminate failure as it has the lowest interlaminar shear failure load. Table 13 shows the predicted failure loads for (0/90/45l-45)s laminate at the $45 interface. Similarly Table 14 shows the predicted failure loads at the i60 interface of (0/60/-60)s laminate. Table 13. Comparison Failure Loads (lb/in) for (0l90/45/-45)§i Laminate Loading Angle Sun and Zhou Hashim-Rotem Present * 7.5 3076 4813 3158 15 2792 4561 2848 22.5 2585 4581 2684 103 Table 13. Comparison Failure Loads (lb/in) for (0l90/4SI-45)s Laminate Loading Angle Sun and Zhou Hashim-Rotem Present 30.0 2768 4561 2706 37.5 2996 4813 2996 Table 14. Comparison of Failure Loads (lb/in) for (0/60/-60)s Laminate Loading Sun and Zhou Hashin-Rotem Present 10 2204 3586 2819 20 1891 3565 1874 30 1822 1741 1822 The failure criterion adopted in this study assumes that the weakest interface controls the free edge failure and also the laminate strength. Thus, for off-axis specimens of the (0/ 90/45l-45)s laminate, the third interface (45/-45) dictates the laminate failure as it has the lowestint‘erlaminar failure load. Similarly for (O/6OI-6O)S laminate specimens, the (60/-60) interface eontrol the interlaminar failure loads. Table 13 and Table 14 show the compari- SOn of the present model results with Sun and Hashin Rotem results. For both off-axis I aminates the cubic zig-zag model shows good agreement with Sun’s results except for the 1 00 case in (O/6Ol-6O)S laminate. In that case, the predicted load was much higher than the pre sent results. This may be a result of matrix cracking which is not accounted for in the p be sent model. 104 Figure 50 and Figure 51 show the variation of the failure load with the loading in angle in (0l90/45/.45)s and (060l-60)S laminates. The strength prediction using cubic zig-zag model shows good agreement with Sun and Zhou experimental results. 5.6 Conclusion Failure in quasi-isotropic (0/90/45l-45)s and (O/60/-60)s with free edges under off-axis loads are investigated based on cubic zig-zag theory. Classical ply failure theory is not adequate in predicting the strength of a laminate because of that it does not include the Interlaminar stresses. A new failure criterion was developed based on free edge interlami- nar stresses. Failure loads predicted by the present model showed good agreement with the experimental results. Quasi-isotropic laminates with free edge effects showed very aniso- tropic strength behavior with respect to the loading direction. Major failure in these lami- nates are interface shear failure which is dominated by the shear stresses at the free edges. 105 0.03030 . . . . . . , . 4 3 0.03025 ~ 4 ._. 3.03020 ~ « C‘. O E O 8 3. 3.03015 r , - 5 Number of element along x ax1s=4 :5 Number of element along 2 axis=3 0.03010 ~ - 0.03005 - - 0.03000 A 1 A 1 A 1 A 1 A 1 A 0.0 5.0 10.0 15.0 20.0 25.0 30.0 Number of elements along y axis Figure 46. Convergence of U-Displacement at Free Edge as a Function of Number of Elements Along y-axis. 106 0.030225 - 1 3.030220 - 3.030215 - ~ 3.030210 - ' U-displacement 0030205 ’ Number of element along y axis=20 A 0.030200 ~ Number of element along 2 axis=3 - 0.030195 ~ - 0.030190 . 1 . 1 . . . 1 0.0 5.0 10.0 15.0 20.0 25.0 Number of elements along x axis Figure 47. Convergence of U-Displacement at Free Edge as a Function of Number of Elements Along x-axis 107 40000 ' I ' I w 1 v ‘ 30000 T 20000 I tn in (ksi) 0 A 1 0.0 0.2 0.4 y/b Figure 48. Interlaminar stress 1 Distribution Along the Width of (0/90/45l-45)s Laminate at (45/45) Interface of 7.50 Off-Axis Loading Case. Tn in (ksi) 108 40000 . I . r . A . 1 30000 I 20000 I I 1 0000 0 A A A A 0.0 0.2 0.4 0.6 0.8 y/b Figure 49.1nterlaminar stress 13x. Distribution Along the Width of (0/6O/-60)s Laminate at (60/-60) Interface of 100 Off-Axis loading Case. 1.0 Failure Loads in (lb/in) 109 4000.0 I 3000.0 l + Experimental results (C.T. Sun and 8.6. Zhou) —— Present 41 2000.0 0.0 Figure 50. Failure Loads of (0/90/45l-45)s Quasi-Isotropic Laminate vs. 5.0 10.0 15.0 20.0 25.0 30.0 Loading Angle in (0) Loading Angle. 35.0 40.0 Failure Loads in (lb/in) 110 4000.0 ' I I I ' l T T ' I ' l ‘ I + Experimental results (C.T Sun and 8.6 Zhou) —- Present 3000.0 ~ _ I + 4 2000.0 - _ 1 J i 1 000.0 A l A l A I A l A l A 1 L l A _ g 0.0 5.0 10.0 15.0 20.0 25.0 30.0 35.0 40.0 Loading Angle in (0) Figure 51. Failure Loads of (0/60/-60)S Quasi-Isotropic Laminate vs. Loading Angle. CHAPTER 6 Summary and Conclusions In this research, the cubic zig-zag theory and associated finite element model is used to study problems such as free edge stresses in resin-rich zones and failure analysis of lami- nated composites. The advantage of using this theory is that it satisfies interlaminar conti- nuity of transverse shear stresses. Whereas, discrete layerwise theories do not generally satisfy this condition explicitly. However, the discrete layer theories will achieve this result with increasing number of sub-divisions. The associated finite element has the topology of an eight-node brick element with seven degrees of freedom at each node. A sublaminate approach is possible with the cubic zig-zag model because of the way the ele- ment was formulated to allow stacking capability through the thickness direction. This is achieved by including shear traction as degrees of freedom in the element. A review of cubic zig-zag theory and the finite element model is discussed in Chapter 2. The finite element model based on cubic zig-zag theory is used to analyze free edge stresses in (45/-45)s angle-ply, (0/90)s cross-ply and (0/90/45l-45)s graphite epoxy lami- nated plate, each plate being subjected to a uniform axial tensile strain. The results of the finite'element model are found be in good agreement with models in the literature as dis- cussed in Chapter 3. A convergence study is performed to show the importance of mesh creation to find a suitable refinement. lll 112 Results have shown that the magnitude of stresses at the free edge is independent of width but dependent on the thickness of the laminates. It is observed that most dominate stress is IJrz which rises from zero in the interior regions to a very high value at the free edge. A very steep gradient is observed in the other component of interlaminar stress at the free edge. So, it is important to capture these local stresses for good design of laminated composites. The important advantages of the present model are evident when compared to other methods utilizing three-dimensional finite element analysis for studying the problem of three-dimensional edge effects. The present model minimizes the number of elements in the thickness direction because the element is quadratic in shear and also has shear stresses degree of freedom at nodes. It also allows for a very fine mesh near the critical stress gradient regions of a composite laminate. Chapter 4 deals with the concept of resin-rich regions and the response of simply sup- ported composite beams with resin and delaminated layers. It is usually found that in spite of careful fabrication of composite materials, there exists a region of resin zone or area of less fibers. These regions at the interlaminar boundaries of plies provide the path for cracks to initiate and grow because of their low strength or toughness. So, it is necessary to include these regions from the analysis point of view and predict the local stresses and strains in these regions accurately and efficiently. Results for four different layups have shown that cubic zig-zag theory is accurate com- pared to first order shear deformation theory (FSDT) in predicting stresses and strains. Cubic zig-zag model is capable of accommodating an extremely thin resin layer between plies without increasing the number of degrees of freedom and without numerical prob- 113 lems. FSDT didn’t account for the presence of thin resin layers at the ply interfaces. Whereas, zig-zag model is able to capture large discontinuous transverse shear strains and continuous shear stresses at the ply interfaces of the laminate. In Chapter 5, based on the stress state in resin layers a failure criterion is developed to study the failure loads due the free edge effects. Quasi-isotropic laminates under off-axis tensile loading are analyzed. Results are shown to compare favorably with those in the lit- erature. These laminates showed anisotropic behavior in strength with loading direction. It is observed that the interface failure is dominated by the large shear stresses at the free edge of composite laminates. A finite element model based on cubic zig-zag theory is developed to predict free edge interlaminar stresses in laminated composites. The results are in good agreement with models in the literature. It has been demonstrated that present model can capture the high interlaminar stresses with fewer elements compared to conventional three-dimensional continuum based elements. Successfully predicted the stresses and strains in the resin-rich regions which exist at the interface of plies of laminated composites and is area of critical importance in the case of delamination. A new failure criterion is developed to predict the failure loads based on the stress state in resin for quasi-isotropic laminates under off-axis loading. 6.1 Recommendation for Future Work Present model can be extended to study the problem such as adhesive joints where again the interlaminar stresses are very important. Because of limited computational resources it is important to have models which are computationally efficient to model pro- 114 gressive failure or crushing of tubes etc. With small and fixed number of degrees of free- dom, present model has necessary kinematic to analyze these problems. This can be implemented by developing a progressive failure algorithm which can be incorporated into finite element model. A failure criterion based on critical energy release rate can be develop to study delamination onset and growth in laminated composites. 10 11 List of References Timoshenko, S. A. and Woinowsky- Krieger, S, “Theory of Plates and Shells,” 2nd ed., McGraw-Hill, Inc., New York (1959) Mindlin, R.D., “Influence of Rotary Inertia and Shear on Flexural Motion of Isotropic Elastic Plates,” J. of Applied Mechanics, Vol. 18, 1951, pp.3l-38 Reissner, E., “The Effect of Transverse Shear Deformation on the Bending of Elastic Plates,” J. of Applied Mechanics, Vol. 12, 1945, pp.69-76 Yang, P. C., and Norris, C. H., and Stavsky, Y., “Elastic Wave Propagation in Hetero- geneous Plates,” Int. J. Sol. Structure, Vol. 2, 1966, pp. 665-684. Lo, K. H., Christensen, R. M., and Wu, E. M., “A High-Order Theory of Plate Defor- mation - Part 1: Homogeneous Plates,” Journal of Applied Mechanics, Vol.44, 1977, pp. 663-668. Reddy, J. N ., “A Simple Higher-order Theory of Laminate Composite Plates,” Journal of Applied Mechanics, Vol. 51, (1984), pp. 745-752 Barbero E.J. and Reddy J .N., “An Accurate determination of stresses in thick lami- nates using a generalized plate theory”, Int. J. Numer. Meths. Engrg., Vol. 29, 1990, pp. 1-14 Tessler, A., “A Higher-Order Plate Theory with Ideal Finite Element Suitability”, Comp. Meth. Appl. Mech. Eng., Vol. 85, 1991, pp. 183-205. Pagano, N. J ., “Exact Solutions for Composite Laminates in Cylindrical Bending,” J. Of Composite Materials, Vol. 3, 1969, pp. 48-111. Reddy, J .N., “Generalization of Two-Dimensional Theories of Laminated Composite Plates”, Comm. Appl. Num. Meths. Vol. 3, 1987, pp.173-180. Toledano, A., and Murakami, H., “A Composite Plate Theory for Arbitrary Laminate Configuration,” Journal of Applied Mechanics, Vol. 54, 1987, pp181-189 115 116 12 Lu, X. and Liu, D., “An Interlaminar Shear Stress Continuity Theory for Both Thin and Thick Composite Laminates”, J. Appl. Mech., Vol.59, 1992, pp.502-509. 13 Lee, CY. and Liu, D., “An Interlaminar Stress Continuity Theory for Laminated com- posite analysis”, Computers & Structures, Vol.42, No.1, 1992, pp. 69-78. 14 DiSciuva, M., “Development of An Anisotropic, Multi-layered, Shear-Deforrnable Rectangular Plate Element,” Computers and Structures, Vol. 21 1985, pp. 789-796. 15 DiSciuva, M., “An Improved Shear-Deforrnation Theory for Moderately Thick Multi- layered Anisotropic Shells and Plates,” J. Appl. Mech., Vol. 54, 1987, pp. 589-596. 16 DiSciuva, M., “A General Quadrilateral Multilayered Plate Element with Continuous Interlaminar Stresses,” Computers and Structures, Vol. 47, 1993, pp. 91-105. 17 Murakami, H., “Laminated Composite Plate Theory with Improved In-Plane Response,”, J. of Applied Mechanics, Vol. 53, 1986, pp. 661-666. 18 Cho, M. And Parmerter, R.R., “Efficient Higher Order Composite Plate Theory for General Lamination Configurations”, AIAA Journal, Vol. 31, 1993, pp. 1299-1306. 19 Averill RC. and Yip, Y. C., “Development of Simple, Robust Finite Elements Based on Refined Theories for Thick Laminated Beams”, Computers and Structures, 1996a 20 Yip, Y. C., “A 3-D Laminate Plate Finite Element Model With Zig-Zag Sublaminate Approximations For Composite and Sandwich Panels”, Doctoral dissertation, Depart- ment of Material Science and Mechanics, Michigan State University. 1996 21 Resisner, E., “On A Certain Mixed Variational Theorem and A Proposed Application,” International Journal for Numerical Methods in Engineering, Vol. 20, 1984, pp. 1366- 1368. 22 Barbero, E. J ., “3-D Finite Element Modeling of Laminated Composites by Layerwise Constant Shear Elements,” Enhancing Analysis Techniques For Composite Materials, Vol. 10, 1991, pp. 229-237 23 Agarwal, B. D. and Broutrnan, L. J., “Analysis and Performance of Fiber Composites,” John Wiley & Sons, Inc. 24 Graber, R. P., “Tensile-Stress Strain Behavior of Graphite/Epoxy Laminate,” NASA CR 3592, NASA, Washington, DC (1982). 25 Puppo, A. H., and Evensen, H. A.m “Interlaminar Shear in Laminated Composites Under Generalized Plane Stresses,” J. Composite of Material, Vol. 4 1970 pp. 204-220 26 Pipes, R. B. and Pagano, N. J., “Interlaminar Stresses in Composites Under Axial 117 Extension,” Journal of Composite Material, Vol. 4 (1970) pp. 538-548. 27 Tang, S. and Levy, A., “A Boundary Layer Theory -Part 11: Extension of Laminate Finite Stripe” Journal of Composite Material, Vol. 9, (1975) pp. 42-52.s 28 Pagano, N. J. and Pipes, R. B., “The Influence of Stacking Sequence on Laminate Strength,” Journal of Composite Material, Vol. 5 (1971) pp. 50-57. 29 Pagano, N. J. and Pipes, R. B., “Some Observations on the Interlaminar Strength of Composite Laminates,” Int. J. Mech. Sci Vol. 15 (1973) pp. 679-688. 30 Wang, S. S. and Choi, 1., “Boundary Layer Effects in Composite Laminates: Part I - Free Edge Stress Singularities,” Journal of Applied Mechanics, Vol. 49 (1982) pp. 541- 548. 31 Wang, S. S. and Choi, 1., “Boundary Layer Effects in Composite Laminates: Part II - Free Edge Stress Solutions and Basic Characteristics,” Journal of Applied Mechanics, Vol. 49 (1982) pp. 549-560 32 Wang, S. S. and Yuan, F. G., “A Singular Hybrid Finite Element Analysis of Boundary Layer Stresses in Composite Laminates,” Int. Journal of Solids & Struct.ures, Vol. 19 (1983) pp. 825-837 33 Whitney, J. M. and Browning, C. E., “Free edge Delamination of Tensile Coupons,” J. Composite of Material, Vol. 6 (1973) pp. 300-303. 34 Wang, A. S. D. and Crossman, F. W., “Some New results on Edge effect in Symmetric Composite Laminates,” Journal of Composite Material, Vol. 11 (1977) pp. 92-106 35 Wang, A. S. D. and Crosswan, F. W., “Calculation of Edge Stresses in Multi-layer Laminates by Sub-structuring,” Journal of Composite Material, Vol. 12, (1978) pp.76- 83.. 36 Pagano, N. J ., “Stress Field in Composite Laminates,” Int. Journal of Solids & Struc- tures, Vol. 14 (1978) pp. 385-400 37 Pagano, N. J., “Free Edge Stress Field in Composite Laminates,” Int. Journal of Solids & Structures. Vol. 14 (1978) p.401. 38 Raju, I. S., and Crews, J. H. Jr., “Interlaminar Stress Similarities at a Straight Free Edge in Composite Laminates,” Computers and Structures, Vol. 14, 1981, pp. 21-28 39 Whitcomb, J. D. and Raju, I. 5., “Analysis of Interlaminar Stresses in Thick Compos- ite Laminates with and without Edge Delamination,” ASTM STP 876, (1985) pp. 69- 94. 118 40 Jones, R., Callinan, R., Teh, T. K. and Brown, K. C., “Analysis of Multi-Layer Larni- nates Using Three-Dimensional Super-elements,” Int. Journal of Numerical Methods in Engineering, Vol. 20, (1984) pp. 583-587. 41 Kim, J. K. and Hong, C. S., “Three-Dimensional Finite Element Analysis of Interlam- inar Stresses in Thick Composite Laminates,” Computers & Structures, Vol. 40, No. 6. 1991. pp. 1395-1404. 42 Tessler A. and Hughes T.J.R., 1983, “An Improved Treatment of Transverse Shear in the Mindlin-Type Four-Node Quadrilateral Element”, Computer Methods in Applied Mechanics and Engineering, Vol. 39, pp. 311-335. 43 Tessler A. and Hughes T.J.R., 1985, “A Three-Node Mindlin Plate Element with Improved Transverse Shear”, Computer Methods in Applied Mechanics and Engineer- ing, Vol. 50. PP. 71-101. 44 Weeknights, S. G., “Theory of Elasticity of an Anisotropic Elastic Body,” Holden-Day (1963) 45 Ting, T. C. T., and Chou, S. C., “Stress Similarities in Laminated Composite,” Proc. Second U.S.A- U.S.S.R symp. on Fracture of Composite Materials, 1981, pp. 265 46 an, W. L., “Free Edge Effects in Laminates Under Extension, Bending and Twisting, Part I: A Stress Function Approach,” Proceedings of the AIAA/ASME/AHS 32th Structures, Structural Dynamics and Material(SDM) Conference, 1991, AIAA Paper No. 91-0959-CP, Part II, pp. 985-995 47 Yin, W. L., “Free Edge Effects in Laminates Under Extension, Bending and Twisting, Part II: Sublaminate/Layer Modeling and Analysis,” Proceedings of the AIAA/ASME/ AHS 33th Structures, Structural Dynamics, and Material(SDM) Conference, 1992, AIAA Paper No. 92-2228-CP, Part 1, pp.48-58. 48 Soni, S. R. and Kim, R. Y., “Delamination of Composite Laminates Stimulated by Interlaminar Shear,” ASTM STP 893, 1986, pp. 286-307. 49 J00, J. W. and Sun, C. T., “A Failure Criterion for Laminates Governed by Free Edge Interlaminar Shear Stress,” J. of Composite Materials, Vol. 26, 1992, pp. 1511-1521. 50 O’Brien, T. K., “Characterization of Delamination Onset and Growth in a Composite Laminate,” Damage of Composite Materials, ASTM, STP 775, 1982, pp. 140-167. 51 O’Brien, T.K., “Analysis of Local Delamination and their Influence on Composite Laminate Behavior,” Delamination and Debonding of Materials, ASTM STP 876, 1985. PP. 282-297. 52 Rybicki, E. F.and Schmueser, D. W. and Fox, J ., “An Energy Release Rate Approach 119 for Stable Crack Growth in Free edge Delamination problem,” J. Composite Materials, Vol. 11, 1977. pp. 471-487 ~ 53 Song, SJ. and Wass, A. M., “A Spring Foundation Model for Mode I Failure of Lami- nated Composites Based on an Energy Criterion,” J. Engng. Tech., Vol. 116, 1994, pp. 512-516 54 Tsai, SW. and V. D. Azzi., “Strength of Laminated Composite Materials,” AIAA Jour- nal, Vol. 4, 1966. Pp. 296-301 55 Cantwell, W. J. and Morton, J ., “The Significance of Damage and Defects and Their Detection in Composite Materials: Review,” J. of Strain Analysis, Vol. 27, No.1, 1992, pp. 29-42. 56 Aitharaju, V. R. and Averill, R.C “An Assessment of Zig-Zag Kinematic Displacement Models for the Analysis of Laminated Composites,” Submitted for Mechanics of Com- posite Materials and Structures, 1997. 57 Reissner, E., “On a Mixed Variational Theorem and a Shear Deformable Plate The- ory,” International Journal for Numerical Methods in Engineering, Vol. 23, 1984, pp. 193-198. 58Rotem, A., “Fatigue Failure Mechanism of Composite Laminates,” Mechanics of Com- posite Materials, 1982, pp. 421-435. 59 Soni, S. R. and Pagano, N. J ., “Elastic Response of Composite Laminates,” Mechanics of Composite Materials, 1982, pp.227-242. 60 Herakovich, C. T., Nagarkar A. and O’Brien, D. A., “Failure Analysis of Composite Laminates with Free edge,” Modem Developments in Composite Materials and Struc- tures, Winter Annual Meeting of ASME, 1979, pp. 53-66 61 Kim, R. Y. and Soni, S. R., “Experimental and Analytical Studies on the Onset of Delamination in Laminated Composites,” J. of Composite Materials, Vol. 18, 1984, pp. 70-80. 62 Zhou, S. G., “Failure of Quasi-Isotropic Composite Materials under Off-Axis Load- ing,” Masters Thesis, Purdue University, 1994. 63 Hashin, Z, “Failure Criteria for Uni-direction Fiber Composite,” J. of Applied Mechanics,” Vol. 47, 1980, pp. 329-335 64 Herakovich, C. T., “Influence of Layer Thickness on the Strength of Angle-Ply Lami- nates,” J. of Composite Materials, Vol. 16, 1982, pp. 216-227. 65 Sun, C. T. and S. G. Zhou., “Failure of Quasi-isotropic Composite Laminates with 120 Free Edges under off-Axis Loading,” J. of Reinforced Plastics and Composites, Vol. 7, 1988, pp. 515-557. 66 Hashin, Z. and Rotem, A., “A Fatigue Failure criterion for Fiber Reinforced Materi- als,” J. of Composite Materials, Vol. 7, 1973, pp. 488-464 67 Rotem, A. and Hashin, 2., “Fatigue of Angle-Ply Laminates,” AIAA Journal, Vol. 14, 1976, pp. 868-872. 68 O’Brien, T. K., “Characterization of Delamination Onset and Growth in a Composite Laminate,” Damage of Composite Materials, ASTM, STP 775, 1982, pp. 140-167. 69 O’Brien, T.K., “Analysis of Local Delamination and their Influence on Composite Laminate Behavior,” Delamination and Debonding of Materials, ASTM STP 876, 1985. Pp. 282-297. 70 Rybicki, E. F.and Schmueser, D. W. and Fox, J ., “An Energy Release Rate Approach for Stable Crack Growth in Free edge Delamination problem,” J. Composite Materials, Vol. 11, 1977. pp. 471-487 71 Wang, A. S. D and Crossman, F. W., “Initiation and Growth of Transverse Cracks and Edge Delamination in Composite Laminates: Part 1. An Energy Method,” J. of Com- posite Materials, Vol. 14, 1980, pp.71-87. 72 Wang, A. S. D and Crossman, F. W., “Initiaion and Growth of Transverse Cracks and Edge Delamination in Composite Laminates: Part 2. An Energy Method,” J. of Com- posite Materials, Vol. 14, 1980, pp.71-87. 73 Rodini, B. T. and Eisenmann, J. R., “An Analytical and Experimental Investigation of Edge Delamination in Composite Laminates,” Proc 4th Conf. Fibrous Composites, San Diego, Cal. 1978. 74 Wu, E. M., “Failure Analysis of Composites with Stress Gradients.” Fracture of Com- posites Materials. 1979, pp.63. 75 Rybicki, E. F. and Kanninen, M. F., “A Finite Element Calculation of Stress Intensity Factors by A Modified Crack Closure Integral,” Eng. Fracture Mechanics, Vol. 9, 1977,pp.470. 76 Irwin, G. R., “Fracture,” Handbuch der Physik, Vol. 6, Springer-Verlag, 1958, pp5.