UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Crack growth and damage modeling of fibre reinforced polymer composites McClennan, Scott Andrew 2004

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2004-0552.pdf [ 13.7MB ]
Metadata
JSON: 831-1.0078773.json
JSON-LD: 831-1.0078773-ld.json
RDF/XML (Pretty): 831-1.0078773-rdf.xml
RDF/JSON: 831-1.0078773-rdf.json
Turtle: 831-1.0078773-turtle.txt
N-Triples: 831-1.0078773-rdf-ntriples.txt
Original Record: 831-1.0078773-source.json
Full Text
831-1.0078773-fulltext.txt
Citation
831-1.0078773.ris

Full Text

CRACK GROWTH AND DAMAGE MODELING OF FIBRE REINFORCED POLYMER COMPOSITES by Scott A n d r e w McClennan B.S. (Aeronaut ical and Astronaut ical Engineer ing) , Un ivers i ty o f I l l ino is , 2001 A THESIS S U B M I T T E D I N P A R T I A L F U L F I L M E N T OF T H E R E Q U I R E M E N T S F O R T H E D E G R E E OF M A S T E R OF A P P L I E D S C I E N C E in T H E F A C U L T Y OF G R A D U A T E S T U D I E S (Department o f Metals and Mater ials Engineer ing) W e accept this thesis as con fo rming to the required standard T H E U N I V E R S I T Y OF B R I T I S H C O L U M B I A October 2004 © Scott A n d r e w McClennan , 2004 THE UNIVERSITY OF BRITISH COLUMBIA FACULTY OF GRADUATE STUDIES Library Authorization 3 In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Name of Author (please print) Date (dd/mm/yyyy) Title of Thesis: Degree: t A c v ^ W r oV / \ ^ V ^ 1 ^ C T ^ ^ C Q - Y E A R : 2.00^ Department of fAoi\<e.n o\<> ~€-The University of British Columbia i i i i l iVancouver, BC Canada grad.ubc.ca/forms/?formlD=THS V ', page 1 of/I . . . . «;.• . ,. . ' . last updated:Q0-Jul-04-Abstract A B S T R A C T Notched tensile strength and the development o f damage in composite laminates are studied in this thesis through exper imental , numer ical and analyt ical methods to develop simple models for predic t ing notched strength using complete ly physical input parameters. A series o f experimental tensile fracture tests using the Overheight Compact Tension ( O C T ) specimen geometry, established by (Kongshavn and Poursart ip, 1999), were conducted to study the development o f the characteristic damage zone. The material used in the tests was a quasi- isotropic carbon f ibre / epoxy. The specimens were mod i f ied so that the notch root radius var ied wh i le mainta in ing a constant crack w i d t h to specimen section rat io. Specimens w i t h notch root radi i less than 16 m m display stable crack g rowth w i t h l i t t le notch sensit iv i ty. Specimens w i t h larger notch root radi i are unstable and display more notch sensit iv i ty. For al l specimens, the height o f the damage zone converges to between 6 m m and 7 m m . The transit ion f r o m stable to unstable behaviour is explained using fracture mechanics equations and a transit ion radius can be determined f r o m three material parameters; elastic modulus, specif ic strain energy and tensile strength. A simple bi l inear cohesive zone model cal led the simple damage model ( S D M ) was developed as a material model i n the A B A Q U S f in i te element code to be used as a demonstrator for mode l ing techniques that can be appl ied to al l strain-softening material models. I n part icular, the relat ionship between input parameters and the element w i d t h was studied. Us ing the same input parameters as the analyt ical model and a suf f ic ient ly ref ined mesh, the numer ical model predicts the peak loads from the series o f experimental O C T tests w e l l . A t ransi t ion f r o m stable to unstable behaviour is predicted at the same radius as seen experimental ly. A method o f determining an appropriate element w i d t h to ensure an accurate predict ion is presented. A mod i f i ed version o f the S D M cal led the adaptive simple damage model ( A S D M ) was also developed, w h i c h automat ical ly scales the input strength to account for the effect o f element w id th . This mod i f i ca t ion al lows larger elements to be used wh i le st i l l obtain ing an accurate solut ion, a useful feature i f large structures are being modeled. i i Table of Contents T A B L E OF CONTENTS Abstract ii List of Tables vi List of figures vii Nomenclature xii Acknowledgements xiv Chapter 1 Introduction 1 Chapter 2 Background and literature review 3 2.1 Int roduct ion 3 2.2 Fracture mechanics 3 2.3 Cohesive crack models 7 2.3.1 Cohesive crack models for intralaminar crack propagation 8 2.3.2 Cohesive zone models for interlaminar crack propagation 11 2.3.3 Mesh Sensitivity 13 2.4 Notched strength predict ion 15 2.4.1 Semi-empirical models 15 2.4.2 Cohesive models 17 2.4.3 Relation of the damage zone size to strength 20 2.5 A specimen for stable damage g rowth 21 2.6 Summary 22 Chapter 3 Experiments 28 3.1 In t roduct ion 28 3.2 Specimens and procedures 28 3.2.1 Material 28 3.2.2 Test equipment 29 3.2.3 OCT tests 29 3.2.4 Four-point bend tests 32 3.3 Results and data reduct ion 33 3.3.1 OCT tests 33 Table of Contents 3.3.2 Four-point bend tests 40 3.4 Summary 41 Chapter 4 Numerical 61 4.1 In t roduct ion 61 4.2 Elastic simulat ions 62 4.2.1 Stress concentrations 62 4.2.2 Stress contours 63 4.2.3 Element size required to resolve the notch tip stress concentration 63 4.3 Development o f the cohesive zone models 64 4.3.1 Simple damage model (SDM) 64 4.3.2 Adaptive simple damage model (ASDM) 65 4.3.3 Algorithm for the SDM and ASDM. 68 4.3.4 How mesh insensitive is the ASDM? 72 4.4 O C T simulat ions 73 4.4.1 The virtual specimens 73 4.4.2 Parametric study on Gc and apeak 75 4.4.3 SDM/ ASDM comparison 75 4.4.4 Simulation of experimental OCT tests 77 4.5 Summary 83 Chapter 5 Discussion and analysis of results 104 5.1 In t roduct ion 104 5.2 Ana ly t i ca l descript ion o f notch root radius sensit iv i ty 104 5.2.1 Regime 1 104 5.2.2 Regime II. 106 5.2.3 Trans ition radius 108 5.2.4 Peak load solution for the OCT specimens 108 5.3 Damage zone size and shape I l l 5.4 Signif icance o f the transi t ion radius 112 5.5 Other material systems and geometries 113 5.6 Choosing an appropriate element w i d t h for strain-softening models 113 5.6.1 Element width case study 114 5.6.2 Guidelines for choosing an appropriate element size 775 i v Table of Contents 5.7 Summary 118 Chapter 6 Conclusions and future work 128 6.1 Conclusions 128 6.2 Future w o r k 130 Appendix A Load versus CMOD plots 132 Appendix B Extra line analysis plots 141 Appendix C Stress distribution and compliance plots 145 Appendix D ASDM UMAT for ABAQUS 148 Appendix E Transition radius and peak load prediction plot for an AS4/3501-6 DCB Specimen 160 E. 1 Specimen 160 E.2 Regime I 160 E.3 Regime I I 161 E.4 Trans i t ion radius 161 E.5 Discussion 161 References 165 v List of Tables LIST OF TABLES Table 3.1 Or ig ina l elastic mater ial properties for CFRP laminate (M i t che l l , 2002) 43 Table 3.2 O C T tests and specimen measurements 44 Table 3.3 Four-point bend tests and specimen measurements 45 Table 3.4 Summary o f O C T tests and analysis 46 Table 3.5 Cr i t ica l fracture energy release rates 47 Table 3.6 Four-point bend test results 48 Table 4.1 Mater ia l properties fo r numer ical O C T tests 85 Table 4.2 Height o f cr i t ica l ly stressed material at peak load 85 Table 4.3 Numer ica l tests 86 Table 4.4 A S D M val idat ion tests 87 Table 5.1 Approx imate material properties for acry l ic 120 Table 5.2 Mesh ing guidel ines based on notch radius 120 Table E. l AS4/3501-6 material properties 163 List of Figures LIST OF FIGURES Figure 2.1 Coordinates at a crack t ip 24 Figure 2.2 A slender notch w i t h a b lunt notch t ip 24 Figure 2.3 Cohesive crack model 24 Figure 2.4 Tract ion - displacement relat ionship for cohesive crack model 25 Figure 2.5 Cohesive crack models for a) duct i le solids and b) quasi-bri t t le solids ( f rom (de Borst , 2003)) 25 Figure 2.6 Damage local isat ion using a cohesive zone model 26 Figure 2.7 Damaged notch specimen o f Kor tschot and Beaumont (Kor tschot and Beaumont, 1990a) 27 Figure 3.1 O C T specimen geometry (thickness B ranged f r o m 7.9 m m to 8.6 m m ) 49 Figure 3.2 O C T specimens w i t h vary ing notch root radi i a) 0.5 m m , b) 2.0 m m , c) 5.0 m m , d) 12.7 m m , e) 19.0 m m and f ) 27.0 m m : 50 Figure 3.3 Scribed lines for l ine analysis (shown for rd4-2) 51 Figure 3.4 O C T test set-up (shown for rd54-2) 51 Figure 3.5 Load-displacement plots for specimen rd32-2 from the machine posi t ion and f r o m the relat ive displacement o f the pins. 52 Figure 3.6 Schematic o f section taken through the damage zone ahead o f the in i t ia l notch t ip. The dashed lines represent the 90° fibres in the outer stacks. The measured damage height d id not include any damage to the outside o f these lines 52 Figure 3.7 Four-point bend test geometry 53 Figure 3.8 Specimen r d l 0 - 2 showing or ientat ion o f crack and extent o f surface damage (scribed lines are 2.5 m m apart) 53 Figure 3.9 Load versus crack mouth opening displacement for specimens o f representative notch root radius. Sol id lines indicate continuous port ions o f the curve wh i le dashed segments indicate segments in w h i c h the crack j umps and the load drops sharply 54 Figure 3.10 Peak loads attained in each O C T test versus notch root radius p. Closed data points denote specimens d isplay ing a plateau on the load versus C M O D curve, open data points denote specimens that showed sharp load drops at the peak load 54 Figure 3.11 Prof i les o f the damage zone along the crack plane 55 List of Figures Figure 3.12 Section o f r d 2 - l a t 0.635 m m ahead o f in i t ia l notch t ip.. . , 55 Figure 3.13 Section o f rd54-2 at 0.635 m m ahead o f in i t ia l notch tip.r. 56 Figure 3.14 Section o f r d 2 - l at 24.13 m m ahead o f in i t ia l notch t ip 56 Figure 3.15 Section o f rd54-2 at 24.13 m m ahead o f in i t ia l notch t ip 57 Figure 3.16 C O D prof i les determined from l ine analysis for rd25-2 57 Figure 3.17 Load versus p in opening displacement for rd25-2 showing locat ion o f photos for l ine analysis 58 Figure 3.18 C O D prof i les determined f r o m l ine analysis for rd28-2 58 Figure 3.19 Load versus p in opening displacement for rd28-2 showing locat ion o f photos for l ine analysis 59 Figure 3.20 Final crack lengths in r d 2 - l and rd54-2 f r o m sect ioning and l ine analysis 59 Figure 3.21 Crack resistance points for rd25-2 and rd28-2 60 Figure 3.22 Four-point bend test load-displacement p lot for b85-90-2 60 Figure 4.1 Elastic stress fields ahead o f notch t ip at P = 15.0 k N 88 Figure 4.2 Elastic stress field for p = 27 m m at exper imental ly determined peak load o f 17.8 k N 88 Figure 4.3 Stress contours ahead o f notch t ip (symmetr ic about the notch plane) for a) p = 1 m m at 15 k N load, b) p = 4 m m at 15 k N load, c) p = 10.5 m m at 15 k N load, d) p = 16 m m at 16.9 k N load and e) p = 27 m m at 17.8 k N load. Elements at the notch t ip are 0.5 m m x 0.5 m m in al l cases. The l ight grey regions indicate a stress o f 460 M P a 89 Figure 4.4 Stress concentrations on the notch plane ahead o f a 4 m m radius notch t ip in an O C T specimen as resolved by di f ferent sized constant stress (cs) and fu l l y integrated quadratic (quad) elements. Stresses are extrapolated to the nodes 90 Figure 4.5 Effect o f element w i d t h on stress calculated at a notch t ip 90 Figure 4.6 Load versus arm displacement for D C B simulat ions showing the effect o f element w i d t h on the peak load predict ion 91 Figure 4.7 Peak loads for D C B simulat ion using element meshes o f vary ing coarseness w i t h and w i thou t scaling the material strength to match the 62.5 p m mesh size 91 Figure 4.8 Strength scaling ratios to match the peak load in the D C B simulat ion using a mesh w i t h 62.5 p m constant stress elements at the crack t ip 92 Figure 4.9 Simple bi l inear damage model 92 Figure 4.10 F l o w chart for S D M / A S D M algor i thm 93 Figure 4.11 Integrat ion point number ing 94 v i i i List of Figures Figure 4.12 Scaled stress-strain curve 94 Figure 4.13 Damage zone size determined by the cohesive zone model 95 Figure 4.14 Example o f energy release determined using elements o f di f ferent size 95 Figure 4.15 V i r tua l O C T specimen mesh for a) fu l l specimen and b) at notch plane 96 Figure 4.16 Parametric study o f G c and a p e a k for rd2 96 Figure 4.17 Reduct ion in e* p e a k at the integrat ion points o f the notch t ip element in rd8 97 Figure 4.18 Peak loads predicted using the A S D M (large diamonds) p lot ted over the experimental results (smal l squares). Open symbols represent cases that were unstable (i.e. numer ica l run crashed or d id not display load plateau exper imental ly) . No te that the ordinte scale begins at 10 k N to magn i fy the results , 97 Figure 4.19 N o r m a l stress contours for rd2 a) at C M O D = 0 . 6 7 9 m m (before peak load), b) at C M O D = 1.70 m m (peak load) and c) at C M O D = 3 . 2 1 m m (after peak load) 98 Figure 4.20 N o r m a l stress contours for rd54 a) at C M O D = 0 . 6 7 4 m m , b) at C M O D = 1 . 7 0 m m and c) at C M O D = 2 . 4 8 m m (peak load) 99 Figure 4.21 N o r m a l stress contours for rd32 a) at C M O D = 0 . 6 7 4 m m (before peak load), b) at C M O D = 2 . 1 1 m m (peak load) and c) at C M O D = 3 . 2 1 m m (after peak load) 100 Figure 4.22 Progression o f the stress d is t r ibut ion along the notch plane for rd2. Each dist r ibut ion is at the C M O D value ( in m m ) noted above i t 101 Figure 4.23 Progression o f the stress d is t r ibut ion along the notch plane for rd28. Each dist r ibut ion is at the C M O D value ( in m m ) noted above i t 101 Figure 4.24 Damage zone w i d t h determinat ion for rd28 (Note that the crack inc lud ing the damage zone is 32.25 m m w ide at this point and the C M O D is 3.9 m m ) . . . 102 Figure 4.25 Compl iance for rd28 before shi f t ing A S D M predic t ion 102 Figure 4.26 Compl iance for rd28 after shi f t ing A S D M predic t ion by 3.25 m m to account for par t ia l ly damaged zone 103 Figure 4.27 Progression o f the crack front i n selected v i r tual O C T specimens 103 Figure 5.1 Stress and energy fai lure cr i ter ia for fracture in regime 1 121 Figure 5.2 Stress and energy fai lure cr i ter ia for fracture in regime I I 121 Figure 5.3 r"1/2 s ingular i ty at the t ip o f a sharp notch in an O C T specimen 122 Figure 5.4 Mesh o f quarter-point elements around a crack t ip ( f r o m (Cook, Ma lkus andPlesha, 1989)) 122 Figure 5.5 Ana ly t i ca l predict ion o f peak loads for O C T specimens 123 Figure 5.6 Cracks in i t ia t ing o f f the centerl ine in a) r d 4 5 - l and b) rd44-2 123 ix List of Figures Figure 5.7 A n ell ipse w i t h a major axis o f 8 m m and a m ino r axis o f 2 m m has the same radius o f curvature at the notch t ip as a circle w i t h a radius o f 16 m m 124 Figure 5.8 Stress d is t r ibut ion ahead o f a notch t ip showing an element o f w i d t h w e that is too large. Mesh must either be ref ined or strength scal ing used to achieve accurate results 124 Figure 5.9 Load versus C M O D results for a sharp-notched O C T specimen modeled w i t h di f ferent size elements 125 Figure 5.10 No tch t ip stress f ields for coarse and f ine meshed O C T specimens w i t h a sharp notch t ip at load equal to 8.3 k N 125 Figure 5.11 Element size contours for G c = 80 k J / m2 126 Figure 5.12 Large constant stress elements used to model the experimental 1 m m notch root radius O C T specimen al l g ive the same result 126 Figure 5.13 Element size contours for G c = 20 k J / m2 127 Figure 5.14 When G c is reduced to 20 kJ /m a dif ference is noted in the peak loads predicted using di f ferent w i d t h elements 127 Figure A . 1 Exper imental l o a d - C M O D results and numerical predic t ion for specimen rd2 133 Figure A . 2 Exper imental l o a d - C M O D results and numerical predic t ion for specimen rd8 133 Figure A .3 Exper imental l o a d - C M O D results and numerical predic t ion for specimen rd21 .' 134 Figure A . 4 Exper imental l o a d - C M O D results and numer ica l predic t ion for specimen rd28 134 Figure A . 5 Exper imental l o a d - C M O D results and numer ica l predic t ion for specimen rd32 135 Figure A . 6 Exper imental l o a d - C M O D results and numerical predic t ion for specimen rd38 135 Figure A . 7 Exper imental l o a d - C M O D results and numerical predic t ion for specimen rd54 136 Figure A .8 Exper imental l o a d - C M O D results for specimen r d l 136 Figure A .9 Exper imental l o a d - C M O D results for specimen rd4 137 Figure A . 10 Exper imental l o a d - C M O D results for specimen rd6 137 Figure A . l 1 Exper imental l o a d - C M O D results for specimen r d l O 138 Figure A . 12 Exper imental l o a d - C M O D results for specimen r d l 7 138 Figure A . 13 Exper imental l o a d - C M O D results for specimen rd25 139 x List of Figures Figure A . 14 Exper imental l o a d - C M O D results for specimen rd35 139 Figure A . 15 Exper imental l o a d - C M O D results for specimen rd44 140 Figure B . l C O D prof i les determined from l ine analysis for r d 2 - l 142 Figure B.2 C O D prof i les determined from l ine analysis for rd2-2 142 Figure B.3 C O D prof i les determined f r o m l ine analysis for rd8-2 143 Figure B.4 C O D prof i les determined from l ine analysis for r d 5 4 - l 143 Figure B.5 C O D prof i les determined from l ine analysis for rd54-2 144 Figure C. 1 Progression o f the stress d is t r ibut ion along the notch plane for rd8. Each dist r ibut ion is at the C M O D value ( in m m ) noted above i t 146 Figure C.2 Progression o f the stress d is t r ibut ion along the notch plane for r d 2 1 . Each dist r ibut ion is at the C M O D value ( in m m ) noted above it 146 Figure C.3 Progression o f the stress dist r ibut ion along the notch plane for rd32. Each dist r ibut ion is at the C M O D value ( in m m ) noted above i t 147 Figure C.4 Compl iance plots from locat ion o f peak stress in numer ical simulat ions. 3.25 m m is subtracted from Aa to account for w i d t h o f damage zone 147 Figure E. l D C B specimen (unidi rect ional f ibres oriented along length o f specimen) 164 Figure E.2 Peak load predict ion for AS4/3501-6 D C B specimen 164 Nomenclature N O M E N C L A T U R E S Strain Parametric element coordinate a Stress e Ang le in radians from notch plane p Notch root radius s Displacement V Parametric element coordinate S*peak Scaled peak strain S*ult Scaled ul t imate strain (7c Unnotched laminate tensile strength (7max M a x i m u m stress ( in stress distr ibut ion) SOpeak In i t ia l characteristic peak strain Yc Cri t ica l specif ic strain energy Yd Dissipated specif ic strain energy Sdam Damage zone strain Speak Current element peak strain Smax M a x i m u m strain in an element ( for scaling) <?N Notched laminate tensile strength Speak Current peak strain CTpeak Peak stress Yr Residual specif ic strain energy Sscale M a x i m u m value o f strain field used for scal ing Suit Ul t imate strain a Crack length ao Mater ia l system characteristic d imension (average stress method) A S D M Adapt ive Simple Damage M o d e l B Specimen thickness C Compl iance C M O D Crack M o u t h Opening Displacement D Depth ( four-po in t bend specimen) do Mater ia l system characteristic d imension (point stress method) D C B Double-cant i lever beam (specimen) E' Effect ive anisotropic modulus E0 In i t ia l modulus Er Residual secant modulus F W o r k per formed by external forces G Strain energy release rate Gc Cri t ica l fracture/ strain energy release rate GF Fracture energy hc Characteristic stable height o f damage zone Nomenclature hd Height o f damage zone he i Element height hs Height o f cr i t ica l ly stressed material Ki Stress intensity factor (mode I ) Kf Stress concentrat ion factor LL Load span ( four-point bend specimen) U Support span ( four-point bend specimen) O C T Over-height Compact Tension (specimen) P Load Pc Peak load (structural strength) P O D Pin Opening Displacement R Crack resistance energy r Distance f r o m notch t ip S D M Simple Damage M o d e l T Interface tract ion t Simulat ion t ime U Elastic strain energy W Section w i d t h o f O C T specimen Went Cr i t i ca l ly stressed w i d t h o f material ahead o f notch t ip VOdam Damage zone w i d t h Wtotal Total w o r k done X Coordinate in direct ion o f crack y Coordinate perpendicular to crack Acknowledgements A C K N O W L E D G E M E N T S First, m y deepest gratitude goes to m y advisor, Dr . Anoush Poursart ip, for g i v ing me the oppor tuni ty to w o r k in the Composites Group at U B C and for the countless hours he has spent keeping me focused in the r ight direct ion. He has the great ab i l i ty to get to the root o f every p rob lem that I b r ing h i m . I am also indebted to m y co-advisor, Dr . Reza V a z i r i , for his assistance and long discussions on the numerical side o f m y thesis. The w o r k that I have done here w o u l d not be possible w i thou t the w o r k o f numerous Composites Group members that have come before me, par t icu lar ly the w o r k o f J. Scott Ferguson, K e v i n Wi l l i ams , Ingr id Kongshavn, Jason M i t c h e l l , K a r i m K a n j i and A n t h o n y F loyd . I w o u l d l ike to thank al l o f the group members that have made m y t ime at U B C as enjoyable as i t was educational. Roger Bennett deserves special recogni t ion for prepar ing the specimens, m o d i f y i n g the testing j i gs and assisting me w i t h any equipment concerns that I had. I am also indebted to A n t h o n y F l o y d for answering every numer ica l question I cou ld come up w i t h and to K a r i m K a n j i for considerable assistance w i t h the experimental tests. A b ig thank y o u goes to m y entire fami l y , especially both o f m y parents, for support ing me throughout m y entire education no matter what I decided to do or what part o f the cont inent I decided to do i t on. Last, thanks Tracy, for being inspirat ional and for being w i t h me throughout, even when you ' re thousands o f miles away. x i v Chapter I Introduction C H A P T E R 1 INTRODUCTION Compared to most structural materials, fibre reinforced po lymer (FRP) composites are very new. They have on ly been used and studied for a l i t t le over h a l f a century, but in that t ime an enormous amount o f l i terature has been produced. Numerous models (empir ica l , analyt ical , and since the rise o f the f in i te element method, numer ical ) have been devised to expla in every aspect o f their behaviour, f r o m processing to damage and fracture. W h i l e the science behind composite materials has been advanced immensely by this research and the development o f these material models, i t presents a daunt ing task for an engineer or designer w h o wishes to use FRP composites in a design. H o w can the engineer narrow d o w n the number o f mode l ing choices and be certain that the model chosen for the design is the best one? In a recent exercise, several leading fai lure theories for composite laminates were tested against fourteen chal lenging test cases (H in ton , Kaddour and Soden, 2002b; H i n t o n and Soden, 1998). The exercise was a massive undertaking and revealed that the di f ferent models had di f ferent benefits, that even the best models could not describe al l test cases, and that the most compl icated models were not always the best models (H in ton , Kaddour and Soden, 2002a). The uncertainty in design using composite materials has been apparent in the aerospace industry. W h i l e the benefits o f using l ightweight , s t i f f and durable FRP composites that can be ta i lored to meet the requirements o f each appl icat ion are obvious, the uncertainty in the behaviour o f new FRPs has kept them f r o m being used more extensively. A l t hough FRPs have been used in re lat ively small quant i ty on most commerc ia l aircraft in the latter ha l f o f the last century, the commerc ia l aerospace industry has been cautious about using them for major load bearing structures unt i l very recently, due to a lack o f understanding o f their behaviour. A h igh number o f f l igh t hours w i t h relat ively l i t t le t ime for inspection has made the uncertainty in composite behaviour too large a factor to overcome. For mi l i ta ry aircraft, in w h i c h there is a much higher ratio o f maintenance to f l igh t hours, composites have always been used more extensively. W i t h the composites knowledge amassed over the past few decades however, the next generation passenger aircraft are using FRP composites in larger quant i ty to realise large 1 Chapter I Introduction gains in efficiency. In order to do so, the materials must go through extensive testing, which requires thousands of hours of experiment and modeling for each new material and application. To be competitive in the commercial aircraft market, the time between conception and roll-out does not allow every newly developed model to be rigorously tested. Complicated models that can explain every micro-mechanical aspect of a material's response under specific conditions may be excellent laboratory tools but are impractical for real applications. Each additional material parameter or variable in a model adds uncertainty. At the end of the day, a designer needs to be confident in his understanding of the physics of the problem and needs to be certain that a design will not fail. Each parameter or variable that is added to a model that is not physically based or that is added for numerical convenience hurts the understanding of the physical processes taking place. The most useful tools for a design engineer are simple models containing few parameters, completely based on physical quantities. Such models can be used to quickly narrow down a list of eligible materials for a given application. They can also be used to quickly test the validity of more complicated models that might be necessary to investigate specialized aspects of a select few prospective materials or structures. Completely new models were not developed in this thesis. Instead, a new way of looking at models that have been around for many years is presented with the goal of providing a clear picture of the mechanisms of fracture in FRP composites. A large number of experimental tensile fracture tests were conducted, which provide insight into the physical processes of fracture. A simple, three-parameter, numerical, cohesive zone model for orthotropic composites is used to help explain the results of the experiments. Although several other researchers have proposed models similar to the one presented here, the interaction of mesh size and input parameters has often necessitated the use of non-physical parameters. A better physical understanding of the problem eliminates these issues. Finally, a simple analytical model, using the same three parameters as the numerical model, is devised from basic fracture mechanics equations. This model provides a designer with a simple tool for evaluating different materials. 2 Chapter 2 Background and Literature Review CHAPTER 2 BACKGROUND AND LITERATURE REVIEW 2.1 Introduction This thesis studies the behaviour of notched composite laminates under applied tensile loads by means of (a) experimental evidence, (b) the development of a numerical composite damage model and, (c) a simple analytical description. The aim is to develop a fully consistent, physically based explanation of the behaviour from all three vantage points. There has been a considerable amount of work on damage and crack propagation in composites over the past several decades with numerous theories and models proposed and developed by a large number of researchers. Some of the models complement one another, some have been modified or improved upon by later researchers and others run parallel to each other. Major contributions to the study of crack and damage propagation from experimental evidence to closed-form solutions to numerical models are discussed here. 2.2 Fracture mechanics The theory of crack propagation in a continuum originated well before the development of advanced fibre-reinforced polymer composites. Griffith established the energy criterion for crack growth in brittle materials (Griffith, 1921) in 1921. Considering the energy required to create the surfaces of a crack, he determined the amount of energy required for a crack to grow. A crack will grow if the energy released during an incremental growth of the crack by length da is greater than or equal to the energy needed for the crack to grow by that amount, expressed as G>R. 2.1 G is called the strain energy release rate defined for unit thickness as the change in energy with crack extension such that G = f c l , 2.2 da 3 Chapter 2 Background and Literature Review where F is the work performed on a plate with free ends by external forces, and U is the stored elastic energy in the plate. R is the "crack resistance" equal to the energy absorbed during crack growth. While this definition was developed for brittle, elastic materials, it is valid for materials with plasticity or damage at the crack tip, although the crack resistance term is then dominated by plastic or damage energy as opposed to only surface energy. Evaluating the elastic energy in terms of the compliance of the plate, a general expression for the energy release rate can be written in terms of load and the compliance, as given by (Broek, 1982) G - * * . 2.3 IB da Alternatively, the condition for fracture can be expressed in terms of the stress or strain field at the tip of a crack. The well-known function for the stress field immediately around a crack tip is attributed to Irwin and given in many textbooks (e.g. (Broek, 1982)) as -JlTtr The stresses ax, ay, and zxy are defined in the directions given in Figure 2.1 in terms of the cylindrical coordinates r and 6. Kj is the mode I1 stress intensity factor describing the magnitude of the stress field. Along the crack plane ahead of the crack tip the function fy(0) is equal to unity. For a center crack in an infinite plate Ki is a function only of the stress in the plate and the length of the crack as given in Equation 2.5. Kj =<j4na~ 2.5 For finite sized structures, K[ also depends on geometry. Solutions for Kj have been compiled for standard test geometries (e.g. (Broek, 1982), (Tada, Paris and Irwin, 1985)) or can be determined through finite element techniques. Again, these equations are elastic solutions, valid if the plastic or damage zone at the crack tip is small relative to other ' There are three modes of cracking. Mode I is crack opening in which the crack surfaces move perpendicular to the crack plane. Modes II and III are shear modes (sliding and tearing) in which the crack surfaces move parallel to the crack plane. 4 Chapter 2 Background and Literature Review specimen dimensions. This means that the damaged material at the crack tip must be constrained by the bulk elastic material of the structure so that stress, strain and crack opening displacements do not differ greatly from the purely elastic case. A critical value of the stress intensity factor Kjc can be determined which corresponds to the stress at failure ac as Kic is called the fracture toughness of the material. If Kjc is known, the far-field stress at which fracture will occur can be determined. In Linear Elastic Fracture Mechanics (LEFM), Kjc is considered to be a material property. It seems that there are two independent approaches for determining the condition for crack initiation. However, the two conditions, energy and stress intensity, are linked. Broek (Broek, 1982) states, "The energy criterion is a necessary criterion for crack extension. It need not be a sufficient criterion. Even if sufficient energy for crack propagation can be provided, the crack will not propagate unless the material at the crack tip is ready to fail: the material should be at the end of its capacity to take load and to undergo further straining. " A relationship can be defined between Ki and Gj, if G/ is written in terms of stress. The full derivation is given by Broek and is only summarised here. The strain energy release rate is the amount of energy released during an infinitesimal length of crack growth. Alternatively, it can be viewed as the energy required to close a crack over an infinitesimal distance given as(Broek, 1982) where S is the infinitesimal length of crack and v is the opening displacement. Substituting Equation 2.4 and crack opening displacement in terms of K/ and E into Equation 2.7 results 2.6 2 r<Jyv G, =lim— dr 2.7 5 Chapter 2 Background and Literature Review 2 in the plane stress L E F M relationship between the stress intensity factor and the strain energy release rate as 2.8 Here E is the modulus for a homogeneous, isotropic material. It is often assumed that the stress criterion for fracture and the energy criterion for fracture are fulfilled simultaneously, e.g. (Broek, 1982; Kanninen and Popelar, 1985). If this is the case then K 2 G c = ^ . 2.9 E Indeed, Kanninen and Popelar (Kanninen and Popelar, 1985) write, "When Ki, for example, attains its critical value, then G must also reach its critical value and [Equation 2.8] implies [Equation 2.9]. Consequently, for linear elastic bodies, the stress intensity factor and the energy balance approaches to fracture are equivalent. " In L E F M , which assumes infinitely sharp crack tips, the simultaneous fulfilment of both conditions is assumed to be true and it is generally accepted for blunter crack tips as well. For blunt notch tips, holes and cut-outs the stress field around the stress concentrator is described by a stress concentration factor Kj instead of a stress intensity factor. While K/ describes the stress field as stresses approach infinity at a sharp crack tip, KT is a factor that relates the finite stress at a notch tip to the far-field stress. The stress field normal to the crack plane around a blunt notch tip can be related to the stress field of an equivalent sharp crack tip a certain distance away from the tip by an equation given by (Creager and Paris, 1967): a = K, p { W'\ K, 0' cos 2Kr' 2r\ 2 ) jj^? 2 + , cos-(, . 9' . W\ 1 + s i n — s i n — I 2 2) 2.10 2 This equation is valid for the case of plane stress. For plane strain the left side of the equation is multiplied by a (1-v2) term. 6 Chapter 2 Background and Literature Review The radius of curvature at the notch tip is p and r' is the distance from the effective notch tip as defined in Figure 2.2. In Equation 2.10, K\ is the stress intensity factor for a sharp crack tip located at a point p/2 behind the actual notch tip. Along the x-axis in front of the notch tip (8' = 0) equation 2.10 reduces to K, p K, cr2 = . ' . 1 . 2.11 As mentioned above, the modulus in Equation 2.8 is valid for isotropic materials. Sih, Paris and Irwin defined an effective modulus E', given by Equation 2.12, in terms of the elastic constants that is valid for anisotropic materials such as F R P composites (Sih, Paris and Irwin, 1965). E'= r 2.12 \aua22 \a22 + 2al2 +a66 V v a n 2au j For plane stress, au=\IEi, a2i=\IE2, an=-vnlEi and a66=VGi2, where Eu E2, V12 and Gn are the orthotropic elastic constants. 2.3 Cohesive crack models Another method of predicting crack initiation and propagation, that is fundamentally different from L E F M , is the cohesive crack model. The concept, originally proposed by Dugdale (Dugdale, 1960) for metals, has given rise to numerous analytical and numerical models for crack growth and damage in composite materials. The cohesive crack model recognises that there is a region of material ahead of a crack tip that is not linear elastic. In F R P composite materials this "damage zone" consists of matrix cracking, fibre breakage, interface separation and fibre pullout. In the cohesive crack approach, a relationship between traction and displacement is defined along a potential crack surface, sometimes called a "fictitious crack", ahead of an actual crack or notch tip as shown schematically in Figure 2.3. This fictitious crack represents the damaged area lumped into a discrete plane. A typical separation relationship is shown in Figure 2.4. The interface is linear-elastic up to a specified critical stress. For strains larger 7 Chapter 2 Background and Literature Review than the critical strain associated with this stress, the material is damaged. The damaged interface in Figure 2.4 obeys a linear relationship until the tractions on the interface are zero at some value of interface separation. The traction-displacement curve is often chosen to be linear for simplicity but is not required to be so. Dugdale used a cohesive crack approach to determine plastic zone sizes in steel panels containing edge slits. The concept was then developed further by Barenblatt (Barenblatt, 1962) to model crack propagation in metal sheets. Since that time, the cohesive approach has been used both as a continuum material model to predict damage growth and propagation, and also as strictly an interface model. When used as a continuum material model (called the cohesive zone model), the traction-displacement relationship becomes a stress-strain relationship. Important considerations have been brought to light from research in both areas and will be discussed in the following sections on interlaminar crack propagation and notched strength prediction. The link between fracture mechanics and cohesive models is the strain energy release rate. The area underneath the traction-displacement curve of a cohesive crack model is an energy release rate; essentially the amount of energy required to completely open the interface for a given increment of crack growth. The condition for complete crack separation in these models is the attainment of either a critical value of the energy release rate Gc or a critical interface displacement Sc. The two criteria are identical i f the bilinear constitutive curve is uniquely defined by one of these two values and by two of the three elastic parameters; EQ, apeak or Speak, where epeak is the strain at the peak stress on the curve. 2.3.1 Cohesive crack models for intralaminar crack propagation In most cases cohesive crack models are implemented in finite element computer codes as zero thickness elements. A common method of implementing the cohesive law is to use a penalty stiffness scaled by a damage parameter co after the peak traction is reached, for example (Alfano and Crisfield, 2001). As illustrated in Figure 2.4 for mode I loading, the basic constitutive relation is V , ifSi <<5,0 i/S,>Sk 2.13 0 8 Chapter 2 Background and Literature Review Here, Tis the interface traction defined in the i'h material direction, where /= 1,2,3. When the interface separation 8 is less than a threshold separation distance So, the interface is elastic. For 8 greater than So but less than a critical separation distance Sc, the interface is damaged by an amount co. After Sc the interface is fully damaged and the crack faces are completely separated. The damage parameter co e[0,l] is a monotonically increasing function, i.e. the amount of damage in the interface can never decrease. The damage parameter is a function of either interface separation as in (Camanho, Davila and de Moura, 2003), (Alfano and Crisfield, 2001), (Geubelle and Baylor, 1998), (Espinosa, Dwivendi and Lu, 2000) or a function of the energy release rate such as in (Borg, Nilsson and Simonsson, 2001; Borg, Nilsson and Simonsson, 2002), (Borg, Nilsson and Simonsson, 2002; Zou, Reid and Li , 2003). Another notable model is that of Xu and Needleman (Xu and Needleman, 1996; Xu and Needleman, 1994). In their model the interface constitutive relations are based on a work potential. The tractions on an interface are determined by differentiating the work potential, which is an exponential function of both normal and shear interface separations. Typically, separate traction-separation relationships are defined for opening (mode I) and shear (mode II and III) modes. Alfano and Crisfield ((Alfano and Crisfield, 2001)) investigated the problem of developing a coupled model for mixed-mode loading, wherein the critical strain energy release rate depends on the proportion of each mode of loading. 2.3.1.1 Shape of the traction-displacement curve Camanho et al. (Camanho, Davila and de Moura, 2003), Alfano and Crisfield (Alfano and Crisfield, 2001) and Geubelle and Baylor (Geubelle and Baylor, 1998) have all used a bilinear traction-displacement relationship. The formulations of Borg (Borg, Nilsson and Simonsson, 2001), Corigliano (Corigliano, Mariani and Pandolfi, 2003) and Xu and Needleman (Xu and Needleman, 1996; Xu and Needleman, 1994) allow for the definition of different curve shapes by varying an exponential parameter in the constitutive relation. The shape of the constitutive curve can play an important role in the response of the model. A curve such as that shown in Figure 2.5 a) can better represent the response of a ductile material while that in Figure 2.5 b) is better for modeling a quasi-brittle material, de Borst ((de Borst, 2003)) observes that for ductile fracture, the most important parameters of the cohesive zone model are the tensile strength and the fracture energy while for quasi-brittle 9 Chapter 2 Background and Literature Review materials, the shape of the stress-separation relation plays a larger role and can be even more important than the value of the tensile strength. He presents the example of a single edge-notched (SEN) tensile specimen modeled with different shaped linear and exponential constitutive softening curves but with the same input value of fracture energy. The resulting global peak loads differed by more than 15%. 2.3.1.2 Input parameter determination For a simple bilinear cohesive crack or cohesive zone model, only a few parameters are needed to define the constitutive behaviour and these are related to physical quantities. A bilinear constitutive curve is completely defined by three parameters: • The initial elastic modulus EQ • Either the peak stress apeak or the strain at peak stress epeak • Either the critical fracture energy Gc or the ultimate strain suu In many of the models reviewed in the literature however, the input parameters are not arrived at independently. Instead they are deduced by fitting the global response of the numerical model, such as the load versus displacement of a notched tensile test, to the results of an experimental test for one or two notch sizes. One of the reasons for this is the dependence of the structural response on mesh discretization when the models are implemented in finite element codes. The peak tensile stress in particular is often adjusted to match experimental results. It is generally agreed that the peak stress does not have a large effect on crack growth (Zou, Reid and L i , 2003) but that proper selection of apeak is important for a good prediction of crack initiation (Alfano and Crisfield, 2001). Geubelle and Baylor noted in their study on impact-induced delamination of composites that the time and location of initial matrix damage depended strongly on apeak (Geubelle and Baylor, 1998). Zou et al. (Zou, Reid and L i , 2003) and Alfano and Crisfield (Alfano and Crisfield, 2001) also noted that the value of the peak stress has a strong influence on the computational efficiency of the finite element simulation, in that the higher apeak, the more refined the mesh must be in the region of the delamination front. Usually, o-peak is chosen to be reasonably close to the physical interlaminar strength of the interface, however, depending on the element size used to mesh the crack plane, an incorrect 10 Chapter 2 Background and Literature Review unstable solution is sometimes obtained when using this value for <Jpeah Therefore, it is sometimes altered so that a stable solution is achieved. Borg et al. (Borg, Nilsson and Simonsson, 2001) noted that when modeling mixed mode delamination in carbon-fibre/ epoxy, if peak stresses equal to the failure stresses in uni-axial tensile tests were used, the peak stresses were of magnitudes that prevented the growth of delamination. Consequently, these researchers used peak stresses substantially lower than the interlaminar strengths of the material. Both Xu and Needleman (Xu and Needleman, 1996; Xu and Needleman, 1994) and Geubelle and Baylor (Geubelle and Baylor, 1998) chose (jpeak as a fraction of the elastic modulus instead of using experimental values. Most of the models are implemented in explicit dynamic finite element codes and several of the models are used to investigate dynamic crack propagation. For most of the models though, the material parameters are independent of strain rate. To model dynamic crack propagation, Geubelle and Baylor (Geubelle and Baylor, 1998) pre-multiplied the input fracture toughness by a strain rate parameter. Corigliano, Mariani and Pandolfi (Corigliano, Mariani and Pandolfi, 2003) included a strain-rate dependent parameter as part of the mode I cohesive law in their model and compared it to a rate-independent model. They concluded that rate-dependency can greatly affect dynamic crack processes but that the effects of rate dependency on the test response during dynamic delamination growth decrease when inertial terms become dominant. 2.3.2 Cohesive zone models for interlaminar crack propagation The concept of the cohesive crack model has also been expanded from an interfacial model to simulate the continuum response of a damaging material. Ffillerborg et al. (Hillerborg, Modeer and Petersson, 1976) appear to be the first to apply the concept of the cohesive crack method of Dugdale and Barrenblatt to a continuum finite element model. Applying the cohesive crack model to a solid continuum introduces certain problems. In the cohesive crack model all damage is restricted to a discrete plane, but the cohesive zone model attempts to represent a volume of damaged material with elements that may not be of the same size as the material's characteristic damage size. The response of the model therefore depends on a length scale. 11 Chapter 2 Background and Literature Review In a cohesive zone model solid elements are used as opposed to zero-thickness interface elements. A stress-strain relation replaces the relationship between traction and displacement. The area beneath the stress-strain curve is a specific energy y (energy per unit volume). The specific energy is multiplied by the length scale to obtain a fracture energy release rate GF (i.e. energy per unit area) that can be related to fracture mechanics and known material properties. In finite element simulations damage tends to localise in a single layer or row of elements. Therefore the element height he is the length scale and GF is calculated as GF=yhe. 2.14 Floyd (Floyd, 2004) has investigated localisation in a composite damage model and reviews the major work done in this area. Damage localisation in a cohesive zone model is demonstrated in Figure 2.6 for a simple three-element system. A l l three elements are defined using the same constitutive model. Initially there is zero strain in each element (a). As the applied displacement is increased, the stress in each element increases along the slope defined by the initial modulus until the strain corresponding to the peak stress <jmax is reached (b). One of the elements will reach this critical strain slightly before the other two and will begin to soften. As softening occurs in the middle element, its residual modulus is reduced, thereby reducing the stress in the entire structure. With the reduction in stress, the load in the other two elements decreases and they shed their elastic energy on to the softening element as they unload along the slope of the initial modulus (c). The unstable system quickly leads to complete damage in the softening element (d). At this point there is zero strength in the softened element and the crack is considered to be fully developed at that location. Although only a three-element example is presented here, the effect is the same for stacks of larger numbers of elements. The implication of localisation is that the damage height effectively becomes the height of a single element and the energy release rate of the material is dependent on that element height. A means of dealing with localisation and the mesh height effect was introduced by Bazant et al. (Bazant, 1984; Bazant and Oh, 1983; Bazant and Oh, 1984) for concrete. In their approach, termed the crack band model, the damaged zone is represented by a band of distributed cracks. The fracture energy release rate GF in a numerical model is a constant related to the height of a layer of elements and the area under the stress-strain curve in the 12 Chapter 2 Background and Literature Review numer ica l model (as in Equat ion 2.14). Physical ly, i t is related to the physical material parameters: the characteristic height o f the damage zone hc for the material system and the specif ic fracture energy yc, therefore GF=ye-he=yc-hc. 2.15 Numer ica l l y , GF can be held constant in a s imulat ion for vary ing he on ly by adjust ing ye, since yc and hc are material parameters. F l o y d (F loyd , 2004) has implemented a crack band scheme in a composite damage model that scales the const i tut ive softening curve for vary ing he wh i l e mainta in ing GF constant. Cohesive zone models appl ied to predict notched strength are discussed in Section 2.4. 2.3.3 Mesh Sensitivity One o f the ma in drawbacks to the cohesive crack model is that the path o f a propagat ing crack is mesh dependent. The crack is forced to g row along the interface where the cohesive law has been def ined. For continuous f ibre re inforced po lymer laminates, this p rob lem does not create large restrict ions because delaminations are k n o w n to f o r m in the inter laminar layers. However , to a l low for delaminat ion to f o r m on any inter laminar layer for a structure in w h i c h the delaminat ion locat ion is not k n o w n a priori, cohesive interface elements should be def ined between al l layers o f the model . I n a cont inuum cohesive zone model , the cohesive law is def ined w i t h i n al l elements a l l ow ing cracks to ini t iate and develop at any locat ion. I t was shown in Section 2.3.2 that for cohesive zone models, energy dissipation depends on the height o f the elements used to discretize the structure for cohesive zone models in w h i c h vo lumetr ic elements are used. Interface cohesive crack models do not have this drawback because the const i tut ive law is def ined in terms o f length-independent quantit ies; t ract ion and displacement. The solut ion for both types o f models, however, is affected by the spatial discret izat ion in the direct ion paral lel to the crack. A s a m i n i m u m requirement, the discret izat ion must be f ine enough to resolve the gross shape o f the stress field around the stress concentrator. The cohesive models account for damage at the notch t ip so there is no stress singular i ty in the stress field as predicted by L E F M . 13 Chapter 2 Background and Literature Review Di f ferent authors have found vary ing degrees o f sensit iv i ty to mesh size along the crack direct ion. A t one extreme, Geubelle and Bay lo r (Geubel le and Bay lor , 1998) conducted a sensit iv i ty study by look ing at the convergence o f the speed o f a rap id ly propagat ing crack in a pre-loaded P M M A strip for various mesh sizes. They found that the fracture speed converged for elements w i t h an approximately 7 pm characteristic d imension. They concluded that the elements needed to be t w o or three t imes smaller than their estimate o f the cohesive zone size for convergence o f the scheme. Conversely, B o r g et al. (Borg , N i lsson and Simonsson, 2002) found that they could use a relat ively coarse mesh and achieve load versus displacement results comparable to those o f a f ine mesh. However , in their model they found i t necessary to decrease the input values o f strength and Gc f r o m actual experimental values by 5 0 % or more. A l f a n o and Cr is f ie ld (A l f ano and Cr is f ie ld , 2001) also conducted a study on the effect o f vary ing mesh size and tensile strength simultaneously. They found that when large tensile strength values (comparable to experimental values) were used, they required a very ref ined mesh to achieve good load versus displacement results for a D C B specimen. They also noted that they cou ld use a coarse mesh to achieve good results i f the model tensile strength was decreased (by up to a factor o f 30). X u and Needleman ( X u and Needleman, 1994) noted that when using a fine mesh, decohesion takes place in a more or less duct i le manner wh i l e larger elements led to a more br i t t le response. Since most FRP composites have damage zones that are several mi l l imet res in height, it is often necessary to use elements that are somewhat smaller than the damage zone height so that stress fields can be proper ly represented. The crack band scaling method developed by F l o y d takes care o f this aspect. For the s imulat ion o f engineering structures however, too fine a mesh can lead to enormous run t imes. A l so in a cohesive zone model , in w h i c h fibre and mat r ix properties are smeared into a cont inuum, it makes sense to use elements that are on the order o f the size o f the damage zone. I f elements are smaller they begin to represent volumes compr is ing on ly fibre or on ly matr ix . A method o f scal ing the material input parameters for re lat ively large element wid ths (on the order o f mi l l imetres or hundreds o f microns as opposed to microns) is presented in Chapter 4. 14 Chapter 2 Background and Literature Review 2.4 N o t c h e d s t r e n g t h p r e d i c t i o n Numerous models based on strength o f materials, fracture mechanics and cohesive methods have been proposed over the years that attempt to predict the strength o f structures containing stress concentrators such as holes and notches. 2.4.1 Semi-empirical models I n the 1970's and early 1980's in part icular, several models were proposed. A w e r b u c h and Madhukar (Awerbuch and Madhukar , 1985) wrote a comprehensive rev iew o f the major models up to 1985. The essential works discussed in their rev iew are: 1) Waddoups, Eisenmann and K a m i n s k i : L E F M approach 2) Wh i tney and Nuismer : Point stress and average stress methods 3) Kar lak : A mod i f ied po in t stress method 4) Pipes, Wetherho ld and Gi l lespie: A mod i f ied po in t stress method 5) M a r and L i n : Composi te stress d is t r ibut ion approach M o s t studies have investigated the response o f plates contain ing either sharp slits or c i rcular holes where the notch size a to specimen w i d t h W was variable. A few studies have looked at the effect o f vary ing the notch t ip radius wh i le mainta in ing a/W as a constant. Results are typ ica l ly presented as the ratio o f notched strength to unnotched strength. The Waddoups, Eisenmann and K a m i n s k i ( W E K ) (Waddoups, Eisenmann and K a m i n s k i , 1971) mode l is a semi-empir ical model in that i t uses L E F M equations along w i t h a characteristic mater ial length that must be determined f r o m f i t t ing curves to experimental data. The characteristic length ac defines what the authors cal l an intense energy region emanating f r o m the edge o f a hole in a plate. This region is considered to be an equivalent crack so that a stress intensity factor for a crack emanating f r o m a hole can be def ined as Here p is the hole radius and for a specimen contain ing no hole the funct ion f(ac/p) is equal to un i ty . The assumption is then made that f(ac/p) is a constant equal to the ratio o f the 2.16 strength o f the control specimen to the strength o f a notched specimen B y 15 Chapter 2 Background and Literature Review conduct ing experimental tests w i t h a notched and an unnotched specimen, the model can be cal ibrated and used to predict strengths for other specimens o f s imi lar geometry. Due to its semi-empir ical nature, good agreement w i t h experimental results was achieved. I n independent tests by A w e r b u c h and Madhukar however, i t was found that the characteristic length was dependent on hole radius in addi t ion to mater ial system, p l y or ientat ion and stacking sequence. F o l l o w i n g the W E K model , Wh i tney and Nu ismer (Wh i tney and Nuismer , 1974) postulated t w o stress cr i ter ia for predic t ing notched strength o f composite laminates that have received considerable attention. The "po in t stress" and "average stress" models also incorporate a characteristic distance at the notch t ip. Fai lure is predicted to occur when the stress at some distance f r o m the notch t ip reaches the unnotched strength o f the mater ial . Instead o f us ing an L E F M solut ion for the stress dist r ibut ion at a sharp notch t ip , the models use the elastic solut ion for the stress d is t r ibut ion at a hole. The theoretical basis o f their approach is that the material must be cr i t ica l ly stressed over a characteristic d imension in order to find a suff ic ient f l aw size that w i l l ini t iate fai lure. I n the point stress model , the stress at a po in t a distance do f r o m the notch t ip along the axis o f the notch is compared to the unnotched strength. The authors consider do to be a material property. The average stress model is the same as the po in t stress model except that the average stress over a distance ao f r o m the notch t ip is compared to the unnotched strength. A s w i t h the W E K model , values for do and ao must be determined by conduct ing experiments and f i t t ing the analyt ical curve to experimental data. I n so do ing they achieve good results but find that the characteristic dimensions vary s l ight ly w i t h hole size. Subsequent investigations o f the model , e .g. (Awerbuch and Madhukar , 1985), ( K i m , K i m and Takeda, 1995), (Tan, 1988), also found that do and ao depend on material system and laminate geometry. Separate research by Kar lak and by Pipes, Wetherho ld and Gi l lespie ( P W G ) (Pipes, Wetherho ld and Gi l lespie, 1979) found that do and ao depend on the material system and also on the notch or hole size in some cases. Kar lak proposed a square-root relat ionship between hole radius and do using a new parameter kg. The determinat ion o f ko must be done exper imental ly for each new material system that is modeled. This adjustment improved the results in most cases. Investigations by P W G led to a more general exponential relat ionship between R and do. 16 Chapter 2 Background and Literature Review Tan (Tan, 1988) developed the poin t stress and average stress models o f Wh i tney and Nu ismer further using classical laminated plate theory. The models are constructed on the basis that fracture in a composite laminate is fibre control led. A n effect ive fibre strength in the region o f a notch t ip containing heavy mat r ix damage is calculated and fai lure is predicted to occur when the stress in any p l y meets this effect ive strength. The models require the normal stress d is t r ibut ion, longi tudinal and compressive strength parameters o f a lamina and, again, a characteristic length that must be determined from an experimental notched specimen. A l l o f the models based on the po in t stress and average stress models general ly produce good results, w h i c h is not unexpected as the models are semi-empir ical and require experimental testing to determine the characteristic dimensions. Add i t i ona l l y , the characteristic dimensions do not represent physical quantit ies and are therefore d i f f i cu l t to jus t i f y . 2.4.2 Cohesive models M u c h o f the recent w o r k on notched strength predict ion has incorporated the use o f cohesive models. Cohesive models are general ly more phys ica l ly based, tak ing into account the development o f a damaged zone (or process zone) ahead o f a notch t ip. I n 1986, Back lund and Aronsson (Back lund and Aronsson, 1986) presented a two-parameter damage zone model ( D Z M ) for notched composite laminates based on the cohesive zone model o f H i l le rborg . I n addi t ion to determining crack in i t ia t ion, the length o f damage zone extension is also obtained. The model is implemented as a material model in a f in i te element code. The const i tut ive curve is considered to be r ig id unt i l crack in i t ia t ion so that the on ly parameters needed to define the l inear stress-displacement curve are the unnotched tensile strength and an apparent fracture energy Gc*, w h i c h is the area beneath the stress-displacement curve. Centre-notched quasi- isotropic carbon/epoxy laminate specimens containing circular and non-ci rcular holes were tested exper imental ly by Back lund and Aronsson and predict ions were made using the D Z M model w i t h 0.25 m m elements and the W E K model . I t should be noted that Gc* (=35 kJ /m ) was chosen to match numerical predict ions to the average experimental fracture load for five notched specimens and not determined independently. The D Z M produced s l ight ly better results than the W E K model w i t h an average deviat ion from the experimental peak loads o f about 3 .3% for twelve 17 Chapter 2 Background and Literature Review specimens o f t w o di f ferent notch sizes. A l t hough i t is not k n o w n h o w close the apparent fracture energy is to the actual strain energy release rate o f the mater ial , Back lund and Aronsson noted that a var iat ion in Gc* o f 10% resulted in a var iat ion in the predicted peak load o f about 3%. Interest ingly, they also note that when model ing square holes that produce a s ingular i ty in the stress field, values o f damage in i t ia t ion stress approach zero as element size is reduced but that the predicted peak load converges for element lengths on the order o f 0.25 to 0.5 m m . Recent ly Shin and W a n g (Shin and W a n g , 2004) have also used a cohesive zone model to predict notched fai lure in open-hole tension tests. Improv ing upon the method o f Back lund and Aronsson they use a three-parameter const i tut ive curve w i t h the effect ive modulus for anisotropic materials as def ined by Sih et al. (S ih, Paris and I r w i n , 1965) in addi t ion to the cr i t ical fracture energy Gc and material tensile strength crc. They also prov ide a good comparison o f notched strength predict ions from the cohesive zone model ( C Z M ) to point stress and average stress model predict ions. The C Z M predicted fai lure loads s l ight ly closer to experimental results. The authors also attempt to measure the w i d t h o f the damage zone ahead o f the notch t ip in experimental specimens and correlate i t to the length o f damaged material predicted by the C Z M before fai lure. The specimens were center-notched tensile specimens (210 m m x 25.4 m m ) w i t h four hole diameters ranging from 3.0 to 12.0 m m . The material systems used were both quasi- isotropic and cross-ply thermoplast ic A S 4 / P E E K and thermosett ing T300/3501 laminates. Unfor tunate ly , the specimen geometry they used is unstable and they were unable to g row a fu l l y developed damage zone. They therefore stopped the tests at 95 % o f the estimated fai lure strength and used radiographic and sectioning techniques to determine damage zone widths in f ront o f a 3 m m diameter notch. For the quasi- isotropic laminates the w i d t h o f the damage zone was approximately 2 m m and for the cross-ply laminates i t was approximately 0.55 m m . For a constant stress level , the w i d t h o f the damage zone is predicted to increase w i t h increasing hole diameter. Since most notched strength studies are conducted w i t h center-notched or edge-notched specimens containing circular cut-outs, the net section area o f a specimen changes w i t h changing hole radius. The fai lure strength o f the specimen is often more sensitive to a reduct ion in the net section than to a change in the notch radius. T o isolate the effect o f 18 Chapter 2 Background and Literature Review notch radius on tensile strength on ly a few researchers have maintained a constant crack length a and constant a/W wh i l e changing the notch t ip radius p (e.g. (H i tchen et al . , 1994), (Hyakutake and Hagio , 1990), ( K a m i y a and Sekine, 1997)). H i tchen et al. (Hi tchen et al . , 1994) used centre-notched tension specimens (160 m m x 25 m m ) but kept a (and therefore a/W) constant. No tch t ip radi i were var ied by d r i l l i ng holes w i t h p o f 0.5 m m , 1.0 m m , 2.0 m m and 5.0 m m at each end o f a central sl it so that 2a was maintained at 10 m m . Hyakutake and Hag io (Hyakutake and Hagio , 1990) conducted experiments w i t h double edge-notched bars where the notch radi i were var ied f r o m 0.08 m m to 10 m m wh i le mainta in ing a constant a o f either 1 m m or 5 m m . K a m i y a and Sekine ( K a m i y a and Sekine, 1997) used a compact tension specimen geometry w i t h values for p o f 0.25 m m , 2.0 m m , 5.0 m m and 10 m m . General ly speaking, a t rend o f increased fai lure load w i t h increased notch radius was observed. Results presented by K a m i y a and Sekine contradict ing this t rend are, as pointed out by the authors, most probably due to holes that were too large compared to the specimen size. I n al l cases the combinat ion o f specimen geometry and material system was unstable so that on ly fracture loads could be determined and not the stable damage zone size. Hyakutake and Hag io (Hyakutake and Hag io , 1990) introduced a concept they cal l " l inear notch mechanics". They propose that the strength o f a notched specimen is related to the stress present at the notch root when the fai lure load is reached. They consider this stress o'max.cip) to be a funct ion on ly o f the notch root radius. Then the nomina l stress in the specimen at fa i lure is where Kj is the stress concentrat ion factor. The stress distr ibut ions in the specimens at the fai lure load were determined using a f in i te element analysis and several experimental notched tests are required to develop a characteristic curve for a material system. This method o f notched strength predict ion seems questionable as i t says the fai lure cr i ter ion for re lat ively sharp notches is based solely on stress at the notch t ip w i t h no consideration g iven to the energy necessary for fracture. 19 Chapter 2 Background and Literature Review 2.4.3 Relation of the damage zone size to strength Other models that attempt to predict the strength o f notched specimens are based on a cr i t ical height o f the damage zone. The cr i t ical height refers to the height o f the damage zone when the fai lure load is reached. Typ ica l l y a shear-lag or statistical analysis is used to determine what the damage zone height w i l l be and then a relat ionship invo lv ing the damage zone height and the fai lure load is def ined. Dharani et al. (Dharani , Jones and Goree, 1983) developed strength solutions for unidi rect ional laminates containing slits or rectangular or circular notches. A shear-lag analysis is used to determine the remote stress when the stress in a fibre at the notch t ip reaches the unnotched strength o f a fibre. K a m i y a and Sekine ( K a m i y a and Sekine, 1997) also used a micro-mechanical analysis to determine the stresses in transverse fibres at a notch t ip . B y mode l ing a single fibre as a p u l l -out p rob lem, they calculate the length o f fibre that w i l l pu l l out f r o m the surrounding mat r ix and then determine the probabi l i ty o f breakage o f a fibre o f that length. The upper bound o f the length o f delaminated area perpendicular to the crack is analogous to the damage zone height. One o f the simplest phys ica l ly based models for strength predic t ion i n notched laminates is that o f Kor tschot and Beaumont (Kor tschot and Beaumont, 1991; Kor tschot , Beaumont and Ashby , 1991; Kor tschot and Beaumont, 1990a; Kor tschot and Beaumont, 1990b). They investigated the size and shape o f the sub-cri t ical damage zone in double-edge notched ( D E N ) tensile specimens and use fracture mechanics equations to predict strength based on the m a x i m u m stress experienced by the 0° plies at the notch t ip. Us ing a radiographic technique, these authors were able to measure the size and shape o f the damage zone dur ing loading, pr ior to the onset o f fast fracture. For the cross-ply graphite/epoxy specimens studied, the damage zone consisted o f spl i t t ing in the 0° pl ies, transverse cracks in the 90° pl ies and a tr iangular delaminat ion area ahead o f the notch t ip as shown in Figure 2.7. They found that increasing damage zone height (i.e. split length) led to higher specimen strength. Damage height decreased w i t h increasing notch w i d t h and therefore strength also decreased w i t h increasing notch w id th . Us ing a static, plane-stress f in i te element model incorporat ing the orthotropic and layered properties o f the laminate to determine the stress concentrat ion, 20 Chapter 2 Background and Literature Review they determined a semi-empir ical relat ionship for KT in terms o f the damage zone height and notch w id th . Notched specimen strength could then be related to the tensile strength o f a 0° p ly . A f te r not ic ing a size effect in specimens o f di f ferent w i d t h , Kor tschot and Beaumont inc luded a W e i b u l l d ist r ibut ion term in the model to account for it. They investigated the effect o f lay-up and notch w i d t h on the notched strength but d id not consider notch root radius. Thei r unstable specimen geometry prevented the development o f a stable damage zone. 2.5 A specimen for stable damage growth A specimen geometry for stable crack g rowth was developed by Kongshavn and Poursartip (Kongshavn and Poursart ip, 1999). They found that by increasing the height o f a compact tension ( C T ) specimen they were able to g row damage stably i n the specimen. The increased size o f the Over-height Compact Tension ( O C T ) specimen ensured that the boundaries d id not interfere w i t h the re lat ively large damage zone in the quasi- isotropic FRP laminate. The stabi l i ty o f the system can be described in terms o f fracture mechanics. The equ i l ib r ium state o f crack propagat ion was g iven by the equal i ty in Equat ion 2 . 1 . Stabi l i ty o f the equ i l ib r ium state is determined by the rate o f increase o f G w i t h respect to a compared to the rate o f increase o f R w i t h respect to a. The three possibi l i t ies are: dG dR + ,. — < — stable da da 3G 95 . . . . . . . — = — cr i t ical stabi l i ty 2.18 da da dG dR — > — unstable da da I f i t is assumed that crack resistance is a material property as in L E F M then R = Gc is a constant and the stabi l i ty condi t ion is ^ < 0 . 2.19 da The change in energy release rate w i t h crack g rowth can be obtained by di f ferent iat ing Equat ion 2.3 y ie ld ing 21 Chapter 2 Background and Literature Review 8G P2 82C — = - . 2.20 da IB da2 For C T and O C T specimens — - is posi t ive, al though the increased height o f the O C T da specimen stiffens the arms, w h i c h decreases the change in compl iance w i t h crack g rowth . Under displacement control the O C T can be stable i f the load drop w i t h a small increment o f crack g rowth is enough to make G < R. For tough material systems and for sharp notches under displacement contro l , load drops w i t h crack extension do b r ing G be low R, causing crack arrest. For br i t t le material systems and other notch geometries the specimen is no longer stable as w i l l be discussed later. Kongshavn and Poursart ip, and later M i t c h e l l (M i t che l l , 2002), used the O C T geometry to study the size and shape o f the characteristic stable damage zone w i t h the intent o f developing a strain-softening damage model for composite laminates. 2.6 Summary • I n fracture mechanics, the onset o f fracture in a cont inuum can be predicted based on an energy cr i ter ion (G > R) or a stress cr i ter ion (K/ > Kic). For elastic materials, the t w o cr i ter ia can be l inked using Equat ion 2.8. • Cohesive crack models are often used to model g rowth o f fracture and damage in composite materials. These models are mesh sensitive. Due to local isat ion and issues w i t h energy dissipation, they do not give converged solutions w i t h increased mesh ref inement. • A l t h o u g h cohesive crack and cohesive zone models use simple, physical input parameters, interact ion o f the input parameters w i t h mesh size is not w e l l understood and leads to parameters being adjusted to f i t experimental results. • Simple, semi-empir ical models for predic t ing notch fai lure have also been proposed, notably the point stress and average stress cr i ter ions, but these models also require the f i t t ing o f the analyt ical curves to one or more sets o f experimental data. • A l l experimental investigations o f notched strength, w i t h the exception o f (Kongshavn and Poursart ip, 1999), that were rev iewed used unstable specimen 22 Chapter 2 Background and Literature Review geometries such as centre-notched tension, double edge-notched tension or compact tension tests. The researchers were therefore unable to make observations on crack growth or the stable height of the damage zone. Damage zone size at failure was predicted to depend on notch size and geometry. 23 Chapter 2 Background and Literature Review Figure 2.3 Cohesive crack model 24 Chapter 2 Background and Literature Review T &max 1 do 8C 5 Figure 2.4 Traction — displacement relationship for cohesive crack model a) b) Figure 2.5 Cohesive crack models for a) ductile solids and b) quasi-brittle solids (from (de Borst, 2003)) 25 Chapter 2 Background and Literature Review (a) (b) (c) (d) Figure 2.6 Damage localisation using a cohesive zone model 26 Chapter 2 Background and Literature Review Chapter 3 Experiments C H A P T E R 3 EXPERIMENTS 3.1 Introduction T w o di f ferent specimen geometries were tested dur ing the course o f this study to investigate the behaviour and characteristics o f the damage zone in a f ibre-re inforced po lymer composite and to obtain material property inputs for the numer ica l model . The t w o geometries tested were a mod i f ied Over-height Compact Tension ( O C T ) specimen and a four-po in t bend specimen. The O C T specimens were mod i f ied to include circular notch t ips o f vary ing radi i . The p r imary objectives in testing the O C T specimens were to: 1) Obtain a relat ionship between peak load and notch t ip radius to val idate the numerical mater ial model 's abi l i ty to cope w i t h stress fields o f vary ing magnitude. 2) Produce stable crack g rowth so that the energy release rate for the structure can be determined. 3) Investigate the size and shape o f the process zone. The four-po in t bend specimens were tested p r imar i l y to obtain the f lexural strength o f the material system. 3.2 Specimens and procedures 3.2.1 Material The material used for both the O C T tests and the four-point bend tests was a quasi- isotropic carbon fibre/ epoxy laminate designated B A - T - 7 and manufactured by The Boe ing Company, Seattle, Washington, U S A . I t is the same material as used b y M i t c h e l l (M i t che l l , 2002) to investigate the notched behaviour o f stitched and unst i tched R F I CFRP laminates. The laminates used in this study are manufactured using a resin film infus ion (RFI ) process and are unreinforced (unsti tched) through the thickness. The fibres are a mix ture o f standard and intermediate modulus f ibres and the epoxy resin is 3501-6. The lay-up orientat ion is [+45/-45/02/90/02/-45/+45J6. The or ig ina l ly determined elastic properties for the material system are g iven i n Table 3 .1 . 28 Chapter 3 Experiments 3.2.2 Test equipment A n M T S hydraul ic uniaxia l testing machine w i t h an Instron control ler (Mode l# 8500R.) was used for both the O C T tests and the four point bend tests. The load cel l ( M T S M o d e l # 661.21A-02) has a 50 k N capacity. Separate testing j igs , designed and bui l t in-house, were used for each type o f test. Mach ine contro l and data acquisi t ion were conducted w i t h a dedicated PC using Instron Wavemaker 7.0 software. 3.2.3 OCT tests The Over-height Compact Tension ( O C T ) specimen geometry developed by Kongshavn and Poursart ip (Kongshavn and Poursart ip, 1999) was chosen to produce stable crack g rowth i n a specimen so that the composite damage zone could be investigated and to a l low for val idat ion o f the numer ical model . Kongshavn demonstrated that the O C T specimen w i t h a sharp notch t ip is stable under displacement contro l and is large enough so that the boundaries do not greatly affect the damage zone size or shape. I n this study the O C T specimen is s l ight ly mod i f ied by cut t ing circular notches w i t h vary ing radi i at the notch t ip . I t w i l l be demonstrated later that this has an effect on the stabi l i ty o f crack g rowth i n the specimen. 3.2.3.1 Specimens The O C T specimen geometry w i t h nomina l measurements is shown in Figure 3 .1 . The notch root radius was var ied f r o m a relat ively sharp radius o f 0.5 m m to a very large radius o f 27 m m . Several o f the specimens are displayed in Figure 3.2 to i l lustrate the dif ference i n the size o f the notches. The panel thickness B had a nomina l value o f 8.5 m m but var ied between 7.9 m m and 9.0 m m due to inconsistencies in the panel manufactur ing. The O C T specimens were cut f r o m a large panel so that the material 0° d i rect ion was paral lel to the notch as shown i n Figure 3 .1 . I t was important to have perfect ly c i rcular and undamaged notch t ips. They were cut at l o w speed w i t h d iamond sintered saw blades to reduce damage and no sanding or f in ish ing was required. A complete l ist o f the O C T tests conducted is g iven in Table 3.2 along w i t h the important specimen dimensions. The specimens are numbered corresponding to the diameter o f the notch t ip cut-out. General ly t w o tests were conducted at each notch root radius. A l l tests were done under displacement control . The first test at each notch root radius was displaced 29 Chapter 3 Experiments at a constant rate to an ul t imate displacement and then unloaded. The second test o f each set included unloads at t w o or three locations along the softening por t ion o f the load-displacement curve. The tests that included per iodic unloads are noted in Table 3.4. 3.2.3.2 OCT test set-up and procedure The O C T test set-up is shown in Figure 3.4. The specimen is loaded in tension through pins located above and be low the notch. To prevent tw is t ing o f the specimen (as noted by M i t c h e l l (M i t che l l , 2002)) , a st i f fening support was attached to the back edge o f the specimen. I t has l i t t le effect on the response o f the specimen and loading was stopped before the crack grew into the region covered by the stiffener. For al l tests an extensometer ( Instron 2620-825 w i t h ±5 m m travel) was f i xed between the upper and lower edges o f the notch mouth as shown in Figure 3.4 to measure crack mouth opening displacement ( C M O D ) . Smal l grooves were cut i n the specimen to prevent s l ipp ing o f the knife-edges. I t was realised after some in i t ia l tests that compl iance in the loading j i g created a dif ference in the displacement o f the machine and the displacement at the loading pins as can be seen in Figure 3.5. The precise displacement o f the pins is necessary in the procedure for determining Gc. Therefore, for the remainder o f the tests, the relat ive displacement o f the pins was also measured using an extensometer attached between grooves on the loading pins. The displacement rate dur ing load and unload for al l tests was 0.508 m m / m i n (0.02 in /min) . The machine cross-head was displaced unt i l approximately 4.5 m m o f C M O D and then displacement was reversed to the po in t o f zero load. Since there was some residual deformat ion, the po in t o f zero load d id not coincide w i t h zero displacement. For certain tests (noted in Table 3.2) part ial unloads were conducted at intervals along the softening por t ion o f the l o a d - C M O D curve. I n these cases the specimens were unloaded to 5 0 % o f the load at the in i t ia t ion o f the unload cycle. For al l tests the values recorded dur ing testing were load at the pins, machine posi t ion and C M O D . For some tests the relat ive p in displacement was also recorded. The data output frequency was 10s" 1 . Photographs o f the specimen were taken dur ing each test using a N i k o n d ig i ta l S L R D 1 0 0 camera w i t h a 105mm macro N i k o r lens mounted on a stationary t r ipod in front o f the test 30 Chapter 3 Experiments set-up. The photographs were taken in large (6.1 mega-pixel ) . jpg format. Photographs were taken at regular intervals dur ing in i t ia l loading before damage, f requent ly near the peak load and f o l l o w i n g each incremental load drop dur ing crack propagat ion. For each photograph the C M O D was noted and, in addi t ion, the internal c lock on the camera was synchronised w i t h the data acquisi t ion t im ing so that each photograph could be placed on the l o a d - C M O D curve dur ing post-processing. For each specimen approximately 25 to 40 photographs were taken depending on the durat ion o f loading. 3.2.3.3 Post-processing 3.2.3.3.1 Line analysis The l ine analysis technique used by Kongshavn and Poursart ip (Kongshavn and Poursart ip, 1999) and M i t c h e l l (M i t che l l , 2002) to obtain the crack opening displacement ( C O D ) pro f i le is also used in this invest igat ion. Before test ing, each specimen was prepared by scr ib ing hor izontal l ines 2.5 m m apart, paral lel to the notch in the region ahead o f the notch t ip as shown in Figure 3.3. The lines are used to provide points o f reference for measuring the crack opening pro f i le f r o m photographs taken dur ing the test. A photograph o f the specimen at a g iven C M O D dur ing the test is compared to a photograph o f the specimen at the start o f the test. The relat ive displacements o f points along the scribed lines are measured. Us ing several photographs throughout the test, C O D prof i les at di f ferent stages o f crack propagat ion are obtained. Software developed in-house was used to make the measurements. Results from the l ine analysis are presented in Section 3.3.1.3. 3.2.3.3.2 Sectioning Cross-sections were cut perpendicular to the crack plane in front o f the in i t ia l notch t ip for certain specimens to produce a pro f i le o f the damage zone along the length o f the crack. A slow-speed d iamond sintered cut t ing wheel w i t h a 0.32 m m (0.013 in) th ick blade was used to produce a very smooth surface that d id not require po l ish ing. Sections were cut every 0.635 m m to 1.905 m m and the remain ing surface was photographed w i t h the dig i ta l camera and a 105 m m macro lens using large . t i f f format for very h igh resolut ion. A 6.35 m m (0.25 in) marker was placed in each photograph so that the images cou ld be calibrated and the height o f the damaged area could be measured using image analysis software. The damage height hc was measured from the t ip o f v is ib le damage on one side o f the crack plane to the t ip o f v is ib le damage on the other side. Damage i n the outer pl ies was not included in the 31 Chapter 3 Experiments measurement o f damage height because i t could be either extensive or not depending on the amount o f resin on the surface. For example, considerably more surface damage was apparent on the top surface o f the laminate (away f r o m the resin f i l m in the R F I process) than on the bot tom surface. For consistency, i t was decided not to include any damage in the pl ies outside o f the 90° p l y in the outer stacks as shown in Figure 3.6. Results f r o m sectioning are presented in Section 3.3.1.2. 3.2.4 Four-point bend tests Four po in t bend tests were conducted in addi t ion to the O C T tests to obtain a value o f material strength for input into the numer ical material model . The A S T M Standard Test M e t h o d for Flexural Properties o f Unre in forced and Reinforced Plastics and Electr ical Insulat ing Mater ia ls by Four-Point Bend ing ( D 6272 - 00) was used as a guidel ine. 3.2.4.1 Specimens The bend specimens were cut from the same panel o f material used in the O C T tests. Six specimens o f three di f ferent depths were cut in each o f the 0° and 90° directions o f the mater ial . Test ing details for each specimen are g iven in Table 3.3. The tests were al l conducted in the plane o f the laminate. 3.2.4.2 Test set-up and procedure A schematic o f the test set-up is shown i n Figure 3.7. A test j i g was designed and bu i l t i n -house f o l l o w i n g the guidelines out l ined in A S T M D 6272 - 00. The t w o lower noses are r ig id l y attached to the crosshead wh i le the t w o upper noses are r ig id l y connected but a l lowed to swive l about the connect ion to the crosshead. This set-up al lows the j i g to sel f-al ign as i t loads. The diameter o f the loading noses is 5 m m , w h i c h is smaller than the recommended diameter g iven in the standard. For the 90° tests, al l o f the specimens fa i led in the region between the loading noses indicat ing tensile fai lure. V isua l inspect ion o f the loading points revealed no damage due to compression under the pins so the diameter o f the loading noses was not a factor in the fai lure. For al l o f the 0° specimens however, compressive fai lure occurred first under one o f the inner loading noses. Here the nose diameter was probably too small and va l id results were not obtained. A n extensometer ( Instron 2620-825 w i t h +5 m m travel) was fixed between a point r i g id l y attached to the base o f the test j i g and a poin t r i g id l y attached to the axis o f the swive l j o i n t 32 Chapter 3 Experiments on the top section o f the test j i g to measure displacement. The lower crosshead was moved up to load the specimen. The displacement rate depended on the specimen dimensions and was determined using guidel ines f r o m A S T M D6272 - 00. See Table 3.3. Photographs were also taken at regular intervals dur ing each test. However , fa i lure was always catastrophic and l i t t le useful in fo rmat ion was obtained f r o m the photographs. 3.2.4.3 Post-processing Load and displacement data were output at a frequency o f 10 s"1 and impor ted to M ic roso f t Excel to create plots. The stress i n the outer fibres o f the beam at the point o f cr i t ical load was calculated using the linear, smal l -def lect ion equation for stress g iven by the A S T M standard ( A S T M , 2000) as ^ Peril Ls ~ ! acnt  = V • 3.1 BD2 Results for the four-point bend tests are presented in Section 3.3.2. 3.3 Results and data reduction 3.3.1 OCT tests Table 3.4 summarizes the analysis that was conducted on each specimen. I n each specimen, a damage zone propagated from the notch t ip rough ly paral lel to the in i t ia l notch. For al l tests, the crack grew paral lel to the material 0° direct ion and there was negl ig ib le damage at the loading pins. Surface damage at the end o f a test is shown for a typ ica l specimen in Figure 3.8. 3.3.1.1 Load - Displacement behaviour F r o m the results o f notched tensile testing rev iewed in the l i terature, it was expected that the peak load attained by each specimen before fai lure w o u l d increase w i t h notch root radius. However , very l i t t le notch sensit iv i ty was observed for specimens containing notch root radi i ranging from 0.5 m m al l the w a y to 12.7 m m . Plots o f load versus C M O D are presented in Append ix A separately for each di f ferent notch radius tested. The results f r o m six representative specimens are also plot ted together in Figure 3.9 so that important features o f the responses can be compared. I n al l o f the experimental load versus displacement plots in 33 Chapter 3 Experiments this thesis sol id lines are used to show continuous port ions o f the curve and dashed segments are used to denote segments in w h i c h the crack has j u m p e d or load has dropped sharply between sampl ing points. I n al l cases the in i t ia l behaviour is l inear elastic. I n this phase, al l o f the strains in the laminate are reversible and no damage has occurred. A n increase in the compl iance o f the specimens is not iced as the radius o f the notch root is increased. This is an effect o f us ing re lat ively small specimens w i t h large notch root radi i . A t some point before the peak load is reached the behaviour becomes s l ight ly non-l inear. I n some cases, usual ly the smaller notch root radius specimens, (e.g. rd2-2, r d 8 - l and rd lO-1) sl ight shifts a long the C M O D axis are not iced indicat ing very smal l amounts o f damage are beginn ing to occur. I n the larger notch root radius specimens (e.g. rd32-2 and rd54-2) the sl ight shifts are not noticeable but the curve softens s l ight ly f r o m a straight l ine. The non-l inearit ies begin to be noticeable at around 11 k N o f load for al l specimens. The most dist inct feature di f ferent iat ing the responses o f the small and large notch root radius cases is the behaviour around peak load. I n Figure 3.9, specimens rd2-2, r d 8 - l , rd lO-1 and r d 2 5 - l are al l p lot ted using th in lines to dist inguish their behaviour from that o f specimens rd32-2 and rd54-2. The onset o f damage in the smaller notch radius specimens occurs at rough ly the peak load but does not result in an immediate load drop. Instead the peak load is maintained for a short te rm wh i le C M O D increases, result ing in a plateau. For the t w o larger notch root radius specimens, p lot ted here using th ick grey l ines, the behaviour is di f ferent. There is no region o f stable damage growth . Instead the peak load occurs r ight at the onset o f damage fo l l owed immediate ly by a large drop in load. The smaller notch root radius specimens appear to exhib i t stable fracture wh i le the larger ones temporar i ly exhib i t unstable fracture. The large load drop in the large notch radius specimens is due to a j u m p in crack length as w i l l be shown Chapter 4. The O C T specimens are large enough to contain the j u m p in crack length. The sharp load drop in the large notch 3 In the literature, the maximum load attained by a specimen is usually called the failure load because the specimens tested are unstable and have no residual strength after the failure load is reached. In the tests done here, the geometry is stable and the specimens are able to sustain load after the peak on the load versus displacement plot. The "failure load" is therefore referred to as the "peak load" in this thesis. 34 Chapter 3 Experiments root radius specimens reduces the stored energy in the specimen to less than the fracture energy and rapid crack g rowth is arrested. A s seen in Figure 3.9, the effect o f the notch root radius becomes negl ig ib le as damage progresses. Upwards o f 3.5 m m o f C M O D , the load versus C M O D curves for al l specimens are approximately the same. I t is remarkable that as notch root radius is increased, the transi t ion between specimens exh ib i t ing stable damage g rowth at peak load and those exh ib i t ing unstable damage g rowth is not gradual, but rather appears to be dist inct ly def ined. A l l o f the specimens w i t h notch root radi i o f 14 m m or less always displayed a sizeable plateau at about the peak load, wh i l e al l o f the specimens w i t h notch root radi i o f 16 m m or greater had load versus C M O D curves that dropped sharply at the onset o f damage. I t appears that at around 16 m m there is a t ransi t ion in behaviour. This transit ion in peak load at about a notch root radius o f 16 m m is more clearly def ined in Figure 3.10 o f peak loads plot ted against the notch root radius p. The peak loads are not normal ised w i t h the thickness o f the specimen. I t is assumed that the tensile fracture load o f a cont inuous fibre re inforced composite laminate is predominant ly dependent on the breaking o f f ibres and not mat r ix cracking. A l l o f the specimens used in this study have the same number and orientat ion o f pl ies and any differences in thickness are due to the amount o f resin and the fibre-volume fract ion. A n alternative is to normal ise the peak loads w i t h respect to the net section area o f each specimen. The stress at the notch t ip o f a C T or O C T specimen is due to a combinat ion o f tension and bending (Tada, Paris and I r w i n , 1985). The fai lure strength o f a specimen is therefore dependent on the net section area being loaded. The specimens were cut to a nomina l w i d t h W o f 81 m m but actual ly var ied f r o m 80.3 m m to 82.5 m m . The nomina l stress in the specimen is expressed as fo l lows a N = a N Tension N Bending f m 6P\ a + W-a P J 3.2 + [W-af W-a 2P(2W + a) [W - of 35 Chapter 3 Experiments W h i l e peak loads are di rect ly related to the dimensions a and W, normal isat ion w i t h respect to (2W+a)l(W-a)2 d id not change any o f the trends in the peak load versus C M O D plo t and i t was decided to leave the p lo t unnormal ised for clar i ty. A s shown in Figure 3.10, peak load is rough ly constant at about 15.5 k N between notch radi i o f 3.0 m m and 12.7 m m . Between p equal to 14.15 m m and 16 m m a transi t ion occurs after w h i c h the peak loads range between about 16.5 and 18.5 k N . 3.3.1.1.1 Variability Greater var iab i l i ty is seen in the peak loads o f the large notch root radius tests as w o u l d be expected. The plateau exhib i ted b y specimens w i t h smaller p indicates the development o f a damage zone ahead o f the notch t ip before a discrete crack is formed. Th is damage zone development absorbs the effects o f any small f laws that may affect the fai lure load o f more "b r i t t l e " materials and geometries. For large p, s low development o f the damage zone does not occur and f laws are more l i ke ly to affect the peak load attained by the specimen. A more detailed discussion o f what causes the differences in behaviour between specimens contain ing notches o f small root radius and those w i t h notches o f large root radius is g iven i n Chapter 5 after al l o f the experimental and numerical results are presented. 3.3.1.2 The size and shape of the damage zone The damaged regions ahead o f the in i t ia l notch t ip o f three specimens were sectioned to obtain a clearer picture o f the processes tak ing place at the notch t ip dur ing fracture. A relat ively sharp notched specimen, r d 2 - l , and the bluntest notch t ip tested, rd54-2, were chosen along w i t h a mid-s ized notch radius specimen, rd8-2. Results o f damage height versus distance from the in i t ia l notch t ip along the crack plane are presented in Figure 3 .11. I t is evident that the damage height at the notch t ip is much larger in specimen rd54-2 when the crack init iates than in specimen r d 2 - l . Photographs o f the sectioned material jus t ahead o f the in i t ia l notch t ips are shown in Figure 3.12 and Figure 3.13. For the large p specimen, the damage is spread out over almost 20 m m and the actual crack locat ion is not we l l def ined. For the small p specimen, damage is concentrated in a very small region on ly 2 to 3 m m in height. I n this case, al l o f the fibre breakage occurs along a wel l -def ined l ine through the thickness o f the sample. 36 Chapter 3 Experiments The extent o f the damaged area in specimen rd54-2 qu ick ly decreases unt i l a stable damage height o f about 7.0 m m is reached 10 m m from the in i t ia l notch t ip. I n specimen r d 2 - l the opposite t rend is observed; the damage height increases from about 2.5 m m to about 6.0 m m by the t ime the crack has g rown about 10 m m from the in i t ia l notch t ip . A l t h o u g h specimen rd8-2 was not sectioned al l the w a y to the back end o f the specimen, the damage zone was seen to commence w i t h a height o f approximately 5.0 m m and progress to around 6.0 m m at 10 m m from the in i t ia l notch t ip. These trends indicate that when fracture init iates from a notch o f re lat ively smal l root radius, the damage zone grows to its stable characteristic height as the damage zone progresses. Therefore, for small notch root rad i i , such as in r d 2 - l , the stress concentrat ion is reduced by the format ion o f the damage zone. A reduct ion in stress concentrat ion to be low its cr i t ical value causes the crack to arrest. Conversely, for a specimen w i t h a re lat ively large notch root radius, the notch ef fect ively "sharpens" to the stable characteristic height o f the damage zone for the material system. The sharpening behaviour o f specimen rd54-2 helps to explain w h y the crack j u m p s unstably when the peak load is reached. A s the damage zone forms at the notch t ip , the stress concentrat ion grows because the notch root radius is ef fect ively becoming smaller. The smal l di f ference in the stable characteristic height displayed by specimens r d 2 - l and rd54-2 is w i t h i n the experimental error o f the method used for measuring the damage zone height. W h i l e the characteristic height determined is on ly an approximate value, the important aspect o f this data is the sharpening versus b lunt ing behaviour that is observed for the large and small notch root radius specimens respectively. I n both the small and large radius specimens the damage zone terminates approximately 36 m m from the in i t ia l notch t ip . This corresponds to a C M O D o f 4.4 m m for r d 2 - l and 4.75 m m for rd54-2. The dif ference is explained by the dif ference in compl iance o f the t w o specimens. Since a large amount o f material is cut out o f rd54-2 i t is more compl iant and has a larger C M O D . A s the damage zone tapers out to v i rg in mater ial , the height o f the damage zone decreases to zero. The distance over w h i c h this decrease takes place gives an estimate o f the w i d t h o f the damage zone. The height o f the damage zone begins to decrease at around 29 to 30 m m for 37 Chapter 3 Experiments both r d 2 - l and rd54-2. The w i d t h o f the damage zone is therefore estimated to be about 6 m m . The decrease in the height o f the damage zone between 30 m m and 36 m m shown in Figure 3.11 comes f r o m the visual height o f damage in the cross sections. I t becomes more d i f f i cu l t to measure the damage height as damage becomes l ighter. Therefore the tapering o f the damage height approaching v i rg in material may not be the actual shape o f the damage zone. I t is more l i ke ly that the height o f the damage zone remains closer to the characteristic damage height unt i l complete ly undamaged material is reached. Results f r o m dep ly ing by Kongshavn and Poursart ip (Kongshavn and Poursart ip, 1999) and M i t c h e l l (M i t che l l , 2002) suggest that this is the case. T o determine the damage zone w i d t h more rel iably, residual strength testing o f the sectioned material should be done. Neg l ig ib le strength w o u l d indicate a f u l l y separated crack wh i l e increasing residual strength w o u l d indicate decreasing amounts o f damage unt i l the strength o f the undamaged material is reached at the end o f the damage zone. Qual i ta t ive ly , i t was noted that the th in sections (less than 1 m m th ick) became more d i f f i cu l t to pu l l apart and contained some unbroken fibres at about 30 m m from the in i t ia l notch t ip. 3.3.1.3 Determining the critical fracture energy Gc One o f the pr inc ipa l benefits o f obtain ing the stable g rowth o f a crack from a stable specimen geometry is that i t a l lows the determinat ion o f the fracture energy o f the material s imp ly by measuring the area under the load-displacement curve. This method is not va l id fo r rapid, unstable crack g rowth because the energy release rate can be larger than the cr i t ical value Gc and the area under the curve w i l l include energy dissipated through v ibrat ion, noise, heat, etc. For stable crack g rowth , the amount o f released energy is equal to Gc. A n average value o f Gc is therefore calculated as the total w o r k done at the pins d iv ided by the area o f crack g rowth as This method assumes a constant crack resistance R, w h i c h is also assumed in the numer ical model . This is a good assumption for this material since i t appears that Gc stabilizes after crack in i t ia t ion. A s seen in Figure 3 .11 , a stable zone o f sel f -s imi lar crack g rowth was observed for both in i t ia l ly b lunt notches and sharp notches. 38 Chapter 3 Experiments For calculat ion o f the w o r k done, it is important to measure the p in opening displacement ( P O D ) and not C M O D , since the load is measured at the pins. Unfor tunate ly , P O D was on ly recorded in the later tests and is not available for the major i t y o f the smaller notch radius cases. Specimens rd25-2 and rd28-2 both displayed a plateau on the load versus P O D curves indicat ing stable development o f the damage zone and they are used here to determine Gc values. Crack lengths are determined f r o m the crack opening prof i les o f the damaged specimens. T o obtain the crack opening displacements as damage progressed in the O C T specimens, the l ine analysis technique described in Section 3.2.3.3 is used. The l ine analysis plots for specimens rd25-2 and rd28-2 are shown in Figure 3.16 and Figure 3.18 respectively. Due to surface damage, i t was not possible to use the l ine closest to the plane o f the notch. Typ ica l l y either the 3 r d or 4 t h l ines on either side o f the notch plane were used to make the measurements. The relat ive separation o f the lines therefore includes both the C O D and the straining o f 7.5 m m to 10.0 m m o f material on either side o f the crack plane leading to a small error in the measurements. A correct ion is made by measuring the relat ive separation in the scribed lines used in the analysis jus t p r io r to the first non- l inear i ty i n the l o a d - p i n displacement curve. Assuming that no damage has occurred at this point , the relative separation at the notch t ip is approximately the amount o f elastic strain in the mater ial . This strain is noted in the f igures as a dashed l ine and the C O D s determined f r o m the analysis are measured from the intersection o f the C O D prof i le w i t h the dashed l ine. L ine analysis was also conducted for specimens r d 2 - l , rd2-2, rd8-2 , r d 5 4 - l and rd54-2 and the plots are presented in Append ix B. The final crack lengths determined by l ine analysis are compared to the values obtained from sect ioning for specimens r d 2 - l and rd54-2 in Figure 3.20. I n both specimens, the l ine analysis value underestimates the total length o f damaged material obtained from sectioning. M i t che l l (M i t che l l , 2002) observed the opposite t rend but included the apparent crack length from elastic straining o f the material between the lines used in the analysis. W i t h o u t inc lud ing the apparent crack length from elastic strains, the l ine analysis gives an approx imat ion to the crack length fa l l ing somewhere in the damage zone. For Gc determinat ion, this approx imat ion is reasonable. I n Equat ion 3.3, a is 39 Chapter 3 Experiments an effect ive through-crack representing the real crack composed o f regions o f f u l l damage and part ial damage. Plots o f load versus P O D for specimens rd25-2 and rd28-2 are g iven in Figure 3.17 and Figure 3.19 respectively. The locations where l ine analysis was conducted are marked and the change in crack length is noted on the plots. The amount o f w o r k done by the system is calculated by assuming an elastic unloading path back through the or ig in . No te that the o r ig in has been shif ted s l ight ly to account for the in i t ia l take-up o f slack in the system. A lso , it can be seen f r o m the f igures that the actual un loading path does not return to the or ig in but shows about 1.75 m m o f residual displacement. I t is assumed that this residual is not due to plastic deformat ion, but rather to the damaged surfaces not being able to f i t back together proper ly w i t h i n the damage zone. The w o r k done and the calculated values o f Gc are tabulated versus change in crack length for both specimens in Table 3.5. Exper imental points on the crack resistance curve are plot ted in Figure 3 .21. There is a large scatter in the results due to the error in approx imat ing the crack length using l ine analysis. A t short lengths o f crack g rowth a small error i n Aa used to calculate Gc has a large effect on the solut ion. A t longer crack lengths the solut ion is less sensitive to small deviations in Aa so that a more rel iable, al though st i l l approximate, solut ion is obtained. Even w i t h the uncertainty in the results at small Aa, a trend o f in i t ia l ly increasing crack resistance is apparent. However , this is not true i?-curve behaviour. The lower Gc values determined at small Aa correspond to the development o f the damage zone. A f te r the damage zone height has stabil ised at its characteristic height (as shown in Figure 3.11), Gc appears to stabilise as w e l l . Us ing the two largest Aa points on the load versus P O D plot o f each specimen to calculate Gc, a value o f about 80 to 85 kJ /m is determined. 3.3.2 F o u r - p o i n t b e n d tests Four-point bend tests were conducted to determine the unnotched tensile strength o f the mater ial . Results f r o m the four-point bend tests are presented in Table 3.6. The load versus displacement p lot for specimen b85-90-2 is g iven in Figure 3.22 as a representative case. The rest o f the plots for the 90° or ientat ion specimens are presented in Append ix C. 40 Chapter 3 Experiments Separate tests were done w i t h the specimens oriented in both in-plane directions. I n the f i rst case, the laminate 90° di rect ion was oriented along the axis o f the beam w i t h the stronger 0° d i rect ion oriented in the transverse direct ion. This or ientat ion corresponds to the or ientat ion o f the O C T specimens and is used to determine the y-d i rect ion strength for the numer ical model o f the O C T specimens. Tests were also conducted w i t h the 0° f ibres oriented a long the axis o f the beam to determine strength in the x-d i rect ion. However , va l id results were on ly obtained for the 90° or ientat ion, w h i c h is the important d i rect ion for numer ical mode l input. I n al l o f the 90° or ientat ion tests fai lure occurred in the centre o f the span due to tension in the outer f ibres. For the specimens oriented in the 0° d i rect ion, fai lure occurred in compression under the loading pins and therefore gave incorrect tensile strength values. The strength values for the 0° or ientat ion specimens are on ly lower bounds on the tensile strength. A s shown by the sharp load drop in Figure 3.22, fai lure was catastrophic, resul t ing in a w e l l -def ined peak load. The stiffness o f the specimen, however, becomes non- l inear at a point before the peak load indicat ing that a small amount o f damage has taken place before fai lure. The dashed l ine in the f igure il lustrates that the curve first becomes non- l inear at around 1 k N o f load. Equat ion 3.1 is used to calculate the stress in the outer fibres at peak load. A s stated in the specif icat ion Equat ion 3.1 is on ly str ict ly va l id for materials that are l inear up to the point o f rupture so a sl ight error w i l l be introduced in using this equation. A n approx imat ion to the tensile fai lure strength o f the material is then taken as the average o f the tests. For the 90° or ientat ion the f lexural strength is taken f r o m the average o f five tests to be approximately 460 M P a . For the 0° or ientat ion the m i n i m u m value o f the f lexura l strength is taken as the m a x i m u m f lexural strength recorded f r o m five tests o f about 1000 M P a . 3.4 S u m m a r y • A transit ion in the load-displacement behaviour o f the O C T specimens was noted at a notch root radius o f about 16 m m . O C T specimens w i t h radi i less than this value displayed stable crack g rowth w i t h l i t t le notch sensit iv i ty. O C T specimens w i t h larger radi i were unstable and showed more notch sensit iv i ty. 41 Chapter 3 Experiments • Post-test sectioning o f the O C T specimens reveals that the damaged height o f mater ial is in i t ia l ly dependent on the radius o f the notch root. A large in i t ia l radius creates a ta l l , but dif fuse damage zone wh i le a sharp notch concentrates the damage in a short (height) damage zone. A f te r stable crack g rowth is achieved, a characteristic damage height o f 6 m m to 7 m m is produced for al l in i t ia l notch sizes. • A cr i t ical fracture energy value o f 80 to 85 kJ /m and a strength o f 460 M P a was determined for the material used in this thesis. 42 Chapter 3 Experiments Table 3.1 Original elastic material properties for CFRP laminate (Mitchell, 2002) P r o p e r t y V a l u e Ex 75 GPa Ey 30 GPa Oxy 0.161 GXy 17.1 GPa Chapter 3 Experiments Table 3.2 OCT tests and specimen measurements. Specimen Notch root at, (mm) ^(mm) B (mm) Identifier radius (mm) r d l - 1 0.5 33.71 81.84 8.98 r d 2 - l 1 33.3 81.71 8.49 rd2-2 1 33.35 81.92 8.93 r d 4 - l 2 33.42 82.52 9.04 r d 4 - 2 . 2 33.43 80.53 8.32 r d 6 - l 3 33.51 82.17 8.55 rd6-2 3 33.07 80.90 8.71 rd6-3 3 48.1 81.90 8.75 r d 8 - l 4 33.35 81.91 8.89 rd8-2 4 33.45 82.10 8.85 r d l 0-1 5 33.03 80.70 8.55 r d l 0-2 5 33.2 81.75 8.57 r d l 7-1 8.5 33.33 81.00 8.43 r d 2 1 - l 10.5 33.43 81.20 8.13 r d 2 5 - l 12.7 33.99 81.91 8.07 rd25-2 12.7 34.1 80.50 8.00 r d 2 8 - l 14.15 34 80.50 8.10 rd28-2 14.15 33.8 80.30 8.50 r d 3 2 - l 15.9 33.71 81.93 8.70 rd32-2 15.8 34.2 81.20 8.60 rd32-3 15.8 34.18 81.32 7.91 r d 3 5 - l 17.5 33.8 81.24 7.96 rd35-2 17.5 34.03 81.41 8.30 r d 3 8 - l 19 33.64 81.96 8.26 rd38-2 19 33.76 81.42 8.35 r d 4 4 - l 22.25 33.46 82.02 8.56 rd44-2 22 33.92 81.06 8.32 r d 5 4 - l 27 32.85 81.3 8.83 rd54-2 27 33.16 81.12 8.14 44 Chapter 3 Experiments Table 3.3 Four-point bend tests and specimen measurements Specimen I d e n t i f i e r O r i e n t -a t ion D ( m m ) B ( m m ) L s ( m m ) L L ( m m ) L s / D L o a d i n g Rate ( m m / m i n ) b56-90-l 90 5.92 7.79 135 45 22.8 5.70 b56-90-2 90 5.81 7.81 90 30 15.5 2.58 b85-90-l 90 8.70 8.21 135 45 15.5 3.88 b85-90-2 90 8.65 8.18 135 45 15.6 3.90 bl50-90-l 90 15.20 8.18 135 45 8.9 2.22 bl50-90-2 90 15.23 7.8 135 45 8.9 2.21 b56-0-l 0 6.00 8.05 90 30 15.0 2.50 b56-0-2 0 6.05 8.03 90 30 14.9 2.48 b85-0-l 0 8.64 8.07 135 45 15.6 3.90 b85-0-2 0 8.67 8.03 135 45 15.6 3.89 bl50-0-l 0 15.18 8.13 135 45 8.9 2.22 bl50-0-2 0 15.14 8.06 135 45 8.9 2.23 45 Chapter 3 Experiments Table 3.4 Summary of OCT tests and analysis Specimen L o a d - L o a d - Per iod ic Cyc l i c L ine Gc calc- Sections I d e n t i f i e r C M O D Pin Un loads U n l o a d i n g Ana lys is u la t ion Disp. r d l - 1 r d 2 - l rd2-2 r d 4 - l rd4-2 >/ r d 6 - l rd6-2 V • rd6-3 r d 8 - l rd8-2 •/ rd lO-1 r d l 0-2 V r d l 7-1 V r d 2 1 - l r d 2 5 - l rd25-2 r d 2 8 - l rd28-2 r d 3 2 - l s rd32-2 •/ rd32-3 s r d 3 5 - l rd35-2 r d 3 8 - l V rd38-2 Y r d 4 4 - l Y rd44-2 S S r d 5 4 - l rd54-2 46 Chapter 3 Experiments Table 3.5 Critical fracture energy release rates rd25-2 Aa (mm) W ( J ) Gc (kJ/m2 ) 8.4 3.2 48.0 25.0 13.0 64.9 27.4 19.0 86.8 32.4 23.3 90.0 rd28-2 Aa (mm) W{A) Gc (kJ/m 2 ) 12.5 5.0 47.4 14.5 11.6 94.0 26.5 19.0 84.2 33.5 23.8 83.6 Chapter 3 Experiments Table 3.6 Four-point bend test results Specimen M a x i m u m Fa i l u re load M a x . stress d isp lacement ( k N ) at f a i l u r e ( m m ) ( M P a ) b56 -90 - l N A N A N A b56-90-2 4.981 1403.9 479.3 b85 -90 - l 7.237 2190.9 476.0 b85-90-2 6.530 1989.0 438.7 b l 5 0 - 9 0 - l 3.819 6394.3 456.8 b 150-90-2 3.881 6089.8 454.4 D56-0-1 3.962 3263.4 1013.5* b56-0-2 4.087 3295.6 1009 .1* b 8 5 - 0 - l 5.395 4085.7 915.6* b85-0-2 5.855 3644.8 815.2* b150-0-1 N A N A N A B150-0-2 3.467 10644.5 777.8* * Strength values for the 0° specimens are the m i n i m u m bound o f strength because fai lure occurred in compression under the loading pins instead o f tension between the pins. 48 Chapter 3 Experiments ^ 106 m m ^ c — 38.7 m m W= 80.6 m m a = 33 m m Figure 3.1 OCT specimen geometry (thickness B ranged from 7.9 mm to 8.6 mm) 49 Chapter 3 Experiments Chapter 3 Experiments Figure 3.3 Scribed lines for line analysis (shown for rd4-2) Figure 3.4 OCT test set-up (shown for rd54-2) 5 1 Chapter 3 Experiments 0 0.5 1 1.5 2 2.5 3 3.5 4 Displacement (mm) Figure 3.5 Load-displacement plots for specimen rd32-2 from the machine position and from the relative displacement of the pins. Figure 3.6 Schematic of section taken through the damage zone ahead of the initial notch tip. The dashed lines represent the 90 °fibres in the outer stacks. The measured damage height did not include any damage to the outside of these lines. 52 Chapter 3 Experiments Spec length B = 8.5 m m . ca cC± L L u Figure 3.7 Four-point bend test geometry Figure 3.8 Specimen rdl 0-2 showing orientation of crack and extent of surface damage (scribed lines are 2.5 mm apart) 53 Chapter 3 Experiments 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 C M O D (mm) Figure 3.9 Load versus crack mouth opening displacement for specimens of representative notch root radius. Solid lines indicate continuous portions of the curve while dashed segments indicate segments in which the crack jumps and the load drops sharply. O o O O o • • t * p (mm) 20 25 30 Figure 3.10 Peak loads attained in each OCT test versus notch root radius p. Closed data points denote specimens displaying a plateau on the load versus CMOD curve, open data points denote specimens that showed sharp load drops at the peak load. 54 Chapter 3 Experiments 20 18 16 14 ? 12 E 10 •a 0 Figure 3.13 rd54-2 Figure 3.15 O | sections w o • r X Figure 3.12 6 4 2 0 rd8-2 Damage zone development 10 Figure 3.14 Stable damage zone 15 20 25 x (mm) Partial damage 30 Figure 3.11 Profiles of the damage zone along the crack plane 35 Figure 3.12 Section of rd2-lat 0.635 mm ahead of initial notch tip 55 Chapter 3 Experiments Figure 3.13 Section of rd54-2 at 0.635 mm ahead of initial notch tip Figure 3.14 Section of rd2-l at 24.13 mm ahead of initial notch tip 5 6 Chapter 3 Experiments Figure 3.15 Section of rd54-2 at 24.13 mm ahead of initial notch tip 1.8 Distance from init ial notch tip (mm) Figure 3.16 COD profiles determinedfrom line analysis for rd25-2 57 Chapter 3 Experiments Aa =0 mm POD (mm) Figure 3.17 Load versus pin opening displacement for rd25-2 showing location ofphotos for line analysis — - image38 POD=l .7 mm -«-image42 P0D=2.1 mm - • - image44 POD=2.5 mm - * - image48 POD=3.2 mm image51 POD=3.8 mm -2 0 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 Distance from initial notch tip (mm) Figure 3.18 COD profiles determinedfrom line analysis for rd28-2 58 Chapter 3 Experiments POD (mm) Figure 3.19 Load versus pin opening displacement for rd28-2 showing location ofphotos for line analysis • Sectioning (full damage height) • Sectioning (partial damage height) 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 Crack length (mm) Figure 3.20 Final crack lengths in rd2-l and rd54-2 from sectioning and line analysis 59 Chapter 3 Experiments 100 90 80 70 S S 60 C | 50 40 30 20 10 0 T 3 « O 2.5 2.0 1.5 1.0 0.5 0.0 -0.5 10 15 20 A a (mm) 25 • rd25-2 • rd28-2 30 Figure 3.21 Crack resistance points for rd25-2 and rd28-2 3 4 P i n displacement (mm) 35 Figure 3.22 Four-point bend test load-displacement plot for b85-90-2 60 Chapter 4 Numerical C H A P T E R 4 N U M E R I C A L 4.1 Introduction A numerical study o f damage development and propagat ion in FRP composite laminates is undertaken in this thesis. I n part icular, the over-height compact tension tests described in Chapter 3 are simulated us ing a three-parameter damage model . There are several goals to this endeavour: 1) T o c lar i fy the state o f stress and strain at the notch t ip at in i t ia t ion o f fracture and dur ing propagat ion using f in i te element numer ical models. 2) T o develop a phys ica l ly based cohesive zone model that is described by independently determined material parameters capable o f match ing experimental results. 3) T o investigate the abi l i ty o f the cohesive zone model to predict the transit ion in structural behaviour w i t h notch root radius observed in the experimental O C T tests. 4) T o investigate the dependence o f the model response on f in i te element mesh size and conf igurat ion and the interact ion o f mesh size and material parameters. Speci f ical ly , the abi l i ty to accurately model crack propagat ion on a large scale using a coarse mesh is explored. The first goal can be par t ly achieved using elastic simulat ions o f the O C T tests at the po in t o f peak load. Here, the size and shape o f the stress fields give clues to the response o f the various sized notches in the O C T tests. T o achieve the remain ing three goals, t w o related user material models ( U M A T s ) have been developed and implemented in the A B A Q U S f in i te element computer program. The first model , cal led the simple damage model ( S D M ) , is a cohesive zone model requi r ing ( for each material direct ion) the in i t ia l elastic modulus Eo, the strain at the onset o f damage epeak and the cr i t ical specif ic fracture energy yc as input parameters. The second model is a mod i f i ca t ion o f the S D M . I t includes the abi l i ty to adapt automat ical ly the material const i tut ive relat ionship for each element based on mesh size and the strain field. 61 Chapter 4 Numerical The development o f these models contr ibuted signi f icant insight into the observed response o f the experimental O C T tests. 4.2 Elastic simulations The adaptive damage model presented in this chapter makes use o f the strain d is t r ibut ion w i t h i n an element. Before developing the model , elastic simulat ions o f the O C T specimens were carried out using the A B A Q U S f in i te element computer program to gain a better understanding o f the stress and strain distr ibut ions at the notch t ip . A n estimate o f the element size required to resolve the stress or strain f ie ld is also made. 4.2.1 Stress concentrations Elastic stress distr ibut ions along the notch plane ahead o f the notch t ip at the poin t o f peak load were obtained from static f in i te element analyses o f the O C T specimens using the orthotropic material properties g iven i n Table 4 . 1 . Note that the modulus E2 was increased s l ight ly from 30 GPa, as g iven by the material supplier, to 32 GPa to better match the compl iance o f the experimental specimens. O n l y the bo t tom h a l f o f the specimens were modeled, w i t h symmetry condit ions appl ied along the notch plane. Two-d imens iona l , f u l l y -integrated, 8-noded, plane stress elements were used w i t h 0.5 m m x 0.5 m m elements i n the region o f the notch plane and stresses were quadrat ical ly extrapolated to the nodes. The computed stress fields for specimens w i t h notch radi i o f 1, 4, 10.5, 16 and 27 m m are presented in Figure 4.1 for a load o f 15.0 k N . Th is load level is approximately the peak load found exper imental ly for specimens displaying a stable plateau on the load versus C M O D plots (i.e. p = 0.5 to 14 m m ) . The linear-elastic material model used does not a l low any damage to occur at the notch t ip . I t was seen from the experiments that damage does occur before the peak load is reached, as evidenced by sl ight non- l inear i ty in the l o a d - C M O D curves. The elastic stress f ields are therefore not complete ly accurate, however they do a l low a reasonable approx imat ion o f the stress states for compar ison between notch sizes. The material tensile strength o f 460 M P a that was determined f r o m the four-po in t bend test and is used later in the damage model is marked in Figure 4.1 w i t h a dashed l ine. Clear ly, a notch root radius o f 27 m m creates a stress concentrat ion at the notch t ip that is less than the tensile strength when the load is at a cr i t ical level for smaller notch sizes. A t a notch root radius o f 16 m m , w h i c h is rough ly the transi t ion point between the stable and unstable fracture 62 Chapter 4 Numerical initiation seen experimentally, the stress is approximately the same as the tensile strength. For the smaller notch root radii the stresses at the notch tip are well above the tensile strength of the material meaning that the stress criterion for fracture has been met. Whether or not fracture initiates depends on whether sufficient energy is stored in the system. Increasing the load to the experimentally observed peak load of 17.8 kN for the 27 mm notch root radius specimen magnifies the elastic stress concentration at the notch tip to slightly above the tensile strength as shown in Figure 4.2. 4.2.2 Stress contours From the elastic simulations it is also possible to estimate the extent of critically stressed material for a given load. The contour plots in Figure 4.3 show the elastic stress fields in the region of the notch tip. Each notch size is shown at approximately its peak load, as determined from the experimental tests. The light grey areas around the notch tips are areas that experience a stress over 460 MPa at the peak load for that notch root radius. Generally, there is an elliptical shape to the stress contours with the critically stressed area extending by about one to two millimetres ahead of the notch tip, with a gradual decrease in the width with increasing notch root radius. Conversely, it can be seen that as the notch root radius is increased, the height of critically stressed material around the notch tip increases. Table 4.2 lists the approximate height of critically stressed material at the experimentally determined peak load for each of the specimens from the elastic analyses. Although the trend is increasing height with increasing notch root radius, the critically stressed height dips to 8.5 mm for p equal to 27 mm from 9.5 mm for p equal to 16 mm suggesting that perhaps experimental specimens rd54-l and rd54-2 failed prematurely. The significance of the critically stressed height is discussed further in Chapter 5. 4.2.3 Element size required to resolve the notch tip stress concentration It was found that relatively large, fully-integrated, 8-noded elements that have quadratic shape functions do a much better job of resolving the stress concentration ahead of a notch tip than smaller, constant stress elements containing a single integration point. Typically, explicit dynamic finite element codes that are used to model dynamic crack growth require 63 Chapter 4 Numerical the use of reduced integration, constant or linear stress elements for efficiency. Static, implicit codes on the other hand are often more efficient using higher order elements. The stress concentration at the tip of an OCT specimen with a four-millimetre radius notch is modeled in Figure 4.4 using four different width constant stress elements ranging from 0.25 mm to 2.0 mm and also using 0.25 mm and 0.5 mm wide quadratic elements. All of the elements match the stress field away from the notch tip well. The notch tip extrapolations of the stress are not very good for the constant stress elements but they do approximate the true solution fairly well at the location of the integration points (e.g. at 0.25 mm along the notch plane for the 0.5 mm wide constant stress element). The higher order elements do a better job of predicting the notch tip stress. There is little difference between the 0.25 mm and 0.5 mm solutions. 4.3 Development of the cohesive zone models 4.3.1 Simple damage model (SDM) A three-parameter cohesive zone model for fracture of composite laminate materials was developed as a UMAT in ABAQUS. The model is intended as a simple representation of a more advanced composite damage model that has been developed by the Composites Group at the University of British Columbia called CODAM (Williams, Vaziri and Poursartip, 2003). CODAM (Composite DAMage Model) uses physical input parameters describing the onset and saturation of damage for fibre and matrix in a composite laminate. A linear relationship between a strain function that incorporates strain interactions from other material directions and a monotonically increasing damage parameter is defined for both fibre and matrix in each direction. Total damage is defined as a combination of fibre and matrix damage depending on fibre and matrix properties and lay-up sequence. Material moduli are then reduced according to a linear relationship with total damage, resulting in a strain-softening stress-strain curve similar to the curves used in cohesive zone models. The model developed here loses some of the generality of CODAM but is still physically-based and is used to demonstrate characteristics of softening models in a simple manner. The model is therefore termed the SDM (Simple Damage Model). It was decided to implement the model in ABAQUS to take advantage of the features of the implicit code. An implicit solution scheme uses an iterative technique such as the Newton-Raphson method to solve non-linear problems and is ideal for structural problems involving 64 Chapter 4 Numerical long time durations since the step size can be relatively large. Explicit, direct integration methods on the other hand require much less computation per time step but require relatively small time steps, making it ideal for wave propagation problems but not long duration structural dynamics problems. The unstable growth of a crack in a continuum would be best modeled using an explicit integration scheme since the crack travels at close to the wave speed of the material and inertia is important in the solution. However, the OCT tests that are modeled in this thesis are quasi-static tests and display, for the most part, stable crack growth. Therefore, an implicit code, which does not include inertia, can give an accurate solution much more efficiently than an explicit dynamic code. In addition, although COD AM is implemented in the LS-DYNA explicit finite element code, higher-order elements are not available in the code.4 Higher-order elements containing a quadratic strain field (available in ABAQUS) are required in the modified SDM, described below, in order to achieve a good description of the strain field in a single element. 4.3.2 Adaptive simple damage model (ASDM) A modified version of the SDM will be described along with the development of the SDM. The modified version attempts to take into account the effect of element width5 that was noted in the literature review. A coarse mesh in the width direction introduces error as demonstrated schematically in Figure 4.5. For a notch-tip element in which stress is calculated at the centre of the element, the stress at the calculation point is less than critical when the stress at the edge of the element has reached the critical level. This is due to the stress distribution at the notch tip. The wider the element, the lower the calculated stress. With the typically used Gauss integration scheme or a constant stress element, an element's calculation points are located at some point within the element and not on the edges. This means that the stress field at a notch tip will not be calculated at the tip, but rather a certain distance ahead of it. That distance depends on the size of the element and the integration scheme that is used. In higher 4 The reason for this is that higher-order elements yield higher maximum frequencies than lower-order elements. In explicit time integration schemes the required step size to achieve a correct solution decreases with the maximum frequency since the step size must be small enough so that information does not propagate across more than one element per time step. 5 Width refers to the element dimension in the direction of crack growth. 65 Chapter 4 Numerical order elements error is reduced because the integration points are located closer to the element edge. The effect of sampling the stress field away from the notch tip is that damage initiation in the element is delayed. When the stress at the notch tip/ element edge is at the critical level for crack initiation, the stress at the calculation point will be lower than critical. When the stress at the calculation point reaches the critical level, the stress at the notch tip may be considerably higher than critical. For a large element there may be enough energy G in the system to propagate the crack, but it cannot initiate until the stress criterion is also met. In that case the applied load must be increased until the stress seen by the element reaches the material strength. In doing so, the stress at the notch tip will increase to above the material strength and G will surpass Gc. With the sudden initiation of the crack, the additional stored elastic energy is released into the crack causing it to propagate unstably until G <GC. The difference between the actual stress at the notch tip when damage initiates and the critical stress depends on the stress concentration as well as the location of the integration point. From Figure 4.1, it can be seen that, using a finite sized element, a much larger error would be incurred in modeling an OCT specimen with a notch tip radius of 1 mm than a notch tip with a radius of 27 mm. The much steeper stress distribution in the sharper notch case means that the difference between the actual notch tip stress and the calculated stress at the integration point is very large. As a notch tip approaches the infinitely sharp case smaller and smaller elements should be required in order to minimize the difference between the actual and calculated stresses. Using the existing composite damage model (CODAM) the mesh width effect was seen in modeling a unidirectional carbon fibre/ epoxy double-cantilever beam (DCB) specimen. The input parameters in the direction perpendicular to the notch (out-of-plane) were determined to match those of a cohesive zone model with a three-parameter softening constitutive curve (initial modulus, strain at onset of damage and an ultimate strain at complete separation of the interface). The region ahead of the notch tip was discretized with four different constant stress element widths; 62.5 pm, 125 pm, 250 pm and 500 pm. The height of the elements was maintained at 250 pm. The arms of the virtual specimens were opened at a quasi-static loading rate past crack initiation. The effect of element width can be seen in Figure 4.6. The 66 Chapter 4 Numerical peak load was strongly dependent on the width of the element with the apparent strength of the specimen increasing with element size. Under the same loading conditions and for the same material properties, peak load varied from 19.5 N for 62.5 pm elements to 29.3 N for 500 pm elements. There are two ways of dealing with this situation. The first is to refine the mesh at the tip of a notch to very small elements to reduce the error. A problem with this method is immediately apparent. As mentioned in Chapter 2, mesh refinement when using a strain-softening material model does not lead to a converged solution due to localisation and the dependence of GF on the height of the localised band of elements. Another option is to adapt the constitutive material model to take into account the distance of the integration point from the notch tip. In this manner, damage would initiate in an element when the stress at the element edge reaches the material strength. The stress-strain curve at the integration point would then need to be adjusted so that the critical stress level is lower while still maintaining the area beneath the curve for matching GC. The second option is the solution technique that has been incorporated into the SDM. Since it involves adapting the constitutive model to the stress concentration and element size, it is called the Adaptive Simple Damage Model (ASDM). Returning to the problem of the DCB specimen modeled with different element widths, the method is demonstrated. When the problem was first modeled with a series of simulations using elements of different width, different global peak loads were obtained (light grey points in Figure 4.7). Scaling the material strength of the larger sized elements against the 62.5 pm elements and re-running the solutions resulted in equivalent peak loads being predicted (black points in Figure 4.7). To determine the scaling ratio for each element size, an elastic simulation of the DCB was done with 62.5 pm elements in the region of the crack plane to obtain the stress distribution ahead of the crack tip as shown in Figure 4.8. The constitutive curves of the 125 pm, 250 pm and 500 pm element width specimens were then scaled by the ratio of stress at the element midpoints to the stress at the midpoint of the 62.5 pm element. In the ASDM, the procedure described above for the DCB specimen is done automatically in the material model. The strain distribution is determined across a second order element and 67 Chapter 4 Numerical the strength at the integration point is scaled by the ratio of the strain at the integration point to the strain at the element edge. 4.3.3 Algorithm for the SDM and ASDM The SDM and A S D M are implemented as UMATs in A B A Q U S for use with plane stress two-dimensional solid elements. The SDM can be used with elements of any order or integration scheme but the A S D M requires fully-integrated, 8-noded elements to function properly. Similar to the more advanced composite damage model C O D A M , the S D M input parameters are derived from the properties of fibre and matrix smeared together to form the orthotropic material properties of each layer of a laminate. These layers are called sub-laminates. In the SDM, the material behaviour in each principal direction for a sub-laminate is defined by three parameters: initial modulus EQ, elemental specific fracture energy ye and characteristic critical strain Sopeak (strain at damage initiation). These values are all invariant with element size. The peak stress apeak (strength) is determined from EQ and £opeak and the ultimate strain Suit (strain at damage saturation) follows from <7peak and ye. ^peak The softening behaviour described by these parameters is shown in Figure 4.9. Damage initiates when s is equal to speak and a linear, negative tangent modulus reduces stress with increasing strain until complete damage is reached at suit. Softening only occurs in tension. In compression the material is assumed to be undamaged (linear elastic). Currently there is no strain interaction between material directions so strains in a given direction only cause damage in the same given direction. At this stage the model is suitable for modeling Mode I crack growth. Further development is necessary to define the constitutive relationships for shear loading as would be necessary to model mode II or mode III crack growth. Before the onset of damage, the model is identical to an orthotropic elastic material model. The constitutive curve is defined separately in each material direction. The model is not limited to modeling only interfaces, but can be used as the material model for an entire structure, E()£peak 2ye peak 4.1 4.2 68 Chapter 4 Numerical allowing the location of stress concentrations to determine the location of crack growth. As an elastic damage model there are no plastic strains before complete damage occurs and unloading takes place on a line through the origin. In this material model the criterion for crack initiation is strain equal to the damage initiation strain epeak and the criterion for complete separation of the crack is the dissipated specific fracture energy equal to the critical specific fracture energy yc. Since damage is an irreversible process, the amount of fracture energy dissipated yd is calculated at each step and the residual fracture energy yr, equal to yc - yd, is stored as a history variable by the code and used to keep track of the current state of the material. The current state is uniquely defined by yr, the monotonically increasing current peak strain6 epeak and the residual secant modulus Er as shown in Figure 4.9. Initially speak is set equal to the strain at onset of damage eopeak and the modulus and fracture energy are set equal to their undamaged values. £ — F peak 0 peak Er=E0 4.3 Yd =7e The algorithm for the SDM and A S D M is shown as a flowchart in Figure 4.10. Both models follow the same basic structure with the A S D M including some additional operations. A double box around a process indicates a calculation or subroutine that is only included in the A S D M . The A S D M routine is briefly described here and the U M A T code is included in Appendix D. The following steps are carried out at each time or load increment and iterated until equilibrium between the internal and external loads is achieved (i.e. convergence). Loop over each element: Step 1. Determine the strain profile in each element The strain profile at a notch tip or damage zone is used to carry out parameter scaling as opposed to the stress profile because, for a softening material model, the strain field is always increasing towards a notch tip while the stress field will begin to decrease towards the notch tip once softening has initiated. For an initially elastic material, as long as the crack is 6 This strain is referred to as the peak strain since it coincides with the peak stress on the constitutive curve. 69 Chapter 4 Numerical small compared to the other planar dimensions, the surrounding elastic material controls the strain in the damaged material at the notch tip. In a fully integrated, 8-noded, plane stress element there are nine integration points. Using a Gauss scheme the integration points are located as shown in Figure 4.11. The x-direction and y-direction strain fields across each row and column of integration points are determined by fitting a second order polynomial to the strains at each point. This is similar to using the element shape functions but allows the x - and y-strain fields to be described separately for each row or column using three polynomial coefficients. For the first time or load step the strains are zero everywhere and subsequently the strains are taken from a global storage matrix, which holds the integration point strains from the previous time step. Using a second order polynomial allows the strain function to closely approximate the r'l/2 distribution found near a sharp notch. In the current methodology, there is no provision to ensure continuity in the strain function at the edges of the element. This can lead to small errors in the magnitude of the strain field at the element edges. To fix the problem the strains in neighbouring elements would need to be taken into consideration, greatly increasing the complexity of the problem. Loop over each integration point within the element: Step 2. Determine the scaling strain First, the integration point numbering sequence within the element is determined with respect to the local material axes to orient the element. There are six strain fields in the x-direction and six strain fields in the y-direction represented by separate polynomials in each element corresponding to the dashed segments shown in Figure 4.11. Each integration point is at the intersection of two of these segments, one segment in the x-direction and one segment in the y-direction. The maximum normal strains Smax, on each of the two segments extrapolated to the element edges, corresponding to the integration point are determined and become the x - and y-scaling strains for that integration point. That is, the x-strain on the y-direction segment and the y-strain on the x-direction segment. If the element is within the strain field caused by a notch or other stress concentrator then the maximum strains usually fall on the element edge. The element edge that has the highest strain is in the direction of the notch tip. 70 Chapter 4 Numerical Step 3. Decide whether or not to scale the constitutive curve The constitutive curve at an integration point is only scaled during the initial elastic loading of the element. The objective of scaling is to modify the peak stress value so it is not necessary to scale the strength after the peak stress has been reached and the element is softening. Also, to ensure the correct fracture energy, the constitutive curve is not scaled when reloading with a residual modulus. Therefore scaling takes place i f 0<£< £peak A N D Speak < SOpeak- 4.4 In the code a slight adjustment to this criterion was necessary so that scaling would not take place when s is close to speak- Since the strains used for scaling are from the previous time steps, it could occur that the current strain had already passed the new strain at onset of damage, causing numerical errors. This adjustment does not change the effect of scaling in any way since the material is linear elastic before the damage onset strain is reached. If scaling does not take place then step 4 is skipped. Step 4. Determine the new constitutive curve If the element is still elastic and the current strain is less than the strain at onset of damage then scaling takes place. The new damage onset strain is calculated as (for each direction): £ * — £ peak 0 peak f £ ^ V^max 7 4.5 Here s is the strain from the previous time step. A new ultimate strain is calculated so that the fracture energy remains constant: ' 0 C peak P * = LS 4 f. Ene* The modified stress-strain curve after scaling for an integration point is shown in Figure 4.12. Step 5. Calculate stress The current peak strain is tracked throughout the solution and is updated in this step. If scaling took place then 71 Chapter 4 Numerical 4.7 otherwise, Speak —max (s, Epeak , £0peak )• 4.8 W i t h the stress-strain curve n o w def ined, e is s imply mapped to the curve to determine the current stress cr. Step 6 . Return the tangent stiffness tensor to the FE program A B A Q U S requires the tangent stiffness tensor and updated stresses to obtain a converged solut ion. The tangent modu l i are therefore determined and the stiffness tensor is assembled and returned to the code. Step 7. Update history variables The increment o f fracture energy dissipated is calculated and the total dissipated fracture energy is updated. The residual secant modulus is then calculated from the residual fracture energy as The strains at each integrat ion point are also stored i n a global storage array so that they can be used to determine the strain prof i les in Step 1. Continue integration point loop. Continue element loop. 4.3.4 How mesh insensitive is the ASDM? The A S D M procedure is intended to reduce the mesh sensit iv i ty o f the model . The question o f h o w coarse the mesh can be and st i l l return accurate results is discussed br ie f ly . The damage zone size Wdam predicted by the numer ica l model is dependent on the values o f Speak and Suit that are used. I n the A S D M these values can also change throughout the course o f the solut ion. Figure 4.13 il lustrates h o w the shape o f the strain field and the strain parameters affect the size o f the damage zone. I f the distance £dam is increased (i.e. speak is r 2yr m i n E, 4.9 ult peak J 72 Chapter 4 Numerical decreased and/or suit is increased), the width of the damage zone increases, resulting in a tougher constitutive curve. If the element width used to discretize the area in front of a notch tip is increased, the A S D M scales Speak, and consequently crpeak, to a lower value. To maintain the critical fracture energy, suit is increased. As element width is increased, Wdam also increases. Thus the material behaviour is slightly changed in using the scaling technique. Element widths should be less than the width of the damage zone in order to accurately represent it. The size of the damage zone in the numerical model is slightly altered by using the A S D M . Another question is whether large elements will be able to accurately represent the energy released in growing a crack part way through an element. Using the example shown in Figure 4.14, it is seen that the energy released in growing a crack a certain length will be constant regardless of the element width used. In Figure 4.14 a) a crack is grown for 2 mm through two 1 mm wide elements. Each element has a constitutive curve as shown in the figure, with a one-dimensional critical fracture energy of 100 mJ/mm. Since each element is 1 mm wide, the energy released per segment is 100 mJ and the total energy released for a 2 mm wide crack is 200 mJ. In Figure 4.14 b) a single 5 mm wide element is used to discretize the entire 5 mm width in front of the crack. The element's constitutive curve also has a fracture energy of 100 mJ/mm, but since it is 5 mm wide, the energy released per segment is 500 mJ. A crack growth of 2 mm is represented by the release of 40% of the segment's critical fracture energy, equal to 200 mJ. Therefore the energy release is the same. The difference between these two situations is that there is a difference in the residual strength of the interface at the location of the initial notch tip. It appears that the main limitation on the width of elements used to discretize the area ahead of a notch is the ability of the elements to resolve the stress and strain fields. 4.4 OCT simulations 4.4.1 T h e v i r t u a l specimens Seven O C T specimens containing notches with root radii of 1 mm, 4 mm, 10.5 mm, 14 mm, 16 mm, 19 mm and 27 mm were modeled. To be consistent with the experimental tests, each 73 Chapter 4 Numerical v i r tua l specimen is referred to using its notch diameter, for example rd2 for the specimen w i t h a notch root radius o f 1 m m . A s w i t h the experimental specimens, the nomina l dimensions o f the v i r tual specimens are the same as g iven in Figure 3.1 but the actual v i r tua l dimensions vary s l ight ly to match the experimental dimensions o f each separate specimen. Depending on the case being run, the mesh size and conf igurat ion var ied. Ha l f - symmet ry versions and fu l l versions o f each specimen were used depending on the test. I n al l cases, the solut ions were for plane stress. The fu l l version o f rd21 is shown in Figure 4.15 as a representative mesh. For al l specimens, a band o f rectangular, u n i f o r m size elements were used to mesh the area around the notch plane ahead o f the notch t ip . Element size was increased away f r o m the notch plane. The dimensions and tests for each v i r tua l specimen are tabulated in Table 4.3. I n the numer ical model the loading pins are modeled w i t h an elastic material model using material properties for steel. Non-penetrat ing contact is def ined between the pins and the holes in the O C T specimens. The st i f fening brace along the back edge o f the specimen that was used in the experiments is not inc luded in the model . For al l specimens the height o f elements in the region o f the notch plane was kept constant at 0.5 m m . A s per F l o y d (F loyd , 2004) the fracture energy is determined as the element height t imes the specif ic fracture energy. A constant element height was used so that the w id thwise discret izat ion could be studied independently. A b r i e f note is made here o f an interesting phenomenon. When constant stress elements were used, it was found that strain w o u l d localise to a band o f elements, a single element in height, in the same manner as noted by F l o y d for the composite damage model implemented in the expl ic i t f in i te element code L S - D Y N A . In the second-order elements used in this study however, there are nine integrat ion points in three rows o f three w i t h i n each element. I t was not certain whether strain w o u l d localise to an entire element (i.e. al l three rows o f integrat ion points) or to on ly one or t w o rows w i t h i n the element. I t turned out that in every case tested in this thesis, strain always localised to t w o rows o f integrat ion points - the midd le r o w and one o f the edge rows. Further invest igat ion is necessary to determine the reason for this pattern. However , the height o f the local isat ion band is therefore not the entire element height, but rather a f ract ion o f it. Consider ing the Gauss integrat ion scheme, 74 Chapter 4 Numerical each integrat ion po in t represents a weighted area o f the element. The integrat ion po in t in the midd le is weighted by 8/18 and integrat ion points on the edges are weighted by 5/18. The element height therefore needs to be mu l t ip l ied by 13/18 to get the actual height o f the crack band. 4.4.2 Parametric study on Gc and apeak The S D M was used to investigate the sensit iv i ty o f the model to the t w o input parameters cont ro l l ing the fracture response. Gc and apeak were var ied independently on a half-symmetry version o f rd2. Elements w i t h a w i d t h o f 0.5 m m were used to discretize the notch plane region. The simulat ions used the basic material properties g iven in Table 4 . 1 . I t was on ly found necessary to s l ight ly increase E2 to 32 GPa ( f rom the or ig inal value o f 30 GPa f r o m the material supplier) to match the elastic compl iance o f the experimental tests. The load versus C M O D plo t shown in Figure 4.16 demonstrates that as Gc is increased, the magnitude o f the entire post-peak por t ion o f the curve is increased. V a r y i n g Opeak does not change the overal l magnitude o f the curve but s l ight ly adjusts the peak load value that is attained. A n increase in apeak leads to an increase i n peak load as w e l l . The overal l sensit iv i ty o f the solut ion to these parameters is not large w i t h i n the range tested. The experimental results seem to be bounded by the 70 k J / m 2 curve and the 85 k J / m 2 . I t was therefore decided to use the parameters corresponding to the rough experimental values determined in Chapter 3 for the rest o f the simulat ions. That is a strength o f 460 M P a and a r) cr i t ical fracture energy o f 80 kJ /m . 4.4.3 SDM/ ASDM comparison T o test the effectiveness o f the A S D M , ha l f -symmetry simulat ions were conducted o f v i r tua l specimen rd8 using 0.5 m m and 1.0 m m w ide elements in the region o f the notch plane. I n both cases the height o f the elements was 0.5 m m . The val idat ion tests are l isted in Table 4.4. First, to examine the effect o f using di f ferent sized elements, the S D M and constant stress elements were used for both element wid ths w i t h the exper imental ly determined propert ies; Gc equal to 80 k J / m 2 and apeak equal to 460 M P a . The l o w end o f the exper imental ly determined Gc values were found to w o r k the best. A stable solut ion was obtained for both mesh sizes and surpr is ingly no dif ference in response was observed between the coarse and 75 Chapter 4 Numerical fine mesh cases. The load versus C M O D plots were ident ical w i t h a peak load o f 16.1 k N in both cases. F r o m the results o f the previous simulat ions on D C B specimens this outcome was not expected. A s imulat ion using the same input parameters was per formed w i t h 0.5 m m wide second order elements. A g a i n the peak load was approximately the same, as noted in Table 4.4. The reason for this outcome can be explained s imply. Referr ing to the elastic stress d is t r ibut ion for a 4 m m notch radius O C T specimen in Figure 4 . 1 , the stresses near the notch t ip are higher than the tensile strength to a distance o f almost 2 m m ahead o f the notch t ip when the load is 15 k N . The stress cr i ter ion for fracture is therefore satisfied in al l elements w i t h i n this d imension and the condi t ion for fracture to progress is the fracture energy equal to the cr i t ical fracture energy. Since al l three meshes have the same cr i t ical fracture energies, there is v i r tua l ly no dif ference in their responses. I t appears that as long as the elements are smaller ( in w i d t h ) than the w i d t h o f cr i t ica l ly stressed material at fracture, the mesh is adequately ref ined and no scaling is necessary. Test ing this theory further, the input Gc was changed to 10 kJ /m wh i le keeping al l other parameters the same. The numerical solut ion became unstable at the peak load indicat ing a transit ion in behaviour from one o f stable crack g rowth (plateau in the load-displacement plots) to one o f temporary instabi l i ty (sharp load drop). The peak load values st i l l permi t ted a comparison to be made. I n this case a dif ference in peak load was seen between the 0.5 m m wide and 1.0 m m wide constant stress elements. The dif ference was approximately 6% w i t h peak loads o f 10.3 and 10.9 k N for the 0.5 m m and 1.0 m m cases respectively. Us ing quadratic elements w i t h the S D M , the peak loads were closer together and between the constant stress values at 10.4 and 10.5 k N for 0.5 m m wide and 1.0 m m wide elements respectively. Nex t the A S D M was used to model the 4 m m notch root radius specimen w i t h Gc equal to 10 k J / m 2 . Figure 4.17 shows the reduct ion in e*peak in the nine integrat ion points o f the element situated on the notch plane immediate ly ahead o f the notch t ip for a mesh o f 1.0 m m wide elements. No te that the midd le co lumn o f integrat ion points have a larger reduct ion in e*peak because they represent a larger area o f the element. For both 0.5 m m wide and 1.0 m m wide elements, the predicted peak load was 10.1 k N , a sl ight reduct ion from the S D M results but i t 76 Chapter 4 Numerical is consistent, showing that there is no dif ference in using 0.5 m m or 1.0 m m w ide elements when the A S D M is used as the material model . I t turns out that i n the case o f the O C T specimens modeled in this thesis, the tough mater ial makes scaling for element w i d t h v i r tua l ly unnecessary i f elements are adequately ref ined to match the stress dist r ibut ion. A n effect o f element w i d t h was on ly not iced when the cr i t ical fracture energy was reduced, mak ing the material more br i t t le. Element w i d t h also becomes more important i f the notch root radius is increased because the structural response becomes more "b r i t t l e " . The w i d t h o f cr i t ica l ly stressed material is less in larger notch root radi i specimens when the energy cr i ter ion is met than in smaller notch root radi i specimens (see Figure 4.1). For notch root radi i larger than the transit ion radius, the cr i t ica l ly stressed w i d t h is close to zero so f in i te sized elements are always larger than the cr i t ica l ly stressed w i d t h . I f the onset o f fracture is control led by the stress cr i ter ion as opposed to the energy cr i ter ion, element w i d t h is very important. Elements larger than the cr i t ica l ly stressed w i d t h o f mater ial when the energy cr i ter ion is met w i l l affect the structural response. The effect is mi t igated, however, by the fact that the stress concentrat ion becomes shal lower as the notch root radius increases. Therefore, structures that are most sensitive to element w i d t h are ones made o f br i t t le material so that fai lure occurs due to satisfaction o f the stress cr i ter ion, but w i t h relat ively sharp notches so that there is a large stress concentrat ion. This observation is consistent w i t h l i terature. Guidel ines for choosing an appropriate element size are discussed in Chapter 5. B o t h the S D M and the A S D M have been shown to w o r k w e l l for s imulat ion o f the O C T specimens and the A S D M is used to per fo rm the s imulat ion o f the experimental tests. 4.4.4 S i m u l a t i o n o f e x p e r i m e n t a l OCT tests The A S D M model was used to simulate the experimental O C T tests using the seven v i r tua l specimens l isted in Table 4.3. A l l o f the specimens were modeled in fu l l as plane stress problems w i t h un i f o rm, 8-noded 0.5 m m x 0.5 m m elements in the region o f the notch plane. The material and cohesive properties g iven in Table 4.1 were used for al l o f the specimens. For an element height o f 0.5 m m and local isat ion in t w o rows o f integrat ion points w i t h i n an element, the ye corresponding to a Gc o f 80 kJ /m is 2.215 x 10 kJ /m . 77 Chapter 4 Numerical In i t ia l l y the simulat ions were conducted using ha l f models o f the specimens w i t h symmetr ic boundary condit ions appl ied along the plane o f the notch. This was found to be an adequate representation for most o f the specimens, but at large notch radi i the crack w o u l d sometimes g row o f f o f the centre l ine causing small errors in the peak load determined. Simulat ions were next conducted on a fu l l mesh b y ref lect ing the mesh for the ha l f -symmetry models about the notch plane. Th is resulted in more consistent predict ions w i t h the crack either propagat ing through the r o w o f elements immediate ly above the notch plane or immedia te ly be low the notch plane. For v i r tua l specimen rd2 however, the h igh strain concentrat ion in the re lat ively large 0.5 m m h igh elements at the notch t ip caused local isat ion to occur simultaneously in both elements on either side o f the notch plane before propagat ing a long on ly one row. This caused an error in the predicted peak load. A l l o f the meshes were therefore changed so that there was a r o w o f elements centred along the-notch plane ahead o f the notch t ip . This conf igurat ion resulted in consistent predict ions from all specimens. A couple o f notes to keep in m i n d when meshing w i t h a strain-softening material model : • Ha l f - symmet ry should not be used so that the damage zone is free to propagate a long the mid-plane or to branch to either side depending on the condit ions • I f a re lat ively sharp notch is being modeled, the notch t ip should not be located at a boundary node (mid-s ide nodes are f ine) since softening can then occur in both elements bounding the notch t ip resul t ing in an energy release rate that is too h igh 4.4.4.1 Load-CMOD behaviour The load versus C M O D predict ions for each o f the v i r tua l specimens are p lot ted over the experimental results i n the plots in Append ix A . The peak loads predicted by each v i r tua l test are plot ted versus notch root radius in Figure 4.18. The numerical predict ions match the experimental results remarkably w e l l at almost al l notch radi i . W h i l e the material parameters were determined f r o m the same set o f specimens as the tests that were modeled, the model was not cal ibrated by f i t t ing the load-displacement p lo t o f a s imulat ion at one notch root radius to experiment and then using the result ing input parameters as is often the case in the l i terature. Compar ing the numer ica l ly predicted peak loads to the experimental results, some trends are apparent. First, the sharpest notch t ip specimens ( for p equal to 0.5 to 2 m m ) have the lowest 78 Chapter 4 Numerical peak loads, increasing f r o m about 14.5 k N to 15.5 k N at p equal to 3 m m . Between notch root radi i o f 3 m m and 14.15 m m there is very l i t t le increase in peak load, w h i c h is matched very closely by the numer ica l predict ions. A t approximately p equal to 16 m m , there is clearly a transit ion po in t at w h i c h peak loads begin to increase w i t h increasing p. This po in t w i l l be referred to as position- A g a i n , the numer ical predict ions also predict this t ransi t ion. There is a deviat ion in the predicted peak load f r o m the experimental results at the largest notch root radi i . For p equal to 27 m m , the predicted peak load is approximately 19.5 k N cont inu ing the t rend o f increasing peak load w i t h p, wh i le the experimental result levels o f f at approximately 17.8 k N . No te that there are t w o experimental points for p equal to 27 m m located together. Invest igat ing the ind iv idual load versus C M O D plots more closely, i t can be seen that the numer ica l ly predicted shape o f the curves is s imi lar to the experimental results. A t small notch root rad i i , the material begins to soften early, causing a w ide , rounded peak on the load versus C M O D plot. This rounded peak in the numer ica l specimens represents the jagged plateau seen exper imental ly in the smal l p specimens. See for example Figure A . l and Figure A . 2 . A s p increases, softening o f the material is delayed to higher loads, creating a more l inear response up to the peak load. A t large p the peak o f the curve is more pronounced. For specimens rd21 and rd28, the peak is s l ight ly more pronounced in the numer ica l specimens than seen exper imental ly. Especial ly for rd21 (Figure A 3 ) , the experimental specimen has a more noticeable plateau w i t h a s l ight ly lower peak load, perhaps indicat ing that the notch t ip o f specimen rd21-1 contained a small f law, causing i t to damage prematurely. In specimen rd32, the sharp load drop, caused by a sudden j u m p in the crack length after the peak load, is w e l l represented by a steep drop in the numer ical curve. For the t w o largest notch root rad i i , 38 m m and 54 m m , shown in Figure A . 6 and Figure A .7 respectively, the experiments showed unstable crack g rowth f o l l o w i n g the peak load w i t h very large load drops before the crack was arrested. The numerical simulat ions were unable to handle this instabi l i ty and the runs crashed at or jus t after the peak loads were reached. In the imp l i c i t integrat ion scheme used, the solut ion is stable as long as the global stiffness matr ix is posi t ive def ini te. Therefore the numerical model can handle localised softening but crashes i f there is global instabi l i ty as we see at large p. 79 Chapter 4 Numerical 4.4.4.2 Stress fields The stress field contours, for stress normal to the crack, for the A S D M simulat ion o f rd2 are shown in Figure 4.19. A s C M O D is increased a) a stress concentrat ion appears at the t ip o f the notch. B y the t ime the peak load is reached b) the elements at the notch t ip have begun softening and the stress concentrat ion has detached f r o m the in i t ia l notch t ip . A considerable amount o f damage is predicted to occur before the peak load is reached. A t c) the C M O D is 3.21 m m and the damage zone has progressed rough ly 30 m m . I n Figure 4.19 b) and c) a " t a i l " o f h igh stress behind the stress concentrat ion is apparent. Th is ta i l is due to the cohesive zone const i tut ive curve. The stresses along the tai l represent par t ia l ly damaged material on the softening por t ion o f the cohesive zone stress-strain curve. For comparison, the predicted stress contours for rd54 are p lot ted in Figure 4.20. A t a) the C M O D o f 0.67 m m is approximately the same as in Figure 4.19 a) for rd2. W h i l e rd2 already shows normal stresses at the notch t ip equal to the material strength, the notch t ip stresses for rd54 are we l l be low cr i t ical . Even at a C M O D o f 1.70 m m as shown in Figure 4.19 b) and Figure 4.20 b) for rd2 and rd54 respectively, the stresses are be low cr i t ical for rd54 wh i le rd2 is at its peak load. Short ly after, rd54 reaches its peak load at a C M O D equal to 2.48 m m , shown in Figure 4.20 c). The stresses at the notch t ip have jus t reached the material tensile strength and the stress concentrat ion has not yet detached from the in i t ia l notch t ip. The s imulat ion crashes at this point. A s hypothesized f r o m the exper imental results, the specimen goes unstable as soon as the notch t ip stresses reach the tensile strength. A l s o , as seen in the elastic s imulat ion results, the height o f h igh ly stressed material jus t pr ior to peak load is much greater in the larger notch root radius specimen. A band approximately 20 m m h igh is stressed at close to the cr i t ical level when v i r tua l specimen rd54 becomes unstable. Conversely, in rd2 the height o f material stressed at close to the tensile strength increases from on ly a couple o f mi l l imetres at the notch t ip to several mi l l imetres at the point o f peak load. I t appears that the height o f the cr i t ical stress zone either grows to the stable propagat ion height for sharp notches (stable behaviour) or decreases to the stable propagat ion height for b lunt notches (unstable behaviour) . For the notch root radius at the transit ion from stable to unstable behaviour, the stress contours should remain rough ly the same size as the crack propagates. The stress contours at three C M O D values are p lo t ted for specimen rd32 (16 m m notch root radius) in Figure 4 .21 . 80 Chapter 4 Numerical The magnitude o f the stress contours does in fact reduce as the crack propagates. This may be due to there being a large cut-out o f material at the in i t ia l notch t ip that alters the gross stress dist r ibut ion away from the notch plane. A s the crack t ip progresses away from the in i t ia l notch t ip , the effect o f the cut-out is diminished. 4.4.4.3 Crack width X Y - p l o t s o f the stress distr ibut ions in the damaging band o f elements (along the notch plane) are useful for look ing at the propagat ion o f the crack front. Stress d is t r ibut ion propagat ion plots for rd2 and rd28 are presented i n Figure 4.22 and Figure 4.23 respectively for increasing values o f C M O D . I t was not possible to create one o f these plots for rd54 because the run crashed after crack in i t ia t ion. Plots for rd8, rd21 and rd32 are i n Append ix C. The posi t ion o f the crack front is determined from the locat ion o f the peak stress for a g iven C M O D or t ime. This posi t ion marks the w i d t h o f material that has damaged, whether f u l l y or part ia l ly . The poin t at w h i c h the t ra i l ing edge o f the stress curve goes to zero marks the w i d t h o f fu l l y cracked material ( interface w i t h zero residual strength). The w i d t h o f the damage zone in the numerical model is the distance between these two points Wdam,e- A s discussed in Section 4.3.4, this damage w i d t h is not the same as the physical w i d t h o f the damage zone because its value is strongly dependent on the shape o f the stress-strain curve and somewhat dependent on mesh size. The shape o f the stress-strain curve determines the value o f suit for g iven y and epeak inputs. The bi l inear stress-strain relat ionship assumed numer ica l ly is on ly an approx imat ion o f the physical stress-strain relat ionship. U n t i l this point , i t has been suff ic ient for produc ing accurate load-displacement predict ions because the two factors that inf luence the predict ions are Gc and apeak, w h i c h are accurately captured by the bi l inear model . However , the w i d t h o f the damage zone depends on euit, wh i ch is not necessarily accurate i n the bi l inear model . I f i t is assumed, for the moment , that the actual const i tut ive behaviour o f a characteristic vo lume o f the material is b i l inear, as in the numerical model , then the characteristic damage zone w i d t h o f the material Wdam,c can be determined from the characteristic suit and the strain f ie ld as was shown in Figure 4.13. A predict ion for Wdam,c is made using the rat io o f element height ( local isat ion band) to the characteristic damage height. Us ing Equat ion 2.14 w i t h the exper imental ly determined 81 Chapter 4 Numerical characteristic damage height o f 6.5 m m and a Gc o f 80 m J / m m 2 , yc is approximately equal to 12.3 m J / m m 3 . The elemental specif ic strain energy ye is 221.54 m J / m m 3 for the 0.5 m m h igh quadratic elements. Us ing Equat ion 4.2, the ul t imate strain at fai lure in the numer ica l model is then determined to be 0.963. The strain pro f i le for rd28 at C M O D equal to 3.9 m m is shown in Figure 4.24. For a strength o f 460 M P a and in i t ia l modulus o f 32 GPa, the locat ion o f the peak and ul t imate strains are also plot ted. The numer ica l w i d t h o f the damage zone Wdame is read as approximately 18 m m . The characteristic w i d t h o f the damage zone is m u c h smaller. I f a bi l inear consti tut ive curve is used as the actual stress-strain response o f the mater ial , an ul t imate strain o f 0.053 is obtained using Equat ion 4.2. A s shown in Figure 4.24, this corresponds to a Wdam,c o f approximately 1 m m . F r o m the specimens that were sectioned, this value seems too small . Sect ioning suggests a damage zone w i d t h o f about 6 m m . Therefore, the physical stress-strain softening behaviour o f the material must not be linear. 4.4.4.4 Compliance plots Us ing the peak stress locat ion to determine the w i d t h o f cracked mater ial , compl iance plots are generated for each specimen. The numerical compl iance, as determined using the A S D M , is plot ted in Figure 4.25 for rd28. The locations o f the peak stresses i n Figure 4.22 and Figure 4.23 g ive the w i d t h o f material that contains any amount o f damage for a given C M O D . The w i d t h o f material that is fu l l y cracked is less. This can be seen in Figure 4.25 by compar ing the A S D M compl iance predict ion to the compl iance p lo t as determined from a series o f static, elastic simulat ions in w h i c h the interface constraints were successively removed. Exper imenta l points f r o m the l ine analysis determinat ion o f crack length are also inc luded on the plot . B y shi f t ing the A S D M compl iance predict ion by 3.25 m m along the C M O D axis, i t is superposed on the elastic predic t ion, as shown in Figure 4.26. S imi lar to what was seen exper imental ly when the length o f damaged material f r o m sectioning was compared to the crack length determined by line, analysis, the length o f the crack inc lud ing part ial damage is approximately equivalent to a shorter, but fu l l y separated crack. Here, shi f t ing the crack w i d t h by 3.25 m m places the t ip o f the equivalent discrete crack in the centre o f the approximately 6 m m w ide damage zone that was determined f r o m sectioning. 82 Chapter 4 Numerical 4.4.4.5 Crack velocity The static, imp l i c i t model used does not a l low quanti tat ive predict ions o f crack ve loc i ty to be made because inert ia is not included in the solut ion. B y p lo t t ing the posi t ion o f the crack front versus a normal ised temporal parameter, however, the sharp j u m p in crack posi t ion at large notch radi i can be seen compared to smaller notch radi i . I n Figure 4.27 the posi t ion o f the crack front is plot ted versus normal ised s imulat ion t ime t/to. The total durat ion o f the s imulat ion is to. V i r tua l specimen rd2 begins to fracture early and the increase in crack length w i t h t ime is gradual. A s notch root radius is increased, the in i t ia t ion o f fracture is delayed. A t p equal to 16 m m (rd32), fracture initiates at a normal ised t ime o f about 0.4 and the posi t ion o f the crack front j umps f r o m 0 m m to 20 m m rapid ly . B y about a normal ised t ime o f 0.6, the crack fronts for al l o f the notch radi i are almost coincident. 4 . 5 Summary • Static, elastic simulat ions o f the O C T specimens help expla in the transi t ion in behaviour seen exper imental ly . For re lat ively sharp notches, the normal stress at the notch t ip attains the tensile material strength pr ior to the onset o f structural softening. O n the other hand, the elastic stress f ie ld solutions reveal that the normal stress at the notch t ip is less than cr i t ical for large radi i when the peak load is attained for a sharp notch. • A three-parameter cohesive model cal led the S D M was developed as a U M A T in A B A Q U S . The S D M was mod i f ied to a l low it to automat ical ly adapt its const i tut ive behaviour based on element size and the strain d is t r ibut ion to account for the effect o f element w i d t h on peak load predict ion. The model is cal led the A S D M . • Us ing the exper imental ly determined Gc and apeak, the A S D M did a good j o b o f predic t ing the experimental O C T results and also displayed a transi t ion in behaviour at a notch root radius o f about 16 m m . Rapid crack g rowth in v i r tual specimens w i t h notch root radi i larger than the transit ion radius was unstable and the simulat ions stopped prematurely. • Scal ing to account for element w i d t h was found to be unnecessary for tough materials as long as the element w i d t h was less than the w i d t h o f cr i t ica l ly stressed material at 83 Chapter 4 Numerical the onset o f softening. Element w i d t h becomes increasingly important for br i t t le materials and sharp notch t ips. • The physical damage zone w i d t h cannot be predicted using the S D M w i thou t assuming a shape for the actual physical stress-strain relat ion. Crack length however can be re l iab ly predicted using the S D M by p lo t t ing the compl iance o f the v i r tua l specimen. • Simulat ions o f the O C T specimens using the A S D M reveal that considerable damage and localised softening takes place ahead o f re lat ively sharp notch tips before discrete crack format ion and the onset o f structural softening. For specimens w i t h notch root radi i larger than the transi t ion radius however, very l i t t le damage occurs pr io r to discrete crack fo rmat ion and the onset o f unstable crack advance. 84 Chapter 4 Numerical Table 4.1 Material properties for numerical OCT tests Property Value Ex 75 GPa Ey 32 GPa Oxy 0.161 Gxy 17.1 GPa @peak 460 M P a GF 80 k J / m 2 Table 4.2 Height of critically stressed material at peak load Notch root radius . (mm) Approximate peak load (kN)* Height of critically stressed material (mm) 1 15 3.5 4 15 4.75 10.5 15 6 16 16.9 9.5 27 17.8 8.5 Taken f r o m O C T experiment results Chapter 4 Numerical T a b l e 4.3 N u m e r i c a l tests V i r t u a l Spec imen a ( m m ) W ( m m ) Mesh Size/ C o n f i g u r a t i o n s Tests rd2 33.0 81.5 H a l f model : 0.5 m m quads Fu l l model : 0.5, 1.0 m m quads 1) Parametric study o f Gc and Gpeak 2) Exper iment s imulat ion rd8 33.0 81.8 H a l f model : 0.5, 1.0 cs 0.5, 1.0 quads Fu l l model : 0.5 m m quads 1) S D M / A S D M compar ison 2) Exper iment s imulat ion rd21 34.0 81.2 Fu l l model : 0.5 m m quads Exper iment s imulat ion rd28 33.8 80.0 Fu l l model : 0.5 m m quads Exper iment s imulat ion rd32 34.0 81.5 Fu l l model : 0.5 m m quads Exper iment s imulat ion rd38 33.0 81.0 Fu l l model : 0.5 m m quads Exper iment s imulat ion rd54 33.0 81.2 Fu l l model : 0.5 m m quads Exper iment s imulat ion 86 Chapter 4 Numerical Table 4 . 4 ASDM validation tests Sim# Model Gc (kJ/m2) Element type Element width (mm) Peak Load (kN) • 1 S D M 80 Const, stress 1.0 16,1 2 S D M 80 Const, stress 0.5 16.1 3 S D M 80 Quadratic 0.5 16.2 4 S D M 10 Const, stress 1.0 10.9 5 S D M 10 Const, stress 0.5 10.3 6 S D M 10 Quadratic 1.0 10.6 7 S D M 10 Quadratic 0.5 10.4 8 A S D M 10 Quadratic 1.0 10.1 9 A S D M 10 Quadratic 0.5 10.1 87 Chapter 4 Numerical 1400 1200 600 -, 550 500 -| ^ .450 g 400 w 3 5 0 Vi 2 300 " I •3 2 5 0 1 I 200 Z 150 \ 100 50 J 0 0 Notch plane 6 8 10 12 14 Distance along notch plane (mm) 16 20 Figure 4.1 Elastic stress fields ahead of notch tip at P — 15.0 kN _Tensile strength 6 8 10 12 14 Distance along notch plane (mm) 16 20 Figure 4.2 Elastic stress fieldfor p — 27 mm at experimentally determined peak load of 17.8 kN 88 Chapter 4 Numerical a) p = 1 mm b) p = 4 mm c) p= 10.5 mm d)p= 16 mm (MPa) '22 +4. S00e+02 + 3. 833e+0Z +3. 067e+02 + 2. 300e+0Z +1. S33e+02 +7. 667e+01 -1. 5ZSe-0S 667e+01 S33e+02 300e+02 067e+02 833e+02 600e+02 988e+02 e)p = 27 mm Figure 4.3 Stress contours ahead of notch tip (symmetric about the notch plane) for a) p = 1 mm at 15 kN load, b) p = 4 mm at 15 kN load, c) p = 10.5 mm at 15 kN load, d) p = 16 mm at 16.9 kN load and e) p = 27 mm at 17.8 kN load. Elements at the notch tip are 0.5 mm x0.5 mm in all cases. The light grey regions indicate a stress of460 MPa. 89 Chapter 4 Numerical 0 4 , , , H : , , H , , , 0 1 2 3 4 5 6 7 8 9 10 Distance along notch plane (mm) Figure 4.4 Stress concentrations on the notch plane ahead of a 4 mm radius notch tip in an OCT specimen as resolved by different sized constant stress (cs) andfully integrated quadratic (quad) elements. Stresses are extrapolated to the nodes. i I I I 1 Figure 4.5 Effect of element width on stress calculated at a notch tip 90 Chapter 4 Numerical 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 8 (mm) Figure 4.6 Load versus arm displacement for DCB simulations showing the effect of element width on the peak load prediction 35 n 30 •a es o -J it « 0-25 20 10 4 0 • ^ Without strength scaling • With strength scaling 0 -I , , , , , , 0 0.1 0.2 0.3 0.4 0.5 0.6 Element width (mm) Figure 4.7 Peak loads for DCB simulation using element meshes of varying coarseness with and without scaling the material strength to match the 62.5 pm mesh size 91 Chapter 4 Numerical 60 55 50 45 H 40 ^ 35 o-S 30 b 25 20 15 10 5 62.5 /um element (100% strength) {--\ 125 /im element (84.5% strength) 250 jum element (64.2 % strength) 500 /um element (42.2 % strength) 62.5 /um solution 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.1 Distance ahead of initial crack tip (mm) 0.9 Figure 4.8 Strength scaling ratios to match the peak load in the DCB simulation using a mesh with 62.5 pm constant stress elements at the crack tip Yc=Yr+ Yd ^Opeak £peak Figure 4.9 Simple bilinear damage model S 92 Chapter 4 Numerical Initialize (for each direction): £peak — ^Opeak Er = E0 Yr=Yz Read user input: ^ 0 1 . £<)2, G ] 2 , V 1 2 , ? c l , 7c2> ^peakb ^peak2 g e t s t a t e v a r i a D | e s Determine the strain distribution in the element (initially all strains are zeroed) Determine the maximum strain, £'ci„, av, across the element No Update the current peak strain: smak = max(£, Calculate the current stress based on the current strain and the constitutive curve Determine the tangent modulus necessary to arrive at current strain Assemble the tangent stiffness tensor Calculate amount of strain energy dissipated Calculate initial sM: 2rc ^0£0 peak Yes ( > £*peak - min(£bp O T f o £ 0 p e a k £ ) V ^elmax ) Update the current peak strain: ^oeak £ peak Loop over integration points Update parameters: Yr=n- Yi E, = m i n ( £ n - — Update history variables: epeak, yn Er Store integration point strains in global matrix for determining element strain gradient Figure 4.10 Flow chart for SDM/ASDM algorithm 93 Chapter 4 Numerical Figure 4.11 Integration point numbering ^ — - ^ S c a l e d c u r v e £ peak Sopeak Figure 4.12 Scaled stress-strain curve 94 Chapter 4 Numerical ^ult \-am Speak am Separated Damaged Elastic x (mm) Figure 4.13 Damage zone size determined by the cohesive zone model Figure 4.14 Example of energy release determined using elements of different size 95 Chapter 4 Numerical a) h) Figure 4.15 Virtual OCT specimen mesh for a) full specimen and b) at notch plane 0 0.5 1 1.5 2 2.5 3 3.5 4 4 CMOD (mm) Figure 4.16 Parametric study of Gc and apeakfor rd2 96 Chapter 4 Numerical 0.016 0.015 ] 0.25 Figure 4.17 Reduction in £*peak at the integration points of the notch tip element in rd8 Z •a « o ct PH 20 -i 19 18 i 17 16 15 14 13 12 J, 11 10 • • # • • B • • O • o 0 15 p (mm) 20 25 30 Figure 4.18 Peak loads predicted using the ASDM (large diamonds) plotted over the experimental results (small squares). Open symbols represent cases that were unstable (i.e. numerical run crashed or did not display load plateau experimentally). Note that the ordinate scale begins at 10 kN to magnify the results. 97 Chapter 4 Numerical Figure 4.19 Normal stress contours for rd2 a) at CMOD=0.679 mm (before peak load), b) at CMOD= 1.70 mm (peak load) and c) at CMOD=3.21 mm (after peak load) 98 Chapter 4 Numerical Figure 4.20 Normal stress contours for rd54 a) at CMOD=0.674 mm, b) at CM0D=1.70 mm and c) at CMOD=2.48 mm (peak load) 99 Chapter 4 Numerical Figure 4.21 Normal stress contours for rd32 a) at CMOD=0.674 mm (before peak load), b) at CMOD=2.11 mm (peak load) and c) at CMOD=3.21 mm (after peak load) 100 Chapter 4 Numerical 500 400 300 200 "3 100 E i-o * 0-| 100 -200 15 20 x (mm) 25 30 35 40 Figure 4.22 Progression of the stress distribution along the notch plane for rd2. Each distribution is at the CMOD value (in mm) noted above it. 3.90 | 4 43^ 4 6 9 4.95 u CO E u o Z -200 15 20 25 x (mm) 30 35 40 45 Figure 4.23 Progression of the stress distribution along the notch plane for rd28. Each distribution is at the CMOD value (in mm) noted above it. 101 Chapter 4 Numerical x (mm) Figure 4.24 Damage zone width determination for rd28 (Note that the crack including the damage zone is.32.25 mm wide at this point and the CMOD is 3.9 mm) 10 15 20 25 A a (mm) 30 35 40 Figure 4.25 Compliance for rd28 before shifting ASDM prediction 102 Chapter 4 Numerical 0.0018 -, 0.0000 -! 1 , , 1 1 , 1 , 0 5 10 15 20 25 30 35 40 A a (mm) Figure 4.26 Compliance for rd28 after shifting ASDM prediction by 3.25 mm to account for partially damaged zone. 45 -| 40 -35 -30 -E 25 -B, « 20 -15 -10 -5 -0 • 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t/t„ Figure 4.27 Progression of the crack front in selected virtual OCT specimens 103 Chapter 5 Discussion and Analysis of Results C H A P T E R 5 DISCUSSION AND ANALYSIS OF RESULTS 5.1 Introduction The experimental and numerical results presented in Chapter 3 and Chapter 4 prov ide evidence o f a transit ion in the load versus displacement behaviour o f a notched structure w i t h increasing notch root radius. A l t hough invest igat ing notch sensit iv i ty in composites (or any mater ial) is not unique, the tests done here have isolated the notch root radius f r o m the notch length and section area in a stable test specimen, w h i c h has prov ided a clear picture o f the mechanisms causing fracture growth . A compl icated model is not needed; in fact the transi t ion in behaviour can be explained s imply using fracture mechanics equations. Trad i t ional ly , the onset o f fracture at a notch t ip is determined using either the energy cr i ter ion for fracture or the stress cr i ter ion as described in Chapter 2. W h i l e both cr i ter ia must be satisfied for fracture to occur, it is usual ly assumed that both cr i ter ia are satisfied simultaneously (Broek, 1982). Under certain condit ions they are not satisfied simultaneously, leading to t w o regimes o f fracture; energy dr iven fracture and stress dr iven fracture. 5.2 Analytical description of notch root radius sensitivity 5.2.1 Regime I Energy dr iven fracture w i l l be cal led regime I. I n this regime the stress cr i ter ion for fracture is satisfied (i.e. a > ac) before the energy cr i ter ion is satisfied (i.e. G > R). This typ ica l l y occurs for relat ively sharp notches, where h igh stress concentrations are present. Since the stress at the t ip has already reached the material strength, the energy in the system controls the onset and propagat ion o f fracture. Under displacement contro l , fracture in this regime is stable because fracture initiates as soon as it is energetical ly favourable. The stored energy in the system never exceeds the crack resistance energy. Specimens r d l through rd28 fa l l into this category. The schematic in Figure 5.1 helps explain the behaviour o f structures in regime I. The stress concentrat ion at a notch t ip increases w i t h appl ied load. Since the notch t ip is re lat ively sharp, stresses become cr i t ical in a small region at the notch t ip at load P' before there is 104 Chapter 5 Discussion and Analysis of Results suff ic ient energy in the global system to cause discrete crack growth . The h igh local energy causes damage (or p last ic i ty in metals), b lun t ing the notch t ip. A s the notch t ip blunts, the notch t ip stress remains constant w i t h increasing load. W h e n G reaches Gc at P!c, the stress cr i ter ion and the energy cr i ter ion are simultaneously satisfied and crack g rowth init iates. W i t h the simultaneous satisfaction o f the stress and energy cr i ter ia, Equat ion 2.8, repeated here as Equat ion 5.1 using the orthotropic modulus, can be used to determine the peak loads in the O C T specimens analyt ical ly. K) = E G, 5.1 For a k n o w n Gc, Equat ion 5.1 is solved for Kjc. Then a relat ionship between load and Kj can be used to determine the peak load for di f ferent geometries. A relat ionship for C T specimens is g iven by (Srawley and Gross, 1972) using a coeff ic ient F2 to a l l ow for vary ing specimen dimensions as in Equat ion 5.2. K, = aNJW-a-F2(a/W) 5.2 Here, <JN is the nomina l stress at the notch t ip , compr is ing stress due to tension and stress due to bending as given in Equat ion 5.3 (Tada, Paris and I r w i n , 1985). °N =  A N + °N Tension Bending f W-a^ 6P a + c r N = ^ — + V t J 5.3 a W-a (W-af _ 2P(2W + a) N~ [W-af The values o f F2 p rov ided in the reference for Equat ion 5.2 are for isotropic materials w i t h standard specimen dimensions. T o determine an F2 value for the orthotropic O C T specimens tested in this thesis, a numerical s imulat ion o f a sharp-notched O C T specimen was conducted in A B A Q U S using an orthotropic elastic mater ial model w i t h plane stress quadratic elements. Tr iangular elements (approx. 50pm) were used at the notch t ip w i t h the mid-s ide nodes moved to the quarter-point o f the sides radiat ing f r o m the t ip. Such a conf igurat ion 1/2 creates an r" s ingular i ty in the stress and strain f ields o f the element to match the expected stress and strain field at the t ip o f a sharp notch (Cook, Ma lkus and Plesha, 1989). The stress 105 Chapter 5 Discussion and Analysis of Results field is displayed in Figure 5.3. The mode I stress intensity factor Kj is determined using a v i r tua l crack closure technique ( V C C T ) f r o m the displacement o f nodes on the quarter-point elements at 9= ±n. The equation for Ki is g iven by (Cook, Ma lkus and Plesha, 1989) as K, = 2G f - V / 2 K + l TC <2lj [ ( 4 v B 2 - v c 2 ) - ( 4 v s l - v c l ) ] . 5.4 where G is the shear modulus, i, B and C are g iven in Figure 5.4 and 3 - v K = - 5.5 1 + v for plane stress. For a g iven load, the value o f K/ f r o m Equat ion 5.4 is substituted into Equat ion 5.2 to determine F2. For the material and specimen geometry studied here a value o f 0.782 was determined. Comb in ing Equations 5 .1 , 5.2 and 5.3 and rearranging, a solut ion for peak load for the O C T specimens in mode I loading is g iven as pl_ {W-a)24WG-c ^ c 2(2W + a)slW -a-F2(a/W) ' The superscript / indicates regime I. Since p does not enter into the equation, Equat ion 5.6 predicts that the peak load is constant over a range o f notch root rad i i . The notch radi i o f all specimens in regime I are small enough that the stress at the notch t ip is above the material strength before the energy cr i ter ion is satisfied. Therefore, the solut ion for a sharp notch t ip should be equivalent to the solut ion for a blunter notch t ip that is st i l l smaller than the transi t ion radius. 5.2.2 Reg ime I I The second regime is the stress dr iven regime, wh ich , for a g iven material system, comprises specimens w i t h notch root radi i larger than those in regime I (e.g. rd32 through rd54). For very b lunt notches or unnotched specimens, the stress concentrat ion is sl ight or non-existent. Therefore, i t is possible to store excess strain energy in a specimen before the stress at any po in t reaches the tensile material strength. Fracture init iates when stress eventual ly reaches 106 Chapter 5 Discussion and Analysis of Results the strength. A s w i t h regime I, the stress and energy cr i ter ia are shown schematical ly for regime I I in Figure 5.2. Due to the blunt notch t ip the stress concentrat ion is less in this regime and the notch t ip stresses therefore increase less rap id ly w i t h load. The energy in the system reaches the cr i t ical level for crack propagat ion at P" before the stresses reach the cr i t ical level. Damage does not ini t iate at this point . The energy in the system continues to increase unt i l F 7 ^ when the stress cr i ter ion is also satisfied. A discrete crack init iates at this point . The release o f the excess energy when cracking initiates drives the crack fo rward i n an unstable manner unt i l the combinat ion o f a and P br ings the energy release rate to less than or equal to the crack resistance. A t this point , the effect ive notch root radius o f the crack has sharpened f r o m its in i t ia l size to correspond to the height o f the mater ia l 's characteristic damage height. A f te r a crack j u m p o f several mi l l imetres in the blunt notched O C T specimens, the large drop in load caused cracking to arrest. Stable crack g rowth then ensued at the characteristic damage height. A n estimate o f the peak load in regime I I for vary ing notch root radi i in the O C T specimens is der ived as fo l lows. The normal stress dist r ibut ion equation o f (Creager and Paris, 1967) for b lunt notch t ips is g iven in Equat ion 2.11 in terms o f the sharp notch stress intensity factor and repeated here as Equat ion 5.7. ^2nr' 2r ^2nr' The effect ive distance f r o m the notch t ip is g iven by r'=r + - p . 5.8 2 The m a x i m u m stress occurs at the notch edge (i.e. r = 0) therefore, ' - • ^ • 0 ) - ^ J i ( f c j + ^ f e 5.9 2K, (7 — max np 107 Chapter 5 Discussion and Analysis of Results Fracture in i t ia t ion occurs when amax is equal to the material strength ac. The cr i t ical stress intensity factor KJc is therefore K i c = — • 5.10 For the O C T specimen, the relat ionship between Ki and P is k n o w n f r o m Equations 5.2 and 5.3 in regime I such that (W - of K * = - £ zr^^W-aF,. 5.11 Equat ing Equations 5.10 and 5.11 yields an estimate for the peak load i n regime I I as g iven in Equat ion 5.12. pu = CTcJ^p~{W-a)2 5 ^ 4(2W + ayjW-a F2 I n this regime the peak load is related to the square root o f the notch root radius. 5.2.3 Transition radius The notch root radius at w h i c h the specimens transi t ion f r o m the stable regime I behaviour to the unstable regime I I behaviour is obtained by equating Equations 5.6 and 5.12 y ie ld ing 4E'GC rtransition ? * 7l(Tc A s w i t h the numerical model , knowledge o f on ly three material parameters; effect ive modulus, cr i t ical fracture energy release rate and tensile strength, is enough to determine the transi t ion radius between stable and unstable behaviour. 5.2.4 Peak load solution for the OCT specimens Us ing the same parameters as were used for the S D M to mode l the O C T specimens tested in Chapter 3, the analyt ical peak load solut ion for vary ing notch root radi i is obtained. N o m i n a l specimen dimensions are given in Figure 3 .1 . However , since the actual specimen dimensions var ied s l ight ly f r o m these nomina l dimensions for each specimen (see Table 3.2), the solut ion is presented as a band encompassing the var iat ion in actual dimensions. The effect ive modulus E' is calculated using Equat ion 2.12 and the modu l i f r o m Table 4 . 1 . The 108 Chapter 5 Discussion and Analysis of Results tensile strength is taken to be the same as the peak stress for the S D M , 460 M P a , and Gc is also the same as the S D M , at 80 k J / m 2 . The analyt ical ly predicted peak loads are presented in Figure 5.5 on top o f the experimental and numerical results. The analyt ical predict ions match the experimental and numer ica l results fa i r ly closely. I n regime I, the experimental and numerical results for the most part fa l l inside the analyt ical predict ion band, however, the lower peak loads at the sharpest notch tips ( for p equal to 0.5 m m to 2 m m ) are not captured by the analyt ical predict ion. The analyt ical model does not expla in this notch sensit iv i ty at re lat ively sharp notch t ips. Theoret ical ly , the development o f a damage zone should el iminate notch sensit iv i ty for al l radi i in regime I. The trend is also predicted by the numer ical model , re in forc ing the experimental results, but the exact reason for the trend has not been determined. A s shown in Figure 4.19, Figure 4.20 and Figure 4 . 2 1 , the shape o f the normal stress field at a notch t ip is disrupted by the presence o f a large radius notch t ip. I t is speculated that the reason for the notch sensit iv i ty at re lat ively sharp notch t ips is due to this dif ference in the shape o f the stress f ields, however further invest igat ion is necessary to c lar i fy the reason for the trend. The presence o f this small effect at sharp notch tips does not affect the overal l a im o f this analyt ical model , wh i ch is to g ive reasonably accurate predict ions o f structural strength and to determine the transit ion w i t h notch root radius from stable to unstable behaviour. Exper imenta l ly and numer ica l ly , the transit ion point was found to be at a radius o f about 16 m m . Ana ly t i ca l l y , the transi t ion notch root radius is calculated to be s l ight ly larger at 17.67 m m . The regime I I predic t ion and locat ion o f ptransition is sensitive to the strength parameter. Increasing the strength by 4 % to 480 M P a yields a ptranSition o f 16.2 m m , closer to the experimental results. A l so the regime I predict ion is very sensitive to Gc. Decreasing Gc 2 2 from 80 kJ /m to 75 kJ /m wh i le keeping al l else constant, changes position to 16.6 m m . A f te r the transi t ion point , the analyt ical peak load predict ion increases along the regime I I curve. Here, the predict ion o f peak load fo l lows more closely the numer ical predict ion using the A S D M than the experimental results, w h i c h actual ly saw a sl ight decrease in peak load at the largest notch root radius. Go ing back to the elastic stress field contours presented in Figure 4.3, an unexpected result was seen. The height o f cr i t ica l ly stressed material at the 109 Chapter 5 Discussion and Analysis of Results experimental peak load o f 17.8 k N for p equal to 27 m m was less than the height o f cr i t ica l ly stressed material at peak load for p equal to 16 m m . N o w , using the peak load o f approximately 19.5 k N , predicted both numer ica l ly and analyt ical ly for a 27 m m notch root radius, the cr i t ica l ly stressed height is 17 m m , considerably larger than the cr i t ica l ly stressed height o f a 16 m m notch root radius. A statistical explanat ion for the experimental decrease in structural strength at the largest notch size is proposed. W e i b u l l theorized in 1939 that a size effect on strength exists due to random f laws in a cont inuum. In regime I I , a structure w i l l fa i l when the stress in a small element o f the structure reaches the cr i t ical strength. For a g iven nomina l stress, the probabi l i ty o f a structure contain ing a small element w i t h a strength equal to the nomina l stress is g iven by the W e i b u l l d is t r ibut ion. The probabi l i ty o f a structure contain ing a small element o f this strength increases w i t h the size o f the structure due to a larger stressed area. This weakest l i nk explanat ion for size effect applies we l l to tensile bars that see a u n i f o r m stress dist r ibut ion. Bazant (Bazant, 2002), however, argues that W e i b u l l theory does not ho ld for notched specimens. Increasing the size o f a notched specimen does not change the magni tude o f the stress concentrat ion at a notch t ip i f the structure is suf f ic ient ly large that the boundary condit ions do not interfere w i t h the stress concentrat ion. Add i t i ona l l y , W e i b u l l theory does not ho ld for structures d isplaying stable crack g rowth . For quasi-bri t t le materials, the format ion o f a damage zone around a notch t ip , pr ior to peak load and the development o f a discrete crack, absorbs the effect o f any f laws in the mater ial . A di f ferent k i n d o f size effect occurs in the O C T specimens though. That is the effect o f increasing the size o f the notch root radius wh i le mainta in ing the size o f the structure. I n regime I, the arguments o f Bazant ho ld because o f the fo rmat ion o f a sizeable damage zone pr io r to the development o f a discrete crack and the peak load. However , in regime I I the specimens fa i l in a br i t t le manner w i thou t the development o f a damage zone pr ior to the peak load. Fai lure in this regime is therefore subject to the occurrence o f f laws w i t h i n the region o f cr i t ica l ly stressed material ahead o f the notch t ip . A s the notch root radius is increased, this cr i t ica l ly stressed region is spread out over a larger height. A s evidence o f this, cracks in the large notch root radius specimens (rd44 and rd54) were seen to init iate s l ight ly o f f centre as shown in Figure 5.6. W i t h more cr i t ica l ly stressed area around the 110 Chapter 5 Discussion and Analysis of Results notch t ip for larger radius notches, there is more l i ke l ihood o f a f l aw being present, leading to lower strengths than expected in very b lunt notches. Further invest igat ion is needed to ver i f y this hypothesis. 5.3 D a m a g e zone size a n d shape I t has been shown that a transit ion in the behaviour o f the O C T specimens occurs at a circular notch root radius o f approximately 16 m m . I t was also shown that the characteristic damage height o f the material is about 7 m m . Is there a l i nk between these t w o values? The height o f the damage zone hd must be related to the extent o f cr i t ica l ly stressed mater ial at peak load. The cr i t ica l ly stressed heights hs for vary ing notch root radi i were presented in Table 4.2. These heights rough ly correspond to the in i t ia l hd that were determined by sect ioning specimens r d 2 - l and rd54-2. For p equal to 1 m m , the cr i t ica l ly stressed height f r o m the elastic s imulat ion was about 3.5 m m . Exper imenta l ly , hd was approximately 3.0 m m . A t the transit ion radius, hs should be equal to hd, w h i c h should be equal to the characteristic stable damage zone height hc, at the in i t ia t ion o f crack g rowth so that the notch neither blunts nor sharpens w i t h the development o f the damage zone. F r o m the elastic s imulat ion for p equal to 16 m m (taken as ptransition), the cr i t ica l ly stressed height is approximately 9.5 m m , s l ight ly larger than, but close to, the stable damage zone height o f 7.0 m m determined from sectioning. Conversely, at the large notch root radius, hs and hd should in i t ia l ly be much higher. F r o m sectioning, the in i t ia l damage height o f rd54-2 was close to 20 m m . A l t h o u g h at the peak load determined exper imental ly , hs was on ly 8.5 m m , there is large var iab i l i ty at this notch root radius. F r o m the numer ica l ly and analyt ical ly determined peak load o f 19.5 k N for p equal to 27 m m , hs is great ly increased to 17 m m . Obv ious ly , the approximate values obtained here cannot be used to develop a quanti tat ive relat ionship between hd and hs, but they do show a trend between the in i t ia l notch size, the extent o f cr i t ica l ly stressed material and the damage zone height. Regardless o f the in i t ia l notch root radius, as a crack progresses, the damage height stabilizes to the characteristic damage zone height hc. Since a t ransi t ion notch root radius creates a 111 Chapter 5 Discussion and Analysis of Results stress d is t r ibut ion that neither blunts nor sharpens w i t h the development and propagat ion o f a damage zone, the question arises o f w h y the plransUion is not equal to hJ2. A s seen in Figure 4.20, a large notch root radius creates a stress f ie ld at the notch t ip that is rather el l ipt ical in shape as opposed to circular. I t is speculated then, that the shape o f the leading edge o f the damage zone is also el l ipt ical . I f this is the case, the radius o f curvature o f the characteristic e l l ipt ical damage zone should be rough ly the same as ptransition- The radius o f curvature o f an ell ipse is g iven by 5.14 Here, A is the semimajor axis and B is the semiminor axis o f the ell ipse. A s shown in Figure 5.7, an ell ipse w i t h a radius o f curvature o f 16 m m , equal to the transit ion notch root radius, can be constructed f r o m a major axis o f 8 m m and a m ino r axis o f 2 m m . These dimensions correspond rough ly to hc and the w i d t h o f cr i t ica l ly stressed material at peak load for a notch root radius o f 16 m m (see Figure 4.3 d). A n ell ipse seems to approximate the shape o f the leading edge o f the damage zone fa i r ly w e l l a l though not exact ly. The radius o f curvature o f an ell ipse, w i t h a major axis the same as the damage zone height, provides a relat ionship to the circular notch root radius at the transi t ion po in t between stable and unstable behaviour. 5.4 S i g n i f i c a n c e o f t h e t r a n s i t i o n r a d i u s There is stronger evidence o f the l i nk between the stable damage height and the transit ion radius. The transit ion radius is the notch root radius at w h i c h the stress cr i ter ion and energy cr i ter ion for fracture are satisfied at exact ly the same t ime. For a notch w i t h radius equal to Ptranshion, the normal stress the notch t ip jus t reaches the material strength as G equal to Gc is also satisfied. Therefore, for all in i t ia l notch root rad i i , the damage zone height converges to the height that gives a notch root radius equal to ptransition so that the stress cr i ter ion and energy cr i ter ion are satisfied simultaneously. I n notches in i t ia l ly blunter than ptransition this results in sharpening o f the notch t ip w i t h g rowth o f the crack and for notches in i t ia l ly sharper than position this leads to b lunt ing o f the notch t ip as damage grows. 112 Chapter 5 Discussion and Analysis of Results 5.5 Other material systems and geometries Throughout this thesis, specimens w i t h notches sharper than position have been cal led stable. Theoret ical ly , this should be true for al l materials. I n br i t t le materials, the size o f ptransitwn decreases, so the l i ke l ihood o f stable fracture is reduced. The transit ion radius is real ly a measure o f toughness o f the mater ial . The larger the transi t ion radius, the tougher the mater ial . The locat ion o f ptransition is, o f course, very sensitive to material properties as given in Equat ion 5.13. For br i t t le materials ( lower Gc) the transit ion radius moves to the left on Figure 5.5. A lso , ptransition is related to one over the strength squared showing a strong inverse relat ionship between strength and toughness. Before the signif icance o f ptransition was realised, four experimental O C T tests were also conducted using an isotropic acryl ic material instead o f the composite. The notch radi i tested were 0.5 m m , 2 m m , 4 m m and 8.5 m m . This material is more br i t t le w i t h approximate material properties g iven in Table 5 .1 . Unstable, fast fracture occurred in al l cases. Us ing these properties and Equat ion 5.13, position is calculated to be 0.174 m m . Theoret ical ly , this material cou ld display stable crack g rowth under displacement control i f the in i t ia l notch root radius was made less than 0.174 m m . The transit ion radius and peak load predict ion p lo t for a unidi rect ional carbon f ibre D C B specimen is der ived as an example in Append ix E. This is an interesting case because i t gives insight into the al lowable thickness o f a Te f lon insert that can be used to st i l l achieve stable fracture. 5.6 Choosing an appropriate element width for strain-softening models A n improved understanding o f the transit ion in behaviour f r o m stable to unstable crack g rowth in notched specimens gives a better understanding o f h o w ref ined a mesh should be in the direct ion o f a propagat ing crack. A s ment ioned in Section 4.4.3, for stable crack g rowth behaviour, the w i d t h o f an element at a notch t ip needs to be smaller than the w i d t h o f cr i t ica l ly stressed material wcrit i f no strength scal ing is used. Elements smaller than wcri, w i l l begin softening pr ior to the satisfaction o f the energy cr i ter ion leading to stable crack growth . I f the notch t ip element is w ider than wcrih as shown in Figure 5.8, damage in the 113 Chapter 5 Discussion and Analysis of Results notch t ip element w i l l be delayed unt i l after fu l f i lment o f the energy cr i ter ion, leading to an incorrect unstable crack g rowth behaviour.. I f the physical crack g rowth behaviour o f the structure being modeled is unstable, there w i l l be an effect o f element w i d t h regardless o f h o w w ide the element is. I n this case it is impract ical to ref ine the mesh to very small w id ths and strength scaling should be used. 5.6.1 E l e m e n t w i d t h case s tudy A n example o f h o w a large element w i d t h , used w i thou t scal ing, can lead to erroneous predict ions is presented. I n the or ig inal simulat ions by F l o y d (F loyd , 2004) o f the experimental O C T tests conducted by M i t che l l (M i t che l l , 2002) , using the C O D A M mater ial model implemented in L S - D Y N A , a relat ively coarse mesh o f 1.25 m m by 1.25 m m constant stress elements was used so that the crack band model cou ld be demonstrated. The specimen geometry, mater ial and orientat ion was the same as for the tests done in this thesis except that the laminate consisted o f f ive stacks instead o f six so that thickness was reduced f r o m 8.5 m m to 7.2 m m and there was loose through-thickness reinforcement in the f o r m o f Kev la r st i tching. The notch t ip had a negl ig ib le radius o f curvature. I n the previous work , damage parameters were chosen to match the numerical s imulat ion to the exper imental ly determined peak load. A strength o f 420 M P a was used and it was found that in order to match the peak load, the input GC needed to be 29 k J / m 2 , considerably lower than the value o f 80 kJ /m determined in the present work . The load versus C M O D result f r o m the C O D A M simulat ion is presented in Figure 5.9. A l t h o u g h a stable response was seen exper imental ly , the numer ical solut ion predicted a br i t t le, unstable response. The fracture energy dissipated dur ing the v i r tual test was determined by measuring the area under the load versus P O D curve (not shown). D i v i d i n g by the area o f crack g rowth gives GF equal to over 40 kJ /m , considerably larger than the input GC. The same mesh and material and damage parameters were used w i t h the S D M model in A B A Q U S and the results are also shown in Figure 5.9. A g a i n , the peak load reasonably matches the experimental result, but an unstable specimen response is predicted. The predic t ion looks wrong and to ver i fy , the analyt ical solut ion in Equat ion 5.6 is used w i t h the same input parameters as used numer ica l ly . Equat ion 5.6 predicts a peak load o f 8.3 k N , considerably di f ferent than the 13.3 k N peak load predicted numer ica l ly in Figure 5.9. The 114 Chapter 5 Discussion and Analysis of Results reason for the error is made clear by p lo t t ing the stress d is t r ibut ion ahead o f the crack t ip using the coarse mesh at the analyt ical ly predicted peak load o f 8.3 k N . Figure 5.10 shows that the numer ica l ly determined stress at the notch t ip has not even reached the tensile strength threshold because the coarse mesh does a poor j o b o f matching the real stress d is t r ibut ion. No te that the downturn near the notch t ip in the coarse mesh stress field predic t ion in Figure 5.10 is not due to damage but rather due to the inabi l i ty o f the mesh to match the real stress dist r ibut ion. A n elastic solut ion o f what the stress d is t r ibut ion should look l ike i f there were no damage is also plot ted in the f igure. 2 2 B y mak ing the fracture energy input 29 kJ /m instead o f 80 kJ /m , the transi t ion radius has shifted to the left. I n addi t ion, the coarse mesh blunts the stress d is t r ibut ion at the notch t ip so that fracture becomes stress contro l led as opposed to energy control led. To achieve a correct predict ion there are two alternatives. The strength o f the const i tut ive model can be scaled d o w n so that damage initiates proper ly . This can be done as a pre-processing step or automat ical ly i n a material mode l such as the A S D M . The other alternative is to ref ine the mesh. Here the mesh is ref ined to second-order 0.5 m m by 0.5 m m elements as used in the rest o f this thesis. The solut ion for the " f i ne mesh" is also g iven in Figure 5.9. A stable response is achieved w i t h a peak load o f about 8.3 k N as predicted analyt ical ly. The stress d is t r ibut ion p lot ted at peak load i n Figure 5.10 shows that damage has occurred in elements up to 3 m m in f ront o f the in i t ia l notch before fai lure. N o w , i t is seen that a Gc input o f 29 k J / m 2 is much too l o w to match the experimental results. Us ing the composite laminate properties used in the rest o f this thesis o f Gc equal to 80 kJ /m and crc equal to 460 M P a , a solut ion very close to the experimental results is predicted as shown in Figure 5.9. The sl ight di f ference may, i n part, be due to the through-thickness st i tching that was not accounted for in the model . 5.6.2 Guidelines for choosing an appropriate element size The element size required to achieve a correct solut ion when using a strain-softening material model is dependent on both material properties and the geometry o f the structure being meshed. For fracture control led by the energy cr i ter ion, the w i d t h o f elements used to discretize the plane ahead o f a notch can be larger than for unstable stress control led fracture. 115 Chapter 5 Discussion and Analysis of Results First, considering regime I where fracture is control led by satisfaction o f the energy cr i ter ion, the mesh size must be smal l enough such that the stress in the notch t ip element cre is higher than the material strength <rc when the energy cr i ter ion is satisfied. This is expressed as cre =aac. 5.15 The factor a has a value greater than or equal to 1.0. Assuming for the moment that the notch t ip is sharp and the element size is small enough to closely approximate the notch t ip stress d is t r ibut ion, Equat ion 5.15 is wr i t ten in terms o f the stress intensity factor as 1 K, ' 5.16 a ^Inr Here, r is replaced by the distance to the calculat ion point o f the notch t ip element. For a constant stress element that distance is ha l f o f the element w i d t h wJ2. A l so , for regime I, the stress intensity factor is replaced w i t h the cr i t ical fracture energy result ing in 1 4^GC <TC=-^J=^. 5.17 v Solv ing for we at the cr i t ical stress gives EG w e = — T T - 5.18 na crc Us ing elements o f size we, the stress cr i ter ion w i l l be satisfied at the same t ime as the energy cr i ter ion for a value o f a equal to 1.0. This is the minimum requirement. The factor a should be a number larger than 1.0 to a l low development o f the damage zone pr ior to the onset o f fracture. A l s o , a should be increased to account for the error in calculat ing the stress field w i t h a coarse mesh, al though this error is small at the locat ion o f the integrat ion point. A s was shown in Figure 4.4 for a 4 m m notch root radius, even fa i r ly coarse constant stress elements are able to match the true solut ion closely at the distance o f the first integrat ion po in t near the notch t ip. The small error in the coarse mesh approx imat ion o f the stress field is incorporated into a. 116 Chapter 5 Discussion and Analysis of Results The above element w i d t h solut ion is for a sharp notch t ip . For a blunter notch w i t h a radius st i l l less than position, smaller elements than g iven by Equat ion 5.18 are required. Figure 4.1 demonstrates that as p is increased, wcrit decreases. A n estimate o f the element size needed for blunter notches is made by using the equation for the stress f ie ld at a b lunt notch t ip o f (Creager and Paris, 1967). For the element calculat ion po in t at we/2 Equat ion 5.17 becomes I f the approx imat ion is made that p/(we+p) « 1, then A s p is increased, element w i d t h should be made smaller. A t a po in t equal to the transi t ion radius, we needs to be equal to zero to achieve the exact solut ion. A s stated in Section 4.3.2, in regime I I when p is greater than ptransition and fracture is stress contro l led, the mesh either needs to be ref ined to small element sizes to reduce the error or strength scaling should be used. O f course as p is increased, the stress concentrat ion is reduced so that the size o f element needed to accurately model the stress f ie ld decreases. A l so , as p approaches in f in i ty , as is the case for an unnotched specimen, the stress concentrat ion approaches 1.0. When there is u n i f o r m stress across the specimen, there is no error in calculat ing the stress at any point w i t h i n the section. Guidel ines for choosing an appropriate element w i d t h are g iven in Table 5.2. Equat ion 5.18 is p lot ted in Figure 5.11 for the material properties used to model the experimental O C T specimens (Table 4.1). Each contour in the p lot represents a di f ferent value o f a. A l ine representing the 460 M P a strength is also plot ted on the f igure, showing that for this tough mater ial , re lat ively large elements can be used w i thout result ing in the stress cr i ter ion being satisfied after the energy cr i ter ion. Re-runn ing v i r tual specimen rd2 (p = 1.0 m m ) w i t h 0.5 m m , 1.0 m m and 2.0 m m w ide constant stress elements (wh i le he maintained constant at 0.5 m m ) , match ing stable solutions are achieved in al l three cases as f 5.19 117 Chapter 5 Discussion and Analysis of Results shown in Figure 5.12. I t is noted that the constant stress elements s l ight ly over-predict the peak load compared to the or ig inal second-order element solut ion. This may be due to the ar t i f ic ia l stiffness that is numer ica l ly added to the constant stress elements to prevent zero-energy hourglass modes. Even the very w ide 2.0 m m elements though give a solut ion s imi lar to the 0.5 m m wide elements. Comparat ive ly , i f Equat ion 5.18 is plot ted for Gc equal to 20 k J / m 2 as in Figure 5.13, the a contours indicate that much more ref ined elements are necessary. Trac ing the 460 M P a strength l ine on the p lot shows that about 0.28 m m wide constant stress elements w o u l d be required to achieve a good solut ion ( for a = 2.0). The s imulat ion o f specimen rd2 was again carr ied out using 0.5, 1.0 and 2.0 m m wide constant stress elements w i t h no strength scal ing to ver i f y the element w i d t h effect. A s shown in Figure 5.14, increasing values o f peak load are predicted w i t h increasing element w id th . The ref ined second-order element mesh produces a softer response w i t h damage occurr ing before the peak load is attained. The s imulat ion st i l l crashed after the peak load was reached because the numerical model cou ld not handle the sharp load drop. The guidelines presented in Table 5.2 can be used to help determine an appropriate element size for a g iven situation. The guidel ines are based on the transit ion radius, w h i c h incorporates material properties and geometry in a single parameter. I n regime I, Equat ion 5.18 can be p lot ted to prov ide a guidel ine for h o w large an element can be, but precise values o f a have not been determined. Setting a equal to 2.0 appears to result in a suf f ic ient ly small element for notches w i t h root radi i less than ptransition-5.7 S u m m a r y • T w o behaviour regimes can be ident i f ied for structures w i t h vary ing notch root radi i . I n the f i rst regime, for sharp notches, the steep stress concentrat ion at the notch t ip ensures that the stress cr i ter ion is satisfied before the energy cr i ter ion. Fracture is therefore contro l led by energy release and is stable under displacement control . I n the second regime, for notches blunter than the transit ion radius, the stress cr i ter ion is satisfied after the energy cr i ter ion is fu l f i l l ed . The stress control led fracture in this regime is unstable. 118 Chapter 5 Discussion and Analysis of Results • A transit ion radius between the t w o regimes can be def ined as the intersection o f the energy and stress cri ter ia. I t is a measure o f the toughness o f the material . • The height o f the characteristic damage zone for a material is related to the t ransi t ion radius. I f the leading edge o f the damage zone is taken to be rough ly e l l ip t ica l , i t has the same radius o f curvature as the transit ion radius. • For notches in i t ia l l y sharper than the transit ion radius, the notch blunts w i t h damage g rowth un t i l the stress cr i ter ion and the energy cr i ter ion are satisfied simultaneously, resul t ing in the characteristic damage height. For notches in i t ia l l y blunter than the transi t ion radius, the notch sharpens w i t h damage unt i l the same condi t ion is met. • Simple guidelines based on the transi t ion radius were presented in Table 5.2. For notches w i t h root radi i less than the transi t ion radius, an element needs to be smaller than the cr i t ica l ly stressed w i d t h o f mater ial at the notch t ip . For notches w i t h root radi i larger than the transi t ion radius, scal ing o f the peak stress in a strain-softening model should be conducted to a l low for larger element sizes. 119 Chapter 5 Discussion and Analysis of Results Table 5.1 Approximate material properties for acrylic P r o p e r t y V a l u e E 3.3 GPa 110 M P a Gc 500 J / m 2 Table 5.2 Meshing guidelines based on notch radius . . . P M e s h i n g a p p r o a c h p ^ Ptransition Use Equat ion 5.18 to determine al lowable element w i d t h w i t h a> 2.0 to a l low for pre-peak softening. A s p —> position, cx should be increased. p ~ Ptransition Elements need to be ref ined to accurately model the notch t ip stress d is t r ibut ion. Scal ing can be used to a l low for larger element sizes. p Ptransition Use strength scaling to a l low for larger element sizes on large structures. Element w i d t h is not important. N o strength scaling is necessary. 120 Chapter 5 Discussion and Analysis of Results 'G. 1.0 Damage development/ notch blunting ^> K{p\ K. 1.0 P' Pc Figure 5.1 Stress and energy failure criteria for fracture in regime I G. 1.0 Excess s * s energy > 1 ' K. 1.0 P" Pe Figure 5.2 Stress and energy failure criteria for fracture in regime II 121 Chapter 5 Discussion and Analysis of Results Chapter 5 Discussion and Analysis of Results Stress Criterion (Regime II) p (mm) Figure 5.5 Analytical prediction ofpeak loads for OCT specimens a) b) Figure 5.6 Cracks initiating off the centerline in a) rd45-l and b) rd44-2 123 Chapter 5 Discussion and Analysis of Results 20 -20 J -20 -15 -10 -5 0 5 x (mm) Figure 5.7 An ellipse with a major axis of 8 mm and a minor axis of 2 mm has the same radius of curvature at the notch tip as a circle with a radius of 16 mm. a Figure 5.8 Stress distribution ahead of a notch tip showing an element of width we that is too large. Mesh must either be refined or strength scaling used to achieve accurate results. 124 Chapter 5 Discussion and Analysis of Results 0 0.5 1 1.5 2 2.5 3 3.5 4 CMOD (mm) Figure 5.9 Load versus CMOD results for a sharp-notched OCT specimen modeled with different size elements 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 x (mm) Figure 5.10 Notch tip stress fields for coarse and fine meshed OCT specimens with a sharp notch tip at load equal to 8.3 kN 125 Chapter 5 Discussion and Analysis of Results 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 C M O D (mm) Figure 5.12 Large constant stress elements used to model the experimental 1 mm notch root radius OCT specimen all give the same result 126 Chapter 5 Discussion and Analysis of Results 0 100 200 300 400 500 600 700 800 900 1000 o-c (MPa) Figure 5.13 Element size contours for Gc = 20 kJ/m2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 CMOD (mm) Figure 5.14 When Gc is reduced to 20 kJ/m a difference is noted in the peak loads predicted using different width elements 127 Chapter 6 Conclusions and Future Work CHAPTER 6 CONCLUSIONS AND FUTURE WORK 6.1 Conclusions Both fracture mechanics and cohesive zone methods have been used in this thesis. The two methods in fact complement each other, providing explanations for different aspects of the other's behaviour. This combined approach has led to a clearer picture of the mechanisms of fracture in FRP composites, as well as, perhaps, other materials. The numerical cohesive model has shown that the LEFM fracture criteria; stress and energy, are not satisfied simultaneously under all conditions. For sharp notch tips, the stress criterion is automatically satisfied for any applied load and the energy criterion drives fracture. For blunt notch tips, the delayed fulfillment of the stress criterion leads to excess energy storage and unstable crack growth. For its part, using the LEFM fracture criteria has provided a means of estimating the width of element necessary to accurately model fracture under different conditions. The development of a simple analytical model, using only three physical input parameters, provides a designer with a simple tool for evaluating candidate materials, designs or other more complicated models. Light has also been shed on the interaction between element size and input parameters for strain-softening finite element material models, which allows greater efficiency as well as certainty in predictions. The following conclusions are made: • From experimental tests and numerical predictions it was found that structural strength is sensitive to notch root radius mainly at large notch root radii. Before a certain transition radius, which is a function only of material properties, the development of a damage zone absorbs the effect of the notch root radius. • For the tough material used for the OCT specimens, the transition radius was found to be about 16 mm. For more brittle materials, the transition radius is less. Theoretically, all materials have a transition radius between stable and unstable behaviour. For very brittle materials in which notch sensitivity is apparent at 128 Chapter 6 Conclusions and Future Work relatively sharp notch tips (e.g. glass), the transition radius is very small and therefore not noticeable in most tests. • The reason for the transition in behaviour is the separate fulfilment of the stress and energy criteria. Fracture can therefore be separated into two regimes. In regime I, the notch root radius is small and the normal stress at the notch tip reaches the material tensile strength prior to satisfaction of the energy criterion. Fracture occurs when the energy in the system is equal to the energy required to form a discrete crack. In regime II, blunter notch tips create smaller stress concentrations. Satisfaction of the stress criterion is therefore delayed until after the fulfilment of the energy criterion. The resulting excess energy in the system when the stress criterion is satisfied drives the crack forward unstably. • A material's transition radius can be related to the characteristic height of the damage zone. If the front of the damage zone is taken to be roughly elliptical like the normal stress field contours, then the radius of curvature of the front of the damage zone is equal to the transition radius. • Sectioning analysis of experimental OCT specimens reveals that the height of damage along a crack plane stabilises to a characteristic height as the crack advances. For notches initially blunter than the transition radius, the damage height is initially larger than the characteristic height and decreases to the characteristic height over a few millimetres of crack growth. For notches initially sharper than the transition radius, the notch tip blunts with the development of the damage zone. Discrete crack formation does not occur until the characteristic damage zone height is fully developed. • A simple numerical cohesive model was developed as a user material model in the ABAQUS finite element code as a demonstrator for modeling techniques that can be applied to any strain-softening material model. Using three experimentally determined parameters, the model closely matches the experimental results. It was determined that the experimental input parameters give a good prediction as long as the element width used to mesh the area around the notch tip is sufficiently refined. In regime I, the element width must be smaller than the width of critically stressed 129 Chapter 6 Conclusions and Future Work material at the notch tip when the energy criterion is satisfied. In regime II, the element width is of critical importance because fracture is driven by satisfaction of the stress criterion and there can be large differences in the stress calculated at an element integration point and the actual notch tip. • An addition was made to the simple damage model (SDM) to allow for larger element widths near the transition radius and in regime II. The adaptive simple damage model (ASDM), automatically scales the input tensile strength while maintaining the fracture energy so that softening initiates in an element when the stress at the element edge reaches the material tensile strength. • Guidelines for selecting an appropriate element width based on the radius of the notch root relative to the transition radius of the material were given. In regime I, the fracture mechanics stress criterion provides an equation for determining an appropriate element size. In regime II, strength scaling should be used if the mesh cannot be adequately refined. As the notch root radius approaches infinity (i.e. an unnotched structure) the stress concentration drops to unity and strength scaling becomes less important. 6.2 Future work There are several areas in which the work in this thesis can be advanced. • First, only one material system, a tough quasi-isotropic composite laminate, was tested in this thesis. It is theorised that the transition radius model presented here is valid for most other materials as well. Validation tests need to be done to prove this. It would be interesting to test the theory further for a more brittle material such as acrylic and also to test it for metals that develop plastic zones before fracture. Testing on other types of specimens and composite laminates would also be useful. • Further investigation is necessary to determine the reason for the initial notch sensitivity in regime I at relatively sharp notch tips. • To be useful, the A S D M technique should be incorporated into a more robust strain-softening composite damage code, such as C O D A M . Doing so would allow it to be tested under more varied conditions. 130 Chapter 6 Conclusions and Future Work • This thesis focused on only mode I fracture. Notch sensitivity in shear mode loading (mode II and mode III) should also be studied. A strain-softening model incorporating interaction between material directions to allow for shear loading should be used to determine the effect of element size under shear loads. It is hypothesized that under shear loading the height of the element would become important as the width is important for mode I loading. • Further investigation is necessary to determine ideal values of the factor a in the element width equation in regime I. 131 Appendix A Load versus CMOD Plots APPENDIX A L O A D VERSUS C M O D PLOTS 132 Appendix A Load versus CMOD Plots 16 n 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 CMOD (mm) Figure A.l Experimental load-CMOD results and numerical prediction for specimen rd2 16 i 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 CMOD (mm) Figure A.2 Experimental load-CMOD results and numerical prediction for specimen rd8 133 Appendix A Load versus CMOD Plots 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 CMOD (mm) Figure A.3 Experimental load-CMOD results and numerical prediction for specimen rd21 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 CMOD (mm) Figure A.4 Experimental load-CMOD results and numerical prediction for specimen rd28 134 Appendix A Load versus CMOD Plots 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 CMOD (mm) Figure A. 5 Experimental load-CMOD results and numerical prediction for specimen rd32 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 CMOD (mm) Figure A. 6 Experimental load-CMOD results and numerical prediction for specimen rd38 135 Appendix A Load versus CMOD Plots 20 n 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 C M O D (mm) Figure A. 7 Experimental load-CMOD results and numerical prediction for specimen rd54 16 i 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 CMOD (mm) Figure A. 8 Experimental load-CMOD results for specimen rdl 136 Appendix A Load versus CMOD Plots -a a o rd4-l rd4-2 Figure A. 9 Experimental load-CMOD results for specimen rd4 73 a o — rd6-l — rd6-2 — rd6-3 Figure A. 10 Experimental load-CMOD results for specimen rd6 137 Appendix A Load versus CMOD Plots 16 i 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 CMOD (mm) Figure A. 11 Experimental load-CMOD results for specimen rdlO 16 , 0.0 0.5 1.0 1.5 . 2.0 2.5 3.0 3.5 4.0 4.5 5.0 CMOD (mm) Figure A. 12 Experimental load-CMOD results for specimen rdl 7 138 Appendix A Load versus CMOD Plots 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 C M O D (mm) Figure A. 13 Experimental load-CMOD results for specimen rd25 C M O D (mm) Figure A. 14 Experimental load-CMOD results for specimen rd35 139 Appendix A Load versus CMOD Plots C M O D (mm) Figure A. 15 Experimental load-CMOD results for specimen rd44 140 Appendix B Extra Line Analysis Plots APPENDIX B E X T R A LINE ANALYSIS PLOTS 141 Appendix B Extra Line Analysis Plots Figure B.l COD profiles determinedfrom line analysis for rd2-l 2.2 i Distance from initial notch tip (mm) Figure B.2 COD profiles determined from line analysis for rd2-2 142 Appendix B Extra Line Analysis Plots 1 .6 n Distance from initial notch tip (mm) Figure B.3 COD profiles determined from line analysis for rd8-2 2 i Distance from intial notch tip (mm) Figure B.4 COD profiles determined from line analysis for rd54-l 143 Appendix B Extra Line Analysis Plots 144 Appendix C Stress Distributions and Compliance Plots APPENDIX C STRESS DISTRIBUTION AND C O M P L I A N C E PLOTS 145 Appendix C Stress Distributions and Compliance Plots 0 5 10 15 20 25 30 35 40 45 JC (mm) Figure C1 Progression of the stress distribution along the notch plane for rd8. Each distribution is at the CMOD value (in mm) noted above it. 0 5 10 15 20 25 30 35 40 x (mm) Figure C.2 Progression of the stress distribution along the notch plane for rd21. Each distribution is at the CMOD value (in mm) noted above it. 146 Appendix C Stress Distributions and Compliance Plots 500 400 ^ 3 0 0 PH 200 100 u o Z, 100 -200 1.76 / 2.01 2.11 2.22 2.51 9 7 Q 3 8 4 4.36 3.05 3.32 3.58 | 4.10 1 4 6 2 10 15 20 25 x (mm) 30 35 40 45 Figure C.3 Progression of the stress distribution along the notch plane for rd32. Each distribution is at the CMOD value (in mm) noted above it. 0.0020 0.0018 H 0.0016 0.0014 | 0.0012 | 0.0010 i .2 a 0.0008 | CJ 0.0006 0.0004 -I 0.0002 0.0000 0 10 15 20 25 A a (mm) 30 35 40 Figure C.4 Compliance plots from location of peak stress in numerical simulations. 3.25 mm is subtracted from Aa to account for width of damage zone. 1 4 7 Appendix D ASDM UMATfor ABAQUS A P P E N D I X D A S D M U M A T F O R A B A Q U S c c UBC Adaptive Simple Damage Model f o r Abaqus v 4.3 c c E l a s t i c - l i n e a r s o f t e n i n g composite damage model c c Scott McClennan, J u l y 2004 c c c (c) Copyright The U n i v e r s i t y of B r i t i s h Columbia, Composites Group c 2004. A l l r i g h t s reserved. c c This o r t h o t r o p i c m a t e r i a l model i s f o r use with second-order plane c s t r e s s elements. I t i s a cohesive zone model r e q u i r i n g the o r t h o t r o p i c c e l a s t i c m a t e r i a l constants and two a d d i t i o n a l damage parameters. I t i s c intended to model mode I crack propagation i n laminate composites. The c model has the b u i l t - i n a b i l i t y to adapt i t s c o n s t i t u t i v e behaviour c based on mesh s i z e and the s t r e s s - f i e l d to allow a coarser mesh than i s c u s u a l l y r equired by a cohesive zone model. c Q* ***************************************************** ***************** C USER SUBROUTINE TO INITIALIZE STATE VARIABLES Q* ********************************************************************* * SUBROUTINE SDVINI(STATEV,COORDS,NSTATV,NCRDS,NOEL, NPT, 1 LAYER,KSPT) c INCLUDE 'ABA_PARAM.INC' c DIMENSION STATEV(NSTATV),COORDS(NCRDS) C C ******* STATE VARIABLES USED *************** c STATEV(1)=currpeakepsnpt(1) c STATEV(2)=currpeakepsnpt(2) c STATEV(3)=currpeakepsnpt(3) c STATEV(4)=Erl(1) c STATEV(5)=Erl(2) c STATEV(6)=Erl(3) c STATEV(7)=gammadis(1) c STATEV(8)=gammadis(2) c STATEV(9)=gammadis(3) c STATEV(10)=peakepsnpt(1) c STATEV(11)=peakepsnpt(2) c STATEV(12)=peakepsnpt(3) c STATEV(13)=ultepsnpt(1) c STATEV(14)=ultepsnpt(2) c STATEV(15)=ultepsnpt(3) c STATEV(16)=scale(1) c STATEV(17)=scale(2) c Q *******s^ATE VARIABLES************** do 11 i=l,17 statev(i)=0.0 11 continue 148 Appendix D ASDM UMATfor ABAQUS c RETURN END c c 0 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C STRAINFIELD SUBROUTINE 0 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * c t h i s subroutine c a l c u l a t e s the polynomial c o e f f i c i e n t s f o r the c s t r e s s d i s t r i b u t i o n matrix and stores them i n array s t f i e l d ( 6 , 2 , 3 ) c subroutine strainfield(polycoeffs,global,noelmax,nptmax,NOEL) c double p r e c i s i o n p o l y c o e f f s , g l o b a l i n t e g e r polynpts(6,3) in t e g e r i , j , k c dimension polycoeffs(6,2,3),global(noelmax,nptmax,2) c c 0 ********** POLYNOMIAL NPT LOOKUP ******************* c Each i n t e g r a t i o n p o i n t i s ass o c i a t e d with a polynomial i n both the c x - d i r e c t i o n and y - d i r e c t i o n c c **** SHELLS **** 0 **** x—DIR **** j = l do 21 i = l , 3 do 22 k=l,3 p o l y n p t s ( i , k ) = j j=j+l 22 continue 21 continue 0 **** Y—DIR **** polynpts(4,1)=1 polynpts(4,2)=4 polynpts(4,3)=7 polynpts(5,1)=2 polynpts(5,2)=5 polynpts(5,3)=8 polynpts(6,1)=3 polynpts(6,2)=6 polynpts(6,3)=9 c c ************** CALCULATE POLYNOMIAL COEFFS *********** c S i m i l a r to the element shape f u n c t i o n s , the 2nd order polynomial c c o e f f i c i e n t s are c a l c u l a t e d based on the p o s i t i o n of the i n t e g r a t i o n c p o i n t s c do 61 i = l , 6 do 62 j=l,2 polycoeffs(i,j,3)=GLOBAL(NOEL,polynpts(i,2),j) polycoeffs(i,j,2)=(GLOBAL(NOEL,polynpts(i,3),j) 1 -GLOBAL(NOEL,polynpts(i,1),j))/(-2.0*Sqrt(0.6) ) polycoeffs(i,j,1)=(GLOBAL(NOEL,polynpts(i,1),j)-Sqrt(0.6) 1 * p o l y c o e f f s ( i , j , 2 ) - p o l y c o e f f s ( i , j , 3 ) ) / 0 . 6 c 149 Appendix D ASDM UMATfor ABAQUS 62 continue 61 continue c r e t u r n end c c 0 ******************************************************** c SCALESTRAIN SUBROUTINE Q * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * c t h i s subroutine c a l c u l a t e s the maximum value of s t r a i n on the c polynomial curve a s s o c i a t e d with each i n t e g r a t i o n p o i n t f o r each c p r i n c i p a l d i r e c t i o n c subroutine scalestrain(maxst,polycoeffs,global,orient,noelmax,NOEL,NPT) c c polynomial c o e f f i c i e n t s double p r e c i s i o n p o l y c o e f f s , g l o b a l , o r i e n t c current s t r a i n s double p r e c i s i o n s t r a i n c i n t e g r a t i o n p o i n t lookup t a b l e double p r e c i s i o n nptlup(9,2) c parametric coordinates double p r e c i s i o n s i , s2 s3, s4, s5, s6, s7 c s c a l i n g s t r a i n array ( p o l y d i r e c t i o n , s t r a i n d i r e c t i o n ) double p r e c i s i o n maxst(2) c l o c a l v a r i a b l e s i n t e g e r i,j,k,m,n double p r e c i s i o n kl,k2 c dimension polycoeffs(6,2,3),global(noelmax,9,2), orient(noelmax,3) c c ************ NPT & EDGE COORDINATES ****************** s7=-1.0 s6=-Sqrt(0.6) s5=-4.0/9.0 s4=0.0 s3=4.0/9.0 s2=Sqrt(0.6) sl=1.0 c Q **************** LOOKUP TABLE ************************ c nptlup(9,2) -> values are p o l y i d x , p o l y i d y f o r each npt c c c *** poly 1,2,3 *** k=l do 110 i=l,7,3 do 111 j=0,2 nptlup(i+j,1)=k 111 continue k=k+l 110 continue c c *** poly 3,4,5 *** do 120 i=l,7,3 do 121 j=0,2 150 Appendix D ASDM UMATfor ABAQUS nptlup(i+j,2)=j+4 121 continue 120 continue c c c The o r i e n t a t i o n of the element r e l a t i v e to the g l o b a l coordinates i s c determined so that the maximum normal s t r a i n s f o r the polynomial i n c each the x- and y - d i r e c t i o n s can be determined f o r each i n t e g r a t i o n c p o i n t c c maxst(l)= max x s t r a i n on y poly c maxst(2)= max y s t r a i n on x poly c c ***** x - d i r max s t r e s s ***** c i f ( a b s ( o r i e n t ( n o e l , 3 ) - o r i e n t ( n o e l , 1 ) ) . g t . 1 a b s ( o r i e n t ( n o e l , 2 ) - o r i e n t ( n o e l , 1 ) ) ) then m=2 n=l e l s e m=l n=2 endif . kl=0.0 k2=0.0 i f (abs(polycoeffs(nptlup(npt,m),1,3)-1 g l o b a l (noel,npt,1)). l t . (0.00001*global(noel,npt,1))) then kl=s3 k2=s5 e l s e i f (abs(polycoeffs(nptlup(npt,m),1,1)*s2**2.0+ 1 polycoeffs(nptlup(npt,m),1,2)*s2+polycoeffs(nptlup(npt,m),1,3)-1 global(noel,npt,1)) . I t . (0.0001*global (noel,npt,1))) then k l = s l k2=s3 el s e kl=s5 k2=s7 endif c maxst(1)=max(0.0, g l o b a l ( n o e l , n p t , 1 ) , 1 polycoeffs(nptlup(npt,m),1,1)*kl**2.0+ 1 polycoeffs(nptlup(npt,m),1,2)*kl+polycoeffs(nptlup(npt,m),1,3), 1 polycoeffs(nptlup(npt,m),1,1)*k2**2.0+ 1 polycoeffs(nptlup(npt,m),1,2)*k2+polycoeffs(nptlup(npt,m),1,3)) c c ***** y - d i r max s t r e s s ***** c kl=0.0 k2 = 0 . 0 i f ( a b s ( p o l y c o e f f s ( n p t l u p ( n p t , n ) , 2 , 3 ) -1 g l o b a l ( n o e l , n p t , 2 ) ) . I t . ( 0 . 0 0 0 1 * g l o b a l ( n o e l , n p t , 2 ) ) ) then kl=s3 k2=s5 e l s e i f (abs(polycoeffs(nptlup(npt,n),2,l)*s2**2.0+ 1 polycoeffs(nptlup(npt,n),2,2)*s2+polycoeffs(nptlup(npt,n),2,3)-1 g l o b a l ( n o e l , n p t , 2 ) ) . l t . ( 0 . 0 0 0 1 * g l o b a l ( n o e l , n p t , 2 ) ) ) then k l = s l 151 Appendix D ASDM UMAT for ABAQUS k2=s3 el s e kl=s5 k2=s7 endif c maxst(2)=max(0.0, g l o b a l ( n o e l , n p t , 2 ) , 1 polycoeffs(nptlup(npt,n),2,1)*kl**2.0+ 1 p o l y c o e f f s ( n p t l u p ( n p t , n ) , 2 , 2 ) * k l + p o l y c o e f f s ( n p t l u p ( n p t , n ) , 2 , 3 ) , 1 polycoeffs(nptlup(npt,n),2,1)*k2**2.0+ 1 polycoeffs(nptlup(npt,n),2,2)*k2+polycoeffs(nptlup(npt,n),2,3)) c c re t u r n end c ********************************************************** c The f o l l o w i n g u s e r - m a t e r i a l subroutine determines the damaging m a t e r i a l c response of an o r t h o t r o p i c m a t e r i a l . Q ******************* *.* ********************************************* C USER SUBROUTINE FOR MATERIALS Q * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE, SPD,SCD, 1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN, 2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI,NSHR,NTENS, 3 NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT, 4 DFGRDO,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC) c INCLUDE 'ABA_PARAM.INC c parameter noelmax=700 parameter nptmax=9 c CHARACTER*80 MATERL C DIMENSION STRESS(NTENS),STATEV(NSTATV), 1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS), 2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1), 3 PROPS(NPROPS),COORDS(3),DROT(3,3), 4 DFGRDO(3,3),DFGRD1(3,3) c c m a t e r i a l constants double p r e c i s i o n E0(3) double p r e c i s i o n nxy, nyx c compliance tensor terms double p r e c i s i o n c, sn(2) c model parameters double p r e c i s ion peakepsch(3) double p r e c i s ion gammac(2) c determined parame t e r s double p r e c i s ion STRANcurr(3) double p r e c i s ion Et(3) double p r e c i s ion E r l ( 3 ) double p r e c i s ion peakepsnpt(3) double p r e c i s ion currpeakepsnpt double p r e c i s ion ultepsch(3) double p r e c i s ion ultepsnpt(3) 152 Appendix D ASDM UMATfor ABAQUS d o u b l e p r e c i s i o n g a m m a r l(3) d o u b l e p r e c i s i o n g a m m a d i s(3) d o u b l e p r e c i s i o n f s d o u b l e p r e c i s i o n d f s d o u b l e p r e c i s i o n f s c u r r d o u b l e p r e c i s i o n s i g m a ( 3 ) c s t r a i n d i s t r i b t i o n m a t r i x c o n t a i n i n g p o l y n o m i a l c o e f f s d o u b l e p r e c i s i o n s t f i e l d ( 6 , 2 , 3 ) c l o c a l s t r a i n s t o r a g e m a t r i x f o r c u r r e n t n p t d o u b l e p r e c i s i o n l o s t r a i n ( 9 , 2 ) c g l o b a l s t r a i n m a t r i x f o r a l l n p t d o u b l e p r e c i s i o n g l s t r a i n ( n o e l m a x , n p t m a x , 2 ) c m a x i m u m e l e m e n t a l s t r a i n s ( x , y ) d o u b l e p r e c i s i o n m a x e p s(2) c l o c a l o r i e n t a t i o n c o o r d i n a t e s s t o r a g e d o u b l e p r e c i s i o n l o c o o r d s ( 3 ) c g l o b a l o r i e n t a t i o n c o o r d i n a t e s s t o r a g e d o u b l e p r e c i s i o n g l c o o r d s ( n o e l m a x , 3 ) c l o c a l v a r i a b l e s a n d f l a g s c DOUBLE PRECISION d e b u g i n t e g e r i , j , k i n t e g e r s t e p i n t e g e r s c a l e (2) c Q * * * * * * * * * * USER INPUT * * * * * * * * * * * * * * * * * * * * * * * * * * * c EO (1) = p r o p s (1) ~-EO(2) = p r o p s ( 9 ) n x y = p r o p s (2) EO(3) = p r o p s ( 3 ) g a m m a c ( l ) = p r o p s (4) gammac(2) = p r o p s ( 1 2 ) p e a k e p s c h ( l ) = p r o p s ( 5 ) p e a k e p s c h ( 2 ) = p r o p s ( 1 3 ) c n y x = n x y n y x = n x y * EO (2)/EO(1) c c c * * * * * * * * * * CALCULATE CHARACTERISTIC ULTIMATE STRAIN * * * * * * * * * * C d o 210 i = 1,NDI u l t e p s c h ( i ) = ( 2 . 0 * g a m m a c ( i ) ) / ( E O ( i ) * p e a k e p s c h ( i ) ) i f ( u l t e p s c h ( i ) . l t . p e a k e p s c h ( i ) ) t h e n w r i t e (*,*) 'WARNING: u l t e p s ( ' , i , ') < p e a k e p s ( ' , i , ') ' e n d i f 210 c o n t i n u e c p e a k e p s c h ( 3 ) = m i n ( p e a k e p s c h ( 1 ) , p e a k e p s c h ( 2 ) ) u l t e p s c h ( 3 ) = m a x ( u l t e p s c h ( 1 ) , u l t e p s c h ( 2 ) ) c c g a m m a c ( 3 ) = 0 . 5 * E 0 ( 3 ) * p e a k e p s c h ( 3 ) * u l t e p s c h (3) c c w r i t e ( 6 , * ) ' e n d c a l c u l a t e u l t i m a t e s t r a i n ' c Q ********** CURRENT STRAINS * * * * * * * * * * * * c 153 Appendix D ASDM UMATfor ABAQUS do i = l , 3 STRANcurr(i)=STRAN(i)+DSTRAN(i) enddo fs=abs(STRAN(3)) dfs=abs(DSTRAN(3)) fscurr=fs+dfs c Q ********** GET STATE VARIABLES *********************** c do 260 i = l , 3 if(STRAN(i).eq.0.0) then E r l ( i ) = E 0 ( i ) gammarl(i)=gammac(i) currpeakepsnpt(i)=peakepsch(i) peakepsnpt(i)=peakepsch(i) u l t e p s n p t ( i ) = u l t e p s c h (i) e l s e c urrpeakepsnpt(i)=statev(i) E r l ( i ) = s t a t e v ( i + 3 ) gammadis(i)=statev(i+6) peakepsnpt(i)=statev(i+9) ultepsnpt(i)=statev(i+12) e n d i f 260 continue scale(1)=statev(16) scale(2)=statev(17) c c write(6,*) 'end s t a t e v a r i a b l e s ' c ******* AUTOMATIC TIME INCREMENTATION **************** c c PNEWDT=0.5 c c ********* DETERMINE POLYNOMIAL COEFFICIENTS ********** C c a l l s t r a i n f i e l d ( s t f i e l d , g i s t r a i n , n o e l m a x , n p t m a x , N O E L ) c c C ********* CALCULATE SCALING STRAIN VALUE ************* c c a l l scalestrain(maxeps,stfield,glstrain,glcoords,noelmax,NOEL,NPT) c c write(6,*) ' s c a l e s t r a i n subroutine f i n i s h e d ' c Q * * * * * * * * * * * * * T O SCALE OR NOT TO SCALE? *************** c Only scale once during the e l a s t i c l oading p o r t i o n of the s t r e s s - s t r a i n c curve c do 270 i=2,2 i f (maxeps(i) .ge.0.0.and.glstrain(NOEL, NPT, i ) .gt.0 . 0 . 1 and.(STRANcurr(i).gt.(0.3*peakepsnpt(i))). 1 and. (STRANcurr(i) . I t . (0.5*peakepsnpt ( i ) ) ) . 1 and.(DSTRAN(i).gt.0.0). 1 and. (scale (i) . ne . 1)) then c s c a l e ( i ) = 1 write(6,*) ' s c a l i n g noel, npt, d i r ' , NOEL, NPT, i write(6,*) ' g l s t r a i n , STRANcurr', glstrain(NOEL,NPT,i), STRANcurr(i) 154 Appendix D ASDM UMATfor ABAQUS i f ( m a x e p s ( i ) . I t . g l s t r a i n ( n o e l , n p t , i ) ) then w r i t e (6,*) 'WARNING: maxeps(',i,') i s l e s s than g l s t r a i n (', noel, ', ',npt, ', ',i, ' ) ' e l s e peakepsnpt(i)=min(peakepsch(i) , 1 p e a k e p s c h ( i ) * ( g l s t r a i n ( n o e l , n p t , i ) / m a x e p s ( i ) ) ) endif i f (STRANcurr(i).gt.peakepsnpt(i)) then w r i t e ( * , * ) 'current s t r a i n i s l a r g e r than scaled peak s t r a i n f o r 1 element, npt, d i r : ' , noel, npt, i w r i t e (*,*) 'DSTRAN=', DSTRAN(i) w r i t e (*,*) 'STRANcurr:', STRANcurr(i) w r i t e (*,*) 'peakepsnpt:', peakepsnpt (i) peakepsnpt(i)=STRANcurr(i) endif u l t e p s n p t ( i ) = (2.0*gammarl(i))/(Erl(i)*peakepsnpt ( i ) ) currpeakepsnpt(i)=peakepsnpt (i) e l s e currpeakepsnpt(i)=max(STRANcurr(i),currpeakepsnpt(i),peakepsnpt(i)) endif 270 continue currpeakepsnpt(3)=max(fscurr,currpeakepsnpt(3),peakepsnpt(3)) c c c w r i t e (6,*) 'end to s c a l e or not' c c Q ***************** p o i S S O N REDUCTIONS ****************** C c s n ( l ) = s n l 2 , sn(2) = sn21 c write(6,*) 'Update Poissons...' s n ( l ) = (Erl(1)/E0(1))*nxy sn(2) = (Erl(2)/E0(2))*nyx c c c = (1. - sn ( 2 )*sn(1)) i f (c . l e . 0.0) then w r i t e ( * , * ) 'Warning: s i n g u l a r s t i f f n e s s tensor' w r i t e ( * , * ) s n ( l ) , sn(2) endif c c write(6,*) 'end Poisson reductions' c ********** DETERMINE UPDATED STRESSES **************** C sigma=0.0 c c do 310 i = l , 2 i f (STRANcurr(i).ge.ultepsnpt(i)) then sigma(i)=0.0 e l s e i f (STRANcurr(i) .ge.currpeakepsnpt ( i ) ) then sigma (i) = ( (ultepsnpt (i) -^STRANcurr (i) ) / (ultepsnpt (i) -1 p e a k e p s n p t ( i ) ) ) * ( E 0 ( i ) * p e a k e p s n p t ( i ) ) endif 310 continue c c write(6,*) 'End s t r e s s update' 155 Appendix D ASDM UMAT for ABAQUS c Q ************* DETERMINE TANGENT MODULI ****************** C i f (STRANcurr(1).ge.ultepsnpt(1)) then Et(l)=0.0 e l s e i f (STRANcurr(1).It.currpeakepsnpt(1).or.DSTRAN(1).le.0.0) then Et(1) = E r l (1) else Et (1) = (sigma (1) -STRESS (1) ) / ( (DSTRAN (1) +sn (2) * DSTRAN (2) ) /c) endif c i f (STRANcurr(2).ge.ultepsnpt(2)) then Et(2)=0.0 e l s e i f (STRANcurr(2).It.currpeakepsnpt(2).or.DSTRAN(2).le.0.0) then Et(2) = E r l ( 2 ) e l s e Et(2)=(sigma(2)-STRESS(2)-(sn(2)*Et (1 )/c)*DSTRAN(1))/(DSTRAN(2)/c) endif c c Et(3)=E0 (3) c c write(6,*) 'End determine tangent moduli' c c ***********ASSEMBLE TANGENT STIFFNESS TENSOR********** C c write(6,*) 'Assemble s t i f f n e s s tensor...' do 430 kl=l,NTENS do 431 k2=l,NTENS DDSDDE(kl,k2)=0.00 431 continue 430 continue c i f (c . l e . 0.0) then w r i t e (*,*) 'Warning: s i n g u l a r s t i f f n e s s tensor' w r i t e ( * , * ) s n ( l ) , sn(2) el s e DDSDDE(1,1) = Et (1) / C DDSDDE(1,2) = (sn(2) * E t ( l ) ) / c DDSDDE(2,1) = DDSDDE(1,2) DDSDDE(2,2) = Et (2) / c DDSDDE (3,3) = Et(3) DDSDDE(1,3) =0.0 DDSDDE (2,3) =0.0 DDSDDE (3,1) =0.0 DDSDDE(3,2) =0.0 endif c c write(6,*) 'end-tangent s t i f f n e s s tensor assembly' c Q ************** STRESS CALCULATION********************** c C do 440 i=l,NTENS do 441 j=l,NTENS STRESS(i)=STRESS(i)+DDSDDE(i,j)*DSTRAN(j) 441 continue 156 Appendix D ASDM UMATfor ABAQUS 440 continue c do 444 i=l,2 i f (STRESS(i) .gt. ( E r l ( i ) * c u r r p e a k e p s n p t ( i ) ) ) then STRESS(i) = E r l ( i ) * S T R A N c u r r ( i ) endif 444 continue c c w r i t e (6,*) 'end s t r e s s c a l c u l a t i o n ' c ****** ELASTIC STRAIN ENERGY CALCULATION ************** C do 450 i=l,NTENS SSE=SSE+STRESS(i)*DSTRAN(i)12 .0 450 continue c c ******* DISSIPATED STRAIN ENERGY CALCULATION ************ c do 455 i=l,2 if(STRAN(i) .ge.peakepsnpt ( i ) ) then gammadis(i)=min(gammac(i), 1 max(gammadis(i), (0.5*E0(i)*peakepsnpt ( i ) * * 2 . 0) + 1 (((E0(i)*peakepsnpt(i))+STRESS(i)) 1 2 . 0)*(STRANcurr(i) 1 peakepsnpt(i))-(0.5*STRESS(i)*STRANcurr(i)))) endif 455 continue c i f ( f s . g t . p e a k e p s n p t (3)) then gammadis(3)=min(gammac(3), 1 max(gammadis(3), (0.5*E0(3)*peakepsnpt(3)**2.0)+ 1 ( ( (E0 (3)*peakepsnpt(3))+abs(STRESS(3)))12 . 0)* 1 (abs(STRANcurr(3))-peakepsnpt ( 3 ) ) -1 (0.5*STRESS(3)*STRANcurr(3)))) endif c c w r i t e (6,*) 'end s t r a i n energy c a l c u l a t i o n ' c Q **************** UPDATE PARAMETERS ******************* c c ** c h a r a c t e r i s t i c parameters ** c c w r i t e (6,*) 'update parameters' do 460 i = l , 3 i f ( c u r r p e a k e p s n p t ( i ) . g e . u l t e p s n p t ( i ) ) then gammarl (i)=0 . 0 el s e gammarl(i)=min(gammac(i),gammac(i)-gammadis(i)) endif 4 60 continue c c ** element parameters ** c do 470 i=l , 2 E r l ( i ) = m i n ( E r l ( i ) , ( 2 . 0 * g a m m a r l ( i ) ) / ( u l t e p s n p t ( i ) * c u r r p e a k e p s n p t ( i ) ) ) 470 continue c E r l ( 3 ) = m i n ( E r l ( 3 ) , (2 . 0*gammarl (3) ) /(ultepsnpt(3)*currpeakepsnpt (3))) c 157 Appendix D ASDM UMATfor ABAQUS c c c c c c c c c w r i t e ( 6 , * ) ' e n d u p d a t e p a r a m e t e r s ' * * * * * * * * * * * * * U P D A T E S T A T E V A R I A B L E S w r i t e ( 6 , * ) ' u p d a t e s t a t e v a r i a b l e s ' S T A T E V ( 1 ) = c u r r p e a k e p s n p t ( 1 ) S T A T E V ( 2 ) = c u r r p e a k e p s n p t ( 2 ) S T A T E V ( 3 ) = c u r r p e a k e p s n p t ( 3 ) S T A T E V ( 4 ) = E r l ( 1 ) S T A T E V ( 5 ) = E r l (2) S T A T E V ( 6 ) = E r l ( 3 ) S T A T E V ( 7 ) = g a m m a d i s ( 1 ) S T A T E V ( 8 ) = g a m m a d i s ( 2 ) S T A T E V ( 9 ) = g a m m a d i s ( 3 ) S T A T E V ( 1 0 ) = p e a k e p s n p t ( 1 ) S T A T E V ( 1 1 ) = p e a k e p s n p t ( 2 ) S T A T E V ( 1 2 ) = p e a k e p s n p t ( 3 ) S T A T E V ( 1 3 ) = u l t e p s n p t ( 1 ) S T A T E V ( 1 4 ) = u l t e p s n p t ( 2 ) S T A T E V ( 1 5 ) = u l t e p s n p t ( 3 ) S T A T E V ( 1 6 ) = s c a l e ( 1 ) S T A T E V ( 1 7 ) = s c a l e (2) * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * D E B U G G I N G * * * * * * * * * * * * * * * * * * * * * * * s t e p = s t e p + l ( s t e p e q w r i t e 6, *) w r i t e 6, *) w r i t e 6, *) w r i t e 6, *) w r i t e 6, *) w r i t e 6, *) w r i t e (6 , *) w r i t e (6 , *) w r i t e (6 , *) w r i t e (6 , *) w r i t e (6 , *) w r i t e (6 , *) w r i t e (6 , *) w r i t e (6 , *) w r i t e (6 , *) w r i t e (6 , *) w r i t e (6 , *) w r i t e (6 , *) w r i t e (6 , *) w r i t e (6 , *) w r i t e (6 , *) w r i t e (6 , *) w r i t e (6 , *) w r i t e (6 , *) 1) t h e n ' e l n u m b e r : ' , N O E L , ' , N P T ' , N P T ' c o o r d s x , y : ' , COORDS ( 1 ) , COORDS (2) ' l o c o o r d s 1 , 2 , 3 : ' , l o c o o r d s ( l ) , l o c o o r d s ( 2 ) , l o c o o r d s ( 3 ) ' g l c o o r d s 1 , 2 , 3 : ' , g l c o o r d s ( n o e l , 1 ) , g l c o o r d s ( n o e l , 2 ) , g l c o o r d s ( n o e l , 3 ) ' p e a k e p s c h ( 1 , 2 , 3 ) ' , p e a k e p s c h ( 1 ) , p e a k e p s c h ( 2 ) , p e a k e p s c h ( 3 ) ' p e a k e p s n p t ( 1 , 2 , 3 ) ' , p e a k e p s n p t ( 1 ) , p e a k e p s n p t ( 2 ) , p e a k e p s n p t ( 3 ) ' c u r r p e a k e p s n p t ( 1 , 2 , 3 ) ' , c u r r p e a k e p s n p t ( 1 ) , c u r r p e a k e p s n p t ( 2 ) , c u r r p e a k e p s n p t ( 3 ) ' u l t e p s c h ( 1 , 2 , 3 ) ' , u l t e p s c h ( 1 ) , u l t e p s c h ( 2 ) , u l t e p s c h ( 3 ) ' u l t e p s n p t ( 1 , 2 , 3 ) ' , u l t e p s n p t ( 1 ) , u l t e p s n p t ( 2 ) , u l t e p s n p t ( 3 ) ' g a m m a c ( 1 , 2 , 3 ) ' , g a m m a c ( 1 ) , g a m m a c ( 2 ) , g a m m a c ( 3 ) ' g a m m a c ( 1 , 2 , 3 ) ' , g a m m a c ( 1 ) , g a m m a c ( 2 ) , g a m m a c ( 3 ) ' g a m m a d i s ( 1 , 2 , 3 ) ' , g a m m a d i s ( 1 ) , g a m m a d i s ( 2 ) , g a m m a d i s ( 3 ) ' g a m m a r l ( 1 , 2 , 3 ) ' , g a m m a r l ( 1 ) , g a m m a r l ( 2 ) , g a m m a r l ( 3 ) ' S T R A N ( 1 , 2 , 3 ) ' , S T R A N ( 1 ) , S T R A N ( 2 ) , S T R A N ( 3 ) ' S T R A N c u r r ( 1 , 2 , 3 ) ' , S T R A N c u r r ( 1 ) , S T R A N c u r r ( 2 ) , S T R A N c u r r (3) ' D S T R A N ( 1 , 2 , 3 ) ' , D S T R A N ( 1 ) , D S T R A N ( 2 ) , D S T R A N ( 3 ) ' S T R E S S ( 1 , 2 , 3 ) ' , S T R E S S ( 1 ) , S T R E S S ( 2 ) , S T R E S S ( 3 ) ' f s ' , f s ' s i g m a ( 1 , 2 , 3 ) : ' , s i g m a ( l ) , s i g m a ( 2 ) , s i g m a ( 3 ) ' s n ( l ) , s n ( 2 ) : ' s n ( l ) , s n ( 2 ) ' E r l ( 1 , 2 , 3 ) ' , E r l ( l ) , E r l ( 2 ) , E r l ( 3 ) ' E t ( l , 2 , 3 ) : ' , E t ( l ) , E t ( 2 ) , E t ( 3 ) ' D D S D D E ( 1 , 1 ) , ( 1 , 2 ) : ' , D D S D D E ( 1 , 1 ) , D D S D D E ( 1 , 2 ) ' D D S D D E ( 2 , 1 ) , ( 2 , 2 ) : ' , D D S D D E ( 2 , 1 ) , D D S D D E ( 2 , 2 ) 158 Appendix D ASDM UMAT for ABAQUS w r i t e ( 6 , * ) 'DDSDDE( 3 , 3) : ', D D S D D E ( 3 , 3 ) w r i t e ( 6 , * ) ' maxeps ( 1 , 2 ) ' . , maxeps ( 1 ) , maxeps ( 2 ) w r i t e ( 6 , * ) ' g l s t r a i n ( n o e l , n p t , 1 ) ', g l s t r a i n ( n o e l , n p t , 1 ) w r i t e ( 6 , * ) ' g l s t r a i n ( n o e l , n p t , 2 ) ' , g l s t r a i n ( n o e l , n p t , 2 ) w r i t e ( 6 , * ) ' s t f i e l d ( 1 , 2 , 1 ) ', s t f i e l d ( 1 , 2 , 1 ) w r i t e ( 6 , * ) ' s t f i e l d ( 1 , 2 , 2 ) ', s t f i e l d ( 1 , 2 , 2 ) w r i t e ( 6 , * ) ' s t f i e l d ( l , 2 , 3 ) ', s t f i e l d ( 1 , 2 , 3 ) w r i t e ( 6 , * ) ' s t e p ' , s t e p s t e p = 0 e n d i f c c *********** MATERIAL MODEL ENDS HERE ***************************** C Q *********** LOCAL STRAIN STORAGE ********************* c s t o r e t h e s t r a i n s i n x , y d i r e c t i o n f o r c u r r e n t NPT i n t h e c t e m p o r a r y s t o r a g e a r r a y l o s t r a i n ( n p t m a x , 3 ) c do 6 1 0 i = l , 2 l o s t r a i n ( n p t , i ) = S T R A N c u r r ( i ) 6 1 0 c o n t i n u e c i f ( n p t . e q . l ) t h e n l o c o o r d s ( 1 ) = c o o r d s ( 1 ) e l s e i f ( n p t . e q . 7 ) t h e n l o c o o r d s ( 2 ) = c o o r d s ( 1 ) e l s e i f ( n p t . e q . 3 ) t h e n l o c o o r d s ( 3 ) = c o o r d s ( 1 ) e n d i f c c i f NPT = NPTMAX t h e n c o p y t h e l o c a l s t r a i n s t o t h e g l o b a l s t r a i n c a r r a y . c i f ( n p t . e q . n p t m a x ) t h e n do 6 2 0 i = l , n p t m a x do 6 2 1 j = l , 2 g l s t r a i n ( n o e l , i , j ) = l o s t r a i n ( i , j ) 6 2 1 c o n t i n u e 6 2 0 c o n t i n u e do 6 2 2 i = l , 3 g l c o o r d s ( n o e l , i ) = l o c o o r d s ( i ) 6 2 2 c o n t i n u e c l o s t r a i n = 0 . 0 l o c o o r d s = 0 . 0 c e n d i f c w r i t e ( 6 , * ) 'end o f m a t e r i a l l o o p ' c RETURN END 159 Appendix E Transition Radius and Peak Load Prediction Plot for a DCB Specimen APPENDIX E TRANSITION RADIUS AND P E A K L O A D PREDICTION P L O T FOR AN AS4/3501-6 DCB SPECIMEN E. l Specimen The double-canti lever beam specimen ( D C B ) modeled here is shown in Figure E. l w i t h dimensions. I t is modeled w i t h the mater ial properties for unid i rect ional AS4/3501-6 carbon fibre/epoxy g iven i n Table E . l . This is the same combinat ion o f specimen and material used by (Kan j i , 2003). E.2 Regime I Energy dr iven fai lure occurs for sharp notches. The energy fai lure cr i ter ion is developed in the same w a y as for the O C T specimen. Since the stress cr i ter ion is satisfied f i rst in regime I, the relat ionship between the cr i t ical stress intensity factor and the cr i t ical energy release rate is g iven by Equat ion 2.9 repeated here as KIC2=E'GJc. E . l The stress intensity factor for a D C B specimen is g iven in terms o f load and specimen dimensions by (Tada, Paris and I r w i n , 1985) as 2V3~(%)b K,=C y ^ - , E.2 h/2 where C = 1.0 plane stress 1 F T, plane strain v2 Combin ing Equat ion E. 1 and Equat ion E.2 yields the predic t ion for peak load in regime I , JE'GeBh' Pc = - T = • E - 4 ClSa For a g iven Gc, the peak load is invariant w i t h notch root radius in this regime. 160 Appendix E Transition Radius and Peak Load Prediction Plot for a DCB Specimen E.3 Regime II The stress fai lure cr i ter ion is used to predict the peak load in regime I I . A s was g iven in Equat ion 5.10, the sharp notch stress intensity factor can be related to the cr i t ical stress in a blunt notch through Equat ion 5.9. The cr i t ical stress intensity factor is then C o m b i n i n g the Ki-P relat ion from Equat ion E.2 w i t h Equat ion E.5, the peak load i n regime H i s I 3 / Pf = ^ H . E . 6 C4V3a A s w i t h the O C T specimen, the peak load is related to the cr i t ical stress and the square root o f the notch root radius. E.4 T r a n s i t i o n r a d i u s The transit ion radius occurs at the point when the peak load in regime I is equal to the peak load in regime I I . C o m b i n i n g Equat ion E.4 and Equat ion E.6, the transit ion radius is 4E'GC o = F 7 r transition 2 " ncjc Since ptransition is a material property, Equat ion E.7 for a D C B specimen is exact ly the same as Equat ion 5.13 for an O C T specimen. E .5 D i s c u s s i o n The peak load predict ion p lot is presented in Figure E.2. The peak load is predicted to be about 16.5 N for sharp notches and for this material the transi t ion radius is predicted to be 0.21 m m . I t w o u l d actual ly be quite d i f f i cu l t to obtain a notch radius o f this size in the delaminat ion specimen because it w o u l d mean that the resin layer between plies on the delaminat ion plane w o u l d have to be over 0.4 m m th ick at the notch t ip. The typ ica l cured p l y thickness for this material ( f ibre and resin) is 0.132 m m (Hexcel Corporat ion, 1998). 161 Appendix E Transition Radius and Peak Load Prediction Plot for a DCB Specimen The interface strength plays a large role in determining ptransition- I f a material system w i t h a stronger more br i t t le mat r ix is considered, ptransition w i l l greatly decrease. For example, i f crc is doubled to 200 M P a , ptransition becomes 0.05 m m . The size o f ptransition i n a D C B specimen is important because notches in D C B specimens are usual ly not perfect ly sharp. A nonadhesive insert such as Te f lon is usual ly inserted between plies dur ing processing to create the in i t ia l notch. A S T M provides guidel ines o f h o w th ick the insert can be wi thout af fect ing the results due to notch sensit iv i ty. I n A S T M D-5528 ( A S T M , 2002) the m a x i m u m thickness o f insert a l lowable is 0.013 m m . The transi t ion radius provides an indicat ion o f h o w th ick an insert can actual ly be w i thou t causing unstable fracture. For unid i rect ional AS4/3501-6 , the A S T M standard provides a safe value. For very strong, br i t t le interfaces, however, it becomes more important to use a very th in nonadhesive insert to create the in i t ia l notch. A theoretical br i t t le material w i t h an inter laminar strength o f 400 M P a and a strain energy release rate o f 60 J / m A 2 w o u l d have a transi t ion radius o f 0.0064 m m , meaning that a 0.013 m m nonadhesive insert w o u l d be th ick enough to make fracture unstable. 162 Appendix E Transition Radius and Peak Load Prediction Plot for a DCB Specimen T a b l e E . l AS4/3501-6 m a t e r i a l p rope r t i es P a r a m e t e r V a l u e Ei ' 102 Gpa' E2 9 Gpa1 Vl2 0.3' G,2 7.1 Gpa1 Gic 125 J / m 2 1 1 100 M P a 1 " ' f r o m ( K a n j i , 2003) I I f r o m Hexce l data sheets (Hexcel Corporat ion, 1998) I I I estimated f r o m values in l i terature, e.g. (Jackson and M a r t i n , 1993; L i f sh i t z and Leber, 1998; W i s n o m , Jones and H i l l , 2001) 163 Appendix E Transition Radius and Peak Load Prediction Plot for a DCB Specimen a P A h thicknessB, a»h Figure E.l DCB specimen (unidirectional fibres oriented along length of specimen) | Stress criterion Energy criterion \ Transition radius / I — 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 p (mm) 0.45 0.5 Figure E.2 Peak load prediction for AS4/3501-6 DCB specimen 164 References R E F E R E N C E S Al fano , G and Cr is f ie ld , M A . "Fin i te element interface models for the delaminat ion analysis o f laminated composites: mechanical and computat ional issues", Internat ional Journal for Numer ica l Methods in Engineer ing, V o l . 50, pp. 1701-1736, 2 0 0 1 . A S T M . D 6272: "Standard test method for f lexural properties o f unreinforced and re inforced plastics and electr ical insulat ing materials by four-po in t bend ing" , pp. 923-931 , 2000. A S T M . D 5528: "Standard Test M e t h o d for M o d e I Inter laminar Fracture Toughness o f Unid i rect ional F iber-Reinforced Polymer M a t r i x Composi tes", pp. 292-302, 2002. Awerbuch , J and Madhukar , M S . "Notched Strength o f Composi te Laminates: Predict ions and Exper iments—a Rev iew" , Journal o f Reinforced Plastics and Composites, V o l . 4, p p . 3 - 1 5 9 , 1985. Back lund , J and Aronsson, C-G. "Tensi le fracture o f laminates w i t h holes", Journal o f Composi te Mater ia ls, V o l . 20, pp. 259-286, 1986. Barenblatt , G I . "The mathematical theory o f equ i l ib r ium cracks in br i t t le f racture", Advances in A p p l i e d Mechanics, V o l . 7, pp. 55-129, 1962. Bazant, Z . "Size Ef fect in B lun t Fracture: Concrete, Rock, M e t a l " , Journal o f Engineer ing Mechanics, V o l . 110, pp. 518-535, 1984. Bazant, Z . Scal ing o f Structural Strength. London , U K : Hermes Penton L td . , 2002. Bazant, Z and O h , B H . "Crack band theory for fracture o f concrete", Mater iaux et Construct ions.Mater ials and Structures, V o l . 16, pp. 155-177, 1983. Bazant, ZP and Oh, B H . "Rock fracture v ia strain-softening f in i te elements", Journal o f Engineer ing Mechanics, V o l . 110, pp. 1015-1035, 1984. B o r g , R, N i lsson, L, and Simonsson, K. "S imula t ion o f delaminat ion in fiber composites w i t h a discrete cohesive fai lure mode l " , Composites Science and Technology, V o l . 6 1 , pp. 667-677, 2 0 0 1 . B o r g , R, N i lsson, L, and Simonsson, K. " M o d e l i n g o f delaminat ion using a discretized cohesive zone and damage fo rmu la t ion" , Composites Science and Technology, V o l . 62, pp. 1299-1344, 2002. Broek, D. Elementary Engineer ing Fracture Mechanics. The Hague: Mar t inus N i j h o f f Publishers, 1982. 165 References Camanho, PP, Dav i la , C G , and de M o u r a , M F . "Numer ica l s imulat ion o f mixed-mode progressive delaminat ion in composite materials", Journal o f Composi te Mater ia ls, V o l . I n press, 2003. Cook, R D , Ma lkus , D S , and Plesha, M E . Concepts and Appl icat ions o f F in i te Element Analys is . N e w Y o r k : John W i l e y & Sons, 1989. Cor ig l iano, A , Mar ian i , S, and Pandol f i , A . "Numer ica l model ing o f rate-dependent debonding processes in composites", Composi te Structures, V o l . 6 1 , pp. 39-50, 2003. Creager, M and Paris, PC. "Elast ic field equations for b lunt cracks w i t h reference to stress corosion crack ing" , International Journal o f Fracture Mechanics, V o l . 3, pp. 247-252, 1967. de Borst , R. "Numer ica l aspects o f cohesive-zone models" , Engineer ing Fracture Mechanics, V o l . 70, pp. 1743-1757, 2003. Dharan i , L R , Jones, W F , and Goree, JG. "Mathemat ica l mode l ing o f damage in unid i rect ional composites", Engineer ing Fracture Mechanics, V o l . 17, pp. 555-573, 1983. Dugdale, DS. " Y i e l d i n g o f steel sheets containing sl i ts", Journal o f Mechanics and Physics o f Sol ids, V o l . 8, pp. 100-104, 1960. Espinosa, H D , D w i v e n d i , S, and L u , H-C. " M o d e l i n g impact induced delaminat ion o f woven fiber re inforced composites w i t h contact/cohesive laws" , Computer Methods in A p p l i e d Mechanics and Engineer ing, V o l . 183, pp. 259-290, 2000. F loyd , A . " A n engineering approach to the s imulat ion o f gross damage development in composite laminates", Ph.D. Thesis. Department o f Metals and Mater ia ls Engineer ing, The Un ivers i ty o f Br i t i sh Co lumbia , 2004. Geubel le, P H and Bay lor , JS. " Impact- induced delaminat ion in composites: a 2 D s imulat ion" , Composites Part B, V o l . 29B , pp. 589-602, 1998. G r i f f i t h , A A . "The phenomena o f rupture and f l o w in sol ids", Phi losophical Transactions o f the Roya l Society o f London. V o l . A 2 2 1 , pp. 163-197, 1921. Hexce l Corporat ion. "Hexcel Composites Product Data Sheet", 1998. H i l le rborg , A , Modeer , M , and Petersson, P-E. "Analys is o f Crack Format ion and Crack G r o w t h in Concrete by Means o f Fracture Mechanics and Fini te Elements", Cement and Concrete Research, V o l . 6, pp. 773-782, 1976. H in ton , M J , Kaddour , A S , and Soden, PD. " A comparison o f the predict ive capabil i t ies o f current fai lure theories for composite laminates, j udged against experimental evidence", Composites Science and Technology, V o l . 62, pp. 1725-1797, 2002a. 166 References H i n t o n , M J , Kaddour , A S , and Soden, P D . "Evaluat ion o f fai lure predict ion in composite laminates: background to 'part B ' o f the exercise", Composites Science and Technology, V o l . 62, pp. 1481-1488, 2002b. H in ton , M J and Soden, PD. "Predict ing fai lure in composite laminates: the background to the exercise", Composites Science and Technology, V o l . 58, pp. 1001-1010, 1998. H i tchen, SA, Og in , SL , Smi th , PA , and Soutis, C. "The effect o f f ibre length on fracture toughness and notched strength o f short carbon f ibre/ epoxy composites", Composites, V o l . 25, pp. 407-413, 1994. Hyakutake, H and Hagio , T. "Fracture o f FRP plates containing notches or a c ircular hole under tension", International Journal o f Pressure Vessels and P ip ing, V o l . 44, pp. 277-290, 1990. Jackson, W C and M a r t i n , R H . " A n Inter laminar tensile strength specimen", pp. 333-354, 1993. K a m i y a , S and Sekine, H . " A discussion on the damage mechanism and the apparent fracture strength o f notched f iber- re inforced composite laminates", Journal o f Composi te Mater ia ls, V o l . 3 1 , pp. 580-595, 1997. K a n j i , K M J . " A n invest igat ion into the m ixed mode delaminat ion behaviour o f br i t t le composite laminates", M.A.Sc. Thesis. Department o f Metals and Mater ia ls Engineer ing, The Un ivers i ty o f Br i t i sh Co lumbia , 2003. Kann inen, M F and Popelar, C H . Advanced Fracture Mechanics. N e w Y o r k : O x f o r d Un ivers i ty Press, 1985. K i m , JK, K i m , D S , and Takeda, N . "Notched Strength and Fracture Cr i ter ion in Fabric Composi te Plates", Journal o f Composi te Mater ia ls, V o l . 29, pp. 982-998, 1995. Kongshavn, I and Poursart ip, A . "Exper imenta l Invest igat ion o f a Strain Softening Approach to Predict ing Fai lure in Notched Fibre Reinforced Composi te Laminates" , Composites Science and Technology, V o l . 59, pp. 29-40, 1999. Kor tschot , M T and Beaumont, P W R . "Damage Mechanics o f Composi te Mater ia ls. I V . The Effect o f L a y - U p on Damage Growth and Notched Strength", Composites Science and Technology (Previously: Fibre Science & Technology) , V o l . 40, pp. 167-179, 1991. Kor tschot , M T , Beaumont, P W R , and Ashby, M F . "Damage Mechanics o f Composi te Mater ia ls. I I I . Predict ion o f Damage G r o w t h and Notched Strength", Composites Science and Technology (Previously: Fibre Science & Technology) , V o l . 40, pp. 147-165, 1991. 167 References Kortschot , M T and Beaumont, PR. "Damage Mechanics o f Composi te Mater ia ls. I. Measurements o f Damage and Strength", Composites Science and Technology, V o l . 39, pp. 289-301 , 1990a. Kor tschot , M T and Beaumont, PR. "Damage Mechanics o f Composi te Mater ia ls. I I . A Damaged-Based Notched Strength M o d e l " , Composites Science and Technology, V o l . 39, pp. 303-326, 1990b. L i fsh i tz , J M and Leber, H . "Response o f f iber-re inforced polymers to h igh strain-rate loading in inter laminar tension and combined tension/shear", Composites Science and Technology, V o l . 58, pp. 987-996, 1998. M i t c h e l l , J. "Notched Behaviour o f Sti tched and Unst i tched R F I CFRP Laminates" , M.A.Sc. Thesis. Department o f Metals and Mater ia ls Engineer ing, The Un ivers i ty o f Br i t i sh Co lumbia , 2002. Pipes, R B , Wetherho ld , R C , and Gi l lespie, JW. "Notched strength o f composite mater ials", Journal o f Composi te Mater ia ls, V o l . 13, pp. 148-160, 1979. Shin, CS and W a n g , C M . " A n improved cohesive zone model for residual notched strength predic t ion o f composite laminates w i t h di f ferent or thotropic lay-ups", Journal o f Composi te Mater ia ls, V o l . 38, pp. 713-736, 2004. Sih, G C , Paris, PC, and I r w i n , GR. " O n cracks in rect i l inear ly anisotropic bodies", Internat ional Journal o f Fracture Mechanics, V o l . 1, pp. 189-202, 1965. Srawley, JE and Gross, B. "Stress intensity factors for bend and compact specimens", Engineer ing Fracture Mechanics, V o l . 4, pp. 587-589, 1972. Tada, H , Paris, P, and I r w i n , GR. The Stress Analys is o f Cracks Handbook. St. Lou is , M issour i : Paris Productions Incorporated and D e l Research Company, 1985. Tan, SC. "Ef fect ive Stress Fracture Mode ls for Unnotched and No tched Mul t id i rec t iona l Laminates" , Journal o f Composi te Mater ia ls, V o l . 22, pp. 322-340, 1988. Waddoups, M E , Eisenmann, JR, and K a m i n s k i , B E . "Macroscopic fracture mechanics o f advanced composite mater ials", Journal o f Composi te Mater ia ls , V o l . 5, pp. 446-454, 1971. Wh i tney , J M and Nuismer, RJ. "Stress Fracture Cr i ter ia for Laminated Composites Conta in ing Stress Concentrat ions", Journal o f Composi te Mater ia ls, V o l . 8, pp. 253-265, 1974. W i l l i a m s , K V , Vaz i r i , R, and Poursart ip, A . " A Physical ly Based Cont inuum Damage Mechanics M o d e l for T h i n Laminated Composi te Structures", International Journal o f Solids and Structures, V o l . 40, pp. 2267-2300, 2003. 168 References W i s n o m , M R , Jones, M I , and H i l l , GFJ. "Inter/laminar tensile strength o f carbon f ibre-epoxy - specimen size, layup and manufactur ing effects", Advanced Composi te Letters, V o l . 10, pp. 171-177, 2 0 0 1 . X u , X - P and Needleman, A . "Numer ica l simulat ions o f dynamic crack g rowth along an interface", Internat ional Journal o f Fracture, V o l . 74, pp. 289-324, 1996. X u , X P and Needleman, A . "Con t inuum mode l l ing o f interfacial decohesion", So l id State Phenomena, V o l . 35-36, pp. 287-302, 1994. Z o u , Z , Reid, SR, and L i , S. " A cont inuum damage model for delaminations in laminated composites", Journal o f Mechanics and Physics o f Sol ids, V o l . 5 1 , pp. 333-356, 2003. 169 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0078773/manifest

Comment

Related Items