Peridynamic analysis of dynamic fracture: influence of peridynamic horizon, dimensionality and specimen size

In peridynamic models for fracture, the dissipated fracture energy is regularized over a non-local region denoted as the peridynamic horizon. This paper investigates the influence of this parameter on the dynamic fracture process in brittle solids, using two as well as three dimensional simulations of dynamic fracture propagation in a notched plate for two loading cases. The predicted crack speed for the various scenarios of the initially stored energy, also known as the velocity toughening behavior as well as characteristics of the crack surface topology obtained in different crack propagation regimes in 3D computational simulations are compared with the experimentally observed crack velocity and fracture surfaces for Polymethyl Methacrylate (PMMA) specimens. In addition, we investigate the influence of the specimen size on the dynamic fracture process using two dimensional peridynamic simulations. The fracture strengths and the velocity toughening relationship obtained from different specimen sizes are compared with the Linear Elastic Fracture Mechanics (LEFM) size effect relationship and with results from experiments, respectively.


Introduction
During dynamic crack propagation, the stress field around the crack tip, the energy absorbed in the fracture process, as well as the topology of the fracture surface changes as a function of the crack velocity. There are bifurcation points in the stress fields at certain levels of crack propagation velocities after which velocity toughening mechanisms, such as microbranching get activated and cause an increased dissipation rate due to additional surface creation [76]. In the following we provide a brief overview of the analytical, experimental as well as computational analyses performed to characterize dynamic crack propagation. The contradictions between the

LEFM and experiments
Linear elastic Fracture mechanics (LEFM) is an established theoretical framework to characterize and predict the initiation as well as the subsequent motion (position and velocity) of a moving crack, which started with the pioneering work of Inglis [29] and Griffith [43]. LEFM predicts that cracks can not propagate arbitrarily fast, there are speed limits that they have to obey. The investigation of the stress field around a dynamically moving brittle crack of fixed length performed by Yoffe [87] suggested an explanation for this speed limit based on the bifurcation of the dynamic stress state around the moving crack tip. She showed that, as the crack speed increases to 60% of the shear wave speed, the maximum circumferential stress shifts from 0 • (i.e. the plane in front of the crack) to ±60 • . However, the resulting dynamic stress intensity factor from Yoffe's solution was independent of the crack tip velocity, which is incorrect. This was later corrected by Broberg [12], who considered a crack that starts propagating at constant velocity from a zero initial length. The resulting ratio of dynamic to static stress intensity factor for the same crack length from Broberg's analysis decreases with increasing crack velocity and vanishes at the Rayleigh wave speed v r , which is the theoretical velocity limit for a mode-I crack. This result was later confirmed by Atkinson and Eshelby [4]. Afterwards, Freund [40,41] analyzed dynamic crack propagation with a constant as well as a non-constant velocity while considering the interactions between the stress waves and cracks. The dynamic stress intensity factor obtained by Freund [40] was in agreement with [12] and [4], i.e. the terminal velocity for a mode-I crack is v r .
One of the contradictions between the predictions of LEFM [40] and the observations from experiments [20,70,76] is related to the terminal velocity of cracks. According to LEFM, a crack should continuously accelerate up to the Rayleigh wave speed v r . However, in experiments triggering a mode-I crack, this value is never realized. The limiting crack velocity observed in experiments for soda-lime glass is about 0.51v r [20], around 0.7v r for PMMA [76] and about 0.45v r for Homalite-100 [70]. This discrepancy is attributed to the assumption of a single sharp crack in the theoretical analyses. In experiments, however a single sharp crack is observed only up to a certain velocity and after that, an increased roughening of the fracture surface accompanied by an increase in the dissipation rate is detected. LEFM predictions [40] for brittle dynamic fracture are shown to be valid [75], as long as the assumption of a single propagating crack is replicated in the experiments. From postmortem fracture surface examinations [39,67,69], it has been shown, that the surface roughness of a propagating crack is an increasing function of the crack tip velocity. This increasing fracture surface roughness is attributed to the different micro-structural processes occurring at the crack tip as well as in the process zone in front of it. Ravi-Chandar and Knauss [69] attributed this fracture surface roughness to three main sources: the material heterogeneity (voids or inclusions), the material micro-structure and the interaction of the stress waves with the propagating crack (leading to the so-called "Wallner lines" [67]). Fineberg et al. [38,39,74] reported a well defined transition in the motion of the crack tip, from a smooth propagation to an oscillatory propagation above a critical velocity. They proposed the existence of a dynamic instability (micro-branching instability) above a critical crack tip velocity, which governs the crack tip oscillations and showed that these oscillations were correlated to periodic micro-structures observed on the fracture surface.

Computational analysis
A number of different classes of continuum models have been put forward to bridge the gap between dynamic fracture experiments and computational simulations. These continuum models are mostly governed by the partial differential equations of continuum mechanics. These partial differential equations are solved using approximation functions and their derivatives, which are constructed either locally using a mesh or using scattered points without any mesh connectivity. The first class of methods is known as mesh-based methods, such as Finite Element Method (FEM), while the second is known as mesh-free methods, such as the Diffuse Element Method (DEM) [61] or the Element Free Galerkin (EFG) method [9]. For a review of mesh-free methods we refer to [31,48,73]. To represent the discontinuity resulting from a crack, discrete as well as smeared representations of cracks as strong and weak discontinuities are employed in these models.
Cohesive zone models (CZM) [28,86] equip interfaces of standard finite elements with a cohesive law to model crack propagation. These models provide autonomous initiation, propagation and branching of a crack [86], however, the crack path is highly dependent on the particular topology of the discretization [37,92]. A framework for interface elements using a variational (phase-field like) formulation was recently presented in [52]. In [86] some features of fast crack propagation were captured, such as crack branching and the existence of limit velocity for cracks. Impact damage in brittle materials was modeled in [28] using a linearly decaying cohesive relation, discussing the existence of an intrinsic time scale related to the critical crack opening displacement and the elastic wave velocity even when using a rate-independent cohesive relation. More recently, [93] presented a rate-dependent cohesive law and compared the crack propagation velocities with a rate-independent cohesive relation which over-estimated the crack propagation velocities. Some critical artifacts of dynamic fracture simulations using cohesive zone models are also presented in [37].
FEM in conjunction with cohesive elements experiences mesh bias when representing a discrete crack (a strong discontinuity). This mesh-dependence can be addressed by using adaptive re-meshing techniques which, however, become particularly challenging in 3D. Circumventing the use of a mesh in mesh-free methods facilitates modeling arbitrary cracks without re-meshing [10,62]. A second option for having a discrete representation of cracks without the need for re-meshing is accomplished by the eXtended FEM (XFEM) [8,60], which allows the crack to propagate through the finite elements, employing enrichments with a strong displacement discontinuity kinematics. This removes the mesh bias and allows the crack to propagate in any direction. However, it requires additional degrees of freedom and suitable crack indicators for the advancement of the crack tip as well as crack tracking which adds additional complexity [68]. In particular, additional criteria for crack branching, coalescence and arrest are required. A loss of hyperbolicity criterion was presented in [11] to predict the onset of branching.
For modeling complex fragmentation processes in brittle materials it is desirable to represent crack propagation independently from the discretization with minimal or no external intervention such as tracking of discontinuities or decisions on the direction of crack growth. A number of numerical models for fracture simulations rely on non-local averaging approaches to regularize the surface energy of the discrete crack over a length in order to avoid mesh dependence. The non-local approaches generally provide an autonomous crack initiation, growth and branching and one does not need to track the crack tip. Some models in this category include the non-local damage model [66], gradient enhanced damage models [65], eigen-erosion models [63], phase-field models [19] and peridynamics [77]. The similarities and differences between the gradient damage and phase-field damage models are analyzed in [34]. Phase-field models for dynamic fracture [18,46] regularize the discontinuity over a characteristic length, which defines the width of the failure zone. Hence, the representation of the discontinuity is not discrete anymore. In [14], micro-branching instability was modeled using a phasefield model. It was reported that the length and the frequency of the micro-branches depends on the phase-field internal length scale. Some of the features of dynamic fracture propagation, such as a limiting crack speed, toughening of the material with increasing crack speed and crack branching, were qualitatively reproduced in [15]. The effect of the internal length scale of the phase-field on the crack propagation velocity and dissipated energy is discussed in [18]. However, the mechanisms which govern the change of the dynamics of cracks with respect to this internal length scale still do not seem to be well understood.
Besides continuum models, atomistic simulations [1,2] have been used to investigate dynamic fracture processes. In [22], it was argued that the large deformations at the crack-tip leads to the transition of the material around the crack tip to a hyper-elastic zone. This hyper-elastic zone sets a local value of the Rayleigh wave velocity around the crack tip which limits the crack propagation velocity.

Dynamic fracture modeling using peridynamics
Peridynamics [77,81] is a non-local re-formulation of partial differential equations of classical continuum mechanics which provides continuum model capabilities to incorporate discontinuities into kinematic fields (i.e. no additional enrichment functions or crack tracking is required). The peridynamic formulation employs spatial integrals to compute the strain energy at a point instead of the spatial derivatives of the classical continuum, which circumvents the problem of computing gradients over discontinuities. The dissipated fracture energy is regularized in peridynamics over a non-local region, known as the horizon, similar to classical continuum methods using non-local regularization as discussed in Sect. 1.2. Unlike other non-local continuum models for fracture, however, peridynamics introduces damage at the level of connections (or bonds) between the material points.
The ("sub-continuum"-level) bond breakage results in the evolution of a discrete fracture surface autonomously, so multiple discrete cracks can initiate, propagate and interact with each other. Fracture modeling with peridynamics combines the properties of a regularized damage approach and a discrete fracture approach in one model. This non-local regularization introduces a new parameter in the formulation, the radius of the peridynamic horizon. The peridynamic horizon, just like the regularization length in other classical continuum methods, can be viewed as a consequence of the mathematical formulation. A clear argument from the physics of the problem for the need of such a parameter does not exist. However, as will be shown later in this paper, it does affect some aspects of the physics of the problem, especially when inertial forces are in play, such as wave propagation and dynamic crack propagation. In some of the non-local regularization methods used in conjunction with (local) classical continuum, this regularization length has a quantitative affect on the quasi-static fracture behavior. An example is the regularization length in phase-field models, which affects the limit tensile stress at which crack nucleation occurs [18]. It will be shown in the current contribution, that this is not the case in peridynamics. The role of the horizon in static analyses and in elastic wave propagation analysis assuming a peridynamic continuum is relatively well explained [6,25,26,77]. However, the effect of the peridynamic horizon on dynamic fracture processes is not fully understood.
A mesh-free method based on bond-based peridynamic continuum presented in [79] was able to reproduce dynamic fracture and fragmentation of brittle solids subjected to impact loads. Bond-based peridynamics utilizes central force potentials between material points to characterize the force interactions. This leads, however, to a fixed Poisson's ratio of 1/4 for 3D 1/3 for 2D problems. State-based peridynamics [81] does not suffer from this limitation. It can also incorporate constitutive relations from the classical continuum theory (i.e. constitutive relations involving a deformation gradient) using the correspondence framework in peridynamics [81]. However, the correspondence framework suffers from zero energy mode instability [7,78], stabilization conditions and stabilizing terms to the strain energy density have been proposed as remedies in [78,83,84]. Bond-based peridynamics uses a critical-stretch based criterion for bond failure [77]. It is noted, that for state-based peridynamics, a critical-energy based failure criterion, which takes into account the total bond energy (i.e. isotropic as well as the deviatoric part of the deformation) was introduced in [49]. A comparison between the critical-stretch and critical-energy failure criteria is presented in [35,91].
Dynamic crack propagation and branching in brittle materials, using a bond-based peridynamic model with criticalstretch failure criterion, was investigated in [88,89]. Crack branching and overall fracture patterns were shown to be in qualitative agreement with the experiments. Using a convergence study with respect to the horizon size, they showed, that the crack velocity profile stays below the Rayleigh wave speed and that it converges to the maximum crack velocity observed in experiments, once the horizon reaches sub-millimeter values. A peridynamics model was also used to analyze the dynamic fracture of shale and rock material in [24,32]. In [3], a comparative study between cohesive zone models, the XFEM and peridynamics was provided and the crack propagation patterns and velocities were compared. The crack velocities obtained using all three methods were overestimated compared to the experiments, but remained below the Rayleigh wave velocity. Dynamic brittle fracture and crack branching was investigated in [16] using bond-based peridynamic simulations performed using various transient boundary conditions. A stress wave pile-up mechanism at the crack tip, which deflects the damage away from the symmetry line in front of a mode-I crack, was investigated as the possible reason for crack branching.

Objectives and organization
The aim of this study is to investigate the dynamic fracture behavior using an ordinary state-based peridynamics model [81] along with a critical-energy based criterion for bond failure [49]. We perform two as well as three dimensional simulations of dynamic crack propagation. The motivation of performing these analyses in two and three dimensions is twofold. Firstly, the peridynamic method suffers from the well-known surface effects [59] caused by the fact, that the material points at the boundaries of the solution domain do not have complete nonlocal neighborhoods. This influences three dimensional simulations differently than it does for the two dimensional case. And secondly, the mechanisms involved in dynamic crack propagation, such as micro-branching, are shown to be dependent on the specimen thickness [14]. Performing both 2D and 3D simulations enables the investigation of the influence of dimensionality on the dynamic crack propagation. With this background, the paper investigates the role of the peridynamic horizon in simulations of dynamic crack propagation in plates made of PMMA for different levels of initial quasi-static loading of the specimen before dynamic crack propagation initiates. Considering four different values for the radius of the peridynamic horizon, the influence of the peridynamic horizon on the dynamic fracture is investigated using delta-convergence analysis, in regards to crack propagation patterns and velocities, the energy threshold for branching and the evolution of kinetic, elastic and dissipative portions of the total energy.
To the best knowledge of the authors, a simulation tool which is able to reproduce and validate the relationship between the crack velocity and the fracture energy release rate observed in brittle dynamic fracture experiments, without incorporating an a priori phenomenological crack velocity or strain rate dependency, does not exist. Such phenomenological rate-dependent criteria have been used to influence the crack velocity in conjunction with interface elements (CZM), (see, e.g. [93]), non-local continuum damage models [85], phase-field models [36] as well as peridynamics [27]. In this study we do not resort to a phenomenological rate-dependent criterion; instead we utilize a simple criticalenergy based damage criterion. Hence, the dependence of the dissipated fracture energy on the crack propagation velocity is obtained purely from the interaction of the crack tip with the simulation domain.
The remainder of the paper is organized as follows: Sect. 2 provides a brief introduction into the ordinary state-based peridynamic formulation and the critical-energy based damage criterion. The simulation setup and the loading conditions are explained in Sect. 3. In Sects. 4 and 5, two and three dimensional simulations of crack propagation in a PMMA plate are presented. Both sections discuss the sensitivity of the fracture process with respect to the peridynamic horizon size for a monotonic loading case as well as a spectrum of quasi-static loading cases. The effect of the peridynamic horizon on the crack propagation speed, the evolution and balance of the elastic, the kinetic and the dissipated energy and the overall dissipated fracture energy -crack velocity relationship as well as the threshold for branching are discussed. The dissipated fracture energy -crack propagation velocity relationship obtained from two as well as three dimensional simulations is compared with experiments. The qualitative differences between the results from two and three dimensional crack propagation simulations and the influence of the specimen size on the fracture process are discussed in Sect. 6 along with comparisons with experimental results. Finally, the findings from the peridynamic simulations are summarized in the concluding remarks (Sect. 7).

Peridynamic formulation
The peridynamic formulation of a continuum is characterized by the direct interaction of a material point x with a set of material points x in a volume defined by a cut-off radius δ, known as the peridynamic horizon H x [81]. Equilibrium at a material point x at time t, subjected to a force density b(x, t), is given by the equation of motion x − x is the force, which a material point x exerts on x . Angular brackets are used to denote quantities that T operates on.

Constitutive modeling for elastic material response
The The scalar force state t[x, t] depends on a scalar stretch-like quantity, denoted as the extension state e[x, t], which characterizes the kinematics of the model. Extension is defined as the difference of the length of the deformed and undeformed relative position vectors, i.e. e ξ ξ ξ = |Y ξ ξ ξ | − |ξ ξ ξ |, where we use x − x = ξ ξ ξ . The extension state e can be further decomposed additively into an isotropic (e i ) and a deviatoric (e d ) extension state. The isotropic extension state can be represented in terms of a scalar-valued volume dilatation θ(x) that is defined to match the volumetric strain of a classical continuum model under isotropic loading conditions. According to [81], θ(x) is defined as: where ω(ξ ξ ξ) is the influence function and m is the weighted volume. The material parameters corresponding to e i and e d are determined by strain energy equivalence with the classical continuum model [81]. This leads to the force state t[x, t] for three dimensional problems given by: where, K and μ are the bulk and shear moduli, respectively. For two dimensional LPS models, we follow the formulation presented in Chapter 6 in [17]. 2D equation of motion is obtained by replacing the volume dV x in Eq. 1 by t 2D d A x , where t 2D is the plane-stress thickness and d A x is the area of the material point x . Dilatation for the planestress case is obtained by scaling θ(x) from Eq. 3 by γ , where

Surface effects
As mentioned earlier, peridynamic material models are defined by strain energy equivalence with classical continuum models via the definition of a peridynamic neighborhood H x . However, as material points located near a boundary do not have a complete non-local neighborhood, the effective material properties near the surface of a peridynamic model are computed to be slightly different from those in the bulk. A number of methods have been recently proposed for correcting this peridynamic surface effect (see for e.g. [55,59]). These correction procedures specify different material properties for the affected points at the boundaries. This is usually achieved by comparing the differences in the internal force for simple loading conditions obtained from peridynamics and classical continuum. The accuracy of such techniques depends on the specific geometry and boundary conditions. Additionally, it is also unclear how these corrections influence the total energy in the system. A consistent solution for these surface effects, i.e. independent of the problem type and geometry as well as without influencing the total energy of the system, is provided by the correspondence model [81]. However, correspondence formulation suffers from zero energy modes [7,78] and stabilization techniques [84] again introduce some fictitious parameter which influences the total energy. Recently, the idea of bond-associated correspondence has been put forward [21,30,57]. These models have enhanced stability along with the benefits of the original (nodal) correspondence model [81], but they do so at an extremely high computational cost. In future research, we plan to explore the applicability of the bond-associated correspondence formulation to model dynamic fracture and fragmentation.

Constitutive modeling for fracture
Crack initiation and propagation is modeled in peridynamics by irreversibly breaking the connections (bonds) between the material points. This results in an autonomous crack propagation at the continuum level without the explicit need of a criterion for the direction or for the length of the propagating crack. A critical energy based bond failure criterion [49] is used to define a critical threshold for the bonds to fail irreversibly. The stored energy density in a bond w ξ can be calculated by projecting the relative force vector on the relative displacement vector as [49]: The integrand of w ξ is called the dual force density and w ξ represents the energy density in a bond due to a relative displacement of two material points x and x from η η η = 0 at time (t = 0) to η η η(t f inal ). The energy density in a bond w ξ is fully recoverable unless it exceeds some critical energy density level w c . w c was defined in [49] in terms of the fracture energy G c , a material property that can be determined experimentally. Figure 1b shows a fracture plane separating the two halves of a body. The bonds between all material points A (along the dashed line 0 < z < δ) and all material points B in the spherical volume of radius δ are bridging a unit fracture surface. It is assumed that, when all these bonds exceed w c , a fracture surface of unit area is created which results in an energy release of G c : The integration domain of Eq. 7 is shown in Fig. 1b. For further details we refer to [44,49]. Eq. 7 is solved for w c , leading to For the two dimensional case, following [35], the energy required to break all the bonds per unit fracture area is given by: which leads to the critical energy density level w c 2D as: In Eqs. 8 and 10, δ acts as a localization limiter and guarantees a constant energy release rate independent of the discretization. Further details are provided in [35]. Both failure criteria were implemented in the open source software Peridigm [54,64]. In a structural simulation, once the condition w ξ ≥ w c is satisfied for a bond, it breaks, i.e. the stiffness of the bond vanishes irreversibly and it does not contribute to the internal force any more. A damage variable d can be defined as the ratio of the number of broken bonds to the total (initial) number of bonds at a node. This damage parameter d will be used for generating contour plots along the computed crack path in the next sections.

Discretization
The widely used peridynamic implementation to date is the mesh-free discretization technique presented in [79]. It uses a set of nodes, where each node is associated with a cell of known volume in the reference configuration. This discretization scheme results in a strong-form collocation method with a relatively low computational expense, but leads to only a first order convergence rate [56,72]. In [82], it was shown that a higher order convergence rate can be obtained using quadrature based finite difference discretizations or by employing finite element discretizations. A Reproducing Kernel (RK) enhanced approximation approach was presented in [56] for bond-based peridynamics and in [45] for the correspondence peridynamics framework. It was shown, that the introduction of RK shape functions in the solution approximation allows for arbitrary smoothness and completeness of the approximation. These higher-order techniques, however, significantly increase the implementation complexity and the computational cost.
In this study, we investigate the ability of the frequently used mesh-free peridynamics framework [79] to model the velocity toughening relationship and the fracture behavior of PMMA for a spectrum of loadings. Higher order discretization techniques are not within the scope of this paper.
The mesh-free discretization of the correspondence model of peridynamics [81] is well investigated for uniform grids and similarities have been drawn with Smoothed Particle Hydrodynamics (SPH) discretization [42] as well as with RKPM discretization with synchronized derivatives [13,48]. The LPS model used in the current study is, however, different from the correspondence models, as the shape tensor (moment matrix in the context of other mesh-free methods) is not computed. For the general state-based peridynamics (i.e. without using the correspondence framework) and its discretizations, especially using non-uniform grids, similar analyses do not seem to exist. In this study, we have employed uniform grids only. Local mesh refinement techniques in the framework of dual horizon peridynamics are addressed in [47,71].

Simulation setup
The rectangular PMMA specimen used in the simulations has a length of 32 mm, a height of 16mm and a thickness of 0.5mm. The crack is supposed to initiate from an initial notch with a length of 3mm, as shown in Fig. 2a. The material parameters for PMMA used in all simulations are specified as follows: Young's modulus E = 3.09 × 10 9 Pa, Poisson's ratio ν = 0.35, Density ρ = 1180K g/m 3 , Fracture energy G c = 300J /m 2 and Rayleigh wave speed v r = 906m/s.
In peridynamic analyses, two kinds of convergence studies can be performed. One is δ− convergence [88], which is performed by keeping the nodal density n (number of discrete points along the radius of the peridynamic horizon) constant, while the horizon size δ is systematically decreased, i.e. δ → 0, the solution in this case approaches the classical continuum solution [80]. δ and n are related to each other by δ = n x, where x is the mesh size. The second is n− convergence [72], where the horizon δ is kept fixed and the nodal density n is increased. As n → ∞ the solution approaches the analytical solution of peridynamics for the particular horizon size. In this study, we focus on δ− convergence, considering four values for the peridynamic horizon (δ = 280, 210, 140, and 70μm) with an approximate mesh size x = 80, 60, 40, and 20μm, respectively. This results in n = 3.5 for all simulations. For n lower than 3, we have observed in some cases, that the cracks starts following the mesh, resulting in staircase fracture patterns. As very high values of n may become computationally expensive, we have chosen n = 3.5. Nonetheless, a n−convergence study will be part of future investigations. This discretization leads to around 0.08, 0.142, 0.32, and 1.28 million material points in two dimensions and 0.48, 1.138, 4.16, and 32.0 million material points in three dimensions for δ = 280, 210, 140, and 70μm, respectively. All simulations are carried out using an extended version of Peridigm [54,64].
In the experiments [38,39,76], the specimens were loaded incrementally every 10 − 20sec at very low strain rates in order to dampen the elastic waves induced by the loading. A similar loading condition in a full dynamic simulation applying such a slow loading procedure would be extremely time consuming. In order to eliminate the effect of elastic waves due to transient loads, we use a combination of an implicit quasi-static and an explicit dynamic solver. In the quasi-static part of the solution, the specimen boundaries are subjected to prescribed displacements such that a desired level of elastic strain energy per unit area is stored in front of the crack tip. Subsequently, the simulation switches to an explicit dynamic solver, while keeping the displacement at the specimen boundaries fixed, and the crack is allowed to propagate. Due to the nonlocal nature of peridynamics, the displacement boundary conditions are prescribed on a finite thickness layer at the boundary [88]. Traction boundary condition in the current formulation can be applied as an equivalent body force acting on this volume. For recent developments in consistent handling of the nonlocal traction boundary condition in peridynamics, we refer to [90]. In all simulations, we use a boundary layer of 240μm thickness to prescribe the displacement boundary condition. This thickness has been chosen such, that there are at least three nodes in the boundary layer for the coarsest simulation. For the explicit dynamic solver, a time step of 1.0 × 10 −9 s is used for all simulations in order to avoid any influence of the temporal discretization. A similar combination of solvers was also used in [15,23,93] for the same purpose.
In each simulation, a displacement U (see Fig. 2a) is applied in y direction on the top and bottom surfaces of the specimen quasi-statically (while keeping the crack stationary), i.e. ignoring the inertia effects in the simulation domain (see Fig. 2b). The specimen is unconstrained in x direction. This loading condition provides a constant value of stored elastic energy G per unit area prior to fracture in front of the initial notch: Once the crack propagates completely through the specimen, G would be approximately the total amount of energy dissipated in the fracture process. After achieving the desired value of elastic energy G by applying U according to Eq. 11, the analysis switches to an explicit dynamic solver without changing the boundary conditions and henceforth, the crack is allowed to propagate dynamically taking inertia effects into account. The crack propagation speed v c =ȧ is computed by tracking the crack tip every 100 solution steps, i.e. the computed crack speed is averaged over 100 time steps, where a(t) is the position of the crack tip. As simulations use a time step size of 1.0 × 10 −9 sec, the crack velocity is computed ten times in one μsec. The crack tip is defined as the most right node in the specimen experiencing a damage level larger than 0.4. In case of macroscopic crack branching, the faster propagating branch is tracked for determination of the crack velocity.
The rate of energy dissipated per unit crack area released during fracture propagation is computed from: as E el (t), E kin (t) and E diss (t) denote the total elastic, kinetic and the dissipated energy as a function of time. C is the initially (quasi-statically) stored elastic strain energy and a(t) is the area created by the propagating crack as a function of time. a is computed by assuming that a single crack is propagating through the medium, i.e. without taking the additional area created by the micro-branches into account. Hence, an increase in is expected when the micro-and macrobranches start to initiate and propagate (crack roughening), see Fig. 4. consider monotonically increasing loads to investigate the influence of the peridynamic horizon on the (quasi-static) fracture strength and the effective stiffness. In a subsequent series of simulations, we use a combination of an implicit quasi-static and an explicit dynamic solver to load the plates to a desired level prior to fracture initiation. Comparisons between cracking patterns, crack speed, evolution of energy dissipation rate as well as the elastic and kinetic energy at different loadings, obtained for varying horizon sizes, are presented.

Crack propagation under monotonically increasing load
The motivation behind these simulations is to investigate the influence of the peridynamic horizon on the fracture strength and effective stiffness. These simulations are performed at a loading rate which gives a straight propagating crack without any branches. The horizontal boundaries of the PMMA plates (Fig. 2a) are loaded monotonically until a crack starts propagating from the initial notch. The displacement is applied at a rate of 3mm/s. The load-displacement relationship obtained at the top boundary (i.e. at y = 16mm) is shown in Fig. 3a.
It can be seen that the effective stiffness is not sensitive with respect to the horizon size. The same holds for the fracture strength, which provides almost identical results for δ = 140 and 70μm and a marginal reduction of 4% for the largest horizon size 280μm. The figure confirms, that by virtue of the determination of the energy density in a bond according to Eq. 10, the fracture energy G c is preserved independently from the horizon size. Figure 3b shows the temporal evolution of the crack tip position. Due to slight differences in the fracture strengths, the cracks initiated at slightly different times in the four specimens. However, the crack propagation velocity is almost same (i.e. v c ≈ 538m/s) for all four sizes of the peridynamic horizons.

Dynamic crack propagation for different levels of quasi-static loadings
In these simulations, the specimens are loaded quasistatically up to different levels while restricting the initial notch (shown in Fig. 2a) from propagating. After reaching the desired loading value, the simulations switch from an implicit to an explicit dynamic solver and the crack is allowed to propagate. In this subsection, we investigate the influence of the peridynamic horizon on the cracking pattern, the crack tip velocity, dissipation rate, energy balance and the crack velocity toughening mechanism. This comparison is performed for four loading levels U = 21μm, U = 25μm, U = 30μm, and U = 40μm. These loading levels are related to four levels of elastic prestressing of the plate, which lead to four different characteristics of crack propagation (smooth crack propagation, rough crack propagation, crack branching and multiple crack branching). Each of these cases has been simulated by peridynamics model using four levels of peridynamic horizons (δ=280 μm, 210 μm, 140 μm, 70 μm).   It should be noted that in our simulations U = 21μm is the smallest load at which cracks propagate, i.e. at U = 20μm cracks do not yet propagate in all specimens. With an increase in load ( U = 25μm), the crack path starts to exhibit slight roughness along the crack path (second row of Fig. 4). At U = 30μm (third row of Fig. 4), macroscopic branching attempts along the main crack can be  seen. And at U = 40μm (fourth row of Fig. 4), two stable daughter cracks propagate. The temporal evolution of the crack tip position and the normalized crack propagation velocity with respect to crack tip position, for all four considered loading cases in Fig. 4, are presented in Fig. 5. Crack propagation velocity is normalized with respect to the Rayleigh wave speed (v r = 906m/s).
For U = 21μm, the crack tip position is almost identical for all four horizon sizes, as can be seen in Fig. 5a. As the crack paths obtained for this loading value were smooth and featureless (first row of Fig. 4), this results in smooth crack propagation velocity with a few fluctuations between 25 to 30mm crack length, as shown in Fig. 5b. These few fluctuations can be attributed to a slight roughness along the end of the crack path, which is hard to see for δ = 20, 40 and 60 μm, but is visible for δ = 280μm. The maximum difference between the average crack propagation velocities obtained from the largest (δ = 280μm) and the smallest (δ = 70μm) horizon size is less than 1% of v r .
A rough crack path is obtained for U = 25μm, as shown in second row of Fig. 4. It can be observed, that the cracks start smooth and develop some roughness afterwards, connected with a slight deviation of the crack path from a straight horizontal line. This smooth crack propagation in the early stage corresponds to the smooth crack tip velocity up to 7mm in Fig. 5d. The maximum difference between the average crack propagation velocities in this rough crack propagation regime predicted by the different horizon sizes is around 6% of v r .
For the loading level U = 30μm, the main crack starts to develop branches along its complete path as presented in the third row of Fig. 4. These branches are not stable and they get arrested after a certain length. The frequency of these branches along the main crack increases with decreasing horizon size. The crack velocities plotted in Fig. 5f show high fluctuations from the beginning, which is a consequence of the non smooth crack propagation. The vertical lines in Fig. 5f indicate crack branching events. The maximum difference between the average crack propagation velocities for the investigated range of peridynamic horizons is around 3% of v r .
When the loading is increased to U = 40μm, two stable branches develop soon after propagation starts (fourth row of Fig. 4). Both stable branches develop further subbranches along their path. In this case, the crack tip position and velocity is tracked for the faster propagating branch in Fig. 5h. The crack velocity has very high fluctuations which is a consequence of the crack roughening and branching. The stable daughter branches are developed for all four horizon sizes at around same crack length except for δ = 140μm which branches at the initial notch. The initial branching angles for all horizon sizes are around 35 • . The maximum difference between the average crack propagation velocities is increasing with increase in U . For U = 40μm, it is around 7% of v r . It should be noted that, the difference between the crack propagation velocities with respect to the horizon size for straight and smooth crack propagation was negligible and this difference increased as the crack path started to roughen and branch.

Dissipation rate and energy balance
The influence of the peridynamic horizon on the evolution of elastic, kinetic and dissipate d energy (Eq. 12) as well as the instantaneous dissipation rate (Eq. 13) over crack length are discussed in this section. Figure 6 presents the evolution of the dissipation rate ( /G c ) and the elastic, kinetic and dissipated energies (normalized with respect to total energy C in Eq. 12) plotted over the crack tip position for the specimens shown in Fig. 4.
By comparing Fig. 6a, c, e, and g, it can be seen that the dissipation rate increases with increasing U . The average crack velocity for δ = 70μm increases from 496m/s (Fig.  5a) at U = 21μm to 561m/s (Fig. 5g) at U = 40μm. As the crack velocity increases, the normalized dissipation rate ( /G c ) increases from around 1.0 (Fig. 6a) to approximately 4.0 (Fig. 6g). This is caused by the additional surface created by crack roughening and branching at higher crack propagation velocities. The increase in dissipation with increasing crack velocity is known as crack velocity toughening, which ultimately sets the speed limit on propagating cracks. The influence of the horizon size on the dissipation for a spectrum of loading levels will be discussed in the next subsection.
A crack propagating smoothly should consume approximately G c amount of energy per unit area, which is the case for the loading level U = 21μm (first row of Fig. 4). As shown in Fig. 6a, the fracture patterns for all four horizon sizes are smooth, which leads to a monotonic growth of /G c up to 1.0 . As the load level is increased to U = 25μm, we observe a rough crack path, however, without yet developing macroscopic branches, as presented in the second row of Fig. 4. It can be seen that /G c in Fig. 6c oscillates between 1.0 and 2.0, but never exceeds 2.0. A crack which forms two new smooth branches would consume 2 G c amount of energy, i.e. the energy required to form two fracture surfaces of unit area. The case U = 25μm is still (sightly) below this state. /G c > 2.0 was also suggested in [15], using a phase-field damage model, as an energetic criterion for macroscopic branching.
For δU = 30 and 40μm, the crack patters show macroscopic branches for all four horizon sizes as shown in the third and fourth row of Fig. 4, respectively. Fig. 6e, g present the normalized dissipation rate for these two loading cases. Vertical lines in these figures correspond to the major branching events during crack propagation. In Fig. 6e ( U = 30μm), the dissipation rate ( /G c ) for all four horizon sizes is oscillating around 2.0. This is a case where one crack tip is not enough to dissipate all the energy and two equivalent macroscopic cracks (i.e. branching) would dissipate too much energy. In the case U = 30μm, the main crack forms short lived branches along the way, which relaxes the main crack tip. According to Fig. 6g, for the case U = 40μm, /G c stays well above 2.0. This indicates, that two stable branches develop, which again develop their own sub branches during propagation.It is noted, that the initial crack branching occurs at almost the same position of the crack tip for all four horizon sizes. Figure 6 shows, that the growth of the dissipated energy is accompanied by the growth of the kinetic energy during the crack propagation. This is because the evolution of the kinetic energy around the tip for a straight propagating crack depends directly on the crack speed. As a mode-I crack separates the material into two pieces during dynamic fracture, the crack faces move with a certain velocity (proportional the to crack propagation velocity) perpendicular to the direction of the propagating crack, this motion contributes to the total kinetic energy of the specimen. In the context of peridynamic simulations of dynamic fracture, this will be explained for the case of a straight and smooth propagating crack, without any branches, as obtained for the loading level U = 21μm. As the crack propagates, a damaged region is created along the crack faces. As the material is less stiff in this damaged zone as compared to the bulk material, a part of the released energy is transformed into kinetic energy in this region. This process is illustrated in Fig. 7 at three different positions of the crack tip for different horizon sizes. Only marginal differences of the size of the concentrated kinetic energy density zone around the crack tip is observed for this loading case, which is confirmed by the small influence of δ on the evolution of the total kinetic energy in Fig. 6b. Nonetheless, there are differences with respect to the wave dispersion properties which result in different patterns of trailing waves for each horizon size [25].
As the loading is increased to 30 and 40μm, according to Fig. 6f, h, slightly larger differences of the plots for the kinetic energy obtained from different horizon sizes can be observed. These differences can be attributed to the fact that there are different number of crack tips with respect to δ (at an arbitrary time instant during crack propagation), and every crack tip has a damaged region associated with it where kinetic energy gets concentrated. This slightly influences the evolution of the total kinetic energy of the specimen. Corresponding with the observations from the left column of Fig.  6, a slight influence of the horizon size on the evolution of the dissipated energy during crack growth is observed, which also manifests itself in the evolution of the decaying elastic energy.

Crack velocity toughening mechanism
In this section, we focus on the velocity toughening behavior, which governs the crack tip speed limit, and the influence of the peridynamic horizon size on it. The velocity toughening mechanism is characterized by the relationship between the energy stored prior to fracture G (from Eq. 11) and the crack propagation speed v c . With increasing G, the crack propagation velocity increases up to an asymptotic value which sets the crack tip speed limit. Crack propagation simulations are performed for all four horizon sizes at different load levels characterized by different initial energies G stored in the specimen (quasi-statically) before a crack starts propagating. Various levels of G are selected from Eq. 11 by choosing equal intervals of U , for G < 2.0K J/m 2 an interval of 5.0μm is selected and for G > 2.0K J/m 2 it is set to 10.0μm.
The average crack velocities v c , computed as explained in Section 3, and normalized with respect to v r for the selected levels of G and for all considered horizon sizes, are plotted in Fig. 8. It is noted, that at higher load levels, and for larger horizon sizes, the cracks simultaneously started propagating from initial notch as well as from the top and bottom corner of the right hand vertical boundary (see bottom right specimen in Fig. 16). These specimens were not considered in Fig. 8. In addition, also experimental data from [76] is included in Fig. 8 in black color for comparison. It is noted, that results from experiments presented in the literature (e.g.  [38,93]) performed on PMMA show velocity toughening characteristics in a similar range.
The fracture patterns obtained from the finest spatial resolution, i.e. for δ = 70μm (in green color in Fig. 8), for nine increasing levels of the initial stored energy G are shown in Fig. 9. A transition is observed from a single smooth propagating crack at G = 170J /m 2 to a rough crack path with short lived branches along the way at G = 348J /m 2 and to stable macro-branching for larger load levels from G = 618J /m 2 onwards. This is the point when a single crack tip is not sufficient to dissipate all the available energy flux and the daughter micro-branches grow into macro-branches. At higher loading levels, the complexity of the fracture process increases, as the generation of daughter micro-branches evolving into new macro-branches repeats itself, eventually leading to a tree-like fracture pattern at δ = 100μm. The initial branching angles obtained at lower loading levels increases from 37 • at U = 30μm to approximately 60 • at U = 70μm. At higher loads, the branching angle almost remains constant and eventually slightly reduces to 57 • at δ = 100μm. Figure 8 shows, that with decreasing size of the peridynamic horizon, the computed crack velocities show a converging behavior (see the results for δ = 70 and 140μm). This convergence is in agreement with [16], who also performed a δ-convergence study with respect to the crack propagation speed using bond-based peridynamics with a critical stretch based fracture criterion, using a plate specimen similar to the one used in this study. They applied transient loads on the specimen boundaries and the crack faces, but without a spectrum of loading levels, which could not have been considered by using only an explicit dynamic solver. They reported, that, as the horizon size is sufficiently reduced, the crack initiation times, fracture patterns, trend and magnitude of the crack propagation speeds become identical ( Fig. 10 in [16]). In our current study, we can compare fracture patterns and crack propagation speeds for different horizon sizes at different loading levels. It can be observed from Fig. 8, that the minimum crack velocity obtained for all four horizon sizes is almost identical (with differences within 1% of v r ). As the loads are increased the fracture patterns stay qualitatively similar as shown in Fig. 4. However, an additional observation from Fig. 8 is, that as the energy stored per unit area G increases, the difference between the average crack velocities for the smallest and the largest horizon size also increases. It is due to the difference in the spatial resolution, as the crack propagating in the finest mesh can generate a larger number of micro-and marco-branches as compared to the coarsest mesh, which is accompanied by a reduction of the average crack propagation speed. It is worth to note, that the macroscopic branches which propagate through the whole length of the specimen start around a similar level of stored energy G = 618J /m 2 for all horizon sizes. This level of G is marked in Fig. 8 by a horizontal gray line. Beyond this loading level, the crack velocities start developing larger differences with respect to the horizon size.
The maximum crack tip velocity observed in experiments [38,76,93] for PMMA is around 70%v r . It must be noted, that the PMMA specimens used in the experiments were approximately ten times larger than the ones used in the current simulations (see Fig. 2a). Performing a convergence study using a full scale specimen would be computationally extremely expensive at this level of spatial resolution and is beyond the context of this study. However, this difference in the specimen size and geometry can cause the cracks to propagate at different speeds, because the stress intensity at the crack tip evolves differently per unit crack length for specimens with different sizes. In addition, the waves reflecting from the boundaries interact more frequently with the crack tip in a smaller specimen as compared to a larger specimen, which also influences the stress intensity at the crack tip. This dependence on the geometry and the size of the specimen for dynamic fracture experiments is also pointed out in [33,51,53]. The terminal velocity for the crack propagation obtained for the smallest horizon size (δ = 70μm) is around 15% ∼ 20% of v r higher than recorded in the experiments. However, the trend of the curves for δ = 70μm and for the experiment in Fig. 8 looks very similar. This quantitative disagreement between the simulations and the experiments [38,76,93] is primarily attributed to the difference in the size and geometry of the specimens. The influence of the specimen size will be discussed in more detail in Sect. 6. Also, small scale material nonlinearities may affect the results as reported in [58]. However, the paper is limited to linear elastic behavior of the material, because the primary aim of this study is to understand the influence of the peridynamic nonlocality and its effects on the dynamic fracture process.

Three dimensional simulations of dynamic crack propagation
This section investigates the influence of three dimensionality on the dynamic fracture process by performing simulations similar to Sect. 4, using, however, now a three dimensional peridynamic model of the PMMA plate shown in Fig. 2a with a thickness of t = 0.5mm. As discussed in Sect. 2, the LPS model suffers from surface effects, since the material points located along the boundaries have incomplete nonlocal neighborhoods. Therefore, for the three dimensional case, the parameter t/δ becomes relevant. For the 3D simulations, it is expected, that the effective stiffness K e f f will depend on the peridynamic horizon because of the considerably larger number of material points along the boundary with incomplete neighborhoods as compared to the 2D case in Sect. 4.1. Hence, one goal of this section is to provide an insight on the influence of these surface effects on the fracture process.

Crack propagation under monotonic load
The load-displacement relationship obtained and the temporal evolution of the crack tip position for the monotonic loading is shown in Fig. 10. As expected, the effective stiffness (slightly) depends on the horizon size (Fig. 10a). For the finest discretization using δ = 70μm, the effective stiffness is almost identical to the K e f f obtained from the two dimensional cases in Sect. 4.1. As regard to the fracture strength, slight increase with decreasing horizon size is observed. When comparing the results for the smallest horizon size, the 3D analysis (F u = 127N ) leads to a slightly smaller ultimate load as compared to the 2D case ((F u = 137N ). Figure 10b shows the temporal evolution of the crack tip position. Due to the differences in the fracture strengths, the cracks start propagating at different times in the four analyses. Unlike in the 2D case presented in Sect. 4.1, the crack propagation velocities for all four peridynamic horizons also show differences in the range of 5% v r . This can be attributed to the differences in the effective elastic properties obtained for different horizon sizes.
For constant density, the stiffness of a material is directly proportional to the volumetric and surface wave speed. For mode I cracks, the latter governs the maximum crack propagation velocity.

Dynamic crack propagation for different levels of quasi-static loadings
In analogy to Sect. 4.2, crack propagation at different levels of stored elastic energy is investigated in this section.

Cracking pattern and crack tip velocity
The cracking patterns obtained from three dimensional simulations are presented in Fig. 11. The fracture pattern is similar to the 2D simulation. As before, a transition from smooth to a rough crack propagation is observed when increasing the load level from U = 21μm to 25μm. When the load level is increased to U = 30μm, macro-branches start developing for the two large horizon sizes only, while a rough single crack with a few small micro-branches is observed for the two small horizon sizes (δ = 140μm and δ = 70μm). For U = 40μm, similar to the 2D simulations, macrobranching is observed for all horizon sizes.
The temporal evolution of the crack tip position and the crack propagation velocity (normalized by v r ) with respect to crack tip position is presented in Fig. 12 for all four loading cases shown in Fig. 11. For U = 21μm, a smooth single crack with slight crack roughening when reaching the right face of the plate is realized, significant differences of the position of the crack tip position are observed in (Fig. 12a). The maximum difference between the average crack propagation velocities obtained from the different horizon sizes is around 9% of v r , while this difference was only 1% for the two dimensional case. For U = 40μm, except for δ = 140μm, all analyses lead to at least two stable branches propagating through the complete length of the specimen, as shown in row 4 of Fig. 11. The initial branching angles, in contrast to the two dimensional case, depend on the horizon size, with the largest branching angle associated with the largest horizon.
In contrast to the 2D simulations, we observe a distinct region of micro-branching along the main crack (see for δ = 70μm at U = 40μm before branching in Fig. 11).
Such micro-branches were not observed in 2D simulations, the cracks rather formed short lived macro-branches which also diverted the crack from propagating in a straight line, as can be seen in Fig. 9. In the region of micro-branching, the crack velocity increases up to 86%v r at 8mm crack length for δ = 70μm, this velocity level is shown by a horizontal dashed green line in Fig. 12h. Around 10mm crack length, two macro-branches start developing, this event is marked by the sudden decrease in crack velocity (around 70%v r ). From this observation we can deduce, that a single crack in the micro-branching region propagates faster than compared to the individual macro-branches if it had branched. As such  micro-branching regime is not available in 2D simulations, this results in the crack propagation velocities obtained from 3D simulations to be slightly higher than the crack velocities obtained in 2D. Figure 13 presents the evolution of the normalized dissipation rate /G c and the normalized elastic, kinetic and dissipated  energies over the crack tip position for all simulations shown in Fig. 11. In contrast to the two dimensional case, where for smooth crack propagation /G c was independent of δ ( Fig.  6a), the three dimensional simulations show a dependence on the peridynamic horizon (see Fig. 13a). This dependence is associated with differences in the effective stiffness and crack propagation velocities with respect to the horizon size. /G c follows a similar trends as K e f f (Fig. 10a), i.e. it increases with increasing δ. This also leads to a slight influence of the peridynamic horizon on the evolution of the portions of energies (Fig. 13b).

Dissipation rate and energy balance
As discussed in Sect. 4.2.2, the kinetic energy around a moving crack tip is directly proportional to the speed of the crack tip. And since the cracks in row 1 of Fig. 11 ( U = 21μm) propagate at different velocities (Fig. 12a), there are differences in the kinetic energy density around the propagating crack tip with respect to δ, as shown in Fig. 14. The region of concentrated kinetic energy around the crack tip for δ = 70μm stays significantly smaller than the one obtained for δ = 280μm during the crack propagation. This leads to the differences observed in the evolution of the total kinetic energy seen in Fig. 13b.
As discussed earlier, the dissipation rate is higher for a larger horizon, this is why we see macro-branching taking place in larger horizon sizes before the smaller ones in Fig.  11. In contrast to the 2D simulations, here we observe that in 3D it is possible to have a single crack propagating with global /G c > 2.0, i.e. without macro-branching, as can be seen for δ = 70 and 140μm in Fig. 13e. At U = 40μm, macro-branching is observed for all four horizons (row 4 of Fig. 11). In the three dimensional case a single crack can sustain significantly higher /G c as compared to the two dimensional case where macro-branching occurs as soon as /G c > 2.0. This high dissipation rate is maintained by the development of these quasi-periodic micro-branches before the development of stable daughter cracks. These microbranches occurring before stable macro-branches can be seen for δ = 70μm in row 4 of Fig. 11.

Crack velocity toughening mechanism
In analogy to the 2D case, the average crack velocities v c normalized with respect to v r are plotted in Fig. 15 for different levels of G stored in the specimen prior to fracture and for all considered horizon sizes. At the load level U = 40μm (i.e. G = 618J /m 2 ), initial cracks developed macroscopic branches soon after initiation, which fully propagated through the whole length of the specimen (except δ = 140μm). This level of G is shown in Fig. 15 using a horizontal gray line. Compared to the 2D case, the crack speed, even for the smallest horizon size, is considerably large and almost approaches 95%v r . As discussed in Sect. 5.2.1, the micro-branching mechanism observed in 3D simulations leads to a faster propagating crack tip in comparison to the speed of the individual branches if it had formed macro-branches (which is the case for the 2D simulations). A detailed discussion comparing 2D and 3D simulations is provided in Sect. 6.
The fracture patterns resulting from nine increasing levels of the initial stored energy G are presented in Fig. 16 for the finest spatial resolution, i.e. for δ = 70μm. The figure shows the transition from a single smooth propagating crack at U = 241J /m 2 , to a rough crack path at G = 348J /m 2  (Experimental data taken from [76]) and to the formation of micro-branches at G = 618J /m 2 , which finally lead to macro-branching. A maximum branching angle of 40 • is obtained for U = 60μm which is smaller than the maximum branching angle for the 2D case recorded as 60 • at U = 70μm. Figure 17 shows a three dimensional view of the fracture patterns obtained from six levels of U , ranging from 25 to 70μm after 10μs. The transition of the surface roughness from a smooth to a rough state for the fracture obtained for G = 348J /m 2 is presented in Fig. 18. It is interesting to note, that the roughness of the simulated fracture surfaces is able to qualitatively reproduce the well known mirror-misthackle transition. Fracture surfaces obtained by experiments in [74] are also presented in the second row of Fig. 18 for a qualitative comparison. A surface, denoted as mirror is characterized by the propagation of a single sharp crack, creating a smooth and featureless crack surface. The misty zone is characterized by the daughter cracks which initiated (locally along the specimen thickness) but got arrested before forming a micro-branch. Finally, the hackle zone is caused by daughter cracks, which successfully form quasi-periodic micro-branches, some of which grow into macro-branches. The periodic features on the crack surface obtained in the hackle zone from the simulations are not in a quantitative agreement with the experiment, nonetheless the periodic topology is qualitatively similar.

Influence of the peridynamic horizon
Considering the load displacement relationship obtained from two and three dimensional crack propagation simula-tions of PMMA plates under monotonic loads presented in Figs. 3a and 10a, it is observed, that the effective stiffness almost shows no dependence on the peridynamic horizon in the 2D case, while a notable dependence is observed for the 3D case. Related to this observation, in contrast to the 2D simulations, the 3D simulations also show an influence of the peridynamic horizon on the crack propagation velocities (compare Figs. 5 and 12 for a loading level of U = 21μm). This is not surprising, as the crack propagation velocity depends on the elastic wave propagation velocities, which are a function of the effective stiffness. The dependence of the effective stiffness on the horizon size for the three dimensional case is a consequence of the surface effects of peridynamic models. In the two dimensional case, the points with incomplete neighborhoods are only present at the vertical and horizontal boundaries. In the three dimensional case, in addition to the vertical and horizontal boundaries, all points along the x-y surfaces (at z = 0 and z = 0.5mm Fig. 2a) also have incomplete neighborhoods. From the relationships between G and v c obtained from two and three dimensional simulations, one can observe that the difference between the crack propagation velocities with respect to δ increases with increasing G (compare Figs. 8 and 15). An increase in G corresponds to crack propagation involving crack roughening, micro-as well as macro-branching. As discussed earlier, the crack velocity reduces as soon as the crack branches. Crack roughening and micro-branching are also localized bifurcation attempts along the crack front which also reduce the main crack tip velocity. A crack tip propagating in a simulation domain with a very fine spatial resolution is able to form these features along the crack path more frequent as compared to a simulation with a coarse spatial resolution. Evidently, in both two and three dimensions, simulations with a high spatial resolution (e.g. x = 20μm) is able to branch more frequent (slightly reducing crack speed at each bifurcation) as compared to simulation with a coarse discretization (e.g.  . 19 Comparison of the stored energy-crack tip velocity relationship computed for 2D and 3D simulations for δ = 70μm. The horizontal gray line marks the threshold for macroscopic branching x = 80μm). Hence, as G is increased, micro-and macrobranching starts, and the frequency at which the crack is able to bifurcate will depend on the spatial resolution. This leads to different crack propagation velocities observed for different peridynamic horizons at increasing loads.

Influence of dimensionality
As the effective stiffness for δ = 70μm in 2D as well as 3D simulations is nearly identical, we can assume that the influence of the surface effects is not significant in this case. However, the crack propagation velocities obtained from two and three dimensional simulations for δ = 70μm differ significantly, as is illustrated in Fig. 19, in which the normalized crack propagation velocities (v c /v r ) obtained from 2D and 3D simulations at different levels of energy per unit area (G) stored in front of the initial notch are plotted. If only the effective stiffness would be governing the crack propagation velocity, one should not observe significant differences in the crack propagation velocities. Comparing again the fracture topologies for both cases, it must be concluded, that one of the major aspects that contributes to the difference of the crack propagation velocities obtained from 2D and 3D simulations is the micro-branching instability. The fracture patterns obtained from the 2D simulations presented in Fig. 9 are qualitatively different than the ones obtained from the three dimensional simulations shown in Fig. 16. The key qualitative differences can be best seen for U = 30μm and U = 40μm (before branching). For the ease of comparison, these two cases are presented side by side in Fig. 20. In the 3D case the macro-crack propagates in a straight line with short lived branches on either sides, while in the 2D case a pronounced crack tip splitting is observed which results in the crack tip being diverted from propagating in a straight horizontal line.
As discussed in Sect. 5.2.1, macro-branching reduces the crack propagation speed of individual branches significantly as compared to a single crack propagating in the microbranching regime. Micro-branching can be interpreted as an unsuccessful attempt by the crack tip to form a stable branch [14,75]. Bifurcation of the crack front occurs locally along the thickness of the specimen, and, given sufficient energy, this local bifurcation will expand along the thickness of the specimen to form a macro-branch. During expansion of the secondary crack along the thickness, the main crack keeps propagating. This reduces the driving force at this newly formed micro-branch and it gets arrested. When these localized bifurcations of the crack front are unable to expand through the specimen thickness, the fracture surface results in the "misty" zone shown in Fig. 18. This process can happen repeatedly in a three dimensional simulation and provides a mechanism by which these simulations can even maintain /G c > 2.0 without forming stable branches (see Fig. 13e, g). In the periodic micro-branching regime, the crack surface generates a pattern known as the "hackle" zone, as shown in Fig. 18. In contrast, two dimensional simulations cannot provide such mechanisms for a propagating crack tip, as any bifurcation at the crack tip will lead to branching.

Influence of the specimen size
The differences between the crack propagation velocities obtained from the peridynamic simulations and experiments (Fig. 19) is attributed mainly to the different sizes of the PMMA specimens used in the experiment and the simulations. In the experiments performed in [38,76], the largest dimension of PMMA specimen used was 400mm. In the experiment reported in [93], the maximum dimension of the specimen is 320mm. In contrast, the simulations presented here use a specimen which has the largest dimension of 32mm.
Such large difference in the specimen size can lead to the differences observed in crack propagation velocity. The reason can be seen in the scaling for the simple case where the energy release rate ∝ σ 2 l (σ is the applied stress and l is the instantaneous crack length). Considering two specimens where one is scaled in all dimensions (including the initial notch) to be significantly smaller than the other. A crack initiates in both specimens when = G c , i.e. when the energy release rate corresponds to the Griffith condition [43]. In the smaller specimen, since the initial crack length l is small, the condition = G c is realized at the cost of very high applied stresses σ . As a crack starts to propagate (assuming constant σ ), the energy release rate will double when l (the initial notch length) doubles, i.e. over a very short propagation distance for the small specimen and vice versa for the large specimen. In other words, the evolution of Fig. 20 Comparison of the fracture patterns obtained from 2D (left column) and 3D simulations (right column) for δ = 70μm over the crack length will differ significantly for different specimen sizes (with higher for smaller specimen) and the corresponding crack tip acceleration and velocity will also differ accordingly [75]. Additionally, depending on the size of the specimen, the elastic waves reflecting from the boundaries interact at a different frequency with the crack tip. This also influences the stress intensity at the crack tip.
The choice of specimen size used in the paper so far is constrained by the computational expense. However, in order to assess the ability of the peridynamics model to correctly reproduce the fracture mechanics size effect [5] for quasistatic loading and the influence of the specimen size on the computed toughness characteristics, a number of additional simulations were performed using additional specimen sizes. The dimensions of the PMMA specimen presented in Fig.  2a are scaled according to a factor λ (except the thickness, which is kept constant at 0.5mm). Assuming λ = 1.0 for the specimen size 16mm × 32mm, we perform peridynamic simulations for different specimen sizes up to λ = 10 (which corresponds to the size of the specimens used in the experiments [38,76,93]). These simulations are performed using a two dimensional formulation since a 3D simulation would be computationally extremely expensive. We use a discretization size of x = 40μm with δ = 140μm, which leads to 8.0 and 32.0 million material points for λ = 5 and λ = 10, respectively. A discretization size of 20μm for λ = 10.0 would lead to 128 million material points in 2D and 3.2 billion material points for 3D analyses. Solving such a large system of equations is out of the context of this study.
Two sets of simulations are performed for various values of the size parameter λ. The first set is concerned with monotonic loading, similar to Sect. 4.1, considering λ = 0.1, 1.0, 5.0, and 10.0. Loading is applied at a rate of 3mm/s. These simulations are performed to investigate the relationship between the failure stress σ and the specimen size λ. According to Linear Elastic Fracture Mechanics, the fracture strength σ for specimens of different sizes λ should follow the scaling σ = c/ √ λ. On a log − log scale, this gives a straight line relationship, with a slope of −1/2, between σ and λ. c is the strength for λ = 1.0, obtained for δ = 140μm as c = 8.535 M Pa (see Fig. 3a). The expected strength scaling of the LEFM is compared with the strengths obtained from the peridynamic simulations (on a log − log scale) in Fig. 21. This figure confirms, that the strength scaling obtained from peridynamic simulations for different λs is able to reproduce the LEFM size effect [5].
In the second set of simulations, we investigate the velocity toughening behavior, i.e. the relationship between G and v c , for λ = 1.0, 5.0, and 10.0. These simulations are performed similar as reported in Sect. 4.2. Fig. 22 compares the computational results for the different specimen sizes with the experimental data from [76] (in black color). It shows a reduction of the crack propagation speed with increasing λ. For λ = 10.0, which agrees with the actual size of the specimen tested in the laboratory, we start to see an agreement between the experimental values [76] and the simulations for higher levels of G. However, for low values of G, the computed crack speed is still higher than observed in the experiment. A possible cause for this deviation is the larger influence of small-scale material nonlinearities [58] at lower values of G, which become insignificant compared to the total fracture induced dissipation at higher values of G. Figure 22 demonstrates the dependence of the relationship between G and v c on the specimen size for a given material. This dependence has been discussed in [33,51,53], however further computational as well as experimental investigation in to the characterization of this size effect is needed.

Conclusions
In this study, the influence of the peridynamic horizon on the velocity toughening relationship and the fracture behavior of PMMA has been investigated for a spectrum of loadings. The simulations were performed using two as well as three dimensional peridynamic formulations and the specimens were loaded using displacement controlled boundary conditions. Two loading scenarios were investigated in both cases: i) crack propagation under monotonically increasing loads and ii) dynamic crack propagation in quasi-statically preloaded plates. From the 2D monotonic loading simulations, which results in a single straight crack propagation path, we found that the effective stiffness as well as the crack propagation velocity is independent of δ and the fracture strengths converge as δ is decreased. In contrast, the 3D simulations show a slight dependence of the effective stiffness and the crack propagation speed on δ.
In the second loading case, the crack tip velocity was computed for different levels of the initially stored energy G in the specimen. Qualitative differences in the crack propagation between the 2D and the 3D simulations have been observed. In 3D simulations the crack front can bifurcate locally (along the specimen thickness) while keeping a relatively straight crack propagation path. A local crack front bifurcation can grow into a micro-branch, which arrests soon or grows into a macro-branch, depending on the available energy and the interaction of this micro-branch with the main crack front. These arrested localized bifurcations and micro-branches leave distinct localized features on the crack surface, known as "mirror", "mist" and "hackle" fracture surfaces. The high resolution of the presented 3D peridynamics simulations is able to reproduce these surface features reported in the literature. The fact, that such mechanisms are not available in a 2D analysis explains the qualitative differences observed between the crack propagation processes in 2D and 3D simulations.
In this study, we found two sources responsible for influencing the crack propagation speed with respect to the size of the peridynamics horizon. The first mechanism is the crack surface roughening and micro-branching, which reduce the main crack propagation speed at each bifurcation attempt. Since a finer spatial resolution, characterized by a smaller horizon, provides more possibilities for the crack tip to bifurcate as compared to a coarse mesh, the average crack propagation speed is lower in a simulation with a smaller discretization size. The differences in crack velocities with respect to the horizon size increase with increasing loading level. This influences 2D and 3D simulations in a similar manner. The second cause is a consequence of the nonlocal nature of the peridynamic formulation and the influence of missing non-local neighborhoods along the surfaces. Since this influence is larger for 3D simulations as compared to 2D analyses, where this effect is negligible, the effective stiffness and, consequently, the crack velocity and the crack branching characteristics obtained from 3D and 2D simulations differ.
The influence of the specimen size on fracture propagation was also analyzed using two dimensional simulations. It was shown, that the peridynamics model is able to exactly reproduce the theoretical fracture mechanics size effect law. As the specimen size in the simulations was increased to match the size of the specimen actually used in the experiment, an agreement of the average crack propagation velocities obtained for higher levels of the initially stored energy G was observed. This result showed, that the crack velocity toughening behavior is actually dependent on the size of the specimen and should not be considered as a material property as usually assumed in numerical models. Further experimen-tal as well as computational studies are needed to characterize this size effect.
However, for lower levels of G the crack propagation velocities obtained from simulations are still higher than the experimental values. We suspect, that this difference at lower levels of stored energies can be attributed to the fact, that PMMA is not an ideal brittle material and that local material nonlinearities at the crack tip contribute substantially to the total energy dissipation for lower levels of G. This relative contribution becomes insignificant as compared to the energy dissipated in the fracture process at higher levels of G.
It is noted, that the present investigation has been restricted to structured grids. In future research, the δ−convergence study presented here needs to be supplemented with a n−convergence study. In addition, besides considering also unstructured grids, the benefit of using discretization techniques with higher order convergence rates [45,56,72,82] could be further topics of research.