

small (250x250 max)
medium (500x500 max)
Large
Extra Large
large ( > 500x500)
Full Resolution


A COUPLED HYGROTHERMAL COHESIVE LAYER MODEL FOR SIMULATING DEBOND GROWTH IN BIMATERIAL INTERFACE By YONG WANG Bachelor of Science Tsinghua University Beijing, China 1986 Master of Science Tsinghua University Beijing, China 1989 Submitted to the Faculty of the Graduate College of the Oklahoma State University in partial fulfillment of the requirements for the Degree of DOCTOR OF PHILOSOPHY July, 2006 ii COPYRIGHT By Yong Wang July, 2006 iii A COUPLED HYGROTHERMAL COHESIVE LAYER MODEL FOR SIMULATING DEBOND GROWTH IN BIMATERIAL INTERFACE Thesis Approved: Thesis Adviser Dean of the Graduate College iv ACKNOWLEDGMENTS I would like to express my sincere thanks to my advisor, Professor Samit Roy, for his supervision, constructive guidance, financial support and inspiration throughout the study. His profound contribution to my development, professional and otherwise, is deeply appreciated. My appreciation extends as well to the members of my supervisory committee Dr. Hongbing Lu and Dr. Andrew S. Arena, Jr. from Department of Mechanical and Aerospace Engineering, Dr. Roger Zierau from Department of Mathematics, and to the faculty of the Department of Mechanical and Aerospace Engineering. I also would like to express my thanks to Professor Kenneth M. Liechti and his group member Dr. Soojae Park and Mr. Dewei Xu from The University of Texas at Austin for providing test data and valuable suggestions for this study. I wish to express my sincere gratitude to my parents for their unconditional love and confidence in me all the time. This dissertation is dedicated to my beloved wife Bo Lan and my daughter Hongfan Wang, without their love and support this work would have not been possible. I would like to thank them for their consistent encouragement and love, and for all that they have given to me in my life. Finally, the financial support of the National Science Foundation (Grant No. CMS0296167) is acknowledged with thanks. v TABLE OF CONTENTS Chapter Page I. INTRODUCTION.……………………………………………………………….……....1 1.1 Applications of Adhesive Material………………………………………..……….1 1.2 Cohesive Zone Model…………………………………………………..….....……3 1.3 Coupled Hygrothermal Effect on Cohesive Layer……………………..….…........5 1.4 Cohesive Layer Model……………………………………………..………………6 1.5 Objective and Contents…………………………………………..………….……..7 II. COHESIVE LAYER MODEL……………………………………………………….11 2.1 Cohesive Layer Configuration………………………..……..………………..…..11 2.2 Cohesive Layer Constitutive Equations……………………………….………….13 2.2.1 Triangular StressStrain TractionSeparation Law………………………......13 2.2.2 2D Cubic StressStrain TractionSeparation Law……………...……….......14 2.2.3 3D Cubic StressStrain TractionSeparation Law………………………..…16 2.3 Cohesive Layer WorkofSeparation…………………..…………………....……17 2.4 Failure Criteria of Mixed Mode Debond………………………………..……..…17 2.4.1 Criterion Based on Prescribed Maximum Strain………………………….....18 2.4.2 Criterion Based on Strain Energy Release Rate…………………..………....18 III. DEVELOPMENT OF AN ANALYTICAL SOLUTION FOR COHESIVE LAYER MODEL AND MODEL VERIFICATION…………………...…………...22 vi 3.1 Analytical Solution from Cohesive Zone Model……………………….……...…22 3.2 Analytical Solution from Cohesive Layer Model…………………………..….…27 3.3 Comparison between Analytical Solution and FEA Results………………..........32 3.4 Conclusions……………………………………………..........................…….…..40 IV. HYGROTHERMAL EFFECT ON COHESIVE LAYER………………………….42 4.1 Introduction……………………………………………………..…..….…………42 4.2 Fickian Diffusion…………………………………………………..………..……43 4.3 NonFickian Diffusion (Strain Assisted Diffusion)………………….………...…44 4.4 Cohesive Layer Diffusion Boundary Conditions……………………………....…51 4.5 Cohesive Layer WorkofSeparation………………………………………..……52 4.6 Cohesive Layer Degradation and FEA Simulation………………………..…..….54 4.6.1 Cohesive Layer Degradation Due to Moisture Concentration………….....…54 4.6.2 FEA Simulation Results…………………………………………….……......55 4.7 Conclusions………………………………………………………………..…...…61 V. NUMERICAL SIMULATION OF PENINSULA BLISTER TEST………………...62 5.1 Introduction……………………………………………………………….…..…..63 5.2 Constitutive Law and Failure Criterion of Cohesive Layer…………………...….67 5.3 Specimen Geometry and Energy Release Rates…………………..…...............…67 5.4 TimeDependent Behavior of Polymeric Thin Film…………………………...…71 5.5 Simulation of Debonding in the Peninsula Blister Test………………………..…72 5.5.1 Simulation under Small Deformation…………………………...……….......72 5.5.2 Peninsula Blister Specimen and Test Results…………………...…………...74 5.5.3 Simulation with Large Deformation and Residual Stresses………..…...…...75 vii 5.5.4 Simulation Including TimeDependent Effect………………………….……81 5.6 Conclusions……………………………………………………….………...…….86 VI. SIMULATION OF TIMEDEPENDENT DEBOND GROWTH…………….……87 6.1 Introduction………………………………………………………………....…….87 6.2 Failure Criterion Based on WorkofSeparation……………………….…........…88 6.3 Nonlinear Viscoelastic Model and Fracture Energy…………………..….…....…90 6.4 Moving Wedge Test……………………………………………………......…..…91 6.5 Numerical Simulation of Wedge Test……………………………………...……..92 6.6 Conclusions……………………………………………………………..…….......98 VII. FRP BONDED STRUCTURE UNDER BLAST LOAD…………….……….….100 7.1 Introduction……………………………………………………….…..…….…...100 7.2 Implicit Integration Methods…..………………………………..…..………..…101 7.2.1 Houbolt Method…………………………………………….………..……..102 7.2.2 Newmark Method………………………………………………………......103 7.2.3 Wilsonθ Method…………………………………………….……….……104 7.2.4 HHT Method………………………………………………………………..104 7.3 Blast Load……………………………………………………….………………106 7.4 Dynamic Response of FRP Bonded Structure………………………….……….109 7.4.1 Short Term and Long Term Responses in Cohesive Layer…..………….…109 7.4.2 Critical Debond Locations of FRP Bonded Beam...…………….……..…...112 7.4.3 Explosive at Different Locations………………..………………………….113 7.4.4 Concrete Beam with Initial Cracks………………..………………………..115 7.5 Conclusions……………………………………………………………………...118 viii VIII. DYNAMIC ANALYSIS WITH MATERIAL AND GEOMETRIC NONLINEARITY………………………………………………………………..119 8.1 Modeling of Material and Geometric Nonlinearity…………………………..…119 8.1.1 Material Nonlinearity………………………………………………….……120 8.1.2 Numerical Formulation for ElastoPlastic Problems……………………….123 8.1.3 Geometric Nonlinearity…………………………………………………….128 8.1.4 HilberHughesTaylor (HHT) method for temporal discretization………...131 8.2 Model Verification………………………………………………………………131 8.2.1 Single Element Extension Verification …………...………………………..132 8.2.2 Simple Shear Verification…………………………………………………..133 8.3 Dynamic Response of a Circular Steel Plate with Coating under Blast Load......137 8.3.1 Modeling of a Clamped Circular Steel Plate with Coating……………..…..137 8.3.2 Numerical Damping of HHT Method………………………………………138 8.3.3 Dynamic Responses under Different Load Levels………………………....141 8.3.4 Effects of Coating on Plate Responses………………………………..……143 8.3.5 TimeDependent Effect of Polymeric Coating……………………………..147 8.4 Conclusions………………………………………………………...……………149 IX. CONCLUSIONS……………………………………………….…………….……151 BIBLIOGRAPHY.………………………………………………………………...…....154 ix LIST OF TABLES Table Page 31 Material properties for concrete and epoxy adhesive………….…………………...33 32 Variation of damage length with different cohesive layer thicknesses………….…39 33 Analytical and numerical solutions of damage length……………………………...40 51 Different element size in peninsula direction………………………………...…….76 52 Effect of residual stresses on debond process as predicted by FEA simulations .....77 53 Material properties of a viscoelastic epoxy …………………………………….….83 54 Peninsula blister test data and FEA predictions with elastic and viscoelastic thin film……………………………………………………………..…85 61 max σ and max ε under different strain rate ε&……………………………………..…94 62 FEA results with various wedge speeds……………………………………………97 71 Impulses of different load distributions...…………………………………………108 72 Material properties of concrete, epoxy adhesive and FRP…………….……....….110 81 Constants defining the yield surface………………………………………………124 82 Material properties of steel plate and polymer coating…………………………...138 x LIST OF FIGURES Figure Page 2.1 Three characteristic zones near debond tip and corresponding stressstrain relations ...………………………………………………………………………..…………..12 22 Normalized triangular stressstrain tractionseparation law for a cohesive layer.....14 23 Cubic stressstrain tractionseparation law for a cohesive layer……………….…..15 24 Cubic stressstrain tractionseparation law and strain energy release rate for a cohesive layer……………………..……..……………………………………19 31 A cohesive layer in a double cantilever beam (DCB)……………………………...23 32 Schematic debonding of a double cantilever beam (DCB) showing cohesive layer with elastic zone and damaged zone ( x1∈[0,l])……….…….….....23 33 Transverse deformation comparison between beam centerline and cohesive interface using FEA…………………………………….…….…………..28 34 Transverse stress distribution in the beam and cohesive layer…………………….29 35 Finite element mesh of a DCB specimen with symmetry boundary conditions…...33 36 Stress distribution with different correction coefficient k …………………………35 37 Reaction force comparison for different normalized cohesive layer thicknesses…36 38 Reaction force comparison for different debond lengths…………………………..37 39 Damage length comparison for different cohesive layer thicknesses………………38 310 Stress distributions for different cohesive layer thicknesses…………………..….39 xi 41 A cohesive layer with moisture diffusion in a DCB beam…………………………45 42 Depiction of the influence of moisture in cohesive layer on workofseparation….55 43 Comparison of predicted moisture concentration profiles.………………………...56 44 Debond growth in a cohesive layer due to moisture degradation……..……………57 45 Transverse mechanical strains and corresponding failure zones at different moments ……………………………………………….………………58 46 Transverse stresses along bond length at different moments………………………59 47 Predicted evolution of debond length with time……………………………………60 48 Reaction force decreases with time………………………………..……………….61 51 Peninsula blister specimen and possible debond sites…………………………...…65 52 Effect of relative peninsula width on the energy release rate………..……………..69 53 Effect of blister region aspect ratio on the energy release rate (Analytical solution is derived from l 2a > 2)…………………………………….70 54 Peninsula blister responses over large debond lengths..……………………………73 55 3D FEA mesh with 20node brick element………………………………………..76 56 FEA simulation of peninsula blister responses (with large deformation and residual stresses)…………..………………….………79 57 Simulated 3D film displacement profiles at different stages of debonding….……80 58 Predicted liquid pressure and film deflection increase when debond approaches the end of the peninsula……………………...…………………….…..81 59 Peninsula blister test data and FEA simulation results (with large deformation and residual stresses)………….…………………….…….84 61 Stressstrain relation of epoxy under different strain rates (test data)………….…..89 xii 62 Specimen of a moving wedge test……………………………………………….…91 63 Debond length vs. time (test result and FEA prediction)…………………………..93 64 Vertical reaction force vs. time (test result and FEA prediction)………………….93 65 Fracture energy 2Γ vs. debond speed……………………………………………..95 66 FEA mesh and contour for Jintegral………………………………………………97 71 A simply supported concrete beam bonded with FRP under blast load…………..107 72 Triangular and exponential distributions of blast load……………………..……..107 73 Schematic FEA mesh and distribution of blast load along the beam……………..109 74 Short term response y σ in the top and bottom cohesive layer at midspan...…….111 75 Long term response x σ in the top and bottom cohesive layer at midspan...…….111 76 Typical distribution of stress y σ in a cohesive layer ………………………….…113 77 Typical distribution of stress x σ in a cohesive layer ………………………….…113 78 FEA mesh of a simply supported beam under off side load………………………115 79 Stress limits of x σ as function of explosive location along beam axis…………...115 710 A simply supported FRP bonded beam with an initial crack in concrete …....….116 711 Comparison of axial stress x σ in the cohesive layer before and after concrete cracking (Cracking occurred at t = 50ms )…………………………….117 712 Comparison of shear stress xy τ in the cohesive layer before and after concrete cracking (Cracking occurred at t = 50ms )………………………….….117 81 Mathematic models for representation of strain hardening behavior………….….121 82 Elastoplastic linear strain hardening behavior for uniaxial case…………..……..122 83 Incremental stress changes at a point in an elastoplastic continuum…………….126 xiii 84 Single element extension (plane strain conditions)……………………………….132 85 Stressstrain relation comparison between NOVA and ABAQUS……………….133 86 Single element under simple shear (displacement control)……………………….134 87 Single element shear stress vs. shear strain for simple shear ………….…………134 88 Applied shear force vs. shear strain (global comparison, 19×19 mesh)…………135 89 Shear stress at the center of the specimen (local comparison, 19×19 mesh).……135 810 Equivalent plastic strain at the center of the specimen (19×19 mesh) ……..…..136 811 Deformation of a simple shear specimen with 200% shear strain (19×19 mesh) …………………………………………...………………………………………136 812 Axisymmetric model of a circular steel plate with coating…………....………...137 813 FEA mesh for simulating steel plate with coating (dimensions not scaled)……..138 814 Numerical damping of HHT method on highfrequency modes……………...…139 815 Numerical damping of HHT method on plate deflection………………………..140 816 Central deflection history under various peak pressure levels.……………..…...142 817 Maximum equivalent plastic strain at plate center under various peak pressure levels……………………………………………………...…..….142 818 Minimum and maximum deflections under various peak pressure levels.……....143 819 Central deflection history of the plate…………………………………………...144 820 Radial stress history in polymer coating…………………………………………145 821 History of equivalent plastic strain in steel plate…………………………….…..145 822 Maximum equivalent plastic strain at the bottom surface of steel plate…..…….146 823 Equivalent plastic strain in steel plate and permanent deformation (not scaled) …………………………………………………………………………..………..146 xiv 824 Central deflection history of the plate…………………………………….……..148 825 Radial stress history in steel plate at the plate center……………………...…….148 826 Radial stress history in polymer coating at the plate center ……………….…....149 xv NOMENCLATURE ij σ Stress component ( Pa ) ij ε Strain component σ Normal stress ( Pa ) τ Shear stress ( Pa ) max σ Maximum normal stress in cohesive layer ( Pa ) max τ Maximum shear stress in cohesive layer ( Pa ) max ε Maximum normal strain in cohesive layer max γ Maximum shear strain in cohesive layer c h Cohesive layer thickness (m ) φ sep Cohesive layer workofseparation ( J /m2 ) G Strain energy release rate ( J /m2 ) E Young’s modulus ( Pa ) υ Poisson’s ratio ρ Material density (Kg /m3 ) ψ Helmholtz free energy ( J /m3 ) m Moisture concentration (Kg /m3 ) xvi Γ Fracture energy ( J /m2 ) R External nodal force vector ( N ) F Internal nodal force vector ( N ) M Mass matrix K Stiffness matrix Δt Time step (sec ) D Elasticity matrix ep D Elastoplastic stiffness matrix ys σ Yield stress ( Pa ) H′ Strain hardening function ( Pa ) p ε Equivalent plastic strain σ Equivalent stress ( Pa ) i J Stress invariant i J ′ Deviatoric stress invariant 1 CHAPTER I INTRODUCTION 1.1 Applications of Adhesive Material The use of structural adhesives is rapidly increasing as they offer distinct advantages over conventional mechanical fastening technique. Laminated composites and thin film structures are some of the most popular applications in industry. Fiber reinforced polymer (FRP) composite is a new application of cohesive material in industry and civil engineering. It is necessary to better understand the nature and the reliability of the bonding between the layers or at the bimaterial interface. Many polymeric materials, including structural adhesives, exhibit nonlinear and timedependent behaviors and quite sensitive to the change of temperature and moisture penetrant concentration. Therefore the loadcarrying requirements of the cohesive layer, timedependent material properties, and coupled hygrothermal effect in the adhesive layer, stress concentration and nonlinear deformation near the crack tip, crack initiation and propagation (debond) within the interface between two materials are some of the essential factors that should be considered to get a reliable solution. Fiber reinforced polymers are a class of advanced composite materials that have been extensively used as lightweight, performanceenhancing materials in aerospace, automobile, and defense industries for quite some time. Over the past few years there has been extensive research into their potential applications in the construction industry. 2 However, the actual application of FRP composite in civil engineering sector has been slow especially as concrete reinforcement. One of the chief reasons for their slow acceptance is because of a lack of reliable predictive models and sound design guidelines for their use in civil infrastructure applications. One promising prospect of the application of FRP composite materials in civil engineering is infrastructure repair and retrofit. FRP materials have been used to strengthen the concrete beam element of buildings and bridges with low cost and high strengthening effect [1]. This allowed increasing the strength and/or ductility of these structures while benefiting from the FRP material advantages including: ease of application, high strengthtoweight ratio, and excellent resistance against corrosion and chemical attacks. Ritchie et al. [2] studied the behavior of concrete beams strengthened by bonding FRP (glass, carbon, and aramid) plates to the tension zone and showed that FRP reinforcement increased beam stiffness by 1779% and beam ultimate strength by 4097%. New uses of FRP sheets to upgrade the resistance of steel structures have recently been studied. A considerable increase in the strength and stiffness of the rehabilitated steel bridge girders was observed [3]. A major concern for such retrofitting is the debonding of polymeric adhesive that could compromise the reinforcing effect of the FRP. When exposed to harsh environments, degradation of the adhesive bond could lead to delamination of the FRP reinforcement that could ultimately lead to catastrophic failure. Thin film structure appears mainly as coating in a wide variety of applications and in multilayer structures in microelectronic devices and package. Since poor bonding results in crack or delamination, fracture mechanics is a natural approach for 3 characterizing the resistance to failure and making durability or reliability prediction. At the same time, several experimental methods, including peel test, blister test, indentation test, and scratch test etc., have been used to determine the interfacial bond strength or toughness. One problem common to all of these methods is that global plastic dissipation makes it difficult to extract the true toughness of the interface. A new test method, peninsula blister test, can minimize the plastic dissipation and hence be an effective approach to measure the fracture toughness of the thin film structures. Laminated composites have the advantage of low weight and high strength compared to the structural metals, and hence obtain an increasing application in aerospace industry. Delamination, which is created when two layers debond from each other, is a common type of failure mode in layered composites without the throughthethickness reinforcement. The initiation of delamination growth is usually controlled by mode I and mode II fracture toughness. Numerical approach with mixed mode failure criterion is an effective way to investigate the delamination of the layered composites. 1.2 Cohesive Zone Model An important issue when considering failure is the observation that most engineering materials are not perfectly brittle in the Griffith sense but display some ductility after reaching the strength limit. In fact, for most engineering materials there exists a small zone in front of the crack tip, in which smallscale yielding, micro cracking and void initiation, growth and coalescence take place. If this fracture process zone is sufficiently small compared to the structural dimensions, linearelastic fracture mechanics concept can apply. However, if this is not the case, the cohesive forces that exist in this 4 fracture process zone must be taken into account. The most powerful and natural way is to use cohesive zone model, which was introduced by Barenblatt [4] and Dugdale [5] for elasticplastic fracture in ductile metals, and for quasibrittle materials by Hillerborg et al. [6] in his socalled fictitious crack model. The fracture process zone approach of Needleman [7, 8] and Tvergaard and Hutchinson [9, 10] involves attributing a prescribed tractionseparation law to the interface and, because it allows crack growth to occur, the associated plastic dissipation from loading and unloading of points that are passed by the crack front is rigorously accounted for. As a result, the selected tractionseparation law determines the workofseparation (or adhesive fracture energy), which is the work required to create a unit area of fully developed crack [11]. In the past two decades or so, cohesive zone models have become very popular and have been recognized to be an important tool for describing fracture in engineering materials. Especially when the crack path is known in advance, either from experimental evidence or because of the structure of the material (such as in laminated composites), cohesive zone model has been used with great success. Song and Waas [12], Shawan and Waas [13], and ElSayed and Sridharan [14] successfully employed cohesive zone model to investigate the fracture properties in laminated composites. In those cases, the finite element mesh was constructed such that the known crack path coincides with the element boundaries. The most common failure form of FRP composite plate bonded concrete structure is the delamination of the FRP plate from the concrete component. The debond process of FRP plate usually initiates and propagates along the adhesiveconcrete or adhesiveFRP 5 interface, which is known in advance. The cohesive zone model is thus a good tool for investigation of local fracture processes in FRP delamination. Fatigue crack growth is traditionally characterized via linear elastic fracture mechanics concepts where crack growth rates are correlated with the change in energy release rate or the maximum value of the energy release rate in a cycle. This approach has worked well for metals and polymers alike, especially in dry, room temperature environments, where conditions are still generally linearly elastic. Correlation between crack growth rates and elastic fracture parameters do become suspect in polymers near their glass transition and when saturated by a solvent. Cohesive zone modeling offers a solution to this difficulty in the sense that if the neartip damage can be accounted for in the tractionseparation law of the interphase then the local nonlinear inelastic behavior of the material can be coupled into any analysis directly [15]. 1.3 Coupled Hygrothermal Effect on Cohesive Layer Moisture can cause a host of reliability problems at interfaces including interface bond degradation and debonding. Two mechanisms can be identified. First, moisture at an interface can reduce the interface bonding strength dramatically by altering the chemical bonds. Second, when an interface with a crack or a cracklike defect is under tensile stress, stress corrosion may allow crack growth at stresses much lower than critical fracture would require [16]. The influence of moisture diffusion on crack growth along an interface is not yet fully understood. Environmental cracking in a polymer typically occurs in the presence of a penetrant, such as moisture, and mechanical strain. It has been postulated that the 6 mechanism involved in environmental crack growth in a polymer involves a small zone of craze formation and/or plasticization at the crack tip. For thermoset resins, such as epoxy, energy absorption at the crack tip is primarily by a shear yielding process and not by crazing. Consequently, for a thermoset epoxy, the zone of plasticization ahead of the crack tip must be determined using a diffusion law for nonporous media, such as Fick’s law. However, quite frequently, polymer composites exhibit deviations from the classical Fickian treatment, termed as anomalous or nonFickian diffusion, especially at elevated temperatures and stress levels, and at high relative humidity. Sophisticated hygrothermal models have been developed and verified by Roy [1720] to account for anomalous diffusion. 1.4 Cohesive Layer Model Cohesive layer model employs a thin layer of material, which is placed between two adjacent layers for laminated structures and multilayer structures, or along the bimaterial interface, or along a predicted cracking path in a single material (e.g. concrete and metal) to simulate the elasticplastic failure in ductile or quasibrittle material in the vicinity of the debond tip. It allows debond (or failure) to initiate and grow in these elements along a prescribed debond path. The vicinity of debond tip can be divided into three characteristic zones: elastic zone, damage zone, and debonded zone. The corresponding stressstrain tractionseparation relation may take different forms in each zone for different materials, different loading and environmental conditions to cope with any particular nonlinear behavior in front of the debond tip. 7 The thickness of cohesive layer is an important parameter in the cohesive layer model. It should be noted that the cohesive layer thickness is not arbitrary in cohesive layer model, and it is related to some characteristic lengthscale of the debond process, such as crack opening displacement (COD). Cohesive layer thickness can be determined from the maximum deformation at the debond tip and the maximum strain max ε in cohesive layer at failure. Environmental degradation is usually found in adhesive material and therefore results in the change of the material properties such as max ε and max σ (maximum stress) When investigating the moisture degradation, twodimensional and threedimensional moisture diffusion in the cohesive layer can be directly simulated in a twodimensional or threedimensional cohesive layer. Mixed mode I and mode II fracture is the common failure form of cohesive layer and sometimes even includes mode III. Mixed mode failure and corresponding failure criteria can also be easily implemented in cohesive layer model to predict the debond process in adhesive layer. 1.5 Objective and Contents The objective of this research is to construct a cohesive layer model from fundamental principles of continuum mechanics and thermodynamics, take into account the strain rate dependent material properties, nonFickian Hygrothermal effects as well as diffusioninduced degradation in the cohesive layer. By means of the cohesive layer model, the effect of ratedependent material properties, environmental degradation of the adhesive material, dynamic response involving material and geometric nonlinearity under blast load, quasistatic debond initiation and propagation of the adhesive layer were 8 studied to provide a better understanding of the strengthening effect and reliability of FRP plated structures. In Chapter I, a brief literature review of the application and research status of structural adhesive and FRP bonded structure, and the motivations of this research are presented. In Chapter II, twodimensional and threedimensional cohesive layer constitutive models with prescribed tractionseparation laws were constructed from fundamental principles of continuum mechanics and thermodynamics. Based on debond tip deformation, workofseparation or strain energy release rate, criteria for mixed mode I and mode II debond (and even includes mode III) were developed to predict the debond initiation and propagation of the cohesive layer. In Chapter III, an analytical solution was derived by introducing a correction term into the original Williams’ solution to predict the transverse stress in a cohesive layer when considering the deformation of a stiff substrate. Implementation of the cohesive layer model into a testbed finite element code was carried out and code verification was performed. Benchmark comparisons of finite element prediction of both global critical load and local stress field with analytical solution for a DCB specimen resulted in good agreement after modifications were made to the original Williams’ solution. A sensitivity study was conducted to evaluate the influence of cohesive layer thickness on local parameters such as damage zone length, and global parameter such as critical force. In Chapter IV, a twodimensional cohesive layer constitutive model involving strain dependent, nonFickian hygrothermal effects as well as diffusion induced degradation in the cohesive layer was constructed. Numerical simulation of a wedgetest 9 including debond growth caused by synergistic interactions between local stress and diffusing moisture was also presented to demonstrate the ability of the cohesive layer model to simulate environmental cracking. In Chapter V, a threedimensional cohesive layer model and corresponding mixed mode failure (debond) criterion were implemented in a testbed finite element code to simulate the full threedimensional peninsula blister test. Issues such as large deformation, timedependent material behavior, and residual stresses in the thin film were considered in the simulation model. Distinctive numerical techniques were successfully employed to simulate the unique liquid loading process. FEA simulation results were also compared with analytical solution and test data. Good agreement was obtained. In Chapter VI, cohesive layer model with strainrate dependent tractionseparation constitutive law was implemented in a testbed FEA code to simulate a moving wedge test. Timedependent material properties of the adhesive material were considered and quasistatic debond growth of the adhesive layer was successfully simulated by this code. Results predicted by the computational model were benchmarked through comparison with analytical solutions and mixed mode fracture tests. In Chapter VII, cohesive layer model was used to study the dynamic response of a FRP bonded concrete beam under blast loading. Implicit HilberHughesTaylor (HHT) method was employed in the model to allow better control of numerical damping. Long term and short term responses were obtained and their effects on the failure of the adhesive layer were investigated. Dynamic responses of the structure with an initial crack and its effect on debond initiation were also studied. 10 In Chapter VIII, a twodimensional implicit dynamic finite element formulation including material and geometric nonlinearity was derived and implemented into a testbed FEA code. Model verification under very large deformation was successfully performed through comparison with ABAQUS FEA predictions. Subsequently, the NOVA3D FEA model was applied to a circular steel plate with a polymer coating subjected to intensive blast loading, and the effect of polymer coating on the nonlinear dynamic response was numerically investigated. In Chapter IX, conclusions are presented based on the study of the cohesive layer model and its applications on various engineering structures. 11 CHAPTER II COHESIVE LAYER MODEL This chapter gives a detailed description of the cohesive layer model for twodimensional and threedimensional cases. It includes the definition of cohesive layer, the constitutive laws for the cohesive layer, the concept of workofseparation, and the mixed mode failure criteria. 2.1 Cohesive Layer Configuration Ductile polymeric adhesive materials usually have a nonlinear normal and shear stressstrain response. In the event of crack initiation and propagation in such polymeric materials, there exists a damage zone ahead of the debond tip, in which, craze and void initiation, growth and coalescence take place. The cohesive forces in the damage region must be taken into account to capture the behavior of the failing material in this zone, especially if the zone size is not sufficiently small compared to characteristic structural dimensions. The vicinity of debond tip can be divided into three common zones: an elastic/viscoelastic zone, a damage zone, and a debonded zone, as depicted in Fig. 21(a). The corresponding stressstrain (or tractionseparation) relation may take different nonlinear forms in each zone for different materials and different loading conditions. 12 (a) Debond tip traction force 0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Normalized strain Normalized stress (b) Tractionseparation law near debond tip Fig. 2.1 Three characteristic zones near debond tip and corresponding stressstrain relations In order to model this nonlinear behavior near the debond tip, a cohesive layer, which is a thin layer of cohesive finite elements, can be placed between two adjacent layers for laminated structures and multilayer structures, or along the bimaterial interface, or along a predicted cracking path in a single material (e.g. concrete and metal). It allows debond (or failure) to initiate and grow in these elements and different stressr z Elastic zone Damage zone Debonded zone 13 strain tractionseparation laws can be selected to cope with any particular nonlinear behavior in front of the debond tip. An example of the tractionseparation law for a cohesive layer is shown schematically in Fig. 21(b). The thickness, c h , of cohesive layer is an important parameter in the cohesive layer model. It is not arbitrary, but is directly related to a characteristic length parameter (δ ), such as crack opening displacement (COD). In this simulation, δ ε max c = h and max ε is the maximum strain that could be reached at debond tip. Environmental degradation in cohesive material is included in the cohesive layer model through the change of the material properties max σ and max ε . When investigating the moisture induced degradation, twodimensional and threedimensional moisture diffusion in the cohesive layer can be directly simulated. Mixed mode failure and corresponding failure criteria can also be easily established to predict the debond process. 2.2 Cohesive Layer Constitutive Equations 2.2.1 Triangular StressStrain TractionSeparation Law Triangular stressstrain tractionseparation law is a simple and commonly used model for cohesive material (Fig. 22), especially in theoretical analysis. Considering the stress and strain of mode I debond (opening mode) in the direction perpendicular to the debond surface, three types of zone are defined as follow: Elastic zone: when the transverse strain max 3 ε ≤ ε , stress linearly increases with strain, stress reaches its maximum value max σ =σ at max 3 ε = ε . 14 Damage zone: when the transverse strain max 3 ε > ε , stress decreases gradually from its maximum to zero as strain approaches max ε . Debond (failure) zone: when the transverse strain max ε >ε , stress remains zero which implies a full debond or separate of the cohesive layer. Fig. 22 Normalized triangular stressstrain tractionseparation law for a cohesive layer The maximum stress max σ and maximum strain max ε , which are material properties, are functions of time, strain rate, temperature and moisture concentration etc., and represent the prescribed maximum stress and strain that could be reached in cohesive layer when the cohesive layer debonds along a specified direction. 2.2.2 2D Cubic StressStrain TractionSeparation Law Based on fundamental principles of continuum mechanics, for twodimensional case, a more accurate cohesive layer constitutive relationship takes the cubic form as 15 employed by Needleman [7] (Fig. 23). In the direction perpendicular to the debond surface, the transverse normal stress is given by, ⎪ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎪ ⎨ ⎧ = < > − + ≤ ≤ = 0 4 27 4 27 0 1 ( 2 ) 0 1 4 27 22 max 22 22 max max 22 22 3 22 2 max 22 22 22 ε σ ε ε ε σ ε σ ε ε ε ε σ (21) 1.4 1.2 1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 1 1.2 0.4 0.2 0 0.2 0.4 0.6 0.8 1 1.2 Fig. 23 Cubic stressstrain tractionseparation law for a cohesive layer In the direction of debond growth, linear elastic response is assumed, and the axial normal stress is given by, 11 max 11 max max 11 ε σ ε ε σ σ = = (22) where, the normalized strain in a given direction is defined as max ε ε ε ii ii = (no sum on i). σ / σmax ε / εmax 16 A similar treatment to the one developed in Eq. (21) is prescribed for the shear 12 σ response, with the proviso that the shear stress is independent of the sign of shear strain. 2.2.3 3D Cubic StressStrain TractionSeparation Law Threedimensional cohesive layer treatment is often necessary for some structures like peninsula blister specimen, which is used to measure the interfacial fracture toughness in thin film structures. A full threedimensional cohesive layer model was constructed in this study to meet this requirement. By extending the above twodimensional cohesive layer model just described, the constitutive law for a threedimensional cohesive layer can be expressed in a similar manner as shown in Eqs. (23) and (24). Again, nonlinear responses are considered for transverse stress components 31 σ , 32 σ , and 33 σ , while for other stress components linear elastic response is assumed. 2 3 max 3 3 3 3 3 3 max 3 3 27 ( 2 ) 0 1 4 0 1 ( 1, 2, 3) 27 0 4 i i i i i i i i i σ ε ε ε ε σ ε σ ε ε ⎧ − + ≤ ≤ ⎪⎪⎪⎪ = = > ⎨⎪⎪⎪ < ⎪⎩ (23) max max max jk jk jk σ σ ε σ ε ε = = ( j, k =1, 2 ) (24) 17 2.3 Cohesive Layer WorkofSeparation For a given stressstrain tractionseparation law, the workofseparation (or strain energy) is the work needed to fully separate a unit area of cohesive layer, which is given by the total area under the prescribed stressstrain curve. Under the triangular stressstrain tractionseparation law as shown in Fig. 22, for pure mode I debond, c c I sep h d h max max 0 2 ( ) 1 max φ σ ε ε σ ε σ = ∫ = (25) and for pure mode II debond, c c II sep h d h max max 0 2 ( ) 1 max φ τ γ γ τ γ γ = ∫ = (26) Similarly, under the cubic stressstrain tractionseparation law as shown in Fig. 2 3, c c I sep h d h max max 0 16 ( ) 9 max φ σ ε ε σ ε σ = ∫ = (27) c c II sep h d h max max 0 16 ( ) 9 max φ τ γ γ τ γ γ = ∫ = (28) where c h is the thickness of the cohesive layer. 2.4 Failure Criteria of Mixed Mode Debond Mixed mode I and mode II debond is the common failure form of cohesive layer, while pure mode I or mode II debond is only a special case under certain conditions. It is 18 necessary to establish a failure criterion such that it contains both the contributions of mode I and mode II debond, and in some cases even includes mode III. 2.4.1 Criterion Based on Prescribed Maximum Strain At the debond tip, when the strains satisfy the following condition, cohesive layer will debond [21], δ γ γ α ε ε α = ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ + ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ 2 max 2 2 max 1 y xy (29) where y ε and xy γ are the transverse normal strain and shear strain respectively, max ε and max γ are the prescribed normal and shear failure strains of the cohesive layer, respectively. 1 α , 2 α , and δ are constants and 1 1 2 α =α =δ = was taken in this study. 2.4.2 Criterion Based on Strain Energy Release Rate The strain energy release rate in the cohesive layer during mixed mode debond, G , due to the tractionseparation force can be partitioned into the opening (mode I) and shear (mode II) components, I G and II G respectively, in such a way that, I II G = G +G (210) Each individual component can be calculated by integrating the mode I and II tractionseparation curves (Fig. 24) = ∫ t G d I ε σ ε ε 0 ( ) (211) 19 = ∫ t G d II γ τ γ γ 0 ( ) (212) where max σ σ = σ , max τ τ = τ , max ε ε = ε , max γ γ = γ are normalized stresses and strains in the specific directions, respectively. 1.0 0.8 0.6 0.4 0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0.2 0 0.2 0.4 0.6 0.8 1 1.2 Normalized strain Normalized stress GI Fig. 24 Cubic stressstrain tractionseparation law and strain energy release rate for a cohesive layer Considering the energy required to separate the cohesive layer, the cohesive layer debond process can be better predicted by means of the following criterion e G G G G n II II m I I C C = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ (213) where, I G and II G are the respective values of the ambient energy release rates given by the area under the corresponding stressstrain curves under a given applied loading. 20 IC G and IIC G are the critical strain energy release rates in pure mode I and mode II debond, respectively. m, n, and e are material constants, Kutlu and Chang [22] found that m = n = e = 1 provided the best fit to their experimental data. When neglecting energy dissipation in the bulk adhesive, critical strain energy release rate of cohesive layer is equal to the workofseparation sep φ (the total area under mode I or mode II tractionseparation curves for 0 ≤ε ≤ 1, Fig. 24), which is the energy necessary to generate unit debond length (2D case) or area (3D case). Thus for the cubic stressstrain tractionseparation law, integrating Eq. (23) over the limits [0, ε ] for the case of pure mode I debond and [0, γ ] for the case of pure mode II debond, ) 4 1 3 2 2 (1 4 ( ) 27 2 3 4 max max 0 σ ε ε σ ε ε ε ε ε I = c ∫ = c − + G h d h (214) ) 4 1 3 2 2 (1 4 ( ) 27 2 3 4 max max 0 τ γ γ τ γ γ γ γ γ II = c ∫ = c − + G h d h (215) And therefore, I Ic G →G as ε →1, and II IIc G →G as γ →1 giving I c c G h d h C max max 0 16 ( ) 9 max σ ε ε σ ε ε = ∫ = (216) II c c G h d h C max max 0 16 ( ) 9 max τ γ γ τ γ γ = ∫ = (217) For triangular stressstrain tractionseparation law, a similar procedure can be applied. Where, c h is the cohesive layer thickness. It should be noted that the cohesive layer thickness is not arbitrary in cohesive layer model, and it is related to some 21 characteristic lengthscale of the debond process, such as COD. Determination of cohesive layer thickness from test data will be discussed in a later section. A phase angle is defined to describe the modemix of the failure (debond) in a cohesive layer ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ = − 2 1 * * tan 1 I II G G φ (218) where * I G and * II G are the strain energy release rates at which debond of cohesive layer initiates. In the case of the presence of antiplane shear stress (mode III), the mixed mode debond criterion can be expressed as + + = 1 IIIc III IIc II Ic I G G G G G G (219) and the phase angle between mode I and mode III is given by, ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ = − 2 1 * * tan 1 I III G G ϕ (220) 22 CHAPTER III DEVELOPMENT OF AN ANALYTICAL SOLUTION FOR COHESIVE LAYER MODEL AND MODEL VERIFICATION The objective of this chapter is to find an analytical solution to the cohesive damage zone at the interface between a fiber reinforced polymer (FRP) plate and concrete substrate. An analytical solution was derived to predict the stress in cohesive layer when considering the deformation in a stiff substrate. A twodimensional cohesive layer model with a prescribed stressstrain tractionseparation law as described in Chapter II was employed in this study. For comparison purpose, the cohesive layer model was implemented into a testbed finite element code (NOVA3D). Detailed benchmark comparisons of analytical results with finite element predictions for a double cantilever beam specimen for model verification were performed and issues related to cohesive layer thickness were investigated. It was observed that the assumption of a rigid substrate in analytical modeling can lead to inaccurate analytical prediction of cohesive damage zone length and stress distribution near debond tip. 3.1 Analytical Solution from Cohesive Zone Model Williams and Hadavinia [23] used a cohesive zone model with various simple forms of cohesive tractionseparation laws to analyze the global features and local stress 23 distribution of a double cantilever beam specimen (DCB, shown in Fig. 31). The DCB specimen is modeled as a cantilever beam with elastic foundation as shown in Fig. 32. The deformation of the beam is given by the equation EI w dx d v = 4 4 (31) y x P P Concrete substrate Cohesive layer Fig. 31 A cohesive layer in a double cantilever beam (DCB) a l Crack Damaged Zone Elastic Zone x1 x2 v Fig. 32 Schematic debonding of a double cantilever beam showing cohesive layer with elastic zone and damage zone ( 1 x ∈[0,l]) P 24 Where E is the Young’s modulus of the beam, v is the deflection of beam at its midplane, b is the width, and w is the distributed load per unit length of the beam that can be related to the stress in the cohesive layer by w = −bσ y . The stress y σ in the cohesive zone is modeled by a triangular elasticlineardamage tractionseparation law referring to the stress in the cohesive zone (as depicted in Fig. 22), where max v is the displacement at final fracture, max σ is the maximum stress in the cohesive zone at 3 max v v = . (a) In the damage zone: ( ) 2 3 max max max v v v b w b y = − = − σ σ (32) Let max 4 max 1 2 3 EIv bσ λ = (33) Thus Eq. (31) becomes ( ) max 4 4 1 4 v v dx d v = λ − (34) and the corresponding solution is, max 1 1 1 2 1 1 1 1 1 2 1 1 v = v + B sinhλ x + B coshλ x + C sinλ x + C cosλ x (35) (b) In the elastic zone: v v w b b y max max 3 σ = − σ = − (36) Let 25 max 4 max 2 4 3 EIv bσ λ = (37) v dx d v 4 4 2 4 = −4λ (38) and the solution is given by, ( sin cos ) 1 2 2 2 2 2 v = e−λ2x2 A λ x + A λ x (39) There are a total six unknown coefficients in Eqs. (35, 39). Along with the unknown damage length l and the critical force cr P , it requires eight boundary conditions to determine the beam deformation and the corresponding stresses in the cohesive layer. Two force boundary conditions are provided by force and moment equilibrium at the crack tip 0 1 x = . Continuity conditions imposed at the boundary of the elastic zone and damage zone yield another four boundary conditions by matching v , dx dv , 2 2 dx d v , 3 3 dx d v at x = l 1 and 0 2 x = . Finally, two displacement boundary conditions are max v = v at 0 1 x = and 3 max v v = at x = l 1 . From 1 max v (0) = v , 0 2 2 B + C = (310) EI P a dx d v (0) = cr 2 1 2 , where a is the crack length 2 1 2 2 2EIλ B + C = Pcra (311) EI P dx d v (0) = cr 3 1 3 26 3 1 1 1 2EIλ B − C = Pcr (312) 3 ( ) max 1 v v l = 3 sinh cosh sin cos max max 1 1 2 1 1 1 2 1 v v + B λ l + B λ l + C λ l + C λ l = (313) 3 (0) max 2 v v = 3 max 2 A = v (314) 1 ( ) 2 (0) dx dv l dx dv = ( cosh sinh cos sin ) ( ) 1 1 1 2 1 1 1 2 1 2 1 2 λ B λ l + B λ l + C λ l − C λ l = λ A − A (315) ( ) (0) 2 2 2 2 1 2 dx d v l dx d v = 1 2 1 1 2 1 1 1 2 1 2 2 1 λ (B sinhλ l + B coshλ l − C sinλ l −C cosλ l) = −2λ A (316) ( ) (0) 3 2 3 3 1 3 dx d v l dx d v = ( cosh sinh cos sin ) 2 ( ) 1 2 3 1 1 2 1 1 1 2 1 2 3 1 λ B λ l + B λ l − C λ l + C λ l = λ A + A (317) Eqs. (310 ~ 317) are nonlinear in terms of damage length l. Solutions are sought by an iterative numerical predictorcorrector method as follows: Damage length l is varied from 0.0 to 1.0 mm with interval of 0.001~0.01mm for the current specimen, and thereby Eqs. (310 ~ 316) become a set of linear equations that can be solved to obtain the constants coefficients 1 A , 2 A , 1 B , 2 B , 1 C , 2 C and the crack initiation load cr P for a specified damage length l. These constants and corresponding damage length l are then substituted into Eq. (317) and, because the equation is not exactly solved, the solution 27 error is numerically estimated. The correct solution for damage length l and the constants is the one that minimizes the error. 3.2 Analytical Solution from Cohesive Layer Model In Williams’ solution, deformation of the beam transverse to beam axis is neglected, therefore the displacement in ydirection (v) at the centerline of the beam is considered as the deformation of the cohesive zone. Maximum deformation max v at the debond tip, which is independent of the geometry and material properties of the beam, is the characteristic length scale of the cohesive zone model. When applying cohesive zone approach to a cohesive layer model, maximum deformation can be expressed as max max ε c v = h , where max ε is the maximum strain in the cohesive layer at failure, which is assumed to be a material property. As a result, the cohesive layer thickness c h is no longer arbitrary, but is uniquely determined by the relation max max ε h v c = . When considering the transverse deformation of the beam, c h is an important factor in evaluating the relative stiffness of beam and the cohesive layer in the transverse direction. The Young’s modulus of the beam (concrete) is usually much higher than that of the cohesive layer. On the other hand, the thickness of the beam is also much greater than that of the cohesive layer. Consequently the transverse deformation of the beam is comparable to the deformation in the cohesive layer and thus cannot be neglected. A small lateral deformation of the beam will greatly change the stress in the cohesive layer. FEA results clearly show the difference between the deflection of the beam at the 28 centerline and the deformation of the cohesive layer (that is, the displacement at the beamcohesive layer interface represents the deformation of the cohesive layer) (Fig. 3 3). As can be seen in Fig. 33 the displacement at the interface is generally smaller than the displacement at the centerline of the beam, the latter being the summation of the transverse deformations of the beam and the cohesive layer. Displacement ( mm ) X coordinate ( mm ) 4.0E04 2.0E04 0.0E+00 2.0E04 4.0E04 6.0E04 8.0E04 1.0E03 1.2E03 0 0.5 1 1.5 2 2.5 3 Interface Beam center line Fig. 33 Transverse deformation comparison between beam centerline and cohesive interface using FEA Thus, the total displacement v in ydirection at the centerline of the beam is composed of two parts: deformation b v of the beam and deformation c v of the cohesive layer, such that b c v = v + v (318) Assume that the transverse stress at the interface and in the cohesive layer is y c σ =σ , and = 0 y σ at the free top surface of the beam. A distribution law for y σ 29 through the height of the beam ( b h ) must be assumed to calculate the deformation of the beam (see Fig. 34). The actual distribution law for transverse stress can be obtained from elastic FEA analysis. Three idealized distribution laws for y σ were evaluated in this study: beam cohesive layer hc hb y stres s σc hb 0 Fig. 34 Transverse stress distribution in the beam and cohesive layer (a) Linear: (h y) h b b c y = − σ σ (b) Quadratic: 2 2 (h y) h b b c y = − σ σ (319) (c) Cubic: 3 3 (h y) h b b c y = − σ σ The corresponding lateral deformations in the lower half of the beam are, respectively: (a) Linear: b h h c b b b b c b h y b y E h y dy h E h dy E v dy b b σ b σ σ ε 8 ( ) 3 2 0 2 0 2 0 = ∫ = ∫ = ∫ − = (b) Quadratic: b h h c b b b b c b h y b y E h y dy h E h dy E v dy b b σ b σ σ ε 24 ( ) 7 2 0 2 0 2 2 2 0 = ∫ = ∫ = ∫ − = (320) 30 (c) Cubic: b h h c b b b b c b h y b y E h y dy h E h dy E v dy b b σ b σ σ ε 64 ( ) 15 2 0 2 0 3 3 2 0 = ∫ = ∫ = ∫ − = where b E is the Young’s modulus of the beam, and b h is the height of the beam. The above results can be generalized as b c b b E h v k σ = , which represents the lateral deformation of the beam, where k is the coefficient determined by the distribution of y σ along the ydirection within the beam given by Eq. (320). From Eq. (320), the values of k are 8 3 , 24 7 , 64 15 for linear, quadratic and cubic distributions, respectively. The analytical solution to the DCB specimen bonded by a cohesive layer with small but finite thickness can now be derived as follows. (a) In the damage zone: From Eq. (32), max max max 3 2 σ σ v v v c y = − (321) Displacement at the centerline of the beam is ⎟ ⎟⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = + = − − b b b c y E v v v v v k h max max max 3 2 σ σ (322) ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ − = − = − b b b b b b y E E I v k h v v b E I b dx d v max max 4 max 4 3 2 ( ) σ σ (323) By defining ⎟ ⎟⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − = b b b b E E I v k h b max max 4 1 3 2 σ λ (324) 31 The governing equation becomes ( ) max 4 4 1 4 v v dx d v =λ − (325) (b) In the elastic zone: ⎟ ⎟⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = + = + b b b c y E v v v v k h max max 3σ σ (326) ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ + = − = − b b b b b y E E I v k h v b E I b dx d v max max 4 4 3 4 4 σ σ (327) Defining, ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ + = b b b b E h k v E I b max max 4 2 3 4 σ λ (328) v dx d v 4 4 2 4 = −4λ (329) Governing Eqs. (325, 329) are of a similar form with the original Williams’ Eqs. (34, 38). But the coefficients 1 λ and 2 λ are different because of the presence of the extra term b b E k h which represents the transverse deformation of the beam (note that k = 0 gives the original Williams’ solution). A solution procedure and boundary conditions similar to the ones used for solving Eqs. (34, 38) can be employed to solve Eqs. (325, 329) to determine the unknown constants. Regarding the boundary conditions for Eqs. (325, 329), the two force boundary conditions and four displacement continuity conditions are same as the original solution 32 because these conditions are directly related to the beam. But the two displacement conditions, which are related to the cohesive layer, are a little different. First, at x = l 1 and 0 2 x = , the deformation of cohesive layer c v (not v at the centerline of the beam) equals 3 max v and the corresponding stress is max σ . 3 max v v v vc b = − = (330) Therefore the corresponding Eqs. (313, 314) become 3 ( ) max 1 v v l c = 3 ( sinh cosh sin cos ) max max max 1 1 2 1 1 1 2 1 v E h v B l B l C l C l k b + + + + − b = σ λ λ λ λ (331) 3 (0) max 2 v v c = 3 max max 2 v E h A k b − b = σ (332) On the other hand, at 0 1 x = , transverse stress = 0 y σ in the concrete beam, so that the lateral deformation of beam = 0 b v , therefore max v v v v v b c c = + = = which is same as Eq. (310). 3.3 Comparison between Analytical Solution and FEA Results A double cantilever beam (DCB) consisting of two substrates bonded with a thin layer of epoxy adhesive is shown in Fig. 31. The entire layer of epoxy is modeled with 33 special cohesive layer elements which obey the prescribed stressstrain tractionseparation law. A cohesive layer thickness h mm c = 0.02 and beam height h mm b =1 were used in this example, and the initial debond length a =15 mm corresponds to the unbonded portion of the beam. Unlike the usual linear elastic fracture mechanics modeling method of a sharp crack tip of zero tip radius, the localized fracture process zone in the current study has a small but finite thickness c h . Substrate beam is modeled as linear elastic, and the material properties used in this analysis are listed in Table 31. Table 31 Material properties for concrete and epoxy adhesive Concrete Epoxy Young’s modulus ( GPa ) 27.5 3.85 Shear modulus ( GPa ) 11.0 1.54 Poisson’s ratio 0.25 0.25 max σ ( MPa )  30.0 max ε  0.0526 Fig. 35 Finite element mesh of a DCB specimen with symmetry boundary conditions 34 Considering the symmetry boundary conditions along the centerline of epoxy adhesive layer, only half of the DCB specimen is modeled. Fig. 35 shows the actual finite element mesh (8node quadrilateral element) with applied symmetry boundary conditions. A very fine mesh is used to model the sharp stress gradients in the damage region. Plane strain conditions are assumed. Convergence study shows that converged results were obtained when the geometries of the elements in the crack tip area are smaller than 10 1 of the length of the damage zone. A constant beam tip displacement (instead of applied force) is specified in the FEA analysis to simulate a wedge test and to ensure that stable debond growth would occur. Debond growth (or cohesive layer element failure) is characterized by the transverse mechanical strain in the cohesive elements exceeding the specified maximum strain max ε in the cohesive layer beyond which the transverse stress goes to zero as defined by the cohesive constitutive law depicted in Fig. 22. Instead of node release and element deletion schemes used in most finite element codes, a failed element remains active in the subsequent analysis while the stiffness of the element is reduced close to zero. Analytical results with different k values are obtained and compared with numerical results from FEA as shown in Fig. 36. It can be seen that the best agreement with the FEA results is reached with the correction coefficient 64 k = 15 , which represents a cubic distribution of transverse stress in the beam. This cubic stress distribution is also verified by FEA simulation of the stress field in the concrete beam. The length of damage zone l decreases as the k increases, which implies a greater transverse deformation in the 35 beam under the same stress at the interface of beam and the cohesive layer. Due to the complexity of the transverse stress distribution, especially near the crack tip, k may take different values under different combination of material properties and geometries of the concrete beam and the cohesive layer. Further, as shown in Fig. 37, the deviation of the FEA results from the analytical solutions within the damage zone is likely attributable to the use of eightnode quadrilateral elements with quadratic interpolation that could result in linear variation of throughthickness strain within the fracture localization zone. Fortunately, the stress distribution in the damage zone does not have a significant effect on the process of debond initiation and propagation in the cohesive layer. 2.E+07 1.E+07 5.E+06 0.E+00 5.E+06 1.E+07 2.E+07 2.E+07 3.E+07 3.E+07 4.E+07 0 0.5 1 1.5 2 Unmodified k=0 Modified k=15/64 Modified k= 7/24 Modified k= 3/8 NOVA3D Fig. 36 Stress distribution with different correction coefficient k X coordinate (mm) = 0.02 b c h h Transverse stress (Pa) 36 Furthermore, different cohesive layer thickness and debond lengths were used to verify the agreement between analytical solutions and FEA simulation results, for both global and local metrics as discussed in the following paragraphs. Fig. 37 Reaction force comparison for different normalized cohesive layer thicknesses To study the effect of cohesive layer thickness on a global metric such as critical force (Pcr), the comparisons are conducted under the assumption that the cohesive layer deformation at the debond tip is equal to the maximum strain max ε of the cohesive layer (critical strain). For most practical cases, the ratio of cohesive layer thickness to beam thickness ( b c h h ) is in the range of 0.02 ~ 0.15. Critical reaction force versus normalized cohesive layer thickness ( b c h h ) is studied and good agreement is obtained between 37 analytical and FEA solutions. As shown in Fig. 37, the critical reaction force by FEA is slightly lower than the analytical predictions. This error is probably due to the fact that at the crack tip ( x1 = 0 ) the transverse displacement at the interface is actually slightly greater than that at the beam centerline (see Fig. 33), thereby reducing the reaction force required. To study the effect of debond length on reaction force, two kinds of boundary conditions are used: (a) Deformation control: under the same critical deformation at debond tip for FEA and analytical solution (b) Displacement control: under the same free end displacement for FEA and analytical solution. Good agreement is observed in both cases as shown in Fig. 38. 200 300 400 500 600 700 800 900 10 15 20 25 Debond length ( mm ) Reaction force ( N ) Analytical Deformation control (FEM) Displacement control (FEM) Fig. 38 Reaction force comparison for different debond lengths = 0.02 b c h h 38 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 0.05 0.1 0.15 0.2 0.25 Williams Modified Williams FEM Fig. 39 Damage length comparison for different cohesive layer thicknesses For local stress verifications, cohesive layer stress distribution in crack tip region is determined by evaluating damage length l. It can be observed from Fig. 39 and from Table 32 that the results from FEA are in good agreement with the modified Williams’ analytical solution when the normalized thickness lies between 0.02 and 0.15. When > 0.20 b c h h , that is, when the thickness of cohesive layer is relatively large, the deformation of the beam is relatively small compared with the deformation of the cohesive layer. As a result, the influence of the correction factor is no longer significant and the FEA solution approaches the original (unmodified) Williams’ solution. Fig. 310 shows how cohesive layer thickness c h influences local transverse stress distribution near crack tip as characterized by the damage length l. Good agreement is observed between analytical prediction and FEA results. Normalized thickness hc/hb Damage length ( mm ) 39 Fig. 310 Stress distributions for different cohesive layer thicknesses Table 32 Variation of damage length with different cohesive layer thicknesses Damage length ( mm ) b c h h Analytical NOVA3D Error ( % ) 0.02 0.186 0.190 2.15 0.05 0.336 0.332 1.19 0.10 0.446 0.437 2.02 0.15 0.513 0.534 4.09 0.20 0.562 0.587 4.45 It was observed that the modification to Williams’ model has little effect on global metrics (i.e. critical force and free end displacement), while for localized damage zone length and local stress distribution the effect of different values of k is significant as 40 depicted in Fig. 36, Fig. 39 and Table 33. It can be concluded that value of k has significant effect on the damage length, and 64 k = 15 results in the best agreement with FEA for the cohesive layer model with normalized thickness ranging from 0.02 to 0.15. For smaller or larger cohesive layer thickness, the relative lateral deformation of cohesive layer and concrete beam will change significantly. Further analytical studies are necessary to address these specific conditions. Table 33 Analytical and numerical solutions of damage length Correction factor k 0 15/64 7/24 3/8 FEA Damage length l ( mm ) 0.328 0.186 0.156 0.116 0.190 3.4 Conclusions An analytical solution was derived by introducing a correction term into the original Williams’ solution to predict the transverse stress in a cohesive layer when considering the deformation in a stiff substrate. Implementation of the cohesive layer model into a testbed finite element code was carried out and code verification was performed. Benchmark comparisons of finite element prediction of both global critical load and local stress field with analytical results for a DCB specimen resulted in good agreement after modifications were made to the original Williams’ solution. A sensitivity study was conducted to evaluate the influence of cohesive layer thickness on local parameters such as damage zone length, and global parameter such as critical force. From 41 the present studies it can be concluded that both local and global cohesive layer parameters are fairly sensitive to the cohesive layer thickness, whereas the correction factor (k) to Williams’ original solution significantly influences the local stress distribution and damage length. 42 CHAPTER IV HYGROTHERMAL EFFECT ON COHESIVE LAYER The objective of this chapter is to model the synergistic bond degradation mechanism that may occur at the interface between a fiber reinforced polymer (FRP) and concrete. For this purpose, a twodimensional cohesive layer involving straindependent, nonFickian hygrothermal effects as well as diffusion induced degradation in the cohesive layer was constructed. The model was implemented in a testbed finite element code (NOVA3D). Results from demonstration cases involving synergistic bond degradation were presented. 4.1 Introduction The influence of moisture diffusion on crack growth along an interface is not yet fully understood. Environmental cracking in a polymer typically occurs in the presence of a penetrant, such as moisture, and stress. It has been postulated that the mechanism involved in environmental crack growth in a polymer involves a small zone of craze formation and/or plasticization at the crack tip. For thermoset resins, such as epoxy, energy absorption at the crack tip is primarily by a shear yielding process and not by crazing. Consequently, for a thermoset epoxy, the zone of plasticization ahead of the crack tip must be determined using a diffusion law for nonporous media, such as Fick’s 43 law. However, quite frequently, polymer composites exhibit deviations from the classical Fickian treatment, termed as anomalous or nonFickian diffusion, especially at elevated temperatures and stress levels, and at high relative humidity. Some researchers have suggested that the deviation can be explained by a twostage Fickian process [24, 25]. Others claim that the diffusion process in a polymer is really nonFickian [26, 27]. Sophisticated hygrothermal models have been developed and verified by Roy [1720] to account for anomalous diffusion. The objective of this chapter is to model the synergistic bond degradation mechanisms that might occur at the interface due to interactions between stress, cohesive damage, and penetrant diffusion. For this purpose, a twodimensional cohesive layer constitutive model with a prescribed tractionseparation law is constructed from basic principles of continuum mechanics and thermodynamics, taking into account nonFickian hygrothermal effects that are likely to occur within the cohesive layer. 4.2 Fickian Diffusion Diffusion is the process by which matter is transported from one part of a system to another as a result of random molecular motions. In 1855, Fick first put diffusion on a quantitative basis by adopting the mathematical equation of heat conduction. The theory of diffusion in isotropic substance is therefore based on the hypothesis that the rate of transfer of diffusing substance through unit area of a section is proportional to the concentration gradient measured normal to the section, i.e. x f D C ∂ = − ∂ (41) 44 where C is the concentration of diffusing substance, D is called diffusion coefficient or diffusivity, and is the function of coordinate x, y, and z (location) and concentration C. Considering the mass balance of diffusing substance, we have = 0 ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ z f y f x f t C x y z (42) ⎟⎠ ⎞ ⎜⎝ ⎛ ∂ ∂ ∂ + ∂ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ ∂ + ∂ ⎟⎠ ⎞ ⎜⎝ ⎛ ∂ ∂ ∂ = ∂ ∂ ∂ z D C y z D C x y D C t x C (43) The fundamental differential equation of diffusion in an isotropic medium (where D is independent of the concentration C and location) can be expressed as ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ∂ + ∂ ∂ + ∂ ∂ = ∂ ∂ ∂ 2 2 2 2 2 2 z C y C x D C t C (44) and this is the socalled Fickian diffusion. 4.3 NonFickian Diffusion (Strain Assisted Diffusion) For a twodimensional cohesive layer of finite thickness c h , under planestrain conditions as shown in Fig. 41, the Helmholtz free energy per unit volume is given by, ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 4 18 12 3 17 22 12 2 12 2 16 22 12 3 15 22 4 14 22 3 13 12 2 12 12 22 12 2 11 22 3 7 11 12 8 11 22 9 12 22 10 22 2 6 12 2 5 22 2 0 1 11 2 22 3 12 4 11 , , , , , , , , , , , , , , , , , , , ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ρψ ε ε ε ε ε C m T C m T C m T C m T C m T C m T C m T C m T C m T C m T C m T C m T C m T C m T C m T C m T C m T C m T C m T + + + + + + + + + + + + + = + + + + + (45) where the mechanical strain components are defined as, 45 ( ) ( ) 11 11 REF REF ε = E −α T − T −β m − m ( ) ( ) 22 22 REF REF ε = E −α T −T −β m − m 12 12 ε = E and, ρ : mass density of material in the cohesive layer 11 ε : mechanical strain component in x direction 22 ε : mechanical strain component normal to crack face (in y direction) 12 ε : shear strain component tangential to crack face ij E : total (kinematic) strain components m : moisture concentration in the cohesive layer at time t REF m : reference moisture concentration T : temperature in the cohesive layer at time t REF T : reference temperature α (T) : linear coefficient of thermal expansion β (T) : linear coefficient of moisture expansion moisture diffusion y x cohesive layer concrete beam Fig. 41. A cohesive layer with moisture diffusion in a DCB beam 46 From reduced entropy inequality, the cohesive stresses are defined by, ( , ) 2 ( , ) ( , ) ( , ) 1 4 11 7 12 8 22 11 11 σ ρ ψ C m T C m T ε C m T ε C m T ε E = + + + ∂ = ∂ (46) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 22 2 5 22 8 11 9 12 10 22 22 2 3 11 22 12 12 12 14 22 2 2 3 15 22 12 16 22 12 17 12 , 2 , , 3 , 2 , , 4 , 3 , 2 , , C mT C C mT C mT C mT E C mT C mT C mT C mT C mT C mT σ ρ ψ ε ε ε ε ε ε ε ε ε ε ε ε ε = ∂ = + + + + ∂ + + + + + + (47) Similarly, ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )3 18 12 2 17 22 12 12 2 16 2 3 15 22 2 12 22 12 13 12 2 3 6 12 7 11 9 22 11 22 12 12 3 , 4 , 2 , 3 , , 2 , , 2 , , , , ε ε ε ε ε ε ε ε ε ε ε ε ε σ ρ ψ C m T C m T C m T C m T C m T C m T C m T C m T C m T C m T C m T E + + + + + + = + + + + ∂ = ∂ (48) Chemical potential of the diffusing vapor is defined by Weistman [28], ∂m μ = ρ ∂ψ or, ⎥⎦ ⎤ ⎢⎣ ⎡ − ∂ ∂ + ⎥⎦ ⎤ ⎢⎣ ⎡ − ∂ ∂ + ⎥⎦ ⎤ ⎢⎣ ⎡ − ∂ ∂ + ⎥⎦ ⎤ ⎢⎣ ⎡ − ∂ ∂ + ∂ ∂ = 3 14 22 4 22 2 14 10 22 3 22 10 5 22 2 22 5 4 11 2 11 0 4 3 ( , ) ( ) 4 ( , ) ( ) 2 ( , ) ( ) 2 ( , ) ( ) ε β ε ε β ε μ ε β ε ε β ε C m T T m C m T T C m C C m T T m C m T T C m C m C (49) From conservation of mass, the governing Equation for twodimensional moisture diffusion is, x y m f f t x y ∂ ⎛ ∂ ∂ ⎞ = −⎜ + ⎟ ∂ ⎝ ∂ ∂ ⎠ (410) 47 where the moisture flux, f = fxnˆx + fynˆy r , in the absence of temperature gradients is given by, ˆ ˆ x y f D x f D y μ μ = − ∂ ∂ = − ∂ ∂ (411) Assuming isotropic material and using the chain rule, 22 22 ˆ ˆ ˆ x f D m D T D m x T x x μ μ μ ε ε ⎛ ∂ ⎞ ∂ ⎛ ∂ ⎞ ∂ ⎛ ∂ ⎞ ∂ = − ⎜ ⎟ − ⎜ ⎟ − ⎜ ⎟ ⎝ ∂ ⎠ ∂ ⎝ ∂ ⎠ ∂ ⎝ ∂ ⎠ ∂ (412) 22 22 ˆ ˆ ˆ y f D m D T D m y T y y μ μ μ ε ε ⎛ ∂ ⎞ ∂ ⎛ ∂ ⎞ ∂ ⎛ ∂ ⎞ ∂ = − ⎜ ⎟ − ⎜ ⎟ − ⎜ ⎟ ⎝ ∂ ⎠ ∂ ⎝ ∂ ⎠ ∂ ⎝ ∂ ⎠ ∂ (413) Assuming isothermal condition and substituting Eqs. (49, 412, 413) into Eq. (4 10), gives 22 22 m m m D m D D m D t x x ε x y y ε y ∂ = ∂ ⎛ ∂ + ∂ε ⎞ + ∂ ⎛ ∂ + ∂ε ⎞ ∂ ∂ ⎜⎝ ∂ ∂ ⎟⎠ ∂ ⎜ ∂ ∂ ⎟ ⎝ ⎠ where ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ∂ = ∂ ⎟⎠ ⎞ ⎜⎝ ⎛ ∂ = ∂ 22 ˆ , ˆ ε μ μ ε and D D m D D m ⎥⎦ + ⎤ ∂ ∂ − ∂ ∂ + + ∂ ∂ − ∂ ∂ + + ∂ ∂ − ⎢⎣ ⎡ ∂ ∂ + + ∂ ∂ − ∂ ∂ + ∂ ∂ = 2 22 2 14 3 22 14 4 2 22 14 2 22 2 10 2 22 3 10 2 22 10 2 2 5 22 2 5 2 22 5 2 2 11 4 2 4 2 11 4 2 2 0 2 8 ( ) 6 ( ) ( ) 2 ( ) 6 ( ) 6 ( ) ( ) ˆ 4 ( ) 2 ( ) 4 ( ) β ε β ε β ε β ε β ε ε ε β ε β ε β ε T C T T m C m T C T T C m C m C T C T m C m C T C T m C m C m C D D m (414a) 48 ⎥⎦ − ⎤ ∂ ∂ + ⎢⎣ ⎡ − ∂ ∂ − + ∂ ∂ = 2 14 22 3 22 14 10 22 2 22 10 22 5 5 4 12 ( ) ( ) ˆ 2 2 ( ) ( ) 3 6 ( ) ( ) ε β ε ε β ε β ε ε T C T m C T C T m C T C T m C D D (414b) If the unknown material coefficients are expanded in a Taylor Series about a reference moisture concentration value REF m , and retaining terms up to second order in change in moisture concentration, ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 2 2 , , ˆ , 1 , , ( 0, 1,......,18) REF REF k k k k REF REF REF m m k REF k REF k REF C m T C m T C m m C m m m m C m T C m T m C m T m k ⎛ ∂ ⎞ ⎛ ∂ ⎞ = + ⎜ ⎟ − + ⎜ ⎟ − ⎝ ∂ ⎠ ⎝ ∂ ⎠ = ⎡⎣ + Δ + Δ ⎤⎦ = % (415) In order to benchmark the present model against an established cohesive zone model, some of the unknown material coefficients were determined by assuming a cubic tractionseparation law similar to the one proposed by Needleman [7] and modified by ElSayed and Sridharan [14] for a finitethickness cohesive layer, giving, ( ) ( )[ ( ) ( ) 2 ] 0 0 0 0 C m,T Cˆ m ,T 1 C m ,T m C~ m ,T m REF REF REF = + Δ + Δ ( , ) 0 1 C m T = ( , ) 0 2 C m T = ( , ) 0 3 C m T = ( ) ( )[ ( ) ( ) 2 ] 4 4 4 C m,T T 1 C m ,T m C~ m ,T m REF REF MAX = MAX + Δ + Δ ε σ ( ) ( )[ ( ) ( ) 2 ] 5 5 5 1 , ~ , 8 C m,T 27 T C m T m C m T m REF REF MAX = MAX + Δ + Δ ε σ 49 ( ) ( )[ ( ) ( ) 2 ] 6 6 6 1 , ~ , 8 C m,T 27 T C m T m C m T m REF REF MAX = MAX + Δ + Δ γ τ ( ) ( ) ( ) ⎪⎭ ⎪⎬ ⎫ = = = , 0 , 0 , 0 9 8 7 C m T C m T C m T assume ν = 0 ( ) ( )[ ( ) ( ) 2 ] 10 2 10 10 1 , ~ , 2 C m,T 9 T C m T m C m T m REF REF MAX = − MAX + Δ + Δ ε σ ( , ) 0 11 C m T = ( , ) 0 12 C m T = ( ) ( )[ ( ) ( ) 2 ] 13 2 13 13 1 , ~ , 2 C m,T 9 T C m T m C m T m REF REF MAX = − MAX + Δ + Δ γ τ ( ) ( )[ ( ) ( ) 2 ] 14 3 14 14 1 , ~ , 16 C m,T 27 T C m T m C m T m REF REF MAX = MAX + Δ + Δ ε σ ( , ) 0 15 C m T = ( , ) 0 16 C m T = ( , ) 0 17 C m T = ( ) ( )[ ( ) ( ) 2 ] 18 3 18 18 1 , ~ , 16 C m,T 27 T C m T m C m T m REF REF MAX = MAX + Δ + Δ γ τ Substituting these definitions in expressions for , , , 11 22 12 σ σ σ (Eqs. 46, 47, and 48) results in, ( ) ( ) ( ) MAX MAX C T E T T ε σ = = 4 2 ˆ gives ( )[ ( ) ( ) ] 11 2 11 4 4 σ σ T 1 C m ,T m C~ m ,T m ε MAX REF REF = + Δ + Δ (416) 50 Similarly, ( )[ ( ) ( ) ] ( )[ ( ) ( ) ] ( )[ ( ) ( ) ] 3 22 2 14 14 14 2 22 2 10 10 10 22 2 22 5 5 5 4 ˆ 1 , ~ , 3 ˆ 1 , ~ , 2 ˆ 1 , ~ , ε ε σ ε C T C m T m C m T m C T C m T m C m T m C T C m T m C m T m REF REF REF REF REF REF + + Δ + Δ + + Δ + Δ = + Δ + Δ or ( ){[ ( ) ( ) ] [ ( ) ( ) ] [ ( ) ( ) ] 3 } 22 2 14 14 2 22 2 10 10 22 2 22 5 5 1 , ~ , 21 , ~ , 1 , ~ , 4 27 ε ε σ σ ε C m T m C m T m C m T m C m T m T C m T m C m T m REF REF REF REF MAX REF REF + + Δ + Δ − + Δ + Δ = + Δ + Δ (417) where the maximum tensile stress in the cohesive layer is MAX σ and characteristic interface length δ is related to the layer thickness c h by, MAX c h ε = δ . And finally the shear stress can be expressed as, ( ){[ ( ) ( ) ] [ ( ) ( ) ] [ ( ) ( ) ] 3 } 12 2 18 18 2 12 2 13 13 12 2 12 6 6 1 , ~ , 21 , ~ , 1 , ~ , 4 27 ε ε σ τ ε C m T m C m T m C m T m C m T m T C m T m C m T m REF REF REF REF MAX REF REF + + Δ + Δ − + Δ + Δ = + Δ + Δ (418) where MAX MAX MAX γ ε ε ε ε ε ε ε ε 12 12 22 22 11 11 = , = , = , the maximum shear stress in the cohesive layer is MAX τ . The consistent diffusivities m D and ε D for the cohesive layer are obtained by substituting the definitions of the material coefficients into Eqs. (414a, 414b), 51 { ( ) ( )[ ] ( )[ ] ( )[ ] ( )[ ] ⎭ ⎬ ⎫ + − + Δ + + Δ + Δ − − + Δ + + Δ + Δ + − + Δ + + Δ + Δ + − + Δ + + Δ + Δ = 2 22 2 14 14 3 2 14 14 22 4 3 14 22 22 2 10 10 2 2 10 10 22 3 2 10 22 2 5 5 2 5 5 22 2 5 22 2 4 4 2 4 4 11 2 4 11 0 0 2 ~ 4 ( 2 ~ ) 2 (1 ~ ) 16 27 2 ~ 4 ( 2 ~ ) 2 (1 ~ ) 2 9 2 ~ 4 ( 2 ~ ) 2 (1 ~ ) 8 27 2 ~ 4 ( 2 ~ ) 2 (1 ~ ) ˆ 2 ˆ ~ ( ) ε β ε β ε ε σ ε β ε β ε ε σ ε β ε β ε σ ε β ε β ε σ T C C C m C m C m T C C C m C m C m T C C C m C m C m T C C C m C m C m D D C T C T MAX MAX MAX MAX MAX MAX MAX MAX m (419) ( )[ ] ( )[ ] ( )[ ] ⎭ ⎬ ⎫ + + Δ − + Δ + Δ − + Δ − + Δ + Δ ⎩ ⎨ ⎧ = + Δ − + Δ + Δ 2 22 2 14 14 3 3 14 14 22 22 2 10 10 2 2 10 10 22 2 5 5 22 5 5 ( 2 ~ ) 3 (1 ~ ) 4 27 ( 2 ~ ) 2 (1 ~ ) 2 27 ( 2 ~ ) (1 ~ ) 4 ˆ 27 ε β ε ε σ ε β ε ε σ ε β ε σ ε T C C m C m C m T C C m C m C m D D T C C m C m C m MAX MAX MAX MAX MAX MAX (420) 4.4 Cohesive Layer Diffusion Boundary Conditions Assuming that the chemical potential of the ambient vapor on the exposed boundary of the cohesive zone remains constant with respect to time [28], the resulting concentration at the boundary of the cohesive zone (crack tip) can be derived as, ( ) ij BOUNDARY b μ T ,m,ε = μ (421) 52 { ( )[ ] ( )[ ] ( )[ ] MAX BOUNDARY MAX MAX MAX MAX MAX b T C C m C m C m T C C m C m C m T C C m C m C m C C C m ⎭ ⎬ ⎫ + + Δ − + Δ + Δ − + Δ − + Δ + Δ + + Δ − + Δ + Δ = + Δ 3 22 2 14 14 4 3 14 14 22 2 22 2 10 10 3 2 10 10 22 22 2 5 5 2 5 5 22 0 0 0 ( 2 ~ ) 4 (1 ~ ) 16 27 ( 2 ~ ) 3 (1 ~ ) 2 9 ( 2 ~ ) 2 (1 ~ ) 8 27 ˆ ( 2 ~ ) ε β ε ε σ ε β ε ε σ ε β ε ε σ μ An equation about the unknown moisture concentration Δm at the boundary can be derived in the form of aΔm2 + bΔm + c = 0 . Therefore the boundary concentration can be solved a m m b b ac b REF 2 = + − ± 2 − 4 (422) 4.5 Cohesive Layer WorkofSeparation For pure Mode I or Mode II fracture, the workofseparation at the cohesive interface (of finite thickness c h ) per unit volume in the presence of moisture concentration (Δm ) is given by ( ) {[ ] [ ] [ ] } ( ) [ ] [ ] [ ] ⎭ ⎬ ⎫ − + Δ + Δ + + Δ + Δ ⎩ ⎨ ⎧ = + Δ + Δ ⎭ ⎬ ⎫ − + Δ + Δ + + Δ + Δ ⎩ ⎨ ⎧ = + Δ + Δ = ⎥ ⎥⎦ ⎤ ⎢ ⎢⎣ ⎡ = ∫ ∫ ∫ ∫ 2 10 10 14 14 2 5 5 22 3 22 2 14 14 2 10 10 22 1 0 22 2 5 5 1 0 22 22 0 22 22 1 ~ 4 1 ~ 1 3 2 1 ~ 2 1 4 27 21 ~ 1 ~ 1 ~ 4 27 max C m C m C m C m T h C m C m C m C m C m C m d h T C m C m d dV h d MAX MAX c c MAX MAX c MAX V I sep σ ε ε ε ε σ ε ε φ σ ε σ ε ε ε (423) 53 ( ) { } } ( ) max 1 12 12 12 12 0 0 1 2 6 6 12 0 2 2 2 3 13 13 12 18 18 12 12 2 6 6 13 13 ) 27 1 4 2 1 1 27 1 1 4 2 2 1 3 II sep c MAX V c MAX MAX MAX MAX c d dV h d h T Cm C m C m C m C m C m d T h C m C m C m C m γ φ σ ε σ γ ε τ γ ε ε ε ε τ γ ⎡ ⎤ = ⎢ ⎥ = ⎢⎣ ⎥⎦ = ⎧ ⎡ + Δ + Δ ⎤ ⎨ ⎣ ⎦ ⎩ − ⎡⎣ + Δ + Δ ⎤⎦ + ⎣ + Δ⎤ + ⎡ Δ ⎦ = ⎧ ⎡ + Δ + Δ ⎤ ⎨ ⎣ ⎦ ⎩ − ⎡⎣ + Δ + Δ ∫ ∫ ∫ ∫ % % % % % 2 18 18 1 1 4 C m C m ⎫ ⎤ ⎡ ⎤ + + Δ + Δ ⎬ ⎦ ⎣ ⎦⎭ % (424) or write in the short form [ 2 ] 1 2 ( ) 1 ( ) ( ) 16 9 T h A T m A T m MAX MAX c I sep φ = σ ε + Δ + Δ (425) [ 2 ] 3 4 ( ) 1 ( ) ( ) 16 9 T h A T m A T m MAX MAX c II sep φ = τ γ + Δ + Δ (426) The thickness of the cohesive layer, c h , is a measure of the fracture localization zone and it is directly related to the characteristic length scale of the fracture process zone controlled by specific fracture mechanism [29]. The area under the normalized stressstrain curve ( max 0 ε 1 ε ≤ ≤ ) shown in Fig. 42 represents the workofseparation of the cohesive layer. From Eqs. (425, 426) and experimental results, it can be observed that the workofseparation decreases at elevated temperature and higher moisture concentration due to the physical or chemical degradation at the interface through a decrease in the maximum peel stress and a corresponding reduction in the area under the curve. When considering the timedependent behavior of the polymeric adhesive, workof separation is also a function of time. Lower workofseparation would imply lower critical strain energy release rate, and therefore, lower resistance to crack growth. 54 4.6 Cohesive Layer Degradation and FEA Simulation The proposed cohesive layer model was implemented in an inhouse finite element code, NOVA3D. In the previous section, the parameters i C and i C~ are material constants to be determined experimentally. Moisture diffusion tests and fracture experiments are necessary to characterize these coefficients. 4.6.1 Cohesive Layer Degradation due to Moisture Concentration When moisture diffuses into the cohesive layer, two effects occur concurrently. First, the cohesive layer begins to swell thereby causing the local stress state to change due to the constraining effect of the surrounding adherends. Secondly, it is likely that moisture will penetrate the bulk cohesive layer to reach the cohesive interface (or interphase), and then rapidly diffuse along the interface. At the interface, water molecules typically react with the chemical bonds across the interface especially in the presence of tensile stress, because stress provides additional driving force for the bond rupturing process. Such chemical reactions transform strong covalent bonds to weak Van der Waals bonds, thereby significantly weakening the interface strength and fracture toughness. The bondstrength degradation could be important even when the change in moisture concentration is relatively small (~10%) [16]. Due to a lack of available bond degradation data from ongoing experiments to allow characterization of material coefficients in Eq. (415) at present time, a simple 55 bond strength degradation scheme (as illustrated in Fig. 42) was employed for the present analysis while preserving the basic framework presented in Eq. (415). Fig. 42. Depiction of the influence of moisture in cohesive layer on workofseparation 4.6.2 FEA Simulation Results Fully coupled stress and diffusion analyses were invoked in this investigation, analogous to simulating a wedgetest under wet conditions. For the DCB specimen, moisture diffusion analysis is activated only in the cohesive layer, with moisture boundary conditions applied at the exposed surface at x = 0 as depicted in Fig. 41. A cohesive layer thickness h mm c = 0.02 and beam height h mm b = 1 were used, Material properties used in this analysis are listed in Tables 31, and the diffusivity of the cohesive layer is 5.22×10−8mm2 / s . With time, moisture gradually diffuses from the exposed debond tip ( x = 0 ) into the originally dry elements in the cohesive layer. Due to the lack of the material 56 coefficients defining nonFickian diffusivities given by Eq. (415), linear Fickian diffusion with constant boundary concentration (m 1 10 8 g /mm3 b = × − ) was modeled in the cohesive layer for this demonstration case. Fig. 43 depicts the moisture concentration profiles plotted along the bond length with the origin at the original location of the debond tip at four different time steps. The corresponding analytical solutions for onedimensional Fick’s law are also plotted in Fig. 43 for diffusion model verification. It is evident that the concentration profiles predicted by the finite element diffusion analysis are in excellent agreement with the analytical solution for the linear case. Incidentally, the horizontal dashed line in Fig. 43 corresponds to 10% of saturation concentration and its purpose will be discussed in the following paragraphs. 0.0E+00 2.0E09 4.0E09 6.0E09 8.0E09 1.0E08 1.2E08 0 0.5 1 1.5 2 8.0E4 sec. 1.0E6 sec. 3.7E6 sec. 1.0E9 sec. Fick's Law Fig. 43 Comparison of predicted moisture concentration profiles A phenomenological stepfunction degradation law is assumed in the analysis such that when the local moisture concentration is greater than or equal to 10% of the X coordinate (mm) Concentration (g/mm3) 57 saturation concentration corresponding to the dashed horizontal line in Fig. 43, the corresponding maximum stress ( max σ ) in the cohesive tractionseparation law is reduced by 10% . For numerical stability, a linear degradation is assumed when the concentration is between 0 and 10%. The length of the moistureinduced cohesive strength degradation zone at various time steps is indicated by the intersection of the concentration profile and the horizontal dashed line in Fig. 43. t = 8 ×104 sec debond damaged undamaged 0.16 mm beam cohesive layer (a) Under dry condition t = 3.7 ×106 sec 0.57 mm (b) Degradation due to moisture concentration Fig. 44 Debond growth in a cohesive layer due to moisture degradation In this demonstration case, a constant beam tip displacement is applied at all times for the DCB, simulating a wedge test. The debond growth predicted at the crack tip by 58 the finite element analysis is 0.16 mm at initial time ( t = 0 ) under dry conditions as shown in Fig. 44(a). Over time, as moisture diffusion and subsequent bond degradation takes place, debond propagation occurs and the failure length increases to 0.57 mm over a period of 3.7×106 seconds, as depicted in Fig. 44(b). Failure (or debond) length is determined by observing if the transverse mechanical strain, 22 ε , has exceeded the prescribed maximum transverse strain, max ε , along the bond length. Fig. 45 shows that the transverse mechanical strain monotonically increases with time due to the formation of cohesive damage, material failure, and resultant debond propagation. The location of the crack tip of the failure zone at various time steps is indicated by the intersection of the mechanical strain and the horizontal dashed line in Fig. 45 indicating failure strain, max ε (debond length increases from 0.16 mm to 0.57 mm due to strength degradation caused by moisture diffusion). 0 0.5 1 1.5 2 2.5 3 3.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 8.0E4 sec. 1.0E6 sec. 3.7E6 sec. Normalized strain ε y /ε max X coordinate ( mm ) Fig. 45 Transverse mechanical strains and corresponding failure zones at different moments 59 The reduction in corresponding transverse stress distribution at different time steps is plotted along bond length in Fig. 46. The progressive reduction in the peak stress magnitude is due to degradation of bond strength caused by moisture ingress. Evolution of debond length with time is plotted in Fig. 47, and the decrease in reaction force due to the degradation of cohesive layer stiffness with time is shown in Fig.48. It is evident that for the present case, debond growth is driven by a synergistic interaction of moisture diffusion and transverse stress near the debond tip as discussed in the following paragraph. 1.0E+07 0.0E+00 1.0E+07 2.0E+07 3.0E+07 0 0.5 1 1.5 2 2.5 3 8.0E4 sec. 1.0E6 sec. 3.7E6 sec. Stress ( Pa ) X coordinate ( mm ) Fig. 46. Transverse stresses along bond length at different moments There are two important milestones during the diffusionassisted debond growth process, 4 t1 = 8×10 second and 6 2 t = 3.7 ×10 second (see Fig. 47). For 1 0 ≤ t < t , the moisture concentration ahead of the crack tip is not high enough to result in significant 60 degradation to the cohesive layer. Therefore, no debond growth occurs during this period. At roughly t = t1 , the 10% concentration front impinges on the debond tip and triggers debond growth as shown in Fig. 46. For 1 2 t ≤ t ≤ t , the moisture concentration ahead of the debond tip becomes high enough (>10%) such that the onset of bond degradation occurs, resulting in steady debond growth as shown in Fig. 47. In this regime, the debond growth is driven by the rate of propagation of the 10% concentration front, and therefore can be said to be diffusioncontrolled. Finally, for 2 t > t , the debond driving force (shown in Fig. 48) falls below a threshold value such that no further synergistic debond growth is possible. In this case, the peak transverse stress in the cohesive layer has decreased to max 0.9σ as shown in Fig. 46 due to the assumed moisture degradation characteristics. At the same time the transverse deformation (Fig. 44) and debond length (Fig. 47) have reached steady state values, respectively. Sqrt( hour ) Debond length and 10% concentration front ( mm ) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0 10 20 30 40 50 60 Debond length 10% concentration front t1 t2 Fig. 47 Predicted evolution of debond length with time 61 Fig. 48 Reaction force decreases with time 4.7 Conclusions A twodimensional cohesive layer constitutive model with a prescribed tractionseparation law was constructed from basic principles of continuum mechanics and thermodynamics, taking into account concentrationdependent and straindependent non Fickian hygrothermal effects that are likely to occur within a cohesive layer. Implementation of the model in a testbed finite element code was carried out and code verification was performed. Numerical simulation of a wedgetest involving debond growth caused by synergistic interactions between local stress and diffusing moisture was also presented to demonstrate the ability of the cohesive layer model to simulate environmental cracking. Reaction force ( N ) Sqrt ( hour ) 62 CHAPTER V NUMERICAL SIMUILATION OF PENINSULA BLISTER TEST The objective of this paper is to examine threedimensional and viscoelastic effects in the peninsula blister test for thin film adhesion, which is considered an effective way to measure the interfacial fracture toughness of bonded thin film structures. As will be demonstrated in this Chapter, analytical solutions are sometimes inadequate for accurately simulating peninsula blister experiments because linear elastic behavior is assumed in the thin film as well as in the debonding process zone at the interface. For this purpose, a threedimensional cohesive layer model and corresponding liquid loading simulation algorithm were developed and implemented into an inhouse testbed finite element analysis (FEA) code (NOVA3D). Features such as threedimensional mixedmode debonding, large displacements and rotations, residual stresses in the thin film, and timedependent (viscoelastic) effects in the thin film were considered. Numerical convergence and a stable debond growth were obtained over a fairly large debond length when suitably refined FEA mesh and liquid volume increment were employed. Detailed benchmark comparisons of the finite element predictions with analytical solutions and experimental results were performed, and good agreement was obtained. 63 5.1 Introduction Laminated thin film structures appear in a wide variety of applications such as multilayer structures, microelectronic devices, and packages. The increasing application of thin films in industry has made it necessary to better understand the nature and the reliability of the bonding between layers or bimaterial interface. Since poor bonding results in a crack or delamination, fracture mechanics is a natural approach for characterizing the resistance to failure and making durability or reliability prediction. For this purpose, several experimental methods, including peel tests, blister tests, indentation tests, scratch tests etc., have been used to determine the interfacial bond strength or toughness. The peel test is a simple and perhaps the most commonly used method for examining the strength of adhesively bonded layers. In a peel test, a thin flexible strip that is bonded to a substrate is pulled off at a certain angle to the underlying substrate. In the absence of any plastic deformation and residual stresses, adhesive fracture energy Γ can be derived directly [30] from the peel force through Γ = (1− cosϕ )P (51) where P is the peel force per unit width of the film and ϕ is the peel angle. In reality, plastic deformation in the peeling arm and crack tip singularity are present in all the peel test specimens, and were investigated by Kim et al. [31], Kinloch et al. [32], and Wei and Hutchinson [33]. It implies that a significant portion of the mechanical energy applied in peeling is dissipated and, if not properly accounted for, will lead to significant overestimates of the adhesive fracture energy. This is particularly true 64 when the yield strength of the peeling film is relatively low and the degree of adhesion is high. The pressurized circular blister test was developed by Williams [34] as an alternative approach to minimize dissipation effects and overcome some of the other drawbacks of the peel test. In many blister tests, a pressurized fluid (usually water is used and considered as incompressible) is injected between the faces of the crack between the film and substrate via a hole through the substrate. A uniform fluid pressure is applied to the thin film and hence the problems associated with mechanical contact and gripping are avoided. Furthermore, the axisymmetry of the circular blister configuration greatly minimize edge effects and effects due to specimen nonuniformity. The standard blister configuration (circular blister) was first introduced by Dannenberg [35] in 1961 to measure the adhesion of thick organic coating to metal. As tougher interface were examined, problems were encountered with tensile failure in blister film prior to debond at the interface. Another disadvantage of circular blister is that the strain energy release rate increases as the fourth power of the debond radius, thus making accurate evaluation of the debond radius essential, and resulting in a very unstable fracture specimen. Several alternative configurations with higher strain energy release rates were developed. These configurations include island blister proposed by Allen and Senturia [36], constrained blister proposed by Napolitano et al. [37] and Chang et al. [38], and peninsula blister proposed by Dillard and Bao [39]. In a constrained blister specimen, even though the strain energy release rate is independent of debond radius, the effect of friction between the delaminating film and the constraint could be significant, raising concern about the difficulty of analyzing contact 65 problems with large deformation. In the island blister specimen, debond growth occurs radially inwards on the island. A moderate increase in compliance is produced by a relative small increase in debond area, thereby giving rise to high strain energy release rates. As the island radius decreases with debonding, the calculated strain energy release rates increase in an unbounded manner. Fig. 51 Peninsula blister specimen and possible debond sites The peninsula blister is an extension of the island blister concept. Its name is derived from the fact that debonding will occur along a narrow “peninsula” which extends into the blister region, as shown in Fig. 51. The peninsula blister configuration not only maintains the high energy release rate of the island blister but also provides a constant energy release rate as the adhesive debond progresses. Added features include the larger debond area and additional data points that can be obtained from a single specimen. A consistent set of analyses by Liechti and Shirani [40] showed that the only configuration among circular, island, and peninsula blister specimen with a homogeneous 2 1 3 L l 2a 2b 66 delaminating layer that did not suffer any yielding was a relatively thick peninsula blister configuration. A subsequent cohesive zone fracture analysis by Shirani and Liechti [41] found that only 6% of the work input was dissipated as global plastic dissipation in a polyimide thin film on aluminum, a relatively tough interface. Disadvantages of the peninsula blister test include the lack of symmetry and difficulties in the fabrication of the specimen for certain material systems. Despite its promising advantages, the peninsula blister configuration has not gained wide acceptance in the adhesion testing community. One problem is the relatively complex nature of the specimen fabrication although this can be alleviated by suitable masks. Another is the relatively complex nature of the analysis particularly when large deformations and rotations and residual stresses are encountered. Both issues have recently been addressed [42]. One issue that is yet to be resolved on the analytical front is an accounting of the viscoelastic nature of many adhesives, particularly under long term loadings [43]. On the other hand, problems were also encountered in numerical simulation of peninsula blister test, such as the loss of axisymmetric geometry, mixed mode I and mode II nonlinear fracture (and even mode III), large deflections and rotations in the film and adhesive layer, and large debond length. Another problem is the special loading methodology of the blister test. It employs a liquidvolumecontrol loading approach which is obviously different from any general applied force or displacement loading methods. The purpose of the work presented here was to consider optimum geometries for peninsula blister specimens while accounting for geometric nonlinearities, viscoelastic effects and mixedmode debonding. A threedimensional cohesive layer model and 67 corresponding liquid loading simulation algorithm were developed and implemented into an inhouse testbed FEA code (NOVA3D). Numerical convergence was achieved with suitable increment of liquid volume and with reasonable FEA mesh refinement. Steadystate debond growth was obtained over a fairly large portion of the specimen. The finite element results were benchmarked against analytical solutions and experimental data. 5.2 Constitutive Law and Failure Criterion of Cohesive Layer Peninsula blister specimen is a full threedimensional configuration. Threedimensional stressstrain tractionseparation law is adopted to simulate the cohesive layer as shown in Eqs. (23, 24). In peninsula blister specimen, shear force perpendicular to the peninsula direction (corresponding to mode III fracture) was observed in peninsula blister test and simulation. A more complex failure criterion that contains the contribution of mode III debond was used in this case. + + = 1 IIIc III IIc II Ic I G G G G G G (52) Detailed description of the 3D cohesive layer model and its failure criterion is given in Chapter II. 5.3 Specimen Geometry and Energy Release Rates Proper choices of the geometry and, in particular, the aspect ratios of the peninsula blister specimen and regions within it are essential in order to ensure that 68 debond initiation and propagation occurs in the 1direction along the narrow peninsula area as shown in Fig. 51. Debonding at sites 2 and 3 is also possible but is undesirable. From the standpoint of linear elastic fracture mechanics, the energy release rate can be used to predict crack initiation and propagation; it can be viewed as the driving force for crack propagation or debond growth. When the deflection of the film is on the order of the film thickness, the theory of linear elastic thin plates can be applied in this case. Dillard and Bao [39] determined the analytical solution of the strain energy release rates for debonding at the three sites under a given liquid pressure p [ 5 5 ] 2 1 16 ( ) 1440 a a b Db G = p − − (53) D G p a b 288 2 ( )4 2 = − (54) D G p a 288 2 4 3 = (55) where D is the bending stiffness of the film given by 12(1 2 ) 3 −υ D = Eh (56) E , υ , and h are its Young’s modulus, Poisson’s ratio, and thickness, respectively. Numerical results from FEA analysis (NOVA3D) are in very good agreement with Dillard’s analytical solutions as shown in Fig. 52 and 53. It can be seen from Fig. 53 that higher energy release rate G can be achieved at site 1 for a small peninsula width. Therefore, when the peninsula width is relatively small (e.g. b 0.2 a < ), debond will 69 occur preferentially at site 1 rather than site 3 (because G1 > G3 ). It should also be noted that 1 G is greater than 2 G for all b a ratios, implying debond of the peninsula would always propagate along the peninsula length rather than across the width. As evident from Fig 52, in order to obtain large energy release rates with the peninsula specimen, the most effective way is to select small values of b a . It is suggested that b 0.1 a = is likely to be the most practical choice, as overly small b a will increase difficulties in fabricating specimen. 0.5 0 0.5 1 1.5 2 2.5 3 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Normalized peninsula width b/a Strain energy release rate G*D/P/P Analytical G1  FEA G2  FEA G3  FEA Fig. 52 Effect of relative peninsula width on the energy release rate Analytical G1 – FEA G2 – FEA G3 – FEA 70 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Aspect ratio of the blister region l /2a Strain energy release rate G*D/P/P FEA b/a = 0.05 FEA b/a = 0.1 FEA b/a = 0.2 Analytical Fig. 53 Effect of blister region aspect ratio on the energy release rate (Analytical solution is derived from l 2a > 2 ) Aspect ratio of the blister region ( 2 l a ) also deserves examination. The derivation of strain energy release rates discussed earlier (Eqs. 53 ~ 55) was based on a plate with infinite length and clamped along the two remaining boundaries. According to Timoshenko and WoinowskyKrieger [44], for an aspect ratio of 2, the deviation in the midplane deflection of a fully clamped plate (all four sides clamped) from the infinite plate solution is only 2%. It can be concluded that when the aspect ratio 2 2 l a > in the wider inflated region, and L l 2 a b − > − at the end of the peninsula (i.e. the narrow inflated regions at the two sides of the peninsula), errors in the above derivation from infinite plate solution are negligible. Furthermore, a constant energy release rate 1 G is obtained when debond occurs within the limits of 4a < l < L − 2(a −b) . FEA simulations point to 71 similar conclusion as shown in Fig. 53. When aspect ratio 1.5 2 l a < , value of energy release rate 1 G increases as the debond grows. For larger blister lengths, energy release rate remains constant and converges to the solution of the infinitely long plate. Additional verifications using numerical techniques will be presented in section 5. 5.4 TimeDependent Behavior of Polymeric Thin Film Polymeric thin films typically display some timedependent behavior, such as, creep and/or relaxation due to molecular motion. A linear viscoelastic model was employed in the current study to simulate the timedependent behavior of the polymeric film. The threedimensional constitutive equations of linear viscoelastic material are, τ τ σ t K t τ θ d t kk ∫ ∂ = − ∂ 0 ( ) 3 ( ) (57) τ τ μ τ d e S t t t ij ij ∫ ∂ ∂ = − 0 ( ) 2 ( ) (58) where ij s and ij e are the stress and mechanical strain deviators, respectively, kk θ =ε is the volumetric strain. The bulk relaxation modulus K(t) and shear relaxation modulus μ (t) can be expressed using Prony Series as, Σ= − ∞ = + n i t i K t K K e i 1 ( ) τ (59) Σ= − ∞ = + n i t i t e i 1 μ ( ) μ μ τ (510) 72 The linear viscoelastic model for the thin film was implemented into the FEA simulation code NOVA3D. FEA simulations of debond in a peninsular blister specimen using elastic and viscoelastic thin film were conducted and compared with test data. 5.5 Simulations of Debonding in the Peninsula Blister Experiment 5.5.1 Simulation under Small Deformation A threedimensional cohesive layer model with 20node quadratic brick element was implemented into a testbed FEA code NOVA3D. Convergence studies showed that stable debonding was obtained within a fairly large debond region if the aspect ratio of the element in the debond tip area and the increment of liquid volume (loading increment) are sufficiently small. As the liquid is injected into the blister at a constant rate, the pressure also increases linearly with time (or liquid volume) before debond occurs. When the failure criterion given by Eq. (52) is satisfied, the cohesive layer elements at the debond tip fail and debonding initiates. Instead of the node release and element deletion schemes used in most finite element codes, a failed element remains active in the subsequent analysis while the stiffness of the element is reduced to approximately zero. As the cohesive layer elements fail, new surfaces are generated and the pressurized liquid occupies the newly debonded volume . In the present work, the original aspect ratio at the beginning of simulation is set as 0.75 2 l a = . During the simulation, as the debond initiates and propagates, the film deflection and aspect ratio increase while liquid pressure decreases, as shown in Fig. 54. 73 When the aspect ratio is greater than 1.5, liquid pressure stabilizes and remains constant, and a constant strain energy release rate is obtained. However, the maximum deflection in the thin film increases continuously until the aspect ratio is greater than 2. As alluded to earlier, this result agrees with Timoshenko’s solution for infinite plate [44]. Debond process near the end of the peninsula was also studied in this paper. When the aspect ratio in the inflated region L l 1.5 a b − < − , pressure will start to increase due to the boundary effects (These simulation results will be discussed in Section 6.3). From both analytical and FEA simulations, it can be concluded that the effective test domain for peninsula blister test is: 1.5 2 l a > and L l 1.5 a b − > − . When debond occurs within this domain, constant strain energy release rates will be obtained. Time (sec) Liquid pressure (KPa) Deflection (w/h), Aspect ratio 0 50 100 150 200 250 0 10 20 30 40 0 0.5 1 1.5 2 2.5 Liquid pressure Film deflection Aspect ratio Fig. 54 Peninsula blister responses over large debond lengths 74 5.5.2 Peninsula Blister Specimen and Test Results The data used here was obtained from a peninsula blister specimen with dimensions L = 35.56mm , l =15.24mm , 2a =10.16mm , 2b =1.27mm ( 2 0.125 2 b a = ) with an initial aspect ratio 1.5 2 l a = (as illustrated in Fig. 51). The epoxy thin film in the current simulation was assumed as linear elastic and its Young’s modulus was E =1.78GPa , Poisson’s ratio was υ = 0.3 , and film thickness was 0.127 fh = mm . Residual stresses R R 5.7 x y σ =σ = MPa were obtained from bulge tests and, for analysis, were applied to the thin film before loading. The aluminum substrate was considered to be rigid both. The critical energy release rates (or fracture energy) extracted from the blister test was 100 ~ 110 / 2 cG = J m (including residual stresses) and 130 ~ 140 / 2 cG = J m (without residual stress). In addition, a critical pressure of 80.9kPa was observed at debond initiation, along with a constant pressure of 86.5kPa during stable debonding. The adhesive material which was used to bond the epoxy thin film to the substrate is Hysol EA9696, with an adhesive strength of max σ = 42.7MPa . The maximum strain in the adhesive layer is assumed to be max ε = 0.1. When applying cohesive layer model to this specimen, a cohesive layer thickness 0.0416 ch = mm was determined from Eq. (2 16) (including residual stresses in the thin film and a fracture energy 100 / 2 cG = J m ). It should be clarified here that in the present case, the thickness of adhesive layer is different from the cohesive layer thickness. The latter is determined by adhesive 75 material properties and corresponding fracture energy or debond energy as presented in Eq. (216). Consequently, only a portion of the adhesive layer is modeled by the special threedimensional cohesive layer elements and the remainder of the adhesive layer was modeled by regular elements. Exact measurements of the adhesive layer thickness in peninsula blister specimen are difficult and the debond process is not sensitive to the adhesive layer thickness. Therefore we can reasonably assume that the adhesive layer thickness is equal to the cohesive layer thickness. In other words, only one layer of cohesive layer elements was employed to represent the entire adhesive layer. 5.5.3 Simulation with Large Deformation and Residual Stresses For peninsula blister specimen, both the thin film and adhesive layer experience large deflections and rotations near the debond tip during debond growth. An updated Lagrangian (UL) formulation [45] was employed to cope with the geometrical nonlinear characteristics in the thin film (nonlinear straindisplacement relations, while the stressstrain relation in the thin film is still linear). The UL formulation was combined with Cauchy stress and Almansi strain tensor as energy conjugates. A typical finite element mesh for the threedimensional peninsula specimen is shown in Fig. 55. Only one layer of threedimensional 20node brick element was used to model both the thin film and a portion of the adhesive layer. In the plane of the film, element dimension along the debond (peninsula) direction is greatly dependent on the film thickness as smaller film thickness will experience larger bending deformation near the debond tip; a fine mesh was used to account for the sharp stress gradients and large curvature near the debond tip. Convergence studies showed that the critical pressure is 76 rather sensitive to the element dimension in the peninsula direction (Table 51). However, when the ratio of element size in the xdirection to film thickness was Δx hf = 0.4 , convergence was achieved. Y Z X Node # = 7389 Ele # = 978 Cohesive Ele # = 186 Fig. 55 3D FEA mesh with 20node brick element (Viewed from the bottom, elements on the top are the cohesive layer element) Table 51 Different element size in peninsula direction Element size ( f Δx h *) 1.6 0.8 0.4 0.2 Critical pressure (kPa) 110.5 100.0 95.5 94.0 Deflection ( f w h ) 3.87 3.70 3.62 3.59 Liquid volume (mm3) 46.6 44.2 43.0 42.6 *: Δx is the element dimension in peninsula direction, f h is the film thickness, w is the maximum deflection in the thin film 77 Residual stresses, due to mismatch of thermal coefficients of expansion are usually present in the thin film following cure. Assuming a thermal expansion coefficient α = 0.0001/o C and a temperature change ΔT = −23oC , resulted in uniform residual stresses R R 5.7 x y σ =σ = MPa throughout the thin film. The effect of residual stresses on critical pressure and maximum deflection were considered (Table 52) for equibiaxial residual stresses ranging from 0 to 15 MPa. The impact on critical pressure and maximum deflection were noticeable. The mode I/II phase angle φ ranged from 38o ~ 53o while the mode I/III phase angle ϕ varied from 9o ~ 13o over the residual stress levels that were considered. They remained nearly constant during debond growth. This indicates that mixedmode I/II debond failure is dominant and the mode III contribution is relatively small, and may be neglected in peninsula blister specimen. Table 52 Effect of residual stresses on debond process as predicted by FEA simulations Residual stress R y R x σ =σ ( MPa ) 0 3.0 5.7 10.0 15.0 Critical pressure ( KPa ) 79.7 88.6 95.5 106.0 116.6 Maximum deflection ( f w h ) 4.01 3.81 3.62 3.35 3.08 Liquid volume (mm3 ) 48.0 45.2 43.0 40.2 36.8 Phase angle (I and II) ( o ) 53.5 49.5 46.1 41.8 37.8 Phase angle (I and III) ( o ) 12.6 11.9 11.3 10.4 9.7 Note: Phase angle is the average value for the elements along the perpendicular direction of the peninsula. 78 Measured and predicted pressure and maximum deflection history are given in Fig. 56. Stable debond pressure predicted by the elastic FEA simulation is P = 96.0KPa , which is approximately 11% higher than the test value of P = 86.5KPa . The error is probably due to the fact that the viscoelastic effects in the thin film were ignored. Another possible source of error is that the simple failure criterion given by Eq. (52) may not exactly reflect the real debond conditions. The initial aspect ratio ( 2 l a ) of current peninsula blister specimen is 1.5, and therefore the pressure would remain constant during debond growth process while maximum deflection increases continuously as the debond propagates until the aspect ratio reaches 2.5. This is in agreement with that derived from the theory of thin plates. Threedimensional film displacement profiles from FEA simulations at different stages of debonding are shown in Fig. 57. The debonding process near the end of the peninsula region was also studied. According to the analytical solution derived from the theory of elastic plates, pressure will start to increase when the aspect ratio in the narrow region L l 1.5 a b − < − . FEA simulation results showed (Fig. 58) that when the debond tip approached the end of the peninsula (an aspect ratio L l 1.0 a b − < − ), the pressure and film deflection started to increase with debond growth. 79 Time (x1000 sec) Pressure (KPa) Debond length (mm) 0 1 2 3 4 0 20 40 60 80 100 0 1 2 3 4 5 6 7 8 9 Liquid pressure Debond length Debond initiation (a) Test data Time (s) Liquid pressure (KPa) Deflection (w/h), Debond length (mm) 0 50 100 150 200 0 20 40 60 80 100 0 1 2 3 4 5 6 7 8 Liquid pressure Film deflection Debond length (b) Simulation results Fig. 56 FEA simulation of peninsula blister responses (with large deformation and residual stresses) 80 (a) Aspect ratio = 1 (b) Aspect ratio = 2 (c) Aspect ratio = 2.5 Fig.57 Simulated 3D film displacement profiles at different stages of debonding 81 Aspect ratio (Ll)/(ab) Liquid pressure (KPa) Filmdeflection (w/h) 1.5 1 0.5 0 40 60 80 100 120 140 160 180 200 3.5 4 4.5 5 5.5 6 Liquid pressure Film deflection Fig. 58 Predicted liquid pressure and film deflection increase when debond approaches the end of the peninsula 5.5.4 Simulation Including TimeDependent Effect In the simulation of geometrically nonlinear, elastic debonding in thin film peninsula blister specimens with an initial aspect ratio of 1.5, critical pressure remains constant once the debond has initiated. However, in the actual peninsula blister test, a moderate increase in pressure was observed in the early stages of debond growth (Fig. 5 6). This suggested that there might be some timedependent effects in the epoxy thin film. Due to unavailability of the timedependent material properties of the epoxy thin film in the current blister test, a demonstration simulation considering linear viscoelastic effects in the thin film was performed to examine the influence of timedependent material properties on the critical pressure and the debond process. 82 For the purpose of highlighting the effect of linear viscoelasticity of the thin film, the material properties for a linear viscoelastic epoxy thin film corresponding to the elastic film were defined and scaled from the tested values as follows: the tensile modulus E0 =1.78GPa and Poisson’s ratio 0 υ = 0.3 of the elastic thin film correspond to the values of the viscoelastic tensile relaxation modulus, E(t) , and Poisson’s ratio, υ (t) , at the time t = 0 , i.e., 0 0 ( ) 1.78 t E Et GPa = = = 0 0 ( ) 0.3 t υ υ t = = = For a typical linear viscoelastic epoxy thin film, material properties such as bulk relaxation modulus K(t) , shear relaxation modulus μ (t) , tensile modulus E(t) , and Poisson’s ratio υ (t) at room temperature can be expressed in the form of a Prony Series. The Prony Series coefficients for the epoxy thin film are listed in Table 53 using data obtained for an epoxy [46]. For a liquid injection rate R = 80mm3 / hour , the simulation accounting for large deformations and residual stresses showed that the liquid pressure remained constant after debond initiation for a linearly elastic thin film. However, for a viscoelastic thin film (Fig. 59a), the liquid pressure did increase at the early stages of debond growth before reaching a steady state value. Stable debond processes with constant liquid pressures were obtained in both cases with constant but slightly different energy release rates. 83 Table 53 Material properties of a viscoelastic epoxy ( ) i Log τ i K (MPa) i μ (MPa) i E (MPa) i υ 1 4.3 24.64 0.22 0.21 3.11 x 104 2 3.3 25.54 0.45 0.11 4.95 x 104 3 2.3 104.24 20.88 68.48 8.31 x 104 4 1.3 29.14 20.32 29.72 1.13 x 103 5 0.3 150.32 22.11 93.33 3.19 x 103 6 0.7 155.61 160.32 387.51 2.40 x 103 7 1.7 440.03 161.77 491.45 5.44 x 103 8 2.7 330.48 310.46 758.45 6.53 x 103 9 3.7 438.32 159.01 552.47 1.19 x 102 10 4.7 226.82 187.85 377.34 1.66 x 102 11 5.7 206.62 71.09 256.68 1.74 x 102 12 6.7 245.83 63.80 137.01 2.59 x 102 13 7.7 49.88 21.35 95.43 6.03 x 102 14 8.7 169.51 20.31 47.03 1.18 x 102 15 9.7 34.86 11.30 53.36 3.00 x 103 16 10.7 82.59 6.45 15.69 3.04 x 103 17 11.7 11.97 3.08 17.01 2.20 x 103 18 12.7 40.06 0.25 4.14 9.92 x 104 19 13.7 8.54 1.29 7.19 1.57 x 104 ∞ 2337.41 12.34 36.91 0.4974 *: At room temperature 23.5°C and relative humidity 35% 84 Time (x1000 sec) Liquid pressure (KPa) 0 1 2 3 4 0 20 40 60 80 100 Test data Elastic film Viscoelastic film (a) Liquid pressure history Time (x1000 sec) Debond length (mm) 0 1 2 3 4 0 1 2 3 4 5 6 7 Test data Elastic film Viscoelastic film (b) Debond length history Fig. 59 Peninsula blister test data and FEA simulation results (with large deformation and residual stresses) 85 Measurements and simulation results for both elastic and viscoelastic thin films are shown in Fig. 59 and tabulated in Table 54 for comparison. Fig. 59(a) compares the pressure history as a function of time, while Fig. 59(b) depicts comparison of the debond growth history. From Fig. 59 and Table 54 it is evident that including viscoelastic timedependence in the thin film response results in a significant reduction in the error between simulations and test data. It should be noted that the slight discrepancy in the time scales of the simulation and test data is due to some leakage of liquid that was observed during the peninsula blister test especially at elevated pressure levels. When a lower liquid injection rate R =10mm3 / hour was employed for the viscoelastic thin film simulation, the debond initiation pressure decreased from 88.5KPa (corresponding to R = 80mm3 / hour ) to 84.0KPa while the peak pressure and stable debond pressure remained unchanged. Thus it is likely that the debond initiation pressure will be closer to the test value of 80.9kPa as the injection rate is decreased. Table 54 Peninsula blister test data and FEA predictions with elastic and viscoelastic thin film Test data Linear elastic thin film Viscoelastic thin film Event Time Pressure Time Pressure Error Time Pressure Error Debond initiation 2250 80.9 1880 94.5 16.8 2210 88.5 9.39 Peak pressure 2980 88.9 1910 96.0 7.99 2340 95.5 7.42 Stable debond 3400 86.5 1910 96.0 11.0 2490 91.5 5.78 Unit: Time  second, Pressure  kPa , Error  % (predicted pressure compared with test data) 86 5.6 Conclusions Peninsula blister specimen is an effective way to measure the interfacial fracture toughness of a variety of adhesive bond systems as it offers very high energy release rate values and maintains constant energy release rates over a large range of debond lengths. A threedimensional cohesive layer model and a corresponding threedimensional mixedmode failure (debond) criterion were developed based on the principles of continuum mechanics and fracture mechanics. It was implemented in a finite element code to simulate quasi static debonding in the peninsula blister test. Numerical convergence and stable debond growth were obtained over a fairly large range of debond lengths. The results from FEA simulations were in reasonable agreement with both an analytical solution and test data. Suitable geometries for the peninsula blister specimen were studied by both analytical and FEA approaches and guidelines were reiterated for the design of peninsula blister specimen. FEA simulation results also showed that large deflections, timedependent material behavior, and residual stresses in the thin film are important factors that should be considered in simulations of the peninsula blister test in order to extract the interfacial fracture toughness of a given adhesive system. 87 CHAPTER VI SIMULATION OF TIMEDEPENDENT DEBOND GROWTH The objective of this chapter is to model the synergistic bond degradation mechanisms that may occur at the interface between a fiber reinforced polymer (FRP) and a substrate. FEA Simulation of a wedge test was conducted, and the timedependent effect in the adhesive layer was involved in the simulation model. The results predicted by the computational model were benchmarked through comparison with analytical solutions and mixed mode fracture tests. Steady debond growth was obtained after the wedge front entered the originally bonded area, which is consistent with the observations from wedge tests. 6.1 Introduction Polymeric thin films usually display some timedependent behaviors due to molecular motion. This effect was also observed in blister test and creep test on a fully clamped polymeric plate. A linear viscoelastic model was introduced to simulate the timedependent behavior of the polymeric film by peninsula blister test (Chapter V). NonFickian hygrothermal effects on cohesive layer are derived in Chapter IV. It should be noted that in the present approach, expansion of Helmholtz free energy in terms of convolution integrals was not carried out to directly include viscoelasticity in the 88 cohesive stressstrain law. This is because the use of convolution integrals in addition to temperature and moisture dependence would render the cohesive stressstrain law intractable as far as characterization of the convolution coefficients is concerned. Therefore, in the interest of tractability, a simplified approach is employed where the ratedependent behavior in the cohesive layer is implemented through the characterization of ratedependence of the maximum stresses and maximum strains in the cohesive layer as presented in Table 61. The remainder of the polymeric adhesive outside the cohesive layer is modeled as a nonlinear viscoelastic continuum with timedependent constitutive behavior. The influence of temperature and moisture concentration on the workofseparation and on crack growth is derived from firstprinciples. The model is implemented in a testbed finite element code NOVA3D. Results predicted by the cohesive layer model are benchmarked through comparison with experimental data from mixedmode fracture experiments performed using a moving wedge test. Ratedependent debond process was also investigated with this model under different debond speeds. 6.2 Failure Criterion Based on WorkofSeparation Several mixed mode failure criteria of the cohesive layer are described in Chapter II. A new failure criterion based on fracture energy was introduced to predict the debond process of a wedge test. Fracture energy 2Γ of the adhesive was extracted from the wedge test. Let the fracture energy 2Γ equal to the workofseparation of cohesive layer with mixed mode I and mode II fracture, gives 2 ( II ) sep I sep sep Γ =ηφ =η αφ +βφ (61) 89 Where I sep φ and II sep φ are the workofseparation of pure mode I and mode II debond, respectively (Chapter II and Chapter IV) For pure mode I fracture under dry conditions (Δm = 0 ): c I sep h 16 max max φ = 9 σ ε For pure mode II fracture under dry conditions (Δm = 0 ): c II sep h 16 max max φ = 9 τ γ sep φ is the total workofseparation of the mixed mode debond. It is easy to derive the relation II sep I sep sep φ =αφ +βφ , where α and β are constants determined by the failure type (i.e. mode mix) of the cohesive layer. Comparing the idealized cohesive stressstrain curve in Fig. 23 with the uniaxial tension test data of the adhesive material in Fig. 61, it was found that the actual workofseparation from the test data is less than the theoretic one predicted by tractionseparation law due to the premature failure of epoxy primer. Thus a correction factor η = 0.75 was introduced to reflect the difference. Fig. 61 Stressstrain relation of epoxy under different strain rates (test data) 90 Therefore cohesive layer thickness c h can be obtained from Eq. (61) if the fracture energy 2Γ or workofseparation sep φ , which is the energy necessary to generate unit debond length, is known from experiment or analytical methods. 6.3 Nonlinear Viscoelastic Model and Fracture Energy A linear viscoelastic model is described in Chapter V, material constitutive equations are listed in Eqs. (57, 58) and the timedependent material properties can be expressed with Prony series as shown in Eqs. (59, 510). A nonlinear viscoelastic model for the bulk adhesive using modified free volume approach [47] was also implemented in the code. A strainbased formulation proposed by Popelar and Liechti [48] entails a nonlinear shiftfactor based on freevolume given by, ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + − ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ + Δ + Δ + Δ + Δ + = − s eff eff s s d v v v v d d f f B f T c T c f B Log a ε ε α β δθ α β δθ 2.303 2.303 ( ) (62) The generalized Jintegral for large deformation with pseudostress model proposed by Schapery [49] is also used for fracture modeling, ∫Γ ⎥⎦ ⎤ ⎢⎣ ⎡ ∂ ∂ = Φ − dL x u J dy T R R i v i (63) For this case, the critical work input required to initiate the crack (or, the fracture initiation energy) can be expressed in terms of the timedependent farfield parameter v J . ∫ ∂ ∂ Γ = − t v i R d J E D t 0 2 ( ) τ τ τ (64) 91 6.4 Moving Wedge Test To study the properties of polymer adhesive and its debond process under different environmental conditions and strainrates,
Click tabs to swap between content that is broken into logical sections.
Rating  
Title  Coupled Hygrothermal Cohesive Layer Model for Simulating Debond Growth in Bimaterial Interface 
Date  20060701 
Author  Wang, Yong 
Department  Mechanical Engineering 
Document Type  
Full Text Type  Open Access 
Abstract  Structural adhesive is extensively used in laminated composite, thin film structure and fiber reinforced polymer (FRP) bonded concrete structure. Debond initiation and propagation of adhesive structures under different environmental and loading conditions is the main concern of this study. Twodimensional and threedimensional cohesive layer models were constructed and implemented in a testbed finite element code. Model verification was successfully performed; at the same time the cohesive layer model was also applied to simulate the debond initiation and growth in various adhesive structures. Through the application of the cohesive layer model to different adhesive applications, such as DCB beam, wedge test specimen, peninsula blister specimen, FRP bonded structures and polymer coating structures, it is evident that the cohesive layer model is a very good tool to simulate the adhesive debond and debond growth in adhesive structures under different load conditions. Different stressstrain tractionseparation laws can be selected to cope with any particular nonlinear behavior in front of the debond tip and moisture diffusion in the cohesive layer can be directly simulated in a 2D or 3D case. The cohesive layer thickness is an important parameter in the cohesive layer model, which can be determined from the maximum deformation at the debond tip or the interfacial bond strength. Factors such as timedependent and strainratedependent material properties, synergistic bond degradation, strain hardening, long term and short term dynamic responses as well as hygrothermal effect on the cohesive layer debond initiation and growth were investigated by means of cohesive layer model. Material and geometric nonlinearity consideration was found to be essential when evaluating the structural dynamic responses under blast load. (Copies available exclusively from MIT Libraries, Rm. 140551, Cambridge, MA 021394307. Ph. 6172535668; Fax 6172531690.) 
Note  Dissertation 
Rights  © Oklahoma Agricultural and Mechanical Board of Regents 
Transcript  A COUPLED HYGROTHERMAL COHESIVE LAYER MODEL FOR SIMULATING DEBOND GROWTH IN BIMATERIAL INTERFACE By YONG WANG Bachelor of Science Tsinghua University Beijing, China 1986 Master of Science Tsinghua University Beijing, China 1989 Submitted to the Faculty of the Graduate College of the Oklahoma State University in partial fulfillment of the requirements for the Degree of DOCTOR OF PHILOSOPHY July, 2006 ii COPYRIGHT By Yong Wang July, 2006 iii A COUPLED HYGROTHERMAL COHESIVE LAYER MODEL FOR SIMULATING DEBOND GROWTH IN BIMATERIAL INTERFACE Thesis Approved: Thesis Adviser Dean of the Graduate College iv ACKNOWLEDGMENTS I would like to express my sincere thanks to my advisor, Professor Samit Roy, for his supervision, constructive guidance, financial support and inspiration throughout the study. His profound contribution to my development, professional and otherwise, is deeply appreciated. My appreciation extends as well to the members of my supervisory committee Dr. Hongbing Lu and Dr. Andrew S. Arena, Jr. from Department of Mechanical and Aerospace Engineering, Dr. Roger Zierau from Department of Mathematics, and to the faculty of the Department of Mechanical and Aerospace Engineering. I also would like to express my thanks to Professor Kenneth M. Liechti and his group member Dr. Soojae Park and Mr. Dewei Xu from The University of Texas at Austin for providing test data and valuable suggestions for this study. I wish to express my sincere gratitude to my parents for their unconditional love and confidence in me all the time. This dissertation is dedicated to my beloved wife Bo Lan and my daughter Hongfan Wang, without their love and support this work would have not been possible. I would like to thank them for their consistent encouragement and love, and for all that they have given to me in my life. Finally, the financial support of the National Science Foundation (Grant No. CMS0296167) is acknowledged with thanks. v TABLE OF CONTENTS Chapter Page I. INTRODUCTION.……………………………………………………………….……....1 1.1 Applications of Adhesive Material………………………………………..……….1 1.2 Cohesive Zone Model…………………………………………………..….....……3 1.3 Coupled Hygrothermal Effect on Cohesive Layer……………………..….…........5 1.4 Cohesive Layer Model……………………………………………..………………6 1.5 Objective and Contents…………………………………………..………….……..7 II. COHESIVE LAYER MODEL……………………………………………………….11 2.1 Cohesive Layer Configuration………………………..……..………………..…..11 2.2 Cohesive Layer Constitutive Equations……………………………….………….13 2.2.1 Triangular StressStrain TractionSeparation Law………………………......13 2.2.2 2D Cubic StressStrain TractionSeparation Law……………...……….......14 2.2.3 3D Cubic StressStrain TractionSeparation Law………………………..…16 2.3 Cohesive Layer WorkofSeparation…………………..…………………....……17 2.4 Failure Criteria of Mixed Mode Debond………………………………..……..…17 2.4.1 Criterion Based on Prescribed Maximum Strain………………………….....18 2.4.2 Criterion Based on Strain Energy Release Rate…………………..………....18 III. DEVELOPMENT OF AN ANALYTICAL SOLUTION FOR COHESIVE LAYER MODEL AND MODEL VERIFICATION…………………...…………...22 vi 3.1 Analytical Solution from Cohesive Zone Model……………………….……...…22 3.2 Analytical Solution from Cohesive Layer Model…………………………..….…27 3.3 Comparison between Analytical Solution and FEA Results………………..........32 3.4 Conclusions……………………………………………..........................…….…..40 IV. HYGROTHERMAL EFFECT ON COHESIVE LAYER………………………….42 4.1 Introduction……………………………………………………..…..….…………42 4.2 Fickian Diffusion…………………………………………………..………..……43 4.3 NonFickian Diffusion (Strain Assisted Diffusion)………………….………...…44 4.4 Cohesive Layer Diffusion Boundary Conditions……………………………....…51 4.5 Cohesive Layer WorkofSeparation………………………………………..……52 4.6 Cohesive Layer Degradation and FEA Simulation………………………..…..….54 4.6.1 Cohesive Layer Degradation Due to Moisture Concentration………….....…54 4.6.2 FEA Simulation Results…………………………………………….……......55 4.7 Conclusions………………………………………………………………..…...…61 V. NUMERICAL SIMULATION OF PENINSULA BLISTER TEST………………...62 5.1 Introduction……………………………………………………………….…..…..63 5.2 Constitutive Law and Failure Criterion of Cohesive Layer…………………...….67 5.3 Specimen Geometry and Energy Release Rates…………………..…...............…67 5.4 TimeDependent Behavior of Polymeric Thin Film…………………………...…71 5.5 Simulation of Debonding in the Peninsula Blister Test………………………..…72 5.5.1 Simulation under Small Deformation…………………………...……….......72 5.5.2 Peninsula Blister Specimen and Test Results…………………...…………...74 5.5.3 Simulation with Large Deformation and Residual Stresses………..…...…...75 vii 5.5.4 Simulation Including TimeDependent Effect………………………….……81 5.6 Conclusions……………………………………………………….………...…….86 VI. SIMULATION OF TIMEDEPENDENT DEBOND GROWTH…………….……87 6.1 Introduction………………………………………………………………....…….87 6.2 Failure Criterion Based on WorkofSeparation……………………….…........…88 6.3 Nonlinear Viscoelastic Model and Fracture Energy…………………..….…....…90 6.4 Moving Wedge Test……………………………………………………......…..…91 6.5 Numerical Simulation of Wedge Test……………………………………...……..92 6.6 Conclusions……………………………………………………………..…….......98 VII. FRP BONDED STRUCTURE UNDER BLAST LOAD…………….……….….100 7.1 Introduction……………………………………………………….…..…….…...100 7.2 Implicit Integration Methods…..………………………………..…..………..…101 7.2.1 Houbolt Method…………………………………………….………..……..102 7.2.2 Newmark Method………………………………………………………......103 7.2.3 Wilsonθ Method…………………………………………….……….……104 7.2.4 HHT Method………………………………………………………………..104 7.3 Blast Load……………………………………………………….………………106 7.4 Dynamic Response of FRP Bonded Structure………………………….……….109 7.4.1 Short Term and Long Term Responses in Cohesive Layer…..………….…109 7.4.2 Critical Debond Locations of FRP Bonded Beam...…………….……..…...112 7.4.3 Explosive at Different Locations………………..………………………….113 7.4.4 Concrete Beam with Initial Cracks………………..………………………..115 7.5 Conclusions……………………………………………………………………...118 viii VIII. DYNAMIC ANALYSIS WITH MATERIAL AND GEOMETRIC NONLINEARITY………………………………………………………………..119 8.1 Modeling of Material and Geometric Nonlinearity…………………………..…119 8.1.1 Material Nonlinearity………………………………………………….……120 8.1.2 Numerical Formulation for ElastoPlastic Problems……………………….123 8.1.3 Geometric Nonlinearity…………………………………………………….128 8.1.4 HilberHughesTaylor (HHT) method for temporal discretization………...131 8.2 Model Verification………………………………………………………………131 8.2.1 Single Element Extension Verification …………...………………………..132 8.2.2 Simple Shear Verification…………………………………………………..133 8.3 Dynamic Response of a Circular Steel Plate with Coating under Blast Load......137 8.3.1 Modeling of a Clamped Circular Steel Plate with Coating……………..…..137 8.3.2 Numerical Damping of HHT Method………………………………………138 8.3.3 Dynamic Responses under Different Load Levels………………………....141 8.3.4 Effects of Coating on Plate Responses………………………………..……143 8.3.5 TimeDependent Effect of Polymeric Coating……………………………..147 8.4 Conclusions………………………………………………………...……………149 IX. CONCLUSIONS……………………………………………….…………….……151 BIBLIOGRAPHY.………………………………………………………………...…....154 ix LIST OF TABLES Table Page 31 Material properties for concrete and epoxy adhesive………….…………………...33 32 Variation of damage length with different cohesive layer thicknesses………….…39 33 Analytical and numerical solutions of damage length……………………………...40 51 Different element size in peninsula direction………………………………...…….76 52 Effect of residual stresses on debond process as predicted by FEA simulations .....77 53 Material properties of a viscoelastic epoxy …………………………………….….83 54 Peninsula blister test data and FEA predictions with elastic and viscoelastic thin film……………………………………………………………..…85 61 max σ and max ε under different strain rate ε&……………………………………..…94 62 FEA results with various wedge speeds……………………………………………97 71 Impulses of different load distributions...…………………………………………108 72 Material properties of concrete, epoxy adhesive and FRP…………….……....….110 81 Constants defining the yield surface………………………………………………124 82 Material properties of steel plate and polymer coating…………………………...138 x LIST OF FIGURES Figure Page 2.1 Three characteristic zones near debond tip and corresponding stressstrain relations ...………………………………………………………………………..…………..12 22 Normalized triangular stressstrain tractionseparation law for a cohesive layer.....14 23 Cubic stressstrain tractionseparation law for a cohesive layer……………….…..15 24 Cubic stressstrain tractionseparation law and strain energy release rate for a cohesive layer……………………..……..……………………………………19 31 A cohesive layer in a double cantilever beam (DCB)……………………………...23 32 Schematic debonding of a double cantilever beam (DCB) showing cohesive layer with elastic zone and damaged zone ( x1∈[0,l])……….…….….....23 33 Transverse deformation comparison between beam centerline and cohesive interface using FEA…………………………………….…….…………..28 34 Transverse stress distribution in the beam and cohesive layer…………………….29 35 Finite element mesh of a DCB specimen with symmetry boundary conditions…...33 36 Stress distribution with different correction coefficient k …………………………35 37 Reaction force comparison for different normalized cohesive layer thicknesses…36 38 Reaction force comparison for different debond lengths…………………………..37 39 Damage length comparison for different cohesive layer thicknesses………………38 310 Stress distributions for different cohesive layer thicknesses…………………..….39 xi 41 A cohesive layer with moisture diffusion in a DCB beam…………………………45 42 Depiction of the influence of moisture in cohesive layer on workofseparation….55 43 Comparison of predicted moisture concentration profiles.………………………...56 44 Debond growth in a cohesive layer due to moisture degradation……..……………57 45 Transverse mechanical strains and corresponding failure zones at different moments ……………………………………………….………………58 46 Transverse stresses along bond length at different moments………………………59 47 Predicted evolution of debond length with time……………………………………60 48 Reaction force decreases with time………………………………..……………….61 51 Peninsula blister specimen and possible debond sites…………………………...…65 52 Effect of relative peninsula width on the energy release rate………..……………..69 53 Effect of blister region aspect ratio on the energy release rate (Analytical solution is derived from l 2a > 2)…………………………………….70 54 Peninsula blister responses over large debond lengths..……………………………73 55 3D FEA mesh with 20node brick element………………………………………..76 56 FEA simulation of peninsula blister responses (with large deformation and residual stresses)…………..………………….………79 57 Simulated 3D film displacement profiles at different stages of debonding….……80 58 Predicted liquid pressure and film deflection increase when debond approaches the end of the peninsula……………………...…………………….…..81 59 Peninsula blister test data and FEA simulation results (with large deformation and residual stresses)………….…………………….…….84 61 Stressstrain relation of epoxy under different strain rates (test data)………….…..89 xii 62 Specimen of a moving wedge test……………………………………………….…91 63 Debond length vs. time (test result and FEA prediction)…………………………..93 64 Vertical reaction force vs. time (test result and FEA prediction)………………….93 65 Fracture energy 2Γ vs. debond speed……………………………………………..95 66 FEA mesh and contour for Jintegral………………………………………………97 71 A simply supported concrete beam bonded with FRP under blast load…………..107 72 Triangular and exponential distributions of blast load……………………..……..107 73 Schematic FEA mesh and distribution of blast load along the beam……………..109 74 Short term response y σ in the top and bottom cohesive layer at midspan...…….111 75 Long term response x σ in the top and bottom cohesive layer at midspan...…….111 76 Typical distribution of stress y σ in a cohesive layer ………………………….…113 77 Typical distribution of stress x σ in a cohesive layer ………………………….…113 78 FEA mesh of a simply supported beam under off side load………………………115 79 Stress limits of x σ as function of explosive location along beam axis…………...115 710 A simply supported FRP bonded beam with an initial crack in concrete …....….116 711 Comparison of axial stress x σ in the cohesive layer before and after concrete cracking (Cracking occurred at t = 50ms )…………………………….117 712 Comparison of shear stress xy τ in the cohesive layer before and after concrete cracking (Cracking occurred at t = 50ms )………………………….….117 81 Mathematic models for representation of strain hardening behavior………….….121 82 Elastoplastic linear strain hardening behavior for uniaxial case…………..……..122 83 Incremental stress changes at a point in an elastoplastic continuum…………….126 xiii 84 Single element extension (plane strain conditions)……………………………….132 85 Stressstrain relation comparison between NOVA and ABAQUS……………….133 86 Single element under simple shear (displacement control)……………………….134 87 Single element shear stress vs. shear strain for simple shear ………….…………134 88 Applied shear force vs. shear strain (global comparison, 19×19 mesh)…………135 89 Shear stress at the center of the specimen (local comparison, 19×19 mesh).……135 810 Equivalent plastic strain at the center of the specimen (19×19 mesh) ……..…..136 811 Deformation of a simple shear specimen with 200% shear strain (19×19 mesh) …………………………………………...………………………………………136 812 Axisymmetric model of a circular steel plate with coating…………....………...137 813 FEA mesh for simulating steel plate with coating (dimensions not scaled)……..138 814 Numerical damping of HHT method on highfrequency modes……………...…139 815 Numerical damping of HHT method on plate deflection………………………..140 816 Central deflection history under various peak pressure levels.……………..…...142 817 Maximum equivalent plastic strain at plate center under various peak pressure levels……………………………………………………...…..….142 818 Minimum and maximum deflections under various peak pressure levels.……....143 819 Central deflection history of the plate…………………………………………...144 820 Radial stress history in polymer coating…………………………………………145 821 History of equivalent plastic strain in steel plate…………………………….…..145 822 Maximum equivalent plastic strain at the bottom surface of steel plate…..…….146 823 Equivalent plastic strain in steel plate and permanent deformation (not scaled) …………………………………………………………………………..………..146 xiv 824 Central deflection history of the plate…………………………………….……..148 825 Radial stress history in steel plate at the plate center……………………...…….148 826 Radial stress history in polymer coating at the plate center ……………….…....149 xv NOMENCLATURE ij σ Stress component ( Pa ) ij ε Strain component σ Normal stress ( Pa ) τ Shear stress ( Pa ) max σ Maximum normal stress in cohesive layer ( Pa ) max τ Maximum shear stress in cohesive layer ( Pa ) max ε Maximum normal strain in cohesive layer max γ Maximum shear strain in cohesive layer c h Cohesive layer thickness (m ) φ sep Cohesive layer workofseparation ( J /m2 ) G Strain energy release rate ( J /m2 ) E Young’s modulus ( Pa ) υ Poisson’s ratio ρ Material density (Kg /m3 ) ψ Helmholtz free energy ( J /m3 ) m Moisture concentration (Kg /m3 ) xvi Γ Fracture energy ( J /m2 ) R External nodal force vector ( N ) F Internal nodal force vector ( N ) M Mass matrix K Stiffness matrix Δt Time step (sec ) D Elasticity matrix ep D Elastoplastic stiffness matrix ys σ Yield stress ( Pa ) H′ Strain hardening function ( Pa ) p ε Equivalent plastic strain σ Equivalent stress ( Pa ) i J Stress invariant i J ′ Deviatoric stress invariant 1 CHAPTER I INTRODUCTION 1.1 Applications of Adhesive Material The use of structural adhesives is rapidly increasing as they offer distinct advantages over conventional mechanical fastening technique. Laminated composites and thin film structures are some of the most popular applications in industry. Fiber reinforced polymer (FRP) composite is a new application of cohesive material in industry and civil engineering. It is necessary to better understand the nature and the reliability of the bonding between the layers or at the bimaterial interface. Many polymeric materials, including structural adhesives, exhibit nonlinear and timedependent behaviors and quite sensitive to the change of temperature and moisture penetrant concentration. Therefore the loadcarrying requirements of the cohesive layer, timedependent material properties, and coupled hygrothermal effect in the adhesive layer, stress concentration and nonlinear deformation near the crack tip, crack initiation and propagation (debond) within the interface between two materials are some of the essential factors that should be considered to get a reliable solution. Fiber reinforced polymers are a class of advanced composite materials that have been extensively used as lightweight, performanceenhancing materials in aerospace, automobile, and defense industries for quite some time. Over the past few years there has been extensive research into their potential applications in the construction industry. 2 However, the actual application of FRP composite in civil engineering sector has been slow especially as concrete reinforcement. One of the chief reasons for their slow acceptance is because of a lack of reliable predictive models and sound design guidelines for their use in civil infrastructure applications. One promising prospect of the application of FRP composite materials in civil engineering is infrastructure repair and retrofit. FRP materials have been used to strengthen the concrete beam element of buildings and bridges with low cost and high strengthening effect [1]. This allowed increasing the strength and/or ductility of these structures while benefiting from the FRP material advantages including: ease of application, high strengthtoweight ratio, and excellent resistance against corrosion and chemical attacks. Ritchie et al. [2] studied the behavior of concrete beams strengthened by bonding FRP (glass, carbon, and aramid) plates to the tension zone and showed that FRP reinforcement increased beam stiffness by 1779% and beam ultimate strength by 4097%. New uses of FRP sheets to upgrade the resistance of steel structures have recently been studied. A considerable increase in the strength and stiffness of the rehabilitated steel bridge girders was observed [3]. A major concern for such retrofitting is the debonding of polymeric adhesive that could compromise the reinforcing effect of the FRP. When exposed to harsh environments, degradation of the adhesive bond could lead to delamination of the FRP reinforcement that could ultimately lead to catastrophic failure. Thin film structure appears mainly as coating in a wide variety of applications and in multilayer structures in microelectronic devices and package. Since poor bonding results in crack or delamination, fracture mechanics is a natural approach for 3 characterizing the resistance to failure and making durability or reliability prediction. At the same time, several experimental methods, including peel test, blister test, indentation test, and scratch test etc., have been used to determine the interfacial bond strength or toughness. One problem common to all of these methods is that global plastic dissipation makes it difficult to extract the true toughness of the interface. A new test method, peninsula blister test, can minimize the plastic dissipation and hence be an effective approach to measure the fracture toughness of the thin film structures. Laminated composites have the advantage of low weight and high strength compared to the structural metals, and hence obtain an increasing application in aerospace industry. Delamination, which is created when two layers debond from each other, is a common type of failure mode in layered composites without the throughthethickness reinforcement. The initiation of delamination growth is usually controlled by mode I and mode II fracture toughness. Numerical approach with mixed mode failure criterion is an effective way to investigate the delamination of the layered composites. 1.2 Cohesive Zone Model An important issue when considering failure is the observation that most engineering materials are not perfectly brittle in the Griffith sense but display some ductility after reaching the strength limit. In fact, for most engineering materials there exists a small zone in front of the crack tip, in which smallscale yielding, micro cracking and void initiation, growth and coalescence take place. If this fracture process zone is sufficiently small compared to the structural dimensions, linearelastic fracture mechanics concept can apply. However, if this is not the case, the cohesive forces that exist in this 4 fracture process zone must be taken into account. The most powerful and natural way is to use cohesive zone model, which was introduced by Barenblatt [4] and Dugdale [5] for elasticplastic fracture in ductile metals, and for quasibrittle materials by Hillerborg et al. [6] in his socalled fictitious crack model. The fracture process zone approach of Needleman [7, 8] and Tvergaard and Hutchinson [9, 10] involves attributing a prescribed tractionseparation law to the interface and, because it allows crack growth to occur, the associated plastic dissipation from loading and unloading of points that are passed by the crack front is rigorously accounted for. As a result, the selected tractionseparation law determines the workofseparation (or adhesive fracture energy), which is the work required to create a unit area of fully developed crack [11]. In the past two decades or so, cohesive zone models have become very popular and have been recognized to be an important tool for describing fracture in engineering materials. Especially when the crack path is known in advance, either from experimental evidence or because of the structure of the material (such as in laminated composites), cohesive zone model has been used with great success. Song and Waas [12], Shawan and Waas [13], and ElSayed and Sridharan [14] successfully employed cohesive zone model to investigate the fracture properties in laminated composites. In those cases, the finite element mesh was constructed such that the known crack path coincides with the element boundaries. The most common failure form of FRP composite plate bonded concrete structure is the delamination of the FRP plate from the concrete component. The debond process of FRP plate usually initiates and propagates along the adhesiveconcrete or adhesiveFRP 5 interface, which is known in advance. The cohesive zone model is thus a good tool for investigation of local fracture processes in FRP delamination. Fatigue crack growth is traditionally characterized via linear elastic fracture mechanics concepts where crack growth rates are correlated with the change in energy release rate or the maximum value of the energy release rate in a cycle. This approach has worked well for metals and polymers alike, especially in dry, room temperature environments, where conditions are still generally linearly elastic. Correlation between crack growth rates and elastic fracture parameters do become suspect in polymers near their glass transition and when saturated by a solvent. Cohesive zone modeling offers a solution to this difficulty in the sense that if the neartip damage can be accounted for in the tractionseparation law of the interphase then the local nonlinear inelastic behavior of the material can be coupled into any analysis directly [15]. 1.3 Coupled Hygrothermal Effect on Cohesive Layer Moisture can cause a host of reliability problems at interfaces including interface bond degradation and debonding. Two mechanisms can be identified. First, moisture at an interface can reduce the interface bonding strength dramatically by altering the chemical bonds. Second, when an interface with a crack or a cracklike defect is under tensile stress, stress corrosion may allow crack growth at stresses much lower than critical fracture would require [16]. The influence of moisture diffusion on crack growth along an interface is not yet fully understood. Environmental cracking in a polymer typically occurs in the presence of a penetrant, such as moisture, and mechanical strain. It has been postulated that the 6 mechanism involved in environmental crack growth in a polymer involves a small zone of craze formation and/or plasticization at the crack tip. For thermoset resins, such as epoxy, energy absorption at the crack tip is primarily by a shear yielding process and not by crazing. Consequently, for a thermoset epoxy, the zone of plasticization ahead of the crack tip must be determined using a diffusion law for nonporous media, such as Fick’s law. However, quite frequently, polymer composites exhibit deviations from the classical Fickian treatment, termed as anomalous or nonFickian diffusion, especially at elevated temperatures and stress levels, and at high relative humidity. Sophisticated hygrothermal models have been developed and verified by Roy [1720] to account for anomalous diffusion. 1.4 Cohesive Layer Model Cohesive layer model employs a thin layer of material, which is placed between two adjacent layers for laminated structures and multilayer structures, or along the bimaterial interface, or along a predicted cracking path in a single material (e.g. concrete and metal) to simulate the elasticplastic failure in ductile or quasibrittle material in the vicinity of the debond tip. It allows debond (or failure) to initiate and grow in these elements along a prescribed debond path. The vicinity of debond tip can be divided into three characteristic zones: elastic zone, damage zone, and debonded zone. The corresponding stressstrain tractionseparation relation may take different forms in each zone for different materials, different loading and environmental conditions to cope with any particular nonlinear behavior in front of the debond tip. 7 The thickness of cohesive layer is an important parameter in the cohesive layer model. It should be noted that the cohesive layer thickness is not arbitrary in cohesive layer model, and it is related to some characteristic lengthscale of the debond process, such as crack opening displacement (COD). Cohesive layer thickness can be determined from the maximum deformation at the debond tip and the maximum strain max ε in cohesive layer at failure. Environmental degradation is usually found in adhesive material and therefore results in the change of the material properties such as max ε and max σ (maximum stress) When investigating the moisture degradation, twodimensional and threedimensional moisture diffusion in the cohesive layer can be directly simulated in a twodimensional or threedimensional cohesive layer. Mixed mode I and mode II fracture is the common failure form of cohesive layer and sometimes even includes mode III. Mixed mode failure and corresponding failure criteria can also be easily implemented in cohesive layer model to predict the debond process in adhesive layer. 1.5 Objective and Contents The objective of this research is to construct a cohesive layer model from fundamental principles of continuum mechanics and thermodynamics, take into account the strain rate dependent material properties, nonFickian Hygrothermal effects as well as diffusioninduced degradation in the cohesive layer. By means of the cohesive layer model, the effect of ratedependent material properties, environmental degradation of the adhesive material, dynamic response involving material and geometric nonlinearity under blast load, quasistatic debond initiation and propagation of the adhesive layer were 8 studied to provide a better understanding of the strengthening effect and reliability of FRP plated structures. In Chapter I, a brief literature review of the application and research status of structural adhesive and FRP bonded structure, and the motivations of this research are presented. In Chapter II, twodimensional and threedimensional cohesive layer constitutive models with prescribed tractionseparation laws were constructed from fundamental principles of continuum mechanics and thermodynamics. Based on debond tip deformation, workofseparation or strain energy release rate, criteria for mixed mode I and mode II debond (and even includes mode III) were developed to predict the debond initiation and propagation of the cohesive layer. In Chapter III, an analytical solution was derived by introducing a correction term into the original Williams’ solution to predict the transverse stress in a cohesive layer when considering the deformation of a stiff substrate. Implementation of the cohesive layer model into a testbed finite element code was carried out and code verification was performed. Benchmark comparisons of finite element prediction of both global critical load and local stress field with analytical solution for a DCB specimen resulted in good agreement after modifications were made to the original Williams’ solution. A sensitivity study was conducted to evaluate the influence of cohesive layer thickness on local parameters such as damage zone length, and global parameter such as critical force. In Chapter IV, a twodimensional cohesive layer constitutive model involving strain dependent, nonFickian hygrothermal effects as well as diffusion induced degradation in the cohesive layer was constructed. Numerical simulation of a wedgetest 9 including debond growth caused by synergistic interactions between local stress and diffusing moisture was also presented to demonstrate the ability of the cohesive layer model to simulate environmental cracking. In Chapter V, a threedimensional cohesive layer model and corresponding mixed mode failure (debond) criterion were implemented in a testbed finite element code to simulate the full threedimensional peninsula blister test. Issues such as large deformation, timedependent material behavior, and residual stresses in the thin film were considered in the simulation model. Distinctive numerical techniques were successfully employed to simulate the unique liquid loading process. FEA simulation results were also compared with analytical solution and test data. Good agreement was obtained. In Chapter VI, cohesive layer model with strainrate dependent tractionseparation constitutive law was implemented in a testbed FEA code to simulate a moving wedge test. Timedependent material properties of the adhesive material were considered and quasistatic debond growth of the adhesive layer was successfully simulated by this code. Results predicted by the computational model were benchmarked through comparison with analytical solutions and mixed mode fracture tests. In Chapter VII, cohesive layer model was used to study the dynamic response of a FRP bonded concrete beam under blast loading. Implicit HilberHughesTaylor (HHT) method was employed in the model to allow better control of numerical damping. Long term and short term responses were obtained and their effects on the failure of the adhesive layer were investigated. Dynamic responses of the structure with an initial crack and its effect on debond initiation were also studied. 10 In Chapter VIII, a twodimensional implicit dynamic finite element formulation including material and geometric nonlinearity was derived and implemented into a testbed FEA code. Model verification under very large deformation was successfully performed through comparison with ABAQUS FEA predictions. Subsequently, the NOVA3D FEA model was applied to a circular steel plate with a polymer coating subjected to intensive blast loading, and the effect of polymer coating on the nonlinear dynamic response was numerically investigated. In Chapter IX, conclusions are presented based on the study of the cohesive layer model and its applications on various engineering structures. 11 CHAPTER II COHESIVE LAYER MODEL This chapter gives a detailed description of the cohesive layer model for twodimensional and threedimensional cases. It includes the definition of cohesive layer, the constitutive laws for the cohesive layer, the concept of workofseparation, and the mixed mode failure criteria. 2.1 Cohesive Layer Configuration Ductile polymeric adhesive materials usually have a nonlinear normal and shear stressstrain response. In the event of crack initiation and propagation in such polymeric materials, there exists a damage zone ahead of the debond tip, in which, craze and void initiation, growth and coalescence take place. The cohesive forces in the damage region must be taken into account to capture the behavior of the failing material in this zone, especially if the zone size is not sufficiently small compared to characteristic structural dimensions. The vicinity of debond tip can be divided into three common zones: an elastic/viscoelastic zone, a damage zone, and a debonded zone, as depicted in Fig. 21(a). The corresponding stressstrain (or tractionseparation) relation may take different nonlinear forms in each zone for different materials and different loading conditions. 12 (a) Debond tip traction force 0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Normalized strain Normalized stress (b) Tractionseparation law near debond tip Fig. 2.1 Three characteristic zones near debond tip and corresponding stressstrain relations In order to model this nonlinear behavior near the debond tip, a cohesive layer, which is a thin layer of cohesive finite elements, can be placed between two adjacent layers for laminated structures and multilayer structures, or along the bimaterial interface, or along a predicted cracking path in a single material (e.g. concrete and metal). It allows debond (or failure) to initiate and grow in these elements and different stressr z Elastic zone Damage zone Debonded zone 13 strain tractionseparation laws can be selected to cope with any particular nonlinear behavior in front of the debond tip. An example of the tractionseparation law for a cohesive layer is shown schematically in Fig. 21(b). The thickness, c h , of cohesive layer is an important parameter in the cohesive layer model. It is not arbitrary, but is directly related to a characteristic length parameter (δ ), such as crack opening displacement (COD). In this simulation, δ ε max c = h and max ε is the maximum strain that could be reached at debond tip. Environmental degradation in cohesive material is included in the cohesive layer model through the change of the material properties max σ and max ε . When investigating the moisture induced degradation, twodimensional and threedimensional moisture diffusion in the cohesive layer can be directly simulated. Mixed mode failure and corresponding failure criteria can also be easily established to predict the debond process. 2.2 Cohesive Layer Constitutive Equations 2.2.1 Triangular StressStrain TractionSeparation Law Triangular stressstrain tractionseparation law is a simple and commonly used model for cohesive material (Fig. 22), especially in theoretical analysis. Considering the stress and strain of mode I debond (opening mode) in the direction perpendicular to the debond surface, three types of zone are defined as follow: Elastic zone: when the transverse strain max 3 ε ≤ ε , stress linearly increases with strain, stress reaches its maximum value max σ =σ at max 3 ε = ε . 14 Damage zone: when the transverse strain max 3 ε > ε , stress decreases gradually from its maximum to zero as strain approaches max ε . Debond (failure) zone: when the transverse strain max ε >ε , stress remains zero which implies a full debond or separate of the cohesive layer. Fig. 22 Normalized triangular stressstrain tractionseparation law for a cohesive layer The maximum stress max σ and maximum strain max ε , which are material properties, are functions of time, strain rate, temperature and moisture concentration etc., and represent the prescribed maximum stress and strain that could be reached in cohesive layer when the cohesive layer debonds along a specified direction. 2.2.2 2D Cubic StressStrain TractionSeparation Law Based on fundamental principles of continuum mechanics, for twodimensional case, a more accurate cohesive layer constitutive relationship takes the cubic form as 15 employed by Needleman [7] (Fig. 23). In the direction perpendicular to the debond surface, the transverse normal stress is given by, ⎪ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎪ ⎨ ⎧ = < > − + ≤ ≤ = 0 4 27 4 27 0 1 ( 2 ) 0 1 4 27 22 max 22 22 max max 22 22 3 22 2 max 22 22 22 ε σ ε ε ε σ ε σ ε ε ε ε σ (21) 1.4 1.2 1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 1 1.2 0.4 0.2 0 0.2 0.4 0.6 0.8 1 1.2 Fig. 23 Cubic stressstrain tractionseparation law for a cohesive layer In the direction of debond growth, linear elastic response is assumed, and the axial normal stress is given by, 11 max 11 max max 11 ε σ ε ε σ σ = = (22) where, the normalized strain in a given direction is defined as max ε ε ε ii ii = (no sum on i). σ / σmax ε / εmax 16 A similar treatment to the one developed in Eq. (21) is prescribed for the shear 12 σ response, with the proviso that the shear stress is independent of the sign of shear strain. 2.2.3 3D Cubic StressStrain TractionSeparation Law Threedimensional cohesive layer treatment is often necessary for some structures like peninsula blister specimen, which is used to measure the interfacial fracture toughness in thin film structures. A full threedimensional cohesive layer model was constructed in this study to meet this requirement. By extending the above twodimensional cohesive layer model just described, the constitutive law for a threedimensional cohesive layer can be expressed in a similar manner as shown in Eqs. (23) and (24). Again, nonlinear responses are considered for transverse stress components 31 σ , 32 σ , and 33 σ , while for other stress components linear elastic response is assumed. 2 3 max 3 3 3 3 3 3 max 3 3 27 ( 2 ) 0 1 4 0 1 ( 1, 2, 3) 27 0 4 i i i i i i i i i σ ε ε ε ε σ ε σ ε ε ⎧ − + ≤ ≤ ⎪⎪⎪⎪ = = > ⎨⎪⎪⎪ < ⎪⎩ (23) max max max jk jk jk σ σ ε σ ε ε = = ( j, k =1, 2 ) (24) 17 2.3 Cohesive Layer WorkofSeparation For a given stressstrain tractionseparation law, the workofseparation (or strain energy) is the work needed to fully separate a unit area of cohesive layer, which is given by the total area under the prescribed stressstrain curve. Under the triangular stressstrain tractionseparation law as shown in Fig. 22, for pure mode I debond, c c I sep h d h max max 0 2 ( ) 1 max φ σ ε ε σ ε σ = ∫ = (25) and for pure mode II debond, c c II sep h d h max max 0 2 ( ) 1 max φ τ γ γ τ γ γ = ∫ = (26) Similarly, under the cubic stressstrain tractionseparation law as shown in Fig. 2 3, c c I sep h d h max max 0 16 ( ) 9 max φ σ ε ε σ ε σ = ∫ = (27) c c II sep h d h max max 0 16 ( ) 9 max φ τ γ γ τ γ γ = ∫ = (28) where c h is the thickness of the cohesive layer. 2.4 Failure Criteria of Mixed Mode Debond Mixed mode I and mode II debond is the common failure form of cohesive layer, while pure mode I or mode II debond is only a special case under certain conditions. It is 18 necessary to establish a failure criterion such that it contains both the contributions of mode I and mode II debond, and in some cases even includes mode III. 2.4.1 Criterion Based on Prescribed Maximum Strain At the debond tip, when the strains satisfy the following condition, cohesive layer will debond [21], δ γ γ α ε ε α = ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ + ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ 2 max 2 2 max 1 y xy (29) where y ε and xy γ are the transverse normal strain and shear strain respectively, max ε and max γ are the prescribed normal and shear failure strains of the cohesive layer, respectively. 1 α , 2 α , and δ are constants and 1 1 2 α =α =δ = was taken in this study. 2.4.2 Criterion Based on Strain Energy Release Rate The strain energy release rate in the cohesive layer during mixed mode debond, G , due to the tractionseparation force can be partitioned into the opening (mode I) and shear (mode II) components, I G and II G respectively, in such a way that, I II G = G +G (210) Each individual component can be calculated by integrating the mode I and II tractionseparation curves (Fig. 24) = ∫ t G d I ε σ ε ε 0 ( ) (211) 19 = ∫ t G d II γ τ γ γ 0 ( ) (212) where max σ σ = σ , max τ τ = τ , max ε ε = ε , max γ γ = γ are normalized stresses and strains in the specific directions, respectively. 1.0 0.8 0.6 0.4 0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0.2 0 0.2 0.4 0.6 0.8 1 1.2 Normalized strain Normalized stress GI Fig. 24 Cubic stressstrain tractionseparation law and strain energy release rate for a cohesive layer Considering the energy required to separate the cohesive layer, the cohesive layer debond process can be better predicted by means of the following criterion e G G G G n II II m I I C C = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ (213) where, I G and II G are the respective values of the ambient energy release rates given by the area under the corresponding stressstrain curves under a given applied loading. 20 IC G and IIC G are the critical strain energy release rates in pure mode I and mode II debond, respectively. m, n, and e are material constants, Kutlu and Chang [22] found that m = n = e = 1 provided the best fit to their experimental data. When neglecting energy dissipation in the bulk adhesive, critical strain energy release rate of cohesive layer is equal to the workofseparation sep φ (the total area under mode I or mode II tractionseparation curves for 0 ≤ε ≤ 1, Fig. 24), which is the energy necessary to generate unit debond length (2D case) or area (3D case). Thus for the cubic stressstrain tractionseparation law, integrating Eq. (23) over the limits [0, ε ] for the case of pure mode I debond and [0, γ ] for the case of pure mode II debond, ) 4 1 3 2 2 (1 4 ( ) 27 2 3 4 max max 0 σ ε ε σ ε ε ε ε ε I = c ∫ = c − + G h d h (214) ) 4 1 3 2 2 (1 4 ( ) 27 2 3 4 max max 0 τ γ γ τ γ γ γ γ γ II = c ∫ = c − + G h d h (215) And therefore, I Ic G →G as ε →1, and II IIc G →G as γ →1 giving I c c G h d h C max max 0 16 ( ) 9 max σ ε ε σ ε ε = ∫ = (216) II c c G h d h C max max 0 16 ( ) 9 max τ γ γ τ γ γ = ∫ = (217) For triangular stressstrain tractionseparation law, a similar procedure can be applied. Where, c h is the cohesive layer thickness. It should be noted that the cohesive layer thickness is not arbitrary in cohesive layer model, and it is related to some 21 characteristic lengthscale of the debond process, such as COD. Determination of cohesive layer thickness from test data will be discussed in a later section. A phase angle is defined to describe the modemix of the failure (debond) in a cohesive layer ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ = − 2 1 * * tan 1 I II G G φ (218) where * I G and * II G are the strain energy release rates at which debond of cohesive layer initiates. In the case of the presence of antiplane shear stress (mode III), the mixed mode debond criterion can be expressed as + + = 1 IIIc III IIc II Ic I G G G G G G (219) and the phase angle between mode I and mode III is given by, ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ = − 2 1 * * tan 1 I III G G ϕ (220) 22 CHAPTER III DEVELOPMENT OF AN ANALYTICAL SOLUTION FOR COHESIVE LAYER MODEL AND MODEL VERIFICATION The objective of this chapter is to find an analytical solution to the cohesive damage zone at the interface between a fiber reinforced polymer (FRP) plate and concrete substrate. An analytical solution was derived to predict the stress in cohesive layer when considering the deformation in a stiff substrate. A twodimensional cohesive layer model with a prescribed stressstrain tractionseparation law as described in Chapter II was employed in this study. For comparison purpose, the cohesive layer model was implemented into a testbed finite element code (NOVA3D). Detailed benchmark comparisons of analytical results with finite element predictions for a double cantilever beam specimen for model verification were performed and issues related to cohesive layer thickness were investigated. It was observed that the assumption of a rigid substrate in analytical modeling can lead to inaccurate analytical prediction of cohesive damage zone length and stress distribution near debond tip. 3.1 Analytical Solution from Cohesive Zone Model Williams and Hadavinia [23] used a cohesive zone model with various simple forms of cohesive tractionseparation laws to analyze the global features and local stress 23 distribution of a double cantilever beam specimen (DCB, shown in Fig. 31). The DCB specimen is modeled as a cantilever beam with elastic foundation as shown in Fig. 32. The deformation of the beam is given by the equation EI w dx d v = 4 4 (31) y x P P Concrete substrate Cohesive layer Fig. 31 A cohesive layer in a double cantilever beam (DCB) a l Crack Damaged Zone Elastic Zone x1 x2 v Fig. 32 Schematic debonding of a double cantilever beam showing cohesive layer with elastic zone and damage zone ( 1 x ∈[0,l]) P 24 Where E is the Young’s modulus of the beam, v is the deflection of beam at its midplane, b is the width, and w is the distributed load per unit length of the beam that can be related to the stress in the cohesive layer by w = −bσ y . The stress y σ in the cohesive zone is modeled by a triangular elasticlineardamage tractionseparation law referring to the stress in the cohesive zone (as depicted in Fig. 22), where max v is the displacement at final fracture, max σ is the maximum stress in the cohesive zone at 3 max v v = . (a) In the damage zone: ( ) 2 3 max max max v v v b w b y = − = − σ σ (32) Let max 4 max 1 2 3 EIv bσ λ = (33) Thus Eq. (31) becomes ( ) max 4 4 1 4 v v dx d v = λ − (34) and the corresponding solution is, max 1 1 1 2 1 1 1 1 1 2 1 1 v = v + B sinhλ x + B coshλ x + C sinλ x + C cosλ x (35) (b) In the elastic zone: v v w b b y max max 3 σ = − σ = − (36) Let 25 max 4 max 2 4 3 EIv bσ λ = (37) v dx d v 4 4 2 4 = −4λ (38) and the solution is given by, ( sin cos ) 1 2 2 2 2 2 v = e−λ2x2 A λ x + A λ x (39) There are a total six unknown coefficients in Eqs. (35, 39). Along with the unknown damage length l and the critical force cr P , it requires eight boundary conditions to determine the beam deformation and the corresponding stresses in the cohesive layer. Two force boundary conditions are provided by force and moment equilibrium at the crack tip 0 1 x = . Continuity conditions imposed at the boundary of the elastic zone and damage zone yield another four boundary conditions by matching v , dx dv , 2 2 dx d v , 3 3 dx d v at x = l 1 and 0 2 x = . Finally, two displacement boundary conditions are max v = v at 0 1 x = and 3 max v v = at x = l 1 . From 1 max v (0) = v , 0 2 2 B + C = (310) EI P a dx d v (0) = cr 2 1 2 , where a is the crack length 2 1 2 2 2EIλ B + C = Pcra (311) EI P dx d v (0) = cr 3 1 3 26 3 1 1 1 2EIλ B − C = Pcr (312) 3 ( ) max 1 v v l = 3 sinh cosh sin cos max max 1 1 2 1 1 1 2 1 v v + B λ l + B λ l + C λ l + C λ l = (313) 3 (0) max 2 v v = 3 max 2 A = v (314) 1 ( ) 2 (0) dx dv l dx dv = ( cosh sinh cos sin ) ( ) 1 1 1 2 1 1 1 2 1 2 1 2 λ B λ l + B λ l + C λ l − C λ l = λ A − A (315) ( ) (0) 2 2 2 2 1 2 dx d v l dx d v = 1 2 1 1 2 1 1 1 2 1 2 2 1 λ (B sinhλ l + B coshλ l − C sinλ l −C cosλ l) = −2λ A (316) ( ) (0) 3 2 3 3 1 3 dx d v l dx d v = ( cosh sinh cos sin ) 2 ( ) 1 2 3 1 1 2 1 1 1 2 1 2 3 1 λ B λ l + B λ l − C λ l + C λ l = λ A + A (317) Eqs. (310 ~ 317) are nonlinear in terms of damage length l. Solutions are sought by an iterative numerical predictorcorrector method as follows: Damage length l is varied from 0.0 to 1.0 mm with interval of 0.001~0.01mm for the current specimen, and thereby Eqs. (310 ~ 316) become a set of linear equations that can be solved to obtain the constants coefficients 1 A , 2 A , 1 B , 2 B , 1 C , 2 C and the crack initiation load cr P for a specified damage length l. These constants and corresponding damage length l are then substituted into Eq. (317) and, because the equation is not exactly solved, the solution 27 error is numerically estimated. The correct solution for damage length l and the constants is the one that minimizes the error. 3.2 Analytical Solution from Cohesive Layer Model In Williams’ solution, deformation of the beam transverse to beam axis is neglected, therefore the displacement in ydirection (v) at the centerline of the beam is considered as the deformation of the cohesive zone. Maximum deformation max v at the debond tip, which is independent of the geometry and material properties of the beam, is the characteristic length scale of the cohesive zone model. When applying cohesive zone approach to a cohesive layer model, maximum deformation can be expressed as max max ε c v = h , where max ε is the maximum strain in the cohesive layer at failure, which is assumed to be a material property. As a result, the cohesive layer thickness c h is no longer arbitrary, but is uniquely determined by the relation max max ε h v c = . When considering the transverse deformation of the beam, c h is an important factor in evaluating the relative stiffness of beam and the cohesive layer in the transverse direction. The Young’s modulus of the beam (concrete) is usually much higher than that of the cohesive layer. On the other hand, the thickness of the beam is also much greater than that of the cohesive layer. Consequently the transverse deformation of the beam is comparable to the deformation in the cohesive layer and thus cannot be neglected. A small lateral deformation of the beam will greatly change the stress in the cohesive layer. FEA results clearly show the difference between the deflection of the beam at the 28 centerline and the deformation of the cohesive layer (that is, the displacement at the beamcohesive layer interface represents the deformation of the cohesive layer) (Fig. 3 3). As can be seen in Fig. 33 the displacement at the interface is generally smaller than the displacement at the centerline of the beam, the latter being the summation of the transverse deformations of the beam and the cohesive layer. Displacement ( mm ) X coordinate ( mm ) 4.0E04 2.0E04 0.0E+00 2.0E04 4.0E04 6.0E04 8.0E04 1.0E03 1.2E03 0 0.5 1 1.5 2 2.5 3 Interface Beam center line Fig. 33 Transverse deformation comparison between beam centerline and cohesive interface using FEA Thus, the total displacement v in ydirection at the centerline of the beam is composed of two parts: deformation b v of the beam and deformation c v of the cohesive layer, such that b c v = v + v (318) Assume that the transverse stress at the interface and in the cohesive layer is y c σ =σ , and = 0 y σ at the free top surface of the beam. A distribution law for y σ 29 through the height of the beam ( b h ) must be assumed to calculate the deformation of the beam (see Fig. 34). The actual distribution law for transverse stress can be obtained from elastic FEA analysis. Three idealized distribution laws for y σ were evaluated in this study: beam cohesive layer hc hb y stres s σc hb 0 Fig. 34 Transverse stress distribution in the beam and cohesive layer (a) Linear: (h y) h b b c y = − σ σ (b) Quadratic: 2 2 (h y) h b b c y = − σ σ (319) (c) Cubic: 3 3 (h y) h b b c y = − σ σ The corresponding lateral deformations in the lower half of the beam are, respectively: (a) Linear: b h h c b b b b c b h y b y E h y dy h E h dy E v dy b b σ b σ σ ε 8 ( ) 3 2 0 2 0 2 0 = ∫ = ∫ = ∫ − = (b) Quadratic: b h h c b b b b c b h y b y E h y dy h E h dy E v dy b b σ b σ σ ε 24 ( ) 7 2 0 2 0 2 2 2 0 = ∫ = ∫ = ∫ − = (320) 30 (c) Cubic: b h h c b b b b c b h y b y E h y dy h E h dy E v dy b b σ b σ σ ε 64 ( ) 15 2 0 2 0 3 3 2 0 = ∫ = ∫ = ∫ − = where b E is the Young’s modulus of the beam, and b h is the height of the beam. The above results can be generalized as b c b b E h v k σ = , which represents the lateral deformation of the beam, where k is the coefficient determined by the distribution of y σ along the ydirection within the beam given by Eq. (320). From Eq. (320), the values of k are 8 3 , 24 7 , 64 15 for linear, quadratic and cubic distributions, respectively. The analytical solution to the DCB specimen bonded by a cohesive layer with small but finite thickness can now be derived as follows. (a) In the damage zone: From Eq. (32), max max max 3 2 σ σ v v v c y = − (321) Displacement at the centerline of the beam is ⎟ ⎟⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = + = − − b b b c y E v v v v v k h max max max 3 2 σ σ (322) ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ − = − = − b b b b b b y E E I v k h v v b E I b dx d v max max 4 max 4 3 2 ( ) σ σ (323) By defining ⎟ ⎟⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − = b b b b E E I v k h b max max 4 1 3 2 σ λ (324) 31 The governing equation becomes ( ) max 4 4 1 4 v v dx d v =λ − (325) (b) In the elastic zone: ⎟ ⎟⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = + = + b b b c y E v v v v k h max max 3σ σ (326) ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ + = − = − b b b b b y E E I v k h v b E I b dx d v max max 4 4 3 4 4 σ σ (327) Defining, ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ + = b b b b E h k v E I b max max 4 2 3 4 σ λ (328) v dx d v 4 4 2 4 = −4λ (329) Governing Eqs. (325, 329) are of a similar form with the original Williams’ Eqs. (34, 38). But the coefficients 1 λ and 2 λ are different because of the presence of the extra term b b E k h which represents the transverse deformation of the beam (note that k = 0 gives the original Williams’ solution). A solution procedure and boundary conditions similar to the ones used for solving Eqs. (34, 38) can be employed to solve Eqs. (325, 329) to determine the unknown constants. Regarding the boundary conditions for Eqs. (325, 329), the two force boundary conditions and four displacement continuity conditions are same as the original solution 32 because these conditions are directly related to the beam. But the two displacement conditions, which are related to the cohesive layer, are a little different. First, at x = l 1 and 0 2 x = , the deformation of cohesive layer c v (not v at the centerline of the beam) equals 3 max v and the corresponding stress is max σ . 3 max v v v vc b = − = (330) Therefore the corresponding Eqs. (313, 314) become 3 ( ) max 1 v v l c = 3 ( sinh cosh sin cos ) max max max 1 1 2 1 1 1 2 1 v E h v B l B l C l C l k b + + + + − b = σ λ λ λ λ (331) 3 (0) max 2 v v c = 3 max max 2 v E h A k b − b = σ (332) On the other hand, at 0 1 x = , transverse stress = 0 y σ in the concrete beam, so that the lateral deformation of beam = 0 b v , therefore max v v v v v b c c = + = = which is same as Eq. (310). 3.3 Comparison between Analytical Solution and FEA Results A double cantilever beam (DCB) consisting of two substrates bonded with a thin layer of epoxy adhesive is shown in Fig. 31. The entire layer of epoxy is modeled with 33 special cohesive layer elements which obey the prescribed stressstrain tractionseparation law. A cohesive layer thickness h mm c = 0.02 and beam height h mm b =1 were used in this example, and the initial debond length a =15 mm corresponds to the unbonded portion of the beam. Unlike the usual linear elastic fracture mechanics modeling method of a sharp crack tip of zero tip radius, the localized fracture process zone in the current study has a small but finite thickness c h . Substrate beam is modeled as linear elastic, and the material properties used in this analysis are listed in Table 31. Table 31 Material properties for concrete and epoxy adhesive Concrete Epoxy Young’s modulus ( GPa ) 27.5 3.85 Shear modulus ( GPa ) 11.0 1.54 Poisson’s ratio 0.25 0.25 max σ ( MPa )  30.0 max ε  0.0526 Fig. 35 Finite element mesh of a DCB specimen with symmetry boundary conditions 34 Considering the symmetry boundary conditions along the centerline of epoxy adhesive layer, only half of the DCB specimen is modeled. Fig. 35 shows the actual finite element mesh (8node quadrilateral element) with applied symmetry boundary conditions. A very fine mesh is used to model the sharp stress gradients in the damage region. Plane strain conditions are assumed. Convergence study shows that converged results were obtained when the geometries of the elements in the crack tip area are smaller than 10 1 of the length of the damage zone. A constant beam tip displacement (instead of applied force) is specified in the FEA analysis to simulate a wedge test and to ensure that stable debond growth would occur. Debond growth (or cohesive layer element failure) is characterized by the transverse mechanical strain in the cohesive elements exceeding the specified maximum strain max ε in the cohesive layer beyond which the transverse stress goes to zero as defined by the cohesive constitutive law depicted in Fig. 22. Instead of node release and element deletion schemes used in most finite element codes, a failed element remains active in the subsequent analysis while the stiffness of the element is reduced close to zero. Analytical results with different k values are obtained and compared with numerical results from FEA as shown in Fig. 36. It can be seen that the best agreement with the FEA results is reached with the correction coefficient 64 k = 15 , which represents a cubic distribution of transverse stress in the beam. This cubic stress distribution is also verified by FEA simulation of the stress field in the concrete beam. The length of damage zone l decreases as the k increases, which implies a greater transverse deformation in the 35 beam under the same stress at the interface of beam and the cohesive layer. Due to the complexity of the transverse stress distribution, especially near the crack tip, k may take different values under different combination of material properties and geometries of the concrete beam and the cohesive layer. Further, as shown in Fig. 37, the deviation of the FEA results from the analytical solutions within the damage zone is likely attributable to the use of eightnode quadrilateral elements with quadratic interpolation that could result in linear variation of throughthickness strain within the fracture localization zone. Fortunately, the stress distribution in the damage zone does not have a significant effect on the process of debond initiation and propagation in the cohesive layer. 2.E+07 1.E+07 5.E+06 0.E+00 5.E+06 1.E+07 2.E+07 2.E+07 3.E+07 3.E+07 4.E+07 0 0.5 1 1.5 2 Unmodified k=0 Modified k=15/64 Modified k= 7/24 Modified k= 3/8 NOVA3D Fig. 36 Stress distribution with different correction coefficient k X coordinate (mm) = 0.02 b c h h Transverse stress (Pa) 36 Furthermore, different cohesive layer thickness and debond lengths were used to verify the agreement between analytical solutions and FEA simulation results, for both global and local metrics as discussed in the following paragraphs. Fig. 37 Reaction force comparison for different normalized cohesive layer thicknesses To study the effect of cohesive layer thickness on a global metric such as critical force (Pcr), the comparisons are conducted under the assumption that the cohesive layer deformation at the debond tip is equal to the maximum strain max ε of the cohesive layer (critical strain). For most practical cases, the ratio of cohesive layer thickness to beam thickness ( b c h h ) is in the range of 0.02 ~ 0.15. Critical reaction force versus normalized cohesive layer thickness ( b c h h ) is studied and good agreement is obtained between 37 analytical and FEA solutions. As shown in Fig. 37, the critical reaction force by FEA is slightly lower than the analytical predictions. This error is probably due to the fact that at the crack tip ( x1 = 0 ) the transverse displacement at the interface is actually slightly greater than that at the beam centerline (see Fig. 33), thereby reducing the reaction force required. To study the effect of debond length on reaction force, two kinds of boundary conditions are used: (a) Deformation control: under the same critical deformation at debond tip for FEA and analytical solution (b) Displacement control: under the same free end displacement for FEA and analytical solution. Good agreement is observed in both cases as shown in Fig. 38. 200 300 400 500 600 700 800 900 10 15 20 25 Debond length ( mm ) Reaction force ( N ) Analytical Deformation control (FEM) Displacement control (FEM) Fig. 38 Reaction force comparison for different debond lengths = 0.02 b c h h 38 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 0.05 0.1 0.15 0.2 0.25 Williams Modified Williams FEM Fig. 39 Damage length comparison for different cohesive layer thicknesses For local stress verifications, cohesive layer stress distribution in crack tip region is determined by evaluating damage length l. It can be observed from Fig. 39 and from Table 32 that the results from FEA are in good agreement with the modified Williams’ analytical solution when the normalized thickness lies between 0.02 and 0.15. When > 0.20 b c h h , that is, when the thickness of cohesive layer is relatively large, the deformation of the beam is relatively small compared with the deformation of the cohesive layer. As a result, the influence of the correction factor is no longer significant and the FEA solution approaches the original (unmodified) Williams’ solution. Fig. 310 shows how cohesive layer thickness c h influences local transverse stress distribution near crack tip as characterized by the damage length l. Good agreement is observed between analytical prediction and FEA results. Normalized thickness hc/hb Damage length ( mm ) 39 Fig. 310 Stress distributions for different cohesive layer thicknesses Table 32 Variation of damage length with different cohesive layer thicknesses Damage length ( mm ) b c h h Analytical NOVA3D Error ( % ) 0.02 0.186 0.190 2.15 0.05 0.336 0.332 1.19 0.10 0.446 0.437 2.02 0.15 0.513 0.534 4.09 0.20 0.562 0.587 4.45 It was observed that the modification to Williams’ model has little effect on global metrics (i.e. critical force and free end displacement), while for localized damage zone length and local stress distribution the effect of different values of k is significant as 40 depicted in Fig. 36, Fig. 39 and Table 33. It can be concluded that value of k has significant effect on the damage length, and 64 k = 15 results in the best agreement with FEA for the cohesive layer model with normalized thickness ranging from 0.02 to 0.15. For smaller or larger cohesive layer thickness, the relative lateral deformation of cohesive layer and concrete beam will change significantly. Further analytical studies are necessary to address these specific conditions. Table 33 Analytical and numerical solutions of damage length Correction factor k 0 15/64 7/24 3/8 FEA Damage length l ( mm ) 0.328 0.186 0.156 0.116 0.190 3.4 Conclusions An analytical solution was derived by introducing a correction term into the original Williams’ solution to predict the transverse stress in a cohesive layer when considering the deformation in a stiff substrate. Implementation of the cohesive layer model into a testbed finite element code was carried out and code verification was performed. Benchmark comparisons of finite element prediction of both global critical load and local stress field with analytical results for a DCB specimen resulted in good agreement after modifications were made to the original Williams’ solution. A sensitivity study was conducted to evaluate the influence of cohesive layer thickness on local parameters such as damage zone length, and global parameter such as critical force. From 41 the present studies it can be concluded that both local and global cohesive layer parameters are fairly sensitive to the cohesive layer thickness, whereas the correction factor (k) to Williams’ original solution significantly influences the local stress distribution and damage length. 42 CHAPTER IV HYGROTHERMAL EFFECT ON COHESIVE LAYER The objective of this chapter is to model the synergistic bond degradation mechanism that may occur at the interface between a fiber reinforced polymer (FRP) and concrete. For this purpose, a twodimensional cohesive layer involving straindependent, nonFickian hygrothermal effects as well as diffusion induced degradation in the cohesive layer was constructed. The model was implemented in a testbed finite element code (NOVA3D). Results from demonstration cases involving synergistic bond degradation were presented. 4.1 Introduction The influence of moisture diffusion on crack growth along an interface is not yet fully understood. Environmental cracking in a polymer typically occurs in the presence of a penetrant, such as moisture, and stress. It has been postulated that the mechanism involved in environmental crack growth in a polymer involves a small zone of craze formation and/or plasticization at the crack tip. For thermoset resins, such as epoxy, energy absorption at the crack tip is primarily by a shear yielding process and not by crazing. Consequently, for a thermoset epoxy, the zone of plasticization ahead of the crack tip must be determined using a diffusion law for nonporous media, such as Fick’s 43 law. However, quite frequently, polymer composites exhibit deviations from the classical Fickian treatment, termed as anomalous or nonFickian diffusion, especially at elevated temperatures and stress levels, and at high relative humidity. Some researchers have suggested that the deviation can be explained by a twostage Fickian process [24, 25]. Others claim that the diffusion process in a polymer is really nonFickian [26, 27]. Sophisticated hygrothermal models have been developed and verified by Roy [1720] to account for anomalous diffusion. The objective of this chapter is to model the synergistic bond degradation mechanisms that might occur at the interface due to interactions between stress, cohesive damage, and penetrant diffusion. For this purpose, a twodimensional cohesive layer constitutive model with a prescribed tractionseparation law is constructed from basic principles of continuum mechanics and thermodynamics, taking into account nonFickian hygrothermal effects that are likely to occur within the cohesive layer. 4.2 Fickian Diffusion Diffusion is the process by which matter is transported from one part of a system to another as a result of random molecular motions. In 1855, Fick first put diffusion on a quantitative basis by adopting the mathematical equation of heat conduction. The theory of diffusion in isotropic substance is therefore based on the hypothesis that the rate of transfer of diffusing substance through unit area of a section is proportional to the concentration gradient measured normal to the section, i.e. x f D C ∂ = − ∂ (41) 44 where C is the concentration of diffusing substance, D is called diffusion coefficient or diffusivity, and is the function of coordinate x, y, and z (location) and concentration C. Considering the mass balance of diffusing substance, we have = 0 ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ z f y f x f t C x y z (42) ⎟⎠ ⎞ ⎜⎝ ⎛ ∂ ∂ ∂ + ∂ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ ∂ + ∂ ⎟⎠ ⎞ ⎜⎝ ⎛ ∂ ∂ ∂ = ∂ ∂ ∂ z D C y z D C x y D C t x C (43) The fundamental differential equation of diffusion in an isotropic medium (where D is independent of the concentration C and location) can be expressed as ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ∂ + ∂ ∂ + ∂ ∂ = ∂ ∂ ∂ 2 2 2 2 2 2 z C y C x D C t C (44) and this is the socalled Fickian diffusion. 4.3 NonFickian Diffusion (Strain Assisted Diffusion) For a twodimensional cohesive layer of finite thickness c h , under planestrain conditions as shown in Fig. 41, the Helmholtz free energy per unit volume is given by, ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 4 18 12 3 17 22 12 2 12 2 16 22 12 3 15 22 4 14 22 3 13 12 2 12 12 22 12 2 11 22 3 7 11 12 8 11 22 9 12 22 10 22 2 6 12 2 5 22 2 0 1 11 2 22 3 12 4 11 , , , , , , , , , , , , , , , , , , , ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε ρψ ε ε ε ε ε C m T C m T C m T C m T C m T C m T C m T C m T C m T C m T C m T C m T C m T C m T C m T C m T C m T C m T C m T + + + + + + + + + + + + + = + + + + + (45) where the mechanical strain components are defined as, 45 ( ) ( ) 11 11 REF REF ε = E −α T − T −β m − m ( ) ( ) 22 22 REF REF ε = E −α T −T −β m − m 12 12 ε = E and, ρ : mass density of material in the cohesive layer 11 ε : mechanical strain component in x direction 22 ε : mechanical strain component normal to crack face (in y direction) 12 ε : shear strain component tangential to crack face ij E : total (kinematic) strain components m : moisture concentration in the cohesive layer at time t REF m : reference moisture concentration T : temperature in the cohesive layer at time t REF T : reference temperature α (T) : linear coefficient of thermal expansion β (T) : linear coefficient of moisture expansion moisture diffusion y x cohesive layer concrete beam Fig. 41. A cohesive layer with moisture diffusion in a DCB beam 46 From reduced entropy inequality, the cohesive stresses are defined by, ( , ) 2 ( , ) ( , ) ( , ) 1 4 11 7 12 8 22 11 11 σ ρ ψ C m T C m T ε C m T ε C m T ε E = + + + ∂ = ∂ (46) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 22 2 5 22 8 11 9 12 10 22 22 2 3 11 22 12 12 12 14 22 2 2 3 15 22 12 16 22 12 17 12 , 2 , , 3 , 2 , , 4 , 3 , 2 , , C mT C C mT C mT C mT E C mT C mT C mT C mT C mT C mT σ ρ ψ ε ε ε ε ε ε ε ε ε ε ε ε ε = ∂ = + + + + ∂ + + + + + + (47) Similarly, ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )3 18 12 2 17 22 12 12 2 16 2 3 15 22 2 12 22 12 13 12 2 3 6 12 7 11 9 22 11 22 12 12 3 , 4 , 2 , 3 , , 2 , , 2 , , , , ε ε ε ε ε ε ε ε ε ε ε ε ε σ ρ ψ C m T C m T C m T C m T C m T C m T C m T C m T C m T C m T C m T E + + + + + + = + + + + ∂ = ∂ (48) Chemical potential of the diffusing vapor is defined by Weistman [28], ∂m μ = ρ ∂ψ or, ⎥⎦ ⎤ ⎢⎣ ⎡ − ∂ ∂ + ⎥⎦ ⎤ ⎢⎣ ⎡ − ∂ ∂ + ⎥⎦ ⎤ ⎢⎣ ⎡ − ∂ ∂ + ⎥⎦ ⎤ ⎢⎣ ⎡ − ∂ ∂ + ∂ ∂ = 3 14 22 4 22 2 14 10 22 3 22 10 5 22 2 22 5 4 11 2 11 0 4 3 ( , ) ( ) 4 ( , ) ( ) 2 ( , ) ( ) 2 ( , ) ( ) ε β ε ε β ε μ ε β ε ε β ε C m T T m C m T T C m C C m T T m C m T T C m C m C (49) From conservation of mass, the governing Equation for twodimensional moisture diffusion is, x y m f f t x y ∂ ⎛ ∂ ∂ ⎞ = −⎜ + ⎟ ∂ ⎝ ∂ ∂ ⎠ (410) 47 where the moisture flux, f = fxnˆx + fynˆy r , in the absence of temperature gradients is given by, ˆ ˆ x y f D x f D y μ μ = − ∂ ∂ = − ∂ ∂ (411) Assuming isotropic material and using the chain rule, 22 22 ˆ ˆ ˆ x f D m D T D m x T x x μ μ μ ε ε ⎛ ∂ ⎞ ∂ ⎛ ∂ ⎞ ∂ ⎛ ∂ ⎞ ∂ = − ⎜ ⎟ − ⎜ ⎟ − ⎜ ⎟ ⎝ ∂ ⎠ ∂ ⎝ ∂ ⎠ ∂ ⎝ ∂ ⎠ ∂ (412) 22 22 ˆ ˆ ˆ y f D m D T D m y T y y μ μ μ ε ε ⎛ ∂ ⎞ ∂ ⎛ ∂ ⎞ ∂ ⎛ ∂ ⎞ ∂ = − ⎜ ⎟ − ⎜ ⎟ − ⎜ ⎟ ⎝ ∂ ⎠ ∂ ⎝ ∂ ⎠ ∂ ⎝ ∂ ⎠ ∂ (413) Assuming isothermal condition and substituting Eqs. (49, 412, 413) into Eq. (4 10), gives 22 22 m m m D m D D m D t x x ε x y y ε y ∂ = ∂ ⎛ ∂ + ∂ε ⎞ + ∂ ⎛ ∂ + ∂ε ⎞ ∂ ∂ ⎜⎝ ∂ ∂ ⎟⎠ ∂ ⎜ ∂ ∂ ⎟ ⎝ ⎠ where ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ ∂ = ∂ ⎟⎠ ⎞ ⎜⎝ ⎛ ∂ = ∂ 22 ˆ , ˆ ε μ μ ε and D D m D D m ⎥⎦ + ⎤ ∂ ∂ − ∂ ∂ + + ∂ ∂ − ∂ ∂ + + ∂ ∂ − ⎢⎣ ⎡ ∂ ∂ + + ∂ ∂ − ∂ ∂ + ∂ ∂ = 2 22 2 14 3 22 14 4 2 22 14 2 22 2 10 2 22 3 10 2 22 10 2 2 5 22 2 5 2 22 5 2 2 11 4 2 4 2 11 4 2 2 0 2 8 ( ) 6 ( ) ( ) 2 ( ) 6 ( ) 6 ( ) ( ) ˆ 4 ( ) 2 ( ) 4 ( ) β ε β ε β ε β ε β ε ε ε β ε β ε β ε T C T T m C m T C T T C m C m C T C T m C m C T C T m C m C m C D D m (414a) 48 ⎥⎦ − ⎤ ∂ ∂ + ⎢⎣ ⎡ − ∂ ∂ − + ∂ ∂ = 2 14 22 3 22 14 10 22 2 22 10 22 5 5 4 12 ( ) ( ) ˆ 2 2 ( ) ( ) 3 6 ( ) ( ) ε β ε ε β ε β ε ε T C T m C T C T m C T C T m C D D (414b) If the unknown material coefficients are expanded in a Taylor Series about a reference moisture concentration value REF m , and retaining terms up to second order in change in moisture concentration, ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 2 2 , , ˆ , 1 , , ( 0, 1,......,18) REF REF k k k k REF REF REF m m k REF k REF k REF C m T C m T C m m C m m m m C m T C m T m C m T m k ⎛ ∂ ⎞ ⎛ ∂ ⎞ = + ⎜ ⎟ − + ⎜ ⎟ − ⎝ ∂ ⎠ ⎝ ∂ ⎠ = ⎡⎣ + Δ + Δ ⎤⎦ = % (415) In order to benchmark the present model against an established cohesive zone model, some of the unknown material coefficients were determined by assuming a cubic tractionseparation law similar to the one proposed by Needleman [7] and modified by ElSayed and Sridharan [14] for a finitethickness cohesive layer, giving, ( ) ( )[ ( ) ( ) 2 ] 0 0 0 0 C m,T Cˆ m ,T 1 C m ,T m C~ m ,T m REF REF REF = + Δ + Δ ( , ) 0 1 C m T = ( , ) 0 2 C m T = ( , ) 0 3 C m T = ( ) ( )[ ( ) ( ) 2 ] 4 4 4 C m,T T 1 C m ,T m C~ m ,T m REF REF MAX = MAX + Δ + Δ ε σ ( ) ( )[ ( ) ( ) 2 ] 5 5 5 1 , ~ , 8 C m,T 27 T C m T m C m T m REF REF MAX = MAX + Δ + Δ ε σ 49 ( ) ( )[ ( ) ( ) 2 ] 6 6 6 1 , ~ , 8 C m,T 27 T C m T m C m T m REF REF MAX = MAX + Δ + Δ γ τ ( ) ( ) ( ) ⎪⎭ ⎪⎬ ⎫ = = = , 0 , 0 , 0 9 8 7 C m T C m T C m T assume ν = 0 ( ) ( )[ ( ) ( ) 2 ] 10 2 10 10 1 , ~ , 2 C m,T 9 T C m T m C m T m REF REF MAX = − MAX + Δ + Δ ε σ ( , ) 0 11 C m T = ( , ) 0 12 C m T = ( ) ( )[ ( ) ( ) 2 ] 13 2 13 13 1 , ~ , 2 C m,T 9 T C m T m C m T m REF REF MAX = − MAX + Δ + Δ γ τ ( ) ( )[ ( ) ( ) 2 ] 14 3 14 14 1 , ~ , 16 C m,T 27 T C m T m C m T m REF REF MAX = MAX + Δ + Δ ε σ ( , ) 0 15 C m T = ( , ) 0 16 C m T = ( , ) 0 17 C m T = ( ) ( )[ ( ) ( ) 2 ] 18 3 18 18 1 , ~ , 16 C m,T 27 T C m T m C m T m REF REF MAX = MAX + Δ + Δ γ τ Substituting these definitions in expressions for , , , 11 22 12 σ σ σ (Eqs. 46, 47, and 48) results in, ( ) ( ) ( ) MAX MAX C T E T T ε σ = = 4 2 ˆ gives ( )[ ( ) ( ) ] 11 2 11 4 4 σ σ T 1 C m ,T m C~ m ,T m ε MAX REF REF = + Δ + Δ (416) 50 Similarly, ( )[ ( ) ( ) ] ( )[ ( ) ( ) ] ( )[ ( ) ( ) ] 3 22 2 14 14 14 2 22 2 10 10 10 22 2 22 5 5 5 4 ˆ 1 , ~ , 3 ˆ 1 , ~ , 2 ˆ 1 , ~ , ε ε σ ε C T C m T m C m T m C T C m T m C m T m C T C m T m C m T m REF REF REF REF REF REF + + Δ + Δ + + Δ + Δ = + Δ + Δ or ( ){[ ( ) ( ) ] [ ( ) ( ) ] [ ( ) ( ) ] 3 } 22 2 14 14 2 22 2 10 10 22 2 22 5 5 1 , ~ , 21 , ~ , 1 , ~ , 4 27 ε ε σ σ ε C m T m C m T m C m T m C m T m T C m T m C m T m REF REF REF REF MAX REF REF + + Δ + Δ − + Δ + Δ = + Δ + Δ (417) where the maximum tensile stress in the cohesive layer is MAX σ and characteristic interface length δ is related to the layer thickness c h by, MAX c h ε = δ . And finally the shear stress can be expressed as, ( ){[ ( ) ( ) ] [ ( ) ( ) ] [ ( ) ( ) ] 3 } 12 2 18 18 2 12 2 13 13 12 2 12 6 6 1 , ~ , 21 , ~ , 1 , ~ , 4 27 ε ε σ τ ε C m T m C m T m C m T m C m T m T C m T m C m T m REF REF REF REF MAX REF REF + + Δ + Δ − + Δ + Δ = + Δ + Δ (418) where MAX MAX MAX γ ε ε ε ε ε ε ε ε 12 12 22 22 11 11 = , = , = , the maximum shear stress in the cohesive layer is MAX τ . The consistent diffusivities m D and ε D for the cohesive layer are obtained by substituting the definitions of the material coefficients into Eqs. (414a, 414b), 51 { ( ) ( )[ ] ( )[ ] ( )[ ] ( )[ ] ⎭ ⎬ ⎫ + − + Δ + + Δ + Δ − − + Δ + + Δ + Δ + − + Δ + + Δ + Δ + − + Δ + + Δ + Δ = 2 22 2 14 14 3 2 14 14 22 4 3 14 22 22 2 10 10 2 2 10 10 22 3 2 10 22 2 5 5 2 5 5 22 2 5 22 2 4 4 2 4 4 11 2 4 11 0 0 2 ~ 4 ( 2 ~ ) 2 (1 ~ ) 16 27 2 ~ 4 ( 2 ~ ) 2 (1 ~ ) 2 9 2 ~ 4 ( 2 ~ ) 2 (1 ~ ) 8 27 2 ~ 4 ( 2 ~ ) 2 (1 ~ ) ˆ 2 ˆ ~ ( ) ε β ε β ε ε σ ε β ε β ε ε σ ε β ε β ε σ ε β ε β ε σ T C C C m C m C m T C C C m C m C m T C C C m C m C m T C C C m C m C m D D C T C T MAX MAX MAX MAX MAX MAX MAX MAX m (419) ( )[ ] ( )[ ] ( )[ ] ⎭ ⎬ ⎫ + + Δ − + Δ + Δ − + Δ − + Δ + Δ ⎩ ⎨ ⎧ = + Δ − + Δ + Δ 2 22 2 14 14 3 3 14 14 22 22 2 10 10 2 2 10 10 22 2 5 5 22 5 5 ( 2 ~ ) 3 (1 ~ ) 4 27 ( 2 ~ ) 2 (1 ~ ) 2 27 ( 2 ~ ) (1 ~ ) 4 ˆ 27 ε β ε ε σ ε β ε ε σ ε β ε σ ε T C C m C m C m T C C m C m C m D D T C C m C m C m MAX MAX MAX MAX MAX MAX (420) 4.4 Cohesive Layer Diffusion Boundary Conditions Assuming that the chemical potential of the ambient vapor on the exposed boundary of the cohesive zone remains constant with respect to time [28], the resulting concentration at the boundary of the cohesive zone (crack tip) can be derived as, ( ) ij BOUNDARY b μ T ,m,ε = μ (421) 52 { ( )[ ] ( )[ ] ( )[ ] MAX BOUNDARY MAX MAX MAX MAX MAX b T C C m C m C m T C C m C m C m T C C m C m C m C C C m ⎭ ⎬ ⎫ + + Δ − + Δ + Δ − + Δ − + Δ + Δ + + Δ − + Δ + Δ = + Δ 3 22 2 14 14 4 3 14 14 22 2 22 2 10 10 3 2 10 10 22 22 2 5 5 2 5 5 22 0 0 0 ( 2 ~ ) 4 (1 ~ ) 16 27 ( 2 ~ ) 3 (1 ~ ) 2 9 ( 2 ~ ) 2 (1 ~ ) 8 27 ˆ ( 2 ~ ) ε β ε ε σ ε β ε ε σ ε β ε ε σ μ An equation about the unknown moisture concentration Δm at the boundary can be derived in the form of aΔm2 + bΔm + c = 0 . Therefore the boundary concentration can be solved a m m b b ac b REF 2 = + − ± 2 − 4 (422) 4.5 Cohesive Layer WorkofSeparation For pure Mode I or Mode II fracture, the workofseparation at the cohesive interface (of finite thickness c h ) per unit volume in the presence of moisture concentration (Δm ) is given by ( ) {[ ] [ ] [ ] } ( ) [ ] [ ] [ ] ⎭ ⎬ ⎫ − + Δ + Δ + + Δ + Δ ⎩ ⎨ ⎧ = + Δ + Δ ⎭ ⎬ ⎫ − + Δ + Δ + + Δ + Δ ⎩ ⎨ ⎧ = + Δ + Δ = ⎥ ⎥⎦ ⎤ ⎢ ⎢⎣ ⎡ = ∫ ∫ ∫ ∫ 2 10 10 14 14 2 5 5 22 3 22 2 14 14 2 10 10 22 1 0 22 2 5 5 1 0 22 22 0 22 22 1 ~ 4 1 ~ 1 3 2 1 ~ 2 1 4 27 21 ~ 1 ~ 1 ~ 4 27 max C m C m C m C m T h C m C m C m C m C m C m d h T C m C m d dV h d MAX MAX c c MAX MAX c MAX V I sep σ ε ε ε ε σ ε ε φ σ ε σ ε ε ε (423) 53 ( ) { } } ( ) max 1 12 12 12 12 0 0 1 2 6 6 12 0 2 2 2 3 13 13 12 18 18 12 12 2 6 6 13 13 ) 27 1 4 2 1 1 27 1 1 4 2 2 1 3 II sep c MAX V c MAX MAX MAX MAX c d dV h d h T Cm C m C m C m C m C m d T h C m C m C m C m γ φ σ ε σ γ ε τ γ ε ε ε ε τ γ ⎡ ⎤ = ⎢ ⎥ = ⎢⎣ ⎥⎦ = ⎧ ⎡ + Δ + Δ ⎤ ⎨ ⎣ ⎦ ⎩ − ⎡⎣ + Δ + Δ ⎤⎦ + ⎣ + Δ⎤ + ⎡ Δ ⎦ = ⎧ ⎡ + Δ + Δ ⎤ ⎨ ⎣ ⎦ ⎩ − ⎡⎣ + Δ + Δ ∫ ∫ ∫ ∫ % % % % % 2 18 18 1 1 4 C m C m ⎫ ⎤ ⎡ ⎤ + + Δ + Δ ⎬ ⎦ ⎣ ⎦⎭ % (424) or write in the short form [ 2 ] 1 2 ( ) 1 ( ) ( ) 16 9 T h A T m A T m MAX MAX c I sep φ = σ ε + Δ + Δ (425) [ 2 ] 3 4 ( ) 1 ( ) ( ) 16 9 T h A T m A T m MAX MAX c II sep φ = τ γ + Δ + Δ (426) The thickness of the cohesive layer, c h , is a measure of the fracture localization zone and it is directly related to the characteristic length scale of the fracture process zone controlled by specific fracture mechanism [29]. The area under the normalized stressstrain curve ( max 0 ε 1 ε ≤ ≤ ) shown in Fig. 42 represents the workofseparation of the cohesive layer. From Eqs. (425, 426) and experimental results, it can be observed that the workofseparation decreases at elevated temperature and higher moisture concentration due to the physical or chemical degradation at the interface through a decrease in the maximum peel stress and a corresponding reduction in the area under the curve. When considering the timedependent behavior of the polymeric adhesive, workof separation is also a function of time. Lower workofseparation would imply lower critical strain energy release rate, and therefore, lower resistance to crack growth. 54 4.6 Cohesive Layer Degradation and FEA Simulation The proposed cohesive layer model was implemented in an inhouse finite element code, NOVA3D. In the previous section, the parameters i C and i C~ are material constants to be determined experimentally. Moisture diffusion tests and fracture experiments are necessary to characterize these coefficients. 4.6.1 Cohesive Layer Degradation due to Moisture Concentration When moisture diffuses into the cohesive layer, two effects occur concurrently. First, the cohesive layer begins to swell thereby causing the local stress state to change due to the constraining effect of the surrounding adherends. Secondly, it is likely that moisture will penetrate the bulk cohesive layer to reach the cohesive interface (or interphase), and then rapidly diffuse along the interface. At the interface, water molecules typically react with the chemical bonds across the interface especially in the presence of tensile stress, because stress provides additional driving force for the bond rupturing process. Such chemical reactions transform strong covalent bonds to weak Van der Waals bonds, thereby significantly weakening the interface strength and fracture toughness. The bondstrength degradation could be important even when the change in moisture concentration is relatively small (~10%) [16]. Due to a lack of available bond degradation data from ongoing experiments to allow characterization of material coefficients in Eq. (415) at present time, a simple 55 bond strength degradation scheme (as illustrated in Fig. 42) was employed for the present analysis while preserving the basic framework presented in Eq. (415). Fig. 42. Depiction of the influence of moisture in cohesive layer on workofseparation 4.6.2 FEA Simulation Results Fully coupled stress and diffusion analyses were invoked in this investigation, analogous to simulating a wedgetest under wet conditions. For the DCB specimen, moisture diffusion analysis is activated only in the cohesive layer, with moisture boundary conditions applied at the exposed surface at x = 0 as depicted in Fig. 41. A cohesive layer thickness h mm c = 0.02 and beam height h mm b = 1 were used, Material properties used in this analysis are listed in Tables 31, and the diffusivity of the cohesive layer is 5.22×10−8mm2 / s . With time, moisture gradually diffuses from the exposed debond tip ( x = 0 ) into the originally dry elements in the cohesive layer. Due to the lack of the material 56 coefficients defining nonFickian diffusivities given by Eq. (415), linear Fickian diffusion with constant boundary concentration (m 1 10 8 g /mm3 b = × − ) was modeled in the cohesive layer for this demonstration case. Fig. 43 depicts the moisture concentration profiles plotted along the bond length with the origin at the original location of the debond tip at four different time steps. The corresponding analytical solutions for onedimensional Fick’s law are also plotted in Fig. 43 for diffusion model verification. It is evident that the concentration profiles predicted by the finite element diffusion analysis are in excellent agreement with the analytical solution for the linear case. Incidentally, the horizontal dashed line in Fig. 43 corresponds to 10% of saturation concentration and its purpose will be discussed in the following paragraphs. 0.0E+00 2.0E09 4.0E09 6.0E09 8.0E09 1.0E08 1.2E08 0 0.5 1 1.5 2 8.0E4 sec. 1.0E6 sec. 3.7E6 sec. 1.0E9 sec. Fick's Law Fig. 43 Comparison of predicted moisture concentration profiles A phenomenological stepfunction degradation law is assumed in the analysis such that when the local moisture concentration is greater than or equal to 10% of the X coordinate (mm) Concentration (g/mm3) 57 saturation concentration corresponding to the dashed horizontal line in Fig. 43, the corresponding maximum stress ( max σ ) in the cohesive tractionseparation law is reduced by 10% . For numerical stability, a linear degradation is assumed when the concentration is between 0 and 10%. The length of the moistureinduced cohesive strength degradation zone at various time steps is indicated by the intersection of the concentration profile and the horizontal dashed line in Fig. 43. t = 8 ×104 sec debond damaged undamaged 0.16 mm beam cohesive layer (a) Under dry condition t = 3.7 ×106 sec 0.57 mm (b) Degradation due to moisture concentration Fig. 44 Debond growth in a cohesive layer due to moisture degradation In this demonstration case, a constant beam tip displacement is applied at all times for the DCB, simulating a wedge test. The debond growth predicted at the crack tip by 58 the finite element analysis is 0.16 mm at initial time ( t = 0 ) under dry conditions as shown in Fig. 44(a). Over time, as moisture diffusion and subsequent bond degradation takes place, debond propagation occurs and the failure length increases to 0.57 mm over a period of 3.7×106 seconds, as depicted in Fig. 44(b). Failure (or debond) length is determined by observing if the transverse mechanical strain, 22 ε , has exceeded the prescribed maximum transverse strain, max ε , along the bond length. Fig. 45 shows that the transverse mechanical strain monotonically increases with time due to the formation of cohesive damage, material failure, and resultant debond propagation. The location of the crack tip of the failure zone at various time steps is indicated by the intersection of the mechanical strain and the horizontal dashed line in Fig. 45 indicating failure strain, max ε (debond length increases from 0.16 mm to 0.57 mm due to strength degradation caused by moisture diffusion). 0 0.5 1 1.5 2 2.5 3 3.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 8.0E4 sec. 1.0E6 sec. 3.7E6 sec. Normalized strain ε y /ε max X coordinate ( mm ) Fig. 45 Transverse mechanical strains and corresponding failure zones at different moments 59 The reduction in corresponding transverse stress distribution at different time steps is plotted along bond length in Fig. 46. The progressive reduction in the peak stress magnitude is due to degradation of bond strength caused by moisture ingress. Evolution of debond length with time is plotted in Fig. 47, and the decrease in reaction force due to the degradation of cohesive layer stiffness with time is shown in Fig.48. It is evident that for the present case, debond growth is driven by a synergistic interaction of moisture diffusion and transverse stress near the debond tip as discussed in the following paragraph. 1.0E+07 0.0E+00 1.0E+07 2.0E+07 3.0E+07 0 0.5 1 1.5 2 2.5 3 8.0E4 sec. 1.0E6 sec. 3.7E6 sec. Stress ( Pa ) X coordinate ( mm ) Fig. 46. Transverse stresses along bond length at different moments There are two important milestones during the diffusionassisted debond growth process, 4 t1 = 8×10 second and 6 2 t = 3.7 ×10 second (see Fig. 47). For 1 0 ≤ t < t , the moisture concentration ahead of the crack tip is not high enough to result in significant 60 degradation to the cohesive layer. Therefore, no debond growth occurs during this period. At roughly t = t1 , the 10% concentration front impinges on the debond tip and triggers debond growth as shown in Fig. 46. For 1 2 t ≤ t ≤ t , the moisture concentration ahead of the debond tip becomes high enough (>10%) such that the onset of bond degradation occurs, resulting in steady debond growth as shown in Fig. 47. In this regime, the debond growth is driven by the rate of propagation of the 10% concentration front, and therefore can be said to be diffusioncontrolled. Finally, for 2 t > t , the debond driving force (shown in Fig. 48) falls below a threshold value such that no further synergistic debond growth is possible. In this case, the peak transverse stress in the cohesive layer has decreased to max 0.9σ as shown in Fig. 46 due to the assumed moisture degradation characteristics. At the same time the transverse deformation (Fig. 44) and debond length (Fig. 47) have reached steady state values, respectively. Sqrt( hour ) Debond length and 10% concentration front ( mm ) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0 10 20 30 40 50 60 Debond length 10% concentration front t1 t2 Fig. 47 Predicted evolution of debond length with time 61 Fig. 48 Reaction force decreases with time 4.7 Conclusions A twodimensional cohesive layer constitutive model with a prescribed tractionseparation law was constructed from basic principles of continuum mechanics and thermodynamics, taking into account concentrationdependent and straindependent non Fickian hygrothermal effects that are likely to occur within a cohesive layer. Implementation of the model in a testbed finite element code was carried out and code verification was performed. Numerical simulation of a wedgetest involving debond growth caused by synergistic interactions between local stress and diffusing moisture was also presented to demonstrate the ability of the cohesive layer model to simulate environmental cracking. Reaction force ( N ) Sqrt ( hour ) 62 CHAPTER V NUMERICAL SIMUILATION OF PENINSULA BLISTER TEST The objective of this paper is to examine threedimensional and viscoelastic effects in the peninsula blister test for thin film adhesion, which is considered an effective way to measure the interfacial fracture toughness of bonded thin film structures. As will be demonstrated in this Chapter, analytical solutions are sometimes inadequate for accurately simulating peninsula blister experiments because linear elastic behavior is assumed in the thin film as well as in the debonding process zone at the interface. For this purpose, a threedimensional cohesive layer model and corresponding liquid loading simulation algorithm were developed and implemented into an inhouse testbed finite element analysis (FEA) code (NOVA3D). Features such as threedimensional mixedmode debonding, large displacements and rotations, residual stresses in the thin film, and timedependent (viscoelastic) effects in the thin film were considered. Numerical convergence and a stable debond growth were obtained over a fairly large debond length when suitably refined FEA mesh and liquid volume increment were employed. Detailed benchmark comparisons of the finite element predictions with analytical solutions and experimental results were performed, and good agreement was obtained. 63 5.1 Introduction Laminated thin film structures appear in a wide variety of applications such as multilayer structures, microelectronic devices, and packages. The increasing application of thin films in industry has made it necessary to better understand the nature and the reliability of the bonding between layers or bimaterial interface. Since poor bonding results in a crack or delamination, fracture mechanics is a natural approach for characterizing the resistance to failure and making durability or reliability prediction. For this purpose, several experimental methods, including peel tests, blister tests, indentation tests, scratch tests etc., have been used to determine the interfacial bond strength or toughness. The peel test is a simple and perhaps the most commonly used method for examining the strength of adhesively bonded layers. In a peel test, a thin flexible strip that is bonded to a substrate is pulled off at a certain angle to the underlying substrate. In the absence of any plastic deformation and residual stresses, adhesive fracture energy Γ can be derived directly [30] from the peel force through Γ = (1− cosϕ )P (51) where P is the peel force per unit width of the film and ϕ is the peel angle. In reality, plastic deformation in the peeling arm and crack tip singularity are present in all the peel test specimens, and were investigated by Kim et al. [31], Kinloch et al. [32], and Wei and Hutchinson [33]. It implies that a significant portion of the mechanical energy applied in peeling is dissipated and, if not properly accounted for, will lead to significant overestimates of the adhesive fracture energy. This is particularly true 64 when the yield strength of the peeling film is relatively low and the degree of adhesion is high. The pressurized circular blister test was developed by Williams [34] as an alternative approach to minimize dissipation effects and overcome some of the other drawbacks of the peel test. In many blister tests, a pressurized fluid (usually water is used and considered as incompressible) is injected between the faces of the crack between the film and substrate via a hole through the substrate. A uniform fluid pressure is applied to the thin film and hence the problems associated with mechanical contact and gripping are avoided. Furthermore, the axisymmetry of the circular blister configuration greatly minimize edge effects and effects due to specimen nonuniformity. The standard blister configuration (circular blister) was first introduced by Dannenberg [35] in 1961 to measure the adhesion of thick organic coating to metal. As tougher interface were examined, problems were encountered with tensile failure in blister film prior to debond at the interface. Another disadvantage of circular blister is that the strain energy release rate increases as the fourth power of the debond radius, thus making accurate evaluation of the debond radius essential, and resulting in a very unstable fracture specimen. Several alternative configurations with higher strain energy release rates were developed. These configurations include island blister proposed by Allen and Senturia [36], constrained blister proposed by Napolitano et al. [37] and Chang et al. [38], and peninsula blister proposed by Dillard and Bao [39]. In a constrained blister specimen, even though the strain energy release rate is independent of debond radius, the effect of friction between the delaminating film and the constraint could be significant, raising concern about the difficulty of analyzing contact 65 problems with large deformation. In the island blister specimen, debond growth occurs radially inwards on the island. A moderate increase in compliance is produced by a relative small increase in debond area, thereby giving rise to high strain energy release rates. As the island radius decreases with debonding, the calculated strain energy release rates increase in an unbounded manner. Fig. 51 Peninsula blister specimen and possible debond sites The peninsula blister is an extension of the island blister concept. Its name is derived from the fact that debonding will occur along a narrow “peninsula” which extends into the blister region, as shown in Fig. 51. The peninsula blister configuration not only maintains the high energy release rate of the island blister but also provides a constant energy release rate as the adhesive debond progresses. Added features include the larger debond area and additional data points that can be obtained from a single specimen. A consistent set of analyses by Liechti and Shirani [40] showed that the only configuration among circular, island, and peninsula blister specimen with a homogeneous 2 1 3 L l 2a 2b 66 delaminating layer that did not suffer any yielding was a relatively thick peninsula blister configuration. A subsequent cohesive zone fracture analysis by Shirani and Liechti [41] found that only 6% of the work input was dissipated as global plastic dissipation in a polyimide thin film on aluminum, a relatively tough interface. Disadvantages of the peninsula blister test include the lack of symmetry and difficulties in the fabrication of the specimen for certain material systems. Despite its promising advantages, the peninsula blister configuration has not gained wide acceptance in the adhesion testing community. One problem is the relatively complex nature of the specimen fabrication although this can be alleviated by suitable masks. Another is the relatively complex nature of the analysis particularly when large deformations and rotations and residual stresses are encountered. Both issues have recently been addressed [42]. One issue that is yet to be resolved on the analytical front is an accounting of the viscoelastic nature of many adhesives, particularly under long term loadings [43]. On the other hand, problems were also encountered in numerical simulation of peninsula blister test, such as the loss of axisymmetric geometry, mixed mode I and mode II nonlinear fracture (and even mode III), large deflections and rotations in the film and adhesive layer, and large debond length. Another problem is the special loading methodology of the blister test. It employs a liquidvolumecontrol loading approach which is obviously different from any general applied force or displacement loading methods. The purpose of the work presented here was to consider optimum geometries for peninsula blister specimens while accounting for geometric nonlinearities, viscoelastic effects and mixedmode debonding. A threedimensional cohesive layer model and 67 corresponding liquid loading simulation algorithm were developed and implemented into an inhouse testbed FEA code (NOVA3D). Numerical convergence was achieved with suitable increment of liquid volume and with reasonable FEA mesh refinement. Steadystate debond growth was obtained over a fairly large portion of the specimen. The finite element results were benchmarked against analytical solutions and experimental data. 5.2 Constitutive Law and Failure Criterion of Cohesive Layer Peninsula blister specimen is a full threedimensional configuration. Threedimensional stressstrain tractionseparation law is adopted to simulate the cohesive layer as shown in Eqs. (23, 24). In peninsula blister specimen, shear force perpendicular to the peninsula direction (corresponding to mode III fracture) was observed in peninsula blister test and simulation. A more complex failure criterion that contains the contribution of mode III debond was used in this case. + + = 1 IIIc III IIc II Ic I G G G G G G (52) Detailed description of the 3D cohesive layer model and its failure criterion is given in Chapter II. 5.3 Specimen Geometry and Energy Release Rates Proper choices of the geometry and, in particular, the aspect ratios of the peninsula blister specimen and regions within it are essential in order to ensure that 68 debond initiation and propagation occurs in the 1direction along the narrow peninsula area as shown in Fig. 51. Debonding at sites 2 and 3 is also possible but is undesirable. From the standpoint of linear elastic fracture mechanics, the energy release rate can be used to predict crack initiation and propagation; it can be viewed as the driving force for crack propagation or debond growth. When the deflection of the film is on the order of the film thickness, the theory of linear elastic thin plates can be applied in this case. Dillard and Bao [39] determined the analytical solution of the strain energy release rates for debonding at the three sites under a given liquid pressure p [ 5 5 ] 2 1 16 ( ) 1440 a a b Db G = p − − (53) D G p a b 288 2 ( )4 2 = − (54) D G p a 288 2 4 3 = (55) where D is the bending stiffness of the film given by 12(1 2 ) 3 −υ D = Eh (56) E , υ , and h are its Young’s modulus, Poisson’s ratio, and thickness, respectively. Numerical results from FEA analysis (NOVA3D) are in very good agreement with Dillard’s analytical solutions as shown in Fig. 52 and 53. It can be seen from Fig. 53 that higher energy release rate G can be achieved at site 1 for a small peninsula width. Therefore, when the peninsula width is relatively small (e.g. b 0.2 a < ), debond will 69 occur preferentially at site 1 rather than site 3 (because G1 > G3 ). It should also be noted that 1 G is greater than 2 G for all b a ratios, implying debond of the peninsula would always propagate along the peninsula length rather than across the width. As evident from Fig 52, in order to obtain large energy release rates with the peninsula specimen, the most effective way is to select small values of b a . It is suggested that b 0.1 a = is likely to be the most practical choice, as overly small b a will increase difficulties in fabricating specimen. 0.5 0 0.5 1 1.5 2 2.5 3 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Normalized peninsula width b/a Strain energy release rate G*D/P/P Analytical G1  FEA G2  FEA G3  FEA Fig. 52 Effect of relative peninsula width on the energy release rate Analytical G1 – FEA G2 – FEA G3 – FEA 70 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Aspect ratio of the blister region l /2a Strain energy release rate G*D/P/P FEA b/a = 0.05 FEA b/a = 0.1 FEA b/a = 0.2 Analytical Fig. 53 Effect of blister region aspect ratio on the energy release rate (Analytical solution is derived from l 2a > 2 ) Aspect ratio of the blister region ( 2 l a ) also deserves examination. The derivation of strain energy release rates discussed earlier (Eqs. 53 ~ 55) was based on a plate with infinite length and clamped along the two remaining boundaries. According to Timoshenko and WoinowskyKrieger [44], for an aspect ratio of 2, the deviation in the midplane deflection of a fully clamped plate (all four sides clamped) from the infinite plate solution is only 2%. It can be concluded that when the aspect ratio 2 2 l a > in the wider inflated region, and L l 2 a b − > − at the end of the peninsula (i.e. the narrow inflated regions at the two sides of the peninsula), errors in the above derivation from infinite plate solution are negligible. Furthermore, a constant energy release rate 1 G is obtained when debond occurs within the limits of 4a < l < L − 2(a −b) . FEA simulations point to 71 similar conclusion as shown in Fig. 53. When aspect ratio 1.5 2 l a < , value of energy release rate 1 G increases as the debond grows. For larger blister lengths, energy release rate remains constant and converges to the solution of the infinitely long plate. Additional verifications using numerical techniques will be presented in section 5. 5.4 TimeDependent Behavior of Polymeric Thin Film Polymeric thin films typically display some timedependent behavior, such as, creep and/or relaxation due to molecular motion. A linear viscoelastic model was employed in the current study to simulate the timedependent behavior of the polymeric film. The threedimensional constitutive equations of linear viscoelastic material are, τ τ σ t K t τ θ d t kk ∫ ∂ = − ∂ 0 ( ) 3 ( ) (57) τ τ μ τ d e S t t t ij ij ∫ ∂ ∂ = − 0 ( ) 2 ( ) (58) where ij s and ij e are the stress and mechanical strain deviators, respectively, kk θ =ε is the volumetric strain. The bulk relaxation modulus K(t) and shear relaxation modulus μ (t) can be expressed using Prony Series as, Σ= − ∞ = + n i t i K t K K e i 1 ( ) τ (59) Σ= − ∞ = + n i t i t e i 1 μ ( ) μ μ τ (510) 72 The linear viscoelastic model for the thin film was implemented into the FEA simulation code NOVA3D. FEA simulations of debond in a peninsular blister specimen using elastic and viscoelastic thin film were conducted and compared with test data. 5.5 Simulations of Debonding in the Peninsula Blister Experiment 5.5.1 Simulation under Small Deformation A threedimensional cohesive layer model with 20node quadratic brick element was implemented into a testbed FEA code NOVA3D. Convergence studies showed that stable debonding was obtained within a fairly large debond region if the aspect ratio of the element in the debond tip area and the increment of liquid volume (loading increment) are sufficiently small. As the liquid is injected into the blister at a constant rate, the pressure also increases linearly with time (or liquid volume) before debond occurs. When the failure criterion given by Eq. (52) is satisfied, the cohesive layer elements at the debond tip fail and debonding initiates. Instead of the node release and element deletion schemes used in most finite element codes, a failed element remains active in the subsequent analysis while the stiffness of the element is reduced to approximately zero. As the cohesive layer elements fail, new surfaces are generated and the pressurized liquid occupies the newly debonded volume . In the present work, the original aspect ratio at the beginning of simulation is set as 0.75 2 l a = . During the simulation, as the debond initiates and propagates, the film deflection and aspect ratio increase while liquid pressure decreases, as shown in Fig. 54. 73 When the aspect ratio is greater than 1.5, liquid pressure stabilizes and remains constant, and a constant strain energy release rate is obtained. However, the maximum deflection in the thin film increases continuously until the aspect ratio is greater than 2. As alluded to earlier, this result agrees with Timoshenko’s solution for infinite plate [44]. Debond process near the end of the peninsula was also studied in this paper. When the aspect ratio in the inflated region L l 1.5 a b − < − , pressure will start to increase due to the boundary effects (These simulation results will be discussed in Section 6.3). From both analytical and FEA simulations, it can be concluded that the effective test domain for peninsula blister test is: 1.5 2 l a > and L l 1.5 a b − > − . When debond occurs within this domain, constant strain energy release rates will be obtained. Time (sec) Liquid pressure (KPa) Deflection (w/h), Aspect ratio 0 50 100 150 200 250 0 10 20 30 40 0 0.5 1 1.5 2 2.5 Liquid pressure Film deflection Aspect ratio Fig. 54 Peninsula blister responses over large debond lengths 74 5.5.2 Peninsula Blister Specimen and Test Results The data used here was obtained from a peninsula blister specimen with dimensions L = 35.56mm , l =15.24mm , 2a =10.16mm , 2b =1.27mm ( 2 0.125 2 b a = ) with an initial aspect ratio 1.5 2 l a = (as illustrated in Fig. 51). The epoxy thin film in the current simulation was assumed as linear elastic and its Young’s modulus was E =1.78GPa , Poisson’s ratio was υ = 0.3 , and film thickness was 0.127 fh = mm . Residual stresses R R 5.7 x y σ =σ = MPa were obtained from bulge tests and, for analysis, were applied to the thin film before loading. The aluminum substrate was considered to be rigid both. The critical energy release rates (or fracture energy) extracted from the blister test was 100 ~ 110 / 2 cG = J m (including residual stresses) and 130 ~ 140 / 2 cG = J m (without residual stress). In addition, a critical pressure of 80.9kPa was observed at debond initiation, along with a constant pressure of 86.5kPa during stable debonding. The adhesive material which was used to bond the epoxy thin film to the substrate is Hysol EA9696, with an adhesive strength of max σ = 42.7MPa . The maximum strain in the adhesive layer is assumed to be max ε = 0.1. When applying cohesive layer model to this specimen, a cohesive layer thickness 0.0416 ch = mm was determined from Eq. (2 16) (including residual stresses in the thin film and a fracture energy 100 / 2 cG = J m ). It should be clarified here that in the present case, the thickness of adhesive layer is different from the cohesive layer thickness. The latter is determined by adhesive 75 material properties and corresponding fracture energy or debond energy as presented in Eq. (216). Consequently, only a portion of the adhesive layer is modeled by the special threedimensional cohesive layer elements and the remainder of the adhesive layer was modeled by regular elements. Exact measurements of the adhesive layer thickness in peninsula blister specimen are difficult and the debond process is not sensitive to the adhesive layer thickness. Therefore we can reasonably assume that the adhesive layer thickness is equal to the cohesive layer thickness. In other words, only one layer of cohesive layer elements was employed to represent the entire adhesive layer. 5.5.3 Simulation with Large Deformation and Residual Stresses For peninsula blister specimen, both the thin film and adhesive layer experience large deflections and rotations near the debond tip during debond growth. An updated Lagrangian (UL) formulation [45] was employed to cope with the geometrical nonlinear characteristics in the thin film (nonlinear straindisplacement relations, while the stressstrain relation in the thin film is still linear). The UL formulation was combined with Cauchy stress and Almansi strain tensor as energy conjugates. A typical finite element mesh for the threedimensional peninsula specimen is shown in Fig. 55. Only one layer of threedimensional 20node brick element was used to model both the thin film and a portion of the adhesive layer. In the plane of the film, element dimension along the debond (peninsula) direction is greatly dependent on the film thickness as smaller film thickness will experience larger bending deformation near the debond tip; a fine mesh was used to account for the sharp stress gradients and large curvature near the debond tip. Convergence studies showed that the critical pressure is 76 rather sensitive to the element dimension in the peninsula direction (Table 51). However, when the ratio of element size in the xdirection to film thickness was Δx hf = 0.4 , convergence was achieved. Y Z X Node # = 7389 Ele # = 978 Cohesive Ele # = 186 Fig. 55 3D FEA mesh with 20node brick element (Viewed from the bottom, elements on the top are the cohesive layer element) Table 51 Different element size in peninsula direction Element size ( f Δx h *) 1.6 0.8 0.4 0.2 Critical pressure (kPa) 110.5 100.0 95.5 94.0 Deflection ( f w h ) 3.87 3.70 3.62 3.59 Liquid volume (mm3) 46.6 44.2 43.0 42.6 *: Δx is the element dimension in peninsula direction, f h is the film thickness, w is the maximum deflection in the thin film 77 Residual stresses, due to mismatch of thermal coefficients of expansion are usually present in the thin film following cure. Assuming a thermal expansion coefficient α = 0.0001/o C and a temperature change ΔT = −23oC , resulted in uniform residual stresses R R 5.7 x y σ =σ = MPa throughout the thin film. The effect of residual stresses on critical pressure and maximum deflection were considered (Table 52) for equibiaxial residual stresses ranging from 0 to 15 MPa. The impact on critical pressure and maximum deflection were noticeable. The mode I/II phase angle φ ranged from 38o ~ 53o while the mode I/III phase angle ϕ varied from 9o ~ 13o over the residual stress levels that were considered. They remained nearly constant during debond growth. This indicates that mixedmode I/II debond failure is dominant and the mode III contribution is relatively small, and may be neglected in peninsula blister specimen. Table 52 Effect of residual stresses on debond process as predicted by FEA simulations Residual stress R y R x σ =σ ( MPa ) 0 3.0 5.7 10.0 15.0 Critical pressure ( KPa ) 79.7 88.6 95.5 106.0 116.6 Maximum deflection ( f w h ) 4.01 3.81 3.62 3.35 3.08 Liquid volume (mm3 ) 48.0 45.2 43.0 40.2 36.8 Phase angle (I and II) ( o ) 53.5 49.5 46.1 41.8 37.8 Phase angle (I and III) ( o ) 12.6 11.9 11.3 10.4 9.7 Note: Phase angle is the average value for the elements along the perpendicular direction of the peninsula. 78 Measured and predicted pressure and maximum deflection history are given in Fig. 56. Stable debond pressure predicted by the elastic FEA simulation is P = 96.0KPa , which is approximately 11% higher than the test value of P = 86.5KPa . The error is probably due to the fact that the viscoelastic effects in the thin film were ignored. Another possible source of error is that the simple failure criterion given by Eq. (52) may not exactly reflect the real debond conditions. The initial aspect ratio ( 2 l a ) of current peninsula blister specimen is 1.5, and therefore the pressure would remain constant during debond growth process while maximum deflection increases continuously as the debond propagates until the aspect ratio reaches 2.5. This is in agreement with that derived from the theory of thin plates. Threedimensional film displacement profiles from FEA simulations at different stages of debonding are shown in Fig. 57. The debonding process near the end of the peninsula region was also studied. According to the analytical solution derived from the theory of elastic plates, pressure will start to increase when the aspect ratio in the narrow region L l 1.5 a b − < − . FEA simulation results showed (Fig. 58) that when the debond tip approached the end of the peninsula (an aspect ratio L l 1.0 a b − < − ), the pressure and film deflection started to increase with debond growth. 79 Time (x1000 sec) Pressure (KPa) Debond length (mm) 0 1 2 3 4 0 20 40 60 80 100 0 1 2 3 4 5 6 7 8 9 Liquid pressure Debond length Debond initiation (a) Test data Time (s) Liquid pressure (KPa) Deflection (w/h), Debond length (mm) 0 50 100 150 200 0 20 40 60 80 100 0 1 2 3 4 5 6 7 8 Liquid pressure Film deflection Debond length (b) Simulation results Fig. 56 FEA simulation of peninsula blister responses (with large deformation and residual stresses) 80 (a) Aspect ratio = 1 (b) Aspect ratio = 2 (c) Aspect ratio = 2.5 Fig.57 Simulated 3D film displacement profiles at different stages of debonding 81 Aspect ratio (Ll)/(ab) Liquid pressure (KPa) Filmdeflection (w/h) 1.5 1 0.5 0 40 60 80 100 120 140 160 180 200 3.5 4 4.5 5 5.5 6 Liquid pressure Film deflection Fig. 58 Predicted liquid pressure and film deflection increase when debond approaches the end of the peninsula 5.5.4 Simulation Including TimeDependent Effect In the simulation of geometrically nonlinear, elastic debonding in thin film peninsula blister specimens with an initial aspect ratio of 1.5, critical pressure remains constant once the debond has initiated. However, in the actual peninsula blister test, a moderate increase in pressure was observed in the early stages of debond growth (Fig. 5 6). This suggested that there might be some timedependent effects in the epoxy thin film. Due to unavailability of the timedependent material properties of the epoxy thin film in the current blister test, a demonstration simulation considering linear viscoelastic effects in the thin film was performed to examine the influence of timedependent material properties on the critical pressure and the debond process. 82 For the purpose of highlighting the effect of linear viscoelasticity of the thin film, the material properties for a linear viscoelastic epoxy thin film corresponding to the elastic film were defined and scaled from the tested values as follows: the tensile modulus E0 =1.78GPa and Poisson’s ratio 0 υ = 0.3 of the elastic thin film correspond to the values of the viscoelastic tensile relaxation modulus, E(t) , and Poisson’s ratio, υ (t) , at the time t = 0 , i.e., 0 0 ( ) 1.78 t E Et GPa = = = 0 0 ( ) 0.3 t υ υ t = = = For a typical linear viscoelastic epoxy thin film, material properties such as bulk relaxation modulus K(t) , shear relaxation modulus μ (t) , tensile modulus E(t) , and Poisson’s ratio υ (t) at room temperature can be expressed in the form of a Prony Series. The Prony Series coefficients for the epoxy thin film are listed in Table 53 using data obtained for an epoxy [46]. For a liquid injection rate R = 80mm3 / hour , the simulation accounting for large deformations and residual stresses showed that the liquid pressure remained constant after debond initiation for a linearly elastic thin film. However, for a viscoelastic thin film (Fig. 59a), the liquid pressure did increase at the early stages of debond growth before reaching a steady state value. Stable debond processes with constant liquid pressures were obtained in both cases with constant but slightly different energy release rates. 83 Table 53 Material properties of a viscoelastic epoxy ( ) i Log τ i K (MPa) i μ (MPa) i E (MPa) i υ 1 4.3 24.64 0.22 0.21 3.11 x 104 2 3.3 25.54 0.45 0.11 4.95 x 104 3 2.3 104.24 20.88 68.48 8.31 x 104 4 1.3 29.14 20.32 29.72 1.13 x 103 5 0.3 150.32 22.11 93.33 3.19 x 103 6 0.7 155.61 160.32 387.51 2.40 x 103 7 1.7 440.03 161.77 491.45 5.44 x 103 8 2.7 330.48 310.46 758.45 6.53 x 103 9 3.7 438.32 159.01 552.47 1.19 x 102 10 4.7 226.82 187.85 377.34 1.66 x 102 11 5.7 206.62 71.09 256.68 1.74 x 102 12 6.7 245.83 63.80 137.01 2.59 x 102 13 7.7 49.88 21.35 95.43 6.03 x 102 14 8.7 169.51 20.31 47.03 1.18 x 102 15 9.7 34.86 11.30 53.36 3.00 x 103 16 10.7 82.59 6.45 15.69 3.04 x 103 17 11.7 11.97 3.08 17.01 2.20 x 103 18 12.7 40.06 0.25 4.14 9.92 x 104 19 13.7 8.54 1.29 7.19 1.57 x 104 ∞ 2337.41 12.34 36.91 0.4974 *: At room temperature 23.5°C and relative humidity 35% 84 Time (x1000 sec) Liquid pressure (KPa) 0 1 2 3 4 0 20 40 60 80 100 Test data Elastic film Viscoelastic film (a) Liquid pressure history Time (x1000 sec) Debond length (mm) 0 1 2 3 4 0 1 2 3 4 5 6 7 Test data Elastic film Viscoelastic film (b) Debond length history Fig. 59 Peninsula blister test data and FEA simulation results (with large deformation and residual stresses) 85 Measurements and simulation results for both elastic and viscoelastic thin films are shown in Fig. 59 and tabulated in Table 54 for comparison. Fig. 59(a) compares the pressure history as a function of time, while Fig. 59(b) depicts comparison of the debond growth history. From Fig. 59 and Table 54 it is evident that including viscoelastic timedependence in the thin film response results in a significant reduction in the error between simulations and test data. It should be noted that the slight discrepancy in the time scales of the simulation and test data is due to some leakage of liquid that was observed during the peninsula blister test especially at elevated pressure levels. When a lower liquid injection rate R =10mm3 / hour was employed for the viscoelastic thin film simulation, the debond initiation pressure decreased from 88.5KPa (corresponding to R = 80mm3 / hour ) to 84.0KPa while the peak pressure and stable debond pressure remained unchanged. Thus it is likely that the debond initiation pressure will be closer to the test value of 80.9kPa as the injection rate is decreased. Table 54 Peninsula blister test data and FEA predictions with elastic and viscoelastic thin film Test data Linear elastic thin film Viscoelastic thin film Event Time Pressure Time Pressure Error Time Pressure Error Debond initiation 2250 80.9 1880 94.5 16.8 2210 88.5 9.39 Peak pressure 2980 88.9 1910 96.0 7.99 2340 95.5 7.42 Stable debond 3400 86.5 1910 96.0 11.0 2490 91.5 5.78 Unit: Time  second, Pressure  kPa , Error  % (predicted pressure compared with test data) 86 5.6 Conclusions Peninsula blister specimen is an effective way to measure the interfacial fracture toughness of a variety of adhesive bond systems as it offers very high energy release rate values and maintains constant energy release rates over a large range of debond lengths. A threedimensional cohesive layer model and a corresponding threedimensional mixedmode failure (debond) criterion were developed based on the principles of continuum mechanics and fracture mechanics. It was implemented in a finite element code to simulate quasi static debonding in the peninsula blister test. Numerical convergence and stable debond growth were obtained over a fairly large range of debond lengths. The results from FEA simulations were in reasonable agreement with both an analytical solution and test data. Suitable geometries for the peninsula blister specimen were studied by both analytical and FEA approaches and guidelines were reiterated for the design of peninsula blister specimen. FEA simulation results also showed that large deflections, timedependent material behavior, and residual stresses in the thin film are important factors that should be considered in simulations of the peninsula blister test in order to extract the interfacial fracture toughness of a given adhesive system. 87 CHAPTER VI SIMULATION OF TIMEDEPENDENT DEBOND GROWTH The objective of this chapter is to model the synergistic bond degradation mechanisms that may occur at the interface between a fiber reinforced polymer (FRP) and a substrate. FEA Simulation of a wedge test was conducted, and the timedependent effect in the adhesive layer was involved in the simulation model. The results predicted by the computational model were benchmarked through comparison with analytical solutions and mixed mode fracture tests. Steady debond growth was obtained after the wedge front entered the originally bonded area, which is consistent with the observations from wedge tests. 6.1 Introduction Polymeric thin films usually display some timedependent behaviors due to molecular motion. This effect was also observed in blister test and creep test on a fully clamped polymeric plate. A linear viscoelastic model was introduced to simulate the timedependent behavior of the polymeric film by peninsula blister test (Chapter V). NonFickian hygrothermal effects on cohesive layer are derived in Chapter IV. It should be noted that in the present approach, expansion of Helmholtz free energy in terms of convolution integrals was not carried out to directly include viscoelasticity in the 88 cohesive stressstrain law. This is because the use of convolution integrals in addition to temperature and moisture dependence would render the cohesive stressstrain law intractable as far as characterization of the convolution coefficients is concerned. Therefore, in the interest of tractability, a simplified approach is employed where the ratedependent behavior in the cohesive layer is implemented through the characterization of ratedependence of the maximum stresses and maximum strains in the cohesive layer as presented in Table 61. The remainder of the polymeric adhesive outside the cohesive layer is modeled as a nonlinear viscoelastic continuum with timedependent constitutive behavior. The influence of temperature and moisture concentration on the workofseparation and on crack growth is derived from firstprinciples. The model is implemented in a testbed finite element code NOVA3D. Results predicted by the cohesive layer model are benchmarked through comparison with experimental data from mixedmode fracture experiments performed using a moving wedge test. Ratedependent debond process was also investigated with this model under different debond speeds. 6.2 Failure Criterion Based on WorkofSeparation Several mixed mode failure criteria of the cohesive layer are described in Chapter II. A new failure criterion based on fracture energy was introduced to predict the debond process of a wedge test. Fracture energy 2Γ of the adhesive was extracted from the wedge test. Let the fracture energy 2Γ equal to the workofseparation of cohesive layer with mixed mode I and mode II fracture, gives 2 ( II ) sep I sep sep Γ =ηφ =η αφ +βφ (61) 89 Where I sep φ and II sep φ are the workofseparation of pure mode I and mode II debond, respectively (Chapter II and Chapter IV) For pure mode I fracture under dry conditions (Δm = 0 ): c I sep h 16 max max φ = 9 σ ε For pure mode II fracture under dry conditions (Δm = 0 ): c II sep h 16 max max φ = 9 τ γ sep φ is the total workofseparation of the mixed mode debond. It is easy to derive the relation II sep I sep sep φ =αφ +βφ , where α and β are constants determined by the failure type (i.e. mode mix) of the cohesive layer. Comparing the idealized cohesive stressstrain curve in Fig. 23 with the uniaxial tension test data of the adhesive material in Fig. 61, it was found that the actual workofseparation from the test data is less than the theoretic one predicted by tractionseparation law due to the premature failure of epoxy primer. Thus a correction factor η = 0.75 was introduced to reflect the difference. Fig. 61 Stressstrain relation of epoxy under different strain rates (test data) 90 Therefore cohesive layer thickness c h can be obtained from Eq. (61) if the fracture energy 2Γ or workofseparation sep φ , which is the energy necessary to generate unit debond length, is known from experiment or analytical methods. 6.3 Nonlinear Viscoelastic Model and Fracture Energy A linear viscoelastic model is described in Chapter V, material constitutive equations are listed in Eqs. (57, 58) and the timedependent material properties can be expressed with Prony series as shown in Eqs. (59, 510). A nonlinear viscoelastic model for the bulk adhesive using modified free volume approach [47] was also implemented in the code. A strainbased formulation proposed by Popelar and Liechti [48] entails a nonlinear shiftfactor based on freevolume given by, ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + − ⎟ ⎟⎠ ⎞ ⎜ ⎜⎝ ⎛ + Δ + Δ + Δ + Δ + = − s eff eff s s d v v v v d d f f B f T c T c f B Log a ε ε α β δθ α β δθ 2.303 2.303 ( ) (62) The generalized Jintegral for large deformation with pseudostress model proposed by Schapery [49] is also used for fracture modeling, ∫Γ ⎥⎦ ⎤ ⎢⎣ ⎡ ∂ ∂ = Φ − dL x u J dy T R R i v i (63) For this case, the critical work input required to initiate the crack (or, the fracture initiation energy) can be expressed in terms of the timedependent farfield parameter v J . ∫ ∂ ∂ Γ = − t v i R d J E D t 0 2 ( ) τ τ τ (64) 91 6.4 Moving Wedge Test To study the properties of polymer adhesive and its debond process under different environmental conditions and strainrates, 



A 

B 

C 

D 

E 

F 

I 

J 

K 

L 

O 

P 

R 

S 

T 

U 

V 

W 


