Peridynamic Simulation of a Mixed-Mode Fracture Experiment in PMMA Utilizing an Adaptive-Time Stepping for an Explicit Solver

In this paper, a benchmark analysis of a peridynamic correspondence energy-based damage model is presented. The benchmark is an experimental setup of a Polymethyl methacrylate (PMMA) plate with a hole. The plate has a minotch and is subject to a compressive load. With increasing loads, a crack initiates at the tip of the notch and continuously grows. The benchmark is modeled utilizing the peridynamic correspondence formulation as a two-dimensional problem. To reduce numerical issues due to bond failure, an adaptive time-stepping method for a Verlet time integration schema is proposed. The method limits the maximum number of broken bonds per material point by adapting the time-step size. This allows the correspondence formulation to be significantly more stable. The benchmark involves a sensitivity analysis based on the Morris method, which is performed in this context. As a result, uncertainties and the impact of geometrical, numerical and material parameters are evaluated and discussed.


Introduction
Wind energy is one of the main pillars for the production of low-carbon emission electricity. To keep electricity affordable, the turbines are becoming increasingly larger and as a result the structural loads are increasing as well. Much more precise methods are needed to design these larger turbines [1][2][3]. In particular, the adhesive joints in the wind rotor blades are a critical component. In rotor blades, usually these adhesive joints are thick and alternative methods to cohesive zone modelling are needed, e.g., the peridynamic theory. Peridynamics was introduced in the early 2000s by Stewart Silling and has been continuously developed since then. The main advantage of the method is that there are no singularities in the region of damage as in classical continuum mechanics [4]. This is achieved by an integral description of the conservation of momentum. The disadvantage is that an integration domain must be defined. The size of the domain in a macro-scale is no material parameter and therefore has to be chosen by the user. For homogeneous stress states without discontinuities and a integration domain size which approaches zero, the solution converges to that of classical continuum mechanics [5][6][7]. For problems with discontinuity, convergence can only be ensured for one domain size. That means the error runs towards zero with a fixed horizon. If the domain size approaches zero, the same problems result as with classical continuum mechanics. In summary, it can be said that the chall of Peridynamics is to clearly define integration domain size rules for damage problems from a mathematical or material modeling point and not from the best practice point of view.
In order to predict failure in a thick adhesive layer, crack propagation in isotropic material must be simulatable. Such a material model has already been implemented and verified for the ordinary state-based formulation [8,9]. The methodology has already been applied at the microscale [10] to study more complex phenomena. Rädel et al. [11] showed that the convergence of the pure ordinary state-based formulation is significantly worse compared to the so-called correspondence formulation. This is due to the bad geometrical approximation. The position-aware linear solid constitutive model solves this issue [12]. However, the complexity of the formulation increases. In perspective, the models to be computed are increasing in size and complexity and the correspondence formulation provides clearly more flexibility in material modelling compared to ordinary-state based formulations. Even if isotropic materials can be modeled by ordinary-state based formulations, the correspondence formulation is used in this paper. The aim of our research is to gradually increase the complexity of various questions and thus to gain confidence in the methodology, e.g., the evaluation of manufacturing deviations of fibre composites [13,14]. The field of Peridynamics being a promising approach is shown by the large number of research questions that have been addressed in recent years [15][16][17].
This paper aims to lay the foundation for complex analyses on bonded fiber composite structures. For this purpose, the damage model developed and verified in [9] is adapted for the correspondence formulation. One assumption is that the peridynamic model needs less parameters to model the same crack propagation behavior compared to fracture mechanical models included in finite element analysis tools. This has been shown by Willberg et al. [9] for an isotropic plate with a hole under tensile loading. This paper elaborates on the behavior and the quality of the peridynamic modelling for a more complex benchmark. The experimental setup of the benchmark is based on the dataset provided by Julien Réthoré [18]. He tested a PMMA (Polymethyl methacrylate or plexiglas) plate with a hole and a notch starting from that hole in an angle of approximately 30 • with respect to the horizontal axis. The dataset has some limitations in terms of documentation. Therefore, one goal of the paper is to show a methodical approach to still generate added value from such a dataset. The reason for using this data set is that there are still to few experimental data sets that represent complex load conditions. Within this paper, a sensitivity analysis is carried out in order to determine the robustness of the proposed modelling approach. For this purpose, the Morris method is used to vary several most relevant parameters. This is to show that the model actually allows a physical description of the problem and does not represent a physically motivated curve fit. Even if the experiment is designed as a quasi-static problem, a Verlet solver is utilized. This is due to the fact that the quasi-static solver was not yet adapted to the energy-based damage model at the time the model was created. The goal of the current presented benchmark is not only to recalculate the results. Gee et al. [19] used the benchmark results to analyze two extended finite element method (XFEM) models. Both models were in good agreement, but needed three parameters to describe the damage initiation and progression. With the peridynamic model we try to use only one parameter to get the same results. The paper is structured as follows. The problem is described, and the requirements for modeling are derived. This includes the consideration of uncertainties. Afterward, the methodology of robustness analysis is presented and justified. Subsequently, the theory of material and damage modeling in peridynamics is presented. One focus lies on the compensation of the so-called zero energy modes and an adaptive time-step size technique is presented and its effect on the solutions is analyzed, in order to compute the models more efficiently. After the problem has been analyzed and the theoretical basics of modeling have been set, the benchmark is calculated and the results are explained and discussed.
Finally, the paper will be concluded and further research on this topic will be discussed.

Experimental and Numerical Setup
The following description is based on the dataset provided by Réthoré [18]. The geometry information is shown in Fig. 1. It shows an image of the plate with a speckle pattern for displacement field measurements. As visible in the figure the precise notch path was determined.
The plate is loaded under pressure with a displacement of 1 × 10 −4 m min −1 . The experiment is designed to initiate crack growth at the tip of the notch. This crack grows steadily as the load increases under mixed mode conditions. The pressure is introduced by an external displacement in the x-direction at the right side at x = 140 × 10 −3 m . The left side at x = 0m is fixed. The notch causes deformation in the x-and y-directions. Stress concentrations occur at the notch tip. This concentration leads to a crack initiation at the notch tip and with further loading the crack grows upwards and to the left. Considering the total  Table 1 Fig. 2 Crack path in digital image correlation technique (DIC) measurement results at the end of the experiment [18] duration of the experiment of 480 s and the calculated maximum crack propagation speed of 9.4 × 10 −4 ms −1 , dynamic behavior is neglected. Moreover, the crack-growth can be considered as quasi-static. Effects such as a heating crack tip described by Mehrmashhadi et al. [20] are therefore also neglected. Figure 2 shows the crack path provided by the benchmark data. The final crack-tip xand y-coordinates are given in the figure. The red lines define the positions of the load introduction and the bearing, respectively. The coordinates of the milled pre-notch and the final crack path are given in Table 1.
Additionally, to the crack pattern the benchmark data provided the load-displacement curve (cf. Fig. 3) and the crack-length plotted against the load steps. For better comparability, the crack-length evolution curve is shown in relation to the external displacement in Fig. 3b. With regard to the material parameters, the benchmark offers only limited information. The Young's modulus given in the dataset ( 5000 × 10 6 Nm −2 ) is overestimated to the authors' knowledge. Based on the load-displacement curve provided, the value should be approximately 3250 × 10 6 Nm −2 . The Poisson's ratio is not provided and is assumed by the authors. Gee et al. [19] provide an energy release rate ranging from 0.23 kJ m −1 to 3.7 kJ m −1 and fitted this value to 0.24 kJ m −1 . In Fig. 3 of the benchmark description [18], the mode I fracture toughness K IC is given with ≈ 1.95 MPa∕ √ m . For plane stress conditions, the energy release rate G IC can be determined as:  with E as the Young's modulus [21], which is much smaller as the value given by Gee et al. An overview over the material data of the PMMA plate is provided in Table 2. The aforementioned uncertainties in the input parameters emphasize the need to perform a sensitivity analysis to estimate the effect of parameter variation.

Sensitivity Analysis
The numerical experiment of the plate made of isotropic material with hole, pre-notch under tensile loading to generate a stable crack propagation, is a great challenge. Unlike the physical experiment, the numerical experiment requires a number of input variables that are intrinsically present in the physical experiment. In other words, the effects, interactions or influences of multiple parameters are summed up in the final experimental result (cf. Eq. (1)). A virtual experiment represented in a numerical model can individually map a large amount of these parameters. This is what enables a systematic investigation into the various sources of uncertainties. However, in Sect. 2 it becomes clear that correct numerical modelling requires knowledge about geometry, material properties, constraints and boundary conditions, etc. Hence, we need to know the uncertainties in the output values of the numerical model that result from uncertainties in the input parameters. In our work, we use sensitivity analysis to investigate different sources of uncertainty in the output values. This approach is motivated by Saltelli [22], who describes the purpose of a sensitivity analysis as follows: The study of how the uncertainty in the output of a model (numerical or otherwise) can be apportioned to different sources of uncertainty in the model input. [23, p. 45] The numerical model takes several hours to run through the analysis. For the computationally time-consuming models, quantitative methods are not appropriate [24]. For that reason, in this work we rely on the Morris method [25] with the extensions proposed by [26] to derive sensitivity indices. The Morris method is computationally cheaper than more sophisticated variance-based sensitivity analysis methods (e.g., FAST and Sobol) [27]. Morris' method can be considered as a screening technique. It varies one-factor-at-a-time and may therefore be referred to as an OAT approach. Given the context of the benchmark experiment, we want to use the Morris method to figure out which "possible settings" yield the "desired properties" in the model. This enables the deduction of the effects of individual parameters on the model behavior.
For our sensitivity analysis, we make use of nine parameters. All of them are given with their mean value, range and a short description in Table 2. Most of them are determined by an analysis of the experimental setting and the given experimental results in Sect. 2. Besides the fact that the mass density is a physical value in our setup, it is also a numerical parameter to increase the time integration step within the used explicit timeintegration solver.
The whole experiment is designed as a quasi-static problem. Therefore, time scaling is assumed as a valid method as long as a steady crack growth occurs. The influence of larger mass densities, and its effect, is analyzed to check this assumption. Two other experimentally not obvious parameters are considered also: the peridynamical horizon and a variation of the load introduction m l . It is assumed that the external displacement is not applied fully horizontally. We realize this by the following function: The Poisson's ratio has been included, because often its effect in fracture mechanics is not taken into account [28]. However, due to the occurrence of shear stresses at the notch tip the initiation as well as the direction of crack propagation is affected by this parameter. With the help of this function, the effect of not ideal boundary conditions is considered.
The Morris method trajectories (sets of parameter combinations) are based on a Latin Hypercube Sampling (LHS) of size 20. This OAT design is defined in the box [0, 1] 9 , with p = 9 model input parameters. We define the number of randomly chosen trajectories to r = 5 . Each input variable is uniformly distributed with bounds listed in Table 2 and illustrated for the geometrical parameter in Fig. 4. With r repetitions per parameter p, the Morris screening requires 50 simulations.
The Morris method uses so-called elementary effects ee i to measure the overall effect * i and the higher-order effect i of each parameter i on the model response. From the sample of size r, the mean and standard deviation for each elementary effect are obtained:

Peridynamic Formulation and Stabilization
The paper follows the assumptions and notations from Silling [4]. Within the neighborhood H , with the volume V , defined by a spherical domain the horizon , the force volume density state for the bond interaction between the positions and ′ is defined as the integral balance of momentum Three variations of the peridynamic model are currently being used, the bond-based as special case of the ordinary state-based, the ordinary state-based and the non-ordinary state-based formulation. In this order, flexibility increases, but so does the complexity of the formulations.

Material Model
A special formulation of the non-ordinary state-based model was introduced by Silling et al. [4] in 2007. This so-called correspondence formulation defines an integral non-local deformation gradient to calculate the bond force vector density states, where is the shape tensor with V defined as the neighborhood volume. For each bond , there is an influence function ⟨ ⟩ , an undeformed vector state ⟨ ⟩ and a deformed vector state ⟨ ⟩ . The shape tensor has to be positive definite and symmetric. The advantage of using the non-local deformation gradient is that classical continuum mechanics constitutive models can be used in Peridynamics. The peridynamic force density vector state is thus The Piola-Kirchhoff stress tensor with respect to an orthonormal basis can be determined as In this benchmark, we deal with a linear elastic material. In this case, the Cauchy stress can be derived using Hook's law as with the fourth-order elasticity tensor and the Green-Lagrange strain tensor defined as For correspondence models, the so called zero-energy modes could occur [29]. These modes are non-physical and lead to unstable or unreasonable solutions. Several stabilization methods were published to overcome this problem [30][31][32][33][34][35]. A promising approach was published by Wan et al. in 2019 [36]. Instead of a bond-based stabilization method proposed by Silling [37], Wan et al. developed a state-based stabilization method. As positive side effect, this method stabilizes the solution for anisotropic material as well. The corrected force density state C with suppression of the zero-energy mode is where is given in Eq. (8). Following Wan et al. [36], the suppression force density state S is with as the non-uniform deformation state caused by the zero-energy mode. If the approximated non-local deformation gradient exactly maps each undeformed bond to the deformed configuration no zero-energy mode occurs. In that case, the non-uniform deformation state is zero and the corrected force density state C is equal to the force density state . To obtain the second-order tensor 1 , the fourth-order elasticity tensor from Eq. (10) maps the second-order shape tensor inverse utilizing the double dot product as

Damage Model
An energy-based damage model has been utilized [8]. This model is superior to the widely used and easy to implement critical stretch model. The reason is that the critical stretch cannot be determined experimentally. The change in distance between two points does not say anything about how it is achieved. This becomes clear when decomposing into an isotropic and a deviatoric part as it is done typically for isotropic materials. The fracture energy of a material is equal. The critical stretch for a shear dominated load case is different from a tension dominated load case. Hence, the same fracture energy can be generated by different critical stretches. As a result, the critical stretch fitted to the energy release rate does not reproduce the characteristic value in a virtual experiment. Following Foster et al. [8] for the correspondence material, the bond energy density is where the projected force density state P is given by projecting the corrected force vector states C onto the relative displacement with state A bond fails if it exceeds a specific critical energy w c . Foster et al. [8] derived a maximum elastic bond potential value w c . This value is based on the energy release rate G 0 and the horizon . The elastic bond potential is given in Eq. (18) (a) for three-dimensional and Eq. (18) (b) for the two-dimensional case with thickness h.
This criterion does not differ between traction and compression. In this paper, the energy-based formulation is adopted. Only bonds with positive stretch � ⟨ ⟩� − � ⟨ ⟩� > 0 are analyzed whether the critical energy w c is reached or not. The whole theory is implemented in the Peridigm framework [38]. The analysis in this paper is conducted using a Verlet time-integration schema. The whole solving process is sketched in Willberg et al. [9] in "Section 4 Implementation", cf. Algorithm 1. The time-integration is a placeholder and can be a static iterative solver as well as dynamic implicit or explicit Verlet solver. Within this study, the latter one was used, because most of the development was done there.

Adaptive Time-Stepping
As already described in Sect. 4.2, damages in the numerical model occur if bonds are deleted once a damage criterion is fulfilled. This might cause breakage of multiple bonds connected to one material point in one time-step. This causes parts of the material's history to be lost, since no new load paths are established. Deleted bonds affect the local stress distribution and can influence the crack initiation and / or propagation. Besides, this modelling shortcoming potentially contributes to other numerical issues. Using a representative model, we found that the calculation of the non-local deformation gradient and the compensation of zero-energy modes no longer work reliably. We determined the reliability threshold for the amount of broken bonds per material point and time step to be seven. As shown in Eq. (14), the non-uniform deformation state ⟨ ⟩ is determined by subtracting the deformed vector state ⟨ ⟩ from the deformed state obtained by the non-local deformation gradient . As shown in Eq. (6), the non-local deformation gradient is determined via the influence function ⟨ ⟩ , including the currently detached bonds. This means a determined state without deleted bonds will be compared with a state involving deleted bonds. This will ultimately lead to numerical issues if the nonuniform deformation state, and thus, the suppression force density state S , see Eq. (13), becomes too large. This section introduces a novel adaptive time-stepping method for a Verlet time-integration schema, which will be utilized to dampen the numerical issues mentioned above. Moreover, the schematic implementation in the Peridigm framework [38] and an exemplary verification will be outlined.
In order to generate a precise crack prediction, every broken bond should initiate a redistribution of the body forces. Otherwise, this behavior could possibly lead to instabilities and inaccurate results. However, the amount of bonds involved in such a model would significantly increase the computational time. A trade-off between efficiency and accuracy can be achieved by implementing a method that controls the maximum number of broken bonds. Ni et al. [39] demonstrated the need for such considerations and proposed a similar approach but for an implicit quasi-static algorithm. They obtained the total number of broken bonds per iteration and compared them to a predefined maximum value N m b . If more bonds are broken than a user-defined value, the peridynamic stiffness matrix will be updated. A similar approach has been done by Zhao et al. [40].
As the solver used in this paper utilizes an explicit analysis, a different concept is required. Typically, the smallest time-step is given by the largest Eigenfrequency of the model [41]. Unfortunately, this time-step may be numerically stable but too large to prevent multiple bond failure.
In contrast to Ni et al. [39], only the number of broken bonds of each node will be considered. Ni et al. allowed much larger time-steps, so the overall model was considered. This does not apply to the explicit Verlet time-integration. Hereby, the speed of sound has to be resolved via the time-integration method. Hence, in each time-step information could only travel from one point to the neighboring point. Thus, the maximum influential sphere of information (i.e., load path redistribution) complies to the smallest distance between two material points in the overall model. Eventually, this significantly increases the level of efficiency and eases the implementation of this adaptive algorithm.
Another unique feature of our proposed concept is that the number of broken bonds will not be controlled directly through an iterative update, but rather by means of a reduced time-step. By repeating the last step with a smaller time value, the applied body forces are reduced and so are the number of broken bonds. The implementation, which is visualized in Fig. 13, requires a number of parameters to be defined: -Maximum number of broken bonds per step and node N m b -Stable number of broken bonds per step and node N s b -Number of steps before potential time raise step d -Time step reduce factor t reduce -Time step raise factor t raise -Time step maximum t max -Time step minimum t min The adaptive time-stepping starts by determining if an iteration is accompanied by a high, normal or low number of broken bonds. To detect this, in every iteration for each node the number of broken bonds will be cumulated. If overall N b ≤ N s b is valid, the timestep is sufficient. Based on the previous step, a decision is made between three different scenarios: 1. The overall damage calculation is sufficient, and the time-step remains unchanged: dt n+1 = dt n 2. The previous iteration revealed a high number of broken bonds, and therefore, the timestep will be reduced: dt n+1 = dt n t reduce 3. The time-step is lower than the maximal time-step and the previous iteration was stable, thus the time-step can be raised: dt n+1 = dt n t raise In order to verify the implementation of the adaptive time-stepping method, a simplified but representative model was used. As visualized in Fig. 5 (left), this model is capable of reproducing aforementioned instabilities more efficiently. The numerical instability increases the local loads and bond energies. Therefore, many bonds break as shown in the figure. The adaptive time-stepping compensates this effect, and the crack propagates without numerical instability, cf. Fig. 5 (right).

Results and Discussions
Initial configurations were calculated to determine a suitable discretization and horizon size. Several convergence criteria were chosen, namely the gradient of the forcedisplacement curve, the force and displacement of crack initiation and the crack path pattern. Convergence was reached if no significant changes between to discretization and their corresponded horizon exist. Figure 7 shows a valid configuration for which the crack length of both the experiment and the peridynamic simulation matches. The pre-notch was implemented utilizing a so-called bond filter. A rectangular plane is introduced at the center of the crack width b c ∕2 . Each bond which would penetrate the plane is not created. The discretization and horizon size are decisive for the accuracy of the geometric approximation of the hole. At coarser discretizations and as a result of larger horizons, the likelihood increases that an irregular crack appears on the right-hand side of the hole. The irregular crack is caused by local stress peaks (a kind of numerical notch effect), which promotes an early failure [42]. This effect can be reduced by increasing the spatial resolution and / or using a conforming grid, meaning that all material points are positioned at the hole's circumference. According to advanced analyses, an ordinary-state-based formulation requires a finer discretization than the correspondence formulation to represent the crack pattern. Figure 6a, b shows the crack path for the ordinary and the correspondence formulation with the same horizon size and discretizations. It can be seen that ordinary state-based formulation does not represent the crack path accurately. A crack occurs at the right side of the hole. This happens also for the correspondence formulation, but for coarser spatial discretization and therefore larger horizon sizes. The force-displacement curve also differs between the ordinary state-based and correspondence formulation. Both are in agreement with Mitchel et al. [12] which show that the pure ordinary state-based formulation has problems with surface effects.
This behavior agrees with the convergence studies of Rädel et al. [11] which showed that the geometric approximation of the correspondence formulation is significantly better than the pure ordinary state-based formulation. It must be noted in the paper that a sensitivity analysis between models was performed. The error due to modeling in less important when the models are compared to each other and the sensitivity of each parameter is considered, since the modeling error exists in all models.

Notes About the Result
Many of the calculations became unstable before the final time step. The reason for this is not always obvious. In order to perform the sensitivity analysis without loss of results, the analyses were used up to an external displacement of 0.59 × 10 −3 m . This value is sufficiently far from the crack initiation of all calculations and the actual experiment. Therefore, a sensitivity analysis after crack initiation is possible. The load-displacement curves and the crack length displacement curves instead show the results up to the end of the respective calculations. Figure 2 shows the crack path in the experiment. The crack starts at the crack tip of the milled notch and grows in a curve to the left. Figure 8 shows all the results of the parameter variation. To show the probability of spatial occurrence of a crack path, the colored pixels of a damaged index area were summed up. A high number means that a large number of cracks run through this position, meaning that the crack path is robust against parameter variations of Table 2. All simulated cracks (red color) start at the crack tip of the initial milled notch and run as a curve to the left, as in the experiment. After the curved shaped path, the crack continues to grow in an almost horizontal direction. The curves differ from the experiment depending on the parameter set in the y-position. The width of this distribution is in the range of 8 × 10 −3 m in y-direction. The peridynamic solution is able to model robustly the experimental crack path initiation. Changes in parameters do not strongly affect the overall pattern of the initiation path, but do influence the crack length. Roughly 50% of all analysis were able to predict the full crack.

Load-Displacement Curve
The experimental load-displacement curve is given in Fig. 3. Basically, several characteristics of this curve challenge the peridynamic simulation: -There is a kink approx. at 0.36 × 10 −3 m displacement, illustrated as a dashed line in Fig. 3a. This phenomenon is not discussed in the benchmark. -The reaction force of the experimental data before that kink shows slightly nonlinear behavior. The crack initiation occurs u x ≈ 0.43 × 10 −3 m . These aspects are not discussed exhaustively in the benchmark. -The experimental curve shows slightly nonlinear behavior until 0.1 × 10 −3 m.
The simulated load-displacement curves are shown in Fig. 9. The noise from the explicit time integration schema was filtered using a low-pass Butterworth filter. The raw data are presented in Fig. 12 in the appendix. The aforementioned kink in the experimental curve was not present in any of the peridynamic simulations. This agrees with the findings of Gee et al. [19]. Thus, a significantly higher final reaction force is derived from the simulation compared to the experiment.
The reaction forces before the kink are slightly higher compared to the experiment ones. A further reduction in the Young's modulus will lead to better results. Assuming that the experimental curve is linear (cf. dashed line through origin in Fig. 9), the simulation fits better to the experiment and the final reaction force of the simulation differs between .
Though the damage initiation is not apparent for both the unfiltered and filtered curves. Since the effect in the experiment is also very small (small discontinuity at u x = 0.43 × 10 −3 m ), this effect is likely to disappear in the simulation due to the applied filter in Fig. 9 or the noise of the numerical time integration in Fig. 12. In summary, the peridynamic simulations reproduce the experimental findings up to damage initiation quite well. However, as by Gee at al. [19] the overall behavior of the curves cannot be reproduced, because the kink was not reproducible. The damage initiation and progression are discussed in the next subsection.  Figure 3b shows the crack propagation against the applied external displacement. It can be seen that the crack does not grow uniformly. The crack propagation might be divided into three phases: After crack initiation, the crack grows quickly, slows down afterwards and finally accelerates again. The resulting crack length-displacement curve appears to be a tangent function.

Crack Length-Displacement Curve
The peridynamic simulations are able to describe the change of crack propagation speed (cf. Fig. 10). For all solutions that form a solid crack, the initial crack propagation is followed by a deceleration and then an acceleration of crack growth. However, the simulations are only able to reproduce this behavior qualitatively. For those simulations which were stable until the full crack forms, the whole crack forming process is faster compared to the experiment. For all results, the crack length-displacement curves are not very smooth. As depicted in Fig. 3b at some point, the crack grows almost instantaneously. The reasons for this are not entirely clear.
It might be due to the explicit time integration method. The loading rate within the peridynamic simulation is significantly higher than in the experiment. This can lead to the fact that the local load redistribution during crack growth does not take place. As a result, the local crack tip stresses are increased compared to the experiment. This effect has to be studied in future work.
It might further be that PMMA in the peridynamic model lacks viscoelastic behavior [43]. Lowering the loading rate would decrease the local stresses due to local creeping [44] at the crack tip resulting in a smoother crack propagation process. However, no data are provided for the viscoelastic behavior of PMMA in the benchmark, because only one loading rate was utilized. Hence, it is difficult to say whether the results are substantially affected by this. Fig. 9 Load-displacement curve of simulation, experimental data and extrapolated experimental data without the kink 1 3

Sensitivity
To obtain more insight into the sources of divergence between numerical results and experiments, a sensitivity analysis was performed (cf. Sect. 3). The Morris method, as all screening methods, provides sensitivity measures that tend to be qualitative, i.e., capable of ranking the input factors in order of importance, but do not attempt to quantify by how much one given factor is more important than another. A quantitative method would provide an estimate, for example, of the exact percentage of total output variance that each factor (or group of factors) accounts for. As discussed in Subsec. 6.2 and visualized in Fig. 8, the overall crack path is in good agreement with the experimental results. The parameters under consideration mainly influence the crack initiation, crack length and reaction force. Therefore, these output parameters will be utilized for the calculation of elementary effects. The mean * i and standard deviation i of elementary effect per considered parameter obtained by Eqs. (3) and (4) can be seen in Fig. 11.
The Morris method provides qualitative sensitivity measures in terms of ranking the input factors in order of importance. Hence, from the results it can be concluded that almost every parameter significantly affects the chosen output parameters. Concerning the simulated crack length, the following conclusions can be made: -The loading condition m l seems to be of low influence. Since measuring this parameter during tests is very challenging, it might be neglected in simulations for the sake of simplicity. -The influence of the Poisson's ratio for this experimental setup is very low. Due to the free edges in the y-direction, the plate is able to deform freely. This results in no additional local loads due to the Poisson effect inhibition. -The most influential parameter is the horizon . This confirms the importance of an exact identification of the horizon. A larger horizon smears local effects. This can be the crack tip or the geometric resolution. I.e., a larger horizon reduces the quality of the geometric modeling of the hole. Therefore, the crack initiates earlier on the right side. At the crack tip, a horizon too large leads to a worse estimation of the stresses at the crack tip. -The remaining parameters, i.e., geometric characteristics ( l c , b c and c ) and material properties (E, G IC and ) are also highly influential and consequently vital for accurate crack prediction.
In terms of the simulated reaction force show, it can be concluded: Fig. 11 Mean * i and standard deviation i of elementary effects determined with the Morris method, input parameter description can be found in Table 2 -The loading condition m l seems to be of low influence. Hence, the crack propagation is almost invariant against variation of the load introduction. -The influence of material properties (E, G IC , and ) slightly exceeds the influence of geometric parameters ( l c , b c and c ). Determination of precise material properties should therefore be prioritized over the determination of accurate geometric properties. -Compared to the crack length effect, the Young's modulus E calculated standard deviation i is much smaller. This implies less interaction with other parameters and less nonlinearity, respectively.
In conclusion, the sensitivity analysis based on the Morris method showed that the influence of loading conditions can be neglected. In addition, it became clear that the determination of precise material parameters and peridynamic parameters such as the horizon is of high importance. Overall, the radar chart in Fig. 11 confirms that the peridynamic model framed in this paper allows a physical description of the problem. Nonetheless, according to Fig. 11 all parameters show large standard deviation indicating significant interaction and nonlinear effects, respectively. This needs to be investigated further in order to provide more qualified conclusions.

Conclusion
In this paper, a benchmark was calculated to verify a correspondence energy-based damage model. An open data set was used for the benchmark. As is often the case, the documentation of open source datasets is not sufficient. In addition, only one experiment was conducted, so scatter could not be measured. One goal was to extract as much value as possible from such a dataset despite these limitations. In order to make reasonable use of the data and to generate added value, a sensitivity analysis was performed. An OAT approach, i.e., the Morris method, was used. It was shown that the modeling approach provides a relatively robust prediction of the experimental result. The crack path could be determined robustly using the present method. No influence of varied parameters on the crack path was found. The simulated crack propagation qualitatively matches the experiment. Although, the crack length plotted against the externally applied displacement changes at different rates. This can be observed in the experiment also. However, the overall damage process of the peridynamic model grows too quickly. The time of crack initiation (external displacement) can be determined within the parameter variation. The associated reaction force, however, is not. This is mainly because of a kink in the experimental load-displacement curve. The source of this kink was not explained in the benchmark and another publication working with this benchmark could not solve this issue [19]. If this kink is ignored, the crack initiation time and the associated force match with the experiment and the publication of Gee at al. [19]. In summary, it can be said that the peridynamic modelling is qualitatively able to describe the benchmark provided by Julien Réthoré [18]. With the proper model fitting, quantitative description might be also possible. This says nothing about the quality of the peridynamic or any other model. There are some lessons learned for a realistic benchmark:

3
-Material data must be determined in independent experiments. In the benchmark presented, at least three material parameters have not been given (Young's modulus, Poisson's ratio, energy release rate). -The experimental design and setup of the experiment are barely described. The realization of the boundary conditions leaves room for interpretation. For example, the load-displacement curve bends before damage sets in. It is not clear why this is the case. A better description of the experimental setup would be helpful. -The result data are not presented in a uniform way. Results are presented in steps and/or displacements. The interpretation of when exactly a crack initiation occurred is up to the person who wants to simulate the experiment. -Typically PMMA shows viscoelastic behavior [43], which may affect the crack growth. To test plastic or viscous effects, the experiment should be unloaded in the loading process to measure residual deformations. -The significant nonlinearity of the force-deformation curve was not explained as a phenomenon in the benchmark.
However, to further qualify the peridynamic model, the experimental basis needs to be improved. A simulation-oriented experimental procedure is essential for a reliable conclusion about the model quality, in terms of narrowing the room for interpretation and increasing comparability between different methods.