Comparison of peridynamic and phase-field models for dynamic brittle fracture in glass

We report computational results obtained with three different models for dynamic brittle fracture. The results are compared against recent experimental tests on dynamic fracture/crack branching in glass induced by impact. Two peridynamic models (one using the meshfree discretization, the other being the LS-DYNA’s discontinuous-Galerkin implementation) and a phase-field model lead to interesting and important differences in terms of reproducing the experimentally observed fracture behavior and crack paths. We monitor the crack branching location, the angle of crack branching, the crack propagation speed, and some particular features seen in the experimental crack paths: small twists/kinks near the far edge of the sample. We discuss the models’ performance and provide possible reasons behind the failure of some of the models to correctly predict the observed behavior.


Introduction
Crack branching is a common occurrence in dynamic brittle fracture and has been studied for many decades, theoretically [1,2], experimentally [3][4][5][6][7], and computationally [7][8][9][10]. Some analytical solutions for crack branching in brittle systems exist, but they are generally limited to a crack running in an infinite domain and at constant speed [11] . For practical engineering applications, in which boundary conditions influence the evolution of dynamic cracks and multiple cracks affect one another, computational methods have been employed to explain dynamic brittle fracture experimental observations [12][13][14][15][16][17][18][19][20]. Invariably, experimental setups designed to study this problem in detail consider a pre-crack (from which the main crack that eventually branches initiates) in order to make tacking of crack path/speed easier, and also makes dynamic crack growth a matter of the fracture energy, not of material's ultimate strength. The cohesive zone model (CZM) [21], many times used in the context of the extended finite element method (XFEM) [18,22,23], are examples of FEM-based models that can simulate crack growth. XFEM uses special enrichment functions (based on analytical results from fracture mechanics) near the crack tip and along the crack faces. The embedded discontinuities in the XFEM displacement field reduce the need for remeshing usually associated with FEM modeling of crack growth [24]. Such models have been used for modeling cracking from low speed impact of brittle materials [25]. The XFEM, as described in [26], requires a damage criterion and tracking of stresses around the crack tip to decide when and how crack branching is to be inserted [24].
The study on dynamic brittle fracture in [27] compared results from three FEM-based methods: XFEM, inter-element crack method, and element deletion method. In modeling of dynamic crack branching in a glass sheet, XFEM did not capture roughening of the crack front (parasitic cracks [28]) observed in experiments to happen before the crack branching event, and crack speeds were significantly above the experimentally measured values. The inter-element and element deletion methods were less accurate in capturing the observed crack trajectories [27].
Although atomistic [29] and particle-type [30] models can simulate crack branching in brittle materials without a specific criterion, these methods also find crack propagation speeds and crack branching angles different from what is observed experimentally [26]. The likely reason is that dynamic crack propagation (including crack branching) can be significantly affected by wave reflections from boundaries, for example [8], requiring the entire structure to be modeled and thus being outside of the reach of atomistic calculations.
Continuum-based approaches like the phase-field model (PFM) and peridynamics (PD) have shown promising results in fracture problems. PFM [31][32][33] has been used to model crack propagation, merging cracks, and crack branching in brittle solids. In PFMs, the crack is considered as diffuse damage instead of a sharp discontinuity. For brittle fracture, this regularization is obtained using a variational approach for the classical Griffith energy balance. The phase-field (damage) variable, introduced as an extra degree of freedom per node in a finite element (FE) discretization, allows one to approximate the fracture paths [34]. By solving an energy minimization problem, damage evolution (and implicitly, cracks) are modeled without the need for remeshing [34]. Quasi-static simulations of brittle fracture with PFM have shown some excellent match with experiments [35,36], but also some crack features (in thermally driven cracks, for example) not seen experimentally [37]. Peridynamics, a nonlocal form of continuum mechanics, was introduced in [38] for modeling damage and fracture [39][40][41][42]. Since then, it has been extended to a variety of other problems in which domain changes/discontinuities are part of the solution [43][44][45][46]. Peridynamics has been successful in modeling crack branching [8,26,47] and has been mathematically proven to recover both static and dynamic Griffith's fracture [48,49]. An extensive review and catalog of available experimental setups that have been used for calibration/validation of numerous peridynamic-based simulations has recently appeared in [50]. An implementation of PD based on the discontinuous-Galerkin method has recently been introduced in the commercial code LS-DYNA for simulating fracture in brittle materials [51].
Soda-lime glass is the prototype brittle material, commonly used in virtually all industries. Brittle cracks in soda-lime glass usually start from pre-existing surface flaws [52]. Crack branching in glass can occur when stress waves interact with a rapidly moving crack or when the strain energy around the fracture process zone leads to a wave "pile-up" in front of the crack tip, leading to a split crack path [8]. Several attempts have been made to use dynamic crack propagation in brittle materials to compare the performance of different approaches [17,[53][54][55]. Validating the results from computational simulations against experiments that provide sufficient amount of details is useful not only for allowing one to select the most appropriate model for the problem at hand but also to improve the models in case they show deficiencies.
In this paper, we test some of the most recent methods that showed promise in simulating dynamic brittle fracture. We evaluate the performance of bond-based PD with the meshfree implementation [38], the LS-DYNA's implementation of PD with discontinuous Galerkin method [56], and the PFM (as described in [20,32]). We test these three methods on some recently published impact-induced crack branching experiments on soda-lime glass [4] that provide detailed information regarding crack propagation speed, branching angle, and point of branching. In this particular experiment, the crack path (post branching) presents some interesting fine details as it approaches the far end of the sample. Such details will be used to see if they are reproduced by the fracture models tested here.
The article is organized as follows: section 2 presents a brief review of peridynamics and discusses the meshfree discretization and the Discontinuous-Galerkin implementation in LS-DYNA; in section 3, a brief review of the phase-field model for fracture is given; in section 4 we review the experimental results that will be used for validating the numerical models; numerical results from the three models are shown in section 5, and these are discussed in section 6; conclusions are given in section 7.

A brief review of peridynamics
In peridynamics [38], each material point is connected through a peridynamic bond to other points within a certain neighborhood region called "the horizon". The peridynamic bonds transfer forces between points and their failure determine the damage information at each point. In peridynamics, one replaces the equation of motion by an integrodifferential equation in which spatial derivatives are eliminated. This allows peridynamics to easily resolve mathematical difficulties and inconsistencies present in the classical theory when damage (diffuse or localized as cracks), for example, develops and evolves in the domain. The equations of motion in PD are: where u is the displacement vector field as a function of point and time t, ρ is the mass density, is the body force density. Also, f is the peridynamic pairwise force function that describes the interaction between material points and ̂. Moreover, the horizon region is the internal sub-region (see Fig. 1), defined as: Let = (̂, ) − ( , ) be the relative displacement and =̂− be the relative position in the reference configuration between two material points of ̂ and . The interaction between material points is set to zero beyond the horizon size δ, such that When the pairwise force derives from a micro-elastic potential , a micro-elastic material is defined as: A linear micro-elastic material is obtained if we consider: where ( ) is called the bond micromodulus function and s is the relative elongation of the bond, or the bond strain: The pairwise force derived from Eq. 4 and Eq. 5 is The material model we use here is one with n=1 in the PD kernel function ( ,̂)/|̂− | [57]. Other options are possible, including the one with n=2 ( [57,58]) which leads to convergence properties to the classical solution that are independent from the discretization size. In this study, we use the plane stress condition, and the corresponding conical micromodulus(see [59]): where E is the elastic Young's modulus, and = 1/3 is the Poisson ratio. At a point x the elastic strain energy density ( ) is obtained by integrating Eq. 5 over the horizon: In peridynamics, damage is modeled using a critical bond strain concept, allowing a bond to break and no longer sustains any force [38,60] once its strain crosses a critical value. In this study, once a peridynamic bond breaks, it remains broken [44,45]. The critical bond strain parameter is obtained from the measured fracture energy, 0 [60]. In 2D, the connection between s0 and G0 is [26] For the conical micromodulus given in Eq. 9, one obtains the critical strain as (see [61])

Brief review of the meshfree PD discretization
The PD equations of motion (see Eq. 1) can be discretized with a variety of methods [42,62]. Here, we use meshfree discretization with one-point Gaussian quadrature over the domain [60] which uses a uniform grid with spacing and nodes centered in the middle of each such cells. We use the partial volume integration algorithm given in [26] for the numerical quadrature. A comparison of different partial volume schemes with an exact computation is given in [63]. Some techniques were introduced in [64] to modify the micromodulus near a boundary to account for the incomplete horizon and reduce the so called "peridynamic surface effect", and, similarly, to modify the critical strain at nodes that have suffered some damage already, into the so-called "damage-dependent material" [62,65]. We do not consider such modifications in this study. The meshfree discretization provides great advantages for problems that involve fracture, crack branching, and fragmentations [8,66]. The discretized peridynamic equations of motion for Equation 1 are: For the explicit time integration, we use the Velocity-Verlet algorithm [26,64]: where , ̇ and, ̈ denote displacement, velocity, and acceleration vectors, respectively.
Like in the regular FEM method based on the continuous Galerkin method, the approximation fields of the solution ( ) and test function ( ), are constructed using shape functions: where denotes the nodal displacement of node .
In the discretization of the continuous FEM, the adjacent elements share nodes whereas each element has its own nodes in the discontinuous-Galerkin FEM domain, i.e., the total nodal number equals the total element number times the node number per element. Equation 14 contains two levels of integration. The first level is over the computational domain, discretized using Gaussian quadrature: where denotes the total number of Gaussian points in the domain. The remaining integration in Equation 16 represents the nonlocal effect at a Gaussian point . Because the first level of integration is carried out through Gaussian points, the bond connectivity is defined as: where ( ) is a bond between Gaussian point and its neighbor Gaussian point ̂ in its horizon.
Therefore, the integration domain in Equation 16 which is a sphere centered at can be discretized by Gaussian points. Finally, the discretized bond-based explicit dynamics peridynamics governing equation becomes: where ̂ is the total number of neighbors of and is the FEM shape function. For each Gaussian point , B loops over all the nodes of its element and will eventually loop all the nodes in the mesh after all Gaussian points ( ) are calculated. Because the approximation field is constructed by regular FEM shape function, the boundary conditions can be enforced as in the standard FEM. The boundary effect of the nonlocal constitutive law is to be corrected in the sense of equivalent elastic energy density throughout the domain, as described in detail in [56].
The horizon size ( = 1.01 ) is defined based on the characteristic length ( ) of the element which is the diagonal length of each element [51,56]. The implies that δ is not constant over the domain if a non-uniform mesh is used. In this paper, since the problem domain is not rectangular, we use non-uniform discretizations, but the departure from non-uniformity is small. LS-DYNA creates peridynamic bonds automatically before the solution starts (see Fig. 2 for a 2D case). A pre-crack can be defined with a discontinuity in the geometry (see Fig. 2). The material model used for PD simulating crack is MAT_ELASTIC_PERI which takes the fracture energy as the input parameter. LS-DYNA uses ELFORM #48 as the peridynamic solid element formulation [67].

A brief review of phase-field modeling
In phase-field modeling (PFM) discontinuities are avoided by introducing a scalar field ( , ) used to regularizes a sharp crack surface topology by diffusive crack zones [20,68,69]. The phase-field depends on position and time and varies in [0,1]. The material is intact when = 0 and is fully broken when = 1. PFM for brittle fracture was proposed in [70], based on Griffith's theory. PFM considers a reference configuration Ω, external Ω, and the internal discontinuity boundary Γ(see Fig. 3). The crack surface density function of the solid is introduced in [32] as follows: where is the phase-field variable, 0 is a length scale parameter handling the transition region of , and ∇ is the spatial gradient of the phase-field variable. The sharp crack topology is obtained as 0 → 0.
Assuming the critical fracture energy density, , as independent of the crack velocity, the surface energy due to the formation of cracks can be approximated as follows: As discussed in [71], the jump of the elastic energy during crack evolution is exactly offset by that of the surface energy. In order to have cracks grow only under tension (not under compression) [31], the elastic energy density, ( ), where is the strain tensor, into tensile ( + ( )) and compressive ( − ( )) parts and the effect of phase-field is only assumed on the tensile elastic energy density: where ( ) is called degradation function representing the change of the elastic properties of the material between the broken and undamaged states [72]. As suggested in [71], ( ) is commonly chosen Also, is a small stability parameter for avoiding numerical singularities as the tends to 1. The total Lagrange energy functional which should be minimized can be expressed in terms of the independent fields ( , t) and ϕ( , t) is as follows: where = is Cauchy stress tensor defined by Both the displacement field ( , t) and phase-field ϕ( , t) can be solved from these equations of motion. Moreover, to prevent crack healing, the irreversibility condition Γ( , s) ∈ Γ( x, t) ( < ) is enforced in the strong-form equations by introducing a strain energy history field, ( , ), which satisfies the Kuhn-Tucker conditions for loading and unloading.
unloading or compression, + − < 0 and ̇( , )= 0. In either case ̇( + − ) = 0. In fact, ( , ) represents the maximum of the tensile part of the elastic energy in the deformation history of point : By replacing + ( , ) by ( , ) in Eq. 23, the strong form is rewritten as: Phase-field formulation can be implemented in any FEM commercial software, like ABAQUS [73][74][75][76] or open source FEM package [77]. Phase-field modeling is inherently a multi-field problem. An implementation of the PFM in COMSOL Multiphysics is presented in [78][79][80]. PFM has been implemented using implicit and recently with explicit methods [81]. We use the same COMSOL code in the present work in which an implicit staggered time integration scheme is adopted to solve the coupled system of equations. In PFM, monolithic or staggered algorithms can be implemented. In monolithic scheme the displacement and the phase-field variables are determined simultaneously using one loop of iterations, which is more efficient. However, the staggered implementation calculates the field variables alternatively leading to a significantly better robustness [32]. Our code uses three main modules: Solid Mechanics Module, History-strain Module, and Phase-field Module which solve the three fields of , and , respectively. Each module is solved independently and sequentially during each iteration. These three modules are used along with the Storage Module to calculate and store the internal variables during each time step. These modules are all written in strong forms and solved based on standard finite element discretization in space domain and finite difference discretization in the time domain [78]. The relationship between the modules, the corresponding governing differential equations, boundary conditions, and initial conditions are shown in Fig. 4. The process of solution in COMSOL is as follows:  The tensile part of elastic energy, + , is input into the History-strain module to update the maximum of the tensile part of the elastic energy in the deformation history, , of the integration point.
 Finally, the phase-field variable, , is obtained from Phase-field module and is imported to Storage Module to calculate the relative error and determine if the convergence criterion is met. b) The "Solid Mechanics Module" solves for the updated displacement field: ( new ).
c) If the convergence criterion (norm-2 relative error for displacements < ) is not satisfied, we repeat step (b) with initial guess fields: ( new , , ). Otherwise, new is considered as +1 , ( + 1)th guess, and continue to step (d).
e) "Phase-field Module" uses +1 and +1 to update the phase-field variable field ( new ).
f) If the convergence criterion (‖ − ‖ 2 /‖ ‖ 2 < ) is not satisfied, we repeat step (e) with initial guess of: ( +1 , +1 , new ). Otherwise, new is considered as +1 , ( + 1)th guess, and we continue to step (g). g) At the end of this procedure, we test the convergence criterion (norm-2 relative error for phase-field, deformation history, and displacement variables field < ). If it is satisfied, we use the results as the initial guess for the next time step ( + 1)th. Otherwise, we use them as a new guess for the current th time step.

4.
Review of the experimental results in [4] For testing the different fracture codes discussed above in terms of their performance on simulating dynamic brittle fracture in glass, we choose the experimental results provided in [4] (see Fig. 6.a), which offers detailed information on the crack branching phenomenon, induced, in this case by impact. The crack propagation speed was monitored. In this experiment, impact on a notch with a pre-crack leads to a crack propagating straight before it branches. Interestingly, just before the two crack branches reach the far end of the sample, some small wiggles are noticed. It will be interesting to see if the computational methods tested in this work do capture such fine dynamic crack path details. The experiment uses the Digital Gradient Sensing (DGS) technique in conjunction with ultrahigh-speed photography to monitor the dynamic crack growth and branching in soda-lime glass. A Hopkinson pressure bar was used for loading a commercially procured 150 mm × 100 mm × 4.65 mm soda-lime glass plate with a 40° V-notch (see Fig.   6.b). The 19 mm deep V-notch was extended by 8 mm using a circular saw of 300 µm thickness. The impulse load on the specimen was generated using a gas-gun launcher aligned co-axially with the longbar. Crack velocities were estimated by the backward difference method from 32 images of crack-tip positions.
Experimental results of three fractured specimens are provided in [4] and are reproduced here in Fig. 7 for convenience. Similar to the observations reported in [82,83] on glass fracture, crack velocities reach maximum values between 1450 m/s and 1550 m/s. The crack-tip velocity was evaluated in [4] within ±100 m/s and ±150 m/s prior to and after branching, respectively. In all three tests presented in [4], crack branching was observed beyond the middle of the specimens (at 53%, 57%, and 58% of the width, measured from the left side of the sample), and occurred between 18-21 µs after the start of crack propagation. The time of the measurement in the experiment is set by the time the striker hit the long bar and a preset delay in the recording camera. The experimental results for crack propagation speed are reported starting from the moment the crack starts to grow. The branching angle were not reported in [4]. We used the photographs shown in Fig. 7 to measure the branching angles at a distance of 2mm (shown by two red lines) away from the estimated branching point. We estimated the branching angles to be approximately 57°, 57°, 55° for three specimens. These estimates are sensitive (±10%) to the placement of the branching point.
Reflecting waves from the glass sides interacting with the advancing cracks post-branching are the likely cause for the small twists/wiggles in the crack paths just before the two crack branches reach the far end of the sample (at 96% of the plate width).

Geometry and loading conditions setup
We simulate the sample geometry shown in Fig. 6. The pressure loading on the V-notch surfaces was not measured in the experiment. Instead, the impactor velocity was provided. Rather than simulating the impactor itself in the computational tests below, we choose to estimate the notch pressure loading in order to reduce the computational cost of the full model (which would require modeling the impact bar).
We have performed a similar calculation for a similar experiment on PMMA reported in [84] using a full LS-DYNA model of the impact bar hitting the notch, to estimate the pressure profile on the V-notch (see [85]) up to fracture initiation. Here we follow the same approach and consider the pressure profile shown  Fig. 8.a for all computational tests below. We apply the pressure profile on the V-notch surfaces as shown in Fig. 8.b to match the actual loading conditions on the notch (see comment in the previous section). Note that if the full-length of the V-notch is loaded, the obtained failure patterns (see Appendix A) differ from the ones produced by the partial loading of the notch surfaces (which is the actual one), pointing to the importance wave propagation in the sample have on fracture evolution. We consider a coefficient of friction of 0.35 between the glass and the long bar [86]. With this friction force, the angle between the normal to the surface and the total force is (20° + tan −1 (0.35)).

Fig. 8. (a) Pressure profile applied on the V-notch surfaces, (b) The pressure loading is applied over the portions in the V-notch in
contact with the impactor (see detail in Fig. 6.a).

Results from the peridynamic meshfree implementation
In general, the nonlocal length scale (the horizon) in PD modeling should be smaller than geometrical features to prevent unrealistic nonlocal interactions [64]. If the model does not have any geometrical features (e.g. notches), an appropriate horizon size would be one in the same range as the material lengthscale exhibited under particular loading conditions, if present [8,39]. We perform a δ-convergence study (horizon size decreases while keeping the number of nodes covered by a node's horizon roughly the same, see [87]) to see if crack propagation features (crack path, speed, angle of branching, location of branching) converge. We use horizon sizes of 2.0, 1.0, and 0.5 mm (given the thickness of the pre-crack) with a horizon factor (the ratio of horizon size to the discretization size) of five. The horizon size of 1.0 mm is in the same range as the "characteristic distance" for dynamic fracture in glass found in [2].
We convert the pressure profile shown in Fig. 8 to the body force components, and apply them on the PD nodes located on the V-notch surface as shown in Fig. 8b. We solve the 2-D bond-based PD model with the rest of the surfaces free. PD simulation results (damage maps and crack velocities) are shown in Fig. 9 and Fig. 10. These results do not vary significantly with changes in the horizon. The damage map in Fig.   9.b (when the horizon is 1 mm) shows a remarkable match with the experimentally observed crack path.
This match may be related to the closeness between the nonlocal region (the horizon size) and the characteristic length in dynamic brittle fracture of glass [2,88,89]. Results obtained with a horizon of 1 mm show that branching happens at around 28.5 μs from the moment the load is applied to the notch surfaces in the simulation, and at location 0.57W away from the left edge of the sample. With the same scheme of measuring the branching angle as in Fig. 7a, we find a branching angle of around 45°. Moreover, small twists/wiggles in the crack branches paths are predicted to happen at 0.96W from the left edge of the sample (see Fig. 9.b). This detail is observed in two of the three experimental samples as well (see Fig.   7a).

Results from LS-DYNA's discontinuous-Galerkin PD implementation
We setup two 3D finite element (FE) models with the non-uniform discretization with average element Damage maps for two mesh sizes are shown in Fig. 11. The location of the crack tips is obtained using the post-processing of the damage maps figures. A first branching event occurs around 29 μs from the moment the load is applied, at a location 0.5W from the left edge of the plate. The branching angle are roughly 40° and 37° for horizon sizes of 4.5 mm and 1.5 mm, respectively. We measure the branching angle using the same approach as shown in Fig. 7. A secondary branching event, not seen in the corresponding experiments, happens at a location around 0.82W from the left edge of the plate. The crack tip velocities are calculated using the backward finite difference and are shown in Fig. 12.

Results from the phase-field simulations
We perform PFM simulations using three length scales of 0.6 mm, 1.2 mm, and 2.4 mm. We model the pre-crack as an induced crack in the phase-field. To this aim, we prescribe an initial strain-history field such that a pre-crack in the phase-field is introduced [20]. The strain-history field is initialized as: where is a line representing the initial crack and ( , ) is the distance from point to the line . is a pre-defined large scalar corresponding to the phase-felid values close to 1 (fully broken state) to artificially generate the initial crack along line .
We show the damage maps with a fixed ratio of 2 of length scale parameter ( 0 ) to mesh size (∆ ) (for the three length scales above) in Fig. 13. Smaller length scale narrows the damage regions. For the model W with the smallest length scale, the crack branches around 29 μs and at a location around 0.57w from the left edge, with a branching angle of 28°. Since the damage is relatively thick here, we measure the branching angle by using the undamaged triangular region between the crack branches (see Fig. 13.c).
The velocity of crack tip (see Fig. 15) is calculated using post-processing the location on the rightmost node with damage of = 0.75.

Discussion
From the start, we should emphasize that uncertainty exists in the experimental tests in terms of the precise impact loading conditions. The partial damage created asymmetrically on the notch surfaces and elsewhere in the sample (see Fig. 7a), is indicative that it might not be realistic to expect a perfectly symmetrical impact. Nevertheless, some crack system features appear to persist between different experimental tests, and these features, most likely, are not overly sensitive to small perturbations in the loading conditions. Therefore, we expect to recover such features from our numerical simulations. Note also that our computational samples are load-free except for the impact loading region. This is different in the experiments, where the samples are held in place by putty layer, which changes slightly the way elastic waves are reflected from those boundaries (compared to the free edges case), implicitly influencing subsequent dynamic crack propagation.
Damage maps obtained by each of the three numerical models show a convergent behavior as the length-scale used in each model is decreased, but they converge to different results. The meshfree PD implementation is the only one that matches the experimentally observed crack patterns well. The histories of crack tip velocities obtained using the smallest length scales for each method are compared in Fig. 15, where we plot them along the data from experiments. All models find a steadily rising crack tip velocity up to a plateauing level, with some small dip(s) before branching takes place, and staying below 0.6 of the Rayleigh waves speeds in glass before and after branching. All three numerical models show crack branching at around 29 ±0.5 s from the moment the loading is applied. Cracks branched at around 30.7 ±1.5 s in experiments, when we match the start of crack propagation with the meshfree PD results.
We conclude, therefore, that the estimated time for crack branching is within the experimental range for all fracture models used here. When a critical stress intensity factor is reached, the crack splits into two or more branches, each propagating with about the same speed as the parent crack but with a smaller process zone [66]. Experiment ran a test with the meshfree implementation of PD using a horizon factor of 3.0 (see Fig. 16); no secondary branching is noticed and the crack pattern (location of branching, angle of branching) is very close to the ones shown in Fig. 9 and observed experimentally. Therefore, something else must cause the secondary branching seen in the LS-DYNA implementation results. We suspect that the way DG discretization models damage evolution might be responsible for this behavior. The crack branching angles of 45 ∘ , 37, and 28 ∘ were obtained with the meshfree PD model, LS-DYNA's DG implementation of PD, and PFM, respectively (see Fig. 17). A crack branching angle of around 55° is measured by using the pictures reported in [4]. Notice similar large discrepancy in the crack branching angles given by the PFM in comparison with experimentally measured values has been reported before in the literature ( [66]). One possible reason for this behavior is the "damage leakage" or "damage bleeding" notable in phase-field simulations of dynamic brittle fracture. After a crack passes a certain location, damage continues to advance into the material, in the direction perpendicular to the crack surface (see Movies 1 and 2). This extension of damage state is non-physical, the observed cracks are sharp, thin (see Fig. 7a). Once damage bleeds into the material, waves propagating through the material reflect differently than they should, leading to a different landscape of waves travelling through the body, reflecting from the computational sample boundaries (external or newly created crack surfaces). The newly created fracture surfaces grow more and more different from the actual crack surfaces of the real sample. Since wave reinforcements (between waves induced upon impact, those generated by the fracture process itself, and those reflecting from boundaries) are a deciding factor influencing the evolution of dynamic brittle cracks, this could explain the discrepancy between the crack paths obtained with the PFM and the experimental results. Notice that this computational phenomenon (unstoppable "thickening" of damage around crack lines) can be observed in many phase-field simulations of dynamic brittle fracture published in the literature (e.g. [90,91], see Fig. 18 For the quasi-static case, other forms of degradation functions in PFM have been recently proposed to increase the accuracy in predicting the critical loads associated with crack nucleation and preserving the linear elastic response in the bulk material prior to fracture [94]. These degradation functions avoid a noticeable localization of the phase field before the nucleation of a crack [95][96][97] or approximate a cohesive-type response for materials that exhibit some ductility after reaching their ultimate strength [98]. We have not investigated whether these newer PF models improve upon the results shown here.
We note that the PFM and the meshfree PD model predict the same branching locations at around 0.57W from the left edge of the sample, close to what experiments show. The PD implementation in LS-DYNA gives the branching location at 0.5W.
The notable crack path turns, as they approach the right side of the sample, are intriguing features of this impact experiment (see Fig. 7). Only the meshfree implementation of PD was able to reproduce these features (see Movie 3). LS-DYNA's PD implementation shows secondary branching instead, (see Fig. 11 and Movie 4), while the phase-field model shows slight undulations in the crack paths before these reach the boundary of the sample (see Movies 1 and 2). The two branches in the meshfree PD solution stop at 55 μs, arrest for 15 μs, then turn and continue growing until they reach the right side of the sample (see Movie 2). The twists/wiggles in the crack paths are likely caused by waves interacting with the advancing cracks. The influence of elastic stress waves bouncing from the boundaries on the behavior of running cracks in brittle materials has also been discussed in [99]. This important aspect shows why, in trying to predict dynamic brittle fracture, it is essential to model the entire sample geometry.
In terms of efficiency, we cannot make a fair comparison between the methods because we ran the three different models on different machines. We only report the computational cost of each method for each of the hardware used. The in-house GPU-enabled meshfree PD code ran on a single CPU with an

Conclusions
We evaluated the performance of three different computational models on dynamic brittle fracture in soda-lime glass induced by impact. We compared the results with some recently published experimental data on crack branching in glass. The three different models (an in-house meshfree implementation of a peridynamic model, a discontinuous-Galerkin implementation of peridynamics available in the commercial software LS-DYNA, and a phase-field model implemented in COMSOL) show some important differences in the fracture behavior they produce. For the angle of crack branching, the meshfree PD gives the closest value to the experimentally measured ones, while the phase-field model shows an angle size half of the observed value. The likely reason for the discrepancy between the phase-field model and experimental results is the "bleeding" of damage into the crack banks after the passing of the crack. The results support the claim that dynamic brittle fracture can be fully understood by using elasticity and a model capable of representing Griffith fracture. They also show that not all discretizations of peridynamic models are able to allow for unrestricted damage evolution, which is a necessary condition in dynamic brittle fracture, given the complexity of possible crack paths. While in quasi-static problems, phase-field models have been shown quite capable of modeling brittle fracture, the same cannot be said about dynamic brittle fracture problems. It remains to be seen if future developments can resolve this issue.

Acknowledgment
The

Appendix A
In this appendix, we show the effect of loading region. We apply the same pressure profile (see Fig. 8.a) on the entire length of the V-notch. As expected, different loading generates different waves in the material. The crack branches later (0.6W) with similar branching angle (see Fig. 19) compared with original loading (see Fig. 9.b). However, here the crack does not rest near the edge and consequently we do not observe wiggles near the edge of the glass (see Movie 5). The change in the load region will change the shape of advancing waves that eventually alter the crack branching location and crack behavior near the edge.