Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Nonlinear XFEM modeling of delamination in fiber reinforced composites considering uncertain fracture… Motamedi, Damoon 2013

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

Item Metadata

Download

Media
24-ubc_2013_fall_damoon_motamedi.pdf [ 4.35MB ]
Metadata
JSON: 24-1.0073939.json
JSON-LD: 24-1.0073939-ld.json
RDF/XML (Pretty): 24-1.0073939-rdf.xml
RDF/JSON: 24-1.0073939-rdf.json
Turtle: 24-1.0073939-turtle.txt
N-Triples: 24-1.0073939-rdf-ntriples.txt
Original Record: 24-1.0073939-source.json
Full Text
24-1.0073939-fulltext.txt
Citation
24-1.0073939.ris

Full Text

NONLINEAR XFEM MODELING OF DELAMINATION IN FIBER REINFORCED COMPOSITES CONSIDERING UNCERTAIN FRACTURE PROPERTIES AND EFFECT OF FIBER BRIDGING by Damoon Motamedi  B.Sc., The University of Tehran, 2006 M.Sc., The University of Tehran, 2008  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF  DOCTOR OF PHILOSOPHY in THE COLLEGE OF GRADUATE STUDIES (Civil Engineering)  THE UNIVERSITY OF BRITISH COLUMBIA (Okanagan)  June 2013  © Damoon Motamedi, 2013  Abstract Initiation and propagation of a crack in composite materials can affect their global mechanical properties severely. From a numerical modeling perspective, most conventional macro-level methods reported for composite laminates are based on the assumption that a Representative Volume Element (RVE) of the material is periodically repeated over the entire sample. However, a considerable amount of spatial non-uniformity in material and geometrical parameters can exist in both unidirectional (UD) and woven fabric composites. The scattered distribution of fibers, fibers penetration between composite layers, voids within the matrix, human errors during sample preparation, and imperfect thickness distribution can be among the most common sources of such non-uniformity. In turn, these non-uniformities can make the numerical simulation of composites under the assumption of a periodic RVE unreliable, and thereby, the stochastic modeling of effective material properties becomes essential for a more precise assessment of composites’ mechanical behaviour. In the present work, a new three-dimensional (3D) stochastic extended finite element method (XFEM) is proposed and implemented to model the delamination surface in composite samples by integrating the capabilities of the finite element method (FEM) commercial software (ABAQUS) into a user-defined FORTRAN code and MATLAB package. XFEM is known to offer significant advantages over conventional FEM by enabling optimal convergence rates in the presence of pronounced discontinuities/singularities such as cracks. The effect of nonlinear modeling parameters such as cohesive zone length, penalty stiffness factor and large deformation are also considered in the proposed approach to add to the accuracy of simulations. The XFEM model is first tested and validated against previously reported data in the literature. Next, a statistical distribution is sought from data nonii  repeatability during a set of double cantilever beam (DCB) and end-notched flexure (ENF) tests conducted on Poly (phenylene Sulfide) PPS/Glass thermoplastic composite samples. Results from the experiments and XFEM are compared and demonstrate the capability of the new numerical approach in capturing non-repeatable material response, often seen during the fracture testing of UD composites to characterize their mode I and mode II fracture properties.  iii  Table of Conttents Absttract .......................................................................................................................................... ii  Tablle of Conten nts .......................................................................................................................... iv  List of o Tables ............................................................................................................................... viii  List of o Figures ............................................................................................................................... ix  List of o Symbols ............................................................................................................................. xv  Ackn nowledgemeents .................................................................................................................... xxv  Dediication .................................................................................................................................. xxvi  1 Chapter: C Inttroduction ............................................................................................................... 1  1.1   Fibers ............................................................................................................................... 1   1.2   Matrix .............................................................................................................................. 2   1.3   Classiffication of FR RP Composiite ................................................................................... 4   1.4   FRP Co omposite Materials App plications ......................................................................... 5   1.5   FRP Co omposite Materials Weaaknesses andd Applicationn Limitationns ......................... 6   1.6   Experim mental and Numerical N Modeling M of F FRP Fracturre Propertiess......................... 10   1.7   Random mness in FR RP Fracture Properties P ...................................................................... 14   1.8   Motivaation and Ob bjectives of the Work ....................................................................... 15   1.8.1   Potentiial for Practiical Applicattions and Exxpected Origginality .................................. 17   1.9   Thesis Outline ................................................................................................................ 19   iv  2 Chapter: C Baackground ............................................................................................................. 21  2.1   Elastic Mechanicall Behaviour of FRP com mposite Mateerials ...................................... 21   2.2   Failuree Modes of FRP F Compon nents ............................................................................. 23   2.3   Fracturre Mechanics........................................................................................................ 26   2.4   Damag ge Mechanics Models........................................................................................... 28   2.5   Compaarison between Fracture Mechanics aand Damagee Mechanics ......................... 31   2.6   Introdu uction of Exttended Finitee Element M Method ....................................................... 32   2.7   Application of XFE EM in Lineaar Elastic Fraacture Mechhanics ..................................... 40   2.8   Application of XFE EM in Elasto o-Plastic Fraacture Mechaanics ..................................... 44   2.9   Application of Coh hesive Interfface in Fractuure Mechaniics ......................................... 47   2.9.1   Cohesiive Zone Mo odel ................................................................................................... 48   2.10   Summaary ........................................................................................................................ 54   3 Chapter: C Extended Finitte Element Method M Impllementation and Validatiion..................... 55  3.1   Large Deformation D n Formulatio on ................................................................................... 56   3.2   Nonlin near Solvers .......................................................................................................... 60   3.3   Modeliing Contact on Material Interfaces U Using XFEM M ............................................ 63   3.4   Implem mentation off the Cohesiv ve Zone Moddel ............................................................. 66   3.5   Other Numerical N Im mplementatiion Details ..................................................................... 69   3.5.1   Node Selection S forr Enrichmen nt .................................................................................... 69   v  3.5.2   Numerrical Integrattion of Disco ontinuous Fiields .......................................................... 71   3.5.3   Implem mentation off the Integrattion Bound A Approach .................................................. 74   3.5.4   Simulaation Algoritthm.................................................................................................... 75   3.6   Numerrical Examplles of Modee I and II Frracture Testss: Validationn of XFEM Code ............................................................................................................................... 78   3.6.1   Numerrical Simulattion of the DCB D Tests ..................................................................... 78   3.6.1.1   Effectss of Differen nt Modeling Approachess ................................................................. 79   3.6.1.2   Effectss of Mesh Siize and Coheesive Zone L Length ....................................................... 82   3.6.1.3   Effectss of Differen nt Penalty Sttiffness Factoors ............................................................. 85   3.6.2   Numerrical Simulattion of the ENF E Tests ...................................................................... 92   3.6.2.1   Effectss of Cohesiv ve Zone Leng gth ................................................................................. 92   3.6.2.2   Effectss of Differen nt Penalty Sttiffness Factoors ............................................................. 94   3.7   Summaary ........................................................................................................................ 96   4  Chapter: C Modeling M Ran ndomness Effect E in UD D Laminates Delamination: A Non-  C  RVE Approach A .............................................................................................................. 97   4.1   Samplee Preparation n: Poly (phen nylene Sulfiide) (PPS)/G Glass FRP ............................ 100   4.2   Elastic Mechanicall Properties of o PPS/Glas s FRP Compposites................................. 103   4.3   Fracturre Tests on th he Fabricateed PPS/Glasss Compositees ......................................... 105   4.3.1   DCB and a ENF Tesst Results......................................................................................... 106   4.4   Stochastic Fracturee Properties .................... . ................................................................ 112  vi  4.5   Numerical Results and Discussions ........................................................................ 115   4.6   Summary ................................................................................................................. 125   5  Chapter: Conclusions and Future Work Recommendation ............................................ 126   5.1   XFEM Model Development ................................................................................... 126   5.2   Performed Deterministic Simulations .................................................................... 126   5.3   Performed Stochastic Simulations .......................................................................... 127   5.4   Potential Future Work............................................................................................. 128   References.............................................................................................................................. 130  Appendices ............................................................................................................................ 143  Appendix A: ABAQUS User-element Subroutine for Nonlinear XFEM Analysis .............. 143  Appendix B: Experimental Calculations According to ASTM D5528-01 [97] .................... 157   vii  List of Tables Table 1-1 Mechanical properties of fibers in commonly used FRPs [3] .................................. 2  Table 1-2 Mechanical properties of matrices in commonly used FRP composites [3] ............ 4  Table 3-1 Mechanical properties of T300/977-2 and AS4/PEEK samples [52] .................... 79  Table 4-1 Mechanical properties extracted from tensile testing ........................................... 104  Table 4-2 Final set of elastic properties of manufactured UD PPS/Glass FRP composites CCCCCC with “1” being the fibers direction; “2” and “3” are perpendicular directions to CCCCCC fibers. .................................................................................................................... 105  Table 4-3 Fracture properties of PPS/Glass samples extracted from DCB and ENF tests... 112  Table 4-4 Employed standard deviation schemes in stochastic simulations of DCB test .... 115   viii  List of Figures Figure 1-1 Sample applications of FRP composite materials in different industries: (a) & CCCCCC  (b) Airbus Military A400M aircraft CFRP wing (Martin Chainey / Airbus  CCCCCC  Military), (c) Cross-sections of glass fiber reinforced polymer (GFRP)  CCCCCC  structural members used, and (d) Pontresina Bridge made of FRP [5] .................. 6   Figure 1-2 FRP composite materials common failure modes: (a) Fiber breakage, (b) CCCCCC  Fiber pull-out, (c) Matrix cracking, and (d) Interlaminar delamination ................. 8   Figure 1-3 Different scales of material modeling during homogenization of material CCCCCC  properties .............................................................................................................. 10   Figure 1-4 Fracture energy toughness measurement fixtures: (a) Double Cantilever CCCCCC  Beam, (b) End-Notched Flexure, (c) Edge Crack Torsion, and (d) Mixed  CCCCCC  Mode Bending tests [8]......................................................................................... 12   Figure 1-5 Forming process of a heated FRP composite laminate: (a) Initiation of the CCCCCC  forming (Step 1), compression and cooling cycle (Step 2), removing the  CCCCCC  male die and possible dimensional distortion of the part (Step 3); and (b)  CCCCCC  Delamination spotted in the corner of a compression moulded component [54] . 18   Figure 2-1 Comparison of tensile stress-strain curve of fiber materials [58] ......................... 24  Figure 2-2 Stress-strain curves for the matrix and fiber materials: (a) Fibers have larger CCCCCC  failure strain than that the matrix, and (b) The matrix has larger failure strain  CCCCCC  than fiber [58] ....................................................................................................... 25   Figure 2-3 Dominant fracture modes of a cracked body: (a) Opening mode (mode I), (b) CCCCCC  Sliding (shearing) mode (mode II), and (c) Tearing mode (mode III) [59] .......... 26  ix  Figure 2-4 An arbitrary finite element mesh with a discontinuity (circles represent the CCCCCC  enriched nodes of the mesh) ................................................................................. 33   Figure 2-5 Unit tangential and normal vectors for the Heaviside function and nearest CCCCCC  point to X on the crack surface; X*........................................................................ 34   Figure 2-6 The influence domain of node J in an arbitrary finite element mesh..................... 38  Figure 2-7 Local crack-tip coordinates and the contour Γ and its interior area, VΓ................ 42  Figure 2-8 Schematic of the cohesive zone in front of crack in a given step of numerical CCCCCC  simulation ............................................................................................................. 52   Figure 3-1 Explicit solver approach and possible drift error in nonlinear problems (dots CCCCCC  show numerical solution steps) [96] ..................................................................... 61   Figure 3-2 Newton-Raphson iterative solver approach in nonlinear FEM problems [96] ..... 62  Figure 3-3 Bilinear traction-separation law for modeling the material degradation .............. 67  Figure 3-4 Variation of φ1 values in an example meshed object (square represents the CCCCCC  positive value and circle denotes the negative value; the rectangle shows the  CCCCCC  crack plane)........................................................................................................... 70   Figure 3-5 Variation of ψ1 values for the front edge of the crack plane in the example CCCCCC  meshed object (square represents the positive value and circle denotes the  CCCCCC  negative value) ...................................................................................................... 71   Figure 3-6 Accepted nodes in the example meshed object based on Equations (3-39) and CCCCCC  (3-41) of level-set criteria ..................................................................................... 71   Figure 3-7 Sub-triangles of a 2D element with third order Gauss quadrature........................ 72  x  Figure 3-8 Sub-tetrahedrals of a 3D element with third order Gauss quadrature................... 73  Figure 3-9 Integration points of a meshed object with enriched elements (red stars are CCCCCC  integration points and blue dots are integration points within the integration  CCCCCC  bound) ................................................................................................................... 75   Figure 3-10 MATLAB-ABAQUS simulation algorithm employed for modeling the CCCCCC  delamination ......................................................................................................... 77   Figure 3-11 A comparison between DCB test results via different methods on T300/977CCCCCC  2 samples .............................................................................................................. 81   Figure 3-12 Load-Displacement DCB test results for the fine mesh (le = 0.4 mm) CCCCCC  simulation with different cohesive zone lengths for T300/977-2 samples ........... 83   Figure 3-13 Load-Displacement DCB test results for the coarse mesh (le = 1.25 mm) CCCCCC  simulation with different cohesive zone lengths for T300/977-2 samples ........... 84   Figure 3-14 The comparison between DCB test load-displacement results of T300/977-2 CCCCCC  samples with different penalty stiffnesses ............................................................ 87   Figure 3-15 The comparison between DCB test load-displacement results of AS4/PEEK CCCCCC  samples with different penalty stiffnesses ............................................................ 88   Figure 3-16 The comparison between DCB test load-displacement results for fine mesh CCCCCC  analysis of AS4/PEEK and previous works ......................................................... 90   Figure 3-17 The comparison between DCB test load-displacement results for coarse CCCCCC  mesh analysis of AS4/PEEK and previous works ................................................ 91   xi  Figure 3-18 The comparison between ENF test load-displacement results for AS4/PEEK CCCCCC  and previous works ............................................................................................... 93   Figure 3-19 Effect of penalty stiffness value on the ENF test load-displacement results CCCCCC  for AS4/PEEK sample .......................................................................................... 95   Figure 4-1 Microscopic images of fibers and matrix distribution of PPS/Glass UD tape: CCCCCC  (a) Corner of the tape, and (b) Middle of the tape [105] .................................... 101   Figure 4-2 (a) Forming cycle used for preparing PPS/Glass test samples using (b) an CCCCCC  automated press apparatus (Wabash MPI 100 ton) ............................................ 103   Figure 4-3 Snapshot of elastic mechanical properties for a typical woven PPS/Glass ply CCCCCC  from material data sheets .................................................................................... 105   Figure 4-4 Experimental test set-ups: (a) DCB, and (b) ENF .............................................. 106  Figure 4-5 The variation of fracture energy toughness versus crack length for the three CCCCCC  tested samples using: (a) Compliance Calibration Method, (b) Modified  CCCCCC  Beam Theory Method, (c) Modified Compliance Calibration Method [97] ...... 108   Figure 4-6 Illustration of the fiber bridging zone (FBZ) during crack propagation; as the CCCCCC  crack length increases, FBZ emerges in the cracked region up to the fiber’s  CCCCCC  rupturing displacement; after fiber breakage, the FBZ effect vanishes from  CCCCCC  the region which has exceeded the failure opening displacement ...................... 109   Figure 4-7 Different images of a DCB test sample: (a) macro scale image of fiber CCCCCC  bridging, (b) X-ray micro-tomography image of fiber bridging along the  CCCCCC  sample thickness, and (c) attenuation of the X-ray reflection due to  CCCCCCC absorption; demonstrating uneven distribution of fibers .................................... 110  xii  Figure 4-8 ENF test repeat results with a constant crack length (43 mm): (a) the CCCCCC  variation of fracture energy toughness versus the mid-span displacement, and  CCCCCC  (b) the variation of bending load versus the mid-span displacement. ................ 111   Figure 4-9 Proposed stochastic bilinear traction-separation behaviour (Rand2 is a CCCCCC  random number taken from a 2-parameter Weibull distribution; GCL and GCH  CCCCCC  correspond to the lower and upper limits of GC via Equation 4-1) .................... 113   Figure 4-10 Comparison of the opening force in stochastic simulations of DCB tests CCCCCC  with experimental data using: (a) constant standard deviation formulation,  CCCCCC  and (b) standard deviation as a function of crack length .................................... 116   Figure 4-11 Comparison of predicted fracture energy toughness via stochastic simulations CCCCCC of DCB tests with experimental data using: (a) constant standard deviation CCCCCC  formulation, and (b) standard deviation as a function of crack length ............... 116   Figure 4-12 Comparison of measured opening force with predicted values in stochastic CCCCCCC and deterministic simulations of continuous DCB test using: (a) fracture CCCCCCC energy/toughness remains equal to the average value of experiments, (b) CCCCCCC fracture toughness only changes with increase in delamination length, (c) CCCCCCC fracture toughness increases with extension of delamination with constant CCCCCCC standard deviation formulation, and (d) fracture toughness increases with CCCCCCC extension of delamination with standard deviation as a function of CCCCCCC delamination length ............................................................................................. 117  Figure 4-13 Evolution of the cohesive zone in front of crack upon loading in a given CCCCCC  simulation step .................................................................................................... 120   xiii  Figure 4-14 Stages of delamination propagation within the DCB numerical model: (a) CCCCCCC Onset of rigid hardening in the process zone, (b) Apex of the bilinear CCCCCCC traction-separation law, and (c) Deterioration of the cohesive stiffness ............. 121  Figure 4-15 Comparison of stochastic measured and predicted force-displacement CCCCCC values in ENF tests on the PPS/Glass samples ................................................... 122  Figure 4-16 Stages of delamination propagation within the ENF numerical model: (a) CCCCCC  Apex of rigid hardening in the process zone, (b) Initial stage of crack  CCCCCC  propagation, and (c) Extensive deterioration of material ................................... 124   xiv  List of Symbols Al arg e  Critical energy release rate large displacement correction factor  a0  Initial delamination/crack length  acr  Crack/delamination length in the current configuration  acr  Crack/delamination increment  aij  Material compliance matrix components  B and Bij  Cartesian derivatives of shape functions matrix and components  B  Transformed Cartesian matrix of derivatives  BCoh  Cohesive/contact matrix of shape functions  bH  Crack surface (Heaviside) enriched degree of freedom  C and Cijkl  Material stiffness matrix and its components  Cij  Plane strain/stress material stiffness matrix components  c tip  Crack-tip enriched degrees of freedom  c ij  Imaginary constants of anisotropic material  Dmg  Damage index  DInterface  Matrix of interface material properties  DSep  Elastic-plastic constitutive matrix  Dij  Imaginary compliance components of anisotropic material  di  Arm of loading on the test sample  EGreen  Green strain tensor xv  Ei  Elastic moduli in different directions  EL  Linear part of Green strain tensor  Elongitudinal  Longitudinal Young’s modulus of anisotropic material  ENL  Non-linear part of Green strain tensor  Etransverse  Transverse Young’s modulus of anisotropic material  en  Unit normal vector of the crack alignment  es  Unit tangential vector of the crack alignment  ez  Unit binormal vector of the crack alignment  F  Resisting force against the deformation of numerical model  F  Deformation gradient in the reference coordinate system  F enr  Crack-tip enrichment functions  Fi  Second order strength tensors  Fij  Fourth order strength tensors  f1 ( X , Y , Z )  Arbitrary continuous function  f1' ( X , Y , Z )  Arbitrary discontinuous function  fb  Vector of body forces  ft  Vector of external forces  G  Total energy release rate  G I , G II and G III  Fracture energy toughness for different fracture modes I, II and III  GIC  Mode I critical fracture energy toughness xvi  GIIC  Mode II critical fracture energy toughness  GC  Total/general critical fracture energy release rate  GC (ave )  Mean critical fracture energy release rate  GC ( std )  Standard deviation of critical fracture energy release rate  GCH  Highest critical fracture energy release rate  GCL  Lowest critical fracture energy release rate  GS  Matrix of shape function derivatives  Gij  Shear moduli in different directions  g i ( )  Component of anisotropic crack-tip enrichment functions  Fext  External load acting on a body  Fint  Internal force of a body  H (x)  Heaviside step function at point x  I 33  Identity matrix  J  J-integral  J Auxiliary  Auxiliary J-integral  J Superimposed  Superimposed J-integral  K0  Initial stiffness matrix  K I , K II and K III  Stress intensity factors  KCoh  Cohesive/contact portion of tangential stiffness matrix  KGeo  Geometrical portion of tangential stiffness matrix xvii  KMat  Material portion of tangential stiffness matrix  KT  Tangential stiffness matrix  K Pen and K ii  Penalty stiffness variable and its matrix form for three dimensional problems  K ij  Stiffness matrix components for linear finite element method  aux K Iaux , K IIaux and K III  Auxiliary stress intensity factors  2 K 1node and K node  Sets of nodes associated with crack-tip regions  L  Test sample length  L(c) and Lij (c)  Orthotropic material fracture properties constant  lch  Cohesive zone length  le  Finite element mesh size and element length  M and Mi  Interaction integral  M Local  Local interaction integral  MS  Re-arranged second Piola-Kirchhoff stress tensor  mLine  Slope of ordinary least square technique line  m poly  Number of sub-polygons in an element for numerical integration  mij  Direction cosines  N ( x)  Conventional finite element method shape function  Nall  Finite element mesh nodes  N enr  Normalized enriched shape function  xviii  Nf  Enriched nodes of the finite element mesh  g N node  Set of nodes corresponding to the crack face  ni  Vector of outward normal direction  nG  Gauss quadrature order  t  t  ( n1 , n2 , n3t )  Unit vector of the crack edge  P and Pij  Nominal stress tensor and its components  PContact  Vector of contact forces  Pgrip  Loading on a the test sample (from the tensile machine crosshead)  pi q'  qi  Anisotropic material compliance constant depending on longitudinal material properties Smoothing function of contour integral Anisotropic material compliance constant depending on transverse material properties  R0  Right hand side vector of external forces  RCOD  Ratio of opening to sliding displacements  RSIF  Ratio of dynamic stress intensity factors  Rs (u h )  Residual forces  Randi  Randomly chosen number from a Weilbull distribution  (r, θ)  Local polar coordinate system at the crack-tip  S  Material compliance matrix  xix  Sij  Shear strengths in different material directions  T  Traction of interfacial material  Tmax  Maximum interfacial strength  Tij  Interlaminar material strength in different directions  tBlock  Loading block thickness  tck  Test sample thickness  ti  Surface traction vector  u enr  Enriched degree of freedom  uk  Displacement field of the integration domain  du k  Displacement increment  uh  General vector containing all the nodal parameters  u ord  Ordinary degrees of freedom  u g ( x)  Displacement field at point x  ui  Displacement field in different directions  ukAuxilairy  Auxiliary displacement field in different directions  tip u xtip , u tip y and u z  Displacement fields near the crack-tip in different directions  Ve  Element volume  V*  Ratio of each sub-polygon volume to element volume  vij  Poisson’s ratios in different directions  W  Strain energy density  xx  W aux  Auxiliary strain energy density  w  Test sample width  wg  Gauss points weights  wp  Modification factor for the weight of each sub-polygon  XC  Compressive strength in fibers direction  Xi  Reference (Lagrangian) coordinates  XT  Tensile strength in fibers’ direction  X*  Point on the crack surface  X’  Arbitrary point in finite element mesh  (X, Y, Z)  Reference coordinates  ( X 1 , Y1 , Z1 )  Arbitrary point 1 coordinates on the crack surface  ( X 2 , Y2 , Z 2 )  Arbitrary point 2 coordinate on the crack edge  *  *  *  ( X 1 , Y1 , Z1 )  Nodal coordinates in the reference configuration  (x, y, z)  Current coordinates  YC  Compressive strength in perpendicular direction to fibers  YT  Tensile strength in perpendicular direction to fibers  YSnormal  Normalized failure threshold  ZC  Compressive strength in perpendicular direction to fibers  ZT  Tensile strength in perpendicular direction to fibers    Penalty stiffness scaling factor  xxi  1  Weibull distribution shape parameter  1  Weibull distribution scale parameter    Far-field surface path  c  Crack surface path  t  Traction surfaces  V  Body domain    Opening/Sliding displacement of the grip load    Relative opening/sliding displacement  0  Relative critical opening/sliding displacement   eff  Effective separation displacement  f  Relative failure opening/sliding displacement   fH  Highest relative failure opening/sliding displacement   fL  Lowest relative failure opening/sliding displacement  n  Normal crack-tip opening displacement   n0  Relative normal crack-tip opening displacement parameter  t  Tangential crack-tip opening displacement  t0  Relative tangential crack-tip opening displacement parameter    Relative crack displacement in the global coordinate system   and  kl  Green strain tensor and its components  xxii   Auxilairy  Auxiliary strain tensor components  εf,ult  Fiber failure strain  εm,ult  Matrix failure strain  ε’f  Fiber strain at matrix failure  ε’m  Matrix strain at fiber failure  0  Crack-tip angle  i  Transformed angles in anisotropic material  ij  i , i ,  xi and  yi  Anisotropic material properties derived from material governing equation  *  Correction factor of relative crack displacements ratio  ( 1 ,  2 ,  3 )  Local coordinates of finite element method  (u k )  Virtual work   and  ij  Second Piola Kirchhoff stress tensor and its components  i  Re-arranged second Piola Kirchhoff stress tensor   eff  Effective traction  σf,ult  Fiber failure stress  σm,ult  Matrix failure stress    Cauchy stress tensor  σ’f  Fiber stress at matrix failure mode  σ’m  Matrix stress at fiber failure mode  xxiii   ijAuxilairy  Auxiliary stress tensor   ( x)  Conventional shape function    n ,  t   Traction-separation model potential function  0  Material fracture energy  1  Element nodal distance value from the crack face  ,   Orthotropic material parameters relating the stress intensity factors to crack opening/sliding   ( x)  General enrichment function  1  Element nodal distance from the edges of the crack plane  e  Element domain  xxiv  Acknowledgements I wish to thank Dr. Abbas S. Milani for his great supervision, mentorship, inspiration, help and constant encouragement during my PhD research and studies at the University of British Columbia. I would also like to acknowledge the valuable assistance of Drs. Martin Bureau, Francis Thibault, David Butcher, Hicham Mir, and Zohir Benrabah from the National Research Council of Canada – Industrial Materials Institute (NRC-IMI) for training and providing access to their lab facilities. Finally, financial support from the Natural Sciences and Engineering Research Council (NSERC) of Canada is greatly acknowledged.  xxv  Dedication To my parents and family for their unconditional love and support.  xxvi  1 Chapter: Introduction Composites are generally defined as multiphase materials made by commingling two or more existing materials to achieve required mechanical, thermal or electrical properties. Fiber reinforced polymer (FRP) can be categorized as a composite material that is composed of a base polymer material called matrix and reinforced by rebars called fibers. A complete bonding between these two material phases can provide sufficient functionality toward application of composites. The effective mechanical properties of FRP composites are of great interest to aerospace and structural engineers and are directly related to the properties of individual components of each given FRP as reviewed below.  1.1  Fibers  Fibers are materials made into long filaments with 10μm diameter. The main duties of fibers in composites consist of carrying the external loads and providing the required stiffness, strength and thermal stability. Fibers should demonstrate the following characteristics to satisfy the desirable mechanical properties of FRP composites [1]: 1- High modulus of elasticity. 2- High ultimate strength. 3- Low variation of strength among their distributions. 4- High stability of their strength during handling. 5- High uniformity of diameter and surface dimension among bundles of fibers.  1  In the present industrial applications, there are mainly three dominating types of fibers: carbon, glass and aramid fibers. Carbon fibers have high elastic modulus and fatigue failure strength in comparison to the other two types of fibers; however, the economic aspects of their application should always be considered [2]. Glass fibers have demonstrated good mechanical strength relative to their low cost, which has increased their application as reinforcing agent in FRP composites in cost-sensitive industries such as construction [1]. In Table 1-1, a summary of mechanical properties of commonly used fibers is presented.  Table 1-1 Mechanical properties of fibers in commonly used FRPs [3]  Longitudinal Density Material  (kg/m3)  Longitudinal  Relative  Strength (MPa)  cost  Modulus of Elasticity (MPa)  E-glass  2540  70000  3450  Low  S-glass  2500  86000  4500  Moderate  Graphite, high modulus  1900  400000  1800  High  Graphite, high strength  1700  240000  2600  High  Kevlar 29  1450  80000  2800  Moderate  Kevlar 49  1450  130000  2800  Moderate  1.2  Matrix  The matrix is the base material of a FRP composite which confines and bonds the fibers together. It is mainly composed of plastic ingredients that are less environmentally hazardous  2  and more corrosion resistant in comparison to fibers. The followings are the main duties and features of matrix materials in FRP composites [1]: 1- Transferring the external loads to the fibers while binding them together. 2- Protecting the fibers from environmental hazards and corrosion. 3- Providing the general shape of a given FRP structure. 4- Defining some effective mechanical properties of the composite (mainly perpendicular to the fiber directions) such as ductility, impact strength, etc. In industrial applications, FRP composite materials are often manufactured with two types of matrices: thermosets and thermoplastics. Thermosets provide higher strength in comparison to thermoplastics; however, once cured, they cannot be reheated and reused or reformed. This is in contrast to thermoplastic matrices which can be re-employed for recycling through melting and solidification cycles. In Table 1-2, mechanical properties of common matrix materials are presented.  3  Table 1-2 Mechanical properties of matrices in commonly used FRP composites [3]  Density  Tensile Modulus  Tensile Strength  (kg/m3)  (MPa)  (MPa)  Epoxy  1200-1400  2500-5000  50-110  Phenolic  1200-1400  2700-4100  35-60  Polyester  1100-1400  1600-4100  35-95  Nylon  1100  1300-3500  55-90  PEEK  1300-1350  3500-4400  100  PPS  1300-1400  3400  80  Polycarbonate  1200  2100-3500  55-70  Polyethylene  900-1000  700-1400  20-35  Teflon  2100-2300  -  10-35  Material  1.3  Classification of FRP Composite  Numerous methods for classification of FRP composites exist in the literature. For instance, the classification can be based on the different matrices: thermosets and thermoplastics. Another approach for classifying the FRP composites is based on fiber material type: carbon, glass or aramid. Fiber length can also be applied to differentiate the FRP composites: shortfiber and long-fiber composites. Short-fiber reinforced composites consist of a dispersed phase of discontinuous fibers (fiber length is less than hundred times diameter) with random or preferred orientation of fibers. On the other hand, long-fiber reinforced composites contain continuous fibers bonded together by the matrix. They can be composed of a unidirectional  4  orientation of fibers or a bidirectional orientation of fibers (e.g., woven composites), or 3D textile preforms.  1.4  FRP Composite Materials Applications  Today fiber Reinforced Polymer (FRP) composite materials are widely used in high tech engineering applications including aeronautical, marine and automotive industries. These materials have high strength-to-weight ratios, good corrosion resistance, superior fracture toughness, and can be engineered based on required strength or performance objectives for a given structure. To give few examples, leading aircraft manufacturers such as Airbus and Boeing have increased the application of FRP composites in their products from 3% to 20% in Airbus A380 and over 50% in Boeing 787 [4]. Figure 1-1 depicts one of the largest aircraft wings, with a 42.4 meter span, made of carbon FRP composite assembled for the Airbus Military A400M transport aircraft. In construction industry, FRP composites have become a good alternative for innovative construction and their applications have already been extended to upgrading and retrofitting the existing structures as well as constructing various types of off-shore platforms, buildings and bridges. Pontresina Bridge (Figure 1-1-d), with a span length of 12.5 m, was constructed in 1997 across the Flanz River in Switzerland. The structural truss system of the bridge was made of FRP composite [5].  5  (a)  (b)  (c)  (d)  Figure 1-1 Sample applications of FRP composite materials in different industries: (a) & (b) Airbus Military A400M aircraft CFRP wing (Martin Chainey / Airbus Military), (c) Cross-sections of glass fiber reinforced polymer (GFRP) structural members used, and (d) Pontresina Bridge made of FRP [5]  1.5  FRP Composite Materials Weaknesses and Application Limitations  Although FRP composite structures have proven to provide numerous advantages, initiation and propagation of cracks in these materials can affect their mechanical properties drastically. The most common FRP composite material failure modes are classified into fiber breakage, fiber-pull out, matrix cracking and interlaminar delamination (Figure 1-2). Fiber breakage 6  failure in FRP composites depends on the fibers strength distribution. This failure mode occurs when the subjected tensile stress in the FRP increases and leads to failure of low strength fibers in the laminate while the high strength fibers are still carrying the external load. For a better understanding of fiber breakage, distribution of the fibers strength in FRP has been investigated by several researchers and it was realized that the strength of glass fibers follows the Rayleigh distribution [6] while the strength of carbon fibers fits the Weibull and Gauss distributions [7]. In the cases where fiber strength is greater than the interface bond between the matrix and the fibers, the fibers can be pulled out of the matrix due to the tensile loading. In several cases, the fiber pull-out occurs at the fibers end, or at a fracture surface in the laminate. Another common failure mode in FRP composites attributed to the matrix material is called the matrix cracking. Since strength in the matrix is commonly lower than the fibers, normally the first stage of failure starts with matrix deterioration. Matrix cracking usually takes place before the entire FRP laminate matches its failure point, and it demonstrates a ductile failure in comparison to the brittle fiber failures. Among other failure modes, delamination is perhaps the most common failure mode and may occur because of a weak bonding between composite layers, existing cracks in the matrix material, broken fibers and fatigue or impact loadings. It can drastically reduce the structural stiffness and weaken the tensile or shear capacity of the FRP structures under service loads [8].  7  σmax  σmax Fiber Pull-out  Fiber Breakage  σmax  σmax (a)  (b)  σmax Matrix Cracking  σmax  σmax  Interlaminar Delamination  σmax (c)  (d)  Figure 1-2 FRP composite materials common failure modes: (a) Fiber breakage, (b) Fiber pull-out, (c) Matrix cracking, and (d) Interlaminar delamination  Mechanical defects are not the only cumbersome issue affecting the FRP performance and applications. Another undesirable feature of FRP composite materials is the complex heterogeneous and often stochastic properties of these materials which result in randomness and uncertainty in their manufacturing processes and material compositions in final products. The deterministic approaches used in many investigations ignored the spatial variability that exists in such material behaviour, especially at the micro-level scale, and this effect can entail errors into larger scale simulations. Traditionally, a common modeling approach for FRP’s is 8  the implementation of a homogenization technique [9] and relating a certain scale of material properties to the larger scale, mostly by averaging measured properties (Figure 1-3). Homogenization techniques can be categorized into micro-level, meso-level and macro-level, where the uncertainty and variability are assumed for different scales of material composition. At the micro scale, microstructure of the material composition is examined with methods such as the moving-window technique (MWT) [10]. This approach requires sophisticated experimental instruments and deals with the complexity of the FRP composite material conformation and may not be suitable for non-research purposes. Meanwhile, at the mesolevel, study of FRP composites has received considerable attention over the past decades and proven to be one of the effective ways to harness uncertainty to homogenization processes [11, 12]. However, difficulties with boundary condition assumptions in modeling still exist in this technique and have made its application for large structures challenging [13]. The largest scale in material composition studies is the macro-scale. Such studies avoid the microscopic complexity of composite materials and the numerous random variables. In macro-level techniques, an existing material’s behaviour and randomness are captured by coupon size test results while relying on continuum-mechanics based structural modeling formulations [14].  9  a1 << a2 << a3  a1 ~10-6 m Micro-scale  a3 ~1 m  a2 ~10-3 m Meso-scale  Macro-scale  Figure 1-3 Different scales of material modeling during homogenization of material properties  1.6  Experimental and Numerical Modeling of FRP Fracture Properties  In order to improve mechanical performance of FRP composite materials in the presence of process-induced or loading-induced cracks, extensive studies on the fracture properties of FRP composites have been performed, both experimentally and numerically. Experimentally, Owen and Bishop [15] applied the double edge-notched tensile test to measure the mode I critical stress intensity factor for varying orientations of UD glass FRP composites (Figure 14). Gaggar and Broutman [16] utilized both single and double edge-notched tensile tests as well as a notched bend test to extract the critical stress intensity factors. Mower and Li [17] summarized the experimental results from previous investigations and concluded that the linear elastic fracture mechanics (LEFM) is not a valid approach for FRP composite materials with long fibers and a nonlinear material constitutive model is required to accurately characterize the fracture energy toughness of FRP composites. The fracture energy toughness of unidirectional FRP composite materials under a double cantilever beam (DCB) test can be calculated using the modified beam theory, compliance calibration and modified compliance calibration methods [18]. The modified beam theory approach is recommended by O’Brien  10  and Martin [19] and has shown the most repeated value for the critical fracture energy toughness. In order to measure the mode II fracture energy toughness of the FRP composites, the end-notched flexure (ENF) test was suggested and employed by Davies et al. [20]. The ENF test is today one of the most recognized testing methods for mode II study of FRP composite materials. However, due to the unstable nature of this test, only the crack initiation values can be extracted from the test results. Edge delamination phenomena is another failure mode in FRP composite laminates studied by Lee [21] using the edge crack torsion (ECT) test. This test can be applied to extract the mode III interlaminar fracture properties of specimens. The extensive investigation on the fracture phenomenon of test samples showed delamination problems mainly consist of mixed mode fracture mechanics characteristics. Subsequently, the mixed mode bending (MMB) test has been designed as one of the recent methods for mixed mode characterization (mode I and II interaction).  11  Pgrip  Pgrip  a0 w  tck  tck L  a0  L  (a)  (b)  d1  Pgrip  Pgrip w L  a0  Pgrip Initial Crack  tck Z  Pgrip  Y  d1 d2  (c)  Pgrip  X  a0  L  (d)  Figure 1-4 Fracture energy toughness measurement fixtures: (a) Double Cantilever Beam, (b) EndNotched Flexure, (c) Edge Crack Torsion, and (d) Mixed Mode Bending tests [8]  Regarding numerical modeling of FRP composites delamination, numerous investigations have been performed over the past few decades. Hillerborg et al. [22] introduced a combination of finite element method and an analytical solution to simulate the crack growth in concrete structures. This approach is often referred to as ‘fictitious crack modeling’, where a traction-separation law instead of the conventional stress-strain relationship is utilized in the  12  crack-tip zone to capture the degradation of the material properties due to damage. Later, Xu and Needleman [23] applied an energy potential function to implement the cohesive zone model (CZM) during the analysis of interface debonding. CZM application was extended to FRP composites by Camacho and Ortiz [24], Camanho et al. [25], Blackman et al. [26], Gao and Bower [27], Segurado and LLorca [28], Cox and Yang [29] and Nishikawa et al. [30] while they have improved the cohesive interface models. Based on these reports, CZM has proven to be capable of modeling the ‘large process zones’—in the present case the FRP composite delamination. Despite its capability to model the progressive delamination, some severe disadvantages of applying large process zone have been noted and need to be tackled. These include numerical instability (elastic snap-back), reduction of stress intensity upon delamination initiation, and the fictitious softening of the original body in the process zone. In other investigations, a newly introduced feature of the finite element method, known as the extended finite element method (XFEM), has been more recently implemented for numerical modeling of progressive delamination in FRP composites. The original XFEM approach was first introduced by Belytschko and Black [31] and enhanced by Moёs et al. [32]. They implemented the concept of a partition of unity method (PUM), introduced earlier by Melenk and Babuška [33], to develop a method to model material discontinuity. In the basic XFEM approach, a Heaviside step function is implemented to model the crack surface by adding extra degrees of freedom to each node of the so called ‘enriched elements’ [32]. They introduced a framework capable of considering cracks and frictional contact with the zerothickness process zone in 2D problems. Xiao et al. [34] utilized this approach combined with a statically admissible stress recovery (SAR) technique in modeling cohesive cracks with a softening law composed of linear segments. Later, the approach was implemented by Unger  13  et al. [35] to model the cohesive crack in concrete specimens. The application of the 2D model was extended to composite materials by Benvenuti [36], who regularized XFEM for embedded cohesive interfaces. Sosa et al. [37] demonstrated the effectiveness of XFEM in 3D modeling of fiber metal laminate delaminations and compared their results with existing CZM results. XFEM and CZM concepts will be discussed in more detail in Chapters 4 and 5.  1.7  Randomness in FRP Fracture Properties  As addressed in Section 1.5, studying composite materials based on the assumption of Representative Volume Element (RVE) and the subsequent homogenization have been the basis of the work of several researchers. Kaw [38] used the RVE approach for unidirectional (UD) composites and Peng and Cao [39] developed a dual homogenization technique for woven fabrics. However, other investigations have shown that there can be considerable spatial non-uniformity both in UD composites [40] and woven fabrics [41, 42], which may hinder the full capability of RVE approach for fracture simulations [43]. There can be a variety of sources for such non-uniformity of material properties in composites. Examples include random distributions of fibers within samples, fiber penetration between layers, existence of voids within the matrix, human error in manufacturing process, uneven heating or cooling of samples during molding. Hence, there is a need for developing new models that can consider the heterogeneousness characteristics of FRP composites and include statistical distributions of their mechanical properties as well as the pre-existing defects in test specimens [44, 45]. Stochastic modeling of the mechanical behaviour of composites can be especially important in predicting critical loads and critical displacement values, as well as crack formation patterns in large structures [46, 47]. Among the most recent works on 14  stochastic modeling of FRP composites, Ashcroft et al. [48] emphasized the effect of material uncertainty and non-uniformity in predicting delamination phenomena.  1.8  Motivation and Objectives of the Work  Based on the above background review, a number of investigations have been performed to study the mechanical properties of FRP composites. In particular, fracture behaviour of these materials is an interesting topic where a variety of experimental tests and numerical procedures have been proposed by different research groups. The heterogeneous nature of FRP fracture properties has made the research in this area a challenging task. With respect to numerical simulations, different numerical methods have been applied to extract the correct behaviour of FRP composites and the results demonstrated the necessity for further investigation in this field for more realistic simulations.  Among the various numerical  methods employed to model the delamination in composite materials, the interface element method combined with a cohesive law has received great attention by numerous researchers to date. Espinosa et al. [49] implemented this method for modeling dynamic delamination of woven GFPR composite materials with an anisotropic visco-plastic material model in conjunction with a cohesive law. Cohesive zone model properties such as maximum interface strength, fracture parameters, penalty stiffness and cohesive zone length were studied by Turon et al. [50] and Harper and Hallett [51] to overcome the existing implementation obstacles of cohesive zone models in numerical simulations. A generalized framework for implementing the cohesive XFEM in modeling delamination was introduced by Benvenuti [36]. In that work, the fundamentals of XFEM with cohesive law characteristics were studied and examples of DCB test were modeled. In addition to FEM-based methods, a mesh-free 15  based method was utilized by Barbieri and Meo [52] to extensively simulate the mode I and II, DCB, ENF and end-loaded split (ELN), tests in 2D with the focus on nonlinear aspects of crack problems. Most of the above mentioned studies, however, are based on deterministic properties of composite materials, especially their fracture properties, while other experimental results demonstrated the variability of the composite fracture phenomenon [48]. Among recent stochastic works, Ashcroft et al. [48] introduced micro-structure randomness in the form of fracture properties of CFRP composite materials into the numerical simulation of DCB tests using interface cohesive elements in the FEM model. Non-uniformity and random distribution of material fracture properties were considered by means of uniform and Weibull distributions and results emphasized the need for further studies on including micro structural randomness for accurate predictions of fracture performance of composite laminates. Therefore, the present thesis is primarily aimed at developing and examining an enhanced numerical approach for simulating the composite fracture tests considering both material and geometric nonlinearities along with stochastic fracture properties. A simplified approach is introduced to implement the cohesive zone model and hence to avoid the numerical softening due to existence of a large process zone. An ABAQUS user-element subroutine is developed and linked to MATLAB to implement the nonlinear XFEM for DCB test simulation (mode I fracture). The cohesive zone model is associated with enriched elements to consider a bilinear traction-separation law in the crack front using the XFEM contact model following the work of Khoei et al. [53]. The model is also implemented to investigate the ENF test for mode II fracture properties. Stochastic distributions are employed to the fracture properties of the material via the bilinear traction-separation law, and results are compared with both available data in the literature and a set of performed tests on Poly Phenylene Sulfide (PPS)/Glass UD 16  comp posites. Finaally, a set of o sensitivity y analyses aare perform med to identtify the effeects of differrent model variables v succh as the messh size, coheesive zone leength and peenalty stiffneess.  1.8.1  Potentia al for Practical Applica ations and E Expected Orriginality  As ad ddressed beffore, composite materialls are rapidlyy replacing m metallic matterials in diffferent indusstries. Moulding process is an esseential part of any compposite manuffacturing linne that dealss with form ming structu ural compon nents for vehicles, airpplanes, boaats, etc. A rising challenge for opttimal implem mentation off FRP compoosites in form ming processses, howeverr, is to underrstand their behaviour and a conformaability againnst heat and pressure varriations as w well as differrent geomettries of the mould. For instance, dduring moulding processs of a thickk FRP comp posite into a curved/ben nt part, it is possible thaat in the corrner of the ppart, delaminnation between layers in nitiates and propagates p and a eventuallly be a reasoon for the reejection of thhe part ure 1-5) [54]]. The samee problem may m happen during serviice under exxcessive loaadings, (Figu even for a part th hat had been n originally manufactureed with no fflaw. Thereffore, it is reqquired for leading indu ustries to em mploy advaanced simulaation tools as part of their design and manu ufacturing prrocesses to predict p damaage phenom mena (in the ppresent casee, delaminatiion) in FRP composite materials as a accurately y as possibble, while nnot ignoringg the unavooidable variaation in material properties of compo osites (in the current casee, the fracturre propertiess).  17  Delamination onset in the bent/corner aarea  Step 1  Step 2  Step 3  (a)  (b)  Figure F 1-5 Fo orming processs of a heated FRP F compositte laminate: (aa) Initiation oof the forming (Step 1), co ompression an nd cooling cyccle (Step 2), reemoving the m male die and p possible dimen nsional distorttion of the part p (Step 3); and (b) Delam mination spottted in the corn ner of a comprression mould ded componen nt [54]  The expected oriiginality of the present work withinn the contexxt of applicaation of advvanced numeerical modelling approacches for com mposite desiggn and manuufacturing caan be summ marized as folllows.   Performin ng mode I and II fraccture tests for a PPS//Glass UD laminate ussed in aerospacee and automo otive industrries.    For the first f time in nvestigating the observeed randomneess (non-reppeatability) in the material response r from m the abovee fracture tests.    Enhancin ng the existiing XFEM approach inn the literatture to moddel the com mposite fracture tests underr uncertainty y (combininng XFEM with contaact and cohhesive  18  modeling capabilities under one framework, along with stochastic fracture properties), hence moving towards more reliable damage prediction tools.   Implementing stochastic characteristics of the tested PPS/Glass UD composites in the XFM simulations and validating with experimental data.    Studying the effect of different modeling factors (such as penalty stiffness) in the enhanced XFEM approach using a series of sensitivity analysis.  1.9  Thesis Outline  This thesis is organized into five chapters. The first chapter, presented above, focused on a general literature review on experimental and numerical modeling of the fracture mechanics of FRP composite materials and also the necessity for considering the variability in microstructure of these materials. In Chapter 2, damage modeling and fracture properties of FRP composites are reviewed. Also in this chapter, the XFEM method and its applications are reviewed and discussed. In the final part of Chapter 2, cohesive zone and contact surface implementations of XFEM are described and finite element formulations are presented. In Chapter 3, results from XFEM simulations are compared to those extracted from other numerical methods and benchmark experimental tests in the literature. Modeling parameters such as penalty stiffness, mesh size and cohesive zone length are studied via a set of sensitivity analysis to examine the presented XFEM approach effectiveness in modeling fracture properties of UD composites.  Chapter 4 presents the fabrication process and  experimental procedures used to obtain elastic moduli and fracture properties of PPS/Glass samples. Following, stochastic features of the measured fracture properties are introduced to the XFEM model. Results are then compared with those from experiments using DCB and 19  ENF test set-ups. Finally, in Chapter 5, the undertaken numerical and experimental procedures as well as the main results are summarized. Future work recommendations are also included in Chapter 5.  20  2 Chapter: Background In this chapter, elastic mechanical behaviour, damage properties and fracture mechanics of FRP composite materials along with different underlying modes of failure are first reviewed. Advantages of damage modeling and fracture mechanics will be compared in modeling the post-failure behaviour of FRP composite materials. Next, basic definition and properties of XFEM will be described and its application in modeling the LEFM and EPFM will be discussed. The implementation of damage mechanics through cohesive zone model (CZM) will also be presented and effective parameters in such a model will be summarized.  2.1  Elastic Mechanical Behaviour of FRP composite Materials  As mentioned in Chapter 1, in FRP the matrix provides the integrity of the composite by holding fibers together. It also has greater elongation characteristic than fibers which forces the fibers to carry the maximum load before the matrix fails. Fibers, on the other hand, provide high strength and stiffness to the material system. Such a material composition will lead to anisotropic material properties and entails an appropriate technique for extracting the global (macro-level) material behaviour. The conventional elastic constitutive relationship between stress ij and strain  kl for a FRP composite material, similar to an orthotropic material, can be written as: 0 0 0   11   11  C1111 C1122 C1133   C 0 0 0    22   22   1122 C2222 C2233  33  C1133 C2233 C3333 0 0 0    33      0 0 C1212 0 0   212   12   0  13   0 0 0 0 C1313 0   213       0 0 0 0 C1323  2 23   23   0  (2-1)  21  where Cijkl (i,j,k,l = 1, 2, 3) are the stiffness matrix, C, components. For easier demonstration of the material parameters, especially for FEM implementation, the following vector-form compliance equation is used to describe the material constitutive behaviour:    S where S  C1  (2-2)  where the compliance matrix, S, is re-arranged as follows:   1  E  1   12  E  1   13  E S 1   0    0    0      21   31     32  0  0  0  0  1 E3  0  0  0  0  1 G12  0  0  0  0  1 G13  0  0  0  0  E2  1 E2      23 E2  E3 E3   0    0    0    0    0   1   G23   (2-3)  Ei is the elastic modulus in direction i, Gij are the shear moduli and vij are the Poisson’s ratios. Generally, the unidirectional composite materials have a transversely isotropic behaviour which yields a relationship between material properties as follows:  E1  E2 , 12 13 , G12  G13  (2-4)  Subsequently, the compliance matrix can be redefined as:  22   1  E  1  v12  E1  v  12 E S 1   0   0    0   v12 E1 1 E1 v  23 E1  v12 E3 v  32 E3 1 E3  0      0  0  0  0  0  0  0  1 G12  0  0  0  0  1 G12  0  0  0  0   0   0    0    0   0   1   G32   (2-5)  The American Society for Testing and Materials (ASTM) standards [55-57] has proposed several experimental tests for extracting the material elastic constants at macro-scale.  2.2  Failure Modes of FRP Components  Most FRP composite materials demonstrate a brittle behaviour when imposed to in-plane loading. Typically, the response curve of such laminated FRP composites under tensile test starts from the origin and increases linearly up to the failure point. At this point, the material faces a form of irreversible damage (e.g., fiber breakage). After this point, it is most likely that the load-deformation curve drops to zero and the FRP composite loses its capacity to carry further load. The stress-strain relationships for several fiber materials are demonstrated in Figure 2-1.  23  6000 E-Glass Aramid-49 Standard Carbon High-Modulus Carbon Ultra High-Modulus Reinforcing Steel  5000  Stress (MPa)  4000 3000 2000 1000 0 0  1  2  3  4  5  Strain (%)  Figure 2-1 Comparison of tensile stress-strain curve of fiber materials [58]  Damage initiation in FRP composite materials has a direct relationship with matrix and the fiber material properties and the existing flaws in the structure. A difference between the elongation limit of the matrix and fibers is a dominant parameter in tensile failure of a laminated FRP composite. In the case where failure strain of fibers is larger than the matrix, the matrix cracking is expected to happen earlier. However, if the fiber ratio of the composite is large enough (greater than about 10%), fibers can continue carrying the load up to their rupture. On the other hand, if the fibers have lower elongation than the matrix, which is the most common case, it will force fibers to carry the maximum load according to their capacity and fail due to breakage. Both cases can theoretically be justified when no flaw exists in the FRP composite material (Figure 2-2). Incomplete bonding at the interface of the matrix and fibers, air entrapment, uneven distribution of fibers within the matrix, and premature cracks 24  between layers of laminate are several common manufacturing defects in FRP composites [58].  Stress, σ  Stress, σ σf,ult  σf,ult  Fiber  Fiber  σ’f σm,ult  σm,ult Strain, ε  Matrix εm,ult  (a)  εf,ult  σ’m  Matrix εf,ult  Strain, ε εm,ult  (b)  Figure 2-2 Stress-strain curves for the matrix and fiber materials: (a) Fibers have larger failure strain than that the matrix, and (b) The matrix has larger failure strain than fiber [58]  Depending on the matrix and fiber material properties and the above mentioned imperfections, matrix cracking, fiber pull-out, fiber breakage or delamination is expected to occur and cause a crack formation in a FRP structure due to the intrinsic brittleness of these materials under extensive loadings. These damage modes will reduce the structure’s capacity to endure extra loadings; however, they may not lead to complete failure and collapse of the structure. Examining the initiation, stability and growth of defects are directly linked to the comprehensive study of fracture mechanics and the theory of plasticity.  25  2.3  Fracture Mechanics  Structural damage associated with the crack causes the local stresses and global deformations to increase for the body near the cracked area. Three independent modes are applied to define any coupled fracture deformation of a given structure (Figure 2-3) [59]: 1. Opening mode (mode I), when two faces of the crack are pulled away in the crack’s plane. 2. Sliding (shearing) mode (mode II), when two faces of the crack are sliding over each other in the crack’s plane. 3. Tearing mode (mode III), when two faces of the crack are taken apart out of the crack’s plane.  (a)  (b)  (c)  Figure 2-3 Dominant fracture modes of a cracked body: (a) Opening mode (mode I), (b) Sliding (shearing) mode (mode II), and (c) Tearing mode (mode III) [59]  26  Degradation and loss of the integrity of a structure can be drastically increased due to the crack formation (and propagation) which entails a careful scrutiny of different fracture modes with respect to the composite material composition and structural loadings. In this regard, the classical fracture mechanics can be categorized into linear elastic and elastic-plastic fracture depending on the type of material under study. Linear elastic fracture mechanics (LEFM) discusses deformation and stress fields around the crack-tip when a small plastified zone forms in the cracked region relative to the crack length. In such a case, singular stress fields emerge at the crack-tip. In LEFM, the stress intensity factors (SIFs) are applied to assess the stability of the crack by comparing it to critical SIFs extracted from experiments. The Westergaard [60] solution is a well-known approach in applying SIFs to estimate the displacement and stress fields near the crack-tip. Despite the capability of LEFM in estimating the fracture properties of brittle materials, with the extension of the plastic zone or the fracture process zone the singular stress fields vanish from the crack-tip region (e.g., crack propagation in steel structures). For such cases, elastic-plastic fracture mechanics (EPFM) proposes more accurate solutions and considers a plastic zone in front of the crack-tip region where extensive plastic deformation is expected to emerge. Wells [61] proposed the crack-tip opening displacement (CTOD) as a failure threshold for a structure. As an extension to the nonlinear analysis of such plastic zones, the J-integral was introduced by Rice [62] to accurately evaluate the energy release rate of materials when the Griffith theory is imprecise. In general, fracture mechanics of a composite material is more challenging in comparison to other homogenous materials. The heterogeneous composition of FRP composite materials entail complex fracture phenomenon such as delamination. As mentioned before, delamination is one of the failure modes in FRP composite materials and can occur because of 27  a weak bonding between composite layers, an existing crack in the matrix material, broken fibers, fatigue, impact loadings, etc. In multi-layer laminated FRP composite structures, fibers in each layer are confided by adjacent layers on the top and bottom faces and the matrix acts as the bonding agent. When a crack forms in the matrix, the brittle nature of the matrix lets the LEFM accurately assess the crack stability and propagation pattern. On the other hand, in unidirectional laminated FRP composite structures, especially when formed using a moulding process with high compression pressure, fibers can penetrate into adjacent layers and the bonding between layers can be affected by fibers penetration (called fiber bridging). In this case, fibers onset an extra resisting force and prevent the crack from opening and the fracture behaviour of the structure demonstrates similar properties to those experienced during an elastic-plastic fracture. This aspect of FRP composites has increased the need for more advanced techniques of modeling their fracture phenomena and the governing failure modes.  2.4  Damage Mechanics Models  Several numerical investigations have been focused on the plasticity and damage modeling to consider defects in FRP composite structures. In these models, some failure criteria are predefined for the FRP composite and, based on the stress fields in the structure, failed elements are identified. Depending on the failure mode and numerical modeling technique, distorted elements may be eliminated from the model, or their stiffness would be degraded to simulate the material softening. Erdogan and Sih [63] introduced the maximum stress criterion which compares the principal stresses in each direction of an element with the material allowable stresses to evaluate the dominant failure mode. Hoffman [64] proposed that when the maximum strain criterion is similar to the maximum stress criteria, the principal strains in the 28  material directions are compared to the material strains capacity. Tsai [65] applied the Hill’s yield surface and redefined the Tsai-Hill theory. The yield surface based on the Tsai-Hill criterion is defined as:   11   T11  2       22   T 22  2       33   T33   1 1 1   2  2  2  T11 T 22 T33  2       12  T12   2       23   T 23  2       13   T13      2     11 22   33  22   11 33   YS normal   (2-6)  where  ij , Tij and YS normal are the stresses in different material directions, material strengths and normalized failure threshold, respectively. Tsai and Wu [66] also presented a failure theory based on the strengths criterion for anisotropic materials. The Tsai-Wu failure criterion could differentiate between material tensile and compressive strength and ignores the interaction between failure modes. Based on this criterion, the failure surface can be expressed as follows:  Fij  i j  Fi i  YS normal  (2-7)  where  i , Fi and Fij are the re-arranged vector of stress tensor, second and fourth order strength tensors, respectively. Specifically for unidirectional (UD) laminated composite materials, Hashin [67] introduced a set of failure modes based on the combination of fibers and matrix strengths. The failure criteria proposed by Hashin [67] can be summarized as follows: 1- Tensile fiber failure for  11 ≥ 0:  29  2  Failure   11   122   132  1      S122  1 No Failure  XT   (2-8)  2- Compressive fiber failure for  11 < 0:    11   XC  2    1     1   Failure  (2-9)  No Failure  3- Tensile matrix failure for  22   33 > 0: 2  2 Failure   22   33   23   22 33  122   132  1       2 S 23 S122  1 No Failure   YT  (2-10)  4- Compressive matrix failure for  22   33 < 0:     Y  2          33   22 2 33  C   1 22 4 S 23  2 S 23   YC    232   22 33 S 232     122   132 S122    2    Failure  1   1 No Failure  (2-11)  5- Interlaminar tensile failure for  33 > 0: 2  Failure   33   1      ZT   1 No Failure  (2-12)  6- Interlaminar compression failure for  33 < 0:  30  2  Failure   33   1      Z C   1 No Failure  (2-13)  where X T , YT , Z T , X C , YC , ZC , S12 , S13 and S 23 are the tensile, compressive and shear strengths in different material directions.  2.5  Comparison between Fracture Mechanics and Damage Mechanics  In general, fracture mechanics focuses on local discontinues where macro-cracks are present. It observes the singular stress and strain fields in front of the crack-tip, especially in brittle materials, and provides an accurate estimation of the damage evolution when a flaw exists in the structure. On the other hand, damage mechanics tends to evaluate the stress and strain state in the structure and assess whether the material faces degradation in each loading step. In addition, it can be applied to study the global behaviour of the structure and predict the failure initiation and expansion trend due to the material deterioration. Hence, depending on the problem encountered, an appropriate approach should be selected to accurately evaluate the structural behaviour of the composite. In the present work, the focus on the local damage modelling of FRP composite materials would prompt the implementation of fracture mechanics due to the local nature of delamination and brittle properties of such materials. However, the large process zone existing in the delamination front will require the implementation of elastic-plastic fracture mechanics to ensure the accurate evaluation of the structural deformation (this will be formulated in detail in Section 3.2).  31  2.6  Introduction of Extended Finite Element Method  Numerical modeling is an important part of most engineering applications. In many cases, a structure’s dimensions and the test set-up configurations cause limitations when performing full-scale experimental studies, increasing the demand for undertaking more numerical analyses. A variety of numerical modeling techniques have been proposed in the past decades and can be categorized into mesh free methods such as Smoothed Particle Hydrodynamics (SPH), Element-Free Galerkin Method (EFGM), Finite Difference Method (FDM) and Meshless Methods as well as mesh-based methods such as Finite Element Method (FEM) and Boundary Element Method (BEM). The FEM can be directly enhanced and used in modeling discontinuities by introducing the partition of unity method (PUM), proposed by Melenk and Babuška [33], to approximation functions. The basis of PUM is similar to regular finite element approximation where the summation of all shape functions at any Gauss point is equal to one. The new method is known as the extended finite element method (XFEM) [31-32]. The XFEM has demonstrated to be a more accurate and convenient solution where the conventional finite element produces rough or highly oscillatory results. In XFEM, similar to conventional FEM, the finite element mesh is generated regardless of discontinuities locations. Then, specific search algorithms such as the level-set or fast marching methods are utilized to identify the location of any discontinuity with respect to the existing mesh and distinguish the different types of required enrichments for the affected mesh elements. Next, additional auxiliary degrees of freedom are added to the conventional FEM approximation in selected nodes around the discontinuity. These degrees of freedom assist the model in capturing the displacement jumps caused by discontinuities. 32  Assume a discontinuity (a crack) within an arbitrary finite element mesh (Figure 2-4). The displacement field of point x, u g (x) , inside the domain can be described in two parts; the conventional finite element approximation, and the XFEM enriched field defining the discontinuity [32]: u g ( x)     I nI N all  I  ( x )u Iord     J n J N  J  ( x ) ( x )u Jenr  (2-14)  f  where  (x) is the conventional shape function (also often shown in the literature by N ( x) ),   (x) is the general enrichment function, Nall is the finite element mesh nodes and N f is the enriched nodes of the mesh, u  ord  is the classic degrees of freedom at each node and u  enr  is  the additional enriched degree of freedom at the J th enriched node.  Crack  Figure 2-4 An arbitrary finite element mesh with a discontinuity (circles represent the enriched nodes of the mesh)  In order to choose the enrichment function, any discontinuous function in the problem domain can be employed to estimate the displacement field approximation in vicinity of the crack. A 33  function that satisfies such a requirement is the Heaviside step function ( H (x) ). It gains a ‘+1’ value on one side of the crack and ‘-1’ on the other side of the crack and can be utilized when the crack propagation is modeled by straight line segments. To find the Heaviside function value at each node of an element, tangential and normal vectors of the crack surface curve should be measured. If X * is the nearest point of a crack to the point X’, Figure 2-5, and en is the unit normal vector of the crack alignment in which es  en  ez ( es is the unit tangential vector), then using a scalar product between the distance vector of the element’s nodes and the normal vector of the crack surface, the Heaviside function value can be calculated.  1 ; if ( X ' X * ).en  0 H (x)     1 ; otherwise  (2-15)  en es X* S X  Figure 2-5 Unit tangential and normal vectors for the Heaviside function and nearest point to X on the crack surface; X*  Another set of functions utilized in XFEM are those extracted from an analytical solution of the near crack-tip displacement fields in LEFM as follows [68]: 34  u xtip   1 2G12  r      K I cos 3  4v12  cos   K II sin 3  4v12  2  cos  2  2 2   (2-16)  u tip y   1 2G12  r      K I sin 3  4v12  sin    K II cos 3  4v12  2  cos  2  2 2   (2-17)  u ztip   2 G12  r  K III sin 2 2  (2-18)  tip where u xtip , u tip y , u z , K I , K II and K III are the crack-tip displacement field and stress intensity  factors (SIFs) in the three standard fracture modes, and (r, θ) are the local polar coordinate system at the crack-tip. Subsequently, the extracted crack-tip enrichment functions were proposed by Moёs et al. [32] as:  F  enr    (r ,  )  4  i 1          r sin , r cos , r sin sin  , r cos sin   2 2 2 2    Among presented functions,  r sin   2  (2-19)  is the only discontinuous function and the remaining  functions are continuous. A more general set of enrichment functions can be achieved by studying the crack asymptotic displacement fields in anisotropic materials. For such a case, with general boundary conditions and the structure subjected to arbitrary forces, the following characteristic equation can be obtained using the equilibrium of forces and compatibility conditions [69]: a11 4  2a16  3   2a12  a66   2  2a26   a22  0  (2-20)  where aij is the material compliance matrix components. 35  Lekhnitskii [69] illustrated that the roots of Equation (2-20) are always complex or purely imaginary    k    kx  i ky , k  1, 2 with the conjugate pairs as 1 , 1 and  2 ,  2 .  Employing the above equations, the displacement fields in the vicinity of the crack-tip were elicited by Sih et al. [70] by means of analytical functions and complex variables as follows:      (2-21)      (2-22)  u xtip  K I   1  2r Re  μ1 p 2 cos θ  μ 2 sin θ  μ1 p1 cos θ  μ1 sin θ  π  μ1  μ 2   u tip y  KI    1 2r Re  μ1 q 2 cos θ  μ 2 sin θ  μ1 q1 cos θ  μ1 sin θ  π   μ1  μ 2  For pure mode II:      (2-23)      (2-24)  u xtip  K II   1  2r p 2 cos θ  μ2 sin θ  p1 cos θ  μ1 sin θ  Re  π  μ1  μ2   u tip y  K II   1  2r q2 cos θ  μ 2 sin θ  q1 cos θ  μ1 sin θ  Re  π  μ1  μ 2   where Re denotes the real part of complex variable. pk and qk can also be defined [70]:  pk  a11k2  a12  a16 k  qk  a12 k   a22  k   a26  (2-25) (2-26)  Using the above equations, the near crack-tip enrichment functions for the crack in an anisotropic material are expressed as [71]:  F  enr      4   (r, )i 1   r cos 1 g1 ( ) , r cos 2 g 2 ( ) , r sin 1 g1 ( ) , r sin 2 g 2 ( )  2 2 2 2    (2-27)  36  where gk ( ) and k , k  1, 2  are defined as: g k ( )   cos     sin     xy sin   2  xy  2  (2-28)   ky sin     cos   kx sin      k  arctan  (2-29)  Moёs et al. [32] substituted the Heaviside and near crack-tip enrichment functions in the XFEM approximation and presented the following equation: u g ( x)    ( x)u I  I        l 2 lK node  J  ( x) H ( x)bJH     F 2  i  enr    k 1 kK node  J g n J N node  I nI N all   l ( x)  cltip  i         k ( x)  cktip i F enr ( x) i    1  i  1     ( x) i      (2-30)  2  g where N node is the set of nodes that have a crack face (but not a crack-tip) in their influence  domain, b H is the vector (in multi-dimensional problems) of additional degrees of freedom which are related to the modeling of crack faces (not crack-tips), c tip is the vector of enr additional degrees of freedom for modeling crack-tips, F ( x) is the crack-tip enrichment 2 are the sets of nodes associated with the crack-tip 1 and 2 in their function and K 1node and K node  influence domains, respectively. Figure 2-6 illustrates a finite element mesh for modeling the existing discontinuity in Figure 2-4. Nodes that need to be enriched with Heaviside and near crack-tip functions are distinguished by circles and rectangles, respectively. Other nodes that are not affected by the crack remain well within the classical finite element framework.  37  Influence domain of node J  J  Crack  Enriched by Heaviside  Enriched by near-tip functions  Figure 2-6 The influence domain of node J in an arbitrary finite element mesh  The XFEM displacement approximation in Equation (2-30) can be implemented in numerical modeling with LEFM to predict the displacement field in a cracked structure. Considering the total potential energy governing the problem, we can write:    V   ij    ij dV   f i b   u k dV   f i t   u k dt V  t  (2-31)  where V , t , f i b , fi t and u k denote the body domain, traction surfaces, body forces, external forces vectors and displacement field, respectively. Discretizing the Equation (2-31) and applying the variational formulation, the following linear algebraic equation can be written:  K 0u h  R 0  (2-32)  38  Where K0, R0 and u h denote the initial stiffness matrix, right hand side vector of external forces and general vector containing all the nodal parameters including the ordinary, crack face and the crack-tip enriched degrees of freedom vectors, respectively:    u h  u ord , b H , c tip    T  (2-33)  In Equation (2-32), the initial stiffness matrix and the right hand side vector of external forces are also defined as:   K iju u  H ord K 0   K ijb u  K ctipu ord  ij  ord ord    R0  R0  u ord  K iju  ord H  b  K ijb K ijc  , R0  bH  K iju  H H  b  b  c tip  c  K ijb  tip H  , R0  ord tip  K ijc  H tip  c  tip tip  c          T  (2-34)  (2-35)  The stiffness arrays K ijrs ( r, s  u ord , b H , c tip ) in Equation (2-34) include the classical, enriched and coupled components of XFEM approximation: K ijrs   e ( Bir ) T DB sj d   ( r, s  u ord , b H , c tip )  (2-36)  where B is the matrix of derivatives of shape functions and is defined as:  Biu  ord   N i, X  0   0   N i ,Y  0   N i , Z  0 N i ,Y 0 N i, X N i ,Z 0  0  0  N i ,Z   0  N i ,Y   N i , X   (2-37)  39  Bib  H   N i H  H i , X  0   0    N i H  H i ,Y  0    N i H  H i ,Z    Bic  tip  N i H  H i , X N i H  H i ,Z   N F enr k  F enr k i  i  0   0  k k  N i F enr  Fi enr  0   k k enr  Fi enr  N i F      0   0  N i H  H i ,Z   0  N i H  H i ,Y   N i H  H i , X   0 N i H  H i ,Y 0  0    ,X  N F  0 enr k  i   N F N F  ,Y  ,Z  k  enr k enr k  0  N F  0 enr k  i   Fi enr  k   Fi enr  k  0    ,Y  0  i  i   Fi enr  (2-38)    N F N F  enr k  i  i  k  0  ,X ,Z   Fi enr  enr k   Fi enr  Fi  k  enr k       ,Z     ,Y   ,X        (2-39)  where H i and Fi enr are the Heaviside and crack-tip enrichment functions at element nodes. Normalizing the enrichment functions by mean of deducting the enrichment functions at element nodes will readily satisfy the PUM fundamental (i.e., the summation of normalized enrichment functions for a given node will be the unity). Also, the k index is used to differentiate between different enrichment functions, respectively.  2.7  Application of XFEM in Linear Elastic Fracture Mechanics  As addressed earlier, the linear elastic fracture mechanics (LEFM) is a suitable tool for the analysis of many materials with a small plastic or process zone in front of the crack-tip. The SIFs are the only required variables for assessing the stress and displacement fields around the crack-tip under different loading and subsequently, measuring the stability of the crack by  40  comparing SIFs to their critical values extracted from the experiments. In terms of numerical modeling of LEFM problems, XFEM provides an accurate estimation of the stress and displacement fields around the crack-tip using the Heaviside and crack-tip enrichment functions. To further the application of XEFM in LEFM, it is also required to introduce a post-processing procedure to evaluate the SIFs. For this purpose, the interaction integral, known as M-integral, was proposed by Moёs et al. [32] and Dolbow et al. [72] for isotropic materials and extended to orthotropic materials by Kim and Paulino [73]. The interaction integral is essentially derived from the J-integral introduced by Rice [72] as in Equation (240). Jk   u j   Wn k  t i X k   c       d   (2-40)  where W  1 2   ij ij is the strain energy density, X i , ti and nk represent Lagrangian coordinates, vectors of surface traction and vectors of outward normal direction, respectively, and the integral paths Γ and Γc denote far-field and crack surface paths, respectively, as shown in Figure 2-7.  41  Г VГ x2 x1 θ  X2 Vε  Гε  Гc X1 Figure 2-7 Local crack-tip coordinates and the contour Γ and its interior area, VΓ  The form of equation (2-40) is not perfectly suited for FEM implementations and an equivalent form can be obtained by applying the divergence theorem and several additional assumptions, as discussed by Kim and Paulino [73]:    2 u j u j  ij    W q',k dV    ij   J k     ij  X k  V  V    X j X k X k    1 Cijlm    q ' dV lm ij  2 X  k   (2-41)  where q ' is a function smoothly changing from q'  0 at the exterior boundary Γ to q '  1 near the crack-tip. In the present study, q ' is assumed to be varying linearly from 1 to 0. Noting that the value of q is constant near the crack-tip area, distinguished by shaded unaffected elements in Figure 2-7, the gradient of q ' vanishes in Equation (2-41). By adopting auxiliary displacement, stress and strain fields ( ukAuxilairy,  ijAuxilairy ,  Auxilairy ) and ij  superimposing them to displacement, stress and strain fields ( uk ,  ij ,  ), the superimposed ij  J-integral will contain the following parts:  42  J Superimposed  J  J Auxiliary  M  (2-42)  where M is the interaction integral in the local crack-tip 2D Polar coordinate and consists of terms involving products of actual and auxiliary fields. The Cartesian coordinate definition of M-integral is then as follows:  u iaux  u M k    ij   ijaux i  W  W aux q ',k dV X k X k  V        2 u aux  ijaux i      ij  X j X k X k V         Cijlm aux     q ' dV  X k ij lm     (2-43) (i, j , k  1, 2, 3)  To extract the M-integral in local crack-tip Polar coordinates, a simple transformation can be employed:  M Local  M1 cos  M 2 sin  (2-44)  Now, if one wants to relate the superimposed J-integral to the energy release rate, G, of the two actual and auxiliary fields, the following relationship between SIFs and M-integral can be derived:  J  G  GI  GII  c11K I2  c12 K I K II  c22 K II2      2    (2-45)        J Superimposed  c11 K I  K Iaux  c12 K I  K Iaux K II  K IIaux  c22 K II  K IIaux        2  J Auxiliary  c11 K Iaux  c12 K Iaux K IIaux  c22 K IIaux      2  (2-46)    (2-47)    (2-48)  2  M Local  2c11K I K Iaux  c12 K Iaux K II  K I K IIaux  2c22 K II K IIaux where c11 , c12 and c 22 are defined as follows:  43  c11    a 22  1   2   Im 2  1  2   (2-49)  c12    a22  1  a11   Im Im1  2  2  1 2  2  (2-50)  c22   a11 Im1   2  2  (2-51)  Finally, by solving the following algebraic equations for the auxiliary field for mode I and II, SIFs can be extracted:  M Local  2c11K I  c12 K II  for K Iaux  1.0 and K IIaux  0.0  (2-52)  M Local  2c11K I  c12 K II  for K Iaux  1.0 and K IIaux  0.0  (2-53)  2.8  Application of XFEM in Elasto-Plastic Fracture Mechanics  In contrast to LEFM, EPFM deals with ductile materials with a relatively large plastic or process zone in front of the crack-tip. Irwin [74] proposed a simple plastic zone correction to SIFs in order to consider the plastification effects. Alternative solutions were also introduced by Dugdale [75] and Barenblatt [76] to correct plastic zone characteristics. Wells [61] offered crack-tip opening displacement (CTOD) as an independent variable to measure the fracture energy toughness of materials. This approach covered not only the LEFM but also established a way to investigate the EPFM in materials. Later, Rice [62] introduced a path independent contour integral, the J-integral, to describe the stress and strain distribution near the crack-tip. It is worth mentioning that the J-integral was earlier extracted by Eshelby [77], however; the application of this method in studying LEFM and EPFM was proposed by Rice [62].  44  As mentioned in Section 4.2, the interaction integral method has been frequently adopted in LEFM crack analyses to evaluate mixed mode SIFs. One disadvantage of this method is related to the dependency of the method on the auxiliary field exact solution. In terms of EPFM, XFEM potential has left behind such limitations and provided a convenient solution. Motamedi and Mohammadi [78 - 79] implemented dynamic J-integral and CTODs to assess the SIFs in orthotropic materials. Definition of Equation (2-41) can be applied to measure the tangential component of the J-integral which corresponds to the rate of changes in the potential energy per unit crack extension, namely, the energy release rate (G):  G  J1 cos 0  J 2 sin 0  (2-54)  where 0 is the crack angle. In order to accurately evaluate the stress and displacement fields around the crack-tip, the component separation method, proposed by Aliabadi and Sollero [80], was employed. They discussed that the following relationship between the stress intensity factors and the CTODs of the crack faces can be obtained:  n   8r   D11 D12   K I             D21 D22   K II   t   (2-55)  with   p  1 p2   p  p2   , D12  I m 1  D11  I m 2 1       1 2 2     1   q  1q 2   q  q2   , D22  I m 1  D21  I m 2 1  1   2   1   2   (2-56)  Where p1, p2, q1 and q2 are defined as follows:  45  p1  a11 1  12   a12    (2-57)    p2  a11 1 22  a12 q1   q2   (2-58)  a12 1  12   a22  (2-59)  1  a12 1  12   a22  (2-60)  2  Another set of relationships between the energy release rate and the SIFs was proposed by Wu [81]:  G  (1/ 2) K T L(c) 1 K  (2-61)  where L(c) was defined by Dongye and Ting [82] for orthotropic materials with the symmetry planes coinciding with the coordinate planes. The non-zero components of L (c ) are: (2-62)  L33  C 55 C 44  C66C22 L11  C11C66 L22      1 2  (2-63)  where Cij is the constitutive/stiffness coefficients and   (C11C 22  C122 ) C 66 C 66  (2-64)    (C66  C11C22 ) 2  (C12  C66 ) 2  (2-65)  Therefore, the ratio of opening to sliding displacements, RCOD , can be written as:  46  RCOD    2 D21K I  D22 K II  1 D11K I  D12 K II  (2-66)  And the ratio of dynamic stress intensity factors, RSIF , can be obtained as:  RSIF   K I RCOD D21  D22  K II D21  RCOD D11  (2-67)  Substitution for K I from Equation (2-61) into Equation (2-67) leads to the following relationship for K II : K II   2G 2 L11 (c ) RSIF  L22 (c )  (2-68)  And for K I , by definition we have:  K I  K II RSIF  2.9  (2-69)  Application of Cohesive Interface in Fracture Mechanics  In several cases where the crack front experiences a large scale processing zone (e.g., large blunting with large scale plasticity, fibrous rupture or ductile rupture), the stress singularity disappears from the crack-tip and traction forces emerge on the surface of the crack to resist the extensive increase of CTOD. This failure behaviour is common in FRP composites, especially when two adjacent layers have a bridged fiber orientation. When premature crack initiates the growth, a damage zone appears in front of the crack tip and dissipates the high stress intensity expected in LEFM. This damage zone occurs in form of a cohesive zone (e.g., fiber bridging) and hinders the identification of the crack-tip using conventional methods. In 47  such problems, a complimen ntary approaach to EPFM M based on ddamage mecchanics is reqquired to pro ovide accuraate results.  2.9.1  Cohesivee Zone Mod del  One of the advan nce methodss for simulatting the largge process zoone with tracction forces is the cohessive zone model (CZM). The basis of o this modeel is depictedd in the worrk of Dugdalle [75] and Barenblatt B [7 76], where th hey improveed the LEFM M by correctiing the plasttic zone in frront of the crack-tip c to include the defect proccess zone annd de-cohesiion. Dugdalle [75] assum med a finitee stress equaal to the yieeld stress of the materiaal, which waas contradictory to the crackopening stress in n brittle mateerials. Baren nblatt [76] fuurther investiigated the trraction existiing on c surfacce and linkeed the tracttion magnituude to the distance froom the cracck-tip. the crack Need dleman [83] implementted a potenttial functionn to characterize the ttraction-sepaaration modeel on the craack surface by considerring a cohessive zone forr ductile intterfaces in m metals. The potential p fun nction was deefined as folllows:            n ,  t    0 1  1    n     n     t2    ex  exp xp p  n0    n 0    t20   (2-70)  wherre 0 , n ,  t ,  n 0 and  t 0 are the material m fraccture energy,, normal annd tangentiall crack opening displaceements, norm mal and tang gential param meter relativve to the moddel. Based oon this poten ntial function n, with coheesive traction ns on the crrack surface,, normal andd shear tracttions (   n and a t ) are functions off crack open ning/sliding displacemennts. The cohhesive normaal and tangeential tractio ons can be extracted e by finding the derivative oof the potenntial functionn over the crrack opening g/sliding displacements as follows: 48  n     n ,  t   n  (2-71)  t     n ,  t   t  (2-72)  Tvergaard and Hutchinson [84] have extended Needleman’s work [83] on CZM to elasticplastic materials, while Cui and Wisnom [85] presented the application of this model for materials with perfectly plastic behavior. Later, Xu and Needleman [23] applied their CZM into the finite element analysis for the first time to simulate the dynamic fracture analysis. Camacho and Ortiz [24] utilized a linear degradable CZM to study the impact damage in brittle materials. Although presented models have been utilized in earlier researches to extract the behaviour of material interface layers, they have had disadvantages such as introducing softening and numerical instability to finite element models, especially during the crack propagation steps. To overcome these problems, a rigid cohesive model was proposed by Geubelle and Baylor [86]. In their method, a high initial stiffness, known as penalty stiffness, was applied to the interface elements for an undamaged region of the material. Material behaviour and the material degradation were assessed by damage indices to reduce the initial stiffness. A bilinear traction-separation law was also employed to model the initiation and propagation of delamination in composite plates. Despite reliable results acquired by this approach, the convergence became dependent on the penalty stiffness. Considering inappropriate penalty stiffness would result in numerical instabilities, especially when the material degradation is commenced in the elements of the front region. It is also worth adding that the above mentioned models are accurate for modeling the pure fracture modes, such as mode I, mode II and mode III, however, for the mixed-mode problems such as mode I/II, some contradictory issues are raised against the fundamental basis of CZM. For instance, 49  because of a large process zone at the crack front for mode I (normal opening), it is not expected to observe any stress singularity. This assumption may be opposed by the mode II (tangential sliding) acting on a same plane as mode I and causes a stress concentration due to the lack of cohesive tractions on the crack surface. To overcome such a difficulty, Ortiz and Pandolfi [87] employed the concept of effective separation,  eff , and effective traction,  eff , as follows: 2   eff   n2   *  t2  (2-73)  2   eff   n2   *  t2  (2-74)  where  * is a dimensionless correction factor relating the crack sliding to the crack opening. These two parameters can be simply related by employing a cohesive potential function as presented in Equation (2-70). Mi et al. [88] and Alfano and Crisfield [89] improved the mixed-mode failure model’s capabilities by developing a damage parameter relating the interfacial material strength to the crack relative displacements. The following damage index formulation in modeling materials degradation, progressive delamination and various other applications are presented by several researchers [25, 51, 90 and 52]: T  K Pen   0 (elastic part )  T  (1  Dmg ) K Pen  0     f ( softening part )   f  (decohision part ) T  0  (2-75)  where T, Dmg ,  ,  0 and  f are, respectively, the current effective traction/strength of interfacial material ( T   eff ), the damage index, relative opening/sliding displacement, relative  critical  opening/sliding  displacement  and  relative  failure  opening/sliding 50  displacement of the crack faces. The damage index for each individual failure mode can be defined as follows:  Dmg   f      0        0   f  (2-76)  The above model assumes the element degradation happens when the crack relative displacement exceeds the critical value, 0 , defined in terms of the penalty stiffness, KPen , and maximum interfacial strength, Tmax , as:  0   Tmax K Pen  (2-77)  According to Equation (2-77), selecting an appropriate value for KPen becomes important to establish a stable rigid cohesive model. While choosing a large value for KPen may help with true estimation of the elements stiffness before the crack initiation, it will reduce the required critical relative displacement value for the crack initiation and may cause numerical instability upon the crack initiation, known as the elastic-snap back [50]. Earlier research have been undertaken to formulate the KPen based on different types of material properties. Turon et al. [50] proposed a simple relationship between the transverse modulus of elasticity, Etransverse , specimen thickness, tck, and penalty stiffness, KPen :  K Pen   Etransverse t ck  (2-78)  where  was proposed to be equal to 50 to prevent the stiffness loss. A wide range of penalty stiffness values have been considered by other researchers for different traction-separation 51  laws, material properties and numerical simulation purposes. In essence, a comparison between numerical analysis and experimental result would be required to identify an optimum value for penalty stiffness in each given application. In addition to the penalty stiffness, the length of cohesive zone is another critical factor of the CZM. As opening or sliding displacements initially increase, elements in the cohesive zone gradually reach the maximum interfacial strength. Upon this point, the element’s stiffness moves into the softening region of traction-separation law and experiences irreversible degradation. The maximum length of cohesive zone occurs when the crack-tip elements are debonded completely (Figure 2-8).  x  T Tmax  T  T  T  x Cohesive zone length  T Tmax T  T  T  x Maximum cohesive zone length  Figure 2-8 Schematic of the cohesive zone in front of crack in a given step of numerical simulation  52  Accordingly, choosing the correct value for the length of cohesive zone is essential in numerical modeling of delamination to prevent numerical difficulties; such as a softening problem due to the implementation of the traction-separation law instead of a conventional constitutive relationship. Earlier works have studied this topic extensively. Hillerborg et al. [22] proposed a characteristic length parameter for isotropic materials as follows:  lch  Elongitudinal  GIC 2 Tmax  (2-79)  where lch , GIC , Elongitudinal are the cohesive zone length, the critical energy release rate and the longitudinal Young modulus of the material, respectively. For various traction-separation laws, Planas and Elices [91] introduced a different equation for isotropic materials. However, in orthotropic materials, like FRP composite laminates, Jin and Sun [92], and Yang et al. [93] demonstrated the effect of longitudinal, transverse and shearing moduli as well as laminate thickness on the cohesive zone length. They suggested a modified formulation for measuring the cohesive zone length in the slender composite laminates as follows: 1/ 4   G  lch   Etransverse 2IC  t ck3 / 4 Tmax    (2-80)  It is necessary to mention that the number of elements within the cohesive zone model is directly related to the cohesive zone length, and for realistic simulations, it is required to have a sufficient number of elements within this region. A range of different values has been proposed in earlier works [25 and 94], and it is clear that it is difficult to estimate an exact value that could work for all simulations.  53  2.10 Summary  In this chapter, elastic mechanical behaviour formulation of FRP composite materials was briefly described. Next, the expected failure modes for different FRP composite compositions were illustrated. General history of the fracture mechanics and damage mechanics were presented for a better understanding of subsequent steps in advanced failure modeling of FRP composite materials. Fundamentals of XFEM were introduced and its application and implementation in modeling LEFM and EPFM problems were presented. Also, CZM modeling of delamination interface was introduced and effective parameters in modeling the failure surface using this method were outlined. The following chapter will focus on implementation of principals of hybrid XFEM and CZM. It will be shown that this hybrid method will notably assist conventional FEM with numerical modeling of different fracture modes in composites.  54  3 Chapter: Extended Finite Element Method Implementation and Validation In today’s modern industries, numerical tools demonstrate a great capability to handle stress analysis problems with cumbersome structural shapes and material nonlinearity. Among different computational methods, FEM has gained a significant attention and is promoted with current technological and industrial needs. Numerous commercial FEM packages are available for performing advanced stress analysis for a wide range of linear and nonlinear materials. In the present research, the ABAQUS package has been implemented as it provides a large elements library as well as different material properties options for different types of analyses. It also contains different options for modeling contact surfaces and adaptive mesh analysis of structures. More importantly, a major interest of employing ABAQUS for modeling the delamination problem with XFEM was its flexibility for linking the userelement subroutines to the FEM solver [95]. In the present work, a new user-defined 3D element has been developed using Lagrangian formulation. Large deformation XFEM has been introduced to the element formulation with extended functionality to model the cohesive zone element properties (e.g., tractionseparation law) and interface contact (the code has been included in the Appendix). In the next sections, the underlying nonlinear FEM formulation and the extension of its application to XFEM modeling of CZM and interface contact in FRP composites will be discussed. The validation of the code using data updated in the literature will be presented.  55  Finally a set of sensitivity analysis on XFEM model parameters will be introduced for the first time.  3.1  Large Deformation Formulation  During FEM simulations, the reliability of results decreases when higher terms of deformation are neglected. Therefore, it is necessary to apply large deformation formulation in the present analysis especially as it is intended to evaluate the geometrical and constitutive material nonlinearities of FRP composites under excessive loadings. Regardless of the state of deformation, the equilibrium between the internal forces and external forces is always established [53]:  Pij X i   f jb  0  (3-1)  where Pij is the nominal stress components. The above differential equation is written in the reference configuration for the Lagrangian description. Based on the above equation, one needs to apply a small or large deformation formulation via the second Piola-Kirchhoff stress and Green strain, EGreen, tensors, and consequently use the constitutive material properties to formulate the normal stresses. In small displacement theory, only the linear term of Green strain tensor is being utilized to calculate the second Piola-Kirchhoff tensor via appropriate constitutive relations. However, in the large displacement theory, the nonlinear portion of Green strain tensor is considered in stress calculations. The Green strain tensor can be written as follows [53]:  56  EGreen  EL  ENL  (3-2)  where the linear part, EL , and the nonlinear part, ENL , are defined as:   u1   X   u  2    Y   u3    EL   Z u2 u1      X Y   u3 u2   Y  Z   u u   1 3  Z X   E NL  (3-3)   1  u1  2 1  u 2  2 1  u3  2             2  X  2  X  2  X    1  u1  2 1  u 2  2 1  u3  2            2 2   Y Y   2  Y      2  1  u1  2 1  u 2  2 1  u3               2  Z  2  Z  2  Z    u1 u1 u 2 u 2 u3 u3      X Y      X Y X Y   u1 u1  u 2 u 2  u3 u3   Y Z Y Z Y Z   u1 u1 u 2 u 2 u3 u3       Z X Z X Z X   (3-4)  To implement the above equations into FEM formulation, it is required to expand the equilibrium equation using the conventional FEM weak form approach. The obtained equilibrium can be written as [53]:   F  V  T     P dV    u k V  T     f b dV    u k t  T  f t dt  0  (3-5)  57  where P and F are the nominal stress tensor and the deformation gradient respectively. Applying the standard FEM Galerkin discretization process and rewriting Equation (3-5) in terms of nodal variables and FEM shape function leads to the following equation:  B  V  T  P dV   N T f b dV   N T f t dt  0 V  (3-6)  t  where B contains the Cartesian derivatives of the shape functions, Bij   N j X i  .  In order to maintain the virtual work principal and to preserve the constitutive material stressstrain relationship, we need to rewrite the above equation using the second order PiolaKirchhoff stress tensor, which is a symmetric stress tensor in contrast to the nominal stress tensor. The transformed form of Equation (3-6) can be expressed as follows [53]: T   (u k )   B  dV   N T f b dV   N T f t dt  0 V  V  (3-7)  t  where B can be defined using the deformation gradient vector, F, as:   x F   X  y X         BFB N  i  Y  N i  Z  N  i  Z  z X  x Y  N i x X X N i x Y Y N i x Z Z x N i  X X x N i  Y Y x N i  X X  y Y  x Y x Z x Z  z Y  N i Y N i Z N i Z  x Z  N i y X X N i y Y Y N i y Z Z y N i  X X y N i  Y Y y N i  X X  y Z  y Y y Z y Z  z  Z   N i Y N i Z N i Z  (3-8) N i z X X N i z Y Y N i z Z Z z N i  X X z N i  Y Y z N i  X X         z   Y  z  Z  z   Z   (3-9)  58  The variation of the discretized FEM Galerkin method with respect to du h is: T  d (u )   B d dV   B T dF  dV  K T du h V  V  (3-10)  where K T is the tangential stiffness matrix. If the nodal vectors are substituted into the above equation, the tangential stiffness matrix can be formed as [53]: T  KT  K Mat  K Geo   B C B dV   Gs M s Gs dV V  T  V  (3-11)  where KMat , KGeo , GS and M S are the material and geometrical portion of tangential stiffness, shape function derivatives matrix and re-arranged second Piola-Kirchhoff stress tensor, respectively. More specifically GS , and M S are defined as follows.  N i  X   0    0  N  i  Y G S   0   0   N i  Z   0   0   0 N i X 0 0 N i Y 0 0 N i Z 0   0   0   N i  X   0    0  N i   Y   0   0   N i  Z   (3-12)  59   11 I 33  12 I 33  13 I 33  M S    22 I 33  23 I 33   sym.  33 I 33   (3-13)  where I 33 is the unity matrix in the three dimensional domain. In the final loading step, based on the total nodal displacement, the Cauchy stress tensor,   , is calculated by transforming the second order Piola-Kirchhoff stress tensor [53]:     F F T det(F )  3.2  Nonlinear Solvers  (3-14)  The nonlinear analysis of structures normally requires an iterative solver to find nodal variables under the equilibrium condition in each step. In nonlinear FEM, purely incremental, known as explicit, and predictive/corrective, known as implicit, solvers have frequently been employed in the literature [96]. The first assumption in both solvers is the equilibrium of acting forces on the body. Thus, one can write the equilibrium of the body according to the internal forces, Fint , and the external loads, Fext :  Fint  Fext  0  (3-15)  In a case where the structure’s response is nonlinear, simply solving the first order linear equation will not satisfy Equation (3-15). The difference appearing in a nonlinear FEM between external loads and internal forces called residual forces, Rsi , in the ith loading step. Accordingly, Equation (3-15) is rewritten in the following format:  60  Fext  Fint  Rsi  (3-16)  In Figure 3-1, the relationship between internal forces, external forces and residual forces are depicted.  Load  Explicit Solver  Drift Error  Actual Structure Response  Displacement Figure 3-1 Explicit solver approach and possible drift error in nonlinear problems (dots show numerical solution steps) [96]  These residual forces correspond to the new structural configuration after experiencing the external load. In purely incremental/explicit solvers, no corrective procedure is applied to diminish the residuals. Hence, in such methods, small increments of external loading should be imposed to the structure to ensure the residual forces in each numerical step remain within an acceptance tolerance [96].  61  Load Fext Rs2  KT(2)  Rs1  KT(1)  Fint  K0  uh0 uh1 uh2  Displacement  Figure 3-2 Newton-Raphson iterative solver approach in nonlinear FEM problems [96]  On the other hand, in corrective/implicit solvers, the residual forces are moderated with an iterative correction method such as the Newton-Raphson or Quasi Newton-Raphson methods. In such solvers, the structure's tangent stiffness, K T , needs to be evaluated in each iteration to extract the displacement correction due to residual forces, which can be derived from the Taylor series as follows: Rs i 1 (u ih1 )  Rs i (u ih )   Rs i (u ih ) h (u i 1  u ih ) h u  (3-17)  Rsi (uih1 ) is the Jacobian Matrix presuming that: u h  KT   Rs i  ( Fext  Fint )  u h u h  (3-18)  Using Equation (3-17), displacement correction can be implemented by solving the following equation in each iteration.  KT(i ) uih  Rsi  (3-19)  62  After the ith iteration, the total displacement is the summation of previous iterations and the new residual forces. Tangent stiffness should be evaluated using the updated displacement (Figure 3-2).  uih  uih1  uih  (3-20)  In terms of convergence criteria, the iterative solver in ABAQUS stops the iteration steps based on two criteria. In the first one, if the residual forces can be negligible at every single degree of freedom in comparison to an overall residual force tolerance, external loads and internal forces are considered to be in equilibrium. The overall tolerance value can be set depending on the user’s demand, and if it remains intact, ABAQUS assumes 0.5% of the average force in the entire structure at the given time step. Another threshold for accepting the solution is based on the last displacement correction. The last correction should be relatively smaller than the fraction of the total incremental displacement (1% by default); otherwise, ABAQUS performs another iteration step [96].  3.3  Modeling Contact on Material Interfaces Using XFEM  In recent investigations, Khoei et al. [53] introduced a new modeling technique to simulate nonlinear 3D contact problems using the large deformation formula and XFEM. The proposed tangential stiffness matrix (merely for the interface material) based on nonlinear XFEM was defined as:  K T  K Mat  K Geo   B T DSep B dV   G s M S Gs dV T  V  V  (3-21)  63  where B and Gs have an enriched part added to the conventional FEM part of nodal vectors and DSep is the elastic-plastic constitutive matrix. These matrices can be redefined as follows:    B  Bu    ord  Gs  Gsu  Bb  ord  H  H  Gsb   H  (3-22)    (3-23) H  where B b and Gsb are defined as (H is the Heaviside function):  Bb  H   N i H  x   X X   N i H  x  Y Y   N i H  x   Z Z   N i H  x   N i H  x   X Y  Y X    N i H  x  N i H  x  Z Y  Y Z    N H  x  N H  x i i   X Z  Z X   N i H  y X X  N i H  y Y Y   N i H  y Z Z  N i H  y  N i H  y  X Y Y X  N i H  y  N i H  y  Z Y Y Z  N i H  y  N i H  y  Z X X Z   N i H  z   X X   N i H  z  Y Y   N i H  z   Z Z (3-24)  N i H  z  N i H  z    X Y  Y X  N i H  z   N i H  z   Z Y Y Z   N i H  z   N i H  z    Z X X Z   64  G sb  H    N i H   X   0    0   N H  i   Y    0   0    N i H   Z   0   0   0  N i H  X 0 0  N i H  Y 0 0  N i H  Z 0     0   N i H  X   0    0  N i H   Y   0   0   N i H  Z  0  (3-25)  As mentioned before, to satisfy the PUM fundamentals, the Heaviside function should be deducted by Heaviside value at each element’s node. The elastic-plastic constitutive matrix is defined as:   K 11 0 0    D   0 K 22 0   0 0 K 33    ep S  (3-26)  where K ii is the penalty (contact) stiffness assigned to the local coordinates of the contact surface. K 11 provides the impenetrable characteristic to the normal direction of the contact plane which follows the Kuhn-Tucker thresholds [73]:   n  0, PContact  0,  n   PContact   0  (3-27)  where PContact contains the vector of contact forces, respectively.  65  The remaining terms in the elastic-plastic constitutive matrix, K 22 and K 33 , create the friction forces and prevent the contact surfaces from abrupt sliding. For these terms, standard static and dynamic friction laws can be applied to perform the analysis [53].  3.4  Implementation of the Cohesive Zone Model  As described in Section 2.9, depending on the FRP composite material lay-up, we may expect to observe a large processing zone during decohesion. Hence, it is required to apply the cohesive crack modeling into the XFEM analysis. For this purpose, a bilinear tractionseparation law is utilized instead of the conventional constitutive relationship for the interface material to embed cohesive behaviour into the crack-tip region by means of re-arranging the nodal displacement components [35]. They also implemented a cohesive/contact transformation matrix ( BCoh ) to re-arrange the nodal degrees of freedom and rewrote crack opening and sliding displacements and tractions on crack faces as follows:    BCoh u k  (3-28)  T  D Interface   (3-29)  where D Interface includes the cohesive interface material properties as described in Equation (2-74). As proposed by Geubelle and Baylor [86], a rigid cohesive zone model is implemented to simulate the crack initiation. This model applies an initial rigid stiffness in enriched elements before damage initiation and provides a good interpretation of material deterioration while the relative crack displacement reaches the failure limit. Also, the application of this model with XFEM improves numerical simulations consistency and  66  reduces difficulties such as snap-back and results fluctuation which are mostly caused by the reduction of stiffness in fully damaged elements to zero [28].  T Tmax  Kpen  δ0  δf  δ  Figure 3-3 Bilinear traction-separation law for modeling the material degradation  The cohesive/contact transformation matrix can be extracted by finding the displacement in an enriched element. The displacement vector of any point in the enriched element, u g  x  , is:    u g ( x)   N i u iord     N j   N k H k  H j b Hj i j   k    N     0  0  u ord    N enr   b H   (3-30)  where,    N enr  N j  Nk H k  H j  j  k   (3-31)  The conventional FE shape function’s value remains constant for different points in the enriched element while the enriched shape function’s value demonstrates an odd function property with respect to the interface position: N enr bottom    N enr top   (3-32)  67  Thus, the global relative crack displacement,  , can be described in the form of the displacement difference between two points above and beneath the crack surface. u ord  N enr  H   N b     u k (top )  u k (bottom )  N   20 N  enr        u ord   N enr  H  b     u ord   H  b   (3-33)  (3-34)  In order to find the relative crack displacements in a global coordinate system, a simple transformation based on the normal and tangential directions, mij , of the crack plane with respect to the global coordinate can be employed:   m11   m21 m31  m12 m22 m32  m13   X   m11    m23   Y   2 m21 m31 m33   Z   m12 m22 m32  m13  u ord  m23  0 N enr  H   BCoh u k b  m33       (3-35)  Consequently, Equation (3-35) can be substituted into Equation (3-29) and used in the tangential stiffness formulation to introduce process zone properties within enriched elements (also compare to Equation 3-21): K T  K Mat  K Geo  K Coh   B  V  T  DSep B dV   G T M S G dV    ( BCoh ) D Interface BCoh dc T  V  (3-36)  c  Finally, in order to evaluate the internal forces, one can simply employ the Equation (3-37) as follows: T  T  Fint   B  dV   N enr f t dC V  ΓC  (3-37)  68  3.5  Other Num merical Imp plementation n Details  The numerical implementaation of th he presenteed approachh not onlyy requires good ntional FEM,, but also enntails the noddal selectionn for the enrriched, underrstanding off the conven enhan nced approaach for num merical integrration over a discontinuuous field aand contact bbound appro oximation which w will be described in n the next suub-sections.  3.5.1  Node Selection for Enrichment E t  In orrder to differrentiate betw ween ordinarry and enricched nodes, all elementss have to unndergo selecction criteria known as th he level-set technique. IIn the first sttep, the distaance of eachh node from the crack face should d be determiined. Then, in each eleement, nodaal distance vvalues shoulld be compaared. Any cracked c elem ment would contain a nnode with a positive distance valuee and a negaative one. Th he following g relationshipp interprets tthe first seleection criteriion for the crracked elem ments:  1maxx (element)  1min (elemeent)  0  (3-38)  wherre  1 is the element e nodaal distance value v from thhe crack facee and is definned as:             1 ( X * )  X 1*  X 1 n1  Y1*  Y1 n2  Z1*  Z1 n3 *  *  (3-39)  *  In Eq quation (3-3 39), ( X1 , Y1 , Z1 ) and (X1, Y1, Z1) are the eleement’s nodde coordinatees and arbitrrary point coordinates c on the crack surface. IIn the seconnd step, thee boundary oof the crack ked element should be identified i ussing another level-set teechnique. Insstead of the crack surfaace, the nodaal distance frrom the edgees of the craack plane is ccalculated annd nodes encclosed by all edges are considered. c The T selection n criterion c an be demonnstrated as fo follows: 69   1 min (element)  0  (3-40)  where  1 is equal to the element nodal distance from the edges of the crack plane defined as:               1 ( X * )  X 1*  X 2 (n2t n3  n3t n2 )  Y1*  Y2 (n3t n1  n1t n3 )  Z1*  Z 2 (n2t n1  n1t n2 )  (3-41)  In Equation (3-41), ( X 1 , Y1 , Z1 ) is an arbitrary point coordinate on the crack edge and ( n1t , n2t , n3t ) is the unit vector of the crack edge. As an example, in Figures 3-4 and 3-5, a meshed object was analyzed using the above mentioned level-set techniques with a crack plane situated in the middle layer of the mesh. These figures depict the level-set variables  1 and   1 , and the resulting enriched nodes as shown in Figure 3-6.  Figure 3-4 Variation of φ1 values in an example meshed object (square represents the positive value and circle denotes the negative value; the rectangle shows the crack plane)  70  Fig gure 3-5 Variiation of ψ1 va alues for the frront edge of th he crack planee in the examp ple meshed ob bject (square represen nts the positivee value and cirrcle denotes th he negative vaalue)  Figu ure 3-6 Accep pted nodes in the t example meshed m object based on Equ uations (3-39) and (3-41) of levelset criteria  3.5.2  Numericcal Integrattion of Disco ontinuous F Fields  In orrder to obtain n a weak fo orm integratiion formulattion (for FEM M) for the sstiffness andd force matriices, numeriical integratiion should be b performedd on discrettized elemennts. Methodss such  71  as Gauss quadrature and Simpson’s rule are well-known for their applications in numerical computations and provide accurate results in continuous fields. On the other hand, employing such numerical integration methods on discontinuous fields will not suffice and lead to pivot points in the corresponding equations’ system. Therefore, to prevent such problems from being introduced into fracture problems through cracked elements, auxiliary sub-triangles should be employed to discretize discontinuous material domains. Then, a numerical integration scheme can be utilized to evaluate the integration over each sub-triangle and consequently over the cracked element (Figure 3-7). For 3D models, a similar approach can be applied and sub-tetrahedral elements replace the sub-triangles in order to deliver sufficient integration points in discontinuous fields (Figure 3-8).  Crack  Figure 3-7 Sub-triangles of a 2D element with third order Gauss quadrature  72  Crack  Figure 3-8 Sub-tetrahedrals of a 3D element with third order Gauss quadrature  The numerical integration of a discontinuous field with subdivided elements is similar but not identical to the continuous case. For instance, the numerical integration of an arbitrary continuous function, f1 ( X , Y , Z ) , over a single element can be summarized as:   e  nG  f 1 ( X , Y , Z ) d  Ve  w g i f 1 ( 1i ,  2i ,  3i )  error i 1  (3-42)  where Ve is the element volume, e is the element domain; nG , wg and ( 1 , 2 , 3 ) are the Gauss quadrature order, points weights and local coordinates, respectively. In order to adapt the above formulation to a discontinuous function, f1' ( X , Y , Z ) , the following modification should be considered:   e  m poly nG  f ' ( X , Y , Z ) d  Ve   w p j w g i f ' ( 1i ,  2i ,  3i )  error j 1 i 1  (3-43)  73  wherre mpoly is th he number of o sub-polygo ons and w p is the modiification facctor for the w weight of eaach sub-poly ygon, which can be evalluated by thhe ratio of eaach sub-polyygon volumee, V * , over the element volume:  w  p  j  3.5.3    V j* Ve  (3-44)  Implemeentation of the t Integrattion Bound Approach  Anotther numericcal aspect off the interfacce contact aand cohesivee zone modeel simulationn with XFEM M is the inteegration (con ntact) bound d. Applicatioon of the inttegration bouund in the viicinity of th he interface can eliminaate the need for definingg new integgration pointts in the intterface surfaace and decreease the com mputational time. t It also prevents thee elastic-plasstic relationsship to be mistakenly m in ntroduced to the integrattion points w which are faar from the interface. D Despite the advantage off this method d with regulaar sub-elemeent methods,, the contactt bound dimeension ot explicit an nd one should perform sensitivity s aanalyses to oobtain an opptimum valuue in a is no given n application n. In Figure 3-9, integraation points in the entirre elements and sub-eleements and within w the integration bound is illustrrated by red stars and bllue dots.  74  Figurre 3-9 Integra ation points of a meshed ob bject with enriiched elements (red stars arre integration points and blue dotts are integrattion points witthin the integrration bound))  3.5.4  Simulatiion Algorith hm  A MA ATLAB cod de in concertt with ABAQ QUS was im mplemented tto perform tthe simulatioons. In the first f step, MATLAB M co ode reads the nodes c oordinates, connectivityy tables, material mech hanical properties, crack k surface an nd fracture innformation from a userr defined daata file (*.DA AT). Using the level set techniqu ue, describeed in Sectioon 3.4, the MATLAB code recog gnizes the ordinary elem ments and en nriched elem ments of thee model. It tthen assembbles an inputt file for AB BAQUS execcution. ABA AQUS utilizees the develloped user element subroutine (UEL L) along with h generated input file (*.INP) to runn the model aand write thee stress, strain and displacement fiellds in the reesults file. At A this stage,, MATLAB code reads these fieldss from d calculates the t energy reelease rate ussing the J-inntegral to evaaluate the staability the laatter file and of thee delaminatiion. If the deelamination is unstable, MATLAB rre-writes a nnew input file with exten nded delamin nation, otheerwise, it only increasess the grip oppening displlacement (w with no delam mination pro ogression). The T iterativee procedure continues uuntil the maxximum num mber of steps (defined by y the user priior to simulaation) is reacched. This nuumber shoulld be large ennough 75  to capture the response of the structure accurately, while being not too large to make the simulations excessively costly; in the present simulation cases it varied between 500 to 1000 steps. The summary of the above implementation procedure is shown as a flow-diagram in Figure 3-10.  76  Start the simulation  MATLAB reads the user input *.DAT file, performs the level-set technique and generates the ABAQUS *.INP file  ABAQUS runs *.INP using UEL subroutine, writes the displacement, strain and stress results on *.DAT file  MATLAB performs the post-processing, evaluates delamination stability, and extends the crack length and rewrites the user input *.DAT file.  No  Does Maximum Number of runs reached?  Yes  End of the simulation  Figure 3-10 MATLAB-ABAQUS simulation algorithm employed for modeling the delamination  77  3.6  Numericall Examples of Mode I and a II Fractture Tests: V Validation oof XFEM C Code  In th he next sections, severall benchmark k examples aare numericcally simulatted using the new XFEM M framewo ork presenteed in the previous p seections and the ABAQ QUS user-defined subro outine provid ded in the Appendix. A Th hese examplles are well ddescribed inn the literaturre and otherr researcherss have simiilarly perforrmed numeriical modelinng and experimental teests to verify fy their propo osed approacches and cod des as addresssed below.  3.6.1  Numericcal Simulatiion of the DCB D Tests  The Double Can ntilever Beaam (DCB) is one of thhe standard tests to evaaluate the m mode I a damage pproperties oof materials. DCB samplles are interllaminar fracture energy toughness and main nly fabricated d based on th he ASTM D5528-01 D [977]. Compositte materials considered in this exam mple are T3 300/977-2 carbon c fiber-reinforcedd epoxy and AS4/PEE EK carbon fiberreinfo forced poly ether ether ketone, wh hich are widdely used inn the aerosppace industrries to manu ufacture airfframe structtures and to o replace steeel componnents especiaally when a high service temperatu ure is requirred [52]. Th he T300/977 -2 models hhad a 150 m mm length, 220 mm h, and 1.98 mm thickn ness for each h arm with 55 mm iniitial crack (F Figure 3-11). For width AS4//PEEK, mod del dimension ns were 102 2 mm long, 225.4 mm widde and 1.56 m mm thick foor each arm, with a 32.9 9 mm initiall crack. Thee material prroperties forr each speciimen are givven in Tablee 3-1 [52].  78  Table 3-1 Mechanical properties of T300/977-2 and AS4/PEEK samples [52]  T300/977-2 Elastic Properties  E11 = 150 GPa E22 = E33 = 11 GPa G12 = G13= 6 GPa G23 = 3 GPa v12 = v13 = 0.25 v23 = 0.5  Fracture Properties  Tmax = 45 MPa GIC = 268 J/m2  AS4/PEEK Elastic Properties  E11 = 122.7 GPa E22 = E33 = 10.1 GPa G12 = G13 = 5.5 GPa G23 = 2.2 GPa v12 = v13 = 0.25 v23 = 0.48  Fracture Properties  Tmax = 80 MPa GIC = 969 J/m2  Previous works on both types of above composite samples have been reported using a cohesive interface layer method via the conventional FEM as well as the mesh-free method by Camanho et al. [25], Turon et al. [50], Barbieri and Meo [52]. In the present study, effects of different important modeling variables such as the interface stiffness (e.g., the penalty factor) and the cohesive region length are assessed via the XFEM model. Results (in the following sub-sections) are compared to the previous standard numerical approaches to provide further understanding of the advantages of the XFEM method in terms of numerical accuracy and stability.  3.6.1.1  Effects of Different Modeling Approaches  Turon et al. [50] investigated the effective cohesive zone length for T300/977-2 specimens. They suggested a cohesive zone length of 0.9 mm from numerical simulations on a very fine mesh (with element length, le, of 0.125 mm). Based on their work, the size of elements in the cohesive zone region should not exceed 0.5 mm, and a minimum of two elements are required in this region for acceptable modeling results. In Figure 3-11, the present XFEM global force-  79  displacement (F-Δ) results with the fine mesh simulation with the CZM penalty stiffness (Kpen) of 1×106 N/mm3 and the cohesive zone length of 1.5 mm are compared to the results from Camanho et al. [25], Turon et al. [50], Barbieri and Meo [52] by means of different numerical approaches.  80  w = 20 mm  L = 150 mm  tck= 3.96 mm  70 60  a0 = 55 mm  F (N)  50 40 30 Present 3D Nonlinear XFEM with the Fine Mesh (Kpen = 10e6 N/mm³) 2D Meshfree (Barberieand Meo [52])  20  Cohesive Zone Model (Turon et al. [50])  10  Decohesion Element Method (Camanho et al. [25]) Experiment (Camanho et al. [25])  0 0  2  4  6  8  10  12  14  16  18  Crack opening (mm)  Figure 3-11 A comparison between DCB test results via different methods on T300/977-2 samples  81  Figure 3-11 shows that all models predict a similar trend of load-displacement during delamination of the sample. The mesh-free method [52] overestimates the stiffness of the material and leads to raising the peak opening force by 5%, while the Turon et al. [50] cohesive finite element approach underestimates the resisting force by 10% in comparison to the experimental data obtained by Camanho et al. [25]. The XFEM estimates the peak opening force with 3% difference from the experiment and, similar to Camanho et al. [25], it provides a more conservative estimation of the fracture behavior of the DCB samples.  3.6.1.2  Effects of Mesh Size and Cohesive Zone Length  DCB tests of T300/977-2 specimens were simulated using two different mesh sizes, namely element lengths of 0.4 mm and 1.25 mm, to demonstrate the effect of coarse and fine mesh on the XFEM results. In addition, in order to illustrate the influence of cohesive zone length in each case, XFEM simulations were re-run (Figures 3-12 and 3-13) with different lch values and the fixed penalty stiffness of 1×106 N/mm3 [50].  82  w = 20 mm  L = 150 mm  tck= 3.96 mm  70  a0 = 55 mm  60  F (N)  50 40 30 Present 3D Nonlinear XFEM with the Fine Mesh (lch = 1.5mm)  20  Present 3D Nonlinear XFEM with the Fine Mesh (lch = 2.5mm) Present 3D Nonlinear XFEM with the Fine Mesh (lch = 3.5mm)  10  Experiment (Camanho et al. [25])  0 0  1  2  3  4  5  6  7  8  9  Crack opening (mm)  Figure 3-12 Load-Displacement DCB test results for the fine mesh (le = 0.4 mm) simulation with different cohesive zone lengths for T300/977-2 samples  83  w = 20 mm  L = 150 mm  tck = 3.96 mm  70 a0 = 55 mm  60  F (N)  50 40 30 Present 3D Nonlinear XFEM with the Coarse Mesh (lch = 1.5mm)  20  Present 3D Nonlinear XFEM with the Coarse Mesh (lch = 2.5mm) Present 3D Nonlinear XFEM with the Coarse Mesh (lch = 3.5mm)  10  Experiment (Camanho et al. [25])  0 0  1  2  3  4  5  Crack opening (mm)  6  7  8  9  Figure 3-13 Load-Displacement DCB test results for the coarse mesh (le = 1.25 mm) simulation with different cohesive zone lengths for T300/977-2 samples  84  Recalling Figure 3-12, in the fine mesh (le = 0.4 mm) model, it is observed that using 3 (lch = 1.5 mm) to 6 (lch = 2.5 mm) elements within the cohesive zone would lead to an accurate estimation of the experimental data, while increasing this critical value to 8 (lch = 3.5 mm) elements would introduce an unrealistic global softening behavior to the model. In the coarse mesh (le = 1.25 mm) runs (Figure 3-13), only in the case with 3 (lch = 3.5 mm) elements the simulation result became relatively agreeable with the experimental values. It is worth adding that in an earlier work by Harper and Hallett [51], they had also obtained load-displacement results using different mesh sizes in the interface elements. Namely, for smoother numerical results, they decreased the elements size to prevent the dynamic effects of larger elements failure such as a sudden drop of the fracture energy release rate. They also introduced a global damping factor of 5% into the simulations to dissipate the oscillation caused by the cohesive element debonding and the loss of stiffness in each step of crack propagation. In the present study, the enriched elements in the cohesive zone have the aggregation of stiffness from XFEM approximation and the traction-separation law. Hence, when complete debonding occurs, the affected elements’ stiffness does not completely disappear by elimination of the cohesive zone stiffness, and the XFEM approximation can inherently prevent the oscillations to a certain degree, without adapting a damping ratio into the model. This can be especially beneficial regarding computational time in the case of explicit analysis.  3.6.1.3  Effects of Different Penalty Stiffness Factors  As discussed in Section 3.4, the accuracy of the bilinear traction-separation law in modeling the process zone is directly dependent on the penalty stiffness value, whose optimum value may change from one material or fracture mode to another. In this section, in order to evaluate 85  the accuracy of XFEM predictions against different penalty stiffness values, a set of simulations with fine mesh were performed with a wide range of penalty stiffnesses, varying from 102 N/mm3 to 105 N/mm3, and a similar cohesive zone length (lch = 2.5 mm). Results were compared to the experimental data for both T300/977-2 and AS4/PEEK specimens, respectively (Figures 3-14 and 3-15).  86  w = 20 mm  L = 150 mm  70 tck = 3.96 mm  60  F (N)  50  a0 = 55 mm  40 Present 3D Nonlinear XFEM with the Fine Mesh (Kpen = 10e5 N/mm³)  30  Present 3D Nonlinear XFEM with the Fine Mesh (Kpen = 10e4 N/mm³)  20  Present 3D Nonlinear XFEM with the Fine Mesh (Kpen = 10e3 N/mm³)  10  Present 3D Nonlinear XFEM with the Fine Mesh (Kpen = 10e2 N/mm³) Experiment (Camanho et al. [25])  0 0  1  2  3  4  5  6  7  8  9  Crack Opening (mm)  Figure 3-14 The comparison between DCB test load-displacement results of T300/977-2 samples with different penalty stiffnesses  87  w = 25.4 mm  180  L = 102 mm  tck = 3.12 mm  160 140  a0 = 32.9 mm  F (N)  120 100  Present 3D Nonlinear XFEM with the Fine Mesh (Kpen = 10e5 N/mm³)  80  Present 3D Nonlinear XFEM with the Fine Mesh (Kpen = 10e4 N/mm³)  60  Present 3D Nonlinear XFEM with the Fine Mesh (Kpen = 10e3 N/mm³)  40  Present 3D Nonlinear XFEM with the Fine Mesh (Kpen = 10e2 N/mm³)  20  Experiment (Barbieri and Meo [52])  0 0  1  2  3  4  Crack opening (mm)  5  6  7  Figure 3-15 The comparison between DCB test load-displacement results of AS4/PEEK samples with different penalty stiffnesses  88  According to Figures 3-14 and 3-15, XFEM results are less sensitive to the larger order of penalty stiffness values (from 103 to 105 N/mm3) in comparison to the conventional finite element method discussed by Turon et al. [50]. Also, within the above recommended K Pen range, two sets of complimentary simulations on AS4/PEEK samples were run to see the effect of interaction between the mesh size and the penalty stiffness. According to the results in Figures 3-16 and 3-17, the mesh sensitivity decreases using lower values of the penalty stiffness, and vice versa. As AS4/PEEK has a higher critical energy release rate (Table 3-1), a larger cohesive zone region is expected in comparison to T300/977-2 samples and, hence, the sensitivity of simulations to the element size is reduced.  89  w = 25.4 mm  L = 102 mm  tck = 3.12 mm  180 160  a0 = 32.9 mm  140  F (N)  120 100 80  Present 3D Nonlinear XFEM with the Fine Mesh (Kpen = 10e5 N/mm³)  60  Present 3D Nonlinear XFEM with the Fine Mesh (Kpen = 10e4 N/mm³)  40  2D Meshfree (Barbieri and Meo [52])  20  Experiment (Barbieri and Meo [52])  0 0  1  2  3  4  5  6  7  8  Crack opening (mm)  Figure 3-16 The comparison between DCB test load-displacement results for fine mesh analysis of AS4/PEEK and previous works  90  w = 25.4 mm  L = 102 mm  tck = 3.12 mm  180 160  a0 = 32.9 mm  140  F (N)  120 100 80  Present 3D Nonlinear XFEM with the Coarse Mesh (Kpen = 10e5 N/mm³)  60  Present 3D Nonlinear XFEM with the Coarse Mesh (Kpen = 10e4 N/mm³)  40  2D Meshfree (Barbieri and Meo [52])  20  Experiment (Barbieri and Meo [52])  0 0  1  2  3  4  5  6  7  8  Crack opening (mm) Figure 3-17 The comparison between DCB test load-displacement results for coarse mesh analysis of AS4/PEEK and previous works  91  3.6.2  Numericcal Simulatiion of the ENF Tests  The End Notch h Flexure (E ENF) test was w also em mployed in this researrch to veriffy the effectiveness of the XFEM model in stu udying the m mode II fraccture properrties of com mposite materrials. The EN NF test conffiguration is similar to thhe three-poinnt bending/D DCB test, the only differrence is that the sample has a pre-assigned crackk in the midddle layer. EN NF samples ccan be prepaared accordiing to ASTM M standard D5528-01 [97]. Similaar to DCB tests, AS4/PEEK (Tablle 3-1) was utilized for ENF test siimulations. T The specimeen is 102 m mm long, 25..4 mm wide and 1.56 mm m thick for each e arm, wiith a 32.9 mm m initial craack. pared their mesh-free m m method resultts with the E ENF experim mental Barbieri and Meo [52] comp valuees. Based on their work, it was illustrrated that thhe mode II faailure has an abrupt natuure and perfo orming its nu umerical sim mulation req quires a relattively largerr cohesive zzone, compaared to DCB B mode, to capture c the entire e failuree in specimeens. In the section to folllow, we exxamine this effect e using the t XFEM model. m  3.6.2.1  Effectts of Cohesiv ve Zone Len ngth  r were considered tto model a ssudden modde II failure w within Diffeerent sizes of cohesive region one failure f step and a the resullts were com mpared with tthose by Barrbieri and M Meo [52] (Figgure 318).  92  w = 25.4 mm  L = 102 mm  900 tck = 3.12 mm  800 700  a0 = 32.9 mm  F (N)  600 500 400  Present 3D Nonlinear XFEM with the Fine Mesh (lch = 15 mm)  300  Present 3D Nonlinear XFEM with the Fine Mesh (lch = 25 mm)  200  Present 3D Nonlinear XFEM with the Fine Mesh (lch = 35 mm) 2D Meshfree (Barbieri and Meo [50])  100  Experiment (Barbieri and Meo [50])  0 0  1  2  3  4  5  6  Middle point deflection (mm)  Figure 3-18 The comparison between ENF test load-displacement results for AS4/PEEK and previous works  93  The results from XFEM illustrate smoother softening behaviour at the peak point of the loaddisplacement curve. The difference between the experimental curve and the numerical simulation after the peak is perhaps related to the slanted behaviour of the ENF loaddisplacement curve after crack opening and may be interpreted as an abrupt crack extension under increasing load under only one-step simulation. Also from Figure 3-18, it is evident that increasing the cohesive zone length leads to material softening and consequently decreases the peak value of the numerical load-displacement curve. Finally, comparing lch values in Figures 3-12 and 3-18 confirms that, for ENF simulation, much larger process zone lengths are required when compared to DCB simulations.  3.6.2.2  Effects of Different Penalty Stiffness Factors  In this section, different penalty stiffness values, KPen, were employed to evaluate the mode II response of AS4/PEEK samples. The length of the cohesive zone was kept constant, lch = 15 mm, and the crack propagation happened within numerous steps. A range between 103 N/mm3 to 106 N/mm3 was considered. The effect of this parameter on simulation results is illustrated in Figure 3-19.  94  w = 25.4 mm  L = 102 mm  1200 tck = 3.12 mm  1000 a0 = 32.9 mm  F (N)  800 Present 3D Nonlinear XFEM with the Fine Mesh (Kpen = 10e3N/m³) Present 3D Nonlinear XFEM with the Fine Mesh (Kpen = 10e4N/m³) Present 3D Nonlinear XFEM with the Fine Mesh (Kpen = 10e5N/m³) Present 3D Nonlinear XFEM with the Fine Mesh (Kpen = 10e6N/m³) 2D Meshfree (Barbieri and Meo [52])  600 400 200  Experiment (Barbieri and Meo [52])  0 0  1  2  3  4  5  6  Middle point deflection (mm)  Figure 3-19  Effect of penalty stiffness value on the ENF test load-displacement results for AS4/PEEK sample  95  It is clear from Figure 3-19 that applying penalty stiffness lower than 104 N/mm3 will result in extensive softening to the model which tremendously reduces the peak load. Such behaviour is derived from lower rigidity in the hardening region of the traction-separation law; however, such dependency of results on Kpen may become useful for adapting different material behaviours in different simulations. On the other hand, extensive hardening in the case with Kpen higher than 105 N/mm3 can be unrealistic as it reduces the real deflection of the specimen and prevents the actual failure crack opening to be modeled. Finally, comparing Figures 3-14 and 3-19, it is evident the ENF test is much more sensitive to the range of Kpen, however, a value of 105N/mm3 seems optimal for both DCB and ENF tests.  3.7  Summary  In this chapter, the nonlinear FEM formulation was summarized for delamination simulations and different solvers for extracting the structures response were introduced. The extension of the nonlinear formulation towards XFEM modeling was then introduced and numerical implementation of nonlinear XFEM in simulating the contact interfaces and CZM was demonstrated. Finally, two (DCB and ENF) example problems from the literature were simulated with the introduced XFEM method and results were compared with the conventional FEM, cohesive interface element and element free models. In the next chapter, the stochastic nature of FRP composite materials will be discussed and the randomness of fracture properties will be introduced to the DCB and ENF simulations of tested PPS/Glass material.  96  4 Chapter: Modeling Randomness Effect in UD Laminates Delamination: A Non-RVE Approach As discussed in Chapter 1, the multi-scale nature of composite materials is tied to the mixture of multiple constituent materials. More specifically, in such materials, comingling of the matrix and fibers results in heterogeneous material characteristics and often cumbersome procedures needed to analyze their mechanical properties at different scales. Studying mechanical properties of FRP composite materials can be classified into three different groups. Macro-scale is the largest scale for analyzing the FRP composites and includes coupon size models for experimental, analytical or numerical analyses. In such a scale, anisotropic or orthotropic material properties may be assigned to the material structure. The other extent of this type of analysis is related to micro-scale behaviours of FRP composites with a focus on the constitutive relationships of individual components of the composite as well as the interaction between the matrix and the fibers. The intermediate level of composite materials study is known as meso-scale which links the micro-scale analysis to the macroscale analysis. Investigation of a laminate lay-up with more concentration on individual plies’ mechanical and geometrical properties and orientation is an example of the meso-level study which can also include fracture phenomena such as matrix cracking and fiber breakage as well as delamination and fiber bridging between layers. For an accurate analysis of structural behaviour of composites, it is required to account for the different scales of material properties and include their effects in simulations. The most basic and yet tangible scale of material analysis in practice is perhaps the macroscale since it can be implemented to determine the effective properties of large composite structures. Most homogenization techniques aim at this level of analysis and employ smaller97  scale representative volume elements (RVE’s) to extract the averaged (effective) mechanical properties of the material at macro-level. The RVE should have specific characteristics; it has to be large enough to contain a sufficient number of heterogeneous characteristics of the composite under study and should assume a periodically distributed properties and boundary conditions in adjacent RVE’s. Such distribution assumptions through the entire structure can make the RVE homogenization technique vulnerable to existing defects in specific plies and/or the random distribution of fibers and their bridging. It also may monitor the damage mechanism and propagation of a crack, as these phenomena are relatively present in a local scale rather than a global scale. As a result, full scale (non-RVE) modeling is required for more effective damage modeling. The earliest attempts to assign statistical (random) properties to the crack location using stochastic modeling of laminates go back to the work of Wang et al. [44] and Fukunaga et al. [98]. More detailed investigation on the heterogeneous nature of laminates due to nonuniform fiber distribution was performed by Baxevanakis et al. [40] where they employed an image analysis technique and demonstrated the unreliable aspects of periodic and quasiperiodic assumptions of fiber distribution material, especially in the case of damage modeling of FRP composites. Further study on RVE assumptions by Bulsara et al. [99] suggested the efficient size effect of RVE’s. Trias et al. [43] compared stress and strain distributions of a periodic (RVE) model with a random model and concluded that a periodic assumption should be employed for extracting the effective properties of structures in global scale, while including the randomness is required for local analysis of structures such as matrix cracking and crack propagation. Silberschmidt [100] connected the microstructural randomness to the macro-level analysis and revealed that the fluctuation of mechanical properties is a result of 98  non-uniform fiber distribution. In the most recent investigation, Ashcroft et al. [48] applied a Weibull distribution to a set of cohesive elements and modeled the damage evolution of CFRP laminates. In their work, fracture energy toughness was defined statistically and exploited from a random distribution. The influence of the randomness of fiber distribution would increase in unidirectional (UD) FRP laminates as the possibility of fibers penetration within different adjacent layers of composite raises and leads to larger process zones in front of the crack. The response of the structure in such a case during the interlaminar crack growth will demonstrate a resistance while damage evolves and forms a resistance curve known as an R-curve. Fiber bridging results to an increase in fracture toughness of the material from initiation to steady-state delamination extension [101]. More studies on R-curves and fiber bridging showed that despite the fact that fiber bridging can directly affect the effective material properties and hence the shape of the R-curve, it can also be the case that the R-curve is affected by specimen geometry [102]. Nairn [103] combined the energy release rate of fracture mechanics with a CZM traction-separation law to represent the fiber bridging and extracted the R-curve of the model process zone. Airoldi and Dávila [104] applied experimental data and FEM to extract the cohesive element parameters and predicted the R-curve of delamination tests with fiber bridging. As also evident from their work, a high scatter of fiber bridging can exist in UD laminates and would have a significant effect on their macro-level performance. In the next sections of this chapter, the elastic mechanical properties of PPS/Glass will be assessed. R-curves of fracture experimental tests (DCB and ENF) with regards to uneven fiber/bridging distribution in the fabricated PPS/Glass UD laminates (Section 2.2) will be  99  studied first. Then, stochastic fracture properties of the material will be adapted into the developed XFEM model and used to capture non-repeatable material response.  4.1  Sample Preparation: Poly (phenylene Sulfide) (PPS)/Glass FRP  In recent decades the fiber reinforced thermoplastic composites have drawn more attention to some high-tech industries as compared to thermoset composites, as they are lighter, tougher, more sustainable, and more cost-effective with the right manufacturing process. Among different polymers used in thermoplastic composite industries, Poly Phenylene Sulfide (PPS), Poly Butylene Terephthalate (PBT) and Poly Ether Ether Ketone (PEEK) have demonstrated strong thermal and mechanical performances [105]. In the present work, PPS resin was chosen as a base material for manufacturing in-house test laminates. PPS is a semi-crystalline polymer and has excellent mechanical, thermal and physical properties. Its strength and affordability can help PPS to fill the gap between the partially crystalline industrial plastics and semi-crystalline resins such as PEEK. In addition, reinforced PPS with glass fibers offers economical advantages in comparison to carbon fibers by delivering high resistance characteristics against chemical and solvents, low moisture absorption, short thermoforming mould cycles with a greater creep strength and durability against temperature changes [105]. Typical applications of PPS/Glass composites in current industries vary from components in construction equipment, pumps parts and impellers to plastic parts in motor vehicles. In Figure 4-1, a microscopic image of UD tape of PPS/Glass, commercialized by TENCATE ADVANCED COMPOSITES [105], is depicted at two different regions of a sample. This particular product has an average void level less than 0.2% with a more even fiber-matrix distribution in the corner of the tape. 100  (a)  (b)  Figure 4-1 Microscopic images of fibers and matrix distribution of PPS/Glass UD tape: (a) Corner of the tape, and (b) Middle of the tape [105]  The common forming procedures used for making PPS/Glass laminates are mainly similar and their minor differences, apart from cost, lie on the range of required temperature or pressure as follows. 1- Press Lamination: In this process, PPS/Glass UD plies can be stacked in any desired orientation in a frame mould. Then, the frame is placed into a heated platen press where the assembly temperature is increased to approximately 330-360oC at a contact pressure until the PPS matrix melting point is reached. Namely, the platen pressure can be raised to 1.7 MPa and the assembly kept under this state for approximately 30 minutes. Then, the assembly fixture temperature is reduced by a cooling cycle flow passing through platens while the pressure is maintained. 2- Autoclave Lamination: Similar to the press lamination process, in this process PPS/Glass laminas are laid-up in any desired orientation and placed in a vacuum bag throughout the entire process. A high temperature resisting material is used for the vacuum bag. The fixture 101  assembly is placed in the autoclave for 30 minutes in a temperature close to 315-330oC. During this time, the pressure on the fixture is increased from ambient to 0.68 MPa. After the heating process, the part is cooled down to a temperature less than 93oC and the pressure can be dropped to ambient while the vacuum is released. In the present work, the press lamination method (using a Wabash 100 Ton Press shown in Figure 4-2) was utilized to form the required test laminates. In the assembly process, 14 layers of PPS/Glass UD lamina were stacked. Polyimide Teflon (0.147 mm) with high temperature performance (melting at 426oC) was placed in the middle of the stacked pile to represent the predefined crack in the specimens. The entire assembly was brought together in a preheated moulding fixture and then placed in a heated platen press where the assembled fixture’s temperature was sustained at approximately 350oC. The heating platen pressure was raised to 0.44 MPa and was kept in this state for approximately 28 minutes. This timing span let the inside of the mould reach the melting point of PPS at 350oC within 23 minutes and permitted the fibers to consolidate over 5 minutes in the molten matrix. Then, the assembly’s temperature was reduced via a cooling flow passing through platens for 20 minutes while the pressure was maintained at 0.44 MPa (Figure 4-2). It is worth adding that, 2D thermoplastic composite plies can also be formed into 3D complex shapes using the above apparatuses and appropriate moulds. For the PPS/Glass laminate, it can be stacked and heated to around 330-345oC using an infrared oven for a maximum of 8 minutes and then transferred to a predesigned core/cavity mould where it can be formed to a 3D shape under a pressure between 0.1-0.4 MPa. Depending on the part shape and dimensions, the production process may take under 10 minutes. This process is often referred to as thermoforming or stamping. 102  23  Cooling  Consolidation  20  Heating  Temperature (oC)  350  28  48  Time (minutes)  Figure 4-2 (a) Forming cycle used for preparing PPS/Glass test samples using (b) an automated press apparatus (Wabash MPI 100 ton)  4.2  Elastic Mechanical Properties of PPS/Glass FRP Composites  As mentioned earlier, the American Society for Testing and Materials (ASTM) standards [5557] can be employed to extract the material elastic constants. In the first step, the longitudinal and transverse elastic moduli of the PPS/Glass samples were measured based on ASTM D 3039/D 3039M [56]. The test specimens were prepared with 250 mm length, 20 mm width, and 3 mm thickness. After the test specimen was firmly aligned and tightened in the tensile machine grips, a transducer was mounted on the mid-span, mid-width of the specimen. Then, the loading cell started to load the specimen with the rate of 2 mm/min. Observed failure was explosive failure in the gage area at the middle of specimen (XGM) for longitudinal fibers, while in the transverse direction, lateral failure happened close to the grip at top or bottom  103  (LAT/LAB) of the specimen. Table 4-1 contains the test results for longitudinal and transverse directions.  Table 4-1 Mechanical properties extracted from tensile testing  Test repeat 1 2 3 4 5 6 7 Average Standard Deviation  Longitudinal (fiber) Direction  5.77 9.68 9.64 8.34 6.84 8.91 6.30 7.93  Modulus of Elasticity (MPa) 45326.07 43303.08 43601.28 50547.07 46821.49 39836.18 38975.62 44058.68  1.40  3460.23  Elongation (%)  Transverse (perpendicular to fiber) Direction Elongation Modulus of Elasticity (%) (MPa) 0.13 3490.58 0.06 1334.61 0.08 1176.69 0.09 1797.79 0.10 1513.28 0.09 1289.65 0.02 2100.46 0.08 1814.72 0.03  697.73  As we needed to acquire the remaining material properties for subsequent numerical simulations, the average Young’s modulus extracted from experimental data above were compared to the supplied material data sheets and brochures (see Figure 4-3 for a snapshot) to find and adjust the closest set of material properties given the actual fiber volume fraction, reinforcement architecture, etc (for confidentiality reasons the material’s full specification and composition have not been disclosed). The final set of elastic properties implemented in the subsequent simulations of this thesis is given in Table 4-2.  104  Figure 4-3 Snapshot of elastic mechanical properties for a typical woven PPS/Glass ply from material data sheets  Table 4-2 Final set of elastic properties of manufactured UD PPS/Glass FRP composites with “1” being the fibers direction; “2” and “3” are perpendicular directions to fibers.  Effective Elastic Properties  4.3  E11 = 44,400 MPa  G12 = 880 MPa  v12 = 0.25  E22 = 1800 MPa  G23 = 660 MPa  v23 = 0.48  E33 = 1800 MPa  G13 = 880 MPa  v13 = 0.25  Fracture Tests on the Fabricated PPS/Glass Composites  Both DCB and ENF tests were conducted on the PPS/Glass UD specimens according to the procedure described in ASTM D5528-01 [97]. The main goal of these tests was to measure the mode I and II fracture energy toughness properties of the material. In DCB test, specimens with pre-inserted delamination were put in the tensile machine and underwent opening displacement with the rate of 2 mm/min on the grips. At the onset of delamination extension, 105  the force fo on the loading celll was record ded while deelamination iis was perm mitted to proppagate for 5 mm, obserrved and co ontrolled vissually, befoore unloadinng the specimen for thee next loadiing step with h extended delamination d n length. In the case of ENF test, a similar speecimen confiiguration with pre-inseerted delamination wass employedd while the test set-upp was coinccided with th he three poin nt bending teest. The loadding on the m mid-span waas increased with a rate of o 2 mm/miin until the delamination commencced in the sppecimen (Figgure 4-4). A At this pointt, the force was w recorded d and used in n energy releease rate calcculation. Thee calculationn steps via reecorded dataa from DCB and ENF tessts have beenn described in Appendixx B.  (a)  (b)  4 Experimen ntal test set-up ps: (a) DCB, an nd (b) ENF Figure 4-4  4.3.1  DCB and ENF Testt Results  In order to measu ure the modee I fracture energy e toughhness, three different caalculation meethods o the DCB B test resultss with three repeats: Com mpliance Caalibration Method, were employed over Modiified Beam Theory T Meth hod and Mo odified Comppliance Caliibration Metthod. For thee ENF 106  test, based on the nature of the mode II failure, crack propagation happens due to excessive shear deformation and has an abrupt nature as addressed in Chapter 3. Such behaviour makes it difficult to control the external load to achieve a target crack extension. Due to this test limitation, in the present study only the first step of crack propagation in ENF tests was considered while the test repeats were carried out with the same pre-assigned crack length (43 mm). Figures 4-5 and 4-6 depict the obtained experimental data from three repeats of the DCB test. According to the observed increasing trend between crack length and the critical energy release rate, it could be concluded that the material experiences fiber bridging during delamination. This was further verified by macro- and micro- imaging of samples (Figure 47). It is worth adding that fiber bridging is more likely to occur in unidirectional laminates than in woven fabric composites, since the layers in UD’s are laid-up in a single orientation and during the moulding phase, as the matrix melts, there is less geometrical confinement to prevent fibers from penetrating into adjacent layers. The clear data scatter through the test repeats in Figure 4-5 also shows that, despite the DCB sample coupons have been cut from the same large (master) plate, the distribution of effective material properties could randomly vary from one coupon to another. This observation may be due to a non-uniform reinforcement pattern, uneven pressure or heating/cooling effect during the compression moulding stage of the master plate. Even within the same plate, more uniform pressure or heat concentration in the middle section of the plate may have been present compared to its edge sections (also recall Figure 4-1).  107  w  L  tck  Specimen No. 1 Specimen No. 2 Specimen No. 3  0.8  0.7  G (mJ/mm2)  G (mJ/mm2)  0.7  Specimen No. 1 Specimen No. 2 Specimen No. 3  a0  0.8  0.6 0.5 0.4 0.3  0.6 0.5 0.4 0.3  0.2 40  50  60  70  80  0.2  90  40  50  Crack Length (mm)  60  80  90  Crack Length (mm)  (a)  (b) Specimen No. 1 Specimen No. 2 Specimen No. 3  0.8 0.7  G (mJ/mm2)  70  0.6 0.5 0.4 0.3 0.2 40  50  60  70  80  90  Crack Length (mm) (c) Figure 4-5 The variation of fracture energy toughness versus crack length for the three tested samples using: (a) Compliance Calibration Method, (b) Modified Beam Theory Method, (c) Modified Compliance Calibration Method [97]  108  Matrix cracking  Ω  Crack tip extends; Matrix cracking  Ω  Further Loading  Γ Γ Fiber bridging  Further Loading  Fiber bridging zone extends  Delamination continues; Matrix cracking  Ω  Γ  Some earlier fiber bridges are broken (d)  Figure 4-6 Illustration of the fiber bridging zone (FBZ) during crack propagation; as the crack length increases, FBZ emerges in the cracked region up to the fiber’s rupturing displacement; after fiber breakage, the FBZ effect vanishes from the region which has exceeded the failure opening displacement  According to ASTM 5528-01 [97], the modified beam theory provides more conservative result in comparison to other methods, and hence it is recommended for design purposes. Finally, in addition to the ASTM D5528-01 [97] discontinuous DCB test, continuous crackopening tests, i.e., without unloading and re-loading for the next initial crack length configuration, were performed to monitor the effect of progressive failure in the specimens (results to be discussed in Section 4.5).  109  5 mm  (a)  (b)  (c)  Figure 4-7 Different images of a DCB test sample: (a) macro scale image of fiber bridging, (b) X-ray micro-tomography image of fiber bridging along the sample thickness, and (c) attenuation of the X-ray reflection due to absorption; demonstrating uneven distribution of fibers  For ENF tests, as mentioned earlier capturing the step-by-step crack propagation was not feasible, therefore only the initial step of crack propagation was used for the analysis (Figure 4-8). Similar to the DCB test, continuous loading was also employed for the ENF test to observe the effect of fiber bridging with extensive crack evolution in the test coupons (results 110  to be discussed in Section 4.5). The summary of fracture properties extracted from DCB and ENF tests is shown in Table 4-3.  1.7 1.6  Load (N)  Fracture Energy Toughness (N.mm/mm2)  1.8  1.5 1.4 1.3 1.2 6.0  6.5 7.0 Critical Opening Displacement (mm)  7.5  400 390 380 370 360 350 340 330 320 310 300 6.0  (a)  6.5 7.0 Critical Opening Displacement (mm)  7.5  (b)  Figure 4-8 ENF test repeat results with a constant crack length (43 mm): (a) the variation of fracture energy toughness versus the mid-span displacement, and (b) the variation of bending load versus the midspan displacement.  111  Table 4-3 Fracture properties of PPS/Glass samples extracted from DCB and ENF tests  Fracture Properties  4.4  Mode I (DCB)  Mode II (ENF)  Tmax = 9.6 MPa  Tmax = 5.5 MPa  GIC(ave) = 0.48 kJ/m2  GIIC(ave) = 1.48 kJ/m2  As a function of delamination  As a function of delamination  length:  length:  GIC = 0.0075acr + 0.04 kJ/m2  GIIC = 0.045acr - 0.5 kJ/m2  (For acr> 40 mm)  (For acr> 40 mm)   0 = 9.6e-05 mm   0 = 5.5e-05 mm  Stochastic Fracture Properties  As addressed in Section 4.1, FRP composite materials demonstrate a large amount of randomness in material properties due to an uneven distribution of fibers in the matrix, and the possible penetration of fibers within different layers of laminate during the forming processes (especially in the case of unidirectional laminates). In the present work, fracture properties of the tested samples are considered to have a stochastic nature based on nonrepeatability observed in repeated experiments, as opposed to deterministic approaches where averaged values of experimental results are assumed for the estimation of material properties. Namely, a random number within the range of experimentally measured fracture energy 112  toughness values was picked to form a stochastic bilinear traction-separation law of enriched elements in the cohesive zone (Figure 4-9): GC  GC ( ave )   1  Rand 1  Rand 2  GC ( std )  (4-1)  where Rand1, GC(ave) and GC(std) are the random integer (odd or even to assign a random sign), average and standard deviations of fracture energy toughness, respectively. .  T Tmax  GCL = GC(ave) – Rand2 × GC(std) GCH = GC(ave) + Rand2 × GC(std)  Kpen GCL<GC<GCH δ0  δfL  δf  δfH  δ  2GCL/Tmax < δf < 2GCH/Tmax Figure 4-9 Proposed stochastic bilinear traction-separation behaviour (Rand2 is a random number taken from a 2-parameter Weibull distribution; GCL and GCH correspond to the lower and upper limits of GC via Equation 4-1)  Following the fiber bridging discussion in Section 4.3 and the observed experimental trend in Figure 4-1, the fracture energy toughness distribution was considered to be a function of the crack length and hence, a linear interpolation was utilized to extract GC(ave) for each specific crack length. For GC(std) , it can be a constant or in a more general form it can scale with GC(ave), which in turn becomes a function of crack length. The second random number (Rand2) 113  corresponding to individual enriched elements in each stage of damage evolution was taken from a uniform distribution with a range of 0 to 1. The randomly selected values were then converted to a Weibull two-parameter distribution between 0 and 1 via: 1  1  1 Rand2 weibull   ln1  Rand2 uniform   1   (4-2)  where 1 > 0 is a shape parameter and 1 > 0 is the scale parameter of distribution and both are considered to be equal to 3. It should be added that, according to conventional bilinear traction-separation behaviour, a direct relationship exists between the critical fracture energy toughness, GC , failure crack opening displacement,  f , and maximum interface strength,  Tmax : GC   Tmax f 2  (4-3)  Therefore, the obtained statistical distribution of the fracture energy toughness can be converted into the variation of failure crack opening and/or maximum interface strength of material via Equation (4-1). Khokhar et al. [106] introduced the randomness into their simulation by implementing a relationship between random fracture energy toughness and the maximum interface strength by keeping the failure crack opening displacement constant. To improve the numerical simulation’s convergence, the present study assumes a constant value for the maximum interface strength (see Table 4-3) while the failure crack opening displacement is randomly varied during damage evolution (see Figure 4-4). This approach relies on constant penalty stiffness and prevents the over-strengthening of elements’ stiffness in the process zone. It also stands with the fact that the crack length extension in test  114  specimens is a function of the crack-tip opening displacement (CTOD) and the energy release rate in front of crack-tip.  4.5  Numerical Results and Discussions  An XFEM model of the PPS/Glass composite samples under DCB tests were established using the developed ABAQUS user-defined element subroutine. To consider the stochastic aspect of fracture properties in conjunction with Equation (4-1), two different approaches were employed. In the first approach, the standard deviation of fracture energy toughness was assumed to be constant during the crack propagation (i.e., equal to the standard deviation of the entire DCB experimental points in Figure (4-5). In the second approach, a linear function was assigned using test data to relate the fracture energy toughness standard deviation to the crack length (Table 4-4).  Table 4-4 Employed standard deviation schemes in stochastic simulations of DCB test  Standard Deviation Constant  Function of delamination size  GC(std) = 0.121 kJ/m2  GC(std) = 0.2125 - 0.0016acr kJ/m2  (For acr > 40 mm)  (For acr > 40 mm)  Figure 4-9 shows a comparison of opening force between measured and predicted values via the above two approaches under the standard (discontinuous) DCB tests [97]. 115  Numerical Specimen No. 1 Specimen No. 2 Specimen No. 3 Estimated Experimental Average  50 45  45 40  F (N)  F (N)  40  Numerical Specimen No. 1 Specimen No. 2 Specimen No. 3 Estimated Experimental Average  50  35  35  30  30  25  25  20 12  14  16  18  20  22  24  26  28  20  30  12  Critical crack opening (mm)  14  16  (a)  18  20  22  24  26  Critical crack opening (mm)  28  30  (b)  Figure 4-10 Comparison of the opening force in stochastic simulations of DCB tests with experimental data using: (a) constant standard deviation formulation, and (b) standard deviation as a function of crack  0.8  0.8  0.7  0.7  0.6  0.6 G (N/mm)  G (N\mm)  length  0.5 0.4  Numerical  0.3  Specimen No. 1  0.2  Specimen No. 2  Linear (Estimated Average)  50  60 70 Crack Length (mm)  80  Numerical Specimen No. 1 Specimen No. 2 Specimen No. 3 Linear (Estimated Average)  0.2  0 40  0.4 0.3  Specimen No. 3  0.1  0.5  90  0.1 0 40  (a)  50  60 70 Crack Length (mm)  80  90  (b)  Figure 4-11 Comparison of predicted fracture energy toughness via stochastic simulations of DCB tests with experimental data using: (a) constant standard deviation formulation, and (b) standard deviation as a function of crack length  116  40  35  35  30  30  25  25  F (N)  F (N)  40  20  20 15  15 Deterministic Simulation  10  Experimental Result No. 1  5  Experimental Result No. 2  5  Deterministic Simulation  10  Experimental Result No. 1  Experimental Result No. 2  0  0  0  0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 Crack Opening (mm)  2  4  6  (b)  40  40  35  35  30  30  25  25  F (N)  F (N)  (a)  20 Stochastic Simulation No. 1 Stochastic Simulation No. 2 Stochastic Simulation No. 3 Stochastic Simulation No. 4 Stochastic Simulation No. 5 Stochastic Simulation No. 6 Experimental Result No. 1 Experimental Result No. 2  15 10 5 0 0  2  4  6  8 10 12 14 16 18 20 22 24 26 28 30 Crack Opening (mm)  8 10 12 14 16 18 20 22 24 26 28 30 Crack opening (mm)  20 Stochastic Simulation No. 1 Stochastic Simulation No. 2 Stochastic Simulation No. 3 Stochastic Simulation No. 4 Stochastic Simulation No. 5 Stochastic Simulation No. 6 Experimental Result No. 1 Experimental Result No. 2  15 10 5 0 0  2  4  6  (c)  8 10 12 14 16 18 20 22 24 26 28 30 Crack opening (mm)  (d)  Figure 4-12 Comparison of measured opening force with predicted values in stochastic and deterministic simulations of continuous DCB test using: (a) fracture energy/toughness remains equal to the average value of experiments, (b) fracture toughness only changes with increase in delamination length, (c) fracture toughness increases with extension of delamination with constant standard deviation formulation, and (d) fracture toughness increases with extension of delamination with standard deviation as a function of delamination length  117  The stochastic fracture energy toughness values (via the XFEM force values in Figure 4-10 and the standard formulas of ASTM D5528-01 [97]) were also calculated and compared to the measured values in Figure 4-11. Similarly, the comparison of results under the continuous test mode, i.e. progressive damage, is depicted in Figure 4-12. According to Figure 4-11, prediction results for the discrete (standard) DCB case with constant standard deviation tend to follow experimental results that are closer to the average data points. This observation has been directly reflected in the continuous-mode tests in Figure 4-12. As depicted in Figure 4-12(a), the XFEM simulation with constant fracture toughness and no stochastic effect (i.e., similar to a deterministic simulation) tends to follow the lower bond of experimental repeats. Choosing the lower bond fracture limit can be quite acceptable in practical applications where the safety has the highest impact on design. In Figure 4-12(b) the XFEM simulation with fracture toughness as a function of delamination length, and with no standard deviation effect, shows an increase in resisting force during crack evolution, and a better overall prediction capability. This response in particular demonstrates that the implemented XFEM model in macro-level simulations can capture the effect of fiber bridging in the meso/micro level by increasing the fracture energy toughness during each step of delamination extension. Studying stochastic prediction cases, Figures 412(c) and (d), it is seen that having the standard deviation varying with the crack length can reproduce a more realistic (wider) range of results by means of a higher variation introduced to fracture properties. It can also be observed from Figures 4-12(c) and (d) that utilizing a constant standard deviation reduces the fluctuation of stochastic simulations and demonstrates a smoother trend, while the variable standard deviation method leads to larger fluctuations during the crack propagation steps. Interestingly, the highest scatter/fluctuation at the opening 118  displacement of 22 mm in Figures 4-12(d) coincides with the high deviation observed in mode I fracture experiments in Figure 4-11(b) when the delamination length is close to 65 mm. The numerical results oscillation in Figure 4-12 would have been caused at each crack propagation stage when the energy release rate reaches the critical fracture energy toughness values. More specifically, when CTODs reaches its onset (i.e., crack starts opening), the corresponding elements in front of the crack are progressed within the rigid (high stiffness) portion of the bilinear traction-separation law up to the apex where the maximum interface strength is reached. After this point, the material faces sudden softening and the opening force reduces to the complete failure of elements, leading to the propagation of the crack into the next element where again a local increase of the global opening load is expected. Also, the delamination propagation inherits different size of kinks/jumps which links to the fact that we are allowing the failure opening displacement to vary for different enriched elements. This assumption causes the numerical model to deviate from a smooth trend observed in experimental results. Figures 4-13 and 4-14 show the schematic of crack evolution based on the bilinear traction-separation law and the XFEM model at different stages of delamination.  119  Onset of rigid hardening:  T  T  δ  δ  Apex of the bilinear  T δ  Further loading  traction-separation law:  T  T  T δ  δ  δ  Further loading  Softening Initiation of cohesive zone:  T  T  Deterioration of  T δ  δ  δ  cohesive zone  Further loading  (progression of damage to next step):  T  T -  δ  T δ  δ  Figure 4-13 Evolution of the cohesive zone in front of crack upon loading in a given simulation step  120  (a)  (b)  (c) Figure 4-14 Stages of delamination propagation within the DCB numerical model: (a) Onset of rigid hardening in the process zone, (b) Apex of the bilinear traction-separation law, and (c) Deterioration of the cohesive stiffness  121  In the case of ENF tests, the experimental set-up was simulated using the same ABAQUS user-element subroutine. Randomness was introduced into the analysis only by means of the constant standard deviation method. The results for continuous delamination are depicted in Figure 4-15.  Stochastic Simulation No. 1  450  Stochastic Simulation No. 2 Stochastic Simulation No. 3  400  Stochastic Simulation No. 4  350  Experimental Results No. 1 Experimental Results No. 2  F (N)  300  Experimental Results No. 3 Experimental Results No. 4  250  Experimental Results No. 5  200 150 100 50 0 0  1  2  3  4  5  6  7  8  9  10  Middle point deflection (mm) Figure 4-15 Comparison of stochastic measured and predicted force-displacement values in ENF tests on the PPS/Glass samples  As illustrated in Figure 4-15, the stochastic simulations of the ENF test have resulted in a great agreement with non-repeatable experimental data. In performed simulations, critical (opening) deflection varied due to the traction-separation law’s dependency on the failure crack sliding displacement and, as illustrated, it also affects the variation of the maximum critical flexural load. Based on the observed trend in simulation data, it can be concluded that 122  the fiber bridging effect in ENF test has a minimal effect in global behaviour of the results. Figure 4-16 shows the XFEM model contours under different stages of delamination in ENF test. It should be added that from an application perspective, for forming processes of composite preparation, such as compression moulding set-up in Figure 1-10, mode II fracture would be more relevant due to sliding between layers of the laminate under the punch load.  123  (a)  (b)  (c) Figure 4-16 Stages of delamination propagation within the ENF numerical model: (a) Apex of rigid hardening in the process zone, (b) Initial stage of crack propagation, and (c) Extensive deterioration of material  124  4.6  Summary  In this chapter, different methods of fabricating PPS/Glass test samples were reviewed and their elastic mechanical properties were extracted. The stochastic nature of DCB and ENF tests on PPS/Glass composite was illustrated. This characteristic was next introduced into the XFEM numerical analysis by means of a random fracture energy toughness distribution. Failure crack opening/sliding displacements were employed to relate the fracture energy toughness randomness to the traction-separation law. In order to demonstrate the capability of the present method in capturing the DCB and ENF test results, several stochastic simulations were performed and the results were compared to experimental data.  125  5 Chapter: Conclusions and Future Work Recommendation 5.1  XFEM Model Development  A new framework was presented to numerically simulate the fracture behaviour of FRP composite materials and, more specifically, the unidirectional PPS/Glass laminates. For this purpose, the ABAQUS finite element package was utilized as a simulation engine. In order to extend the capability of ABAQUS in modeling crack and delamination contact interfaces in large delaminations, a user-element subroutine was developed to introduce nonlinear XFEM element properties including CZM and contact. In these elements, the right hand side vector and the stiffness matrix were defined and, according to the employed degrees of freedom, were assembled into the global system of equations. In addition, the stochastic fracture properties of the composite samples were adapted into the code to capture the randomness seen through non-repeatable test results. The model may be used with both implicit and explicit nonlinear solvers for both deterministic and stochastic simulations.  5.2  Performed Deterministic Simulations  Following earlier works in the literature, the performed benchmark deterministic simulations demonstrated the effectiveness of the combined XFEM - cohesive zone model (CZM) and contact interface modeling approach in 3D numerical analysis of mode I and II fracture mechanics of fiber reinforced composites in the presence of large deformation effects and interface material nonlinearity. Sets of sensitivity analysis were also performed to evaluate the effect of modeling parameters on the XFEM numerical predictions. For more reliable simulations, it was found that a minimum of two elements is required within the cohesive zone region (regardless of critical length) in front of the crack tip. On the other hand, 126  considering a very long cohesive zone would introduce a global softening of material into simulations and can lead to an underestimation of the peak opening force. A maximum of six elements with a fine mesh is recommended as the limit within the cohesive zone region for mode I fracture analysis of the studied unidirectional composites. It was also observed that reducing the penalty stiffness value in the traction-separation law improves the convergence of numerical simulations and reduces the mesh size sensitivity. However, using conventional FEM this can again cause a softening problem and reduce the peak opening force. The XFEM approach with embedded CZM was found to be less sensitive to the aforementioned effects, particularly when the penalty stiffness value is chosen arbitrarily within the range of transverse and longitudinal moduli of the composite. In the case of mode II fracture analysis, an increase in the penalty stiffness could cause extensive flexural deformation without any crack formation while a small value for the penalty stiffness may lead to a lower critical flexural force.  5.3  Performed Stochastic Simulations  For stochastic analysis, the prediction of fracture behaviour of fabricated unidirectional PPS/Glass composites was presented. DCB tests, both in the standard and continuous modes, were conducted for extracting the experimental material fracture properties along with their variation through test non-repeatabilities. Based on the experimental data and obtained x-ray images, it was concluded that fiber bridging was present in the specimens during delamination and the energy release rate can be a function of the crack length. In order to reproduce the experimental results with numerical simulations, under the stochastic behaviour of the material the bilinear traction-separation cohesive behaviour was applied within the framework 127  of nonlinear XFEM. It was found that the approach is capable of predicting delamination surfaces with the traction due to possible fiber bridging effects. To take the present stochastic material properties effect into account, even under a given crack length, the fracture energy toughness value was randomized using two different standard deviation methods. Specifically, applying the constant standard deviation method (i.e., using the overall standard deviation of the entire data sets) demonstrated a low variation of predicted/model values, but induced fewer numerical fluctuations in the continuous test mode. In contrast, considering a standard deviation as a function of crack length captured a wider range of experimental data points by increasing randomness effects in the fracture properties, though it induced some fluctuations in the numerical curves. For the ENF test, performed to study the mode II fracture energy toughness of the same PPS/Glass samples, due to the abrupt nature of shear failure, capturing the dependency of the mode II fracture energy toughness on the crack length was not possible. Such failure behaviour, however, demonstrated a lower randomness as represented by more repeatable continuous mode test data. Therefore, the average fracture energy toughness in this mode was assumed to be independent of the crack length and was only affected by the standard deviation of the test results. The numerical simulation results were well-agreeing with the experimental ones and the variation of the failure deflection was captured by the stochastic XEM model.  5.4  Potential Future Work  Future works may be defined based on three different areas. First, improvements can be made to the FEM shell analysis by introducing a similar nonlinear XFEM into the standard shell 128  element formulation, employing the mass and damping matrices for dynamic analysis of fracture tests and applying different traction-separation laws for modeling different material types. Second, both DCB and ENF tests can be repeated with larger size specimens to capture the stabilizing regions of R-curves observed in experiments with longer delamination sizes. Such tests can help to scrutinize the stable and unstable nonlinear regions of fiber bridging. Third, other fracture tests such as the end loaded split (ELS), mixed-mode bending (MMB) and edge crack torsion (ECT) can be studied to verify the advantages of the XFEM method in modeling more complex fracture tests and potentially improving the traction-separation law application for such tests. Another significant aspect of stochastic modeling is to understand the effect of number of random simulations. A larger number of simulations would produce a more distinguished region of predictions between upper bond and lower bonds of data. Finally, more focus can be put on increasing the number of experimental repeats to gain a better estimation of deviation of data from mean values and hence to improve the reliability of the simulations. This can also include the introduction of a more precise random distribution of fracture properties, specifically by including a non-liner dependency of the fracture energy toughness, mean and standard deviation on the crack length via further experimentation and statistical analysis. Also, an integration of a microscopic analysis with a signal-to-noise (GC(ave)/GC(std)) based stochastic XFEM modeling framework can be worthwhile.  129  References [1] Tuakta, Ch (2005) Use of fiber reinforced polymer composite in bridge structures. University of Massachusetts Institute of Technology. USA. [2] Meier, U (2000) Advanced solution with composites in construction. Proceeding of the international Conference on Composites in Construction, (Figueiras, J, Juvandes, L, Faria, R, & Eds.), pp. 3-7. [3] Mazumdar, S (2002) Composites manufacturing: materials, products and process engineering. CRC Press. [4] Potluri, P, & Manan, A (2007) Mechanics of non-orthogonally interlaced textile composites. Composites Part A: Applied Science and Manufacturing 38:1216-1226. [5] Keller, T (2001) Recent all-composites and hybrid fiber-reinforced polymer bridges and buildings. Progress in Structural Engineering Materials 3(2):132-140. [6] Durham, SD, & Pagett, WJ (1997) Cumulative damage models for system failure with application to carbon fibers and composites. American Statistical Association and the American Society for Quality Control 39(1):34-44. [7] Zhang, TF, Yang, XH, Sun, WSh, & Cai, ZJ (2011) A two-parameter model of FRP laminates stiffness reduction. Advanced Materials Research 236-238:1187-1194. [8] Prasad, MS, Venkatesha, CS, & Jayaraju, T (2011) Experimental methods of determining fracture toughness of fiber reinforced polymer composites under various loading conditions. Journal of Minerals & Materials Characterization & Engineering 10(13):12631275.  130  [9]  Sriramula,  S,  &  Chryssanthopoulos,  MK  (2009)  Quantification  of  uncertainty modelling in stochastic analysis of FRP composites. Composites Part A: Applied Science and Manufacturing 40(11):1673–1684. [10] Baxter, SC, & Graham, LL (2000) Characterization of random composites using moving-window technique. Journal of Engineering Mechanics 126(4):389-397. [11] Guilleminot, J, Soize, C, Kondo, D, & Binetruy, C (2008) Theoretical framework and experimental procedure for modelling mesoscopic volume fraction stochastic fluctuations in fiber reinforced composites. International Journal of Solids and Structures 45(21):5567– 5583. [12] Soize, C (2008) Tensor-valued random fields for meso-scale stochastic model of anisotropic elastic microstructure and probabilistic analysis of representative volume element size. Probabilistic Engineering Mechanics 23(2–3):307–323. [13] Whitcomb, J, Srirengan, K, & Chapman, CS (1994) Evaluation of homogenization for global/local stress analysis of textile composites. Collection of Technical Papers-AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference 3:1649-1663. [14] Mehrez, L, Moens, D, & Vandepitte, D (2012) Stochastic identification of composite material properties from limited experimental databases, part I: Experimental database construction. Mechanical Systems and Signal Processing 27:471–483. [15] Owen, MJ, Bishop, PT (1973) Critical stress intensity factors applied to glass reinforced polyester resins. Journal of Composite Materials 7:146-159.  131  [16] Gaggar, S, Broutman, LJ (1977) Fracture toughness of random glass fiber epoxy composites: an experimental investigation. Flaw Growth and Fracture ASTM STP 631:310330. [17] Mower, TM, Li, VC (1987) Fracture characterization of short fiber reinforced thermoset resin composites. Engineering Fracture Mechanics 26(4):593-603. [18] ASTM-International (2001) Standard Test Method for Mode I, D 5528 01: interlaminar fracture toughness of unidirectional fiber-reinforced polymer matrix composites. West Conshohocken: American Society for Testing and Materials. [19] O’Brien, TK (1998) Interlaminar fracture toughness: the long and winding road to standardization. Composites Part B: Engineering 29B(1): 57–62. [20] Davies, P (1993) ESIS protocols for interlaminar fracture testing of composites. France, IFREMER brochure. [21] Lee, SM (1993) An edge crack torsion method for mode III delamination fracture testing. Journal of Composite Technology and Research 15(3): 193-201. [22] Hillerborg, A, Modeer, M, & Petersson, PE (1976) Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cement and Concrete Research 6:773-782. [23] Xu, XP, & Needleman, A (1994) Numerical simulations of fast crack growth in brittle solids. Journal of the Mechanics and Physics of Solids 42:1397-1434. [24] Camacho, GT, & Ortiz, M, (1996) Computational modelling of impact damage in brittle materials. International Journal of Solids and Structures 33:2899-2938.  132  [25] Camanho, PP, Davila, CG, & De Moura, MF (2003) Numerical simulation of mixed-mode progressive delamination in composite materials. Journal of Composite Materials, 37:1415-24. [26] Blackman, BRK, Hadavinia, H, Kinloch, AJ, & Williams, JG (2003) The use of a cohesive zone model to study the fracture of fiber composites and adhesively-bonded joints. International Journal of Fracture 119:25-46. [27] Gao, YF, & Bower, AF (2004) A simple technique for avoiding convergence problems in finite element simulations of crack nucleation and growth on cohesive interfaces. Modelling and Simulation in Materials Science and Engineering 12:453-463. [28] Segurado, TMJ, & LLorca, CTJ (2004) A new three-dimensional interface finite element to simulate fracture in composites. International Journal of Solids and Structures 41:2977-2993. [29] Cox, B, & Yang, Q (2005) Cohesive models for damage evolution in laminated composites. International Journal of Fracture133:107-37. [30] Nishikawa, M, Okabe, T, & Takeda, N (2007) Numerical simulation of interlaminar damage propagation in CFRP cross-ply laminates under transverse loading. International Journal of Solids and Structures 44:3101-3113. [31] Belytschko, T, & Black, T (1999) Elastic crack growth in finite elements with minimal remeshing. International Journal for Numerical Methods in Engineering 45: 601620.  133  [32] Moёs, N, Dolbow, J, & Belytschko, T (1999) A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering 46: 131-150. [33] Melenk, JM, & Babuška, I (1996) The partition of unity finite element method: basic theory and applications. Computer Methods in Applied Mechanics and Engineering 139: 289-314. [34] Xiao, QZ, Karihaloo, BL, & Liu, XY (2006) Incremental-secant modulus iteration scheme and stress recovery for simulating cracking process in quasi-brittle materials using XFEM. International Journal for Numerical Methods in Engineering 69(12):26062635. [35] Unger, JF, Eckardta, S, & Könke, C (2007) Modelling of cohesive crack growth in concrete structures with the extended finite element method. Computer Methods in Applied Mechanics and Engineering 196(41):4087-4100. [36] Benvenuti, E (2008) A regularized XFEM framework for embedded cohesive interfaces. Computer Methods in Applied Mechanics and Engineering 197:4367-78. [37] Sosa, JLC, & Karapurath, N (2012) Delamination modelling of GLARE using the extended finite element method. Composites Science and Technology 72:788–791. [38] Kaw, AK (1997) Mechanics of composite materials. New York: CRC Press. [39] Peng, X, & Cao, J (2002) A dual homogenization and finite element approach for material characterization of textile composites. Composites Part B: Engineering 33: 45–56. [40] Baxevanakis, C, Jeulin, D, & Renard, J (1995) Fracture statistics of a unidirectional composite. International Journal of Fracture 73(2):149–181. 134  [41] Elleithy, R (2000) The hierarchical structure and flexure behavior of woven carbon fiber epoxy composite. Polymer Composites 27(5):716-723. [42] Sejnoha, M, & Zeman, J (2008) Micromechanical modeling of imperfect textile composites. International Journal of Engineering Science 46(6):513-526. [43] Trias, D, Costa, J, Mayugo, J, & Hurtado, J (2006) Random models versus periodic models for fiber reinforced composites. Computational Materials Science 38(2):316324. [44] Wang, ASD, Chou, PC, & Lei, SC (1984) A Stochastic Model for the Growth of Matrix Cracks in Composite Laminates. Journal of Composite Materials 18(3):239-254. [45] Yushanov, S, & Bogdanovich, A (1998) Stochastic theory of composite materials with random waviness of the reinforcements. International Journal of Solids and Structures 35(22):2901-2930. [46] Silberschmidt, VV (2005) Matrix cracking in cross-ply laminates: effect of randomness. Composites Part A: Applied Science and Manufacturing 36(2):129-135. [47] Skordos, A, & Sutcliffe, M (2008) Stochastic simulation of woven composites forming. Composites Science and Technology 68(1):283-296. [48] Ashcroft, IA, Khokar, ZR, & Silberschmidt VV (2012) Modelling the effect of micro-structural randomness on the mechanical response of composite laminates through the application of stochastic cohesive zone elements. Computational Materials Science 52(1):95– 100.  135  [49] Espinosa, HD, Dwivedi, S, & Lu, HC (2000) Modeling impact induced delamination of woven fiber reinforced composites with contact/cohesive laws. Computer Methods in Applied Mechanics and Engineering 183:259-290. [50] Turon, A, Da’vila, CG, Camanho, PP, & Costa, J (2007) An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models. Engineering Fracture Mechanics, 74:1665-1682. [51] Harper, P, & Hallett, SR (2008) Cohesive zone length in numerical simulations of composite delamination. Engineering Fracture Mechanics, 75(16):4774-4792. [52] Barbieri, E, & Meo, M (2009) A mesh-free penalty-based approach to delamination in composites. Composites Science and Technology 69(13):2169-2177. [53] Khoei, AR, Biabanaki, SOR, & Anahid, M (2009) A Lagrangian extended finite element method in modeling large plasticity deformations and contact problems. International Journal of Mechanical Sciences 51(5):384-401. [54] Bhattacharyya, D, Bowis, M, & Jayaraman, K (2003) Thermoforming woodfibre–polypropylene composite sheets. Composites Science and Technology 63(3):353365. [55] ASTM-International (2011) Test Method D 4762 11A: Standard Guide for Testing Polymer Matrix Composite Materials. West Conshohocken: American Society for Testing and Materials. [56] ASTM-International (2011) Test Method D 3039/D 3039M - 00: Standard Test Method for Tensile Properties of Polymer Matrix Composite. West Conshohocken: American Society for Testing and Materials. 136  [57] ASTM-International (2011) Test Method D 7617/7617M 11: Standard Test Method for Transverse Shear Strength of Fiber-reinforced Polymer Matrix Composite Bars. West Conshohocken: American Society for Testing and Materials. [58] Bisby, LA, & Fitzwilliam, J (2006) An introduction to FRP composites for construction. ISIS, A Canadian Network of Centres of Excellence. [59] Hellen, K (1984) Introduction to Fracture Mechanics. McGraw-Hill, New York. [60] Westergaard, HM (1939) Bearing pressures and cracks. Journal of Applied Mechanics. [61] Wells, AA (1961) Unstable crack propagation in metals: cleavage and fast fracture. Proceedings of the Crack Propagation Symposium, Paper 84. [62] Rice, JR (1968) Path-independent integral and the approximate analysis of strain concentration by notches and cracks. Journal of Applied Mechanics, Transactions ASME 35 (2):379–386. [63] Erdogan, F, & Sih, GC (1963) On the crack extension in plane loading and transverse shear. Journal of Basic Engineering 85:519-527. [64] Hoffman, O (1967) The brittle strength of orthotropic materials. Journal of Composite Materials 1(2):200-206. [65] Tsai, SW (1965) Strength characteristics of composite materials. NASA CR-224. [66] Tsai, SW, & Wu, EM (1971) A general theory of strength for anisotropic materials. Journal of Composite Materials 5(1):58-80.  137  [67] Hashin, Z (1980) Failure criteria for unidirectional fiber composites. Journal of Applied Mechanics 47(2):329-334. [68] Gdoutos, EE (2005) Fracture mechanics: An introduction. Springer, Netherlands. [69] Lekhnitskii, SG (1963) Theory of an anisotropic elastic body. San Francisco, Holden-Day. [70] Sih, GC, Paris, PC, & Irwin, GR (1965) On cracks in rectilinearly anisotropic bodies. International Journal of Fracture Mechanics. 1:189–203. [71] Asadpoure, A, & Mohammadi, S (2007) Developing new enrichment functions for crack simulation in orthotropic media by the extended finite element method. International Journal for Numerical Methods in Engineering. 69:2150-2172. [72] Dolbow, J, Moёs, N, & Belytschko, T (2001) An extended finite element method for modeling crack growth with frictional contact. Computational Methods in Applied Mechanics and Engineering, 190:6825-6846. [73] Kim, JH, & Paulino, GH (2003) The interaction integral for fracture of orthotropic functionally graded materials: evaluation of stress intensity factors. International Journal of Solids and Structures, 40(1):3967-4001. [74] Irwin, GR (1961) Plastic zone near a crack and fracture toughness. Proceedings, Seventh Sagamore Ordnance Materials Conference, Syracuse University Research Institute, 1:4-63. [75] Dugdale, DS (1960) Yielding in steel sheets containing slits. Journal of the Mechanics and Physics of Solids, 8(2):100-104.  138  [76] Barenblatt, GI (1962) The mathematical theory of equilibrium cracks in brittle fracture. Advances in Applied Mechanics, 7(1):55-129. [77] Eshelby, JD (1956) The continuum theory of lattice defects. Solid State Physics, 3(1):79-144. [78] Motamedi, D, & Mohammadi, S (2010) Dynamic crack propagation analysis of orthotropic media by the extended finite element method. International Journal of Fracture, 161(1):21-39. [79] Motamedi, D, & Mohammadi, S (2011) Fracture analysis of composites by time independent moving-crack orthotropic XFEM. International Journal of Mechanical Sciences, 54(1):20–37. [80] Aliabadia, MH, & Sollero, P (1998) Crack growth analysis in homogeneous orthotropic laminates. Composites Science and Technology, 58(10):1697–1703. [81] Wu, KC (2000) Dynamic crack growth in anisotropic material. International Journal of Fracture, 106:1-12. [82] Dongye, C, & Ting, TCT (1989) Explicit expressions of Barnett-Lothe tensors and their associated tensors for orthotropic materials. Quarterly of Applied Mathematics, 47: 723-734. [83] Needleman, A (1987) A continuum model for void nucleation by inclusion debonding. Journal of Applied Mechanics, 54(3):525-531. [84] Tvergaard, V, & Hutchinson, JW (1992) The relation between crack growth resistance and fracture process parameters in elastic–plastic solids. Journal of the Mechanics and Physics of Solids, 40(6):1377-1397. 139  [85] Cui, W, & Wisnom, MR (1993) A combined stress-based and fracture-mechanics based model for predicting delamination in composites. Composites, 24(6):467-474. [86] Geubelle, PH, & Baylor, JS (1998) The impact-induced delamination of laminated composites: A 2D simulation. Composites Part B: Engineering, 29(5):589-602. [87] Ortiz, M, & Pandolfi, A (1999) Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis. International Journal of Numerical Methods and Engineering, 44:1267-1282. [88] Mi, U, Crisfield, MA, & Davies, GAO (1998) Progressive delamination using interface elements. Journal of Composite Materials, 32(14):1246-1272. [89] Alfano, G, & Crisfield, MA (2001) Finite element interface models for the delamination analysis of laminated composites: mechanical and computational issues. International Journal for Numerical Methods in Engineering, 50:1701-1736. [90] Fan, C, Ben Jar, PY, & Cheng ,JJR (2008) Cohesive zone with continuum damage properties for simulation of delamination development in fiber composites and failure of adhesive joints. Engineering Fracture Mechanics, 75:3866-3880. [91] Planas, J, & Elices, M (1991) Nonlinear fracture of cohesive materials. International Journal of Fracture, 51:139-57. [92] Jin, ZH, & Sun, CT (2005) Cohesive zone modeling of interface fracture in elastic bi-materials. Engineering Fracture Mechanics, 72(12):1805-1817. [93] Yang, QD, Cox, BN, Nalla, RK, & Ritchie, RO (2006) Fracture length scales in human cortical bone: the necessity of nonlinear fracture models. Biomaterials, 27:2095-113.  140  [94] Carpinteri, A, Cornetti, P, Barpi, F, & Valente, S (2003) Cohesive crack model description of ductile to brittle size-scale transition: dimensional analysis vs. renormalization group theory. Engineering Fracture Mechanics, 70:1809-939. [95] Giner, E, Sukumar, N, Tarancon, JE, & Fuenmayor, FJ (2009) An ABAQUS implementation of the extended finite element method. Engineering Fracture Mechanics, 76(3):347–68. [96] SIMULIA (2008) ABAQUS 6.8-1 Documentation, Rhode Island, USA. [97] ASTM-International (2001) Standard Test Method for Mode I, D 5528 01: interlaminar fracture toughness of unidirectional fiber-reinforced polymer matrix composites. West Conshohocken: American Society for Testing and Materials. [98] Fukunaga, H, Tsu-Wei, C, Peters, PWM, & Schulte, K (1984) Probabilistic failure strength analysis of graphite epoxy cross ply laminates. Journal of Composite Materials, 18(4):339-56. [99] Bulsara, VN, Talreja, R, & Qu, J (1999) Damage initiation under transverse loading of unidirectional composites with arbitrarily distributed fibers. Composites Science and Technology, 59(5):673-82. [100] Silberschmidt, VV (2006) Effect of micro-randomness on macroscopic properties and fracture of laminates. Journal of Material Science, 41(20):6768-776. [101] Spearing, SM, & Evans, AG (1992) The role of fiber bridging in the delamination resistance of fiber-reinforced composites. Acta Metallurgica et Materialia, 40(9):2191-9.  141  [102] Tamuzs, V, Tarasovs, S, & Vilks, U (2001) Progressive delamination and fiber bridging modeling in double cantilever beam composite specimens. Engineering Fracture Mechanics, 68(5):513–25. [103] Nairn, J (2009) Analytical and numerical modeling of R curves for cracks with bridging zones. International Journal of Fracture, 155(2):167–81. [104] Airoldi, A, & Dávila, CG (2012) Identification of material parameters for modeling delamination in presence of fiber bridging. Composite Structures, 94:3240-9. [105] TENCATE ADVANCED COMPOSITES USA, INC. (www.tencate.com) [106] Khokhar, ZR, Ashcroft, IA, & Silberschmidt, VV (2009) Simulations of delamination in CFRP laminates: Effect of microstructural randomness. Computational Materials Science, 46(2):607-13.  142  Appendices Appendix A: ABAQUS User-element Subroutine for Nonlinear XFEM Analysis     SUBROUTINE UEL(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS, 1 PROPS,NPROPS,COORDS,MCRD,NNODE,U,DU,V,A,JTYPE,TIME,DTIME, 2 KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,NPREDF, 3 LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,NJPROP,PERIOD) C INCLUDE 'ABA_PARAM.INC' +  PARAMETER ( ZERO = 0.D0, HALF = 0.5D0, ONE = 1.D0, SEVEN=7.0D0, EIGHT=8.0D0 )  C DIMENSION RHS(NDOFEL,*),AMATRX(NDOFEL,NDOFEL),PROPS(*), 1 SVARS(*),ENERGY(8),COORDS(MCRD,NNODE),U(NDOFEL), 2 DU(NDOFEL,*),V(NDOFEL),A(NDOFEL),TIME(2),PARAMS(*), 3 JDLTYP(MDLOAD,*),ADLMAG(MDLOAD,*),DDLMAG(MDLOAD,*), 4 PREDEF(2,NPREDF,NNODE),LFLAGS(*),JPROPS(*) C C  C C C C  REAL*8 GAUSSPOINT(140,MCRD), + XCR(MCRD+1,MCRD), + dNdx(NNODE,MCRD+1), + U_MIDPNT(NNODE,MCRD), + B(MCRD*MCRD, NDOFEL), + DNDXi(NNODE,MCRD), + C_COORDS(NNODE,MCRD), + Gs(2*MCRD,MCRD*MCRD), + U_Ms(MCRD*MCRD,MCRD*MCRD), + C(2*MCRD,2*MCRD), + Dep(MCRD,MCRD), + QT(MCRD,MCRD), + QTT(MCRD,MCRD), + QR(MCRD,MCRD), + QRT(MCRD,MCRD), + DGep(MCRD,MCRD), + FB(MCRD,MCRD), + FBT(MCRD,MCRD), + STRGLB(MCRD,MCRD), + STRGLOB(MCRD,MCRD), + dUdx(MCRD,MCRD), + Nx(1,NNODE), + GWEIGHT(140), + FF(MCRD*MCRD), + STRSG(MCRD*MCRD), + COOR(MCRD), + DLTU(MCRD), + dV(MCRD) USER-ELEMNT ARRAYS GENERAL ELEMENT VALUES GAUSS INTEGRATION VARIABLES (3 INTEG POINT) ARRAYS FOR 3D ELEMENT REAL*8 U_PHI(NNODE) REAL*8 UNeg(MCRD), + UPos(MCRD),  143  + UNg(MCRD), + UPs(MCRD), + U1(MCRD), + UBAS(NDOFEL), + EC(2*MCRD), + STRS(2*MCRD), + BTSTRS(NDOFEL), + Trct(MCRD), + LHSC(NDOFEL), + RHSC(NDOFEL) REAL*8 BT(NDOFEL,MCRD*MCRD), + GsT(MCRD*MCRD,2*MCRD) REAL*8 JACB(MCRD,MCRD), + INVJACB(MCRD,MCRD), + INVFB(MCRD,MCRD), + INVFBT(MCRD,MCRD), + IUNIT(MCRD,MCRD), + DGBc(MCRD,NDOFEL), + BTGT(NDOFEL,2*MCRD), + GB(2*MCRD,NDOFEL), + MsB(MCRD*MCRD,NDOFEL), + BTMsB(NDOFEL,NDOFEL), + BcTDGBc(NDOFEL,NDOFEL), + KMAT(NDOFEL,NDOFEL), + KGEM(NDOFEL,NDOFEL), + KCON(NDOFEL,NDOFEL) REAL*8 Bc(MCRD,NDOFEL), + BcT(NDOFEL,MCRD), + CSTRN(MCRD,MCRD), + EEC(MCRD,MCRD), + BcRes(NDOFEL,MCRD) REAL*8 Xi, Yi, Zi, WEIGHT, ENRCOH, ENRJC, + E11, E22, E33, G12, G23, G31, + Nu12, Nu21, Nu23, Nu32, Nu31, Nu13, + DETFB, DETJ, nu, uacrt, Jsdv, + HPOINT, ENRELM, NNINT, KPEN,TMAX, + Jcrt, Dmg1, Dmg2, Dmg3 C INTEGER COUNT1, COUNT2, INTP, IINTP, I, J, K, L, M, N C C ************************************************************************************* C *** INITIALISATION: IMPORTANT!! FORTRAN DOES NOT PUT ZEROS IN THERE AUTOMATICALLY *** C ************************************************************************************* C OPEN LEVEL-SET VALUES FROM PRE-PROCESSING OPEN(16, FILE='C:/ABAQUS-Matlab/IOstat +/philvlst.dat',STATUS='OLD') C NINTP = 24 ! Number of integration points uacrt = PROPS(1) ! Crack Lenght dV(1) = PROPS(2) ! Material Initiation Opening Failure dV(2) = PROPS(3) ! Material Final Opening Failure E11 = PROPS(4) ! Material Constant E22 = PROPS(5) ! Material Constant E33 = PROPS(6) ! Material Constant G12 = PROPS(7) ! Material Constant G23 = PROPS(8) ! Material Constant G31 = PROPS(9) ! Material Constant Nu12 = PROPS(10) ! Material Constant  144  Nu23 = PROPS(11) ! Material Constant Nu31 = PROPS(12) ! Material Constant KPEN = PROPS(13) ! Material Constant Jsdv = PROPS(14) ! G Standard Deviation C C  CALCULATING THE TRACTION-SEPARATION CONSTANT TMAX = KPEN*dV(1) Jcrt = TMAX*dV(2)/2 WIDTH = 1 ! Element Tickness C ********************************** C *** ZERO THE REQUIRED MATRICES *** C ********************************** ENRELM = 0.0 C DO I = 1, NNODE Nx(1,I) = 0.0 U_PHI(I) = 0.0 DO J = 1, MCRD+1 dNdx(I,J) = 0.0 ENDDO ENDDO C DO J = 1, MCRD COOR(J) = 0.0 ENDDO C DO I = 1, 2*MCRD DO J = 1, 2*MCRD C(I,J) = 0.0 ENDDO ENDDO C DO I = 1, MCRD*MCRD DO J = 1, NDOFEL B(I,J) = 0.0 BT(J,I) = 0.0 ENDDO ENDDO C DO I = 1, 2*MCRD DO J = 1, MCRD*MCRD Gs(I,J) = 0.0 GsT(J,I) = 0.0 ENDDO ENDDO C DO I = 1, NNODE DO J = 1, MCRD C_COORDS(I,J) = 0.0 ENDDO ENDDO C DO I = 1, MCRD DO J = 1, MCRD IUNIT(I,J) = 0.0 JACB(I,J) = 0.0 INVJACB(I,J) = 0.0 Dep(I,J) = 0.0 DGep(I,J) = 0.0  145  QT(I,J) = 0.0 QTT(I,J) = 0.0 QR(I,J) = 0.0 QRT(I,J) = 0.0 ENDDO ENDDO C DO I = 1, NNODE DO J = 1, MCRD U_MIDPNT(I,J) = 0.0 ENDDO ENDDO C DO I = 1, NSVARS SVARS(I) = 0.0 ENDDO C DO I = 1, NDOFEL DO J = 1, NDOFEL AMATRX(I,J) = 0.0 KGEM(I,J) = 0.0 KMAT(I,J) = 0.0 KCON(I,J) = 0.0 ENDDO ENDDO C  C C C  C C C  DO I = 1, NDOFEL RHS(I,1) = 0.0 RHSC(I) = 0.0 LHSC(I) = 0.0 UBAS(I) = 0.0 ENDDO ******************** *** DUMMY ARRAYS *** ******************** DO I = 1, MCRD IUNIT(I,I) = 1.0 ENDDO INTS = 1 ! Integration point scheme (1: gauss) STYPE = 1 ! Element type (1: B8, 2: T4) **************** *** HOOK LAW *** **************** Nu21 = Nu12*E22/E11 Nu32 = Nu23*E33/E22 Nu13 = Nu31*E11/E33 Delt = (1-Nu12*Nu21-Nu23*Nu32-Nu13*Nu31-2*Nu12*Nu23*Nu31)  C C(1,1) = E11*(1-Nu32*Nu23)/Delt C(1,2) = E11*(Nu21+Nu31*Nu23)/Delt C(1,3) = E11*(Nu31+Nu21*Nu32)/Delt C(2,1) = E22*(Nu12+Nu13*Nu23)/Delt C(2,2) = E22*(1-Nu13*Nu31)/Delt C(2,3) = E22*(Nu32+Nu31*Nu12)/Delt C(3,1) = E33*(Nu13+Nu12*Nu23)/Delt C(3,2) = E33*(Nu23+Nu13*Nu21)/Delt C(3,3) = E33*(1-Nu12*Nu21)/Delt C(4,4) = G12/2 C(5,5) = G23/2  146  C C C  C C C  C C C C C C  C(6,6) = G31/2 ******************************************* *** FINDING DEFORM SHAPE OF COORDINATES *** ******************************************* DO I = 1, NNODE DO J = 1, MCRD C_COORDS(I,J) = COORDS(J,I) + U(2*MCRD*(I-1)+J) ENDDO ENDDO ************************************************ *** FINDING REFERENCE COORDINATE DEFORMATION *** ************************************************ DO I = 1, NDOFEL U(I) = U(I) - DU(I,1) ENDDO ************************* *** CALLING LEVEL SET *** ************************* CALL LVLSETRDR(U_PHI,JELEM,NNODE) ! READ THE Level-Set FROM MATLAB ******************************************************* *** CALLING THE LOCAL CRACK'S PLANE IN EACH ELEMENT *** ******************************************************* CALL MIDPLNFIND(NNINT,U_MIDPNT,QT,U_PHI,COORDS,NNODE,MCRD)  C  C C C C C C  DO I = 1, MCRD DO J = 1, MCRD QTT(I,J) = QT(J,I) ENDDO ENDDO ********************************* *** CALLING GAUSS COORDINATES *** ********************************* CALL SUBTGAUSS(GAUSSPOINT,GWEIGHT,NNODE,MCRD) **************************************** *** LOOKING FOR OPENING DISPLACEMENT *** **************************************** COUNT1 = 0 COUNT2 = 0  C DO I = 1, MCRD UPos(I) = 0.0 UNeg(I) = 0.0 UPs(I) = 0.0 UNg(I) = 0.0 ENDDO C ************************************ C *** LOOP OVER INTEGRATION POINTS *** C ************************************ DO INTP = 1, NINTP HPOINT = 0.0 DO I = 1, MCRD U1(I) = 0.0 DLTU(I) = 0.0 ENDDO C COOR(1) = GAUSSPOINT(INTP,1) COOR(2) = GAUSSPOINT(INTP,2) COOR(3) = GAUSSPOINT(INTP,3) WEIGHT = GWEIGHT(INTP)  147  C C C C C C  C C C  *********************************** *** CALLING THE SHAPE FUNCTIONS *** *********************************** CALL LAGRANGEBASIS(COOR,Nx,dNdx,NNODE,MCRD) ********************************* *** CALCULATING Gpt-LEVEL SET *** ********************************* DO I = 1, NNODE HPOINT = HPOINT + U_PHI(I)*dNdx(I,4)/ABS(U_PHI(I)) ENDDO IF (HPOINT .GT. 0.0) THEN HPOINT = 1.0 ELSEIF (HPOINT .LT. 0.0) THEN HPOINT = -1.0 ELSE HPOINT = 0.01 WEIGHT = 0.0 ENDIF ******************************************** *** CHECKING FOR Gpt WITHIN CONTACTBOUND *** ******************************************** DO I = 1, NNODE DO J = 1, MCRD U1(J) = U1(J) + dNdx(I,4)*U(6*(I-1)+J) + + (HPOINT - ABS(U_PHI(I))/U_PHI(I))* + (dNdx(I,4)*U(6*(I)+J-3)) ENDDO ENDDO  C IF (HPOINT .LT. 0.0)THEN COUNT1 = COUNT1 + 1 UNeg(1) = UNeg(1) + U1(1) UNeg(2) = UNeg(2) + U1(2) UNeg(3) = UNeg(3) + U1(3) ELSE COUNT2 = COUNT2 + 1 UPos(1) = UPos(1) + U1(1) UPos(2) = UPos(2) + U1(2) UPos(3) = UPos(3) + U1(3) ENDIF C ENDDO DO I = 1, MCRD DO J = 1, MCRD UPs(I) = UPs(I) + QT(I,J)*UPos(J)/COUNT2 UNg(I) = UNg(I) + QT(I,J)*UNeg(J)/COUNT1 ENDDO ENDDO C DLTU(1) = UPs(1) - UNg(1) DLTU(2) = UPs(2) - UNg(2) DLTU(3) = UPs(3) - UNg(3) c C C C C C C  ********************************************* *** FORMING STIFFNESS AND RESIDUAL MATRIX *** ********************************************* ******************************* *** CALLING ELAST-PLAST RELATIONSHIP ***  148  C  ******************************* DO I = 1,MCRD DO J = 1,MCRD DGep(I,J) = 0.0 ENDDO ENDDO  C CALL ELSPLC(Dep,DLTU,dV,uacrt,C,MCRD,ENRELM,JTYPE,KPEN, 1 Dmg1, Dmg2, Dmg3) C C C  ASSIGNING DAMAGE INDECES AND CRACK OPENING DISPLACEMENT TO USER-VARIABLES SVARS(1)=Dmg1 SVARS(2)=Dmg2 SVARS(3)=Dmg3 SVARS(4)=DLTU(1) SVARS(5)=DLTU(2) SVARS(6)=DLTU(3)  C C C  TRANSFORIMG THE ELASTIC-PLASTIC RELATIONSHIP INTO LOCAL CRACK PLANE DO I = 1, MCRD DO J = 1, MCRD DO K = 1, MCRD DO L = 1, MCRD DGep(I,J) = DGep(I,J) + QTT(I,K)*Dep(K,L)*QT(L,J) ENDDO ENDDO ENDDO ENDDO  C C C C  ************************************ *** LOOP OVER INTEGRATION POINTS *** ************************************ DO I = 1, NDOFEL DO J = 1, NDOFEL BTMsB(I,J) = 0.0 ENDDO ENDDO  C  C C C  C C  DO IINTP = 1, NINTP COOR(1) = GAUSSPOINT(IINTP,1) COOR(2) = GAUSSPOINT(IINTP,2) COOR(3) = GAUSSPOINT(IINTP,3) WEIGHT = GWEIGHT(IINTP) *********************************** *** CALLING THE SHAPE FUNCTIONS *** *********************************** CALL LAGRANGEBASIS(COOR,Nx,dNdx,NNODE,MCRD) Xi = COOR(1) Yi = COOR(2) Zi = COOR(3) DO I = 1, MCRD DO J = 1, MCRD JACB(I,J) = 0.0 ENDDO ENDDO ******************************* *** FORMING JACOBIAN MATRIX ***  149  C  ******************************* DO I = 1, MCRD DO J = 1, MCRD DO K = 1, NNODE JACB(I,J) = JACB(I,J) + COORDS(I,K)*dNdx(K,J) ENDDO ENDDO ENDDO  C + +  C C C  DETJ = JACB(1,1)*(JACB(2,2)*JACB(3,3) - JACB(2,3)*JACB(3,2)) JACB(1,2)*(JACB(2,1)*JACB(3,3) - JACB(3,1)*JACB(2,3)) + JACB(1,3)*(JACB(2,1)*JACB(3,2) - JACB(2,2)*JACB(3,1)) IF (DETJ .LT. 0.0) THEN DETJ = (-1)*DETJ ENDIF INVJACB(1,1) = (JACB(2,2)*JACB(3,3) - JACB(2,3)*JACB(3,2))/DETJ INVJACB(1,2) = (JACB(1,3)*JACB(3,2) - JACB(1,2)*JACB(3,3))/DETJ INVJACB(1,3) = (JACB(1,2)*JACB(2,3) - JACB(1,3)*JACB(2,2))/DETJ INVJACB(2,1) = (JACB(2,3)*JACB(3,1) - JACB(2,1)*JACB(3,3))/DETJ INVJACB(2,2) = (JACB(1,1)*JACB(3,3) - JACB(1,3)*JACB(3,1))/DETJ INVJACB(2,3) = (JACB(1,3)*JACB(2,1) - JACB(1,1)*JACB(2,3))/DETJ INVJACB(3,1) = (JACB(2,1)*JACB(3,2) - JACB(2,3)*JACB(3,1))/DETJ INVJACB(3,2) = (JACB(1,2)*JACB(3,1) - JACB(1,1)*JACB(3,2))/DETJ INVJACB(3,3) = (JACB(2,2)*JACB(1,1) - JACB(2,1)*JACB(1,2))/DETJ ******************************** *** FORMING DERIVATIES dN/dx *** ******************************** DO I = 1, NNODE DO J = 1, MCRD DNDXi(I,J) = 0.0 ENDDO ENDDO  C  C C C  DO I = 1, NNODE DO J = 1, MCRD DO K = 1, MCRD DNDXi(I,J) = DNDXi(I,J) + dNdx(I,K)*INVJACB(K,J) ENDDO ENDDO ENDDO ******************************* *** CALCULATING B-LEVEL SET *** ******************************* HPOINT = 0.0 DO I=1,MCRD DO K=1,MCRD dUdx(I,K)=0.0 ENDDO ENDDO  C  C C  DO I=1,MCRD DO K=1,MCRD DO J=1,NNODE dUdx(I,K)=dUdx(I,K)+DNDXi(J,K)*U(6*(J-1)+I) ENDDO ENDDO ENDDO ********************************* *** CALCULATING Gpt-LEVEL SET ***  150  C  ********************************* DO I = 1, NNODE HPOINT = HPOINT + U_PHI(I)*dNdx(I,4)/ABS(U_PHI(I)) ENDDO IF (HPOINT .GT. 0.0) THEN HPOINT = 1.0 ELSEIF (HPOINT .LT. 0.0) THEN HPOINT = -1.0 ELSE HPOINT = 0.01 WEIGHT = 0.0 ENDIF DO I=1,MCRD DO K=1,MCRD DO J=1,NNODE dUdx(I,K) = dUdx(I,K)+(HPOINT-ABS(U_PHI(I))/U_PHI(I))* + DNDXi(J,K)*U(6*(J-1)+I+MCRD) ENDDO ENDDO ENDDO C ************************ C *** FORMING B MATRIX *** C ************************ DO I = 1, MCRD*MCRD DO J = 1, NDOFEL B(I,J) = 0.0 BT(J,I) = 0.0 ENDDO ENDDO C CALL GRADB(B,DNDXi,HPOINT,U_PHI,NNODE,MCRD,NDOFEL) C BT = TRANSPOSE(B) C DO I = 1, MCRD*MCRD FF(I) = 0.0 ENDDO DO I = 1, 2*MCRD EC(I)=0.0 ENDDO C DO I = 1, MCRD DO J = 1, MCRD FB(I,J) = 0.0 FBT(I,J) = 0.0 ENDDO ENDDO C FF = MATMUL(B,U) C dUdx(1,1) = FF(1) dUdx(2,1) = FF(2) dUdx(3,1) = FF(3) dUdx(1,2) = FF(4) dUdx(2,2) = FF(5) dUdx(3,2) = FF(6) dUdx(1,3) = FF(7) dUdx(2,3) = FF(8) dUdx(3,3) = FF(9)  151  C C C C C  ************************************************** *** FORMING GREEN-LAGRANGIAN & LARGE STRAINS *** ************************************************** DO I = 1, MCRD*MCRD IF (I .EQ. 1 .OR. I .EQ. 5 .OR. I .EQ. 9) THEN FF(I) = FF(I) + 1.0 ENDIF ENDDO  C FB(1,1) = FF(1) FB(2,1) = FF(2) FB(3,1) = FF(3) FB(1,2) = FF(4) FB(2,2) = FF(5) FB(3,2) = FF(6) FB(1,3) = FF(7) FB(2,3) = FF(8) FB(3,3) = FF(9) C ******************************* C *** FORMING INVF & GRADF *** C ******************************* DO I = 1, MCRD DO J = 1, MCRD INVFB(I,J) = 0.0 ENDDO ENDDO DETFB = FB(1,1)*(FB(2,2)*FB(3,3) - FB(2,3)*FB(3,2)) + FB(1,2)*(FB(2,1)*FB(3,3) - FB(3,1)*FB(2,3)) + + FB(1,3)*(FB(2,1)*FB(3,2) - FB(2,2)*FB(3,1)) IF (DETFB .LT. 0.0) THEN DETFB = (-1)*DETFB ENDIF INVFB(1,1) = (FB(2,2)*FB(3,3) - FB(2,3)*FB(3,2))/DETFB INVFB(1,2) = (FB(1,3)*FB(3,2) - FB(1,2)*FB(3,3))/DETFB INVFB(1,3) = (FB(1,2)*FB(2,3) - FB(1,3)*FB(2,2))/DETFB INVFB(2,1) = (FB(2,3)*FB(3,1) - FB(2,1)*FB(3,3))/DETFB INVFB(2,2) = (FB(1,1)*FB(3,3) - FB(1,3)*FB(3,1))/DETFB INVFB(2,3) = (FB(1,3)*FB(2,1) - FB(1,1)*FB(2,3))/DETFB INVFB(3,1) = (FB(2,1)*FB(3,2) - FB(2,3)*FB(3,1))/DETFB INVFB(3,2) = (FB(1,2)*FB(3,1) - FB(1,1)*FB(3,2))/DETFB INVFB(3,3) = (FB(2,2)*FB(1,1) - FB(2,1)*FB(1,2))/DETFB C FBT = TRANSPOSE(FB) INVFBT = TRANSPOSE(INVFB) C DO I = 1, 2*MCRD DO J = 1, MCRD*MCRD Gs(I,J) = 0.0 GsT(J,I) = 0.0 ENDDO ENDDO C CALL TANGSMS(Gs,FF,MCRD) C GsT = TRANSPOSE(Gs) GB = MATMUL(Gs,B) EC = MATMUL(GB,U)  152  C C C  ************************** *** STRESS CALCULATION *** ************************** DO I = 1, 2*MCRD STRS(I) = 0.0 ENDDO DO I = 1, MCRD*MCRD STRSG(I) = 0.0 ENDDO DO I = 1, MCRD DO J = 1, MCRD STRGLB(I,J) = 0.0 STRGLOB(I,J) = 0.0 ENDDO ENDDO  C STRS = MATMUL(C,EC) C STRGLOB(1,1) = STRS(1) STRGLOB(2,2) = STRS(2) STRGLOB(3,3) = STRS(3) STRGLOB(2,1) = STRS(4) STRGLOB(1,2) = STRS(4) STRGLOB(1,3) = STRS(5) STRGLOB(3,1) = STRS(5) STRGLOB(2,3) = STRS(6) STRGLOB(3,2) = STRS(6) C  C C C  STRGLB = MATMUL(STRGLOB,IUNIT) STRSG(1) = STRGLB(1,1) STRSG(2) = STRGLB(1,2) STRSG(3) = STRGLB(1,3) STRSG(4) = STRGLB(2,1) STRSG(5) = STRGLB(2,2) STRSG(6) = STRGLB(2,3) STRSG(7) = STRGLB(3,1) STRSG(8) = STRGLB(3,2) STRSG(9) = STRGLB(3,3) ************************************ **** FORMING Gs & Ms MATRICES **** ************************************ DO I = 1, MCRD**2 DO J = 1, MCRD**2 U_Ms(I,J) = 0.0 ENDDO ENDDO  C DO I = 1, MCRD U_Ms(I,I) = STRSG(1) U_Ms(I+MCRD,I+MCRD) = STRSG(5) U_Ms(I+2*MCRD,I+2*MCRD) = STRSG(9) U_Ms(I+MCRD,I) = STRSG(4) U_Ms(I,I+MCRD) = STRSG(2) U_Ms(I+2*MCRD,I) = STRSG(7) U_Ms(I,I+2*MCRD) = STRSG(3) U_Ms(I+2*MCRD,I+MCRD) = STRSG(8) U_Ms(I+MCRD,I+2*MCRD) = STRSG(6) ENDDO C ********************************************  153  C C  *** FORMING COHESIVE & CONTACT Bc MATRIX *** ******************************************** DO I = 1, MCRD DO J = 1, NDOFEL Bc(I,J) = 0.0 BcT(J,I) = 0.0 BcRes(J,I) = 0.0 ENDDO ENDDO DO I = 1, NNODE BcRes(6*I-5,1) = 0.0 BcRes(6*I-4,2) = 0.0 BcRes(6*I-3,3) = 0.0 BcRes(6*I-2,1) = dNdx(I,4)*(HPOINT-ABS(U_PHI(I))/U_PHI(I)) BcRes(6*I-1,2) = dNdx(I,4)*(HPOINT-ABS(U_PHI(I))/U_PHI(I)) BcRes(6*I,3) = dNdx(I,4)*(HPOINT-ABS(U_PHI(I))/U_PHI(I)) Bc(1,6*I-2) = -2*dNdx(I,4)*(HPOINT-ABS(U_PHI(I))/U_PHI(I)) Bc(2,6*I-1) = -2*dNdx(I,4)*(HPOINT-ABS(U_PHI(I))/U_PHI(I)) Bc(3,6*I) = -2*dNdx(I,4)*(HPOINT-ABS(U_PHI(I))/U_PHI(I)) ENDDO  C C C C  BcT = TRANSPOSE(Bc) ******************************** *** FORMING STIFFNESS MATRIX *** ******************************** DO I = 1, 2*MCRD DO J = 1, NDOFEL GB(I,J) = 0.0 BTGT(J,I) = 0.0 ENDDO ENDDO  C GB = MATMUL(Gs,B) BTGT = TRANSPOSE(GB) C DO I=1, NDOFEL DO J=1, NDOFEL DO K=1, 2*MCRD DO L=1, 2*MCRD KMAT(I,J) = KMAT(I,J) + + BTGT(I,K)*C(K,L)*GB(L,J)*WEIGHT*DETJ*WIDTH ENDDO ENDDO ENDDO ENDDO C MsB = MATMUL(U_Ms,B) BTMsB = MATMUL(BT,MsB) C DO I = 1, NDOFEL DO J = 1, NDOFEL KGEM(I,J) = KGEM(I,J) + BTMsB(I,J)*WEIGHT*DETJ*WIDTH ENDDO ENDDO C DGBc = MATMUL(DGep,Bc) BcTDGBc = MATMUL(BcT,DGBc) C DO I = 1, NDOFEL  154  DO J = 1, NDOFEL KCON(I,J) = KCON(I,J) + BcTDGBc(I,J)*WEIGHT*DETJ*WIDTH ENDDO ENDDO C DO I = 1, MCRD Trct(I) = 0.0 ENDDO C  C C C  DO I = 1, MCRD DO J = 1, NDOFEL Trct(I) = TrcT(I) + DGBc(I,J)*(U(J)+DU(J,1)) ENDDO ENDDO **************************** *** RESIDUAL CALCULATION *** **************************** DO I = 1, MCRD*MCRD FF(I) = 0.0 ENDDO DO I = 1, 2*MCRD EC(I)=0.0 ENDDO  C  C C C  DO I = 1, MCRD*MCRD DO J = 1, NDOFEL FF(I) = FF(I) + B(I,J)*(U(J)+DU(J,1)) ENDDO ENDDO ************************************************* *** FORMING GREEN-LAGRANGIAN & LARGE STRAIN *** ************************************************* DO I = 1, MCRD*MCRD IF (I .EQ. 1 .OR. I .EQ. 5 .OR. I .EQ. 9) THEN FF(I) = FF(I) + 1.0 ENDIF ENDDO  C FB(1,1) = FF(1) FB(2,1) = FF(2) FB(3,1) = FF(3) FB(1,2) = FF(4) FB(2,2) = FF(5) FB(3,2) = FF(6) FB(1,3) = FF(7) FB(2,3) = FF(8) FB(3,3) = FF(9) C DO I = 1, 2*MCRD DO J = 1, MCRD*MCRD Gs(I,J) = 0.0 ENDDO ENDDO C CALL TANGSMS(Gs,FF,MCRD) C GB = MATMUL(Gs,B) BTGT = TRANSPOSE(GB) C  155  DO I=1, MCRD DO J=1, MCRD CSTRN(I,J) = 0.0 EEC(I,J) = 0.0 ENDDO ENDDO C DO I=1, MCRD DO J=1, MCRD DO K=1, MCRD CSTRN(I,J) = FB(K,I)*FB(K,J)+CSTRN(I,J) ENDDO ENDDO ENDDO DO I=1, MCRD DO J=1, MCRD IF (I .EQ. J) THEN EEC(I,J) = 0.5*CSTRN(I,J)-0.5 ELSE EEC(I,J) = 0.5*CSTRN(I,J) ENDIF ENDDO ENDDO C EC(1) = EEC(1,1) EC(2) = EEC(2,2) EC(3) = EEC(3,3) EC(4) = EEC(1,2) EC(5) = EEC(2,3) EC(6) = EEC(3,1) C STRS = MATMUL(C,EC) C C C  XFEM ELEMENTS RESIDUAL FORCES DUE TO LARGE DEFORMATION DO I=1, NDOFEL DO J=1, 2*MCRD RHS(I,1) = RHS(I,1) + BTGT(I,J)*STRS(J)*WEIGHT*DETJ*WIDTH ENDDO ENDDO  C C C  XFEM RESIDUAL FORCES DUE TO COHESIVE REGION OR CONTACT INTERFACE  DO I=1, NDOFEL DO J=1, MCRD RHS(I,1) = RHS(I,1) + + BcRes(I,J)*Trct(J)*WEIGHT*DETJ*WIDTH RHSC(I) = RHSC(I) + + BcRes(I,J)*Trct(J)*WEIGHT*DETJ*WIDTH ENDDO ENDDO ENDDO ! END IF LOOP OVER GAUSS POINTS C ************************************** C *** FORMING RIGHT-HAND-SIDE MATRIX *** C ************************************** AMATRX = KMAT + KGEM + KCON RETURN END  156  Appendix B: Experimental Calculations According to ASTM D5528-01 [97]  Double Cantilever Beam (DCB) Test: I.  The Modified Beam Theory (MBT) Method:  Based on the beam theory expressions, the critical energy release rate for the DCB test is written as follows [97] (see also Figure 1-4):  GIC   3Pgrip  (B-1)  2wacr  where Pgrip,  , w and acr are the applied load, the displacement of the load, specimen width and the delamination length, respectively. Due to possible rotation in the crack tip front, which is neglected in standard beam theory, a correction factor is considered by treating the test with a longer crack length, acr+ acr . To determine the crack length increase, acr , experimentally, the cube root ratio of the opening displacement over the applied load, known as beam compliance, must be plotted versus delamination length, acr, at onset of all delamination propagations. The least square root method should then be utilized to find the yintercept which is equal to the crack length increase. By substituting the new crack length into the above beam theory equation, one obtains:  GIC   3Pgrip   (B-2)  2w(acr  acr )  157  II.  The Compliance Calibration (CC) Method:  In the compliance calibration (CC) method, the least squares of log(   ) versus log(acr ) at Pgrip  onset of all delamination propagations must be plotted and the slope of this line, n, is implemented to correct the energy release rate:  GIC  III.  nPgrip   (B-3)  2wacr  The Modified Compliance Calibration (MCC) Method:  In the modified compliance calibration (CC) method, the normalized delamination length over the thickness versus the cube root ratio of opening displacement over the applied load at onset of all delamination propagations must be plotted and the slope of this line, mLine , is implemented to calculate the energy release rate: 2  3Pgrip ( GIC    2/3 ) Pgrip  (B-4)  2mLine w t ck  where tck is the specimen thickness. Finally, in order to account for shortening of arm of moment and the rotation of the loading block, a large displacement correction factor should be multiplied into the energy release value calculated from either of the above methods. This correction factor is calculated as follows:  158  Al arg e  1   3   10  acr      2    t ck t Block    3  4 2  2 acr           (B-5)  where t Block is the loading block thickness. End Notch Flexure (ENF) Experimental Test: Similar to the DCB test, for ENF test data the compliance method is employed to calculate the critical energy release rate [19]: 2  GIIC   9acr Pgrip  2  2 E Longitudinal w 2 t ck  3  (B-6)  159  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

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

Comment

Related Items