Applied Science, Faculty of
Materials Engineering, Department of
DSpace
UBCV
Li, Zhuo
Master of Applied Science - MASc
University of British Columbia
Brittle fracture has been investigated for more than a century. Various models managed to capture the fracture phenomena, among which the phase-field models gained popularity in the mechanics and the physics community in the past two decades. However, due to distinct backgrounds, detailed comparisons between the phase-field models in the two communities are missing. Therefore, to understand the strengths and limitations of these models, this study aims to compare them in the same settings. Quasi-static and dynamic brittle fractures are investigated in the current work, of which the main difference is the loading speed. Since numerous phase-field models have been developed in both communities, in the present work, only the most representative models are selected, i.e., Miehe et al.'s model and the Karma-Kessler-Levine (KKL) model. Through comparisons with the Griffith theory, in quasi-static simulations, the unphysical dependence of the energy release rate on the load in the KKL model was observed. Meanwhile, in dynamic fracture, both models showed good agreement with the theories by Eshelby and Adda-Bedia in some aspects while deviations in others. The results suggest that for quasi-static fracture, Miehe et al.’s model is a reasonable option. While for dynamic fracture, comparisons with experiments are necessary to evaluate the accuracy of the models.
Phase-Field Modelling of CrackPropagationbyZhuo LiA THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Materials Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)June 2021© Zhuo Li 2021The following individuals certify that they have read, and recommend to the Faculty ofGraduate and Postdoctoral Studies for acceptance, the thesis entitled:Phase-Field Modelling of Crack Propagationsubmitted by Zhuo Li in partial fulfillment of the requirements for the degree of Masterof Applied Science in Materials Engineering.Examining Committee:Matthias Militzer, Professor, Materials Engineering Department, UBCSupervisorWarren J. Poole, Professor, Materials Engineering Department, UBCSupervisory Committee MemberChadwick Sinclair, Professor, Materials Engineering Department, UBCSupervisory Committee MemberiiAbstractBrittle fracture has been investigated for more than a century. Various models managed tocapture the fracture phenomena, among which the phase-field models gained popularity inthe mechanics and the physics community in the past two decades. However, due to distinctbackgrounds, detailed comparisons between the phase-field models in the two communitiesare missing. Therefore, to understand the strengths and limitations of these models, thisstudy aims to compare them in the same settings. Quasi-static and dynamic brittle fracturesare investigated in the current work, of which the main difference is the loading speed.Since numerous phase-field models have been developed in both communities, in thepresent work, only the most representative models are selected, i.e., Miehe et al.’s model[1] and the Karma-Kessler-Levine (KKL) model [2]. Through comparisons with the Griffiththeory [3], in quasi-static simulations, the unphysical dependence of the energy release rateon the load in the KKL model was observed. Meanwhile, in dynamic fracture, both modelsshowed good agreement with the theories by Eshelby [4] and Adda-Bedia [5] in some aspectswhile deviations in others.The results suggest that for quasi-static fracture, Miehe et al.’s model is a reasonableoption. While for dynamic fracture, comparisons with experiments are necessary to evaluatethe accuracy of the models.iiiLay SummaryA seemingly unimportant crack can lead to disasters like airplane crashes, pipeline explosions,etc. Millions of dollars are lost each year because of structural failures. Fracture mechanics,as a subject dealing with crack evolution is, therefore, extremely important in engineeringdesigns. However, after several decades of development, there are still challenges in under-standing fracture mechanisms. As a reliable and inexpensive alternative for experiments,computational modelling quickly becomes an indispensable part of fracture mechanics. Inrecent years, phase-field models stood out among fracture models and gained popularity inthe mechanics and the physics community. In detail, both communities use different formu-lations of these models to capture fracture phenomena, but the comparison between them ismissing in the literature. Therefore, the present work aims to compare these models so thatresearchers who are new to the field can decide which model to use.ivPrefaceThe purpose of the Master’s thesis is to have detailed comparisons between two popularphase-field fracture models by the mechanics and the physics community. The thesis waswritten to fulfill the graduation requirements of the Master of Applied Science Programat the Department of Materials Engineering, The University of British Columbia, PointGrey campus. I was responsible for developing the computer programs, implementing thesimulations, analyzing the data, and writing the thesis.The research question was brought out by my supervisor, Prof. Matthias Militzer, whoalso provided assistance and guidance throughout every stage of the research. This thesis iswritten by the author with discussions and editing assistance from Prof. Matthias Militzer.In Chapter 2, Figures 2.7 and 2.9 have been taken with permission from the cited re-sources.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xivList of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32.2 Crack Propagation in Quasi-static Brittle Fracture . . . . . . . . . . . . . . 32.3 Basic Phenomena in Dynamic Brittle Fracture . . . . . . . . . . . . . . . . 72.4 Brief Overview of the Fracture Models . . . . . . . . . . . . . . . . . . . . . 82.5 The Phase-Field Fracture Models in the Mechanics Community . . . . . . . 112.5.1 Basics of the Mechanical Models . . . . . . . . . . . . . . . . . . . . 112.5.2 Interpretations of the Internal Length Scale (l) . . . . . . . . . . . . 152.5.3 Model Variations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.5.4 Model Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.6 The Phase-Field Fracture Models in the Physics Community . . . . . . . . . 25vi2.6.1 General Overview of the Phase-Field Fracture Models in the PhysicsCommunity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 252.6.2 Basics of the KKL Model . . . . . . . . . . . . . . . . . . . . . . . . 262.6.3 Applications of the KKL Model . . . . . . . . . . . . . . . . . . . . . 283 Scope and Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 314 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 324.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 324.2 Energy Functional and Governing Equations . . . . . . . . . . . . . . . . . 334.2.1 The Mechanical Model . . . . . . . . . . . . . . . . . . . . . . . . . . 334.2.2 The KKL Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 344.3 Numerical Techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 354.4 Initial and Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . 364.4.1 Initial Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 364.4.2 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . 375 Benchmarking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 405.1 Sensitivity Analysis of the Threshold Values in the Iteration Methods . . . . 405.2 Sensitivity Analysis of the Mesh Size . . . . . . . . . . . . . . . . . . . . . . 425.3 Sensitivity Analysis of the Time Step Size . . . . . . . . . . . . . . . . . . . 455.4 Sensitivity Analysis of the Crack Front Position . . . . . . . . . . . . . . . . 475.4.1 φfront in the Quasi-static Setting . . . . . . . . . . . . . . . . . . . . 475.4.2 φfront in the Dynamic Setting . . . . . . . . . . . . . . . . . . . . . . 505.5 Sensitivity Analysis of the Geometric Configuration in the Quasi-static Setting 545.6 Sensitivity Analysis of the Geometric Configuration in the Dynamic Setting 565.6.1 Initial Crack Length . . . . . . . . . . . . . . . . . . . . . . . . . . . 565.6.2 Domain Length . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 575.6.3 Domain Width . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 586 Simulations of Quasi-static Crack Propagation . . . . . . . . . . . . . . . . 616.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 616.2 General Comparisons between the Models . . . . . . . . . . . . . . . . . . . 616.2.1 Crack Profile . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 616.2.2 Stress and Elastic Energy Profiles in the Crack Plane . . . . . . . . . 63vii6.2.3 Energy Release Rate and Critical Energy Release Rate . . . . . . . . 666.3 Influence of the Irreversibility Criterion . . . . . . . . . . . . . . . . . . . . 676.4 Influence of the Residual Stiffness . . . . . . . . . . . . . . . . . . . . . . . . 726.5 Influence of the Crack Tip Curvature . . . . . . . . . . . . . . . . . . . . . . 746.6 Influence of the Input Elastic Energy . . . . . . . . . . . . . . . . . . . . . . 766.7 Interaction between the Crack and a Grain Boundary . . . . . . . . . . . . 787 Simulations of Dynamic Crack Propagation . . . . . . . . . . . . . . . . . 827.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 827.2 Investigations of the Branching Phenomena . . . . . . . . . . . . . . . . . . 837.2.1 Branching Velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . . 837.2.2 Branching Angle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 867.2.3 Branching Energy Release Rate . . . . . . . . . . . . . . . . . . . . . 867.3 Influence of the Input Elastic Energy . . . . . . . . . . . . . . . . . . . . . . 888 Summary and Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 928.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 928.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94AppendicesA Software Source Codes for the Mechanical Model . . . . . . . . . . . . . . 105A.1 The Mechanical Model for 2D Quasi-static Fracture . . . . . . . . . . . . . 105A.2 The Mechanical Model for 2D Dynamic Fracture . . . . . . . . . . . . . . . 123B Software Source Codes for the KKL Model . . . . . . . . . . . . . . . . . . 127B.1 The KKL Model for 2D Quasi-static Fracture . . . . . . . . . . . . . . . . . 127B.2 The KKL Model for 2D Dynamic Fracture . . . . . . . . . . . . . . . . . . . 134viiiList of Tables4.1 Parameters in 2D quasi-static simulations . . . . . . . . . . . . . . . . . . . 385.1 Geometry parameters in 2D dynamic simulations . . . . . . . . . . . . . . . 515.2 Geometry parameters for three cases . . . . . . . . . . . . . . . . . . . . . . 597.1 Parameters in 2D dynamic simulations . . . . . . . . . . . . . . . . . . . . . 83ixList of Figures2.1 Schematics of three fracture modes: (a) mode I, (b) mode II and (c) mode III. 42.2 The schematic description of the maximum tensile stress criterion. . . . . . . 62.3 Schematic profiles for (a) the sharp-interface model and (b) the diffuse-interfacemodel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.4 1D analytical crack (φ) profiles with different values of the internal length scale(l) generated from Equation (2.10). . . . . . . . . . . . . . . . . . . . . . . . 132.5 1D homogeneous solutions for (a) the phase-field parameter (φ) and (b) thestress (σ) generated from Equations (2.18) and (2.19). . . . . . . . . . . . . . 162.6 1D homogeneous solutions for (a) the phase-field parameter (φ), and (b) thecorresponding energy density (E) with cubic and quartic degradation function. 182.7 (a) Crack deflection and penetration behaviours at the grain boundary (blackline) with different grain boundary angles (ϕint). (b) Schematics of the influ-ence of the bulk toughness (Gbulkc ) on the interface toughness (Gintc ). . . . . . 232.8 The 1D analytical phase-field profile in the KKL model generated from Equa-tion (2.33). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 282.9 (a) The simulation setup in which a boundary with an angle θ separating twomaterials with I being the weak material with toughness ΓI and II being thetough material with toughness ΓII . (b) Comparison between the critical load(∆y/∆yc) for continuous crack propagation at the boundary from phase-fieldsimulations and that from the theoretical predictions with different boundaryangles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 304.1 Simulation setups for (a) the snap-back of a 1D elastic band and (b) crackpropagation in a 2D rectangular sheet. . . . . . . . . . . . . . . . . . . . . . 334.2 Schematics of the initial guess (solid line) and the analytical solution (dottedline) for the φ profile at interface. . . . . . . . . . . . . . . . . . . . . . . . . 374.3 φ profiles with two boundary conditions in both models. . . . . . . . . . . . 39x5.1 Crack interface profiles and stress distributions in 1D simulations with differentφthreshold and uthreshold in both models. . . . . . . . . . . . . . . . . . . . . . 415.2 (a) Critical ∆φthreshold and (b) critical ∆uthreshold for the results to convergeunder different mesh sizes in the mechanical and the KKL model. . . . . . . 425.3 Crack interface profiles and stress distributions in 1D simulations with differentmesh sizes in both models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 435.4 The schematic representation of the simulation domain in 2D. . . . . . . . . 445.5 φ profiles and elastic energy profiles in the crack plane with different meshsizes in both models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 455.6 φ profiles and elastic energy profiles in the crack plane with different time stepsizes in both models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 465.7 Elastic energy release and surface energy increment with respect to the crackpropagation length for (a) the mechanical model and (b) the KKL model withdifferent φfront. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 485.8 φ profiles (a) at the right half of the interface and (b) in the crack plane inthe first few computational loops in the mechanical model. . . . . . . . . . . 495.9 φ profiles (a) at the right half of the interface and (b) in the crack plane inthe first few computational loops in the KKL model. . . . . . . . . . . . . . 505.10 Crack velocity with respect to the propagation length with different φfront for(a) the mechanical model and (c) the KKL model. Phase-field profiles in thecrack plane at different times for (b) the mechanical model and (d) the KKLmodel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 525.11 The contour line of φ = 0.5 in the KKL model. . . . . . . . . . . . . . . . . . 535.12 Contour lines of φ = φfront for (a) the mechanical model and (b) the KKLmodel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 535.13 Elastic energy release and surface energy increment with respect to the crackpropagation length for (a) the mechanical model and (b) the KKL model withdifferent domain sizes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 555.14 Crack velocity with respect to the propagation length with different initialcrack lengths for (a) the mechanical model and (b) the KKL model. . . . . . 575.15 Crack velocity with respect to the propagation length with different domainlengths for (a) the mechanical model and (b) the KKL model. . . . . . . . . 585.16 Crack velocity with respect to the propagation lengths with different domainwidths for (a) the mechanical model and (b) the KKL model. . . . . . . . . 59xi5.17 Contour lines of φ = 0.5 with different domain widths for (a) the mechanicalmodel and (b) the KKL model. . . . . . . . . . . . . . . . . . . . . . . . . . 606.1 (a) Right half of the crack interface profiles and (b) crack profiles in the crackplane for both the mechanical model and the KKL model. . . . . . . . . . . 626.2 Elastic energy profiles in the crack plane for both the mechanical model andthe KKL model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 636.3 The shear stress in the y-direction (σzy) in the crack plane for both the me-chanical model and the KKL model. . . . . . . . . . . . . . . . . . . . . . . 646.4 Comparisons of the elastic energy profiles in the crack plane with the analyticalsolution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 666.5 (a) Elastic energy release and (b) surface energy increment with respect to thecrack propagation length for the mechanical model and the KKL model. . . 676.6 φ and elastic energy profiles in the crack plane for both models with differentirreversibility criteria. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 696.7 Right half of the crack interface profiles for (a) the mechanical model and (b)the KKL model with different irreversibility criteria. . . . . . . . . . . . . . . 706.8 Elastic energy and surface energy with respect to the propagation length for(a) the mechanical model and (b) the KKL model with two irreversibility criteria. 716.9 φ and elastic energy profiles in the crack plane along the propagation lengthfor both models with different residual stiffness. . . . . . . . . . . . . . . . . 736.10 Average number of iterations solving the mechanical equilibrium equation withrespect to the residual stiffness in the mechanical model and KKl model. . . 746.11 The schematics of the crack tip shape. . . . . . . . . . . . . . . . . . . . . . 756.12 Elastic energy profiles in the crack plane in the beginning and shortly afterpropagation with different κcur for (a) the mechanical and (b) KKL model. . 756.13 Contour lines of φ = 0 with different κcur in (a) the mechanical model and (b)the KKL model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 766.14 Elastic energy release and surface energy increment with respect to the cracktip position. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 776.15 Crack path at the boundary with different Gintc to Gbulkc ratios. . . . . . . . . 797.1 (a) The crack velocity with respect to the crack tip position and (b) the contourof φ profile with level set to 0.5 for both models. . . . . . . . . . . . . . . . . 85xii7.2 (a) The elastic energy profiles and (b) φ profiles in the crack plane for themechanical model and the KKL model. . . . . . . . . . . . . . . . . . . . . . 857.3 Contour lines of φ = 0.5 in the mechanical model and the KKL model. . . . 867.4 The critical energy release rate with respect to the crack tip position. . . . . 887.5 The crack velocity with respect to the crack tip position in the (a) mechanicalmodel and the (b) KKL model with different Ein. . . . . . . . . . . . . . . . 897.6 Contour lines of φ = 0.5 with different Ein in the (a) mechanical model and(b) the KKL model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 907.7 The critical energy release rate with respect to the crack tip position withdifferent Ein for the (a) mechanical model and the (b) KKL model. . . . . . 90xiiiList of SymbolsLatin LettersA Cross-sectional area of an elastic banda Crack extension lengthC Stiffness tensorC Young’s modulus in 1D and shear modulus in 2DCr Rayleigh wave speedCs Shear wave speedd Interpolation function for the energy related to strain hardeningdx Mesh sizedt Time step sizeE Young’s modulusEsurface Surface energyEelastic Elastic energyE ′ Parameter that is related to Young’s modulus and Poisson’s ratioEin Elastic energy input in the systemEc Critical elastic energy density to break the materialF Free energy functionalfinterfacial Density functional of the interfacial energyfelastic Density functional of the elastic energyfplastic Density functional of the energy related to strain hardeningfdissipation Density functional of the dissipated energy because of plasticityG Energy release rateGM Energy release rate in the mechanical modelGK Energy release rate in the KKL modelGc Critical energy release rateGMc Critical energy release rate in the mechanical modelGKc Critical energy release rate in the KKL modelxivG1.5 Energy release rate when Ein = 1.5Gc.G2.0 Energy release rate when Ein = 2.0Gc.G1.5c Critical energy release rate when Ein = 1.5Gc.G2.0c Critical energy release rate when Ein = 2.0Gc.Gbulkc Critical energy release rate in the bulkGintc Critical energy release rate in the grain boundaryGdc Critical energy release rate in dynamic fractureGθ Energy release rate as a function of kinking angle θg Degradation functionh Geometrical functionK Stress intensity factorKc Fracture toughnessKI Stress intensity factor for mode I fractureKII Stress intensity factor for mode II fractureKIII Stress intensity factor for mode III fractureL Length of an elastic bandLx Length of the 2D sheetLy Width of the 2D sheetl Internal length scalelc Initial crack lengthM Mobilitym Factor in the double-well functionn Normal vectorp Interpolation function for the dissipated energy because of plasticityr Distance from the crack tiprtip Crack tip radiusS Strain-energy density factort TimeU Elastic energyu Displacement fieldV Double-well functionx PositionxvGreek LettersΓl Numerical crack surface area∆ Boundary displacementδ Crack interface thicknessε Strain tensorεc Critical strain for crack initiationη Residual stiffnessθ Kinking angleκ Coefficient before the gradient term in the KKL modelκcurv Crack tip curvatureλ Branching angleµ Shear modulusν Poisson’s ratioρ Material densityσ Stress tensorσi,j Stress at (i, j) in the polar coordinate systemσtensile Tensile stressσc Failure stressφ Phase-field parameterφfront Definition of the crack front positionϕint Grain boundary anglexviList of AbbreviationsBCL Battelle Columbus LaboratoriesFDM Finite difference methodGDP Gross national productKKL Karma-Kessler-LevinLEFM Linear elastic fracture mechanicsMD Molecular dynamicsxviiAcknowledgementsI cannot begin to express my gratitude to my supervisor Prof. Matthias Militzer for givingme this exciting project and providing me the encouragement and guidance throughout theproject. I would also like to extend my sincere thanks to Dr. Michael Greenwood fromCanmetMATERIALS for his constructive advice and invaluable discussions. I am alwaysgrateful to have such a group of supportive co-workers in the microstructure group. Dr.Hocine Lebed never wavered in his support whenever I have problems in the research. AyushSuhane and Hariharan Umashankar are always willing to have discussions with me to clear mydoubts. I am also thankful to Tianbi Zhang who taught me to use the departmental cluster,MICRESS, Thermo-calc, and kindly provided me the template for this thesis. Special thanksto Sabyasachi Roy who gave me valuable suggestions about thesis writing.Finally, I am deeply indebted to my parents for their love and support.xviiiChapter 1IntroductionAccording to a 1983 report by the Battelle Columbus Laboratories (BCL) [6], the annualcost related to fracture in the US is $119 billion, which is about 4% of the GDP (grossnational product). A tremendous amount of money is spent each year on safety inspection,maintenance, repair, and replacement of the structures. Apart from the economic impact,fracture also leads to numerous catastrophes, e.g., the breakage of the liberty ships in theNorth Atlantic, the crashes of the first commercial jet airliner (de Havilland Comet), thefailure of the Aloha Airlines Flight 243, etc, which result in the loss of human lives. Theseevents of failure happen due to various reasons, e.g., overloading, defects in materials, lackof maintenance, deficiencies in structural design, and environmental effects. Thanks to theefforts of engineers and researchers, people nowadays have a much deeper understanding offracture mechanisms, which enables safer designs with fewer failures. Meanwhile, fracturemechanics becomes one of the most important subjects in mechanics with billions of dollarsput into research each year [6]. However, there are still challenges in fracture mechanics,which require collaboration between researchers in multiple fields, i.e., physics, mechanics,and materials.Understanding the fracture mechanism has always been a challenge as different materialsoften exhibit distinct fracture phenomena. Depending on whether the material goes throughplastic deformation, fracture can be categorized into brittle fracture with little or no plasticdeformation and ductile fracture with plastic deformation. Brittle fracture is often observedin ceramics. Meanwhile, ductile fracture can be found in most metals. Compared with duc-tile fracture, brittle fracture is relatively well-understood, although the crack nucleation inbrittle fracture is still a challenging topic. The existing theories can predict some importantfracture phenomena, e.g., critical condition for propagation, crack path, and velocity. How-ever, these theories often utilize assumptions that may not be applicable in reality. Moreover,the simplified settings adopted in theories lead to difficulty in predicting fracture phenomenain practice, as the settings are normally more complex. As a result, computational modellingas an alternative way to study fractures has gained popularity in recent years. Since some1of the computational models are based on theories, they are the perfect tools to verify theo-ries and understand the fracture mechanisms. Besides, the flexibility of applying the modelsto arbitrary settings enables comparisons with experiments and the predictions of fracturephenomena.The phase-field model as a popular model for simulating phase transformation has beenapplied to solve fracture problems for two decades in both the mechanics and the physicscommunity. In both communities, the phase-field models were successfully applied to brittlefracture and showed some promising results. Since these models appear to be similar withonly slight differences in formulations, it would be interesting to know the differences betweenthem. However, there have not been any detailed comparisons between these models in theliterature. Therefore, the main objective of the work is to compare phase-field fracture modelsin the mechanics and the physics community to understand their strengths and limitations,so that researchers who are new to the field can have a better idea of which model to use.This thesis is organized as follows: Chapter 2 provides a literature review, in which thebasics of the fracture theories and some fracture models are briefly covered, along with a de-tailed review of the phase-field fracture models in the mechanics and the physics community.Chapter 3 presents the scope and the objectives of the work. In Chapter 4, the methodol-ogy including the equations and simulation settings is given. To derive the proper rangesof the model parameters, sensitivity analyses are presented in Chapter 5. Chapters 6 and 7focus on the model applications, in which the results are compared with fracture theories. InChapter 8, a summary and an outlook into future work are presented.2Chapter 2Literature Review2.1 IntroductionIn this chapter, some relevant fracture theories and models are reviewed. The first part of thechapter introduces the fracture theories for both quasi-static and dynamic fracture. Then,some of the popular fracture models are briefly overviewed. The latter part of the chapteris dedicated to an extensive review of the phase-field fracture models in two communities:the mechanics and the physics community. Although the chapter focuses on brittle fracture,phase-field fracture models with plasticity are also reviewed as pure brittle fracture is hardlyseen in engineering materials.2.2 Crack Propagation in Quasi-static Brittle FractureCrack propagation is a complex phenomenon that exhibits distinct patterns under differentloading conditions. In general, the loading direction can be arbitrary and induces complexstress states in the body, i.e., tensile and compressive stress and in-plane and out-of-planeshear stress, which are difficult to analyze. Therefore, to simplify the analysis, fracture isnormally identified by three modes depending on the loading direction:• Mode I or the opening mode where the load is normal to the crack plane (Figure 2.1a).• Mode II or the sliding mode where the load is parallel to the crack plane and normalto the crack front (Figure 2.1b).• Mode III or the tearing mode where the load is parallel to the crack plane and thecrack front (Figure 2.1c).From the definitions of fracture mode, it is obvious that fracture is a 3D problem. However,because of the complexity in 3D stress analysis, most fracture theories and simulations arebased on simplified 2D scenarios with some assumptions. Among these assumptions, threeof them are most popular:3(a) (b) (c)Figure 2.1: Schematics of three fracture modes: (a) mode I, (b) mode II and (c) modeIII. The black arrow indicates the loading direction.• Plane strain condition where the plate is thick such that the strain is assumed to bezero in the thickness direction.• Plane stress condition where the plate is thin such that the stress is assumed to be zeroin the thickness direction.• Antiplane shear condition where only the displacement in the thickness direction isnon-zero and does not vary within the thickness direction.Although distinct from other conditions, the antiplane shear condition is rarely observedin reality, it was originally designed and widely used to study the 2D mode III fracture intheories.Apart from the loading direction, fracture can also be categorized based on the loadingspeed, i.e., quasi-static and dynamic fracture. In quasi-static fracture, the loading rate isextremely slow such that the stress is assumed to be in equilibrium. Meanwhile, the crackfront moves slowly, and the velocity is assumed to be zero. When the actual loading rate andcrack velocity are zero instead of being close to zero as in the quasi-static case, it becomes thestatic fracture. Although these two types of fracture are different in experiments, they aretreated the same way in theories and simulations, as they are all governed by the mechanicalequilibrium equation. Therefore, quasi-static and static fractures are not distinguished inthis thesis. In dynamic fracture, however, stress waves are utilized to characterize the rapidloading, and the inertia effect due to the fast crack movement is taken into consideration [7].4The fundamental theories for quasi-static brittle fracture can be dated back to the seminalstudies by Griffith [3] and Irwin [8], which are considered as the foundations of the linearelastic fracture mechanics (LEFM). Griffith assumed that the elastic energy release balancesthe surface energy increment at the crack tip. To identify the amount of elastic energy releaseper unit area of the newly created crack surface, the energy release rate (G) was introducedand is given as [3]:G = −∂U∂a(2.1)where U is the elastic energy, and a is the crack surface area which is the crack length in 2D.The unit of G is J/m2. When G is larger than the critical value (Gc) known as the criticalenergy release rate, the crack will propagate. To enable more thorough investigations intobrittle fracture, Irwin derived the asymptotic stress solution near the crack tip. In the polarcoordinate, the near-tip stress reads [9]:σij =K√2pirfij(θ) (2.2)where σij is the stress, K is the stress intensity factor, r is the distance from the crack tip,fij(θ) is a dimensionless quantity related to the load and geometry. The unit of K is Pa√m.Based on the stress solution, Irwin suggested that the crack propagates when K reachesthe critical value (Kc) known as the fracture toughness. Although the Griffith theory is anenergetic approach, while the stress intensity factor is related to the near-tip stress field, theyare equivalent. The relation of G and K in 2D is given by [9]:G =K2IE ′+K2IIE ′+K2III2µ(2.3)with:E ′ =E plane stressE1−ν2 plane strain(2.4)where KI , KII and KIII are the stress intensity factors in three fracture modes, E′ is relatedto Young’s modulus (E) and Poisson’s ratio (ν), and µ is the shear modulus.Both the Griffith theory and Irwin’s formulations are powerful tools in terms of predict-ing the start of crack propagation. However, they are unable to describe the crack path.Therefore, many theories [10–12] were developed to predict the crack propagation directionby extending either the Griffith theory or Irwin’s formulations. These criteria can be dividedinto two categories: the local criterion [10, 11] in which the propagation is determined by5the crack tip fields, and the global criterion [12] which is based on the energy balance of thesystem.As for the local criterion, either the stress or elastic energy field at the crack tip regionis evaluated. By comparing the stress at the crack tip, Erdogan and Sih [10] introducedthe maximum tensile stress criterion in which they stated that the crack propagates intothe direction normal to the maximum tensile stress. Due to its simplicity, this criterion waswidely used in the mechanics community [13–15]. In numerical applications, the tensile stressis evaluated along a circle centred at the crack tip with a radius r (Figure 2.2) [16]. Since theaccuracy of the kinking angle depends on the number of points along the circle, this methodis sensitive to the mesh size at the crack tip. Besides, to deal with the stress singularityat the crack tip, great efforts have been made in the past few decades to develop numericaltools, e.g., the re-meshing technique [15], the extended finite element method [13, 14], etc.σtensilercrackFigure 2.2: The schematic description of the maximum tensile stress criterion. Thesolid arrow indicates the direction of the maximum tensile stress (σtensile) alongthe circle. The dashed arrow indicates the crack propagation direction which isperpendicular to the direction of the maximum tensile stress.Arguing that the maximum tensile stress criterion becomes invalid when r = 0, as thesingular stresses cannot be compared, Sih and Macdonald [11] introduced the minimumstrain energy density criterion. They defined a term called the strain-energy density factor(S), which is the elastic energy density multiplied by the distance to the crack tip (r). Sincethe stress has 1/√r singularity (Equation (2.2)) and the energy density has 1/r singularity,by multiplying the energy density with r, S does not have a singularity even at the cracktip. Then they claimed that the crack should propagate in the direction with the minimumstrain energy density factor. However, since the criterion is implemented in the same wayas the maximum tensile stress criterion and the stress singularity still exists, it also requiresusing complicated numerical methods.The mesh dependency is the main problem in the models utilizing the local criteria. Themaximum strain energy release rate criterion as a global approach, however, does not have6such a problem. With this criterion, the crack is supposed to propagate in the direction withthe maximum strain energy release rate. Simple as it may seem, the numerical implementa-tion is complicated. The most straightforward approach is the real crack extension method[13]. In this method, the crack is manually extended into different directions with small dis-tances, after which the energy release rates in those directions are evaluated. Although thismethod is easy to implement, it suffers from low computational efficiency as the mechanicalequations are calculated many times in different directions for each small crack extension[16]. To address this problem, several methods have been developed, e.g., the virtual crackextension method and the Gθ method. In the virtual crack extension method [17–19], theanalytical solutions for the energy release rate in different directions are approximated withthe finite element method. Therefore, instead of numerically solving the mechanical equa-tions in all directions, the direction with the maximum energy release rate can be determinedfrom the analytical solutions. Similarly, in the Gθ method [20], the energy release rate is for-mulated as a function of the kinking angle (θ) with certain assumptions, and θ is derivedby minimizing the function Gθ. Although these methods are much more efficient, they areconceptually more complicated than the real crack extension method. Besides, the analyticalrelations between the energy release rate and the kinking angle may not be accurate.2.3 Basic Phenomena in Dynamic Brittle FractureDistinct from quasi-static fracture, the crack velocity is an important feature in dynamicfracture. The first calculation of the crack velocity based on continuum mechanics was givenby Freund [21] who studied the propagation of a semi-infinite crack in an infinitely largedomain. The result suggested that the limiting velocity for dynamic crack propagation is theRayleigh wave speed (Cr). However, various experiments showed that the limiting velocityis around 0.5Cr, and the value depends on the materials and the loading condition [7].The discrepancy in the limiting velocity may come from the straight crack path assumptionmade in Freund’s analysis, as the crack branching can be observed in dynamic fracture inmany materials. Therefore, one possible reason for the low limiting velocity observed inthe experiments is the breakdown of Freund’s theory. However, even today, a satisfyingexplanation for the branching behaviour is still lacking.Yoffe [22] was the first to perform the continuum mechanics analysis on the branchingphenomenon. In the study, Yoffe analyzed the stress distribution of a crack translating witha constant velocity in an infinitely large medium and adopted the maximum tensile stress7criterion. From the analytical solution, it was observed that when the crack speed reaches0.6Cr, the rearrangement of the stress at the crack front leads to the crack bifurcation with thebranching angle equal to 60°. Eshelby [4] then argued that since the surface energy doublesafter crack bifurcation, to maintain the total energy of the system, the speed after branchingshould decrease. Thus, Eshelby assumed that the velocity is zero after branching. Then,from the energy balance between the energy flux into the crack tip and the surface energy,the critical branching velocity was found to be around 0.6Cr in the limiting case, where thebranching angle is zero. However, Eshelby’s result was considered as a rough estimate as nosolution was provided for arbitrary branching angles [5]. By adopting Eshelby’s approach,Adda-Bedia [5] then studied the bifurcation of a semi-infinite crack in an infinitely long stripunder the antiplane shear condition. The branching angle and critical branching velocitywere shown to be 79.2° and 0.39Cs, where Cs is the shear wave speed.The branching angle and velocity by Adda-Bedia [5] agree qualitatively well with exper-iments. However, Ravi-Chandar [7] argued that the crack speed after branching does notdecrease significantly from experimental observations, which suggests that the zero-velocityassumption by Adda-Bedia and Eshelby is inappropriate. From the micrographs of the crackfronts, Ravi-Chandar suggested that the crack branching is a result of the interaction ofmicrocracks in a small zone ahead of the crack tip (process zone), which cannot be explainedproperly from the continuum mechanics scale. Since the LEFM is the foundation of thetheories and models in this project, the continuum mechanics theories are the focuses of thecurrent work. A detailed explanation of other possible branching mechanisms can be foundin [7]. Although understanding the branching mechanism remains a challenge, experimentaland simulation investigations will continuously help unveil the mystery behind it.2.4 Brief Overview of the Fracture ModelsExperiments play a critical role in understanding fracture mechanisms. However, the dy-namics of the fracture process is challenging to investigate experimentally, since the crackpropagates very fast, i.e., several thousands of meters per second [23]. Therefore, theoriesand models are required to have a thorough understanding of the fracture mechanisms. Sincethe macrocrack originates from the coalescence of microcracks, or the separation of atoms atan even smaller scale, the atomic-scale models that can deal with the motion of atoms seemto be ideal to study the mechanisms. Molecular dynamics (MD), as one of the most popularatomistic models, utilizes Newton’s equations of motion to describe the motion of atoms and8has successfully been used to reproduce many fracture phenomena, e.g., interaction of thecrack with dislocation, atom configuration, void, etc. Even though with the cheaper andmore powerful computational resources nowadays, MD simulations can deal with systemsof up to several million atoms, they are still too small to have enough complexity in thestructure [23].Therefore, atomistic models have to be coupled with models in larger scales, i.e., mesoscale(100nm to 100µm) and macroscale (> 0.5mm). Depending on the treatment of the crackinterface, these models can be divided into two categories: sharp-interface model where theproperties, e.g., mass, density, etc, go through an abrupt transition at the crack interface(Figure 2.3a) and diffuse-interface model where the properties vary smoothly at the inter-face (Figure 2.3b). Most of the sharp-interface models rely on various propagation criteriadiscussed in Section 2.2. Because of the stress singularity at the crack tip, they suffer fromthe mesh sensitivity issue. As for the diffuse-interface models, the smooth interface removesthe singularity at the crack tip. As a result, they are relatively insensitive to the mesh size.Although the diffuse-interface crack description and the non-singular crack tip stress deviatefrom the LEFM, the diffuse-interface models have great advantages in capturing complexcrack patterns, especially in 3D.(a) (b)Figure 2.3: Schematic profiles for (a) the sharp-interface model and (b) the diffuse-interface model.Among the diffuse-interface models, the phase-field fracture models have gained popu-larity in both the mechanics and the physics community over the past two decades. Thephase-field model was originally designed to describe solidification. The first successful ap-9plication of the phase-field concept was made by Kobayashi [24] who studied dendrite growthin solidification. The ability to reproduce the complex morphology of dendrites with simpleformulations and low computational cost, then made the phase-field model emerge to be apopular computational tool to solve problems in different fields, i.e., solidification [24–27],recrystallization [28–30], solid-state phase transition [31–34], crack propagation [2, 35–37],etc.The beauty of the phase-field model lies within its simple and elegant formulations. In themodel, a phase-field variable (φ) is related to physical properties, e.g., density, concentration,grain, phase, etc, which vary smoothly at the interface. Then, the interface migration isgoverned by the partial differential equations referred to as the phase-field equations. Fornon-conserved variables, e.g., grain, phase, etc, the phase-field equation is described by theAllen-Cahn equation as [38]:∂φ∂t= −MδFδφ(2.5)where M is the mobility related to the interface migration speed, and F is the free energyfunctional which contains the interfacial energy, bulk free energy, elastic energy, etc. As forthe conserved variables, e.g., concentration, density, etc, the phase-field equation is describedby the Cahn-Hilliard equation as [38]:∂φ∂t= ∇ ·M∇δFδφ(2.6)By coupling different fields and adding the corresponding energy terms into the free energyfunctional, the phase-field model can be applied to a wide variety of problems.To describe brittle fracture, the phase-field parameter is used to identify the state of thematerial, i.e., intact body or crack. Depending on the models, φ can either be conserved ornon-conserved. Therefore, both the Allen-Cahn and the Cahn-Hilliard equation were used inthe literature.102.5 The Phase-Field Fracture Models in theMechanics Community2.5.1 Basics of the Mechanical ModelsVarious criteria presented in Section 2.2 were developed to predict the crack path, and theywere widely used in finite element simulations. However, most of them need to be appliedrepetitively at every small crack increment, and the accuracy largely depends on the meshsize.With the global energy minimization assumption, Francfort and Marigo [39] proposed thevariational formulation for quasi-static crack propagation in brittle materials. They showedthat by minimizing the potential of the system, the crack path can be determined withoutadditional criteria. Besides, compared with the damage models [40–42], the crack initiationcan be captured by the model without having empirical relations. As demonstrated by variousnumerical examples [39], the model can predict crack initiation and propagation. However,it cannot simulate crack propagation under arbitrary loadings [39]. Besides, the definition ofthe crack as a discontinuity adds complexity to the numerical implementation.In the image segmentation model, Ambrosio and Tortorelli [43] approximated the originaldiscrete functional by Mumford and Shah [44] with a continuous energy functional. Since thevariational formulation closely resembles the Mumford-Shah energy functional, inspired byAmbrosio and Tortorelli [43], Bourdin et al. [45] replaced the discontinuous crack with thediffuse-interface crack in the variational formulation. In the approximated formulation, anorder parameter φ was introduced to represent the crack with φ = 0 inside the crack and φ = 1in the bulk. With the Γ-convergence technique developed by Ambrosio and Tortorelli [43],the continuous energy functional was shown to converge to the variational formulation whenthe diffuse-interface crack converges to the discontinuous crack with a vanishing interfacethickness. The continuous energy functional by Bourdin et al. is given by:F =∫Ω(finterfacial + felastic) dΩ (2.7)where finterfacial and felastic are the density functionals of the interfacial energy and the elasticenergy. The interfacial energy density is given as [45]:finterfacial = Gc(l(∇φ)2 + h(φ)l)(2.8)11with [45]:h(φ) =(1− φ)24(2.9)where l is the internal length scale that controls the crack interface thickness, and h(φ) isthe geometrical function that determines the shape of the crack interface. The whole termnext to Gc in Equation (2.8) is called the crack surface density function. With this form ofthe crack surface density function, the 1D crack interface profile is given by [46]:φ = 1− e− |x|2l (2.10)As shown in Figure 2.4, with l approaching 0, the diffuse-interface crack converges to thediscontinuous crack. The elastic energy density reads [45]:felastic = (g(φ) + η)12ε(u) : [Cε(u)] (2.11)where u is the displacement field, C is the stiffness tensor, and ε(u) is the infinitesimal straintensor given by:ε(u) =12(∇u+ (∇u)T ) (2.12)where T denotes the transpose of a tensor. 12ε(u) : [Cε(u)] is the elastic energy density inthe intact body, and g(φ) is the degradation function interpolating the elastic energy densityin the crack (φ = 0) and the bulk (φ = 1). The degradation function g(φ) based on theAmbrosio-Tortorelli functional is defined as [45]:g(φ) = φ2 (2.13)η next to the degradation function in Equation (2.11) is the residual stiffness, which is a non-zero constant to ensure numerical stability [45]. To avoid the overestimation of the elasticenergy, the value of it should be sufficiently small, e.g., less than 10−6. Through optimizingthe displacement field and the crack field to minimize the energy functional in Equation(2.7), both the stress distribution and the crack path are determined. The optimization wasachieved utilizing the alternative minimization algorithm [45] in which at each time, oneparameter is optimized while the other parameter is fixed. To prevent the unphysical self-healing behaviour of the crack, φ was not allowed to evolve once it becomes 0 [45]. Comparedwith the variational formulation, Bourdin et al.’s model [45] is easier to implement and hashigher computational efficiency.Kuhn and Mu¨ller [36, 47] then adopted the energy functional by Bourdin et al. [45].Instead of using the alternative minimization algorithm, the crack field was governed by the12Figure 2.4: 1D analytical crack (φ) profiles with different values of the internal lengthscale (l) generated from Equation (2.10).Ginzburg-Landau evolution equation. Combining Equations (2.7), (2.8), (2.11), the phase-field equation reads [36]:1M∂φ∂t= −δFδφ= 2Gcl∇2φ+ Gc2l(1− φ)− φε : [Cε] (2.14)Since the stress is assumed to be in the equilibrium state in quasi-static fracture, the displace-ment field is governed by the mechanical equilibrium equation. By ignoring the body force,which is the force acting on the volume of the body, e.g., gravitational force, electromagneticforce, etc, the equation reads:∇ · σ = 0 (2.15)where σ is the stress derived from the elastic energy density in Equation (2.11) given as [36]:σ =∂felastic∂ε= (φ2 + η)Cε (2.16)Then for the first time, the model was called a phase-field model for its close resemblanceto the phase-field models used to describe phase transformation. The order parameter forthe crack, therefore, was referred to as the phase-field parameter. However, since the energydissipation is assumed to happen instantaneously in brittle fracture, the time-dependent term1/Mdφ/dt could delay the dissipation process and lead to the deviation from the Griffiththeory. Kuhn and Mu¨ller [36] then showed that when the mobility is sufficiently large,the solutions are in good agreement with the Griffith theory. The non-conserved governing13equation for the phase-field parameter (Equation (2.14)) also makes it clear that the crackwas treated as a non-conserved field in the models by Kuhn and Mu¨ller [36] and Bourdinet al. [45]. As a result, the volume of the body decreases during crack evolution. However,since the crack volume is much smaller than the domain size in simulations, the volume lossduring crack evolution remains negligible.Based on continuum mechanics and thermodynamics, Miehe et al. [1, 48] proposed analternative continuous model in which the energy functional and the displacement equationare the same as Equations (2.7) and (2.15). Two types of governing equations for the phase-field parameter were introduced: the rate-dependent equation identical to Equation (2.14)and the rate-independent equation by setting M to infinity such that [1]:2Gcl∇2φ+ Gc2l(1− φ)− φε : [Cε] = 0 (2.17)Their model in essence is the same as those by Kuhn and Mu¨ller [36] and Bourdin et al. [45].However, the phase-field variable was viewed as a damage variable with φ = 0 as the fullydamaged state or the crack state and with φ < 1 as the partially damaged state containingmicrocracks. Therefore, instead of only constraining the evolution of φ = 0, the phase-fieldparameter is forced to decrease (dφ/dt ≤ 0) so that both the crack and the partially damagedarea cannot heal [1].Different irreversibility criteria, i.e., φ = 0 and dφ/dt ≤ 0, could potentially lead todifferent model responses. Therefore, Linse et al. [49] compared these two criteria by studyingthe uniaxial tension of a 1D elastic band. The accuracy of the results was evaluated bycomparing the crack surface area with the analytical value derived from the surface densityfunction in Equation (2.8). With φ ≤ 0, the crack surface area agrees well with the analyticalvalue when the internal length scale (l) to domain size ratio is less than 0.05. With dφ/dt ≤ 0,however, the error surprisingly increases to over 100% with the decrease of the l to domainsize ratio. They pointed out that the difference occurs only after the peak load since thepartially damaged area recovers to the intact state (φ = 1) when only constraining φ = 0,while with the other criterion, φ keeps decreasing. The different predictions of the surfacearea also result in a slightly different load-displacement response after peak load [49]. Basedon the accuracy of the surface area predictions, they then suggested that Miehe et al.’sirreversibility criterion, i.e., dφ/dt ≤ 0, could only be used before the peak load. However,the accuracy of the model cannot simply be determined from the 1D crack profile, since it isonly a numerical feature. The influence of the irreversibility criteria should be investigatedwith 2D simulations in which the energy release rate and the crack pattern can be evaluated.142.5.2 Interpretations of the Internal Length Scale (l)Internal Length Scale as a Regularization ParameterThe internal length scale (l) was interpreted by Bourdin et al. [45] as a regularizationparameter, which controls the convergence of the diffuse-interface crack to the discontinuouscrack. Therefore, theoretically, l should be zero such that the energy functional can convergeto the variational formulation. However, the smaller l requires a finer mesh to solve theequation, which leads to a higher computational cost.Internal Length Scale as a Fitting ParameterEver since Amor et al. [50] illustrated that l can be related to the failure strength of thematerial, it has been more accepted as a fitting parameter rather than a regularizationparameter. Amor et al. [50], Borden [51], and Pham and Marigo [52] studied the 1D uniaxialtension of an elastic band and provided the analytical stress-strain relation by assuming thehomogeneous damage (∇φ = 0) in the band. Combining Equations (2.15) and (2.17), thehomogeneous solutions for the phase-field parameter and the stress read [51]:φ =(2lEGcε2 + 1)−1(2.18)σ =(2lEGcε2 + 1)−2Eε (2.19)The phase-field parameter and the stress with respect to the strain are shown in Figure 2.5.The peak stress and the corresponding strain in Figure 2.5b are the failure strain and stress,which are given by [51]:εc =√Gc6lE, σc =916√EGc6l(2.20)The homogeneous solutions are stable only before the critical load as the damaged arealocalizes, and ∇φ = 0 is no longer valid. Since σc agrees perfectly well with the failurestrength obtained from both 1D and 2D uniaxial tension simulations [52–54], it can be usedto determine the value of l. Besides, Equation (2.20) shows that σc becomes infinite whenl approaches zero, which is consistent with the infinite strength predicted by the Griffiththeory in the defect-less material. Various studies [55–58] have shown that the interpretationof l as a fitting parameter enables the accurate representation of the crack path and theload-displacement response. For example, Nguyen et al. [55] performed 3D three-point15bending simulations with l calculated from the failure strength reported from experiments,and good agreement in the load-displacement curve was obtained compared with experiments.However, even viewing l as a fitting parameter, the l to domain size ratio should still besufficiently small to maintain accuracy. Zhang et al. [56] showed that although the load-displacement curve is accurate using a small domain size, the crack path deviates fromexperiments. Good accuracy of both the load-displacement curve and the crack path can beachieved only with a small internal length scale to domain size ratio.(a) (b)Figure 2.5: 1D homogeneous solutions for (a) the phase-field parameter (φ) and (b)the stress (σ) generated from Equations (2.18) and (2.19). εc and σc are thecritical strain and stress for the breakage of the 1D band.Until now, the interpretation of l as a fitting parameter is more popular, as it doesnot have to be arbitrarily small to converge to the discontinuous crack and can be relatedto experimental results. However, the connection of l with both the interface thicknessand the failure strength leads to the dependence of the failure strength on the interfacethickness. As pointed out by Zhang et al. [56], the l to domain size ratio still has to bevery small to converge to the variational formulation, regardless of whether l is interpretedas a regularization parameter or a fitting parameter. Therefore, Wu [59] developed a modelin which the failure strength is separated from l. Although the model has a similar energyfunctional as those proposed by Bourdin et al. [45] and Miehe et al. [1], the geometricalfunction (h(φ)) and the degradation function (g(φ)) are completely different. The modelremoves the dependence of the failure strength on l. As a result, the failure strength doesnot have to change, while the numerical convergence to the discontinuous crack can still be16guaranteed [59]. However, the model has not been evaluated extensively, and the limitationof the model is still unknown.2.5.3 Model VariationsThe models by Miehe et al. [1] and Bourdin et al. [45] share the same energy functionalwith only minor differences in the way they implement the model, e.g., different governingequations and irreversibility criteria. However, the models utilizing that energy functionalall showed the same problems. For example, the stress-strain response before the peak loadis nonlinear as shown in Figure 2.5b. Therefore, alternative forms of the energy functionalwere developed.Alternative Forms of the Geometrical Function (h(φ))The nonlinear stress-strain response is caused by the decrease of φ with the non-zero strainshown in Figure 2.5a, which is referred to as the damage accumulation in the literature. Thisunphysical phenomenon can also lead to strain energy release far from the crack [51]. Insteadof adopting the quadratic function (Equation (2.9)), Pham et al. [52] suggested using thelinear function, i.e., 9(1 − φ)/64, as the geometrical function. They showed that the stress-strain relation becomes perfectly linear, and φ = 1 before the failure stress. Besides, thecrack interface profile is described by the quadratic function, which does not spread to theentire domain like the exponential function in Equation (2.10). Therefore, the boundaryeffect can be less significant. However, the linear function leads to the numerical issue as φcan exceed the 0 to 1 interval [60]. As a result, φ has to be bounded manually. Borden [51]then suggested using the cubic function, which results in the similar linear elastic behaviourbefore the failure stress and the bounded phase-field parameter. However, since an additionalparameter is introduced in the cubic function, the failure stress becomes dependent on thatparameter [51].Alternative Forms of the Degradation Function (g(φ))An alternative way to improve the linear elastic behaviour is to modify the degradationfunction since it also enters in the homogeneous solution. The degradation function shouldbe constructed in a way such that g(0) = 0 and g(1) = 1. Meanwhile, the derivative of g(φ)with φ should be zero at φ = 0 to avoid the elastic driving term inside the crack as shown inEquation (2.14) or (2.17) [61]. Kuhn et al. [61] then compared three degradation functions17that satisfy these requirements: the quadratic function, the cubic function, and the quarticfunction.Homogeneous solutions with quartic and cubic function were given (Figure 2.6) in whichthe perfectly linear elastic behaviour and φ = 1 were observed before the failure stress.However, with the cubic and the quartic function, they found that the crack nucleationhappens only when the numerical perturbation is applied on the phase-field parameter. Thisis because there are two solutions: φ = 1 and the homogeneous solution, and for bothsolutions, φ and the gradient of E with φ at the critical strain (εcc or εqc) are the same asshown in Figure 2.6 [61]. Therefore, the system will stay in φ = 1 without perturbation.However, the influence of the perturbation was not evaluated.(a) (b)Figure 2.6: 1D homogeneous solutions for (a) the phase-field parameter (φ), and (b)the corresponding energy density (E) with cubic and quartic degradation func-tion. εcc and εqc are the critical strains for the breakage of the 1D band with thecubic and the quartic degradation function. Adapted from [61].Higher-order Gradient Term of Phase-Field ParameterTo improve the convergence rate of the diffuse interface crack to the discontinuous crack,Borden [51] introduced the fourth-order phase-field term into the surface density function.The interfacial energy density then becomes [51]:finterfacial = Gc(l34(∆φ)2 +l2(∇φ)2 + (1− φ)24l)(2.21)18At the same time, the 1D crack interface profile reads [51]:φ = 1−(1 +|x|l)e−|x|l (2.22)Compared with the original equation with the second-order phase-field term, the faster con-vergence rate of the interfacial energy to the analytical value (Gc) with the decrease of l wasobserved in 1D simulations. Besides, the high-order model gets rid of the residual stiffness[51]. However, the interface profile is still described by an exponential function that spreadsto the entire domain. As a result, although the convergence rate to the analytical value isfaster, the accuracy of the interfacial energy is still prone to the boundary effect. Besides,the Γ-convergence of the model was not available [51].2.5.4 Model ApplicationsCrack Initiation and Quasi-static Crack PropagationAlthough the interpretations of the parameters and the forms of the model can be different,these models are the same in essence and are capable of describing fracture phenomena.Among various fracture phenomena, crack nucleation has long been a difficult problem inthe mechanics community. Same as the variational formulation, the phase-field model isalso capable of simulating nucleation. The 1D homogeneous solution suggests that the cracknucleates when reaching the failure stress (Figure 2.5b). Meanwhile, φ decreases to 0.75during the nucleation (Figure 2.5a). Kuhn and Mu¨ller [54] studied the crack nucleation froma circular notch in a 2D sample. Compared with the homogeneous solution, they foundthat φ decreases to 0.8 during nucleation, which is quite close to 0.75. Besides, the criticalapplied load for crack nucleation observed in 2D simulations is consistent with the failurestress in the homogeneous solutions. Therefore, the critical condition for crack nucleation in2D phase-field simulations, i.e., critical load and decrease in φ, can be correctly describedby the 1D homogeneous solutions. This finding was further proved by Tanne´ et al. [62]with more common experimental setups in which the crack nucleates from U and V-notches.However, they also found that the damage accumulation far from the crack could potentiallylead to deviations of the critical load from experiments, which can be avoided by using thelinear geometrical function.Because of the simplified elastic energy expression in Equation (2.11), as pointed outby Amor et al. [50], the model is unable to distinguish between compression and tension.Therefore, Amor et al. [50] developed a so-called ‘anisotropic’ model [63] which accounts for19the asymmetric behaviour in tension and compression by decomposing the strain into thetensile strain, the compressive strain, and the shear strain, of which only the tensile strainand the shear strain contribute to crack propagation. Later on, similar ‘anisotropic’ models[64, 65] were proposed, of which the main difference is the way they decompose the strain.However, anisotropy normally refers to the direction-dependent mechanical and materialproperties. Since the anisotropic elasticity is often less significant than the anisotropic failurestrength [66], and it can be incorporated into the stiffness tensor in the constitutive equation,the anisotropic surface energy is worth more attention. By making the crack interfacialenergy direction-dependent, Li et al. [67, 68], Teichtmeister et al. [69], and Nguyen et al.[66] proposed anisotropic models where interesting crack patterns, i.e., zigzag crack path[66–69], and good agreement with experiments [66] were observed.Dynamic Crack PropagationAll those models were originally designed to study quasi-static crack propagation. However,with only a few changes, they can be extended to the dynamic case. In dynamic fracture, therelease of elastic energy also leads to an increase in kinetic energy. The energy functional,therefore, becomes [70]:F =∫Ω(12ρ(∂u∂t)2+Gc(l(∇φ)2 + h(φ)l)+ (g(φ) + η)12ε : [Cε])dΩ (2.23)where ρ is the material density, and the first term is the kinetic energy. Since ρ is the same inthe bulk and the crack, similar to the non-conversed volume discussed in Section 2.5.1, it alsodeviates from the reality. From the conservation of the linear momentum, the displacement-field equation reads [70]:ρ∂2u∂t2= ∇ · σ = ∇ · ((φ2 + η)Cε) (2.24)As discussed in Section 2.3, the mechanisms behind the limiting velocity, branching veloc-ity, and branching angle are still not well-understood. Therefore, to gain insights into themechanisms, many researchers tried to reproduce the experimental results with phase-fieldsimulations. For the limiting velocity, various phase-field studies [70–72] have reported it tobe around 0.6Cr, which agrees well with the commonly observed experimental results. Bleyeret al. [73] then showed that the crack accelerates to 0.9Cr when suppressing the branchingwith a tough material surrounding the crack path or with circular holes in the crack path.This finding suggests that the limiting velocity is the Rayleigh wave speed and is influenced20by the branching behaviour. The branching velocity from phase-field simulations, however,has not been reported in the mechanics community. Instead of focusing on the branchingvelocity, Bleyer et al. [73] showed that the crack starts branching when the energy releaserate reaches 2Gc, which is consistent with Eshelby’s energetic theory [4]. As for the branchingangle, Borden et al. [70] and Hofacker and Miehe [74] adopted the setting of the dynamicshear loading experiment by Kalthoff and Winkler [75]. Utilizing the phase-field model, theysuccessfully reproduced the crack path observed in the Kalthoff-Winkler experiment.Crack Propagation in Heterogeneous MaterialsMost applications of the phase-field fracture models were in the macroscale, in which thematerials are often assumed to be homogeneous with the same material properties acrossthe domain. In the mesoscale, however, the material properties can vary within the domainbecause of the complex microstructures. Therefore, to study the interaction between crackpropagation and microstructures, the critical energy release rate and the stiffness can nolonger be homogeneous.Hossain et al. [76] were the first to study crack propagation in a heterogeneous material.The model is the same as the homogeneous one, but the critical energy release rate and thestiffness become the functions of the position. Toughening of the material was observed witheither heterogeneous elasticity or critical energy release rate in their simulations with simplelaminar structures. Using the idea of Hossain et al. [76], Nguyen et al. [77, 78] then appliedthe model to more complex microstructures derived from experiments. In the simulations,different elastic properties were assigned to different phases. However, the critical energyrelease rate is assumed to be homogeneous in the domain. Through 2D and 3D simulations,they observed the inter-granular and trans-granular crack propagation behaviours.Schneider et al. [79] extended the multi-phase-field model by Nestler et al. [80] forphase transformation to brittle fracture. In the energy functional, additional terms wereadded to account for the interfacial energy of the phase boundaries. Since there is no phasetransformation in the process, only the crack interface was allowed to evolve [79]. For thefirst time, heterogeneous critical energy release rate, elasticity, and grain boundary weakeningeffect were all taken into account.However, all these simulations lack quantitative comparisons with either experiments orfracture theories. The first quantitative comparison was made by Hansen-Do¨rr et al. [81]who compared the results with the crack deflection theory by He and Hutchinson [82]. In21the theory, the crack deflects or penetrates the boundary depending on the grain boundaryangle and the interface toughness (Gintc ) to the bulk toughness (Gbulkc ) ratio as shown inFigure 2.7a. In Hansen-Do¨rr et al.’s model [81], the grain boundary toughness is lower thanthe bulk toughness and kept constant across the grain boundary. However, they pointedout that the energy release rate at the grain boundary is not equal to the grain boundarytoughness because of the interaction between the crack interface and the grain boundary.As shown in Figure 2.7b, since the crack interface spreads to the bulk, the bulk toughnessalso contributes to the energy release rate. Therefore, with a wider grain boundary or anarrower crack interface, the influence of the bulk toughness can be reduced. The model bySchneider et al. [79] can also capture the grain boundary effect. However, the grain boundaryenergy is implicitly incorporated into the model through the interfacial energy term for thephase boundary. As a result, the actual critical energy release rate at the grain boundary isunknown. Therefore, the comparison with the crack deflection theory becomes much moretricky to do. Chen et al. [83] also developed a phase-field model accounting for the grainboundary effect, which showed good agreement with the crack deflection theory. However,instead of having a constant critical energy release rate at the grain boundary, the grainboundary toughness was defined as the lowest value in the middle of the boundary, whichvaries smoothly to the bulk. As a result, the actual energy release rate at the grain boundaryis not the grain boundary toughness they claimed. Besides, the results seem questionableas the dynamic simulation was compared with the crack deflection theory, which is for thequasi-static case.22(a)(b)Figure 2.7: (a) Crack deflection and penetration behaviours at the grain boundary(black line) with different grain boundary angles (ϕint). (b) Schematics of theinfluence of the bulk toughness (Gbulkc ) on the interface toughness (Gintc ). In (b),the black solid curve is the crack surface profile, the red solid line indicates thewidth of the crack interface, and the red dotted line indicates the thickness of thegrain boundary. Adapted from [81]. Reprinted with permission from Elsevier.23Crack Propagation with Plastic DeformationThe phase-field fracture model is very sophisticated in terms of describing perfect brittlefracture. However, for most materials, plastic deformation accompanies fracture. Unlike thelinear stress-strain relation for brittle materials, beyond the yield stress, the stress becomeseither constant in a perfect plastic material or increases linearly or nonlinearly with strainbecause of strain hardening. Upon unloading, not all the strain is recovered due to thepermanent plastic deformation. Therefore, the strain can be decomposed into the recoverableelastic strain and the unrecoverable plastic strain. The potential of the system with plasticdeformation then reads [84]:F =∫Ω(felastic + fplastic + fdissipation) dΩ (2.25)where fplastic is the plastic free energy, and fdissipation is the dissipated energy caused byplasticity. Combining the energy functionals for fracture (Equation (2.7)) and plastic defor-mation (Equation (2.25)), several researchers [85–92] then proposed the phase-field modelsfor ductile fracture. The energy functionals used in these models can be generalized as [84]:F =∫Ω(fintefacial + felastic + p(φ)fplastic + d(φ)fdissipation) dΩ (2.26)where p(φ) and d(φ) are the coupling functions between fracture and plastic deformation.Most of the models considered a linear hardening behaviour, while Miehe et al. [91] gavea more general formulation for both linear hardening and nonlinear hardening. However, themain difference is how they couple plasticity.Distinct from all the other models, there is no coupling between the plastic deformationand crack propagation in the model by Duda et al. [86], which means that p(φ) and d(φ) areboth equal to one. As a result, crack propagation is only driven by the elastic energy, whichleads to the same crack propagation behaviour like that in brittle fracture.Similar to the model by Duda et al. [86], another unique model was proposed by Ambati etal. [87] in which p(φ) and d(φ) are equal to one. However, instead of employing the originaldegradation function (g(φ)) for brittle fracture like Duda et al. [86], g(φ) is modified tocouple plasticity. In the degradation function, they introduced a threshold value identifyingthe critical plastic deformation, and fracture is triggered when plastic deformation exceedsthe critical value.All the other models coupled the plastic deformation through p(φ) and d(φ) while g(φ)was kept the same as that in brittle fracture. Borden et al. [89], Kuhn et al. [92], and Miehe24et al. [91] suggested using the same quadratic function for g(φ), p(φ), and d(φ), while variousforms of the coupling functions were used in other studies [85–88, 90].The different choices of the degradation function lead to different responses of the model,e.g., plastic strain evolution, stress-strain diagram, etc [84]. In general, the coupling of crackpropagation with plasticity remains a challenge, and most of the models are in the theoreticalstage with no comparison with experiments. The only exception is the work by Ambati et al.[90], in which the necking behaviour was captured, and the load-displacement curves were ingood agreement with experiments.2.6 The Phase-Field Fracture Models in the PhysicsCommunity2.6.1 General Overview of the Phase-Field Fracture Models inthe Physics CommunityTo understand dynamic fracture phenomena, phase-field models have also become popular inthe physics community. Although the models share the same name as those in the mechanicscommunity, they are different in many aspects.The first phase-field fracture model in the physics community was developed by Aransonet al. [93], where the crack was identified as the normalized mass density of the material with0 in the crack and 1 in the bulk. The displacement equation is similar to Equation (2.24)except that the density is replaced by the normalized mass density. Besides, to interpolatethe stiffness in the crack interface, the stiffness was assumed to be linearly dependent on thenormalized density. The evolution of the normalized density is governed by the Ginzburg-Landau type evolution equation with a very simple energy functional only consisting of thegradient energy term and the double-well term. Utilizing the model, they then captured thebranching behaviour and the velocity oscillation phenomenon. However, the oversimplifiedenergy functional leads to the unphysical logarithmic dependence of the crack thickness onthe domain size [93].To address this problem, Karma et al. [2] proposed a phase-field model later on known asthe Karma-Kessler-Levine (KKL) model. Although the energy functional in the KKL modelis unique, the governing equations are the same as those in the mechanical models.Based on the well-defined sharp-interface equations for the interface migration in the phase25transformation problem, Brener and Spatschek [94] and Spatschek et al. [95] developed aphase-field fracture model. In the model, they assumed that the crack interface migrationcan be described by the same set of sharp-interface equations. Utilizing the model, they thenobserved the branching instability and showed that the critical energy for crack propagationagrees well with the Griffith theory. More importantly, instead of ignoring the crack tipshape like all the other phase-field fracture models, the crack tip curvature can be chosenfreely. However, the sharp-interface equations describing phase transformation may not beapplicable for crack propagation. Therefore, quantitative comparisons with experiments andtheories are required to verify the model.Almost all the phase-field fracture models are non-conversed. The only exception is themodel by Eastgate et al. [96] in which the phase-field parameter is defined as the normalizedmass density and governed by the diffusion equation. Unfortunately, the phase-field param-eter inside the crack cannot go to 0, and the value is proportional to the strain in the crack[97].Although many models showed promising results, most of them lack quantitative com-parisons with experiments or theories. The KKL model, as an exception, is the most well-developed model, and has reproduced various dynamic fracture phenomena, e.g., dynamiccrack branching [35, 98, 99], crack path oscillation [100], lance-shaped structure [101], etc.Therefore, only the KKL will be discussed in this context.2.6.2 Basics of the KKL ModelAs in other phase-field fracture models, the crack is identified with an order parameter φ inthe KKL model. With the phase-field parameter, the energy functional reads [98]:F =∫Ω(12ρ(∂u∂t)2+κ2(∇φ)2 + V (φ) + g(φ)(12ε : [Cε]− Ec))dΩ (2.27)where κ is the coefficient related to the surface energy, Ec is the critical elastic energy densityto break the material, and V (φ) is the double-well function given as [2]:V (φ) = mφ2(1− φ)2 (2.28)where m determines the height of the double-well barrier. As demonstrated by Karma andLobkovsky [35], in mode III fracture, when m equals 0, the critical branching velocity is theclosest to the theoretical prediction by Adda-Bedia [5]. However, in other fracture modes,the proper range of m was not given. The first term in Equation (2.27) is the kinetic energy,26and the rest of the terms are the interfacial and elastic energy. The elastic energy expression,i.e., g(φ)12ε : [Cε], is almost the same as those in the mechanical models, except that thedegradation function is different. According to Karma et al. [2], the order of the functionshould be greater than 3 so that the strain can be released through the crack. Therefore, thedegradation function was defined as [2]:g(φ) = 4φ3 − 3φ4 (2.29)The interfacial energy density, however, is different from those in the mechanical models andis given as [2]:finterfacial =κ2(∇φ)2 + V (φ) + (1− g(φ)) Ec (2.30)As in the models by Kuhn and Mu¨ller [36], and Miehe et al. [48], the phase-field parameteris governed by the Ginzburg-Landau evolution equation given as [2]:1M∂φ∂t= −δFδφ= κ∇2φ− 2mφ(1− φ)(1− 2φ)− 12(φ2 − φ3)(12ε : [Cε]− Ec)(2.31)The mobility (M), similar to that in the mechanical models, delays the energy dissipation,and therefore should be sufficiently large to obtain good agreement with the Griffith theory[35]. As in the mechanical models, to prevent the unphysical recovery of the crack, twokinds of irreversibility criteria were introduced, i.e., φ = 0 [35] and dφ/dt ≤ 0 [98]. Thecomparison between two irreversibility criteria, however, is missing in the literature. Basedon the conservation of the linear momentum, the displacement equation is given by:ρ∂2u∂t2= ∇ · ((4φ3 − 3φ4)Cε) (2.32)By solving Equation (2.31) with the displacement set to 0, the crack interface profile in1D (Figure 2.8) reads [102]:∂φ∂x=√2Ecκ(1− g(φ) + V (φ)Ec)(2.33)Combining Equations (2.30) and (2.33), the critical energy release rate reads [35]:Gc = 2√2κEc∫ 10√1− g(φ) + V (φ)Ec dφ (2.34)As in the mechanical models, the diffuse-interface approximation of the crack leads to thedeviation from the LEFM. Therefore, the interface thickness to domain size ratio has to besufficiently small. As demonstrated by Karma and Lobkovsky [35], the critical branchingvelocity showed good agreement with the fracture theory when the interface thickness todomain size ratio approaches 0.27Figure 2.8: The 1D analytical phase-field profile in the KKL model generated fromEquation (2.33).2.6.3 Applications of the KKL ModelDynamic Crack PropagationDistinct from the mechanical models, all the applications of the KKL model are limited to dy-namic brittle fracture with no application for quasi-static fracture. However, compared withthe mechanical models, the KKL model reproduced broader dynamic fracture phenomenaand enables deeper insight into the dynamic fracture mechanisms.Utilizing the basic KKL model, Karma and Lobkovsky [35] studied 2D branching phe-nomena in mode III fracture and found that the critical branching velocity is 0.41Cs and thebranching angle is 70° which are close to Adda-Bedia’s predictions, i.e., 0.39Cs and 79.2°.However, the model cannot distinguish between compression and tension because of thesimplified elastic energy expression, which limits its application to more general loading con-ditions. Henry and Levine [100], then decomposed the elastic energy to guarantee that thematerial cannot break under compression. Utilizing the ‘anisotropic’ model, they observedthe oscillated crack path under biaxial strain, which agrees well with the experiments.Henry [98] then studied the 2D branching phenomena in mode I fracture. Although thebranching velocity is in good agreement with the theoretical prediction, it is much higherthan the experimental observations. Henry then suggested that the discrepancy comes fromthe 2D simplification of 3D experiments. To demonstrate that point, they did a series of 3D28simulations. Starting with the 3D straight crack propagation under mode I loading, Henry[103] observed that the crack front in the thickness direction is elliptical rather than straight asassumed in 2D. Then Henry and Adda-Bedia [99] studied the 3D branching phenomena underthe same loading condition. They found that with a thin plate, the branch starts locally andspreads to the whole thickness of the plate, and eventually looks like the pattern observedin 2D simulations. However, the crack pattern becomes very complex with a thick plate,which agrees well with the experiments. Meanwhile, the crack velocity is lower in the thickplate because more elastic energy is consumed to create the complex crack pattern leavingless energy for the kinetic energy [99]. Therefore, this finding gives one plausible explanationfor the higher branching velocity predictions in 2D simulations compared with experiments.To understand the lance-shaped structure under mixed-mode (I and III) loading conditionobserved in experiments, Pons and Karma [101] managed to reproduce the phenomenon with3D simulations and provided insights into the possible mechanisms of the formation of thestructure.The anisotropic model with the direction-dependent fracture energy was proposed byHakim and Karma [104], in which the anisotropy was added through the gradient term(κ/2(∇φ)2) in the interfacial energy in Equation (2.30). Mesgarnejad et al. [105] thenapplied the model to study crack propagation in orientationally ordered composites, and thecrack path was in good agreement with experiments.Crack Propagation in Two-Phase MaterialsApart from the numerous applications in isotropic and anisotropic fracture, Henry [106] alsoapplied the model to study crack propagation in two-phase materials. In the study, onlythe heterogeneous critical energy release rate was considered and was by modifying Gc inEquation (2.34). With the model, Henry studied crack propagation in a two-phase system, inwhich a titled boundary separates the materials with different values of Gc (Figure 2.9a). Aninitial crack was then placed in the weak material. By applying the load less than the amountallowing crack propagation in the strong material but higher than that for crack propagationin the weak material, the crack was supposed to either propagate along the boundary orstop at the boundary. The critical energy for the continuous crack propagation along theboundary at various grain boundary angles was compared with the theory. With Gc in thestrong material much larger than that of the weak material, the results agreed qualitativelywell with the theory (Figure 2.9b). However, when Gc for the strong material becomes twice29that of the weak material, the unphysical penetration of the grain boundary was observed.Therefore, Henry suggested that the material heterogeneity cannot be incorporated simply bymodifying Gc as the interaction between the crack interface and the heterogeneous criticalenergy release rate leads to the unphysical crack penetration behaviour. However, moreinvestigations into the unphysical behaviour are required as the interaction between the crackinterface and the heterogeneity also exists in the mechanical models, while the mechanicalmodels showed good agreement with theories.(a) (b)Figure 2.9: (a) The simulation setup in which a boundary with an angle θ separatingtwo materials with I being the weak material with toughness ΓI and II being thetough material with toughness ΓII . (b) Comparison between the critical load(∆y/∆yc) for continuous crack propagation at the boundary from phase-fieldsimulations and that from the theoretical predictions with different boundaryangles. Green and red markers are from phase-field simulations with differentdomain sizes. Adapted from [106]. Reprinted with permission from Elsevier.30Chapter 3Scope and ObjectivesThe work aims to compare phase-field fracture models developed by the mechanics and thephysics community. In particular, two fracture models are examined: the phase-field fracturemodel by Miehe et al. [1] shortened as the mechanical model and the KKL model, whichis the most well-developed model in the physics community. Since the standard mechanicalmodel and the KKL model are for brittle fracture in isotropic materials, most simulationsare based on that setting. The overall goal can be achieved through a sequence of objectivesgiven below:• The mechanical model and the KKL model will be implemented in the quasi-staticsetting. For the KKL model, it will be the first time in the literature to have quasi-staticapplications. First, parameter analysis will be carried out to derive the proper rangesof the model parameters. Then, simulation results will be compared with the Griffiththeory to evaluate their abilities in capturing the quasi-static fracture phenomena.The fracture phenomena in heterogeneous materials will also be briefly covered byinvestigating crack propagation in a bi-grain system.• The dynamic fracture phenomena will be studied utilizing the mechanical model andthe KKL model. First, parameter analysis will be performed to determine the properranges of the parameters. Then, the theories by Eshelby [4] and Adda-Bedia [5] will beutilized to evaluate the accuracy of the models.31Chapter 4Methodology4.1 IntroductionTwo types of phase-field fracture models, i.e., the mechanical model and the KKL model, wereemployed in this study to assess their abilities in capturing quasi-static and dynamic fracturephenomena. In the present work, two scenarios were studied, i.e., the snap-back of a 1D elasticband (Figure 4.1a) and crack propagation in a 2D rectangular sheet (Figure 4.1b). Becausethis project does not intend to study crack initiation, in both scenarios, an initial crack withφ = 0 is manually placed in the domain. Since the computational cost in 2D is much higherthan that in 1D, the 1D simulations are to determine the ranges of the numerical parameters.Distinct from the 1D case where the loading direction is constrained to the x-axis, there aremultiple loading directions in 2D. Among all the fracture modes, the mode III fracture wasinvestigated in 2D as other loading conditions, e.g., mode I, II, and mixed-mode, wouldrequire the use of the ‘anisotropic’ models to distinguish between compression and tension.Besides, in both communities, there are different types of ‘anisotropic’ models, which wouldadd extra complexity to the analysis. Therefore, although the mode III fracture is of littlepractical relevance, it is still investigated in this project. Meanwhile, the energy functionalsand governing equations are briefly reviewed, since they are more specific rather than genericas those in Chapter 2. Then along with the initial and boundary conditions defined based onthe simulation setting, various numerical techniques were utilized to translate the models intocomputer programs. Since this work does not intend to compare results with experiments,the parameters used in the simulations are dimensionless.320 x∆−∆L(a)xyz∆−∆LxLylc(b)Figure 4.1: Simulation setups for (a) the snap-back of a 1D elastic band and (b) crackpropagation in a 2D rectangular sheet.4.2 Energy Functional and Governing Equations4.2.1 The Mechanical ModelAlthough numerous attempts have been made to replace the energy functional proposed byBourdin et al. [45], and some showed promising results, serious drawbacks of these modelsmake them less popular. Therefore, Bourdin et al.’s energy functional [45], as the mostpopular one in the mechanics community, was adopted in this project. In the 1D setting, thestiffness tensor (C) is reduced to Young’s modulus (E), and the 1D strain tensor (ε) can besimplified to ∇u. Similarly, in 2D mode III fracture, with the antiplane shear condition, C isreduced to shear modulus (µ). Meanwhile, the displacement field becomes a scalar field withnon-zero values only in the thickness direction (z-direction). Therefore, the energy functionalin Equation (2.12) is simplified to:F =∫Ω(Gc(l(∇φ)2 + (1− φ)24l)+(φ2 + η) 12C(∇u)2)dΩ (4.1)where C represents Young’s modulus (E) in 1D and the shear modulus (µ) in 2D. Instead ofderiving the phase-field and displacement field solution by alternative minimization method,the rate-independent phase-field equation and the displacement equations proposed by Mieheet al. [1] were implemented since they are much more popular in the mechanics community.33Later on, the mechanical model would refer to Miehe et al.’s model. Based on Equation(4.1), the governing equations read:2Gcl∇2φ+ Gc2l(1− φ)− Cφ (∇u)2 = 0 (4.2)∇ · σ = C∇ · ((φ2 + η)∇u) =0 quasi-staticρ∂2u∂t2dynamic(4.3)Distinct from the displacement equations, extra constraints are required when solving thephase-field equation such that the crack (φ = 0) cannot recover to the intact state (φ = 1).Unless explicitly mentioned, dφ/dt ≤ 0 was employed in the simulations.4.2.2 The KKL ModelAs for the KKL model, since the branching velocity in mode III fracture agrees well withAdda-Bedia’s theory [5] when neglecting the double-well function according to Karma andLobkovsky [35], the double-well term was removed in the simulations as well. Therefore, theenergy functional can be simplified to:F =∫Ω(κ2(∇φ)2 + (4φ3 − 3φ4)(12C(∇u)2 − Ec))dΩ (4.4)Compared with the rate-dependent phase-field equation used in the literature, the rate-independent equation has higher accuracy. Therefore, to be consistent with the mechanicalmodel, the rate-independent equation was adopted for the KKL model. By setting M toinfinity in Equation (2.31), the phase-field equation becomes:κ∇2φ− 12 (φ2 − φ3)(12C (∇u)2 − Ec)= 0 (4.5)As in the mechanical model, dφ/dt ≤ 0 was enforced to prevent the recovery of the crack. Bysubstituting the simplified energy functional (Equation (4.4)) into the mechanical equilibriumequation (∇ · σ = 0) and Equation (2.32), the displacement equations in quasi-static anddynamic fracture read:∇ · σ = C∇ · ((4φ3 − 3φ4)∇u) =0 quasi-staticρ∂2u∂t2dynamic(4.6)344.3 Numerical TechniquesSince the phase-field equations and the displacement equations are non-linear partial dif-ferential equations, they were discretized by the explicit finite difference method (FDM).With FDM, the computational domain was first discretized into linear grids in 1D or squareCartesian grids in 2D. Then, the space derivatives were approximated with the central finitedifference method. For example, the second derivative of φ at the grid point (i, j) in 2Dreads:∇2φ = ∂2φ∂x2+∂2φ∂y2≈ φi+1,j − 2φi,j + φi−1,jdx2+φi,j+1 − 2φi,j + φi,j−1dx2(4.7)where dx is the mesh size.After discretization, the phase-field and mechanical equilibrium equation were solvedby iteration methods. In the iteration methods, the discretized equation is reformulatedsuch that the grid-point value is calculated by its nearest neighbours. For example, in themechanical model, φ at the grid point (i) in 1D is given as:φi = ϕ(φi+1, φi−1) =2Gcl (φi+1 + φi−1) + dx2(Gc2l+ C (∇u)2)4Gcl + dx2(Gc2l+ C (∇u)2) , i = 1, 2, . . . n (4.8)where ϕ is the function of the nearest neighbours of φi. Then an initial guess (φ01, φ02, . . . φ0n)is made based on the initial condition. By substituting the initial guess into the righthand-side of Equation (4.8), the first approximation (φ11, φ12, . . . φ1n) after the first iteration can bederived. The iteration then continues until the largest relative difference between the currentand previous iteration, i.e., maxi=1,2,...n(|(φki − φk−1i )/φki |), becomes smaller than an imposedthreshold value.Both Gauss-Seidel iteration and Jacobi iteration are based on the same idea. However,the rule updating the grid-point value is different. With the Jacobi iteration, the updaterule is expressed as: φki = ϕ(φk−1i−1 , φk−1i+1 ), and the points are not updated until the end of oneiteration. With Gauss-Seidel iteration, the updated value enters into the calculation of thenext grid point immediately. For example, after calculating φ12 with ϕ(φ01, φ03), φ12 is used tocalculate the next grid point such that φ13 = ϕ(φ12, φ04) and so on.Since the current and previous iterations appear on different sides of the update rule in theJacobi iteration (φki = ϕ(φk−1i−1 , φk−1i+1 )), the parallel programming is accessible as the domaincan be divided into subdomains, and the grid points in those subdomains can be updatedparallelly. Although the Gauss-Seidel iteration cannot be parallelized, the convergence rateis much faster than that of the Jacobi iteration. After numerous tests, it was observed that35the computational cost is the lowest when solving the phase-field equation with the Gauss-Seidel iteration and the mechanical equilibrium equation with the paralleled code by Jacobiiteration. Therefore, the phase-field equations were solved by the Gauss-Seidel iteration, andthe mechanical equilibrium equations were solved by the Jacobi iteration and parallelizedwith the OpenMP library.The only time-dependent equations are the displacement equations in the dynamic casewhere the time derivative is approximated by the explicit Euler scheme given as:∂2u∂t2≈ ut+dt − 2ut + ut−dtdt2(4.9)where dt is the time step, ut−dt, ut and ut+dt are displacements at previous, current and nexttime step. To obtain the suitable ranges of numerical parameters, i.e., the threshold valuesin the iteration methods, the mesh size, the time step size, etc, a series of sensitivity analyseswere carried out in Chapter 5.4.4 Initial and Boundary Conditions4.4.1 Initial ConditionsAfter discretizing the domain into Cartesian grids, the initial and boundary condition shouldbe assigned based on the simulation settings. In 1D, an initial crack (φ = 0) is placed in themiddle of the band with two ends stretched as shown in Figure 4.1a. The crack interface isinitialized by the analytical phase-field profile (Equation (2.10) or (2.33)).For 2D simulations, the setting is shown in Figure 4.1b: an edge crack is positionedin the middle of a rectangular sheet with the top and bottom boundary stretched in thethickness direction (z-direction). Since the 1D analytical crack profile cannot be assigned toa 2D domain directly. The initial phase-field profile is determined by solving the phase-fieldequation with the displacement (u) set to 0. In this case, the initial guess for the phase-field parameter has to be set in the iteration methods. Because dφ/dt ≤ 0 was enforced asthe irreversibility criterion, the initial guess of φ cannot be chosen arbitrarily. Instead, itshould be narrower than the analytical profile such that φ can decrease to reach that profile.Figure 4.2 shows the initial guess of the phase-field parameter in the simulations. Apart fromthe crack where φ equals 0, φ at other points is set to 1. The crack interface then becomesa linear line with a slope equal to 1/dx. For the displacement field, the domain is elasticallystretched initially such that the strain is distributed uniformly.36φ = 0φ = 1initial guessanalytical profiledxFigure 4.2: Schematics of the initial guess (solid line) and the analytical solution (dot-ted line) for the φ profile at interface. The circular dots represent the meshes.4.4.2 Boundary ConditionsAfter setting the initial conditions, proper boundary conditions should be employed through-out the simulations. As for the displacement field, the boundary condition depends on theloading condition. In 1D simulations, since the two ends are stretched, the Dirichlet bound-ary condition was applied where the boundary displacements were fixed and set to ∆ on theright end and −∆ on the left end (Figure 4.1a). Neglecting the initial crack, the band iselastically stretched, such that the input elastic energy (Ein) in the system reads:Ein =12C(2∆L)2L =2C∆2L(4.10)where L is the length of the elastic band.Similarly, the Dirichlet boundary condition was applied at the top and bottom boundaryin 2D simulations where u equals ∆ at the bottom and −∆ at the top (Figure 4.1b). Theinput elastic energy is defined as the elastic energy along the width direction (y-direction)per unit length in the length direction (x-direction) since it has the same unit as the criticalenergy release rate. Therefore, the input elastic energy reads:Ein =12C(2∆Ly)2Ly =2C∆2Ly(4.11)where Ly is the width of the sheet. Since the left and right boundary are free from theexternal load, the boundary condition reads:σ · n = C (∇u) · n = C(∂u∂x,∂u∂y)· n = C∂u∂x= 0 (4.12)where n is the vector normal to the boundary. The gradient of displacement in the x-directionequals zero, which corresponds to the zero-flux boundary condition. In the traditional fracture37models, i.e., sharp interface models, the zero-flux boundary condition should also be imposedon the crack surface since it is free from the load. However, in both the mechanical modeland the KKL model, this condition is implicitly incorporated into the equation as inside thecrack, i.e., φ = 0, the stress is zero in both models (Equations (4.3) and (4.6)).Unlike the boundary conditions for the displacement field, which have clear physicalmeanings, the phase-field boundaries are more related to the numerical outcomes. In themechanics community, the zero-flux boundary condition was applied at all the phase-fieldboundaries. However, in the KKL model, the Dirichlet boundary condition (φ = 1) wasadopted at the top and bottom boundary to avoid the breakage of the loading boundaries.To compare these boundary conditions, 2D quasi-static simulations were carried out. In thesimulations, each type of boundary condition, i.e., Dirichlet condition or zero-flux condition,was applied at all phase-field boundaries at a time in both models. The parameters forthe simulations are listed in Table 4.1. To translate the models into computer programs,C language and the parallel programming library (OpenMP) were used. The codes can befound in Appendices A and B. With the output files from the program, the results werevisualized with Python utilizing the Matplotlib library.Table 4.1: Parameters in 2D quasi-static simulationsSymbol Parameter ValueLx Length of the 2D sheet 50Ly Width of the 2D sheet 50lc Length of the initial crack 15l Internal length scale in the mechanical model 1C Shear modulus 1Gc Critical energy release rate 0.5Ein Input elastic energy 1δ Crack interface thickness 4κ Coefficient before the gradient term in the KKL model 1.26Ec Critical elastic energy to break the material 0.05η Residual stiffness 10−10dx Mesh size 0.1∆φthreshold Threshold value for φ in the iteration methods 10−5∆uthreshold Threshold value for u in the iteration methods 10−738As shown in Figures 4.3a and 4.3c, the Dirichlet condition leads to the distortion of φat the left boundary for both models, and the non-zero gradient of φ at the rest of theboundaries for the mechanical model. The zero-flux condition, on the contrary, works wellon all boundaries for both models (Figures 4.3b and 4.3d). The non-zero gradient of φ atboundaries can affect the stress distribution and alter the crack propagation behaviour. Sincethe zero-flux boundary condition was used in the mechanical model in the literature, and itis identical to the Dirichlet condition at all the boundaries except the left one in the KKLmodel, it was adopted in both models in the following simulations.(a) (b)(c) (d)Figure 4.3: φ profiles with two boundary conditions in both models. (a) and (b) arefor Dirichlet condition and zero-flux condition in the mechanical model. (c) and(d) are for Dirichlet condition and zero-flux condition in the KKL model.39Chapter 5Benchmarking5.1 Sensitivity Analysis of the Threshold Values inthe Iteration MethodsTwo threshold values, i.e., the threshold values for the phase-field parameter (∆φthreshold) andthe displacement field (∆uthreshold), are introduced in the iteration methods. To determinetheir values, 1D quasi-static simulations were performed. In the simulations, C, which isYoung’s modulus in 1D, was set to 1, and the length of the elastic band (L) was set to 50.Three mesh sizes were investigated, i.e., 0.01, 0.05 and 0.1. The values of the rest of theparameters are given in Table 4.1.At first, simulations were carried out with the mesh size set to 0.1. To determine thethreshold values more efficiently, one threshold value was evaluated at a time by only solvingone of the governing equations. For ∆φthreshold, only the phase-field equation was solved,and the displacement was set to 0. By comparing the phase-field profiles with the analyticalsolutions (Equations (2.10) and (2.33)), ∆φthreshold can be determined. As shown in Figures5.1a and 5.1c, the phase-field profiles converge to the analytical solutions when ∆φthreshold isbelow 10−5 for both models. Similarly, by only solving the mechanical equilibrium equationwhile fixing the phase-field profile, the appropriate value of ∆uthreshold is determined. Asshown in Figures 5.1b and 5.1d, the stress curves converge when ∆uthreshold is below 10−7 forboth models.Then the threshold values were evaluated with different mesh sizes. As shown in Fig-ure 5.2, the threshold values decrease when reducing the mesh size. One possible explanationis that by decreasing the mesh size, the resolution increases. Hence, much higher accuracyis required to match all the grid-point values between the current and previous iteration.40(a) (b)(c) (d)Figure 5.1: Crack interface profiles and stress distributions in 1D simulations with dif-ferent φthreshold and uthreshold in both models. (a) and (b) are for the mechanicalmodel. (c) and (d) are for the KKL model. All the profiles are plotted based onthe right half of the 1D elastic band, so x = 0 indicates the middle of the band.41(a) (b)Figure 5.2: (a) Critical ∆φthreshold and (b) critical ∆uthreshold for the results to convergeunder different mesh sizes in the mechanical and the KKL model. Markersrepresent the data points.5.2 Sensitivity Analysis of the Mesh SizeThe mesh size does not only affect the threshold values in the iteration methods but also theaccuracy of the numerical solutions. Although a few grid points inside the crack interfaceare sufficient to have a smooth interface profile, much finer meshes are essential to accuratelycapture the discontinuous displacement inside the crack. To determine the mesh size, 1Dquasi-static simulations with various mesh sizes were carried out. According to Bourdin etal. [45], the non-zero mesh size results in the overestimation of the critical energy release rateby a factor of 1 +dx/4l in the mechanical model. Therefore, since l was 1 in the simulations,only the mesh sizes below 0.1 were investigated to ensure the accuracy of the critical energyrelease rate. Figure 5.3 compares the phase-field and stress profiles for different mesh sizesin both models. The analytical stresses for both models are 0 since all the tension is releasedafter the breakage of the elastic band. With the decrease of the mesh size, the phase-fieldand the stress profiles become closer to the analytical solutions. Even though the stresses donot become 0, they are below 0.07% of the initial value, which is negligible. The mesh sizeseems to have a small influence on 1D solutions. However, it might not be the case in 2D as42the stress ahead of the crack tip becomes higher.(a) (b)(c) (d)Figure 5.3: Crack interface profiles and stress distributions in 1D simulations withdifferent mesh sizes in both models. (a) and (b) are for the mechanical model.(c) and (d) are for the KKL model. All the profiles are plotted based on theright half of the 1D elastic band.Due to the high computational cost, only two mesh sizes, i.e., dx = 0.08 and dx = 0.1,were investigated in 2D quasi-static simulations. The values of the rest of the parameters43can be found in Table 4.1. Since there are two shear stresses, i.e., ∂u/∂x and ∂u/∂y, in 2D,to simplify the analysis, the elastic energy was analyzed instead. The phase-field and elasticenergy profiles along the dotted line in Figure 5.4 for two mesh sizes were then compared.Since the profiles along the dotted line will be studied extensively in the later chapters,they are referred to as the profiles in the crack plane. As shown in Figure 5.5, negligibledifferences were observed when varying the mesh size. However, the computational costbecomes significantly higher with dx equal to 0.08, which hinders investigations into themesh size below that value. Therefore, dx was set to 0.1 in the quasi-static simulations.xyFigure 5.4: The schematic representation of the simulation domain in 2D. The dottedline lies in the crack plane along the x-direction.44(a) (b)(c) (d)Figure 5.5: φ profiles and elastic energy profiles in the crack plane with different meshsizes in both models. (a) and (b) are for the mechanical model. (c) and (d) arefor the KKL model.5.3 Sensitivity Analysis of the Time Step SizeSimilar to the non-zero mesh size, the non-zero time step size (dt) in the displacement equa-tion also leads to the discretization error. Therefore, 2D dynamic simulations with varioustime step sizes were performed. Apart from the parameters used in the 2D quasi-static sim-45ulations (Table 4.1), the density of the material (ρ) was set to 1. As shown in Figure 5.6, inboth models, the crack profiles converge when dt is below 0.06, and the elastic energy profilesin the crack plane overlap when dt is below 0.02. Therefore, dt was set to 0.02 in dynamicsimulations.(a) (b)(c) (d)Figure 5.6: φ profiles and elastic energy profiles in the crack plane with different timestep sizes in both models. (a) and (b) are for the mechanical model. (c) and (d)are for the KKL model.465.4 Sensitivity Analysis of the Crack Front PositionDistinct from the sharp interface models where the interface position is definite, in phase-field models, the crack interface distributes over a narrow band of which the width is theinterface thickness (δ). As a result, the crack front position varies with the value of φ, which isreferred to as φfront in the current work. The uncertainty of the crack front position leads todifficulty in capturing the energy release rate, crack velocity, and branching angle. Therefore,2D quasi-static and dynamic simulations were carried out to determine the value of φfront.Since φ = 0 represents the crack, φfront should be close to 0 to represent the position of thecrack tip. Therefore, only the φfront below 0.5 was analyzed.5.4.1 φfront in the Quasi-static SettingAccording to Equation (2.1), the energy release rate can be calculated from the elastic energyof the system and the crack propagation length. Similarly, the surface energy of the systemand the propagation length are required to determine the critical energy release rate. Sinceφfront affects the propagation length, it could potentially influence these two parameters.Therefore, 2D quasi-static simulations were performed to evaluate the influence of φfront onthe energy release rate and the critical energy release rate. The parameters are given inTable 4.1.Figure 5.7 compares the energy of the system with respect to the crack propagation lengthwith different φfront. In almost all the cases except the one with φfront = 0.1 in the KKLmodel, the slopes of the elastic energy and surface energy curves change abruptly in thebeginning.For the mechanical model, the abrupt transition can be attributed to two factors. Thefirst is that the initial phase-field profile, which is the analytical solution, is unstable. Asshown in Figure 5.8a, the crack interface profile quickly stabilizes and deviates from theanalytical solution after the first computational loop. Meanwhile, the surface energy andthe elastic energy transit to the stable regime where the slopes of the curves stabilize. Theother reason is that the position of φfront = 0.1 lags behind that of φfront = 0.5 as shownin Figure 5.8b, which also leads to the larger surface energy increment and elastic energyrelease in the case with φfront = 0.1 (Figure 5.7a). Especially in the first few loops, after thesame amount of elastic energy release and surface energy increment, with φfront set to 0.1,the almost zero propagation length leads to larger energy gradients at the beginning.As in the mechanical model, in the KKL model, with φfront set to 0.1, the abrupt transition47(a)(b)Figure 5.7: Elastic energy release and surface energy increment with respect to thecrack propagation length for (a) the mechanical model and (b) the KKL modelwith different φfront.48(a) (b)Figure 5.8: φ profiles (a) at the right half of the interface and (b) in the crack planein the first few computational loops in the mechanical model. The blue and thered dotted line in (b) correspond to the crack front position with φfront set to0.5 and 0.1.in the energy curves occurs because the position of φfront = 0.1 lags behind that of φfront = 0.5(Figure 5.9b). For the case with φfront = 0.5, the energy curves are linear (Figure 5.7b)because the crack interface profile does not change abruptly in the first few computationalloops as shown in Figure 5.9a.In both models, the initially unstable regime is sensitive to φfront. As a result, at thesame crack propagation length, the energy varies with φfront. However, in both models, theslopes of the energy curves in the stable regime do not change with φfront (Figure 5.7). Sincethe slopes of the elastic energy and surface energy curve represent the energy release rateand the critical energy release rate, they are much more important than the total amount ofenergy in the system. By neglecting the initially unstable regime and only focusing on theslopes at the stable regime, the results are insensitive to φfront. Therefore, φfront was set to0.5 in quasi-static simulations as it results in linear energy curves in the KKL model.49(a) (b)Figure 5.9: φ profiles (a) at the right half of the interface and (b) in the crack planein the first few computational loops in the KKL model. The blue and the reddotted line in (b) correspond to the crack front position with φfront set to 0.5and 0.1.5.4.2 φfront in the Dynamic SettingAs in the quasi-static case, φfront can also affect the dynamic fracture phenomena. In theliterature, some people argued that the crack velocity and branching angle cannot be deter-mined accurately because of the diffuse interface in the phase-field models. However, φfrontwas set to 0.5 by Karma and Lobkovsky [35] to calculate the velocity and branching angle.Therefore, to examine the influence of φfront on the crack velocity and branching angle, 2Ddynamic simulations were carried out. The geometry parameters are listed in Table 5.1. Thedensity (ρ) was set to 1, the time step size (dt) was set to 0.02, and the rest of the materialparameters can be found in Table 4.1.Figures 5.10a and 5.10c compare the crack velocity with respect to the crack tip positionwith different φfront. The unit of the crack velocity is the shear wave speed (Cs), whichequals√µ/ρ. All the curves predict the same trend: the crack velocity increases thendecreases during blunting and becomes unstable after branching. The velocity oscillationafter branching is commonly observed in experiments, which is related to the interaction ofthe microcracks at the crack tip. With the decrease of φfront, the velocity becomes lower in50Table 5.1: Geometry parameters in 2D dynamic simulationsSymbol Parameter ValueLx Length of the 2D sheet 120Ly Width of the 2D sheet 120lc Length of the initial crack 50both models, which is because the slope of the line at the crack front (ab) decreases duringcrack propagation as shown in Figures 5.10b and 5.10d. Therefore, the rotation of the lineab adds the additional angular velocity. Meanwhile, the velocity curves shift leftwards whendecreasing φfront as the position of the smaller φfront lags behind that of the larger φfront.As for the branching velocity, which is indicated by the dotted lines in Figures 5.10a and5.10c, when varying φfront from 0.5 to 0.1, it decreases by 20% in the mechanical modeland 9% in the KKL model. Compared with the KKL model, the branching velocity in themechanical model is sensitive to φfront. Therefore, when analyzing the branching velocity inthe mechanical model, the 20% error induced by φfront should be taken into consideration.Another noticeable feature in the velocity curves is that the velocity increases abruptlyimmediately after branching, which is because of the way the velocity is calculated. Here,the crack tip position is defined as the point with the maximum x-axis value in the contourline of φ = φfront. With this definition, the critical time for crack bifurcation is when thenumber of points with the maximum x-axis value changes from one to two. However, withthis method, immediately after branching, significant errors occur in the velocity calculation.As shown in Figure 5.11, at the moment of branching, the crack tip position changes from ato b or c. As a result, the sudden increase in the y coordinate of the crack tip position leadsto a major increase in the velocity, which is a numerical artifact.The diffuse interface also leads to the problem of defining the branching angle. Figure 5.12compares the contour lines of φ = φfront. Since the branches are not thin lines, the actualbranching angle cannot be extracted from the contour directly. The most straightforwardway is to select lines inside the branches and view the angle between the selected lines as thebranching angle. Here, the line is selected visually such that the branch is almost symmetricwith it. With this method, the branching angle increases 8° in the mechanical model anddecreases 2° in the KKL model when varying φfront from 0.1 to 0.5. Considering the errorintroduced when selecting the lines, the branching angle in both models is insensitive to51(a) (b)(c) (d)Figure 5.10: Crack velocity with respect to the propagation length with different φfrontfor (a) the mechanical model and (c) the KKL model. Phase-field profiles inthe crack plane at different times for (b) the mechanical model and (d) theKKL model. The dotted lines in (a) and (c) indicate the branching velocities.φfront. To be consistent with the work by Karma and Lobkovsky [35], φfront was set to 0.5in dynamic simulations.52Figure 5.11: The contour line of φ = 0.5 in the KKL model. a, b and c are the pointsindicating the crack tip.(a) (b)Figure 5.12: Contour lines of φ = φfront for (a) the mechanical model and (b) the KKLmodel.535.5 Sensitivity Analysis of the GeometricConfiguration in the Quasi-static SettingThe geometric configuration refers to the crack length, domain width, and domain length.In dynamic fracture, the radiated elastic wave from the crack tip can be reflected at theboundaries and interact with the crack. Therefore, to avoid the interaction between thecrack and the boundary, the crack length has to be very long. In the quasi-static case,however, this situation does not exist with the mechanical equilibrium condition. Therefore,the crack length was not analyzed in the quasi-static setting.The domain size, however, does influence the accuracy of the simulations. As mentionedin Chapter 2, in both the mechanical model and the KKL model, the interface thickness todomain size ratio has to be very small to obtain good agreement with either the Griffith theoryor the experimental results. Since varying the interface thickness requires the adjustment ofthe mesh size, only the domain size was analyzed. Therefore, to determine the domain size,2D quasi-static simulations were performed. Square domains were used in the simulations,and two domain lengths, i.e., 50 and 100, were investigated. In case the initial crack lengthinfluences the results, the initial crack length to domain size ratio was kept constant whenvarying the domain size. The rest of the parameters are listed in Table 4.1.Figure 5.13 compares the elastic energy and surface energy with respect to the propagationlength with different domain sizes. In both models, with the increase of the domain size, theelastic energy and the surface energy curves shift upwards. For the elastic energy, that iscaused by the difference in the initial elastic energy in the system. Although Ein is thesame with different domain sizes, according to Equation (4.11), it only accounts for the inputenergy in the y-direction (width direction). Therefore, considering that the length of thedomain is different, the initial elastic energy in the system is higher with a larger domainsize. The difference in the initial surface energy is because the initial crack length also varieswith the domain size. However, neither the elastic energy nor the surface energy affects crackpropagation.On the contrary, the energy release rate and critical energy release rate, which can becalculated from the slopes of the elastic energy and surface energy curves in the stable regime,determine the crack propagation behaviour. In the mechanical model, when reducing thedomain size by half, the energy release rate and critical energy release rate decrease by 7.7%and 9.6% (Figure 5.13a). While in the KKL model, these two parameters do not vary withthe domain size (Figure 5.13b). In the mechanical model, ideally, larger domain sizes should54(a)(b)Figure 5.13: Elastic energy release and surface energy increment with respect to thecrack propagation length for (a) the mechanical model and (b) the KKL modelwith different domain sizes.55be investigated until the results converge. However, since the computational cost becomessignificantly higher even by setting the domain size to 100, the smaller domain size, i.e., 50,was adopted in the following quasi-static simulations.5.6 Sensitivity Analysis of the GeometricConfiguration in the Dynamic SettingIn dynamic fracture, the geometric configuration is even more important than that in quasi-static fracture as the elastic wave can interact with the crack and lead to distinct fracturephenomena. In fracture theories, since it is impossible to consider all the scenarios, moregeneralized and simplified settings are usually adopted. In these settings, some assumptions,e.g., infinite large domain and semi-infinite crack, seem impractical in experiments. Undercertain conditions, however, the theoretical predictions can be treated as the asymptoticsolutions of the experimental results.There are two settings commonly adopted in fracture theories: a crack in an infinitedomain where the boundaries do not affect the crack tip fields, and a semi-infinite crack inan infinitely long strip where only the boundaries in the width direction affect the crack fields.Adda-Bedia [5] studied the dynamic fracture phenomena in 2D considering a semi-infinitecrack in an infinitely long strip under antiplane shear condition. Therefore, 2D dynamicsimulations were carried out to make the simulation setting equivalent to that in Adda-Bedia’s theory. In the simulations, the density (ρ) was set to 1, the time step size (dt) wasset to 0.02, and the rest of the material parameters can be found in Table 4.1.5.6.1 Initial Crack LengthSince having a semi-infinite crack is impossible in the simulations, analysis has to be madebefore selecting the initial crack length. In the simulations, different crack lengths wereanalyzed. The domain width and length were set to 120. Figure 5.14 compares the crackvelocity with different initial crack lengths in both models. All the curves converge beforethe propagation length is 20. However, the curves diverge in the later stage as the cracktip fields start interacting with the right boundary. Since the velocities are the same withdifferent lc in the early stage, lc was set to 50.56(a) (b)Figure 5.14: Crack velocity with respect to the propagation length with different initialcrack lengths for (a) the mechanical model and (b) the KKL model.5.6.2 Domain LengthIn Adda-Bedia’s theory, the domain length was assumed to be infinitely long, which is im-practical in the phase-field simulations. Therefore, to determine the proper range of thedomain length, 2D dynamic simulations were carried out. In the simulations, the domainwidth and the initial crack length were set to 120 and 50. Figure 5.15 compares the crackvelocity with different domain lengths in both models. In all the cases, the spike in thevelocity curves is the numerical artifact caused by the method calculating the velocity, whichhas been demonstrated in Section 5.4. When the domain length is above 120, the velocitycurves converge. Therefore, Lx was set to 160.57(a) (b)Figure 5.15: Crack velocity with respect to the propagation length with different do-main lengths for (a) the mechanical model and (b) the KKL model.5.6.3 Domain WidthThe last geometry parameter to be determined is the domain width. Distinct from the othertwo parameters, the domain width is unspecified in the fracture theory. However, both themechanical model and the KKL model require the use of a large domain width to increasethe accuracy. Therefore, 2D dynamic simulations were carried out to determine the domainwidth. According to Knauss [107], in experiments and simulations, the relative size of theinitial crack length with respect to the domain width has to be large to approximate the semi-infinite crack in theories. Therefore, to avoid varying the relative size of the crack lengthwhen increasing the domain width, the ratio of the initial crack length to domain width todomain length was kept constant. The geometry parameters for the simulations are listed inTable 5.2.Figure 5.16 compares the crack velocity with different values of the domain width. Withthe increase of the domain width, the crack velocity decreases in both models, which is con-sistent with the observations by Henry [98]. When varying Ly from 120 to 160, in maximum,the velocity decreases by 5.7% in the mechanical model and 5.8% in the KKL model. Con-sidering that the velocity could vary by up to 20% with different φfront as demonstrated inSection 5.4, the change in the crack velocity with different domain widths is negligible.58Table 5.2: Geometry parameters for three casesCase Domain length (Lx) Domain width (Ly) Initial crack length (lc)1 160 120 502 186 140 58.33 213 160 66.7(a) (b)Figure 5.16: Crack velocity with respect to the propagation lengths with different do-main widths for (a) the mechanical model and (b) the KKL model.Figure 5.17 compares the branching angle with different domain widths in both models.When increasing the domain width from 140 to 160, the branching angle only increases 5°in the mechanical model. Considering the error induced from selecting the dotted lines inFigure 5.17a, the increase in the branching angle is negligible. In the KKL model, the contourlines converge when the domain width is above 140 as shown in Figure 5.17b. Therefore, tohave an accurate representation of the crack shape with relatively lower computational cost,Ly was set to 140, and the geometry parameters of case 2 in Table 5.2 were used in dynamicsimulations.59(a) (b)Figure 5.17: Contour lines of φ = 0.5 with different domain widths for (a) the mechan-ical model and (b) the KKL model. The ticks of the x-axis and y-axis are notgiven as they are different for three cases in order to match the crack shapes.The dotted lines in (a) indicate the branching angles.60Chapter 6Simulations of Quasi-static CrackPropagation6.1 IntroductionThis chapter focuses on the quasi-static applications of the mechanical model and the KKLmodel. Utilizing the computer programs in Appendices A and B, general comparisons be-tween the models were first carried out. Then to assess the robustness of the models, someof the material parameters were evaluated. Since crack propagation in heterogeneous ma-terials is a challenging topic, these models were also extended to study crack propagationin a bi-grain system to clarify their abilities in describing the role of the grain boundaryduring fracture. All the simulations were in 2D, and the values of the parameters used in thesimulations can be found in Table 4.1 unless explicitly stated.6.2 General Comparisons between the Models6.2.1 Crack ProfileAs shown in Figure 6.1a, the major discrepancy in the interface profiles between both modelsis that φ in the bulk is less than 1 in the mechanical model. This is one of the side effectsof using the quadratic degradation function (Equation (2.9)), as φ can decrease from 1 withany non-zero strain according to Equation (4.2). The phenomenon was referred to as damageaccumulation in the literature, which has been discussed in Chapter 2. With higher stress atthe crack front region, the effect of damage accumulation becomes more severe. As a result,in the mechanical model, φ in the crack front is even lower than that in the bulk as indicatedby the red dotted line in Figure 6.1b.From Figure 6.1b, it is also noticeable that inside the crack, φ is non-zero in the KKLmodel, because of the degradation function (Equation 2.29). The non-zero φ also leads to61(a) (b)Figure 6.1: (a) Right half of the crack interface profiles and (b) crack profiles in thecrack plane for both the mechanical model and the KKL model.a significant amount of residual elastic energy inside the crack since the elastic energy issensitive to φ. As shown in Figure 6.2, the residual elastic energy in the KKL model isaround 30% of the corresponding peak elastic energy. Although there is also residual elasticenergy in the mechanical model, but because of the close-to-zero φ inside the crack, i.e.,around 0.002, the magnitude is negligible. In the literature, some people in the mechanicscommunity suggested setting φ to 0 when it is below an imposed threshold value. However,since that approach is not well accepted and may also affect other model outcomes, it was notincluded in the present work. In conclusion, the KKL model is more accurate in describingthe phase-field profiles since it does not have the damage accumulation issue. However,the non-zero φ in the KKL model leads to a large residual elastic energy inside the newlydeveloped crack.62Figure 6.2: Elastic energy profiles in the crack plane for both the mechanical modeland the KKL model. The blue dotted lines indicate the residual elastic energyinside the crack.6.2.2 Stress and Elastic Energy Profiles in the Crack PlaneThe crack profile is a numerical feature of the model which may not influence crack propaga-tion. However, it can affect physical features, e.g., stress and elastic energy distribution. Inthis section, the stress and elastic energy in both models are compared with the analyticalsolutions.In 2D mode III fracture, considering a crack in an infinitely large domain, the analyticalstress ahead of the crack tip in the crack plane reads [108]:(σzxσzy)= 0σ∞zy√1−(a/x)2 (6.1)where σzx and σzy are the shear stresses in x direction (propagation direction) and y direction(normal to the propagation direction), a is the crack length, x is the x-coordinate of the stress,and σ∞zy is the stress far away from the crack tip. The analytical stress in other directionscan be found in [108]. Since the general stress solution is much more complicated, only thestress distribution in the crack plane was evaluated. In the simulations, σ∞zy is 0.173, and the63crack length (a) was determined from the φ profiles.From simulations, it was observed that σzx in both models is within a margin of theorder of 10−5, which agrees well with the analytical solution, i.e., zero. For σzy, however,large deviations can be observed in both models as shown in Figure 6.3. Close to the cracktip, because of the diffuse interface, both models cannot capture the stress singularity aspredicted by the analytical solution. However, compared with the mechanical model, theKKL model agrees better with the analytical solution near the crack tip. In the KKL model,the small deviation from the analytical solution ahead of the peak stress is because of theboundary effect. Far from the crack tip, the KKL model is almost overlapping with theanalytical solution. For the mechanical model, as in the near-tip stress field, the stress islower than the analytical solution. The lower stress is because φ is less than 1 ahead of thecrack tip, and the stiffness is a function of φ in the phase-field model, i.e., g(φ)µ. The loss ofstiffness in the mechanical model, as a result, underestimates the stress ahead of the cracktip. Therefore, the KKL model is more accurate in describing the stress distribution far fromthe crack tip.Figure 6.3: The shear stress in the y-direction (σzy) in the crack plane for both themechanical model and the KKL model.The near-tip stress can be captured accurately by neither of the models. The elastic64energy profile in both models, on the contrary, exhibits much better agreement with theanalytical solution. With the analytical stress in Equation (6.1), the analytical elastic energyahead of the crack tip in the crack plane reads:Eelastic =12µ(σ2zx + σ2zy)=(σ∞zy)22µ(1− (a/x)2) (6.2)As shown in Figure 6.4, in both models, compared with the stress distribution, the peakelastic energy is much higher than the far-field elastic energy. However, contrary to the trendobserved in the stress profiles, the peak elastic energy in the mechanical model is higher thanthat in the KKL model. It is because the stress and elastic energy at the interface have thefollowing relation:Eelastic =12µg(φ)(σ2zx + σ2zy)(6.3)Because of the dependence of the stiffness (µ) on the degradation function (g(φ)), the elasticenergy at the interface also relies on the φ profile. As a result, in the mechanical model, eventhough the peak stress is lower, the peak elastic energy is higher compared with that in theKKL model. Far from the crack tip, the elastic energy in the KKL model still agrees betterwith the analytical solution due to the damage accumulation issue in the mechanical model.In conclusion, both models are not accurate in describing the near-tip stress and elasticenergy distribution because of the diffuse interface. As for the stress and elastic energy farfrom the crack tip, the KKL model is more consistent with the analytical solution. Since theelastic energy contains all the stresses, i.e., σzx and σzy, and can also provide the informationon the energy state of the system, the elastic energy is analyzed in the following sections.65Figure 6.4: Comparisons of the elastic energy profiles in the crack plane with the an-alytical solution.6.2.3 Energy Release Rate and Critical Energy Release RateAlthough the stress, elastic energy, and crack profile provide some information about theaccuracy of the models, the only way to quantify the accuracy is by examining the energyrelease rate and the critical energy release rate. Figure 6.5a compares the elastic energy withrespect to the propagation length for both models. The slope of the curve indicates the energyrelease rate. Theoretically, the energy release rate always equals the critical energy releaserate. In the mechanical model, the energy release rate (GM) is 1.2Gc, while GK is 1.5Gcin the KKL model. Compared with the analytical value, i.e., Gc, both models overestimatethe energy release rate, and the mechanical model is more accurate. For the critical energyrelease rate, it can be calculated from the slope of the surface energy curve in Figure 6.5b.In both models, the critical energy release rate is slightly overestimated but the mechanicalmodel is more accurate.In both models, the energy release rate is larger than the critical energy release rate,which suggests that some elastic energy is lost not because of crack propagation. However,it is not surprising as the energy minimization equation (phase-field equation) can inevitablylead to the decrease of the total energy of the system, even though that results in deviations66from the LEFM. However, interestingly, in the KKL model, the energy release rate happensto equal the input elastic energy, which is also 1.5Gc. It suggests that the energy releaserate depends on the applied load, i.e., input elastic energy, instead of the material property,i.e., critical energy release rate. However, it might also simply be a coincidence and causedby some inappropriate choices of the material parameters. To examine the possible reasons,more detailed analyses are given in the latter part of this chapter.In conclusion, both the mechanical model and the KKL model can accurately describethe critical energy release rate. However, both models overestimate the energy release rate,and whether the energy release rate in the KKL model depends on the input elastic energyrequires further investigations.(a) (b)Figure 6.5: (a) Elastic energy release and (b) surface energy increment with respectto the crack propagation length for the mechanical model and the KKL model.GM and GK are the energy release rates for the mechanical model and the KKLmodel. GMc and GKc are the critical energy release ratees for the mechanicalmodel and the KKL model.6.3 Influence of the Irreversibility CriterionAs mentioned in Chapter 2, the irreversibility criterion should be included in both the me-chanical model and the KKL model to avoid the crack recovery, i.e., φ increases from 0 to671. Two types of irreversibility criteria were introduced in the literature, i.e., force φ = 0 atthe crack and dφ/dt ≤ 0. Although some 1D comparisons of these two criteria were carriedout in the literature, they focused on crack initiation while neglecting the influence on crackpropagation. Therefore, in this section, these irreversibility criteria were studied with 2Dcrack propagation simulations.Figure 6.6 compares the phase-field and the elastic energy profiles in the crack planefor two irreversibility criteria in both models. When the crack tips reach the same position(Figures 6.6a and 6.6c), the differences in the peak elastic energies are negligible as shownin Figures 6.6b and 6.6d. In atomic scale, the peak elastic energy is directly related to theatomic bound ruptures. While in the phase-field models, it controls the reduction of φ atthe crack front. Therefore, the peak elastic energy determines the crack motion in the nextmoment. Since the peak elastic energies are close, the tendency of crack propagation isinsensitive to the irreversibility criterion. However, before the peak elastic energy, the elasticenergy profiles are less smooth for φ = 0 than those with dφ/dt ≤ 0 in both models. Thepossible reason is that inside the newly developed crack, φ is slightly larger than 0, whichmakes the criterion φ = 0 unable to constrain the evolution of φ inside the crack. Since theelastic energy profile is quite sensitive to φ, the unconstrained φ inside the crack leads to thedistortion in the elastic energy profiles.Distinct from the phase-field profiles in the crack plane, in both models, the crack interfaceprofiles are quite different as shown in Figure 6.7. After crack propagation when little elasticenergy is left, the interface profiles recover to the analytical solutions with φ = 0, while theprofiles deviate from the analytical solutions with dφ/dt ≤ 0. However, even though theinterface profile is more accurate with φ = 0, it does not necessarily mean that the criticalenergy release rate is also more accurate.Therefore, to quantify the effect of the irreversibility criterion on crack propagation, theenergy release rate and critical energy release rate were evaluated. Figure 6.8 compares theenergy with respect to the crack propagation length with different irreversibility criteria.The energy release rate and the critical energy release rate were obtained from the slopesof the elastic energy and surface energy curves in the stable regime. In both models, theenergy release rate is insensitive to the irreversibility criterion. The critical energy releaserate, however, varies with the irreversibility criterion in both models and is more sensitiveto it in the mechanical model. It is not surprising as the interface profile in the mechanicalmodel is more sensitive to the irreversibility criterion. In the mechanical model, the criticalenergy release rate is 1.03Gc with dφ/dt ≤ 0 and 0.81Gc with φ = 0. Meanwhile, in the KKL68(a) (b)(c) (d)Figure 6.6: φ and elastic energy profiles in the crack plane for both models with dif-ferent irreversibility criteria. (a) and (b) are for the mechanical model. (c) and(d) are for the KKL model.model, it is 1.06Gc with dφ/dt ≤ 0 and 0.98Gc with φ = 0. In both models, the cases withφ = 0 underestimate the critical energy release rate, while the cases with the other criterionoverestimate the critical energy release rate. Surprisingly, in the mechanical model, with φ =0, the critical energy release rate is underestimated for around 20%, even though the interface69(a) (b)Figure 6.7: Right half of the crack interface profiles for (a) the mechanical model and(b) the KKL model with different irreversibility criteria.profile is quite accurate compared with the analytical solution (Figure 6.7a). However, it isstill reasonable as the 1D interface profile cannot accurately represent the critical energyrelease rate, which is an outcome of the surface energy change in 2D. Compared with themechanical model, the change of critical energy release rate with different irreversibilitycriteria is negligible in the KKL model.In conclusion, in both models, with dφ/dt ≤ 0, the elastic energy profile is smooth, andthe critical energy release rate is more accurate with a maximum 6% overestimation comparedwith Gc. Therefore, dφ/dt ≤ 0 is a better option at least in this specific setting.70(a)(b)Figure 6.8: Elastic energy and surface energy with respect to the propagation lengthfor (a) the mechanical model and (b) the KKL model with two irreversibilitycriteria.716.4 Influence of the Residual StiffnessTo prevent numerical instability, residual stiffness was introduced in the mechanical model.The value of it should be small, e.g., less than 10−6, to avoid the overestimation of theelastic energy. In the KKL model, however, residual stiffness has not been introduced in theliterature. Therefore, to understand whether residual stiffness is required in the KKL model,2D simulations were carried out in this section.First, the initial crack thickness, where φ = 0, was set to 1dx. In this situation, in bothmodels, simulations were stable even without the residual stiffness. However, with the crackthickness larger than 2dx, the simulations became unstable in both models unless the residualstiffness is non-zero. Therefore, the residual stiffness is required even in the KKL model whenthe crack thickness is larger than 2dx. In the present work, the residual stiffness was addedto the KKL model in the same way as the mechanical model. By adding η to Equation (4.4),the energy functional in the KKL model reads:F =∫Ω(κ2(∇φ)2 + (4φ3 − 3φ4 + η)(12C(∇u)2 − Ec))dΩ (6.4)To understand the role of the residual stiffness, the crack thickness in the simulationswas set to 10dx. Since the time scale does not exist in the quasi-static case, the number ofloops solving the phase-field equation and the mechanical equilibrium equation was treatedas a numerical time scale. Figure 6.9 compares the phase-field profiles and the elastic energyprofiles in the crack plane after the same number of loops. Almost no difference can beobserved in the phase-field profiles when varying η in both models. However, as expected,the elastic energy inside the crack is overestimated with a large residual stiffness in bothmodels. With the residual stiffness less than 10−10, the elastic energy profiles converge inboth models.Meanwhile, when solving the mechanical equilibrium equation with the Jacobi iterationmethod, the average numbers of iterations in each loop become close when η ≤ 10−10 asshown in Figure 6.10. Since the mechanical equilibrium equation is much more costly to solvethan the phase-field equation, the number of iterations solving the mechanical equilibriumequation reflects the computational time of the simulation. Therefore, in both models, withthe decrease of the residual stiffness, the computational time will eventually converge. It isalso noticeable that the computational cost of the mechanical model is significantly lowerthan that of the KKL model, which is possibly due to the use of the different degradationfunctions (Equations (2.13) and (2.29)).72(a) (b)(c) (d)Figure 6.9: φ and elastic energy profiles in the crack plane along the propagation lengthfor both models with different residual stiffness. (a) and (b) are for the mechan-ical model. (c) and (d) are for the KKL model.73Figure 6.10: Average number of iterations solving the mechanical equilibrium equationwith respect to the residual stiffness in the mechanical model and KKl model.The markers represent the data points.6.5 Influence of the Crack Tip CurvatureAlthough the crack is assumed to be a sharp notch in the LEFM, it can have different shapesin reality, which can affect crack propagation. For example, the crack tip curvature (κcur),which equals the inverse of the crack tip radius (rtip) (Figure 6.11), can alter the expressionof the crack tip fields, and therefore changes the stress intensity factor. However, in theliterature, the curvature effect was investigated by neither the mechanical model nor theKKL model. To understand the abilities of these models in capturing the curvature effect,2D simulations with various crack tip curvatures were carried out in this section. In thesimulations, the shape of the crack tip, where φ = 0, was described by parabolas. The crackthickness was set to 10dx.Figure 6.12 compares the phase-field and elastic energy profiles in the crack plane withvarious crack tip curvatures. In the beginning, the elastic energy at the crack tip becomeshigher when increasing the crack tip curvature. Then, shortly after the crack starts propa-gating, the energy profiles overlap. Therefore, the crack tip curvature only affects the initial74xyrtipFigure 6.11: The schematics of the crack tip shape. The crack tip radius (rtip) is theradius of the circle that can approximate the curve of the crack tip.short period of crack propagation, and its influence on the overall process is negligible.This finding can be further demonstrated by the contours of the phase-field profiles.Figure 6.13 compares the contour lines of φ = 0 with different curvatures. Shortly after thecrack starts propagating, the crack becomes a sharp notch in both models. Therefore, thecrack shape effect can be captured by neither of the models.(a) (b)Figure 6.12: Elastic energy profiles in the crack plane in the beginning and shortly afterpropagation with different κcur for (a) the mechanical and (b) KKL model.75(a) (b)Figure 6.13: Contour lines of φ = 0 with different κcur in (a) the mechanical modeland (b) the KKL model.6.6 Influence of the Input Elastic EnergyIn fracture theories, the energy release rate and critical energy release rate should be inde-pendent of the load in quasi-static brittle fracture. Therefore, ideally, both models shouldbe insensitive to the input elastic energy. However, as shown in Figures 6.14a and 6.14c, theenergy release rate in both models are sensitive to the input elastic energy and become closerto the analytical value, i.e., Gc, with the decrease of Ein. In the KKL model, the energyrelease rate equals the input elastic energy in both cases. Considering all the analyses inthe previous sections, it is safe to say that all the parameters are chosen properly, and thedependence of the energy release rate on Ein is not a coincidence. This unphysical behaviourcan also be related to the sensitiveness of the critical energy release rate with the input elasticenergy in the KKL model (Figure 6.14d). For the mechanical model, on the contrary, thecritical energy release rate is relatively insensitive to the input elastic energy (Figure 6.14b).Therefore, the KKL model is not appropriate to simulate quasi-static fracture as the en-ergy release rate equals the input elastic energy. Although never mentioned in the literature,this conclusion could be the reason that nobody has applied the KKL model to quasi-staticfracture.76(a) (b)(c) (d)Figure 6.14: Elastic energy release and surface energy increment with respect to thecrack tip position. (a) and (b) are for the mechanical model. (c) and (d) arefor the KKL model. G1.5 and G2.0 are the energy release rates with Ein equalto 1.5Gc and 2.0Gc. G1.5c and G2.0c are the critical energy release ratees withEin equal to 1.5Gc and 2.0Gc.776.7 Interaction between the Crack and a GrainBoundaryFracture in heterogeneous materials has always been a challenging topic in phase-field fracturemodelling. In the mechanics community, various attempts have been made to incorporatematerial heterogeneity, and some of the results [81] even showed quantitative agreement withthe theory. On the contrary, in the physics community, the only application to considermaterial heterogeneity is by Henry [106], who demonstrated that the KKL model is notsuitable for simulating fracture in heterogeneous materials. Despite the similarities betweenthe mechanical model and the KKL model, both communities have different opinions aboutwhether the phase-field model can account for material heterogeneity. Therefore, to under-stand the ability of the phase-field models in capturing fracture phenomena in heterogeneousmaterials, the mechanical model and the KKL model were adapted for studying the role ofa grain boundary during crack propagation.The simulation setting is the same as the one used by Hansen-Do¨rr et al. [81]. A 45°tilted grain boundary separates two grains, and two types of critical energy release rate wereassigned to the bulk (Gbulkc ) and the grain boundary (Gintc ). Since the grain boundary widthhas to be wide to have accurate results, the grain boundary width was set to 10, which is overtwice the crack interface thickness. In the simulations, Gbulkc was set to 0.5, and various Gintcwere tested. To be consistent with the other simulations in this thesis, the mode III loadingcondition was studied. Since the energy release rate is overestimated in both models becauseof its sensitiveness with Ein, Ein was set to the lowest value to allow the crack propagationin the bulk so that the energy release rate is the closest to the analytical value. Therefore,Ein is 1.33Gbulkc in the mechanical model and 1.03Gbulkc in the KKL model.Figure 6.15 compares the crack behaviours with different grain boundary to bulk tough-ness ratios in both models. Both models correctly predict the crack deflection and penetrationbehaviour. The critical Gintc to Gbulkc ratio for crack deflection is less than 0.4 in the mechan-ical model and less than 0.3 in the KKL model. Compared with the theoretical ratio fromthe theory by He and Hutchinson [82], which is 0.85, both models deviate from the theory.The discrepancy between the simulations and the theoretical prediction is caused by theinteraction between the crack interface and the grain boundary. Because of that, there arealso some unphysical phenomena in the process.The whole process contains three stages. The first stage is when the crack interface startstouching the grain boundary. In both models, the crack bends slightly downwards to meet78(a) (b)(c) (d)Figure 6.15: Crack path at the boundary with different Gintc to Gbulkc ratios. (a) and(b) are for the mechanical model with ratios equal to 0.4 and 0.7. (c) and (d)are for the KKL model with ratios equal to 0.3 and 0.6. The transparent grayband represents the grain boundary.the grain boundary, which is because the diffuse crack interface spreads to the grain boundaryeven before φ = 0 touching it. In the mechanical model, because φ ahead of the crack frontis less than 1, the crack interface spreads to the grain boundary even earlier than that inthe KKL model. As a result, the contribution of Gintc decreases the crack surface energy. To79further reduce the surface energy, the crack tends to propagate towards the grain boundaryto maximize its interface area inside the boundary. In reality, however, it should not happenas the crack path should remain unchanged before touching the grain boundary.The second stage is when most of the crack interface is inside the boundary. When thecrack propagates along the x-axis inside the grain boundary, the crack surface energy, i.e.,critical energy release rate, should be Gintc . Since the critical energy release rate varies withthe propagation direction, when the crack deflects at the grain boundary, Gintc increases andbecomes a function of the deflection angle. Therefore, theoretically, the crack deflectionoccurs when the critical energy release rate deflecting at the boundary is lower than thatpropagating straight in the bulk. However, as discussed in Section 6.6, the critical energyrelease rate is overestimated in both models, and the overestimation becomes larger whenincreasing the load. The overestimation of the critical energy release rate inside the grainboundary can be indicated from the crack interface thickness in Figures 6.15a and 6.15c.Since the KKL model is more sensitive to the load, the crack interface widens significantlyinside the crack as shown in Figure 6.15c. As a result, the critical energy release rate insidethe grain boundary is overestimated in both models but more in the KKL model, which isthe main reason for the deviations from the theory by He and Hutchinson. Apart from that,in the phase-field models, the grain boundary is defined as a thick band. Because of it, thecrack deflection angle tends to be less than the grain boundary angle as the critical energyrelease rate is lower with a smaller angle. Strictly speaking, during deflection, the deflectionangle equals the grain boundary angle. Here, it is assumed that the crack deflects at theboundary when the difference between the deflection angle and the grain boundary angle iswithin 10°.In the end, even when the crack deflects at the grain boundary, it will finally leave theboundary as shown in Figures 6.15a and 6.15c. However, theoretically, it should not happen.For a sharp interface crack, once the crack propagates along the boundary, it will never gointo the bulk since the critical energy release rate along the boundary is lower than that in thebulk. This unphysical penetration of the crack into the bulk is similar to the finding by Henry[106]. As in the first stage, this unphysical phenomenon is caused by the interaction betweenthe crack interface and the grain boundary. Since the crack deflection angle is smaller thanthe grain boundary angle, as the crack propagates, part of the crack interface will touch thebulk. As a result, the crack surface energy increases because of the contribution of Gbulkc . Asthe crack propagates, more of the crack interface goes into the bulk, and the crack surfaceenergy keeps increasing. When the surface energy reaches the value allowing the crack to80propagate straight inside the bulk, the crack will leave the grain boundary. In conclusion,because of the large deviations from the theory and the unphysical phenomena, the sharpinterface models might be better options than the mechanical model and the KKL model forsimulating the interaction between the crack and the grain boundary.81Chapter 7Simulations of Dynamic CrackPropagation7.1 IntroductionThe previous chapter demonstrates that the KKL model exhibits unphysical phenomena inquasi-static simulations. In dynamic fracture, both models have been successfully appliedto reproduce experimental results and theoretical predictions. To compare their abilitiesin capturing the branching phenomenon in dynamic fracture, 2D dynamic simulations werecarried out. As mentioned in Chapter 2, some of the assumptions adopted in 2D simulationsare not valid in experiments, which are in 3D. Therefore, 2D dynamic branching theoriesespecially those based on the LEFM would be ideal to be compared with the phase-fieldsimulations. Among the theories, Adda-Bedia’s theory [5] is based on the LEFM and con-siders the antiplane shear loading condition, which is also the loading condition used in thecurrent work. Besides, it was employed by Karma and Lobkovsky [35] to verify the resultsof the KKL model. Adda-Bedia’s theory, therefore, was adopted to evaluate the results ofboth models. In practice, branching happens either when the elastic wave interacts with thecrack tip or when the crack velocity reaches a fraction of the material wave speed [109]. Inthe simulations, since the loading rate is zero, the only elastic wave is that radiated from thecrack tip, which can interact with the crack only after being reflected from the boundary.Since Adda-Bedia did not consider the branching induced by the interaction with the elasticwave, all the analyses were carried out before the radiated elastic wave is reflected from theboundary. The parameters used in the simulations were determined in Chapter 5 and can befound in Table 7.1 unless stated explicitly.82Table 7.1: Parameters in 2D dynamic simulationsSymbol Parameter ValueLx Domain length 186Ly Domain width 140lc Initial crack length 53.8C Shear modulus 1Gc Critical energy release rate 0.5Ein Input elastic energy 1δ Crack interface thickness 4l Internal length scale in the mechanical model 1κ Coefficient before the gradient term in the KKL model 1.26Ec Critical elastic energy to break the material 0.05η Residual stiffness 10−10ρ Material density 1dx Mesh size 0.1dt Time step size 0.02∆φthreshold Threshold value for φ in the iteration methods 10−5∆uthreshold Threshold value for u in the iteration methods 10−77.2 Investigations of the Branching Phenomena7.2.1 Branching VelocityIn quasi-static fracture, crack propagation is extremely slow such that the velocity is zero.In dynamic fracture, however, the crack velocity is one of the most important features.During dynamic crack propagation, the crack first accelerates as the elastic energy flows intothe kinetic energy. After reaching the critical velocity, crack bifurcation happens. In themeantime, the energy that flows into the crack tip is enough to create the surface energy oftwo branches. Since the crack tip becomes blunt before bifurcation, more energy flows intothe surface energy. As a result, less energy transforms into kinetic energy, and the crackvelocity decreases before branching. After branching, the microcracks developed at the cracktip region lead to velocity oscillation.Figure 7.1a compares the crack tip velocity with respect to the crack tip position for both83models. Both models predict the general trend of dynamic crack propagation. The dottedlines in Figures 7.1a and 7.1b indicate the branching positions. Before branching, the crackaccelerates, and then the velocity decreases during crack tip blunting. After branching, thevelocity becomes unstable. However, compared with the mechanical model, the KKL modelpredicts higher crack velocity. Naturally, one might think that it is because of the differentmagnitudes of the peak elastic energy at the crack tip region since the peak elastic energyacts as the ‘force’ in crack propagation. Figure 7.2a compares the elastic energy in the crackplane for both models. The lower peak elastic energy in the KKL model suggests that itcannot be the reason for the higher crack velocity. The other possible explanation is that thephase-field profiles behave differently even with the same peak elastic energy because of thedifferent phase-field formulations. Apart from the peak elastic energy, it is also noticeablethat in the KKL model, the elastic energy has a secondary spike inside the crack, which iscaused by the delay in the movement of the smaller φ as shown in Figure 7.2b. Since thesmaller φ that lags behind is negligible in magnitude and moves with the crack front, thesecondary spike in the elastic energy is negligible. As in the quasi-static case, because of thedamage accumulation, φ ahead of the crack front in the mechanical model is less than 1 asshown in Figure 7.2b. As a result, the elastic energy far from the crack tip is slightly lowerthan that in the KKL model as shown in Figure 7.2a.Apart from the crack velocity, because of the different phase-field formulations, thebranching position is also quite different as shown in Figure 7.1b. However, in Adda-Bedia’stheory, the crack is assumed to be propagating at a constant velocity before branching, andthe acceleration process is ignored. As a result, neither the actual propagation velocity northe branching position can be known from the theory. Therefore, the critical branchingvelocity is the only feature that can be compared with Adda-Bedia’s theory.In Adda-Bedia’s theory, the critical branching velocity is 0.39Cs. Before comparing withthe phase-field simulations, it has to be mentioned that in the theory, the critical branchingvelocity is derived based on a limiting case where the velocity immediately after branchingis assumed to be zero. In the phase-field simulations, however, that is not true as shown inFigure 7.1a. According to Adda-Bedia, with the increase of the velocity after branching, thecritical branching velocity should become higher. Therefore, the actual branching velocityshould be above 0.39Cs. As shown in Figure 7.1a, the branching velocity is 0.28Cs in themechanical model and 0.38Cs in the KKL model. As mentioned in Section 5.4, the branchingvelocity in the mechanical model is sensitive to φfront. Considering the 20% error induced byφfront, the branching velocity in the mechanical model should range from 0.22Cs to 0.28Cs.84(a) (b)Figure 7.1: (a) The crack velocity with respect to the crack tip position and (b) thecontour of φ profile with level set to 0.5 both models. The dotted lines in both(a) and (b) correspond to the branching positions.(a) (b)Figure 7.2: (a) The elastic energy profiles and (b) φ profiles in the crack plane for themechanical model and the KKL model.85Compared with the theoretical value, both models underestimate the branching velocity, butthe KKL model is more consistent with Adda-Bedia’s theory.7.2.2 Branching AngleIn Adda-Bedia’s theory, the branching angle is 79.2°, which is also derived based on a limitingcase. Figure 7.3 compares the crack shapes of both models. The branching angle is 47.3° in themechanical model and 64.5° in the KKL model. Compared with the theoretical prediction,i.e., 79.2°, both models underestimate the branching angle, but the KKL model is moreconsistent with the theory.Figure 7.3: Contour lines of φ = 0.5 in the mechanical model and the KKL model.The dotted lines are used to determine the branching angles.7.2.3 Branching Energy Release RateAn alternative theory to be compared with is the one by Eshelby [4]. Eshelby suggestedthat the crack bifurcates when the energy release rate becomes twice the critical energyrelease rate of a single crack. Then by studying the limiting case, where the branchingangle and the velocity after branching are zero, the branching velocity was obtained. In the86limiting case, the energy release rate for branching is 2Gc. Although based on the same idea,in Adda-Bedia’s theory, the branching angle is no longer assumed to be zero. Therefore,the energy release rate for branching becomes a function of the branching angle (λ) givenas: 2 ((360− λ)/λ)λ/360Gc [5]. Although the branching velocity in Eshelby’s theory is lessaccurate compared with that in Adda-Bedia’s theory, the fundamental ideas of both theoriesare the same as they are all based on the critical energy condition for branching.Unlike the energy release rate in quasi-static fracture, in dynamic fracture, it cannotbe calculated directly from the elastic energy as the kinetic energy is also involved in thepropagation process. However, since the energy release rate equals the critical energy releaserate according to the Griffith theory, the critical energy release rate can be evaluated instead.As in the quasi-static case, the critical energy release rate can be calculated easily from theslope of the surface energy curves.Figure 7.4 compares the critical energy release rate or equivalently the energy release ratewith respect to the crack tip position in both models. In the beginning, the crack velocity isclose to 0. Therefore, the critical energy release rate theoretically should be Gc. However, inthe mechanical model, the initial critical energy release rate is overestimated and is around1.5Gc. On the contrary, the initial value in the KKL model is quite accurate. Because ofthe overestimation of the initial critical energy release rate in the mechanical model, thecurve is always higher compared with that in the KKL model as shown in Figure 7.4. Withthe initial critical energy release rate obtained from Figure 7.4, the branching angle, andthe relation between the branching angle and the critical energy release rate, the theoreticalenergy release rate for branching then becomes 3.84Gc for the mechanical model and 2.63Gcfor the KKL model. With the branching positions determined from Figure 7.1a and thecritical energy release rate curves in Figure 7.4, the critical energy release rate for branchingin the phase-field simulations is derived, which is 3.68Gc in the mechanical model and 3.9Gcfor the KKL model. Compared with the theoretical values, i.e., 3.84Gc and 2.63Gc, themechanical model underestimates the energy by 4%, while the KKL model overestimates thevalue by 48%. The significant overestimation of the critical branching energy in the KKLmodel suggests that the mechanical model is more consistent with Eshelby’s theory.87Figure 7.4: The critical energy release rate with respect to the crack tip position. Themarkers indicate the critical value for branching. The red dotted lines indicatethe positions of the curves.7.3 Influence of the Input Elastic EnergyTo understand the robustness of the models, parameter analysis was carried out in thissection. In the theories by Eshelby and Adda-Bedia, the critical branching conditions areindependent of the load, which is Ein in the current work. However, as demonstrated inChapter 6, both models are sensitive to Ein in the quasi-static setting. Therefore, it is worthknowing whether these models are also sensitive to Ein in dynamic fracture. Two input elasticenergies, i.e., 2.0Gc and 2.5Gc, were investigated in this section. The rest of the parameterscan be found in Table 7.1.Figure 7.5 compares the velocity with respect to the crack tip position for both modelswith different Ein. In both models, the velocity increases slightly with the increase of theinput elastic energy, because more energy flows into the kinetic energy. For the branchingvelocity, with the increase of Ein, it only increases by 3.4% in the mechanical model anddecreases by 5.3% in the KKL model as indicated by the dotted lines. Therefore, consideringthe errors induced by the selections of φfront and the domain width, the branching velocityis insensitive to Ein in both models.As for the branching angle, in both models, it increases slightly, i.e., less than 10°, with88(a) (b)Figure 7.5: The crack velocity with respect to the crack tip position in the (a) mechan-ical model and the (b) KKL model with different Ein. The dotted lines indicatethe values of the critical branching velocities.the increase of Ein as shown in Figure 7.6. The inner part of the branch, however, does notmove when increasing Ein. Therefore, the increase in the branching angle with Ein is causedby the thickening of the crack tip. Since the variance of the branching angle with Ein isnegligible compared with the error introduced when selecting the dotted lines, the branchingangle is insensitive to Ein.At last, the influence of Ein on the critical energy release rate was evaluated. Figure 7.7compares the critical energy release rate with respect to the crack tip position with differentEin. In both models, the initial critical energy release rate is insensitive to the input elasticenergy. However, the critical energy release rate increases with the increase of Ein becausemore elastic energy transforms into the crack surface energy. As a result, the crack tipthickens. The critical energy release rates for branching are shown in the dotted lines inFigure 7.7. Since the branching angle is insensitive to Ein, the critical energy release ratefor branching based on Eshelby’s theory does not change, which is 3.84Gc for the mechanicalmodel and 2.63Gc for the KKL model. In the mechanical model, with the increase of Ein,the discrepancy with the theory decreases from 4% to 2%. Meanwhile, in the KKL model,the discrepancy increases from 48% to 69%. Therefore, the critical energy release rate forbranching is sensitive to Ein in the KKL model.89(a) (b)Figure 7.6: Contour lines of φ = 0.5 with different Ein in the (a) mechanical modeland (b) the KKL model.(a) (b)Figure 7.7: The critical energy release rate with respect to the crack tip position withdifferent Ein for the (a) mechanical model and the (b) KKL model.In conclusion, both the mechanical model and the KKL model agree well with dynamicfracture theories. In terms of predicting the branching phenomena by Adda-Bedia, the KKL90model is more accurate. Compared with the energetic theory by Eshelby, the mechanicalmodel showed a better agreement. As mentioned earlier, these theories are for the limitingcase in which some assumptions are not valid in the phase-field simulations. Therefore,comparisons with experiments are required to evaluate the accuracy of the models.91Chapter 8Summary and Future Work8.1 SummaryTwo representative phase-field fracture models in the mechanics and the physics community,i.e., Miehe et al.’s model [1] and the KKL model, were compared in the quasi-static andthe dynamic setting. The simulation settings and the numerical techniques for the modelimplementation were illustrated in Chapter 4. Through various comparisons in Chapters 6and 7, the strengths and limitations of these models were elucidated. The main findings aresummarized below.In the quasi-static case, both models were compared with the Griffith theory. It wasfound that the KKL model cannot be applied to quasi-static fracture as the energy releaserate depends on the load (Ein) instead of the material parameter (Gc). The mechanicalmodel, on the other hand, is more appropriate for simulating quasi-static fracture. However,the mechanical model also has some limitations. The crack tip curvature effect, although,exists in reality, cannot be captured. The sensitiveness of the phase-field profile with the loadleads to the overestimation of the energy release rate. For crack propagation in a bi-grainsystem, because of the interaction between the crack interface and the grain boundary, bothmodels deviate from the theory and exhibit unphysical phenomena.As for the dynamic fracture, the KKL model agrees well with the branching angle andvelocity in Adda-Bedia’s theory [5]. The mechanical model, however, underestimates thebranching angle and velocity by 40% and 28%. Apart from Adda-Bedia’s theory, Eshelby’stheory [4] was also adopted to identify the critical branching condition. It was observedthat the mechanical model agrees well with the theory, while the KKL model overestimatesthe theoretical value by 48%. However, since these theories are for the limiting case wheresome assumptions are not valid in phase-field simulations, comparisons with experiments arenecessary to evaluate the accuracy of the models.Therefore, this work illustrated that for quasi-static simulations, the mechanical model isa better choice. For the dynamic case, more investigations into both models are necessary.928.2 Future WorkThe following extensions of the study are suggested:• The selection of the domain size is restricted by the computational cost in the quasi-static simulations. To enable the use of larger domain size, the adaptive meshingtechnique could be used, where only the interface has the fine mesh.• The current work is restricted to the mode III fracture. However, in other fracturemodes, the findings in this project may not be applicable. Therefore, simulations withother fracture modes are required.• In dynamic simulations, the fracture phenomena could be compared with experimentsto evaluate the accuracy of the model predictions. Besides, to compare with experi-ments, these models should be implemented in 3D.• Instead of using φfront to identify the crack front position, alternative methods can beused such that the velocity and crack shape could be captured with higher accuracy.For example, the crack front can be identified by the position with the peak elasticenergy.93Bibliography[1] C. Miehe, M. Hofacker, and F. Welschinger, “A phase field model forrate-independent crack propagation: Robust algorithmic implementation based onoperator splits,” Computer Methods in Applied Mechanics and Engineering, vol. 199,no. 45-48, pp. 2765–2778, 2010.[2] A. Karma, D. A. Kessler, and H. Levine, “Phase-field model of mode iii dynamicfracture,” Physical Review Letters, vol. 87, no. 4, p. 045501, 2001.[3] A. A. Griffith, “Vi. the phenomena of rupture and flow in solids,” PhilosophicalTransactions of the Royal Society of London. Series A, Containing Papers of aMathematical or Physical Character, vol. 221, no. 582-593, pp. 163–198, 1921.[4] J. D. Eshelby, “Energy relations and the energy-momentum tensor in continuummechanics,” in Fundamental Contributions to the Continuum Theory of EvolvingPhase Interfaces in Solids, pp. 82–119, Springer, 1999.[5] M. Adda-Bedia, “Brittle fracture dynamics with arbitrary paths. ii. dynamic crackbranching under general antiplane loading,” Journal of the Mechanics and Physics ofSolids, vol. 52, no. 6, pp. 1407–1420, 2004.[6] R. P. Reed, The economic effects of fracture in the United States, vol. 647. USDepartment of Commerce, National Bureau of Standards, 1983.[7] K. Ravi-Chandar, Dynamic fracture. Elsevier, 2004.[8] G. R. Irwin, “Analysis of stresses and strains near the end of a crack transversing aplate,” Journal of Aplied Mechanics, Transactions ASME, vol. 24, pp. 361–364, 1957.[9] B. Lawn, Fracture of brittle solids. Cambridge University Press, 1993.[10] F. Erdogan and G. C. Sih, “On the crack extension in plates under plane loading andtransverse shear,” Journal of Basic Engineering, vol. 85, no. 4, pp. 519–525, 1963.94[11] G. C. Sih and B. Macdonald, “Fracture mechanics applied to engineeringproblems-strain energy density fracture criterion,” Engineering Fracture Mechanics,vol. 6, no. 2, pp. 361–386, 1974.[12] M. A. Hussain, S. L. Pu, and J. Underwood, “Strain energy release rate for a crackunder combined mode i and mode ii,” in Fracture Analysis: Proceedings of the 1973National Symposium on Fracture Mechanics, Part II, ASTM International, 1974.[13] N. Moe¨s, J. Dolbow, and T. Belytschko, “A finite element method for crack growthwithout remeshing,” International Journal for Numerical Methods in Engineering,vol. 46, no. 1, pp. 131–150, 1999.[14] T. Belytschko and T. Black, “Elastic crack growth in finite elements with minimalremeshing,” International Journal for Numerical Methods in Engineering, vol. 45,no. 5, pp. 601–620, 1999.[15] P.-O. Bouchard, F. Bay, Y. Chastel, and I. Tovena, “Crack propagation modellingusing an advanced remeshing technique,” Computer Methods in Applied Mechanicsand Engineering, vol. 189, no. 3, pp. 723–742, 2000.[16] P.-O. Bouchard, F. Bay, and Y. Chastel, “Numerical modelling of crack propagation:automatic remeshing and comparison of different criteria,” Computer Methods inApplied Mechanics and Engineering, vol. 192, no. 35-36, pp. 3887–3908, 2003.[17] T. K. Hellen, “The finite element calculation of stress intensity factors using energytechniques,” in The Second International Conference on Structural Mechanics inReactor Technology, (Berin, Germany), IASMiRT, 1973.[18] D. M. Parks, “A stiffness derivative finite element technique for determination ofcrack tip stress intensity factors,” International Journal of Fracture, vol. 10, no. 4,pp. 487–502, 1974.[19] T. K. Hellen, “On the method of virtual crack extensions,” International Journal forNumerical Methods in Engineering, vol. 9, no. 1, pp. 187–207, 1975.[20] P. Destuynder, M. Djaoua, and S. Lescure, “Quelques remarques sur la me´canique dela rupture e´lastique,” Journal de Me´canique The´orique et Applique´e, vol. 2, no. 1,pp. 113–115, 1983.95[21] L. B. Freund, Dynamic fracture mechanics. Cambridge University Press, 1998.[22] E. H. Yoffe, “Lxxv. the moving griffith crack,” The London, Edinburgh, and DublinPhilosophical Magazine and Journal of Science, vol. 42, no. 330, pp. 739–750, 1951.[23] S. P. Patil and Y. Heider, “A review on brittle fracture nanomechanics by all-atomsimulations,” Nanomaterials, vol. 9, no. 7, p. 1050, 2019.[24] R. Kobayashi, “Modeling and numerical simulations of dendritic crystal growth,”Physica D: Nonlinear Phenomena, vol. 63, no. 3-4, pp. 410–423, 1993.[25] W. J. Boettinger, J. A. Warren, C. Beckermann, and A. Karma, “Phase-fieldsimulation of solidification,” Annual Review of Materials Research, vol. 32, no. 1,pp. 163–194, 2002.[26] W. J. Boettinger, A. A. Wheeler, B. T. Murray, and G. B. McFadden, “Prediction ofsolute trapping at high solidification rates using a diffuse interface phase-field theoryof alloy solidification,” Materials Science and Engineering: A, vol. 178, no. 1-2,pp. 217–223, 1994.[27] B. Bo¨ttger, J. Eiken, and I. Steinbach, “Phase field simulation of equiaxedsolidification in technical alloys,” Acta Materialia, vol. 54, no. 10, pp. 2697–2704,2006.[28] B. Zhu and M. Militzer, “3d phase field modelling of recrystallization in a low-carbonsteel,” Modelling and Simulation in Materials Science and Engineering, vol. 20, no. 8,p. 085011, 2012.[29] T. Takaki, Y. Hisakuni, T. Hirouchi, A. Yamanaka, and Y. Tomita, “Multi-phase-fieldsimulations for dynamic recrystallization,” Computational Materials Science, vol. 45,no. 4, pp. 881–888, 2009.[30] Y. Suwa, Y. Saito, and H. Onodera, “Phase field simulation of stored energy driveninterface migration at a recrystallization front,” Materials Science and Engineering:A, vol. 457, no. 1-2, pp. 132–138, 2007.[31] D.-H. Yeon, P.-R. Cha, and J.-K. Yoon, “A phase field study for ferrite–austenitetransitions under para-equilibrium,” Scripta Materialia, vol. 45, no. 6, pp. 661–668,2001.96[32] R. G. Thiessen, J. Sietsma, T. A. Palmer, J. W. Elmer, and I. M. Richardson,“Phase-field modelling and synchrotron validation of phase transformations inmartensitic dual-phase steel,” Acta Materialia, vol. 55, no. 2, pp. 601–614, 2007.[33] I. Steinbach and M. Apel, “The influence of lattice strain on pearlite formation infe–c,” Acta Materialia, vol. 55, no. 14, pp. 4817–4822, 2007.[34] M. Militzer, M. G. Mecozzi, J. Sietsma, and S. Van der Zwaag, “Three-dimensionalphase field modelling of the austenite-to-ferrite transformation,” Acta Materialia,vol. 54, no. 15, pp. 3961–3972, 2006.[35] A. Karma and A. E. Lobkovsky, “Unsteady crack motion and branching in aphase-field model of brittle fracture,” Physical Review Letters, vol. 92, no. 24,p. 245510, 2004.[36] C. Kuhn and R. Mu¨ller, “A continuum phase field model for fracture,” EngineeringRracture Mechanics, vol. 77, no. 18, pp. 3625–3634, 2010.[37] A. Schlu¨ter, A. Willenbu¨cher, C. Kuhn, and R. Mu¨ller, “Phase field approximation ofdynamic brittle fracture,” Computational Mechanics, vol. 54, no. 5, pp. 1141–1161,2014.[38] L.-Q. Chen, “Phase-field models for microstructure evolution,” Annual Review ofMaterials Research, vol. 32, no. 1, pp. 113–140, 2002.[39] G. A. Francfort and J.-J. Marigo, “Revisiting brittle fracture as an energyminimization problem,” Journal of the Mechanics and Physics of Solids, vol. 46,no. 8, pp. 1319–1342, 1998.[40] L. Kachanov, Introduction to continuum damage mechanics, vol. 10. Springer Science& Business Media, 1986.[41] J. Mazars and G. Pijaudier-Cabot, “Continuum damage theory—application toconcrete,” Journal of Engineering Mechanics, vol. 115, no. 2, pp. 345–365, 1989.[42] Z. P. Bazant, “Mechanics of distributed cracking,” Applied Mechanics Reviews,vol. 39, no. 5, pp. 675–705, 1986.97[43] L. Ambrosio and V. M. Tortorelli, “Approximation of functional depending on jumpsby elliptic functional via γ-convergence,” Communications on Pure and AppliedMathematics, vol. 43, no. 8, pp. 999–1036, 1990.[44] D. B. Mumford and J. Shah, “Optimal approximations by piecewise smooth functionsand associated variational problems,” Communications on Pure and AppliedMathematics, vol. 42, no. 5, pp. 577–685, 1989.[45] B. Bourdin, G. A. Francfort, and J.-J. Marigo, “Numerical experiments in revisitedbrittle fracture,” Journal of the Mechanics and Physics of Solids, vol. 48, no. 4,pp. 797–826, 2000.[46] B. Bourdin, G. A. Francfort, and J.-J. Marigo, “The variational approach tofracture,” Journal of Elasticity, vol. 91, no. 1-3, pp. 5–148, 2008.[47] C. Kuhn and R. Mu¨ller, “A phase field model for fracture,” Proceedings in Appliedmathematics and Mechanics, vol. 8, no. 1, pp. 10223–10224, 2008.[48] C. Miehe, F. Welschinger, and M. Hofacker, “Thermodynamically consistentphase-field models of fracture: Variational principles and multi-field feimplementations,” International Journal for Numerical Methods in Engineering,vol. 83, no. 10, pp. 1273–1311, 2010.[49] T. Linse, P. Hennig, M. Ka¨stner, and R. de Borst, “A convergence study ofphase-field models for brittle fracture,” Engineering Fracture Mechanics, vol. 184,pp. 307–318, 2017.[50] H. Amor, J.-J. Marigo, and C. Maurini, “Regularized formulation of the variationalbrittle fracture with unilateral contact: Numerical experiments,” Journal of theMechanics and Physics of Solids, vol. 57, no. 8, pp. 1209–1229, 2009.[51] M. J. Borden, Isogeometric analysis of phase-field models for dynamic brittle andductile fracture. PhD thesis, The University of Texas at Austin, 2012.[52] K. Pham, H. Amor, J.-J. Marigo, and C. Maurini, “Gradient damage models andtheir use to approximate brittle fracture,” International Journal of DamageMechanics, vol. 20, no. 4, pp. 618–652, 2011.98[53] C. Kuhn and R. Mu¨ller, “Interpretation of parameters in phase field models forfracture,” The Proceedings in Applied Mathematics and Mechanics, vol. 12, no. 1,pp. 161–162, 2012.[54] C. Kuhn and R. Mu¨ller, “Crack nucleation in phase field fracture models,” in The13th International Conference on Fracture, vol. 579, (Beijing, China), InternationalCongress on Fracture, Citeseer, 2013.[55] T. T. Nguyen, J. Yvonnet, M. Bornert, C. Chateau, K. Sab, R. Romani, andR. Le Roy, “On the choice of parameters in the phase field method for simulatingcrack initiation with experimental validation,” International Journal of Fracture,vol. 197, no. 2, pp. 213–226, 2016.[56] X. Zhang, C. Vignes, S. W. Sloan, and D. Sheng, “Numerical evaluation of thephase-field model for brittle fracture with emphasis on the length scale,”Computational Mechanics, vol. 59, no. 5, pp. 737–752, 2017.[57] K. H. Pham, K. Ravi-Chandar, and C. M. Landis, “Experimental validation of aphase-field model for fracture,” International Journal of Fracture, vol. 205, no. 1,pp. 83–101, 2017.[58] B. B. Yin and L. W. Zhang, “Phase field method for simulating the brittle fracture offiber reinforced composites,” Engineering Fracture Mechanics, vol. 211, pp. 321–340,2019.[59] J.-Y. Wu, “A unified phase-field theory for the mechanics of damage and quasi-brittlefailure,” Journal of the Mechanics and Physics of Solids, vol. 103, pp. 72–99, 2017.[60] J.-Y. Wu, V. P. Nguyen, C. T. Nguyen, D. Sutula, S. Bordas, and S. Sinaie, “Phasefield modeling of fracture,” Advances in Applied Mechancis: Multi-scale Theory andComputation, vol. 52, 2018.[61] C. Kuhn, A. Schlu¨ter, and R. Mu¨ller, “On degradation functions in phase fieldfracture models,” Computational Materials Science, vol. 108, pp. 374–384, 2015.[62] E. Tanne´, T. Li, B. Bourdin, J.-J. Marigo, and C. Maurini, “Crack nucleation invariational phase-field models of brittle fracture,” Journal of the Mechanics andPhysics of Solids, vol. 110, pp. 80–99, 2018.99[63] M. Ambati, T. Gerasimov, and L. De Lorenzis, “A review on phase-field models ofbrittle fracture and a new fast hybrid formulation,” Computational Mechanics,vol. 55, no. 2, pp. 383–405, 2015.[64] G. Lancioni and G. Royer-Carfagni, “The variational approach to fracture mechanics.a practical application to the french panthe´on in paris,” Journal of Elasticity, vol. 95,no. 1-2, pp. 1–30, 2009.[65] F. Freddi and G. Royer-Carfagni, “Regularized variational theories of fracture: aunified approach,” Journal of the Mechanics and Physics of Solids, vol. 58, no. 8,pp. 1154–1174, 2010.[66] T. T. Nguyen, J. Re´thore´, and M.-C. Baietto, “Phase field modelling of anisotropiccrack propagation,” European Journal of Mechanics-A/Solids, vol. 65, pp. 279–288,2017.[67] B. Li, C. Peco, D. Milla´n, I. Arias, and M. Arroyo, “Phase-field modeling andsimulation of fracture in brittle materials with strongly anisotropic surface energy,”International Journal for Numerical Methods in Engineering, vol. 102, no. 3-4,pp. 711–727, 2015.[68] B. Li and C. Maurini, “Crack kinking in a variational phase-field model of brittlefracture with strongly anisotropic surface energy,” Journal of the Mechanics andPhysics of Solids, vol. 125, pp. 502–522, 2019.[69] S. Teichtmeister, D. Kienle, F. Aldakheel, and M.-A. Keip, “Phase field modeling offracture in anisotropic brittle solids,” International Journal of Non-Linear Mechanics,vol. 97, pp. 1–21, 2017.[70] M. J. Borden, C. V. Verhoosel, M. A. Scott, J. Hughes, T., and C. M. Landis, “Aphase-field description of dynamic brittle fracture,” Computer Methods in AppliedMechanics and Engineering, vol. 217, pp. 77–95, 2012.[71] C. Kuhn, A. Schlu¨ter, and R. Mu¨ller, “A phase field approach for dynamic fracture,”Proceedings in Applied Mathematics and Mechanics, vol. 13, no. 1, pp. 87–88, 2013.[72] T. Li, J.-J. Marigo, D. Guilbaud, and S. Potapov, “Gradient damage modeling ofbrittle fracture in an explicit dynamics context,” International Journal for NumericalMethods in Engineering, vol. 108, no. 11, pp. 1381–1405, 2016.100[73] J. Bleyer, C. Roux-Langlois, and J.-F. Molinari, “Dynamic crack propagation with avariational phase-field model: limiting speed, crack branching and velocity-tougheningmechanisms,” International Journal of Fracture, vol. 204, no. 1, pp. 79–100, 2017.[74] M. Hofacker and C. Miehe, “Continuum phase field modeling of dynamic fracture:variational principles and staggered fe implementation,” International Journal ofFracture, vol. 178, no. 1-2, pp. 113–129, 2012.[75] J. F. Kalthoff and S. Winkler, “Failure mode transition at high rates of shearloading,” DGM Informationsgesellschaft mbH, Impact Loading and Dynamic Behaviorof Materials, vol. 1, pp. 185–195, 1988.[76] M. Z. Hossain, C.-J. Hsueh, B. Bourdin, and K. Bhattacharya, “Effective toughnessof heterogeneous media,” Journal of the Mechanics and Physics of Solids, vol. 71,pp. 15–32, 2014.[77] T. T. Nguyen, J. Yvonnet, Q.-Z. Zhu, M. Bornert, and C. Chateau, “A phase fieldmethod to simulate crack nucleation and propagation in strongly heterogeneousmaterials from direct imaging of their microstructure,” Engineering FractureMechanics, vol. 139, pp. 18–39, 2015.[78] T. T. Nguyen, J. Yvonnet, M. Bornert, C. Chateau, F. Bilteryst, and E. Steib,“Large-scale simulations of quasi-brittle microcracking in realistic highlyheterogeneous microstructures obtained from micro ct imaging,” Extreme MechanicsLetters, vol. 17, pp. 50–55, 2017.[79] D. Schneider, E. Schoof, Y. Huang, M. Selzer, and B. Nestler, “Phase-field modelingof crack propagation in multiphase systems,” Computer Methods in AppliedMechanics and Engineering, vol. 312, pp. 186–195, 2016.[80] B. Nestler, H. Garcke, and B. Stinner, “Multicomponent alloy solidification:phase-field modeling and simulations,” Physical Review E, vol. 71, no. 4, p. 041609,2005.[81] A. C. Hansen-Do¨rr, R. de Borst, P. Hennig, and M. Ka¨stner, “Phase-field modellingof interface failure in brittle materials,” Computer Methods in Applied Mechanics andEngineering, vol. 346, pp. 25–42, 2019.101[82] H. Ming-Yuan and J. W. Hutchinson, “Crack deflection at an interface betweendissimilar elastic materials,” International Journal of Solids and Structures, vol. 25,no. 9, pp. 1053–1067, 1989.[83] H. Chen, C. Zhang, Q. Lu, H. Chen, Z. Yang, Y. Wen, S. Hu, and L. Chen, “Atwo-set order parameters phase-field modeling of crack deflection/penetration in aheterogeneous microstructure,” Computer Methods in Applied Mechanics andEngineering, vol. 347, pp. 1085–1104, 2019.[84] R. Alessi, M. Ambati, T. Gerasimov, S. Vidoli, and L. De Lorenzis, “Comparison ofphase-field models of fracture coupled with plasticity,” in Advances in ComputationalPlasticity, pp. 1–21, Springer, 2018.[85] R. Alessi, J.-J. Marigo, and S. Vidoli, “Gradient damage models coupled withplasticity and nucleation of cohesive cracks,” Archive for Rational Mechanics andAnalysis, vol. 214, no. 2, pp. 575–615, 2014.[86] F. P. Duda, A. Ciarbonetti, P. J. Sa´nchez, and A. E. Huespe, “A phase-field/gradientdamage model for brittle fracture in elastic–plastic solids,” International Journal ofPlasticity, vol. 65, pp. 269–296, 2015.[87] M. Ambati, T. Gerasimov, and L. De Lorenzis, “Phase-field modeling of ductilefracture,” Computational Mechanics, vol. 55, no. 5, pp. 1017–1040, 2015.[88] R. Alessi, J.-J. Marigo, and S. Vidoli, “Gradient damage models coupled withplasticity: variational formulation and main properties,” Mechanics of Materials,vol. 80, pp. 351–367, 2015.[89] M. J. Borden, J. Hughes, T., C. M. Landis, A. Anvari, and I. J. Lee, “A phase-fieldformulation for fracture in ductile materials: Finite deformation balance lawderivation, plastic degradation, and stress triaxiality effects,” Computer Methods inApplied Mechanics and Engineering, vol. 312, pp. 130–166, 2016.[90] M. Ambati, R. Kruse, and L. De Lorenzis, “A phase-field model for ductile fracture atfinite strains and its experimental verification,” Computational Mechanics, vol. 57,no. 1, pp. 149–167, 2016.[91] C. Miehe, S. Teichtmeister, and F. Aldakheel, “Phase-field modelling of ductilefracture: a variational gradient-extended plasticity-damage theory and its102micromorphic regularization,” Philosophical Transactions of the Royal Society A:Mathematical, Physical and Engineering Sciences, vol. 374, no. 2066, p. 20150170,2016.[92] C. Kuhn, T. Noll, and R. Mu¨ller, “On phase field modeling of ductile fracture,”GAMM-Mitteilungen, vol. 39, no. 1, pp. 35–54, 2016.[93] I. S. Aranson, V. A. Kalatsky, and V. M. Vinokur, “Continuum field description ofcrack propagation,” Physical Review Letters, vol. 85, no. 1, p. 118, 2000.[94] E. A. Brener and R. Spatschek, “Fast crack propagation by surface diffusion,”Physical Review E, vol. 67, no. 1, p. 016112, 2003.[95] R. Spatschek, M. Hartmann, E. Brener, H. Mu¨ller-Krumbhaar, and K. Kassner,“Phase field modeling of fast crack propagation,” Physical Review Letters, vol. 96,no. 1, p. 015502, 2006.[96] L. O. Eastgate, J. P. Sethna, M. Rauscher, T. Cretegny, C.-S. Chen, and C. R. Myers,“Fracture in mode i using a conserved phase-field model,” Physical Review E, vol. 65,no. 3, p. 036117, 2002.[97] R. Spatschek, E. Brener, and A. Karma, “Phase field modeling of crack propagation,”Philosophical Magazine, vol. 91, no. 1, pp. 75–95, 2011.[98] H. Henry, “Study of the branching instability using a phase field model of inplanecrack propagation,” Europhysics Letters, vol. 83, no. 1, p. 16004, 2008.[99] H. Henry and M. Adda-Bedia, “Fractographic aspects of crack branching instabilityusing a phase-field model,” Physical Review E, vol. 88, no. 6, p. 060401, 2013.[100] H. Henry and H. Levine, “Dynamic instabilities of fracture under biaxial strain usinga phase field model,” Physical Review Letters, vol. 93, no. 10, p. 105504, 2004.[101] A. J. Pons and A. Karma, “Helical crack-front instability in mixed-mode fracture,”Nature, vol. 464, no. 7285, pp. 85–89, 2010.[102] V. Hakim and A. Karma, “Laws of crack motion and phase-field models of fracture,”Journal of the Mechanics and Physics of Solids, vol. 57, no. 2, pp. 342–368, 2009.103[103] H. Henry, “Study of three-dimensional crack fronts under plane stress using a phasefield model,” Europhysics Letters, vol. 92, no. 4, p. 46002, 2010.[104] V. Hakim and A. Karma, “Crack path prediction in anisotropic brittle materials,”Physical Review Letters, vol. 95, no. 23, p. 235501, 2005.[105] A. Mesgarnejad, C. Pan, R. M. Erb, S. J. Shefelbine, and A. Karma, “Crack pathselection in orientationally ordered composites,” Physical Review E, vol. 102, no. 1,p. 013004, 2020.[106] H. Henry, “Limitations of the modelling of crack propagating through heterogeneousmaterial using a phase field approach,” Theoretical and Applied Fracture Mechanics,vol. 104, p. 102384, 2019.[107] W. G. Knauss, “Stresses in an infinite strip containing a semi-infinite crack,” Journalof Applied Mechanics, vol. 33, no. 2, pp. 356–362, 1966.[108] A. T. Zehnder, “Lecture notes on fracture mechanics,” 2007.[109] J. E. Field, “Brittle fracture: its study and application,” Contemporary Physics,vol. 12, no. 1, pp. 1–31, 1971.104Appendix ASoftware Source Codes for theMechanical ModelA.1 The Mechanical Model for 2D Quasi-staticFracture1) main.c/* main loop for mechanical model*/#include <stdio.h>#include <stdlib.h>#include "initial.h"#include "pfsolver.h"#include "displsolver.h"#include "output.h"#include "matrix_vector.h"// current time t, output countint t = 0, count = 0;// total number of loops tot_T, output results every interval_tconst int tot_T = 1000, interval_T = 10;// M: width of the domain (row), N: length (column)const int M = 501, N = 501;// grid sizeconst double dx = 0.1;// threshold values for iteration methodsconst double err_dis = 1.e-8, err_pf = 1.e-7;// crack width, position relative to the middle, and length105const int init_crack_width = 1, init_crack_po = 0, init_crack_length = 15;// shear modulus mu, length scale epsi, residual stiffness etaconst double mu = 1, epsi = 1, eta = 1.e-10;// input elastic energy G(unit of Gc), fracture toughness Gcconst double G = 1.5, Gc = 0.5;// phase-field parameter phi and displacement field udouble **phi, **u;int main(int argc, char const *argv[]){// initialize phase-field parameter and displacement fieldinitial();for (t = 0; t <= tot_T; t++){// solve phase field equationpfsolver();// solve displacement equationdisplsolver();// output to filesoutput();}// free the pointersfree_dmatrix(phi, M);free_dmatrix(u, M);return 0;}2) matrix vector.h#ifndef MATRIX_VECTOR_H#define MATRIX_VECTOR_Hvoid nrerror(char error_text[]);double **dmatrix(int nrow, int ncol);double ***d3array(int n1, int n2, int n3)void free_dmatrix(double **v, int nrow);106void free_d3array(double ***v, int n1, int n2);#endif3) matrix vector.c/* functions for creating and freeing dynamic array */#include <stdio.h>#include <stddef.h>#include <stdlib.h>void nrerror(char error_text[])/* output error message */{fprintf(stderr, "%s\n", error_text);exit(1);}double **dmatrix(int nrow, int ncol)/* allocate a double matrix with nr rows and nc columns */{double **m;/* allocate pointers to rows */m = (double **) malloc((size_t) (nrow * sizeof(double*)));if (!m){nrerror("allocation failure 1 in matrix()");}/* allocate columns */for (int i = 0; i < nrow; i++){m[i] = (double *) malloc((size_t) (ncol * sizeof(double)));if (!m[i] && i == 0){nrerror("allocation failure 2 in matrix()");}}return m;107}double ***d3array(int n1, int n2, int n3)/* allocate a three dimensional double array */{double ***m;m = (double ***) malloc((size_t) (n1 * sizeof(double**)));if (!m){nrerror("allocation failure 1 in 3array()");}for (int i = 0; i < n1; i++){m[i] = (double **) malloc((size_t) (n2 * sizeof(double*)));if (!m[i] && i == 0){nrerror("allocation failure 2 in 3array()");}for (int j = 0; j < n2; j++){m[i][j] = (double *) malloc((size_t) (n3 * sizeof(double)));if (!m[0][j] && j == 0){nrerror("allocation failure 2 in 3array()");}}}return m;}void free_dmatrix(double **v, int nrow)/* free a double matrix */{for (int i = 0; i < nrow; i++){108free(v[i]);}free(v);}void free_d3array(double ***v, int n1, int n2)/* free a double 3d array */{for (int i = 0; i < n1; i++){for (int j = 0; j < n2; j++){free(v[i][j]);}free(v[i]);}free(v);}4) initial.h#ifndef INITIAL_H#define INITIAL_Hvoid initial();double Delta();#endif5) initial.c/* initial conditions for phase field and displacement field */#include <math.h>#include "matrix_vector.h"#include "initial.h"extern const int M;extern const int N;extern const double dx;109extern const int init_crack_width;extern const int init_crack_po;extern const int init_crack_length;extern const double mu;extern const double G;extern const double Gc;extern double **phi;extern double **u;void initial(){// calculate boundary displacement delta (dimensionless)double delta = Delta();// create (M*N) 2d array for phi and uphi = dmatrix(M, N);u = dmatrix(M, N);for (int i = 0; i < M; i++){for (int j = 0; j < N; j++){// initialize phase fieldif ((i >= (M-1)/2 - init_crack_width/2 && i <= (M-1)/2 +init_crack_width/2) && (j >= 0 && j <= init_crack_length/dx)){phi[i][j] = 0;}else{phi[i][j] = 1;}// initialize displacement fieldu[i][j] = 2.*delta/(M-1)*(i - (M-1)/2);}}110}// displacement at boundarydouble Delta(){// real value of input elastic energydouble G_r = G*Gc;// width of the domaindouble w = dx*M;return sqrt(G_r*w/(mu*2.));}6) pfsolver.h#ifndef PFSOLVER_H#define PFSOLVER_Hvoid pfsolver();void pf_eqsolver();double grad2_u(int i, int j);#endif7) pfsolver.c/* solve phase field equation */#include <math.h>#include <stdio.h>#include "pfsolver.h"#include "gphifunc.h"extern const int M;extern const int N;extern int t;extern double Gc;extern const double epsi;extern const double mu;extern const double dx;extern const double err_pf;111extern double **phi;extern double **u;void pfsolver(){pf_eqsolver();}// mechanical equilibrium solvervoid pf_eqsolver(){// Gauss-Seidel iterationint loop = 0;double error = 1;for (loop = 0; loop < 1000000000 && error > err_pf; loop++){error = 0;for (int i = 0; i < M; i++){for (int j = 0; j < N; j++){double a, b, c;a = 4./(dx*dx) + 1./(4.*epsi*epsi)+ mu*grad2_u(i, j)/(2.*Gc*epsi);// zero-flux boundary conditionif (i == 0){if (j == 0){b = 1./(dx*dx)*(phi[i + 1][j] + phi[i + 1][j]+ phi[i][j + 1] + phi[i][j + 1]);}else if (j == N - 1){112b = 1./(dx*dx)*(phi[i + 1][j] + phi[i + 1][j]+ phi[i][j - 1] + phi[i][j - 1]);}else{b = 1./(dx*dx)*(phi[i + 1][j] + phi[i + 1][j]+ phi[i][j + 1] + phi[i][j - 1]);}}else if (i == M - 1){if (j == 0){b = 1./(dx*dx)*(phi[i - 1][j] + phi[i - 1][j]+ phi[i][j + 1] + phi[i][j + 1]);}else if (j == N - 1){b = 1./(dx*dx)*(phi[i - 1][j] + phi[i - 1][j]+ phi[i][j - 1] + phi[i][j - 1]);}else{b = 1./(dx*dx)*(phi[i - 1][j] + phi[i - 1][j]+ phi[i][j + 1] + phi[i][j - 1]);}}else if (j == 0){b = 1./(dx*dx)*(phi[i + 1][j] + phi[i - 1][j]+ phi[i][j + 1] + phi[i][j + 1]);}else if (j == N - 1){113b = 1./(dx*dx)*(phi[i + 1][j] + phi[i - 1][j]+ phi[i][j - 1] + phi[i][j - 1]);}else{b = 1./(dx*dx)*(phi[i + 1][j] + phi[i - 1][j]+ phi[i][j + 1] + phi[i][j - 1]);}c = 1./(4.*epsi*epsi);// store the current phidouble tmp = phi[i][j];// update new phiphi[i][j] = 1./a * (b + c);// force the phi not to increaseif (phi[i][j] > tmp){phi[i][j] = tmp;}// relative difference between current and previous phiif (fabs((phi[i][j] - tmp)/tmp) >= error&& phi[i][j] > 1.e-12){error = fabs((phi[i][j] - tmp)/tmp);}}}}}// second derivative of udouble grad2_u(int i, int j){double u_x, u_y;// Zero flux on left and right, Dirichlet on top and bottomif (i == 0)114{if (j == 0 || j == N - 1){u_y = 0;}else{u_y = (u[i][j + 1] - u[i][j - 1])/(2.*dx);}u_x = (u[i + 1][j] - u[i][j])/(dx);}else if (i == M - 1){if (j == 0 || j == N - 1){u_y = 0;}else{u_y = (u[i][j + 1] - u[i][j - 1])/(2.*dx);}u_x = (u[i][j] - u[i - 1][j])/(dx);}else if (j == 0 || j == N - 1){u_x = (u[i + 1][j] - u[i - 1][j])/(2.*dx);u_y = 0;}else{u_x = (u[i + 1][j] - u[i - 1][j])/(2.*dx);u_y = (u[i][j + 1] - u[i][j - 1])/(2.*dx);}// calculate the stationary phi without strain in first step115if (t == 0){u_x = 0;u_y = 0;}return u_x*u_x + u_y*u_y;}8) displsolver.h#ifndef DISPLSOLVER_H#define DISPLSOLVER_Hvoid displsolver();void eqsolver();double f_u(int i, int j);#endif9) displsolver.c/* solve displacement equation*/#include <math.h>#include <stdio.h>#include <omp.h>#include "gphifunc.h"#include "displsolver.h"#include "matrix_vector.h"extern const int M;extern const int N;extern const int t;extern const double err_dis;extern const double mu;extern const double dx;extern double **phi;extern double **u;116void displsolver(){eqsolver();}// mechanical equilibrium solvervoid eqsolver(){// the temporary u for new step u_tmp_newdouble **u_tmp_new;u_tmp_new = dmatrix(M, N);// use Jacobi iteration to find equilibrium displacement fieldint loop = 0;double error = 1;for (loop = 0; loop < 10000000000 && error > err_dis; loop++){error = 0;#pragma omp parallel for schedule(static) reduction(max:error)for (int i = 1; i < M - 1; i++){for (int j = 0; j < N; j++){double a, b, c, d, e;b = gphi(phi[i][j]) + gphi(phi[i + 1][j]);c = gphi(phi[i][j]) + gphi(phi[i - 1][j]);// zero flux boundary at j = 0 and j = N - 1if (j == 0 && j != N - 1){a = 4.*gphi(phi[i][j]) + gphi(phi[i + 1][j])+ gphi(phi[i - 1][j]) + gphi(phi[i][j + 1])+ gphi(phi[i][1]);d = gphi(phi[i][j]) + gphi(phi[i][j + 1]);e = gphi(phi[i][j]) + gphi(phi[i][1]);u_tmp_new[i][j] = 1./a * (b*u[i + 1][j]117+ c*u[i - 1][j] + d*u[i][j + 1] + e*u[i][1]);}else if (j != 0 && j == N - 1){a = 4.*gphi(phi[i][j]) + gphi(phi[i + 1][j])+ gphi(phi[i - 1][j]) + gphi(phi[i][N - 2])+ gphi(phi[i][j - 1]);d = gphi(phi[i][j]) + gphi(phi[i][N - 2]);e = gphi(phi[i][j]) + gphi(phi[i][j - 1]);u_tmp_new[i][j] = 1./a * (b*u[i + 1][j]+ c*u[i - 1][j] + d*u[i][N - 2] + e*u[i][j - 1]);}else{a = 4.*gphi(phi[i][j]) + gphi(phi[i + 1][j])+ gphi(phi[i - 1][j]) + gphi(phi[i][j + 1])+ gphi(phi[i][j - 1]);d = gphi(phi[i][j]) + gphi(phi[i][j + 1]);e = gphi(phi[i][j]) + gphi(phi[i][j - 1]);u_tmp_new[i][j] = 1./a * (b*u[i + 1][j]+ c*u[i - 1][j] + d*u[i][j + 1] + e*u[i][j - 1]);}// compute the relative differenceif (fabs((u_tmp_new[i][j] - u[i][j])/u[i][j]) >= error&& fabs(u[i][j]) > 1.e-6 && phi[i][j] > 1.e-6){error = fabs((u_tmp_new[i][j] - u[i][j])/u[i][j]);}}}// exchange old and new u#pragma omp parallel for schedule(static)for (int i = 1; i < M - 1; i++){118for (int j = 0; j < N; j++){u[i][j] = u_tmp_new[i][j];}}}free_dmatrix(u_tmp_new, M);}10) output.h#ifndef OUTPUT_H#define OUTPUT_Hvoid output();double eepsilon(int i, int j, int k);double ssigma(int i, int j, int k);double elastic(int i, int j);#endif11) output.c/* output the stress, phase field, elastic energy every interval_T */#include <stdio.h>#include <math.h>#include "gphifunc.h"#include "output.h"extern int count;extern const int M;extern const int N;extern const int t;extern const int interval_T;extern const double mu;extern const double dx;extern double **phi;extern double **u;119void output(){if (t % interval_T == 0){printf("t = %d\n", t);FILE *fp1, *fp2, *fp3, *fp4;char name1[150], name2[150], name3[150], name4[150];sprintf(name1, "Data/phi_%d.txt", count);sprintf(name2, "Data/sigmazy_%d.txt", count);sprintf(name3, "Data/sigmazx_%d.txt", count);sprintf(name4, "Data/energy_%d.txt", count);count ++;fp1 = fopen(name1, "w");fp2 = fopen(name2, "w");fp3 = fopen(name3, "w");fp4 = fopen(name4, "w");for (int i = 1; i < M - 1; i++){for (int j = 1; j < N - 1; j++){fprintf(fp1, "%f\t", phi[i][j]);fprintf(fp2, "%f\t", ssigma(i, j, 0));fprintf(fp3, "%f\t", ssigma(i, j, 1));fprintf(fp4, "%f\t", elastic(i, j));}fprintf(fp1, "\n");fprintf(fp2, "\n");fprintf(fp3, "\n");fprintf(fp4, "\n");}if(fp1 != NULL)fclose(fp1);120if(fp4 != NULL)fclose(fp2);if(fp5 != NULL)fclose(fp3);if(fp6 != NULL)fclose(fp4);}}// calculate the straindouble eepsilon(int i, int j, int k){double ux, uy, epsilonxy;ux = (u[i + 1][j] - u[i - 1][j])/(2.*dx);uy = (u[i][j + 1] - u[i][j - 1])/(2.*dx);if(k == 0){epsilonxy = ux;}else if(k == 1){epsilonxy = uy;}return epsilonxy;}// stressdouble ssigma(int i, int j, int k){return mu*gphi(phi[i][j])*eepsilon(i, j, k);}// elastic energydouble elastic(int i, int j){return gphi(phi[i][j])*mu/2.*(eepsilon(i, j, 0)*eepsilon(i, j, 0)+ eepsilon(i, j, 1)*eepsilon(i, j ,1));121}12) gphifunc.h#ifndef FUNCTION_H#define FUNCTION_Hdouble gphi(double phi);double gphi_prime(double phi);#endif13) gphifunc.c/* interpolation functions */extern const double eta;// elastic energy interpolation function g(phi)double gphi(double phi){return phi*phi + eta;}// derivative of g(phi) with respect to phidouble gphi_prime(double phi){return 2.*phi;}14) makefileFLAG = -O2 -lm -fopenmpOBJ = main.o matrix_vector.o gphifunc.o initial.o pfsolver.odisplsolver.o output.oHEAD1 = matrix_vector.hHEAD2 = gphifunc.hHEAD3 = initial.hHEAD4 = pfsolver.hHEAD5 = displsolver.hHEAD6 = output.h122run : $(OBJ)gcc $(OBJ) -o run $(FLAG)main.o : main.c $(HEAD1) $(HEAD2) $(HEAD3) $(HEAD4) $(HEAD5) $(HEAD6)gcc -c -std=c99 main.cgphifunc.o : gphifunc.cgcc -c -std=c99 gphifunc.cmatrix_vector.o : matrix_vector.cgcc -c -std=c99 matrix_vector.cinitial.o : initial.c $(HEAD1) $(HEAD3)gcc -c -std=c99 initial.cpfsolver.o : pfsolver.c $(HEAD2) $(HEAD4)gcc -c -std=c99 pfsolver.c -fopenmpdisplsolver.o : displsolver.c $(HEAD1) $(HEAD2) $(HEAD5)gcc -c -std=c99 displsolver.c -fopenmpoutput.o : output.c $(HEAD2) $(HEAD6)gcc -c -std=c99 output.cclean :rm -r $(OBJ) rundatacl :rm Data/*A.2 The Mechanical Model for 2D Dynamic FractureThe code for the dynamic crack propagation is almost the same as that for the quasi-staticcase except that the displacement equation solver (displsolver.c) is now for the dynamicdisplacement equation and two additional parameters are introduced, i.e., the time step(dt) and the density (ρ). The new displacement solver is given below:1) displsolver.c/* dynamic displacement solver*/#include <math.h>123#include <stdio.h>#include <omp.h>#include "gphifunc.h"#include "displsolver.h"#include "matrix_vector.h"extern const int M;extern const int N;extern int t;extern const double err_dis;extern const double mu;extern const double rho;extern const double dx;extern const double dt;extern double **phi;extern double ***u;void displsolver(){// solve damped equation in the first stepif (t < 1){#pragma omp parallel for schedule(static)for (int i = 1; i < M - 1; i++){for (int j = 1; j < N - 1; j++){u[i][j][2] = u[i][j][1] + dt*f_u(i, j);}}}else{#pragma omp parallel for schedule(static)124for (int i = 1; i < M - 1; i++){for (int j = 1; j < N - 1; j++){// solve full displacement equationu[i][j][2] = -u[i][j][0] + 2.*u[i][j][1] + dt*dt*f_u(i, j);}}}// exchange old and new ufor (int j = 0; j < N; j++){u[i][j][0] = u[i][j][1];u[i][j][1] = u[i][j][2];}// zero flux boundary on left and rightu[i][0][1] = u[i][1][1];u[i][N - 1][1] = u[i][N - 2][1];u[i][0][0] = u[i][1][0];u[i][N - 1][0] = u[i][N - 2][0];}// elastic energy discretizationdouble f_u(int i, int j){double f1, f2, f3, f4, f_tot;f1 = (gphi(phi[i][j]) + gphi(phi[i + 1][j]))* (u[i + 1][j][1] - u[i][j][1]);f2 = (gphi(phi[i - 1][j]) + gphi(phi[i][j]))* (u[i][j][1] - u[i - 1][j][1]);f3 = (gphi(phi[i][j]) + gphi(phi[i][j + 1]))* (u[i][j + 1][1] - u[i][j][1]);f4 = (gphi(phi[i][j - 1]) + gphi(phi[i][j]))* (u[i][j][1] - u[i][j - 1][1]);125f_tot = 1./rho * mu/(2.*dx*dx) * (f1 - f2 + f3 - f4);return f_tot;}126Appendix BSoftware Source Codes for the KKLModelB.1 The KKL Model for 2D Quasi-static FractureApart from main.c, pfsolver.c and gphifunc.c, the code for the KKL model in thequasi-static case is the same as that in Section A.1. The files that are different from thosein the mechanical model are given below:1) main.c/* main loop for the KKL model*/#include <stdio.h>#include <stdlib.h>#include "initial.h"#include "pfsolver.h"#include "displsolver.h"#include "output.h"#include "matrix_vector.h"// current time t, output countint t = 0, count = 0;// total number of loops tot_T, output results every interval_tconst int tot_T = 1000, interval_T = 10;// M: width of the domain (row), N: length (column)const int M = 501, N = 501;// grid sizeconst double dx = 0.1;// threshold values for iteration methodsconst double err_dis = 1.e-8, err_pf = 1.e-7;127// crack width, position relative to the middle, and lengthconst int init_crack_width = 1, init_crack_po = 0, init_crack_length = 15;// shear modulus mu, coefficient before gradient phi kappa, critical// energy for breakage Ec, residual stiffness etaconst double mu = 1, kappa = 1.26, Ec = 0.05, eta = 1.e-10;// input elastic energy G(unit of Gc), fracture toughness Gcconst double G = 1.5, Gc = 0.5;// phase-field parameter phi and displacement field udouble **phi, **u;int main(int argc, char const *argv[]){// initialize phase-field parameter and displacement fieldinitial();for (t = 0; t <= tot_T; t++){// solve phase field equationpfsolver();// solve displacement equationdisplsolver();// output to filesoutput();}// free the pointersfree_dmatrix(phi, M);free_dmatrix(u, M);return 0;}2) gphifunc.c/* interpolation functions */extern const double eta;// elastic energy interpolation function g(phi)double gphi(double phi)128{return 4.*phi*phi*phi - 3.*phi*phi*phi*phi + eta;}// derivative of g(phi) with respect to phidouble gphi_prime(double phi){return 12.*(phi*phi - phi*phi*phi);}3) pfsolver.c/* solve phase field equation */#include <math.h>#include <stdio.h>#include "pfsolver.h"#include "function.h"extern const int M;extern const int N;extern int t;extern double Ec;extern const double kappa;extern const double mu;extern const double dx;extern const double err_pf;extern double **phi;extern double **u;void pfsolver(){// equilibrium solverpf_eqsolver();}// mechanical equilibrium solvervoid pf_eqsolver()129{// Gauss-Seidel iterationint loop = 0;double error = 1;for (loop = 0; loop < 1000000000 && error > err_pf; loop++){error = 0;for (int i = 0; i < M; i++){for (int j = 0; j < N; j++){double a, b;if (i == 0){// zero-flux boundary conditionif (j == 0){a = 1./4.*(phi[i + 1][j] + phi[i + 1][j]+ phi[i][j + 1] + phi[i][j + 1]);}else if (j == N - 1){a = 1./4.*(phi[i + 1][j] + phi[i + 1][j]+ phi[i][j - 1] + phi[i][j - 1]);}else{a = 1./4.*(phi[i + 1][j] + phi[i + 1][j]+ phi[i][j + 1] + phi[i][j - 1]);}}else if (i == M - 1){if (j == 0)130{a = 1./4.*(phi[i - 1][j] + phi[i - 1][j]+ phi[i][j + 1] + phi[i][j + 1]);}else if (j == N - 1){a = 1./4.*(phi[i - 1][j] + phi[i - 1][j]+ phi[i][j - 1] + phi[i][j - 1]);}else{a = 1./4.*(phi[i - 1][j] + phi[i - 1][j]+ phi[i][j + 1] + phi[i][j - 1]);}}else if (j == 0){a = 1./4.*(phi[i + 1][j] + phi[i - 1][j]+ phi[i][j + 1] + phi[i][j + 1]);}else if (j == N - 1){a = 1./4.*(phi[i + 1][j] + phi[i - 1][j]+ phi[i][j - 1] + phi[i][j - 1]);}else{a = 1./4.*(phi[i + 1][j] + phi[i - 1][j]+ phi[i][j + 1] + phi[i][j - 1]);}b = 1./4.*dx*dx/kappa*gphi_prime(phi[i][j])*(mu/2.*grad2_u(i, j) - Ec);// store the current phidouble tmp = phi[i][j];131// update new phiphi[i][j] = a - b;// force the phi not to increaseif (phi[i][j] > tmp){phi[i][j] = tmp;}// relative difference between current and previous phiif (fabs((phi[i][j] - tmp)/tmp) >= error&& phi[i][j] > 1.e-12){error = fabs((phi[i][j] - tmp)/tmp);}}}}}// second derivative of udouble grad2_u(int i, int j){double u_x, u_y;if (i == 0){if (j == 0 || j == N - 1){u_y = 0;}else{u_y = (u[i][j + 1] - u[i][j - 1])/(2.*dx);}u_x = (u[i + 1][j] - u[i][j])/(dx);}132else if (i == M - 1){if (j == 0 || j == N - 1){u_y = 0;}else{u_y = (u[i][j + 1] - u[i][j - 1])/(2.*dx);}u_x = (u[i][j] - u[i - 1][j])/(dx);}else if (j == 0 || j == N - 1){u_x = (u[i + 1][j] - u[i - 1][j])/(2.*dx);u_y = 0;}else{u_x = (u[i + 1][j] - u[i - 1][j])/(2.*dx);u_y = (u[i][j + 1] - u[i][j - 1])/(2.*dx);}// calculate the stationary phi without strain in the first stepif (t == 0){u_x = 0;u_y = 0;}return u_x*u_x + u_y*u_y;}133B.2 The KKL Model for 2D Dynamic FractureIn the 2D dynamic case for the KKL model, displsolver.c is the same as that in Section A.2and the rest of the code is the same as that in Section B.1.134
Thesis/Dissertation
eng
Materials Engineering
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Militzer, Matthias
Phase-field modelling of crack propagation
Text