7. wék‘i‘" {rL».n§£ .\.A ,fl‘ w-" _ .~~"":. ‘ \ - {iyfi ' $ t 7i. '11 5 , 29$?ng i."* s” ijllJv This is to certify that the thesis entitled A PROGRESSIVE DAMAGE ANALYSIS INCLUDING DELAMINATION FOR FIBER REINFORCED COMPOSITES presented by Sivom Manchiraju has been accepted towards fulfillment of the requirements for the MS. degree in Mechanical Engineering “04%;: OZ; Major Professor’s Signature 9/1/04 Date MSU is an Affirmative Action/Equal Opportunity Institution - - _.--_<:g.—.-.-.-.-.- - LIBRARY Michigan State University PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE MAY“9 151320083 6/01 cJClRC/DatoDuepss-p. 15 A PROGRSSIVE DAMAGE ANALYSIS INCLUDING DELAMINATION FOR FIBER REINFORCED COMPOSITES By SIVOM MANCHIRAJU A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Mechanical Engineering 2004 ABSTRACT A PROGRESSIVE DAMAGE MODEL INCLUDING DELAMINATION FOR FIBER REINFORCED COMPOSITES By Sivom Manchiraju A finite element model for progressive in-plane damage and interfacial delamination was presented for cross-ply composite laminates. The commercial finite element software ABAQUS/STANDARD was used for the thesis research. Each composite laminae was modeled as layers of shell elements based on Mindlin’s plate theory, thereby achieving computational efficiency. Displacement continuity was established across the layers of shell elements. Interlaminar stresses were recovered at the interface with the use of multi- point constraints or an interface layer. Several material constitutive relations were developed to simulate the in-plane damage, such as fiber breakage, matrix cracking and in-plane shear failure. Interlaminar delamination was simulated by removing the displacement constraints or by reducing the stiffness of the interface layer. The delamination scheme was validated for a composite beam subjected to sinusoidal loading by an elasticity solution extended from the Pagano’s problem. The finite element analysis results matched well with the elasticity solution for different interface bonding conditions. Composite laminates with different stacking sequences were subjected to impact loading. The results from the finite element analysis compared well with the experimental results qualitatively. The typical penny-shaped delamination area was observed at the interface and matrix cracking was observed along the fiber direction. However, the delamination size was about 50 % lower than the experimental results. To my parents and brother ACKNOWLEDGMENTS First of all I would like to thank my advisor, Dr. Dahsin Liu for his advice and guidance through out my study. He not only gave me the inputs but also provided support and encouragement all along the study. I would also like to thank Dr. Alfred Loos and Dr. Ronald Averill for serving in this thesis defense committee. I would like to thank my colleagues who provided me the experimental results in order to validate my computational studies. I would also like to thank all the professors who gave me the required background in mechanical engineering. Finally I would like to thank my family and friends who always supported and motivated me. TABLE OF CONTENTS List of Figures ................................................................................................................... vii 1. INTRODUCTION .......................................................................................................... 1 1.1 Background and Motivation ..................................................................................... 1 1.2 Damage Modeling: Background Theory and Literature Survey ............................... 2 1.2.1 Failure Criteria ................................................................................................... 3 1.2.2 Damage Models ................................................................................................. 4 1.2.3 Delamination Models ......................................................................................... 6 1.3 Organization of Thesis .............................................................................................. 9 2. ELASTICITY STUDIES OF DELAMINATION ........................................................ 11 2.1 Extended Pagano’s Problem ................................................................................... 11 2.1.1 Fundamental Equations .................................................................................... 12 2.1.2 Boundary Conditions ....................................................................................... 13 2.1.3 Interfacial Conditions ....................................................................................... 14 2.1.4 Elasticity Solution ............................................................................................ 16 2.2 Finite Element Modeling of Delamination ............................................................. 17 2.2.1 Multi-Point Constraint Method ........................................................................ 18 2.2.2 Interfacial Layer Method ................................................................................. 21 2.3 Comparison between Elasticity and Finite Element Analysis ................................ 25 2.3.1 Perfectly Bonded Composite Beam ................................................................. 25 2.3.2 Imperfectly Bonded Composite Beam ............................................................. 28 2.3.3 Further Elasticity Analysis ............................................................................... 32 3. In-Plane Failure ............................................................................................................. 38 3.1 Basic Assumptions .................................................................................................. 38 3.2 Damage Modes and Stiffiiess Degradations ........................................................... 39 3.2.1 Fiber Damage ................................................................................................... 40 3.2.2 Matrix Damage ................................................................................................ 43 3.3 Validation of Dynamic Analysis ............................................................................. 45 3.4 Simulation of [0/0] Laminates ................................................................................ 47 3.4.1 Matrix Failure Only ......................................................................................... 48 3.4.2 Brittle Fiber Failure .......................................................................................... 49 3.4.3 Progressive Fiber Failure ................................................................................. 49 4. Delamination Failure .................................................................................................... 54 4.1 Delamination Modeling and Criteria ...................................................................... 54 4.1.1 Delamination Failure Only .............................................................................. 56 4.1.2 In-Plane Failure and Delamination Failure ...................................................... 58 4.1.3 Modified Delamination Failure ........................................................................ 59 4.2 Investigation of Delamination ................................................................................. 61 4.2.1 Delamination Shape ......................................................................................... 61 4.4.2 Delamination Size ............................................................................................ 63 5. Simulations of [0/90/0] Laminates ................................................................................ 66 6. Conclusions and Recommendations ............................................................................. 74 6.1 Conclusions ............................................................................................................. 74 6.2 Recommendations ................................................................................................... 76 Appendix A ....................................................................................................................... 77 Matlab code to for the analysis of imperfectly bonded composite laminate Appendix B ....................................................................................................................... 83 Subroutine MPC to enforce the displacement continuity conditions Appendix C ....................................................................................................................... 86 Subroutine URDFIL to calculate the interlaminar stresses from the nodal forces Appendix D ....................................................................................................................... 88 Subroutine USDFLD for changing the stiffness of the interface layer for the delamination model with interface layer Appendix E ....................................................................................................................... 9O UMAT subroutine for the material in-plane damage BIBLIOGRAPHY ............................................................................................................. 95 vi Fig. 2.1 Fig. 2.2 Fig. 2.3 Fig. 2.4 Fig. 2.5 Fig. 2.6 Fig. 2.7 Fig. 2.8 Fig. 2.9 Fig. 2.10 Fig. 2.11 Fig. 2.12 Fig. 2.13 Fig. 2.14 Fig. 2.15 Fig. 2.16 List of Figures A [90/0] beam subjected to sinusoidal loading on the top surface ......... 11 Two shell elements with multi-point constraints to enforce Interfacial displacement continuity ............................................. 19 Two shell elements with springs to enforce the interfacial displacement continuity .......................................................... 20 Finite element mesh showing interfacial area at a node, used to calculate the interlaminar stresses from the nodal forces .................... 20 Composite laminate with a thin interfacial layer .............................. 22 Normalized in-plane stress at the mid-span through the thickness of [90/0] beam ..................................................................... 26 Normalized shear stress along the length of [90/0] beam .................... 27 Normalized normal stress along the length of [90/0] beam .................. 28 Displacement mismatch at the interface at different bonding conditions ........................................................................... 29 Normalized deflections of [90/0] beam for different shear bonding cond1t10ns30 Normalized deflections of [90/0] beam for different bonding conditions ........................................................................... 31 Interlaminar shear stress at the edge for different shear bonding conditions ................................................................ 31 Interlaminar normal stress at the center for different shear bonding conditions ................................................................ 32 Transverse shear stress with different shear slip coefficients and perfect normal bonding ...................................................... 33 Transverse shear stress at the edge for large shear slip coefficient and different normal separation coefficients .................... 34 Interlaminar shear stress distribution along the length of the beam for different shear slip coefficients with perfect normal bonding. ......... 35 vii Fig. 2.17 Fig. 2.18 Fig. 2.19 Fig. 3.1 Fig. 3.2 Fig. 3.3 Fig. 3.4 Fig. 3.5 Fig. 3.6 Fig. 3.7 Fig. 4.1 Fig. 4.2 Fig. 4.3 Fig. 4.4 Fig. 4.5 Fig. 4.6 Fig. 4.7 Normalized transverse normal stress distribution for different normal separation coefficients with no shear slip ............................. 35 Normalized transverse normal stress distribution for different normal separation coefficients with large shear slip ......................... 36 Normalized interlaminar normal stress along the length of the beam for different normal separation coefficients with perfect shear bonding ...................................................................... 37 Uniaxial stress-strain curve along the fiber direction ........................ 41 Force and deflection history from experiment and simulation .............. 47 Force-deflection curves for different impact velocities based on material model with matrix failure only ....................................... 48 force-deflection curve based on a material with brittle fiber failure ....... 49 Force-deflection curves for different impact velocities based on progressive fiber failure model ...................................... 50 Comparison of force-deflection curves at three different impact velocities .................................................................. 52 Absorbed energy vs. impact energy ........................................... 53 Delamination at the bottom interface of [0/90/0] laminate with no in-plane failure ................................................................. 57 Delamination at the bottom interface of [0/90/0] laminate with in-plane failure criterion ................................................... 58 Delamination at the bottom interface of [0/90/0] laminate based on modified delamination failure criterion ............................ 61 Axis of delamination rotates with the fiber orientation of the bottom lamina ................................................................. 62 Delamination area vs. impact energy .......................................... 63 Delamination at the bottom interface of [0/90/0] laminate under various impact energies ................................................... 64 Delamination area vs. fiber angle between laminae .......................... 65 viii Fig. 5.1 Fig. 5.2 Fig. 5.3 Fig. 5.4 Fig. 5.5 Fig. 5.6 Fig. 5.7 Force-deflection curves for [0/90/0] laminate for different impact velocities .................................................................. 67 Comparison of force-deflection curves at three different impact velocities ........................................................ 68 Delamination shape at the bottom interface of [0/90/0] laminate at three different impact velocities .............................................. 69 Delamination area vs. impact energy ............................................ 7O Absorbed energy vs. impact energy ............................................ 71 F orce-deflection curve for [0/90/0] laminate with and without delamination damage .................................................... 72 F orce-deflection curves for different composite laminate lay-ups .............................................................................. 73 1. INTRODUCTION 1.1 Background and Motivation The high strength-to-weight and high stiffness-to-weight ratios of laminated composites have made them popular in a wide variety of engineering applications, especially in high-performance structures. However, laminated composite structures can take maximum advantage of these properties only when their designs do not overlook the composite weaknesses due to various damage modes. Once the damage modes are understood, a better composite structure design can be achieved. When used in defense or aerospace applications, laminated composites are frequently subjected to impact loading. They are quite susceptible to impact loading and their properties are severely degraded upon impact. Because of their heterogeneity and anisotropy within individual laminae and through the laminate thickness, they undergo complex damage process when subjected to impact loading. Because the damage process is critically important to the success of laminated composites, there is a great need to understand the details of it. The studies of the damage process of laminated composites subjected to impact loading have been performed both experimentally and computationally. Experimental studies involve preparation of test specimens and impact loading the specimens. Drop weight test or gas gun test are usually used for impact investigations. Based on the impact energy and absorbed energy and the inspection of the damaged specimens, the damage process is studied. For the computational studies of the damage process, finite element method is frequently used. Computational studies are inexpensive and provide a lot of scope for parametric studies leading to a better understanding of the damage process. But, it is very difficult to develop a computational model which is able to capture the damage process and at the same time is computationally efficient. This thesis research attempts to develop such a finite element procedure for the study of the damage process of laminated composites. 1.2 Damage Modeling: Background Theory and Literature Survey The finite element method is an effective and efficient tool for computational studies. Its use in fiber-reinforced laminated composites has been an active area of research for quite some time. Because of the complexity of the damage mechanisms that can occur in laminated composites, the overall mechanical behavior of laminated composites is highly nonlinear and is extremely difficult to model. In order to model the nonlinear response of laminated composites, continuum models are used. Representation of individual composite laminae by a homogenized continuum model permits the overall simulations of lamina performance and the damage process. A composite laminate can be modeled by using two-dimensional shell elements which use a laminate theory such as the classical laminate theory, Mindlin’s plate theory, etc. A composite laminate can also be modeled by using three-dimensional brick elements. Three-dimensional elements are able to capture the entire three-dimensional stress distribution (at least theoretically) but they are computationally expensive. Thus, shell elements are the most commonly used elements for the study of composite laminates. A commercial finite element software ABAQUS was used for this thesis research because user defined subroutines could be added to the software with ease. ABAQUS uses shell elements based on the first-order shear deformation theory in which the transverse shear strain is assumed to be constant through the thickness of the shell. A number of shell elements have been presented, based on higher-order shear deformation theories which could calculate the transverse stresses more accurately for the study of laminated composites [1, 2]. The damage process in composite laminates is a three-dimensional process. There are intralaminar damage in the form of fiber damage, matrix damage and fiber-matrix debonding and interlaminar damage in the form of delamination. 1.2.1 Failure Criteria In the analysis of composite laminates, a wide range of failure criteria have been used. A failure criterion may combine all the possible damage modes into one equation, e.g. Tsai-Hill failure criterion, Tsai-Wu failure criterion, etc. Failure criteria may also consist of a set of equations; each representing a distinct damage mode. The first two major failure criteria were proposed by Hashin [3, 4]; one for fiber damage and the other for matrix damage. The failure criteria were quadratic in nature since a linear criterion usually over predicted the damage while a polynomial criterion was too complicated. The fiber damage due to tension depends upon the normal stress in the fiber direction and the shear along the fiber. The matrix damage depends upon the normal stress in the matrix direction and the shear stress along the fiber. Yamanda and Sun [5] proposed a failure criterion of a lamina similar to Hashin’s fiber failure criterion but used in-situ strengths in the criterion. Chang and Chang [6] proposed phenomological failure criteria which included criterion for matrix cracking, fiber damage and fiber-matrix shear failure. The fiber damage was assumed to be caused by the stress in the fiber direction and the fiber-matrix shear failure was assumed to be due to the normal stress in the fiber direction and the shear stress. Their failure criterion for fiber-matrix shear failure was identical to Hashin’s failure criterion for fiber damage while their matrix failure criterion was also similar to Hashin’s. Chang and Lessard [7] proposed another failure criteria accounting for the non-linear nature of shear stress-strain relation in the composite materials. They also used in-situ strengths, similar to Yamanda and Sun, in their criteria. Shahid and Chang [8] later on presented a finite element code for the analysis of composite materials based on this failure criterion. Many other failure criteria were available in the literature. Some of them were based on the failure of isotropic materials and then extended to orthotropic behavior by superimposing several cutoff lines that took into account the experimental results [9] while some used statistical technique to generate failure envelopes [10]. 1.2.2 Damage Models Once a failure criterion was met, the material underwent progressive damage. Different damage models were proposed to describe the post-failure behavior of the composite laminates. One of the simplest and most widely used methods for the post- failure modeling was the brittle failure model. Chang and Chang [6] proposed to reduce the elastic constants to zero in order to simulate the post-failure response of composites. This model was called a brittle failure model as the material was assumed to be perfectly brittle and there was a sharp and steep drop in the strength of a damaged composite material. This model was used extensively for the analysis of composite materials [7, 8, 11, 12, 29-34]. It was also used for the damage in the matrix direction in this thesis research. More recently, progressive failure models with post-failure material softening received a great deal of attention. In general, a macroscopic description of damage could be captured by a constitutive model that exhibited a decrease of stress with increasing strain, i.e. strain softening, beyond the point of failure. This strain softening was based on continuum damage mechanics (CDM) theories [13] in which the net effect of fracture was idealized as a degradation of the elastic moduli of the composite materials. All CDM models provided a mathematical description for the dependence of the elastic moduli on the damage state and the evolution of the damage with respect to the loading and unloading behavior of the composite materials. Typically, CDM models used a simple predefined stiffness reduction function, for example, a linear softening relation, an exponential softening relation, or more complex damage growth laws such as Weibull’s function. Bazant and coworkers presented a series of crack band models based on CDM [14-16]. They were extensively used for quasi-brittle materials such as fiber composites and concrete. Johnson et a1. [17, 18] proposed a CDM model with linear softening for laminated composites. They tried to correlate the area under the softening zone with the material fracture toughness. Matzenmiller et al. [19] proposed a general constitutive model based on CDM for brittle materials. The damage function could be of any form. They adopted a damage function based on Weibull’s statistical function for illustration purpose. Williams et a1. [20] implemented the constitutive model proposed by Matzenmiller in commercial finite element code LS-DYNA and studied impact response of laminated composites. They also compared this strain softening model with the brittle failure model proposed by Chang [7]. They concluded that in the post-failure region, the strain softening model was in a better agreement with the experiments when compared to the brittle failure model. Schapery and Sicking [21] used a thermodynamically based progressive damage formation and growth model. They used a polynomial evolution law for the damage in the matrix direction. Basu and Waas [22] extended the work done by Schapery and Sicking to include the fiber failure. They assumed that the failure in the fiber is brittle. Upon failure, the stress in the fiber direction dropped to a plateau which was maintained upon further loading. This fiber failure model would be used in this thesis research. Bazant [15] showed strain softening causes strain localization and thus the strain softening models are not objective. He proposed a non-local continuum damage model for the analysis of brittle materials to overcome strain localization [16]. Kennedy and Nahan [23] implemented a non-local continuum damage model based on exponential softening in LS-DYNA for the analysis of composite materials. 1.2.3 Delamination Models The computational modeling of interlaminar failure, i.e. delamination, was equally challenging. A composite laminate is formed by stacking together individual laminae and these laminae can have different fiber orientations, resulting in non-uniform material properties through the laminate thickness. The mismatch of the material properties between the adjacent laminae could cause the laminae to separate from each other when subjected to loading. This kind of failure was called delamination. Liu [24] proposed a bending stiffness mismatch theory to predict the delamination response of the laminated composites. It stated that the amount of delamination was dependent upon the mismatch of bending stiffness between adjacent laminae. Delamination modeling posed a unique challenge in the composite laminate analysis as it was a three-dimensional damage thus required stress and displacement components in all three dimensions. Delamination modeling could be divided into two groups based on the types of failure criteria being used for the prediction of delamination: strength of materials based approach and fracture mechanics based approach. In both approaches the three-dimensional state of the composite laminate had to be known. Brewer and Lagace [25] proposed a quadratic criterion based on the interlaminar stresses for delamination prediction. Zhou and Sun [26] formulated the same interlaminar failure criterion utilizing interface stresses. Kim and Soni [27] proposed averaging the interlaminar stresses over an arbitrary critical length of one ply thickness. Fenske and Vizzini [28] proposed a failure criterion which included the in-plane stresses at the interface in the quadratic delamination criterion. Three-dimensional brick elements were used in order to calculate the interlaminar stresses. De Moura and Maques [29] used the interlaminar stresses calculated by a shear flexible shell element to predict the delamination in composite laminates subjected to impact loading. However, their study was limited to predicting the delamination shape and not to simulate delamination. Lou et a1. [30] used three-dimensional brick elements and used Chang brittle failure criteria to study impact response of laminated composites. The quadratic failure criterion was used for delamination. Once the criterion was met, the element’s stress bearing capacity in the out-of-plane direction was reduced to zero. Other researchers have also used a similar finite element model to predict the impact response of laminated composites [31-34]. These models were computationally expensive owing to the use of three-dimensional brick elements. Moreover, a very refined mesh was required in order to calculate the inlerlaminar stresses accurately when using a brick element. Li et a1. [35] proposed an interface element which joined the brick elements on either side of the potential delamination interface. Once the delamination took place, the stiffness of the interface element was changed to make the displacements discontinuous across the interface. Sprenger et a1. [36] used a similar interface element to simulate delamination. Because of high computational cost involved in using brick elements, focus was shifted to model delamination by using computationally efficient shell elements. Multi- layered shell elements were used to model the composite laminates. Composite laminae on either side of the potential delamination surface were modeled by shell elements. Displacement continuity conditions were enforced between the shell elements. Zheng and Sun [37] imposed constraints between appropriate nodes of the shell elements to enforce displacement continuity. The strain energy release rate was calculated at the interface by the virtual crack closure method and the constrained conditions between the nodes were removed when the strain energy release rate was higher than the critical value. Klug and Sun [38] also used this model for the study of delamination. Similar constrained conditions were used in this thesis research. Wisnom and Chang [39] developed three-dimensional spring element to enforce the displacement continuity between the corresponding nodes. Their study, however, was based on plain strain elements. The area under the force-deflection curve of the spring was used to calculate the strain energy release rate. Once the delamination criterion was met, the stiffiiess of the spring was reduced to simulate displacement discontinuity across the interface. Other researchers used specially developed interface elements to model delamination [40-46]. Interface elements satisfy the displacement continuity conditions in a penalty sense, i.e., the stiffness of the interface element represented the penalty parameter. In contrast to continuum elements where the stress-strain relations were used, the relations between stresses and associated displacements of the interface governed the behavior of the interface elements. With the use of force and displacement components across the interface, the strain energy release rate across the interface could be calculated and the stiffness of the interface element was reduced to reflect the displacement discontinuous. With the use of interface element, the finite element modeling (preprocessing) was more convenient than the constrained conditions for the displacement continuity. Moreover, with the interface element, the continuity conditions were removed gradually, by changing the stiffness (penalty parameter) of the interface element gradually, resulting in a better representation of the damage process. 1.3 Organization of Thesis This thesis is divided into six chapters. In Chapter Two, an elasticity analysis for a simply supported [90/0] beam with different interfacial bonding conditions is presented. Two finite element models for delamination are presented. The results from finite element simulations are compared with the elasticity solution. A material constitutive model for in-plane damage of the composite laminates is presented in Chapter Three. In Chapter Four, strength of materials based delamination failure criterion is presented and the delamination shape and size are compared with the experimental results. Finite element simulations of composite laminates with a stacking sequence of [0/90/0] subjected to impact loading is presented in Chapter Five. The results from the simulations are compared with the experimental results. Chapter Six summarizes the conclusions of the present study and ends with recommendations for future studies. The Appendices follows Chapter Six and includes the computer programs developed for this thesis research. 10 2. ELASTICITY STUDIES OF DELAMINATION 2.1 Extended Pagano’s Problem In order to validate delamination simulations based on computational schemes such as finite element analysis, an elasticity solution is desired for justifying the computational approaches. Among the very few elasticity problems available for laminated composite analysis, Pagano’s problem [47] offers such an opportunity. Pagano investigated stresses and deformations of laminated composite beams simply supported at both ends and subjected to sinusoidal loading on the top surface, as shown in Fig. 2.1. Since perfect bonding conditions were assumed on the laminate interfaces, Pagano’s solutions could not be used for composite laminates with delamination. In order to account for delamination on laminate interfaces, Pagano’s problem must be modified to include shear slip and normal separation, i.e. delamination, on the laminate interfaces. The extended Pagano’s problem is explained henceforth. q(x) Z LX 90° h“ 0° h’ Fig. 2.1- A [90/0] beam subjected to sinusoidal loading on the top surface. 11 2.1.1 Fundamental Equations Similar to the original Pagano’s problem, the following equilibrium equations are used as the governing equations for the composite beams investigated in the extended Pagano’s problem. 0x,x + 022,2 = 0 02,2 +0xz,x = 0 (2-1) For plane strain condition, the constitutive relations given below are used 6.: =R110x + R120: 82 = R1202: + R220.2 7x2 = R660xz (22) where R). are the reduced compliance coefficients. Besides, the linear strain- displacement relations given below are also used for linear elastic analysis of composite beams. ex = u,x 5y = u,y (2.3) 7x2 = ux,z + uz,x 12 2.1.2 Boundary Conditions The following boundary conditions are necessary if the exact solution to the extended Pagano’s problem is sought. They are identical to those used in the Pagano’s problem. A. The condition on the top surface of the composite beams are given by h . 0. (x5) = qo sm px h an (x, 2) — 0 (2.4) where h is the total thickness, 1 is the length of the composite beams, qo is a constant and It . . p = 7. Theses two boundary condrtrons reveal that the normal stress on the top surface of the composite beams is equal to the prescribed sinusoidal loading while the shear stress vanishes. B. The bottom surface is traction free, hence both the normal stress and the shear stress vanish, i.e. 0.. (x,— g) = a... (x: g) = 0 (2.5) C. The following simply supported boundary conditions are imposed at both ends of the composite beams, i.e. x = 0 and x = I 0.40.2) =a.(1,z)=0 (2.6) uz(0,z) = uz(1,z) = o 13 2.1.3 Interfacial Conditions The deviation of the extended Pagano’s problem from the original Pagano’s problem is primarily based on the interfacial conditions, or the continuity conditions on the laminate interfaces. In the original Pagano’s problem, perfect bonding conditions are assumed on the laminate interfaces, hence both stress and displacement components are continuous through the interfaces of composite beams. However, in the extended Pagano’s problem, the laminate interfaces are allowed to be poorly bonded or completely delaminated. A. The stress continuity conditions shown below remain to be true for composite interfaces with various bonding conditions as they are essentially based on Newton’s third law. In the following equations, a local coordinate system based on the mid-plane of each composite layer is used. The thickness of each layer is represented by h with superscripts u and I denoting the upper and lower layers above and below the laminate interface of interest, respectively. hu 1 h] 0': (x9——2_) = Oz (x27) (27) u h“ 1 h’ 0x2 (LT—2”) : 0x2 (xi-é”) (2.8) When delamination takes place, the continuity of interfacial displacement components across the laminate interface will be lost. The interfacial displacement continuity equations used in the original Pagano’s problem must be modified to include slippage or mismatch on the laminate interfaces for simulations of poor bonding as well as delamination. 14 B. Linear Shear Slip Theory - A linear shear slip theory is used to correlate the interlaminar shear stress with the mismatch in the in-plane displacement between the components above and below the laminate interface, i.e. u I I u —h 1 h , h h 2 )— “x (xi?) = [Jo-x2 (L?) : 1”ng (x,-—2—) (2'9) u¥(x, When the shear slip coefficient ,u vanishes, Eq. (2.9) reduces to ll I h] —u x,-—— =0, 2 ) x( 2) u§(x, i.e. a perfect bonding condition as the in-plane displacement component on the bottom surface of the upper layer is equal to that on the top surface of the lower layer. As 11 increases, the mismatch of the in-plane displacement components between the two layers increases. For a completely delaminated composite interface, 11 becomes very large and the interlaminar shear stress on, becomes negligibly small. That it, in order to have a finite value of displacement mismatch, y should approach infinity. Accordingly, any condition between perfect bonding and complete delamination can be simulated with a value of 11 between zero and infinity. C. Linear Normal Separation Theory - In a similar fashion, imperfect normal bonding conditions can be modeled by using a linear normal separation theory. The mismatch in the out-of-plane displacement between the components above and below the laminate interface is related to the interlaminar normal stress by the following equation ”I!“ h] I I?! ha u;‘(x.—2—)—ué(x,—2—) =ko. (xiv/ca: (1%?) (2.10) 15 Again, different bonding conditions on the laminate interface can be simulated by using different values of k. That is, a laminate interfaces with a perfect bonding condition can be represented by k = 0, that with a complete delamination condition can be represented by k with an infinitely large value while that with an intermediate bonding condition can be represented by a finite k value. 2.1.4 Elasticity Solution The solution proposed by Pagano for perfectly bonded composite beams can be used to solve the extended Pagano’s problem. The following stress equations are assumed, which satisfy the governing equations and the boundary conditions. 4 . ' 2 019) = srn prAJ-imji exp(mj,-z,-) j=1 . 4 0'9) =-pzsinprAJ-i exp(mj,-z,-) (2.11) j=l . 4 0,92) = -p cos prAJ-imfl- exp(mj,-z,-) j=1 The following displacement equations are obtained from substituting the strain- displacement relations and constitutive equations into Eq. (2.11) 4 ' COS x ' ' up = ___pp 2A ji(R1('2) p2 — Rf‘l)m},.)exp(m 1,2,.) (2.12) i=1 16 . 4 . R“) ug‘) =sinprAJ-1(R1('2)mji- "122 P2)exp(mjizi) j=1 11 where superscript (i) represents for the layer number. Hence, the solution process for a two-layer composite beam is reduced to determine eight coefficients A1,, associated with the boundary conditions given by Eqs. (2.4) and (2.5) and the interfacial conditions given by Eqs. (2.7-2.10). Once the eight coefficients are identified, the stress and displacement distributions of the composite beam can be found from corresponding equations, i.e. Eqs. (2.11) and (2.12). With different interfacial shear slip coefficients, ,1: and different interfacial normal separation coefficients, k, different interfacial bonding conditions can be simulated for composite beams. Results from these elasticity studies can be used to justify delamination simulations based on computational schemes such as finite element analysis. The MATLAB program written for this elasticity study is included in Appendix A. 2.2 Finite Element Modeling of Delamination In finite element simulations, each composite laminate was modeled by multiple shell elements with one shell element representing one composite lamina. Displacement continuity conditions were enforced on the laminate interfaces between shell elements such that the composite laminate behaved like a single plate. The shell element was based first order shear deformation theory. It had four nodes and each node had six degrees of freedom. 17 2.2.1 Multi-Point Constraint Method In this thesis research, a two-lamina composite beam was investigated. Each lamina was modeled by a shell element. The following multi-point constraints were imposed on the corresponding nodes of the two shell elements to enforce the interfacial displacement continuity. “2’ =13“ —(Iz“/2)w;‘ —(h’/2)w.’. (2.13) at] =u3" —(h“/2)w;—(h’/2)w; (2.14) .1; =ug (2.15) vi =12: (2.16) vi w; (2.17) where the superscripts u and I denoted the lamina above and below the laminate interface, respectively. u x, u y and u z are displacement components in x-, y- and z-directions, respectively, y/ x , 1,11 y and l/lz are rotation components with respect to y-, x- and 2- directions, respectively. ug, u? and u? are mid-plane displacement components. The degree of freedom for rotation about 2 axis 11/ 2 was left unconstrained as its continuity did not affect the response of the beam. (The continuity conditions for the other two rotational degrees of freedom affected the interlaminar normal stress recovery and are discussed later.) A FORTRAN subroutine USER MPC for the multi-point constraints was 18 written based on these equations and was added to the commercial finite element software ABAQUS [49] for numerical simulations. The enforcement of the displacement continuity on the laminate interface warranted the laminae to behave like a single plate. The computational subroutine is given in Appendix B. Fig. 2.2 depicts the finite element model. It also shows the enforcement of the displacement continuity on the corresponding nodes. shell element at mid-plane h“ ?——0—.'——.———.‘— of upper lamina I [A A A A l\' L __________ _ _ _ _ _.. _____ .J multi-point constraints between corresponding nodes to enforce ‘— displacement continuity r __________ _ _ _ _ _ _____ ”1 (Eqs.2.14-2.18) .V v v v \ l l . . . of lower lamina I I Fig. 2.2-Two shell elements with multi-point constraints to enforce interfacial displacement continuity. Instead of the multi-point constraints, springs could be used to enforce the interfacial displacement continuity. For this, the shell elements had to be offset. The top shell element had to be offset to the bottom surface of the top lamina and the bottom shell element had to be offset to the top surface of the bottom lamina, as shown in Fig. 2.3. The springs should be able to account for the continuity conditions given by Eqs. (2.13-2.17) and the displacement continuity can be altered by reducing the stiffness of the springs to simulate conditions other than perfect bonding. 19 II; o I , . ‘ . ‘._ ‘ f ' l- ‘1 x g r ‘- p ' . ' 9. '7 n ' 1* However, such kinds of springs are not presently available in any commercial finite _ element code and must be developed by users. In this thesis research, only the multi-point constraint method was used for modeling delamination. :— ——————————————————— 1 I I i , : shell element offset from h" "_ _ _ _ _ _ _ _ Midiplflnfi _______ : mid-plane to bottom surface I I I I springs to enforce interfacial <— displacement continuity conditions I I 1 , I hI L. _ _ _ _ _ _ _ Mig'fllinf. _______ : shell element offset from : : mid-plane to top surface 1 I L .............................. J Fig. 2.3- Two shell elements with springs to enforce the interfacial displacement continuity. The interlaminar stresses can be obtained from the nodal forces involved in the constraints. Since all the finite elements used in this study were four-node rectangular elements, the interfacial area of each node can be taken as the sum of four areas surrounding the nodes in which each area is equal to one-quarter of one of the corresponding finite element. Fig. 2.4 shows the details. Fig. 2.4- Finite element mesh showing interfacial area at a node, used to calculate the interlaminar stresses from the nodal forces. 20 Thus, if Fx, Fy and F2 are the nodal forces in x-, y- and z-direction respectively, then the corresponding interlaminar stresses on, 6,2 and (32 at that node can be calculated by dividing these nodal forces by the corresponding interfacial area for the node. The details of the subroutine used for recovering interlaminar stresses are included in Appendix C. Once delamination took place, all the continuity conditions between the corresponding nodes were removed, making the finite element model to behave like separated laminae. In order to prevent possible mutual penetration between the laminae, a contact condition given below was imposed between the separated nodes. ug—ufz 20 ‘ (2'18) This condition was imposed in ABAQUS by using the one-dimensionalgap element /‘~'~--..,g . . GAPUNI‘bfgtween the nodes. The gap elements were used to connect the corresponding e \ ‘Nmm n....--"" nodes in the upper and lower layers. When the layers were perfectly bonded, all the five continuity conditions, Eqs. i(2:13-2.17), and the contact condition, Eq. (2.18), were enforced. It should be noted that when Eq. (2.15), the out-of-plane displacement continuity condition was enforced, the contact condition, Eq. (2.18), became redundant. When delamination took place, the continuity conditions were removed and the contact condition prevented penetration from happening. 2.2.2 Interfacial Layer Method This delamination model involved using an additional thin interfacial layer between the laminae. In the multi-point constraint (MPC) method, only perfect bonding and complete debonding on the laminate interfaces could be modeled as the displacement 21 continuity was either enforced or removed. In order to study a gradual transition from perfect bonding to complete debonding, the interfacial displacement had to be gradually altered. For this purpose, an interfacial layer was introduced between the two laminae. The composite laminate was modeled to have three layers as shown in Fig. 2.5. The composite laminae were modeled as shell elements and an interfacial layer with very small thickness was introduced between them. The interfacial layer had properties identical to the matrix material. Delamination was considered to take place when the interfacial layer was damaged. I I : : composite layer (shell element) multi-point constraints between 4— corresponding nodes of top lamina and interfacial layer F—-—‘F—-—-F———s-—-—-w . ._.———.-———o-———._——.-.¢— interface layer (shell element) I l <— multi-point constraints between — ———————————————————— - -_ corresponding nodes of top lamina and bottom lamina composite layer (shell element) Fig. 2.5- Composite laminate with a thin interfacial layer. The continuity conditions similar to those used in the MPC method were used to enforce the displacement continuity between the top lamina and the interfacial layer and between the interfacial layer and the bottom lamina to bond the three layers perfectly as a composite laminate. In the MPC method, all the continuity conditions were removed to simulate delamination. In the interfacial layer method, only the rotational displacement 22 V29. 1101 dis he {0? continuity conditions were removed and the translational displacement continuity conditions were left untouched; instead the elastic constants of the interfacial layer were reduced to simulate delamination. The displacement continuity conditions enforced at the corresponding nodes of the three shell elements are described henceforth. A. To enforce the displacement continuity between the top lamina and the interfacial layer, the corresponding nodes of the interfacial layer were forced to have the following continuity equations: uB'" =u2“ —(h“/2)w;‘ ~(h’"/2>w;" (2.19) ufi'" =u3" -(h“/2)w‘y‘ —(h’"/2)w;" (2.20) ug’" :4)“ (2.21) vi" =v§§ (2.22) W]? my; (2.23) With the imposition of the above continuity equations, the degrees-of-freedom of the nodes” of the interfacial layer were lost and could not be used to enforce the continuity of displacement with the bottom lamina. Thus, the bottom lamina also used the degrees-of- freedom of the nodes of the top lamina. The nodes of the bottom lamina have the following relations with those of the top lamina: 23 I u u? =14? -(h"/2)w§.‘ -(h'/2)w§-h’"w;" (2.24) ug’ = 1.3“ —(h./2)w; mil/22111; —h,..w;" (225) u3’ = 1.2“ (2.26) vi = v? (2.27) V’i = i”; (2.28) Eqs. (2.24) and (2.25) differ from Eqs. (2.13) and (2.14) by a factorh'" 111’" .The interfacial layer was assumed to have the properties of the matrix. When the composite plate was perfectly bonded and the loading was small, the rotation components Win and tug," in Eqs. (2.24) and (2.25) were small. Since hm , the thickness of the interfacial layer, was also small, the factor hm I/meecame negligibly small. Hence, Eqs. (2.24-2.28) reduced to Eqs. (2.13-2.17) and the introduction of the interfacial layer did not affect the g.-5~ I 3.. _ .‘f” . a g III-K ”£1 “,1: response of the composite beam. Based on the nodal forces involved in the continuity equations, Eqs. (224-228), the interlaminar stresses were calculated. B. For the transition from a perfectly bonded composite beam to a composite beam with poor bonding, the rotational continuity equations for wx and w y , between the top lamina and the interface layer and between the top lamina and the bottom lamina, i.e. Eqs. (2.22), (2.23), (2.27) and (2.28), were removed and the Young’s modulus and Poisson’s ratio of the matrix interface were reduced gradually. This reduction of stiffiiess 24 resulted in shear slippage or the loss of in-plane displacement continuity between the top and the bottom laminae. When the interfacial layer was damaged, the rotations of the interfacial layer were no longer small and the factor hm 111’" became significant. The Eqs. (2.24-2.28) no longer reduced to Eqs. (2.13-2.17). Instead, the displacement continuity between the top and bottom laminae was broken and a displacement mismatch equal to hm tum occurred. Hence, there was a shear slip on the laminate interface. Because the‘continuity of the out- of-plane displacement u 2 was still enforced, there was no problem of penetration. 2.3 Comparison between Elasticity and Finite Element Analysis 2.3.1 Perfectly Bonded Composite Beam The two-lamina composite beam subjected to sinusoidal loading, as shown in Fig. 2.1, was investigated. The two delamination models, the MPC model and the interfacial layer model, were compared with the elasticity solution for the composite beams with perfect bonding on the laminate interface. The normalized deflection in the mid-span of the composite beam from elasticity analysis was 2.7027. Both finite element models gave a normalized deflection of 2.6950. The stresses from elasticity analysis and finite element simulations also compared well. Hence, only the interfacial layer method was used in the later studies. Fig. 2.6 shows the normalized in-plane stress (as defined by Pagano) at the mid-span of the [90/0] composite beam as a function of thickness. Both elasticity solutions (solid line) and finite element solutions (solid square) are presented. They agree with each other very well. 25 NI 0.5 - + elasticity I FEM I I V I0.“ T I I I -700 -500 -300 -100 - 100 30 500 700 -0.5 r Fig. 2.6- Normalized in-plane stress at the mid-span through the thickness of [90/0] beam. Fig 2.7 shows the normalized interlaminar shear stress at thelaminate interface along the half span of the [90/0] beam. Stresses from elasticity solution and FEM analysis (recovered from the nodal forces as explained in section 2.2.2) are presented. Finite“ element simulations with only the rotational constraint for y/ x , with rotational constraints for both I/lx andt/Iy, and with no rotational constraints for wx and VI), gave identical : results. 26 + elasticity 4 I FEM with rotational constarint for ‘I’y A FEM with rotation contraints for ‘Px and ‘I’y 3.5 r X FEM with no rotation constraints 3 a . q " m? 2.5 r i 1 4, i, on 2 _4 1.5 i 1 i 0.5 r 0 f 0 0.1 0.2 0.3 0.4 0.5 (edge) x/l (center) Fig. 2.7- Normalized shear stress along the length of [90/0] beam. Fig 2.8 shows the normalized interlaminar normal stress at the laminate interface along the length of the [90/0] beam. Results from elasticity analysis and those recovered from finite element simulations are shown for comparison. It also shows the effect of the rotational constraints. For the simulation without rotational constraints for W1; and l/ly, the normal stresses calculated at the interface were not accurate. The simulation with only the rotational constraint for 11/ y had good agreement with the elasticity solution. However, the simulation with the rotational constraints for both 1,1! x and W}, yielded small discrepancy with the elasticity solution. In the finite element analysis, the transverse stresses could only be recovered at the laminate interface. Hence, the distribution of 27 transverse shear stress (5x2 and transverse normal stress (52 through the thickness of the beam were not presented. + elasticity 1-6 7 I FEM with rotational constraint for \yy A FEM with rotation constraints for ‘Px and ‘l’y 1 4 4 X FEM with no rotation constraints 0 X X X 1.2 —« X X 0, 1 If x 0.8 4 X 0.6 - X 0.4 ~ 0.2 4 0 I I I 0 0.1 0.2 0.3 0.4 0.5 (edge) x/l (center) Fig. 2.8- Normalized normal stress along the length of [90/0] beam. 2.3.2 Imperfectly Bonded Composite Beam Finite element analysis with an interfacial layer, presented in section 2.2.2, was used to simulate different bonding conditions at the interface of the two-lamina composite beam. The interfacial shear slip was introduced by the factorhm y/m , i.e. the multiplication of thickness of the interfacial layer hm and the rotation component 111’" . When the elastic modulus (E) of the interfacial layer was reduced, the factorh'" y/m increased, resulting in a poor shear-bonding condition. In the elasticity solution, the factor ,u was increased to simulate the poor shear-bonding condition. There was no 28 direct correlation between the two factors, E and 11 that introduced the shear slip in the finite element modeling and elasticity analysis, respectively. However, both of them could introduce the shear slip gradually. Hence, they were compared to study the behavior of the [90/0] beam with different shear-bonding conditions, i.e. p was changed from 10'7 to unity and E was changed fromlO7 to unity. Hence, log ,1: and -log B were chosen as the scales for comparison. Fig. 2.9 shows the interfacial shear slip (mismatch of the displacement along the beam axis) for different values of ,u and E. 3.00E—03 - +elasticity l FEM 2.50E-03 i 2.00E-03 i I 1.50E-03 - 1.00E-03 ~ displacement mismatch, mm 5.00E-04 r 0.00E+00 I | I I -7 -6 -5 .4 -3 -2 -1 0 log it or -log E Fig. 2.9- Displacement mismatch at the interface at different bonding conditions. Different interfacial shear slip conditions simulated different bonding conditions on the laminate interface. It was found that the maximum deflection of the composite beam at very highp and k=0, i.e. poor shear bonding but perfect normal bonding, was 1.44 times the deflection of a perfectly bonded beam (,u=0, k=0). Fig. 2.10 shows the normalized deflection at the mid-span of the [90/0] composite beam for different shear 29 bonding conditions obtained from elasticity analysis and finite element simulations. Good comparisons can be observed between them. Fig. 2.11 shows the normalized deflections at the mid-span of the composite beam at different normal separation and shear slip conditions. The results were obtained from the elasticity analysis. 5 n 4 _ _ 3 - +elasticity W ‘ F I FEM 2 4 1 _ 0 T I I T I f I -7 -6 -5 -4 -3 -2 -l 0 log 11 or - log E Fig. 2.10 - Normalized deflections of [90/0] beam for different shear bonding conditions. The deflection of the composite beam with very high values of ,u and k, i.e. with poor shear bonding and poor normal bonding, was found to be 36.47 times that of a perfectly bonded beam. The finite element simulation gave a ratio of 36.48. This result clearly - indicates the importance of delamination modeling in the analysis of composite laminates. It may also be pointed out that the deflection of each bonding condition converged to a constant. The interlaminar shear stresses at the beam end for different shear bonding conditions obtained from the elasticity analysis and the finite element simulation are shown in Fig.2.12. Because of the symmetry of the problem, the interlaminar shear stress distribution was symmetric on either side of the mid-span of the 30 beam. As shown in the diagram, the interlaminar shear stress gradually decreased and reduced to zero at very poor bonding condition. 100 so 60 W 40 20 J -—< .4 +k=0.01 + k=0.l +k=l +k=10 +k=100000 Fig. 2.11- Normalized deflections of [90/0] beam for different bonding conditions. + elasticity 4 I FEM , T I j -7 -6 -5 -4 -3 -2 -1 0 log it or -log E Fig. 2.12 - Interlaminar shear stress at the edge for different shear bonding conditions. 31 The interlaminar normal stresses at the center of the beam for different shear bonding conditions obtained from the elasticity analysis and the finite element simulation are shown in F ig.2.13. As shown in the diagram, the interlaminar normal stress increases as the shear bonding decreases. 1 I A 3 I I I t—O/‘S/: 0.8 _ 0-5 ‘ +elasticity oz I FEM 0.4 ~ 0.2 a o I I I I I I I I -7 -6 -5 -4 -3 -2 -1 0 1 log p or -log E Fig. 2.13 - Interlaminar normal stress at the center for different shear bonding conditions. 2.3.3 Further Elasticity Analysis The through-the-thickness distributions of transverse shear stresses at the beam end WW“... for different shear slip coefficients and perfect normal bonding conditions are shown in Fig. 2.14. As the shear slip coefficient increased, the shear bonding became poorer and the interlaminar shear stress became smaller. At high shear slip coefficients, there was no interlaminar shear stress on the laminate interface. 32 N | . 0.5 ~\ k=0 0375 - \ +u=0 +p=l.00E-03 -\ 0.25 — +u=10 ; -. i ' ”1’13 It 1 U 0.125 4 04 ~--— \M 4 i . on It s ' ~ - i 15 20 -0.125 - \\ . -0.25 - '0375 T / -0.5 ‘1 I '4 . .v '1 ‘ . , ”"16“, yr“ Fig. 2.14- Transverse shear stress with different shear slip coefficients and perfect normal bonding. Fig. 2.15 shows the through-the-thickness distribution of transverse shear stress at the beam end for a poor shear slip coefficient ( ,u =1000) with different normal separation coefficients. The interlaminar shear stress was always continuous for all the bonding conditions although the value was zero due to the poor shear bonding condition. It can also be observed from Fig. 2.15, in a totally delaminated beam, (with high shear slip and normal separation coefficients) the bottom layer had no response to the tensile loading applied on the top lamina. The two laminae seemed to behave completely independent. From Fig. 2.14 and 2.15, it could be concluded that the transverse shear stress was more affected by shear slip coefficient ( p ) than by the normal separation coefficient (k) 33 0.375 r 0.25 r 0.125 — 0 1 0x1 25 -0.125 -0.25 +k=0 “0.375 +k=l +k=1000 '05 11:1000 Fig. 2.15- Transverse shear stress at the edge for large shear slip coefficient and different normal separation coefficients. The interlaminar shear stress distributions along the length of the beam for different shear bonding conditions and perfect normal bonding condition are shown in Fig. 2.16. The through-the-thickness distributions of normalized transverse normal stress at the mid-span of the beam are plotted for different bonding conditions. Fig. 2.17 shows the cases of different normal separation coefficients with a perfect shear bonding condition. As the normal separation coefficient increased, the normal bonding condition became poor and the interlaminar normal stress became smaller. 34 =0 + u=0 -—I-- p=0.001 + p=10 0 0.1 0.2 0.3 0.4 0.5 (edge) 11/] (center) Fig. 2.16- lnterlaminar shear stress distribution along the length of the beam for different shear slip coefficients with perfect normal bonding. i 0.5 — +k=0 ~1I— k=0.01 +k=10 -05 iF-‘O Fig. 2.17- Normalized transverse normal stress distribution for different normal separation coefficients with no shear slip. 35 Fig. 2.18 shows the cases of a high shear slip coefficient with different normal separation conditions. From Fig.2.l7 and 2.18, it could be concluded that the transverse normal stress was more affected by the normal separation coefficient k than by the shear slip coefficient ,u . i 0.5 T 0.375 - i/ -/ /,.Q 0.25 ~ ,I ,/ ,fi 0.125 - /,,n/ a, 0 m/ _ I 0'z 0.25 [In 0.5 0.75 1 -0.125 /,n’ /" /Dl +k=l '0375 +k=1000 -0.5 ' 11=1000 Fig. 2.18- Normalized transverse normal stress distribution for different normal separation coefficients with large shear slip. The interlaminar normal stress along the length of the beam for different normal separation conditions and perfect shear bonding is shown in Fig. 2.19. As the normal separation coefficient increased, the normal bonding became poorer and the interlaminar normal stress became very smaller. 36 0.9 1 0 0.1 0.2 0.3 0.4 0.5 (edge) x/l (center) Fig. 2.19- Normalized interlaminar normal stress along the length of the beam for different normal separation coefficients with perfect shear bonding. 37 3. In-Plane Failure A wide range of constitutive models have been proposed in the literature for laminated composite materials. Most of the models assume a linear elastic response till the point of failure. They are, however, different in failure criteria for determining the point of failure and in progressive damage models for describing the post-failure behavior. Commonly used failure criteria are Tsai-Hill’s failure criteria, Tsai-Wu’s failure criteria, Hashin’s failure criteria, Chang’s failure criteria, etc., while examples of post-failure models include brittle failure, linear softening, exponential softening, softening based on Weibull’s distribution, etc. A constitutive model is a phenomenological description of material behavior. In establishing a constitutive model, it is necessary to characterize the behavior of the composite material of interest by testing methods. Once the mechanical responses and damage modes of the composite material are identified, mathematical formulations can be established or selected to simulate the composite behavior. The discussions henceforth describe the composite responses and the incorporation of damage modes into a constitutive model. 3.1 Basic Assumptions The following assumptions were made in developing a constitutive model for a composite lamina. A. A composite lamina was assumed to be homogeneous. Its elastic constants were averaged properties obtained from material tests. 38 B. Each composite lamina was modeled by a two-dimensional orthotmpic continuum. The orthotropic nature of the composite material was maintained throughout the damage process. C. Although the stress-strain curves for a composite lamina showed some degrees of nonlinearity, especially in shear, linear elastic relations were assumed till the point of failure. .‘ TD\.)The damage process of a composite lamina was considered to be linearly brittle. ‘. ‘‘‘‘‘‘‘ This was due to the fact that very small permanent deformations remained in the composite lamina afier unloading, as long as the lamina was not very close to complete failure. 3.2 Damage Modes and Stiffness Degradations As mentioned above, a two-dimensional linear elastic orthotropic model was assumed to represent the behavior of a composite lamina till the point of failure. From Hooke’s law, the following stress-strain relation was given - I r0'” Ell EZZVIZ 0 £11 I " V12V21 1 " V12V21 022 = EIIVZI E22 0 £22 (3_1) 1 ’ V12V21 1 ‘ V12V21 012 0 0 612 712 where a is stress, a and y are strains, E is Young’s modulus, v is Poisson’s ratio, G is shear modulus and the subscripts 1 and 2 represent the directions of fiber and matrix, 39 respectively. When a composite lamina undergoes damage, the damage may manifest itself in different forms. They can be divided into fiber damage and matrix damage. 3.2.1 Fiber Damage The stress in the fiber direction of a composite lamina is predominantly carried by the fibers because of their high stiffness in comparison with matrix material. Thus, the stress bearing capability in the fiber direction is not significantly affected by the damage in matrix. Fibers may fail because of tensile breakage. Broken fibers may debond from matrix and cause cavities. Fibers may also buckle or kink due to compressive loading. These fiber failures are catastrophic and can strongly affect the stress bearing capability of the composite lamina to which the fibers belong. The following failure criteria proposed by Hashin [3-4] were used to predict the fiber failure in this thesis research. a. fiber failure in tension: 0'll 20 2 0'11 2 2 0 failed : -— —l 3.2 ef ( X t ) {< O elastic ( ) b. fiber failure in compression: 0"ll < O efz :(ZH)2_1 C {2 0 failed (3 3) < O elastic where X t and X C are lamina strengths of the composite lamina along the fiber direction in tension and in compression, respectively. Fiber damage could also affect the integrity of matrix and hence all stiffness components in the constitutive model should be 40 degraded when fiber damage takes place. As an example, Fig. 3.1 shows a uniaxial loading-unloading-reloading curve. As can be seen, it is an elastic-brittle response with a linear stress-strain relation till the point of failure. At the point of failure, the stress value drops to a low level and the composite lamina continues to maintain the same stress level on further loading. If unloading takes place, the material unloads with a secant modulus. When the stress level becomes zero, there is no permanent deformation in the composite lamina. loading 0 ‘ > +7 unloading/reloading 8 Fig. 3.1- Uniaxial stress-strain curve along the fiber direction. When a composite lamina is damaged, the stiffness matrix given in Eq. (3.1) took the following form: — — * t * E11 E221’12 0 * * * * 1-V12V21 1'V13V21 * * EIIVZI E22 0 (3.4) * * ¥ * l—V12V21 1-V12V21 a: 0 0 G12 41 where the quantities in * represent the corresponding degraded elastic constants. They were defined as follows. E11.1:0'ill/“311 Viz =(V12 /511)Ei1 E52 = 1322/1000 (3.5) V21: Viz/522 /Ei1 Gf'z = 012/1000 It should be noted that of] :01? /5 where ofl is the stress to cause fiber failure. The degradation factor 1/ 1000 was chosen for the property degradation because if the property was degraded further, it gave convergence problems, while a degradation factor greater than 1/ 1000 would result in large post failure stress at large strains. For the degradation of E 1 1, a degradation factor US was chosen as it gave good match with the experimental results. Degradation factors of 1/3, 1/ 10 and 1/ 100 were also investigated. The degradation factors 1/3 and 1/10 gave almost similar results as the factor US but for the degradation factor of 1/ 100, the damage was too severe and the impactor penetrated the laminate at much lesser energy than in the experiments. 42 3.2.2 Matrix Damage Shear stress and transverse normal stress can transfer to both fibers and matrix. But their effects on composite damage are predominantly restricted to matrix damage and fiber-matrix debonding. In a composite lamina, an initial crack in the matrix can propagate in any direction. However, when the crack reaches a fiber-matrix interface, it changes its course to propagate along the fiber direction without crossing into the fiber. Thus, the damage in matrix is aligned in the fiber direction. This damage of matrix depends upon the matrix strength in both normal direction and shearing direction, if debonding between fiber and matrix is not of concern. Being different from matrix cracking that is observed due to transverse tensile stress and shear stress, matrix crushing is observed due to transverse compressive stress. A number of failure criteria were proposed for matrix damage. In this thesis research, failure criteria proposed by Hashin and modified by Chang and Lessard were used. It contained the following basic items. A. matrix cracking due to shear stress and transverse tensile stress: 0'22 2 0 2 022 2 012 2 ZOfailed = — + — ‘1 3.6 em ( Yr ) ( S ) {<0elastic ( ) where Y, and S are the transverse tensile strength and shear strength of the composite lamina, respectively. The elastic constants were degraded to 0.1% of their undamaged values to simulate the brittle failure. The damaged stiffness matrix was defined as 43 En ESZVIZ I a O l-v12v21 I-VIZVZI El 1V2: 1532 . 0 (3.7) 1—v12v21 l—v12v21 0 0 of} _ i where 1532 = 522/1000 * * V12 =(V12 ”510511 and of} = 612/1000 B. matrix crushing due to transverse compressive stress: 0'22 <0 0’ a Z Ofailed em2 = (#12 Hi? —1 I . (3.8) YC S < 0 elastic where Yc is the transverse compressive strength of the composite lamina. The stiffness matrix for the damaged composite was modified as E1 1 ESZVIZ 0 l—vlzvgl 1-v12v31 5116: 1532 . 0 (3.9) 1—v12v21 1’V12V21 0 0 G12 _ i 44 The E32 and v3, are defined in the same manner as in matrix cracking, Eqn. (3.7). The constitutive model consisting of the failure criteria and degraded stiffness matrices given above was incorporated into the commercial finite element code ABAQUS by programming it into a user’s material subroutine UMAT for the analysis of the composite material. The computer subroutine is given in Appendix E. As nonlinear equations are solved in increments, the solution process was incremental. For the 11m increment, the strain increment for the nth increment and the stresses and strains of the n-lth increment were passed on to the subroutine by the software. The solution of nonlinear equations requires the tangent stiffness matrix which depends upon the material constitutive relations. Thus in the subroutine UMAT, the material constitutive relations to calculate the stresses at the end of the nth increment and the tangent stiffness matrix for the nth increment were defined. 3.3 Validation of Dynamic Analysis The commercial finite element package ABAQUS was used in this thesis research. In order to validate the contact algorithm and the dynamic analysis procedures of the finite element code, a plate impact problem studied by Chen and Sun [50] was re-examined. The force and deflection histories were of the primary interests. In this study, a 200 mm x 200 mm graphite/epoxy plate with a stacking sequence of [90/0/90/0/90]s was impacted by a 12.7 mm diameter and 8.537 g spherical projectile at an impact velocity of 3 m/s. The material properties for the graphite/epoxy composite are listed below: 1311:1412 GPa, E22=9.72 GPa, V12=0.3, V23=0.3, G|2=5.53 GPa and G23=3.74 GPa. 45 The thickness of each composite ply was 0.269 mm and the plate was assumed to be simply supported along all four edges. In the finite element analysis, a multiple-shell element model (MPC model) as presented in Chapter 2 was used. Fig.3.2 (a) shows the force history and Fig. 3.2 (b) shows the displacement history obtained fiom the finite element analysis and experiments of Ref. [50]. As can be seen, both the displacement and force histories compare well with the corresponding experimental results. 350 _ Expenment — - — Simulation 300 - 250 - 200 4 force, N 150] 100 4 50~ time, msec (a) 46 0.2 “ Experiment 0.18 - 0.16 — 0.14 - x - 0.12 ~ ‘ \ - - - Simulation tmm / P -h 1 0.08 ~ 0.06 ~ 0.04 ~ 0.02 a o I I I I 1 0 1 00 200 300 400 500 time, msec displacemn (b) Fig. 3.2- Force and deflection history from experiment and simulation. 3.4 Simulation of [0/0] Laminates Having established the finite element procedure, the next step was to investigate the impact tests for unidirectional composite laminates with a stacking sequence of [0/0]. It should be noted that there was no delamination failure in these composite laminates as they had identical material properties through thickness. The studies were performed on a glass/epoxy composite material. The composite laminates had dimensions of 125mm x 125mm and a thickness of 3.8lmm. All four edges were fully constrained. The specimens were impacted with a spherical impactor of 12.88 kg mass and 12.7mm diameter at different impact velocities. The material properties used in the studies are listed below: E1 1:390 GPa, E22 =8.6 GPa, v12 =0.28, G12 =3.8 GPa, G13 =3.8 GPa, G23 =2.9 GPa 47 The material strengths used in the simulations are given as follows: Xt= 1080 MPa, X¢= 620 MPa, Yi= 39 MPa, Yc= 128 MPa, S= 71 MPa In the finite element simulation, the impactor was modeled as a rigid body (analytical rigid body in ABAQUS) and the shell element (S4R) of ABAQUS was used for modeling the composite laminates. 3.4.1 Matrix Failure Only As a beginning simulation, a material constitutive model with only matrix failure was assumed, i.e. fiber failure was not allowed to take place. The force-deflection curves for different impact velocities are shown in Fig. 3.3. A peak force around 8 kN was obtained for the highest velocity, though experimental results showed a peak force of 4.6 kN. The model based on matrix failure only seemed to be too rigid. Inclusion of fiber failure in the constitutive model was deemed necessary. 9 a 8 _ 7 . —v=1.7 m/s -~ v=2.0 m/s 5 * -—~——v=2.2 m/s 5 5 _ / -—v=2.3 m/s 3' , / —v=2.5 m/s E 4 i /,/ —v=2.9 m/s 3 a J / 2 - ,x/ J 1 - /,- o :I/ v i 1 l I T r I 1 Y r T I 1 0 5 10 15 Deflection, mm Fig. 3.3- Force-deflection curves for different impact velocities based on material model with matrix failure only. 48 3.4.2 Brittle Fiber Failure In addition to the matrix failure, a brittle failure model was used for simulating fiber failure. There was no strain softening after the point of failure. When the failure was reached, the elastic constants were reduced to 0.1% of their undamaged levels. It was observed from the simulations that the material stiffness degradation was too severe and the impactor penetrated the specimen at a velocity much lower than that observed in the experiments as evident from Fig.3.4. 5- z X 111“ § "' — Exp. v=2.02 m/s ---- Sim. v=2.02 m/s 0 ‘I I j 0 5 10 15 deflection, mm Fig. 3.4— force-deflection curve based on a material with brittle fiber failure. 3.4.3 Progressive Fiber Failure A progressive fiber failure model was used in conjunction with the matrix failure model. The stress-strain curve is shown in Fig. 3.1. Upon failure, the stiffness in the fiber direction is degraded such that the lamina can only support only 20% of the initial strength. The stress is maintained on further loading. Upon unloading the lamina unloads 49 with the secant modulus. The force-deflection curves resemble that obtained from experiments up to some extent. Fig. 3.5 shows the entire force-deflection curves obtained from simulations at different impact velocities. Figs. 3.6(a), 3.6(b), and 3.6(c) show the comparison of force-deflection curves between experiments and finite element simulations. ———v=1 .2 m/s -—--v=1.7 m/s —v=2.02 m/s —v=2.2 m/s ——-v=2.3 m/s -- ~-~v=2.5 m/s ——v=2.9 m/s Force, KN 4 8 12 16 20 Deflection, mm Fig. 3.5- Force—deflection curves for different impact velocities based on progressive fiber failure model. 50 — Exp. v=1.7 mls Deflection, mm (b) 51 z x -—-— Sim. v=1.7 mls e“ 2 0 IL 0 I I I W I O 8 12 16 20 Deflection, mm (a) 5 . |/rv z 3 r -— Exp. v=2.3 mls x " A ,J 7 --—-Sim. v=2.3 mls r ,7 u? 2 - /1‘V 1 1 // /. - // o 1 I I T l 0 8 12 16 20 Force, KN Fig. 3.6- Comparison of force-deflection curves at three different impact velocities. Fig 3.7 shows the absorbed energy vs. the impact energy from both experimental and simulation results. The impact energy is based on the kinetic energy of the impactor while the absorbed energy is obtained from the integration of the force-deflection curves. As can be observed, the absorbed energy is smaller than the impact energy for small impact energies. The absorbed energy increases with the increase of impact energy. This trend is observed till a point where both the absorbed energy and the impact energy are about equal to each other for the first time. This point is called the point of penetration. The absorbed energy is equal to the impact energy from this point of penetration to the point of perforation (where the absorbed energy and the impact energy are about equal to each other for the last time). Upon perforation, the absorbed energy again becomes ff I I // | \ -/ —-— Exp. v=2.5 mls / ————- Sim. v=2.5 mls 8 12 16 20 Deflection, mm (C) smaller than the impact energy. 52 8 + Experiment -I-— Simulation ~~ Equal Energy Line Absorbed Energy ,J on O o I T I I j ‘1 0 1O 20 30 40 50 60 Impact Energy ,J Fig. 3.7- Absorbed energy vs. impact energy With these simulations, it was concluded that the material constitutive model was able to capture the damage model to a reasonable extent. A material model with linear or exponential strain softening could capture the post-failure behavior in a better way; however, all the attempts to develop such a material model were unsuccessful owing to convergence problems encountered in the finite element simulations. The next step for finite element simulations of composite laminates subject to impact loading should be to incorporate delamination modeling into the computational scheme to capture the out-of- plane damage through the composite thickness. 53 4. Delamination Failure Composite laminates are heterogeneous and anisotropic. Their properties vary from one constituent to another within one lamina and change from one lamina to another through the laminate thickness. Because of the mismatch in material properties, the nonuniform stress distribution through the laminate thickness can cause delamination in composite laminates. As an example, a [0/90] laminate has the highest mismatch of material properties among the [0/0] laminates, and is the one that is most susceptible to delamination. Delamination is a major form of damage as it severely degrades the properties of a composite laminate and renders the composite laminate prone to further damage. In fact, once delamination takes place in a composite laminate subjected to impact loading, it is very easy for the impactor to run through the composite laminate by cracking the matrix and pushing away the fibers in individual laminae. Thus, delamination is an important characteristic of composite damage and composite simulation should include delamination modeling. 4.1 Delamination Modeling and Criteria In Chapter 2, two delamination models were presented, namely the multi-point constraint model and the interfacial layer model. In principle, both were based on using multi-point constraints to enforce the continuity of displacement across the laminate interfaces and to calculate the interlaminar stresses using the constraint forces. Hence, both models gave similar delamination simulations. However, the two models differed from each other on how the continuity of displacements across the laminate interfaces was removed when the criterion for delamination was met. In the multi-point constraint model, the both rotational and translational displacement continuity conditions at the 54 nodes were removed in order to allow displacement discontinuity. In the interfacial layer model, the stiffness of the interfacial element was reduced in order to introduce the displacements discontinuity. Since an additional shell element layer was used as the interfacial layer, the interfacial layer model was computationally more expensive than the multi-point constraint model. Thus, the multi-point constraint mmidwwas...used-in--this"""" study. a” A number of delamination models can be found in the literature. These criteria can be divided into two groups, the fracture-mechanics based models and the strength-of- materials based models. In the fracture-mechanics based models, an initial delamination front must be assumed along with a delamination failure criterion to calculate the onset of delamination. The propagation of delamination front, for example, can be determined by the strain energy release rate obtained fi'om the virtual crack closure method i.e., if the strain energy release rate exceeds a critical level, the delamination propagation will take place. In the strength-of-materials based failure models, it is assumed that delamination takes place due to the presence of large interlaminar stresses. Failure criterion based on interlaminar stresses and interlaminar strengths are used to predict the delamination onset. In this thesis research, strength-of—materials based models were used, [As-[explained in Chapter 2, the interlaminar stresses can be calculated accurately with the proposed delamination modeling. To begin with, a quadratic failure criterion was used to predict the onset of delamination. This criterion has the following form: For 0'33 2 O: 0' 0' o- 2 Odelamination e2=(_1_3_)2+(_23_)2+(_3_3_)2_1 { (4.1) d S S Y, < 0 no delamination 55 where S is the interlaminar shear strength and Y, is the interlaminar transverse strength in tension. It is assumed that delamination takes place only if the interlaminar normal stress . A.n-..,‘~' .rfi-n. 033 is in tension. This criterion was used by many researchers [25-34] and was also used —,. ~_ in this thesis research for simulations of composite laminates subjected to impact loading. 4.1.1 Delamination Failure Only A 125 mm x 125 mm composite laminate with a stacking sequence of [0/90/0] and thickness of 3.81 mm was impacted with an impactor of 12.88 kg mass at a velocity of 3.0 m/s. The composite laminate was_rno.deled_by _tw0§11§]1”§I§IBQD.t§.-..The top shell element represented the 0/90 laminae while the bottom shell element the bottom 0 lamina. Thus, delamination was simulated to take place on the interface between the 0/90 lamina and the bottom 0 lamina. The use of a single interface in delamination simulations is due to the limitation of the commercial finite element code used in this thesis research. _...- ”4.4., -- For a laminate with two interfaces, three shell elements will have to 66115621; i.e. each shell element represents a lamina. However, the process of assembling the three laminae to become a composite laminate will render all the nodes but one in the laminae to become dependent variables. The nodes of the top lamina must be independent for the contact algorithm to work in ABAQUS Standard. Thus, the displacement continuity conditions are applied to the corresponding nodes of the middle lamina and makes the degrees of freedom of the nodes of the middle lamina dependent upon the nodes of the top lamina. Hence, the nodes of the middle lamina are unavailable for imposing the displacement continuity conditions between the middle and the bottom laminae. The nodes of the top lamina will have to be used to enforce the continuity conditions. Accordingly the nodal forces in the bottom lamina will have the combined effect of the 56 top and the middle laminae and thus, cannot be used to calculate the interlaminar stresses at the bottom interface. Besides, it has been noted that the bottom interface has the major delamination damage when subjected to impact loading. It was with these reasons, only two shell elements were used to model the three-lamina composite laminate. When the delamination failure criterion, Eq. (4.1), was used for modeling composite laminates without an in-plane failure criterion, it predicted delamination shape well, as shown in Fig. 4.1. The black elements in the center of the composite laminate shown in Fig. 4.1 (a), collectively represent the delamination area. They were identified from the delaminated nodes shown in Fig. 4.1 (b), to which they attached. As can be observed from Fig. 4.1, the typical peanut-shaped delamination is observed with the major axis aligned in the fiber direction of the bottom 00 lamina. (a) (b) Fig. 4.1 - Delamination at the bottom interface of l0/90/0] laminate with no in-plane failure. 57 4.1.2 In-Plane Failure and Delamination Failure When the delamination failure criterion was used in conjunction with an in-plane failure criterion, it was observed that the delamination size was greatly reduced as shown in Fig. 4.2. This was contrary to what was expected. The in-plane damage modes of the composite material such as matrix damage should have led to increased delamination. In fact, it was observed that the interlaminar shear stresses decreased and the interlaminar normal stress increased when the material stiffness was degraded upon matrix failure. The combination of the changes of interlaminar stresses resulted in the reduction of delamination area. Finite element simulations revealed that the delamination area based on the delamination failure criterion and an in-plane failure criterion proposed in section 3.43 was found to be only 2.73 cm2 while that based on the delamination failure criterion only was 4.29 cmz. I IF] Fig. 4.2 - Delamination at the bottom interface of [0/90/0] laminate with in-plane failure criterion. 58 4.1.3 Modified Delamination Failure It became clear that the foregoing delamination failure criterion must be modified. Fenske and Vizzini [28] proposed a delamination failure criterion with in-plane stress components. They claimed that delamination usually took place in thin matrix rich layers on laminate interfaces. Hence, all six stress components should be included in the justification of delamination onset. The foregoing studies seemed to support their theory. Thus, the delamination failure criterion was modified to include the in-plane stresses in this thesis research. The original delamination criterion proposed by Fenske and Vizzini takes the following form: 2 2 0’ —0'0' +0 Y —Y 0' +0 0 0' 82: I I II [[+(c i)( I 11)+(13)2+(23)2 ‘1 131’, 13x, 5 s 0i; + (Zc - Z, )0'33 _ I Z 0 delaminaton (42) Z ,Z 0 Z ,Z G <0 no delaminatbn where 0', and a” are in-plane principal stresses, Y, and Y C are matrix strengths in tension and compression, respectively, S is interlaminar shear strength and Z, and Zc are interlaminar normal strengths in tension and in compression, respectively. It should be noted that both in-plane stresses and interlaminar stresses are included in the delamination failure criterion. However, Eq. (4.2), the distinction between the tensile interlaminar normal stress and the compressive interlaminar normal stress is not very strong. Although it is well known that the compressive interlaminar normal stress can act against delamination. The fact that delamination has a peanut shape consisting two discontinuous lobes upon impact seems to support this argument. Hence, the two 59 delamination lobes may not be discontinuous at the point of impact if the delamination failure criterion does not distinguish between interlaminar normal compression and interlaminar normal tension sufficiently. The Fenske-Vizzini delaminaiton failure criterion is modified to enhance the role of the compressive interlaminar normal stress. It is proposed that delamination takes place only when the interlaminar normal stress is in tension. Thus, the delamination criterion takes the following form: For 033 20: 2 2 0 —00 +0 Y —Y) 0 +0 ) 0 0 65:1 111 11+(c 1(1 11+,13)2+(23)2 141', 131/, s S (4.3) 1 . . 0 20d] t +(_§i)2_1{ eammazon Y, < 0 no delamination It may be noted that Eq. (4.3) is an extension of Eq. (4.1), with additional terms accounting for in-plane stresses. To calculate the in-plane stress on the laminate interface with the multi-point constraint model, an additional layer of matrix must be introduced in the finite element analysis. Thus, the [0/90/0] laminate which was modeled by two shell elements, the top element representing the 0/90 lamina while the bottom element the bottom 0° lamina, is remodeled as [0/90/M/0]. The top shell element represented the 0/90/M lamina and the bottom element the bottom 0 lamina. The additional matrix layer was very thin with a thickness of 0.01mm and did not have a significant effect on the global response of the composite laminate. With this additional matrix interfacial layer, it was possible to obtain the in-plane stresses on the laminate interface for delamination analysis. 60 When the impact simulation similar to the previous cases was performed with the modified delamination failure criterion, a delamination area of 7.03 cm2 was obtained. "pr3“? I!” thfizmxbwfln ’ I Fig. 4.3- Delamination at the bottom interface of [0/90/0] laminate based on modified delamination failure criterion. 4.2 Investigation of Delamination Both delamination shape and size from finite element simulations require validation. In these studies, both static loading and impact loading were used for delamination validations. 4.2.1 Delamination Shape Composite laminates with dimensions of 125 mm x 125 mm and thickness of 2.6 mm were loaded by a point force at the laminate center. The material in-plane damage was ignored for the simulations. Stacking sequences of [90/0] and [0/0] where 0=15, 30, 45, 60, 75 and 90, were investigated. The delamination shapes at the interfaces of [90/0], [0/90], [0/45] and [0/ 15] are shown in Figs. 4.4(a), 4.4(b), 4.4(c) and 4.4(d). As can be 61 seen from the diagrams, typical peanut-shaped delaminations are obtained with their major axises aligned in the fiber direction of the lamina below the laminate interfaces. ' A ' a. II I ‘II ' 1- II . - Anni..- . All lug. (a) [90/0] (b) [0/90] "LA.-‘i '?'~. allaflllla. .1 .‘H.I IIBHUL ‘~ lillfllil lilflfi ml”..- II!” All Ifiw'qw (6) [0/451 (d) [0/151 Fig. 4.4- Axis of delamination rotates with the fiber orientation of the bottom lamina. 62 4.4.2 Delamination Size The delamination sizes obtained from the simulations were compared with the experimental results. Impact tests were performed on a [0/90/0] composite laminate with dimensions of 125 mm x 125 mm and thickness of 3.81mm. A cylindrical impactor of 12.7 mm diameter and 24 g mass was used. Composite laminates were fixed at all four edges and were impacted at various velocities. The impact tests were considered to be low mass but high velocity with impact energy levels of 13 J, 17 J and 26 J. The simulations did not take into account the material in-plane damage. Figure 4.5 shows the delamination area versus the impact energy from both experiments and finite element simulations. Figs. 4.6(a), 4.6(b) and 4.6(c) show the delamination area from the simulations. It clearly shows that the modified delamination failure criterion is able to capture the delamination shape, though the delamination sizes are not very accurate. Errors of 56.6 %, 55% and 38.6%, are found in the. three simulations. Further improvement of delamination failure criterion is required. 25 ~ ,5, 20 4 0 Experiment 0 g I Simulation 5 15 - .5 o I .E g 10 1 . o o a g 5 , a I- 0 i i 1r a 1 i O 10 20 30 Impact Energy, J Fig. 4 5 - Delamination area vs. impact energy. 63 I=Ilnll¥i .- I. ISI- 5 “QB EIRREWIIIIIIIIIHUUERIIIR fr amnxaaan ‘EAIIDE? "finfllflflflfl (a) 13 J (b) 17 J MINI 'Illlllfi I... II llfllhfi =- II. Cl llglnllflw 'IIDIDRUHDKMEEHBSRI lllbé" null. IHVK‘HW 1:2MRSF“ vauuolll: :l Imxyu» a} «I Il.lll::wnfl::::lflllh '? I IIHIIBII:II Ifllllfla‘l IDII ll I'l .- Ilgfll (c) 26 J Fig. 4.6 - Delamination at the bottom interface of [0/90/0] laminate under various impact energies. Another validation on delamination size was based [O/O] laminates subjected to static point loading. The result is shown in Fig. 4.7. It can be seen that the delamination area increases steadily with the increase in the mismatch of fiber angle between laminae. 3 i 254 0 ° 0 "e 0 § 2‘ ° < j 9 S15 :3 o '5 J E 1 L“ 8 0.5~ o T T T T T 1 0 15 30 45 60 75 90 Angle, 0 Fig. 4.7 - Delamination area vs. fiber angle between laminae. 65 5. Simulations of [0/90/0] Laminates Having proposed a constitutive model to express in-plane damage (including fiber failure and matrix damage) and out-of-plane damage (such as delamination damage), composite laminates with a stacking sequence of [0/90/0] were studied. Because of the large mismatch in material properties through the laminate thickness, delamination was expected to be severe for the cross ply laminate. Once again, composite laminates with dimensions of 125 mm x 125 mm and thickness of 3.81 mm were investigated. Each of them was completely constrained at all four edges and was impacted with a spherical impactor of 12.7 mm diameter and 12.88 kg mass at different impact velocities. For the finite element simulations, a two-shell-element model was used. As the delamination model could handle only one interface, the delamination on the top ‘ h-h‘v -. interface of the 0/96wl‘aiiiinawwas neglected, i.e. the 0/90 laminae were modeled by one shell element. The second shell element modeled the bottom 0 lamina. The composite constitutive model of progressive fiber damage and brittle matrix damage presented in Chapter 3 was used to model the in-plane damage of the composite laminates. The delamination modeling presented in section 4.1.3 was used for modeling the out-of-plane damage. The impactor was modeled as an analytical rigid body and the shell element S4R of ABAQUS was used for the composite lamianae. Fig. 5.1 shows the force-deflection curves from the impact simulations for different impact velocities. There were some convergence problems in the simulations at high impact velocities. The simulations stopped near the maximum load. At High impact velocities, the simulations were unable to capture the penetration process. Figs. 5.2(a), 66 5.2 (b) and 5.2 (c) compare the load-deflection curves from those from the experiments with the simulations. —v=4.0 mls —--—-v=3.6 mls -- v=3.0 mls ._ v=2.7 mls — v=2.4 mls —v=1.9 mls -—--v=1.7 mls —v=1.2 mls force, KN o ' . . . 0 4 8 12 16 deflection, mm Fig. 5.1- Force-deflection curves for [0/90/0] laminate for different impact velocities. — Exp. v=3.6 mls —--- Sim. v=3.6 mls force, KN deflection, mm (a) 67 -— Exp. v=2.4 mls —-- Sim. v=2.4 mls 16 deflection, mm (b) 7 _ 6 a z a; — Exp. v=1.7 m/s 2 -— Sim. v=1.7 mls .2 0 4 8 12 16 deflection, mm (c) Fig. 5.2- Comparison of force-deflection curves at three different impact velocities. 68 The delamination shape obtained from the simulations for the lower interface of the [0/90/0] at three different impact velocities, 3.6 m/s, 2.4 m/s and 1.7m/s are shown in the Fig. 5.3 (a), 5.3 (b) and 5.3 (c), respectively. 'JIIJI lljllllJ 1 11' n. 1' (a) v=3.6 mls (b) v=2.4 mls hi IIIIIIIH' IIIMIIIIIR‘ IIIIDFWWW (c) v=l.7 m/s Fig. 5.3— Delamination shape at the bottom interface of [0/90/0] laminate at three different impact velocities. 69 The delamination area vs. Impact energy from both experiments and the simulations are shown in Fig. 5.4. Clearly the simulations under predicts the delamination area. From Fig. 5.3 and 5.4 it can be concluded that the delamination model is only able to capture delamination damage qualitatively but is not able to match the experimental results quantitatively. 90 1 0 Experiments . 80 a I Simulations "s 70 - °_ 0 g 60 4 e e N c 50 d -.9. 4o .5 o 30 - 3 20 a I 10 q ° u I o T T l l T l 0 20 40 60 80 100 120 Impact energy, J Fig. 5.4- Delamination area at different impact energies. Fig. 5.5 shows the absorbed energy vs. the impact energy from both experimental and simulation results. The impact energy is based on the kinetic energy of the impactor while the absorbed energy is obtained from the integration of the force-deflection curves. 70 120 - a 100 — >7 +Experiment 2’ so , . . g +Slmulatlon I: 60 . ——-Equalenergylme o E o 40 ~ m .n < 20 - 0 4 A i i i r I O 20 4O 60 80 100 120 Impact Energy, J Fig. 5.5- Absorbed energy vs. impact energy To understand the effect of delamination damage on the response of composite laminates, two different simulations were performed. Both simulations were performed for an impact velocity of 2.4 m/s. The first simulation allowed only in-plane damage without delamination damage. The second simulation included both in-plane damage and delamination damage. The force-deflection curves from the two simulations are shown in Fig. 5.6. The energy absorbed by the composite laminate for the two cases were 6.41 J and 9.91 J. The force-deflection curves and the energy absorbed clearly indicate the role of delamination damage. Delamination caused an increase of about 49 % in the energy absorption. 71 .5 force, KN — Experiment, v=2.4 m/s 0) 1 ~—- In-plane damage but no delamination damage. v=2.4 m/s -- ln-plane damage and delamination damage, v=2.4 mls l I 0 5 10 15 deflection, mm Fig. 5.6- Force—deflection curve for [0/90/0] laminate with and without delamination damage. When finite element simulations were performed for composite laminates with a stacking sequence of [0/30/0] and [0/45/0], there were severe convergence problems. Simulations with only in-plane damage and with only out-of-plane damage converged but simulations including both damage models did not converge. Figure 5.7 compares the force-deflection curves of [0/90/0] with and without delamination, [0/45/0] with out delamination and [0/30/0] without delamination. All the simulations were performed at an impact velocity of 2.4 m/s. It can be noted that for a hypothetical case of only in-plane damage and no delamination damage, the [0/90/0] laminate was the most rigid one followed by [0/45/0] and then by [0/30/0] which had maximum in-plane damage. The energy absorption was 6.41 J, 8.23 J and 9.59 J respectively. This was due to the fact that in a [0/90/0], the 90-lamina arrested the cracks in the 0 lamina due to in-plane damage to propagate. 72 However, it may be noted that a [0/90/0] laminate has the maximum mismatch of material properties through the laminate thickness and thus would have the maximum delamination. From the earlier discussions, it was known that delamination degraded the laminate drastically. Thus, it was expected that [0/90/0] composite laminate would perform the worst for a real case with where both material in-plane damage and delamination damage took place. The results shown in Fig. 5.7 seem to point into this direction. Unfortunately, the simulations for [0/45/0] and [0/30/0] composite laminates, did not converge thus making the comparison with [0/90/0] impossible. - - - [0/90/0] with delamination —— [0/90/0] with no delamination ----- [0/45/0] with no delamination —— [0/30/0] with no delamination o l' I, T i l 0 5 10 15 deflection, mm Fig. 5.7- Force-deflection curves for different composite laminate lay-ups. 73 6. Conclusions and Recommendations 6.1 Conclusions A finite element model was presented to study the damage process of composite laminates. The commercial finite element code ABAQUS-Standard was used for the computational studies. The finite element model used two-dimensional shell elements, thus making the model computationally efficient. A homogenized continuum model was used to characterize the composite laminates. Material constitutive relations were formulated into a user subroutine UMAT in ABAQUS for in—plane analysis. The failure criteria proposed by Hashin were used. A brittle damage model was used for the post- failure behavior of the matrix while a progressive damage model was used for the fiber damage. A [O] composite laminate was investigated. The composite laminate was impacted at different impact velocities and the results from the finite element simulations were compared with the experimental results. The force deflection curves from the simulations compared well with the ones from the experiments. An elasticity analysis was performed to study the response of the interface of [90/0] composite beam under different bonding conditions. There was no elasticity analysis available in the literature for investigating delamination in composite laminates resulting in the uniqueness of this study. The well known Pagano’s cylindrical bending problem was extended for an imperfectly bonded composite beam. The [90/0] composite beam was simply supported at the ends and was subjected to a sinusoidal tensile loading at the top surface. A linear shear slip theory and a linear normal separation theory were used for the interfacial displacement to represent different bonding conditions at the interface. 74 For the finite element simulation of the delamination damage, two models were presented, the multi-point constraint model (MPC) and the interface layer model. Both models used shell elements to represent composite laminate. Individual laminae were modeled by shell elements. A user subroutine MPC was written to impose the displacement continuity conditions at the interface. In the interface model, an additional matrix layer of very small thickness was introduced between the laminae. The nodal forces were used to calculate the interlaminar stresses and a stress based failure criterion was used for the delamination damage. Once the delamination was detected, the displacement continuity conditions were removed in the MPC model while for the interface layer model the stiffness of the interface layer was reduced to make the displacement discontinuous. The two delamination models were used to analyze the [0/90] composite beam and the results were compared with the elasticity analysis. The delamination shape and size were validated with different composite layups. A typical penny shaped delamination with the major axis along the direction of the fiber direction of the bottom lamina was obtained but the size predicted from the present model had an error of about 50 % from the experimental results. The two damage models, material in-plane damage model and the out-of-plane damage model for delamination were used to investigate the response of [0/90/0] composite laminate subjected to impact loading at different impact velocities. The results from the simulations were compared with those from the experiments. The simulation results were in good agreement with the experiments. 75 6.2 Recommendations The material constitutive relations developed were not able to capture the penetration process. A more physically based post-failure damage response would be able to represent the laminate response in a better way. All strain softening models had strain localization, making them loose objectivity. The results of all strain softening models depended upon the finite element mesh being used. A non-local continuum damage model would be able to overcome such problems and thus should be investigated. The delamination model used in the present study was based on strength of material based failure criterion. Instead of it, a fracture mechanics based criterion could also be used with the same model as the nodal forces and the displacements across the interface are already available. The delamination size predicted from the current model has an error of about 50 % from the experimental results. The stress field at the delamination front is singular and thus a fracture mechanics based approach might help in reducing this error. A lot of convergence problems were encountered in the thesis research. This might be because of unstable material damage model and point constraints at the nodes across the delamination interface and a sudden removal of these constraint conditions. Moreover, an implicit solver was used for impact simulations. An explicit solver may help achieving convergence. Development of new interface elements for the delamination damage model might also help improving the convergence. 76 Appendix A Matlab code to for the analysis of imperfectly bonded composite laminate. % 90-0 plain strain composite beam under sinusoidal loading % Material properties- l-longitudinal; t-transverse el=25e6; et=le6; glt=0.5e6 ; vlt=.25; vtt=0.25; % 90 degree elu=et; e2u=et; e3u=el; v12u=vtt; v32u=vlt; g12u=gtt; vl3u= v31u*(e1u/e3u); v23u=v13u; r11u=(l/elu)—(e3u*(v13u/elu)"2) ; r12u=-(v12u/elu)-(e3u*vl3u*v23u/(elu*e2u)); r22u=( l/e2u)-(e3u*(v23u/e2u)"2); r21u=rlZu; r66u= 1/g12u; % 0 degree eld=el; e2d=et; e3d=et; v12d=.25; v13d=.25; g12d=glt; r11d=(l/e1d)-(e3d*(vl3d/eld)"2) ; r12d=-(v12d/eld)-(e3d*vl3d*v23d/(e1d*e2d)); r22d=( l/e2d)-(e3d*(v23d/e2d)"2); r21d=rl 2d; r66d=1/g12d; % length and height. h of the beam length=100; p=pi/length; h=2.5; z=2*h; % Constants for imperfect bonding. j-shear slip,'k- normal separation % higher the value of j and k poor is the bonding j=XX; =XX; % Constants mij for Pagano ’3 solution au=r66u+(2*r1 2u); 77 gtt=0.2e6; v3 1u=vlt; v23d=.25; ad=r66d+(2*r12d); bu=((au"2)-(4*rl lu*r22u))"0.5; bd=((ad"2)-(4*r1 1d*r22d))"0.5; ; cu=2*r11u; cd=2*rl 1d; % m for top layer mlu= p*(((au+bu)/cu)"0.5); m2u=-m1u; m3u= p*(((au-bu)/cu)"0.5) ; m4u=-m3u; % m for bottom layer m1d=p*(((ad+bd)/cd)"0.5); m2d=-mld; m3d=p*(((ad-bd)/cd)"0.5); m4d=-m3d; %normal stress on top surface row1= [ -exp(mlu*h/2) -exp(m2u*h/2) -exp(m3u*h/2) -exp(m4u*h/2) O O O 0]; %normal stress on bottom surface row2=[ O 0 O O exp(-mld*h/2) exp(—m2d*h/2) exp(-m3d*h/2) exp(-m4d*h/2)]; % transverse shear stress on the top surface row3=[ mlu*exp(m1u*h/2) m2u*exp(m2u*h/2) m3u*exp(m3u*h/2) m4u*exp(m4u*h/2) O O O O]; % transverse shear stress on the bottom surface row4= [O O 0 O m1d*exp(-mld*h/2)m2d*exp(-m2d*h/2) m3d*exp(-m3d*h/2) m4d*exp(-m4d*h/2)]; % transverse normal stress continuity at the interface row5=[exp(-mlu*h/2) exp(-m2u*h/2) exp(-m3u*h/2) exp(-m4u*h/2) -exp(mld*h/2) - exp(m2d*h/2) -exp(m3d*h/2) -exp(m4d*h/2)]; % transverse shear stress continuity at the interface row6=[mlu*exp(-mlu*h/2) m2u*exp(-m2u*h/2) m3u*exp(-m3u*h/2) m4u*exp(- m4u*h/2) -m1d*exp(mld*h/2) -m2d*exp(m2d*h/2) -m3d*exp(m3d*h/2) - m4d*exp(m4d*h/2)]; % linear shear slip theory for the shear slip at the interface 78 row7=[ (((r12u*p"2)-(rl lu*mlu"2))*exp(-mlu*h/2)) + (j*(pA2)*mlu*exp(-m1u*h/2)), (((r12u*p"2)-(rl 1u*m2u"2))*exp(-m2u*h/2))+(j*(pA2)*m2u*exp(-m2u*h/2)), (((r12u*p"2)-(rl1u*m3u"2))*exp(-m3u*h/2))+ (j*(pA2)*m3u*exp(-m3u*h/2)), (((r12u*p"2)-(rl lu*m4u"2))*exp(-m4u*h/2))+(j*(pA2)*m4u*exp(-m4u*h/2)) -((r12d*p"2)-(r1 ld*mld"2))*exp(mld*h/2) -((r12d*p"2)-(rl ld*m2d"2))*exp(m2d*h/2) -((r12d*p"2)-(rl ld*m3d"2))*exp(m3d*h/2) -((r12d*p"2)-(rl ld*m4d"2))*exp(m4d*h/2)]; % linear normal separation theory for the normal separation at the interface row8=[(((rl2u*mlu)-((r22u/mlu)*p"2))*exp(-mlu*h/2))+(k*(p"2)*exp(-mlu*h/2)), (((rl2u*m2u)-((r22u/m2u)*p"2))*exp(-m2u*h/2))+(k*(p"2)*exp(-m2u*h/2)),... (((rl2u*m3u)—((r22u/m3u)*p"2))*exp(-m3u*h/2))+(k*(pA2)*exp(-m3u*h/2)), (((rl2u*m4u)-((r22u/m4u)*p"2))*exp(-m4u*h/2))+(k*(p"2)*exp(-m4u*h/2)),... -((r12d*m1d)-((r22d/mld)*p"2))*exp(mld*h/2) -((rl2d*m2d)-((r22d/m2d)*p"2))*exp(m2d*h/2) -((rl2d*m3d)-((r22d/m3d)*p"2))*exp(m3d*b/2) -((rl2d*m4d)-((r22d/m4d)*p"2))*exp(m4d*h/2)]; % the LHS, stiffness M= [rowl; row2; row3; row4; rowS; row6; row7; row8]; % the RHS, loading rhS=[1/(p"2); 0; 0; 0; 0; 0; 0; 0;]; % solution A=inv(M)* rhs; 79 a1u=A(l); a2u=A(2); a3u=A(3); a4u=A(4); a1d=A(5); a2d=A(6); a3d=A(7); a4d=A(8); % normalized deflection at the mid-span def2= (((rl2u*m1u)-((r22u/mlu)*p"2))*A(l)*exp(-m1u*h/2))+ (((rl2u*m2u)-((r22u/m2u)*p"2))*A(2)*exp(-m2u*h/2))+ (((rl2u*m3u)—((r22u/m3u)*p"2))*A(3)*exp(-m3u*h/2))+ (((rl2u*m4u)-((r22u/m4u)*p"2))*A(4)*exp(-m4u*h/2)) normaldef= (l00*et*(z"3)*def2)/(length"4) % transverse shear stress through the thickness tstress(l)= -p*((alu*mlu*exp(m1u*h/2)) + (a2u*m2u*exp(m2u*h/2))+ (a3u*m3u*exp(m3u*h/2))+(a4u*m4u*exp(m4u*h/2))); tstress(2)= -p*( (alu*m1u*exp(-m1u*h/2))+(a2u*m2u*exp(-m2u*h/2))+ (a3u*m3u*exp(-m3u*h/2))+(a4u*m4u*exp(-m4u*h/2))); tstress(3)= -p*( (a1d*mld*exp(mld*h/2))+(a2d*m2d*exp(m2d*h/2))+ (a3d*m3d*exp(m3d*h/2))+(a4d*m4d*exp(m4d*h/2))); tstress(4)= -p*( (a1d*mld*exp(-mld*h/2))+(a2d*m2d*exp(-m2d*h/2))+ (a3d*m3d*exp(-m3d*h/2))+(a4d*m4d*exp(-m4d*h/2))); % transverse normal stress though the thickness nstress( 1 )= -(p"2)*( (a1u*exp(m1u*h/2))+(a2u*exp(m2u*h/2))+ (a3u*exp(m3u*h/2))+(a4u*exp(m4u*h/2))); nstress(2)= -(p"2)*( (a1u*exp(-mlu*h/2))+(a2u*exp(-m2u*h/2))+ (a3u*exp(-m3u*h/2))+(a4u*exp(-m4u*h/2))); nstress(3)= -(p"2)*( (a1d*exp(mld*h/2))+(a2d*exp(m2d*h/2))+ (a3d*exp(m3d*h/2))+(a4d*exp(m4d*h/2))); nstress(4)= -(p"2)"‘( (ald*exp(-mld*h/2))+(a2d*exp(-m2d*h/2))+ (a3d*exp(-m3d*h/2))+(a4d*exp(-m4d*h/2))); % inplane stress through the thickness inplanestress(l)= (a1u*(mlu"2)*exp(mlu*h/2)) +(a2u*(m2u"2)*exp(m2u*h/2)) + (a3u*(m3u"2)*exp(m3u*h/2)) +(a4u*(m4u"2)*exp(m4u*h/2)); inplanestress(2)= (alu*(m1u"2)*exp(-m1u*h/2))+(a2u*(m2u"2)*exp(-m2u*h/2))+ 80 (a3u*(m3u"2)*exp(-m3u*h/2))+(a4u*(m4u"2)*exp(-m4u*h/2)); inplanestress(3)= (a1d*(m1d"2)*exp(mld*h/2)) +(a2d*(m2d"2)*exp(m2d*h/2)) + (a3d* (m3d"2)*exp(m3d*h/2)) + (a4d*(m4d"2)*exp(m4d*h/2)); inplanestress(4)= (ald*(m1d"2)*exp(-mld*h/2)) +(a2d*(m2d"2)*exp(-m2d*h/2)) + (a3d*(m3d"2)*exp(-m3d*h/2)) +(a4d*(m4d"2)*exp(-m4d*h/2)); % out of plane de ection throu h the thickness g outplanedef( l )= (((r12u*mlu)—((r22u/ml u)*p"2))*a1u*exp(m lu*h/2)+ ((rl2u*m2u)-((r22u/m2u)*p"2))*a2u*exp(m2u*h/2)+ ((rl 2u*m3u)-((r22u/m3u)*p"2))*a3u*exp(m3u*h/2)+... ((rl2u*m4u)-((r22u/m4u)*p"2))*a4u*exp(m4u*h/2)); outplanedef(2)=(((r12u*m lu)-((r22u/m lu)*p"2))*alu*exp(-m1u*h/2)+ ((rl 2u*m2u)-((r22u/m2u)*p"2))*a2u*exp(-m2u*h/2)+ ((rl 2u*m3u)-((r22u/m3u)*p"2))*a3u*exp(-m3u*h/2)+... ((rl 2u*m4u)-((r22u/m4u)*p"2))*a4u*exp(-m4u*h/2)); outplanedef(3)=(((rl2d*m1d)-((r22d/mld)*p"2))*ald*exp(mld*h/2)+ ((rl2d*m2d)-((r22d/m2d)*p"2))*a2d*exp(m2d*h/2)+ ((rl2d*m3d)-((r22d/m3d)*p"2))*a3d*exp(m3d*h/2)+... ((rl 2d*m4d)-((r22d/m4d)*p"2))*a4d*exp(m4d*h/2)); outplanedef(4)=(((rl 2d*m1d)-((r22d/m1d)*p"2))*a1d*exp(—mld*h/2)+ ((rl 2d*m2d)-((r22d/m2d)*p"2))*a2d*exp(-m2d*h/2)+ ((rl2d*m3d)-((r22d/m3d)*p"2))*a3d*exp(-m3d*h/2)+... ((rl2d*m4d)-((r22d/m4d)*p"2))*a4d*exp(-m4d*h/2)); outplanedef=(l00*et*(2A3)‘outplanedeO/(lengthM); % in plane deflection through the thickness inplanedef( l)=( l/p)*(((r12u*p"2)-(rl 1u*m lu"2))*al u*exp(m 1u*h/2)+ ((r12u*p"2)-(r1 1u*m2u"2))*a2u*exp(m2u*h/2)+... ((r12u*p"2)-(rl lu*m3u"2))*a3u*exp(m3u*h/2)+ ((r12u*p"2)-(rl lu*m4u"2))*a4u*exp(m4u*h/2)); inplanedef(2)=(l/p)*(((r12u*p"2)-(rl 1u*m lu"2))*a l u*exp(-m l u*h/2)+ ((r12u*p"2)-(rl 1u*m2u"2))*a2u*exp(-m2u*h/2)+... ((r12u*p"2)-(rl 1u*m3u"2))*a3u*exp(-m3u*h/2)+ ((r12u*p"2)-(r1 1u*m4u"2))*a4u*exp(-m4u*h/2)); inplanedef(3)=( 1/p)*(((r1 2d*p"2)-(rl 1d*m 1 d"2))*a 1 d*exp(m l d*h/2)+ ((r12d*p"2)-(r1 ld*m2d"2))*a2d*exp(m2d*h/2)+... 81 ((r12d*p"2)-(rl 1d"m3d"2))*a3d*exp(m3d*h/2)+ ((r12d*p"2)-(r1 1d*m4d"2))*a4d*exp(m4d*h/2)); inplanedef(4)=( l /p)*(((r1 2d*p"2)-(rl 1d*mld"2))*ald*exp(-m1d*h/2)+ ((r12d*p"2)-(r1 ld*m2d"2))*a2d*exp(-m2d*h/2)+... ((r12d*p"2)-(r1 ld*m3d"2))*a3d*exp(-m3d*h/2)+ ((r12d*p"2)-(rl ld*m4d"2))*a4d*exp(-m4d*h/2)); 82 Appendix B Subroutine MPC to enforce the displacement continuity conditions SUBROUTINE MPC(UE,A,JDOF,MDOF,N,JTYPE,X,U,UINIT,MAXDOF, * LMPC,KSTEP,KINC,TIME,NT,NF,TEMP,FIELD,LTRAN,TRAN) C INCLUDE 'ABA_PARAM.INC' C DIMENSION A(N),JDOF(N),X(6,N),U(MAXDOF,N),UINIT(MAXDOF,N), * TIME(2),TEMPO\IT,N),FIELD(NF,NT,N),LTRAN(N),TRAN(3,3,N) COMMON DELPARAM(168l),REACTIONX(168l),REACTIONY(l681), l RZ(1681) C OPEN(UNIT= l 7,FILE='debug.dat') if(jtype.eq.100) then a(l)=l a(2)=-l a(3)=-l jdof(1)=3 jdof(2)=3 jdof(3)=3 ue=u(3,2)+u(3,3) C DELAMINATION CHECK C node number passed in as field variable jnodeno=int(field(l , l ,2)) if(kinc.eq. 1 )then delparamGnodeno)=O endif if(delparamO nodeno). ge. 1 )then lmpc=0 endif endif if(jtype.eq.200) then a(l)=l a(2)=0.635 a(3)=-l a(4)=1.27 a(5)=-1 jdof(1)=l jdof(2)=5 jdof(3)=l jdof(4)=5 jdof(5)=l ue=(-O.635*u(5,2))+u(1,3)+(-l.27*u(5,4))+u(1,5) 83 C DELAMINATION CHECK jnodeno=int(field( 1 , l ,3)) if(kinc.eq. 1)then delparamOnodeno)=O endif if(delparam(j nodeno). ge. 1)then lmpc=0 endif endif if(jtype.eq.300) then a(l)=l a(2)=-O.635 a(3)=-l a(4)=-1.27 a(5)=-l jdof(1)=2 jdof(2)=4 jdof(3)=2 jdof(4)=4 jdof(5)=2 ue=(0.635 *u(4,2))+u(2,3)+( l .27*u(4,4))+u(2,5) jnodeno=int(field( l , l ,3)) if(kinc.eq. 1)then delparamfinodeno)=0 endif if(delparam(j nodeno). ge. 1)then lmpc=0 endif endif if(jtype.eq.400) then a(1)=l a(2)=-l jdof(1)=4 jdof(2)=4 ue=u(4,2) jnodeno=int(field( l , l ,2)) if(kinc.eq. 1)then delparamfinodeno)=0 endif if(delparam(jnodeno).ge. 1)then lmpc=0 endif endif if(jtype.eq.500) then a(1)=1 a(2)=-1 84 jdof( l)=5 jdof(2)=5 ue=u(5,2) jnodeno=int(field( l , l ,2)) if(kinc.eq.1)then delparamG nodeno)=0 endif if(delparam(j nodeno). ge. 1)then lmpc=0 endif endif RETURN END 85 Appendix C Subroutine URDFIL to calculate the interlaminar stresses from the nodal forces SUBROUTINE URDFIL(LSTOP,LOVRWRT,KSTEP,KINC,DTIME,TIME) INCLUDE 'ABA_PARAM.INC' DIMENSION ARRAY(513),JRRAY(NPRECD,513),TIME(2), l JNODENUM(1681),VMS(1681) COMMON DELPARAM(1681),REACTIONX(1681),REACTIONY(1681), 1 RZ(1681) REAL NF ORCE,NSTRENGTH,NFAILF ORCE,INPLANE,AINT,BINT,CINT, l DINT,A,D EQUIVALENCE (ARRAY( l ),JRRAY(] , 1)) C OPEN(UNIT=17,FILE='debug.dat') C FIND IN -PLANE STRESS AT THE INTERFACE RECORD -12 call posfil(kstep,kinc,array,jrcd) do k1=l,999999 call dbfile(0,array,jrcd) if (jer .ne. 0) go to 110 key=jrray( l ,2) if(key.eq. I )then jn=irray(1,3) if(jn.le. 168 1)then j=jn j8p=jrray(l,5) endif endif if(key.eq. 12)then if(j sp.eq. 1)then vms(j)=array(3) endif endif C FIND THE NODE NUMBER AND REACTION FORCE FROM RECORD 104 if(key.eq. 1 04) then jn0=jrfay(1,3) if(jno.gt.4000)then if(jno.lt.5682)then i=jno-4000 jnodenum(i)=jno reactionx(i)=array(4) reactiony(i)=array(5) rz(i)=array(6) 86 elarea=9.76 fx=reactionx(i)/elarea fy=reactiony(i)/elarea fz=rz(i)/elarea sforce=((fx**2)+(fy**2)) sstrength=70.0d3 nstrength=39.0d3 nforce=(fz**2) sfailforce=(sstrength**2) nfailforce=(nstrength**2) if(kinc.eq. 1)then delparam(i)=0 endif if(kinc. gt. 1)then if(delparam(i).lt. 1 .O)then if(fz.le.O.dO)then delparam(i)=0.0 endif if(fz.gt.0.0)then inplane=sqrt((vms(i)/(70.65d3))**2) delparam(i)=(sqrt(sforce/sfailforce))+outplane+inplane endif endif endif endif endif END IF END DO 110 CONTINUE RETURN END 87 Appendix D Subroutine USDFLD for changing the stiffness of the interface layer for the delamination model with interface layer. SUBROUTINE USDFLD(FIELD,STATEV,PNEWDT,DIRECT,T,CELENT,TIME, l DTIME ,CMNAME,ORNAME,NFIELD,NSTATV,NOEL,NPT,LAYER,KSPT, 2 KSTEP,KINC,NDI,nshr,coord,jmachtyp,matlayo,laccflg) C INCLUDE 'ABA_PARAM.INC' C PARAMETER(YT=38.0D3,YC=128.0D3,SL=70.0D3) CHARACTER*8O CMNAME,ORNAME CHARACTER*8 F LGRAY(15) DIMENSION FIELD(NFIELD),STATEV(NSTATV),DIRECT(3,3),T(3,3), * TIME(2),coord(*),jmac(*),jmtyp(*) DIMENSION ARRAY(15),JARRAY(15) COMMON sh]3(1681),sh23(l681),SN(1681) INTEGER ELNO,N l ,N2,N3,N4,RELNODE,NSP REAL DELAM,SIG l 3,SIG23,SIG33,VS,OUTPLAN E c OPEN(UNIT= 1 7,FILE='debug’) c ELEMENT NUMBER elno=noel-2000 nsp=5 C DELAMINATION STATE delam=statev( 1 ) if(kinc.ge.2)then c RECOVER NODE NUMBER FROM ELEMENT NUMBER if(elno.lt. 1 600)then relnode=int(elno/40) nl=elno+relnode n2=nl+1 n3=n1+41 n4=n1+42 C CHECK FOR THE CORRESPONDING INTERGATTION POINT if(npt.eq. 1)then sigl3=shl3(nl) sig23=sh23(nl) sig33=SN(nl) endif if(npt.eq.2)then sig13=shl3(n2) 88 sig23=sh23(n2) sig33=SN(n2) endif if(npt.eq.3)then sig13=sh13(n4) sig23=sh23(n4) sig33=SN(n4) endif if(npt.eq.4)then sigl3=shl3(n3) sig23=sh23(n3) sig33=SN(n3) endif c GET THE STRESES CALL GETVRM('SINV',ARRAY,JARRAY,FLGRAY,jrcd, $ jmac, jmtyp, matlayo, laccflg) vs = array(l) if(delam.lt. l .d0)then if(sig33.le.0.d0)then delam=0 else outplane=((sig13/sl)**2)+((sig23/sl)**2)+((sig33/yt)**2) delam= ((vs**2)/(yt*yc))+outplane endif endif if(delam.ge. 1 .d0)then delam=999.0 endif endif endif C UPDATING THE FIELD VARIABLES C if(delam.lt. 1 .0)then field(1)=0 statev(l)=0 endif if(delam. ge. 1.0)then field(l)=l.0 statev( 1):] .0 endif RETURN END 89 Appendix E UMAT subroutine for the material in-plane damage C OOOOOOOOOOOOOOOOO SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, l RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN, 2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI, 3 NSHR,NTENS, 4 NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT, 5 DFGRDO,DFGRD1 ,NOEL,NPT,KS LAY,KSPT,KSTEP,KINC) INCLUDE 'ABA_PARAM.INC' CHARACTER*80 MATERL DIMENSION STRESS(NTENS),STATEV(NSTATV), l DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS), 2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(l), 3 PROPS(NPROPS),COORDS(3),DROT(3,3), 4 DFGRDO(3,3),DFGRD1(3,3), DPRED(1), 5 DELGl1(3),DELGZZ(3),DDSDD(NTENS,NTENS) REAL EMODI 1,EMOD22,GMOD12,XNU,SIG1,SIG2,SIG12, 1 ST] ,ST2,ST3, loadl l,load22, MATERIAL STRENGTH PARAMETER PARAMETER (XT=1080.D3,XC=620.0D3,YT=39.0D3, l YC=128.0D3,SC=89.0D3) OPEN(UNIT= 1 7,FILE='debug.dat') UMAT FOR COMPOSITE MATERIAL SHOULD BE USED WITH PLAIN STRESS CASES- ONLY SIGMAl l,SIGMA22 and SIGMAlZ PROPS(1) - E11 PROPS(2) - E22 PROPS(3) - NU12 PROPS(4) - G12 BRITTLE MATRIX FAILURE AND PROGRESSIVE FIBER FAILURE EF-FIBER FAILURE EM-MATRIX FAILURE EFS- FIBER-MATRD( SHEAR FAILURE ELASTIC PROPERTIES 90 emodl l=props( l) emod22=props(2) xnu12=props(3) gmod12=props(4) xnu2 l=(xnu12*emod22)/emodl l ef=statev( 1) em=statev(2) efs=statev(3) sigld=statev(4) emodl lold=statev(5) stlold=statev(6) if(em.ge. l .dO)then emod22=8.6e3 ‘ if(stress(2).gt.0.d0)then “ gmod12=3.8e3 endif xnu2l=0.062*emod22/(86000O0.0) endif if(efs. ge. l .dO)then gmod12=3.8e3 endif if(efge. l .dO)then emod22=8.6e3 gmod12=3.8e3 emodl l=emodl lold xnu12=0.28*emodl l/(39000OO0.0) C JACOBIAN MATRIX if fiber fails ddsdde(1,1)=-emodl 1*xnu21 *xnu12/ l ((l-(xnu12*xnu21))**2) ddsdde( l ,2)=emodl l *xnu12*(xnu21 * *2)/ l ((l-(xnu12*xnu21))**2) ddsdde(] ,3)=O ddsdde(2, l )=ddsdde( 1,2) ddsdde(2,2)=emod22/( 1 -(xnu l 2 *xnu2 l )) ddsdde(2,3)=0 ddsdde(3 , 1 )=0 ddsdde(3,2)=0 ddsdde(3 ,3)=gmod1 2 endif C JACOBIAN MATRIX IF NO FIBER FAILS if(ef.lt. l .dO)then ddsdde( l , l )=emodl l/(l-(xnu12*xnu2 1)) ddsdde( l ,2)=xnu12 *emod22/( l-(xnu12 *xnu2 1 )) 91 ddsdde( l ,3 )=0 ddsdde(2, l )=ddsdde( l ,2) ddsdde(2,2)=emod22/(1-(xnu1 2 *xnu2 l )) ddsdde(2,3)=0 ddsdde(3 , 1 )=0 ddsdde(3,2)=0 ddsdde(3 ,3 )=gmod1 2 endif C ELASTIC STRESS E stl=stran(l)+ dstran(1) st2=stran(2)+ dstran(2) st3=stran(3)+ dstran(3) sigl=(ddsdde(l,1)*st1)+(ddsdde(l,2)*st2) sig2=(ddsdde(2, 1 )*st 1 )+(ddsdde(2,2)*st2) sig12=ddsdde(3,3)*st3 if(sigl .lt.O.dO)then x=xc else x=xt endif if(sig2.lt.0.d0)then FYC else y=yt endif c FIBER FAILURE if(ef.lt. l .dO)then ef=sqrt((sigl **2)/(x**2)) endif if(ef. ge. 1 .dO)then if(ef.ne.999.d0)then sigld=abs(sig1*0.2) endif endif if(ef.ge. 1 .dO)then ef=999.d0 endif C MATRIX FAILURE if(em.lt. l .dO)then em=sqrt(((si32**2)/(y**2))+((Sig12**2)/(SC**2))) else em=l .dO endif C FIBER MATRIX SHEAR FAILURE if(efs.lt. l .dO)then 92 efs=sqrt(((si g1 * *2)/(x**2))+((si g l 2**2)/(sc* *2))) else efs=l.d0 endif C STIFFNESS DEGRADATION- C MATRD( FAILURE if(em. ge. 1 .dO)then emod22=8.6e3 if(sig2.gt.0.d0)then gmod12=3.8e3 endif xnu2l=0.062*emod22/(8600000.0) endif C FIBER MATRD( SHEAR FAILURE if(efs. ge. l .dO)then gmod12=3.8e3 endif C FIBER FAILURE if(ef. ge. 1 .dO)then emod22=8.6e3 gmod12=3.8e3 C CHECK FOR LAODING/UNLOADING if(abs(st1).gt.stlold)then c loading emodl l=abs(sigld/st1) else c unloading emod11=emod1 lold endif xnu12=O.28*emodl 1/(39000OO0.0) C JACOBIAN MATRIX ddsdde( 1 , l )=-emod1 1 *xnu21 *xnu12/ 1 ((1-(xnu12*xnu21))**2) ddsdde(1,2)=emodl l *xnu12*(xnu2 1 * *2)/ l ((1-(xnu12*xnu21))**2) ddsdde( l ,3 )=0 ddsdde(2, l )=ddsdde( l ,2) ddsdde(2,2)=emod22/(l -(xnu12*xnu21)) ddsdde(2,3)=0 ddsdde(3 , 1 )=0 ddsdde(3,2)=0 ddsdde(3 ,3)=gmod 1 2 endif IF (EF .LT. 1 .DO)THEN ddsdde(] ,1)=emodl l/(l-(xnu12*xnu21)) 93 ddsdde( l ,2)=xnu12 *emod22/( l -(xnu 12 *xnu2 1)) ddsdde( l ,3)=O ddsdde(2, l )=ddsdde( 1 ,2) ddsdde(2,2)=emod22/( 1 -(xnu12 *xnu2 1 )) ddsdde(2,3)=0 ddsdde(3 , 1 )=0 ddsdde(3 ,2)=O ddsdde(3 ,3)=gmod12 endif cl l=emod1 1/(1-(xnu12*xnu21)) c12=xnu12*emod22/(l-(xnu12*xnu21)) c22=emod22/(l-(xnu12*xnu21)) c66=gmod12 stress(l)=(c1 1*stl)+(c12*st2) stress(2)=(C 12* ST 1 )+(C22* ST2) stress(3 )=c66*st3 statev( 1 )=ef statev(2)=em statev(3)=efs statev(4)=si g 1 d statev(5)=emod1 l statev(6)=stl RETURN END 94 BIBLIOGRAPHY 1Reddy, 1. N., “A Simple Higher Order Theory for Laminated Composite Plates,” Journal of composite materials, Vol.51, 1984, pp.745-752. 2Liu, D. and Li, X., “An Simple Overall View of Laminate Theories Based on Displacement Hypothesis,” Journal of composite materials, Vol.30, 1996, pp.1539-1561. 3‘Hashin, Z. and Rotem, A., “A Fatigue Failure Criterion for Fiber Reinforced Materials,” Journal of composite materials, Vol.7, 1973, pp.448-464. 4Hashin, Z., “Failure Criteria for Unidirectional Fiber Composites,” Journal of Applied Mechanics, Vol.47, 1980, pp.329-334. 5Yamanda, S. E. and Sun, C. T., “Analysis of Laminate Strength and Its Distribution,” Journal of composite materials, Vol.12, 1978, pp.275-284. 6Chang, F.-K. and Chang, K. Y., “A Progressive Damage Model for Laminated Composites Containing Stress Concentrations,” Journal of Composite Materials, Vol.21, Sept.l987, pp. 834-855. 7Chang, F.-K. and Lessard, L. B., “Damage Tolerance of Laminated Composites Containing An Open Hole And Subjected To Compressive Loadings: Part I Analysis,” Journal of Composite Materials, Vol.25, Sept.l99l, pp. 1-42. 8Shahid, I. and Chang, F.-K., “An Accumulative Damage Model for Tensile and Shear Failures of Laminated Composite Plates,” Journal of Composite Materials, Vol.29, 1995, pp. 926-981. 9Hart-Smith, “A New Approach to Fibrous Composite Laminate Strength Prediction,” Eighth DOD/NASA/FAA Conference on Fibrous Composites in Structural Design, NASA CP-3087, Part 2, 1989, pp. 663-693. 10Kroll, L. and Hufenbach, W., “Physically Based Failure Criterion for Dimensioning of Thick-Walled Laminates,” Applied composite materials, Vol.4, 1997, pp. 321-332. llAl-Bastaki, N. M. 8., “Design of Fiber Reinforced Composite Structures Subjected to High Strain Rates Using Finite Element Analysis,” Applied Composite Materials, Vol.5, 1998, pp. 223-236. l2Huang, C. H., and Lee, Y. 1., “Experiments and Simulation of the Static Contact Crush of Composite Laminated Plates,” Composite Structures, Vol.61, 2003, pp. 265- 270. 95 13Krajcinovic, D., Damage Mechanics. Elsevier Science: Amsterdam, 1996. l4Bazant, Z. P., Belytschko, T. B. and Chang, T. P., “Continuum Model for Strain Softening,” Journal of Engineering Mechanics ASCE, 110, 1984, pp. 1666-1692. lsBazant, Z. P., “Why Continuum Damage is Nonlocal: Micromechanics Arguments,” Journal of Engineering Mechanics ASCE, l 17, 1991, pp. 1070-1087. 16Bazant, Z. P. and Lin, F. B., “Nonlocal Smeared Crack Model for Concrete fracture,” Journal of Structural Engineering, 114, 1988, pp. 2493-2510. l7Johnson, A. F., “Modelling Fabric Reinforced Composites Under Impact Loads,” Composites : Part A, Vol.32, 2001, pp. 1197-2007. l8Johnson, A. F., Pickett, A. K. and Rozycki, P., “Computational Methods for Predicting Impact Damage in The Composite Structures,” Composites Science and Technology, Vol.61, 2001, pp. 2183-2192. l9Matzenmiller, A., Lubliner, J., and Taylor, R. L., “A Constitutive Model for Anisotropic Damage in Fiber-Composites,” Mechanics of Materials, Vol.20, 1995, pp. 125-152. 20Williams, K. V. and Vaziri, R., “Application of a Damage Mechanics Model for Predictions The Impact Response of Composite Materials,” Composites and Structures, Vol.79, 2001, pp. 997-1011. 2|Schapery, R. A. and Sicking, D. L., “On Nonlocal Constitutive Equations for Elastic and Viscoelastic Composites with Growing Damage,” Mechanical Behavior of Materials and Structures, Vol.79, 2001, pp. 997-1011. 22Basu, S and Waas, A. M., “Computational Modeling of Damage Growth in Composite Laminates,” AIAA, Vol.41, 2003, pp. 1158-1166. 23Kennedy, T. C. and Nahan, M. F., “A Simple Nonlocal Damage Model For Predicting Failure of Notched Laminates,” Composite Structures, Vol.35, 1996, pp. 229- 336. 24Liu, D., “Impact- Induced Delamination- A View of Bending Stiffness Mismatching,” Journal of Composite Materials, Vol.22, 1988, pp. 674-691. 25Brewer, J .C., and Lagace, P.A., “Quadratic Stress Criterion for Initiation of Delamination,” Journal of Composite Materials, Vol.22, 1988, pp. 1141-1155. 96 26 Zhou, S. G. and Sun, C. T., “Failure Analysis of Composite Laminates with Free Edge,” Journal of Composites Technology and Research, 12, 1990, pp. 91-97. 27Kim, R. Y., and Soni, S. R., “Experimental and Analytical Studies on the Onset of Delamination in Laminated Composites,” Journal of Composite Materials, Vol. 18, 1984, pp.70—80. 28Fenske, M. T. and Vizzini, A. J., “The Inclusion of In—Plane Stresses in Delamination Criteria,” Journal of Composite Materials, Vol. 35, 2001, pp. 1325-1342. 29 de Moura, M. F. S. F. and Marques, A. T., “Prediction of Low Velocity Impact Damage in Carbon-Epoxy Laminates,” Composites Part A: Applied Science and Manufacturing, Vol.33, pp. 361-368. 30Lou, R. K., Green, E., R. and Morrison, C., “Impact damage analysis of composite plates,” International Journal of Impact Engineering, Vol.22, 1999, pp. 435-447. 31Choi, Y. and Chang, F.-K., “A Model for Predicting Damage in Graphite/Epoxy Laminate Composites Resulting from Low Velocity Impact,” Journal of Composite Materials, Vol.26, Sept.l992, pp. 2134-67. 1\fjiFinn, S. R. and Springer, G. S., “Delaminations in Composite Plates Under Transverse Static or Impact Loads,” Composite Structures, Vol.23, 1993, pp. 177-190. 33Hou, J. P., Petn'nic, N. and Ruiz. C., “Prediction of Impact Damage in Composite Plates,” Composites Science and Technology, Vol.60, 2000, pp. 273-281. 34Mukherjee, Y. X., Gulrajani, S. N. and Mukherjee, S., “A Numerical and Experimental Study of Delaminated Layered Composites,” Journal of Composite Materials, Vol.28, No. 9, 1994, pp. 837-870. 35Lia, C. F., Hub, N., Yina, Y. J., Sekinec, H. and Fukunagac, H., “Low-Velocity Irnpact-Induced Damage of Continuous F iber-Reinforced Composite Laminates. Part I. An FEM Numerical Model,”Composites Part A: Applied science and manufacturing, Vol.33, 2002, pp. 1055-1062. 36 Sprenger, W., Gruttmann, F. and Wagner, W., “Delamination Growth Analysis in Laminated Structures with Continuum-Based 3d-Shell Elements and a Viscoplastic Softening Model,” Computer Methods in Applied Mechanics and Engineering, 2000, pp. 123-139. 37Zheng, S., and Sun, C. T., “A Double Plate FEM Model for the Impact Induced Delamination Problem,” Composites Science and Technology, Vol.53, 1995, pp. 111-118. 97 38Klug, J., Wu, X. X. and Sun, C. T., “Efficient Modeling of Postbuckling Delamination Growth in Composite Laminates Using Plate Elements,” AIAA, Vol.34, No. 1,1996, pp. 178-184. 39Wisnom, M. R. and Chang, F.-K., “Modeling of Splitting and Delamination in Notched Cross-Ply Laminates,” Composites Science and Technology, Vol.60, 2000, pp. 2849-2856. 40Allix. O. and Ladeveze, P., “Damage Analysis of Interlaminar Fracture Specimens,” Composite Structures, Vol.31, 1995, pp. 61-74. “Roddy, E. D. Jr., Mello, F. J., and Guess, T. R., “Modeling and Initiation of Growth of Delaminations in Composite Structures,” Journal of Composite Materials, Vol.31, No. 8, 1997, pp. 812-821. 42Gonzalves, J.P.M., de Moura, M.F.S.F., de Castro, P.M.S.T., and Marques, A.T., “Interface Element Including Point-to-Surface Constraints for Three-Dimensional Problems with Damage Propagation,” Engineering Computations,Vol.17, 2000, pp. 28- 47. 43Zou, Z., Reid, S. R., Soden, P. D. and Li, S., “Mode Separation of Energy Release Rate for Delamination in Composite Laminates using Sublaminates,” International Journal of Solids and Structures, Vol.38, 2001, pp. 2597-2613. 44Zou, 2., Reid, S. R., Soden, P. D. and Li, S., “Application of Delamination Model to Laminated Composite Structures,” Composite Structures, Vol.56, 2002, pp. 375-389. 145-Liu, M. L. and Yu, J., “Finite Element Modeling of Delamination by Layerwise Shell Element Allowing for Interlaminar Displacements,” Composites Science and Technology, Vol.63, 2003, pp. 517-529. {,46Pantano, A., Averill, R.C., “A Mesh-Independent Interface Technology for Simulation of Mixed-Mode Delamination Growth,” International Journal Of Solids And Structures , Vol. 41, 2004, pp. 3809-3831. 47N. J. Pagano, “Exact Solutions for Composite Laminates in Cylindrical Bending,” Journal of Composite Materials, Vol. 3, 1969, pp. 398-411. 48Liu, D., Xu, L. and Lu, X., “Stress Analysis of Imperfect Composite Laminates with an Interlaminar Bonding Theory,” International Methods In Numerical Methods In Engineering, Vol. 37, 1994, pp. 2819-2839. 49ABAQUS-Theory manual and Standard Users Manual- Version 6.3, Habbit. Karlson and Sorensen. 98 50Sun, C. T. and Chen, J. K., “On the impact of initially stresses laminates,” Journal of Composite Materials, Vol.19, No. 11, 1985, pp. 490-504. 99 IES 1111111111111111111 . *‘m"" ~.-.-_ __ _ , ,