A comparative accuracy and convergence study of eigenerosion and phase-field models of fracture
aa r X i v : . [ m a t h . NA ] J a n A COMPARATIVE ACCURACY AND CONVERGENCESTUDY OF EIGENEROSION AND PHASE-FIELD MODELSOF FRACTURE
A. PANDOLFI , K. WEINBERG AND M. ORTIZ Abstract.
We compare the accuracy, convergence rate and compu-tational cost of eigenerosion (EE) and phase-field (PF) methods. Forpurposes of comparison, we specifically consider the standard test caseof a center-crack panel loaded in biaxial tension and assess the conver-gence of the energy error as the length scale parameter and mesh sizetend to zero simultaneously. The panel is discretized by means of aregular mesh consisting of standard bilinear or Q Introduction
The tracking of crack growth in solids is a free-discontinuity problem in-volving the formation of new internal surfaces. Modelling crack propagationremains a challenging problem in computational mechanics. In recent years,the phase-field (PF) method has gained in popularity, especially for crackpropagation and tracking problems, cf., e. g., [1, 2, 3, 4, 5, 6, 7, 8]. Pro-posed originally by Ambrosio and Tortorelli [9, 10] as a regularization of theMumford-Shah image segmentation functional, the PF method introducesan auxiliary continuous field, the phase field , as a means of representing thestate of the material in vicinity of the crack. In effect, the phase field smoothsthe sharp surfaces of the cracks over a neighboring volume of finite thickness.In this way, the problem is reformulated in terms of displacement and phasefields defined over the entire volume of the domain. The governing equa-tions define a system of second-order partial differential equations, thus inprinciple eschewing the difficulties inherent to evolving boundaries and dis-continuities. However, the convenience of the initial implementation comesat the price of exceedingly fine discretization requirements in the vicinityof the cracks, a doubling of degrees of freedom, a sensitive dependence onthe choice of intrinsic length scale, strongly non-linear and possibly unstabledynamics, difficulties enforcing strict irreversibility, no-healing and positive dissipation, difficulties enforcing crack closure, difficulties enforcing mode-mixity dependent fracture criteria, onerous computing time requirementsand other difficulties, which need to be carefully addressed and assessed.Element erosion (ER) methods [11, 12, 13, 14], consisting of approxi-mating cracks as notches of small but finite width, supply another well-established class of computational methods for simulating crack growthwhich has been extensively used to simulate fracture in a number of areasof application, including fragmentation and terminal ballistics. In seminalwork, Negri [15] noted that some of the early versions of element erosionfail to converge, or converge to the wrong limit, due to mesh-dependency ofthe crack path, and provided a remedy based on the use of local averagesover intermediate scales [16, 17]. Subsequent enhancements of element ero-sion incorporating such local averaging [18, 19, 20, 21, 22, 23, 24, 25, 26]are provably convergent to Griffith fracture in the limit of vanishingly smallmesh sizes [18].Phase-field and element erosion methods have a common variational struc-ture: i) an elastic energy-release mechanism, namely, progressive damage inthe case of PF and abrupt damage in the case of ER; and ii) an energy cost ofdamage, derived from the phase-field and its gradients in the case of PF andfrom an estimate of the fracture area in the case of ER. In both cases, thestatic equilibrium configurations of the solid follow from global energy mini-mization. In addition, crack propagation is modeled in both cases by meansof a rate-independent gradient flow that balances elastic energy-release rateand dissipation.We begin by formalizing the commonalities between PF and ER methodsand showing how they are special cases of a common variational structurebased on the general notion of eigendeformation . Eigendeformations arewidely used in mechanics to describe deformation modes that cost no lo-cal energy, cf., e. g., [27]. In some sense, eigendeformations provide themost general representation of elastic energy-release mechanisms. Not sur-prisingly, therefore, the method of eigendeformations provides a commonframework for both PF and ER methods. Within this unified view, PF andER methods simply correspond to particular choices of restricted classes ofallowable eigendeformations and cost functions thereof.Whereas the convergence properties of finite-element approximations ofEE and PF is well established mathematically [28, 18], a direct quantitativecomparison and assessment of both methods appears to have been missing.In this work, we endeavor to fill that gap by means of selected numericaltests. By convergence we specifically understand convergence of the EE andPF solutions to the Griffith solution as the length parameter ǫ and the meshsize h both tend to zero. Thus, we regard the Griffith solution as exact andthe EE and PF solutions as approximations thereof. Given the variationaland energy minimization principles at work for both EE and PF, it is naturalto measure errors and convergence rates in terms of energy and endeavor toascertain the rate at which the EE and PF energies approach the limiting IGENEROSION
VS.
PHASE-FIELD 3
Griffith energy. We also recall that, for propagating cracks, the energyrelease rate supplies the requisite driving force for crack advance. Therefore,accuracy and convergence of the energy is a sine qua non prerequisite forthe accuracy and convergence of the crack tracking problem.It bears emphasis that we seek to characterize the convergence of solu-tions with respect to two parameters simultaneously, namely, ǫ and h . Thisdouble limit raises the fundamental question of the relative rates at which ǫ and h should be reduced to zero. Mathematical analysis [28, 18] showsthat convergence requires ǫ to decrease to zero more slowly than h , i. e., ǫ must be chosen on an scale intermediate between h and the size of thedomain. Remarkably, this requirement is contrary to rules of thumb oftenused in practice that recommend setting ǫ to a fixed multiple of h . Here, weinstead seek to optimize ǫ for given h as part of the approximation scheme.Specifically, for a given mesh size h we determine the optimal value ǫ h of ǫ by recourse to energy minimization, in the spirit of variational adaption.We show that ǫ h indeed yields the energy closest to the Griffith limit andthus minimizes the energy error for given mesh size h . The resulting mesh-size convergence plots for EE and PF may therefore be viewed as the bestpossible for each method, which makes the method comparison fair.We specifically consider the standard test case of a center-crack panelloaded in biaxial tension as a means of assessing the accuracy, convergencerate and computational cost performance of EE and PF. The panel is dis-cretized by means of a regular square mesh consisting of standard bilinearor Q Multi-field models of brittle fracture
According to Griffith’s criterion for fracture, in a brittle material crackgrowth is results from the competition between elastic energy minimizationand the fracture energy cost of creating new surface. Assuming rate inde-pendence, crack growth in a solid occupying a domain Ω ⊂ R is governedby the potential energyΠ( u ) = E ( u ) + (forcing terms) , (1) A. PANDOLFI, K. WEINBERG, M. ORTIZ where E ( u ) = Z Ω \ J u W ( ε ( u )) dx + G c | J u | , (2)is the total energy, including the elastic energy of the solid and the energycost of fracture, W ( ε ( u )) denotes the strain energy density, ε ( u ) = sym ∇ u the linearized strain tensor, u ( x ) the displacement field, dx the element ofvolume and the forcing terms (not spelled out for brevity) include bodyforces, boundary tractions and prescribed displacements. The jump set J u collects the cracks across which the displacement u may jump discontinu-ously and | J u | denotes the crack surface area. The material-specific param-eter G c is the specific fracture energy density per unit area and measuresthe fracture strength of the solid.The central and all-encompassing governing principle of energy minimiza-tion posits that the displacement field u at any given time is expected to minimize the potential energy Π( u ) subject to monotonicity of the jump set J u , i. e., to the constraint that J u must contain all prior jump sets, andto crack closure constraints. In this manner, the problem of crack trackingis reduced to a pseudo-elastic problem, with monotonicity and closure con-straints, for every state of loading. Such pseudo-elastic problems arise gen-erally for rate-independent inelastic solids under monotonicity constraints(cf., e. g., [30] for a rigorous derivation) and were initially formulated inconnection with deformation theory of plasticity [31].The problem thus defined is a free-discontinuity problem in the sense thatthe displacement field u is allowed to be discontinuous and the discontinuityor jump set J u itself, i. e., the crack surface in the present application, isan unknown of the problem. The existence and approximation propertiesof such problems have been extensively investigated in the mathematicalliterature (cf., e. g., [32] for a review). Free-discontinuity problems arenotoriously difficult to solve computationally, which has spurred the searchfor sundry regularizations of the problem that relax, to good computationaladvantage, the sharpness of the discontinuities. In the present work, wespecifically focus on two such regularizations, eigenfracture and phase-fieldmodels, which we briefly summarize next.2.1. Eigenfracture.
The method of eigenfracture (EF) is an approxima-tion scheme for generalized Griffith models based on the notion of eigende-formation [18]. The approximating energy functional is assumed to be ofthe form E ǫ ( u, ε ∗ ) = Z Ω W ( ε ( u ) − ε ∗ ) dx + G c ǫ |{ ε ∗ = 0 } ǫ | = E e ( u, ε ∗ ) + E iǫ ( ε ∗ )(3)where ε ∗ is the eigendeformation field that accounts for fracture, E e ( u, ε ∗ )is the elastic energy, E iǫ ( ε ∗ ) is the energy cost of the eigendeformation, or in-elastic energy, and ǫ is a small length parameter. The elastic energy E e ( u, ε ∗ ) IGENEROSION
VS.
PHASE-FIELD 5 follows as the integral over the entire domain of the strain energy density W as a function of the total strain ε ( u ) reduced by the eigenstrain ε ∗ . In thismanner, eigendeformations allow the displacement field to develop jumps atno cost in elastic energy.This local relaxation comes at the expense of a certain amount of fractureenergy. The challenge in regularized models of fracture is to estimate the in-elastic fracture energy E iǫ ( ε ∗ ) in a manner that converges properly as ǫ → ǫ -neighborhood { ε ∗ = 0 } ǫ of the support { ε ∗ = 0 } of the eigendeforma-tions scaled by 1 /ǫ , cf. Fig. 3b. Specifically, in this construction { ε ∗ = 0 } is the set of points where the eigendeformations differ from zero, { ε ∗ = 0 } ǫ is the ǫ -neighborhood of { ε ∗ = 0 } , i. e., the set of points at a distance to { ε ∗ = 0 } less or equal to ǫ , and |{ ε ∗ = 0 } ǫ | is the volume of { ε ∗ = 0 } ǫ .Remarkably, the method of eigenfracture is provably convergent [18], inthe sense that the total energy (3) Γ-converges to the Griffith energy (2) inthe limit of ǫ →
0. This convergence property shows that the eigenfracturemethod is indeed physically and mathematically sound. We recall that Γ-convergence of the energy functionals in turn implies convergence of thesolutions as ǫ →
0, i. e., the eigenfracture solutions converge to the solutionsof Griffith fracture in the limit of vanishingly small length parameter ǫ .2.2. Phase-field models of fracture.
In the PF approximation of Griffithfracture, the state of the material is characterized by an additional continu-ous field v ( x ) taking values in the interval [0 ,
1] and v = 0 at the crack. Thecrack set J u is then approximated as a diffuse interface where v = 1. Thecorresponding fracture model traces back to the pioneering work of Ambro-sio and Tortorelli [10], who showed that a two-field functional Γ-converges tothe Mumford-Shah functional of image segmentation. Generalized to three-dimensional elasticity, the two-field functional of Ambrosio and Tortorelliassumes the form E ǫ ( u, v ) = Z Ω (cid:16) ( v + o ( ǫ )) W ( ε ( u )) + G c (cid:16) (1 − v ) ǫ + ǫ |∇ v | (cid:17)(cid:17) dx = E eǫ ( u, v ) + E iǫ ( v ) , (4)where ǫ is a small length parameter and o ( ǫ ) stands in for a positive functionthat decreases to zero faster than the small parameter ǫ . The work of Ambro-sio and Tortorelli, and other similar works [33, 32], subsequently spawned nu-merous variants, extensions and implementations (e. g., [34, 35, 1, 15, 16, 36],but the differential structure of the fracture energy E iǫ ( v ) in (4) has remainedessentially unchanged in the later works. A. PANDOLFI, K. WEINBERG, M. ORTIZ
Eigenerosion.
Eigenerosion (EE) supplies an efficient implementationof the eigenfracture model [19]. To establish the connection between eigen-fracture, eq. (3), and eigenerosion, assume that W ( ε ) is quadratic and re-strict eigendeformations to the particular form(5) ε ∗ = ε ( u ) − ( w + o ( ǫ )) / ε ( u ) , with w taking the values 0 or 1, i. e., w ( x ) ∈ { , } . Inserted into eq. (3)this gives the EE functional E ǫ ( u, w ) = Z Ω ( w + o ( ǫ )) W ( ε ( u )) dx + G c ǫ |{ w = 0 } ǫ | = E eǫ ( u, w ) + E iǫ ( w ) . (6)By Jensen’s inequality and properties of extreme points [37], it follows thatthe range of w can be extended to the entire interval [0 , ≤ w ( x ) ≤ v = √ w and a fracture energy computed by the ǫ -neighborhood construction. Con-versely, PF models of fracture may be viewed as special cases of EE, andhence eigenfracture, where the fracture energy is of the Ambrosio-Tortorellitype.The great advantage of the EE model (6) vs . the conventional Ambrosio-Tortorelli-type phase-field model (4) is that in the former, eq. (6), the phase-field is undifferentiated and evaluates the fracture energy through an inte-gral expression, whereas the latter, eq. (4), requires the phase-field to bedifferentiated. Differentiation in turn requires regularity and conforming in-terpolation, e. g., by the finite-element method. By contrast, the integralform of the fracture energy in (6) allows the phase-field to be approximated,e. g., as piecewise constant 0 or 1, which leads to a considerable increase inimplementational simplicity and robustness [19, 20].2.4. Non-local fracture as an Artificial Neural Network.
Artificialneural networks provide a compelling interpretation of the ǫ -neighborhoodconstruction that suggests an entire class of extensions thereof. To makethis connection, we introduce the mollifier(8) ϕ ǫ ( x ) = (cid:26) / (4 πǫ / , | x | < ǫ, , otherwise , and the activation function(9) f ( w ) = (cid:26) , w ≤ , , otherwise . IGENEROSION
VS.
PHASE-FIELD 7 (cid:272) (cid:396) (cid:258) (cid:272) (cid:364) (cid:2659)(cid:882)(cid:374)(cid:286)(cid:349)(cid:336)(cid:346)(cid:271)(cid:381)(cid:396)(cid:346)(cid:381)(cid:381)(cid:282) (cid:272) (cid:396) (cid:258) (cid:272) (cid:364) (cid:272)(cid:381)(cid:374)(cid:448)(cid:381)(cid:367)(cid:437)(cid:410)(cid:349)(cid:381)(cid:374) (cid:258)(cid:272)(cid:410)(cid:349)(cid:448)(cid:258)(cid:410)(cid:349)(cid:381)(cid:374) (cid:349)(cid:374)(cid:410)(cid:286)(cid:336)(cid:396)(cid:258)(cid:410)(cid:349)(cid:381)(cid:374)
Figure 1.
Artificial Neural Network representation of the ǫ -neighborhood construction for the computation of the frac-ture energy.Next we note that the function(10) w ǫ ( x ) = Z Ω ϕ ǫ ( x − y ) w ( y ) d y = ( ϕ ǫ ∗ w )( x ) , obtained by taking the convolution of ϕ ǫ and w , is positive in the ǫ -neigh-borhood { w = 0 } ǫ of the crack set and vanishes elsewhere. Therefore, thefiltered function f ( w ǫ ( x )) is 1 in the ǫ -neighborhood { w = 0 } ǫ and van-ishes elsewhere, i. e. it is the characteristic function of the ǫ -neighborhood.Finally, we have(11) E iǫ ( w ) = G c ǫ Z Ω f ( w ǫ ( x )) dx = G c ǫ Z Ω f (cid:0) ( ϕ ǫ ∗ w )( x ) (cid:1) dx, which supplies an integral representation of the ǫ -neighborhood construc-tion.In (11), we immediately recognize the structure characteristic of an artifi-cial neural network (cf., e. g., [38]), Fig. 1. Thus, the fracture energy E iǫ ( w ),which is the output of the network, follows from the input w through thecomposition of three operations. The first operation is a convolution ϕ ǫ ∗ w ,which defines a linear neural network. The outcome w ǫ of this operation isfiltered locally by means of the activation function f in the form of a bi-nary switch, with the result that f ( w ǫ ) is the characteristic function of the ǫ -neighborhood of the crack set. Finally, the fracture energy follows as anintegral of f ( w ǫ ), suitably scaled by G c / ǫ . A. PANDOLFI, K. WEINBERG, M. ORTIZ
The artificial neural network interpretation suggests an extension of eigen-fracture to a more general class of fracture models where the fracture en-ergy is of the form (11) but ϕ ǫ and f and a general mollifier and activationfunction. This generalization immediately raises the question of how theaccuracy and convergence properties of the model depend on the choice of ϕ ǫ and f .3. Test case: Slit crack under all-around tension
We proceed to assess the accuracy and convergence of PF and EE ap-proaches by recourse to the standard test case of a slit crack in an infinitesolid subjected to all around tension, Fig. 2. We regard the Griffith solutionas exact and the EE and PF solutions as approximations thereof. Since boththe PF and EE problems are energy minimization problems, we specificallymonitor energy errors and seek to characterize the convergence with respectto the length parameter ǫ and mesh size h simultaneously. xx σ σ σ σ
02 1 (a) r θ r x θ x ra a θ (b) Figure 2. (a) Slit crack in infinite plate under all around tension σ . (b) Definition of the geometrical coordinates defining the stressstate around the crack. Exact reference results.
We specifically consider an infinite solidcontaining a straight crack of length 2 a deforming in plane strain under theaction of equibiaxial stress σ at infinity, see Fig. 2(a). The crack-tip stressfield is mode I since the loads are symmetric with respect to the crack line. IGENEROSION
VS.
PHASE-FIELD 9
Conveniently, the problem has an exact solution [39], namely, σ = σ r √ r r (cid:20) cos (cid:18) θ − θ − θ (cid:19) − a r r sin θ sin 32 ( θ + θ ) (cid:21) , (12a) σ = σ r √ r r (cid:20) cos (cid:18) θ − θ − θ (cid:19) + a r r sin θ sin 32 ( θ + θ ) (cid:21) , (12b) σ = σ r √ r r (cid:20) a r r sin θ cos 32 ( θ + θ ) (cid:21) , (12c)where the coordinates θ , θ , θ , r , r , and r are defined in Fig. 2(b) andthe lower indices correspond to cartesian coordinates ( x , x ) aligned andcentered at the crack. Assuming plane strain conditions, the strains followfrom Hooke’s law as ε = 1 − ν E σ − ν (1 + ν ) E σ , (13a) ε = 1 − ν E σ − ν (1 + ν ) E σ , (13b) γ = 2(1 + ν ) E σ , (13c)where E denotes the Young modulus and ν the Poisson’s ratio. The dis-placements u can then be computed by integrating the relations(14) ε = ∂u ∂x , ε = ∂u ∂x , γ = 2 ε = ∂u ∂x + ∂u ∂x , using Cesaro’s method. Finally, we recall that the strain-energy density ofthe solid is(15) W ( ε ) = λ ε + ε ) + 2 µ ( ε ) , λ = Eν (1 + ν )(1 − ν ) , µ = E ν ) . For a crack-free solid, the stresses reduce to(16) σ = σ = σ , σ = 0 , the strains to(17) ε = ε = (1 − ν )(1 + ν ) E σ , γ = 0 , and the strain-energy density to(18) W = (1 − ν )(1 + ν ) E σ . In order to facilitate numerical calculations, we restrict the analysis to abounded domain Ω surrounding the crack. To that end, we begin by notingthat the restriction of the infinite-body displacement field u to Ω minimizesthe total potential energy(19) Π( u ) = Z Ω W ( ε ( x )) dx dx − Z ∂ Ω σ ij n j u i ds where ε ( x ) are the strains attendant to the trial displacements u ( x ), ∂ Ωdenotes the boundary of Ω, n its outward unit normal and ds is the elementof arclength over ∂ Ω. The potential energy (19) represents a free-standingbody occupying the domain Ω deforming under the action of tractions σ ij n j corresponding to the stress field (12).In the absence of the crack, i. e., for crack length a = 0, an application ofClapeyron’s theorem gives(20) Π( u ) = − W | Ω | , where | Ω | is the area of Ω and W is given by (18). For a finite crack, theminimum value of the potential energy follows directly from an applicationof Rice’s J -integral [40]. Specifically, the energy-release rate is given by(21) − Π( u ) ∂a = Z Γ (cid:16) W ( ε ) n − σ ij n j u i, (cid:17) ds, where Γ denotes a counter-clockwise closed contour surrounding the crackand contained in Ω, n is its outward unit normal, and ds is the element ofarc-length on Γ. Choosing Γ to coincide with the flanks of the crack, togetherwith small loops at the tips, and using the asymptotic K -field gives(22) − Π( u ) ∂a = 2 1 − ν E K , where K I is the mode I stress-intensity factor, which can be computed di-rectly from the stress field (12), with the result(23) K I = σ √ πa . Inserting (23) into (22), integrating with respect to a , and using (20) with(18) as initial condition gives(24) Π( u ) = − (1 − ν )(1 + ν ) E σ | Ω | − − ν E πa σ , In addition, according to the Griffith model (2), the inelastic or fractureenergy expended in the extension of the crack is(25) E i ( u ) = G c a. The exact results (24) and (25) are subsequently taken as a convenient basisfor the analysis of the accuracy and convergence of EE and PF approxima-tion schemes. 4.
Discretization and implementation
Next we describe the discretization of the EE and PF models used incalculations. All calculations are performed on a square domain Ω of size D ≫ a . IGENEROSION
VS.
PHASE-FIELD 11 (a) ε (b) Figure 3. (a) Eigenerosion discretization of a slit crack. Thestring of shaded elements containing the crack are disabled. (b) ǫ -neighborhood construction for the calculation of the crack length. Eigenerosion model.
In the implementation of the EE model, thepotential energy(26) Π ǫ ( u, w ) = E ǫ ( u, w ) − Z ∂ Ω σ ij n j u i ds , with the energy E ǫ ( u, w ) as in (6), is discretized by means of a regularmesh consisting of standard bilinear or Q h . The exactstresses σ ij ( x ) in (26) are taken directly from the analytical solution (12).In order to disentangle the convergence with respect to the mesh size h fromquadrature error, in all the calculations we evaluate all element integralsover the interior and the boundary of the domain exactly using the symboliccomputation program Mathematica [29].The implementation of the EE fracture energy is illustrated in Fig. 3.The crack is represented as a string of missing, or ‘eroded’, elements where w = 0, dark gray elements in Fig. 3(a), with w = 1 elsewhere. The erodedelements approximate the crack geometry and the corresponding inelasticor fracture energy is approximated as(27) E i EE ,ǫ,h = G c ǫ A ǫ,h , where ǫ is a length parameter and A ǫ,h is the area of the ǫ -neighborhood ofthe eroded elements, light gray area in Fig. 3(b), i. e., the set of points at adistance smaller or equal to ǫ from the eroded elements. Examples concernedwith crack growth through slanted meshes have been presented in [19]. Fora slit crack centered and aligned with the mesh, the eroded elements areprecisely those which contain the crack, and A ǫ,h follows exactly as(28) A ǫ,h = ⌈ a/h ⌉ h + 2( ⌈ a/h ⌉ + 1) hǫ + πǫ . In this expression, ⌈ a/h ⌉ is the number of eroded elements. The ceilingfunction ⌈ x ⌉ equals the smallest integer larger or equal x . Inserting (28)into (27), we obtain(29) E i EE ,ǫ,h = G c ǫ (cid:16) ⌈ a/h ⌉ h + 2( ⌈ a/h ⌉ + 1) hǫ + πǫ (cid:17) . As expected, the EE approximation of the inelastic energy E i EE ,ǫ,h dependson the choice of length parameter ǫ . From a variational viewpoint, theoptimal value ǫ h of ǫ is that which minimizes E i EE ,ǫ,h , namely,(30) ǫ h = h r ⌈ a/h ⌉ π . The corresponding optimal inelastic energy is(31) E i EE ,ǫ h ,h ≡ E i EE ,h = G c h (1 + ⌈ a/h ⌉ + p π ⌈ a/h ⌉ ) . This energy is the minimum inelastic energy that can be attained for fixed h . Thus, For ǫ ≫ ǫ h , the error incurred by the ǫ -neighborhood construc-tion becomes dominant, causing E i EE ,ǫ,h to increase, whereas for ǫ ≪ ǫ h the under-resolution of the ǫ -neighborhood by the mesh size dominates andcauses E i EE ,ǫ,h to again increase. With(32) ⌈ a/h ⌉ = 2 ah + 2 δ, ≤ δ < , an asymptotic expansion of (30) gives in the limit of h/a ≪ ǫ h = r ahπ + O( h / ) . A similar asymptotic expansion of the optimal inelastic energy (31) likewisegives(34) E i EE ,h = G c a + G c √ πah + O( h ) . We observe from (33) that, to leading order, the optimal length parameter ǫ h scales as √ ah , i. e., as the geometrical mean of the crack length andthe mesh size. Thus, the optimal length parameter depends not only onthe mesh size but also on the geometry of the crack and, in general, of thedomain. We note that ǫ h → h →
0, albeit at a slower rate, as requiredby convergence [18].We also observe from (34) that, to leading order, the fracture energy erroris of order O( h / ). This rate of convergence is slower than the O( h ) rateof convergence of the elastic energy and, therefore, dominates the overallenergy error. However, this loss of convergence can be remedied simplyby recourse to a standard Richardson extrapolation technique (cf., e. g.,[41, 42]). To this end, we note from (34) that the fracture energy attendantto a mesh of size 2 h is(35) E i EE , h = G c a + 2 G c √ πah + O( h ) . IGENEROSION
VS.
PHASE-FIELD 13
We may replace E i EE ,h by the weighted sum(36) E i EE+RE ,h = λE i EE ,h + (1 − λ ) E i EE , h without disturbing the limit. We then choose the weight λ so as to cancelthe O( h / ) term, with the result(37) λ = √ √ − . From (34) and (36) it then follows that(38) E i EE+RE ,h = G c a + O( h ) , as desired. Inserting (31) and (37) into (36), we find E i EE+RE ,h = √ √ − G c h (1 + ⌈ a/h ⌉ + p π ⌈ a/h ⌉ ) − (1 + √ G c h (1 + ⌈ a/h ⌉ + p π ⌈ a/h ⌉ ) . (39)explicitly. This simple Richardson extrapolation effectively eliminates thelow-order accuracy of the original ǫ -neighborhood construction and restoresthe full order of convergence expected of the finite element method.4.2. Phase-field model.
In order to have a fair comparison with EE, wediscretize the potential energy(40) Π ǫ ( u, v ) = E ǫ ( u, v ) − Z ∂ Ω σ ij n j u i ds with the energy E ǫ ( u, v ) as in (4), by means of a regular mesh likewiseconsisting of standard bilinear or Q h . Conforming inter-polation is used for both the displacement and the phase fields. In orderto separate interpolation errors from numerical quadrature errors, we againevaluate all element integrals over the interior and boundary of the domainexactly using the symbolic computation program Mathematica [29]. Thephase field is unconstrained and, therefore, satisfies free Neumann bound-ary conditions on the outer boundary of the domain.The unknown fields ( u h , v h ) are solved iteratively by the method of al-ternating directions [43], i. e., by successively fixing one of the fields andsolving for the other. Conveniently, the scheme reduces the solution to asequence of linear problems (cf., e. g., [44, 45]). The iteration is primed bysetting, as initial condition for the iteration, v h = 0 on the nodes lying onthe crack and v h = 1 elsewhere. The iteration may then be expected toapproximate the solution for a crack of length 2 a provided that the appliedstress σ equals the critical stress for crack extension, i. e., if(41) G c = 1 − ν E K I , K I = σ √ πa , as in this case the exact elasticity solution (12) and the crack length 2 a jointly minimize the Griffith potential energy (1). The choice of the length parameter ǫ , and its relation to the mesh size h and geometrical features of the domain, is known to have a strong effecton the accuracy and convergence properties of PF approximations [28]. Inparticular, convergence requires h to decrease to zero faster than ǫ ([28],Theorems 4.1 and 5.1). We expect the exact potential energy (24) to beapproached by the PF potential energy from above (cf., e. g., Fig. 4b). Forfixed h , we also expect the PF potential energy to exhibit a minimum ata certain value ǫ h of ǫ . Thus, for small ǫ the mesh size h is unable toresolve the width of the crack, resulting in an overly stiff response and highPF potential energy. Contrariwise, since the exact potential energy (24) isattained from above for ǫ →
0, the PF potential energy diverges from theexact value for large ǫ . As a consequence of these opposing trends, for fixed h the PF potential energy attains a minimum at a certain ǫ h , as surmised(cf., e. g., Fig. 4b). Since, as already noted, the exact potential energy (24)is approached by the PF potential energy from above, the energy minimizing ǫ h results in the least energy error for given h and is therefore optimal fromthe standpoint of convergence.Evidently, ǫ h depends on the geometry of the crack and of the body,the state of damage and the discretization. In calculations, we determine ǫ h numerically by computing the PF minimum potential energy for given h over a range of equally-spaced values of ǫ , interpolating the computedenergies as a function of ǫ and computing the minimum of the interpolatedfunction (cf., e. g., Fig. 4b).5. Numerical Results E ν G c σ a −
10 0.403125
Table 1.
Parameters used in the numerical calculations.We proceed to investigate the relative accuracy and convergence rates ofthe EE and PF methods by way of numerical testing. To this end, we fixthe domain size at D = 5, the crack length at 2 a = 0 . h/D = 0 .
02, 0 .
01, 0 . . . Optimal choice of length parameter.
Fig. 4 shows the computeddependence of the minimum potential energy Π ǫ on ǫ at fixed h for both EE,eq. (26), and PF, eq. (40). As surmised, at fixed h the minimum potentialenergy Π ǫ attains a minimum at a well-defined optimal value ǫ h for both EEand PF. For EE, ǫ h is given analytically and in close form by (30). For PF, ǫ h is determined numerically by interpolating the results shown in Fig. 4b. IGENEROSION
VS.
PHASE-FIELD 15 M i n i m u m po t en t i a l ene r g y h/D = 0.02h/D = 0.01h/D = 0.005h/D = 0.0025h/D = 0.00125Exact (a) EE M i n i m u m po t en t i a l ene r g y h/D = 0.02h/D = 0.01h/D = 0.005h/D = 0.0025h/D = 0.00125Exact (b) PF Figure 4.
Dependence of the minimum potential energy on thelength parameter ǫ for various mesh sizes h . The limiting value ofthe minimum potential energy, (24), for ǫ → As is evident from the figure, the optimized values ǫ h of the length parameterminimize the potential energy error for fixed h and thus result in the bestpossible rate of energy convergence with respect to mesh size h . Mesh size h O p t i m a l l eng t h sc a l e -2 -1 NumericalEEEE+RE (a) EE
Mesh size h O p t i m a l l eng t h sc a l e -4 -3 -2 (b) PF Figure 5.
Optimal value ǫ h of the length parameter ǫ as a func-tion of mesh size h . (a) Eigenerosion, eqs. (30) and (33). (b)Phase-field values obtained from Fig. 4b. The optimal ǫ h values for EE and PF, i. e., the minimizers of the curvesdisplayed in Fig. 4, are shown in Fig. 5 as a function of the mesh size h . ForEE, Fig. 5(a) simply displays the theoretical optimal values, eqs. (30) and h / 2a P o t en t i a l E ne r g y E rr o r -3 -2 -1 PFEEEE+RE (a) h C P U T i m e [ s ] -2 -1 EEPF (b)
Figure 6. a) Energy convergence plots showing normalized po-tential energy errors vs. normalized mesh size h/ a for eigenerosion(EE), eigenerosion with Richardson extrapolation (EE+RE) andphase-field (PF). b) Execution times as a function of mesh size h for eigenerosion and phase-field (33) for purposes of comparison. In all cases, the dependence of the optimal ǫ h on the mesh size h is strongly suggestive of a power-law scaling(42) ǫ h ∼ h / . As expected [28, 18], the optimal length parameter ǫ h converges to zeromore slowly than the mesh size h . The scaling law (42) follows heuristicallyif we assume that near the crack the phase field varies on the scale of ǫ /L ,where L is a characteristic size (intrinsic geometric feature, e. g., domainsize, crack size, ligament size). In order to resolve this variation, we mustchoose h ∼ ǫ /L , whence (42) follows.The square-root nonlinearity of the scaling law (42), or equivalently ǫ h ∼√ Lh , results in optimal values of the length parameter that may be stronglyincommensurate with h , contrary to standard computational practice. Thus,for fine meshes, h ≪ L it follows that ǫ h ≫ h , i. e., ǫ h lags behind h and h resolves ǫ h finely. Contrariwise, for coarse meshes, h ≫ L , we have ǫ h ≪ h ,i. e., ǫ h runs ahead of h and is unresolved by h . This latter regime is clearlyvisible in Fig. 4b, e. g., at h = 0 .
1, for which the corresponding ǫ h is overone order of magnitude smaller.5.2. Energy convergence.
Potential energy convergence plots for eigen-erosion (EE), eigenerosion with Richardson extrapolation (EE+RE) andphase-field (PF) are shown in Fig. 6a. The plot displays least-square fitsof the data by functions of the form(43) inf Π ǫ h = Π + Ch α , IGENEROSION
VS.
PHASE-FIELD 17 with respect to the parameters C and α . The limiting potential energy Π in (43) is given by (24). The exponent α is the rate of convergence and theconstant C shifts the error line vertically in log-log convergence plots. Allenergy errors are computed using the optimal ǫ h corresponding to each meshsize, computed as described earlier. The figures visualize the least error inminimum potential energy as a function h for each method. All terms arenormalized: the mesh size h is normalized by the crack length, 2 a , and thepotential energies are normalized by Π .As may be seen from Fig. 6a, both PF and EE exhibit sublinear (square-root) energy convergence, though the constant for EE is nearly one order ofmagnitude smaller (better) than that for PF, indicating superior accuracyof EE over PF under the conditions of the test. By far the best accuracyand convergence rate is obtained for EE+RE, which exhibits linear energyconvergence, or twice the rate of convergence of EE and EF, and a betterconstant than that of both EE and PF.5.3. Computational cost.
Fig. 6b shows a comparison between executiontimes for EE and PF. We recall that the preceding convergence calculationsare carried out using exact integration of the element energy and nodalforces. This procedure has the advantage of singling out interpolation errorsand deconvolving them from other sources of error such as numerical quad-rature. However, the valuation of the exact integrals is costly and deviatesfrom common practice, which invariable relies on numerical quadrature as afurther approximation. Therefore, in order to obtain practical estimates ofexecution times, we repeat all calculations using standard finite element nu-merical quadrature. Specifically, we use four-point Gaussian quadrature forthe element integrals and two-point Gaussian quadrature for the boundary-edge integrals. All calculations are performed using a sparse linear solver [46]in memory-shared configuration, using a single node, 16-core Intel Skylake(2.1 GHz), with 192 GB Memory and 2666 MT/s speed.Fig. 6b shows that, under the conditions of the test and for the particularcomputer architecture used in the calculations, EE is about one order ofmagnitude faster than PF. The higher efficiency of EE over PF is in fact ex-pected, since EE entails displacement degrees of freedom and a no-iterationdirect solution, whereas PF entails both displacement and phase-field de-grees of freedom and an iterative solution.6.
Summary and concluding remarks
We have presented a comparison of the accuracy, convergence and com-putational cost performance of the eigenerosion (EE) and phase-field (PF)methods for brittle fracture. Both approaches operate on the principle ofminimization of a potential energy functional that accounts for the elasticenergy of the system, the inelastic energy attendant to crack growth and thework of the applied loads. Both approaches can be derived as special casesof the general method of eigendeformations by effecting particular choices of the eigendeformation field. The energy functionals incorporate an intrinsiclength scale ǫ representing an effective crack width. The solution for Griffithfracture is obtained in the limit of ǫ → Acknowledgements
MO gratefully acknowledges funding by the Deutsche Forschungsgemein-schaft (DFG, German Research Foundation) via project 211504053 - SFB1060 and project 390685813 - GZ 2047/1 - HCM.
References [1] B. Bourdin, G. A. Francfort, and J.-J. Marigo. Numerical experiments in revisitedbrittle fracture.
Journal of the Mechanics and Physics of Solids , 48:797–826, 2000.[2] B. Bourdin and A. Chambolle. Implementation of an adaptive finite element approxi-mation of the Mumford-Shah functional.
Numerische Mathematik , 85:609–646, 2000.[3] A. Karma, D. A. Kessler, and H. Levine. Phase-field model of mode III dynamicfracture.
Physical Review Letters , 81:045501, 2001.[4] C. Miehe, M. Hofacker, and F. Welschinger. A phase field model for rate-independentcrack propagation: Robust algorithmic implementation based on operator splits.
Computer Methods in Applied Mechanics and Engineering , 199:2765–2778, 2010.[5] M. J. Borden, C. V. Verhoosel, M. A. Scott, T. J. R. Hughes, and C. M. Landis.A phase-field description of dynamic brittle fracture.
Computer Methods in AppliedMechanics and Engineering , 217–220:77–95, 2012.[6] C. V. Verhoosel and R. de Borst. A phase-field model for cohesive fracture.
Interna-tional Journal for Numerical Methods in Engineering , 96(1):43–62, 2013.
IGENEROSION
VS.
PHASE-FIELD 19 [7] M. Ambati, T. Gerasimov, and L. De Lorenzis. A review on phase-field models of brit-tle fracture and a new fast hybrid formulation.
Computational Mechanics , 55(2):383–405, 2015.[8] C. Bilgen and K. Weinberg. On the crack-driving force of phase-field models in lin-earized and finite elasticity.
Computer Methods in Applied Mechanics and Engineer-ing , 353:348–372, 2019.[9] L. Ambrosio and V. M. Tortorelli. Approximation of functionals depending on jumpsby elliptic functionals via Γ-convergence.
Communications on Pure and Applied Math-ematics , 43:999–1036, 1990.[10] L. Ambrosio and V. M. Tortorelli. On the approximation of free discontinuity prob-lems.
Bollettino dell’Unione Matematica Italiana 6-B , 7:105–123, 1992.[11] G. R. Johnson and R. A. Stryk. Eroding interface and improved tetrahedral elementalgorithms for high velocity impacts in three dimensions.
International Journal ofImpact Engineering , 5:414–427, 1987.[12] T. Belytschko and J. Lin. A three-dimensional impact-penetration algorithm witherosion.
International Journal of Impact Engineering , 5(1–4):111–127, 1988.[13] M. Ortiz and A. E. Giannakopoulos. Crack propagation in monolithic ceramics undermixed mode loading.
International Journal of Fracture , 44:233–258, 1990.[14] T. Borvik, O.S. Hopperstad, and K.O. Pedersen. Quasi-brittle fracture during struc-tural impact of AA7075-T651 aluminum plates.
International Journal of Impact En-gineering , 37:537–551, 2008.[15] M. Negri. Finite element approximation of the Griffith’s model in fracture mechanics.
Numerische Mathematik , 95(4):653–687, 2003.[16] M. Negri. A discontinuous finite element approximation of free discontinuity probems.
Advances in Mathematical Sciences and Applications , 15:283–306, 2005.[17] M. Negri. A non-local approximation of free discontinuity problems in
SBV and
SBD . Calculs of Variations and Partial Differential Equations , 25:33–62, 2005.[18] B. Schmidt, F. Fraternali, and M. Ortiz. Eigenfracture: an eigendeformation ap-proach to variational fracture.
SIAM Multiscale Modeling & Simulation , 7(3):1237–1266, 2009.[19] A. Pandolfi and M. Ortiz. An eigenerosion approach to brittle fracture.
InternationalJournal for Numerical Methods in Engineering , 92(8):694–714, 2012.[20] A. Pandolfi, B. Li, and M. Ortiz. Modeling fracture by material-point erosion.
Inter-national Journal of Fracture , 184(1-2):3–16, 2013.[21] L. Bichet, F. Dubois, Y. Monerie, C. P´elissou, and F. Perales. An eigenerosion methodfor heterogeneous media.
Materiaux et Techniques , 103(3), 2015.[22] F. Stochino, A. Qinami, and M. Kaliske. Eigenerosion for static and dynamic brittlefracture.
Engineering Fracture Mechanics , 182:537–551, 2017.[23] P. Navas, R.C. Yu, B. Li, and G. Ruiz. Modeling the dynamic fracture in concrete:an eigensoftening meshfree approach.
International Journal of Impact Engineering ,113:9–20, 2018.[24] A. Qinami, E.C. Bryant, W.C. Sun, and M. Kaliske. Circumventing mesh bias byr- and h-adaptive techniques for variational eigenfracture.
International Journal ofFracture , 220(2):129–142, 2019.[25] K. Zhang, S.-L. Shen, and A. Zhou. Dynamic brittle fracture with eigenerosion en-hanced material point method.
International Journal for Numerical Methods in En-gineering , 121(17):3768–3794, 2020.[26] A. Qinami, A. Pandolfi, and M. Kaliske. Variational eigenerosion for rate-dependentplasticity in concrete modeling at small strain.
International Journal for NumericalMethods in Engineering , 121(7):1388–409, 2020.[27] T. Mura.
Mechanics of Defects in Solids . Martinus Nijhoff, Dordrecht, 1987.[28] G. Bellettini and A. Coscia. Discrete approximation of a free discontinuity problem.
Numerical Functional Analysis and Optimization , 15(3-4):201–224, 1994. [29] Wolfram Research, Inc. Mathematica, Version 12.1. Champaign, IL, 2020.[30] L. Fokoua, S. Conti, and M. Ortiz. Optimal scaling in solids undergoing ductilefracture by void sheet formation.
Archive for Rational Mechanics and Analysis ,212(1):331–357, 2014.[31] J. B. Martin.
Plasticity : fundamentals and general results . MIT Press, Cambridge,Mass, 1975.[32] L. Ambrosio, N. Fusco, and D. Pallara.
Functions of Bounded Variation and FreeDiscontinuity Problems . Oxford University Press, Oxford – New York, 2000.[33] A. Braides and G. Dal Maso. Nonlocal approximation of the Mumford-Shah func-tional.
Calculus of Variations and Partial Differential Equations , 5:293–322, 1997.[34] G. Cortesani and R. Toader. Implementation of an adaptive finite element approx-imation of the Mumford-Shah functional.
Numerical Functional Analysis and Opti-mization, , 18:921–940, 1997.[35] A. Chambolle and G. Dal Maso. Discrete approximation of the Mumford-Shah func-tional in dimension two.
ESAIM: Mathematical Modelling and Numerical Analysis ,33:651–672, 1999.[36] S. Conti, M. Focardi, and F. Iurlano. Phase field approximation of cohesive fracturemodels.
Annales de l’Institut Henri Poincar´e C, Analyse non lin´eaire , 33(4):1033–1067, 2016.[37] A. J. Larsen. Local minimality and crack prediction in quasi-static griffith fractureevolution.
Discrete & Continuous Dynamical Systems-S , 6(1):121, 2013.[38] M. H. Hassoun.
Fundamentals of artificial neural networks . MIT press, 1995.[39] A. T. Zehnder.
Fracture Mechanics . Springer, 2012.[40] J. R. Rice. A path independent integral and the approximate analysis of strain con-centration by notches and cracks.
Journal of Applied Mechanics , 35:379–386, 1968.[41] C. Brezinski and M. Redivo-Zaglia. Extrapolation methods.
Applied Numerical Math-ematics , 15(2):123–131, 1994.[42] Z. Zlatev, I. Dimov, I. Farag´o, and A. Havasi.
Richardson Extrapolation: PracticalAspects and Applications . De Gruyter, Berlin/Boston, 2017.[43] D. W. Peaceman and H. H. Rachford, Jr. The numerical solution of parabolic andelliptic differential equations.
Journal of the Society for Industrial and Applied Math-ematics , 3(1):28–41, 1955.[44] C. Bilgen, A. Kopaniˇc´akov´a, R. Krause, and K. Weinberg. A phase-field approach toconchoidal fracture.
Meccanica , 53(6):1203–1219, 2018.[45] D. Knees and M. Negri. Convergence of alternate minimization schemes for phase-field fracture and damage.
Mathematical Models and Methods in Applied Sciences ,27(09):1743–1794, 2017.[46] X.S. Li, J.W. Demmel, J.R. Gilbert, L. Grigori, M. Shao, and I. Yamazaki. SuperLUUsers’ Guide. Technical Report LBNL-44289, Lawrence Berkeley National Labora-tory, September 1999. Dipartimento di Ingegneria Civile e Ambientale, Politecnico di Milano,Piazza Leonardo da Vinci 32, 20133 Milano, Italy. Department of Mechanical Engineering, Universit¨at Siegen, 57068 Siegen,Germany. Division of Engineering and Applied Science, California Institute of Tech-nology, 1200 E. California Blvd., Pasadena, CA 91125, USA.
Email address ::