UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations  Vendor Digitized Files   Help 
Material Information
Subjects
Notes
Record Information

Full Text 
FINITE ELEMENT MODELING OF STATIC INDENTATION DAMAGE IN LAMINATED COMPOSITES BY HAN SOO JUNG A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1991 DEDICATION This thesis is dedicated to my fatherinlaw, Dong Seung Kim, who could not wait to see this thesis. He passed away in August 1989. Also this thesis is dedicated to my younger brotherin law, strong and courageous seafarer, Captain Ja Sup Kim, who must have perished with his crew sometime during September 1990 on his last voyage from Chile. ACKNOWLEDGEMENTS First of all, acknowledgement must go to the professors of my supervisory committee. Without the care and guidance of Dr. Bhavani V. Sankar, chairman of the committee, this work could not have been possible. Dr. Lawrence E. Malvern, who also was a member of the supervisory committee for my master program, encouraged me to continue my graduate study to the doctoral program. Dr. C. T. Sun introduced me to composite materials, about which I had no experience before. Dr. Ibrahim Ebcioglu taught me about plate theories. His encouragement always helped me when I was dispirited. Dr. John J. Mecholsky has widened my vision about composite materials. His many suggestions greatly helped this research. I also heartily appreciate the financial assistance I received during the doctoral program. My special appreciation extends to Dr. C. A. Ross and Dr. J. E. Milton at the Graduate Center of the University of Florida at the Eglin Air Force Base in Fort Walton Beach, Florida. Both professors were members of my supervisory committee for the master's degree. Without their help, I could not have begun my graduate education. So many fond memories are there with them. iii I would like to express my appreciation for the financial support of the United States Air Force for my master's program. I appreciate Dr. David Ebeoglu, Mr. Robert J. Arnold and Dr. David J. Butler of the Armament Division, Eglin AFB, Florida. Dr. Ebeouglu was always there whenever I need help. Mr. Arnold was very kind to me and to my family. The advice and cheerfulness of Dr. Butler helped me all the time. I sincerely appreciate the Agency for Defense Development where I devoted fourteen years of my young life. I extend special thanks to Dr. J. H. Hong, Dr. Y. K. Yoon, Dr. Y. H. Kang, Dr. K. S. Chon, Dr. S. H. Do and many other bosses and colleagues: I wish I could name them all. Finally, I am very grateful to my wife Ki Young and two children, En Jae and En Kyum, who must hold the traveling mileage records in their schools, for their endurance and support. And special appreciations must go to my parents, brothers and sisters for their prayers and encouragement. TABLE OF CONTENTS PAGE ACKNOWLEDGEMENTS ....................................... iii LIST OF TABLES ............. .............................. vii LIST OF FIGURES ....................................... viii ABSTRACT ................................................. x CHAPTERS 1 INTRODUCTION .......... ........................ .. 1 1.1 Introduction ...................................... 1 1.2 Experimental Background ........................ 5 1.3 Objective and Scope ............................. 7 2 FINITE ELEMENT PROGRAMMING OF THE CONTACT PROCESS ................................ 15 2.1 Introduction .... ...................... ........ 15 2.2 Axisymmetric Formulation .................... 17 2.3 Constitutive Relations for an Axisymmetric Laminated Plate ............................... 19 2.4 Algorithm for the Contact Program .............. 25 2.5 Axisymmetric Finite Element ................... 28 3 RESULTS OF FINITE ELEMENT ANALYSIS ............. 44 3.1 Introduction ............................... 44 3.2 Evaluation of the Program ................... 45 3.3 Comparison with Test Results ................ 48 3.4 Comparison to Contact Laws .................. 49 3.4.1 Indentation Law .. ...................... 50 3.4.2 Contact Pressure Distribution .......... 54 3.4.3 Contact Radius ........................... 54 3.4.4 Center Deflection ....................... 55 3.5 Friction Effect ............................ 57 4 DAMAGE ANALYSIS ............................... 78 4.1 Failure Criteria ............................... 78 4.2 Implementing Delamination into Finite Element Program ................................ 88 4.2.1 Introduction ............................ 88 4.2.2 Numerical Experiment ................... 89 4.3 Implementing Fiber Failure and Matrix Failure into Finite Element Program ................... 90 4.4 Incremental Finite Element Method for Damage Analysis ................................ 91 5 CONCLUSIONS .................................. 106 5.1 Conclusion ................................... 106 5.2 Comments for the Future Research .............. 107 REFERENCES ............................................ 108 APPENDIX ................................................. 112 BIOGRAPHICAL SKETCH .................................... 134 LIST OF TABLES TABLE PAGE 2.1 Material properties of AS4/35016 ........... 41 3.1 Material properties of the plates used in. the examples ....................... 60 3.2 Material properties of the numerical applications ...................... 67 4.1 Ultimate strength of AS4/35016 lamina ...... 87 vii LIST OF FIGURES FIGURE PAGE 1.1 Impact test scheme ............................ 10 1.2 Test configuration ............................ 11 1.3 Contact force vs. center deflection curve ...... 12 1.4 Cscan picture of impacted composite plate ..... 13 1.5 A photomicrograph of the impacted plate ....... 14 2.1 The continuity of the stress through the thickness of the composite plate ............... 38 2.2 The continuity of the strain through the thickness of the composite plate ............... 39 2.3 Stacking of plies to one element ............... 40 2.4 Indentation of a plate ........................ 42 2.5 Axisymmetric element ........................... 43 3.1 Plate dimensions in numerical examples ......... 59 3.2 Contact force vs. contact radius relationship .. 61 3.3 Interlaminar shear stresses in the graphite/epoxy laminates ....................... 62 3.4 Interlaminar shear stresses in the hybrid laminates ........................ 63 3.5 Contact force vs. center deflection curve ...... 64 3.6 Contact force vs. center deflection curve, elastic range ................................. 65 3.7 Configurations for numerical examples .......... 66 3.8 Contact force vs. indentation curve, isotropic material ............................. 68 viii 3.9 Contact force vs. indentation curve, composite plate (4 mm thickness) ............... 69 3.10 Effect of plate radius on contact force vs. center deflection curve ....................... 70 3.11 Effect of plate radius on contact force vs. indentation curve ............................. 71 3.12 Effect of indenter radius on contact force vs. indentation curve .................... 72 3.13 Contact pressure distribution .................. 73 3.14 Contact force vs. contact radius ............... 74 3.15 Contact force vs. contact radius, comparison ... 75 3.16 Contact force vs. center deflection............. 76 3.17 Friction effect ................................ 77 4.1 Partial failure patterns ...................... 86 4.2 Delamination simulation ....................... 97 4.3 Numerical application of delamination effect ........................... 98 4.4 Delamination effect on contact force vs. center deflection curve ............................... 99 4.5 Geometry of the contact surface ............... 100 4.6 Contact force vs. center deflection curve ..... 101 4.7 Contact force vs. center deflection curve ..... 102 4.8 Damage pattern of the composite plate ......... 103 4.9 Damage pattern of the composite plate ......... 104 4.10 Damage pattern of the composite plate ......... 105 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy FINITE ELEMENT MODELING OF STATIC INDENTATION DAMAGE IN LAMINATED COMPOSITES BY HAN SOO JUNG AUGUST 1991 Chairman : Bhavani V. Sankar Major Department: Aerospace Engineering, Mechanics and Engineering Science The objective of this research is to study progressive failure of circular laminated brittle matrix composite plates due to static indentation loading using finite elements. The plates are assumed to be simply supported along the edge, and are indented by a rigid indenter with a hemispherical nose. A quasiisotropic homogenization process is used to make the problem axisymmetric. Four noded, isoparametric, quadrilateral, axisymmetric solid elements are used to model the plate. A new algorithm is used to analyze the progressive contact between the indenter and the plate. Hashin's failure criteria are modified so that they can be applied to axisymmetric, transversely isotropic laminates. Stiffness matrices of failed elements are modified suitably. New slip surface is introduced when delamination occurs. Very stiff springs are attached to nodes on delamination surfaces to prevent interpenetration of the nodes. The results from incremental damage analysis agree well with experimental results. This program can be applied to analyze impact problems involving large masses impacting at very low velocities. CHAPTER 1 INTRODUCTION 1.1 Introduction Currently composite materials are used in primary structures of aircraft, particularly aircraft wings and frames, because significant weight savings can be achieved due to their high strength to weight ratio. Also, composite materials can easily be tailored to obtain desired properties in different directions. For example, more 00 layers can be used for longitudinal loads, 900 layers for transverse loads and 450 layers for shear loads. Also, there are other advantages such as resistance to damage by fatigue loading, immunity to corrosion, ease of processing, and low coefficient of thermal expansion. Structural forms that are otherwise inconvenient or impossible to manufacture, e.g., more complex airfoils or helicopter blades, can be easily fabricated. Hence, the utilization of composite materials in aerospace and transportation industries is continuously increasing [1,2]. However, it is well known that brittle matrix composites such as graphite epoxy are very susceptible to low velocity impact, which may be caused by dropping tools, runway debris, hailstones, bird strikes and ground service vehicles hitting 2 aircraft structures. High and medium levels of impact energy cause surface damage that is relatively easily detected, and therefore, repairs may be undertaken. However, low velocity impacts which are at quite low energy levels can seriously damage laminated composite structures internally so that frequently it is very hard to detect damage through visual inspections. The damage can be detected only by using nondestructive inspection (NDI) techniques. For example, low velocity impact to an aluminum plate may cause denting damage and the structure does not lose very much strength. However, inside of an impacted composite plate there may be fiber breakage, matrix cracks, and delaminations, while the impacted surface may not show any damage. Commonly it is termed as "Barely Visible Impact Damage" or BVID [2]. After impact, composite materials lose considerable strength because of internal damage and delaminations. [27] It becomes very important to understand the response of composites to the impact of rigid objects and to predict damage in order to make composite materials more reliable. The study of the effect of BVID on static and fatigue strengths of graphite/epoxy panels shows that impact degradation is of more concern in compression loading than in tension because of delaminations which occur during impact. Importance of the impact study is emphasized by the fact that upper surfaces of the aircraft are more likely to receive impacts from tool 3 drops than lower surfaces, and upper surfaces operate in a compressiondominated load regime [2]. Because of the extremely complicated nature of impact on the composite materials, empirical methods have been widely employed. Generally, empirical methods have several distinct limitations. First, it is very expensive and also time consuming to conduct testing and to process the data. Second, it is practically difficult to test all the combinations of material properties, boundary conditions of the composite structures, and impactor velocities and weight variations. This problem can be critical, if the composite material is in a developing stage, and only limited quantities are available. Third, it is very difficult to comprehend all the stress and strain distributions, and damage progression inside composites by empirical methods. And also there are difficulties in matching results from specimens to the real structures. [6,7] Therefore, numerous analytical methods have been developed to overcome the shortcomings of empirical methods. However, so far, because of the extremely complicated nature of impact problem, there are no complete threedimensional models available which can simulate the entire impact process and predict damage progression. There is an urgent need to develop a three dimensional progressive finite element program which calculates contact forces, contact displacements, stresses and strains inside the composite plate. During the contact process, the finite element program must be able to 4 check all the failure criteria for fiber breakage, matrix cracking and delamination. If failure occurs, then the program could modify the stiffness matrix to account for failures . Fortunately, low velocity impact due to large masses can be treated as a static indentation problem, because the impact duration is much longer than the time required by the propagating waves to travel from the impact site to the supports or free edges, and the test results show good agreement with this assumption. Strictly speaking, the terminology of low velocity impact is not clearly defined. Rather low impact velocity is generally defined as the impact velocity which is low enough to justify a static analysis of the response of the structure, especially in the vicinity of impact [3,4,8,9]. This velocity depends upon the mass, shape and the material property of the impactor and also upon the impacted structure boundary conditions and its natural frequencies. Generally it ranges 13 meter/second with impactor mass less than 10 kg [10]. The advantage of static analysis is that the time variable is eliminated. Hence, we can closely look into progressive damage mechanisms with a small amount of computational resources. If time effects are considered, it will be a formidable task to obtain all the meaningful stress and strain distributions inside the plate and to map the damage pattern and its progression because of computational time and storage limitations of the computer system. If 5 coarser elements are employed, it will be very hard to obtain accurate data because meaningful data must be obtained in the small confined area under the indenter. The experimentally determined contact laws account for local damage in an empirical manner, and hence cannot be used to study the local damage. A fully three dimensional formulation is required to study the progressive damage. In the present study a finite element program performs an incremental analysis of the contact problem, and assesses the damage due to indentation. 1.2 Experimental Background In this section some results from an experimental study [11] are presented to bring out the importance and necessity of the finite element modeling of the damage process due to static indentation. Test specimens were made of Hercules carbon prepreg tape AS4/35016, which is an aminecured epoxy resin reinforced with unidirectional carbon fibers. All specimens were cut from a 30.48 cm (12 inch) by 30.48 cm, [00, +450, 450, 90]8s composite plate layup. In order to obtain the specimens with the quality suggested by the manufacturer, the curing instructions was carefully followed as in Ref. [12]. An autoclave (Baron Model BAC24) was used to fabricate all the specimens of composite laminates used in this study. Detailed 6 procedure of fabrication can be found in Ref. [11]. Material property of the 00 lamina can be found in Table 2.1. Figure 1.1 shows the pendulum impact test facility in which all impact tests were conducted. The pendulum includes the support strings connecting the tup to the laboratory ceiling, and the tup, which is a combination of weights, a load cell, and an indenter with one inch diameter hemispherical head. The total mass of the tup is 13.84 kg. Static indentation tests were done in a standard MTS material testing machine. The composite laminate was placed on a circular steel ring of 5.08 cm (2 inch) diameter to simulate a circular plate with simply supported boundary condition. The test configuration is shown in Figure 1.2. The test results of contact force vs. center deflection of the circular composite plate for the static indentation and impact tests are shown in Figure 1.3. This graph suggests that static indentation tests can simulate low velocity impact tests. Also it suggests that indenting process is highly nonlinear due to internal damages. A sudden drop of the contact force is a unique feature of indentation of composite plates. Figure 1.4 shows the picture of the Cscan of the impacted plate. The circular damage pattern justifies the axisymmetric formulation used in the present study. Figure 1.5 shows the picture of the photomicrograph of the polished 7 cross section of the indented plate. This shows matrix cracks due to shear failure and delaminations. However, these data are not sufficient to understand the entire impact or indentation process. They do not explain how damage initiates and propagates. Also, it is very hard to quantify this damage. 1.3 Objective and Scope The major objective of this research is to develop a finite element program which models the static response of laminated composite plates indented by a rigid spherical indenter and to predict progression of damages inside the plate. The results are applicable to low velocity impact of the composite plates. A circular plate and a rigid hemispherical indenter were chosen in modeling because they correspond to many practical impact situations. Quadrilateral, isoparametric, axisymmetric finite elements were employed in this research. By employing axisymmetric elements, the three dimensional problem can be treated as two dimensional by eliminating dependency on 6. In Chapter 2, the finite element formulation and the incremental displacement method are introduced. The input to the finite element program is the displacements necessary to close the gap between the indenter and candidate contact nodes of the plate. The laminated plate consists of many plies that are orthotropic and have different angles of fiber 8 orientation. As axisymmetric finite elements were employed in this research, there must be special considerations of the material properties and failure criteria to account for this discrepancy. The inplane strains, exx, e and yx are continuous at the interface between two adjacent plies whereas the stresses, a a and are not. On the other hand, zz, cxz and yz are continuous through thickness and are better represented in a homogeneous plate model than ezz, Yxz and yyz, which can be discontinuous through the thickness. A set of transversely isotropic elastic constants is derived for a group of plies that can be considered as a quasiisotropic laminate. In Chapter 3, the accuracy and efficiency of the present method were evaluated by comparing the results with some early work [10]. The contact force vs. center deflection curve of a AS4/35016 graphite/epoxy plate was compared with experimental results in the elastic range. Then, some widely used contact laws were compared with the present finite element results. The effect of friction was studied by assuming stick friction conditions between the indenter and the top surface of the plate. In Chapter 4, the failure criteria are introduced. Among many interactive and independent failure modes, three dimensional failure criteria introduced by Hashin [13,14] were employed in the present studies. These criteria have the same quadratic form as TsaiHill's criterion, but they can distinguish different failure modes such as fiber breakage, matrix cracking and delamination. Since we use an axisymmetric formulation, the damage pattern is also going to be axisymmetric, which does not strictly represent the damage pattern in a composite laminate. However, experimental impact studies [15] have indicated that the actual damage pattern is like a spiralling stair case, and the projection of these delaminations can be considered as circular and hence axisymmetric. This was also confirmed by studies in the laboratory [11]. The extent of delamination, fiber breakage and matrix failure were mapped during the contact process. After failures are detected, the total stiffness matrix of the plate must be modified properly depending upon the failure modes. Matrix and fiber failure can be implemented by reducing or eliminating element stiffness matrices. Implementing delamination into the finite element program requires creation of new free surfaces, consequently new nodes. Rigid springs are inserted across the newly created free surfaces to prevent interpenetration of the nodes. Delamination can be assumed as Mode II (Shear Mode) crack propagation. Incremental damage analysis results are compared with experimental results. A summary of conclusions and some suggestions for future research are listed in Chapter 5. Figure 1.1 Impact test Scheme 11 E E E 2 e o ....... ...... I .......... ... . R = 12.7 mm 7 4 mm Figure 1.2 Test configuration 0 2 4 6 Center Deflection (mm) Figure 1.3 Contact Force vs. center deflection Curve Figure 1.4 CScan picture of impacted composite plate Figure 1.5 A photomicrograph of the impacted plate CHAPTER 2 FINITE ELEMENT PROGRAMMING OF THE CONTACT PROCESS 2.1 Introduction The finite element method has been established as a powerful and versatile tool of analysis for contact problems involving isotropic bodies. In the present study, a method for solving the problem of indentation of composite laminates is developed. The contact problem of quadratic surfaces of elastic bodies was first solved analytically by Hertz in 1881. Since then, numerous contact problems have been solved by formulating the problem as an integral equation. The solution is either obtained by solving the equation numerically or, for simple configurations, in closed form. The Hertzian type problem has three general assumptions. First, the contacting solids are homogeneous, isotropic and linear elastic. Second, the contact areas are essentially flat and relatively small compared to the various radii of curvature of the undeformed bodies in the vicinity of the contact interfaces. Third, the contacting solids are perfectly smooth and only normal pressures are considered [6,16]. 16 Hertzian type solutions are not suitable for indentation of composite laminates which are inhomogeneous, anisotropic and also involve significant damages at low load levels. Further, most laminated composites in use can not be adequately represented as a halfspace as in Hertzian contact problem. Because of these complexities, three dimensional analysis, especially three dimensional finite element methods, are required in the present study. In the last two decades, many researchers thought that contact problems could be considered as a special case of constrained minimization of either total or complementary potential energy. The minimization is formulated as a mathematical programming problem, and the solutions are obtained by either using incremental linear programming, quadratic programming or conjugate (modified) gradient methods. The minimization problem can also be formulated as a variational inequality using penalty methods or Lagrange multipliers [1721]. However, these methods have certain disadvantages such as creation of flexibility matrix, the awkwardness involved in introducing additional unknowns, such as slack variables and Lagrange parameters. Nothing can be better than direct formulation of finite element method as a simple design tool, which calculates the detailed stress field in various types of laminates for various boundary conditions accounting for the damages inside the plate. The goal is to use simple and conventional elements so that the method can be 17 implemented in commercial finite element programs. The novel incremental contact algorithm developed in this study is explained in Section 2.4. 2.2 Axisymmetric Formulation The benefit of employing an axisymmetric element is that the three dimensional problem can be treated as two dimensional. Especially in the case of rigid ball indentation, where the contact surface is circular, axisymmetric formulation is a very efficient method to treat the indenting process in a computer program, because upper surface nodes of the plate come into contact successively and the contact area is readily determined. If we use three dimensional solid elements, it will be a formidable task to determine the contact nodes for the case of a spherical indenter, although the employment of three dimensional solid elements has a distinct advantage in that we can use orthotropic material properties and also the failure criteria for orthotropic laminates directly without any modifications. Complexity of the contact algorithm demands employment of axisymmetric elements, and consequently there must be some considerations to convert the material properties of a laminated plate, which are basically orthotropic, into equivalent axisymmetric properties. The circular damage distribution pattern of the C scan pictures of the damaged plates (Figure 1.4) justifies the employment of axisymmetric elements to solve the problem. 18 Stressstrain relationships of the axisymmetric element are discussed extensively in Section 2.3. In the present study u(r,z) and w(r,z) represent the radial and transverse displacements, respectively. The axisymmetric stress and strain components are { o } = { ra zz 'rz ee } (2.1) (e} ={ ez ezz Yrz eee } (2.2) where, uw err r zz W (2.3) au aw Y + v = 8z + Br a" z ar considered as elastic until failure, material nonlinearities As graphite/epoxy composite is very brittle and can be considered as elastic until failure, material nonlinearities are ignored [2]. As will be seen later, the deflections considered are smaller than the plate thickness, and hence geometric nonlinearities are also not considered. 19 2.3 Constitutive Relations for an Axisymmetric Laminated Plate A composite laminate can be considered as having several sublaminates. If each sublaminate is quasiisotropic, then each sublaminate can be idealized as a homogeneous transversely isotropic material. Then the composite laminate can be considered as a laminate made up of several transversely isotropic layers. Hence the axisymmetric formulation can be employed if the external loading and boundary conditions are also axisymmetric. The stressstrain relationship in the Cartesian coordinates are {e} =[S]{o (2.4) where ( o } and ( e ) are, { e } = { e eyy Y~ ezz Yzx Y~ } T (2.5) { a } = { Oxx ayy, xy Ozz Tzx rzy }T (2.6) For the ply which has fiber orientation angle 8 from the x coordinate, lamina constitutive relations can be obtained from {e} = [ T][S][ T {o }e (2.7) where [ S ] is 1 V21 0 v3 0 0 V12 1 0 320 0 E1 E2 E3 1 0 0 0 0 0 C12 0 o o 3 230 1 0 0 E E2 EE 1 0 0 0 0 0 GM o o 0 0 0 G23 and [ To ], [ Te ] are given by, cos20 sin28 2cos6sinO 0 0 0 sin28 cos2 2cosisinO 0 0 0 cos6sinO cosOsinO cos28sin2e 0 0 0 0 0 0 1 0 0 0 0 0 0 cos6 sin6 0 0 0 0 sine cos6 cos2) sin28 cosfsin6 0 0 0 sin28 cos28 cosssine 0 0 0 2cos6sin6 2cos6sinO cos28sin28 0 0 0 0 0 0 1 0 0 0 0 0 0 cosO sine 0 0 0 0 sin6 cosO For the fiber orientation angle 6, { e = [ S]e {o }e [ S ]e = ] [ s ]e[ S]T ]T Averaging [ S ]e for the number of plies from this form is not suitable because inplane stresses such as axx, ao and Txy are not continuous through thickness of the element while exx, e and Yy are continuous. And, also, outofplane stresses ozz, zx and Tzy are continuous through thickness while outofplane strains ezz, yzx and yzy are not continuous. These are well illustrated in Figure 2.1 and Figure 2.2, which Wu and Springer have obtained [23]. These figures show the continuity of the stresses and strains through the thickness of the impacted plate. The best way to average the modulus of the element can be done after denoting [ S ], as, P q [S], = (2.8) r s where p, q, r and s are 3 x 3 matrices. Then the stressstrain relationship can be written as (xx exx ayy yy x p q (2.9) oz [ r s ezz tzx zx tzy zy The stressstrain relationship can be rearranged such that all the components of stress and strain that are continuous through the thickness are on one side and the remaining components are on the other side as shown in Figure 2.9. ax e a e ezz rp srpq zz Szy zy Let us denote [ D ]e [D ], = (2.11) Srp srpqj Then an average [ D ] matrix for a group of N plies can be obtained as, [D]av= [D]e, i=l (2.12) RS Now the elastic matrix [ C ], can be recovered from [ D ]av. The [ C ]av can be expressed in terms of P, Q, R and S as, PQS1R QS1 (2.13) [ e S R S1 The stressstrain relations take the form { 0 }av = [ C ] { e }a For widely used cross ply or quasiisotropic laminates, such as [ 0 90 ]s or [ 00 ,+ 450 900 ]i laminates, [ C ]av becomes transversely isotropic. For axisymmetric problems, where are = 16e = re = e3 = 0 (2.14) 24 "Crr Ce 0 Crz C]e Ce Cee 0 C (2.15) 0 0 C o rz cz C, 0 Cz Further, the present homogenization scheme yielded a symmetric stiffness matrix, i.e., Cij = Cji, and also c'r = ce Crz = Cez (2.16) Sometimes, five engineering elastic constants, such as Er, Ez, Grz, Vre, Vrz can describe the material properties of a transversely isotropic laminate plate [10]. For this case, the stressstrain relationship can be obtained as follows. 1 v 1 1 Vre 0 0 Z Vre 1 0 Vzz E, E E [C]=r E E r (2.17) 1rz V V 1 0z rz 0 Er E, E, r r 25 Table 2.1 shows the material properties of the AS4/35016 graphite/epoxy lamina. Also, the calculated transversely isotropic element material property matrix [ C ],, for a plate, which consists of four plies in [ 00 450 900 ] orientation is given in Table 2.1. This configuration is chosen because of its wide application in aircraft structures. For example, the main wing box for a straight wing aircraft, 00 plies carry the spanwise direct stresses induced by wing bending like a cantilever beam bending, 450 plies carry shear stresses caused by wing twisting, and finally 900 plies carry the chordwise direct stresses due to fore and aft bending of the aircraft wing [2]. 2.4 Algorithm for the Contact Problem Unlike conventional finite element problems, contact problems are mixed boundary value problems where nonzero displacements are specified at parts of the boundary and zero tractions are specified in the rest. Further, the displacements in the contact surface have to satisfy certain conditions specified by the geometry of the indenter. Sankar [10] used a contact algorithm in conjunction with a finite difference method. In the present study, the algorithm has been modified to suit the finite element method. Unlike the finite difference method, the finite element method has the advantage of modeling complex geometries and multilayer media. 26 For the frictionless contact problem, the vertical degrees of freedom of the top surface nodes of the plate become contact candidates. In Equation (2.18), the displacements of these contact candidate degrees of freedom are notated as du and other displacements as dr. The dr degrees of freedom can be condensed to obtain a smaller stiffness matrix as shown below. [K. Kur du Fu (2.18) .ru Krr dr 0 { d } = [ Kr ]1 [ Kru ] { du } (2.19) After condensation, the equilibrium equations take the form [ Kc ] { d = { Fu } (2.20) where [ KC ] = [ K,] [ K ] [ K ][ Kru (2.21) Since we are dealing with an axisymmetric problem, the contact area is always circular. Let us assume that the contact has progressed to Node k1 (see Figure 2.4). The 27 procedure for determining the incremental load required to bring the kth Node into contact is as follows. Equation (2.20) can be partitioned as K.= Kct w. Fe (2.22) Kfe Kff, wf 0 J where wc are the z displacements of nodes within the contact region including k1 and wf are the remaining degrees of freedom for which the forces are zero. Let us use gi to denote the gap between the indenter and nodes outside the contact region, i.e., i 2 k. The first step is to determine the incremental displacement of the indenter to close the gap gk. We first apply a unit displacement to the indenter, i.e., Wi = 1, i = 1,2, ,kl (2.23) The equation (2.22) is solved for Fc and wf. Now the new gap g'i becomes g'i = gi 1 + w., i t k (2.24) Now we can calculate the displacement 6 of the indenter to make the gap g'k identically equal to zero. 6= gkk gk g k (2.25) gk 1 Wk ) After determining 8, the nodal force FC can be computed. By integrating all the nodal forces, we obtain the incremental contact force, AP. The procedure is repeated for the next contact node k+l. When frictional effects are included, radial forces at nodes in the contact region can not be considered as zero. The procedure has to be slightly modified, as is explained in Section 3.5. The damage alters the stiffness matrix. We assume that the indentation process is under displacement control, and hence the wdisplacements of the contact nodes are kept the same when passing from one iteration to the next. The changes in [ K ] cause a change in the contact load (usually reduction) for the same indentation displacement. Detailed damage analyses are presented in Chapter 4. 2.5 Axisymmetric Finite Element In this section, the derivation of the stiffness matrix of a quadrilateral, isoparametric, axisymmetric element is given for the sake of completeness. 29 In isoparametric elements, the element coordinates and the element displacements are expressed with the same shape functions associated with nodal coordinates and displacements respectively as shown below. n r=C N= (Ej') ** i=l n (2.26) u = N (,1) ui i=l n w = Ni (in) w_ i=1 where r and z are element coordinates, & and n are the natural coordinates and u and w are element displacements as shown in Figure 2.5. Subscript i signifies node number. Shape function or interpolation function Ni are expressed in the natural coordinate system, which has variable & and t that each vary from 1 to +1 as shown in Figure 2.5. The characteristic of the shape function Ni is that its value is unity at node i and is zero at all other nodes. For the quadrilateral element, the Ni are as follows. IN, = ( 1 + E ) ( 1 + r ) 4 N2 = ( 1 ) ( + (2.27) V = (1 + ) (1 ) 4 To evaluate the stiffness matrix of the element, we need to calculate the straindisplacement transformation matrix for the small strain case as, Su err = ar 8w Su aw ee = r As the interpolation functions, N,, are defined in the natural coordinate system by Equation (2.27), we need the relationship between r and z derivatives and t and q derivatives. By using chain rules, the following expressions can be obtained. au zar z1 au au I a aj tW I= 1 at (2.28) au ar az au au z 9 z] owI ar} = as, asF an = as (2.29) aw ar az aw aw BT =W RW 1 a1 (2.29) aw r az aw aw Fz _F _T1 anrl In the above equations, J is the Jacobian operator relating the natural coordinate derivatives to the local coordinate derivatives. The inverse of Jacobian matrix must exist in order for the isoparametric element to be used. This means that there is a onetoone correspondence between the natural and the local coordinates of the elements. In actual computations, the sign of the Jacobian determinant is monitored at the Gaussian integration points to ensure the existence of the inverse Jacobian. Generally, a zero or negative value of the Jacobian determinant indicates input data error or an overly distorted element, and the computations are terminated automatically for the purpose of node renumbering [22]. Equations (2.28) and (2.29) can be calculated as follows. ar ( 1 1 1 aE 4 4 4 4 r _1( 1 1 1 ar 1(l+)r1 (IF) r2 (l) r3 (l+) r, 4 4 4 4 (2.30) 1 (1+ ) z (1+1) Z2 (1n) z3+! (11) Z4 z 1 (+)1 1 +1  47 4 4 4 4 (2.31) u= _1 (+()u 1 1 u+ au +4 1 (1, u1 (1+ ) u2 1 (11)U3 (1n) u4 au 1 (1+4) u1 + (1&) u2 (l) u3 (+ ) U4 4 4 4 4 (2.32) aw = (1+ ) w 1 (1+11 2 (1 )3 a 4 4 4 (111) w4 1 1 1 aw= (1+) w+ (1) (1 W3 (1 ) w4 l 4 4 4 4 (2.33) Eventually, the unknown variables of the system will be the nodal displacements { U } { U} = { UI w 2 W2 u3 w3 u4 w4 }T (2.34) Then, f 9U. , Iau u 4 1+ o 0 i o (+) o ((1+ o0 0 Oz au 7 + 0 o 0 (1d) 0 (1+n) 0 (U aw r P_ 1j 01 + 1 (1+l) 0 (1) 0 lrT aw 4 0 1+t 0 1 0 (1t) 0 (1+) az Strain can be expressed with nodal variables { U ), and also stress can be obtained by multiplying the strain vector by the elasticity matrix [ C ]e { e = [B] { U} (2.35) { o } = [ ]e [B]e { U} In the axisymmetric case, unlike other isoparametric elements, the straindisplacement matrix [ B ]e has terms containing r. Thus, for two elements of identical shapes and material properties but at different r locations, the [ K ]e will be different. Also, a problem arises when r is zero because some terms have r in the denominator. A proper numerical scheme, such as Gaussian quadrature, can avoid this problem [24]. In the axisymmetric case, { e ) becomes { e } = {err ezz Yrz eee}T (2.36) and [ B ]e becomes 4 by 8 matrix as follows. 1+8 0 1++ r 0 1i+ 1+1 0 (0 1 & (1&) (1+11) (1+ ) o i + r 1 0 1( (i+n) 0 (2.37) 4 _(i_ ) 0 ( S) 1 ) ( + ) r o (1S) (1a) 0 1n 0 (1+0) r 0 (1+) 1r 0 The element stiffness matrix can be formulated by using the virtual work principle. The principle requires that for a body under equilibrium, f 8eTodV = f 8aTfBdV+ f a'fsdS + 8aTfN (2.38) vs v, S, where a ; actual stress be ; virtual strain that corresponds to the imposed virtual displacement 8a ; virtual displacement fa ; actual body force fs ; actual surface force f ; actual nodal force for the axisymmetric element, the volume integral of any function A becomes f A dV= f A rdrdzdO (2.39) Volume integration can be written in the natural coordinate system. Notice that r, which is a global coordinate, appears in the integration, which is a unique feature of the axisymmetric element. f A d = 2nff A r det [J] dd4 (2.40) If there are no body forces, or initial stresses, Equation (2.38) becomes 6aT [ B ]e[ C ][ B ]e dV { U} = aT fX (2.41) Ve by applying { 8e } = [B ]e { U} { ao = [ C]e { = [ C ] [ C] [B ] { U} (2.42) {6a } = [ N] {BU} The stiffness matrix is obtained as, [K] = [B ][ C][B ]e dV (2.43) V. It is very easy and convenient to use GaussLegendre numerical quadrature integration formula to evaluate Equation (2.43). The GaussLegendre integration points Ei and ij and corresponding weighing factor aij are found in Ref. [24]. Care must be taken to substitute r in Equations (2.37) and (2.43). From Equation (2.26), r can be obtained as 1 1 r = (1+Ei ) (1+Il)r. + (1+) (1+,).r2 (2.44) 1 *+ (iS) (1 .)i r3+ (1+i) (j ) r4 4 4 After substituting Ei and mj into E and at the previous equations, [ K ]e can be obtained numerically as 37 [K ]e = [ B ]i[ Ce[ B aijdet[ J ]i,r i,j = 1 (2.45) The assembly of element stiffness matrices, applying boundary conditions and the solution of resulting linear equations follow the standard procedures that can be found in Ref. [24,25]. 20 .. 20 .0.04 0 0.04 x, COORDINATE (in) x, COORDINATE (in) Figure 2.1 The continuity of the stresses through the thickness of the composite plate, Ref. [23] 0 0.04 0.001 'KI 0 0.001 0.001 0.001 0 0.001 .0.04 0 0.04 x COORDINATE (in) O.OOS i 0.005 O.OOS 1' ;  . , .OO 0.005 0 0.005 0.005 0 0.005 0.04 0 0.04 x, COORDINATE ( in ) Figure 2.2 The continuity of the strains through the thickness of the composite plate, Ref. [23] F SCF One Element Consists of Four Plies 0 Degree 45 Degree 45 Degree 90 Degree Laminated Plate Figure 2.3 Stacking of plies to one element Table 2.1 Material Properties of AS4/35016 (unit; Gpa) E, 144.8 E2=E3 9.7 G12=G13 6.0 G23 3.6 V12= 13 0.3 23 0.34 * Courtesy of Hercules Advanced Materials and Systems Company Axisymmetric Element Property Matrix, for [ 00, 450, 900 ] [C ], (unit ; Gpa) 63.593 4.178 4.178 11.105 0 20.198 0 4.178 0 4.500 4.178 0 63.593 20.198 \F a; Radius of the plate c; Contact radius R; radius of the indenter Gk; Gap between indenter r,u k c k k IGk N Figure 2.4 Indentation of a plate .h.( ,r, u One Element Z,W r 2 =1 1 3 4 1 (= 1 2 1 3 =1 Local Coordinate z Global Coordinate Axisymmetric Element Figure 2.5 CHAPTER 3 RESULTS OF FINITE ELEMENT ANALYSIS 3.1 Introduction In this chapter, the accuracy of the finite element program is evaluated by comparing with the results of Sankar [10], who employed a finite difference method to solve the equilibrium equations. The results obtained using only 384 finite elements (24 x 16) in the present study compared very well with that of Ref.[10]. In Section 3.2, available methods of increasing toughness and damage tolerance are briefly discussed. In Section 3.3, the finite element program was applied to the experimental test configuration shown in Figure 1.2 to obtain contact force vs. center deflection. Without modification of the total stiffness of the plate due to internal damages, finite element results seem to deviate from the experimental results. These deviations can be attributed to various damages as will be discussed in Chapter 4. However, in the elastic range, where damage is minimal, finite element results match the test curves well (Figure 3.6). In Section 3.4, various empirical equations including the Hertzian indentation law are briefly introduced and compared with finite element results in the elastic range. In Section 3.4.1, 45 the finite element program is used to study the indentation of solid cylinders and plates to compare with Hertzian indentation law. Effects of plate radius and indenter radius on indentation and deflection of the plate are studied. In Section 3.4.2, Hertzian contact pressure distribution, which is known to be elliptic, was compared with finite element results. In Section 3.4.3, the relation between contact force and contact radius was studied. It is known that in a half space contact force varies as the cube of the contact radius. Finite element results matched very well with the available empirical equations. In Section 3.4.4, the relation between contact force and center deflection was studied. In Section 3.5, the effect of friction between indenter and plate was considered. 3.2 Evaluation of the Program Hybridization, such as mixing the fiber types, thus utilizing the advantageous properties of one material to overcome deficiencies of another, is one method to deal with impact damages. The higher strain to failure of glass fiber or Aramid fibers (Kevlar) led to substantial increases in threshold energies for impact damage, with the energy level required to produce delamination being increased. The Aramid/graphite/epoxy hybridization is promising because Aramid material is relatively inexpensive and it increases impact resistance while graphite fibers provide compressive 46 strength, which Aramid fibers lack. Tests showed that doubling of residual strength after impact was possible using Aramid/graphite hybridization compared to an all graphite specimen [2]. It is well known that tough resins can significantly reduce the damage caused by impact and substantially improve the residual strength following impact. This is because the onset of delamination and the delamination growth in a composite depends on its interlaminar fracture toughness, which in turn depends on the toughness of the matrix material. Widely used epoxy resins are brittle, but crystalline thermoplastics such as PEEK (Polyetheretherketone) can be used in place of epoxy to increase toughness. Mode I fracture toughness of the graphite/epoxy ranges 80200 J/m2 while the fracture toughness of graphite/PEEK is about 2500 J/m2 [2] Interleaving is a process in which thin layers of thermoplastics are placed in between layers of thermosetting prepreg materials to form a tough interlayer. It is useful if the failure occurs through delamination. It is not effective if the transverse shear is high and, for this case, an improvement of shear modulus of the resin is more important [2]. Through thickness reinforcements, such as a three dimensional braiding process or throughthethickness stitching, can improve impact resistance. However, inplane 47 strength will be reduced because of wavy fibers and irregular fiber orientations. Sankar [10] performed the analysis of indentation of hybrid and interleaved plate problems by using a finite difference method starting from the axisymmetric equilibrium equations. C11 ( ,rr + r1U, r r2 ) + C44 U,z + ( C13 + C44 ) W, = 0 ( C13 + C44 ) ( u,Irz + r1u,z ) + C33W, z + C44 ( w,I + r1W, r ) = 0 (3.1) (3.2) where {rr C11 (aee = C12 ozz C13 C13 rr C13 ee C33 ezz trz = C44 Yrz (3.4) The results obtained by the present finite element program coincide very well with results in Ref. [10]. Figure 3.2 shows contact force vs. contact radius relationship for (3.3) 48 the three different laminated plates. Material properties of the composite laminates are listed in Table 3.1. Figure 3.3 shows interlaminar shear stress in the homogeneous graphite/epoxy laminate. Figure 3.4 shows interlaminar shear stress in the hybrid laminate. Employing the finite element program is very economical because only 384 (24 x 16) elements were used, and computing time was about 30 minutes in a Vax main frame computer. Although this analysis does not account for the internal damages of the plate, it provides information about the possible onset of delamination and its subsequent propagation. As delamination failure occurs where shear stress is high, shear stress distributions such as Figures 3.3 and 3.4 suggest locations where delamination failures possibly may initiate. By obtaining other strain and stress distributions inside plate, possible damage pattern, such as matrix crack and fiber breakage, and its extent can be predicted conservatively and less expensively before detailed damage analysis. Strengthening by inserting tough materials or hybrid laminates to minimize delaminations can be studied by using the finite element program. Detail damage analyses are presented in Chapter 4. 3.3 Comparison with Test Results The finite element program was used to analyze some of the specimens used in the indentation tests discussed in Section 1.2. The contact force vs. center deflection curve was 49 obtained without considering damages inside the plate and was compared with the test results in Figure 3.5. Even in the presumably elastic range, finite element results make the plate appear stiffer. However, from Hashin's failure criteria, which will be discussed in Section 4.1, the damages begin to accumulate at about 140 N contact force or 0.33 mm center deflection. Hence the actual elastic range is very small. Comparing the results only in this elastic range, results of finite element and experimental data showed fairly good agreement. A total of 512 elements (32 x 16) were employed. Fine mesh was used near the center of the plate, and coarser elements were employed outside the contact region. Actually, finite element results form an upper bound for the test results. It can be easily explained that these early deviations of the stiffness of the test results are due to the inherent imperfections in the test specimens such as voids, debonds and nonuniform fiber orientations. 3.4 Comparison of Contact Behavior In this section, various results concerning indentation, contact pressure distribution, contact radius and center deflection obtained using the present finite element method are compared with available results. 3.4.1 Hertzian contact law Because of the extremely difficult nature of analytical solutions of contact problem, which involves permanent deformations and damage such as matrix cracks and delamination, an experimental approach has been taken to determine the contact law of composite materials. Naturally, loading and unloading curves (contact forceindentation relation) are different because damages occur even at low load levels. This is because composite materials, unlike monolithic materials, have very low resistance to the damage due to loads transverse to the fiber direction. For the loading curve of contact between an elastic sphere and elastic half space, the Hertzian contact law [6] is 3 F =K a2 (3.5) where Kc is 4( 1 vs )2 ( 1 v )2 Kc I[ E + E1 (3.6) Subscript s denotes the indenter and t denotes the target. In Equation (3.6), RS, v and Es are the radius, the Poisson's ratio and the Young's modulus of the indenter, respectively. When a metal indenter like a steel ball indents a polymer 51 composite, the indenter can be assumed rigid because Young's modulus of steel is much higher than Ez, the modulus of elasticity of the composite in the transverse direction. For loading, Sun and Yang [6] suggested the contact force indentation relationship in the form = a 2 (3.7) A more complex form of contact law obtained by Conway [26] for an elastic spherical impactor striking a transversely isotropic composite medium has only 3 percent difference from Equation (3.7) in the example problem to be presented in this section. In Equation (3.7), C44 of [ C ] matrix can be replaced by Ez without significant error (4 percent difference in this example problem). The unloading formula suggested by Sun and Yang [6] was K F= c (aao)25 (3.8) cr where a0 = a [1( acr)] am > ac am (3.9) a = 0 am acr 52 In Equation (3.9), am is the indentation corresponding to the maximum contact force before unloading. It is assumed that the critical indentation, ac, depends only on the material properties. Contact force vs. indentation graphs were obtained by the finite element program for the clamped solid cylinder of 20 mm radius and 20 mm thickness of a hypothetical isotropic material, and also for AS4/35016 composite laminate plate of 25.4 mm radius. The material properties are listed in Table 3.1. A total of 512 elements (32 x 16) were employed. Inside the contact region, fine mesh is used, while coarser mesh size is assigned outside the contact region. Each computer run took 75 minutes on the average in the Vax main frame computer. In the case of indentation of the thick cylinder, the center displacement at the top becomes indentation directly. For the plate, indentation is defined as displacement difference between top and bottom surface of the plate at the center, a = w(0) w(h) (3.10) The results for contact forceindentation are presented in Figures 3.8 and 3.9. For the solid cylinder of isotropic material, the finite element results and Hertzian indentation law of Equation (3.6) agree very well as seen in Figure 3.8, although errors were expected because we have not solved the half space contact problem, and the indenter can not be 53 simulated exactly as a sphere in a finite element program. For the plate of isotropic material, contact forces predicted by the finite element model are larger than forces calculated by Hertzian contact law when the indentation is large (Figure 3.8). Sankar [27] noticed this difference in the contact coefficients of Equation (3.5) between contacts of a sphere with a half space and with a plate. He suggested that differences of the definitions of indentation between half space and plate may be the cause of the problem. The finite element contact forces of the composite plate (4 mm thickness) are shown in Figure 3.9. To check the influence of the radius of the plate, 25.4 mm, 38.1 mm, and 50.8 mm radius plates are compared in Figures 3.10 and 3.11 (the indenter radius is 12.7 mm). Although the contact force vs. center deflection curve shows a big difference due to the bending effect of the plate (Figure 3.10), contact force vs. indentation curve does not vary so much (Figure 3.11). At small indentation, contact force seems to be independent of thickness and radius of the target. To verify that contact force is proportional to the square root of indenter radius, two different indenter radii, 6.35 mm and 12.7 mm, were used to indent a 4 mm composite plate (Figure 3.12). The contact force curve obtained using the finite element program for the indenter radius 6.35 mm was denoted as F1 and the curve for 12.7 mm indenter as F2. It may 54 be seen that for the small indentation, the contact force is proportional to the square root of indenter radius. 3.4.2 Contact pressure distribution of the plate The contact pressure distribution for halfspace indentation is as follows [27]. S=3 P (r/C)2 (3.11) 2 3t C2 where P is contact force, and C is contact radius. Finite element results were compared with these equations in Figure 3.13. Obtaining contact stress distribution by experimental test is not easy and no data is available for the tests done as in Section 1.2. Contact pressure distribution calculated by the finite element method in the 4 mm plate (Figure 3.13) is elliptic. Contact pressure of the finite element method was computed from the average vertical stress oz of the elements in contact with the indenter. Although calculated contact forces were larger than predicted by the halfspace formula, Equation (3.11), the contact pressure distribution is still elliptical. 3.4.3 Contact radius of the plate In a halfspace the contact force is proportional to the cube of the contact radius C. By inserting an empirically 55 obtained relationship between indentation and contact radius for the indenter radius as in Equation (3.12) into Hertzian contact law, Equation (3.7), we can obtain Equation (3.13) [28]. C2 a (3.12) Rs = z C3 (3.13) 3RS Finite element results for a 4 mm plate are compared with the above equation in Figure 3.14. In Figure 3.15, contact force vs. contact radius curves of plates of different thicknesses are compared. As plate thickness becomes larger, more contact force is needed to maintain the same radius of contact, as expected. However, for the small contact radius, contact force can be assumed to be independent of thickness of the plate as suggested in Equation (3.13). 3.4.4 Center deflection of the plate The center deflection 5 of the impacted surface of the plate can be expressed as the sum of indentation and of displacements due to bending and shear. In the case of small deflection, membrane effects can be neglected. Indentation was 56 obtained from Equation (3.7). Plate bending stiffness, Kp, can be expressed as [29], 16x (l+v)D (3+v)a2 (3.14) for the axisymmetric plate with thickness h and radius a. In Equation (3.14), D is (3.15) Eh3 D 12 2) 12 (1v2) where, E can be replaced by C1 of the matrix [ C ]av and v, Poisson ratio, can be replaced by v12. Shear deflection can be obtained from aw P r = G  zz ar 2trh Plate shear stiffness, Ks, can be obtained as . 2x Gh [ Ln a Ln C] and C can be obtained from Equation (3.12). Then, the center deflection 6 is expressed as (3.16) (3.17) S1.5 + P + P (3.18) KK, K, Finite element results of a 4 mm plate are compared with Equation (3.18) in Figure 3.16. It is shown that Equation (3.18) can predict good results for the elastic response of the plate. This result suggests that center deflection is a more reliable parameter than indentation to obtain the contact force. 3.5 Friction Effect The effect of friction is studied by assuming that the coefficient of friction is very large. There are only compressive forces between the upper surface of the plate and the indenter. Stick condition requires that there be no relative motion during a load increment between indenter and plate at the contact nodes. This condition is satisfied by constraining the radial displacements of the nodes after they come into contact with the indenter [3031]. This can be easily implemented in the finite element program by modifying the contact algorithm which was developed in Section 2.4. Now, displacements in Equation (2.22) in Section 2.4 must include radial displacements of the upper nodes. In Equation (2.22), wC is replaced by { wk } = { u1, w,, u1 ,2, .... Ui, W1} T (3.19) 58 Treating the vertical displacements in the same way as in Chapter 2.4, we now put zero for the radial displacements in Equation (3.19) during contact. Radial force is very large compared to vertical force because of the high modulus of radial and circumferential stiffness due to the presence of fibers. Hence, matrix failure may occur virtually immediately after contact starts. Calculated contact force vs. center deflection for both frictional and frictionless contact are shown in Figure 3.17. The effect of friction is significant at higher contact forces. However, the radial stresses developed may cause failure even at smaller loads, and the result shown in Figure 3.17 may not be applicable in practical situations. F Rs h C z Rs; Radius of indenter, 10 mm h; Thickness of plate, 2 mm C ; Contact radius a ; Radius of plate, 50 mm Figure 3.1 Plate dimensions in numerical examples Table 3.1 Material properties of the plates used in the examples [10], (unit, GPa) Material Er Ez Grz vrz V r 1 50 7 3.8 0.285 0.265 2 20 8 3.6 0.270 0.255 3 3 3 1.1 0.350 0.350 1. Homogeneous laminate plate material 1, thickness 2 mm 2. Hybrid laminate plate material 1, thickness 0.25 mm material 2, thickness 0.25 mm material 1, thickness 1.5 mm 3. Interleaf laminate plate material 1, thickness 0.25 mm material 3, thickness 0.25 mm material 1, thickness 1.5 mm 3E01 5E01 7E01 9E01 Contact Radius (mm) Figure 3.2 Contact force vs. contact radius relationship 1E.00 9E01 8E01 7E01 6E01 5E01 4E01 3E01 2E01 1E01 OEO00 62 180 170 160 150  S 140 130 120 S 110 L 100 c/h 0.5 L 90 0 80 (0 70  To c/h 0,376 60 60  c/h 0.25 40  30  c/h 0.125 20  10 "0 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Z Coordinate (mm) Figure 3.3 Interlaminar shear stresses in the graphite/epoxy laminates 170  160 150  S, c/h0.5 0 140  a. 130  P 120  S 110 L c/h0.375 S 100  L 0 03 80 ) 70  60  c/h0.26 60 40  30  c/h0.125 20  10  0  0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 Z Coordinate (mm) Figure 3.4 Interlaminar shear stresses in the hybrid laminates 0 0.2 0.4 0.6 0.8 1 Center Deflection (mm) Figure 3.5 Contact force vs. center deflection curve S i / r / .u *.L // .'I Ur~" 0.02 0.04 0.04 I I 0.08 * F.E.M. ..... Impact Test(IMPCB) ....... Static Test(STC1O)  Static Test(STC 13)  Statcl Test(STC 14) S .08 0.08 0I 0.1 Center Deflection (mm) Figure 3.6 Contact force vs. center deflection curve, elastic range I I / ., ,  260 240 220  200  180  180 140  120  100  / / I / 66 Rs: Radius of Indenter ^ a .\ r Solid Cylinder Rs = 12.7 mm a = 20 mm z SPlate r St Figure 3.7 Configurations for numerical examples 67 Table 3.2 Material properties of the Numerical applications Isotropic Material E = 2.0 GPa v = 0.25 G= 0.8 S E(1v) S (1+v) (12v) 2.4 0.8 [ C] = 0.8 v 1 lv v 1v 0 0 12v 2 (1v) v v 1v 1v 0.8 0 2.4 0 0 0.8 0.8 0 0.8 0.8 0 2.4 Composite Material Er = 56.427 GPa Ez = 10.688 GPa Grz = 4.5 GPa rz = 0.263 v = 0.3 63.593 4.178 [C]av= 0 20.198 4.178 0 20.198 11.105 0 4.178 0 4.5 0 4.178 0 63.593 GPa v 1v v 1v 0 1 150 140  130 0 F.E.M.(4 mm Plate) S120 0 F.E.M.(Solid Cylinder) c B S 110  Eq. 3.6 Z 100  90  o 90 o 80 I. S 70  o0 50  40  30  20  0 0 0.02 0.04 Indentation (mm) Figure 3.8. Contact force vs. indentation curve, isotropic material 69 150 140 130 N F.E.M. 120 Eq. 3.7 O 110  100  Z .90 70 o U 80 480  30 270 x ,6 80  4 0 E W 10 NW 20 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 Indentation (mm) Figure 3.9. Contact force vs. indentation curve, composite plate (4 mm thickness) 70 400 350  300 z / 0 250 ' O IL S 200 0 o 0 150  a = 25.4 mm / a = 50.8 mm so Zr  50i ^ ~ i i i i  i i i  0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 Center Deflection (mm) Figure 3.10. Effect of plate radius on contact force vs. center deflection curve 400 350  / S a = 25.4 mm / / o 300  ... a = 38.1 mm 5  a = 50.8 mm / 260  O /. o // U. /," LL 200  0 /SOI 150 50 0  I I I ti I 0 0.004 0.008 0.012 0.016 0.02 0.024 Indentation (mm) Figure 3.11. Effect of plate radius on contact force vs. indentation curve 400 350  S F1 (R = 6.35 mm) 3 ........ F2(R= 12.7 mm) 300 / z F2F1 S 250 0 / ,( 200 / r /" 150 100 0 0.004 0.008 0.012 0.016 0.02 0.024 Indentation (mm) Figure 3.12. Effect of indenter radius on contact force vs. indentation curve 0 .0.1 0.2 0.3 .0.4 Radial Coordinate r (mm) Figure 3.13. Contact pressure distribution 0 0.2 0.4 Contact Radius (mm) Figure 3.14. Contact force vs. contact radius 150 140 130 120 110 100 90 80 70 60  50  40  30  20  10  0 i, , I!  4 mm plate //  3 mm plate /  2 mm plate 1 mm Plate / i,' i Eq. 3.13 // / :'//, ,I '// / / / //// /w ,1// ./// '/./ / 4. oo I W I l I I 0 0.2 0.4 Contact Radius (mm) Figure 3.15. Contact force vs. contact radius, comparison I 1 0 0.01 0.02 0.03 Center Deflection (mm) Figure 3.16. Contact force vs. center deflection 80  / / / 70 _____________ 2  Vertical Force (No Slip) 0o Vertical Force (No Friction) 60  Sso 40 U O 40 a / 20 S40 / 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 Center Deflection (mm) Center Deflection (mm) Figure 3.17. Friction effect CHAPTER 4 DAMAGE ANALYSIS 4.1 Failure Criteria The calculated stresses and strains from the axisymmetric elements are based on the assumption that the plate is transversely isotropic. However, the plate consists of multiple plies where each ply is orthotropic and has different fiber orientation. Hence, direct application of any existing failure criteria which are based on stresses in the material coordinates is not possible. An estimate of stresses in the material coordinates has to be obtained from the axisymmetric results. The existing failure criteria for composite materials can be divided into two types. One is a interactive failure criterion and the other is a set of independent failure criteria for different failure modes. The interactive failure criterion is basically an empirical technique, and determination of the components in the polynomials of the criterion requires difficult and expensive biaxial tests. The most important criterion in this category is the TsaiWu theory, which uses a general quadratic polynomial including 79 linear terms. This theory says that a failure surface in the stress space exists in the form as Fjoj + Fj = 1 ij = 1,.. 6 (4.1) where Fi and F.. are second and fourth order strength tensors. This criterion can not be widely used because of difficulties in obtaining Fij. These criterion only predicts the onset of failure and not the failure modes. However, if Fi and Fj can be obtained correctly, this criterion may be satisfactory, and it is simple to use [13,14]. Independent failure criteria, typically maximum stress theory or maximum strain theory, can predict the onset of failure and also failure modes. Distinguishing failure modes is very important for the design of composite materials and also for the progressive finite element program which must modify its total stiffness matrix depending on the failure modes. However, it lacks accounting for the interaction of various components of stresses or strains in the failure process. Good failure criteria must predict the principal failure modes separately, such as fiber breakages, matrix cracking and delamination, and each failure mode criterion should be interactive such as the quadratic formula while they do not require biaxial tests to determine the coefficients. The quadratic failure criteria must take into account of the 80 compressive or tensile stresses. The work of Hashin, which deals with the threedimensional stresses fits this category and is given in Ref. [13]. It leads to six independent criteria, given in Equations (4.2) to (4.7). As input to the failure criteria, the following material strength parameters are needed. X, ; fiber direction tensile strength Xc ; fiber direction compressive strength YT ; transverse tensile strength Yc ; transverse compressive strength S12 ; shear strength in the 12 plane S1, ; shear strength in the 13 plane S23 ; shear strength in the 23 plane X, ; interlaminar (throughthethickness) tensile stress S1 ; interlaminar shear strength Fiber failure criteria have two modes; compressive fiber buckling mode and tensile fiber breakage mode. For tensile fiber direction stress, ao > 0, Hashin's fiber failure criterion is S + 12 + 13 1 (4.2) T S 12 S13 For the fiber buckling mode, oa < 0, the criterion is S 1 (4.3) There are also two types of failure modes for the matrix cracking, depending upon o2 + 03, such that when (02 + 03) > 0, the matrix cracking criterion is (2 + (Ti2 2 ( 1(434) 23 + 230203 ) + 12 =1 (4.4) YT S23S12 S13 and when (02 + 03) < 0, the matrix cracking criterion is 1 1 (o02 +) + 2( 1 (2+o 3)2 Yc[ 2S23 4S 23 (4.5) T23 02 03 12 3 )2 S23 S2 3 S13 There are two basic approaches to predict delamination initiation [32]. One is the mechanics of materials approach, which essentially compares the local state of stress in the interply matrix layer where delamination occurs with relevant strength parameters. The other approach is based on the application of fracture mechanics principles and will not be 82 discussed here. However, once the delamination initiates, the concept of a critical strain energy release rate, Gc, will be applied to predict delamination growth. In this study, delamination was predicted by the mechanics of materials with a simple quadratic interactive formula which depends upon 03, such that when 03 < 0, only shear stresses induce delamination such as Mode II type, and the criterion is 2 2 T13 +2 =2 1 (4.6) Si. when o3 > 0, it behaves like a combination of Mode I and Mode II type crack initiation, and the criterion is S)2 T13 + 23 =1 (4.7) x, ) SI Here XI and S, can be assumed as the tensile strength and pure shear strength of the matrix. This assumption is justified by the experimental observations that toughened matrix materials apparently reduce the onset of delaminations. Assume damage of a particular nature occurs at a location given by a in the 00 layer. Then the same type and extent of damage will occur in the 450 layer at (450 + a) and in the 450 83 layer at (450 + a) and so on. That means in a quasiisotropic laminate the damage pattern will look like a staircase in the form of a helix. This phenomenon has been observed in impact experiments also [15,33]. At a particular segment where the fiber direction is at angle 0 to the x axis, the stress components in material coordinates can be obtained as follows 01 cos2 sin28 2cos6sin8 0 0 0 o 02 sin28 cos26 2cos0sin8 0 0 0 T12 cosSsinO cosOsinO cos20sin2 0 0 0 on 03 0 0 0 1 0 0 Ozz T13 0 0 0 0 cos8 sin0 oa 3 0 0 0 0 sinO cosO oez (4.8) because oae and Oez are zero for the axisymmetric analysis, stresses can be easily obtained as in Equation (4.9). o0 = cos2 o r + sin2 8 Oe 02 = sin2 08 a + cos2 8 ge tz2 = cos 0 sin 8 (or Ge) (4.9) 03 = zz T13 = COS 0 orz Tz3 = sin 8 orz 84 By combining Equations (4.2) through (4.7) and Equation (4.9), the failure criteria for implementation into the finite element program are obtained as follows. Fiber failure (i) cos28 or + sin28 ao > 0 cos2 o0 + sin2eo) O' cos) sinO (oOr) 2 Cos0 rz2 XT S12 S13 ( ii ) cos2 or + sin2e a < 0 cos28 r + sin28e Xc Matrix failure (i) sin28 or + cos28 o + ozz > 0 1 (sin2 z sinr or zzCS2 0 o) S23 cos9 sin9 (a0 0a) 2 cos rz \2 S12 ( 13) + sin2 or + cos2 O + azz YT 1 (ii) sin28 ar + cos2e o + azz < 0 (sin2e rzsin OzzcOs2e 09 ozz) Scos9 sine (aeo,) cose arz S12 S 13 + 4 (sin28 O +cos28 ae+6zz)2 = 1 4 53' Delamination (i) ozz < 0 (a2 (ii) ozz > 0 orz)2 z)2 SI I = ) 1 Note that the delamination failure criteria are not functions of 8. This fact can also be verified from the Cscan results shown in Figure 1.5. 0 0 Ply One Element 0 45 Ply 45 Ply 45 Ply 0 Ply 90 Ply Figure 4.1 Partial failure patterns Table 4.1. Ultimate strength of AS4/35016 lamina unit (MPa) XT ; 2172 Xc ; 1724 YT ; 53.8 Yc ; 205.5 S12 ; 120.7 S 3 ; 120.7 S23 ; 89.3 X ; 56.5 S, ; 89.3 * Courtesy of Hercules Aerospace Company 88 4.2 Implementing Delamination into the Finite Element Program 4.2.1 Introduction Delamination initiates from where the shear stress is high, such as above the midsurface of the plate close to the contact surface. Then delaminations start to expand radially, and more delaminations appear in the bottom plies. Delaminations can be simulated in the finite element program by creating slip surfaces between elements as seen in Figure 4.2. New nodes are created, and rigid elements (spring elements with very large spring constants K) are added between two nodes across slip surfaces when the normal force is compressive to prevent interpenetration of the nodes. Adding new surfaces and consequently adding new nodes makes a cumbersome procedure, which is a weakness of the finite element method. The total stiffness matrix of the plate must be recalculated after delaminations occur or grow. Also the dimension of the total matrix is increased due to the introduction of new nodes. If tensile stress exists at the tip of delaminations, the spring element must be eliminated to allow Mode I type (Opening Mode) delamination propagation. Then, iteration of the program is necessary to assure that the rigid elements are placed only where compressive stresses exist. However, unless the contact force is high, delaminations seemed to appear only at the upper portion of the plate where throughthethickness stresses oz are compressive. The present finite element program monitors the 89 signs of all oa across the delamination surfaces while implementing Hashin's failure criteria, and it was found that all az's are compressive in the present study. It is safe to assume that in this study the delamination growth is Mode II type (Shear Mode) and to put rigid elements at all the nodes across delamination surfaces. When delaminations become extensive, especially beyond the midsurface of the plate, then the tip of the delaminations could have high tensile stresses. This Mode I crack growth combined with Mode II type will accelerate delamination growth faster than by Mode II delamination alone. 4.2.2 Numerical Experiment A simple numerical experiment was done to check the effect of delaminations on the contact force vs. center deflection curve and also to see the frictional effect across the delamination surface. Assume that the experimental test configuration shown in Figure 1.2 has preexisting delaminations of 12.7 mm diameter uniformly distributed through the thickness. As 16 elements were employed through the thickness, 15 delaminations were assumed in the finite element program. An incremental point force was applied at the center of the plate. Rigid elements were placed across all the delamination surface nodes as shown in Figure 4.3. Figure 4.4 shows contact force vs. center deflection curve of the plate with preexisting delaminations compared 