High-accuracy phase-field models for brittle fracture based on a new family of degradation functions

ContentslistsavailableatScienceDirect

Journal of the Mechanics and Physics of Solids

journalhomepage:www.elsevier.com/locate/jmps

High-accuracy phase-field models for brittle fracture based on a new family of degradation functions

Juan Michael Sargado

^a,∗

, Eirik Keilegavlen

^a

, Inga Berre

^a,b

, Jan Martin Nordbotten

^a,c

aDepartment of Mathematics, University of Bergen, Allégaten 41, Bergen 5007, Norway

bChristian Michelsen Research, Fantoftvegen 38, Bergen 5072, Norway

cDepartment of Civil and Environmental Engineering, Princeton University, E-208 E-Quad, Princeton, NJ 08544, USA

a rt i c l e i n f o

Article history:

Received 8 May 2017 Revised 25 October 2017 Accepted 25 October 2017 Available online 26 October 2017 Keywords:

Fracture Phase-field Degradation function Damage

a b s t ra c t

Phase-field approaches to fracture based onenergy minimization principles have been rapidlygainingpopularityinrecentyears,andareparticularlywell-suitedforsimulating crackinitiationandgrowthincomplexfracturenetworks.Inthephase-field framework, thesurfaceenergyassociatedwithcrackformationiscalculatedbyevaluatingafunctional definedintermsofascalarorderparameteranditsgradients.Theseinturndescribethe fracturesinadiffusesensefollowingaprescribedregularizationlengthscale.Imposingsta- tionarityofthetotalenergyleadstoacoupledsystemofpartialdifferentialequationsthat enforcestressequilibriumandgovernphase-fieldevolution. Theseequationsarecoupled throughanenergydegradationfunctionthatmodelsthelossofstiffnessinthebulkmate- rialasitundergoesdamage.Inthepresentwork,weintroduceanewparametricfamilyof degradationfunctionsaimedatincreasingtheaccuracyofphase-fieldmodelsinpredicting criticalloadsassociatedwithcracknucleationaswellasthepropagationofexistingfrac- tures.Anadditionalgoalisthepreservationoflinearelasticresponseinthebulkmaterial priortofracture.Throughtheanalysisofseveralnumericalexamples,wedemonstratethe superiorityoftheproposedfamilyoffunctionstotheclassicalquadraticdegradationfunc- tionthatisusedmostoftenintheliterature.

© 2017TheAuthors.PublishedbyElsevierLtd.

ThisisanopenaccessarticleundertheCCBYlicense.

(http://creativecommons.org/licenses/by/4.0/)

1. Introduction

Theaccuratesimulationoffractureevolutioninsolidsisamajorchallengeforcomputationalalgorithms,inlargepartdue tocrackpathsthataregenerallyunknownapriori.Inthisregard,phase-fieldapproacheshaveshowngreatpotential with their abilityto automatically determinethe directionofcrackpropagationthrough minimization ofan energyfunctional.

Thephase-fieldframeworknaturallyhandlestheemergenceofphenomenasuchascracknucleationandbranchingwithout the need to introduceadditional criteria. In particular, formulationsderived fromthe variationaltheory of Francfort and Marigo (1998)have received a lotof attentionfrom the applied mechanics communitydue tothe latter’s strong tiesto Griffith’stheory forbrittle fracture.Phase-fieldmodelsbelong tothecategoryofcontinuum approachesforfractureprop-

∗ Corresponding author.

E-mail addresses: [email protected] (J.M. Sargado), [email protected] (E. Keilegavlen), [email protected] (I. Berre), [email protected] (J.M.

Nordbotten).

https://doi.org/10.1016/j.jmps.2017.10.015

0022-5096/© 2017 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license.

( http://creativecommons.org/licenses/by/4.0/ )

agation,utilizingadiffuserepresentationofcracksinplaceofactualdiscontinuities.Theamountofcrackregularizationis controlledviaaprescribedlengthscalethatconstitutesanadditionalparameterofthemodel.

The aimof the presentwork is to addresstwo longstanding issues associatedwith the phase-fieldformulationthat ariseinconjunctionwithuseofthenow-classicalquadraticdegradationfunction.Thefirsthastodowithprematurestiff- ness degradation stemming from theevolution of damage around regions ofstress concentration. The second and more serious issuedealswiththe observed dependenceofsimulated failure loadson the phase-fieldregularizationparameter, a phenomenon that has largely goneunexplored inthe literature until very recently. The problemis mostnoticeable in casesofcrack growthemanatingfroman explicitly meshed initialfracture andundermines the usefulnessofphase-field approachesin solving problems that involve crack initiationat a priori unknown locations (Klinsmann et al., 2015). The parameterwasinitially introducedby Bourdinetal.(2000) asa purelymathematicalconstructthat allowsfortheGrif- fithenergycorresponding toadiscretecracktoberecoveredinthelimitasgoestozero,inthesenseof-convergence (Braides,2006).Ofthetwoaforementionedissues,thefirstmayberemediedbymakinguseofalternativeformulationsas discussedin Pham etal.(2011).On theother hand,the dependency ofmechanicalresponse on isnot yet fullyunder- stood.Ithasbeensuggestedrecentlythatthelattershouldbeviewedasamaterialparameterthatiscloselyconnectedto thecracknucleationstress(Bourdinetal.,2014;Mesgarnejadetal.,2015;Nguyenetal.,2016).Whileweconsiderthetwo pointsraisedaboveasdistinctissues,wenonethelessrecognizethattheyarealsocloselyinter-related,inparticularbecause thefirstoftenexacerbatesthesecond.

Ourcontributioninthefollowingworkistwofold.First,weprovideaconceptualexplanationofhowthechoiceoflength scalecanresulttoeitherdelayoraccelerationoffailureunderquasi-staticconditions.Secondly,weintroduceanewfamily ofdegradationfunctionsthatallows forcorrectlyreproducingtheonsetoffailureforreasonablychosenarbitraryvaluesof theregularizationparameter.Thelatterpointisimportantsincetheproblemofregularization-dependentmaterialresponse isnot solvedinthealternativeformulationspreviouslymentioned,andfurthermoreisnotonlyconfinedtobrittlefracture asdemonstratedin thenumerical resultsofAreias etal.(2016) on crackingin elastoplastic materials.Enthusiasm in the relatively new phase-fieldparadigm has led to its application in a number ofdiverse problems, which include cracking in piezoelectric solids (Miehe et al., 2010b), fluid-driven fracture propagation(Miehe etal., 2015b;Mikeli´c et al., 2015), thermalshock-inducedcracks (Bourdinetal., 2014) andfragmentationofbatteryelectrode particles(Mieheetal., 2015a).

This underscoresthe need forquantitative accuracy withregard to the fracture model, particularly in the caseof crack nucleationthatisoftenthecriticalfailuremechanismformanysuchapplications.

Theremainderofthispaperisstructuredasfollows:WebegininSection2withadiscussionofimportantfundamental conceptsunderlyingthephase-fieldapproachaswellasourmotivationforpursuingthecurrentresearchdirection.Specifics regardingtheformulationusedinthepresentworkare givenin Section3whichalsoincludes importantdetailswithre- gardtonumerics.Inparticularforthecasewherethephase-fieldevolutionequationisnonlinear,weoutlinealinearization schemebasedon atruncatedTaylorseriesapproximation thatavoids theimplementationofnestedloopsinthe solution scheme.Thefollowingtwosectionscontainthemainnoveltiesofthecurrentwork:Section4beginswithanumericalex- ampledemonstratingtheapparent contradictionthat isoftenseenbetweensimulationresultsandwhatisexpectedfrom the-convergencepropertyofthefracturephase-fieldmodel.Thisisfollowedbyadiscussionwhichaimstoexplainwhy thelatteraloneisnotsufficienttoensuretheproperbehaviorofthemodel.InSection5,weproposeanewparametricfam- ilyofdegradationfunctionsthataimstoincreasetheaccuracyofphase-fieldsimulationsbyaddressingkeyissuesdiscussed intheprevioussection.Superiorityoftheresultingformulationovertheclassicalmodelemployingquadraticdegradationis demonstratedviaseveralnumericalexamplesinSection6.Finally,concludingremarksandoutlookaregiveninSection7. 2. Theoreticalaspectsofphase-fieldmodelling

Thephase-fieldframework wasfirstintroduced byFix(1983)andLanger(1986)formodelingphase transitionsinma- terials, andlater extended to free discontinuity problems by Ambrosio andTortorelli (1990) who worked onimage seg- mentation. Its specific applicationto crack propagation in solidsis much more recent, andis the result of independent work byresearchers comingfrom thefields ofphysics (Aranson et al., 2000; Karma etal., 2001) andapplied mechanics (Bourdinetal., 2000).Weadapt thelatterperspective inthisstudy,andfurthermore notethat whilethe originalformu- lationintroducedby Bourdinetal.hasremainedvirtuallyunchangedincurrentusage,theargumenton whatconstitutes propersolutionstotheresultingequationsaswellasthemeaningofkeyquantitiesisfarfromresolved.Inviewofthis,we beginwithashortreviewoftheorypertainingtothephase-fieldformulationforbrittlefracturealongwithadiscussionof significantdevelopmentsinthefieldinordertoprovidecontextforthepresentwork.

2.1. Brittlefracture:fromGriffithtoFrancfort–Marigo

Griffith(1921)canbecreditedasbeingthefirsttoformallystatethethermodynamicprinciplesgoverningthepropaga- tionoffracturesinbrittlematerialsthathasbecomethefoundationofmodernlinearelasticfracturemechanics.According toGriffith’stheory,anexistingcrackwillpropagatewhentherateofenergyreleaseGassociatedwithcrackextensionex- ceedsacriticalvalueequaltothematerialfracturetoughness,Gc.ThiscanbeexpressedviathefollowingsetofKuhn-Tucker conditions(NegriandOrtner,2008):

G−Gc≤0 (1a)

˙

a≥0 (1b)

(

^G−Gc

)

^{a˙=0 (1c)}

wherea˙denotestherateofcracklengthincrease.ThefirstinequalityprecludesthecaseofunstablecrackingwhereG>Gc, whilethesecond isanirreversibilityconstraintthatpreventsunphysicalhealingoffractures.Finally,condition(1c) implies thatGmustbeequaltoGcwhenthecrackisgrowing,andconverselythatacrackcannotextendwhenG<Gc.Animportant shortcomingofGriffith’stheoryisitsinabilitytoaccommodatecracknucleationorpredictthebranchingoffractures.Inan efforttoovercomethislimitation,anextensionoftheframeworkwasdevelopedbyFrancfortandMarigo(1998)intheform ofavariationaltheoryoffracture,whichisbasedonanenergyminimization paradigm.Itstipulatesthatthetotalpotential energycorresponding to alinear elasticbody ^{containinga setof crackpoints can bewritten asa sumofbulk and} surfaceterms:

(

^u,

)

⁼

\ 1

2

(

^u

)

^:C^e:

(

^u

)

^d

+GcHⁿ⁻¹

( )

⁽²⁾

where

(^u)=¹₂[

∇

u+u

∇

^{] isthe symmetric} small-strain tensor,C^e is thestandard linear isotropic elasticitytensor, andHⁿ⁻¹ isthe(ⁿ−1)-dimensionalHausdorff measuregivingthesurfaceareaassociatedwith^{.Eq.(2)isreferred toas} theGriffithfunctional,anditisassumedthat forsomegivenboundaryconditionson^,theunknowndisplacementsu as wellasthecrackset^{canbeobtainedviaaglobal}minimizationofsaidfunctionalsubjecttotheirreversibilitycondition

t+t⊇

t (3)

thatiscomparableto(1b).Incontrast,Griffithrequiresonlystationarityof(2).Furthermoreinthevariationaltheory,^is notrestrictedtoconsistofasinglecrack,andhenceabodythatisinitiallywithoutflaw(=∅)maynucleateacrackifthe resultingconfigurationhaslowertotalenergycomparedtoonewherenocrackforms.Similarly,acrackisallowedtobranch ifthisleadstoa lower potentialenergythan simpleextension. Thedirectionsof advancearenaturallyobtainedasthose leadingtominimumincreasein(2),sothatintheorytheenergyminimizationframeworkisabletohandlecrackinitiation andbranchingwithouttheneedtointroduceadditionalcriteria.

2.2. Phase-fieldandgradientdamagemodels

The main difficulty in performing a direct minimization of the Griffith functional is that u is generally discontinu- ous across ^{, so that (2) contains a locus of jump sets whose locations are a priori unknown andwhich generallydo} not align with the predefined domain discretization used for numerical calculations. To render the problem tractable, Bourdinetal.(2000)adoptedastrategy whereintheminimization isinsteadperformedonanapproximation oftheGrif- fithfunctionalhavingregularizedjumpsetssothatuiscontinuousovertheentiredomain.Thiswasinspiredbytheearlier workofAmbrosioandTortorelli(1990);1992)whosolvedasimilarprobleminimagesegmentationinvolvingthefunctional ofMumfordandShah(1989).Ascalarorderparameterknownasthecrackphase-fieldisintroducedtointerpolatebetween fracturedandintactregions,withacharacteristiclength parametercontrollingtheamountofregularization.Thevalidity ofsuchasstrategyrestsonwhethertheregularizedapproximationtends towardtheoriginal functionalasgoestozero, inthe sense of-convergence(Braides,2006).Althoughan additionalequation governingthe phase-fieldevolutionmust nowbesolvedalongwiththelinearmomentumequation,themainadvantageofthisapproachisthatnumericalsolutions maybeobtainedviaclassicalfiniteelementalgorithmsasbothuandthephase-fieldarecontinuous.Bourdinetal.(2000)’s regularizationof(2)hastheform

(

^u,

φ )

⁼

1 2

(

¹⁻

φ )

²⁺

κ

ε (

^u

)

^:Ce:

ε (

^u

)

^d

+Gc

₁

2

φ

²⁺₂

∇φ

²

d

, (4)

where

φ

^{is the phase-fieldthat takes on values between0 and1,} corresponding respectively to fully intactand broken states.Ontheother hand,

κ

^{isasmallpositive constantmeant toensurepositivityofthebulk energyas}

φ

→1.Thetwo important features of the above expression are (a) the replacement of Hⁿ⁻¹() ^{by an elliptic functional that calculates} the combinedlength ofall thecracks, and(b)thecoefficient (¹−

φ

)²+

κ

^{knownasthe energy}degradation functionthat penalizesthematerialstiffnessaccordingtothevalueof

φ

^.Eq.4isessentiallyadirectadaptationoftheearlierfunctional ofAmbrosioandTortorelli(1992),forwhichaproofof-convergencewassubsequentlygivenbyChambolle(2004).

2.2.1. Alternativevariationalproblems

It waslater suggested by Miehe etal.(2010c) that developmentofelliptic approximationsto Hⁿ⁻¹() ^{could also be} motivatedfromamorephysicalstandpointbyconsideringthe1-dimensionalcaseofaninfinitelylongbarofuniformcross section.Assumingthatthebarisalignedwiththex-axisandthatasinglecrackfullycutsthebaratx=a,thephase-field profilecorrespondingtothissharpcrackisnoneotherthanthediscontinuousscalarfunction

φ (

^x

)

⁼

1, x=a

0, otherwise. (5)

Fig. 1. Phase-field regularizations of a sharp crack at x = a using various candidate functions.

Aregularizedapproximationoftheabovecanthenbemadeviaafunction

φ

∈ ,where

=

φ

^∈H1

(

^R;[0,1]

)

φ (

^a

)

⁼¹

φ (

^x

)

^ismonotonicallydecreasingawayfroma

φ (

^x

)

→0asx→±∞

. (6)

Somecandidatefunctionsare

φ (

^x

)

=exp

−

|

^x−a

|

(7)

φ (

^x

)

⁼

⎧ ⎨

⎩

1−

|

^x_√^−a

|

2

2

, x∈

a−√

2,a+√ 2

0 otherwise

(8)

φ (

^x

)

⁼

1+

|

^x−a

|

exp

−

|

^x−a

|

. (9)

Wenotethatthefunction(7)utilizedbyMieheetal.(2010c)isinfactthesolutionobtainedbyminimizingthefunctionalof Bourdinetal.(2000)inonedimension.Ontheotherhand,(8)isacompactlysupportedfunctionthatleadstoavariational inequalityproblem(Phametal.,2011),while(9)isobtainedbysolvinga4thordergoverningequationandwasintroduced byBorden etal.(2012) withtheaimofloosening themesh sizerequirementswithrespect toatthesame timetaking advantageofnumericalmethodsthatcanadequatelymodelC¹-continuoussolutions.Fig.1showsacomparisonofthethree functionsmentioned above. It can be observed that fora given value of , the amountof crack diffusionis additionally dependentonthespecificformof

φ

^(x).

Higher order phase-fieldformulations forfracture haveso far not achieved the samepopularity astheir lower order counterparts.Thisismainlyduetothehigherorderofcontinuitythattheyrequireinconjunctionwithnumericalsolutions, whichisexpensivetoobtain withtraditionalframeworkssuchasfiniteelements. Inaddition,-convergenceisyettobe proven fortheseformulations.Onthe other hand,Lietal.(2015) point out that theincorporation ofgeneralanisotropic effectsrelatedtothesurfaceenergyrequiresaformulationthatisatleast4thorder.

2.2.2. Damage

Theconnectionbetweenphase-fieldapproachesandnonlocaldamagemodelswasexplored byPhametal.(2011),who notedthat elliptic functionals approximating (2)can be seen asspecific cases of theintegral ofa generalstate function pertainingtoagradientdamagemodel:

W

( (

^u

)

^,

φ

^,

∇φ )

⁼¹₂

^:C

( φ )

^:

^+w

( φ )

⁺¹₂^w12

∇ φ

^·

∇ φ

⁽¹⁰⁾

wherew(

φ

^)isamonotonicallyincreasingfunctionintheinterval[0,1]withw(⁰)=0andw(¹)=w₁.Theysuggestusing the linearform w(

φ

)=w₁

φ

^{which leadsto models having a real elastic phase withno premature decrease in material}

stiffness,atthe costof solvinga variationalinequality problemforthedamage evolution.Incontrast,the quadraticform of w(

φ

^{) employed in(4) leads to material behavior having no real elastic phase, with damage already occurringat the}

onsetofloading(Amoretal.,2009).Bulkenergyreleaseresultingfromevolutionofthephase-fieldisfacilitatedthrougha

damage-dependent elasticitytensorC(

φ

) ^{asseenin(10).The simplestformwhich leadsto isotropicbehavior consistsof} themultiplicativeansatz

C

( φ )

^=g

( φ )

^{Ce, (11)}

inwhichg(

φ

^)istheenergydegradationfunctionmentionedpreviouslythatisnon-negativeintheinterval[0,1]withg(⁰)= 1 and g(¹)=g(¹)=0. Such form howeverhas limited applicability, since it allows for unphysical compressive cracking based on Gc.Lancioni andRoyer-Carfagni (2009) adoptedthe above model to shearcracking by having g(

φ

^{) act onlyon}

thedeviatoricportionofthestrain. Thiswassubsequentlyimprovedupon byAmoretal.(2009)whoproposedthatcrack growth be driven also by the spherical part of the energywhen the volumetric strain is positive. This results in more realisticanisotropicbehaviorwherethematerialisallowedtocrackinvolumetricexpansionandshear,butnotinvolumetric compression.AnalternativeformulationwasintroducedbyMieheetal.(2010a;2010c)inwhichdegradationoccursonlyon tensilecomponentsoftheprincipalstraintensor,leadingtopuremode-Icracks.

2.3. GoingbacktoGriffith

The functional ^(u,

φ

^{) in (4) is neither linearnor convex, which makes the taskof finding globalminimizers non-}

trivial. Bourdin et al. (2000) proposed an alternate minimization algorithm that takes advantage of the convexity of withrespecttoeitheruor

φ

^{whentheotherisheldconstant.}Nonetheless,solutionschemesbasedondescentalgorithms only convergeto localminimizersor saddlepointsandare by themselves inadequatefor obtaining globalminimizers. It wasfoundthatnaiveapplicationofthealternateminimizationoftenresultedinsolutionsthatexhibitedunphysicaldipsin thetotalenergy,particularlywhencrackgrowthisbrutal.Inan attempttoremedythisbehavior,a heuristicbacktracking schemewasdevelopedbyBourdinetal.(2008)(withsubsequentimprovementsbyMesgarnejadetal.(2015))basedonan additionaloptimalityconditionthatenforcesmonotonic evolutionofthetotalenergywhen theloadisalsomonotonically increasing.Fromaphenomenologicalstandpointhowever,themainobjectiontousingglobalenergyminimizationisthatit allowsforevolutionswherethecurrentconfigurationjumpsoverarbitrarilylargeenergybarriersinordertoreachthenew configurationcorrespondingtotheglobalminimizer(NegriandOrtner,2008).Whilethisenablesthestrict preservationof energyconservation,itmayalsoresultinunphysical responsewherecracks propagateatlower energyreleaseratesthan Gc whichviolatesGriffith’scriterion.Recently,Larsen(2010)introducedthenotionof

^{-stabilityasasteppingstonetoward}

formulationsthatcanpredictcrackpathsbasedonlocalminimality,whichisinturnclosertoGriffith’soriginalidea.Aswith Griffith’smodel,suchsolutionswillalsoexhibitdissipationinthetotalenergyincaseswherecrackpropogationoccursin abrutalmanner.However,aspointedout byNegriandOrtner(2008),brutalcrackingisprimarilyadynamicphenomenon whichexplainswhythetotalenergycannotbecompletelyaccountedforinaquasi-staticframework.

2.4. Onthetreatmentof

If we settle forinvoking localversus global minimality, then classical solutionschemes that were previously deemed inadequatearenowrobustwithouttheneedtoperformbacktracking.Weareleftwiththenon-convexityofGriffith’sfunc- tional,butthisiseasily dealtwithvia alternate minimizationasearliermentioned. Ontheother hand,we endup again with Griffith’soriginal conundrum concerning crack initiation.It turns out that the saving grace is none other than the regularization ofthefunctional,whichaswill bediscussed inthelater sectionsactuallyallows forcrackinitiationinthe absenceofstress singularitiesprovidedthatthecharacteristiclengthassociatedwiththeregularizationisfinite.Thisalso brings usmore inlinewithnonlocaldamage theory,inwhich thethicknessofthe localizationzone isusually aconstant parameterassociatedwithsomephysicalinternallength.

The notion of an intrinsic length scale relating Gc to the material strength

σ

c was originally introduced by Irwin (1958) based on the assumption that a plastic zone develops in front of a propagating crack. An estimate of the lengthofthiszoneaheadofthecracktipisgivenby

lc=EGc

σ

c²

. (12)

However,BažantandPijaudier-Cabot(1989)pointout thattheabovequantityisdistinctfromthecharacteristiclengthof thenonlocalcontinuum;thelatterisdefinedas

=

α

^EG

σ

c^2c

(13) anddescribesthewidthofthesofteningzoneanalogoustotheminimumspacingofcracksinadiscrete-crackmodel,with

α

^{beingaconstantthatdependsontheparticularnonlocal}formulation.Thecrucialdifferenceisthatwhilelccanberightly viewedasa materialproperty, isincontrasta model parameterdueto

α

^.Nevertheless withinthe contextofa specific model,the valueof

α

^isalreadypredeterminedso thatagainacts asamaterial parameterdependent onE,Gc and

σ

c. Ontheotherhandithasbeenobserved(e.g.Klinsmannetal.,2015) thatforcracksmodeledasexplicitboundariesinthe mesh,numericallypredictedcriticalenergyreleaseratesexhibitadependenceonwhenthelatterisnotsufficientlysmall.

Itcanthushappenforcertainsituationsthat

σ

^candGcimposeconflictingrequirementson.Furthermoretheassumption

ofmaterial-dependent regularizationparameters meansthat inheterogeneous media,the characteristiclength takesona differentvalueforeachmediumsothatthediffusecrackbecomesthickerorthinnerasitpassesfromonematerialtothe next.

Letusnowconsideranonlocal formulationthat isdependentonan additionalsetofparameters representedbysome vector

η

^sothat

α

=

α

(

η

)^{.Ifweassumethatthefunctionis}invertible,thenwecanrewrite(13)as

η

⁼

α

⁻¹

σ

c²

EGc

(14)

implyingthatcannowbechosenindependentlyofE,G^cand

σ

^{c,providedthat}

η

^{exists.From(10)and(11)wecanseethat}

themoststraightforwardwaytointroducesuchparametersintothemodelisthroughthedegradationfunction,anditisin factthisrealizationthathasmotivatedthepresentwork.

3. Governingequationsandnumericalimplementation

Forthe remainderofthisstudy,wehavechosen toadoptthefunctionalofBourdinetal.(2000)by reasonofitssim- plicity.Incorporatingtheworkdonebyexternalforces,theregularizedtotalpotentialenergyforagivenbody^subjectto boundaryconditionsisgivenby

=

−W=

g

( φ ) ψ ( )

^+Gc

₁

2

φ

²⁺₂

∇ φ

^·

∇ φ

d

−

b·ud

−

∂^tt·udS (15) where

ψ

(

)=¹2

^:C^e:

^{istheHelmholtzfreeenergydensityand}

∂

^{tdenotesthepartoftheboundaryforwhichNeumann} (i.e.traction)conditionsareprescribed.Thequantitybrepresentsthebodyforce,whiletisthevectorofprescribedtractions actingontheNeumannboundary.Byimposingthestationarityof^weobtainthevariationalequation

δ

=

g

( φ ) ∂ψ

∂

^:

δ

^d

−

b·

δ

^ud

−

∂^tt·

δ

^udS

+

g

( φ ) ψ ( ) δφ

^d

+Gc

₁

φ δφ

+

∇φ

·

δ∇φ

d

=0. (16)

Theaboveequalitymustholdforarbitraryvaluesof

δ

^uand

δφ

^{,implyingthat}

g

( φ ) ∂ψ

∂

^:

δ

^d

=

b·

δ

^ud

+

∂^tt·

δ

^{udS (17a)}

g

( φ ) ψ ( ) δφ

^d

^+Gc

₁

φδφ

⁺

∇φ

^·

∇δφ

d

^{=0. (17b)}

Theseconstitutetheweak formofthegoverningequations.Notingthat

σ

=

∂ ψ

/

∂

,theequivalentstrongformulation maybeobtainedbyapplyingGauss’divergencetheoremyieldingthefollowingcoupledsystem:

∇

·[g

( φ ) σ

^]+b=0on

^(18a)

g

( φ ) σ

^{·n=t on}

∂

^{t (18b)}

u=u¯ on

∂

^{u (18c)}

Gc

∇

²

φ

−Gc

φ

=g

( φ ) ψ ( ε )

^on

^(18d)

∇φ

·n=0on

∂

. (18e)

Eq.(18a)–(18c) comprisethelinearmomentumequation andits correspondingboundaryconditions,while(18d)isthe phase-fieldevolutionequation withthe associatedboundary conditiongivenby (18e).As mentionedearlier,

φ

^{should go}

tozeroawayfromthecrack.Thus itistacitlyassumedthatthedomainisofsufficientsizetoprovideadequateseparation betweenthe regularized crack and the boundary, allowing the phase-field to decay to values that are small enough to approximatethiscondition.

Extensionoftheirreversibilitycondition(3)totheregularized caseisnot immediatelyobvious,sincetheintermediate states0<

φ

<1 donot haveastraightforward physicalinterpretation.The naturalcoursefromthe perspectiveof damage mechanicsistoenforcethecondition

φ (

^x

)

t+t≥

φ (

^x

)

t

∀

^x∈

. (19)

Thiscanbeimposedvia anadditionalpenaltyterminthephase-fieldevolutionequation(Mieheetal.,2010c),oralterna- tivelythroughahistoryvariableHthatreplacesthequantity

ψ

⁽

^{)in(18d),definedas(Mieheetal.,2010a)}

H

(

^x,t

)

⁼_s^max

∈[0,t]

ψ ( (

^x,s

) )

⁽²⁰⁾

A closerlookat thephase-field localizationprocess howeverreveals that (19)maynot be the best extension of(3). For exampleinthe1-dimensionalcase,Kuhnetal.(2015)demonstratethatdamagelocalizationcorrespondingtoadiffusecrack involvesnotonlythegrowthof

φ

^{nearthecracktipbutalsoadecreaseinadjacentregionswhichenablesthephase-field}

tocorrectlysettletotheexponential profilegivenby (7).Based onthisobservedbehavior, astrictimpositionof(19)may leadtoanoverestimateofthecracklength.Consequently,onecanuseamodifiedversionof(20)inwhichirreversibilityis imposedonlywhen

φ

^{exceedsacertainthreshold,i.e.}

H

(

^x,t

)

=

max

s∈[0,t]

ψ ( (

^x,s

) )

^if

φ

>

φ

^c

ψ ( (

^x,t

) )

^{otherwise. (21)}

Theparameter

φ

^{crepresentsthemaximumvalueofdamagethatisallowedtohealduringunloading.Foramaterialpoint}

undergoingdamage

φ

≤

φ

^{c,theresulting}stress-straincurveswillbenonlinear,withtheamountofdeparturefromlinearity dependentonthespecificformofthedegradationfunction.Nonethelesshaving

φ

^c>0allowsthestresspathsforsubsequent unloading andreloading to coincide withthe initialloading curve, which cannot be achievedotherwise. Forthe present work,wehavechosentoset

φ

^c=0.5inallofoursimulations.

Thecoupledsystemdescribedin(18)isimplementedinaclassicalfiniteelementframework,withtheprimaryunknowns beingthedisplacementuandphase-field

φ

^.Ina2-dimensionalsetting,theseareexpressedintermsofthecorresponding nodaldegreesoffreedomas

u= m

I=1

N_I^uuI and

φ

= m

I=1

NI

φ

I (22)

wherein NI=

NI 0 0 NI

(23)

withN_I=N_I(^x)^{denotingtheshapefunctionassociatedwithnodeI,andu}Iand

φ

I therespectivedisplacementandphase- fielddegreesoffreedom.Thestrainandphase-fieldgradientaregivenby

^=m

I=1

B^u_IuI and

∇φ

^=m

I=1

B^φ_I

φ

^{I (24)}

inwhich B^u_I =

_N

I,x 0 0 NI,y

N_I,y N_I,x

and B^φ_I =

NI,x

N_I,y

. (25)

The former is the symmetrized gradient matrix associated with the Voigt form of the strain tensor. The test functions andtheir correspondingderivativescan be obtainedfromthe aboveexpressions byreplacing u_I and

φ

Iwith

δ

^uIand

δφ

I

respectively.Notingthatthelattertwoquantitiesmustbearbitrary,numericalapproximationoftheweakformin(17)yields thefollowingnonlinearsystemofequationsateachnode:

r^u_I =

g

( φ )

^BuI^T

σ

^d

−

N_I^uTbd

−

∂N_I^uTtdS=0 (26a)

r^φ_I =

GcB^φ_I^T

∇φ

^+GcNI

φ

d

⁺

NIg

( φ )

^Hd

^{=0. (26b)}

Duetothenon-convexityof(15),wesolvetheabovecoupledsystembymeansofthealternateminimizationalgorithm outlined inBourdinetal.(2000).Thisinvolvescyclingbetween(26a)and(26b):thelinearmomentum equationissolved firstusingvaluesof

φ

^{fromtheprevious iteration,afterwhichtheupdatedvaluesofuareused todeterminethecurrent}

history field according to (21). The phase-field evolution equation is then solved utilizing the updated valuesof H. The iterationsarecarriedoutrepeatedly untiltheprescribedcriteriaonthesizeofresidualsandinter-iterationcorrectionson theunknownsaremet.

Fordegradationfunctionsinwhichthederivativeg(

φ

^{)isnonlinear,thesubsystem}representedby(26b)isalsononlin- ear resultingin theneed toperform nestediterations:an outer loop that cyclesbetweensubsystems, andan inner loop implementingtheNewton-Raphsonmethodforthenonlinearsubsystem.Incontrast,anaive implementationofthe alter- nate minimization algorithm whereinthe linearizedphase-field equation issolved only once beforegoing back tolinear momentumisgenerallyunstableandleadstoincorrectresults.Thatis,inthesubsystem

{ φ}

^mi⁺¹=

{ φ}

^mi −

K^φφ

φ

i^m

−1

r^φ

m

i (27)

correspondingtothe(^m+1)^{thiterationwithintheithtimestep,useoftheexactJacobiangivenby} K_IJ^φφ

φ

i^m

=

GcB^φ_I^TB^φ_J +

_Gc

+Hg

φ

i^m

NINJ

d

⁽²⁸⁾

oftenproducesan incorrectevolution ofthephase-fieldandeventual blow-up.However, wehavefound thatsuch anap- proachcanbemadetoworkbyreplacingg

φ

_i^m

withanapproximateexpressionderivedfromlowordertermsinaTaylor expansion.Assumingthatthedegradationfunctionissmooth,g(

φ

1)canbe obtainedasaninfinitesumoftermsinvolving higherorderderivativesofgevaluatedat

φ

2∈[0,1]:

g

( φ

1

)

=g

( φ

2

)

+g

( φ

2

) ( φ

1−

φ

2

)

+g

( φ

2

)

2

( φ

1−

φ

2

)

²+... (29)

Nowlet

φ

1=1sothatg(

φ

1)=0.Droppinghigherordertermsaswellasthesubscripton

φ

2,weobtain

0≈g

( φ )

^+g

( φ ) (

¹⁻

φ )

⁽³⁰⁾

whichthengivesusourapproximationforthe2ndderivativeofg: g_app

( φ )

^=−g

( φ )

1−

φ

^{. (31)}

Theabovetechniqueallowsforthestraightforwardimplementationofthenumericalmodelwithinmorerigidlystructured softwareframeworks (e.g.commercialprograms) that mightnot havethe capabilityto implementa full Newton scheme separatelyoneachsubsystem.

4. Acuriouscaseofcracknucleation

Pre-existingcracksmaybeaccountedforintwoways:(a)throughinitializationofthephase-fielditself,and(b)bymod- elingthecrackfacesdirectlyasinternalboundariesinthediscretizedgeometry.SicsicandMarigo(2013)provideanalytical resultsshowingthatthe growthoffullydevelopedfracturesdescribedviathephase-fieldobeysGriffith’slawasthereg- ularizationparametergoesto zero.Ontheother handit hasbeenshowninnumerical experimentsthat thesameis not generallytruefortheextensionofcracksbuiltintothemesh.Recentstudies(e.gKlinsmannetal.,2015;Nguyenetal.,2016) haveobservedthat thesimulatedcritical energyreleaserate(oranalogously,the peakload) overshootsthecorrectvalue forsufficiently small.This inconsistencybecomesmore understandableupon the realizationthat propagationof mesh- describedcracks hasmore to dowithnucleation rather than extension inthe context ofphase-fieldapproaches.Similar behaviorcan be observedinthe caseofcracknucleation atanotch orreentrantcorner;the extensionof mesh-modeled crackscanbeseenasalimitingcaseoftheformerwhenthenotchanglegoestozero.Onecanarguethattherearethree distincttypesofsimulatedmaterialresponseinconnectionwiththephase-fieldmodelforbrittlefracture:(a)propagation ofphase-field-describedcracks thatisrelativelywell-understood,(b)cracknucleationintheabsenceofstress singularities whichwillbe discussed inSection 5.4,and(c)cracknucleation inthe presenceofsingularities.The last categorycan be understoodtoincludebothweakstress singularitiesthatariseinconnection withV-notches, andstrongsingularitiesthat areassociatedwithsharpcracks.

4.1. Preliminarynumericalexample

Toillustratethedependenceofthematerialresponseonthephase-fieldlengthscale,wesimulatefracturepropagation inahomogeneous specimencontaininga centercrackandsubjectedtotensileloadingasshowninFig.2a.Assumingthat Histakenlargeenoughsuch thatthetensilestressesattheboundaryareacceptablyuniform,themode-Istressintensity factorcanbecomputedforfinitevaluesoftheratioa/bas

K_I=

σ

^√

π

^aF

(

^a/b

)

⁽³²⁾

whereF(a/b)isashapefactorgivenby F

(

^a/b

)

⁼

1−0.025

_a

b

2

+0.06

_a

b

4

sec

π

^a

2b. (33)

Theabove formulahasa reportedaccuracyof0.1%orbetterforanya/b(Tada etal., 2000).Wenote thatthe casewhere a/b=0andH=∞correspondstotheoriginalfractureproblemofGriffith(1921),forwhichF(⁰)=1.Fortheplane strain case,thecriticalstressintensityfactorandstrainenergyreleaseratearerelatedby

Gc=K_Ic²

1−

ν

²

E

. (34)

Combiningtheabovewith(32),weobtainthefollowingexpressionforthecriticalfailureload:

Pc= b F

(

^a/b

)

EGc

1−

ν

²

π

^a

1/2

. (35)

Fig. 2. Problem of a center-cracked specimen subjected to uniform tension at the far-field, showing (a) the full specimen geometry, (b) computational domain and boundary conditions, and (c) load-displacement curves for different values of and mesh refinement.

TheactualcomputationaldomainisshowninFig.2balongwiththerelevantboundaryconditions;duetosymmetry,only half the geometry needs to be considered. The initial crack ^{of length a in the} computational domain is modeled as partofthe geometry,andno initializationofthe phase-fieldisperformedto accountforits presence.Actual dimensions usedarea=10mm,b=2aandH=10a.Thusa/b=0.5leadingtoF(^a/b)=1.1862.ThematerialconstantsareE=70,000 MPa,

ν

=0.22andGc=0.007N/mm.Fromtheprecedingequation,weobtain thecriticalfailure loadasP_c=68.26N.The simulation is carried out using monotonic displacement control withthe specimen gradually stretched inincrements of ^U=2.5×10⁻⁴ mmuntil failureoccursintheformofbrutal cracking.Inordertoobtain aprecisedeterminationofthe failure load,^{Uisreducedto2}.5×10⁻⁵mmasfailureisapproached.Fig.2cshowsthedependenceofsimulationresults onthephase-fieldcharacteristiclengthaswell astherelativemeshrefinement,/h^e.Itcanbeseenthat largervaluesof leadtoanincreaseindeviationfromlinearbehavior,whereasasmallerdrivesthepeakloadupwards.Moreimportantly, theresultsshow thatthepeakloaddoesnotconvergetotheanalyticalvalueastheregularizationparameterisdecreased, butratherovershootsthetruesolutionforasisreducedpastacertainthreshold.Ascanbeobserved,theseverityofthis phenomenon isalso influencedby themesh refinement,andinparticular isgreater forcoarsermeshesrelative to.We havefoundthatdifferentsetsofmaterialparametersgivequalitativelythesamebehavioraswhatwehaveshown.

Decreasingthevalueevenfurther(uptotwoordersofmagnitude)fortheaboveproblemdoesnotappeartoimprove accuracywithrespecttoPc.AsshowninFig.3,theerrorinthesimulatedfailureloadremains above10%forthesmallest value of(=0.0025mm)used.However, theneed toproperlyresolve thereducedlength scale inthemesh meansthat thecomputationalsizeoftheproblemincreasessubstantiallyforverysmallvaluesof.Wehavefoundthatforthequasi- static caseasemployed inthe presentwork, the average crack advance betweeniterationswithin thecritical time step hasmore orlessa linearrelationwith.Inthe casewhere=0.0025mm forexample,thecracktip movesan average distanceofaround 0.0003mmperiterationsothat onewouldneedaround theorderof33,000iterationsforthecrackto fullyextendthroughthespecimen.Usingameshrefinementratioof/h=8,theresultingdiscretizationcontainedatotal numberof ≈5.97millionDOFs,andasingleiterationofthealternateminimizationalgorithmtookaround35stocomplete onaworkstationequippedwitha quad-coreprocessorrunningat2.8GHzwithhyper-threading,andusingasparsedirect solverwithbothequation assemblyandsolutionfullyparallelized.¹Thisimpliesarunningtime ofmorethan 13daysfor thecriticaltimestep.TosavetimeingeneratingthedataforFig.3,weterminated thesimulationswhenitwasclearfrom inter-iterationoutputthatthecrackhadstartedtopropagate.

4.2. Exploringtheovershootphenomenon

The apparent lackof convergence in thematerial response withrespect to inthe above numerical example seems tocontradictthe-convergencepropertyof(4),howeverthiscanbeexplainedbythefactthatGriffith’scriteriondoesnot

1We exploit the fact that sparsity profiles for the global tangent matrices do not change in the absence of adaptive mesh refinement, so that the symbolic factorization step can be done once in the beginning and then skipped in subsequent iterations. Without such optimization, the running time for each iteration increases to around 70 s.

Fig. 3. Dependence of peak load on for the center-cracked specimen, obtained using quadratic degradation and a mesh refinement ratio of /h = 8 . For comparison, the semi-analytical peak load is P c= 68 . 26 N. Also shown are the typical speeds of crack advance during the critical time step. Note that the left axis uses a standard linear scale for the peak load, while the bottom and right axes have logarithmic scaling. For the above plots, ∈ {0.0 025, 0.0 05, 0.01, 0.02, 0.04, 0.08, 0.16}.

actuallyinvolvetheenergyfunctionaldirectlybutratheritsgradients(Fréchetderivatives).Thisnuancehasnocorresponding counterpartinimagesegmentation,andmakesphase-fieldsimulationofbrittlefractureafundamentally differentproblem fromtheformer,despite thesimilarityofthe GriffithenergytotheMumford–Shahfunctional.Thus, while-convergence oftheregularized approximationto thesharp-boundaryfunctionalis byitself sufficienttoproduce physicallymeaningful resultsinanimagingcontext,thisisnolongerthecaseforbrittlefracture.

Atpresent,ourunderstandingoftheabovephenomenonreliesonnumericalevidenceobtainedfromanalyzingproblems such asthe one presentedin Section 4.1.To elucidate further, recallthat for a material fracturing accordingto Griffith’s theoryassummarizedin(1),thefollowinginequalityapplieswithregardtoenergyincrements:

−

δ

b^e≤Gc

δ

⁽³⁶⁾

forsomearbitrarysmallcrackextension

δ

>0,where_b^{e denotestheexactelasticbulkenergy}correspondingtothepres- enceofadiscrete crack^.Now −

δ

_b^e=G

δ

^whereGis the energy release rate at the cracktip, so thestrict inequality

−

δ

_b^e<Gc

δ

means that the crackmustbe stationarydueto (1c).If−

δ

^e_b=Gc

δ

,thenapositive

δ

^{isadmissibleand} thecrackcan propagate stably.Onthe other hand,thereverse inequality−

δ

_b^e>Gc

δ

^{isgenerallyunderstoodascorre-} spondingtobrutalcracking.Inaquasi-staticframeworkwheredynamiceffectsaredisregarded,fracturepropagationsimply continuesuntilastateisreachedwhereinthecondition

δ

b^e<Gc

δ

^{isonceagainsatisfied,resultinginarrestofthecrack.}

Itcanbeshownthatinmanycases,suchaconditioncannotbesatisfiedforanylength ofcrackadvancewhichresultsin thefracturecuttingthroughtheentirewidthofthedomain.

Similar behavior ismanifested by the evolution of

φ

^{in the} diffuse-crack modelduringbrutal crack propagation, and canobservedbyscrutinizingsuccessiveiterationswithintherelevanttimestep.Incontrasttotheoriginaltheoryhowever, thephase-fieldmodel containsonlythe equalitypartofGriffith’s criterion inthephase-fieldevolution equation. Thatis, (17b)canbewrittenas

−

δ

_b^app=Gc

δ

⁽³⁷⁾

wherein

δ

_b^app=

g

( φ ) ψ ( ) δφ

^d

⁽³⁸⁾

δ

=

₁

φ δφ

+

∇φ

·

∇δφ

d

⁽³⁹⁾

forsomepositive

δ

^{thatarisesfromanarbitrary}incrementalevolutionofthephase-field,denotedby

δφ

^.Theimplications ofthis are immediatelyobvious when one looks atthe strong form ofthe phase-fieldequation in(18d): assuming that g(

φ

⁾<0everywhereexceptat

φ

=1(andthisisinfactnecessaryfordamagetoevolveatall),thenitisclearthat

φ

^must

beginmovingawayfromitsinitialvalueof0fromthemomentthatnonzerostressisinducedinthematerial.Furthermore,

let

η

bdenotetheerrorarisingfromusing_b^appinplaceof_b^e,i.e.

η

b=

_b^app−

b^e. (40)

Pluggingtheaboveinto(37)andwriting

δ

_b^{e intermsofG,weobtainthefollowingrelation:}

δη

b=

(

^G−Gc

) δ

. (41)

From the previous numericalexample, we caninfer that atsome critical loadingUs brutalpropagation ofthe crackwill occur,presumablybecausenow−_b^app>Gc

δ

^forany

δφ

^{.Theequivalentconditionintermsof}

η

bandGisgivenby

δη

b<

(

^G−Gc

) δ

. (42)

The different curvesin Fig.2c demonstratehow theactual value ofU_s dependson . The key ideais that both

δ

^and

δη

b are influenced by , but in varying degreesfrom one another. In particular, it is no longer justthe quantity G−Gc

thatdeterminestheonsetofbrutalcracking;ascanbeobservedfromFig.2c,bothundershootandovershootofthecorrect failureloadarepossible.Thechallengeistohave(42)occurattheprecisemomentthatGexceedsGc,sothatbrutalfracture occursatthecorrectmagnitudeofloading.

ForthenumericalexampleinSection 4.1,theabovescenario canbe achievedvia aproper selectionoftheregulariza- tionparameter.However,theneedtospecify(andobviously>0)bringsintoquestionthebenefitofhavingregularized approximations^{-converge to theGriffith energyatall as} goesto zero.One can arguethat the removalof such are- quirementisnotadisadvantagesinceitlendsmoreflexibilitytothephase-fieldframeworkandlikewiseopensthedoorto otherinteresting andmoreexoticapproximationsof(2),suchasthehigherorderformulationsbyBordenetal.(2014)and Lietal.(2015)thathavesofarnotbeenproventobe-convergenttoGriffith’senergy.Moreimportantly,themainproblem withrelyingon calibratinginorder toobtainthecorrectinstanceoffailure isthat sucha strategyisnot guaranteedto succeedinall possiblecases,inparticular whenthesetup isvery differentfromtheone analyzedabove.This isdemon- stratedinSection6,wherewestudyaproblemforwhichtheaforementionedtechniquedoesnotworkatall,atleastwithin practicallimitations.

4.3. Preservinglinearityinthematerialresponse

An importantconsequence of(41)is that thematerial response ofthe regularized modelinevitablydrifts fromlinear elasticbehaviorpriortofracture,sincegrowthofGasaresultofincreasingUmustbematchedbyacorrespondingincrease intheincrementalerrorterm

δη

b.Since

η

brepresentsthediscrepancybetweenapproximateandtheexactbulkenergies, an ever-increasing increment in theerror termmeans that the simulated material behavior deviatesfurther andfurther fromlinearelasticitywithincreasing Uasevident inFig.2c.Somecontrolon

δη

bcanbeexercisedthroughthefactor

δ

^, i.e.we keep

δη

bsmallbykeeping

δ

^{smallaswell.Howeversince}

δφ

^{isarbitrary,wecanaccomplishthisonlyby careful}

constructionofeitherthedegradationfunction(whichaffectsthebulkenergy),thecracklengthfunctional,orpossiblyboth.

5. Anewfamilyofdegradationfunctions

TheideaspresentedinSections4.2and4.3canbecombinedtogive asetofpropertiesforwhatwewouldconsideran accuratephase-fieldmodelwithregardtotheextensionofmesh-describedcracks:

(a) Thesimulatedcriticaldisplacementshouldpreferablybeclosetothecorrectvalue,and (b) Theaccumulatederror

η

bshouldbekeptsmallpriortotheoccurrenceofbrutalfracture.

Item(b)isquitestraightforward,andisachievedbyhavingbrutal fractureoccuratlowvaluesofthephase-field.Such behaviorisreadilyobservedwiththealternativestoquadraticdegradationthathaveappearedintheliterature,suchasthe quarticfunction

g4

( φ )

⁼⁴

(

¹⁻

φ )

³⁻³

(

¹⁻

φ )

^{4 (43)}

utilizedbyKarmaetal.(2001)inconjunctionwiththeirownphase-fieldtheory,andthecubicfunction g3

( φ )

^=s

(

¹⁻

φ )

³⁻

(

¹⁻

φ )

²

+3

(

¹⁻

φ )

²⁻²

(

¹⁻

φ )

^{3 (44)}

analyzed byBorden (2012),inwhichthe quantity scontrolstheslope ofthedegradation functionatthe unbrokenstate.

Kuhnetal.(2015) haveshownthat allthree functionshave similarpost-failurebehavior instablecrackgrowth,i.e.their differenceslie primarilyinthepredictionofthelevelofstrainorstressatwhichcrackpropagationoccurs,andalsointhe amountofstiffnessreductionobservedpriortotheonsetofcracking.

Item (a)is moredifficult tosatisfy,inparticular sincequantities pertaining tothe bulkenergy arealso dependenton material properties.The degradation function mustthen be parametric, in orderto havethe means ofcompensating for different valuesof theseproperties.We can seethat noneofthe differentfunctions mentioned above possessthelatter property, so thatone is insteadforcedto relyon tweaking asisdone withthe quadraticdegradation function.From a conceptualstandpoint thisisnotentirelysatisfactory,sinceasa parameterbelongstothecrackfunctionaltermandnot the bulk energy.Furthermore, a changein the regularizationparameter leadsto corresponding changes in boththe bulk