Quinn“. .3 !- mus»: mine-3M ‘_ 1" {fix «a , \ . 3 .u.!?... _ dirt .3 a ‘ 9! wawfifls... . “as; .6 v. 3.....- 1.. i!»- (11 ,. ; human? I I , id. fwma :2 1.. In? tow-5.}. tel... . to. .. 5 . $73.9“. 19‘ an? 4, . t . s . In.) (5....» ha... .032 :‘t t «.3»... dang... ‘34.. .h. ‘ fV ~7 x: to z {40... . 3:15 Rwufli3a...... 2. a.» . sin.) figufifi. hmhur...m.m.w.rh.§ue 3.. mi}. .31}... .3, 306 p t ‘3 . . alga-Ha: 1.3. if... rfimufiaiéx}. .3 flu. in «2&1! and! 3.... p.139. v . WO “W ”a 3CD" LIBRARY Michigan State University This is to certify that the thesis entitled ANALYSIS OF DAMAGE IN LAMINATED COMPOSITES SUBJECTED TO BLAST LOADING presented by ANOOP GOYAL has been accepted towards fulfillment of the requirements for the MS. degree in MECHANICAL ENGINEERING ' nature Major Professor's 5723/0 7 Date MSU is an affirmative-action, equal-opportunity employer -<—.--o--¢—- . -n-I-I-O-O-O‘I-O-O-I-0-0-l-O-I-i-O-I-n-nnn .- 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 6/07 p:/ClRC/DateDue.indd-p.1 ANALYSIS OF DAMAGE IN LAMINATED COMPOSITE PLATES SUBJECTED TO BLAST LOADING By Anoop Goyal A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Department of Mechanical Engineering 2007 ABSTRACT ANALYSIS OF DAMAGE IN LAMINATED COMPOSITES SUBJECTED To BLAST LOADING By Anoop Goyal A finite element code based on interlaminar shear stress continuity and interlaminar displacement discontinuity conditions and presented by Lee was utilized in the thesis research to analyze the in-plane damage and interfacial condition of a composite laminate subjected to blast loading. The displacement discontinuity conditions between different composite plies were enforced by using both shear slip approach and embedded layer approach. The delamination scheme was validated for a composite plate subjected to sinusoidal loading by an elasticity solution extended from the Pagano’s problem. For perfectly bonded composites, the finite element solutions matched well with the elasticity solutions given by Pagano. However, for imperfect bonding conditions, the finite element model could only account for the shear slip across the interfaces, but not the normal separation. To understand the nature of blast loading, a brief explanation about the pressure profile and other parameters involved in blast waves was elaborated. Numerical results obtained by simulating the blast loading on composite laminate were presented. Various in-plane and delamination failure criteria were applied using these results. Although the numerical results for interlaminar analysis agreed with that of experiment, the numerical results for in-plane analysis did not match with that of experiment. Absence of geometric non-linearity in the laminate theory and coupling between in—plane stretching and out-of-plane deflection were considered to be the sources of discrepancy. To my dear friends iii ACKNOWLEDGEMENTS First of all, I would like to thank my advisor Dr. Dahsin Liu, who gave his lot of time in guiding me. Although I did many mistakes, however, he went on patiently encouraging me. In other words, he was involved with me in my MS thesis. I would also like to thank Dr. Lee, who helped me Significantly in understanding the Finite Element formulation and giving valuable inputs. 1 would like to Sincerely thank Professor Krishnan who by his example, students and words gave me tremendous boost to carry on with my MS thesis in an enthusiastic way. Lastly I would like to thank Dr. Anshu Manik for his wonderful example and encouragement. iv TABLE OF CONTENTS List of Tables ....................................................................................... viii List of Figures ........................................................................................ ix 1. INTRODUCTION ................................................................................. 1 1.1 Background and Motivation ................................................................. 1 1.2 Organization of Thesis .............................................................................................. 2 2. LITERATURE REVIEW ............................................................................................... 3 2.1 Blast theories and Phenomenological models .......................................................... 3 2.2 Interfacial Damage .................................................................................................... 6 2.3 Blast Simulations ...................................................................................................... 9 3. ELASTICITY SOLUTION FOR BI-DIRECTIONAL COMPOSITES ...................... 13 3.1 Modified Pagano’s Problem ................................................................................... 13 3.1.1 Fundamental Equations .................................................................................... 14 3.1.2 Boundary Conditions ....................................................................................... 14 3.1.3 Interfacial Conditions ...................................................................................... 16 3.1.4 Elasticity Solution ............................................................................................ 17 3.1.5 Embedded layer Approach ............................................................................... 21 3.1.5.1 Determination of critical thickness value ................................................. 21 3.1.4.2 Alteration in properties of embedded layers ............................................. 25 3.1.6 Slip Approach .................................................................................................. 30 3.2 Relationship between the Embedded Layer and the Slip Approach ....................... 35 4. FTNITE ELEMENT SOLUTIONS ............................................................................... 39 4.1 Lee’s Laminate theory ............................................................................................ 39 4.1.1 Displacement Field .......................................................................................... 39 4.1.2 Fundamental Equations .................................................................................... 41 4.1.3 Boundary Conditions ....................................................................................... 43 4.1.4 Interfacial Conditions ...................................................................................... 43 4.1.4.1 Linear Shear Slip Theory ......................................................... 43 4.1.5 Generalized Strain Solutions ........................................................................... 45 4.1.4.2 Embedded Layer Approach ...................................................... 5 3 5. VALIDATION OF FORTRAN CODE ........................................................................ 56 5.1 Perfect Bonding Conditions .................................................................................... 57 5.2 Imperfect Bonding Conditions ................................................................................ 63 5.3 Validation by Pagano’s Slip Approach ................................................................... 67 6. BLAST SIMULATIONS .............................................................................................. 72 6. 1 Experimental Work [39] ........................................................................................ 72 6.2 Analysis of blast propagation ................................................................................. 76 6.2.1 Equations for air blast transmission ................................................................. 77 6.2.2 Pressure profile ................................................................................................ 81 6.2.2.1 Pressure-time profile ................................................................................. 81 6.2.2.2. Pressure-distance profile .......................................................................... 88 6.3 Finite Element Analysis .......................................................................................... 92 6.4 Results and Discussion ........................................................................................... 93 6.4.1 In-plane failure ................................................................................................. 93 6.4.2 Strain Analysis ........................................................................ 101 6.4.2 Delarnination Analysis ................................................................................... 105 7. CONCLUSIONS AND RECOMMENDATIONS ..................................................... 109 7.1 Conclusions ........................................................................................................... 109 7.2 Recommendations ................................................................................................. 110 Appendix A ......................................................................................... 11] Figures of stresses and displacements in order to determine critical thickness of the embedded layer Appendix B ......................................................................................... 114 Matlab code for the analysis of imperfectly bonded composite laminate using the Embedded Layer Approach Appendix C ......................................................................................... 127 Figures of stresses and displacements obtained using the Embedded Layer approach in the Pagano’s solutions Appendix D .......................................................................................... 1.30 Matlab code for the analysis of imperfectly bonded composite laminate using the Slip Approach Appendix E ......................................................................................... 142 Figures of stresses and displacements obtained using the linear shear slip theory and linear normal separation theory in the Pagano’s solutions Appendix F ......................................................................................... 145 Figures of stresses obtained by using the linear shear slip theory in ISSCT Appendix G ........................................................................................ 147 Figures of stresses and displacements obtained using the Embedded Layer approach in ISSCT Appendix H ........................................................................................ 149 Figures of stresses and displacements obtained by comparing ISSCT with Pagano for imperfect bonding condition vi Appendix I ......................................................................................... 151 Figures of stresses obtained by comparing ISSCT-Shear Slip only with Pagano for imperfect bonding conditions BIBLIOGRAPHY ............................................................................... l 52 vii Table 6.1 Table 6.2 Table 6.3 Table 6.4 Table 6.5 LIST OF TABLES Maximum in plane stress values for plyl and ply2 .......................... 98 Result from various failure criterions .......................................... 99 Maximum In-plane strains for plyl and ply 2 ................................ 105 Values of shear stresses across each interface at different instants of time ........................................................................... 104 Results from delamination criterion ........................................... 108 viii Fig. 2.1 Fig. 3.1 Fig. 3.2 Fig. 3.3 Fig. 3.4 Fig. 3.5 Fig. 3.6 Fig. 3.7 Fig. 3.8 Fig. 3.9 Fig. 3.10 Fig. 3.11 Fig 3.12 Fig. 3.13 Fig 3.14 LIST OF FIGURES Sources and Mechanism of Blast waves ..................................... 4 A [0/90/0] plate subjected to sinusoidal loading on the top surface ..................................................................... l3 Normalized shear stress along the height of a [0/0/90/0/0] composite plate for different ELT values ................................. 22 Normalized shear stress along the height of a [0/0/90/0/0] composite plate for different ELT values ................................. 23 Normalized displacement along the height of a [0/0/90/0/0] composite plate for different ELT values ................................. 23 Normalized displacement along the height of a [0/0/90/0/0] composite plate for different ELT values ................................. 24 Error calculation of the stresses and displacements between perfectly bonded [0/0/90/0/0] and [0/90/0] composite .................. 25 Normalized Transverse Shear stress distribution obtained by varying the material properties of the EL .................................. 27 Normalized Transverse shear stress distribution obtained by varying the material properties of the EL .................................. 27 Normalized in-plane displacement distribution obtained by varying the material properties of the EL .................................. 28 Normalized in-plane displacement distribution obtained by varying the material properties of the EL .................................. 29 Normalized shear stress along [0/90/0] composite plate obtained by different Shear slip constants and normal slip constant .............. 32 Normalized shear stress along [0/90/0] composite plate obtained by different shear slip constants and normal slip constant .............. 33 Normalized displacement along [0/90/0] composite plate obtained by different Shear slip constants and normal slip constant .............. 34 Normalized displacement along [0/90/0] composite plate obtained by different shear slip constants and normal slip constant .............. 34 ix Fig. 3.15 Fig. 3.16 Fig. 4.1 Fig. 4.2 Fig. 4.3 Fig. 4.4 Fig. 4.5 Fig. 4.6 Fig. 4.7 Fig. 4.8 Fig. 4.9 Fig. 4.10 Fig. 5.1 Fig. 5.2 Fig. 5.3 Fig. 5.4 Fig. 5.5 Pure Shear deformation in embedded layers ............................... 36 Pure normal deformation in embedded layers ............................. 37 A composite plate subjected to loading on the top surface .............. 39 Nodal Variables for the composite laminate ............................... 41 Normalized shear stress distributions along the height of [0/90/0] composite plate by different shear slip constants D” and Dsy ........ 50 Normalized shear stress distributions along the height of [0/90/0] composite plate by different shear slip constants D” and D”, ........ 51 Normalized displacement distributions along the height of [0/90/0] composite plate by different shear slip constants D” and 13,, . . . ......52 Normalized displacement distributions along the height of [0/90/0] composite plate by different shear slip constants D“ and 13,, ........ 52 Normalized shear stress distribution by varying the material properties of the embedded layers of a [0/1/90/1/0] composite plate...53 Normalized shear stress distribution by varying the material properties of the embedded layers of a [0/1/90/1/0] composite plate...54 Normalized displacement distribution by varying the material properties of the embedded layers of a [0/1/90/1/0] composite plate...55 Normalized displacement distribution by varying the material properties of the embedded layers of a [0/1/90/1/0] composite plate...55 Approaches to analyze interfacial damage in a composite .............. 56 Comparison of Normalized in-plane stress between Lee and Pagano of a perfectly bonded [0/90/0] composite plate ........................... 57 Comparison of Normalized in-plane stress between Lee and Pagano of a perfectly bonded [0/90/0] composite plate ........................... 58 Comparison of Normalized shear stress between Lee and Pagano of a perfectly bonded [0/90/0] composite plate ........................... 59 Comparison of Normalized shear stress between Lee and Pagano of a perfectly bonded [0/90/0] composite plate ........................... 59 Fig. 5.6 Fig. 5.7 Fig. 5.8 Fig. 5.9 Fig. 5.10 Fig. 5.11 Fig. 5.12 Fig. 5.13 Fig. 5.14 Fig. 5.15 Fig. 5.16 Fig. 5.17 Fig. 5.18 Fig. 5.19 Fig 6.1 Fig. 6.2 Comparison of Normalized displacement between Lee and Pagano of a perfectly bonded [0/90/0] composite plate .................. 60 Comparison of Normalized displacement between Lee and Pagano of a perfectly bonded [0/90/0] composite plate .................. 60 Comparison of normalized transverse stresses obtained by 36 elements vs. 100 elements ............................................ 62 Comparison of normalized transverse stresses obtained by 36 elements vs. 100 elements ............................................ 62 Comparison of normalized transverse shear stress along the height of a delaminated [0/90/0] plate between IS SCT and Pagano ................. 64 Comparison of normalized transverse shear stress along the height of a ’ delaminated [0/90/0] plate between ISSCT and Pagano ................. 64 Normalized deflection calculated at the center of the top layer of the [0/90/0] composite plate for different bonding conditions ............... 65 Comparison of normalized in-plane displacement along the height of a delaminated [0/90/0] plate between Lee and Pagano .................... 66 Comparison of normalized in-plane displacement along the height of a delaminated [0/90/0] plate between Lee and Pagano ..................... 66 Comparison of normalized transverse shear stress along the height of a [0/90/0] composite plate for 2 different bonding conditions ............. 69 Comparison of normalized transverse shear stress along the height of a [0/90/0] composite plate for 2 different bonding conditions ............. 69 Comparison of normalized in-plane deflection along the height of a [0/90/0] composite plate for 2 different bonding conditions ............. 70 Comparison of normalized in-plane deflection along the height of a [0/90/0] composite plate for 2 different bonding conditions ............. 70 Normalized deflection calculated at the center of the top layer of the [0/90/0] composite plate for different bonding conditions ............... 7] Shock wave generator at MSU ............................................... 73 A funnel-type wave transformer ............................................. 74 xi Fig. 6.3 Fig. 6.4 Fig. 6.5 Fig. 6.6 Fig. 6.7 Fig. 6.8 Fig. 6.9 Fig. 6.10 Fig. 6.11 Fig. 6.12 Fig. 6.13 Fig. 6.14 Fig. 6.15 Fig. 6.16 Fig. 6.17 Fig. 6.18 Fig. 6.19 Fig. 6.20 Location of pressure transducer ............................................... 75 Loading history for finite element simulation .............................. 75 Blast Wave ....................................................................... 76 Pressure Time History ......................................................... 81 Pressure time profiles given by different models ........................... 88 Schematic diagram of a laminate subjected to blast loading .............. 89 Finite element discretization of the [0/90/90/0] composite plate ......... 93 Variation of transverse normal stress at the interface across the plane of the quarter [0/90/90/0] composite plate ...................................... 94 Variation of transverse normal stress at the interface across the plane of the quarter [0/90/90/0] composite plate ...................................... 95 Variation of transverse shear stress at the interface across the plane of the quarter [0/90/90/0] composite plate .......................................... 95 Variation of transverse normal stress at the center across the height of the [0/90/90/0] composite plate .................................................... 96 Variation of transverse normal stress at the center across the height of the [0/90/90/0] composite plate .................................................... 97 Variation of transverse normal stress at the center across the height of the [0/90/90/0] composite plate .................................................... 97 Experimental Results .......................................................... 100 Variation of in-plane normal strain at the interface across the plane of the quarter [0/90/90/0] composite plate ................................. 101 Variation of in-plane shear strain at the interface across the plane of the quarter [0/90/90/0] composite plate ................................. 102 Variation of in-plane normal strain at the center across the height of the [0/90/90/0] composite plate .......................................... 103 Variation of in-plane normal strain at the center across the height of the [0/90/90/0] composite plate ............................................ 103 xii Fig 6.21 Fig. 6.22 Fig 6.23 Fig. A] Fig. A.2 Fig. A.3 Fig. A.4 Fig. A.5 Fig. C] Fig. C.2 Fig. C.3 Fig. C.4 Fig. 0.5 Fig. E] Fig. E.2 Fig. E.3 Variation of in-plane shear strain across the height of the [0/90/90/0] composite plate ................................................. ..104 Variation of transverse normal stress at the interface across the plane of the quarter [0/90/90/0] composite plate .................................... .106 Variation of transverse normal stress at the interface across the plane of the quarter [0/ 90/90/0] composite plate ................................... Normalized in-plane normal stress distribution for different values of BLT with perfect bonding conditions Normalized in-plane normal stress distribution for different values of BLT with perfect bonding conditions Normalized transverse displacement for different values of BLT with perfect bonding conditions .......................... Normalized transverse normal stress distribution along the height for different values of ELT with perfect bonding conditions .......... Normalized in-plane Shear stress distribution along the height for different values of BLT with perfect bonding conditions. . . . . . . . . . Normalized transverse displacement along the height obtained by varying EL properties of a [0/1/90/1/0] composite plate [I-isotropic]. .. Normalized in-plane normal stress along the height obtained by varying EL properties of a [0/1/90/1/0] composite plate [I-isotropic] . .. Normalized in-plane normal stress along the height obtained by varying EL properties of a [0/ I/90/I/0] composite plate [I-isotropic] . .. Normalized transverse normal stress along the height obtained by varying EL properties of a [0/1/90/1/0] composite plate [I-isotropic]. .. Normalized in-plane Shear stress along the height obtained by varying EL properties of a [0/1/90/1/0] composite plate [I-isotropic]. .. Normalized in-plane normal stress distribution along the height by varying shear slip constants D” and Dsy and normal slip constant k. .. Normalized in-plane normal stress distribution along the height by varying shear slip constants D” and D” and normal slip constant k.. Normalized transverse normal stress distribution along the height by xiii .. 107 111 111 ..112 ..112 113 127 127 128 128 129 142 .142 Fig. E.4 Fig. E.5 Fig. F.l Fig. F.2 Fig. F.3 Fig. G.1 Fig. G.2 Fig G.3 Fig. H.l Fig. H.2 Fig. H.3 Fig. 1.1 Fig. 1.2 varying shear slip constants D“ and D3), and normal slip constant k. .. Normalized in-plane shear stress distribution along the height by varying shear slip constants D“ and D”, and normal Slip constant k. .. Normalized transverse displacement distribution along the height by varying shear Slip constants D” and D3}. and normal slip constant k. .. Normalized in-plane normal stress distribution along the height obtained by varying the Shear slip constants D” and By . . . . . . Normalized in—plane normal stress distribution along the height obtained by varying the shear Slip constants D” and D3,, ............... Normalized in-plane shearstress distribution along the height obtained by varying the shear Slip constants Der and 13,, ............... Normalized in-plane normal stress along the height obtained by varying EL properties of a [0/1/90/1/0] composite plate [I-isotropic]... Normalized in-plane normal stress along the height obtained by varying EL properties of a [0/1/90/1/0] composite plate [I-isotropic].. Normalized in-plane shear stress along the height obtained by varying EL properties of a [0/1/90/1/0] composite plate [I-isotropic]... Comparison of normalized in-plane normal stress along the height of a delaminated [0/90/0] plate between Lee and Pagano [sz=Dsy=0.1, k=0.0417, E=0.024 psi,G=0.01 psi]. . . . . . . . Comparison of normalized in-plane normal stress along the height of a delaminated [0/90/0] plate between Lee and Pagano [sz=Dsy=0.1, k=0.0417, E=0.024 psi,G=0.01 psi] ........... Comparison of normalized in-plane shear stress along the height of a delaminated [0/90/0] plate between Lee and Pagano [sz=Dsy=0.l, k=0.0417, E=0.024 psi,G=0.01 psi]... . . . . Comparison of normalized in plane normal stress along the height of a [0/90/0] composite plate for two different bonding conditions. Comparison of normalized in plane normal stress along the height of a [0/90/0] composite plate for two different bonding conditions... xiv .143 .143 .144 145 .145 .146 .147 ..147 .148 149 .149 150 151 151 1. INTRODUCTION 1.1 Background and Motivation Explosions can be both awesome and devastating. They have their utilities in mining, building demolition, pyrotechnics and even construction. At the same time, they are of interest to those who are involved with disasters such as explosion and earthquakes, civil defense precautions and to those concerned with defense against terrorists [2]. Explosives when detonated result in powerful blast waves which can destroy the whole structure and the property and personnel along with it. Hence it is very important to understand the mechanism behind these blast waves, so that we can design our structures properly. Among all kinds of materials, composite are one of the most widely used. They find their usage in a wide variety of engineering applications ranging from aircrafts and submarines to pressure vessels and automotive parts. This is because of their high strength-to-weight and high stiffness-to-weight ratios, which make them light but at the same time very strong. However at the same time, these materials have their own limitations. Hence, when various types of loading come upon them during manufacturing or service conditions, these materials may result in various complex damage mechanisms. Delarnination is one of the most common types of damage in composites. When composites are subjected to loading conditions, their anisotropy and lack of homogeneity may cause mismatch of properties across different laminas, which in turn may lead them to separate from each other, or in other words delamination. This makes the whole structure very weak and easily penetrable. A finite element model which can capture the mechanism of delamination in composite materials and can simulate it for any given requirement is very significant. This will not only be inexpensive and saves time but also provide a lot of scope for parametric studies leading to a better understanding of the damage process. This thesis work is an attempt to study the mechanism behind blast loading and then incorporating it with in a suitable finite element model to analyze interfacial damage of laminated composites. 1.2 Organization of Thesis This thesis is divided into seven chapters. In Chapter Two, literature survey concerning the various aspects of blast loading, laminate theories dealing with interfacial damage and blast simulations are discussed. Chapter Three elaborates on the elasticity analysis for a simply supported [0/90/0] plate with different interfacial bonding conditions. Two analytical models — (1) Embedded layer approach and (2) Slip approach -— are presented and used for the analysis. Chapter Four gives us an explanation about Lee’s laminate theory. Using the numerical code based on Lee’s laminate theory [1], embedded layer approach and shear slip approach are individually utilized to analyze simply supported [0/90/0] plate with different interfacial bonding conditions. In Chapter Five, the results from finite element simulations in Chapter Four are compared with the elasticity solutions found in Chapter Three in details. Finite element simulations of composite laminates with a stacking sequence of [0/90/90/0] subjected to blast loading is presented in Chapter Six. The results from the simulations are compared with the experimental results. Chapter Seven summarizes the conclusions of the present study and ends with recommendations for future studies. The Appendices follows Chapter Seven and includes various figures and the computer programs developed for this thesis research. 2. LITERATURE REVIEW This chapter gives a detailed account of literatures surveyed for blast theories, various laminate theories dealing with interfacial damage and blast simulations. 2.1 Blast theories and Phenomenological models A blast wave is a pressure wave of finite amplitude generated due to a rapid release of energy in the medium [3]. The study of blast waves generally involves mechanisms and sources for its generation, analysis of pressure profile which includes variation of pressure with respect to time and distance and structural encounter of blast waves. There are many possible mechanisms by which blast waves can get generated. The two main mechanisms are detonation and deflagration. Depending on the amount of energy provided and the nature of the explosive, the explosive may undergo decomposition or oxidation by either deflagration or detonation, differentiated by the speed of the decomposition of the explosive. In the case of deflagration, the speed of decomposition of the explosive is way below the speed of sound, for instance, a typical gun propellant in a gun barrel burn at a rate of 400 mm/s. On the other hand, detonation is the form of the reaction of an explosive where the decomposition rate of the explosive goes in the range of 1500 — 9000 m/s, way above the speed of the sound. This reaction is accompanied by large pressure and temperature gradients at the shock wave front and the reaction is initiated instantaneously. One example of detonation is TNT which is having a detonation velocity of 6,940 m/s [4] There could be different sources for explosion. Depending on the medium and the surroundings, they will be producing different blast effects. The figure below gives a description about that. The nuclear and chemical explosives undergo the detonation mechanism to generate blast waves. Dust explosions, on the other hand are accidental Deflagration [ l l l , lNuclear “Chemical I Gas and Vapor L— Fuel Dust Cloud Gas Explosion ]Surfacel [Underwater ] lUnderground I Figure 2.1 Sources and Mechanism of Blast waves deflagrations which occur in places like flour mills, grain silos etc. due to high levels of suspending combustible dust particles in the air space during various operation [5]. Gas and vapor cloud explosives undergoes both detonation and deflagration depending on the initial conditions and the external environmental factors. For instance, vapor cloud explosive like hydrocarbon when blown in a confined environment with certain boundary conditions can undergo detonation and can lead to intense damaging consequences. On the other hand, the same explosive when blown in an unconfined environment like open air may not produce any damaging effect at all [6]. This thesis research will mainly be focusing on chemical explosives. The pressure waves generated afier the explosion vary as time and distance change and hence are a function of both. Over the period of many years, researchers tried formulating various models for different blast parameters numerically as well as by using available extensive data obtained by testing of explosives [7]. Many authors tried solving the equations for air blast transmission coupled with jump conditions, known as the Rankine—Hugoniot conditions [7] both numerically and analytically and proposed different solutions. Brode [8] was the first one to propose a solution in 1955. Henrych also proposed a solution [4] for the variation of the peak static overpressure with respect to distance from the source of explosion. It was found that the accuracy of predictions and measurements by these models in the near field was somewhat lower than in the medium to far field. This discrepancy was attributed to complexity of the flow processes involved in forming the blast wave close to the charge where the influence of explosive gas was difficult to quantify. Beshara [9] gave a detailed summary of many available models till 1992 for chemical explosions, nuclear explosions as well as unconfined vapor cloud explosions. This paper mainly dealt with external blast loading on aboveground structures and included nmnerical and empirical models for blast parameters like overpressure, time duration, arrival time, reflected and dynamic overpressure. In another paper, Beshara [10] gave description of different models concerning confined blast loading, explosion induced ground shocks and dust explosions. Graham and Kinney in their book [11] gave empirical models for various blast parameters like peak static overpressure, arrival time, duration of blast wave and intensity. Saleh and Adeli also cited another model [12] for overpressure, arrival time and pressure variation w.r.t distance. As the blast waves travel, the pressure profile also varies with respect to the time. When the blast waves reaches at a particular location, the pressure at that point suddenly jumps from atmospheric pressure to the pressure at the shock fiont and as time passes, the pressure at that location decays to atmospheric value, then drops to a partial vacuum and eventually returns to ambient pressure [13]. In order to explain this type of variation of pressure with respect to time, various authors proposed different phenomenological models. Flynn [7] in considering blast loading of structures assumed a linear decay of pressure and presented a two parameter model. This form, admittedly oversimplified, is ofien adequate for response calculations. Ethridge [7] gave an exponential two parameter model, which was more accurate than Flynn. Friedlander [7] gave a three parameter model which is most widely used as it is optimally accurate and also computationally inexpensive. Ethridge [14] presented a four parameter model which were more accurate but computationally more rigorous also. Brode [8, 15] finally gave two most complex five parameters models. All of these models will be explained in details in chapter six. 2.2 Interfacial Damage Although composite laminates find a wide variety of usages in the engineering structures, interfacial damage may take place in them during manufacturing operations or in service. This may reduce the load bearing capacity of the laminate. There are two types of interfacial damage which are observed in composite laminates. The first one is called the weak bonding and the second one is the delamination [16]. In the case of weak bonding there are displacement discontinuities, or displacement jumps between the layers despite non-zero continuous interfacial stresses. In the case of delamination, interfacial stresses vanish on those layers. It may be possible that the weak bonding may lead to delamination. Lot of literature is available on delamination. Two types of delamination namely edge delamination and central delamination [17] have been widely investigated. Both of them can be viewed as a result of interlaminar stress concentration caused by material property variation in the thickness direction and high external loading. The various analytical and numerical methods available for analyzing the delamination in composite plates can be broadly classified as region approach and layerwise models [18]. In a region approach, the delaminated laminate is divided into sublaminates or segments and the continuity conditions are imposed at the delamination junctions. Each of these sublaminates is analyzed using what is known as the equivalent single layer theories (ESL). The ESL plate theories are derived from the three dimensional elasticity theory by making suitable assumptions regarding the stress state or kinematics of deformation, which allow the reduction of a three dimensional problem into two dimensional problem. The ESL theories are firrther classified into classical laminated plate theory; first order shear deformation theory and third order shear deformation theory. A fuller perspective can be found in Reddy [19], Ghughal and Shimpi [20]. Roy and Chakraborty [21] developed a finite element model to study the chances of delamination at the interfaces of hybrid F RP composite laminates using three dimensional eight noded isoparametric solid elements. The delamination initiation at the interface was assessed using the criteria of Choi et.al [22]. In a layerwise model, modeling of the laminate is done using the layerwise theories which are based on piecewise layer-by-layer approximations of the response quantities in the thickness direction. This is further divided into full layerwise theories and the partial layerwise theories. The delamination can be modeled as an embedded layer or by introducing discontinuity functions in the displacement fields [18]. Detailed description can be found in Reddy [19], Ghughal and Shimpi [20] and Carrera [21]. Ghosal et. al. [24] studied the transient analysis of delaminated composite and smart composite plates using the improved layerwise theory or the partial layerwise theory, extended to include large deformation and the interlaminar contact in the delaminated interface. They used the first- order shear deformation-based displacement field to address the overall response of the laminate and layerwise functions to accommodate the complexity of zigzag-like in-plane deformation through the laminate thickness. In order to model the multiple, discrete delamination, the assumed displacement field was supplemented with Heaviside unit step functions, which allows discontinuity in the displacement field. Ghosal et a1 [25] used the F ermi-Dirac distribution fimction in place of the Heaviside unit step fimction to model a smooth transition in the displacement and the strain fields of the delaminated interfaces. The improved layer wise theory was incorporated into Fermi- Dirac distribution function to account transverse shear affects of anisotropic laminated composites. As far as the weak bonding is concerned, its study is again divided into two parts - shear slip and general weak bonding [16]. Shear slip implies no opening up of layers and accounts for shearing mode weak bonding. Liu et a1 [26] analyzed shear slip between the layers by using the interlayer shear slip theory (ISST). ISST implies that shear stresses are continuous across the layers; however the variation of interlaminar transverse displacements is linearly proportional to the interlaminar shear stress. The numerical version of the ISST theory is used in this thesis research. General weak bonding accounts for opening up of layers and hence opening mode weak bonding. Liu et. a1. [27] presented an interlaminar bonding theory, which was an updated version of ISST, in order to also account for general weak bonding. Interlaminar bonding theory, in addition to ISST, implies that interlaminar normal stresses are continuous and the interlaminar normal separation is linearly proportional to the interlaminar normal stresses. ISST and its updated version interlaminar bonding theory make use of layer wise displacement fields. Soldatos et al studied shear slip [28] and updated it to include the general weak bonding [29] using the global or smeared displacement theories rather then the layerwise theories. However, in all of the above studies, weak bonding along the whole interface was considered. Thus in order to account for local weak bonding which occurs on part of an interface, the above study [28, 29] were further extended by Shu [16] to develop a generalized laminate model which accounted for both local weak bonding and local delamination. The method used [28, 29, 16] is based on the successful combination of the stress analysis information obtained from the three-dimensional elasticity solutions of simply supported plates [30] with smeared or global laminate theories. The main advantage with this method was the use of small number of degrees of freedom, irrespective of the number of layers involved, whereas it can accommodate the boundary conditions imposed at the edges. 2.3 Blast Simulations The design of structures subject to blast loading has been a rapidly growing area of interest in the last few years. Lot of work has been done for explosive loading of materials such as concrete [31] and steel [32]. There has been very limited data available on polymer composite structure response to blast loading. Librescu and Nosier [33] investigated the response of laminated rectangular flat panels subjected to explosive blasts and sonic boom loadings. Only the theoretical analysis was being presented. F riedlander’s equation [7] was used to define the pressure time history. It was assumed that the plate dimensions are too small as compared to blast and sonic boom wave front and hence pressure at any instant is uniformly distributed over the plate. The response of the composite plate was characterized in terms of different parameters being found out using high order shear deformation theory and the results were being compared with their counterparts obtained by first order shear deformation theory and classical plate theory. It was found that working in the fiamework of higher order shear deformation theory was much better than the first order shear deformation theory and classical laminate theory. This was because the first order shear deformation theory required incorporation of a shear correction factor, largely dependent on the laminate sequence, relative anisotropy of the layers etc. and the claasical laminate theory gives inaccurate results for thick plates. Turkrnen and Mecitoglu [34] did the experimental, analytical and numerical study of the laminated composite plates subjected to blast loading. The pressure variation with respect to time was defined by using the Friedlander’s decay function while the spatial distribution was approximated by multiplying the Friedlander’s decay function with sine functions. The finite element modeling was done in ANSYS. The finite element model of the plate comprised of an assembly of two-dimensional eight noded shell elements with seven layers in the transverse direction and no interfacial conditions were assumed between the layers. On the theoretical side, dynamic equation of plates was solved by considering a new displacement function. The experiments were carried out on the laminated composite 10 plates with clamped edges for two different blast loads. Theoretical and finite element results were found to be in good agreement with each other. Qualitative agreement was found between the analyses and experimental in the first load for the first loading, however the results did not matched for the second one. Turkmen and Mecitoglu [35] did experimental and numerical investigation of a stiffened composite laminate plate subjected to blast loading. The pressure variation with respect to the time produced by the blast loading was defined using the Friedlander decay function. The laminated composite plate and the stiffner were modeled using the two dimensional, eight-noded, shell elements named SHELL91 in ANSYS. The element had six degrees of freedom. It was assumed that there was no slippage between the element layers and normals to the center plane before deformation remain straight after deformation. Good agreement was found between the peak strains given by the experimental and numerical analysis in both linear and non-linear ranges. Discrepancy was dound on the stiffner, however it was attributed to the presence of the adhesive layer between the plate and the stiffner. Zyskowski et.al. [36] used the numerical code AUTODYNE in order to simulate the structural response following the detonation caused by an explosion cloud in an unvented structure. Initial phases of the blast wave propagation were calculated one dimensionally using quadrilateral elements while the oblique wave reflections from the structure were modeled three dimensionally using brick-shaped elements. It was found that the numerical results calculated by AUTODYNE and the experimental results had a good correlation. Hence it was concluded that structural response caused due to internal blast loading can be properly addressed by the AUTODYNE code. Trabia et.al. [37] investigated different methods using LS-DYNA to simulate both the explosive loading 11 and the structural response of a cylinder subjected to it. One of the methods used for modeling the explosive loading involved using the Lagrangian elements along with the CONWEP air blast function Another two methods used the Arbitrary-Lagrangian- Eulerian (ALE) coupling procedure along with the Jones-Wilkins-Lee-Baker (JWLB) equations of state and Jones-Wilkins-Lee (JWL) equation of states for the explosive detonation products. The cylinder geometry was modeled using Lagrangian and Eulerian solid elements for composite layers. The numerical values for strain matched well with that of the experimental for two out of three cylindrical test cases. 12 3. ELASTICITY SOLUTION FOR BI-DIRECTION AL COMPOSITES 3.1 Modified 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 approach. Among the very few elasticity problems available for laminated composite analysis, Pagano’s problem [30] offers such an opportunity. Pagano had found out a three dimensional elasticity solution for stress and displacement fields of a composite laminate with simply supported boundary conditions subjected to sinusoidal loading on the top surface, as shown in Fig. 3.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 the interfacial damage between composite laminae, a modified Pagano’s problem was necessary [3 8]. Figure3.l. - A [0/90/0] plate subjected to sinusoidal loading on the top surface 13 3.1.1 Fundamental Equations Consider a laminate composed of N orthotropic layers. Figure 3.1 shows such a composite laminate with N=3. The body is simply supported and a normal traction 0', = q0 (x, y) is applied on the top surface. The constitutive equations for any orthotropic layer can be expressed by: 0x C11 C12 C13 8x 0y = C2, C22 C23 8), (3.1) 0'2 C31 C32 C33 ‘9 T}: = 44yyz’ sz : C55712’ Try = Cbbyxjr (32) where Cij are the compliance coefficients. While the governing field equations can be written in terms of the displacement components as, Cllu’xx +C66u’}:l' +C55u’zz +(C‘IZ + C66 )v’xj' +(C13 + C55 )w’n = 0 (C12 + C66 )um +C66v,“ +C,,v,_,,. +C,,v,,, +(C23 + C4, )w = 0 (3.3) ’ yz +C33w =0 ’22 (C13 + C55 )u’xz +(C23 + C44 )V’yz +C55 w,” +C44W ’ y)“ 3.1.2 Boundary Conditions The following boundary conditions are required if the exact solution to the modified Pagano’s problem is sought. They are identical to those used in the Pagano’s problem. A. The conditions on the top surface of the plate, 2 = + h/ 2 are given by 14 0. (Ly-’25) =qo(x, y) = min pr sin qy p, (x, y,+ h/2) = 0 (3.4) T}: (x, y,+ h/ 2) 0 where h is the total thickness, 0 is the length of the composite plates, b is the width of the . . 7t 7t . . composrte plates, qo rs a constant, p = — and q = Z . Theses three boundary condrtrons a reveal that the normal stress on the top surface of the composite beams is equal to the prescribed sinusoidal loading while the shear stresses vanish. B. The bottom surface is traction free, hence both the normal stress and the shear stress vanish, i.e 0'2[x,y,— 121-] = O Tn(x,y,- ’2/2) 0 (3.5) T}: (x9 y9_ h/Z) : 0 C. The following simply supported boundary conditions are imposed at both edges of the composite plates, i.e. Atx=0,a: 0' =v=w=0 (3.6) Aty=0,b: 0'},=u=w=0 Here u, v and w are the displacements in the x, y and 2 directions, respectively. 15 3.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 plates. However, in the modified 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 i and 1’ +1 denoting the adjacent layers corresponding to the laminate interface of interest, respectively. azi(x,y,—h':/2)= azi+l(x,y, Izzy/2) szi(x,ya—hi/2) = szi+l (x,y, hiH/z) (3.7) ryzi(xaya‘hi/2)= Tyzi+1(x’y’ hM/z) i=l,N—-l B. For the case of perfectly bonded composite laminate, the displacement will also be continuous across the interface of two consecutive layers and hence we have 16 ui(x,y,— hi/Z) ui+l (x, y, hH/Z) vi(x,y,— hi/Z) vi+1 (x, y, hi+1/2) (3.8) wi(x,y,- hi/2)= wi+l(x,y,hi+l/2) i=l,N—l However, 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. There are two ways to handle this situation -— the first approach is called the embedded layer approach which will be discussed in section 3.1.5 and the second approach make use of the normal separation theory and shear slip theory which will be explained in detail in section 3.1.6. We will call it the Slip Approach. Both of these approaches were used to get the elasticity solutions for a [0/90/0] composite plate subjected to sinusoidal loading and correlation between them was derived in section 3.1.7 3.1.4 Elasticity Solution The solution proposed by Pagano for perfectly bonded composite plates can be used to solve the modified Pagano’s problem. The following displacement equations are assumed for each single layer, which satisfy the simply support edge conditions (3.6), 17 u = U (z)cos px sin qy v = V(z)sin px cos qy (3.9) w = W(z)sin px sin qy A. If the material of the layer is non-isotropic then equations (3.9) and the governing field equations (3.3) on simplification yield the following stress equations, 3 0-,. = sin pxsin quMijUJ-(z), (i = 1,2,3) j=1 . 3 ryz = C44 srn pxcos qu(mJ-Lj + qu )Vj(z) Fl 3 sz =C55 COSPXSinqul(mj+PRj)I/j(z) (3.10) J: 3 rxy = C66 cospxcosqyzl(q + pLj Ill-(2) J: where 0'l , 0'2 , and 0'3 stands for 0'x , a", , and oz While the displacement equations are given by, 3 U (2 )= ZlUj (Z) J: V(z)=ZLjUj(z) (3.11) 18 U j(z)= chj(z)+ 0152(2) (j=l,2, 3) (3.12) W]. (z) = G jC j(z)+ a 1.1332(2) Hence, the solution process for N layers composite plate is reduced to determine 6N constants ij, ij (where an additional subscript k is introduced to identify the various layers), 6 constants for each layer, associated with the boundary conditions given by Equations (3.4) and (3.5) and the modified interfacial conditions given by either the embedded layer approach or the slip approach. Once these coefficients are identified, the stress and displacement distributions of the composite beam can be found from corresponding equations, i.e. Equations (3.10) and (3.9). B. In the event that the material of a given layer is isotropic, the stress equations (3.10) for that particular layer must be replaced by the following equations, 0x = l— pCHVl +C,2(-qV2 + V3')Jsin pxsin qy 0y = l—anV. +C.2(—pV. + fillsinpxsinqy 0'2 =[C11V3."C12(PV1 +qV2)]Sin szm qy (3.13) 2ryz = (CH —C12XV2' +qV3)sin pxcos qy 27x2 = (CH - C12 XVI + PV3 )COS PX Sin ‘0’ 22'er =(Cn _C22Xqu +pV2)cospxcosqy l9 Where Vi refer to the displacements and primes denote derivatives with respect to z. The displacement function given by Equations (3.11) and (3.12) for an isotropic layer must be replaced by the following equations (3.14) and (3.15), V.(z)=(a1i +a3iz+a5izzjexp.(cz)+(a2i +a4iiz+a6z 2)exp .(—cz) (i=1, 2, 3) (3.14) where V1 , V2, V3 are used in place of U, V, W and aj, are constants related as, =a .=0 a 61 51' “31 = Pa32’qa41 = ”42’ “’31 = pa33’ca41 = 'Pa43 (3.15) qa +6., ..._,,a +£Cl_2if_11 22 23 21 p C12+C11a41 C 3C W12 ”“13 =P"11 p[ (.3122 +C1111]“ “31 Hence the 6 independent constants for an isotropic layer to be determined using the interfacial conditions and boundary conditions are an, an, a3], a.“ and one term from each of the pairs (an, a13) and (an, a23). If the laminate contains one or more isotropic layers, we must replace equations (3.11), (3.12) and (3.10) with equations (3.14), (3.15) and (3.13) for those layers. In this way, we can handle any situation involving a laminate having a combination of isotropic and non-isotropic layers. 20 3.1.5 Embedded layer Approach One way to study the imperfectly bonded interface is by inserting extremely thin layers between the adjacent layers of the composite plate [26]. These embedded layers have a critical thickness value and their properties are isotropic in nature. In order to simulate all the possible bonding conditions from perfect to imperfect to delamination, the material properties (For example the Young’s Modulus, Shear Modulus etc.) of these interfacial layers are varied. Perfect bonding conditions can be simulated by selecting the properties of these embedded layers to be very close to matrix properties. Delamination can be simulated by selecting these properties to be near zero value and any bonding condition in between the perfect bonding and delamination can be simulated by varying the properties of these embedded layers. In conclusion, two steps are involved in embedded layer approach. Firstly, the critical thickness value of the embedded layer is required and secondly, the properties of the embedded layers are altered in order to simulate different bonding conditions. 3.1.5.1 Determination of critical thickness value In order to determine this, we use Pagano’s solution assuming perfect bonding conditions (which was mentioned in section 3.1.4). We will take a [0/90/0] composite as given in Figure 3.1 and introduce two embedded layers such that the modified composite is [0/0/90/0/0], assuming that the properties of these embedded layers is same as that of top layer with orientation 0°. This was done for the purpose of making the analysis easier. It is observed that as the thickness of the embedded layer reduces and approaches a critical value, the three-dimensional elasticity solutions for the composite with embedded layers, 21 i.e. [0/0/90/0/0] approaches the composite without interfacial layers, i.e. [0/90/0] for the case of perfect bonding conditions between the interfacial layers. Following plots clearly show this comparison between various stress and displacement values. 0.2 — -+— ELT-0.1 . 9‘1”" ~ +ELT.o.o7 ....-- —a—ELr-o.05 3,4,9; +ELT-0.03 1‘1?! 0 P —-v— ELT-0.01 19:91:; _ —+—ELT-o.oo1 1!... ~ -+—-ELT-0.0001 183;, +no EL .3“. 0 0.05 0.1 0.15 0 2 0.25 0.3 0.35 Figure3.2- Normalized Transverse shear stress distribution along the height of a [0/0/90/0/0] composite plate for different values of embedded layer thickness (ELT) with perfect bonding conditions 22 0.2 — ’33 O F -~— ELT-0.1 -B- ELT-0.07 —e— ELT-0.05 —€r— ELT-0.03 + ELT-0.01 + ELT-0.001 -—'-- ELT-0.0001 —6-no EL l -0.8 l l l 0 0.05 0.1 0.15 r:(a/ 2 ,0, z) 0.2 0.25 Figure3.3- - Normalized Transverse shear stress distribution along the height of a [0/0/90/0/0] composite plate for different values of embedded layer thickness (ELT) with perfect bonding conditions 0‘8 I I I I I I —¢— ELT-0.1 —B— ELT-0.07 —a— ELT-0.05 —0— ELT-0.03 + ELT-0.01 + ELT-0.001 ~+- ELT-0.0001 +110 EL 1 0.002 _ a U(0,b/2,z) _08 1 1 1 1 - .01 -0.008 -0.006 -0.004 -0.002 Figure3.4- - Normalized in-plane displacement variation along the height of a [0/0/90/0/0] composite 1 1 0.006 0.008 0.01 plate for different values of embedded layer thickness (ELT) with perfect bonding conditions 23 0.6 4 —*— ELT-0.1 ~B- ELTo0.07 —a— ELT-0.05 —e— ELT-0.03 O" + ELT-O .01 H —*— ELT-0.001 ~—+— ELT-0.0001 +110 EL 0.2 — — ° - a N -0.2 - 4 -O.4 - — '0.8 ’- 1 -0 8 1 4 1 1 1 1 1 L 1 -0.025 -0.02 -0.01 5 -0.01 -0.005 0 0.005 0.01 0.01 5 0.02 0.025 Via/2 .0. 21 Figure3.5— - Normalized in-plane displacement variation along the height of a [0/0/90/0/0] composite plate for different values of embedded layer thickness (ELT) with perfect bonding conditions Figures (3.2) and (3.3) show the variation of ‘L'yz and In with respect to the height of the laminate. From these figures it is very clear that as thickness of the embedded layers approaches towards zero, the solution of the [0/0/90/0/0] approaches the solution of [0/90/0]. Figures (3.4) and (3.5) display the variation of displacement across the thickness of the composite for different thickness of the embedded layers. Appendix A gives a list of all the remaining figures for the stresses and displacements. Fig 3.6 below gives a description of the percent error in the stresses and displacements between the [0/0/90/0/0] composite with different ELT (embedded layer thickness) and the [0/90/0] composite with no embedded layer. We can easily see that, when the thickness of the embedded layer is either 0.0001 inch or 0.001 inch, which are 0.011% and 0.11% of the total thickness, the error is very small. Moreover, for all practical purposes the thickness of the matrix joining the adjacent composite layers is generally of 24 the order of 0.001 in and hence we will select the critical thickness of the embedded layer to be equal to 0.001 in [26]. 100 50 - f _/ 0 . _—:—"— #1 A L 1 ._\ .50 \ 15 +TauYZ t m -100 — -I— Sigmax :3 + SlgmaY .150 . *Sigmaz '200 ‘ +TauXZ _250 j +TauXY —+—— Dlst .300 - -350 32 :2 a~° :2 a! :3 ..\° (0 ('3 2 O N ('3 n o n _ o on Q n O O M M N M ' n c“ «5 Thickness maqub‘added layefitop layer) Figure3.6- Error calculation of the stresses and displacements between perfectly bonded [0/0/90/0/0] composite plate and perfectly bonded [0/90/0] composite Thus we find from the above plots that the critical thickness value of the embedded layer is 0.001 in, which is 0.11% of the total thickness of the composite laminate. At this particular thickness, the stress and the displacement values are matching with the least amount of error. For instance, the maximum error for oz is 1.87%, while for 1,2 and In it comes out to be 1.15% and 1.34%. 3. 1.4.2 Alteration in properties of embedded layers After determining the critical thickness value, we make the embedded layers to be isotropic. The thickness of these isotropic embedded layers is 0.001 inch, or 0.1 1% of the 25 total thickness of the composite laminate. In order to simulate all the bonding conditions from perfect to imperfect the shear modulus and young’s modulus of these layers are varied from the order of 106 psi to 0 psi. The solutions are obtained by determining the unknown coefficients using the boundary conditions given by equations (3.4), (3.5) and (3.6) and the interfacial conditions given by equations (3.8) and (3.9). Once these unknown coefficients are evaluated, the stress and displacement distribution of the composite plate can be found from corresponding Pagano’s three dimensional elasticity solutions (given before in section 3.1.4). For the isotropic layers in the composite plate, the elasticity solutions are found out by using equation (3.13), (3.14) and (3.15). While, for the non-isotropic layers, the elasticity solutions are derived by using the equations (3.10), (3.12) and (3.13) 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 B. Based on the elasticity solutions, following plots for stress and displacements were obtained 26 0 02 04 °"’ ace/312) —+— E=2.4E6 & G=1E6 + E=2.4E4 a G=1 E4 + E=2.4E2 a G=1E2 + E=2.4EO & G=1E0 -e—No Embedded Layers 12 14 Figure3.7- Normalized Transverse shear stress distribution for different bonding conditions obtained by varying the material properties of the embedded layers of a [0/1/90/1/0] composite plate. [I — isotropic] 05 -¢—E=2.4E6 & G=1 E6 -8- E=2.4E4 & G=1 E4 + E=2.4E2 & G=1 E2 +E=2.4E0 & G=1E0 +No Embedded Layers .0 055 OJ 015 ___0& 1y: (a/2,0,z) 1 1 025 03 035 Figure3.8- Normalized Transverse shear stress distribution along the height of the laminate for different bonding conditions obtained by varying the material properties of the embedded layers of a [0/1/90/1/0] composite plate. [I— isotropic] 27 04 The through-the-thickness distributions of transverse shear stresses at the center of the [0/90/0] composite plate for different values of properties of embedded layers are shown in figures (3.7) and (3.8). As against section 3.1.4.1, the embedded layer properties are isotropic in nature and hence their Young’s Modulus and Shear Modulus are related by poison’s ratio. Thus the variations of G and E are proportionate to each other. As we decrease the values of E and G of the embedded layers, the interlaminar shear stresses also kept decreasing. We can observe that when the E and G of the embedded layers is of the same order as that of the other layers of composite, the shear stresses values coincide with the case when there are no embedded layers. In other words, this is a perfect bonding condition. Moreover, we can also notice that when the Young’s Modulus and shear Modulus is of the order of 100 then the shear stress vanishes in the embedded layers. This happen only when delamination takes place [16] and all the conditions in between represent imperfect bonding conditions. 0-5 I T I I I I I —1 1 + E=2.4E6 a. G=1 E6 0 4 _ —a— 522.454 8. G=1 E4 ' + E=2.4E2 & G=1E2 4‘ + E=2.4EO & G-1 130 0.3 — . 0.2 — _ 01 1— _. N O '- . ~‘ - —1 -0.1 .. 1h; -' a r _ _7 Iii . __ -0.3 ~ ‘ ‘\ i -0.4 ~ - ‘0 I l 1 l l 1 l l l 13.1 -0.08 -0.06 -0.04 -002 o 0.02 0.04 0.06 0.08 0.1 I7(0,b/2,z) Figure3.9- Normalized in-plane displacement along the height of laminate for different bonding conditions obtained by varying the material properties of the embedded layers of a [0/1/90/1/0] composite plate. [I - isotropic] 28 0-5 fi fi 1 l 1 fl + E=2.4E6 81 G=1 E6 0 4 p + E=2.4E4 & G=1 E4 L ' -v- E=2.4EZ 81 G=1 E2 + E=2.4E0 8. G=1 E0 0.3 — +No Embedded Layers ' 0.2 - A 0.1 — ‘ n N 0 — — -0.1 — - -O.2 - — -0.3 - — -0.4 ~ — -0 1 1 1 h n 1 1 3.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 V(a/ 2 ,0, z) Figure3.10— Normalized in-plane displacement along the height of laminate for different bonding conditions obtained by varying the material properties of the embedded layers of a [0/1/90/1/0] composite plate. [I - isotropic] Figures (3.9) and (3.10) show how the transverse displacements at the center of the [0/90/0] composite plate vary as the bonding condition changes from perfect to imperfect. We can see from the above figures that as the Young’s Modulus and Shear Modulus starts to go below the order of 104, the variation of displacement across each embedded layer is not zero. In other words, the displacement across the interfacial layers is discontinuous. This supports our observation from the figures (3.7) and (3.8). Appendix C covers all the remaining figures for the stresses and displacements. We shall now discuss the second approach to deal with slippage or discontinuity on the laminate interfaces to simulate poor bonding as well as delamination. 29 3.1.6 Slip Approach When interfacial damage takes place across the interfacial layers, the continuity of the interfacial displacements will be lost. The interfacial displacement continuity equations used in the original Pagano’s problem can also be modified using the linear shear slip theory and linear normal separation theory to include slippage or mismatch on the laminate interfaces for simulations of poor bonding as well as delamination. A. Linear Shear Slip Theory — If a composite interface is non-rigid or imperfect, the in- plane displacements across the adjacent layers are going to be different. This kind of displacement discontinuity is also known as shear slip which may lead to shearing mode delamination [26]. 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"(x i)—u’(x h—l)-D z" (x h—l)-D r"(x ..Ii) (316) 3)” 2 9y, 2 5x R 9y, 2 5x xz 3y) 2 ’ _ u I l u v" (x,y.—:—) — v’zi The superscript (i) represents for the layer number, i.e., the ith layer of the composite laminate, and h,- is the thickness of the layer. As shown in figure 4.2 below, (A/21 and IQ.- are the node values of the displacements U and V at the point (x, y, 2,.) in the layer (i), while IA! 21+: and 1921+: are the node values at the same point but in layer (i+1). Moreover, since the layers are imperfectly bonded, the displacements values are not continuous A A A A across the layers and hence U21and Vziare not same as U21+1 and V21+1. However, 8’3 and T’s are the first derivatives of U and V with respect to z axis, respectively. More specifically, 52.- and f2.- 40 ‘1 73--» ..l- Figure 4.2 Nodal Variables for the composite laminate represents the node values of t? and g at the point (x, y, 2,) in layer (i) while S2111 and Z Z T 21+1 are at the same point but in layer (1' +1). Moreover, although the interlaminar shear stresses are continuous across the interface, the interlaminar shear strains are not, unless the material properties in layer (i) are identical to those in layer (i+1). Consequently, A A A A S 2.- and T 2.- are not the same as $21.1 and T 21.1. 4.1.2 Fundamental Equations The linear strain-displacements equations using equation (4.1) are written as follows: 41 - 6u 6021-1 ally—1 5&2: alim- 8(1)=_=___ +__ +__ +— x 6x 6x ¢1 6x ¢2 6x ¢3 6x ¢4 E'(v1)=fl=aV21—1¢l+6521—1¢2+6V21¢3+6521¢4 6y 6y 6y 6y 6y s“’=0 J(c,;’)=fl_*_6_u=6Uzi-1(151+6T2H¢Z+6U21¢3+6th(154 5x 6y 5y 5y 6y 6y 61321-. 6521—1 6192.- 53..- +—— +——— +— +— 6x ¢‘ 6x ¢2 6x ¢3 6x ¢4 7:;)=iw_+@:aV2i-l ¢l+aSZI-l¢2+aV2i ¢3+§fi¢4+fl 6y 62 62 62 62 62 6y - au 6w 6&21-1 5&21-1 50/21 55:21 "’=—+—= +— +-——— +——— .2 62 6x 62 ¢‘ 62 ¢2 62 ¢3 62 For any composite laminate, the constitutive equations are given by: {a} = [618} r a (1‘) -— (i) r N (i) I Q11 912 913 ém— 8" 0' 0y] _ 512 Q22 Q23 Q26 8.” } 0' _ _. _ _ z 913 923 _Q_33 Q36 82 [Try _Q16 Q26 Q36 Q4151 We 1' (1') E Q (i) 7 (i) {In} =|:§45 @55] {7x2} ¢4+— 6w 6x (4.3) (4.4) (4.5) (4.6) (4-7) (4.8) (4.9) In the above notations, the superscript (i) pertains to layer (i). [Q] refers to the stiffiress matrix and it relates the stress and strain matrices. The various coefficients Q, forming the stiffness matrix are known as the stiffness matrix coefficients. 42 4.1.3 Boundary Conditions The following boundary conditions are used in arriving at the solution. A. The conditions on the top surface of the composite plate are given by h 0. (an?) = q(x, y, t) (4.10) h h 1,, (x,y,-2-) = 5.0.13.5) = 0 (4.11) B. The bottom surface is traction free, hence both the normal stress and the shear stress vanish,i.e. h h h 0' , ,—— =2}, , ,—— =rvz , ,—— =0 4.12 .(xy 2) (xy 2) ,(xy 2) ( ) 4.1.4 Interfacial conditions The shear stress is assumed to be continuous through out the composite laminate with various bonding conditions and hence the stress continuity equations remain the same as equations (3.7) mentioned before. Due to delamination, the continuity of interfacial displacement components is lost and hence in order to take that into consideration two different approaches were taken which were similar to the ones taken to solve the modified Pagano’s problem. The first approach makes use of the linear shear slip theory and the second one is the embedded layer approach 4.1.4.1 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 43 laminate interface. The displacement field as mentioned in equations (4.1) assumes that the normal displacement remains constant through the thickness at a particular location in the x-y plane of the laminate. This is because it is expected that oz is extremely small and in most of the cases, it can be neglected. Hence, normal displacement separation theory was not included in the finite element analysis. Using the linear shear slip theory we get, [12141—1921 = 1331'}: i: 1,2,...,n—1 (4.13) A A V21+1— V21 = 03%;? This can also be written as, [12M = (321+ D3132) 1: l,2,...,n —1 (4.14) A A V21+1 = V21+ DASH? Here, (A1214: and Uzi represents the node value of U in the layers (1‘) and (i+1). Similarly, (i) and 7(1) 12 )2 V2.41 and V2.- represents the node value of V. t represents the shear stress in (1+1) :2 the layer (i) and are of same value as r and If” because of the shear stress continuity. BL? and Di: represents the shear slip constants in the layer (i). When these constants vanish, equation (4.14) reduces to, U21+1 = U21 , A A l=1,2,...,n-l (4.15) V21+1 = V21 which is perfect bonding conditions. As the value of these shear slip constants increases, the discontinuity between the displacement values across the adjacent layers of the composite also increases and in this way, one can simulate the various degrees of imperfect bonding by the variation of these shear slip constants. Finally, when these constants approaches infinity, delamination takes place. Moreover, by visualizing the equations (4.14), we can also conclude that for each nodal point, we require only one node value for the displacement in place of two as they are both related to each other by the shear slip constants and hence we can denote them as follows: (921' = (All A A i=l,2,...,n—l V21 = V1 (4.16) U1 =Uo U2n=Un A A ’ A V1=Vo V211 =Vn 4.1.5 Generalized Strain Solutions Considering layer (i) and (1' +1), if we combine equations (4.7), (4.8) and (4.9), we will get: 1' (1') a“) EU) 521‘ + '2:— EUH) @041) 521+] + ‘31:- if) = [44) 4;) ] A = [441) 41-111] A 6w (4.17) x2 Q45 Q55 T21+ _ Q45 55 T21+1+ _ 6x x Expressing equations (4. 17) in another way, we get, A '1 1') 6w “(1') —(i) ‘ — Szi _ Q44 Q45 Ty: _ 6y 418 A " —(1) —(1) 6w ( ' ) T2,. Q45 Q55 TU _ 6x 45 " —(1+1) —(1+1) “ (i) — S 21+1 _ Q 44 Q 45 T}: _ 6y ’ —(.-+1) —(i+11 (4-19) A 6w T 21+1 45 55 "Z — 6x On further simplification, we can get, A —1 1') 6w —<1) —(1) 1 — S21. = 1 [Q55) -Q(45:| {T12} _ 6y A —111—111 —1112 —(1 --11 T21 QSSQM-Q45 _Q45 Q44 sz _a_W 6x 4.20 6,, ( ) A A (11 r (1) a =[ ” '2] — i=1,2...n—l A12 A22 7,2 % 6x A . (.) aw. 51+ A A (1+1) TV, I — A“ =[ H 12] {} - 6y i=1,2...n—l (4.21) T21+1 A12 A22 TE 2“: 6x Now, using the boundary conditions (4.12) we get, 1§|+—— 0 i )zhmll 6y (4.22) 0 A 6w T1+— 6x This gives us, 311-3“: 63:: (4.23) T1 = —— 6x Similarly, on using the boundary condition (4.1 1) we will get, 46 [£271 =—ia"u‘) , 5y f‘Zn =—@ i 6x (4.24) If we combine the equations (4.1) of the displacement fields with equations (4.21), (4.23) and (4.24), along with the interfacial conditions expressed by equations (4. 14) and (4.16), we will get, i (1') _ (1—1) (1—1) (1) (H) (1) (1—1) W u —(sz Tn +U1—1 1+(A12712 +’422 7.12 _—]¢2 (1) (1) (1) (1) +U,.¢3 +(A12 r), +A22rn ——6x (154 (1) _ (1—1) (1-1) (1) (1-1) (1) (1-1) V — (D5? TYZ + VP] #1 + [All Tyz + A12 sz _ 6y ¢2 + 142.9341) + 414: g). i = 2,3,...n —1 A 1 6w 6w 7;). + 11.7. + [24:51:12 + 424;) - 5;)». [um = U0¢1 +( 5y v“) = Vo¢1 + [- %)¢2 + V1¢3 + (441(97):) + A372) - QEJW L f 6w (n)_ (n-l) (n-I) (n) (n—l) (n) (n-l)__ u - (D31 TE +Un-lyl +(A12 Ty: +A22 sz 6x ¢2 6w + Un¢3 + (— EJ¢4 (n) _ (n—l) (n-l) (:1) (71-1) (n) (n-l) _6W V - (D3,? sz + VII-l ”I + [All Tyz + A12 Tn 6y J¢2 6w V3 "_ 4 +.+¢[ ,y) L 47 (4.25) (4.26) (4.27) Equations (4.25) can be written in a further simpler form as follows f (I) _ (I) (I) (H) (I) (I-1) (I) (I) (I) (I) (I) (I) u —Al U,._,+A2 Ix: +A3 7.1: +A4 U,+A5 rm +A6 Ty, +A7 5; 4.. .. .. . .. .. .aw(4.28) v“) = A5014, + Agwrg-n + 4071;” + Am + Agar; + (”13+ .49) — l ' ' 5y The various constants A,,A,,A3,A,,,A5,A6,A7,AS,A9 are definedasfollows: A1“) = ¢1 i: l,2,3...n w: 0 i=1 ’ ng”¢,+A§;’¢, i=2,3...n m: 0 i=1 3 AW, i=2,3...n Ag" =¢3 i=l,2,3...n m_ Ago, i =1,2,3...n—l A5 - . 0 1=n AW— Ag)», i=l,2,3...n—l 6 _ O i=n Ail) = “(¢2 + ¢4) A: {0 i=1 D$”¢,+Al‘{’ 2 i=2,3...n (.,_ AW, i=1,2,3...n-l (4.29) A9 - . 0 I=n Thus, by using the displacement formulations given in equations (4.28) and (4.29), we can obtain the generalized strain formulations as follows: 48 1- (1'1) sf>_ 2:: AII) 51%. +Am 5212 ” +A§I>__.a;-Vz +49% x x x 62' (i) 51": - 62w + (I')__ +A(I’)_ + (I)_ 6x 6x A7 6x2 1 “—1) 8:3.)— :yv_ _ A10) a_gyi__—1+Ali) 6% l) + A8“) 62-6: + A4106??- "I . 61"? . 2 +A"’a—; +A.§"—a)':‘ +A§"§-—:v— 8:” = 0 . . . . (H) . “‘1’ 6r“) 7g): au+_a_v_:Al(I) a(]1—1+6V1-1 +Aé1)arxz +243?) arxz + )7 + 6y 6 6y 6x 6y 6x 6y (I) at (i) (1) 2 A"’[a——"U+ 6K]+ "’ a —+Ag" "" ——-+a-T—’“ +2A"’— 6w 6y 6x 6y 6y 6x 6x6y . 6r‘"” . 61‘” + (I)__.l_'“£_+ (I)_,1_'z_ 6x A9 6x 2) =_a_li+§r}_._—_Ui_l 6A1“) +122- 1) 6AZ i) "—"+T (.i_ 1) 6A“) —+Ui 6A4)” 62 6x 62 62 62 62 TU) a__A5i)+ z.(I’) 6A“) _+_a___A‘(li) 2W_+ 6W 62 ‘2 62 62 6x 6x (430) y(I)_ 6v+ 6w —V. 6A,("+T(H, 614i, +1.01) a_Asi)+ Via/1m )2 62+ 6y 1—1 62 .12 52 )2 62 62 51;» + a__4" 44;” _@+§3 “2 62 "7 62 62 6y 6y Using the above equations, total potential energy, kinetic energy and the interfacial energy of the laminate were calculated. Then using the Hamilton principle of minimum energy, finite element equations were derived by Lee. Hence, the thru-the thickness assembly was performed by hermite cubic functions as shown in the above equations 4.1- 4.30. The in-plane assembly was performed by using both C0 continuity and Cl continuity interpolation functions. For in-plane variables ,U,.,U. ,._ _1, V,,V,._ 1,1" ”,r“"’ r“) r“) bilinear interpolation functions were .112 9 xz’ yz’ 49 6w 62w 62w 62w _9 2 3 2 and " ' — , 12 telIIlS 6x 6x 6y 6x6y used while for the transverse variables w, interpolation functions were used. These finite element formulations were then implemented in a FORTRAN code. Lee’s FORTRAN code was used to analyze interfacial damage for [0/90/0] composite laminate. We used the same loading conditions and the simply support conditions given by equations (3.4) and (3.6) to obtain the stress and displacement variations along the height of the composite laminate. In order to simulate imperfect bonding of different layers, the values of the shear slip constants sz and D5). were varied from 0.0 to 1.0 and values of stresses and displacements were calculated. 0.9 T... . 7.. ._ I -. I I l . . i - fih ‘_ 0.7— ' ' , x. — . _ _ ,... --’ : " i +sz,Dsy-1E-11 . —.- "1 . _ " +sz,Dsy—1E—9 0.6 - r , 6“" _ ‘ —+—sz.Dsy-1 E-7 " ._ ”3.. + sz,Dsy-1E-5 O 5 _ . +sz,Dsy-1E—3 _ ' —e— sz,Dsy-1E-1 N 5'; -€—-sz,Dsy-1 E+1 0.4 " :' —< 0.3 L 5 i 1 f I' ‘7. ~ ,. -'_ . __ 1,1,. _ _ 0.1 — i i — ‘7‘?" . r -. H~'$fi’w :34.“ " _:L ’-~_—._ i l l l 3.5 | 0.5 2 2.5 {SIM/2.2) Figure4.3- Normalized transverse shear stress distribution along the height of [0/90/0] composite plate for different bonding conditions obtained by different shear slip constants D” and D”. Figure 4.3 and 4.4 show the variation of transverse shear stresses at the center of the edge of the [0/90/0] laminate as a unction of thickness for different bonding conditions. As we 50 start to increase the value of shear slip coefficients, the interlaminar shear stress goes on decreasing until it becomes zero when the values of D5, and D5), are 0.1 or 1.0. This is the point of delamination. As we go from perfect bonding conditions (sz and D3), =10'9) to delamination (sz and D5}, =1), transverse stresses vary significantly. 0.9 .. l l I r +sz,Dsy-1E-11 \; —¢—sz, Dsy-1 E-9 0.8 — . —l—sz,Dsy-1E-7 “ +szDsy—1E-S 0 7 _ +sz,Dsy-1E-3 _ ' +sz,Dsy-1E-1 . .. —e—sz,Dsy-1E+1 0.6 _ :2“ '1. ‘_ __ .— 0'5 ' i W 3‘ - . if .~ _ N E i 0.4 - ,2"-="'" ._' ' V _. ‘ — 7 = 6.; 4" wk 1 I 0.3 - f - 0.2 — _ 0.1 — I/ V — l l J l l 3.5 I 0.5 1 1.5 2 2.5 Tyz (a/ 2,0, 2) Figure4.4 -Normalized transverse shear stress distribution along the height of [0/90/0] composite plate for different bonding conditions obtained by different shear slip constants D“ and Dy), Figure 4.5 and 4.6 show the variation of transverse displacements at the center of the edge of the [0/90/0] laminate as a function of thickness for different bonding conditions the We can also see from figures 4.5 and 4.6 that as we go towards imperfect bonding conditions, the discontinuity of the interlaminar displacements goes on increasing. Appendix F gives the list of all the remaining figures for stresses. 51 0.9 0.8- 0.6— 0.4—- 0.3- 0.2— 0.1— l l 1 I .‘N {3‘23wa -'F*-- sz, Dsy—1 E1 1 ‘03: -¢- 05):, Dsy-1 E-9 ' —+—- sz, Dsy-1 E-7 - —v—— sz. Dsy— 1 E-S N + sz. Dsy—1 E—3 ‘ ,2. -6- sz. Dsy-1 E-1 ”‘33:. ——e— sz. Dsy—‘I E+1 - VT *3" ‘ ‘3'fiv‘r. m N V‘SV It? - _ tag: _ ‘~.‘;‘prw "'- .‘O‘WW, “ r', o 3;: — ‘ ‘\'r'. ~. “v.3 N — V‘vli‘so‘ Wm , Vz‘fi‘fiw «o v . 35:5;- i ' . 'svuw ._ l 92 ~05 6 5 1 m 15 2 U—(O,bf2,z) ' x‘IO'6 Figure4.5-Normalized in-plane displacement variation along the height of [0/90/0] composite plate obtained by varying the shear slip constants D”( and By 0.9 0.8— O.7~ 0.6~ 0.5~ 0.4- 0.3— 0.2» 1 l l I —B— szDsy-1E-11 —¢— DmDsy-1E—9 —+- 05x. Dsy-1 E-7 _ -—9— sz, Dsy-1 E-S + szDsy-1 E-3 49- sz, Dsy—1 E-1 —+— sz, Dsy-1 E+1 Figure4.6-Normalized in-plane displacement variation along the height of [0/90/0] composite plate for different bonding conditions obtained by varying the shear slip constants sz and D3), 52 4.1.4.2 Embedded layer Approach Another approach to take care of the discontinuities of displacements across the laminate interfaces is the embedded layer approach which we discussed in section 3.1.5. Embedded layers were introduced in the composite laminate. Each embedded layer was 0.11% of the total thickness of the composite laminate and its properties were isotropic in nature. By varying the Young’s Modulus and Shear Modulus of these layers, all the bonding conditions from perfect to imperfect were simulated. The different values of Young’s Modulus and Shear Modulus were modified in accordance with the relationships derived in section 3.1.7. The shear slip constants D” and D”, were set to zero. Following plots of displacement and stresses were obtained. 1 l 1 fl T 0.9 — ‘ 0.8 — ‘ 0.7 — ' 0.6 — + E=2.4E6 a. G=1E6 ‘ —e— E=2.4E4 a. G=1 E4 ~ 0.5 - + E=2.4E2 & G=1E2 — + E=2.4E0 & G=1E0 0.4 _ —<—>— E=2.4E-2 & G=1E-2 2 —a— E=2.4E-4 & G=1E-4 0.3 — ‘ 0.2 - ‘ 0.1 — l 3.5 f“ ' ' 0.5 1.5 2 1 2' x, (0, b / 2 , z) Figure4.7- Normalized transverse shear stress distribution for different bonding conditions obtained by varying the material properties of the embedded layers of a [0/1/90/1/0] composite plate. [I — isotropic] 53 Figures 4.6 and 4.7 show the variation of transverse shear stress across the height of the composite laminate. If we compare them with Figures 4.3 and 4.4, we find that they are very similar. For perfect bonding, the maximum percent error for In and In come out to be 1.5% and 0.44 % respectively. Similarly for the case of delamination, maximum the percent error for In and ryz is 2.15% and 2.16%. 1 I I I j +E=2.4E6 & G=1 E6 0 9 _ +E=2.4E4 & G=1 E4 L ° —-v— E=2.4E2 & G=1 E2 -+- E=2.4E0 8. G=1 E0 0.8 — -e-E=2.4E-2 & G=1E—2 ‘ +E=2.4E-4 & G=1E-4 0.7 — - 0.6 — — N 0.5 — ~ 0.4 ~ — 0.3 — — 0.2 — — 0.1 — _ l I l 0.5 1 .5 2 .85 5(a/ 2 ,10, 2) Figure4.8- Normalized transverse shear stress distribution for different bonding conditions obtained by varying the material properties of the embedded layers of a [0/1/90/1/0] composite plate. [I — isotropic] Figures 4.9 and 4.10 explain the variation of the transverse displacement across the height of the composite laminate. If we compare them with figures 4.5 and 4.6, we find that the maximum percent error in case of perfect bonding come out to be 2.14% and 0.52% for U and V respectively. Similarly for the case of delamination, the percent errors are 3.37% for both U and V, respectively. Hence we can see that both these techniques are equivalent to each other. Appendix G covers all the remaining plots for stresses 54 0.9 0.8 0.7 0.6I N 0.5 0.4 0.3 0.2 0.1 L +E=2.4E6 & G=1E6 -e— E=2.4E4 8. G=1E4 +E=2.4E2 & G=1E2 +E=2.4E0 8. G=1E0 —e-E=2.4E-2 8. G=1E—2 +E=2.4E—4 8. G=1E-4 -9. -1 -0.5 171b, b/2,z) x 10' 1.5 e Figure4.9- Normalized in-plane displacement distribution for different bonding conditions obtained by varying the material properties of the embedded layers of a [0/1/90/1/0] composite plate. [I — isotropic] 1 0.9 0.8 0.7 0.6 N 0.5 0.4 0.3 0.2 0.1 l -0.5 J l 0 0.5 V(a/2,0,z) —¢— E=2.4E6 & G=1 E6 —~a— E=2.4E4 & G=1 E4 —v— E=2.4E2 & G=1 E2 --+- E=2.4E0 & G=1EO + E=2.4E—2 & G=1 E-2 + E=2.4E-4 & G=1 E-4 .. '91 A Figure4.10- Normalized in-plane displacement distribution for different bonding conditions obtained by varying the material properties of the embedded layers of a [0/I/90/I/0] composite plate. [I - isotropic] 55 5. VALIDATION OF FORTRAN CODE In this thesis, four different approaches are pursued to analyze the interfacial damage in a composite laminate plate. Two of them are analytical approaches and make use of Pagano’s solution and the other two are numerical approaches and make use of Lee’s Laminate theory. This is clearly shown in the following figure. Interfacial damage Analysis I f | Modified Pagano’s Laminate Theory Solution l—L—I l——;l Slip Approach Embedded Layer Embedded Layer Shear Slip Approach Approach Approach Figure 5.1 Approaches to analyze interfacial damage in a composite The solutions proposed by Pagano using the slip approach or the embedded layer approach are exact solutions and hence can be used to validate the numerical solutions using the laminate theory. In order to achieve this end, firstly, the solutions given by Pagano for perfect bonding conditions were used to validate the FORTRAN code. Secondly, solutions in case of imperfect bonding obtained by the modified Pagano‘s problems were matched against the solutions by laminate theory. Based on the comparisons, appropriate conclusions regarding the validity of the Laminate theory were derived. Also different situations for the utility of the Embedded layer vs. the slip approach were analyzed. 56 5.1 Perfect Bonding Conditions Consider a [0/90/0] composite laminate as shown in figure 3.1. In order to attain perfect bonding conditions using the slip approach, the shear slip coefficients, Der and Dy), and the normal separation coefficient, k were set to zero appropriately. In the case of embedded layer approach, the young’s modulus E and the shear modulus G of the embedded layers were set equal to that of the matrix. Thus our composite becomes [0/1/90/1/0], where ‘1’ stands for isotropic embedded layers. The thickness of the embedded layers was set equal to 0.001 in, which was 0.11% of the total thickness. The loading conditions and the simply support conditions were set to be same as that of Pagano’s problem given by equations (3.4), (3.5) and (3.6). Using these 4 approaches, all stresses and displacements were obtained. For each particular stress and displacement, all the four cases were plotted on the same figure as shown in the figures below. 0-5 1 r mi 1 1 0.4 — 0.3 — —e— Pagano-Embedded Layer + Pagano-Shear Slip 0.2 ~ +issct—Shenr Slip +issct-Emboddod layor -o.3 — a .0.4 -— .4 l 1 0.4 0.6 0.8 - l l 1 053.8 -O.8 -O.4 -0.2 .5— {2 0x (a/2,b/02,z) Figure 5.2-Comparison of Normalized in-plane stress along the height of laminate between Lee’s solutions and Pagano’s solutions of a perfectly bonded |0/90/0] composite plate 57 0.5 r I l 1 —-o— Pagano- Emboddod Layer -0-— Pagano-Shear Slip 0,4 _— +issct-Shoar Slip +issct-Emboddod layer 0.3 >— J 0.2 '- -I 0.1 >— E N 0 I— .. -0_2 r— _. -o3— — -oa~ ~ r 0.4 0.6 _ L 1 l 03.8 -O .8 -O.4 -0 .2 ay(cf/2,b/2,z) °'2 Figure 5.3-Comparison of Normalized in-plane stress along the height of laminate between Lee’s solutions and Pagano’s solutions of a perfectly bonded [0/90/0] composite plate Figures 5.2 and 5.3 show the normalized in-plane stress at the center of the [0/90/0] composite plate as a function of thickness. The solutions obtained via all the 4 approaches are shown. We can notice that the accuracy of the Lee’s laminate theory, using either the shear slip or the embedded layer method, to predict the normal stresses is. good. The maximum error for ex and 0,, between the Lee’s laminate theory and Pagano’s solutions come around 10 % and 4 %. Figures 5.4 and 5.5 below give us the values of normalized transverse shear stresses at the center of the edge of the [0/90/0] composite plate across the height of the laminate. It can be seen that except at few places, the shear stresses predicted by the laminate theory matches well with that of the exact solutions. In case of normalized shear stress, ryz, the maximum error occurs at the center of the middle layer while in case of normalized shear stress, rxz, it occurs at the center of the top and bottom layers. The value of the normalized shear stress, ryz and rxz, predicted by laminate theory is 0.293 and for Pagano 58 it is 0.277. On the interface, the value of the normalized shear stress, ryz and In, predicted by Lee’s laminate theory is 0.258 and for Pagano, it is 0.252 0.5 T fi l I I I k 0.4 — — 0.3 h e 0.2 — ' ' 4 0 1 _ +Pagano-Embodded Layer _ ' + Pagano-Shear Slip +issct—Shear Slip N 0 ~ + issct—Embedded layer n -0.1 ~ ~ -0.2 — a -0.3 — - —0.4 — — .0 I 1 l l I I 4?.05 0 0.05 0.1 0.15 0.2 0.25 0.3 2'12 (0, b/ 2 , 2) Figure 5.4-Comparison of Normalized transverse shear stress along the height of laminate between Lee’s solutions and Pagano’s solutions of a perfectly bonded [0/90/0] composite plate 0.5 l l l l I -6- Pagano-Embedded Layer 0 4F + Pagano-Shoo: Slip ' ‘ +issct-Shoar Slip + issct-Embedded layer 0.3 - - 0.2 — i 0.1 ~ 2 N 0 — .. -0.1 I— - -02 — —I -0.3 — ~ -0.4 - -l _0 l l l l l 1 $05 0 0.05 0. 0.15 0.2 0.25 0.3 1 __ 1y, (a/ 2 ,0, 2) Figure 5.5-Comparison of Normalized transverse shear stress along the height of laminate between Lee’s solutions and Pagano’s solutions of a perfectly bonded |0/90/0] composite plate 59 0-5 I I I l I I T I 1 —e— Pagano-Embedded Layer 0 4 _ + Pagano-Shear Slip H ‘ -e—issct-Shear Slip + issct-Embedded layer 0.3 — 0.2 — 0.1 — N 0 — -0.1 — -O.2 — -0.3 — -0.4 — _o I l l l l ' l l l l 3.01 -0.008 -0.006 -0.004 -0.002 0 0.002 0.004 0.006 0.008 0.01 U(O, b/ 2 , 2) Figure 5.6-Comparison of normalized in-plane displacement along the height of laminate between Lee’s solutions and Pagano’s solutions of a perfectly bonded |0/90/0] composite plate 0-5 l l l f l l T I I -0- Pagano-Embedded Layer 0 4 _ + Pagano-Shear Slip 1 ' +issct-Shea Slip + issct-Embedded layer 0.3 F 0.2 — 0.1 — N 0 — -0.1 — -02 _ -o_3 _ -0.4 — _o 5 l l l l l l l l I -0.025 -0.02 -0.015 -0.01 -0.005 _0 0.005 0.01 0.015 0.02 0.025 V(a/ 2 ,0, 2) Figure 5.7-Comparison of normalized in-plane displacement along the height of laminate between Lee’s solutions and Pagano’s solutions of a perfectly bonded |0/90/0] composite plate Figures (5.6) and (5.7) display the variation of normalized in-plane displacements at the center of the edge of the [0/90/0] composite plate across the height of the composite. The 60 maximum error for displacement U between the Lee’s laminate theory and Pagano’s solutions occur at the center of the middle layer, while for displacement it is in the top layer. The value of the normalized displacements U and V, predicted by Lee’s laminate theory is 0.0028 and -0.0159 while for Pagano, it is 0.00225 and -0.015 Another important thing to notice from all the above figures is that the solutions proposed by the Pagano- Embedded layer method match exactly with that of the Pagano-Slip approach. The same is true with ISSCT-Embedded layer method and the ISSCT-Shear slip approach. Since the laminate theory assumes normal stress to be zero and normal deflection to be constant, there are no plots shown for them. The error between the Pagano’s solutions and Lee’s laminate theory’s solutions could be attributed to two factors. Firstly, Lee’s laminate theory assumes that the normal displacement W (z) is constant across the height of the laminate based on the assumption that the normal stress in the thickness direction is negligible compared to other stresses. Secondly, it was also observed that as we refine the meshing of the plate in the finite element solution by increasing the no of elements, the accuracy of the solution increases. Figures 5.8 and 5.9 display the variation of normalized transverse shear stresses at the I center of the edge of the [0/90/0] composite plate using the ISSCT-Shear Slip Approach with different no of elements. We can clearly see that the solution accuracy improves with using 100 elements as compared to 36 elements. In the case of shear stress, ryz, the maximum percentage error by 100 elements is 4.41% as against 6.47 % with 36 elements. . In the case of shear stress, In, the maximum percentage error by 100 elements is 5.47% as against 9 % with 36 elements. 61 0-5 1 T 1 T I I j wk : - . +Pagano 0.4 ~ ' ' 3 3 : ; . _ +100 elements _ —l—36 Elements 0.3 — _,_ . .- _ 0.1 - _ N 0 — _ -o.1 - g -0.2 — - I} J -o.3 - " . ' _ -o.4 — _ '96, 05 0 o be 011 o is 012 o 125 013 o 35 5(0, b/2 , 2) Figure 5.8 Comparison of normalized transverse stresses obtained by 36 elements vs. 100 elements 0.5 l l l l l I -O- Pagano —e— 100 elements . -|—36 Elements 0.4 — 0.3 — 0.2 — 0.1 — 4305 l I l l l - .05 0 0.05 0.1 015 0.2 0.25 0.3 {la/2,0,2) Figure 5.9 Comparison of normalized transverse stresses obtained by 36 elements vs. 100 elements 62 5.2 Imperfect Bonding Conditions In order to validate the Lee’s laminate theory’s solutions with that of Pagano for imperfect bonding conditions, we considered one particular case of delamination. In order to simulate this particular case using the embedded layer approach, the shear modulus G of the embedded layer was set equal to 0.01 psi and the corresponding Young’s modulus E was set equal to 0.024 psi (assuming the poison’s ratio of the embedded layer to be equal to 0.2). In order to simulate the same condition using the slip (or shear slip) approach, the shear slip coefficients, D” and D”, were set equal to 0.1 and the normal separation coefficient, k was set to 0.0417 appropriately using the relationships defined in section 3.2. Figures 5.10 and 5.11 show the normalized transverse stresses at the center of the edge of a [0/90/0] laminate as a function of thickness. We can notice that for imperfect bonding conditions the results calculated by the laminate theory do not match well with that of the Pagano’s theory. This is expected because for the laminate theory, the normal deflection is assumed constant through the thickness. But for Pagano's solution, the variation of normal deflection W with respect to thickness is always considered. As a result of this, there is a loss of contact between the layers. As we see from the figures 5.10 and 5.11 for Pagano's solution, only the top layer takes the load for poor bonding. The other two layers are basically free of stress. That is the reason we have much bigger normal deflections as it can be seen in figure 5.12. On the other hand, the ISSCT keeps the interfaces in contact, although shear slip might occur. The load passes thru all the three layers. 63 0.5 0.4 0.3 0.2 0.1 '°:3. —l— Pagano-ShearSL'p 8. NormalSlip + Pagano-Embedded Layer ‘ -0— lssct-ShearSlip +Issct-Embedded Layer E l J l l l l l 0 0.2 0.4 0.6 0.8 1 1.2 1 4 1,, (0, b/2 , 2) Figure 5.10 Comparison of normalized transverse shear stress along the height of a delaminated [0/90/0] plate between ISSCT and Pagano [sz=Dsy=0.l, k=0.0417, E=0.024 psi, G=0.01 psi] 0.5 l I I l 1 I —|'- Pagano-ShearSLip 8. NormalSlip 0 4 +Pagano-Embedded Layer ' -O- Issct-ShearSlip + lssct-Embedded Layer 0.3 _ 0.2 1 0.1 ~ N O - -0.1 — -0.2 - -0.3 .. -0.4 _ -0. 1 1 4 I J 1 3. O 0.1 0.2_ 0.3 0.4 0.5 0.6 2'},z (a/ 2 ,0, 2) Figure 5.11 Comparison of normalized transverse shear stress along the height of a delaminated [0/90/0] composite plate between Lee and Pagano [sz=Dsy=0.1, k=0.0417, E=0.024 psi,G=0.01 psi] 64 16 T l l l I l I I l 14~ -t- Pagano-Shear Slip and Normal Slip -*- Pagano- Embedded Layer -e— ISSCT-Shear Slip + lSSCT-Embedded Layer SI 10— J l l J -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 log(sz )or log(D_,_,.) Figure 5.12- Normalized deflection calculated at the center of top layer of |0/90/0] composite plate for different bonding conditions. Although these differences exist between the Pagano’s solutions and the Laminate theory approach, there are certain similarities existing between them which make us utilize the laminate theory. We can notice from the figures 5.10 and 5.11 that the stresses at the interface of the adjacent layers in case of slip (or shear slip) approach and stresses in the embedded layers in case of embedded layer approach vanish. Thus the laminate theory predicts the presence of delamination accurately. This same thing is confirmed in the figures 5.13 and 5.14 where the laminate theory clearly shows that the in-plane displacements across the layers become discontinuous in case of shear slip approach or vary significantly across the thickness of the embedded layer, thus predicting delamination. 65 0.5 T 1 1 l l . —1— Pagano-ShearSLip 8. NormalSlip 0 4 _ -' ‘ —*- Pagano-Embedded Layer I ' ' -0- lssct-ShearSlip + lssct—Embedded Layer 0.2 — \ _ -03 l— 4 —0.4 - e _0 n l J l 1 l 1 4?.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 (7(0, b/2 , z) Figure 5.13-Comparison of normalized in-plane displacement along the height of a delaminated [0/90/0] plate between Lee and Pagano [sz=Dsy=0.l, k=0.0417, E=0.024 psi,G=0.0l psi] 0-5 I I r 1 I I I -t— Pagano-ShearSLip 8. NormalSlip O 4 +_ -“— Pagano-Embedded Layer J ' -6- Issct-ShearSlip + Issct-Embedded Layer 0.3 — _. 0.2 — — 0.1 — 1: _ I» is N O — a It: «I l—— } -—i -0.1 F -o.2 ~ 1: a 0 1t ~O.3 — — 1r * Ir .0 4 _ I\ 3 . 0 n -0 l l 1 l l 1 l 4?.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 V(a/2 ,0, 2) Figure 5.14- Comparison of normalized in-plane displacement along the height of a delaminated [0/90/0] plate between Lee and Pagano [sz=Dsy=0.l, k=0.0417, E=0.024 psi,G=0.01 psi] Another important thing to conclude is that although the values of the stresses predicted by the laminate theory do not match exactly with that of the Pagano’s solutions, they are 66 of the same order and always lesser then the maximum value of the stress calculated by Pagano. Also we can observe that the solutions obtained by the Pagano-slip approach matches exactly with that of the Pagano-embedded layer approach. Same is true in the case of laminate theory. Hence we can conclude that slip (or shear slip) approach is equivalent to embedded layer approach in many cases. However, as we will discuss in section 5.3, when we have situations which do not give the provision for normal slip, the embedded layer approach cannot be used. Hence, in order to validate the shear slip model given by the laminate theory, we cannot use the embedded layer approach. This is because the inherent assumptions of constant normal displacement in the laminate theory only allow modeling the shear slip and not the normal slip. Thus we will be using the Pagano-slip approach to validate the shear slip modeling of laminate theory is section 5.4. Appendix H covers all the remaining figures for stresses and displacements for section 5.2 5.3 Validation by Pagano’s Slip Approach In order to simulate only the shear slip conditions for perfect and imperfect bonding conditions using the Pagano’s slip approach, the shear slip coefficients, D” and D8. can are varied in between 0.0 and 10.0 while the normal separation coefficient, k is set to zero appropriately. We noticed before in figures 3.10 and 3.11, that for [0/90/0] composite laminate as soon as the shears slip coefficients, D” and D53, have a value of 1.0, delamination takes place. In creasing the value of the shear slip coefficients, D3,. and 12,, above 1.0 will only be simulating delamination. Similar thing can be done in case of Lee’s laminate theory to simulate the perfect and imperfect bonding conditions. In the 67 following figures, two particular cases of imperfect bonding were analyzed and the results of the Pagano-shear slip approach are compared with the results predicted by the laminate theory. These two cases of imperfect bonding were simulated by setting the shear slip coefficients D” and D3}, to 1) 10'7 2) 10. Figures 5.15 and 5.16 show the variation of the transverse stresses at the center of the edge of a [0/90/0] composite plate as a fimction of thickness. We can see that the results match fairly well. In case of normalized shear stress, ryz, the maximum error occurs at the center of the middle layer while in case of normalized shear stress, In, it occurs at the center of the top and bottom layers. When the shears slip coefficients, D” and DH are 10.0, the value of the normalized shear stress, ryz and In, predicted by laminate theory is 0.465 and for Pagano it is 0.425. On comparing figures 5.15 and 5.16 with figures 5.10 and 5.11, we find that the stress variation across the layers changes significantly when we have normal separation across the layers. Due to normal separation, the load is not allowed to be evenly distributed across all the layers. Figures 5.17 and 5.18 below show the variation of the normalized in-plane displacements at the center of the edge of a [0/90/0] composite plate as a function of thickness. 68 0.5 0.4 0.3 0.2 0.1 Figure 5.15 Comparison of normalized transverse shear stress along the height of a [0/90/0] composite plate for 2 different bonding conditions 0.5 0.4 0.3 0.2 0.1 Figure 5.16 Comparison of normalized transverse shear stress along the height of a |0/90/0] L —o— Pagam-sz=Dsy=1E-7 - —e— lSSCT-sz=Dsy=1 E-7 _ -1—Pagmo-sz=Dsy=1 E1 A + lSSCT-sz=Dsy=1 E1 — —1 — 1 l l l l J l .1 0 0.1 0.2 0.3 0.4 0 5 1,, (0, b/2 , z) 0.6 composite plate for 2 different bonding conditions 69 j l l l l I +Exact-sz=Dsy=1E-7 _ +lssct—sz=0sy=1E-7 . -l—Exact-sz=Dsy=1 E1 +lssct-sz=Dsy=1 E1 _ I l l l l l l .1 0 0.1 0.2 0.3 0.4 0.5 0.6 zy, (a/2 ,0, z) 0.5 I | I l l " _ +Exact-sz=0sy=1E-7 0.4 — . +lssct-sz=Dsy=1E-7 — +Exact-sz=Dsy=1 E1 +lssct-sz=0sy=1 E1 0.3 ‘ .-\ .1 02— \gp _ _ \ . -r N» » ‘ we" " 'V 0.1 ~ " .4, — .1. I, l N 0 - I a .4',’ I!- -0.1 - ”1' - ._ ,i' \I -0.2 - ‘ -— ~0.3 - ~ -0.4 — d _0 1 l J I l 3.03 -0.02 -0.01 0 0.01 0.02 0.03 U (0, b/ 2 , 2) Figure 5.17 Comparison of normalized in-plane deflection along the height of a [0/90/0] composite plate for 2 different bonding conditions 0.5 T I I I I " xx \ +Exact-sz=Dsy=1E-7 0.4 — ' ‘\.- . \w‘ +lssct-sz=Dsy=1E—7 - \~ K, ‘ \“Kw -l—Exact-sz=Dsy=.1E1 0 3 h ‘\.,, xxx +lssct-sz=DsF1El _ ' s‘ \‘ \ ..‘\ \\ 0.2 “ ‘ ‘\\ . xxx _ i \: .‘ \" or- ‘ - N 0 ’— ‘.\ _ -0.1 - — -0.2 - - -0.3 - - -0.4 F —~ -0 l I l l 1 05.03 -0.02 -0.01 0 0.01 0.02 0.03 V(a/2 ,0, 2) Figure 5.18 Comparison of normalized in-plane deflection along the height of a [0/90/0] composite plate for 2 different bonding conditions 70 Figure 5.19 below shows the variation of the normalized normal deflection as a function of the shear slip coefficients D” and D3)” We can see that the results match quite well. It was found that the percentage error comes out to be around 2.7%. Appendix I covers all the remaining figures for stresses and displacements for section 5.2 5-5 l l l l l I I I l +Pagano-Shear Slip + ISSCT-Shear Slip 4.5 l § | 3.5 2.5 l l ?9 -8 -7 -6 l l l J J -4 l “(I -3) -2 -1 0 1 log D” )or log(Dsy Figure 5.19 Normalized deflection calculated at the center of the top layer of the [0/90/0] composite plate for different bonding conditions 71 6. BLAST SIMULATIONS Advanced composites have found tremendous amount of utility from aircrafis and submarines to pressure vessels and automotive parts. In all of these applications, the composite structures are subjected to different kinds of loading. Amongst all of them, blast loading is one of the most significant types of loading. Especially, when it comes to safety of people involved, it is of critical importance. In order to study interfacial damage of the laminated composite plates subjected to blast loading, we incorporated the blast loading in Lee’s FORTRAN code. In order to do that we need to have a very clear understanding of the pressure profile created due to blast wave. This will be discussed in detail in section 6.2 6. 1 Experimental Work [39] In order to simulate the air-blast loading in the laboratory, a shock tube was used. The shock tube consists of combining a shock wave generator and a wave transformer and is capable of providing hi gh-pressure, high-temperature and hi gh-velocity shock waves. Figure 6.1 shows a photograph of a shock tube housed at Michigan State University. The shock wave generator has a constant cross-sectional area. The outer diameter of the tube is 120 mm while the inner diameter 80 mm. The total length of the shock wave generator is 6 m. The left section, 2 m long, stores a relatively high pressure gas and is called the high-pressure chamber. The right section, 4 m long (two sections, each 2 m, joined by flanges), stores a relatively low pressure gas and is called the low-pressure chamber. 72 Figure 6.1 Shock wave generator at MSU [39] In order to produce a shock wave pressure, the diaphragms between the high-pressure chamber and the low-pressure chamber need to be removed instantaneously. The diaphragm chamber is about 40 mm long. It is bounded by two diaphragms, one at each end, and is used to store a gas with a pressure equal to the average of the high-pressure gas and the low-pressure gas. The pressure waves produced from the shock wave generator are of one-dimensional nature and hence a pressure transformer is used to convert these one-dimensional shock waves into spherical or hemispherical blast waves to simulate real blasts. Figure 6.2 shows a schematic diagram of the shock tube which depicts changing the one- dirnensional shock wave into a hemispherical blast wave. A [0/90/90/0] composite plate of 2.4mm thickness and 127m in diameter was mounted inside the Protection vessel as we can see in figure 6.2. The mounting was designed to provide fully clamped boundary conditions. Hence there are no displacements in the in- plane and out-of-plane directions on the boundaries in the tests. 73 Specimen Protection Protection Holder cart vessel High Double- Low Second preaaure diaphragm pressure diaphragm compartment compartment compartment 8. Nozzle Cable connector 1 Signal measurement ‘ so -on I Preseure control 8. _. a.‘ Measurement section I . - WE ....-;. , Control 8. Air compressor Display box 8. Tanks Figure 6.2 - A funnel-type wave transformer [39] The [0/90/90/0] composite plate had the following properties E1] = 38.6 GPa, E22 = 8.27GPa, E33 = 8.27GPa, Gl2 = 4.14GPa, G13: 4.14GPa, G23 = 2GPa v12 = v23 = v13 2 0.26 (6.1) X , = 1062 MPa, X c = 610MPa, S = 72MPa Y,=3lMPa YC=118MPa Here X,, X0 and Y, , Yc refers to the strength of the material in tension or compression in the direction of material coordinates. S is the interlaminar shear strength of the material. The air blast pressure was measured before the test at regular intervals of time using a pressure transducer, P6, fixed at the same position as the composite plate and is shown in figure 6.3 below. The signals obtained from the pressure transducer were amplified by 74 using a charge amplifier and the pressure time profile of figure 6.4 was obtained. This experimental pressure-time 4 '(‘t' I ~ "1— \\\\\‘ 7]. Iv. ' III/III” \\\'\l\“ /////////////A II; Figure 6.3 Location of pressure transducer [39] profile was modeled using various phenomenological models proposed by several authors and was incorporated into the FORTRAN code. That will be shown later in section 6.2.2. The nozzle of the shock tube, as can be seen in figure 6.2, is having a diameter of 12.7mm and hence the blast loading was concentrated only on that area of the plate at the center of the plate. P5 & P6 at '1500-110'psi 0.233, 32.2 P ' P 1.73, 25.65 -a-—- redlcted 6 Pressure(MPa, P N o L; -10-8-6-4-202468101214161820 time(ms) Figure 6.4 Loading history for fmite element simulation [39] 75 6.2 Analysis of blast propagation A blast wave is a pressure wave of finite amplitude generated due to a rapid release of energy in the medium. These waves are accompanied with a transient change in the gas density, pressure and velocity of air surrounding the explosion point. Regardless of the source of the initial pressure disturbance, the properties of the medium as a compressible gas will cause the front of this disturbance to steepen as it passes through the air, until it exhibits nearly discontinuous increase in pressure, density and temperature. The resulting shock front moves supersonically faster than sound speed in the air ahead of it. Figure 6.5 below shows that more clearly. Compressed Medium Explosive gas Wave front Figure6.5 Blast Wave When these blast waves impinges on any surface, the magnitude of loading on the surface at any instant depends on the type of the explosive, weight of the explosive, the distance traveled by the blast waves from the point of explosion to the surface and the time 76 elapsed from the moment when the blast waves first came in contact with the surface. They also depend on other factors like the medium of propagation. In this report we have considered air as the medium. There are basic equations which govern the transmission of blast waves through air. Solving these equations either numerically or deriving a solution empirically gives us the variation of pressure as a function of distance. In order to appreciate the solutions to these equations properly, two fundamental concepts of blast phenomenon are required to be understood - Equivalent TNT weight and scaling laws. Along with this, there are various phenomenological models explaining the pressure time variation of the blast waves, at a particular location, some distance away from the explosive source. 6.2.1 Equations for air blast transmission The air blast transmission equations are generally described in one of the three one dimensional cases, that is, cases in which the shock and flow fields are described in terms of a single spatial variable and time. This is because of the ease with which these equations can be solved. These cases are linear flow, spherically symmetric flow and cylindrically symmetric flow. Of these three cases, the one most applicable to blast waves in air is spherically symmetric. This case applies both to a spherical source far from any reflecting surface and to a hemispherical source located on a perfectly rigid reflector, both of which approximate a number of real blast sources. In this case, all quantities depend only on time t and the distance r from the origin of coordinates. All flows are radial with a single velocity component, u. the fundamental equations are: 77 Conservation. of momentum au 6a + 1 6p (6.2a) at 6r ,0 6r where p is fluid density and p is absolute pressure Conservation of mass (6.2b) a—’0+u(’3—’0+,o@+2£=0 6t 6r 6r r Conservation of energy as uaS (62C) —— + — = 0 6t 6r where S is entropy of a fluid element An equation of state is needed to complete the set of equations and is given by, p =f(p.S) (62d) In the steep gradients within shock fi'onts, the above equations are not all valid because heat conduction and viscosity effects become important. Hence, in blast theory, certain jump conditions are formulated called the Rankine-Hugoniot equations. These equations for a coordinate system moving with a discontinuity are given by: “IPI =u2p2 p1 + plulz = P2 + pzuz2 (6.3) el +-I—)'—+-1—u,2 =e2 +&+-1—u22 .01 2 ,02 2 Here subscripts 1 or 2 denote one side or the other of the discontinuity. 78 Many analytical and numerical solutions are being proposed to these equations. Some yield high accuracy. In order to appreciate those solutions better we will now discuss two very fundamental concepts which are very frequently used in Blast theory. One of them is Equivalent TNT weight and the other is Scaling law. A. Equivalent TNT weight From the time of Second World War, many researchers and scientists have done extensive testing for some of the explosives like TNT and pentolite and have amassed huge datasheets detailing the variations of different blast parameters with respect to time and distance. In order to avoid following the same process for all kinds of explosives, which will be very time consuming and expensive, they have defined the concept of equivalent TNT weight. Given the weight of the explosive we can find an equivalent TNT weight using the heat of combustion of the explosive. The formula for equivalent TNT weight is given by: WTNT : wexp (6.4) Where, Hexp is the heat of combustion of the explosive and Hnw are the heat of combustion of TNT. wexp and Wm, are the weight of explosives and equivalent weight of TNT. Using this equivalent TNT weight one can easily found out the values of all air blast parameters at a particular distance and time from the point of explosion. B. Hopkinson Scaling Law Experimental studies of blast wave phenomenology are often quite difficult and expensive, particularly when conducted on a large scale. Hence various investigators 79 have attempted to generate model or scaling laws that would widen the applicability of their experiment or analysis. The most common form of scaling law is Hopkinson or “cube root” scaling. This law states that self-similar blast waves are produced at identical scaled distances when two explosive charges of similar geometry and the same explosive, but of different size are detonated in the same atmosphere. It is customary to use as the scaled distance a dimensional parameter, Z Z - R (scaled distance) — Wl/3 t r = — scaled time 6.5 W ( ) ( ) 1/3 1/3 I =— scaledim ulse 4 w ( p ) Here, R is distance fi'om center of explosive source, I is the positive impulse of the blast wave, t is time and W is energy of the explosive. Z refers to scaled distance, 2' refers to scaled time and 9‘ refers to scaled impulse. In simple terms, this law states that pressure, velocity, scaled impulse and scaled time are unique functions of scaled distance. t I R Hence, given any explosive and its weight, we can find its equivalent TNT weight using equation (6.4). And then using equations (6.5) and (6.6) we can find everything about the various parameters of blast loading for that particular explosive. 80 Nearly all the data which are reported in the literatures uses the above two concepts. For instance, in their book [5], Smith and Hetherington gives us a plot of side on blast parameters for spherical charges of TNT. We find that the variation of all the blast parameters including pressure, velocities etc are given in terms of Hopkinson’s scaling law. 6.2.2 Pressure profile 6. 2. 2. 1 Pressure-time profile As the blast waves travel, the pressure profile varies as a function of the distance covered from the point of explosion as well as the time incurred. If we just consider the pressure time profile, the pressure at a particular location, exponentially decreases as time passes. This is depicted in the figure below: A p (0 APS i - I--------‘ Time ta+td++td- Q“.- N i Fig 6.6 Pressure Time History [7] Consider that an explosion has occurred at time t = 0. A pressure transducer was fixed at some distance from the point of explosion and it recorded the time history of absolute pressure. The record produced by such a gage is depicted in figure 6.6. For some time 81 after the explosion, the gage records atmospheric pressure p0. At the arrival time, ta, the pressure at that point suddenly jumps from atmospheric pressure, p0 to pressure at the shock front, P5. As time passes fi'om ta to ta+td+, the pressure decays to atmospheric value, p0 , then drops to a partial vacuum and eventually returns to ambient pressure. (P,- p0) or AP, is known as peak side-on overpressure or merely the peak overpressure and to is called the arrival time. The portion of the time history above initial ambient pressure, tf is called the positive phase of time duration, while the portion below ambient pressure, td' is known as the negative phase of the time duration. Similarly, positive and negative impulses are defined as, ta +1; 1:: flp(t)- Poidt ta (6.7) ta H; H; 1.? = flPo -p(t)]dt ta H; In other words, the area under the positive phase and negative phase of the pressure time curve is known as the positive and negative impulse. The slope of the pressure time curve at time t=ta is called the initial decay rate. We can notice from figure 6.4 that the pressure time profile measured in shock tube is little different from the usual pressure time profile obtained fi'om the regular blast explosions. First of all, the pressure rise is not instantaneous but it rises linearly from 20.18MPa to 32.2MPa. The equation governing this pressure rise can be written as, p(t) = 51.588t + 20.18 (6.8) Here, the unit of pressure is in MPa and the unit of time is in ms. 82 Also as time progresses, the pressure do not decay to ambient value of 0.1MPa but it decays to a constant value of 1.37MPa. Hence while modeling we need to take into considerations these factors. Moreover, in order to find out the various constants involved in the phenomenological models, we also need to use the initial decay rate and the impulse calculated from the experimental data of figure 6.4. For our case, the approximate value of initial decay rate and positive impulse can be calculated to be equal t0, d—p = —4.37542 MPa (6.9) dt ms This slope was obtained using the two point formula. Here the initial peak overpressure at time t=0.18ms was taken as the first point, while the overpressure at t=1. 73 ms was taken as the second point. I; = 103.0529 MPams (6.10) The area under the pressure time curve was calculated using the trapezoidal rule. In order to model the exponential pressure variation at a given location with progressive time, a number of authors had proposed different phenomenological models. Using each of these models we will try to model the pressure time profile of figure 6.4. Flynn had proposed a simple linear decay model for the pressure variation as follows: p(t)=p0+(PS—po{l—t—t7] 010atmos b = 'AP,(0.88 +0.072APS) AP, <10atmos _AP,(1.67 —O.011APS) AP, >10atmos c = 8.71+O.1843AP9 — 104 AP, +10 86 (6.21) (6.22) (6.23) (6.24) In our particular case, since the pressure decays to 1.37 MPa in place of ambient pressure, hence we will consider that value in place of atmospheric pressure. Hence, we can easily calculate the values of all the parameters and can obtain the equations as follows; 41.031 1 t _ 1+ll.845 —’ 0-0153 p(t) =1.37 + 3083(1— 0 0’1 58}: [ [001658)] (6.25) Here, t is in s and p(t) is in MPa. The corresponding curve is shown in figure 6.7 Later on Brode proposed another 5 parameters model in 1956, which is as follows: p(t): p. {(2. —p.)[1-;‘.][ae‘“/‘5 +(1—a)e”"/’3]] (6.261 d where, a: APS +0.5 P0 APs P0 [H016 A135] P0 N’s P0 fl = 70+10 (6.27) a: 1+ Hence, we will get the equation as follows, —23.003t —295.037t p(t)=l.37+30.83 1— ’ 0.196e O~0158+0804e 4.0158 (6.28) 0.0158 87 u (h .1 -°— Experimental + Brode(1 956) + F n + Friedlander + Ethridge-Zparameters + Brode(1 955) —-— Emndngammetars 30- 25 Pressure(MPa) N O .6 0| tameims) Figure 6.7 Pressure time profiles given by different models If we closely examine the pressure profiles in the above figure then we can conclude that the pressure time history given by Ethridge has the best fit for the experimental results. Hence we decided to use this pressure-time profile in our analysis. 6.2. 2. 2. Pressure-distance profile The shock or blast wave is produced when the atmosphere surrounding the explosion is forcibly pushed back by the hot gases produced from the explosion source. The front of the wave, called the shock front, has an overpressure much greater then that in the region behind it. This peak overpressure decreases as the shock is propagated outward. For instance, in figure 6.8, it is shown that an explosive was detonated in air at a height of c from the laminated composite plate. For any particular point (x, y) on the composite, the pressure profile p(r, t) is going to be a function of distance r between the explosive and the point (x, y) on the plate and time. Many authors have done analytical and 88 experimental work to derive the relationships between various blast loading parameters and the distance. Saleh and Adili [12] in their work presented the following relationships to determine the detonation pressure. a Figure 6.8 Schematic diagram of a laminate subjected to blast loading The blast pressure at the point of denotation, p (t), is given by equation (6.21) where the Overpressure AP, as a function of detonation velocity can be written as AP, = 4.18x10‘77v2(1 + 0.87). (6.29) In the above equation, AP, is in the unit of kbars, 7 is the specific gravity of the explosive, and v is the denotation velocity in ft/s. On the surface of the laminate, the distributive pressure loading is controlled by the traveling shock wave front. The pressure can be written as 0, t < L 190. t) = ,‘f (6.30) t e_r, t 2 — p( ) v The distance from the location of the blast to the laminate is determined as, if the blast at the center of the laminate, 89 r=\/(x—-‘i)2 +(y—2)2+cz (6.31) 2 2 Where a and b are the length and width of the rectangular plate shown in figure 6.8, while c is the perpendicular distance between the explosive and the plate. Graham and Kinney [11] in their book modeled the experimental data given in the figure 6.6 for various parameters of pressure profiles and gave following models for chemical explosives. For the peak overpressure, the following relationship explains as to how it varies with sealed distance, Z 808[1+(—‘—Z—)2] (632) p 4.5 5 _. p0 {[1+{O§48)][1+[5€3ll[1+lé]l} Here, the proximity factor, Z is the actual distance scaled to an energy release of 1 kg of TNT in the standard atmosphere. Similarly, the travel time, ta, required for the shock front to travel out to various distances is given by, - -l/2 R 6.33 to = —1- 16p dr ( ) do ’6' 1+ S _ 7P0 3 Where, a0 is the speed of sound in the undisturbed atmosphere The duration of the blast wave is one aspect of its ability to cause damage, for the damage inflicted depends in part on how long the damaging forces are applied. Because the positive pressure phase of the blast wave is the more damaging one, the positive phase duration can be taken as an index of the time duration of the entire blast wave system. . . . + . . This time duration, td , lS grven as, 90 10 980 1+(L) t + 0.54 d _ — , (6.34) l 3 l, 2 W (213 (216 1le 1+ —— 1+ — 1+ — 0.02 0.74 6.9 Where, the duration for positive overpressure is in ms. In a similar fashion, positive impulse per unit area of the blast wave was given by, 0.006721+ Z 0.23 4 I / A = ‘/ ( / ) (6.35) 223([1 +(Z/1.55)3 Here Z is the actual distance scaled to an energy release of 1 kg of TNT in the standard atmosphere. Hence using all these parameters, one can easily derive the whole pressure profiles. This pressure profile can then be used for the purpose of simulations. In our particular experiment, due to the relative small dimensions of the loading area when compared to blast wave front it was assumed that the pressure is constant all over the loading area and hence the pressure profile mentioned in figure 6.4 and later on modeled in figure 6.7 can be used safely for the purpose of simulations at all the points in the loading area. In other words, the pressure profile in our case can be written as, r>R 0, P(r,t) ={p(t). r sR (6.36) Here, p(t) is given by equation (6.21). R is the diameter of the shock tube nozzle and is equal to 0.5 inch. The distance of any point on the laminate from the center is determined as r and is given by, r =\/(x--;-)2 +(y-g)2 (6.37) Here a is defined as the width of the entire plate 91 6.3 Finite Element Analysis The finite element model for the plate is shown in the figure 6.9 below. The square plate shown below is having the width equal to that of the diameter of the circular plate used in the experiment. Simulation of the circular plate with rectangular elements was achieved by selecting a rectangular plate of the same dimension as the diameter of the circular plate, which consisted of an assembly of three-dimensional quasi rectangular four-noded, hundred elements. The degrees of freedom per node are dependent on the number of layers involved in the composite material. For four layers, each node had nineteen degrees of freedom. All edges of the plate are modeled by using the clamped boundary conditions. Moreover, in order to simulate the response of the plate more precisely, all the nodes between the four lines shown below and the vertices were also fully clamped. These four lines were drawn on the basis that they act as tangent to the circular plate used in our experiments. Since the shock tube nozzle diameter is only 0.5”, hence the pressure came only on that area of the plate at the center. In other words, a circular area with diameter 0.5” at the center of the plate was loaded with blast pressure as a function of time. For the linear transient analysis of the plate, time integration was done using the Newmark method. Time increment was taken to 0.04ms. The shear slip constants, D5, and Dsy were set to zero. However, the FORTRAN program had the capability to account for delamination. Whenever, the interfacial shear stress, In or ryz, for any node would exceed the interfacial shear strength, the programme would automatically increment the shear slip constants, D5,, or D5,, from 0 to 0.1 in order to account for delamination and reiterate all the solutions for displacements and stresses. 92 11 121 10 100 l 7 120 9 9 / \ 119 V \ 8 118 7 7 117 6 6 116 5 5 115 4 4 \ I 114 3 \ 113 2 2 112 1 11 \ 91 l ]_1_]_> X Fig 6.9 Finite Element discretization of the plate 6.4 Results and Discussion In this study, the air blast loading is obtained using the shock tube. In order to find out the occurrence of failure or delamination, different failure criteria and delamination criteria were used. In order to use these failure criteria, knowledge regarding maximum stress. values is essential. 6.4.1 ln-plane failure Figure 6.10, 6.11 and 6.12 show the variation of the inplane normal stress and shear stress across the plane of the quarter part of the composite plate. The stress was measured at the interface between 00 and 900 layer and at time 0.04 ms and 0.24 ms after the blast loading impact on the plate. We can clearly see that the transverse normal stresses 0,, and c, are having maximum values at node 61, which is the center of the plate, irrespective of any time. In the case of shear stress, oxy, the stress is not having maximum value at the 93 center of the plate, but it is at node 26. The coordinates of that node are (2a/10, 3a/10), where a refers to the width of the square plate. At these node points, through-thickness distribution of stresses were plotted in order to determine the z coordinate at which they have maximum values. Selection of time t equal to 0.04ms, 0.24ms were randomly chosen in order to study the in-plane and through- thickness distribution of stresses. 600 I I I I 1 +StressX-t=0.24ms +StressX-t=0.04ms 400 P " ' - 200— ,I . , - .’ _ I II \ . ‘ 2'“... . .. 0', \\/ " \ -200 e " . .\ _ -400- - -800 1 l _ I I I 1 800O 10 20 30 4O 50 60 70 Node Number Figure 6.10 Variation of transverse normal stress at the interface across the plane of the quarter [0/90/90/0] composite plate 94 100 50- -50 .— -100 I l -150 -200— +StressY-t=0.24ms + StressY-t=0. 04ms l l J l I -2500 10 20 30 Node Number l 40 50 60 70 Figure 6.11 Variation of transverse normal stress at the interface across the plane of the quarter [0/90/90/0] composite plate 45 40 35 3O 25 20 15 10 -5 I r A +StressXY-t=0.24 ms +StressXY-b004ms O 10 20 l l 30 4O Node Number 50 6O 70 Figure 6.12 Variation of transverse shear stress at the interface across the plane of the quarter [0/90/90/0] composite plate 95 On observing the figures 6.13, 6.14 and 6.15, we can clearly notice that the values of stresses linearly increase or decrease along the z coordinate of each ply. Moreover, for each ply they possess highest magnitude on the surface of the ply. This is true at any instant of time. We can also see that the stresses are symmetric across the mid-plane of the composite. Hence in order to analyze the in-plane failure of the composite laminate we will be only focusing on the ply 1 and ply 2, lying above the mid-plane. Based on these observations table 6.1 gives us the value of maximum in-plane stresses for ply l and ply 2 at each instant of time. 1.5 1 I I f I +StressX-t=0.04 ms + StressX-t=0. 24ms 1 — A 0.5 ~ — N 0 -' -1 -o.5 — — -1 - A _1 l l l l l .1100 -1 000 -500 500 1000 1 500 0x (a/ 2, a/ 2, 2) Figure 6.13 Variation of transverse normal stress at the center across the height of the [0/90/90/0] composite plate 96 1.5 l I I 1 I T 1 —e— StressY—t=0. 04ms + StressY-t=0. 24ms 1 I- ‘i 0.5 - _ N 0 — - -0,5 »- _ -1 h —-t _1 I l m 1 I 1 1 :300 -600 4.00 -200 o 400 600 800 cry (af 2, a/ 2, 270 Figure 6.14 Variation of transverse normal stress at the center across the height of the [0/90/90/0] composite plate 1 .5 1 1 1 1 1 1 I r I —-e— StressXY-t=0. 04ms -6- StressXY-t=0. 24ms 1 — _ 0.5 ~ ‘ _. _ N 0 ~— —~ -o_5 _ _. -1 ._ ‘ A. — _1 1 1 l 1 1 1 £300 -90 -60 4o -20 60 so 1 00 ‘ ' 4‘0 030120710, 3afi0, 2) Figure 6.15 Variation of transverse normal stress at the center across the height of the [0/90/90/0] composite plate We can observe in table 6.1 that the magnitude of stresses increases as time progresses. This is expected because the blast loading pressure applied on the plate increases from time t = 1 second to time t = 0.233 second and after that it starts decreasing exponentially, 97 as shown in figure 6.4. Based on the maximum stress values of table 6.1, failures criterions were applied in order to analyze the in-plane failure of composite laminate. Table 6.1 Maximum in plane stress values for plyl and ply2. time (MS) 0'1 0'2 012 MPa MPa MPa Ply 1 Ply2 Ply 1 Ply2 Ply1 Ply2 0 815.51724 512.61379 276.71034 117.98621 54.348276 26.79931 0.04 898.89655 565.02759 305.0069 1 30.04828 59.905517 29.54 0.08 982.27586 617.44828 333.30345 142.11034 65.463448 32.28069 0.12 982 .27586 617.44828 333.30345 142.11034 65.463448 32.28069 0.16 1149.1034 722.27586 389.89655 166.24138 76.57931 37.761379 0.2 1232,4828 774.68966 418.1931 1 78.30345 82.137931 40.501379 0.24 1257.7931 790.62069 426.7931 181 .97241 83.827586 41.334483 The following failure criterions were used in this report [40] 1. Norris Failure Criteria Norris postulated [41] that failure would occur under plane stress following equations is satisfied (3.11 2 g 11 -932 XY 2 +[q—2] 21, S 2 2 [fl)210r(-Ol) 21 X Y 2. Tsai Hill failure Criteria if any one of the Tsai Hill [42] gave a similar failure criterion, which can be stated as follows (111 2 2:; Y) 3. Fischer Theory 00' 2 [12] 21 S Another orthotropic strength criterion given by Fischer 9 [43] is as follows 98 K: (31162—123 Ell (1+V21)+E22 (1+V12) 0'1 02 XY 2712,13,, (1+ v2, )(1+ 14,) where 01 and 02 represents the in-plane principle stresses. 2",, represents the in-plane principle shear stress. X and Y are either tensile or compressive strengths consistent with S 2 )21 the sign of GI and 02, S is the interlaminar shear strength. Based on table 6.1 and equations (6.38), (6.39) and (6.40), failure criterions were applied on a ply by ply basis. A laminate is assumed to fail analytically if the strength criterion of any one of its laminae is reached. In reality, load distribution usually occurs with in a laminate upon actual failure of an individual ply and hence failure of an individual ply need not necessarily cascade into total fracture of the structure [40]. However, such kind of analysis, often known as the progressive damage analysis is out of the scope of this study. Table 6.2 summarizes the result achieved by applying the failure criterions. The values shown are the calculated left hand side values obtained from equations (6.3 8), (6.39) and (6.40), Table 6.2 Result from various failure criterions. time (ms) Norris Hill Fischer ply 1 ply 2 ply 1 ply 2 ply 1 ply 2 0 74.1363 13.0223 80.8016 14.8073 72.4775 12.5780 0.04 89.8803 15.8057 97.9658 17.9716 87.868 15.2667 0.08 107.3390 18.8598 116.7818I 21.4434 104.7392 18.2168 0.12 107.3390 18.8598 116.7818l 21.4434 104.7392 18.2168 0.16 146.9540 25.7716 160.1794 29.3059 143.6625 24.8920 0.2 168.8395 29.6286 184 33.6965 165.0663 28.6162 0.24 176.1484 0.9792 192.0137 35.2244 172.1999 29.9227 99 Table 6.2 very clearly shows us that the composite laminated plate is supposed to fail immediately as soon as it comes in contact with the blast loading. However, the experimental result shows us that this was not the case. The plate did not undergo any failure with the blast loading shown in figure 6.4. Only when the maximum pressure of the blast wave reached 36.5 MPa did the plate failed as shown by the following experimental results 16.8 MPa 26.6 MP8 29.4 MPa 32.2 MPa 34.3 MPa “ ‘ .5“? 7" 1”“ FM“; 36.5MPa Figure 6.16 Experimental Results In order to explain this discrepancy, it was suggested that the strain displacement relations employed in Lee’s Laminate theory as given by equations (4.3) — (4.8) are linear in nature and hence they may over-predict the stresses. This claim was supported by the fact that the through thickness distribution of stresses given by figures (6.13), (6.14) and (6.15) are linear in nature, which indicates the absence of geometric non-linearity in the laminate theory. Moreover, in order to study this more deeply, strain analysis was also performed as follows. 100 6.4.2 Strain Analysis Similar to stress analysis, the variation of in-plane strains are shown in figures 6.17 and 6.18. These figures show us variation of the in-plane normal stress and shear stress across the plane of the quarter part of the composite plate. The strain was measured at the interface between 00 and 900 layer and at time 0.08 ms after the blast loading impact on the plate. We can clearly see that the transverse normal stresses 8. and 8y are having maximum values at node 61, which is the center of the plate, irrespective of any time. In the case of shear stress, 8,3,, the stress is not having maximum value at the center of the plate, but it is at node 38. The coordinates of that node are (3a/10, 4a/10), where a refers to the width of the square plate 0 (16 .v‘ I I I ‘l l ' +Stralnx +StrainY 0.015 — ~0.01 -0.015 r -0.025 — - £33 1 1 4L l I J_ 0 1 o 20 4o 50 60 7o Node Number Figure 6.17 Variation of in-plane normal strain at the interface across the plane of the quarter [0/90/90/0] composite plate 101 0.016 r r I I I +StrainXY 0.014 '- ‘ " , -1 I 0.012 — 4 l 0.01 — _ xy 0.008 — A 0.006 — _ 0.004 ~ ‘ - 0.002- , ' A 0 l :: I _ l 1 1 ‘ 1 l 0 10 20 30 40 50 80 70 Node Number Figure 6.18 Variation of in-plane shear strain at the interface across the plane of the quarter l0/90/90/0] composite plate At these node points, through-thickness distribution of strains were plotted in order to determine the z coordinate at which they have maximum values. Selection of time t equal to 0.04ms, 0.08ms and 0.24ms were randomly chosen in order to study the in-plane and through-thickness distribution of strains. 102 1.5 I I if I I +Strainx - time =0.24 ms +Strainx - t = 0.04 ms .. ’ 1 — _ 0.5 — ,_ " _ N 0 - - -o.5 - — -1 _ , ‘ _. _1 1 l I 1 I 3.03 -0.02 -o.o1 o 8 0.01 0.02 0.03 x Figure 6.19 Variation of in-plane normal strain at the center across the height of the [0/90/90/0] composite plate 1-5 I l l l I l 1 +StrainY - time =0.24 ms +StrainY - t = 0.04 ms , 1 — _ 0.5 - . " _, N 0 - a -o.5 — — -1 I. ,l " _ _1 I I 1 1 1 I I 45.04 -o.03 -o.02 -0.01 0 8 0.01 0.02 0.03 0.04 Figure 6.20 Variation of in-plane normal strain at the center across the height of the [0/90/90/0] composite plate 103 1.5 1 F i I I j —:—Strainle - time =$.24 ms +StrainXY- t = 0.04 ms -0_5 — _ _1 5 l l l l l l l l 1 -0.025 -0.02 -0.01 5 -0.01 —0.005 O 0.005 0.01 0.01 5 0.02 0.025 gxy Figure 6.21 Variation of in-plane shear strain across the height of the [0/90/90/0] composite plate On observing the figures 6.19, 6.20 and 6.21, we can clearly notice that the values of strains linearly increase or decrease along the z coordinate of each ply. Moreover, for each ply they possess highest magnitude on the surface of the ply. This is true at any instant of time. We can also see that the strains are symmetric across the mid-plane of the composite. Hence in order to analyze the in-plane failure of the composite laminate we will be only focusing on the ply 1 and ply 2, lying above the mid-plane. Based on these observations table 6.3 gives us the value of maximum in-plane strains for ply 1 and ply 2 at each instant of time. We can observe that the magnitude of strains increases as time progresses. This is expected because the blast loading pressure applied on the plate increases from time t = 1 second to time t = 0.233 second and after that it starts decreasing exponentially, as shown in figure 6.4. Also, we can see that the as soon as the blast loading come in touch with the plate, the maximum strain values at the center of the plate goes up to 1.85 and 2.43%. These values then increases to 2.85 and 3.74%. Figure 6.16 shows us that during the loading process; internal damage in the specimen 104 should cause the material to behave non-linearly (as the tested specimen showed permanent deformation). Hence the calculated stresses and strains are all larger. 6.3 Maximum In-plane strains for plyl and ply 2 3:?) 31 3 2 312 Ply1 Ply2 Ply1 Ply2 Ply1 Ply2 0.0000 0.0185 0.0090 0.0243 0.0121 0.0139 0.0069 0.0400 0.0204 0.0133 0.0268 0.0099 0.0153 0.0076 0.0800 0.0223 0.0146 0.0292 0.0108 0.0168 0.0083 0.1200 0.0223 0.0146 0.0292 0.0108 0.0168 0.0083 0.1600 0.0260 0.0170 0.0342 0.0126 0.0196 0.0097 0.2000 0.0279 0. 0183 0.0367 0.0135 0.0210 0.0822 0.2400 0.0285 0.0186 0.0374 0.0138 0.0215 0.0106 Another possible reason which could have caused the stress and strain values to increase was the coupling between the in-plane stretching and the out—of-plane deflection. Since the composite was fixed around all the edges, the in-plane constraint prevented a large transverse deflection from happening. In the meanwhile, it also caused the in-plane stresses to increase. This was not found in the experiments because in reality we cannot have fully clamped situations and hence when the blast loading comes on the plate, the plate can stretch in the x-y plane. 6.4.3 Delamination Analysis Composite laminates are heterogeneous and anisotropic. Their properties vary with in a lamina and also from one lamina to another through the laminate thickness. Because of the mismatch in material properties, the non-uniform stress distribution through the laminate thickness may cause delamination in the composite laminate. Once delamination takes place with in a composite laminate, it is very easy for the impact loading to destroy the composite by cracking the matrix and moving aside or destroying the fibers in 105 individual laminae. Hence, analysis of delamination is very critical for composite analysis. Figure 6.22 and 6.23 shows us the variations of shear stress across the plane of the quarter composite laminate. The stress was measured at the interface between 00 and 90° layer and at time 0.04 ms and 0.24 ms after the blast loading impact on the plate. We can observe that the shear stress, Tn, possess maximum values at node 50, while the shear stress, 1),, , is having maximum value at node 4. Hence at these nodes, we can find the maximum value of shear stresses at all the interfaces at each instant of time. . +StressXZ-t=0.24ms + StressXZ—t=0. 04ms —-1 in C3 an I l l l 1 4L t 1 0 2o 30 40 Node Number Figure 6.22 Variation of transverse normal stress at the interface across the plane of the quarter [0/90/90/0] composite plate 106 30 1 I l i T —a— StressYZ-t=0. 24ms + StressYZ—t=0. 04ms 20 10 -10 -20 -30 -40 1O 20 I 30 1 40 Node Number 50 60 70 Figure 6.23 Variation of transverse normal stress at the interface across the plane of the quarter [0/90/90/0] composite plate Table 6.4 below gives us the value of the shear stresses at different instants of time across all the interfaces. We can see that just like table 6.1 the magnitude of shear stresses at each interface goes on increasing as the magnitude of pressure loading keeps on increasing. Table 6.4 Values of shear stresses across each interface at different instants of time time (ms) T)" TH MPa Interface Interface Interface Interface Interface Interface 1 2 3 1 2 3 0.0000 30.1828 33.3924 30.1841 -20.0214 -30.6234 -20.0214 0.0400 33.2690 36.8069 33.2703 -22.0683 -33.7545 -22.0683 0.0800 36.3552 40.2214 36.3572 -24.1 159 -36.8862 -24.1 159 0.1200 36.3572 40.2214 36.3572 24.1159 -36.8862 -24.1159 0.1600 42.5283 47.0503 42.5283 -28.2103 43.1490 -28.2103 0.2000 45.6145 50.4648 45.6145 -30.2579 -46.2807 -30.2572 0.2400 46.5524 51 .5028 46.5524 -30.8800 47.2324 -30.8800 107 The values obtained from table 6.4 were used in order to find out if delamination took place in the composite plate. The delamination failure criterion proposed by Sun and Zhou [44] was used in this particular analysis. The criterion states that in order for delamination to take place across any interface, (6.41) where S is the interlaminar shear strength, while TR and 1),: are the maximum values of shear stresses on a particular interface. Using this criterion following results were generated as shown in table 6.5. Table 6.5 Results from delamination criterion Time (ms) lnterface1 Interface2 Interface 3 0 0.253 0.396 0.253 0.04 0.307 0.481 0.307 0.08 0.367 0.575 0.367 0.12 0.367 0.575 0.367 0.16 0.502 0.786 0.502 0.2 0.578 0.904 0.578 0.24 0.602 0.942 0.602 Hence we find that, with the blast loading profile given in figure 6.4, delamination cannot take place in the composite plate. We can see in the table 6.4 for time t =0.24ms and for ply 2 or the 90° ply, the criterion indicates that at time t = 0.24 ms, the value of the left hand side of equation (6.41) is equal to 0.94 which is very near to 1. In other words, if we increase the loading, then delamination will take place. 108 7. CONCLUSIONS AND RECOMMENDATIONS 7.1 Conclusions An elasticity analysis was performed to study the response of the interface of [0/90/0] composite plate under different bonding conditions. The well known Pagano’s cylindrical bending problem was extended for an imperfectly bonded composite plate. The [0/90/0] composite plate was simply supported at the ends and was subjected to a sinusoidal tensile loading at the top surface. In order to represent different bonding conditions at the interface, two different approaches were taken, the embedded layer approach and the slip approach. In the embedded layer approach, an additional matrix layer of very small thickness was introduced between the laminae and the different bonding conditions were simulated by varying the Young’s Modulus and Shear Modulus of the embedded layers. In the slip approach, a linear shear slip theory and a linear normal separation theory were used to account for the variation in interfacial displacement For the finite element simulation of the interfacial damage, again two approaches were presented, the embedded layer approach and the shear slip approach. In the embedded layer approach, the different bonding conditions were simulated by varying the Young’s Modulus and Shear Modulus of the embedded layers. In the shear slip approach, only the linear shear slip theory was used for the interfacial displacement. The two delamination models given by the finite element simulation were used to analyze the [0/90/0] composite plate and the results were compared with the two approaches used in elasticity analysis. Results matched well for the case of perfectly bonded composite laminated plates. For the case of imperfectly bonded composite, the 109 finite element model qualitatively predicted the presence of delamination, however only shear slip on the interfacial layers could be modeled by it. A brief explanation about the nature of pressure profile and various other parameters involved in blast loading was given. The finite element model was used to investigate the response of [0/90/90/0] composite laminate subjected to blast loading. The results from the simulations were compared with those from the experiments. The simulation results were not in agreement with the experiments. The discrepancy was attributed to the absence of geometric non-linearity in the laminate theory and also due to the coupling between the in-plane stretching and the out-of-plane deflection. 7.2 Recommendations The finite element code used was unable to capture the stress variation caused due to the sudden blast loading on the composite plate. By incorporating non-linearity in the strain displacement relations of the plate and including the coupling effects between the in-plane stretching and the out-of-plane deflection, this defect is expected to be resolved. Hence an improved laminate theory including the non-linearity of strain displacements and coupling effects should be investigated. The laminate theory used in the present study incorporates displacement discontinuity across the interface using the shear slip theory. However, it does not include the normal separation theory to account for the normal separation between the interfaces. Hence, we find that for imperfectly bonded composite, the finite elements results do not match well with that of the modified Pagano’s problem results in case of tensile loading. Thus a revised laminate theory incorporating the normal separation theory would be very useful. 110 Appendix A Figures of stresses and displacements in order to determine critical thickness of the embedded layer 0'6 I I fl I T I 1 l 0.4 — ._ _ Jr J" " -~—— ELT-0.1 + ELT-0.07 0.2 — ——a— ELT-0.05 A , + ELT-0.03 II + ELT-0.01 :: _.__ ELT-0.001 n -—+—— ELT-0.0001 0 '- II ..I 1' +710 EL II N :' I .. l 1 l 0:03 -0.6 -O.4 -0.2 _3 01.2 0.4 0.6 0.0 1 a. (a/2 . b/z , 2) FigureAJ- Normalized in-plane normal stress distribution along the height of a [0/0/90/0/0] composite plate for different values of embedded layer thickness (ELT) with perfect bonding conditions 0.6 I fir 1 I I I I —+— ELT-0.1 + ELT-0.07 —e— ELT-0.05 0.4 ~ + ELT-0.03 “ + ELT-0.01 _.... ELT-0.001 —~—- ELT-0.0001 0.2 — -0—no EL _ o _ —: N -0.2 — e -o.4 ~ — -O.6 — e _ 1 L l l l I 0:88 -o.e -o.4 - .2 0.4 0.6 ’ 8:(a/2,b/2,z) FigureA.2- Normalized in-plane normal stress distribution along the height of a [0/0/90/0/0] composite plate for different values of embedded layer thickness (ELT) with perfect bonding conditions 111 0.6 I I I I -~— ELT-0.1 + ELT-0.07 —a— ELT-0.05 O 4 __ + ELT-0.03 . ' + ELT-0.01 —+- ELT-0.001 ..__,_. ELT-0.0001 +110 EL 0.2 — 0 ~ .. N -0.2 e — -o..L . -0.6 — — _0 l l l l 3.9 1 95 2 1 2.15 2 2 W(a/2,2 135/ 2, z) FigureA.3- Normalized transverse displacement variation along the height of a |0/0/90/0/0] composite plate for different values of embedded layer thickness (ELT) with perfect bonding conditions 0.6 I I I I I I I 0.2 — —B— ELT-0.07 —-a— ELT-0.05 —0~— ELT-0.03 + ELT-0.01 -+— ELT-0.001 —‘-— ELT-0.0001 +no EL —OJ.5 01.6 017 018 Gig 1 0'2 (a/2,b/2,z) F igureA.4 Normalized transvserse normal stress distribution along the height of a [0/0/90/0/0] composite plate for different values of embedded layer thickness (ELT) with perfect bonding conditions 1 l l 0.2 0.3 0.4 112 0.6 l I l l e-\. —+— ELT-0.1 I: \\ + ELT-0.07 0 4 _ \xi§.-~\ —a— enoos ' \—\E_<‘ :. 0 ELT-0.03 . , + ELT-0.01 “as; 1., -¢—ELT-0.001 0.2 - A —‘+— ELT-0.0001 ~ -e—no EL 0 - ._ N -0.2 - - -0.4 — ‘ -O.6 ~ — -O 1 I I I 1 47.06 -0.04 -0.02 0.02 0.04 0.06 {910,01 2) FigureA.5 — Normalized in-plane shear stress distribution along the height of a [0/0/90/0/0] composite plate for different values of embedded layer thickness (ELT) with perfect bonding conditions 113 F Appendix B Matlab code for the analysis of imperfectly bonded composite laminate using the Embedded Layer Approach clc clear all hold on N=5; %No. of Layer point=10; %Calculate 10 points per Layer toplayer l; midlayer 3; botlayer N; tt=0; t(l)=0.3; t(2)=0.001; t(3)=0.3; t(4)=0.001; t(5)=0.3; for k = 1:N tt = tt+t(k); end %tt=l; a=tt*4; %tt*(a/h) P=Pi/a; b=1*a; %rectangular proportion Q=Pi/b; Sigma=-1; x=a/2; Y=b/2; C=((p‘2)+(q‘2))‘o.5; YM = input('Enter the value of Youngs Modulus, E = '); MR = input('Enter the value of Shear Modulus, G = '); o=tt/2; S=a/tt; h(1)=o; for k = 1:N h(k+1)=h(k)-t(k); end for k = 1:N 114 E(1,l)=25€6; E(l,2)=186; E(1,3)=186; G(1,1)=.5e6; G(1,2)=.2e6; G(1,3)=.5e6; V(l,1)=.25; v(1,2)=.25; V(1,3)=.25; E(2,l)=YM; E(2,2)=YM; E(2,3)=YM; G(2,1)=MR; G(2,2)=MR; G(2,3)=MR; V12.1)=(YM/(2*MR))-1i v(2,2)=(YM/(2*MR) ) -1; V(2.3)=(YM/(2*MR))-1; E(3,l)=1€6; E(3,2)=25€6; E(3,3)=le6; G(3,1)=.Se6; G(3,2)=.5€6; G(3,3)=.286; v(3,1)=.01; V(3,2)=.25; V(3,3)=.25; E(4,l)=YM; E(4,2)=YM; E(4,3)=YM; G(4,1)=MR; G(4,2)=MR; G(4,3)=MR; v(4,l)=(YM/(2*MR) ) -1; v(4,2)=(YM/(2*MR))-l; V(4,3)=(YM/(2*MR))-l; E(5,l)=25€6; E(5,2)=l€6; E(5,3)=1€6; G(5.l)=.5e6; G(5,2)=.286; Gl5,3)=.5€6; V(5,l)=.25; v(5,2)=.25; v(5,3)=.25; Layer=k; v(k,4) = v(k,l)*E(k,2)/E(k,1); %V21 v(k,5) = v(k,2)*E(k,3)/E(k,2); %v32 115 v(k,6) deltalk) = (1-(v(k,1)*v(k,4))-(v(k,2)*v(k,5))-(v(k,3)*v(k,6))... = V(k,3)*E(k,3)/E(k,l); %v3l -(2*v(k,1)*v(k,2)*v(k,6)))/(E(k,1)*E(k,2)*E(k,3)); C(k,1) C(k,2) c(k.3) C(k,4) C(k,5) C(k,6) C(k,7) C(k,8) C(k,9) A(k) B(k) cc(k) (1-V(k,2)*v(k,5))/(E(k,2)*E(k,3)*delta(k)l; %Cll (l-v(k,3)*v(k,6))/(E(k,3)*E(k,1)*delta(k)); %C22 (1-v(k,1)*v(k,4))/(E(k,1)*E(k,2)*delta(k)); %C33 G(k,2); %C44 G(k,3); %C55 G(k,l); %C66 C(k,3)*c0 CU(k.j) CM(k,j) cosh(m(k,j)*h(k)); %Cj(z) at h/2 116 (v(k,4)+v(k,6)*v(k,2))/(E(k,2)*E(k,3)*delta(k)l; %C12.C2l (v(k,5)+v(k,6)*v(k,l))/(E(k,3)*E(k,l)*delta(k)li %C23,C32 (v(k,6)+v(k,4)*v(k,5))/(E(k,2)*E(k,3)*delta(k)); %Cl3,C3l cosh(m(k,j)*(h(k)+h(k+l))/2); %Cj(z) at midlayer cosh(m(k,j)*h(k+1)); %Cj(z) at ~h/2 SU(k,j) = sinh(m(k,j)*h(k)); %Sj(z) at h/2 SM(k,j) = sinh(m(k,j)*(h(k)+h(k+l))/2); %Sj(z) at midlayer SB(k,j) = sinh(m(k,j)*h(k+1)); %Cj(z) at -h/2 alfa(k,j) = 1; for g = 1:point+1 hhlk,g) = h(k)-(g-1)*step; CZ(k,g,j) = cosh(m(k,j)*hh(k,g)); SZ(k.g,j) = sinh(m(k,j)*hh(k,g)). end elseif (gamma(k,j)+(B(k)/(3*A(k))))<0 CU(k,j) = cos(m(k,j)*h(k)); %Cj(z) at h/2 CM(k,j) = cos(m(k,j)*(h(k)+h(k+l))/2); %Cj(z) at midlayer CB(k,j) = cos(m(k,j)*h(k+1)); %Cj(z) at —h/2 SU(k,j) = sin(m(k,j)*h(k)); %Sj(z) at h/2 SM(k,j) = sin(m(k,j)*(h(k)+h(k+1))/2); %Sj(z) at midlayer SB(k,j) = sin(m(k,j)*h(k+l)); %Cj(z) at ~h/2 alfa(k,j) = -1; for g = lzpoint+1 end else hh(k,g) = h(k)-(g-1)*step; CZ(k,g,j) = cos(m(k.j)*hh(k.g)); SZ(k,g,j) = sin(m(k.j)*hh(k.g)); ‘Error' end J(k,j) le,j) C(k,3)*C(k,7)*m(k,j)*4 +alfa(k,j)*(m(k,j)‘2)*(-(p02)*(C(k,7)*C(k,8)... +C(k,3)*C(k,9))+q*2*(C(k,5)*2-CIk,2)*CIk,3)... +2*C(k,5)*C(k,7)))+(C(k,9)*p‘2 ... +cIk,2)*q*2)*(C(k,8)*p‘2+C(k,7)*q‘2); =p*q/J(k.j)*(alfa(k.j)*(m(k.j)‘2)*(c(k,3)*(C(k,4)+c(k,9))... ~(Clk.5)+C(k.7))*(C(k,6)+C(k,8)))... —(c1k,4)+c*L(k.j); vzz(k,g.j)= SZ(k.g.j)*L(k,j); Wiz(k,g.j)= alfa(k,j)*SZ(k,g,j)*R(k,j); WZZ(k,g,j)= CZ(k,g,j)*R(k,j); SigmaZlZ(k,g,j)=CZ(k,g,j)*M(k,j,3); SigmaZ2Z(k,g,j)=SZ(k,g,j)*M(k,j,3); TauXZ1Z(k,g,j)=C(k,8)*(alfa(k,j)*SZ(k,g,j)*(m(k,j)+p*R(k,j))); Tauxzzz(k,g,j)=C(k,8)*(CZ(k,g,j)*(m(k,j)+p*R(k.j))); TauYZiZIk,g,j)=C(k,7)*(a1fa(k,j)*SZ(k,g,j)*(mIk,j)*L(k,j)+q*R(k.j))); Tauvz2ZIk,g,j)=C(k,7)*(czIk,g,j)*(m(k,j)*L(k.j)+q*R(k.j))); TauXYiZIk,g,j)=C(k,9)*(q+(p*L(k,j)))*CZ(k.g.j); 118 TauXY2Z(k,g,j)=C(k,9)*(q+(p*L(k,j)))*SZ(k,g,j); SigmaX1Z(k,g,j)=M(k,j,1)*CZ(k,g,j); sigmaxzzIk,g,j)=M = TauXZ(k.9)/(Sigma*S); %plot(TauXZ(k,g),hh(k,g),'*b'); end end %TauYZ at (a/2,0,0)' TauYZmid=0; n =mid1ayer—1; for j = 1:3 TauYYZZmid(j) = C(2,7)*(m(2,j)*L(2,j)+q*R(2,j))*(alfa(2,j)... *(FG(n*6+j)*SM(2,j)+FG(n*6+3+j)*CM(2,j))); TauYZmid = TauYZmid + TauYYZZmid(j); end TauYZmid=TauYZmid/(Sigma*S); 124 %TauYZ Plot for k = 1:N n = k-l; for g = 1:point+1 TauYZ(k,g)=0; for j = 1:3 TauYYZZ(k,g,j) = (FG(n*6+j)*TauYZlZ(k,g,j)+FG(n*6+3+j)*TauYZ2Z(k,g,j)); TauYZ(k,g) = TauYZ(k,g) + TauYYZZ(k,g,j); end TauYZ(k,g) = TauYZ(k,g)/(Sigma*S); %plot(TauYZ(k,g),hh(k,g),'*b'); end end %TauXY at the top surface' TauXYtop=0; n =toplayer-1; for j = 1:3 TauXXYYtop(j) = C(l,9)*(q+p*L(1,j))*(FG(n*6+j)*CU(l,j)... +FG(n*6+3+j)*SU(l,j)); TauXYtop = TauXYtop + TauXXYYtop(j); end TauXYtop=Tauxvtop/(Sigma*s*2); %TauXY at the bottom surface' TauXYbot=O; n =botlayer-1; for j = 1:3 TauXXYYbot(j) = C(N,9)*((q+p*L(N,j))*(FG(n*6+j)*CB(N,j)... +FG(n*6+3+j)*SB(N.j))); TauXYbot = TauXYbot + TauXXYYbot(j); end . TauXYbot=TauXYbot/(Sigma*S‘2); %TauXY Plot for k = 1:N n = k-l; for g = 1:point+1 TauXY(k,g)=0; for j = 1:3 TauXXYY(k,g,j) = (FG(n*6+j)*TauXY1Z(k,g,j)... +FG(n*6+3+j)*TauXY2Z(k,g,j)); TauXY(k,g) = TauXY(k,g) + TauXXYY(k,g.j); end TauXY(k,g) = TauXY(k,g)/(sigma*s*2); %plot(TauXY(k,g),hh(k,g),‘*b'); end end %X Displacement (U) Plot for k = 1:N n = k-l; for g = 1:point+1 DiSpX(k,g)=O; for j = 1:3 DistX(k,g,j) = (FG(n*6+j)*CZ(k,g,j)... 125 +FG(n*6+3+j)*SZ(k,g,j)); Dist(k,g) = Dist(k,g) + Dispxx(k,g,j); end Dist(k,g) = E(1,2)*Dispx(k,g)/(Sigma*tt*S‘3); %plot(Dist(k,g),hh(k,g),'*b'); end end % Y Displacement (V) Plot for k = 1:N n = k-l; for g = 1:point+1 DispY(k,g)=O; for j = 1:3 DispYY(k,g,j) = (FG(n*6+j)*V12(k,g,j)... +FG(n*6+3+j)*V2Z(k,g,j)); DispY(k,g) = DispY(k,g) + DispYY(k,g,j); end DispY(k,g) = E(1,2)*DispY(k,g)/(Sigma*tt*s‘3); %plot(DispY(k,g),hh(k,g),'*b‘); end end %Z Displacement at middle height(W) Plot k = midlayer; n = k-l; Dispz=0; for j = 1:3 DispZZ(k,g,j) = R(k,j)*(FG(n*6+3+j)*CM(k,j)... +alfa(k,j)*FG(n*6+j)*SM(k,j)); Dispz = DispZ + DispZZ(j); end DiSpZ; DispZ = 100*E(1,2)*Dispz/(Sigma*tt*s*4); % plot(S,DispZ,'m.'); %X Displacement (W) Plot for k = 1:N n = k-l; for g = 1:point+1 DiSpZ(k,g)=0; for j = 1:3 DispZZ(k,g,j) = (FG(n*6+j)*W1Z(k.g,j)--- +FG(n*6+3+j)*W2Z(k,g,j)); DispZ(k,g) = Dispz(k,g) + DispZZ(k,g,j); end DispZ(k,g) = 100*E(1,2)*DispZ(k,g)/(Sigma*tt*s‘4); plot(DispZ(k,g),hh(k,g),'*b'); end end 126 Appendix C Figures of stresses and displacements obtained using the Embedded Layer approach in the Pagano’s solutions 0.5 I I I I I j I . + E=2.4E6 s. G=1 56 0.4 _ 1 +E=2.4E4 a. 01-154 _ +E=2.4E2 a. 0:152 . +E=2.4E0 a. G=1 150 0.3 — '. —e—~o Embedded Layers ‘ 0.2 — _ 0.1 — 1, _ v ‘7 N 0 e v d t P 1 r '0.‘ h : —4 V -0. _, -o, _ -04 - _ _ 1 1 l 1 l l 1 0.50 4 10 12 14 16 e W(0, b/ 2 , z) Figure C.l-Normalized transverse displacement along the height of laminate for different bonding conditions obtained by varying the material properties of the embedded layers of a [0/1/90/1/0] composite plate. [I-isotropic] 0.5 I I I 1 I j T I I 0.4 — 4 ' .. ’ ~ 0.3 ~ , “' - . .2" o 2 _ . -¢— E=2.4E6 & G=1Ee fl ' . ' _ ' +E-Z.4E4 a. on E4 . +E=2.4E2 a. GM 52 J 0-1 “ +E=2.4Eo a G=1 so -o—NoEnmukhqumm N 0 — — -o.1 — j s e , ‘ -o.2 — t, , A .0.3 >— y . .4 / -Q4~ , _ - 1 I 1 1 l 1 1 1 1 0.65 _4 3 2 1 2 3 4 5 Z(a/2.b/2.z) Figure C.2-Normalized in-plane normal stress along the height of laminate for different bonding conditions obtained by varying the material properties of the embedded layers of a [0/I/90/I/0] composite plate. [I-isotropic] 127 0.5 1 1 fl 1 1 _ + E=2.4E6 81 G=1E6 0.4 — " —a— E=2.4E4 & G=1 E4 4 + E=2.4EZ & G=1 E2 0.3 — +E=2.4Eo s. G=1EO 1 -o-No Embedded Layers 0.2 - —- 0.1 - _ N 0 — a -0.1 — - -0.2 - — -0.3 r- - -0.4 — d l _ 1 1 1 1 1 03 .5 -1 -O.5 0.5 1 1 .5 0 O'V(a/2,b/2,z) Figure C.3-Normalized in-plane normal stress along the height of laminate for different bonding conditions obtained by varying the material properties of the embedded layers of a [0/1/90/1/0] composite plate. [I-isotropic] 0.5 0.4 0.3 0.2 + E=2.4E6 & G=1 E6 0“ +E=2.4E4 & G=1 E4 ‘ +E=2.4E2 a. G=1E2 ... +E=2.4E0 & G=1E0 - +No Embedded Layers -o.1 _ -o. 4 -o. _ -o.4 - _0_5 1 1 1 1 1 1 o 0.1 0.2 0.3 o 4 0.7 0.8 0.9 1 l . _ 0.5 0.6 a, (a/2,b/2 , z) Figure C.4-Normalized transverse normal stress along the height of laminate for different bonding conditions obtained by varying the material properties of the embedded layers of a [0/1/90/1/0] composite plate. [I-isotropic] 128 0.5 1 1 1 r 1 1 1 1 I ' ‘ . ‘ . + E=2.4E6 & G=1 E6 0.4 r .- .. . +E=2.4E4 8. G=1 E4 * ‘ + E=2.4EZ & G=1 E2 0.3 L '41- E=2.4E0 & G=1EO _ +No Embedded Layers 0.2 ~ , _ r \ 0.1 — — N O — 4 .01 _ .. -02 — .. -03 _ _ -0.4 - — _0 l l l l l l 1 l l 3.25 -0.2 -0.15 -0.1 -0.05 G— 0.0 0.1 0.15 0.2 0.25 r)? (0,0, 2 Figure C.5-Normalized in-plane shear stress along the height of laminate for different bonding conditions obtained by varying the material properties of the embedded layers of a |0/I/90/I/0] composite plate. [I-isotropic] 129 Appendix D Matlab code for the analysis of imperfectly bonded composite laminate using the Slip Approach clc clear all hold on N=3; %No. of Layer point=10; %Calculate 10 points per Layer toplayer = 1; midlayer = 2; botlayer = N; tt=0; sz = input('Enter the value of shear slip constant in x direction, DSX='); Dsy = input('Enter the value of shear slip constant in x direction, Dsy= '); K = input('Enter the value of shear slip constant in X direction, K ='); t(l)=0.3; t(2)=0.3; t(3)=0.3; % Shear Constants Values for k = 1:N tt = tt+t(k); end a=tt*4; %tt*(a/h) P=Pi/a; b=1*a; %rectangular proportion Q=Pi/b; Sigma=1; x=a/2; Y=b/2; c=1(p‘2)+(q‘2))*o.s; O=tt/2; S=a/tt; h(1)=O; for k = 1:N h(k+l)=h(k)-t(k); end for k = 1:N E(l,l)=25€6; E(1,2)=1€6; E(1,3)=1e6; G(l,1)=.5€6; G(1,2)=.2e6; 130 G(l,3)=.5e6; V(1,l)=.25; V(l,2)=.25; V(l,3)=.25; E(2,1)=1e6; E(2,2)=25e6; E(2,3)=1e6; G(2,1)=.5e6; G(2,2)=.5€6; G(2,3)=.2€6; v(2,1)=.01; V(2,2)=.25; V(2,3)=.25; E(3,1)=25€6; E(3,2)=1e6; E(3,3)=le6; G(3,l)=.5e6; G(3,2)=.2e6; G(3,3)=.Se6; V(3,l)=.25; V(3,2)=.25; V(3,3)=.25; Layer=k; v(k,4) = v(k,1)*E(k,2)/E(k,1); %v21 v(k,5) = v(k,2)*E(k,3)/E(k,2); %v32 v(k,6) = v(k,3)*E(k,3)/E(k,l); %v31 delta(k) = (l-(v(k,l)*v(k,4))-(v(k,2)*v(k,5))-(V(k,3)*v(k,6))... -(2*V(k,1)*V(k.2)*V(k.6)))/(E(k.1)*E(k.2)*E(k,3)); C(k,l) = (1-v(k,2)*v(k,5))/(E(k,2)*E(k,3)*delta(k)); %c11 C(k,2) = (l-v(k,3)*v(k,6))/(E(k,3)*E(k,l)*delta(k)); %c22 C(k,3) = (1-v(k,1)*v(k,4))/(E(k,l)*E(k,2)*delta(k)); %c33 C(k,4) = (v(k,4)+v(k,6)*v(k,2))/(E(k,2)*E(k,3)*delta(k)); %c12,c21 C(k,5) = (v(k,5)+v(k,6)*v(k,1))/(E(k,3)*E(k,1)*delta(k)); %c23,c32 C(k,6) = (v(k,6)+v(k,4)*v(k,5))/(E(k,2)*E(k,3)*delta(k)); %Cl3,C31 C(k,7) = G(k,2); %C44 C(k,8) = G(k,3); %css C(k,9) = G(k,l); %C66 A(k) = C(k,3)*C(k,7)*C(k,8); B = (p‘2)*(C(k,7)*(C(k,1)*C(k,3)-(C(k,6))*2)... +C(k,8)*(C(k,3)*C(k,9)-2*C(k,6)*C(k,7)))... +(q‘2)*(C(k,8)*(C(k,2)*C(k,3)-(c(k,5))*2)... +C(k.7)*(C(k.3)*C(k,9)-2*C(k.5)*c(k,8))); CC(k) = (-p‘4)*(c1k,9)*(C(k,1)*c1k,3)-(C(k,6))*2)... +C(k,8) * (C(k,1)*C(k,7) -2*C(k,6) *C(k,9) ) ) . . . +(p‘2)*(q‘2)*(-C(k,1)*(C(k,2)*C(k,3)-(C(k,5>)‘21... -2*(C(k,4)+C(k,9))*(C(k,6)+C(k,8))*(C(k,5)+C(k,7))... -2*C(k,7)*C(k,8)*C(k,9)+2*C(k,1)*C(k,5)*C(k,7)... +C(k,4)*C(k,3)*(C(k,4)+2*C(k,9))... +C(k,6)*C(k,2)*(C(k,6)+2*C(k,8)))- (q‘4)*(C(k,9)*(C(k,2)*C(k,3)... -(C(k,5))*2)+C(k,7)*(C(k,2)*c(k,8)-2*C(k,5)*c1k,9))); 131 D(k) (p‘6)*C(k.1)*C(k,8)*C(k.9)+(p‘4)*(q‘2)*(C(k.8)*(C(k,1)*C(k.2)... -C(k,4))+C(k,9)*(C(k,1)*C(k,7)—2*C(k,4)*C(k,8))>... +(p‘2)*(q‘4)*(C(k,7)*(C(k.1)*C(k,2)-C(k.4)‘2)... +C(k,9)*(C(k.2)*C(k,8)-2*C(k,4)*C(k,7)))... +(q‘6)*c0 'Error' end step (3*CC(k)*A(k)+B(k)*2)/(-3*A(k)*2); (2*(B(k))‘3+9*A(k)*B(k)*CC(k)+27*D(k)*(A(k))‘2)/(- ((f(k))‘2)/4+((d(k))‘3)/27; (h(k)-h(k+1))/point; % Whether the layer is isotropic or not kind: input('Is this an isotropic layer? Enter 0 or 1(0- isotropic/1-non-isotropic): if kind ==l; '); phi(k) = acos(-f(k)*sqrt(27)/(2*(—d(k))‘1.5)); for j = 1:3 gamma(k,j) = 2*sqrt(-d(k)/3)*cos(1/3*(phi(k)+2*(j-1)*pi)); m(k,j) = sqrt(abs(gamma(k,j)+(B(k)/(3*A(k))))); if (gamma(k,j)+(B(k)/(3*A(k))))>0 CU(k,j) = cosh(m(k,j)*h(k)); %Cj(z) at h/2 CM(k,j) = cosh(m(k,j)*(h(k)+h(k+1))/2); %Cj(z) at midlayer CB(k,j) = cosh(m(k,j)*h(k+1)); %Cj(z) at —h/2 SU(k,j) = sinh(m(k,j)*h(k)); %Sj(z) at h/2 SM(k,j) = sinh(m(k,j)*(h(k)+h(k+1))/2); %Sj(z) at midlayer SB(k,j) = sinh(m(k,j)*h(k+1)); %Cj(z) at —h/2 alfa(k,j) = l; for g = 1:point+1 hh(k,g) = h(k)-(g-1)*step; CZ(k,g,j) = cosh(m(k,j)*hh(k,g)). SZ(k.g,j) = sinh(m(k,j)*hh(k,g)); end elseif (gamma(k,j)+(B(k)/(3*Alk))))<0 CU(k,j) = cos(m(k,j)*h(k)); %Cj(z) at h/2 CM(k,j) = cos(m(k,j)*(h(k)+h(k+1))/2); %Cj(z) at midlayer CB(k,j) = cos(m(k,j)*h(k+1)); %Cj(z) at —h/2 SU(k,j) = sin(m(k,j)*h(k)); %Sj(z) at h/2 SM(k,j) = sin(m(k,j)*(h(k)+h(k+l))/2); %Sj(z) at midlayer SB(k,j) = sin(m(k,j)*h(k+l)); %Cj(z) at -h/2 alfa(k,j) = -1; for g = 1:point+1 hh(k,g) = h(k)-(g-l)*step; CZ(k,g,j) = cos(m(k,j)*hh(k,g)); SZ(k.9,j) = sin(m(k,j)*hh(k,g)); end 132 else 'Error‘ end J(k,j) = C(k,3)*c1k,7)*m(k,j)*4 +alfa(k,j)*(m(k,j)‘2)*(-(p‘2)*(C(k,7)*C(k,8)... +C(k,3)*c1k,9))+q*2*(c*c1k,3)... +2*C(k,5)*C(k,7)))+(C(k,9)*p‘2 ... +C(k,2)*q*2)*(C(k,8)*p*2+C(k,7)*q*2); L(k,j)=p*q/J(k,j)*(alfa(k,j)*(m(k,j)*2)*(C(k,3)*(C(k,4)+C(k,9))... -(C(k.5)+C(k.7))*(C(k.6)+C(k.8)))... -(C(k,4)+C(k,9))*(C(k,8)*p‘2+c(k,7)*q‘2)); R(k,j) =p*m(k,j)/J(k.j)*(alfa(k.j)*(m(k,j)A2)*C(k.7)*(C(k.6)... +c(k,8))-(C(k,6)+C(k,8))*(C(k,9)*p‘2+c(k,2)*q‘2)... +; M*CU(1,j>--. +FG(n*6+3+j)*SU(l,j)); SigmaXtop = SigmaXtop + SigmaXXtop(j); end SigmaXtop=SigmaXtop/(Sigma*S‘2); %'SigmaX at the bottom surface‘ SigmaXbot=0; n =botlayer-1; for j = 1:3 SigmaXXbot(j) = M(3,j,1)*i TauXY(k,g) = TauXY(k,g) + TauXXYY(k,g,j); end 140 T! '- TauXY(k,g) = TauXY(k,g)/(Sigma*s‘2); end end %X Displacement (U) Plot for k = 1:N n = k-l; for g = 1:point+1 Dist(k,g)=0; for j = 1:3 DistX(k,g,j) = (FG(n*6+j)*CZ(k,g,j)... +FG(n*6+3+j)*SZ(k,g,j)); Dist(k,g) = Dist(k,g) + DistX(k,g,j); end Dist(k,g) = E(1,2)*Dist(k,g)/(Sigma*tt*s‘3); end end % Y Displacement (V) Plot for k = 1:N n = k-l; for g = 1:point+1 DispY(k,g):O; for j = 1:3 DispYY(k,g,j) = (FG(n*6+j)*V1Z(k,g.j)... +FG(n*6+3+j)*VZZ(k,g,j)); DispY(k,g) = DispY(k,g) + DispYY(k,g,j); end DispY(k,g) = E(l,2)*DispY(k,g)/(Sigma*tt*S‘3); end end %Z Displacement at middle height(W) Plot k = midlayer; n = k-l; DispZ=O; for j = 1:3 DispZZ(k,g,j) = R(k,j)*(FG(n*6+3+j)*CM(k,j)... +a1fa(k,j)*FG(n*6+j)*SM(k,j)); Dispz = Dispz + DispZZ(j); end DispZ; DispZ = 100*E(1,2)*DispZ/(Sigma*tt*SA4); %X Displacement (W) Plot for k = 1:N n = k-l; for g = 1:point+1 DispZ(k,g)=O; for j = 1:3 DispZZ(k,g,j) = (FG(n*6+j)*WIZ(k,g,j)... +FG(n*6+3+j)*W2Z(k,g,j)); DispZ(k,g) = DispZ(k,g) + DispZZ(k,g,j); end DispZ(k,g) = 100*E(1,2)*Dispz(k,g)/(Sigma*tt*s*4); %plot(DispZ(k,g),hh(k,g),‘*b'); end end 141 Appendix E Figures of stresses and displacements obtained using the linear shear slip theory and linear normal separation theory in the Pagano’s solutions 0.5 1 r 1 l 1 ' . -.- ,_,....-»-‘ 0.4 ~ . 7 ""' - 0.3 - ' - ’ +sz=Dsy=1E-11 and K=4.17E-12 0 2 _ +sz=Dsy-1E-9 and K84.17E-1O _ ‘ if-.. - : -_ +sz-Dsy=1E-7 and K-4.17E-8 - +sz=Dsy=1E-5 and K=4.17E-6 0.1 _ +sz=Dsy=1E-3 and K=4.17E-4 ‘ +sz=Dsy=1 El and Kfl4.17E-2 N 0 _ +sz=Dsy=1E+1 and K=4.17E0 _ -0.1 — a at . ‘ _O 2 _ r _ '1’ -o 3 _ 'f _ -o.4 ~ , ' ' — - 1 1 1 1 1 Mia .4 -2 4 6 — 2 axfa/2,b/2,z) FigureE.l-Normalized in-plane normal stress distribution along the height of laminate by varying the shear slip constants D“ and D3y and normal slip constant k of a [0/90/0] composite plate 0.5 r .‘ ' A , +osx=osy=1 E-11 and K=4.17E-12 f . +sz=osy=1E-9 and 104.1 7E-10 0‘4 ” - +sz=DsF1E-7 and K=4.17E-8 ‘ ’ —0-— sz=Dsy=1E-5 and K=4.17E-6 0.3 - " +sz=Dsy=1E-3 and K=4J7E-4 j ,' +sz=Dsy=1 E1 and K=4.17E-2 0 2 _ , , +sz=Dsy=1 E+1 and K=4.17EO H / : I _- A N 0 '- ‘ 7’1 d .0.1 _ _ v v ‘. -o.2 — u _ O 0 I1 -0.3 - ,3 —‘ -o.4 — -« - 1 1 1 1 J 0'3 .5 -‘| -O.5 O 0.5 1 1 .5 ay (a/Z ,b/Z , z) FigureE.2-Normalized in-plane normal stress distribution along the height of a [0/90/0] composite plate by varying the shear slip constants D” and D”, and normal slip constant k 142 0.5 i l I r l I 0.4 — e 0.3 — — 0.2 — — 0.1 — + sz=Dsy=1E~11 and K=4.17E-12 a +sz=Dsy=1E-9 and Kd.17E-1O N O __ + sz=Dsy=1E-7 md KI4.17E-8 .1 —e— sz=Dsy=1E-5 md Kati 7E-6 + sz=Dsy=1E-3 and K=4.17E-4 -0-1 r- —e— sz=Dsy=1E-1 and K=4.17E-2 ‘ —+— sz=Dsy=1E+1 and K=4.17E0 -0.2 — « -O.3 - - -0.4 - a _ 1 1 1 1 1 1 038.2 O 0.2 0.4 0.8 1 1 .2 _ 0.8 az(a/2,b/2, z) FigureE.3-Normali2ed transverse normal stress distribution along the height of a [0/90/0] composite plate by varying the shear slip constants D” and Dsy and normal slip constant k 0.5 1 1 T l 1 1 1 l r ‘ ( +sz=Dsy=1E-11 and K=4.17E-12 \ ‘_ I +sz=Dsy=1E-9 and K=4.17E-1O 0,4 - \ ,_ p ‘ ‘ +sz-Dsy-1 E-7 and K=4.17E-8 1 k ‘ —e— sz=Dsy=1E-5 and K=4.17E-6 t, + DSXBDsyI1 E-3 Ind K34.17E-4 0.3 — I + sz=Dsy=1 E1 and K=4.17E—2 _ + sz=Dsy=1E+1 and K=4.17E0 o 2 — \ _ 0.1 — _. N 0 — _. -o.1 — — -O.2 - - -O.3 - .4 -0.4 — x ‘ -0 l l l J I l l l l - .25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 z:(0,0, z) FigureE.4-Normalized in-pla'ne shear stress distribution along the height of a [0/90/0] composite plate obtained by varying the shear slip constants D“ and D3), and normal slip constant k 143 0.5 1 | I l l 1 l I +sz=Dsy=1E-11 and K=4.17E-12 0.4 _ : +sz=osy=1E-9 and K=4.17E-10 — 1 I +sz=Dsy=1E-7 and K=4.17E-8 o 3 _ “ ' +sz=osy=1E-5 and K=4.17E-6 _ ‘ +sz=Dsy=1E-3 and K=4.17E-4 1 +sz=Dsy=1E-1 and K=4.17E-2 0.2 ~ 1‘ +sz=DsF1E+1 and K=4.17E0 - 11 1 v v v 0-111' 1 ‘ 11 1 41 t N on» 1 _ 11 1 11 1 11_ 1 _ -0.1“ ‘ 1 1 Mi I l 1 1 r 1 I V 1 -0.3 r 1 _ 1 r 1 1 r 1 I _ 1! _ 0'41' 1 qr _o 5 | l l l I l I 0 4 6 _ 10 12 14 16 W(0,bf2 ,2) FigureE.5—Normalized transverse displacement variation along the height of a [0/90/0] composite plate obtained by varying the shear slip constants D“ and D:y and normal slip constant k Appendix F Figures of stresses obtained by using the linear shear slip theory in ISSCT 0.9 1 T 1 1 . ~ , ”paw .— " ‘1 '1 L If J 0.8 — _ ._ - , ' _ . , ' —B- sz. Dsy-1 E-11 Q. _. . ~ ‘ +sz,Dsy-1 E-9 0.7 ~ W ._ ' —+— sz.Dsy—1E~7 _. Wfi' .-. -—v—— sz.Dsy-1E-5 ”fies-‘7' _.- —1-— szDsy-1E-3 W. HIE? - _.a _ 0-6 " v v ' E -e— szDsy-l E-1 ‘ T: —€— 091. Day-1 E+1 05 L— ' .. N 0.4 - ‘ — 0'3 _ I: 123' gar-9W ‘ 0.2 — - M a 0.1 " ~ —.- _-- .03 . V. —1 ‘ .. .7 0“. ‘ 1 ' 1 1 1 go ~20 -10 10 20 30 —— o 0x(a/2,b/2,z) FigureF.l-Normalized in-plane normal stress distribution along the height of a [0/90/0] composite laminate obtained by varying the shear slip constants D111 and Dyy 0-9 1 1 1 - 1 I ; +sz,Dsy-1E-11 —e—sz,Dsy-1 E-9 0.3 — —1—sz.Dsy-1 E-7 * . .- +szDsy-1E-5 0 7 — +sz,Dsy-1E-3 ' —e—os1gosy-1 E-‘l [ +sz.Dsy-1E+1 0.6 - _, —. .- - 0.5 — ._ N 0.4 h M 1: i 0.3 — as”??? - 0.2 - a 0.1 - — l J J l l 330 -20 -10 10 20 30 E(a/2,b/2,z) FigureF.2-Normalized in-plane normal stress distribution along the height of a [0/90/0] composite laminate obtained by varying the shear slip constants D” and Dsy 145 0.9 W i 1 I T 1 ~ -. -1: ,1 ._ +sz,Dsy-1E-11 ' i 1": ; ; . +sz,Dsy-1E-9 03- ' ijns, -d—oupwds7* ' " ' . +szDsy-1E-5 07__ , ‘ ' i121, +sz.Dsy-1E-3 . ' ‘ 7 3's; ;- ; , —e-sz.Dsy-1E-1 1' 1..“ i ?' :- _ +sz,Dsy-1E+1 1 _a .4 0.6 os— “¢§§h,i _ I g; I o4— ‘1«$;-_ _ 0.3 " Olly. ; V ., 5 ‘Tr‘ -‘ 0.2 — -‘ 01— ‘~. g - A pf?! .— 1 1 1 1 I ' ‘ . . R .5 -1 -0.5 0.5 1 1 .5 ti;801lzj FigureF.3-Normalized in-plane shear stress distribution along the height of a [0/90/0] composite plate obtained by varying the shear slip constants D” and Dsy 146 Appendix G Figures of stresses and displacements obtained using the Embedded Layer approach in ISSCT 1 l l l l l 0.9” .‘_l .- - 7:17, _3 : _1 . -—. .'-.' 1'" '3’. I ‘1..- - - . OB" ‘1').{'_-~"' 3"- _1 0.7 ~ . .1»... a“? ‘ + E=2.4Ee & G=1 E6 — _ . _ gr .; *- ' 11:6 +E=2.4E4&G=1E4 o s — 1 = = ‘ . + E=2.4E2 a. G=1 E2 j - E —1— E=2.4Eo a. G=1EO —e—E=2.4E-2 a G=1E-2 ~ 0.5 — +E=2.4E-4 a G=1E—4 ~ 0.4 — — 0.3 " I "5". - - _; a . _ ., ~: —1 ‘1' 7 _'M . 02 _ - -. ’ _ .1 1:" 7",.“2‘. . ... 0.1 ~ , 1 _- «1+ ‘ ”91.2.5“: :1: _' .. ‘ I 1‘ I I ,: _.'_—. —. ' :1"! 1 1 1 310 -20 -1o _0 1o 20 30 O'x(a/2,b/2,z) FigureG.1-Normalized in-plane normal stress distribution along the height of a [0/1/90/1/0] composite plate obtained by varying the material properties of the embedded layers. [1 — isotropic] 1 I l W l i + E=2.4E6 & G=1 E6 0 9 __ + E=2.4E4 81 G=1E4 . ' + E=2.4E2 81 G=1 E2 -+— E=2.4E0 81 G=1E0 0-8 ‘ ‘9“ E=2.4E-2 & G=1E-2 ‘ —B— E=2.4E—4 81 G=1E—4 0.7 - .. 0.6_ ‘3-2 ._ _ 1,“. . 2. : -1 _. - 7. v5.32, :- N 0.5 - . ' _'-—, -= 1"" a 0.4— r .s-‘éfz'fir'??-i:..fl( —J o.3~ : '- ’ " ' '-"‘ — 0.2 *- ~ 0.1 *- A l l l l l 350 -20 -10 0 10 20 30 ay (a/Z ,b/Z , z) FigureG.2-Normalized in-plane normal stress distribution along the height of a [0/1/90/1/0] composite plate obtained by varying the material properties of the embedded layers. [1 — isotropic] 147 0.9 0.8 0.7 0.6 N 0.5 0.4 0.3 0.2 0.1 l l l + E=2.4E6 a G=1 E6 —a— E=2.4E4 a G=1 E4 + E=2.4E2 a G=1 E2 “ -—1— E=2.4E0 & G=1 E0 -6— E=2.4E-2 a. G=1E-2 ~ —a— E=2.4E-4 a. G=1E-4 _- 1— .0.. Txy (0,0, 22) FigureG.3-Normalized in-plane shear stress distribution along the height of a [0/1/90/1/0] composite plate obtained by varying the material properties of the embedded layers. [1 -— isotropic] 148 Appendix H Figures of stresses and displacements obtained by comparing ISSCT with Pagano for imperfect bonding condition 0.5 1 1 1 1 1 0.4 - , _ 0.3 - j 0.2 — '- ' —l— Pagano-ShearSLip & NormalSlip _ , + Pagano- Embedded Layer ' _'0- lssct-ShearSlip 0'1 ~ " +Issct-Embedded Layer ‘ N 0 L -o.1 — — -o 2 E b r i _ -o.3 ~ " _ -o.4 — _ -0536 ‘1‘ I l 4 e _ 6 o;(a/2,b/2,z) Figure [1.1-Comparison of normalized in-plane normal stress along the height of a delaminated [0/90/0] plate between Lee and Pagano [sz=Dsy=0.1, k=0.0417, E=0.024 psi,G=0.01 psi] 0.5 I fl 1 1 1 1 I —t— Pagano-ShearSLip 8. NormalSlip O 4 __ +Pagano-Embedded Layer H ‘ +lssct—ShaarSlip +Issct—Embedded Layer 0.3 — _ 0.2 - _. 0.1 — - N o r— J -O.1 '- -1 -02 _ .. -O.3 — f -O.4 _ 1 l l l l l -O.§2 _1_5 -1 -0.5 1 .5 2 3(a/2,b73,z) Figure [1.1-Comparison of normalized in-plane normal stress along the height of a delaminated [0/90/0] plate between Lee and Pagano [sz=Dsy=0.1, k=0.0417, E=0.024 psi,G=0.01 psi] 149 0'5 l l l l T l , . -t- Pagano-ShearSLip & NormalSlip 0 4 _ rs; +Pagano-Embedded Layer 1 ' "“x ‘ -O- lssct-ShearSlip + lssct-Embedded Layer 0.3 - — 0.2 - - 0.1 — — N 0 - '1 -0.1 — A .02 _ —1 -0.3 — — -0.4 - - -0 1 1 L 1 1 1 1 1 1 3.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 02 0.25 E(Ofi, 2) Figure 113-Comparison of normalized in-plane shear stress along the height of a delaminated [0/90/0] plate between Lee and Pagano [sz=Dsy=0.l, k=0.0417, E=0.024 psi,G=0.01 psi] 150 Appendix I Figures of stresses obtained by comparing ISSCT-Shear Slip only with Pagano for imperfect bonding conditions 0-5 1 1 1 1 1 1 l 0.4 — -1 0.3 — - ~ 0.2 ~ , —-o— Paganosz-Dsy-fl E-7 — r - +ISSCT-sz=Dsy=1E-7 O 1 _ I —t-Pagano-sz=Dsy=1 E1 _4 ' +ISSCT-sz-Dsy-1E1 N o 1— —1 1' no 1 P ‘ --1 0.2 - — -0.3 - i .04 .. _ _ 1 1 1 1 1 1 1 0.52 -1 .5 -1 -0.5 O 0.5 1 1.5 2 ax(a/2,b/2,z) Figure 1.1 Comparison of normalized in plane normal stress along the height of a [0/90/0] composite plate for two different bonding conditions 0.5 1 1 1 1 1 -o— Emt-szstyfi B7 0 4 +Issct-szsosw1E-7 ‘ F -0— Exact-Dwosw1 E1 _ + Issct-sz=Dsy=1 E1 0.3 '— -—1 0.2 - — 0.1 - - N 0 — — -O.1 — — 0.2 — - 0.3 — _ -O.4 e - _ 1 1 1 1 1 1 1 O.§2 -1.5 -1 -O.5 0 0.5 1 1 .5 2 0'), (a/2,b/2, 2) Figure 1.2 Comparison of normalized in plane normal stress along the height of a [0/90/0] composite plate for two different bonding conditions 151 BIBLIOGRAPHY [1] S.P.Chang, “Analysis of composite laminates using interlaminar stress continuity theory with interfacial slip”, MS. Thesis, Da—Yeh University, Taiwan, 2005 [2] The world of Explosives, "The Many Uses of EXplosives", http://www.explosives.org/Useshtm, downloaded 4/28/07, [3] Baker, W.E., et.al. “Explosive hazards and evaluation”, Amsterdam; New York: Elsevier Scientific Pub. Co., 1983 [4] Brode, H.L., “Numerical solutions of spherical blast waves”, Journal of Applied Physics, 26(6), pp. 766-775, 1955 [5] Smith, RD. and Hetherington, J.G., “Blast and Ballistic Loading of Structures”, Butterworth Heinemann, 2004 , [6] Ronald T. Noyes, “OSU Current Report”, http://osuextra.okstate.edu/pdfs/CR- 1737web.pdf. downloaded 5/01/2007 [7] Ccps, “Guidelines for use of vapor cloud dispersion models”, Wiley-AIChE, 2nd Edition, 1996 [8]Baker, W.E., “Explosions in Air”, University of Texas Press, Austin, 1973 [9]Beshara, F .B.H., “Modeling of blast loading on above ground structures-I. General Phenomenology and External Blast”, Computers and Structures, 51(5), pp. 585-596, 1994 [10] Beshara, F.B.H., “Modelling of blast loading on above ground structures-II. Internal blast and ground shock”, Computers and Structures, 51(5), pp. 597-606, 1994 [11]Kinney, G. F. And Graham, K. J ., “Explosive shocks in air”, Springer Verlag, Berlin, New York, 1985 [12]Saleh A. and Adeli H., “Optimal Control of Adaptive Building Structures under Blast Loading,” Mechatronics, 8, pp. 821-844, 1998 [13] Fung, TC. and Chow, S.K., “Responses of Blast Loading by Complex Time Step Method,” Journal of Sound and Vibration, 223 (1), pp. 23-48, 1999 [14] Ethridge, N.H., “A procedure for reading and smoothing pressure time datafrom HE. and Nuclear Explosions”, BRL Memo Report No. 1691, Aberdeen Proving Ground, Md. 152 [15] Brode, H.L., “Point source explosion in air”, RM-l 824-AEC, the Rand corporation, Santa Monica, California [16] Shu X., “A generalized model of laminated composite plates with interfacial damage”, Composite Structures, 74, pp. 237-246, 2006 [17] Lu, X. and Liu, D., “An Interlaminar shear stress continuity theory for both thin and thick composite laminates,” Journal of Applied Mechanics, 59, pp. 502-509, 1992 [18] Della, C.N., Shu, D., “Vibration of delaminated composite laminates: A review”, Applied Mechanics Review, 60, pp., 1-20, 2007 [19] Reddy, J.R., Mechanics of Laminated Composite Plate and Shells: Second Edition, Chapter 3, CRC Press, Boca Raton, FL, 2004 [20] Ghughal, Y.M. and Shimpi, R.P., “A review of refined shear deformation theories of isotropic and anisotropic laminated plates”, Journal of Reinforced Plastics and Composites, 21 (9), pp. 775-819, 2002 [21] Roy, Tarapada and Chakraborty, Debabrata, “Delamination in Hybrid FRP Laminates under low velocity impact”, Journal of Reinforced Plastics and Composites, 25 (18), 2006 [22] Choi, H. Y., His-Young, T. W. and Chang, F. K, “A New Approach towards Understanding Damage Mechanism and Mechanics of Laminated Composites Due to Low Velocity Impact: Part II — Analysis”, Journal of Composite Materials, 25, pp.1012— 1038, 1991 [23] Carrera E., “Historical Review of zig-zag theories for multilayered plates and shells”, Applied Mechanics Review, 56 (3), pp., 287-308, 2003 [24] Ghosal, A., Kim,H.S., Chattopadhyay, A., Prosser, W.H., “ Effect of delamination on transient history of smart composite plates”, Finite Element in Analysis and Design, 41, 850-874, 2005 [25] Ghoshal, A., Kim,H.S., Kim,J., Choi S.B., Prosser, W.H., Tai H., “ Modeling delamination in composite structures by incorporating the F ermi-Dirac distribution function and hybrid damage indicators”, Finite Element in Analysis and Design, 42, 715- 725, 2006 [26] Lu, X. and Liu, D., “Interlayer shear slip theory for cross-ply laminated with nonrigid interfaces”, AMA, 60 (4), pp. 1063-1073, 1992 [27] Liu, D., Xu, L. and Lu, X., “Stress Analysis of Imperfect composite laminates with an interlaminar bonding theory”, International Journal for Numerical Methods in Engineering, 37, pp. 2819-2839, 1994 153 [28] Shu, X., Soldatos, K.P, “An accurate stress analysis model for angle-ply laminates with weakly bonded layers”, Acta Mechanica, 150, pp. 161-178, 2001 [29] Shu, X., Soldatos, K.P, “An accurate de—lamination model for weakly bonded laminates subjected to different sets of edge boundary conditions”, International Journal of Mechanical Sciences, 43, pp. 935-959, 2001 [30] N.J.Pagano, “Exact solutions for rectangular bidirectional composites and sandwich plates”, Journal of Composite Materials, 4, 1970 [31] Beshara, F .B.A., “Nonlinear finite element analysis of reinforced concrete structures subjected to blast loading”, Ph.D. thesis, City University, London, 1991 [32] Balden, V.H. and Nurick, G.N., “Numerical simulation of the post-failure motion of steel plates subjected to blast loading”, International Journal of Impact Engineering, 32, pp. 14-34, 2005 [33] Librescu, L. and Nosier, A., “Response of laminated composite flat panels to sonic boon and explosive blast loadings”, AIAA, 28 (2), pp. 345-3 52, 1990 [34] Turkmen, HS. and Mecitoglu, Z, “Nonlinear structural response of laminated composite plates subjected to blast loading”, AIAA, 37( 12), pp. 1639-1647, 1999 [35] Turkmen, HS. and Mecitoglu, 2., “Dynamic response of a stiffened laminated composite plate subjected to blast load”, Journal of Sound and Vibration, 221(3), pp. 371-389, 1999 [36] Zyskowski, A., So'chet, I., Mavrot, G., Bailly, P., Renard, J ., “ Study of the explosion process in a small scale experiment-structural loading”, Journal of Loss Prevention in the Process Industries, 17, pp. 291-299, 2004 [37] O’Toole, B., Trabia M. B., Thota J., Wilcox, T., Nakelswamy, K.K., “ Structural response of blast loaded composite containment vessels”, SAMPE Journal, 42 (4), pp. 6- 13, 2006 [38] S. Manchiraju, “A progressive damage analysis including delamination for fiber reinforced composites”, MS. Thesis, Michigan State University, Michigan, 2004 [39] Li, G., Li, Q., Liu, D., Raju, BB. and Templeton, D.W., “Designing Composite Vehicles against Blast Attack,” Society of Automotive Engineering, Paper 2007-01-0137, 2007 [40] Sih, G.C. and Skudra, A.M., “Failure Mechanics of Composites”, Handbook of Composites, Volume 3, Elsevier Science Publishers, New York, NY, 1985 154 [41] Norris, C.B., “Strength of orthotropic material subjected to combined stress”, US Forest Products Laboratory Report # 1816, 1950 [42] Azzi, V.D. and Tsai, S.W., “Anisotropic Strength of Composites”, Experimental Mechanics, 5(9), pp. 283-288, 1965 [43] Fischer, L., “Optimization of orthotropic laminates”, Engineering for Industry, 89 (3), Series B, pp. 399-402, 1967 [44] Zhou, S. G. and Sun, C.T., “Failure analysis of composite laminates with free edge”, Journal of Composite Technology and Research, 12, 1990, pp. 91-97 155 11111111111111111