A Non-Local Approach to Simulation of Damage in Composite Structures by Alireza Forghani B.Sc., Sharif University of Technology, 2003 M.Sc., Sharif University of Technology, 2005 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Civil Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) October 2011 © Alireza Forghani, 2011 ii ABSTRACT Fibre reinforced polymers (FRP) are one of the fastest growing materials in advanced structural applications. Due to the complexity and existence of multiple and interactive modes of failure, there is a lack of comprehensive theory that describes the damage behaviour of these materials. The UBC Composite Damage Model (CODAM) is a sub-laminate based model that is designed to simulate the behaviour of laminated composites at the macro (structural) scale. Physical basis and numerical robustness are the main objectives of CODAM development. The original formulation of CODAM that was developed in the mid 1990’s suffers from material objectivity problem that results in spurious dependency of numerical predictions on the choice of the coordinate system. In this thesis, the objectivity issue of the original CODAM is identified and addressed through a new formulation called CODAM2. The new formulation is capable of predicting damage in various laminate lay-ups from quasi-isotropic to unidirectional. It is also capable of simulating the damage-induced orthotropy in the laminate. An orthotropic non-local averaging scheme is developed for CODAM2 to address the localization issue. Compared to isotropic averaging, the orthotropic scheme significantly improves the predicted crack direction and damage pattern in composite laminates. Unlike the conventional isotropic non-local averaging that performs the averaging over a circular (spherical in 3D) zone, in the orthotropic scheme the averaging is performed over an elliptical zone which requires the introduction of two length scales. The performance of CODAM2 equipped with orthotropic averaging is demonstrated through numerical examples. It is shown that the developed model is capable of accurately predicting the damage behaviour in various specimen geometries from sharp-notched to blunt-notched and open-hole specimens. The predictions of the model in terms of load-displacement, crack- tip position, damage height and crack directions agree well with experimental observations and measurements. iii CODAM2 provides a promising numerical tool to simulate the effect of damage on the behaviour of structures made of laminated composites. This model is computationally efficient and yet relatively simple to understand, calibrate and use in practical applications. iv PREFACE This thesis titled “A non-local approach to simulation of damage in laminated composites” presents the research performed by Alireza Forghani. The research was supervised by Dr. Reza Vaziri and co-supervised by Drs. Anoush Poursartip, Fernand Ellyin and Dominik Schotzau at the University of British Columbia. Figures 3-5, 3-13 and 3-14 were published in (Forghani and Vaziri 2009) “Computational Modelling of Damage Development in Composite Laminates Subjected to Transverse Dynamic Loading”, J. Appl. Mech. Trans. ASME, 76 (5), 051304. These figures reflect the work of the thesis author on simulation of impact and blast on composite laminates. A version of contents of section 5-5 is accepted for publication (Forghani et al. 2011) “A Structural Modelling Framework for Prediction of Damage Development and Failure of Composite Laminates”, Composites Sci. Technol. This section includes the original work of the thesis author on development of CODAM2. A version of the results presented as the numerical example in section 9.5 is accepted for publication (Forghani et al. 2011) “A Structural Modelling Framework for Prediction of Damage Development and Failure of Composite Laminates”, Composites Sci. Technol. These results were produced using the code that was developed by the author of this thesis. v TABLE OF CONTENTS ABSTRACT ............................................................................................................................. ii PREFACE ............................................................................................................................... iv TABLE OF CONTENTS ........................................................................................................v LIST OF TABLES ...................................................................................................................x LIST OF FIGURES .............................................................................................................. xii LIST OF SYMBOLS ......................................................................................................... xviii ACKNOWLEDGEMENTS ..................................................................................................xx 1. INTRODUCTION..........................................................................................................1 2. NUMERICAL SIMULATION OF DAMAGE AND FRACTURE...........................3 2.1 Introduction ......................................................................................................................3 2.2 Material Damage, Fracture and Structural Failure ..........................................................3 2.3 Approaches to Simulation of Fracture and Damage ........................................................4 2.3.1 Continuum Approach ................................................................................................4 2.3.2 Discrete Approach .....................................................................................................5 2.3.3 Methods of Embedded Discontinuity and X-FEM ...................................................7 2.4 Comparison of Approaches..............................................................................................8 2.5 Summary ..........................................................................................................................9 3. LOCALIZATION PROBLEM AND NON-LOCAL REGULARIZATION .........14 3.1 Introduction ....................................................................................................................14 3.2 Strain-Softening Behaviour ...........................................................................................14 3.3 Localization Problem .....................................................................................................15 3.4 Remedies to the Localization Problem ..........................................................................16 3.4.1 Crack-band Approach .............................................................................................17 3.4.2 Mesh Orientation Dependency of the Local Damage Models ................................18 3.4.3 Rate Dependent Damage Models ............................................................................19 3.4.4 Non-Local Methods.................................................................................................20 3.4.5 Non-local Averaging Formulation (Isotropic Averaging) ......................................20 3.4.6 Gradient Non-Local Formulations ..........................................................................22 3.4.7 Numerical Examples with Isotropic Averaging ......................................................23 3.5 Summary ........................................................................................................................24 vi 4. SIMULATION OF DAMAGE IN LAMINATED COMPOSITES.........................39 4.1 Introduction ....................................................................................................................39 4.2 Damage Modes ..............................................................................................................39 4.3 Simulation of Damage at Various Scales.......................................................................40 4.3.1 Micro-scale Modelling ............................................................................................40 4.3.2 Meso-scale Modelling .............................................................................................41 4.3.3 Macro-scale Modelling ...........................................................................................42 4.3.4 Multi-scale Modelling .............................................................................................42 4.4 Simulation of Delamination Damage.............................................................................43 4.5 Simulation of Intra-Laminar Damage ............................................................................44 4.6 Continuum Damage Mechanics for Simulation of Intra-laminar Damage ....................45 4.6.1 Determining the Stiffness of a Damaged Material ..................................................45 4.6.2 Damage Evolution Laws .........................................................................................46 4.7 Sub-laminate Approach .................................................................................................47 4.8 Summary ........................................................................................................................48 5. UBC COMPOSITE DAMAGE MODEL (CODAM) ...............................................50 5.1 Introduction ....................................................................................................................50 5.2 Original CODAM ..........................................................................................................50 5.2.1 Background .............................................................................................................50 5.2.2 Formulation .............................................................................................................51 5.3 Objectivity Issue ............................................................................................................52 5.4 Solutions to the Objectivity Issue ..................................................................................53 5.5 CODAM2 .......................................................................................................................54 5.5.1 Background .............................................................................................................54 5.5.2 Formulation .............................................................................................................54 5.5.3 Damage Evolution ...................................................................................................55 5.5.4 Thermodynamic Aspects .........................................................................................59 5.5.5 Stress Locking and Transition to Scalar Damage ...................................................62 5.5.6 Edge Delamination ..................................................................................................63 5.6 Summary ........................................................................................................................65 6. ORTHOTROPIC NON-LOCAL AVERAGING ......................................................73 6.1 Introduction ....................................................................................................................73 6.2 Motivation ......................................................................................................................73 vii 6.3 Formulation ....................................................................................................................74 6.4 Implementation ..............................................................................................................75 6.5 Treatment of Boundaries................................................................................................75 6.6 Prediction of Damage/Crack Patterns ............................................................................76 6.6.1 Unidirectional Laminates ........................................................................................76 6.6.2 Multi-directional Laminates ....................................................................................77 6.6.3 Relation between Fracture Energy and Non-local Averaging Radii .......................78 6.7 Summary ........................................................................................................................80 7. VERIFICATION..........................................................................................................96 7.1 Introduction ....................................................................................................................96 7.2 Simulating Single Elements ...........................................................................................98 7.3 Verification of the Averaging Procedure .......................................................................98 7.4 Comparison between OOFEM and ABAQUS Predictions ...........................................98 7.5 Frame Independency ......................................................................................................99 7.6 Summary ........................................................................................................................99 8. MODEL CALIBRATION .........................................................................................110 8.1 Introduction ..................................................................................................................110 8.2 OCT Test ......................................................................................................................110 8.3 Evaluation of the Non-local Radii ...............................................................................111 8.3.1 Longitudinal Averaging Radius, r1 .......................................................................111 8.3.2 Transverse Averaging Radius, r2 ..........................................................................111 8.4 Calibration of Initiation and Saturation Strains ...........................................................112 8.4.1 Initiation Strains ....................................................................................................112 8.4.2 Saturation Strains ..................................................................................................112 8.4.3 Anti-Stress Locking Mechanism Parameters ........................................................115 8.4.4 Delamination Parameters ......................................................................................115 8.5 Model Calibration: an Example ...................................................................................116 8.5.1 Non-Local Radii ....................................................................................................116 8.5.2 Initiation Strains ....................................................................................................116 8.5.3 Saturation Strains ..................................................................................................117 8.6 Summary ......................................................................................................................118 9. MODEL VALIDATION AND NUMERICAL EXAMPLES ................................129 9.1 Introduction ..................................................................................................................129 viii 9.2 Validation .....................................................................................................................130 9.3 Scope, Applicability and Limitations ..........................................................................130 9.4 Simulation of the Modified OCT Tests with VarIous Notch Radii .............................132 9.4.1 Test Specifications ................................................................................................132 9.4.2 Model Parameters ..................................................................................................132 9.4.3 Predictions .............................................................................................................133 9.4.4 Comparison between the Non-local and Local FE predictions .............................136 9.5 Open Hole Tension ......................................................................................................138 9.5.1 Material Parameters...............................................................................................138 9.5.2 Numerical Predictions ...........................................................................................140 10. SUMMARY, CONCLUSIONS AND FUTURE WORK ........................................174 10.1 Summary ...................................................................................................................174 10.2 Concluding Remarks ................................................................................................176 10.3 Future Work ..............................................................................................................176 REFERENCES .....................................................................................................................178 APPENDICES ......................................................................................................................185 A. A GUIDELINE FOR USERS ...................................................................................185 A.1 General Notes ..................................................................................................................185 A.2 Compiling and Running the Code ...................................................................................185 A.3 Code Calibration .............................................................................................................187 B. DEMONSTRATION OF THE LOCALIZATION PHENOMENON IN A 1D PROBLEM ...........................................................................................................................188 B.1 Quasi-Static Problem.......................................................................................................188 B.2 Dynamic Problem ............................................................................................................189 B.3 Numerical Example .........................................................................................................190 C. FUNDAMENTALS OF CONTINUUM DAMAGE MECHANICS ...........................200 C.1 Reducing the Stiffness .....................................................................................................200 C.1.1 Strain Equivalence Hypothesis .................................................................................201 C.1.2 Energy Equivalence Hypothesis ...............................................................................203 C.1.3 Other Approaches .....................................................................................................204 C.2 Damage Evolution Laws .................................................................................................205 C.2.1 Maximum Dissipation Hypothesis ............................................................................205 C.2.2 Dissipation Potential .................................................................................................206 ix C.2.3 Equivalent Strain Functions ......................................................................................206 D. CODAM2 IMPLEMENTATION ..................................................................................208 D.1 Algorithm ........................................................................................................................208 D.2 Pseudo Code ....................................................................................................................209 D.2.1 Equivalent Strains .....................................................................................................209 D.2.2 Non-local Averaging ................................................................................................210 D.2.3 Calculating the Damage Parameters .........................................................................210 D.2.4 Calculating the Reduction Coefficients ....................................................................211 D.2.5 Calculating the Damaged Stiffness Matrix ...............................................................211 D.2.6 Calculating the Stress Vector....................................................................................211 D.3 Implementation ...............................................................................................................212 D.3.1 OrthoDamage2 ..........................................................................................................212 D.3.2 CODAM2 .................................................................................................................212 D.3.3 StructuralOrthoNonLocalMaterialExtensionInterface .............................................212 D.3.4 CODAM2NL ............................................................................................................213 x LIST OF TABLES Table 2-1. Capabilities and limitations of approaches for simulating fracture and damage. ......................................................................................................................... 10 Table 6-1. Input parameters used in the simulations of the notched specimen presented in this chapter. ..................................................................................................... 82 Table 6-2. The predicted fracture energy from the simulations of a notched specimen made from a unidirectional 90º laminate with different transverse non-local radii. The effective damage heights and their normalization with respect to r2 are presented. .................................................................................................. 82 Table 6-3. The predicted fracture energy from the simulations of a notched specimen made of a quasi-isotropic laminate with different non-local radii. The effective damage heights and their normalization with respect to r1 are presented. ...... 83 Table 7-1. Elastic properties of the material used in verification of a single element. .. 100 Table 7-2. Damage properties of the material used in verification of a single element. 100 Table 8-1. Summary of calibration procedure for non-local CODAM2 parameters. .... 119 Table 8-2. Elastic Properties of IM7-8552 unidirectional lamina reported by the manufacturer (Hexcel Company). ................................................................. 119 Table 8-3. Material parameters for IM7-8552 quasi-isotropic laminate extracted from OCT tests performed by N. Zobeiry (Zobeiry 2010). ................................... 120 Table 8-4. The values of fracture energy density and expected fracture energy calculated for various values of saturation strain ........................................................... 120 Table 8-5. Damage model parameters used in calibration of the model. ...................... 121 Table 9-1. Elastic properties of the laminate used in the modified OCT simulations (McClennan 2004). ....................................................................................... 142 Table 9-2. Damage model parameters obtained based on sharp-notched OCT test data. ....................................................................................................................... 143 Table 9-3. Labelling and dimensions of the simulated specimens considered for the open- hole tension specimens (WWFE-III’s Test Case 12). ................................... 144 Table 9-4. Estimated values of processing-induced residual stresses in various layers due to the layer-wise CTE mismatch in quasi-isotropic laminates considered in Test Case 12. ................................................................................................. 145 Table 9-5. Calculated range of critical strain for onset of delamination in quasi-isotropic laminates considered in Test Case 12. .......................................................... 146 Table 9-6. Values of the input parameters used in numerical predictions of Test Case 12. ....................................................................................................................... 147 Table 9-7. Details of material parameters used in the numerical simulation of Test Case 12................................................................................................................... 148 xi Table 9-8. Summary of predictions of Test Case 12. ..................................................... 149 xii LIST OF FIGURES Figure 2-1. A schematic representation of the finite element field using the methods of embedded discontinuity: (a) two neighbouring elements with a single weak discontinuity in each form a localized band; (b) the localized band is simulated as two weak discontinuities inside one element; and (c) a strong discontinuity inside an element. This figure is reproduced from Figure 1 in (Jirasek 2000). ................................................................................................. 12 Figure 2-2. A schematic chart showing some of the current approaches to model fracture and damage in solids. ...................................................................................... 13 Figure 3-1 (a) A typical Overheight Compact Tension (OCT) specimen; and b) A typical load-POD curve extracted from test on an OCT specimen............................. 26 Figure 3-2. A schematic, showing a finite element that uses strain-softening as an input and predicts load-POD response. .................................................................... 27 Figure 3-3. In crack-band method, the softening part of the stress-strain curve is scaled according to the height of the mesh. ............................................................... 28 Figure 3-4. FE prediction of crack growth in the simulation of a notched specimen made from an isotropic material using MAT_81 local in LS-DYNA. ..................... 28 Figure 3-5. Predicted damage patterns using the crack-band approach on (a) a plate under lateral impact; and (b) a plate subjected to blast load. .................................. 29 Figure 3-7. The averaging zone ΩX around the point X and the bell-shaped weight function w employed for averaging the equivalent strains. ............................ 31 Figure 3-8, (a) Dimensions of the over-height compact tension specimen (OCT) (Kongshavn and Poursartip 1999); and (b) the applied force and the pin- opening displacement...................................................................................... 32 Figure 3-9. Force versus pin-opening displacement (POD) observed from simulations of an over height compact tension (OCT) specimen with different mesh sizes when the equivalent strain ( eq ) is chosen as the non-local averaging parameter......................................................................................................... 33 Figure 3-10. Force-POD predictions observed from OCT simulations with different mesh sizes when Damage Parameter (ω) is chosen as the non-local averaging variable. The erosion of elements with large strains was prevented. .............. 34 Figure 3-11. Force-POD predictions observed from OCT simulations with different mesh sizes when Damage Parameter is chosen as the non-local averaging variable. The erosion strain was set to 150% in these simulations. ............................... 35 Figure 3-12. Predicted crack pattern in an OCT specimen with (a) local crack band approach and (b) non-local averaging model. ................................................. 36 Figure 3-13. Predicted damage pattern in a plate subjected to lateral impact; (a) and (c) using local crack-band model; and (b) and (d) using non-local averaging. .... 37 xiii Figure 3-14. Predicted damage pattern in a plate subjected to lateral blast load using (a) crack-band approach; and (b) non-local averaging. ........................................ 38 Figure 4-1. A schematic showing various scales in studying and simulating the mechanical behaviour of laminated composites. ............................................ 49 Figure 5-1. (a) Schematic showing a RVE of a cross-ply laminate at three stages: (1) before the damage starts (2) formation of microcracks mainly in the matrix (3) microcracks merge and form a through crack in the laminate. (b) schematic nominal stress vs. nominal strain behaviour of the RVE under applied deformation. .................................................................................................... 66 Figure 5-2. Predicted response of an OCT specimen with a quasi-isotropic lay-up. In one case, the local axes are aligned with the crack and in the other one they make a 45º angle. ...................................................................................................... 67 Figure 5-3. An element under unia ial load assuming that the material’s local CS is (a) aligned and ( ) makes angle. .................................................................... 68 Figure 5-4. A schematic showing (a) the effective strain-softening behaviour of layers in the laminate being combined to result in (b) the effective strain-softening behaviour of laminate. .................................................................................... 69 Figure 5-5. (a) A typical damage parameter defined as a hyperbolic function of the corresponding non-local equivalent strain, and (b) the reduction parameter as a linearly decreasing function of the damage parameter. ............................... 70 Figure 5-6. Schematic showing a stress locking situation in a widely opened crack. The two elements shown are on the crack path and subjected to large opening. The assumed strain states are shown on the corresponding Mohr circles. The undamaged fibres in these neighbouring elements form a load carrying path. ......................................................................................................................... 71 Figure 5-7. Schematic showing (a) reduction of laminate stiffness due to delamination, and (b) the gradual drop in the laminate stiffness as implemented numerically. 72 Figure 6-1. Schematic showing two different situations for a material point A where the line connecting the existing damage zones (D1 and D2) to the point A is (a) parallel to the fibres and (b) perpendicular to the fibres. ................................ 84 Figure 6-2. Schematic showing the influence area of (a) a classical isotropic averaging scheme and (b) the proposed orthotropic non-local averaging scheme. The parallel lines show the direction of the fibres. The shaded zones (D1 and D2) represent existing damage zones in the neighbourhood of point A. ............... 85 Figure 6-3. Defining d1 and d2 parameters as the components of the AB vector in the fi re’ direction and transverse direction respectively. .................................... 86 Figure 6-4. The orthotropic weight function plotted in the local coordinate system of the ply (x1 and x2). In this example, r1=2 and r2=1. .............................................. 87 xiv Figure 6-5. Predicted crack growth patterns obtained from the simulation of a notched specimen for a 0º unidirectional laminate using various values of the non-local radii. The dashed lines indicate the fibre orientations. ................................... 88 Figure 6-6. Predicted damage growth patterns in a 60º unidirectional notched specimen (a) using a local damage model (crack-band), (b) an isotropic averaging with (r1=r2=2mm) and (c) an orthotropic averaging with (r1=2mm and r2=0.5mm) . ......................................................................................................................... 88 Figure 6-7. Schematic showing the averaging zone around a typical point in (a) an isotropic, and (b) an orthotropic averaging scheme. ....................................... 89 Figure 6-8. Predicted height of damage in (a) a unidirectional laminate and (b) a multi- directional laminate. ........................................................................................ 90 Figure 6-9. Predicted damage patterns in a notched specimen for a quasi-isotropic laminate with (a) an isotropic weight function with a radius of 2mm, (b) an orthotropic weight function with r1=2mm and r2=1.0mm, and (c) an orthotropic weight function with r1=2mm and r2=0.5mm. ............................. 91 Figure 6-10. Splitting crack forms parallel to the fibres and along the initial notch in the unidirectional 90º laminate. ............................................................................ 92 Figure 6-11. A schematic showing the distinction between the visual and effective damage heights in a section (C-C) through the predicted damage profile. .................. 93 Figure 6-12. Predicted Load vs POD curves for an OCT simulation of a quasi-isotropic laminate. .......................................................................................................... 94 Figure 6-13. An schematic demonstrating the scaling of the averaging weight function near the boundaries of the structure. ....................................................................... 95 Figure 7-2. A single element used for verification of CODAM2, (a) under uniaxial strain along the direction of the fibre, and (b) under proportionally increasing combined transverse and shear strains. ......................................................... 102 Figure 7-3. Predicted stress-strain curves compared to the input data for cases (a) and (b) in Figure 7-2. ................................................................................................. 103 Figure 7-4. A square specimen subjected to uniform uniaxial displacement used in verification of non-local averaging. The dashed line shows the position of the cross-section at which the local and averaged equivalent strains are evaluated and shown in Figure 7-5. .............................................................................. 104 Figure 7-5. Distributions of the predicted local equivalent strain compared to the averaged equivalent strain. ........................................................................................... 105 Figure 7-6. Dimensions of the open-hole specimen simulated for verification of the predicted stress field in the elastic regime. ................................................... 106 Figure 7-7. OOFEM’s predicted stress distri ution over a horizontal cross-section passing through the centre of the open-hole specimen compared with the corresponding predictions using ABAQUS. ................................................. 107 xv Figure 7-8. Schematic showing the choice of the material coordinate system, (a) x1 axis is chosen to be parallel to the notch, and (b) x1 axis makes a 45º angle with the notch. ............................................................................................................. 108 Figure 7-9. Predicted load-POD curves for the OCT specimens with two different choices for the material coordinate systems. ............................................................. 109 Figure 8-1. Schematic showing the OCT specimen geometry and its typical dimensions. ....................................................................................................................... 122 Figure 8-2. Increasing the transverse non-local radius, r2, while keeping the other parameters constant results in decreasing saturation strain for the matrix damage. ......................................................................................................... 123 Figure 8-3. The iterative procedure that is used to calibrate the saturation strains given the height of damage and fracture energy. .......................................................... 124 Figure 8-4. The predicted Load-POD curves from simulation of OCT specimen for various values of the saturation strains. CTA1 and CTA2 represent the experimentally obtained curves by (Zobeiry 2010). ..................................... 125 Figure 8-5. Predicted Load-POD curve for the calibrated model compared to the measured data................................................................................................ 126 Figure 8-6. The predicted position of damage initiation and saturation fronts as a function of POD compared to the experimental measurements. ................................. 127 Figure 8-7. The predicted position of damage initiation and saturation fronts. ............... 128 Figure 9-1. An schematic of a damage zone surrounded by undamaged material. .......... 150 Figure 9-2 (a) Dimensions of the OCT specimen and (b) modified OCT specimen. ..... 151 Figure 9-3. Modified OCT specimens with varying notch radii (a) 0.5mm, (b) 2mm, (c) 5mm, (d) 12.75mm, (e) 19mm and (f) 27mm.(McClennan 2004) ............... 152 Figure 9-4. Predicted load-CMOD curves compared to the experimental data for the sharp-notched OCT specimen (McClennan 2004). ...................................... 153 Figure 9-5. Predicted load-CMOD curves compared to the experimental data for the modified OCT with ρ=12.7 mm (McClennan 2004). ................................... 154 Figure 9-6. Predicted load-CMOD curves compared to the experimental data for the modified OCT with ρ=16mm (McClennan 2004). ....................................... 155 Figure 9-7. Comparison between the predicted load-displacements for the sharp-notched specimen and the blunt-notched specimen with ρ=16mm (McClennan 200 ). 156 Figure 9-8. Predicted peak loads versus notch root radius compared to measured test data (McClennan 2004). ....................................................................................... 157 Figure 9-9. Predicted length of damage zone and length of zone with saturated damage shown over the load-CMOD for the sharp-notched specimen...................... 158 xvi Figure 9-10. Predicted length of damage zone and length of zone with saturated damage shown over the load-CMOD for the blunt-notched specimen with ρ=12.7 mm ....................................................................................................................... 159 Figure 9-11. Predicted length of damage zone and length of zone with saturated damage shown over the load-CMOD for the blunt-notched specimen with ρ=22. mm. 160 Figure 9-12. Predicted matrix damage pattern in the 90º layer for specimens with notch radii of (a) 0.75mm, (b) 3.2mm, (c) 12.7mm and (d) 22.5mm. .................... 161 Figure 9-13. Predicted damage height evolution in front of the notch for specimens with the very sharp and the very blunt notch root radii compared to experimental measurements (McClennan 2004). ............................................................... 162 Figure 9-14. The position of the cross-sections C1 and C2 at which damage profile is plotted, (b) Predicted damage parameter profile corresponding to matrix cracking for specimens with different notch root radii at cross-section C1 immediately in front of the notch and (c) at cross section C2, 20mm away from the original notch tip. ........................................................................... 163 Figure 9-16. (a) and (b) Predicted matrix damage patterns in the 0º layer using local and non-local models, respectively. (c) and (d) predicted longitudinal stress, σy distribution in front of the crack using local and non-local damage models. 165 Figure 9-17. Schematic of the open-hole specimen geometries: (a) square and (b) tall specimen. ...................................................................................................... 166 Figure 9-18. Assumed R-Curve behaviour for the progression of delamination in the open- hole tensile specimen. The solid lines show the estimated value of Gc in each zone and the dashed line shows the increasing trend of the R-Curve. The right hand side axis shows the delamination initiation strain corresponding to the fracture energy on the left axis The correspondence is through the non-linear relation shown in Equation (5.25). ................................................................ 167 Figure 9-19. Predicted nominal stress vs nominal strain curves for the open hole tensile specimens with four different hole sizes. The predicted curves for both square and tall specimen geometries are shown in each case. ................................. 168 Figure 9-20. Predicted nominal stress vs nominal strain curves for the open hole tensile specimens with the hole radius of 3.175mm. The predictions are compared to the experiments performed by Hallett et al. (Hallett et al. 2009). ................ 169 Figure 9-21. Predicted progression of damage in the smallest specimen (square specimen with a central hole size of 3.175mm) ............................................................ 170 Figure 9-22. Example of a predicted stress-strain curve (for the D3.175S specimen) showing the failure strength (defined as the instant when the laminate secant stiffness drops by 20%) and the ultimate strength. ....................................... 171 Figure 9-23. Predicted Failure and Peak Stresses versus Hole Size. ................................. 172 Figure 9-24. Predicted delamination area in different open-hole tensile specimen. ......... 173 xvii Figure D-1. Flowchart of CODAM2 calculations ............................................................ 214 Figure D-2. CODAM2 and CODAM2NL class hierarchy ............................................... 215 xviii LIST OF SYMBOLS A Laminate’s in-plane stiffness matrix Δa Crack length advancement C Material compliance matrix d1, d2 Components of the vector connecting the source point to a generic point in the orthotropic non-local averaging D Hole diameter in open hole specimens D Material stiffness matrix E 0 Initial modulus of the laminate E * Modulus of the delaminated laminate E1 , E2 In-plane longitudinal and transverse modulus of a ply gf Fracture energy density (area under the strain-softening curve) Gc Delamination energy Gf Laminate’s in-plane fracture energy m fG Fracture energy of the matrix FE fG Predicted fracture energy from simulation of OCT specimen used in model calibration G12 In-plane shear modulus of a ply hc Characteristic damage height hd Predicted damage height (or height of one element) in crack band modelling H Height of open hole specimens M Damage effect tensor n Total number of layers in a sub-laminate nD Number of delaminated interfaces through-thickness Q The plane stress in-plane stiffness matrix of a lamina in the principal orthotropic plane r1, r2 Longitudinal and transverse non-local averaging radii R Stiffness reduction coefficient R Vector containing the stiffness reduction coefficients T Transformation matrix for strain t Laminate thickness, also used as current time depending on the context tk k th Lamina thickness xix ΔU Change in the strain energy w Non-local averaging weight function dW Rate of energy dissipation in damage process ΔW Change in the external work x, X Position vectors Y Work conjugate of the state parameters Y Vector containing the Y parameters εc Critical delamination strain εeq Equivalent strain eq Averaged (non-local) equivalent strain γ Fracture energy density λ Eigenvalues of the ply’s stiffness matri ψ Elastic (recoverable) strain energy density ρeff Effective radius used in orthotropic averaging υ12 , υ21 Major and minor Poisson’s ratio of a ply ω Damage parameter Rate of change of damage parameter ω Vector containing damage parameters ΩX Non-local averaging zone (volume) around point X. Subscript m Parameter used to indicate number of plies in a ply-scaled laminate k Generic layer of a sub-laminate m f L D Denotes the mode of damage: fibre ( f ), matrix ( m ), locking ( L ) and delamination ( D ) Superscript i Initiation s Saturation d Damaged 0 Undamaged (or initial) xx ACKNOWLEDGEMENTS This work would have not been possible without the unconditional love and support of my family: my love, Tahoura, my parents, Mohammad Taghi and Soheila to whom I dedicate this thesis and my brothers, Mohammad Javad and Erfan. I would like to express my sincere thanks to my supervisor Professor Reza Vaziri, who gave me the chance for being involved in interesting and challenging research projects and spent hundreds of hours of his time discussing various aspects of my work. His warm attitude has always inspired and encouraged me through these years. I would also like to thank my research advisors Professors Anoush Poursartip, Fernand Ellyin and Dominik Schoetzau for the fruitful discussions and helpful hints. I would like to express my appreciation to Professor Borek Patzak from the Czech Technical University in Prague for his helpful advises regarding the various aspects of numerical implementation of my code. In the past few years I have been a member of the UBC Composites Group. I would like to thank former and current members of this group, Navid Zobeiry, Carla McGregor, Mehdi Haghshenas, Kamyar Gordnian, Sardar Malek Mohammadi, Bryan Louis and Armin Bebamzadeh for their helpful discussions and joyful company during these years. I would like to express my gratitude to my dear friends in Vancouver, Ali Baghani, Alireza Ahmadnia, Mahdi Beheshtian, Mojtaba Mahsuli, Vahid Shah Mansouri, Mehran Shaghaghi, Pirooz Darabi, Saghi Kowsari, Amir Nejat and Ali Rasekh who have been my second family and provided the great support. The financial support for my study was provided by a combination of the UBC University Graduate Fellowship (UGF), and research grants from the Natural Sciences and Engineering Research Council of Canada (NSERC) and the Defence Research and Development of Canada (DRDC). 1 1. INTRODUCTION Fibre-reinforced polymers (FRPs) are one of the fastest growing materials in advanced structural applications. These materials have been in use in military aircraft for a long time but in the past decade, there has been a movement among major aviation companies to make commercial aircraft from laminated FRP. Application of fibre reinforced polymer composites is not limited to the aerospace industry. These materials are being used these days in cars, storage tanks and pressure vessels, sport equipment such as lightweight bikes and skies. They also have applications in civil engineering mainly to retrofit and repair existing structures. Based on the data published on the website of the Toray Company (Toray 2010), a major producer of car on fi res that provides the car on fi re for production of Boeing’s 787 aircraft, as of March 2010, 44% of the carbon fibre made by Toray is used in aerospace applications, 37% is used in industrial and civil engineering applications and 19% in sports applications (Toray 2010). Fibre reinforced polymers are often used in the form of laminates. One of the objectives in designing a laminate is to choose the number of layers and angle of each ply in order to minimize the weight while achieving the required mechanical properties. The elastic properties of a laminate can be relatively easily estimated based on laminate beam and plate theories. By increasing the load level, laminated composites exhibit complex behaviour. This is due to the various modes of failure that emerge and the extreme orthotropy of the laminate constituents, i.e. plies. Unlike conventional engineering materials, there is no comprehensive theory that can describe and predict the behaviour of a laminated composite subjected to ultimate loading. The last round of the Worldwide Failure Exercise (Hinton et al. 2002; Hinton et al. 2004; Kaddour et al. 2004; Soden et al. 1998), showed that the existing failure theories are not capable of predicting the failure of even simple coupons made of laminated composites under general loading scenarios. Furthermore there is a lack of models that are practically applicable to real- size composite structures. 2 The UBC Composites Group has a long history of research on simulation of composite laminates and predicting their damage behaviour. Being physically based and numerically robust have always been the main objectives of the damage models developed in our group. This thesis includes the author’s contribution to the ongoing development of damage models in the UBC Composites Group. The non-locally regularized damage model proposed in this thesis aims to bridge the gap between the academic research and engineering practice by developing a model that is physically-based and computationally efficient and yet relatively simple to understand, calibrate and use. The thesis consists of ten chapters and four appendices. The next chapter (Chapter 2) is a review of the existing approaches to numerical simulation of fracture and damage. Chapter 3 discusses the localization problem and the non-local averaging as a solution to this issue. Chapter 4 presents a review of the current methods to simulate damage in laminated composites. Chapter 5 discusses the sub-laminate approach and the UBC Composite Damage Model (CODAM). The second generation of CODAM, i.e. CODAM2, is introduced in this chapter that addresses material objectivity issues that the original CODAM suffered from. Chapter 6 introduces the orthotropic averaging scheme that is specifically designed for highly orthotropic materials. The performance of the orthotropic non-local model to predict damage pattern and direction of crack growth is demonstrated in this chapter. Chapter 7 is on verification of the implementation of CODAM2 and its non-local averaging scheme. Chapter 8 is on characterization and calibration of the parameters in CODAM2 and the associated orthotropic averaging. Chapter 9 presents some numerical examples to demonstrate the performance of the proposed model to predict the damage behaviour of the specimens made of laminated CFRP. The scope and limitations of the model are discussed in this chapter and the numerical examples are used to validate the model. The conclusions and future works are presented in Chapter 10. Appendix A presents a general guideline for the users of CODAM2 and OOFEM. Appendix B demonstrates the localization problem using a in a 1D example. Appendix C is on various aspects of continuum damage model such as commonly used damage formulations and evolution laws. Appendix D describes the implementation of CODAM2 in the open-source finite element code, OOFEM. 3 2. NUMERICAL SIMULATION OF DAMAGE AND FRACTURE 2.1 INTRODUCTION All solid materials, when subjected to excessive loading, lose their load carrying capacity. This phenomenon is seen in various engineering materials such as metals, concrete, polymers, ceramics, composites and wood. Depending on the microstructure of the material, this may be due to dislocation movements (slip), formation of micro-cracks and cracks, de-bonding or micro-buckling. Simulation of the effect of material damage and fracture on the behaviour of structures has always been a challenge for researchers and engineers. The multi-scale nature of damage that initiates with micro-cracks and turns into macroscopic cracks makes it difficult to come up with a model that can accurately predict the initiation of damage as well as its growth. It should be noted that crack propagation is a structural problem and not a local material problem. In other words knowing the local state of stress/strain at a material point is not adequate to estimate when the crack would start growing. The energy stored in the structure that is released by growth of crack acts as a driving force to grow the crack further. Therefore, evolution of fracture and damage cannot be studied as a material behaviour without including the structure. In addition to physical objectivity and applicability issues, when it comes to simulation of damage and softening behaviour, finite element models face various numerical issues such as mesh size and orientation dependencies. This chapter presents an overview of the existing approaches towards simulation of damage and fracture in solids. 2.2 MATERIAL DAMAGE, FRACTURE AND STRUCTURAL FAILURE Failure of a structure is subjective and has to be defined based on the type of the structure and its application. In some structures such as pipelines and pressure vessels, formation of a small crack may lead to leakage and cause failure. In other structures, failure may be defined as large deformation of the structure or its inability to carry further load. Formation of damage and fracture in materials is one of the major causes of structural failure. However it should be 4 noted that not all cases of structural failure are caused by damage in the material (e.g. elastic buckling). 2.3 APPROACHES TO SIMULATION OF FRACTURE AND DAMAGE One can classify the current approaches to simulation of damage and fracture of the material into two major categories: (i) continuum and (ii) discontinuous approaches. Continuum approach is formulated based on the continuum damage mechanics and the discontinuous approach is defined based on fracture mechanics theories. Here the advantages and shortcomings of each of these approaches are discussed in detail. 2.3.1 Continuum Approach Damage and formation of crack in engineering materials always start with changes in the micro-structure of the material. Formation of slip planes and dislocations in metals, micro- cracks in concrete and matrix of polymer-based composites are examples of such changes. The goal of continuum approaches is to estimate the effect of material degradation on the effective (macroscopic) properties of the material such as its stiffness. Continuum models are based on smearing and homogenization of the effect of cracks and degraded material over a representative volume. By defining the damage parameters as state parameters, these models attempt to express the relation between the effective (nominal) stresses and strains. In the classical continuum damage theories, the damage parameters are often related to the density of cracks in the representative volume. The representative volume element (RVE) in the context of damage refers to the zone that experiences strain localization. It should be noted that the RVE in the context of damage mechanics is fundamentally different from the RVE that is assumed in micromechanics to estimate properties of the undamaged material. The size of the representative volume is related to the size (height) of damage observed in experimental tests and depends on the material’s microstructure and its constituents. In concrete it mainly depends on the size of the aggregate. In laminated composites, the size of damage depends on the lay-up of the laminate, the architecture of the fibre bed (woven or unidirectional) and also the loading mode (e.g. being tensile or compressive). In typical FRP laminated composites, the size of the damage zone is of the order of a few millimetres. 5 A macroscopic crack forms as an ultimate condition when the damage is saturated and the material completely loses its load carrying capability. The continuum damage mechanics approach offers a relatively good representation of the behaviour of the quasi-brittle materials that exhibit formation of micro cracks and material that are spatially distributed over a non-negligible fracture process zone. In a continuum-based damage formulation, damage initiation and evolution laws are often written as functions of state of strain (and strain rate in rate dependent materials). The secant stiffness of the material is then written in terms of the state of damage. The straightforward formulation and relatively simple implementation in commercial FE codes are among the advantages of continuum approaches. The models based on this approach receive the state of strain at each integration point and calculate the corresponding state of stress. Damage parameters are typically stored at the integration points as history parameters. The scalar damage models (crack-band as an example (Bazant and Oh 1983)), rotating crack models (Jirasek and Zimmermann 1998) and the micro-plane models (Bazant and Oh 1985) are among the continuum-based approaches that are developed to simulate damage in originally isotropic materials such as concrete. Some of these models such as the rotating crack and micro-plane models are capable of simulating the damage-induced orthotropy in an originally isotropic material. Although writing the formulation to describe the behaviour of damaged anisotropic materials faces some complexity, the continuum damage approach is very popular in modelling stiffness degradation in laminated composites. The continuum damage models in their classical local form suffer from the well-known localization problem that is discussed in more detail in Chapter 4. 2.3.2 Discrete Approach Instead of employing continuum damage mechanics theory, the discrete methods follow the fracture mechanics approach where the localized damage zone is replaced with a crack line (or surface in 3-D models). This approach offers a good representation of physical reality when damage forms in a narrow localized band. 6 The discrete approach to modelling fracture has been used in different forms. Some of the discrete models simulate the crack as separation between boundaries of elements. In these models, remeshing of the system is needed to simulate crack growth (see (Carpinteri 1989), (Bouchard et al. 2000) and (Réthoré et al. 2004) for instance). (Movahhedy et al. 2000) and (Gadala et al. 2002) have used the Arbitrary Lagrangian-Eulerian method (ALE) and remeshing to simulate severe changes in the structural geometry due to metal cutting. The other popular approach within the discrete context is introduction of interface elements or tie-break type contacts between all or some of the elements as the crack potential (e.g. (Camacho and Ortiz 1996) and (Yang and Chen 2004)). The criteria governing the crack evolution in a discrete model can be based on either LEFM or Cohesive Crack models. When LEFM condition is assumed, the fracture energy release rate is numerically calculated using either the J-Integral method or the Virtual Crack Closure Technique (VCCT) (e.g (Krueger 2004)). These methods have been implemented in some commercial FE software packages such as ABAQUS. The LEFM-based models require an initial crack and therefore they are unable to predict the initiation of crack in an undamaged material. The cohesive crack model is the approach that was proposed to simulate the effect of the non- linear zone ahead of the crack tip observed in metals and quasi-brittle materials. The cohesive crack method defines a softening traction-opening law that governs the behaviour of fracture process zone. A comprehensive overview of cohesive crack model can be found in (Elices et al. 2002) . In the finite element simulation of the cohesive cracks, the traction-opening law is applied to the interface elements or tie-break type contacts (e.g. (Camacho and Ortiz 1996) and (Forghani and Vaziri 2009)). The discrete approach is a popular method for simulation of delamination in laminated composites. This is due to the fact that delamination is often formed as localized cracks and therefore the cohesive crack is a suitable representation of delamination cracks (see (Alfano and Crisfield 2001) and (Forghani et al. 2007) for instance). The cohesive crack is also theoretically capable of predicting mixed-mode cracking conditions. 7 Commercial FE codes like LS-DYNA and ABAQUS, have the capability of simulating cohesive cracks through their interface elements and contact formulations where a softening law governs the traction-opening of the cohesive interface. From the computational modelling viewpoint, the major issue that discrete models face is that the crack is modeled between the boundaries of elements. Therefore these methods are effective when the potential direction of crack growth is known in advance (that is the case in delamination cracks). Otherwise, one has to either put cohesive interface between every two neighbouring elements (e.g. (Camacho and Ortiz 1996)) or an adaptive re-meshing is required to allow propagation of the crack in the structure (see (Carpinteri 1989) for instance). The discrete approach equipped with a cohesive law works best for notched specimens under quasi-static conditions when a localized crack is formed. In dynamic events where cracks are formed in a rather spatially distributed and diffused form, employing the cohesive model leads to a so-called inverse mesh size dependency. In this situation, refining the mesh results in further branching of the crack and increasing the total length of cracks which would result in an increase in the dissipated energy (see (Ruiz et al. 2000) for instance). 2.3.3 Methods of Embedded Discontinuity and X-FEM The methods of embedded discontinuity and the X-FEM are not independent from continuum or discrete approaches. These methods enable simulation of a localization band (continuum) or a discontinuity line (discrete) inside the elements. The idea of introducing discontinuity into the FE domain was first proposed by Ortiz and co- workers (Ortiz et al. 1987) in which a weak discontinuity (a jump in the strain field) was introduced to the element. In this approach, two neighbouring elements with a weak discontinuity form a localized band as schematically shown in Figure 2-1 (a). In a later study, Belytschko et al. (Belytschko et al. 1988) introduced two lines of weak discontinuity into the approximation field in an element. This allowed modelling of a localized band within an element as shown schematically in Figure 2-1 (b). And finally some element formulations were proposed with embedded strong discontinuity, i.e. jump in the displacement field (e.g. (Dvorkin et al. 1990)). In all these methods, the embedded discontinuities are introduced through internal degrees of freedom in each element. Because of the degrees of freedom being 8 internal, some issues regarding the compatibility of displacement field arise that needs to be carefully addressed (see (Jirasek 2000) for more details). The recently introduced mathematical theory of Partition of Unity Method (PUM) (Babuska and Melenk 1997) proposes a proper framework to enrich the finite element space. Using this approach, the Extended Finite Element Method (X-FEM) was developed by Belytschko and co-workers that simulates the crack by enriching the finite element space with jump functions (Belytschko and Black 1999). The degrees of freedom related to the jump in the displacement field are global resulting in a compatible displacement field and a continuous crack profile. The X-FEM approach has been employed successfully to simulate various types of cracking and discontinuity such as cohesive cracks (Forghani 2005; Moes and Belytschko 2002) and cracks with frictional surfaces (Dolbow et al. 2001). The PUM approach also provides the opportunity of simulating notch-tip singular stress/strain fields by employing special shape functions. The enriched PUM-based finite element method has also been employed in simulation of cracking in laminated composites. It has been especially used to develop solid elements that are capable of simulating through-thickness discontinuity to model delamination (see (de Borst and Remmers 2006; Wells et al. 2002) for instance). The important advantage of the enriched finite element formulation is that it eliminates the need to remesh when the crack propagates. By propagation of crack, the elements that are on the crack path get enriched with the discontinuous shape functions. The main shortcoming of this approach is the difficulty with its implementation in commercial FE codes. This is due to the requirement for introduction of new global degrees of freedom to the structure when the crack propagates. The method also requires information about the neighbouring elements for evaluation of the crack evolution law. 2.4 COMPARISON OF APPROACHES Table 2-1 presents a comprehensive comparison between the local continuum damage models, the continuum damage model enhanced with non-local averaging and the X-FEM approach. One can look at crack as the extreme case of a localized damage zone. It has been reported in the literature that the damage-mechanics-based continuum approach (if treated appropriately) 9 and fracture-mechanics-based discrete approaches lead to similar predictions when it comes to simulation of large-scale structures with localized damage (See (Mazars and Pijaudier-Cabot 1996) for instance). 2.5 SUMMARY An overview of the current approaches in simulating fracture in damage in solids, namely the discrete approach, continuum approach and embedded crack approach was presented in this chapter. Figure 2-2 presents a schematic chart that shows the various approaches to simulation of fracture and damage in solids. The numerical objectivity issues associated with local continuum damage models and the non- local averaging as a solution to these issues are discussed in detail in the next chapter. 10 Table 2-1. Capabilities and limitations of approaches for simulating fracture and damage. Method Feature Local Smeared Crack Methods (Classical Damage Theory) Non-local Damage Models The Extended Finite Element Method (X- FEM) and other PUM-based enriched FE formulations Type of Elements Local damage models work best with structured mesh and four-noded quadrilateral elements. No theoretical limitation on types of elements. No theoretical limitation on types of elements. Shape of Elements and Mesh Crack path and damage pattern are highly influenced by the orientation of the mesh. Therefore the mesh needs to be aligned with the crack direction if the direction is known in advance. If the mesh is fine enough, the predictions of non-local models in terms of damage direction, pattern and size are independent of mesh size and orientation. X-FEM is shown to be almost independent of mesh geometry. Limitations on mesh size There are limitations on size of elements. Since all the softening behaviour localizes in one row of elements, large strains will form in that row. The smaller the elements, the larger the strain in them. This sudden jump of strain in non-damaged elements and damaged row may cause instability in numerical solution. Employing the crack band model also puts a limit on the maximum element size to avoid local snap-back in the element. There is an upper limit on the mesh size dictated by non-local averaging method. In order to make the averaging effective, at least a few elements have to be within the averaging region. The mesh should be fine enough to capture the fracture process zone (FPZ) appropriately. The maximum mesh size would depend on the size of the FPZ. The length scale related to the length of the FPZ is implicitly included in the cohesive law that governs the behaviour of the embedded crack. Numerical Instability Due to softening nature of crack growth problem and chance of snap-back, all the methods may encounter numerical instabilities. The dynamic simulation with explicit time integration scheme is the method of choice for highly non-linear problems that include softening. Special solution methods (such as crack length control scheme) can be used to capture the global snap-back behaviour of the structure. See (Carpinteri 1989) and (Moes and Belytschko 2002) for further information. 11 Table 2-1. Capabilities and limitations of approaches for simulating fracture and damage., (cont.) Method Feature Local Smeared Crack Methods (Classical Damage Theory) Non-local Damage Models The Extended Finite Element Method (X- FEM) and other PUM-based enriched FE formulations Ability to model crack growth Crack growth can be modeled when its direction is known in advance. Generally first mode cracks or mixed mode cracks with restrictions on crack path (like delamination problem) can be modeled with smeared crack models. A damage model that is enhanced with a non-local regularization is capable of capturing the crack path and patterns. Since crack can be placed in any position in the element and inter-element compatibility conditions are met in X-FEM, crack growth can be modeled using this method without remeshing. Implementation in commercial codes. Can be incorporated into any code which accepts user-defined softening materials. User defined material models, generally get state of strain as an input and pass the state of stress as an output. The smeared crack material models are built in some codes such as LS-DYNA, ABAQUS and ANSYS. The non-local averaging requires passing of information between the elements and therefore the capability is not typically available for the user-defined materials in the commercial FE codes. LS-DYNA provides non-local averaging for some of its built-in damage and plasticity models. Since some new global degrees are freedom are introduced by X-FEM, and will be added dynamically during analysis, also since information is required from neighbour elements, implementation of X-FEM-based models in the commercial codes faces some difficulties. There has been some recent efforts to include some built-in X-FEM capabilities in codes such as ABAQUS and LS-DYNA. 12 Figure 2-1. A schematic representation of the finite element field using the methods of embedded discontinuity: (a) two neighbouring elements with a single weak discontinuity in each form a localized band; (b) the localized band is simulated as two weak discontinuities inside one element; and (c) a strong discontinuity inside an element. This figure is reproduced from Figure 1 in (Jirasek 2000). (a) (b) (c) 13 Figure 2-2. A schematic chart showing some of the current approaches to model fracture and damage in solids. Continuum Models Scalar Damage Discrete Models Remeshing Cohesive Interface Rotating Crack Microplane Fracture and Damage Isotropic MaterialsOrthotropic Materials Generalized Smeared Cracking Embedded Cracks Partition of Unity Method (PUM) X-FEM Combined Continuum- Discontinuous Approach Embedded Discontinuity Weak Discontinuity Strong Discontinuity VCCT 14 3. LOCALIZATION PROBLEM AND NON-LOCAL REGULARIZATION 3.1 INTRODUCTION Loading of notched-specimens made of quasi-brittle materials usually lead to formation of damage over a zone with a finite height. Test observations show that the height of damage formed in a self-similar damage growth in such materials is a material property and is almost independent of size and geometry of the structure. The numerical simulation of damage growth using classical continuum damage theory is unable to predict such damage patterns. These damage models predict formation of damage in a row of elements. Therefore finite element predictions using the continuum damage models show strong dependency on the size of the mesh. This so-called localization problem is not due to errors in the numerical model but is the manifestation of a fundamental flaw in the mathematical equations that govern the behaviour of the material experiencing damage. This chapter presents an overview of the localization problem and its remedies. 3.2 STRAIN-SOFTENING BEHAVIOUR Laboratory tests on well-designed notched specimens made of quasi-brittle materials (e.g. concrete or laminated composites), show gradual loss of the load carrying capacity of the specimen by increasing the crack opening. Figure 3-1 shows a typical load versus pin opening displacement (POD) curve measured from a tensile test on a notched specimen made from CFRP laminated composite. Numerical simulation of such phenomena in a continuum framework has always been desirable since it is more convenient to implement in FE codes and unlike the discrete approach, it does not require remeshing. The goal is to apply a suitable stress-strain law with softening behaviour at the material level such that when it is employed in a finite element simulation, the predicted structural response agrees with the experimental observations (Figure 3-2). 15 Material models with hardening behaviour or perfectly plastic behaviour are not capable of predicting a softening load-displacement response. It is only the material models with strain- softening behaviour that can lead to prediction of such a structural response. Historically, there are two distinct approaches that both lead to strain-softening models. The first one is the smeared crack approach where the effect of crack is smeared over a region with certain height and the cohesive traction-separation law is translated into a softening stress- strain relation. The family of crack-band models (e.g. (Bazant and Oh 1983)), fixed and rotating crack models (e.g. (Jirasek and Zimmermann 1998)) and the micro-plane models (e.g. (Bazant and Oh 1985)) are within this group. The second approach is the continuum damage mechanics approach that aims to capture the effect of material degradation on its mechanical properties such as its stiffness. By initiation and evolution of damage in the material, its stiffness gradually reduces until the material loses its entire load carrying capacity. These two approaches have a common outcome: they both assign a strain-softening behaviour to the material that is undergoing damage and fracture. 3.3 LOCALIZATION PROBLEM Computational modelling of structures with materials that exhibit softening behaviour faces a fundamental objectivity issue. By changing the size and orientation of the mesh, predictions of the finite element simulation changes significantly. In this chapter, the reasons behind this problem and the possible solutions are discussed. It is widely acknowledged that the continuum damage models lead to the localization phenomenon. The solution for a structure made of a progressively damaging elastic solid predicts localization of damage over a surface of zero thickness leading to the energy dissipated by the damage process being equal to zero. This would happen regardless of the shape of the softening part of the stress-strain curve. The causes and effects of localization have been discussed in detail in several articles pu lished in the 80’s on this subject (see (Belytschko and Bazant 1984), (Bazant et al. 1984) and (Bazant and Belytschko 1985) for instance). 16 From a mathematical point of view, when damage initiates in a strain-softening material, the partial differential equation governing the equilibrium of the quasi-static system loses its ellipticity, becomes ill-posed and the boundary value problem ceases to have a unique solution. In dynamic problems, the initial value problem of stress wave propagation loses its hyperbolicity when the material starts to soften. This implies that the stress waves cannot be transferred through the damaged zone. Additionally, within the damaged zone the solution to the dynamic problem becomes exponential leading to an unbounded velocity of the material points in this region (see (Sluys 1992) among others). The above problem manifests itself as spurious mesh-size dependency in the finite element solution. In such cases damage/crack pattern tends to localize to the smallest length scale of the model that is the height of one element. Therefore, with successive mesh refinement, damage localizes into a zone of zero volume and the numerical predictions fail to converge to a unique solution. Consequently, the global response of the system shows a strong dependency on the spatial discretization (see (Bazant et al. 1984; Belytschko and Lasry 1989; de Borst et al. 1993; Lasry and Belytschko 1988) among others). 3.4 REMEDIES TO THE LOCALIZATION PROBLEM To address the localization issue, various remedies have been proposed. Some of the remedies such as the Crack Band method offer a technique that can be used in the finite element solution to avoid the mesh-size dependency problem but they do not solve the underlying mathematical problem. Other remedies, such as the family of non-local methods, propose a more fundamental solution that addresses the localization issue in the governing equations. These methods prevent the localization by, implicitly or explicitly, introducing a length scale to the governing equations. These so-called localization limiters (Belytschko and Lasry 1989; Jirasek 1998; Lasry and Belytschko 1988) force the damage to grow in a zone with a finite width that is independent of spatial finite element discretization. Non-local averaging, explicit and implicit gradient and finally rate dependent damage models are among the localization limiters. 17 It was shown by Sluys in his PhD thesis (Sluys 1992) that employing any of the localization limiters preserves the hyperbolicity of the wave equation and thus prevents the localization problem. 3.4.1 Crack-band Approach Crack-band approach proposes the simplest remedy to the localization problem in the finite element analysis. Although crack-band does not solve the mathematical issue behind the localization problem, it can be used to prevent the mesh-size dependency in the finite element simulations. In this approach that was originally proposed by Bazant and Oh (Bazant and Oh 1983), the softening part of the stress-strain relation is scaled according to the height of the element to keep the energy dissipated during the damage progression consistent and independent of element height. Assuming that gf is the area under the stress-strain curve known as the fracture energy density, the fracture energy, Gf, can be written in the following form: dff hgG (3.1) where hd is the height of the damaged zone. Since in a local finite element analysis, damage localizes in a single row of elements, hd, would be equal to the height of that element, he. By scaling the fracture energy density, gf, according to the height of element as shown in Equation (3.2), the crack-band method makes the Gf fairly consistent and independent of the mesh size. Figure 3-3 shows a typical linear softening stress-strain curve and the scaling of the softening part of the curve according to the mesh height. fFE f e G g h (3.2) The crack-band model can be seen as a representation of cohesive crack in a smeared sense where the traction-separation law is translated into a stress-strain law. Opening of the crack faces is smeared over the height of an element and is represented by strain. While simplicity and ease of implementation in commercial codes are the advantages of the crack-band approach, this method faces many shortcomings as listed below. 18 - Crack-band method works best in structured meshes with rectangular elements. Triangular elements and unstructured quadrilateral elements bring complications to the concept of height of element in the analysis. - The crack-band method is not compatible with the higher order elements. Using fully- integrated higher order elements, one or more rows of integration points remain undamaged and unload (see (McClennan 2004) for instance). This not only complicates the calculation of fracture energy but may also result in stress locking in the element. - The damage pattern predicted using the crack-band approach does not represent the physical reality as its height is equal to the height of one element. Therefore, damage height becomes mesh-size dependent. This affects the predicted distribution of stress ahead of the damage tip and consequently affects the structural response. - Similar to cohesive crack models, this method is applicable to cases where in reality there is a distinct localization of damage and formation of a crack. Growth of damage and crack in a notched specimen under quasi-static loading is an example of such situation. In very fast dynamic events such as blast and high-velocity impact loading where the damage takes a more spatially distributed (diffuse) form and thus a well- defined path for the crack propagation does not exist, the applicability of the crack band concept becomes questionable (Bazant and Jirasek 2002). - The crack-band method cannot solve the mesh-orientation dependency problem associated with the local finite element solution as discussed in the next section. 3.4.2 Mesh Orientation Dependency of the Local Damage Models The damage and crack pattern predicted by the finite element analysis using classical continuum damage models, including the crack-band method, tend to follow the mesh path. This is due to the fact that predicted local strains and stresses at the integration points depend on the mesh geometry and orientation of the mesh. An example of the so-called orientation-dependency phenomenon is shown in Figure 3-4. In this example, a notched specimen made of an isotropic material is simulated using the LS- DYNA software. MAT_PLASTICITY_WITH_DAMAGE that is a built-in material model 19 with isotropic damage capabilities is employed in this simulation (see (Zobeiry et al. 2008) for more details). It can be seen that the crack pattern, instead of following the expected growth path along the direction of the notch, follows the mesh pattern. The mesh orientation dependency problem is also observed in simulation of dynamic problems. Figure 3-5 shows the damage pattern in plates subjected to lateral impact and blast load. The plates are assumed to be made of an isotropic quasi-brittle material and the same material model in LS-DYNA has been employed for the simulations. It can be seen that in both cases, the predicted damage follows the mesh direction that results in unrealistic damage patterns (see (Forghani and Vaziri 2009) and (Forghani et al. 2007) for more details). It can be concluded that the application of local damage models including the crack-band method is limited to quasi-static analysis where the crack path is known in advance and the mesh pattern is designed according to the expected crack pattern. 3.4.3 Rate Dependent Damage Models Introduction of strain-rate terms into the damage formulation can also overcome the localization issue. Rate dependency delays the evolution of damage at a material point and allows the neighbour points experience high enough stresses that diverts the growth of damage to those points. Needleman (Needleman 1988) showed that the introduction of rate dependency solves the issue of localization in shear band problem. Sluys and de Borst (Sluys and de Borst 1992) introduced a rate-dependent formulation for smeared cracking model to address the localization problem. It was shown in these articles that for the rate-dependent softening solids, the equilibrium equation for quasi-static loading remains elliptic. In dynamic cases, the wave speed remains real and hyperbolicity of the problem is preserved. It was also shown that the rate dependency, implicitly (and through introduction of a time scale), introduces a length scale to the governing equations. Another approach within the same family was proposed by Ladeveze and co-workers by introducing a delay to the damage evolution law that prevents instantaneous damage growth (e.g. (Ladeveze et al. 2000)). 20 Although the rate dependent damage models propose a relatively simple solution to localization problem, they are only applicable to materials that exhibit rate dependent behaviour. In polymer-based laminated composites, although the resin matrix shows rate dependent behaviour, the effect of the strain rate on fibres is typically negligible. According to (McClennan et al. 2003) there is no agreement within the research community on how extensive the rate effect is on the behaviour of laminated composites. The effect of inertia should not be confused with rate dependency. Due to the inertia effect, materials always exhibit higher effective stiffness when subjected to dynamic loads. 3.4.4 Non-Local Methods In a damage formulation enhanced with a non-local scheme, unlike the local formulation, the state of stress at a point does not only depend on the state of strain and history parameters at that point, but also depends on the state of these parameters in a finite neighbourhood of that point. There are two families of non-local formulations, namely, integral (averaging) formulations (Bazant and Pijaudier-Cabot 1988) and gradient-based formulations (Peerlings et al. 1998; Peerlings et al. 2002). Both non-local averaging and gradient-based non-local approaches, introduce a length scale to the mathematical representation of the physical system and thus prevent the localization problem. 3.4.5 Non-local Averaging Formulation (Isotropic Averaging) Non-local averaging is a popular approach among the non-local methods. In this method, a suitable parameter, α is averaged over a neighbourhood of the point as shown in Equation (3.3). ( ) ( )w d X X x X x (3.3) In this equation, is the average of and w is a weighting function such that ( ) 1.0w d X X x . Averaging is usually performed over a finite neighbourhood of the original point using a bell-shaped weight function as shown schematically in Figure 3-7. The 21 Gaussian distribution function or the inverse of the polynomials are often chosen as the weight function. It is critical to choose an appropriate parameter for averaging in a non-local formulation. By writing the stress-strain relation in the following general form, we can identify the various parameters that can be chosen for non-local averaging: (3.4) In this equation, the stiffness of the damaged material is written in terms of a single damage parameter. The discussion presented here is equally applicable to damage models with multiple damage parameters. The initiation and evolution of the damage parameter is governed by equivalent strain parameter (parameter #3 in Equation (3.4)) that is a function of the strain. Following are some of the requirements for an appropriate averaging scheme: - The averaging procedure should not affect the behaviour of the system in the undamaged elastic regime. Consequently it can only be applied to parameters that are involved in damaged stiffness calculations. As a result of this condition, the averaging cannot be applied to the stress and strain fields (parameters #1 and 5 in Equation (3.4)). - The averaging procedure should prevent the localization and make the predictions independent of the mesh size. It can be shown that the damage parameter(s) ω (parameter #2) is not a suitable variable for non-local averaging. The problem with averaging the damage parameter is due to its limited range 0 to 1. An element with a very high local strain level and local damage parameter of 1 with neighbouring elements that are not fully damaged, will experience an averaged damage parameter of less than 1. This will cause further transfer of stress through that fully damaged element and thus result in stress locking. The built-in non-local capability of LS-DYNA code has been used here to demonstrate this point when the damage parameter is chosen for averaging. In this example, a notched specimen is simulated using three different mesh sizes in front of the initial notch: 1mm, ( )eq σ D ε ε 1 2 3 4 5 22 0.5mm and 0.25mm. The built-in non-local capability of LS-DYNA in MAT_PLASTICITY_WITH_DAMAGE (MAT_81), gives the user the option of choosing between equivalent strains (parameter #3 in Equation (3.4)) and damage parameter (parameter #2) as the averaging parameter. MAT_81 is a built-in combined damage and plasticity model with non-local capabilities. In this numerical example, a linear softening stress-strain model is simulated using MAT_81. The simulations show that when equivalent strain is used as the averaging variable, the predictions in terms of load-displacement curves are almost independent of mesh size (see Figure 3-9). When the damage parameter, ω, is chosen for averaging and with the erosion option in LS- DYNA turned off, the structure does not exhibit any softening behaviour as shown in Figure 3-10 due to stress locking. By turning the erosion option on, LS-DYNA deletes the elements that are experiencing large strains. In this case the predictions show softening but the response depends on spatial discretization as well as the value chosen as erosion strain as shown in Figure 3-11. Sudden deletion of elements that are carrying stress despite of their large local strain causes abrupt jumps and noise in the predicted response of the structure. The strain parameter itself (parameter #4) may be a suitable parameter for the non-local averaging. However in the cases where different damage mechanisms are present (as is the case in laminated composites), the size and shape of the averaging zone may be different for each damage mechanism and therefore the equivalent strains (parameter #3) would be the better choice. 3.4.6 Gradient Non-Local Formulations The explicit and implicit gradient formulations are the two forms of the gradient-based non- local approach. Equations (3.5) and (3.6) show the explicit and implicit non-local equations respectively where the non-local equivalent strain is related to the local equivalent strain through partial differential equations (See (Peerlings et al. 2002) for instance). 2eq eq eqc (3.5) 2eq eq eqc (3.6) 23 In these equations eq is the local equivalent strain parameter and eq represents the non-local equivalent strain. In both equations c has a length squared unit and is related to the material characteristic length (damage height). In the explicit formulation as shown in Equation (3.5), the non-local equivalent strain is explicitly written in terms of local equivalent strain. The finite element model enhanced with the explicit gradient formulation requires the calculation of the Laplacian of the equivalent strain. Therefore it requires enforcement of higher order continuity of the displacement field in the finite element space. As it is shown in (Peerlings et al. 2001) and (Jirasek 2007), in the explicit gradient formulation, the non-local equivalent strain depends on the state of strain in the immediate neighbourhood of a point. Therefore unlike the implicit gradient and the non- local averaging formulations, the explicit gradient formulation only introduces a weak non- locality to the problem (Bazant and Jirasek 2002). It is shown in (Peerlings et al. 2001) that in some cases, this formulation leads to unrealistic predictions. In the implicit formulation (Equation (3.6)), a Helmholtz equation implicitly relates the non- local equivalent strain to the local equivalent strain. In the finite element implementation, the non-local equivalent strain is often treated as an independent variable. This variable is evaluated by transforming the Equation (3.6) into weak form and solving the coupled elasticity-non-local problem. It has been shown that the implicit non-local formulation is a special case of the non-local averaging. The implicit gradient formulation is equivalent with the integral form of non-local when Green’s function of the Equation (3.6) is chosen as the weight function, w. See (Peerlings et al. 2001) and (Peerlings et al. 2002) for further details. 3.4.7 Numerical Examples with Isotropic Averaging In this thesis, the non-local averaging has been employed to regularize the damage model and address the localization problem. The averaging formulation is chosen since it is more convenient for implementation in the existing finite element codes. Non-local averaging not only addresses the localization problem, but also helps to improve the predicted damage pattern. The numerical simulations presented in Chapter 3.4.2 are re- analyzed with non-local averaging enhancement and the predicted damage pattern is compared to the predictions of the crack-band method. The predictions presented in this chapter are all performed with LS-DYNA software using MAT_81 and its non-local capability. 24 Figure 3-12 shows the predicted crack pattern from simulation of a notched specimen under quasi-static loading. The mesh pattern was deli erately inclined in front of the notch. It’s shown in Figure 3-12 (a) that instead of growing horizontally, crack has followed the mesh pattern when local crack band model is employed. Applying the non-local averaging improves the predicted crack pattern as it’s shown in Figure 3-12 (b) and makes the crack path to follow the expected horizontal direction. As was shown in Section 3.4.2, the mesh orientation dependency problem is not limited to quasi-static problems. This problem can be observed when the local damage models are used to simulate behaviour of specimens under dynamic loading such as impact and blast load. Figure 3-13 shows the predictions of damage pattern on the surface of a plate that is subjected to low-velocity lateral impact. In these simulations, two different mesh patterns and sizes where used and local crack-band and non-local averaging models were employed in the simulations. Figure 3-13 (a) and (c) shows the predicted damage patterns when local crack band model is used. The predicted damage patterns in these cases are in the form of two localized band of damage (like folding lines) where orientation of these folding lines depends on the mesh orientation. The predicted damage patterns in these cases exhibit unrealistic dependency on the mesh size and do not agree with experimental observations. The predictions of the model enhanced with non-local regularization show similar damage pattern in the two different mesh patterns and are in a better overall agreement with the experimental observations. See (Forghani and Vaziri 2009) for more details on the specifications of the tests and the specimen and discussion about the numerical predictions. The predicted damage pattern of a plate subjected to lateral blast loading is presented in Figure 3-14. The predictions of crack-band model (Figure 3-14 (a)) show an unrealistic crack pattern that has followed the mesh direction. Non-local averaging has improved the predictions and resulted in a more realistic crack pattern as shown in (Figure 3-14 (b)). 3.5 SUMMARY A brief description of the localization issue was presented in this chapter. An overview of the so-called localization limiters was also presented and the non-local averaging method was discussed in detail. 25 In Chapter 6 of this thesis, an orthotropic averaging scheme that is especially designed for modelling damage in orthotropic material is proposed. 26 Figure 3-1 (a) A typical Overheight Compact Tension (OCT) specimen; and b) A typical load-POD curve extracted from test on an OCT specimen. 0 2 4 6 8 10 12 0.0 1.0 2.0 3.0 POD (mm) L o ad ( k N ) 27 Figure 3-2. A schematic, showing a finite element that uses strain-softening as an input and predicts load-POD response. 0 2 4 6 8 10 12 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 POD (mm) L o ad ( k N ) S tr e s s Strain 28 Figure 3-3. In crack-band method, the softening part of the stress-strain curve is scaled according to the height of the mesh. Figure 3-4. FE prediction of crack growth in the simulation of a notched specimen made from an isotropic material using MAT_81 local in LS-DYNA. fin er m esh ε σ gf 29 (a) (b) Figure 3-5. Predicted damage patterns using the crack-band approach on (a) a plate under lateral impact; and (b) a plate subjected to blast load. 30 Figure 3-6. A chart showing the various so-called localization limiters. Localization Limiters Non-local Approach Rate Dependent Models Cosserat Continuum Non-local Integral Formulation (Averaging) Gradient Non-local Formulations Explicit Gradient Formulation Implicit Gradient Formulation 31 Figure 3-7. The averaging zone ΩX around the point X and the bell-shaped weight function w employed for averaging the equivalent strains. 0 0.2 0.4 0.6 0.8 1 -0.5 -0.3 -0.1 0.1 0.3 0.5 X ΩX w 32 Figure 3-8, (a) Dimensions of the over-height compact tension specimen (OCT) (Kongshavn and Poursartip 1999); and (b) the applied force and the pin-opening displacement. a = 33 mm h = 1 0 4 m m d = 19.1 mm W = 80.6 mm 106 mm 8 4 .6 5 m m c = 38.7 mm Initial Notch F F POD (a) (b) 33 Figure 3-9. Force versus pin-opening displacement (POD) observed from simulations of an over height compact tension (OCT) specimen with different mesh sizes when the equivalent strain ( eq ) is chosen as the non-local averaging parameter. 0 4000 8000 12000 0.0 0.5 1.0 1.5 2.0 2.5 POD (mm) F o rc e ( N ) 1mm Eq. Strain 0.5mm- Eq. Strain 0.25mm- Eq. Strain 34 Figure 3-10. Force-POD predictions observed from OCT simulations with different mesh sizes when Damage Parameter (ω) is chosen as the non-local averaging variable. The erosion of elements with large strains was prevented. 0 4000 8000 12000 0.0 0.5 1.0 1.5 2.0 2.5 POD (mm) F o rc e ( N ) 1mm-Damage- NoErosion 0.5mm-Damage- NoErosion 0.25mm-Damage- NoErosion m . m . 5 m 35 Figure 3-11. Force-POD predictions observed from OCT simulations with different mesh sizes when Damage Parameter is chosen as the non-local averaging variable. The erosion strain was set to 150% in these simulations. 0 4000 8000 12000 0.0 0.5 1.0 1.5 2.0 2.5 POD (mm) F o rc e ( N ) 1mm 0.5mm 0.25mm 36 (a) (b) Figure 3-12. Predicted crack pattern in an OCT specimen with (a) local crack band approach and (b) non-local averaging model. 37 Figure 3-13. Predicted damage pattern in a plate subjected to lateral impact; (a) and (c) using local crack-band model; and (b) and (d) using non-local averaging. (a) (b) (c) (d) 38 (a) (b) Figure 3-14. Predicted damage pattern in a plate subjected to lateral blast load using (a) crack-band approach; and (b) non-local averaging. 39 4. SIMULATION OF DAMAGE IN LAMINATED COMPOSITES 4.1 INTRODUCTION One of the most difficult subjects in the numerical simulation of laminated composites is to predict the behaviour of structural components made of these materials under extreme loading conditions. This is due to the fact that the constituents of composite laminates, i.e. fibres and matrix, exhibit different mechanical behaviours. The difference between the behaviours of fibres and matrix and the inherent heterogeneity of the material, lead to complex and interactive modes of damage. Failure of structures made from laminated composites is a multi-scale event in nature. It starts with formation of micro-cracks in the matrix and leads to formation of meso-scale cracks that finally result in structural instability and failure. Unlike conventional construction materials such as concrete or metals with established theories and analytical methods for predicting their behaviour under ultimate loads, there is no comprehensive model for laminated composites that can predict the behaviour under various loading scenarios. The outcomes of the recent world-wide failure exercises (WWFE) showed that the existing failure theories face many difficulties predicting the damage behaviour of CFRP laminated composites (Hinton et al. 2002; Hinton et al. 2004; Kaddour et al. 2004; Soden et al. 1998). In this chapter an overview on computational modelling of damage and fracture in laminated composites is presented. 4.2 DAMAGE MODES Damage can be defined as permanent and irreversible changes in the microstructure of the material that degrades its mechanical properties such as its stiffness and strength. These micro-structural changes in laminated composites can be in the form of cracks in the composite constituents in various directions and orientations (Talreja 1985). Composite laminates are made by stacking unidirectional or woven layers on each other. Due to this structural configuration, two major damage modes, namely intra-laminar and inter- laminar damage modes exist in laminated composites. 40 Damage and cracks that are formed inside the plies are categorized as intra-laminar damage. Fibre breakage, matrix cracking and failure of fibre-matrix interface within the lamina are among the intra-laminar damage mechanisms. Inter-laminar damage or delamination on the other hand is the cracking between the layers of laminated composites. While continuum damage approaches are often chosen to simulate intra-laminar damage modes, delamination is usually simulated by employing the discrete crack and the cohesive crack concepts. 4.3 SIMULATION OF DAMAGE AT VARIOUS SCALES As mentioned earlier, formation and evolution of damage in laminated composites is multi- scale event. It initiates with micro-scale cracking of matrix that leads to formation of meso- scale delamination cracks and fibre breakage. Finally, by increasing the load, macro-scale cracks are formed. The existing damage models usually target one of these scales and study the damage behaviour of the material/structure at that scale. Depending on the target scale of the model, they are designed to capture the behaviour of the material (or material + structure) at different levels of details. Various scales in studying and simulating the mechanical behaviour of laminated composites are schematically shown in Figure 4-1. 4.3.1 Micro-scale Modelling The micro-scale models study the details of the interaction between the fibre and matrix at the scales up to a few hundred microns (see (Grufman and Ellyin 2008; Xia et al. 2000; Zhang et al. 2005) among others). Micromechanical models aim to predict the behaviour of the laminae based on the behaviour of its constituents (fibre and matrix) under various loading combinations. Micro-scale models have been successfully utilized to predict the mechanical properties of undamaged laminae such as the elasticity moduli and coefficients of thermal expansion, etc. When it comes to modelling damage at the micro-scale, two major approaches have been employed. In the first approach, that is the most popular one, a periodical distribution of the reinforcements is assumed. In this approach, the building block of the model is a repeated unit cell which consists of a fibre at its centre surrounded by matrix. Analytical closed form solutions or finite element simulation are employed to predict the behaviour of the unit-cell 41 under various loading scenarios. Fibre, matrix and the fibre-matrix interface are the constituents of a micro-mechanical unit-cell model (see (Xia et al. 2000), (Zhang et al. 2004) and (Xia et al. 2006) for instance). In the second approach, the morphology and distribution of the fibres in the laminae are into account. The size of the RVE in this approach depends mainly on the distribution of the fibres and can be up to a few hundred microns (see (Grufman and Ellyin 2007) and (Grufman and Ellyin 2008) for instance). Although the micro-scale models have been shown to be very successful in predicting the properties of the undamaged material and also in predicting the onset of damage, they are unable to predict damage evolution. This is due to the fact that the representative volume assumed for predicting the behaviour of the virgin material in a micro-scale model cannot be used as a representative volume for the extensively damaged material. In other words, in a damaged material, smearing cannot be performed at the same scale that is performed in the undamaged material. 4.3.2 Meso-scale Modelling At the next scale, the meso-mechanical models study the behaviour of the composite laminate at the ply-level (e.g. (Camanho et al. 2003; Ladeveze et al. 2000; Maimi et al. 2007; Talreja 1985)). At this scale, plies are the building blocks of the model. If the intention is to simulate delamination, cohesive interfaces will be placed between the neighbouring layers to connect them. In those cases where effect of delamination is ignored, the layers are assumed to be perfectly bonded together. The behaviour of individual plies and the interface are assumed to be the inputs for the ply- based models and the response of the laminate is predicted based on these given behaviours. The damage behaviour of the individual plies is either extracted from the tests performed on unidirectional laminates or from an underlying micro-mechanical theory. By using the properties extracted from tests on unidirectional laminates or extracted from micromechanical simulation of individual plies, the ply-based models implicitly assume that the behaviour of a ply in a laminate is independent of its neighbours. 42 Considering the fact that the standard tests suggested for characterization of unidirectional laminates often exhibit instability, the behaviour observed in such a test would not represent the behaviour of a ply in a multi-directional laminate. In a multi-directional laminate, neighbouring plies provide structural support to the damaged layers and bring alternative load paths to the structure. Therefore plies exhibit a different effective behaviour when placed in a multi-directional laminate. The other problem associated with the ply-based models is their computational cost. Since each and every layer of the laminate is simulated, the excessive computational cost (and time) makes the ply-based models impractical for simulation of real size structures. 4.3.3 Macro-scale Modelling Macro-scale models are the only models that can be used practically for simulation of large- scale structures. These models are designed to simulate the overall response of the laminate. The macro-scale models are typically not capable of predicting the details of the damage events in the layers. Numerical efficiency and also accurately predicting the behaviour of the structure in terms of macro-scale parameters are the main objectives of the macro-scale models. 4.3.4 Multi-scale Modelling The multi-scale nature of damage and failure in composite structures has been the motivation for the recently proposed multi-scale models. In the multi-scale models the simulation consists of modelling the behaviour of the structure/material at different scales (see (Chamis et al. 1996), (Chamis et al. 2000), (Ladeveze 2004), (Talreja 2006), (Gonzalez and LLorca 2006) and (Cox and Yang 2006) for instance). A multi-scale model consists of a structural-scale system that calls the smaller scale models to calculate the response of the material to the given macroscopic displacement/strain field. A major issue with the existing multi-scale models is that they are unable to predict the inherent characteristic length scales (like damage height and length of the fracture process zone) and there is a need for explicit incorporation of the length scale into these models (Bazant 2010). 43 4.4 SIMULATION OF DELAMINATION DAMAGE Delamination is defined as cracking along the interface between the layers leading to separation of adjacent layers. Delamination cracks form due to existence of inter-laminar stresses in laminated composites. Any form of discontinuity in the laminate such as holes, free edges, ply drops, etc. lead to formation of inter-laminar stresses. Minimizing the inter-laminar stresses by choosing an optimized stacking sequence to prevent delamination is an important objective in design of laminated composites. Delamination cracks may also occur as a consequence of matrix cracking (see (Delfosse and Poursartip 1997) and (Kashtalyan and Soutis 2000) for instance). Although there is a close correlation between delamination and matrix cracking, meso- and macro-scale numerical models usually treat them independently. Various approaches have been taken to simulate the delamination cracks in laminated composites. Among them are the discrete simulation of crack, crack-band approach and the embedded discontinuity approach based on Partition of Unity Method (PUM). Although some researchers have employed the continuum approach to simulate delamination ((Schipperen and de Borst 2001) and (Ladeveze et al. 2000) for instance), discrete approach is the most popular approach as it offers a better representation of the physical reality (See (Alfano and Crisfield 2001), (Borg et al. 2001), (Yang and Cox 2005) and (Forghani and Vaziri 2009)). Since the path of delamination growth is known in advance, the remeshing problem associated with the discrete approach is eliminated. The Evolution of delamination crack can be simulated using either linear elastic fracture mechanics (LEFM) or cohesive crack assumptions. When LEFM condition is assumed, the fracture energy release rate is numerically calculated using either the J-Integral method or the Virtual Crack Closure Technique (VCCT) (see (Krueger 2004) for instance). These methods have been implemented in some commercial FE codes such as ABAQUS. There are two major problems associated with the LEFM assumption in simulation of delamination cracking. Firstly, it requires an existing crack/notch and therefore LEFM is not capable of predicting the delamination onset in an undamaged material. The method also ignores the existence of the fracture process zone that exists in front of the delamination crack. 44 The cohesive crack method offers a more realistic representation of delamination cracks. The cohesive crack concept can be applied to interface elements (e.g. (Camanho et al. 2003)), cohesive contact formulations ((Espinosa et al. 2000) and (Forghani and Vaziri 2009) for instance) and enriched FE models ((de Borst and Remmers 2006; Remmers et al. 2003) among others). Interface elements with a softening traction-opening behaviour are often used in commercial finite element codes to simulate cohesive cracks. Some FE codes like LS-DYNA, offer a more general technique to simulate the cohesive cracks through the tie-break contacts. Forghani and Vaziri (Forghani and Vaziri 2009) employed the tie-break contact in LS-DYNA to simulate delamination in composite plates subjected to lateral impact and blast loading. They showed that the tie-break contact, if properly calibrated, is capable of simulating delamination cracks. Unlike the LEFM-based approaches, the cohesive crack method does not need a pre-existing crack or notch to determine the delamination initiation and growth condition. In the cohesive crack approach, initiation of the crack is usually determined using a stress-based criterion. Due to its narrow height, the cohesive interface has a very high initial stiffness. The high initial stiffness would place a very tight limit on the maximum time-step size in an explicit dynamic analysis (CFL condition) that would increase the computational cost. Local mass scaling of the interface can be employed to solve the time-step size problem. Another issue associated with the high stiffness is the introduction of local noise and oscillation in the FE solution that may lead to early initiation of delamination cracks. See (Forghani and Vaziri 2009) and (Forghani et al. 2007) for more details. 4.5 SIMULATION OF INTRA-LAMINAR DAMAGE Compared to delamination, the intra-laminar modes of damage are more complex in nature and often more challenging to simulate. This is due to the extreme orthotropy of the plies and the significant mismatch between the ply constituents, i.e. fibre and matrix. The significant difference between the behaviours of the fibres and the matrix and the inherent heterogeneity of the material leads to formation of various intra-laminar damage modes such as matrix cracking, interface failure, fibre breakage and fibre pullout. Formation of any of these damage modes depends on the state of stress/strain that the layer is experiencing. It also very much depends on the status of the neighbouring layers. Matrix cracking for example can cause 45 failure of fibre-matrix interface, leading to delamination and cause matrix cracking in the neighbouring layers. In spite of the fact that formation of damage in a layer is significantly influenced by its neighbours, the current ply-based models tend to simulate layers independently. In this approach, the behaviour of an isolated unidirectional laminate is taken as a representative of the behaviour of a layer in a multi-directional laminate (e.g. (Camanho et al. 2006; Ladeveze et al. 2000)). The main challenge in the ply-based framework is developing a model that is capable of predicting the damage behaviour of the lamina under different loading scenarios. Unlike the unidirectional laminates where intra-laminar damage take the form of localized cracks, in a multi-directional laminate intra-laminar damage forms in a more widely spread and diffused pattern. Therefore the continuum approach offers a better physical representation of intra-laminar damage compared to the discrete approach. 4.6 CONTINUUM DAMAGE MECHANICS FOR SIMULATION OF INTRA- LAMINAR DAMAGE Continuum damage mechanics has been widely employed to simulate the damage behaviour of plies. In this approach, the secant stiffness matrix of the damaged material is written in terms of a few damage parameters. A strain-based evolution law usually governs the damage growth in these models. See (Talreja 1985), (Talreja 1989), (Ladeveze 1995), (Lemaitre et al. 2000) and (Maimi et al. 2007) for instance. When it comes to orthotropic materials, the continuum damage approach faces two major obstacles: first one is the description of the stiffness of damaged material in terms of a set of state parameters (e.g. damage parameters) and the other one is formulation of the evolution law for those damage parameters. 4.6.1 Determining the Stiffness of a Damaged Material In establishing a damage model for laminated composites, the first challenge is to describe the stiffness of the damaged material in terms of the damage parameters. Extending the concepts defined in isotropic damage theory to an orthotropic material is not straightforward. As an example, employing the concept of strain equivalence for an orthotropic material leads to an asymmetric secant stiffness matrix. The secant stiffness describes the elastic behaviour of the damaged material and an asymmetric secant stiffness does not comply with the Betti- 46 Maxwell’s reciprocity law for elastic systems (see (Matzenmiller et al. 1995) and (Jirasek 2007) for instance). Matzenmiller and co-workers modified the asymmetric material compliance matrix by applying reduction to the Poisson’s ratios to come up with a symmetric compliance matrix (Matzenmiller et al. 1995). An alternative way to derive the stiffness matrix of the damaged material is to employ the so- called principle of energy equivalence that leads to a symmetric matrix (see (Chow and Wang 1987) and (Kennedy and Nahan 1996) for instance). However, there is a flaw in the assumptions behind this hypothesis as it does not take the unloaded portion of the material into account. See Appendix C for a more detailed discussion on this subject. Inspired by micromechanics of defects in solids, some researchers have chosen a different approach to define the damage parameters. In this approach, damage parameters are defined based on the size and orientation of the micro-scale cracks distributed over an RVE. See (Krajcinovic 1985) for the formulation proposed for initially isotropic materials and (Talreja 1985), (Talreja 1989) and (Lacy et al. 1997) for laminated composites. Micromechanical calculations are employed in this approach to determine the effect of distributed microcracks on the characteristics of the RVE such as its stiffness. This approach assumes that non- interactive micro-cracks are uniformly spread over the RVE in the damaged material. Such an assumption is valid only for early stages of damage development. At the later stages of damage development, when micro-cracks start to interact, coalesce and form meso-scale cracks, the assumptions are invalid and the approach would not be applicable (Lacy et al. 1997). In another attempt, which is different from the classical damage mechanics approach, Cusatis and co-workers have generalized Bazant’s microplane theory to simulate damage in the laminae (Beghini et al. 2008; Cusatis et al. 2008). 4.6.2 Damage Evolution Laws The second challenge in the continuum damage approach is writing the evolution law for the multiple damage parameters. Evolution laws that govern the condition for initiation and growth of damage parameters are usually written in terms of the strain state in the layer. Due to existence of various damage modes and interactions among them, defining the evolution laws has always been a challenge for continuum damage models. 47 A popular approach to formulate the evolution law is based on the hypothesis of maximum damage dissipation that is analogous with the associate flow rule in plasticity. In this approach, the evolution of each damage parameter is governed by its work conjugate that are called the energy release rate density parameters. Ladeveze and co-workers have employed this hypothesis to derive the evolution law for their damage model for intralaminar damage (Ladeveze 1995; Ladeveze et al. 2000). Although the hypothesis of maximum damage energy offers a straightforward framework for formulating the evolution law, it is a pure mathematical hypothesis that does not necessarily represent the physics behind the damage evolution. A more general approach motivated by the non-associated flow rule in plasticity defines a potential of dissipation function and damage evolution law is written in terms of the conjugates of each damage parameter with respect to the dissipation potential function. The choice of the dissipation potential function depends on the available test data and the purpose of the damage model (see (Krajcinovic 1985), (Lemaitre et al. 2000) and (Lemaitre and Desmorat 2005) for instance). The hypothesis of maximum damage dissipation can be seen as a special case where the Helmholtz Free Energy is chosen as the dissipation potential. 4.7 SUB-LAMINATE APPROACH As discussed before, the ply-based approach encounters two problems: high computational cost and the assumption of the independent behaviour of layers. The sub-laminate based approach that was first introduced by Williams et al (Williams et al. 2003), takes the sub- laminate as the building block of the laminated composite structure. This approach considers the sub-laminates as the representative volume and as the base level for constructing the damage model. The goal of the sub-laminate approach is to predict the overall behaviour of a laminate under a given loading condition. In this approach, in contrast with the ply-based approach, an element in the numerical simulation represents a sub-laminate. The material model assigned to the element is responsible for simulating the overall behaviour of the sub-laminate. Therefore there would be no need for through-thickness discretization and computational cost will be reduced significantly. The other advantage is that the interactions between the layers within the sub-laminate are implicitly contributing to the overall behaviour of the sub-laminate. This 48 approach does not intend and is not able to predict the details of events in the layers of the laminate. But the important objective of the sub-laminate approach is to capture the effect of damage in a sub-laminate on the macroscopic behaviour of the structure. This approach aims to accurately predict the mechanical features of a damaged structure such as its stiffness, load carrying capacity, stability and post-peak behaviour. Since the sub-laminate as a whole is assumed as the building block in this model, delamination within the sub-laminate cannot be explicitly simulated. If the delamination occurs in a small region within the representative volume that is surrounded by undamaged material, the effect of delamination can be implicitly accounted for in the overall behaviour of the sub-laminate. In situations where delamination is the dominant mode of failure, due to significant separation of the layers, the sub-laminate cannot be assumed as the building block and sub-laminate approach would not be applicable. It should be noted that in design of advanced composite structures, designers try to prevent large delaminations by choosing proper stacking sequences, avoiding thick plies and tapering at the free edges. Therefore large delamination cracks are not expected to form under in-plane loading conditions in well- designed composite laminates. As the damage behaviour of the sub-laminate depends on orientation, sequence and thickness of its constituent layers, each variant of the stacking sequence must be regarded as a distinct material for which the material model assigned to the sub-laminate needs to be re-calibrated. 4.8 SUMMARY A brief review of various approaches for simulating damage in laminated composites was presented in this chapter. It was noted that damage can be simulated at different scales. The current approaches to simulation of various damage modes were briefly presented. The ply-based approach and its advantages and shortcomings were discussed in detail. In the next chapter, the UBC COmposite DAmage Model (CODAM) that is a sub-laminate-based model will be discussed in detail. 49 Figure 4-1. A schematic showing various scales in studying and simulating the mechanical behaviour of laminated composites. Micromechanical • Fibre and Matrix Mesomechanical • Plies and interface between them Macromechanical • Sub-laminate • Laminate • Structure Courtesy of Leila Farhang Courtesy of Navid Zobeiry Unit Cell Ply-Based Sub-laminate Homogenized RVE σ ε FE model of the Structure 50 5. UBC COMPOSITE DAMAGE MODEL (CODAM) 5.1 INTRODUCTION The first generation of the UBC COmpoiste DAmage Model (CODAM) was introduced in the mid 90’s when there was a lack of a physically-based and yet computationally efficient model to simulate damage in laminated composites. CODAM was introduced as a damage model intended to be used in simulating large-scale composite structures. The original version of CODAM has been implemented as a user material model in the commercial explicit finite element code, LS-DYNA, and is currently being incorporated in LS-DYNA as a built-in material model. The ability of CODAM to predict the response of composites under a variety of loading scenarios has been demonstrated elsewhere (e.g. quasi-static loading of notched laminates (Floyd 2004; McGregor et al. 2008; Williams 1998; Zobeiry 2004); transverse impact of composite plates (Floyd et al. 2001; Williams 1998; Williams et al. 2003); axial impact of composite tubes (McGregor 2005; McGregor et al. 2007; McGregor et al. 2008)). A simplified version of CODAM has also been successfully implemented in the commercial implicit finite element code, ABAQUS. As discussed in detail in this chapter, the original formulation of CODAM faces two major problems. Firstly, it shows an unrealistic dependency on the choice of the coordinate system. Secondly, like other local damage models it suffers from the localization problem. In this chapter the problems associated with the original CODAM are discussed briefly. Then the second generation of CODAM, CODAM2, is introduced that addresses the objectivity issues of the original CODAM. 5.2 ORIGINAL CODAM 5.2.1 Background CODAM (Williams et al. 2003) is designed to predict the damage behaviour at the sub- laminate level. It accounts for the effects of matrix cracking and fibre breakage on the degradation of laminate’s stiffness. The loss of the laminate’s stiffness is written in terms of three damage parameters that represent the degradation of effective modulus in two perpendicular directions (in-plane principal directions) and the shear modulus in the same coordinate system. Figure 5-1 schematically shows an RVE of a cross-ply laminate subjected 51 to uniaxial strain. Figure 5-1(b) shows typical nominal stress-strain behaviour of the sub- laminate RVE under the applied deformation. 5.2.2 Formulation CODAM defines the stiffness matrix of the damaged laminate based on the three reduction coefficients as defined in Equation (5.1). 0 00 0 0 0 0 0 0 0 0 0 1 1 0 1 Ex Ey x xyEx x Ex Ey xy yx Ex Ey xy yx Ey yd Ex Ey xy yx Gxy xy R R ER E R R R R R E R R SYM R G Q (5.1) The reduction coefficients are defined typically as linear functions of the damage parameters. The formulation of CODAM to predict the stiffness of the damaged laminate is based on the formulation that (Matzenmiller et al. 1995) proposed for calculation of damaged stiffness in a lamina. The damage evolution law in the original CODAM is written based on the equivalent strain functions that are defined in the following generic form (Williams et al. 2003). 2 2 2 2 2 y y xy yzx x zx i i i i i i i i F K K L L S T U (5.2) The coefficients K, L, S, T and U are scalars to account for the relative effect of each strain component on the driving force for damage growth. The damage parameters (ωx and ωy) are then calculated as a function of corresponding equivalent strain function. In this model a general multi-linear relation between the damage parameter and the equivalent strain function can be defined (see (Williams et al. 2003)). The damage parameter for the in-plane shear is defined based on the damage parameters in x and y directions as shown in Equation (5.3). 52 2 2 s x y x y (5.3) shows a typical example of predicted stress versus strain for a classical progressive damage. By changing the ωvs. F and R vs. ωcurves one can get a variety of stress-strain behaviours such as instantaneous failure and brittle fracture. 5.3 OBJECTIVITY ISSUE The original formulation of CODAM shows unrealistic dependency on the choice of coordinate system and it is not capable of predicting the damage-induced orthotropy. Damage parameters and reduction coefficients in original CODAM are written in terms of strain components in the two perpendicular directions x and y that are aligned with the in- plane principal directions of the laminate (see (Williams et al. 2003)). As a result these two directions become the preferred directions in the damage model. The orientation dependency emerges when CODAM is used to simulate a quasi-isotropic laminate. Since in the elastic regime, the in-plane behaviour of a quasi-isotropic laminate does not exhibit any preferred direction, the choice of the material directions becomes arbitrary. Although the predicted elastic response of the laminate is independent of the choice of the material axes (x and y), its damage behaviour highly depends on the arbitrarily chosen directions. Figure 5-2 shows an example of the predicted load-displacement curves for a notched specimen made from a quasi-isotropic laminate. In one case, the material coordinate x is chosen to be parallel to the initial notch and in the second case it makes a 45º angle. It can be seen that there is a significant difference between the damage behaviours of these two systems. In the first model, the element in front of the notch experiences a large normal strain in its own material coordinate system whereas the second model experiences a combination of normal and shear strains as shown in Figure 5-3. Since damage formulation in CODAM is written in terms of the components of strain in the local coordinate system, the different local strain components lead to different equivalent strains and consequently different damage behaviours. 53 5.4 SOLUTIONS TO THE OBJECTIVITY ISSUE There are two major ways to solve the objectivity issue. The first option is writing the damage formulation in terms of strain-based parameters that are independent of co-ordinate system. Strain invariants and principal strains are examples of such parameters. This method is used mainly for simulation of damage in initially isotropic materials like concrete. The fixed crack (de Borst and Nauta 1985), rotating crack (Jirasek and Zimmermann 1998) and microplane (Bazant and Oh 1985; Bazant and Caner 2005) approaches are examples of these methods. Since these methods do not account for any preferred material orientation, they cannot predict the damage in composite laminates that have preferred directions. There have been some models that adopt the concept of microplane method in simulating damage in an orthotropic media (unidirectional ply) (see (Beghini et al. 2008; Cusatis et al. 2008)). This work is based on modal decomposition of the strain field into orthogonal components with respect to orthotropic stiffness of the lamina. This method extends the original decomposition of strain into volumetric and deviatoric components for isotropic materials that is used in the original microplane method. This method belongs to the category of ply-based approaches. The second option would be formulating the damage model in terms of strain components in the directions that play important roles during the damage process. In laminated composites, fibres are mainly responsible for carrying the loads and therefore the direction of the fibres are the material’s preferred directions. It should e noted that the principal axes in an undamaged laminate do not necessarily coincide with fibre directions (take the [0/45/90]s laminate as an example). A damage model that bases its formulation on the decomposed strain components in the direction of the fibres can potentially address the material objectivity issue that original CODAM formulation suffers from and would also be capable of capturing the damage- induced orthotropy in the laminate. The second generation of CODAM, called CODAM2, acknowledges the existence of some preferred directions in the sub-laminate and writes the damage formulation based on the components of strain in those important directions. Similar to any other sub-laminate based method, the sub-laminate is simulated as a whole and the basic iso-strain condition (i.e. all the layers experience a compatible strain field) is assumed here. 54 5.5 CODAM2 5.5.1 Background As it was discussed in the previous section, bringing the fibre orientations into the damage formulation provides the opportunity to address the dependency of the model on the choice of material’s coordinate system. CODAM2 acknowledges the fact that directions of the layers play an important role in the damage behaviour of the laminate. In CODAM2, the response of the laminate is determined through the combination of the effective behaviour of the constituent layers. Figure 5-4 schematically shows the effective behaviour of layers in a [0/45/-45/90]s quasi-isotropic laminate leading to an overall damage behaviour of the laminate as shown in Figure 5-4(b) The key points along the stress-strain curve shown in this figure are labelled numerically from 1 to 5 to help describe the sequence of damage events. Point #1 corresponds to matrix damage initiation in the 90º layer while point #2 corresponds to fibre damage initiation in the 0º layer. At point #3, the fibres of the ±45º layers start to degrade and at point #4, fibre damage in the 0º layer saturates. The degradation of overall stiffness to mitigate stress locking also coincides with point #4. Point #5 corresponds to the instant when damage is completely saturated and the sub-laminate loses its load carrying capacity. It should be noted that CODAM2 does not intend to predict the details of damage in the layers of the laminate at smaller scales. By assuming appropriate effective softening response for layers, the goal is to obtain the overall softening response of the laminate that represents its damage behaviour in a large structure. 5.5.2 Formulation The in-plane secant stiffness of the damaged laminate, A d , is written as the summation of the effective contributions of the layers in the laminate as shown in Equation (5.4). 1 n d T d k k k k k t A Τ Q Τ (5.4) where kΤ is the transformation matrix for the strain vector and d kQ is the in-plane secant stiffness of k th layer in the principal orthotropic plane. And tk is the thickness of the k th layer. 55 As discussed in Chapter 4, the majority of the methods that are proposed to describe the damaged stiffness of an orthotropic material are derived based on two hypotheses, namely the principle of strain equivalence and principle of energy equivalence. Instead of employing the above mentioned hypotheses, a physically-based and yet simple approach has been employed here to derive the damaged stiffness matrix. Two reduction coefficients, fR and mR , that represent the reduction of stiffness in the longitudinal (fibre) and transverse (matrix) directions have been employed. The shear modulus has also been reduced by the matrix reduction parameter. Following the approach proposed by (Matzenmiller et al. 1995), to achieve a symmetric secant stiffness matri , the major and minor Poisson’s ratios have been reduced by fR and mR respectively. This would lead to an effective stiffness matrix d kQ as shown in Equation (5.5). The reduction coefficients are equal to 1 in the undamaged condition and gradually decrease to 0 for a saturated damage condition. 1 12 2 12 21 12 21 2 12 21 12 0 1 1 0 1 f f m f m f m d m k f m m k R E R R E R R R R R E R R SYM R G Q (5.5) The stiffness reduction parameters are defined as linear functions of damage parameters as shown in Figure 5-5(b). 5.5.3 Damage Evolution The damage parameters are explicitly defined in terms of the maximum non-local equivalent strains that the material has experienced. It is common in the literature to adopt the hypothesis of maximum damage dissipation and write the damage parameters as a function of their work conjugates, i.e. energy release rates (e.g. (Ladeveze et al. 2000)). Parameters Yf,k and Ym,k defined later in Equations (5.22) and (5.23) are the work conjugates of damage parameters related to longitudinal and transverse directions in layer k, respectively, and represent the energy dissipation per unit volume due to the advancement of damage. 56 Although from a mathematical point of view employing this hypothesis is desirable, it does not represent the physics of damage growth properly as the work conjugates of damage parameters (Yf,k and Ym,k ) do not have a physical meaning. Due to the interactive form of the driving forces and their dependency on other damage parameters as shown in Equations (5.22) and (5.23), calibration of a damage model written in terms of Yf,k and Ym,k is also a complicated task. The hypothesis of maximum dissipation also leads to the path independency of the dissipated energy implying that the energy dissipated to arrive at a state of damage would be independent of the path taken to get to that damage state. The path independency of damage is not physically justifiable for a material as complex as a laminated composite. The model proposed here does not follow the hypothesis of maximum damage dissipation. Instead, the driving forces for damage are defined in a simple, yet physically based form. The equivalent strain function that governs the fibre stiffness reduction parameter is written in terms of the longitudinal normal strains as shown in Equation (5.6). The equivalent strain function that governs the matrix stiffness reduction parameter is written in an interactive form in terms of the transverse and shear components of the local strain as shown in Equation (5.7). , 11, ; 1... eq f k k k n (5.6) 2 2 12, , 22, 22,( ) ; 1... 2 keq m k k ksign k n (5.7) The sign of the transverse normal strain plays a very important role in initiation and growth of damage since it indicates the compressive or tensile nature of the transverse stress. Therefore, Equation (5.7) is written such that the equivalent strain for matrix damage carries the sign of the transverse normal strain. Within the framework of non-local strain-softening formulations adopted here, all damage modes, be it intra-laminar (i.e. fibre and matrix damage) or overall sub-laminate modes (locking or delamination discussed in Sections 5.5.5 and 5.5.6) are considered to be a function of the non-local (averaged) equivalent strain defined as: 57 ( ) ( ) ( )eq eq w d X X x X x (5.8) where the subscript denotes the mode of damage: fibre ( f ) and matrix ( m ) damage in each layer, k, within the sub-laminate or associated with the overall sub-laminate, namely, locking ( L ) and delamination ( D ). Thus, for a given sub-laminate with n layers, eq and eq are vectors of size 2n + 2. In Equations (5.8), X represents the position vector of the original point of interest and x denotes the position vector of all other points (Gauss points) in the averaging zone denoted by . In classical isotropic non-local averaging approach, this zone is taken to be spherical (or circular in 2D) with a radius of r. The parameter, r, which affects the size of the averaging zone, introduces a length scale into the model that is linked directly to the predicted size of the damage zone. In the orthotropic averaging scheme proposed in Chapter 6, the averaging is performed over an ellipse with a longitudinal radius of r1 and transverse radius of r2. The non-local averaging is performed on the equivalent strain functions related to each damage parameter. Therefore, separate averaging is applied to each damage mechanism with a bell-shaped weight function, wα. The damage parameters, , are calculated as a function of the corresponding averaged equivalent strains. In CODAM2 the damage parameters are assumed to grow as a hyperbolic function (Figure 5-5(a)) of the damage potential (non-local equivalent strains) such that when used in conjunction with stiffness reduction factors that vary linearly with the damage parameters (Figure 5-5(b)), they result in a linear strain-softening response (or a bilinear stress-strain curve) for each mode of damage. The explicit form of this relation is shown in Equation (5.9) below: ; for 0 eq i s eq i eqs i (5.9) where superscripts i and s denote, respectively, the damage initiation and saturation values of the strain quantities to which they are attached. 58 Damage is considered to be a monotonically increasing function of time, t, such that max t (5.10) where t is computed from Equation (5.9) for the current time (load state), and represents the state of damage at previous times t . The stiffness reduction coefficient or normalized modulus (ratio of the damaged to undamaged modulus of the material), R, is defined here as a linear function of the damage parameters: 1 ; for = , , 1 ; for = f m L R D (5.11) where 1 is a constant that defines the rate of modulus loss with the delamination damage, D and effectively determines the maximum reduction of the target stiffness due to formation of delamination, i.e. 1 1DR (see Figure 5-5(b)). See Section 5.5.6 for how the parameter can be estimated. 59 5.5.4 Thermodynamic Aspects Damage models intend to simulate the behaviour of the material subjected to irreversible degradation. Any damage model has to be thermodynamically objective and comply with the energy-related conditions: firstly the recoverable strain energy has to be non-negative because the material cannot store negative potential energy. And secondly the rate of energy dissipation has to be non-negative. 5.5.4.1 Recoverable Strain Energy The density of recoverable strain energy, ψ, in a damaged laminate can be written as shown in Equation (5.12). 1 1 1 1 : 2 2 2 n T d T d k k k k k t t t σ ε ε Α ε ε Q ε (5.12) where t is the total thickness of the sub-laminate. The free energy is written in terms of strain and reduction parameters (R parameters) as the state variables. One can alternatively choose strain and damage parameters as the state parameters. It should be noted that in this chapter, the material is assumed to be under isothermal conditions where the change in the temperature is assumed to be negligible. The strain energy stored in each layer has to be non-negative to ensure that the laminate’s reversible strain energy is non-negative for any given strain field. By Transforming strain into the eigenvectors space, one can re-write the reversible strain energy in each layer as 3 2 1 1 , 1.. 2 k i i i k k n (5.13) where λi represents the i th eigenvalue of the secant stiffness matrix of layer k (i.e. d kQ ), and i ~ is the component of strain vector in the direction of the i th eigenvector defined as :i i ε φ where φi is the normalized i th eigenvector and the operator “:” represents the inner product of the vectors/tensors. As long as all the eigenvalues of the secant stiffness matrix are non- negative, the reversible strain energy remains non-negative for any given state of strain. Determinant of d kQ can be written as a function of the reduction coefficients as shown in the Equation (5.14). 60 2, 12 , , 1 2 , , 12 22 12 211 m kd k f k m k f k m k f m R G R R E E R R E R R Q (5.14) Knowing that the reduction coefficients are bounded between zero and one ( 10 fR and 10 mR ), the determinant of d kQ is non-negative for any given values of Rf and Rm. The eigenvalues of d kQ can be written in terms of the reduction coefficients as follows. 12,66,,1 GRQR kmkmk (5.15) 2 2 , 1 , 2 , 1 , 2 , , 12 2 2, , , 12 21 4 2 1 f k m k f k m k f k m k k f k m k R E R E R E R E R R E R R (5.16) 2 2 , 1 , 2 , 1 , 2 , , 12 2 3, , , 12 21 4 2 1 f k m k f k m k f k m k k f k m k R E R E R E R E R R E R R (5.17) Since the determinant of d kQ is non-negative and the first two eigenvalues can be shown to be non-negative, the third eigenvalue has to be non-negative. The non-negativeness of the eigenvalues is a sufficient condition to make the reversible part of the strain energy remain non-negative for any given state of strain. 5.5.4.2 Rate of Energy Dissipation In the absence of temperature change, the rate of mechanical energy dissipation in the damage process, d , is written as the difference between the rate of the external work and the rate of change of the recoverable strain energy as shown in Equation (5.18). d , : : : : ε R σ ε σ ε ε R ε R (5.18) where R is a vector that contains the reduction parameters related to matrix and fibre damage for all the layers and it’s written as 61 ,1 ,1 ,2 ,2 , , ... ... m f m f m n f n R R R R R R R (5.19) Stress is defined as the work conjugate of strain ( ε σ ), and therefore the rate of energy dissipation can be simplified as follows: d , , 1 , , : n f k m k k f k m k R R R R R R (5.20) By defining the energy release rates, Yj, as the work conjugates of reduction parameters ( j j Y R ), Equation (5.20) can be re-written as: d , , , , 1 n f k f k m k m k k Y R Y R (5.21) In compliance with the second law of thermodynamics, it is necessary for a damage model to exhibit a non-negative rate of mechanical energy dissipation i.e. d 0 (e.g. (Lacy et al. 1997)). This condition is called Clausius-Duhem inequality or Kelvin inequality. See (Lemaitre and Desmorat 2005) and (Jirasek 2007) for more information on thermodynamics of irreversible processes. Since in an irreversible damage process, the values of reduction coefficients do not increase, the rate of change of the reduction parameters is non-positive (i.e. 0iR ). Thus having the energy release rates ( iY ) positive for all the damage mechanisms would be a sufficient condition to meet the Clausius-Duhem inequality. The energy release rates corresponding to fibre damage and matrix damage of the k th layer can be written as shown in Equations (5.22) and (5.23) respectively. 62 2 1 1 , 21 2 , 2 , , , 12 122 1 m k f k f k f k m k E R Y R R R (5.22) 2 2 2 12 1 2 , 12 122 , , , 12 12 1 22 1 f m k m k f k m k E R Y G R R R (5.23) According to these equations, the rates of energy dissipation during the damage growth for both matrix and fibre damage mechanisms are non-negative, therefore the proposed damage formulation satisfies the Clausius-Duhem inequality. The total energy dissipation rate in the structure is written as shown in Equation (5.24) below. : d d d dW d d R R ` (5.24) where Ωd represents the damaged area. It should be noted that although the density of the energy dissipation rate was proven to be non-negative, due to localization problem that was discussed in Chapter 3, in a finite element simulation using local damage model, the size of the damaged zone depends on the finite element discretization and therefore the predicted dissipated energy would be mesh-size dependent. By introducing length scales through a properly formulated non-local damage model, the localization problem can be overcome thus resulting in an energy dissipation rate that is independent of FE discretization and does not violate the thermodynamics principles. 5.5.5 Stress Locking and Transition to Scalar Damage In the context of this thesis, stress locking refers to the situation where a significantly damaged material is still capable of transferring stress. This issue has been observed in damage models that are capable of simulating damage-induced orthotropy such as the rotating crack method (see (Jirasek and Zimmermann 1998) for instance). In the proposed damage model, there are probable situations where a widely open crack is capable of transferring traction. An example of stress locking situation is demonstrated here. In Figure 5-6, two neighbouring elements on the crack path passing through a quasi-isotropic [45/90/-45/0]s laminate are schematically shown. The two Mohr Circles show the 63 corresponding strain states in these elements. In the left element, 0º and -45º fibres experience large strains that lead to damage of these fibres. The 90º and 45º fibres do not experience high strains and they remain undamaged. In the right element on the other hand, due to its strain state, the 0º and 45º fibres are extensively damaged and 90º and -45º are remained undamaged. The existence of undamaged -45º in the right and 45º in the left elements and the 90º in both forms a truss that is capable of carrying the load in the y direction. To overcome the stress locking problem, an overall damage mechanism has been employed to degrade the load carrying capacity of material in all the directions. This mechanism is activated when the material is subjected to high level of strain. The overall damage mechanism acts on all the layers in the laminate and degrades all the components of stiffness of the laminate simultaneously. This damage mechanism is written as a function of maximum principal strain that is an indicator of the crack opening regardless of its direction. 5.5.6 Edge Delamination In the sub-laminate approach, due to the inherent assumption that a sub-laminate as a whole is considered as the representative volume, delamination cannot be simulated explicitly by introduction of cracks between the layers. Rather, only the overall effect of delamination on the in-plane stiffness (and/or bending stiffness) of the laminate can be taken into account. There have been many attempts to predict the effective stiffness of the laminate in the presence of delamination. O’Brien (O'Brien 1981) proposed a method based on linear elastic fracture mechanics (LEFM) to predict the onset of edge delamination and its effect on the laminate’s stiffness. Kashtalyan et al. (Kashtalyan and Soutis 2000) introduced an analytical method to estimate the stiffness of a cross-ply laminate in the presence of matrix cracking and delamination strips. CODAM2 has been employed to predict two of the test cases for the third round of World Wide Failure Exercise (WWFE-III). In these test cases, as discussed in detail in Chapter 9.5, the edge delamination plays an important role and therefore O’Brien’s method to simulate the effect of delamination on the in-plane laminate stiffness has been added to CODAM2. O’Brien (O'Brien 1981) uses LEFM energy release rate arguments to derive the following equation for the critical in-plane strain corresponding to the onset of delamination: 64 D 0 * c c n G t E E (5.25) where nD is the number of delaminations initiated at the edge of the laminate, t is the laminate thickness, E 0 is the initial effective stiffness of the sub-laminate (parallel to the edge) and E * is the sub-laminate’s effective stiffness when the delaminated layers are completely separated, and Gc is the critical fracture energy for edge delamination. Depending on the fracture mode (Modes I, II or mixed-mode), Gc varies between GIc and GIIc. Once again, in the numerical analysis, the maximum principal strain is used as an overall measure of the magnitude of the in-plane strain experienced by the laminate and compared with the threshold value given by Equation (5.25). Near the free edges, the maximum principal strain represents the component of strain parallel (or tangential) to the edge and when it reaches the critical value, the stiffness of the laminate is reduced to E * by setting the parameter β in Equation (5.11) equal to * 1 E E .The resulting stiffness reduction factor, RD = E * / E 0 , is then applied to all components of the sub-laminate’s in-plane stiffness matrix. Figure 5-7(a) shows schematically the reduction of stiffness due to delamination in a one- dimensional loading situation. A sudden drop in the stiffness of the laminate may cause numerical difficulties, therefore in the numerical implementation the stiffness drop occurs gradually as shown in Figure 5-7(b). The equivalent strain that drives the delamination damage, eq D is taken to be equal to the maximum principal strain in the laminate. When overall stiffness reductions RL (for controlling the numerically-induced stress locking problem) and RD (for simulating the effect of free-edge delaminations) are considered simultaneously, the final form of the reduced in-plane stiffness matrix of the damaged sub- laminate (Equation (5.4)) can be written as: 1 n d T d L D k k k k k R R t A Τ Q Τ (5.25) It should be noted that the method presented here does not offer a general solution for simulation of delaminated composites and it is only applicable to loading scenarios where overall laminate stiffness reductions stemming from free edge delaminations occurs. 65 5.6 SUMMARY In this chapter the original CODAM and its shortcomings were discussed briefly. The second generation of the CODAM (i.e. CODAM2) was introduced in detail. Thermodynamic aspects of the proposed damage model were discussed. The next chapter introduces an orthotropic non-local averaging formulation that is specially formulated for orthotropic damage models. 66 (a) (b) Figure 5-1. (a) Schematic showing a RVE of a cross-ply laminate at three stages: (1) before the damage starts (2) formation of microcracks mainly in the matrix (3) microcracks merge and form a through crack in the laminate. (b) schematic nominal stress vs. nominal strain behaviour of the RVE under applied deformation. x1 1 2 3 xE x2 x3 Matrix cracks, delamination and some fibre breakage Complete failure RVE 2 x2 3 x3 Elastic 1 thickness w x1 hc strain-softening model 67 Figure 5-2. Predicted response of an OCT specimen with a quasi-isotropic lay-up. In one case, the local axes are aligned with the crack and in the other one they make a 45º angle. 0.00 5.00 10.00 15.00 20.00 25.00 30.00 0.00 0.10 0.20 0.30 0.40 Opening F o rc e 1 2 68 Figure 5-3. An element under uniaxial load assuming that the material’slocalCSis(a) aligned and (b)ma es angle. 1 2 1 2 Model 1 Model 2 69 (a) (b) Figure 5-4. A schematic showing (a) the effective strain-softening behaviour of layers in the laminate being combined to result in (b) the effective strain-softening behaviour of laminate. 0 200 400 600 800 0.00 0.02 0.04 0.06 0.08 Strain S tr e s s ( M P a ) 1 2 3 4 5 dg f S tr e s s ( M P a ) σ ε σ ε σ ε σ ε σ ε σ ε 70 (a) (b) Figure 5-5. (a) A typical damage parameter defined as a hyperbolic function of the corresponding non-local equivalent strain, and (b) the reduction parameter as a linearly decreasing function of the damage parameter. R i s eq 1 11 71 Figure 5-6. Schematic showing a stress locking situation in a widely opened crack. The two elements shown are on the crack path and subjected to large opening. The assumed strain states are shown on the corresponding Mohr circles. The undamaged fibres in these neighbouring elements form a load carrying path. ε11 x y 0 -4545 90 0 -45 45 90 ε11 0 45 -45 90 12 2 12 2 72 (a) (b) Figure 5-7. Schematic showing (a) reduction of laminate stiffness due to delamination, and (b) the gradual drop in the laminate stiffness as implemented numerically. E 0 E * E 0 E * s D i D c 73 6. ORTHOTROPIC NON-LOCAL AVERAGING 6.1 INTRODUCTION Composite laminates consist of strongly orthotropic constituents (plies) with very different behaviours along the direction of the fibres and in the transverse direction. The non-local averaging schemes were originally introduced to simulate the damage in isotropic materials such as concrete. When it comes to simulation of damage in laminated composites, although using an isotropic averaging scheme addresses the numerical objectivity issues, associated with the mesh size and mesh orientation dependencies, predictions of an orthotropic damage model with isotropic non-local averaging does not properly represent the physical reality. As shown in this chapter, CODAM2 equipped with isotropic averaging cannot accurately predict the pattern and direction of crack growth in highly orthotropic laminates. To address the problems, an orthoptropic averaging scheme is introduced in this chapter. The idea behind the orthotropic averaging is motivated by the physics of damage growth in a fibre reinforced medium. In this chapter, the formulation of the proposed orthotropic averaging is presented in detail. Numerical examples are presented to demonstrate the capability of the new formulation in predicting damage pattern and direction of crack growth in highly orthotropic media. 6.2 MOTIVATION The following hypothetical example demonstrates the need for a non-local scheme that acknowledges the existence of reinforcement in a certain direction in the material. Consider a fibre-reinforced composite material. Assume two situations where a material point (A) is influenced by a remote damaged area exhibiting large strains due to cracking as shown in Figure 6-1. In the first case (Figure 6-1(a)), the line connecting the point A to the damaged area (marked by D1) is parallel to the direction of fibres and in the second case (Figure 6-1(b)), the damage has to cross the fibres to reach the point A. The presence of the fibres, as the reinforcement, makes these two situations very different. The damage in case (a) has more influence on the behaviour of the material point A than the damage in case (b) does. Figure 6-2(a) schematically shows the averaging area when a classical isotropic averaging is used. In this case, both damaged zones D1 and D2 that have similar relative distances to point 74 A, will make similar contributions to the averaged effective strain at point A. Therefore by employing the isotropic averaging scheme, initiation and growth of damage at point A would be equally influenced by these damage zones. On the other hand, an averaging scheme with an orthotropic weight function, as shown in Figure 6-2(b), could assign more weight to D1 than D2. 6.3 FORMULATION In the formulation proposed here, an orthotropic weight distribution is employed to perform the averaging on the equivalent strains corresponding to matrix cracking and fibre breakage. Since the layers in the laminated composite are aligned in different directions, a different weight function is needed to perform the non-local integration in each layer. An isotropic averaging scheme is still used for the overall damage mechanism of the laminate. Here, the weight function for orthotropic averaging is written in terms of the effective distance (radius) as defined in Equation (6.1). 2 2 1 2 1 1 2 :eff d d r r r (6.1) where d1 and d2 are the components of the vector that connects the material point A to its neighbour B in the x1 and x2 directions respectively as shown in Figure 6-3. In this equation, r1 is the non-local radius along the fibres direction and r2 is the non-local radius in the transverse direction and r1>r2. The weight function is written in terms of the effective radius (ρ eff ) as shown in Equation (6.2). 2 2 1 1 eff effw r (6.2) An example of the distribution of weight function is plotted in Figure 6-4. In this example, r1 is chosen to be twice r2 and therefore the weight function decays twice faster in the x2 direction. Unlike the isotropic averaging scheme, the method proposed here introduces two length scales, r1 and r2, to the system. The predicted height of damage in the laminate would depend 75 on these two averaging radii as well as the lay-up and the loading configuration. The effects of lay-up and contributions of the averaging radii on the predicted damage height are discussed in the numerical examples presented in the Section 6.6. 6.4 IMPLEMENTATION CODAM2 has been implemented in the finite element code, OOFEM (Patzak and Bittnar 2001). OOFEM is an open-source code distributed under GNU General Public License. This highly modular code includes the essential finite element procedures and provides an excellent platform for implementation of advanced material models. Unlike the commercial FE codes, the open-source nature of OOFEM allows the developer to make enhancements to specialized elements and constitutive theories that require access to parts of the code other than element and material routines. The built-in isotropic averaging capability has been extended to perform orthotropic averaging schemes required in CODAM2. The details of the algorithm and the objects and classes added to the object oriented code are presented in Appendix D. The visualization of the results is performed using the open-source visualization code ParaView . All the numerical examples of CODAM2 and the orthotropic non-local averaging presented here were analyzed using OOFEM. 6.5 TREATMENT OF BOUNDARIES One of the subjects on non-local averaging that still is controversial is how to apply averaging for elements that are close to the boundaries of the structure. There are different suggestions to treat the averaging near the boundaries that among them are: averaging using the original weight function, re-scaling the weight function to make sure that the integral of the weight function over the averaging zone is equal to 1 and also averaging with local correction. (Jirasek et al. 2004). The current implementation of non-local in OOFEM takes the scaling the weight function approach to handle the averaging near the boundaries. By dividing the averaged equivalent strain by the integral of weight function as shown below, the averaging is re-normalized near the boundaries. 1 ( ) ( ) ( )eq eq w d W X X x X x (6.3) 76 where ( ) ( )W w d X X X x (6.4) Figure 6-13 schematically shows the scaling of the weight function near the boundaries in a 1D case. It can be seen that the weight function is scaled up and normalized to make its integral equal to 1. 6.6 PREDICTION OF DAMAGE/CRACK PATTERNS 6.6.1 Unidirectional Laminates The proposed orthotropic averaging scheme improves the prediction of the size and pattern of damage and size significantly. Following numerical examples demonstrate the improvements that the orthotropic averaging can make in predicting the crack pattern in unidirectional laminates. As the first example, consider a notched specimen made from a unidirectional laminate. The dominating damage mode in this case is formation of splitting cracks parallel to the fibres. The numerical simulation of a notched specimen made from a unidirectional laminate with a 0º angle (i.e. fibres aligned with the loading direction and perpendicular to initial notch) has been carried out using both isotropic and orthotropic averaging schemes. Figure 6-5 shows the predicted crack patterns. In these analyses, the averaging radius along the fibres, r1, was kept constant and equal to 2mm. Three different values of transverse averaging radius, r2, were chosen to examine the effect of this parameter on the predicted crack path. It can be seen from this figure that reducing r2 improves the predicted crack pattern significantly and enables the crack grow along the direction of the fibres as expected from a physical standpoint. The second example looks at the predicted damage pattern in a notched specimen made from an off-axis unidirectional laminate. In this case the fibres make a 60 º angle with the direction of loading. The damage here consists of matrix cracks that form parallel to the direction of the fibres. Figure 6-6(a-c) shows the predicted crack patterns using a local analysis, non-local analysis with isotropic averaging and non-local analysis with orthotropic averaging 77 respectively. In the latter case the transverse averaging radius was chosen to be a quarter of the longitudinal averaging radius. It can be seen that the predicted local damage pattern is completely influenced by the mesh orientation and grows along the direction of the mesh (which in this case is horizontal in front of the crack tip). Clearly the predicted damage pattern obtained from local analysis is physically unacceptable and is another evidence for inapplicability of local analyses. The classical non-local damage with isotropic averaging scheme results in a significant improvement in the predicted crack pattern as shown in Figure 6-6(b). The predicted direction of damage growth is almost parallel to the direction of the fibres. The problem that the isotropic averaging faces is the predicted height of damage that is unrealistically large. This is contrary to the fact that the matrix cracking parallel to the fibres is expected to propagate in a very narrow and localized band. Employing a relatively small transverse radius in the orthotropic averaging scheme can help to simulate a narrow matrix cracking pattern. Figure 6-6(c) shows the predicted damage pattern when an orthotropic averaging scheme is employed. It can be seen that in this case, with a proper choice of the averaging radii, a more realistic prediction of damage growth in terms of both direction and size (height) can be observed. As shown schematically in Figure 6-7(a), when the isotropic averaging scheme is employed, the predicted damage height depends on the non-local radius as it is the only length scale introduced to the system. In the orthotropic averaging scheme, on the other hand, depending on the lay-up and crack growth direction, the height of damage may be controlled by the appropriate selection of the major and minor averaging radii. When simulating unidirectional laminates in which matrix damage propagates parallel to the direction of the fibres, it is the minor averaging radius that determines the height of damage as shown schematically in Figure 6-7(b). 6.6.2 Multi-directional Laminates Figure 6-8 schematically shows the averaging regions that lead to different damage heights, h, in unidirectional and multi-directional laminates. Figure 6-8(a) shows the predicted damage pattern for an off-axis unidirectional laminate where matrix cracks grow parallel to the fibres 78 with the damage height mainly depending on the transverse averaging radius, r2. Figure 6-8 (b) shows the predicted damage pattern of a notched specimen made of a [0/45/90/-45]s quasi- isotropic laminate. In this case, both averaging radii contribute to the damage pattern but the height of damage is dominated by the longitudinal averaging radius, r1. To investigate the effect of the minor non-local radius on the predicted damage patterns in multi-directional laminates, the simulation of crack propagation in a quasi-isotropic laminate is carried out using three different values for r2 while keeping r1 constant. Figure 6-9 shows the predicted damage pattern with (a) isotropic averaging (r1=r2=2mm), (b) orthotropic averaging with r1=2mm and r2=1mm, and (c) orthotropic averaging with r1=2mm and r2=0.5mm. It can be seen that all three predicted damage patterns are quite similar and that they are not significantly affected by the transverse averaging radius. 6.6.3 Relation between Fracture Energy and Non-local Averaging Radii Similar to the damage pattern, the energy dissipated in the damage process depends on the averaging radii as well as lay-up and direction of the crack growth. In this section, the relation between the fracture energy and non-local radii is discussed by way of two examples. The properties of the material used in these examples are presented in Table 6-1. In the first example, a notched specimen of a 90º unidirectional laminate (where fibres are parallel to the initial notch) simulated. Similar to the other example with the unidirectional laminates, the dominant damage mechanism is matrix cracking parallel to the direction of the fibres. As shown in Figure 6-10, the predicted damage pattern is in the form of splitting cracks parallel to the direction of the fibres. The fracture energy can be calculated based on its fundamental definition, i.e. the change in the dissipated energy divided by the change in crack length. The advancement of the damage front represents the change in the crack length here. In these simulations r1 was kept constant (equal to 2mm) while different values of r2 were used for each case. Table 6-2 shows the predicted values for fracture energy, FE fG , the visual height of damage, FE visualh , and the so-called effective damage height, FE effectiveh . The effective height of damage is defined as the ratio of the predicted fracture energy and the fracture energy density as shown in Equation (6.5). 79 : FE fFE effective f G h g (6.5) The predicted damage profile at a cross-section on the crack path has a typical bell-shaped distribution as shown in Figure 6-11. The visual damage height, FE visualh , represents the height of the damage zone while effective damage height, FE effectiveh , is defined as the height of the equivalent fully damaged zone with the same fracture energy. The ratio of the effective damage height and the assumed transverse radius is also presented in this table. The effective damage height defined here can be used to calibrate the damage model as discussed later in Chapter 8. The fracture energy, FE fG , can be calculated using Equation (6.6). FE f W U G b a (6.6) where ΔW is the e ternal work done and ΔU is the change in the strain energy of the system when the crack length advances by an amount Δa and b is the width of the specimen. Instead of a discrete crack with a defined length, here we are dealing with damage. The damage length can be defined based on the zone that is fully saturated or one can use the damage front as the crack tip. As shown later in Chapter 8, when process zone is fully formed, both the damage front and damage saturated zones grow at similar rates. Therefore whether we use the saturated zone or damage front to calculate the change in the crack length does not make any significant difference in calculation of the fracture energy. The finite element simulation of a single element under uniaxial strain has been used to estimate the fracture energy density, gf. It should be noted that the material experiences a bi- axial strain condition in front of the notch that is different from the uniaxial strain condition assumed for estimation of gf. As discussed later in Chapter 8, the estimated value of gf can only be used as a starting point for model calibration. According to Table 6-2, the height of damage and also the energy dissipated in the process of transverse matrix cracking mainly depend on the transverse non-local radius, r2. Reducing the value of r2 in the numerical simulation to the value of 0.5mm has caused the crack to grow 80 unstably. By reducing r2, the fracture energy, Gf, is reduced which at some point leads to an equilibrium path that exhibits snap-back that causes unstable crack growth in a displacement– controlled simulation. See (Carpinteri and Colombo 1989) for further information on snap- back phenomenon. In the second example, the notched specimen has a [0/45/90/-45]s quasi-isotropic lay-up. The results of the numerical simulations carried out for models with different longitudinal and transverse averaging radii are shown in Table 6-3. Similar to the table presented for the unidirectional laminate, Table 6-3 includes the predicted fracture energy, the visual height of damage, the effective height of damage and the effective height of damage to r1 ratio. It can be seen that in the case of quasi-isotropic lay-up, the predicted damage characteristics such as height of damage and the fracture energy are governed mainly by the longitudinal averaging radius, r1, and the change in the transverse averaging radius, r2, does not influence the predictions significantly. As an example, changing r2 from 2mm to 1mm results in only a 5.6% change in the predicted fracture energy. Whereas the same change in r1 results in roughly 50% reduction in the fracture energy. Figure 6-12 shows the predicted Force-Displacement (pin opening displacement, POD) graphs extracted from the same simulations. It can be seen that the overall response of the structure is only slightly affected by changing the transverse averaging radius, r2. Based on these graphs it can be seen that reduction of the transverse radius helps eliminate the delayed initiation of damage and thus the sudden jump in the load seen in the numerical predictions using isotropic averaging. This figure also includes the Force-POD curve for the case where r1 is reduced to 1.0mm. It can be seen that unlike r2, changing r1 significantly affects the global behaviour of the system and its energy dissipation. 6.7 SUMMARY An orthotropic non-local averaging scheme was introduced in this chapter. It was shown that in orthotropic materials, the isotropic averaging is unable to accurately predict the direction of damage growth and its pattern. In the orthotropic averaging, two length scales are introduced and the averaging area shapes an ellipse in 2-D. 81 In this chapter, it was shown that the orthotropic averaging scheme significantly improves the predicted crack direction and leads to a more realistic prediction of damage pattern. 82 Table 6-1. Input parameters used in the simulations of the notched specimen presented in this chapter. Damage mode Parameters Matrix εim = 0.8%, εsm = 8.0% Fibre εif = 1.1%, εsf = 8.0% Overall (Anti-locking) εiL = 8.0%, εsL = 10.0% Table 6-2. The predicted fracture energy from the simulations of a notched specimen made from a unidirectional 90º laminate with different transverse non-local radii. The effective damage heights and their normalization with respect to r2 are presented. Radii Predicted FE fG (MPa) Visual damage height FE visualh (mm) Effective damage height FE effectiveh (mm) 2 FE visualh r 2r hFEeffective r1 =2.0mm, r2 =2.0mm 4.1 4 3.3 2 1.65 r1 =2.0mm, r2 =1.0mm 2.0 2.5 1.6 2.5 1.6 r1 =2.0mm, r2 =0.5mm Unstable 2 N/A 4 N/A Local 0.8 0.5 0.65 N/A N/A 83 Table 6-3. The predicted fracture energy from the simulations of a notched specimen made of a quasi-isotropic laminate with different non-local radii. The effective damage heights and their normalization with respect to r1 are presented. Radii Predicted FE fG (MPa) Visual damage height FE visualh (mm) Effective damage height FE effectiveh (mm) 1 FE visualh r 1r hFEeffective r1 =2.0mm, r2 =2.0mm 80.4 5 3.2 2.5 1.6 r1 =2.0mm, r2 =1.0mm 75.9 5 3.0 2.5 1.5 r1 =2.0mm, r2 =0.5mm 74.0 5 2.96 2.5 1.48 r1 =1.0mm, r2 =1.0mm 41.3 3 1.65 3.0 1.65 84 (a) (b) Figure 6-1. Schematic showing two different situations for a material point A where the line connecting the existing damage zones (D1 and D2) to the point A is (a) parallel to the fibres and (b) perpendicular to the fibres. A A D1 D2 1 2 85 Figure 6-2. Schematic showing the influence area of (a) a classical isotropic averaging scheme and (b) the proposed orthotropic non-local averaging scheme. The parallel lines show the direction of the fibres. The shaded zones (D1 and D2) represent existing damage zones in the neighbourhood of point A. A x1 x2 r1 r2 A r1 D1 D2 D1 D2 86 Figure 6-3. Defining d1 and d2 parameters as the components of the AB vector in the fibre’ direction and transverse direction respectively. A B d2 d1 x1 x2 87 Figure 6-4. The orthotropic weight function plotted in the local coordinate system of the ply (x1 and x2). In this example, r1=2 and r2=1. w ei g h t x2 x1 88 Figure 6-5. Predicted crack growth patterns obtained from the simulation of a notched specimen for a 0º unidirectional laminate using various values of the non-local radii. The dashed lines indicate the fibre orientations. Figure 6-6. Predicted damage growth patterns in a 60º unidirectional notched specimen (a) using a local damage model (crack- band), (b) an isotropic averaging with (r1=r2=2mm) and (c) an orthotropic averaging with (r1=2mm and r2=0.5mm) . r1 = 2.0mm r2 = 2.0mm (a) (b) (c) r1 = 2.0mm r2 = 1.0mm r1 = 2.0mm r2 = 0.5mm (a) local (b) r1=2mm, r2=2mm (c) r1=2mm, r2=0.5mm 89 Figure 6-7. Schematic showing the averaging zone around a typical point in (a) an isotropic, and (b) an orthotropic averaging scheme. (a) (b) 90 Figure 6-8. Predicted height of damage in (a) a unidirectional laminate and (b) a multi-directional laminate. h ≈ 2 .5 r 1 h≈ 2. 5 r 2 r 1 r 2 (a) (b) 91 Figure 6-9. Predicted damage patterns in a notched specimen for a quasi-isotropic laminate with (a) an isotropic weight function with a radius of 2mm, (b) an orthotropic weight function with r1=2mm and r2=1.0mm, and (c) an orthotropic weight function with r1=2mm and r2=0.5mm. (a) r1=2mm, r2=2mm (b) r1=2mm, r2=1mm (c) r1=2mm, r2=0.5mm 92 Figure 6-10. Splitting crack forms parallel to the fibres and along the initial notch in the unidirectional 90º laminate. Direction of Fibres Splitting Crack 93 Figure 6-11. A schematic showing the distinction between the visual and effective damage heights in a section (C-C) through the predicted damage profile. 0 0.2 0.4 0.6 0.8 1 -2 -1 0 1 2 (Distance from damage centre)/r D a m a g e P a ra m e te r FE effectiveh FE visualh C C y/r D a m a g e P a ra m et er ( ω ) y 94 Figure 6-12. Predicted Load vs POD curves for an OCT simulation of a quasi-isotropic laminate. 0 4000 8000 12000 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 POD (mm) F o rc e (N ) r1=2.0-r2=2.0 r1=2.0-r2=1.0 r1=2.0-r2=0.5 r1=1.0-r2=1.0 1 22.0, 2.0r r 1 22.0, 1.0r r 1 22.0, 0.5r r 1 21.0, 1.0r r 95 Figure 6-13. An schematic demonstrating the scaling of the averaging weight function near the boundaries of the structure. 0 0.2 0.4 0.6 0.8 1 1.2 1.4 -1.5 -1 -0.5 0 0.5 1 1.5 w Scaled w x r w W w w W Scaled Weight Function Original Weight Function 96 7. VERIFICATION 7.1 INTRODUCTION Verification and Validation (V&V) are procedures that make sure a numerical code is implemented properly (verification) and represents the physics of the problem that it intends to solve (validation). Verification deals with implementation of the mathematical model whereas Validation makes sure the model and its solution are properly representing the physics of the problem considering its intended use (Figure 7-1). V&V are necessary components of any numerical modelling procedure that help to provide confidence in the predictions of the code. This chapter focuses on the verification of the implemented code and validation is discussed in Section 9.2. The recommendations of ASME’s verification and validation guideline (ASME V&V-10 2006) are used here to perform verification and validation procedures. Verification is to determine if the computational model accurately represents the underlying mathematical model and is properly solving the numerical problem. To do so, predictions of the numerical model has to be compared with closed-form solutions or other reliable numerical solutions of the mathematical model available. Since analytical closed-form solutions are only available for problems with simple geometries and loading conditions this comparison may not be enough to verify all aspects of the implemented code. Debugging the code, going through the code step-by-step and reporting intermediate variables are the complementary ways to check the implementation. CODAM2 is implemented in OOFEM (Patzak and Bittnar 2001) as a material model. OOFEM is an object-oriented code and material models are implemented in the form of material objects with some pre-defined functionalities. The main functions of the material model from the point of view of structural analysis are calculations of the stiffness matrix and stresses at the integration points given the strain tensor. The source code for OOFEM is available and the user has to compile the code on the local machine. There are a couple of numerical benchmarks included in OOFEM that can be used to make sure that the core code is compiled properly and is working as expected. The benchmarks are designed to check various parts of the code such as element formulation, calculation of the stiffness matrices and load vectors and finally solving the problem. 97 After running the benchmarks and gaining confidence about the overall performance of the core code, one can verify the performance of the specific implementation. This includes the following steps: a) Debugging and tracking the code and checking that the proper object members are called. This is extremely important in object-oriented programming where inheritance and over-riding concepts exist. Due to probable mistakes in the coding, a function call may be re-directed to the definition of the same function in the parent class. This would not cause any runtime error and can only be identified by tracking the code. b) Comparing the predictions of the code with analytical solutions for simple models. Followings are the examples of test models used to verify the implementation of CODAM2: o Comparison between predictions of the implemented code and a benchmark software. A comparison between the stress fields predicted by two codes would be helpful to verify the accuracy of the implemented material model in the elastic regime. These analyses would help to verify the calculation of the stiffness matrix, load vector and overall compatibility of the material model when used in a structure. o Performing the analysis on a single square element subjected to a uniform and controlled displacement (strain) field and comparing the predictions to the analytical solution. This would verify the elastic response as well as the local damage behaviour of the material model. o To verify the non-local averaging scheme, a square specimen with dimensions much greater than the averaging radii that is subjected to a uniform displacement field is used. In this situation due to the uniform distribution of strain, the averaged equivalent strain has to be equal to the local strain. Comparing the averaged strain with the local strain would verify that the averaging integration scheme is performed properly. 98 7.2 SIMULATING SINGLE ELEMENTS The behaviour of the local formulation of CODAM2 can be verified by simulating a single element undergoing controlled displacement. The elastic properties and damage parameters used in this simulation are presented in Table 7-1 and Table 7-2 respectively. The single element made of a unidirectional laminate is once subjected to uniaxial displacement along the direction of the fibre and once subjected to a combined transverse and shear strains as shown in Figure 7-2. Figure 7-3 compares the predictions of the numerical code with the input data provided. It can be seen that the predictions completely agree with the input data given. This shows that the reading of the input data and the procedures of calculating equivalent strain, damage parameters and reduction coefficients are functioning as expected in the code. 7.3 VERIFICATION OF THE AVERAGING PROCEDURE Writing a new class for orthotropic averaging is part of the implementation of CODAM2 in OOFEM. To verify the performance of the averaging procedure, a square specimen subjected to a uniform displacement field has been simulated (Figure 7-4). If the averaging procedure works properly, due to uniformity of the strain field, the distribution of the non-local equivalent strains is expected to be very close to the distribution of the local strain field. As shown in Figure 7-5, the local and averaged equivalent strain distributions are very close. In the example shown here the averaging radii are chosen to be r1=2.0mm and r2=1.0mm. 7.4 COMPARISON BETWEEN OOFEM AND ABAQUS PREDICTIONS In this example, the prediction of the implemented code is compared to predictions of ABAQUS. This numerical example intends to verify the performance of the implemented model in elastic regimes. It also verifies the interaction between the implemented part and the core of the code. The open-hole geometry subjected to a far-field uniform displacement is used in this verification example. The numerical simulation is performed on a typical open-hole specimen as shown schematically in Figure 7-6. Although the meshes used in the OOFEM and ABAQUS simulations are different, the mesh densities near the hole were kept the same. 99 Figure 7-7 shows the predicted stress profiles at the centre of the specimen as obtained from OOFEM and ABAQUS. It can be seen that the predictions of the implemented model is in good agreement with the predictions of the benchmark code. This example verifies the calculation of the stiffness matrix as well as stresses. It also shows that the communication between the implemented material model and the core code is functioning properly. 7.5 FRAME INDEPENDENCY In this section, an example will be shown to demonstrate the frame independency of the proposed damage model. As it was discussed in Section 5.3, the original formulation of CODAM suffers from an unrealistic dependency on the choice of the co-ordinate system. This example is designed to show that the predictions of CODAM2 are independent of the assigned material coordinate system. Here an OCT specimen made from a [0/45/-45/90]s quasi-isotropic material is simulated. Since this laminate exhibits an isotropic behaviour under in-plane loading in the elastic regime, the choice of the material coordinate system is arbitrary. In this example the simulation has been performed twice as shown in Figure 7-8: once by choosing the x1 direction to be parallel to the initial notch and then by choosing the x1 direction to make a 45º angle to the initial notch. Figure 7-9 shows the predicted load versus the pin opening displacement (POD) curves for these two simulations. It can be seen that the two predictions are identical to each other as expected. It should be noted that in this example, the material was kept in the same orientation and the local coordinate system attached to the material was rotated. 7.6 SUMMARY In this chapter a few numerical examples were presented to verify the performance of the implemented code in OOFEM. It was shown that the predictions of the implemented code in the elastic regime, both in the simple case of one element and averaging for a specimen under a uniformly distributed strain state are as expected. The next two chapters discuss the characterization procedure and validation of the model. 100 Table 7-1. Elastic properties of the material used in verification of a single element. E11 E22 G12 υ12 132 (GPa) 8.5 (GPa) 6.8 (GPa) 0.3 Table 7-2. Damage properties of the material used in verification of a single element. m i m s f i f s L i L s D i D s 0.007 0.2 0.0105 0.2 Deactivated Deactivated 101 Figure 7-1. A schematic showing various stages in numerical modelling of a physical system and roles of verification and validation. Physical system (e.g. structure) Mathematical Model Assumptions and Simplifications Analysis of the Model Interpretation of Results V e ri fi c a ti o n V a li d a ti o n V e ri fi c a ti o n V a li d a ti o n 102 (a) (b) Figure 7-2. A single element used for verification of CODAM2, (a) under uniaxial strain along the direction of the fibre, and (b) under proportionally increasing combined transverse and shear strains. 1 1 2 2 s s 2 s 1 1 103 Figure 7-3. Predicted stress-strain curves compared to the input data for cases (a) and (b) in Figure 7-2. 0 20 40 60 80 100 0.00 0.05 0.10 0.15 0.20 0.25 FE Input 0 5 10 15 20 0.00 0.05 0.10 0.15 0.20 0.25 FE Input 0 20 40 60 80 100 0.00 0.05 0.10 0.15 0.20 0.25 FE Input 0 200 400 600 800 1000 0.00 0.05 0.10 0.15 0.20 0.25 FE Input 11( )MPa 22 ( )MPa 12 ( )MPa 22 ( )MPa 11 22 12 11 22 12 (a) (b) 104 Figure 7-4. A square specimen subjected to uniform uniaxial displacement used in verification of non-local averaging. The dashed line shows the position of the cross-section at which the local and averaged equivalent strains are evaluated and shown in Figure 7-5. L L 105 Figure 7-5. Distributions of the predicted local equivalent strain compared to the averaged equivalent strain. 0 0.05 0.1 0.15 0.2 0.25 0.3 -0.5 0 0.5 S tr ai n ( % ) x/L Series1 Series2 0.1905 0.1907 0.1909 0.1911 0.1913 0.1915 0.1917 0.1919 -0.5 0 0.5 S tr ai n ( % ) x/L eq eq 106 Figure 7-6. Dimensions of the open-hole specimen simulated for verification of the predicted stress field in the elastic regime. 2ρ W=10ρ H = 4 0 ρ σ∞ x y σ∞ 107 Figure 7-7. OOFEM’s predicted stress distribution over a horizontal cross-section passing through the centre of the open-hole specimen compared with the corresponding predictions using ABAQUS. 0 1 2 3 4 -6 -4 -2 0 2 4 6 OOFEM ABAQUS x yy 108 (a) (b) Figure 7-8. Schematic showing the choice of the material coordinate system, (a) x1 axis is chosen to be parallel to the notch, and (b) x1 axis makes a 45º angle with the notch. 1x 2x 1x 2x 109 Figure 7-9. Predicted load-POD curves for the OCT specimens with two different choices for the material coordinate systems. 0 2 4 6 8 10 12 0.0 0.5 1.0 1.5 2.0 2.5 POD (mm) F o rc e (k N ) a b a b 110 8. MODEL CALIBRATION 8.1 INTRODUCTION CODAM2 is a progressive damage model that is designed to simulate the initiation and propagation of damage in composite laminates. As an input, CODAM2 requires the state of strain at which each damage mechanism initiates. It also requires the saturation strain for each damage mode. The non-local radii are the other parameters that need to be provided as input to the code. CODAM2 takes the sub-laminate-based approach where the sub-laminate is the building block of the model. This model aims to capture the effective behaviour of the sub-laminate as a whole. Therefore the proper way of calibrating CODAM2 is using the data from the tests performed on the same laminate that is being modeled. This chapter offers methods to characterize CODAM2 and calibrate its required input parameters. The proposed characterization procedure is mainly based on the results of tensile tests on notched specimens made of the laminated composite being modelled. 8.2 OCT TEST The Overheight Compact Tension test (OCT) designed by Kongshavn and Poursartip (Kongshavn and Poursartip 1999) is widely being used at UBC’s Composites Group to characterize the damage/fracture behaviour of laminated composites. The geometry of the OCT specimen and its typical dimensions are shown in Figure 8-1. In the design of the OCT specimen the height of the standard Compact Tension (CT) specimen was increased to prevent the interference of the boundary with the damage zone and to provide the condition for a stable and progressive crack growth. The line analysis technique was the traditional method used to measure the fracture energy from the OCT specimen. Recently, Zobeiry (Zobeiry 2010) developed an image-analysis- based technique to measure and extract the material parameters from the OCT test. This technique has been used in the characterization example presented in this chapter. 111 8.3 EVALUATION OF THE NON-LOCAL RADII The non-local radii (r1 and r2) play important roles in the proposed non-local damage model. The global response of the structure undergoing damage as well as the local damage characteristics such as the height of damage directly depend on the choice of the non-local radii. 8.3.1 Longitudinal Averaging Radius, r1 Numerical predictions presented in Section 6.6.3 show that in the multi-directional laminates (such as the quasi-isotropic laminate), the behaviour of the structure in terms of energy absorption and predicted height of damage mainly depend on the longitudinal radius, r1. In such laminates, fibres are the main source of energy absorption and their non-local damage behaviour is mainly governed by r1. Therefore r1 can be determined based on the observed height of damage in notched specimens of multi-directional laminates. Based on the data listed earlier in Table 6-3, the predicted height of damage in a quasi-isotropic laminate is approximately 2.5 times r1. 8.3.2 Transverse Averaging Radius, r2 The transverse radius, r2, on the other hand plays the main role in the behaviour of unidirectional laminates. By changing this parameter one can improve the crack pattern predictions as well as the structural behaviour of unidirectional laminates. In such laminates, matrix damage forms in a localized manner (crack). Therefore the transverse radius, r2 cannot be determined based on the observed damage height. In assuming a value for r2 the following points should be considered. - To make the averaging scheme effective, at least a few elements should be included within the averaging range. Therefore r2 cannot be smaller than the element size. - Since the fracture energy associated with the matrix cracking is relatively small, r2 should be chosen such that it inhibits the local snap-back. This is discussed in Section 8.4.2.2. - Considering the discussion presented in Chapter 6, the r1/r2 ratio should be high enough to ensure the proper crack propagation direction in unidirectional laminates. 112 8.4 CALIBRATION OF INITIATION AND SATURATION STRAINS CODAM2 requires the parameters that define the state of damage as a function of the state of equivalent strain for each damage mode. In the current implementation, the damage versus equivalent strain is defined by a hyperbolic function as shown in Figure 5-5. 8.4.1 Initiation Strains The method proposed by Zoberity (Zobeiry 2010) that uses image processing technique to detect the damage accompanied by sectioning and/or de-plying the laminate can be used to determine the initiation strains of matrix and fibre damage. In this method, photos are taken from the surfaces of the specimen and the image processing technique is used to calculate the surface displacement field and derive the strain field at the various stages of the test. The strain states at which damage initiates in the material can be determined by analyzing the measured displacement field. See (Zobeiry 2010) for more details. Alternatively, the failure strains reported in the material datasheets for the longitudinal and transverse directions can be employed here as initiation strains. These values are often extracted from tensile tests on un-notched unidirectional coupons and structural failure of such coupons would occur right after initiation of damage. 8.4.2 Saturation Strains Saturation strains on the other hand are parameters that are not easy to define. At the stage when the material has lost its load carrying capacity, crack openings are large and the definition of strain becomes questionable. Even in an average and smeared sense, saturation strains are subjective parameters that depend on the judgment and interpretation of the analyst. So there is no definitive way of measuring saturation strains. The fracture energy is a parameter that plays an important role in the behaviour of a material in a structure undergoing damage. This parameter can be used to estimate saturation strains for CODAM2. The goal is to set the saturation strains to achieve a resulting fracture energy that is within the desired range. 113 8.4.2.1 Fibre Saturation Strain The behaviour of the fibres and the way they crack depends on the lay-up and configuration of the layers in a laminate. Tests on notched specimens made of unidirectional 0º layers would not show any fibre fracture. Rather, splitting cracks parallel to the fibres are seen in such tests. A centre-notched test on a 0º unidirectional laminate would result in an unstable growth of fibre breakage in the material. In a multi-directional laminate on the other hand, the neighbouring layers provide support and alternative load-carrying paths. A well-designed notched specimen (such as OCT specimen) made of multi-directional laminate would exhibit progressive stable crack growth. Therefore the data and observations from tests on unidirectional laminates cannot be used to calibrate the behaviour of a layer in a multi-directional laminate. CODAM2 aims to simulate the effective behaviour of a laminate therefore tests performed on the same laminate needs to be used to calibrate the model. Fracture energy is a parameter that implicitly includes contributions from all the damage mechanisms activated. As an example, the fracture energy extracted from an OCT test on a quasi-isotropic laminate includes the energy absorbed in damage of the matrix and fibres in various layers. The idea is to estimate values for the saturation strains such that the resulting energy release rate agrees with the experimental measurements. Since in the case of multi- directional composite laminates, there are many damage modes involved, there is not a direct way to estimate the saturation strains based on the fracture energy. An iterative procedure as shown schematically in Figure 8-3 can be used to estimate the fibre saturation strain that includes following steps. In a multi-directional laminate, the saturation strain for the matrix can be assumed to be equal to the fibre saturation strain. 1- Use a trial value for the fibre saturation strain. 2- Simulate a single element subjected to uniaxial strain and calculate the resulting area under stress-strain curve (which represents the fracture energy density, * fg ). 3- As shown in Chapter 6, in the orthotropic non-local scheme proposed here the energy- wise effective height of damage for the fibre breakage is about 1.6 times the 114 longitudinal averaging radius, r1. Therefore the resulting fracture energy can be estimated as * * 11.6f fG r g (see Table 6-3). 4- Compare the estimated fracture energy, * fG , with the fracture energy measured from the OCT test, OCT fG . If the difference is not acceptable go to step 1. This iterative procedure is continued until the saturation strain results in an acceptable value of fracture energy. This procedure results in a rough estimation of the fracture energy and saturation strain since the condition of the elements in front of the crack tip is not necessarily uniaxial strain. After estimating the fibre saturation strain using the procedure mentioned above, the material model should be used to simulate the same experiment that the material characterization is based on (OCT test here). Various predictions obtained from the numerical simulation, such as the force-displacement curve, fracture energy, damage pattern and crack length versus POD can be compared to the experimental measurements in order to finalize the value of the saturation strain. 8.4.2.2 Matrix Saturation Strain In the multi-directional laminates, matrix makes a minor contribution to the dissipated energy; therefore for simplicity one can set the value of the saturation strain for the matrix equal to the fibre saturation strain. In unidirectional laminates the matrix saturation strain can be estimated based on the fracture energy of due to transverse cracking. This can be found in the datasheets provided by the material manufacturer. Considering the r2 parameter determined in Section 8.3 and the discussion in Chapter 6, one can use the following equation to determine the saturation strain of the matrix. 2 2 2 1.6 m fm s m i G E r (8.1) The coefficient of 1.6 in the denominator is determined based on Table 6-2. As discussed in Chapter 8.3, the value of r2 should not be too high to cause local snap-back. As it can be seen 115 in the above equation and illustrated in Figure 8-2, increasing r2 reduces the value of the saturation strain. To avoid the local snap-back and consequent numerical complications, r2 should be chosen such that: 2 2 2 2 1.6 m f m i G r E (8.2) The saturation strain chosen here should be employed to simulate an OCT specimen made of a unidirectional 90º laminate. The predicted fracture energy, FE fG , from this simulation should then be compared to the prescribed fracture energy, m fG . A good agreement shows that the chosen parameters are suitable. Otherwise the saturation strain has to be modified until a good agreement is achieved and FE m f fG G . 8.4.3 Anti-Stress Locking Mechanism Parameters The laminate-level degradation that is implemented to avoid the stress-locking mechanism is activated when the material is experiencing relatively large strains. To avoid the interference of the anti-locking degradation with fibre damage, the initiation strain for this degradation mechanism should be chosen to be greater than the fibre saturation strain. The saturation strain for this degradation mechanism should be chosen to be large enough to allow for a gradual loss of stiffness in the laminate. When the anti-locking degradation is saturated, the laminate would completely lose its load carrying capacity in all directions. 8.4.4 Delamination Parameters As discussed before, the sub-laminate approach in its current form is unable to predict the onset of delamination and to simulate the effect of delamination in its general form. For the special case of edge delamination, the approach introduced y O’Brien (O'Brien 1981) is implemented here to simulate the effect of edge delamination on effective stiffness of the laminate. To calculate the strain state at which the edge delamination initiates based on Equation (5.25), the delamination fracture energy, Gc, and the number and location of the delamination cracks are required. The saturation strain for this damage mode is assumed to be slightly greater than 116 its initiation strain. See Section 5.5.6 for further details on the way edge delamination is treated in CODAM2. 8.5 MODEL CALIBRATION: AN EXAMPLE In this section, the calibration procedure for CODAM2 with anisotropic averaging will be demonstrated. The calibration is based on the material characterization tests performed by Zobeiry in his doctoral thesis (Zobeiry 2010). The OCT test specimens employed for material characterization had a quasi-isotropic lay-up [90/45/0/-45]4s made of IM7-8552 CFRP prepregs. The elastic properties of the IM7/8552 unidirectional lamina are listed in Table 8-2. 8.5.1 Non-Local Radii The longitudinal averaging radius is estimated based on the visual height of damage. In the OCT tests, the observed height of damage was about 5-6mm. Based on the discussion presented in Section 8.3.1, to achieve a visual height of damage of 5-6mm in the numerical simulation, r1, should be chosen to be about 2mm. The transverse non-local radius, r2, does not have a significant effect on the visual height of damage nor does it affect the fracture energy of the laminate. Since the mesh size in the numerical simulations presented in this chapter is about 0.5mm and greater, r2 is chosen to be 1mm so that the non-local averaging can be efficient (i.e. cover more than one element). 8.5.2 Initiation Strains The value of the initiation strain for fibre damage is estimated based on the value given in (Zobeiry 2010) as the damage initiation for the laminate of about 1.1%. Damage in the matrix initiates at lower strain levels but since the change in the laminate stiffness due to matrix cracking is very minor, the image processing technique used in (Zobeiry 2010) is unable to capture it. A value between 0.6% and 0.8% was estimated for the matrix initiation for this material in (Forghani et al. 2011). The value of 0.8% is reported as the transverse failure strain for the unidirectional laminate. The value of 0.6% was estimated based on a simple analysis on the residual process-induced tensile stresses in the matrix in this laminate that result in initiation of matrix cracking in a strain level lower than 0.8% (See (Forghani et al. 2011) for further details). 117 8.5.3 Saturation Strains The saturation strain for fibre and matrix damage can be estimated based on matching the fracture energy of the laminate with the measured Gf. As shown in Chapter 6, the predicted fracture energy is about 1.6 times the fracture energy density times the longitudinal non-local radius. Based on the simulation of a single element under uniaxial strain, the fracture energy density can be calculated and an estimation of the expected Gf can be made on that basis. Table 8-4 shows the values of expected fracture energy estimated for different values of saturation strain. The final step in the calibration of the model is to simulate the OCT specimen and choose the saturation parameter that leads to the closest prediction of the Load-POD curve when compared to the experimental measurements. All the values of saturated strain listed in Table 8-4 resulted in estimated fracture energies that are within the range observed in the experiments. Therefore the simulation of the OCT specimen is performed three times to try all these values. The values of the damage model parameters employed in the calibration simulation are listed in Table 8-5. Figure 8-4 shows the predicted load-POD curves of the OCT simulations using various values of saturation strains listed in Table 8-4. The predictions are compared to the two sets of measured data from the experiments. The curves that correspond to saturation strains of 20% and 22% show the best agreement with the experimental observations. The one with 22% saturation strain was chosen as the best fit and the rest of results presented here are based on this value. Figure 8-5 shows the predicted Load versus POD for the model using a value of 22% for the saturation strain. It can be seen that the predicted response agrees very well with the measured test data in terms of peak force and post-peak softening behaviour. Figure 8-6 shows the predicted length of the damage zone compared to the measured damage length. In this figure, the lengths of predicted damage front and fully saturated damage zone are shown. It can be seen that the experimentally observed damage length lies within the predicted process zone. 118 Figure 8-7 shows the predicted Load-POD curve overlaid on the damage-tip position graph. It can be seen that up to a POD of about 1.5mm, the process zone in front of the crack is being formed. At this stage, almost all the parameters in the damage model contribute to the response of the structure. After this point, the process zone starts to move forward and the specimen exhibits a self-similar crack growth. The behaviour of the model after this point is mainly determined by the fracture energy, Gf. 8.6 SUMMARY In this chapter, the procedure for calibrating CODAM2 was presented. It was shown how the available data from a sharp-notched specimen along with the information from the elastic properties of the material can be used to calibrate the model parameters. The next step is to validate the model by comparing the numerical predictions with experimental data other than the one used to calibrate the model. This is the subject of next chapter. 119 Table 8-1. Summary of calibration procedure for non-local CODAM2 parameters. Parameter Calibration Procedure Longitudinal averaging radius: r1 Estimated based on the observed height of damage in the notched specimen Transverse averaging radius: r2 Assumed to be a fraction of r1 (mesh size limits the lower bound for r2) Damage Initiation Parameters Damage detection technique or the material datasheets Damage Saturation Estimated based on the Fracture Energy (Gf) of the laminate and the r1 Table 8-2. Elastic Properties of IM7-8552 unidirectional lamina reported by the manufacturer (Hexcel Company). E11 E22 G12 υ12 150 (GPa) 11 (GPa) 4.6 (GPa) 0.3 120 Table 8-3. Material parameters for IM7-8552 quasi-isotropic laminate extracted from OCT tests performed by N. Zobeiry (Zobeiry 2010). Parameter Value Initiation Strain for Fibre Breakage 1.1% Fracture Energy (Gf) 95-115 kJ/m 2 Visual Height of Damage of the Laminate (hc) 5mm Table 8-4. The values of fracture energy density and expected fracture energy calculated for various values of saturation strain Matrix Initiation Saturation Strain Fracture Energy Density, * fg (MPa) Estimated Fracture Energy, * * 1( 1.6 )f fG r g (kJ/m 2 ) 0.6% 18% 29.3 93 0.6% 20% 34.2 110 0.6% 22% 39.1 126 121 Table 8-5. Damage model parameters used in calibration of the model. Parameter Value Longitudinal Non-local Radius (r1) 2.0mm Transverse Non-local Radius(r2) 1.0mm Initiation Strain for Matrix Damage 0.6% Initiation Strain for Fibre Damage 1.1% Saturation Strain for Fibre Breakage and Matrix Cracking 18% - 22% 122 Figure 8-1. Schematic showing the OCT specimen geometry and its typical dimensions. a = 33 mm h = 1 0 4 m m d = 19.1 mm W = 80.6 mm 106 mm 8 4 .6 5 m m c = 38.7 mm Initial Notch 123 Figure 8-2. Increasing the transverse non-local radius, r2, while keeping the other parameters constant results in decreasing saturation strain for the matrix damage. ε2 σ2 124 Figure 8-3. The iterative procedure that is used to calibrate the saturation strains given the height of damage and fracture energy. Guess an initial value for the saturation strain (εs)* Simulate a single element with the material properties of the laminate subjected to uniaxial strain Calculate the fracture energy density, gf fg d εLaminate response σ Estimate the resulting fracture energy, Gf * = 1.6 r1 gf Compare Gf * to the fracture energy measured from the test ( ) If Gf * and are not close enough, adjust the value of the saturation strain accordingly If Gf * and agree, the estimated saturation strain will be used for simulation of the notched specimens OCT fG OCT fG OCT fG 125 Figure 8-4. The predicted Load-POD curves from simulation of OCT specimen for various values of the saturation strains. CTA1 and CTA2 represent the experimentally obtained curves by (Zobeiry 2010). 0 2 4 6 8 10 12 0.0 0.5 1.0 1.5 2.0 L o ad ( k N ) POD (mm) CTA1 CTA2 Sat 18% Sat 20% Sat 22% 126 Figure 8-5. Predicted Load-POD curve for the calibrated model compared to the measured data. 0 2 4 6 8 10 12 0.0 0.5 1.0 1.5 2.0 L o ad ( k N ) POD (mm) CTA1 CTA2 Calibrated Model 127 Figure 8-6. The predicted position of damage initiation and saturation fronts as a function of POD compared to the experimental measurements. 0 5 10 15 20 25 0.5 1.0 1.5 2.0 POD (mm) C ra ck -t ip p o si ti o n ( m m ) CTA1 CTA2 Damage Tip Damage Sat Ave. Damage Length i (Initiation) aturation 128 Figure 8-7. The predicted position of damage initiation and saturation fronts. 0 2 4 6 8 10 12 0.0 0.5 1.0 1.5 2.0 POD (mm) F o rc e (k N ) 0 5 10 15 20 25 C ra ck -t ip p o si ti o n ( m m ) FE, Force FE Damage Tip FE Damage Sat L o ad ( k N ) Load-POD Curve Damage Initiation Damage Saturation 129 9. MODEL VALIDATION AND NUMERICAL EXAMPLES 9.1 INTRODUCTION To demonstrate the performance of CODAM2 in simulating damage initiation and growth in laminated composites a few numerical examples are presented in this chapter. The first example that is also used to validate the model is on simulation of damage behaviour of a series of notched specimen varying from sharp notch to blunt-notch. The second example is on simulation of open-hole tension specimens used as test cases in the third round of the worldwide failure exercise (WWFE-III) (Kaddour et al. 2009). All the numerical simulations presented in this chapter are modelled as dynamic problems with the explicit time integration scheme employed to obtain the solution. The main reason for adopting the explicit dynamic solution is the highly non-linear nature of the structural response and static instabilities that make it impossible for the implicit methods to converge to a satisfactory solution. Modelling quasi-static events using explicit time integration can be very time consuming and computationally expensive. On one hand the rate of application of load/displacement has to be slow enough to represent the quasi-static condition and minimize the effect of inertial forces. On the other hand, the size of time steps is governed by the CFL condition that depends on the finest part of the mesh and the time step can be quite small where fine meshes are used in the discretization. Mass scaling technique that is usually employed to reduce the computation time is not helpful here. Increasing the mass would highlight the inertia effect and requires slower loading rate, and decreasing it would affect the maximum time step size. Mass scaling is usually helpful if it is only applied to part of the structure but not the whole structure. The typical solution time for each of the numerical examples presented in this chapter is about a few days (40-100hrs) on an Intel ® Core2Quad ® 3.2 GHz CPU. Due to the time limitations the loading rates applied were not slow enough to eliminate the inertia effect. The waviness in the predicted load-displacement curves presented in this chapter is a result of the higher frequency stress waves that contribute to the response of the structure. 130 Lack of an appropriate damping in OOFEM also contributes to the waviness of the response. The only damping option available in the version of OOFEM used to implement CODAM2 (OOFEM-1.9) is mass proportional damping. The mass proportional damping is only effective on low-frequency modes and the rigid body motion. This type of damping cannot help to reduce the higher frequency noise in the response. 9.2 VALIDATION Validation is the procedure of ensuring that a model is an accurate representation of the physical reality (from the perspective of intended use of the model). To validate a model, its intended use, the scope of its application and its limitations need to be identified first. The steps to validate a material model are as follows: 1) Verification 2) Calibration 3) Use the already verified and calibrated material model to numerically simulate some test cases 4) Compare the predictions with the test data It should be noted that the geometry and loading condition for the tests that are used for validation have to be different from the test that is used for calibration of the model. For example, if the material model is calibrated using an OCT test, the simulation of the same OCT test cannot be used for model validation purposes. In the following sections, the scope of CODAM2 is discussed and then numerical examples that include both calibration and validation of the model are presented. 9.3 SCOPE, APPLICABILITY AND LIMITATIONS CODAM2, similar to its predecessor CODAM, is a sub-laminate based model where the sub- laminate is assumed to be the representative volume. In this approach it is assumed that the layers in the sub-laminate experience compatible displacement and strain fields in an average sense. When delamination is not widely spread and the damaged zone is surrounded by an undamaged elastic region (as shown in Figure 9-1), at the borders of damage the iso- displacement condition is enforced by undamaged region. Therefore one can assume that all the layers experience similar average strains within the damage zone. 131 It is worth noting that the smearing is applied over the RVE and a sub-laminate-based model is not intended to capture the details of damage events inside the damage zone. The common practice in the design of composite structures is to avoid delamination. The practical laminates used in these structures are made by stacking thin plies and designers try to minimize the inter-laminar stresses by optimizing the stacking sequence. The free edges are also avoided or treated in the design of composite structures. Therefore, it is unlikely for those laminates to delaminate under in-plane loading conditions. Lateral low velocity impact may cause delamination in such laminates. As long as delaminations formed by the lateral impact are limited to small regions and are surrounded by undamaged material, it would not affect the intra-laminar behaviour of the laminate significantly. When delamination is spread over a large region (comparing to the size of the structure) and delamination is a major player in the response of the structure, the iso-strain assumption would not be valid anymore. In such situations, the notion that the original sub-laminate is still the representative volume would not be accurate. Independent displacement fields need to be assigned to the layers that are separated from each other by delamination. The method employed here to determine the effect of delamination on stiffness of the laminate is applicable only to edge delamination. Test observations show that the laminated composites subjected to in-plane compressive load exhibit residual deformation after complete removal of the load. This is mainly due to friction and interlocking of layers. The current version of CODAM2 is an elastic damage model and is not capable of simulating residual deformation and frictional/plastic behaviours. 132 9.4 SIMULATION OF THE MODIFIED OCT TESTS WITH VARIOUS NOTCH RADII 9.4.1 Test Specifications In this example used for validation purposes, the non-local damage model, CODAM2 is employed to simulate the damage behaviour of modified OCT specimens made of CFRP laminates. Tests that are modelled here were designed and conducted by Scott McClennan, a former member of UBC Composites Group. McClennan (McClennan 2004) modified the original design of OCT by blunting the notch (Figure 9-2(b)) to observe the damage behaviour of material when the stress concentration at the notch tip changes. A wide range of specimens with notch root radii of 2mm to 27mm were tested (Figure 9-3). The line analysis technique (Kongshavn and Poursartip 1999) was employed in those tests to track the damage tip. Some of the specimens were also sectioned afterwards to measure the damage height. Blunting the notch eliminates the high stress concentration at the notch tip leading to a stress field with milder gradients. By decreasing the concentration of stress at the notch tip, blunting the notch delays the initiation of damage and eventually leads to an abrupt and unstable crack growth in the specimens with larger notch root radii. In this example the damage model is calibrated using the data from the sharp-notched OCT specimen and then the calibrated model is employed to simulate the specimens with various blunt notches. 9.4.2 Model Parameters The specimens were cut from laminates with [+45/-45/902/0/902/-45/45]6 lay-up manufactured using resin infusion technique by the Boeing Company. The fibres were a mixture of standard and intermediate modulus Hercules AS fibres and the Hercules 3501-6 epoxy resin was used as the matrix. The elastic properties of the laminate are reported in Table 9-1. The parameters used in the damage model such as the non-local radii and initiation and saturation strains were calibrated based on the data from the test conducted on the sharp-notched specimen. The height of damage was measured by sectioning the specimens and reported in (McClennan 2004). The measured damage heights of the sharp notch specimen show a height of about 133 6mm at sections that are far enough from the initial notch. As discussed in Chapter 6, the observed overall height of damage in multi-directional laminates, depend mainly on the longitudinal averaging radius, r1. The numerical simulations showed that a longitudinal radius of 2mm would result in a damage height of about 6mm. The transverse averaging radius, r2, does not have a major contribution to the observed damage height. Based on the discussions in Chapter 6 and considering the mesh size that was about 0.5mm-0.75mm in the area with the potential of damage growth, r2 is chosen to be 1mm that is half of the value chosen for r1. A first estimation for initiation strains for matrix and fibre damages was done based on the failure strain data available for the unidirectional laminate. A failure strain of about 1.05% along the fibres and about 0.6% in the transverse direction is reported for the unidirectional laminate. It should be noted that since the laminate is made by resin infusion through fibres, the available data reported for the unidirectional laminates made from prepreg can only be used as an initial estimation. The first estimation of the saturation strain is performed based on the reported Mode I fracture energy of the laminate. The fracture energy for the laminate as reported in (McClennan 2004) is around 70-85 kJ/m 2 . Based on the discussion in Chapter 6, the fracture energy density of the material model is estimated to be around 22-26MPa. Numerical simulation on a single element under uniaxial strain, suggests that the saturation strain should be about 18%-22%. The numerical simulation of the sharp-notched specimen is then performed using a series of values of initiation and saturation strains around the initial estimates and the predictions in terms of load versus crack mouth opening displacement (CMOD) curves are compared to the experimental data. The parameters that lead to the best agreement with the experimentally measured load-CMOD curve are presented in Table 9-2. 9.4.3 Predictions The calibrated model is employed to predict the response of the blunt-notched specimens. The predictions of the FE model in terms of load-CMOD curves, peak forces and damage patterns are compared to the test measurements. 134 9.4.3.1 Load-CMOD Curves Figure 9-4 shows the predicted load-CMOD curve of the calibrated model for the sharp- notched specimen compared to the test data. Since the model was calibrated with this specific test, there is no surprise that the predictions are in a very good agreement with the measured data. Figure 9-5 and Figure 9-6 show the predicted load-CMOD curves for the blunt-notched specimens with notch root radii of 12.7mm and 16mm. The agreement between the predicted curves and test data in terms of peak load and the post-peak behaviour is very good showing that the proposed non-local damage model is capable of accurately capturing the initiation of damage as well as its propagation. Figure 9-7 compares the predicted load-CMODs of the sharp-notched specimen and the blunt- notched specimen with radius of 16mm. The sharp-notched specimen as expected has a higher initial stiffness. In the sharp-notched specimen, the non-linear behaviour starts earlier and leads to a lower value of the peak force. This is because the stress at the tip of the notch at the same level of external force is lower in the blunt-notched specimen which leads to a later initiation of damage in the material. Unlike the sharp-notched specimen that shows a very gradual and steady damage growth and stiffness reduction, in the specimen with blunt notch, an abrupt drop in the stiffness is seen that is due to a sudden growth of damage in the specimen. The jump in the crack length and, as a result, the drop in the stiffness of the specimen is more significant in the specimens with larger notch root radius. In the specimens with the largest notch root radius, the sudden crack growth leads to instability of the structure under the displacement control condition that is a sign of snap-back in the equilibrium path. This behaviour is in agreement with the test observations. The sharp-notched specimen and the specimens with the smaller radius exhibit stable response with gradual crack growth and reduction in the load carrying whereas the specimens with larger notch radius exhibit a more brittle behaviour with sudden drop in the load and jump in the crack length. A comparison of the predicted peak loads and test data is presented in Figure 9-8. It can be seen that the predictions of the numerical model agree well with the experiments and follow the trend seen in the tests. Starting from the sharp-notched specimen, the peak force increases by increasing of the notch radius. From the root radius of 3.1mm to 16mm, the predicted peak 135 force is almost constant. In specimens with notch root radius greater than 16mm, the predicted peak force starts increasing again. The latter group of specimens are the ones that exhibit an unstable crack growth. Figure 9-9, Figure 9-10 and Figure 9-11 show the predicted position of the damage tip and damage saturation tip for the sharp-notched specimen and blunt-notched specimens with root radii of 16mm and 22.5mm, respectively. These crack position curves are overlaid on load- CMOD curves to demonstrate the correlation between the crack growth and the load- displacement curves. It can be seen from these figures that the sharp-notched specimen exhibits a smooth growth of the damage. By increasing the radius of the notch, the specimens show a more abrupt change in the crack position just after reaching the peak load which coincides with the sharp decrease in the secant stiffness of the specimen. 9.4.3.2 Damage Pattern In addition to accurate prediction of the load-displacement curves, the proposed non-local model is also capable of accurately predicting the damage pattern and height of damage. In these experiments, the height of damage has been measured by sectioning the specimens and looking at the distribution of the cracks over the cross-section (McClennan 2004). The test observations show that the initial height of damage depends on the radius of the notch. The specimens with larger notches show a greater initial height of damage. This is due to the fact that in the specimens with larger notches, the strain distribution in front of the notch is more uniform and therefore a larger region experiences high strains and become damaged. As damage progresses in the material, all the specimens show convergence of the damage zone to a similar height. This is an evidence that the height of damage, when damage is far from the boundaries of the structure, is a characteristic property of the material. The same trend is also seen in the predicted damage patterns. Figure 9-12 shows the predicted pattern of the matrix damage in the 90º layer. It can be seen that the predicted initial damage height depends on the notch radius and specimens with larger notches show a larger predicted damage zone. As we move away from the original notch tip, the predicted damaged area converges to zones with similar heights in different specimens. 136 The predicted patterns of the damage in the blunt-notched specimen with the largest notch radius (ρ=27mm) and the sharp-notched specimen are compared to the measured damage heights in Figure 9-13 .It can be seen that the predicted damage pattern is in a fairly good agreement with the measured data. Figure 9-14 shows the predicted profiles of matrix cracking in the 90º layer at two different cross-sections. In Figure 9-14(a) the cross-section is in the vicinity of the original notch. It can be seen that in the specimens with larger notch radius, near the notch, the damage is more spread and therefore a larger damage height is observed. In all the specimen sizes, although the damage is initiated in areas with the sizes that depend on the notch radius, the height of the zone that reaches the damage saturation is similar for all the specimens. The saturated damage zone represents the area in which macroscopic crack is formed. Figure 9-14 (b) shows the predicted damage profiles at a cross-section that is 20mm away from the original notch tip. It can be seen that at this distance from the notch, all the specimens regardless of their notch root radius have a similar damage profile and similar heights of damage. 9.4.4 Comparison between the Non-local and Local FE predictions In the case of OCT specimen, since the crack direction is known in advance, in a local simulation, the mesh can be designed to make the crack grow in the expected direction. This has been the case in the OCT simulations presented in (McClennan 2004) and (Zobeiry 2010) among others. The height of damage in this case depends on the mesh size. Unlike the non- local model that predicts blunting the crack tip by progression of damage and an increase in the height of damage until reaching the characteristic damage height, the local model localizes the crack into a small height that is representative of a sharp crack. Due to misrepresentation of the crack-tip bluntness, the local models are unable to predict the stress distribution in front of the crack tip. In a local OCT simulation with a mesh size of say 1.0mm the predicted damage pattern would be a relatively sharp crack with height of 1.0mm (Figure 9-15(a)). Whereas in reality the damaged zone is blunt with a height of about 6mm that leads to a very different near tip stress 137 distribution. A non-local model on the other hand is capable of simulating the bluntness at the notch-tip (Figure 9-15(b)) and therefore the non-local model leads to a more realistic representation of the near-tip stress/strain fields. Figure 9-16 (a) and (b) shows the predicted matrix damage pattern of the 90º layer using local and non-local simulations respectively. Figure 9-16 (c) and (d) show the predicted longitudinal stress distribution (σyy) obtained using the local and non-local simulations, respectively. It can be seen that the non-local model predicts a smoother stress distribution. 138 9.5 OPEN HOLE TENSION The open-hole tension simulations presented here are part of the contribution of the UBC Composites Group to the 3 rd round of the Worldwide Failure Exercise (WWFE-III) of composite materials (Forghani et al. 2011). In this exercise, the leading research groups in the field of damage modelling of composites were challenged to perform blind predictions of the behaviour of composite laminates with various geometries under various loading conditions. In this section a summary of the predictions made using CODAM2 for the open-hole tension is presented. The details of the analysis can be found in (Forghani et al. 2011). Geometrical specifications and labelling used to identify the simulated open hole tensile test specimens are presented in Table 9-3. As the heights of the specimens were not specified in the WWFE-III instruction, for each hole size (diameter), two different heights are assumed. Both square specimens (height to width ratio of 1) and tall rectangular specimens with height to width ratio of 4 are considered. The geometries of these specimens are shown schematically in Figure 9-17. The open-hole specimens that are subject of this example are cut form [454/904/− 4/04]s quasi-isotropic laminate made from IM7/8552 prepreg. The lay-up is made by stacking thick layers with thickness of 0.5mm that is not a practical lay-up due to the potential chance of delamination. In this test case, due to the lay-up of the laminate that consists of thick plies, delamination around the hole and also at the edges of specimen become important. The approach presented in Section 5.5.6 of this thesis has been employed here to simulate the effect of delamination on the in-plane stiffness of the laminate. 9.5.1 Material Parameters The input parameters were partly provided by the WWFE-III instruction. The remaining required parameters were estimated based on the available experimental data at UBC or extracted from the literature. With respect to damage behaviour of the material, WWFE-III instruction only provided failure strains of the unidirectional laminates. The initiation strain parameters required in CODAM2 are estimated based on the provided failure strain values given in the WWFE-III instruction. The matrix initiation strain has been 139 modified based on an estimation of the residual thermal stresses in the laminate due to CTE mismatch as presented in Table 9-4. Due to lack of information on the height of damage in this lay-up, the averaging radius is chosen to be 2mm that would lead to a height of damage of about 5-7mm that is the height of damage observer in quasi-isotropic laminates made of the same material but with thinner plies. The mode I fracture energy, Gf, for this laminate was similarly estimated based on the values measured for the same material but in a laminate with thin layers. 9.5.1.1 Delamination The approach proposed y O’Brien (O'Brien 1981), as described in Section 5.5.6 of this thesis, is employed to estimate the stiffness drop due to delamination and to estimate the strain level at which delamination damage occurs. Delamination is likely to initiate either between the 45º and 90º layers or between the 90º and -45º layers. Considering each scenario separately, the in-plane strain at which delamination would initiate, εc, is calculated based on the approach presented in Section 5.5.6 and in Table 9-5. Since delamination generally occurs under a mixed-mode condition, the delamination fracture energy, Gc, is assumed to vary between the two limits of 0.2 and 0.8 kJ/m 2 corresponding to GIc and GIIc values, respectively (see (Kaddour et al. 2009)). These two values of Gc would provide lower and upper bounds for εc. As shown in Table 9-5, the estimated values of the effective modulus of the delaminated laminate, E * , are very close for the two scenarios considered. A simple analysis that considers the equilibrium of forces at the free edge shows that the value of the bending moment (caused by inter-laminar stresses) at the 45º/90º interface is just slightly higher than the one at the 90º/- 45º interface. Due to the very similar stress states existing at these two locations, the delamination would likely initiate at either interface. A value of εc = 0.4% (this is within the range given inTable 9-5) has been chosen here as the in-plane strain at which delamination initiates. 9.5.1.2 R-Curve Behaviour in Delamination Delamination tests performed on laminated composites show an increase in Gc as the delamination crack propagates from the edges and therefore, the material exhibits an R-Curve type behaviour with increasing resistance. In the predictions presented in this paper, the 140 increase of Gc is simulated by dividing the area in the vicinity of the free edge into various regions (strips), as shown in Figure 9-18, and assigning increasing values of critical strain for damage (delamination) initiation. Figure 9-18 also shows the R-Curve assumed for the predictions of this test case. A summary of the input parameters used in the numerical calculations are presented in Table 9-6. Details of the input parameters used in the finite element simulation including initiation and saturation strains for various modes of damage are provided in Table 9-7. 9.5.2 Numerical Predictions 9.5.2.1 Stress-Strain Curve and Damage Pattern The predicted overall (nominal) stress versus strain curves of specimens with four different hole sizes (D = 3.18mm, 6.35mm, 12.7mm and 25.4mm) and proportionately different widths are shown in Figure 9-19. Each graph shows the predictions for both tall and square specimen geometries. The sequence of events is common between the specimens with different sizes. Damage starts with delaminations emanating from the hole’s oundary and free edges of the specimen. With increasing strain level, delamination grows towards the centre of the specimen and is accompanied by matrix cracking within the layers. Further increase in the applied displacement results in initiation of fibre damage. Fibre breakage starts in the 0º layer near the hole’s periphery and continues to grow until the ultimate load is reached. Figure 9-21 shows the various stages of damage progression in the smallest simulated specimen (D3.175S). According to Figure 9-19, in the smaller specimens (the ones with 3.175mm and 6.35mm hole diameters), the damage process is gradual and occurs in a stable manner. In contrast to the smaller specimens, the larger ones show a more brittle response and fail at a lower value of applied strain. In these specimens, the initiation of fibre breakage results in structural instability and ultimate collapse of the laminate. Availability of excessive amount of strain energy in the structure and the inadequate amount of energy absorbed by the fibres during their breakage is the reason for the unstable crack growth and the resulting brittle behaviour in the larger specimens. 141 9.5.2.2 Failure Stress, Delamination Level and Crack Density Table 9-8 presents a summary of the numerical predictions for the open-hole tension specimen. Included in this table are the predicted values of the failure stress, ultimate stress, and maximum matrix crack density as well as the delamination level and location. Failure strength is defined as the stress level at which the secant stiffness of the laminate drops by 20% while ultimate strength corresponds to the maximum (peak) stress that the specimen reaches in the numerical simulation. Figure 9-22 shows an example where the failure and ultimate strengths are labelled on the predicted stress-strain curve for the D3.175S specimen. Figure 9-23 shows the predicted nominal failure and ultimate stresses versus the hole size. Figure 9-24 shows the predicted delamination pattern for the tall specimens at the instant corresponding to the ultimate load. The unsymmetrical predicted delamination pattern in the specimen with 25.4mm hole size (largest specimen) is due to the instability of the structure just after the ultimate load. In unstable conditions, structures with symmetrical geometry and loading condition may exhibit unsymmetrical damage/cracking patterns. 142 Table 9-1. Elastic properties of the laminate used in the modified OCT simulations (McClennan 2004). Parameter Value E1 75 (GPa) E2 32 (GPa) υ12 0.16 G12 17 (GPa) 143 Table 9-2. Damage model parameters obtained based on sharp-notched OCT test data. Parameter Value Matrix damage initiation strain (tensile) i m 0.8% Matrix damage saturation strain (tensile) s m 20% Fibre damage initiation strain (tensile) i f 1.1% Fibre damage saturation strain (tensile) s f 20% Longitudinal radius used for non-local averaging (mm) r1 2mm Transverse radius used for non-local averaging (mm) r2 1mm 144 Table 9-3. Labelling and dimensions of the simulated specimens considered for the open-hole tension specimens (WWFE-III’sTest Case 12). Label Hole Size (mm) Width (mm) Height (mm) Thickness (mm) D3.175S 3.175 15.875 15.875 4 D3.175T 3.175 15.87 63.5 4 D6.35S 6.35 31.75 31.75 4 D6.35T 6.35 31.75 127.0 4 D12.7S 12.7 63.5 63.5 4 D12.7T 12.7 63.5 254.0 4 D25.4S 25.4 127.0 127.0 4 D25.4T 25.4 127.0 508.0 4 145 Table 9-4. Estimated values of processing-induced residual stresses in various layers due to the layer-wise CTE mismatch in quasi-isotropic laminates considered in Test Case 12. Layer Residual Stresses Equivalent Mechanical Strains 45 σ11 = -23.5MPa σ22 = 23.5MPa σ12 = 0 ε11= -0.014% ε22= 0.26% 12= 0 90 σ11 = -23.5MPa σ22 = 23.5MPa σ12 = 0 ε11= -0.014% ε 22= 0.26% 12= 0 -45 σ11 = -23.5MPa σ22 = 23.5MPa σ12 = 0 ε11= -0.014% ε22= 0.26% 12= 0 0 σ11 = -23.5MPa σ22 = 23.5MPa σ12 = 0 ε11= -0.014% ε 22= 0.26% 12= 0 146 Table 9-5. Calculated range of critical strain for onset of delamination in quasi-isotropic laminates considered in Test Case 12. Delamination Interface E * (GPa) εc 45º / 90º 52.35 0.32% - 0.56% 90º / -45º 51.86 0.32% - 0.56% 147 Table 9-6. Values of the input parameters used in numerical predictions of Test Case 12. Parameter Value provided by WWFE-III instruction Value used in the simulation * Comments and Source Matrix damage initiation strain (tensile) i m 0.81% 0.6 % Estimated considering the thermal and cure shrinkage effects. See (Forghani et al. 2011) for details. Fibre damage initiation strain (tensile) i f 1.50% 1.50% Edge delamination critical strain (tensile) c N/A 0.40% See Section 9.5.1.1 for details Delamination energy (kJ/m 2 ) Gc Mode I: 0.2 Mode II: 0.8 0.35 (0.2 - 0.8) See Section 9.5.1.2 for details of the R- Curve behaviour Fracture energy (kJ/m 2 ) Gf N/A 105 (95 - 115) Based on UBC experiments (Zobeiry 2010) Characteristic damage height (mm) hc N/A 5.0 -7.0 Based on UBC experiments (Zobeiry 2010) Radius used for non-local averaging (mm) r N/A 2.0 * Values in parentheses indicate the possible range while the value outside the parentheses is the average value used. 148 Table 9-7. Details of material parameters used in the numerical simulation of Test Case 12. Damage Mechanism () Zone Initiation and Saturation Strains (%) Stiffness reduction rate (β) Delamination (D) Zone 1 i D = 0.40 ; s D = 0.42 0.2 Zone 2 i D = 0.45 ; s D = 0.47 Zone 3 i D = 0.50 ; s D = 0.52 Zone 4 i D = 0.55 ; s D = 0.57 Zone 5 i D = 0.60 ; s D = 0.62 Fibre Damage (f ) All zones i f = 1.5 ; s f = 5.0* 1.0 Matrix Damage (m) All zones i m = 0.6 ; s m = 5.0* 1.0 Stress Locking (L) All zones i L = 5.0 ; s L = 7.0 1.0 * Represents the nominal value for damage saturation strain used in the analyses and corresponds to the fracture energy of Gf = 105 kJ/m 2 . The two other Gf values listed in Table 9-6, namely, 95 and 115 kJ/m 2 result in saturation strains of 4.5% and 5.5%, respectively. 149 Table 9-8. Summary of predictions of Test Case 12. Hole Size (mm) Predicted Failure Stress (MPa) Predicted Ultimate Stress (MPa)* Maximum Crack Density Delamination level / location 3.175 250 381 (378 - 390) 1/0.5mm The whole specimen is delaminated 6.35 295 353 (338 - 361) 1/0.5mm Near the edges and a triangular shape starting from the hole 12.7 305 308 (303 - 312) 1/0.5mm Near the edges and a triangular shape starting from the hole 25.4 276 277 (273 - 281) 1/0.5mm A triangular shape starting from the hole * The values shown here correspond to the nominal saturation strains for fibre and matrix damage, namely 5.0%, while the values in parentheses represent the predictions obtained for the two bounds on the saturation strain, namely, 4.5% and 5.5% corresponding to the two limits of fracture energy, 95 and 115 kJ/m 2 listed in Table 9-6. 150 Figure 9-1. An schematic of a damage zone surrounded by undamaged material. Damage Zone 151 (a) (b) Figure 9-2 (a) Dimensions of the OCT specimen and (b) modified OCT specimen. a = 33 mm h = 1 0 4 m m d = 19.1 mm W = 80.6 mm 106 mm 8 4 .6 5 m m c = 38.7 mm Initial Notch blunted notch tip notch root radii (ρ) a = 33 mm 152 Figure 9-3. Modified OCT specimens with varying notch radii (a) 0.5mm, (b) 2mm, (c) 5mm, (d) 12.75mm, (e) 19mm and (f) 27mm.(McClennan 2004) a) d) c) f) e) b) 153 Figure 9-4. Predicted load-CMOD curves compared to the experimental data for the sharp-notched OCT specimen (McClennan 2004). 0 4 8 12 16 20 0 1 2 3 4 5 P ( k N ) CMOD (mm) Exp FE FE Experiment 154 Figure 9-5. Predicted load-CMOD curves compared to the experimental data for the modifiedOCTwithρ=12.7mm (McClennan 2004). 0 4 8 12 16 20 0 1 2 3 4 5 CMOD (mm) P ( k N ) Exp FE FE Experiment 155 Figure 9-6. Predicted load-CMOD curves compared to the experimental data for the modified OCT with ρ=16mm (McClennan 2004). (d) (e) 0 4 8 12 16 20 0 1 2 3 4 5 CMOD (mm) P ( k N ) Exp FE FE Experiment 156 Figure 9-7. Comparison between the predicted load-displacements for the sharp-notched specimen and the blunt-notched specimen withρ=16mm (McClennan 2004). 0 4 8 12 16 20 0 1 2 3 4 5 P ( k N ) CMOD (mm) Exp-Sharp FE-Sharp Exp-16mm FE-16mm Sharp Experiment Sharp FE Blunt Experiment (ρ=16.0mm) Blunt FE (ρ=16.0mm) 157 Figure 9-8. Predicted peak loads versus notch root radius compared to measured test data (McClennan 2004). 10 12 14 16 18 20 0 5 10 15 20 25 30 P ea k L o ad ( k N ) (mm) FE Test 158 Figure 9-9. Predicted length of damage zone and length of zone with saturated damage shown over the load-CMOD for the sharp- notched specimen. 0 4 8 12 16 0 0.5 1 1.5 2 2.5 3 3.5 CMOD (mm) P ( k N ) 0 5 10 15 20 25 30 D am ag e L en g th ( m m ) FE Initiation Saturation 159 Figure 9-10. Predicted length of damage zone and length of zone with saturated damage shown over the load-CMOD for the blunt- notchedspecimenwithρ=12.7 mm. 0 4 8 12 16 20 0 0.5 1 1.5 2 2.5 3 3.5 CMOD (mm) P ( k N ) 0 5 10 15 20 25 30 D am ag e L en g th ( m m ) FE Initiation Saturation 160 Figure 9-11. Predicted length of damage zone and length of zone with saturated damage shown over the load-CMOD for the blunt- notchedspecimenwithρ=22. mm. 0 4 8 12 16 20 0 0.5 1 1.5 2 2.5 3 3.5 CMOD (mm) P ( k N ) 0 5 10 15 20 25 30 D am ag e L en g th ( m m ). FE Initiation Saturation 161 Figure 9-12. Predicted matrix damage pattern in the 90º layer for specimens with notch radii of (a) 0.75mm, (b) 3.2mm, (c) 12.7mm and (d) 22.5mm. 1 0 m m (a) (b) (c) (d) 162 Figure 9-13. Predicted damage height evolution in front of the notch for specimens with the very sharp and the very blunt notch root radii compared to experimental measurements (McClennan 2004). 0 4 8 12 16 20 24 0 5 10 15 20 25 30 Distance from notch tip (mm) D am ag e H ei g h t (m m ) FE-Sharp FE-27 Exp-Sharp Exp-27 163 Figure 9-14. The position of the cross-sections C1 and C2 at which damage profile is plotted, (b) Predicted damage parameter profile corresponding to matrix cracking for specimens with different notch root radii at cross-section C1 immediately in front of the notch and (c) at cross section C2, 20mm away from the original notch tip. 0 0.2 0.4 0.6 0.8 1 -15 -10 -5 0 5 10 15 Position(mm) D a m a g e FE-6 FE-25 FE-32 FE-45 0 0.2 0.4 0.6 0.8 1 -15 -10 -5 0 5 10 15 Position(mm) D a m a g e FE-6 FE-25 FE-32 FE-45 C1 C2 (a) (b) (c) ρ=3mm ρ=12.5mm ρ=16mm ρ=27.5mm ρ=3mm ρ=12.5mm ρ=16mm ρ=27.5mm 164 (a) (b) Figure 9-15. Predicted matri damage patterns from simulation of a specimen with ρ=27. mmusing(a) local and (b) non-local damage models. 165 Figure 9-16. (a) and (b) Predicted matrix damage patterns in the 0º layer using local and non-local models, respectively. (c) and (d) predicted longitudinal stress,σy distribution in front of the crack using local and non-local damage models. (a) (b) (c) (d) 166 Figure 9-17. Schematic of the open-hole specimen geometries: (a) square and (b) tall specimen. σ σ W=5D D H = 5 D D W=5D H = 2 0 D σ σ (a) (b) 167 Figure 9-18. Assumed R-Curve behaviour for the progression of delamination in the open-hole tensile specimen. The solid lines show the estimated value of Gc in each zone and the dashed line shows the increasing trend of the R-Curve. The right hand side axis shows the delamination initiation strain corresponding to the fracture energy on the left axis The correspondence is through the non-linear relation shown in Equation (5.25). Z o n e 1 Z o n e 2 Z o n e 3 Z o n e 4 Z o n e 5E d g e 0 200 400 600 800 0 2 4 6 Crack Length, a (mm) G c (J /m 2 ) 0.30 0.43 0.52 0.60 εi D (% ) 168 Figure 9-19. Predicted nominal stress vs nominal strain curves for the open hole tensile specimens with four different hole sizes. The predicted curves for both square and tall specimen geometries are shown in each case. 0 100 200 300 400 500 0 0.003 0.006 0.009 0.012 0.015 25.4 Tall 25.4 Square 0 100 200 300 400 500 0 0.003 0.006 0.009 0.012 0.015 3.175 Tall 3.175 Square 0 100 200 300 400 500 0 0.003 0.006 0.009 0.012 0.015 12.7 Tall 12.7 Square 0 100 200 300 400 500 0 0.003 0.006 0.009 0.012 0.015 6.35 Tall 6.35 Square N o m in a l S tr e s s ( M P a ) Nominal Strain 169 Figure 9-20. Predicted nominal stress vs nominal strain curves for the open hole tensile specimens with the hole radius of 3.175mm. The predictions are compared to the experiments performed by Hallett et al. (Hallett et al. 2009). 0 100 200 300 400 500 0 0.003 0.006 0.009 0.012 0.015 Nominal Strain N o m in al S tr es s (M P a) 3.175 Experiment 3.175 Square 3.175 Tall 170 Figure 9-21. Predicted progression of damage in the smallest specimen (square specimen with a central hole size of 3.175mm) 0 100 200 300 400 500 0 0.003 0.006 0.009 0.012 0.015 Nominal Strain N o m in a l S tr e s s ( M P a ) Nominal Strain=0.36% Nominal Strain=0.57% Delamination Delamination 90º Matrix Damage 90º Matrix Damage Nominal Strain=1.25% 0º Fibre Breakage 90º Matrix Damage N o m in a l S tr e s s ( M P a ) 171 Figure 9-22. Example of a predicted stress-strain curve (for the D3.175S specimen) showing the failure strength (defined as the instant when the laminate secant stiffness drops by 20%) and the ultimate strength. 0 20 40 60 80 0 0.002 0.004 0.006 0.008 0.01 Nominal Strain S e c a n t M o d u lu s ( G P a ) 20% drop in modulus S e c a n t M o d u lu s ( G P a ) 0 100 200 300 400 500 0 0.003 0.006 0.009 0.012 Nominal Strain N o m in a l S tr e s s ( M P a ) Ultimate strength Delamination growth Failure strength at 20% drop in modulus N o m in a l S tr e s s ( M P a ) 172 Figure 9-23. Predicted Failure and Peak Stresses versus Hole Size. 0 100 200 300 400 500 0 10 20 30 Hole Size (mm) S tr e s s ( M P a ) Failure Stress Peak Stress 173 Figure 9-24. Predicted delamination area in different open-hole tensile specimen. Hole sizes: 3.175mm 6.35mm 12.7mm 25.4mm Delamination 174 10. SUMMARY, CONCLUSIONS AND FUTURE WORK 10.1 SUMMARY Computational modelling of damage and fracture in laminated composite structures is very challenging. Pioneering researchers around the world have suggested various approaches to simulate the details of the complex damage process in laminated composites. Regardless of the validity and success of such model in their predictions, they are still not being used by the engineers in practice due to their increasing level of complexity and high computational cost. The non-local damage model proposed in this thesis aims to address this issue by developing a model that is physically-based and computationally efficient and yet relatively simple to understand, calibrate and use. Chapter 2 in this thesis reviewed various approaches to simulation of damage and fracture in materials/structures. Advantages and shortcomings of the Continuum, Discrete and Enriched FE formulations were discussed. Chapter 3 took a close look at the localization problem in strain softening models. The so- called localization limiter methods were discussed in detail in this chapter. Chapter 4 presented an overview of the current approaches to damage simulation in laminated composites. This chapter looked at modelling of damage in laminated composites from both physical and numerical perspectives. Chapter 5 was devoted to presentation of the UBC COmposite DAmage Model (CODAM). The sub-laminate approach and the first generation of CODAM were introduced. The origins of the objectivity problems associated with the formulation of the original CODAM were discussed in detail and examples were presented to demonstrate the problem. The second generation of CODAM, i.e. CODAM2 that has been developed by the author was then introduced. The formulation of CODAM2 was presented and its approach towards formulating the damaged stiffness and the damage evolution law were discussed in detail. It was also shown that the formulation CODAM2 is thermodynamically objective. Chapter 6 introduced the orthotropic averaging scheme developed by the author. This chapter discussed the motivations behind an orthotropic averaging scheme for materials with highly 175 orthotropic constituents (like laminated composites) from a physical point of view. The formulations for the orthotropic averaging scheme which uses two non-local radii were presented in this chapter. The effect of the non-local radii on the predicted damage pattern and structural behaviour of the laminated composites were discussed. It was shown that the orthotropic averaging scheme offers a significant improvement in predicted crack direction and damage pattern in specimens made from unidirectional laminates. Chapter 7 provided verifications for the implementation of CODAM2 in the open-source FE code, OOFEM. Examples were shown to verify the performance of the implemented code. These examples included comparison between the predictions of CODAM2 in elastic regime with predictions of ABAQUS, simulation of a single element (and therefore predicted stress- strain behaviour at an integration point) for various loading scenarios. Performance of the non- local averaging scheme was also examined in an example with uniform strain distribution. Using another example, it was shown that CODAM2 successfully predicts the traction free condition for a widely opened crack that verifies the performance of the anti-stress locking mechanism in this model. Chapter 8 discussed the calibration procedure for parameters required in CODAM2. The initiation and saturation strains for fibre and matrix damage mechanisms as well as the two non-local radii are the parameters used in this model. It was shown how the measurements of a notched specimen that exhibits stable and self-similar crack growth, like OCT, can be used to calibrate CODAM2. Chapter 9 presented some numerical examples to demonstrate the capabilities of CODAM2 in predicting damage behaviour in various specimens. The performance of CODAM2 was validated by comparing the model predictions with experimental measurements for modified OCT specimens with a range of sharp to blunt notches. This chapter also included the numerical simulation of open-hole tension specimens with various sizes. The latter example was part of the contribution of UBC Composites Group to the third round of Worldwide Failure Exercise for composites (WWFE-III). 176 10.2 CONCLUDING REMARKS This thesis presented the latest revision of the UBC COmposite DAmage Model (CODAM). The newly developed CODAM2, improves the previously available models in the following ways: - Formulation It was shown that CODAM2 addresses the spurious dependency of the original CODAM on the choice of the coordinate system by explicit incorporation of the directions of the plies into the formulation. This model is capable of simulating formation and propagation of damage in quasi-isotropic as well as orthotropic laminates. CODAM2 is also capable of modelling the damage-induced orthotropy in laminated composites. The relatively simple formulation of CODAM2 and its straight-forward calibration procedure makes it convenient for engineers to use in practical applications. - Orthotropic Averaging Scheme The orthotropic non-local averaging has been incorporated into CODAM2 to address the localization problem and improve the predicted pattern of damage and direction of crack growth. It was shown that CODAM2 equipped with orthotropic non-local averaging is capable of accurately predicting the damage pattern and overall load-displacement behaviour in a variety of specimen geometries and loading conditions. Successful prediction of the behaviour of modified OCT specimens with various notch root radii showed that the developed model has the capability of predicting the whole spectrum of energy-driven damage growth (sharp- notched specimen) to stress-driven catastrophic damage propagation (blunt-notched specimen). 10.3 FUTURE WORK To make the model capable of simulating a wider range of damage propagation problems in laminated composites, the model needs to be improved further. The following are some of the features that will enhance the model capabilities. 177 - Permanent Deformations One of the enhancements would be to add the capability of simulating the plateau behaviour and residual deformation observed in compressive tests on notched specimens made from laminated composites (McGregor et al. 2008; Zobeiry 2004). Residual deformation was also observed in low velocity lateral impact tests on composite panels (Delfosse et al. 1995; Delfosse and Poursartip 1997). To be able to predict the residual deformations, a coupled damage-plasticity formulation needs to be employed for modelling the behaviour of the material under compressive loading. The damage with plasticity model employed in (Forghani and Vaziri 2009) for an isotropic material was shown to be successful in predicting the unloading behaviour of the impacted plate as well as predicting the residual deformation. The challenge here is to adopt an appropriate plasticity model that is compatible with the fundamentals of CODAM2. - Delamination Another piece of the puzzle that needs to be worked on is to develop a more general solution for simulation of delamination within the framework of the sub-laminate-based approach. As discussed earlier, the basic assumption of compatibility of displacement/strain fields in the sub-laminate approach is not applicable to delamination-dominated responses. One of the possible solutions to this problem would be to adopt the combined continuum- discontinuous PUM-based approach (e.g. (Remmers et al. 2003)) into the sub-laminate-based framework. Although this solution can potentially solve the delamination problem, it will add complications to the model and make it less practical. - Change in the Laminate Lay-up One of the advantages of laminated composites is their tailorability. The new manufacturing techniques such as tape laying method enables the designers to perform an extremely optimized design by dynamically changing the lay-up and orientation of the layers. The challenge here is how to take these changes in the number of layers and their orientations into account in a damage model that intends to predict the ultimate behaviour of such structures. It should be noted that even the ply-based models, in spite of their high level of detail and through-thickness discretization, are incapable of simulating the effect of dynamic change in the orientation of the layers in the laminate. 178 REFERENCES Alfano, G., and Crisfield, M. A. (2001). "Finite Element Interface Models for the Delamination Analysis of Laminated Composites: Mechanical and Computational Issues." Int J Numer Methods Eng, 50(7), 1701-36. ASME V&V-10. (2006). Guide for Verification and Validation in Computational Solid Mechanics, 1st Ed., American Society of Mechanical Engineers, New York. Babuska, I., and Melenk, J. M. (1997). "The Partition of Unity Method." Int J Numer Methods Eng, 40(4), 758. Bazant, Z. P. (2010). "Can Multiscale-Multiphysics Methods Predict Softening Damage and Structural Failure?" Int J Mult Comp Eng, 8(1), 61-67. Bazant, Z. P., and Caner, F. C. (2005). "Microplane Model M5 with Kinematic and Static Constraints for Concrete Fracture and Anelasticity. I: Theory." J. Eng. Mech., 131(1), 31-40. Bazant, Z. P., and Jirasek, M. (2002). "Nonlocal Integral Formulations of Plasticity and Damage: Survey of Progress." J. Eng. Mech., 128(11), 1119-1149. Bazant, Z. P., and Pijaudier-Cabot, G. (1988). "Nonlocal continuum damage, localization instability and convergence." Proc., Preprint - American Society of Mechanical Engineers, Publ by American Soc of Mechanical Engineers (ASME), , 18-7. Bazant, Z. P., and Belytschko, T. B. (1985). "Wave Propagation in a Strain-Softening Bar: Exact Solution." J. Eng. Mech., 111(3), 381-389. Bazant, Z. P., and Oh, B. H. (1985). "Microplane Model for Progressive Fracture of Concrete and Rock." J. Eng. Mech., 111(4), 559-582. Bazant, Z. P., Belytschko, T. B., Chang, T. (1984). "Continuum Theory for Strain-Softening." J. Eng. Mech., 110(12), 1666-1692. Bazant, Z. P., and Oh, B. H. (1983). "Crack Band Theory for Fracture of Concrete." Materiaux Et Constructions, Materials and Structures, 16(93), 155-177. Beghini, A., Cusatis, G., Bazant, Z. P. (2008). "Spectral Stiffness Microplane Model for Quasibrittle Composite Laminates - Part II: Calibration and Validation." J. Appl. Mech. Trans. ASME, 75(2), 0210101-0210106. Belytschko, T., and Black, T. (1999). "Elastic Crack Growth in Finite Elements with Minimal Remeshing." Int J Numer Methods Eng, 45(5), 601-620. Belytschko, T., and Bazant, Z. P. (1984). "Strain-softening materials and finite element solutions." Proc., Constitutive Equations: Macro and Computational Aspects. Presented at the Winter Annual Meeting of the American Society of Mechanical Engineers. ASME, , 253-272. Belytschko, T., and Lasry, D. (1989). "Study of Localization Limiters for Strain-Softening in Statics and Dynamics." Comput. Struct., 33(3), 707-715. 179 Belytschko, T., Fish, J., Engelmann, B. E. (1988). "Finite Element with Embedded Localization Zones." Comput. Methods Appl. Mech. Eng., 70(1), 59-89. Borg, R., Nilsson, L., Simonsson, K. (2001). "Simulation of Delamination in Fiber Composites with a Discrete Cohesive Failure Model." Compos. Sci. Technol., 61(5), 667-677. Bouchard, P. O., Bay, F., Chastel, Y., Tovena, I. (2000). "Crack Propagation Modelling using an Advanced Remeshing Technique." Comput. Methods Appl. Mech. Eng., 189(3), 723-742. Brenner, S. C., and Scott, L. R. (2008). Mathematical Theory of Finite Element Methods, Springer, New York, NY, USA. Camacho, G. T., and Ortiz, M. (1996). "Computational Modelling of Impact Damage in Brittle Materials ." International Journal of Solids and Structures, 33(20-22), 2899. Camanho, P. P., Davila, C. G., Pinho, S. T., Iannucci, L., Robinson, P. (2006). "Prediction of in Situ Strengths and Matrix Cracking in Composites Under Transverse Tension and in-Plane Shear." Compos Part A Appl Sci Manuf, 37(2), 165-76. Camanho, P. P., Davila, C. G., de Moura, M. F. (2003). "Numerical Simulation of Mixed- Mode Progressive Delamination in Composite Materials." J Compos Mater, 37(16), 1415- 1438. Carpinteri, A. (1989). "Post-Peak and Post-Bifurcation Analysis of Cohesive Crack Propagation." Eng. Fract. Mech., 32(2), 265-278. Carpinteri, A., and Colombo, G. (1989). "Numerical Analysis of Catastrophic Softening Behaviour (Snap-Back Instability)." Comput. Struct., 31(4), 607-636. Chamis, C. C., Murthy, P. L. N., Gotsis, P. K., Mital, S. K. (2000). "Telescoping Composite Mechanics for Composite Behavior Simulation." Comput. Methods Appl. Mech. Eng., 185(2- 4), 399-411. Chamis, C. C., Murthy, P. L. N., Minnetyan, L. (1996). "Progressive Fracture of Polymer Matrix Composite Structures." Theor. Appl. Fract. Mech., 25(1), 1-15. Chow, C. L., and Wang, J. (1987). "ANISOTROPIC THEORY OF ELASTICITY FOR CONTINUUM DAMAGE MECHANICS." Int. J. Fract., 33(1), 3-16. Cox, B., and Yang, Q. (2006). "In Quest of Virtual Tests for Structural Composites." Science, 314(5802), 1102-1107. Cusatis, G., Beghini, A., Bazant, Z. P. (2008). "Spectral Stiffness Microplane Model for Quasibrittle Composite Laminates - Part I: Theory." J. Appl. Mech. Trans. ASME, 75(2), 0210091-0210099. de Borst, R., Sluys, L. J., Muhlhaus, H. -., Pamin, J. (1993). "Fundamental Issues in Finite Element Analyses of Localization of Deformation." Eng. Comput., 10(2), 99-121. de Borst, R., and Nauta, P. (1985). "Non-Orthogonal Cracks in Smeared Finite Element Model." Eng. Comput., 2(1), 35-46. de Borst, R., and Remmers, J. J. C. (2006). "Computational Modelling of Delamination." Composites Sci. Technol., 66(6), 713-722. 180 Delfosse, D., and Poursartip, A. (1997). "Energy-Based Approach to Impact Damage in CFRP Laminates." Compos Part A Appl Sci Manuf, 28(7), 647-655. Delfosse, D., Poursartip, A., Coxon, B. R., Dost, E. F. (1995). "Non-penetrating impact behavior of CFRP at low and intermediate velocities." Proc., Proceedings of the 5th Symposium on Composite Materials: Fatigue and Fracture, ASTM, Philadelphia, PA, USA, , 333-350. Dolbow, J., Moes, N., Belytschko, T. (2001). "An Extended Finite Element Method for Modeling Crack Growth with Frictional Contact." Comput. Methods Appl. Mech. Eng., 190(51-52), 6825-6846. Dvorkin, E. N., Cuitino, A. M., Gioia, G. (1990). "Finite Elements with Displacement Interpolated Embedded Localization Lines Insensitive to Mesh Size and Distortions." Int J Numer Methods Eng, 30(3), 541-564. Elices, M., Guinea, G. V., Gomez, J., Planas, J. (2002). "The Cohesive Zone Model: Advantages, Limitations and Challenges." Eng. Fract. Mech., 69(2), 137-163. Espinosa, H. D., Dwivedi, S., Lu, H. -. (2000). "Modeling Impact Induced Delamination of Woven Fiber Reinforced Composites with contact/cohesive Laws." Comput. Methods Appl. Mech. Eng., 183(3-4), 259-290. Floyd, A. (2004). "An Engineering Approach to the Simulation of Gross Damage Development in Composite Laminates." PhD Thesis, Deparment of Civil Engineering, the University of British Columbia, Vancouver, BC, Canada, . Floyd, A., Cronin, D., Worswick, M., Poursartip, A., Vaziri, R., Plumtree, A. (2001). Three- Dimensional Numerical Simulation of the Ballistic Response of Laminated Composite Plates, Report submitted to Defence Research Establishment Valcartier (DREV) under Contract No. W7701-9-3228/001/XSK, Canada. Forghani, A., Zobeiry, N., Poursartip, A., Vaziri, R. (2011). "A Structural Modeling Framework for Prediction of Damage Development and Failure of Composite Laminates." Composites Sci. Technol., . Forghani, A., and Vaziri, R. (2009). "Computational Modeling of Damage Development in Composite Laminates Subjected to Transverse Dynamic Loading." J. Appl. Mech. Trans. ASME, 76(5), 051304. Forghani, A., McGregor, C., Zobeiry, N., Vaziri, R., Ellyin, F., Poursartip, A. (2007). Modeling of Effect of Blast Load on Composite Panels, 1st Ed., UBC-DRDC Contract No. W7701-062619/001/QCL, Canada. Forghani, A. (2005). "Application of the Extended Finite Element Method (X-FEM) in Modelling Cohesive Cracks." MSc Thesis, Department of Civil Engineering, Sharif University of Technology, Tehran, Iran, . Gadala, M. S., Movahhedy, M. R., Wang, J. (2002). "On the Mesh Motion for ALE Modeling of Metal Forming Processes." Finite Elements Anal. Des., 38(5), 435-459. Gonzalez, C., and LLorca, J. (2006). "Multiscale Modeling of Fracture in Fiber-Reinforced Composites." Acta Materialia, 54(16), 4171-4181. 181 Grufman, C., and Ellyin, F. (2008). "Numerical Modelling of Damage Susceptibility of an Inhomogeneous Representative Material Volume Element of Polymer Composites." Composites Sci. Technol., 68(3-4), 650-657. Grufman, C., and Ellyin, F. (2007). "Determining a Representative Volume Element Capturing the Morphology of Fibre Reinforced Polymer Composites." Composites Sci. Technol., 67(3-4), 766-775. Hallett, S. R., Green, B. G., Jiang, W. G., Wisnom, M. R. (2009). "An Experimental and Numerical Investigation into the Damage Mechanisms in Notched Composites." Compos Part A Appl Sci Manuf, 40(5), 613-624. Hinton, M. J., Kaddour, A. S., Soden, P. D. (2004). "A further Assessment of the Predictive Capabilities of Current Failure Theories for Composite Laminates: Comparison with Experimental Evidence." Composites Sci. Technol., 64(3-4), 549-588. Hinton, M. J., Kaddour, A. S., Soden, P. D. (2002). "A Comparison of the Predictive Capabilities of Current Failure Theories for Composite Laminates, Judged Against Experimental Evidence." Composites Sci. Technol., 62(12-13), 1725-1797. Jirasek, M. (2007). "Course Notes on Modeling of Localized Elastic Deformation." . Jirasek, M., Rolshoven, S., Grassl, P. (2004). "Size Effect on Fracture Energy Induced by Non-Locality." Int. J. Numer. Anal. Methods Geomech., 28(7), 653-670. Jirasek, M., and Zimmermann, T. (1998). "Rotating Crack Model with Transition to Scalar Damage." J ENG MECH-ASCE, 124(3), 277-284. Jirasek, M. (2000). "Comparative Study on Finite Elements with Embedded Discontinuities." Comput. Methods Appl. Mech. Eng., 188(1), 307-330. Jirasek, M. (1998). "Nonlocal Models for Damage and Fracture: Comparison of Approaches." Int J Solids Struct, 35(31), 4133-4145. Kaddour, A. S., Hinton, M. J., Li, S., Smith, P. (2009). "Instructions to Contributors of the Third World-Wide Failure Exercise (WWFE-III): Part (A)." . Kaddour, A. S., Hinton, M. J., Soden, P. D. (2004). "A Comparison of the Predictive Capabilities of Current Failure Theories for Composite Laminates: Additional Contributions." Compos. Sci. Technol., 64(3-4), 449-476. Kashtalyan, M., and Soutis, C. (2000). "Effect of Delaminations Induced by Transverse Cracks and Splits on Stiffness Properties of Composite Laminates." Compos Part A Appl Sci Manuf, 31(2), 107-119. Kennedy, T. C., and Nahan, M. F. (1996). "Simple Nonlocal Damage Model for Predicting Failure of Notched Laminates." Composite Structures, 35(2), 229-236. Kongshavn, I., and Poursartip, A. (1999). "Experimental Investigation of a Strain-Softening Approach to Predicting Failure in Notched Fibre-Reinforced Composite Laminates." Compos. Sci. Technol., 59(1), 29-40. Krajcinovic, D. (1985). "Continuous damage mechanics revisited: Basic concepts and definitions." Proc., Winter Annual Meeting - American Society of Mechanical Engineers. ASME, , ASME, New York, NY, USA. 182 Krueger, R. (2004). "Virtual Crack Closure Technique: History, Approach, and Applications." Appl. Mech. Rev., 57(1-6), 109-143. Lacy, T. E., McDowell, D. L., Willice, P. A., Talreja, R. (1997). "On Representation of Damage Evolution in Continuum Damage Mechanics." Int J Damage Mech, 6(1), 62-95. Ladeveze, P. (1995). "Damage Computational Approach for Composites: Basic Aspects and Micromechanical Relations." Comput. Mech., 17(1-2), 142-150. Ladeveze, P. (2004). "Multiscale Modeling and Computational Strategies for Composites." Int. J. Numer. Meth. Engng, 60, 233-253. Ladeveze, P., Allix, O., Deu, J., Leveque, D. (2000). "Mesomodel for Localisation and Damage Computation in Laminates." Comput. Methods Appl. Mech. Eng., 183(1), 105-122. Lasry, D., and Belytschko, T. (1988). "Localization Limiters in Transient Problems." Int J Solids Struct, 24(6), 581-597. Lemaitre, J., and Desmorat, R. (2005). Engineering Damage Mechanics : Ductile, Creep, Fatigue and Brittle Failures, Springer-Verlag, Heidelberg, , DEU. Lemaitre, J., Desmorat, R., Sauzay, M. (2000). "Anisotropic Damage Law of Evolution." Eur J Mech A Solids, 19(2), 187-208. Maimi, P., Camanho, P. P., Mayugo, J. A., Davila, C. G. (2007). "A Continuum Damage Model for Composite Laminates: Part I - Constitutive Model." Mech. Mater., 39(10), 897-908. Matzenmiller, A., Lubliner, J., Taylor, R. L. (1995). "Constitutive Model for Anisotropic Damage in Fiber-Composites." Mech. Mater., 20(2), 125-125. Mazars, J., and Pijaudier-Cabot, G. (1996). "From Damage to Fracture Mechanics and Conversely: A Combined Approach." Int J Solids Struct, 33(20-22), 3327-3342. McClennan, S., Vaziri, R., Poursartip, A. (2003). Incorporation of Strain-Rate Modelling Capability in a Damage Mechanics Based Constitutive Model for Composites, 1st Ed., Report submitted to GM Company, . McClennan, S. A. (2004). "Crack Growth and Damage Modeling of Fibre Reinforced Polymer Composites." MAc Thesis, Departments of Materials Engineering, the University of British Columbia, Vancouver, BC, Canada, . McGregor, C. (2005). "Simulation of Progressive Damage Development in Braided Composite Tubes Under Axial Compression." M. A. Sc. Thesis, Department of Civil Engineering, the University of British Columbia, Vancouver, BC, Canada, . McGregor, C., Zobeiry, N., Vaziri, R., Poursartip, A. (2008). "A Constitutive Model for Progressive Compressive Failure of Composites." J. Composite Mater., 42(25), 2687-2716. McGregor, C. J., Vaziri, R., Poursartip, A., Xiao, X. (2007). "Simulation of Progressive Damage Development in Braided Composite Tubes Under Axial Compression." Composites Part A: Applied Science and Manufacturing, 38(11), 2247-2259. Moes, N., and Belytschko, T. (2002). "Extended Finite Element Method for Cohesive Crack Growth." Eng. Fract. Mech., 69(7), 813-833. 183 Movahhedy, M., Gadala, M. S., Altintas, Y. (2000). "Simulation of the Orthogonal Metal Cutting Process using an Arbitrary Lagrangian–Eulerian Finite-Element Method." J. Mater. Process. Technol., 103(2), 267-275. Needleman, A. (1988). "Material Rate Dependence and Mesh Sensitivity in Localization Problems." Comput. Methods Appl. Mech. Eng., 67(1), 69-85. O'Brien, T. K. (1981). Characterization of Delamination Onset and Growth in a Composite Laminate, NASA Technical Memorandum 81940, . Ortiz, M., Leroy, Y., Needleman, A. (1987). "Finite Element Method for Localized Failure Analysis." Comput. Methods Appl. Mech. Eng., 61(2), 189-214. Patzak, B., and Bittnar, Z. (2001). "Design of Object Oriented Finite Element Code." Adv. Eng. Software, 32(10-11), 759-767. Peerlings, R. H. J., de Borst, R., Brekelmans, W. A. M., Geers, M. G. D. (2002). "Localisation Issues in Local and Nonlocal Continuum Approaches to Fracture." Eur J Mech A Solids, 21(2), 175-189. Peerlings, R. H. J., Geers, M. G. D., de Borst, R., Brekelmans, W. A. M. (2001). "A Critical Comparison of Nonlocal and Gradient-Enhanced Softening Continua." Int. J. Solids Structures, 38(44), 7723-7746. Peerlings, R. H. J., de Borst, R., Brekelmans, W. A. M., Geers, M. G. D. (1998). "Wave Propagation and Localisation in Nonlocal and Gradient-Enhanced Damage Models." Journal De Physique. IV : JP, 8(8), 293-300. Remmers, J. J., Wells, G. N., de Borst, R. (2003). "A Solid-Like Shell Element Allowing for Arbitrary Delaminations." Int J Numer Methods Eng, 58(13), 2013-40. Réthoré, J., Gravouil, A., Combescure, A. (2004). "A Stable Numerical Scheme for the Finite Element Simulation of Dynamic Crack Propagation with Remeshing." Comput. Methods Appl. Mech. Eng., 193(42-44), 4493-4510. Ruiz, G., Ortiz, M., Pandolfi, A. (2000). "Three-Dimensional Finite-Element Simulation of the Dynamic Brazilian Tests on Concrete Cylinders." Int J Numer Methods Eng, 48(7), 963-994. Schipperen, J. H. A., and de Borst, R. (2001). "A Numerical Analysis of Mixed-Mode Delamination in Carbon-Epoxy Prepregs." Compos. Struct., 54(4), 445-451. Sluys, L. J. (1992). "Wave Propagation, Localization and Dispersion in Softening Solids." PhD Thesis, Delft University of Technology, Netherlands, . Sluys, L. J., and de Borst, R. (1992). "Wave Propagation and Localization in a Rate- Dependent Cracked Medium - Model Formulation and One-Dimensional Examples." Int. J. Solids Structures, 29(23), 2945-2958. Soden, P. D., Hinton, M. J., Kaddour, A. S. (1998). "Comparison of the Predictive Capabilities of Current Failure Theories for Composite Laminates." Composites Sci. Technol., 58(7), 1225- 1254. Talreja, R. (1989). "Damage Development in Composites: Mechanisms and Modelling." J. Strain Anal. Eng. Des., 24(4), 215-222. 184 Talreja, R. (1985). "Continuum Mechanics Characterization of Damage in Composite Materials." Proc. R. Soc. Lond. , Ser. A: Math. Phys. Sci., 399(1817), 195-216. Talreja, R. (2006). "Multi-Scale Modeling in Damage Mechanics of Composite Materials." J. Mater. Sci., 41(20), 6800-6812. Toray. (2010). "Toray group net sales ratio by application." <http://www.toray.com/ir/individual/ind_012.html> (03/03, 2011). Wells, G. N., de Borst, R., Sluys, L. J. (2002). "A Consistent Geometrically Non-Linear Approach for Delamination." Int J Numer Methods Eng, 54(9), 1333-55. Williams, K. V., Vaziri, R., Poursartip, A. (2003). "A Physically Based Continuum Damage Mechanics Model for Thin Laminated Composite Structures." Int J Solids Struct, 40(9), 2267- 300. Williams, K. V. (1998). "A Physically-Based Continuum Damage Mechanics Model for Numerical Prediction of Damage Growth." PhD Thesis, Department of Materials Engineering, the University of British Columbia, Vancouver, BC, Canada, . Xia, Z., Zhou, C., Yong, Q., Wang, X. (2006). "On Selection of Repeated Unit Cell Model and Application of Unified Periodic Boundary Conditions in Micro-Mechanical Analysis of Composites." Int J Solids Struct, 43(2), 266-278. Xia, Z., Chen, Y., Ellyin, F. (2000). "Meso/micro-Mechanical Model for Damage Progression in Glass-fiber/epoxy Cross-Ply Laminates by Finite-Element Analysis." Composites Sci. Technol., 60(8), 1171-1179. Yang, Q., and Cox, B. (2005). "Cohesive Models for Damage Evolution in Laminated Composites." Int. J. Fract., 133(2), 107-137. Yang, Z., and Chen, J. (2004). "Fully Automatic Modelling of Cohesive Discrete Crack Propagation in Concrete Beams using Local Arc-Length Methods." Int. J. Solids Structures, 41(3-4), 801-826. Zhang, Y., Xia, Z., Ellyin, F. (2005). "Viscoelastic and Damage Analyses of Fibrous Polymer Laminates by micro/meso-Mechanical Modeling." J. Composite Mater., 39(22), 2001-2022. Zhang, Y., Xia, Z., Ellyin, F. (2004). "Evolution and Influence of Residual stresses/strains of Fiber Reinforced Laminates." Composites Sci. Technol., 64(10-11), 1613-1621. Zobeiry, N. (2010). "Extracting the Strain Softening Response of Composites using Full-Field Displacement Measurement." PhD Thesis, the University of British Columbia, . Zobeiry, N., Forghani, A., McGregor, C., Vaziri, R., Poursartip, A. (2008). "Progressive Damage Modeling of Composite Materials Under both Tensile and Compressive Loading Regimes." MechanicalResponseofComposites’,Series:ComputationalMethods inApplied Sciences, Vol. 10, Camanho, P.P. Davila, C.G., Pinho, S.T.; Remmers, J.J.C. (Eds.), Porto, Portugal Ed., Springer, 179-195. Zobeiry, N. (2004). "Progressive Damage Modeling of Composite Materials Under Compressive Loads." M. A. Sc. Thesis, Department of Civil Engineering, the University of British Columbia, Vancouver, BC, Canada, . 185 APPENDICES A. A GUIDELINE FOR USERS A.1 GENERAL NOTES The user of this code should be aware of the scope of the developed model. CODAM2 local and non-local are designed to capture the growth of damage and its progress leading to ultimate structural failure. The model is applicable to cases where delamination is not widely extended and is not the dominant mode of damage. CODAM2 does not intend to predict the details of damage formation and propagation in the layers of the laminate. The model in its current version is not capable of simulating the residual strain that is observed in compression. See Section 9.3 for details on Scope, Applicability and limitations of the presented non-local damage model. The methodology for calibration is described in Chapter 8. A.2 COMPILING AND RUNNING THE CODE Since the source code for OOFEM written in C++ is provided, the code is theoretically platform independent. See the instruction provided on the OOFEM website: www.oofem.org for information on how to compile the code on various platforms. The command line format for executing OOFEM after compilation is: oofem -f InputFileName (-rn) By applying the optional parameter –rn, OOFEM performs a renumbering on all the degrees of freedom to optimize the computation time. OOFEM website providing the guidelines on how to prepare an input file is: http://www.oofem.org/resources/doc/oofemInput/html/oofemInput.html. CODAM2 and CODAM2NL are the keywords to call the local and non-local CODAM2 material models, respectively. Here are examples of the command line to call CODAM2 and CODAM2NL that needs to be inserted in the input file. CODAM2 1 d 0.000002 E 00.1 n 0.28 Ex 150000.0 Ey 11000.0 Ez 0.0001 NYyz 0.3 NYxz 0.0 NYxy 0.3 Gyz 0.0001 Gxz 0.0001 Gxy 4600.0 tAlphax 0.0 tAlphay 0.0 186 tAlphaz 0.0 lcs 6 0.0 1.0 0.0 1.0 0.0 0.0 paramFile sspar1.dat angles 4 0 90 45 -45 talpha 0.00 CODAM2NL 1 d 0.000002 E 00.1 n 0.28 Ex 150000.0 Ey 11000.0 Ez 0.0001 NYyz 0.3 NYxz 0.0 NYxy 0.3 Gyz 0.0001 Gxz 0.0001 Gxy 4600.0 tAlphax 0.0 tAlphay 0.0 tAlphaz 0.0 lcs 6 0.0 1.0 0.0 1.0 0.0 0.0 paramFile sspar1.dat angles 4 0 90 45 -45 talpha 0.00 r1 2.0 r2 1.0 strratio 0.0 regionmap 2 0 1 Parameters: d: Density, E: Elastic modulus of isotropic matrix, n: Poisson’s ratio of the isotropic matrix. These parameters belong to an earlier version of CODAM2 and are there only for debugging purposes. Ex, Ey, Ez: Young’s modulus of the anisotropic material. NYyz, NYxz, NYxy: Major Poisson’s ratios. Gyz, Gxz, Gxy: Shear Moduli. tAlphax, tAlphay, tAlphaz: Thermal expansion coefficients. lcs: Local Coordinate System: the six numbers represent two vectors that define the 1 and 2 coordinates in the physical coordinate system. paramFile: Name of the file that contains initiation and saturation strains. angles: The number of layers and their angle of orientation with respect to the first vector given in lcs. talpha: Thermal expansion coefficient for the isotropic part (CODAM2.4) r1: Longitudinal non-local radius. r2(<= r1): Transverse non-local radius (use r2=r1 to enable an isotropic averaging scheme) strratio: Use Zero. For debugging purposes. regionmap: List of the cross-section numbers that are included or excluded from averaging. Use 0 for the cross-sections that are excluded and 1 for the ones that are included. 187 The parameters file consists of 5 columns of numbers representing the parameters for delamination, overall, fibre and matrix damage as follows: , D i t , L i t , f i t , m i t , D s t , L s t , f s t , m s t , D i c , L i c , f i c , m i c , D s c , L s c , f s c , m s c D L f m Example: 0.0060 0.080 0.011 0.0080 0.1000 0.100 0.080 0.0800 0.0200 0.035 0.010 0.0080 0.5000 0.080 0.035 0.0500 0.2000 1.000 1.000 1.0000 A.3 CODE CALIBRATION The parameters required for running the model needs to be obtained based on the discussion earlier in Chapter 8. The elastic parameters and the damage initiation strains can be obtained from the datasheets provided by the material manufacturers. A test like OCT that exhibits progressive and stable damage growth can be employed to calibrate the non-local radii and damage saturation parameters. 188 B. DEMONSTRATION OF THE LOCALIZATION PHENOMENON IN A 1D PROBLEM The localization phenomenon under quasi-static and dynamic conditions is demonstrated here using a one-dimensional example. A finite element simulation of a 1-D bar is also employed to demonstrate the mesh size dependency. B.1 QUASI-STATIC PROBLEM Consider an axial bar with a uniform cross-section subjected to essential (Dirichlet) boundary conditions on both ends as shown in Figure B-1. The weak form of the equilibrium equation (in the absence of body force) can be written as 0t L du dv E dx dx dx (B.1) Where u is the incremental displacement along the bar, v is an admissible displacement field and E t is the tangent stiffness of the material. Suppose that part of the bar (L d ) is damaged and exhibits negative tangent stiffness and part of it has remained undamaged (L u ) with positive tangent stiffness. The weak form can be re- written as 0 0 u d d L L du dv du dv E dx E dx dx dx dx dx (B.2) where E 0 is the modulus of the undamaged material that has a positive value and E d represents the tangent stiffness of the material in the softening zone. Since the tangent stiffness of the damaged part is determined based on the loading/unloading condition, the sign and value of E d is a function of the strain rate, ( du dx ) and can be positive or negative. In the weak form of this boundary value problem, the bi-linear form A(u ,v) is defined as 0( , ) u d d L L du dv du dv A u v E dx E dx dx dx dx dx (B.3) 189 Based on the Lax-Milgiram theorem, the coercivity condition that is an essential condition for well-posedness of the mathematical representation and uniqueness of the solution is satisfied if 2 ( , ) FA u u C u (B.4) where represents the H 1 norm and CF is a positive constant. See (Brenner and Scott 2008) for more details. Since the tangent stiffness of the material is partly positive and partly negative over the length of the bar, the bi-linear form defined by (B.3) does not satisfy the coercivity condition. Therefore the weak form does not necessarily have a unique solution and the problem described by the weak form is ill-posed. B.2 DYNAMIC PROBLEM The differential equation of the dynamic equilibrium for the 1-D problem in the absence of rate-dependent terms and body force is written as 2 2 2 2 u u x t t (B.5) where is the density of the material. The stress gradient can be written as the product of the tangent stiffness and strain gradient as tE x x (B.6) Substituting (B.6) into (B.5) we obtain the differential equation of the dynamic equilibrium in terms of strain 2 2 2 2 t u uE x t (B.7) Equation (B.7) is the governing equation of an initial value problem that describes the propagation of strain in a one dimensional bar. In an undamaged material, where the sign of the tangent stiffness is positive, this equation is a hyperbolic equation that represents wave propagation with a propagation speed of 190 0 0 Ec (B.8) In a damaged material with a softening behaviour, the sign of the tangent stiffness E t becomes negative and Equation (B.7) loses its hyperbolicity. Under this condition, wave speed becomes imaginary and this equation no longer represents wave propagation in a solid material. As shown here, the continuum damage theory in its local form predicts that the damaged material, even just after the initiation of damage, loses its capability to transfer stress waves. B.3 NUMERICAL EXAMPLE To demonstrate the mesh-size dependency problem, a bar of length 1.0 is simulated with four different spatial discretization using 1, 2, 4, and 8 elements. An explicit dynamic analysis has been performed on this structure. A linear softening model with damage initiation strain of 1.0% and saturation strain of 4.9% has been employed. The simulations have been performed under displacement controlled condition where one end of the bar was pulled under a constant rate. Fibure B-3 shows the predicted force versus elongation curves for all the four mesh sizes. It can be seen that by refining the mesh, the response of the structure changes significantly. When only one element is employed, the structure loses its load carrying capability at an elongation of 0.049. By employing two elements, this value reduces to 0.0245. The model with four elements predicts a failure elongation of 0.01225. The model with eight elements fails abruptly when it reaches the peak force and exhibits a brittle behaviour. In Figure B-3, the response of the model with eight elements was truncated after its failure. Figure B-4 shows the complete response of this eight-element model. It can be seen that after the failure, the system continues to oscillate which indicates that part of the strain energy that was stored in the elastic system is converted into kinetic energy. In the FE models with more than one element, simulations predicted localization of strain in only one of the elements as depicted in Figure B-5. The rest of elements remained undamaged and unloaded after reaching the peak load. Figures B-6 and B-7 show the predicted values of external work and dissipated energy as a function of bar elongation. It can be seen that the predicted patterns for both external work and 191 dissipated energy strongly depend on the number of elements. By increasing the number of elements, both of these energy quantities decrease. The values of external work and dissipated energy calculated at the end of the numerical simulation for different models with different number of elements are shown in Figure B-8. This figure shows reduction of both of these parameters by increasing the number of elements. In the model with eight elements, there is a difference between the predicted external work and dissipated energy. This shows that the available strain energy in this case has been greater than the capability of the model to dissipate energy. The difference has turned into kinetic energy and resulted in oscillation of the structure as shown in Figure B-4. The quasi-static equilibrium path for the model with eight elements exhibits the snap-back phenomenon. In an explicit dynamic simulation under displacement controlled condition, as performed in this example, the snap-back behaviour cannot be captured. In this situation, the dynamic simulation shows an abrupt jump in the force and oscillation after the peak (Figure B-9). 192 Figure B-1. A one dimensional bar that is partly damaged and partly undamaged supported at one end and subjected to a displacement increment at the other end. Figure B-2. A typical linear strain softening curve showing (a) the initiation and saturation strains and (b) the loading and unloading paths. Lu DamagedUndamaged Displacement increment εs σ εi ε σu σ ε Loading Un loa din g (a) (b) 193 Figure B-3. Predicted Force-Elongation responses for a bar subjected to an axial load (constant velocity at one end). The bar was simulated with five different numbers of elements: 1, 2, 4 and 8 elements. 0 0.4 0.8 1.2 0 0.01 0.02 0.03 0.04 0.05 0.06 Elongation F o rc e 1 Element 2 Elements 4 Elements 8 E le m e n ts 194 Figure B-4. Predicted Force-Elongation response of the bar subjected to axial tension for the model with eight elements. 0 0.4 0.8 1.2 0 0.01 0.02 0.03 0.04 0.05 0.06 Elongation F o rc e 195 Figure B-5. Predicted pattern of the elements that unloaded (solid shade) and damaged (dashed shade) for all the four mesh sizes. 1 Element 2 Elements 4 Elements 8 Elements 196 Figure B-6. Predicted external work versus elongation for the FE models with different number of elements. 0 0.005 0.01 0.015 0.02 0.025 0 0.01 0.02 0.03 0.04 0.05 0.06 Elongation E x te rn al W o rk 1 Element 2 Elements 4 Elements 8 Elements 197 Figure B-7. Predicted dissipated energy versus elongation for the FE models with different number of elements. 0.00 0.01 0.02 0.03 0.00 0.01 0.02 0.03 0.04 0.05 0.06 Elongation D is si p at ed E n er g y 1 Element 2 Elements 4 Elements 8 Elements 198 Figure B-8. Predicted final values of the external work and dissipated energy for the FE models with different number of elements. 0.00 0.01 0.02 0.03 0 2 4 6 8 Number of Elements E n er g y External Work Dissipated Energy 199 Figure B-9. Predicted force-elongation curve for the FE model with 8 elements shown with the solid line compared to the quasi-static equilibrium path shown by dashed line. 0 0.4 0.8 1.2 0 0.005 0.01 0.015 0.02 Elongation F o rc e Dynamic response under displacement controlled condition Equilibrium Path 0.0 06 12 5 200 C. FUNDAMENTALS OF CONTINUUM DAMAGE MECHANICS Continuum damage mechanics aims to simulate the effect of damage and material degradation in a smeared manner by reducing the material’s stiffness. Assuming no permanent deformation in the material, the relation between the nominal stress and strains of a representative volume can be written in the following form: dσ = D (ω)ε (C.1) where σand εare the nominal stress and strain on the RVE and Dd is the stiffness matrix of the damaged material. ω are the parameters that describe the state of damage in the material known as damage parameters. The main challenges in the continuum damage mechanics are firstly writing the D d in terms of damage parameters and secondly writing the evolution laws for the damage parameters. In this Appendix, some of the common approaches used to formulate the damage models are discussed. C.1 REDUCING THE STIFFNESS Various approaches have been proposed in the literature to describe the stiffness reduction of the damaged material. Some of these approaches intend to simulate the effect of damage in originally isotropic materials such as concrete and some others are proposed to simulate the effect of damage in orthotropic materials. Scalar damage model is perhaps the simplest form of continuum damage mechanics models where all the components of the stiffness are reduced by the same amount. In this model, the state of damage is described with only one damage parameter as shown below. 0(1 )σ = D ε (C.2) The scalar damage models are meant to simulate damage in isotropic materials and are unable to predict damage-induced orthotropy. 201 In order to be able to simulate the effect of damage in orthotropic materials and also capture the damage induced orthotropy, the damage model needs to describe the damaged material in terms of a set of damage parameters. In an initially isotropic material damage parameters describe the state of damage in various directions and thus the damage-induced orthotropy can be simulated. In orthotropic materials made from different constituents, where multiple damage mechanisms exist, each damage parameter describes the state of a distinct damage mechanism. In this section the hypotheses of strain equivalence and energy equivalence that are employed by many researchers to formulate the damaged stiffness matrix are discussed in detail. C.1.1 Strain Equivalence Hypothesis The hypothesis of strain equivalence assumes that the strain experienced by the undamaged portion of the material in an RVE is equal to the nominal strain. To demonstrate the method a uniaxial example is presented here. Consider a RVE that is subjected to a nominal stress of σ and has a nominal strain of ε. The relation between the effective strain and stress in the undamaged portion of the material is governed by the stiffness of the virgin material: 0 ˆˆ D = (C.3) The strain equivalence principle assumes that both intact and damaged portions of the material experience the same strain field (i.e. ̂ ). Assuming the strain equivalence principle, the strain in the undamaged material is equal to the nominal strain therefore 0ˆ D (C.4) Assuming that the damage parameter ω represents the ratio of the damaged cross-section of the material to the undamaged cross-section, the nominal stress can be written in terms of the stress in the undamaged material as: ˆ(1 ) (C.5) Substituting (C.5) into (C.4) leads to a relation between the nominal stress and nominal strain: 0(1 )D (C.6) 202 Generalization of this approach to a 3D model for an anisotropic material faces a couple of difficulties. One of them is writing the relation between the nominal shear stress and strains considering the damage parameters. The second problem that is a more fundamental one is that using this approach for a 2D or 3D anisotropic material results in an asymmetric stiffness matrix. Assume that in general, the relation between the nominal stress and the stress that undamaged material experiences is written in following form: ˆσ =M(ω)σ (C.7) where M is a diagonal matrix whose components are written in terms of damage parameters. Equation (C.8) shows the M matrix for a typical 2D condition with three damage parameters. 1 2 12 1 0 0 0 1 0 0 0 1 M(ω) (C.8) Thus the relation between the nominal stress and strain can be written as shown below 0σ =M(ω)D ε (C.9) The stiffness of the damaged material would then be d e D (ω) =M(ω)D (C.10) The damaged stiffness calculated from Equation (C.10) is not symmetrical in general. As an example, in 2D case, assuming three damage parameters (ω1, ω2 and ω12) that correspond to the three components of stress (σ1, σ2 and τ12), the relation between the nominal stress and the stress in the undamaged portion of the material is written as 1 1 1 2 2 2 12 12 12 ˆ1 0 0 ˆ0 1 0 ˆ0 0 1 (C.11) Now under plane stress conditions, the resulting damaged stiffness matrix can be written as: 203 1 1 1 12 2 12 21 12 21 2 21 1 2 2 12 21 12 21 12 12 (1 ) (1 ) 0 1 1 (1 ) (1 ) 0 1 1 0 0 (1 ) d E E E E G D = (C.12) It can be seen that the stiffness matrix shown in Equation (C.12) is not symmetrical. The secant stiffness matrix that relates the total strain to the total stress governs the elastic behaviour of the damaged material and has to comply with the Betty-Ma well’s reciprocity law. The reciprocity law imposes the symmetry condition on the secant stiffness matrix. There have been some remedies suggested by researchers to this problem. For instance, Matzenmiller and co-workers suggested some modification to the strain equivalence approach to accommodate reduction of the Poisson’s ratios (Matzenmiller et al. 1995). Their suggested formulation was backed by experimental data on unidirectional laminates. The approach taken in this thesis to write the damaged stiffness is similar to the approach taken by Matzenmillar et al. except that they employed three damage parameters in their model whereas in this thesis only two damage parameters were introduced. C.1.2 Energy Equivalence Hypothesis The Energy Equivalence hypothesis provides an alternative way to formulate the stiffness matrix of a damaged material. This principle in its dual form assumes that the density of the complementary energy in the damaged material subjected to nominal stress is equal to the complementary energy density in the undamaged portion of the material subjected to effective stress (see (Chow and Wang 1987), (Kennedy and Nahan 1996) and (Jirasek 2007) for instance). This assumption would provide a relation between the macroscopic stress and the effective stress in the intact portion of the material. 01 1 ˆ ˆ 2 2 T d Tσ C σ = σ C σ (C.13) where σ is the nominal stress acting on the RVE and σ̂ represents the stress that is applied to the undamaged portion of the material within the RVE. 204 Similar to the strain equivalence principle, the damage parameters are employed to express the nominal macroscopic stress in terms of the effective stress in the intact material. Substituting (C.7) into Equation (C.13) and simplifying yields 0 T dM (ω)C M(ω) C or 0d TD M(ω)D M (ω) (C-14) Since M is a diagonal matrix, the damaged stiffness matrix derived based on this principle will be automatically symmetric. The author believes that there is a flaw in the basic assumption made for the energy equivalence principle. Assuming the equality of strain (or complementary) energy densities in the intact material within the RVE and the damaged RVE is not a suitable assumption since the volume of the intact material is not equal to the gross volume of the RVE. In other words, in a damaged material not all the gross volume of the RVE contributes to storing the strain energy. Therefore this assumption should be re-stated as the density of the strain energy in the damaged material subjected to nominal strain times the volume of the RVE is equal to the strain energy density in the undamaged portion of the material subjected to the effective strain times the volume of the intact material. C.1.3 Other Approaches Some researchers have taken other paths to derive the stiffness matrix of the damaged material. For example, (Ladeveze 1995) and (Lemaitre et al. 2000) start with writing an expression for the free energy in terms of the damage parameters and then derive the constitutive relationship from it. In another approach, Talreja (Talreja 1985; Talreja 1989) uses the micromechanical assumptions to define the damage parameters and derive the stiffness of the damaged material. 205 C.2 DAMAGE EVOLUTION LAWS Various approaches have been taken to derive the evolution laws for damage in materials. In this section, the methods based on maximum dissipation hypothesis, general dissipation potential and equivalent strain functions are discussed. C.2.1 Maximum Dissipation Hypothesis The maximum dissipation hypothesis is one of the popular approaches taken to derive the evolution laws for damage parameters. This approach that is analogous to the associated flow rule in plasticity, defines the evolution laws based on the work-conjugates of reduction parameters (or damage parameters). Suppose that is the strain energy density of the damaged material that is written in terms of state of strain and reduction parameters as internal state parameters. ( ) ε,R (C.15) The work conjugates of a reduction parameter (also known as damage driving force) is defined as i i Y R (C.16) The rate of energy dissipation density, d , is defined as the difference between the rates of the external work and the recoverable strain energy : :d σ :ε σ :ε ε R ε R (C.17) Knowing that stress is defined as the work conjugate of strain, the rate of dissipated energy density can be simplified as : :d R Y R R (C.18) Based on Equation (C.18), the maximum rate of energy dissipation is achieved if the rate of damage evolution is parallel to the damage driving forces. R Y (C.19) 206 The hypothesis of maximum energy dissipation writes the evolution law based on Equation (C.19) where the rate of evolution of each damage parameter is a function of its own work conjugate. A consequence of this hypothesis is the path-independency of the dissipated energy meaning that the amount of energy dissipated to get from damage state 1ω to 2ω is independent of the load path taken. Although the hypothesis of maximum energy dissipation provides a straight-forward approach to derive an evolution law, it is a purely mathematical hypothesis and does not necessarily agree with the physical reality. C.2.2 Dissipation Potential The concept of dissipation potential provides a general framework for definition of an evolution law. In this approach that is analogous to the non-associated flow rule in plasticity, a dissipation potential function, F d , is defined and the rate of damage evolution is derived based on this function as follows: dF R R (C.20) The dissipation potential function can be defined by fitting a function to the available experimental data. The maximum energy dissipation hypothesis can be seen as a special case of the dissipation potential when the strain energy density is chosen as the dissipation potential function. C.2.3 Equivalent Strain Functions Some researchers have chosen the concept of equivalent strain in defining the evolution law. In this approach evolution law for each damage parameter is written in terms of an equivalent strain function. The equivalent strains are defined as a function of the components of strain. Williams et al. in (Williams et al. 2003) proposed a general expression for the equivalent strain functions. The problem with the equivalent strain functions defined in (Williams et al. 2003) is introduction of numerous model parameters that needs to be calibrated based on the available test data. 207 In the formulation presented in this thesis, a rather simple approach has been taken to define the equivalent strain functions. The equivalent strain function that governs the fibre damage depends on the longitudinal strain in the fibre direction as shown in Equation (5.6) and an interactive expression is chosen to define the equivalent strain function for the matrix damage as shown earlier in Equation (5.7). 208 D. CODAM2 IMPLEMENTATION CODAM2 is implemented as a new material model in OOFEM. In this section, the algorithm, pseudo code and coding overview are presented. D.1 ALGORITHM The algorithm for calculation of stress state given the strain state in an Explicit Dynamic analysis: • Receive the total strain vector • Calculate Equivalent Strain Functions: – Transform the strain to each layer according to its orientation and calculate the equivalent strains for each layer: • Longitudinal strain in the direction of fibres that governs fibre damage propagation (sign of the strain is included. Tensile: positive and Compressive: negative) • The combined transverse and shear strains that governs the matrix damage propagation. (sign of this equivalent strain is determined based on the – Calculate the maximum (absolute) principal strain that is used to calculate the delamination and overall damages • Average the equivalent strains: – Calculate the effective distance between the origin and neighbour points based on Equation (6.1). – Calculate the weight function based on Equation (6.2). – Calculate the averaged Equivalent Strain parameters. • Calculate the Damage Parameters and check damage growth conditions: – In this part of the code the temporary damage parameters are calculated based on the current (temporary) equivalent strain functions 209 – The temporary damage parameters are compared to the damage parameter history and thus damage growth condition is checked. • Calculate the reduction coefficients – The reduction coefficients are calculated as a function of the recently updated damage parameters • Calculate the Damaged Stiffness Matrix (Secant): – Apply the reduction coefficients and calculate the damaged stiffness matrix for each layer. – Rotate and assemble the stiffness of layers – Apply the delamination and overall reductions • Calculate the Stress Vector: – Multiplying the reduced secant stiffness to the total strain vector and calculating the Stress Vector D.2 PSEUDO CODE D.2.1 Equivalent Strains • Equivalent strains used for delamination and overall damage: – Calculating the principal strains (εpr) – Finding the principal strain with the maximum absolute value εprmax – Reporting this as the equivalent strain used for delamination and overall damages – The sign of the strain is included in the reported equivalent strain, i.e. a positive equivalent strain shows governing tensile situation and a negative one is a sign of a governing compressive condition. • Equivalent strain for fibre breakage and matrix damage: – A loop on number of layers (i) 210 • Transforming strains to the local CS of each layer. The essential tools are provided by OOFEM. εTε ii • Calculating the equivalent strain for fibre breakage based on Equation (5.6). • Calculating the equivalent strain for matrix cracking based on Equation (5.7). D.2.2 Non-local Averaging • The non-local averaging is performed over the equivalent strain vector calculated in the previous step. jj eqeq wεε • Averaging is performed on a circle-shaped region. The weight coefficients are calculated based on the orthotropic averaging function proposed in Chapter 6. • OOFEM stores the location and weight of each target integration point in an array attached to each integration point object D.2.3 Calculating the Damage Parameters • Temporary damage parameters are calculated based on temp averaged equivalent strains according to Equation (5.9). • Note that the initiation and saturation strains are different for each of the 4 mechanisms and they also depend on the sign of the equivalent strain (being tensile or compressive). • Damage growth condition: If Damage Growth :temp temp • If the damage growth condition is satisfied, the damage history parameters stored in each Gaussian point will be updated. 211 D.2.4 Calculating the Reduction Coefficients • Reduction coefficients are calculated based on the damage history parameters: 1iR • Where is the maximum reduction of the stiffness term due to damage mechanism α. This parameter is read from the input file. D.2.5 Calculating the Damaged Stiffness Matrix • Matrix and fibre reduction coefficients are applied to the stiffness of each layer and thus the damaged stiffness of the layer is calculated: 1 12 2 12 21 12 21 2 12 21 12 0 1 1 0 1 f f m f m f m d m k f m m k R E R R E R R R R R E R R SYM R G Q • The contributions of layers in the stiffness of the laminate are added up and the reductions due to delamination and overall damage are applied at the laminate level: 1 n d T d k k k k k t A Τ Q Τ D.2.6 Calculating the Stress Vector • The reduced in-plane secant stiffness of the laminate is multiplied by the total strain vector to calculate the stress: 1 1 d n kt σ A ε 212 D.3 IMPLEMENTATION To implement the CODAM2 in OOFEM platform, a couple of new classes have been implemented as shown in the code hierarchy (Figure D-2). D.3.1 OrthoDamage2 OrthoDamage2 is the parent class that includes the common functions that a typical orthotopic damage model would need. All the state variables and history parameters (including Equivalent Strains, Damage Parameters and Reduction Coefficients) are declared in the status class linked to OrthoDamage2. By instantiating the OrthoDamage2, an OrthotropicLinearElasticMaterial is instantiated that provides the calculation for elastic stiffness as well as strain transformation (rotation) calculations. Calculation of the damaged stiffness matrix based on the calculated reduction coefficients is performed in this class. D.3.2 CODAM2 The formulation of the CODAM2 is implemented in the CODAM2 class. This class is called when the local version of CODAM2 material is used in the input file. Elastic material properties and number and orientation of the plies is read from the input line. The damage parameters such as initiation and saturation strains are also read form the parameters file. Calculation of equivalent strains, damage parameters and reduction coefficients is performed in this class. So basically the details and formulations of CODAM2 are in this class. D.3.3 StructuralOrthoNonLocalMaterialExtensionInterface This class is responsible for providing the essential functionalities for orthotropic averaging. Calculation of weight functions for orthotropic averaging is defined and performed in this class. 213 D.3.4 CODAM2NL This class is instantiated when CODAM2NL material is called in the input file. The non-local radii are read from the input file and used in this class. The CODAM2 class is called to perform the functionalities that are common between the local and non-local implementations. CODAM2NL is responsible for calculation of averaged equivalent strain. In this class the computeEquivalentStrain is overwritten on the CODAM2’s function. This function in CODAM2NL class performs the averaging on the equivalent strains. To perform the non-local averaging, the code loops over all the integration points and asks them to calculate the local equivalent strain based on the given local strain. When the local equivalent strains are calculated and stored, the loop over the integration points restart and this time averaging is performed over the neighbourhood and non-local strain is calculated. 214 Figure D-1. Flowchart of CODAM2 calculations receiving total strain vector Calculating the equivalent strains Averaging on equivalent strains Calculating the temp damage parameters Checking the damage growth condition Calculating the reduction coefficients Calculating the reduced secant stiffness Calculating and passing total stress vector to the code 215 Figure D-2. CODAM2 and CODAM2NL class hierarchy Structural Material OrthoDamage2 CODAM2 StrucruralNonLocalMaterial ExtensionInterface CODAM2NL OrthotropicLinear ElasticMaterial StrucruralOrthoNonLocalMa terialExtensionInterface
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A non-local approach to simulation of damage in composite...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A non-local approach to simulation of damage in composite structures Forghani, Alireza 2011
pdf
Page Metadata
Item Metadata
Title | A non-local approach to simulation of damage in composite structures |
Creator |
Forghani, Alireza |
Publisher | University of British Columbia |
Date Issued | 2011 |
Description | Fibre reinforced polymers (FRP) are one of the fastest growing materials in advanced structural applications. Due to the complexity and existence of multiple and interactive modes of failure, there is a lack of comprehensive theory that describes the damage behaviour of these materials. The UBC Composite Damage Model (CODAM) is a sub-laminate based model that is designed to simulate the behaviour of laminated composites at the macro (structural) scale. Physical basis and numerical robustness are the main objectives of CODAM development. The original formulation of CODAM that was developed in the mid 1990’s suffers from material objectivity problem that results in spurious dependency of numerical predictions on the choice of the coordinate system. In this thesis, the objectivity issue of the original CODAM is identified and addressed through a new formulation called CODAM2. The new formulation is capable of predicting damage in various laminate lay-ups from quasi-isotropic to unidirectional. It is also capable of simulating the damage-induced orthotropy in the laminate. An orthotropic non-local averaging scheme is developed for CODAM2 to address the localization issue. Compared to isotropic averaging, the orthotropic scheme significantly improves the predicted crack direction and damage pattern in composite laminates. Unlike the conventional isotropic non-local averaging that performs the averaging over a circular (spherical in 3D) zone, in the orthotropic scheme the averaging is performed over an elliptical zone which requires the introduction of two length scales. The performance of CODAM2 equipped with orthotropic averaging is demonstrated through numerical examples. It is shown that the developed model is capable of accurately predicting the damage behaviour in various specimen geometries from sharp-notched to blunt-notched and open-hole specimens. The predictions of the model in terms of load-displacement, crack-tip position, damage height and crack directions agree well with experimental observations and measurements. CODAM2 provides a promising numerical tool to simulate the effect of damage on the behaviour of structures made of laminated composites. This model is computationally efficient and yet relatively simple to understand, calibrate and use in practical applications. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-10-21 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0050751 |
URI | http://hdl.handle.net/2429/38154 |
Degree |
Doctor of Philosophy - PhD |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2011-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2011_fall_forghani_alireza.pdf [ 3.97MB ]
- Metadata
- JSON: 24-1.0050751.json
- JSON-LD: 24-1.0050751-ld.json
- RDF/XML (Pretty): 24-1.0050751-rdf.xml
- RDF/JSON: 24-1.0050751-rdf.json
- Turtle: 24-1.0050751-turtle.txt
- N-Triples: 24-1.0050751-rdf-ntriples.txt
- Original Record: 24-1.0050751-source.json
- Full Text
- 24-1.0050751-fulltext.txt
- Citation
- 24-1.0050751.ris
Full Text
Cite
Citation Scheme:
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>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0050751/manifest