A Review of Benchmark Experiments for the Validation of Peridynamics Models

Peridynamics (PD), a non-local generalization of classical continuum mechanics (CCM) allowing for discontinuous displacement fields, provides an attractive framework for the modeling and simulation of fracture mechanics applications. However, PD introduces new model parameters, such as the so-called horizon parameter. The length scale of the horizon is a priori unknown and need to be identified. Moreover, the treatment of the boundary conditions is also problematic due to the non-local nature of PD models. It has thus become crucial to calibrate the new PD parameters and assess the model adequacy based on experimental observations. The objective of the present paper is to review and catalog available experimental setups that have been used to date for the calibration and validation of peridynamics. We have identified and analyzed a total of 39 publications that compare PD-based simulation results with experimental data. In particular, we have systematically reported, whenever possible, either the relative error or the R-squared coefficient. The best correlations were obtained in the case of experiments involving aluminum and steel materials. Experiments based on imaging techniques were also considered. However, images provide large amounts of information and their comparison with simulations is in that case far from trivial. A total of six publications have been identified and summarized that introduce numerical techniques for extracting additional attributes from peridynamics simulations in order to facilitate the comparison against image-based data.


Introduction
Peridynamics [113] provides a novel framework for the modeling and simulation of fracture mechanics applications. Peridynamics (PD) can be viewed as a non-local generalization of classical continuum mechanics (CCM), in which the mechanical behavior is modeled in terms of an integral equation of the displacement field rather than by a partial differential equation. The general framework for peridynamics can be further classified into two approaches: (1) bond-based PD [113,117] and (2) state-based PD [118].
In all cases, the PD formulation allows for discontinuous displacement fields, which may be seen as an attractive feature of the method for the treatment of problems dealing with fracture mechanics. Another advantage of PD is that failure phenomena, like in molecular dynamics or atomistic lattice models, can be explicitly included in the constitutive models. In other words, no other external criterion is needed within the peridynamics framework for modeling the initiation and growth of cracks, in contrast to continuum models. This implies that PD is well suited for the simulation of fracture mechanics applications initially devoid of cracks. We refer to [10,87] for a general description of PD and to [116] for a review of its advantages.
Although the PD framework appears to be a promising approach for many problems in solid mechanics, a collective effort remains to be done to validate the bond-based and state-based constitutive models in order to gain confidence in the approach's predictability and reliability. In particular, PD introduces new parameters, such as the so-called horizon, that are a priori unknown and need to be identified. Both calibration and validation processes require experimental data to assess model's adequacy and predictability [41,42]. Another challenge is the treatment of boundary conditions in the case of non-local models. For example, boundary traction does not naturally appear in the non-local formulation [88], which implies that external loads should be applied in a manner different from that in classical continuum mechanics.
The main objective of this survey paper is to review and catalog the available experimental data that have already been used in the literature for the validation of peridynamics models and simulations. Our goals are essentially to overview the experimental setups and observables that have been achieved so far and to summarize the validation results presented by the authors, if available.
The experiments are classified here into three main categories: (1) the first category is concerned with dynamical experiments that involve wave propagation phenomena; (2) the second category lists the experiments that focus on characterizing the initiation and growth of cracks for various materials; (3) the last category gathers the experiments whose outputs are digital images of observed phenomena. This classification was arbitrarily chosen to simplify the presentation of the experiments.
The experimental data has been classified by type: (1) scalar values of observables; (2) scalar functions of input parameters; (3) and digital images. We report the relative error between experiments and predictions for scalar observables. We estimate the so-called coefficient of determination R-squared [26], also referred to as the Pearson correlation coefficient R, between the data and model predictions for scalar functions, like series of points. Data stored in the form of images provide a large amount of information, which makes comparison with simulations less trivial. It may be necessary to develop advanced visualization methods [60,61,140] to take full advantage of the richness of imaging techniques and be able to fully compare the results from non-local model simulations such as those obtained from peridynamics or molecular dynamics. We think that these methods could be useful to the community and therefore report on existing methods, e.g., extraction of fragments or crack surfaces.
We hope that this survey will provide a useful resource for those interested in the validation of PD models. We also believe that such a contribution could lead to the definition of benchmark problems and experiments for the analysis of PD simulations. This has been in fact the source of motivation for writing such a paper.
The paper is structured as follows: Section 2 presents some preliminaries, in particular, the approach that we followed to collect the publications relevant to the survey and a brief description of the method to estimate the coefficient of determination R-squared. Section 3 presents the list of dynamical experiments based on wave propagation phenomena for the validation of crack propagation and Section 4 reviews the experiments used for the validation of crack initiation/propagation. Section 5 reports in a concise fashion the relative error or the coefficient of determination for each of the previous experiments that provide a confidence level in the validity of the peridynamics models. Section 6 describes advanced visualization techniques and methods to extract additional attributes for comparison against experimental data. Section 7 provides conclusions and perspectives on the state-of-the-art in the validation of peridynamics models.

Literature Search Strategy
Our objectives in this paper were to find published articles dealing with the validation of peridynamics models based on experimental data, either in the form of observational data or digital images. References for this survey were thus collected using the two databases Google Scholar and Web of Science TM based on the following search strategy:

A search with keywords "peridynamics + experiment"
and "peridynamics + benchmark" resulted in a total of 53 references, of which 39 papers directly dealt with comparison of simulation results against experimental data. Note that we chose to include a search based on the keyword "benchmark" to find works in which authors used experimental data to verify the implementation of their codes. These papers do not necessarily address model validation. Although we will concentrate on the results of the 39 papers most relevant to this study, the remaining 14 articles are also included in the list of references. 2. A search with keywords "peridynamics + computer graphics" and "peridynamics + visualization" was also carried out and resulted in five papers, among which three utilize techniques from computer graphics for the comparison with experimental data and two use peridynamics for physics-based rendering.

Choice of Validation Metric
Validation processes involve the comparison of metrics and threshold values between experimental data sets and predictions. Many metrics, either statistical or deterministic, have been introduced in the past depending on the nature and availability of data. For consistency throughout the paper, we have chosen to consider two simple metrics. Wave dispersion and propagation [133,143] [ 18] Let us define the relative error between an observable y and the predicted valueŷ as, Let us consider a series of observables y i , i = 1, . . . , n, each associated with predicted valuesŷ i that can be computed in terms of a model f (x i ), where x i denotes the predictor variables. The coefficient of determination R 2 is defined as: where the error sum of squares, SSE, and the total sum of square deviations from the averageȳ, SST, are given by: R 2 is a statistical measure ranging from 0 to 1 that indicates how the predictionsŷ i approximate the data points y i . In particular, if R 2 = 1, i.e. SSE = 0, then the data points and predictions fit perfectly. On the other hand, if R 2 = 0, i.e., SSE = SST, then all predictions lie on a perfectly horizontal line. In other words, the value of R 2 provides a measure of the correlation between the data points y i and the predictionsŷ i . Some of the references did not provide any R-square coefficient. In that case, we chose to extract the values from the experimental and simulation data by using the web-based application WebPlotDigitizer 1 and computed the R 2 coefficient using the SciPy package. 2 Note that the computed R 2 value depends on the extraction procedure used to obtain the data from the published figures. The procedure that we followed here was first to take a screen shot of the figures and upload them to the WebPlotDigitizer.
Then, we selected a range on the curve with tick marks on the x-axis and the y-axis, and provided the corresponding values on the axes. Following this calibration process, we defined uniformly distributed points along the two curves obtained from the experiment and simulation. We acknowledge that the choice of the calibration points and the equidistant points may influence the R 2 value. For a given data set, we repeated the extraction process three times and observed that the variations would be on the order of 2%. We deemed this difference as acceptable.
We wish to emphasize at this point that accessing raw published experimental data in experimental fracture mechanics is difficult. A study of the long-term availability of raw experimental data was carried out in [35]. The main result of the study was that out of the 187 articles published from 2000 to 2016, only 11 data sets were still available and 42% of the authors who were contacted did actually reply. This situation explains why we had to resort to alternatives means to analyze the published data. Table 1 summarizes the experiments involving wave propagation phenomena found in the collected publications. The experiments are grouped by bond-based peridynamics (B) and state-based peridynamics (S) to denote the type of PD formulation that attempted to reproduce the experiment. The type of material used in the experiment and the references for the simulations and experiments are listed as well. The missing R 2 coefficient and relative errors were also evaluated when applicable. Most of the figures in this section, have been adapted and extended from those in the corresponding articles. To distinguish between the measurements and the applied load, the applied load was colored in blue.

Propagation of Stress Waves in a Half-Plane
Nishawala et al. [93] studied the propagation of stress waves in an allyl diglycol carbonate (CR-39) half-plane. The sample's length, the height, and the thickness were 1000 mm, 500 mm, and 6.655 mm, respectively (Fig. 1). Fig. 1 Geometry of the half-plane [25] where a triangular pulse load p(t) of 20 μs with maximum amplitude 20.7 × 10 3 N was applied to the top and at the middle of the plate thickness 6.655 mm. The plate's material was allyl diglycol carbonate or CR-39 An triangular impulse load of 20 μs with maximum amplitude of 20.7 × 10 3 N was applied to the top edge using a semi-cylindrical pit. Nishawala et al. [93] compared experimentally measured displacements [25] in x-direction (m) as a function of the position (m), for three different times. The final simulation time was 208 μs. The R 2 correlation is 0.74 at 60 μs, 0.35 at 92 μs, and 0.15 at 139 μs. Nishawala et al. [93] concluded that additional studies for the relation between the horizon and wave speed were needed.

Edge-on Impact Experiment
The edge-on impact (EOI) experiment is designed to visualize dynamic fracture in brittle materials using a highspeed photographic technique [122]. In this experiment, a cylinder hits a square plate on one edge, at a high velocity. Figure 2 shows the dimensions of the experimental setup. The wave front velocity and the impact damage visualization for aluminum oxynitride (ALON) and polymethyl methacrylate (PMMA) [90,[121][122][123] were addressed in Diehl et al. [34] and Zhang [141].
Zhang [141] compared the simulated and experimentally measured average propagation speed of the primary damage front [90]. The average propagation speed in the simulation was estimated at about 8.7 km s −1 while the corresponding value measured in experiments was about 8.4 km s −1 . The propagation speed of individual cracks obtained from PD simulations was about 4.8 km s −1 and about 4.6 km s −1 for cracks growing in y-direction. In comparison, the value obtained from experiments was about 4.4 km s −1 .
Diehl et al. [34] compared the predicted wave speed of the damage front with that measured in [121][122][123] and reported a relative error e r of 10%. Diehl et al. [34] studied various convergence scenarios considered in [38], as shown in Fig. 3. In this work, the best results with respect to the relative error were obtained using a horizon of three times the nodal spacing, as suggested in [115]. The horizon in [115] was adjusted such that the peridynamics simulations produce the same dispersion curves (frequency (Hz) vs phase velocity of the longitudinal waves (m s −1 )) as those measured in experiments for a given material.

Split-Hopkinson Pressure Bar
PD simulations of wave propagation phenomena in a Split-Hopkinson pressure bar (SHPB) [14,23] were performed by Jia [64] in the case of steel. SHPB is usually used for determining the dynamic stress-strain response of metals [57,76]. Figure 4 shows a diagram of the SHPB experimental setup, which consists of three bars: the striker bar, the incident bar, and a transmission bar. The specimen is placed between the incident bar and the transmission bar. The striker bar impacted the incident bar at a velocity of 17.5 ms −1 and a wave propagated through the bars. The black dots represent strain gages used to measure the strain variations produced by traveling waves. Jia [64] compared the predicted strain as a function of time with experimental results up to 2 ms. Using the experimental data provided in [23], the R 2 coefficient is evaluated at 0.99, which is a very good fit to the experimental time series.  Foster et al. [44] also studied the Hopkinson pressure bar experiment [46]. Figure 4 depicts the sketch of the experimental setup. Their objective was to compare the predicted true strain -stress curves with the experimentally measured curves, for various strain rates. The R 2 coefficient was computed here as 1.00 and 0.97 in the cases of strain rates of 1150 s −1 and 2900 s −1 , respectively.

Wave Dispersion in Bera and Massilion (Dry) Sandstone
Wave dispersion (frequency (Hz) vs phase velocity of the longitudinal waves (m s −1 )) and propagation in two different types of sandstone (Berea and Massilion) [133,143] was studied by Butt et al. [18] for different levels of confining pressures (10, 20, and 40 MPa). Figure 5 shows the geometry of the square plate (1500 mm × 1500 mm) with a thickness of 2 mm. A time-dependent load in displacement propagating in y-direction is applied to the center (x = y = 750 mm) of the square plate.
Butt et al. [18] compared the predicted and experimental dispersion curves. For both sand stones, the experimental dispersion data measured at a confining pressure of 20 MPa [133] was used to fit the peridynamic model parameters. The R 2 coefficient between the fitted material parameters and the experimental data used for the fitting was estimated at 1.00 and 0.93 for the Massilion sandstone and Berea sandstone, respectively.
The calibrated model was subsequently used to validate the predictions made for confining pressures of 10 MPa and 40 MPa. Table 2 shows the R 2 correlation for the experimental data obtained for different confining pressures. Butt et al. [18] reported excellent agreement when the calibrated material properties for 20 MPa were used at various confining pressures. Table 3 lists the mechanical tests found in the literature that dealt with crack initiation and propagation for the validation of peridynamics models. We grouped the publications with respect to the type of materials involved, namely, composite, steel/aluminum, glass, concrete, and other materials. Table 3 also lists the type of experiments, the PD models used in the simulations, where B refers again to bond-based peridynamics and S to state-based peridynamics, and the corresponding references. The missing values of the R 2 coefficient and relative error where also evaluated, when applicable. Fig. 4 Sketch of the experimental setup of the Split-Hopkinson pressure bar (SHPB), also referred to as the Kolsky's bar. The SHPB experimental setup which consists of three bars: the striker bar, the incident bar, and a transmission bar. The striker bar impacts the incident bar at a velocity of 17.5 ms −1 , which induces a wave that propagates through the bars. The black dots represent strain gages used to measure the strain variations produced by traveling waves. In [46] the SHPB bar was extended by a pulse shaper

Dynamic crack growth around stiff inclusion
The dynamic growth of a crack around a stiff inclusion was numerically studied with peridynamics by Agwai et al. [2]. In [75], an 140 mm × 42 mm × 8 mm epoxy, embedded with a central stiff glass inclusion (fiber) of diameter d = 4 mm and featuring an initial crack, is impacted at a velocity of 5.3 m s −1 using a pneumatic hammer with a hemispherical tip. Figure 6 depicts the setup. Experiments were conducted for specimens with weakly and strongly bonded fibers. Agwai et al. [2] qualitatively compared the predicted crack paths for weak and strong interfaces in the visualized simulation results with those observed in the images during the experiments. The objective of the study was primarily to investigate the influence of the bonding interface on the crack path. The simulations qualitatively captured the different failure modes that were observed in the experiment.

Damage growth in a double-lap joint test using a six-bolt specimen
Oterkus et al.
[96] simulated the initiation and propagation of damage in a six-bolt carbon fiber/epoxy system (HTA7/6376) specimen [120] during a fatigue experiment. The fatigue behavior of a double-lap bolted joint with multiple fasteners was induced by subjecting the specimen to a cyclic tensile load σ min / σ max = −1 until failure. Figure 7 depicts the sample's geometry.
[96] quantitatively compared the simulated and predicted failure modes using photo micro graphs. Failure essentially occurred in the matrix near the bolt holes. The failure modes predicted in the simulations were consistent with the experimental observations in the images.

Progressive damage in center-cracked laminates
Kilic et al. [70] simulated the crack initiation and growth in center-cracked laminates using peridynamics and compared their results with the experimental observations in [6,69,134]. Figure 8 shows the rectangular plate utilized in the simulations of dimension 10.16 mm×5.08 mm having an initial crack of length 1.27 mm. The plate is made of four plies of thickness 0.0413 mm each so that its total thickness is 0.1651 mm. The uniform tension applied in experiment [134] is gradually applied by prescribing the velocity boundary conditions v = ±1.27 × 10 −7 mm per time step at the nodes located within a distance of 0.097 mm from both vertical ends. Note that they used a 15 times smaller inplane dimension than the one in the experiment [106], due to limitation in the computational resources. For lamina the elastic properties of T800/3900-2 pre-preg tape reported by [106] were chosen. Lamina orientations of fibers are 0 • fibers are running parallel to the x-axis and 90 • fibers are running parallel to the y-axis.
were considered for three-ply laminates. Kilic et al. [70] obtained an asymmetric delamination pattern due to the presence of a 45 • ply.
Moreover, the multiple splitting around the crack tip was also observed in [6]. Finally, Kilic et al. [70] [12,65] was considered in Hu et al. [58]. Figure 9 shows a diagram of the 200 mm×100 mm plate with an initial crack of length 25 mm subjected Fig. 6 Sketch of the geometry with a width of 8 mm used in [75]. A pneumatic hammer with a hemispherical tip is used to induce an impact at a velocity of 5.2 m·s −1 in the center of the top of the plate to a uniform tensile load σ of ±40 Pa. The simulations showed a symmetric path of splitting fracture mode and bond breaking only in the matrix, which is consistent with the experiments [12].

Kalthoff-Winkler experiment
The Kalthoff-Winkler experiment consists in striking a plate featuring two initial cracks of length 50 mm in the middle with a steel impactor, as sketched in Fig. 10. The observable in this experiment is the crack angle, along with the fact that the cracks propagate through the free surface [66][67][68]. The crack propagates from the initial crack tip at an angle around 68 • with respect to the initial crack direction.
Silling [114] correctly predicted the crack angle with respect to that measured in the experiment and observed that the cracks propagated all the way to the free surface. Amani et al. [3] and Zhou et al. [144] reported that the crack angle predicted by the state-based PD model matched the angle observed in the experiments. Gu et al. [52] also showed that the crack angle and final crack path were qualitatively in good agreement with the experimental observations, even if the crack initiation time obtained in the simulations, 28 μs, was Sketch of the geometry (top view (a) and front view (b)) of six-bolt composite specimen [120]. The specimen plates were joined by six titanium bolts with protruding hexagon heads slightly lower than that measured in the experiments, 29 μs.

Compact tension experiments
Compact tension (CT) tests are used to study fracture in metals [89,91,131] and are standardized by ASTM E647-00 [5] and ISO 7539-6 [62]. Figure 11a shows a plate specimen of dimension 126 mm×121 mm with four circular cutouts. For this configuration, Yolum et al. [135] simulated the crack mouth opening (in mm) with respect to the applied force (in kN) and compared their results with the experimental results of [77,89]. For aluminum alloy (D16AT) CT specimens with four and eight circular cutoffs, the R 2 coefficient was computed as 1.00 in both cases. For a plate devoid of cutoff, the R 2 coefficient, using the experimental results provided in [131], is 1.00. Zhang [141,142] simulated the crack paths observed in [91] on a modified CT test [142], whose configuration is shown in Fig. 11b. The dimensions of the plate are in this case  [70] which is smaller than in the experiments due to computational limitations. The uniform tension applied in experiment [134] is gradually applied by prescribing the velocity boundary conditions v = ±1.27 × 10 −7 mm per time step at the nodes located within a distance of 0.097 mm from both vertical ends 40 mm×40 mm×8 mm. Subsequently, Zhang [141] simulated the position (X, Y ) (in mm) of the crack path and the fatigue life for three values of the PD horizon, i.e., δ = {0.6, 1.2, 2.4} mm. In the case of the modified CT test, the R 2 coefficient for the crack path positions and fatigue life is given by R 2 = 0.99 and R 2 = 0.99 for δ = 0.6, and R 2 = 0.97 and R 2 = 0.99 for δ = 1.2, respectively. Zhang reported that the case δ = 2.4 produces a short line in the plots due to a different crack path. Therefore, this case is not considered for the R 2 correlation.

Taylor impact experiments
A Taylor impact test [4,21], in which a cylindrical projectile is shot on a hard target, was considered by Amani et al. [3] and by Foster et al. [43,45]. The experiment is designed for low velocity impacts in order to induce relatively small deformations. Amani et al. [3] qualitatively compared the post test deformation using dynamic structural light (DSL) images from [4] with the deformation obtained from their simulations. Moreover, Foster et al. [45] compared the normalized length and diameter of the damage as a function of the impact velocity for the same experiment.  The coefficient of determination was estimated as R 2 = 1.00 and R 2 = 0.99 for the normalized length and normalized diameter, respectively. The true stress, as a function of the Lagrangian strain, was also estimated and compared to the experimental results of [21]. The R 2 coefficient is in this case equal to 0.96.

Ballistic impact test
A ballistic impact test was simulated by Tupek et al. [127] using state-based peridynamics and compared to the experimental results of [132], in which a hard-steel spherical projectile of diameter 13.97 mm, impacted an extruded 6061-T6 aluminum sandwich panel (Fig. 12) with a velocity of 370 m s −1 , 530 m s −1 , or 900 m s −1 . The observable in the experiment was the residual velocity versus the impact velocity. The predicted values of the observable lead to a R 2 coefficient of 0.99.

Crack paths in a quenched glass plate
Experiments with single crack path [13,103,136] and multiple crack paths [102,137] in a quenched glass plate were numerically studied by Kilic et al. [71]. In the experimental setup, a thermal stress was applied to the plate by placing it between two heat reservoirs. Depending on the specimen size and temperature difference, three types of crack propagation path were observed in the experiments: straight crack, oscillating crack, and branched crack as shown in Fig. 13. The rectangular thin plate of dimension 24 mm×64 mm×0.4 mm featuring a single initial crack, similar to the one used in the experiment [136], is shown in Fig. 14a. Kilic et al. [71] predicted straight cracks for small temperature differences and observed that the crack would oscillate first and then propagate in a straight direction for intermediate temperatures, and that crack branching would occur at higher temperatures. In other words, the PD simulations were able to capture most of the observed crack behaviors with respect to the temperature differences. Figure 14b shows the plate of dimension 48 mm×64 mm×0.4 mm for the setup with multiple cracks. Kilic et al. [71] numerically studied a plate with 2, 3, 5, 6, 7, 10, and 16 initial cracks and observed that there is an upper limit on the number of cracks that can propagate inside a specific plate. This is consistent with the experimental results in [102].

Crack propagation speed in a pre-cracked glass plate
Gu et al. [52] used the bond-based PD model to study the propagation speed of a crack in a pre-cracked glass plate subjected to a step tensile loading. Gu et al. [52] considered the plate of dimension 100 mm×40 mm with an initial crack of length 50 mm, as shown in Fig.  15. Gu et al. [52] computed the crack propagation speed (in ms −1 ) as a function of time (μs) and compared it against the maximal fracture velocity v F = 1580 m s −1 obtained in the experiment in [15]. The maximal propagation speeds obtained in the simulations were as Fig. 11 Sketches of the compact tension (CT) tests used to study fracture in metals [89,91,131] and standardized by ASTM E647-00 [5] and ISO 7539-6 [62]. a CT tension test with four holes and b modified geometry

Crack growth in a pre-notched glass sheet
Crack growth in a pre-notched glass sheet [100] under tensile loading was studied by Agwai et al. [2]. Figure 15 depicts the pre-notched plate of dimension 100 mm×40 mm and the initial crack of length 50 mm. The PD simulations and experiments both put in evidence that the crack splits into two branches, which grow until they eventually reach the right boundary.

Fast crack growth in an edge-cracked plate under tension
Fast growth of a crack in an edge-cracked PMMA plate under tension [39] was considered by Agwai et al. [2]. The plate had dimension of 380 mm×440 mm and the initial crack was 4 mm long, as shown in Fig. 16. The crack velocity predicted over the time period up to 80 μs was compared with the velocity obtained in the experiment, ranging from 0 to 1000 m s −1 . The corresponding R 2 coefficient was evaluated in this case as R 2 = 0.72. 5. Impact damage on a thin glass plate with polycarbonate backing Impact damage on a thin glass plate with poly carbonate backing [8,20,40] was studied by Hu et al. [59]. They qualitatively compared the damage patterns in the glass layer under three projectile speeds. Hu et al. [59] reported that the main fracture patterns observed experimentally were correctly captured.

Tensile test
A tensile test [124] was simulated by Diehl et al. [32] using a bond-based PD model. The objective was to predict the Poisson ratio of a PMMA sample under a time-dependent loading for different values of the horizon parameter m ∈ {4, 5, . . . , 12, 13}. Thus, the neighborhood B δ (x i ) of a discrete PD node x i contains 2m + 1 nodes in [x i − δ, x i + δ] in each direction [11]. The R 2 coefficient of the Poisson ratio as a function of time was obtained as 0.65. Note that the parameters of the exponential model considered in this study depended on the energy release rate and shear modulus, and not on the horizon, as in many other models. The horizon in the model was determined by the failure strain, which can be measured during the tensile test. Fig. 13 Schematics of the crack growth depending on the specimen size and temperature difference: straight crack (a), oscillating crack (b), and branched crack (c) observed in experiment [136] Fig. 14 Sketches of the single crack (a) and multiple crack glass plate (b). In the experimental setup, a thermal stress was applied to the plate by placing it between two heat reservoirs to initiate the crack growth

Lap-splice experiment
Gerstle et al. [48] presented simulation and laboratory results of a reinforced concrete lap-splice benchmark problem. Gerstle et al. [48] quantitatively compared the crack pattern obtained in the laboratory experiment and in the simulations. The crack pattern observed in the experiment consisted of two or three longitudinal splitting cracks, each one being approximately co-planar with the axis of the concrete cylinder. However, the crack patterns predicted by the PD simulations were quite different and the authors suggested that such a difference could be explained by variations in the loading rates.

3-point bending beam experiment
Gu et al. [51] performed some simulations of a 3-point bending beam experiment [63]. Figure 17 illustrates the plate (320 mm×70 mm) used in the experiment. The velocity applied at the top of the beam Fig. 15 Sketch of the geometry of pre-cracked glass plate under step tensile loading [15,100] was v = 0.075 mm min −1 . The objectives of the study were twofold: (1) to accurately predict the crack pattern observed in the experiment [63], which was actually the case; (2) to predict the crack mouth opening displacement with respect to the load. The experimental and simulation data give in this case a coefficient of determination of R 2 = 0.61. Aziz [7] performed simulations of the crack mouth opening displacement versus the applied load (CMODload) for three different configurations (D3, D6, and D9) as shown in Fig. 18 and Table 4 [19]. The values of the R 2 coefficient are collected in Table 5. Here, we observe that the values of R 2 vary depending   [19,63] with an initial crack. The velocity applied at the top of the beam was v = 0.075 mm min −1 on the length and position of the notch. The crack mouth opening displacement is usually insensitive to concrete damage at the loaded point and support. Alternatively, the load point displacement, a quantity that is very sensitive to the damage in concrete at the load points, was also computed with respect to the applied load (LPD-load). We observe in Table 5 that the coefficient R 2 suffers a non-negligible drop for D9. In a different exercise, the damage model was calibrated with respect to the experimental CMOD-load diagrams and simulations for the CMOD-load were repeated. The new values of R 2 are also reported in Table 5.

Failure in a Brazilian disk
Failure in a Brazilian disk in compression [54] was considered by Gu et al. [51]. Figure 19 shows the diagram of a Brazilian disk of diameter d =100 mm with an initial crack of 30 mm at the center of the disk subjected to a compressive velocity boundary condition v = 0.05 m s −1 . The objective of this problem was to reproduce the crack pattern observed in the experiment when the crack propagates through the disk. Gu et al. [51] reported that the numerical results were similar to the experimental ones. They also concluded that nonordinary state-based PD could be important in the future for the prediction of damage processes.

Anchor bolt pullout experiment
The Anchor bolt pullout experiment [128], shown in Fig. 20, was simulated by Lu et al. [83]. In this  Table 4  The width 152.4 mm and the notch width 2.54 mm were the same for all three configurations experiment, the concrete structure is fixed with two anchors and the bolt is pulled out. Measurements considered for three configurations [128] are shown in Table 6. For these three configurations, the peak loads were predicted and compared to the values obtained in the experiments (see Table 7). For one of these configurations, the load-displacement response (in kN versus μm) was recorded and simulated: the calculation of the corresponding coefficient of determination yielded R 2 = 0.92. Finally, the failure mode was quantitatively assessed with respect to the experimental observations and the conclusion was that the crack direction and crack branching obtained in the PD simulations were in good agreement with those observed in the experiments.

Rupture phenomena in bio-membranes
The rupture phenomenon in bio-membranes [49] was considered by Taylor et al. [125]. The objective in this experiment was to compare the fractal dimension of ruptures in actual membranes with that extracted from PD simulations. The fractional dimension was postprocessed from the computer simulations via image processing techniques and the reticular cell counting (box counting) method. The predicted fractional In the case of the CMOD-load, the damage model was also calibrated against the experimental CMOD-load data and simulations were repeated Fig. 19 Experimental setup of a Brazilian disk in compression [54]. A compression load of ± 0.05 m s −1 was applied at the lower and upper tips of the disk dimensions, namely 1.7, 1.63, and 1.56, computed for different shear moduli, were close to the fractional dimension, 1.7, obtained in the experiment.

Crushing of brittle ice against vertical structures
The PD simulation of brittle ice [119] crushed by a vertical structure was performed by Liu et al. [81]. In the experiment, a thin layer of ice of dimension 1000 mm×1000 mm×60 mm is cut in the middle by a rotating cylindrical structure. The mean force and the peak force in the ice obtained in the simulations were compared with those measured in the experiments for several values of the penetration velocity. Table 8 lists these values. We observe that the relative error strongly depends on the velocity of the cylinder.

Cross-sectional nanoindentation
Oterkus et al. [95] studied the crack patterns for cross-sectional nanoindentation (CSN) [105]. Figure 21 sketches the rectangular silicon substrate specimen with three layers of copper metallization embedded in inter layer dielectric (ILD) and etch stops (ES) [94]. In the Fig. 20 Experimental setup of the anchor bolt pullout [128] for the three different configurations shown in Table 6. In this experiment, the concrete structure is fixed by using two anchors and the bolt is pulled out Table 6 Measurements for the three different configurations [83] for the sketch shown in Fig. 20 Configuration W (mm) L (mm) a (mm) d (mm) b (mm) t (mm)   1  300  300  100  50  15  5  2  600  600  200  100  30  10  3  900  900  300  150  45  15 experiment, the silicon substrate was indented by a Berkovich diamond, which induced a V-shaped crack. Note that in the simulation, an initial 40 • V-shaped crack was introduced similarly as in the experiment. Along the normal to the crack faces, a velocity of ± 0.5 m s −1 was applied. The crack paths observed in the simulation and experiment showed that the main crack propagated through the silicon and the first layer of copper metallization embedded in interlayer dielectric.

Dynamic crack propagation in functionally graded materials
Experiments of dynamic crack propagation in functionally graded materials (epoxy/soda-lime glass) (FGM) [1,73,74] were used for comparison with PD simulations by Cheng et al [24]. The experimental setup consisted of a plate of dimension 152 mm×43 mm with an initial crack of length 8.6 mm subjected to 3-point loading, as shown in Fig. 22. First, the crack pattern observed in the experiment [74] was quantitatively compared with that obtained in the simulation. The conclusion of this study was that, despite the fact that the dynamic loading used in the simulations was different from that in the experiments, a close resemblance between the predicted and experimental crack paths could still be observed. Second, the evolution of the crack length over time [73] was considered for two cases: (1) pre-crack on the stiffer side (E 1 > E 2 ) and (2) pre-crack on the more compliant side (E 1 < E 2 ). When using the linear curve fit of elastic modulus and density to compute the peridynamics micro modulus, the R 2 coefficient is 0.99 for E 1 < E 2 and 0.99 for E 1 < E 2 . When using the piece-wise linear variation measured from the experiments, the R 2 coefficient is now given by 0.98 for E 1 < E 2 and 0.99 for E 1 < E 2 . Table 7 Relative error in the predicted peak loads [83] with respect to those obtained in the experiments [128] Configuration Simulation (kN) Experiment (kN) Relative error (%) 1 17

Summary of Validation Results for Peridynamics Modeling
We reported the R 2 values with respect to the Young's modulus of the materials to create a visual representation of the results. We emphasize that the relative error and coefficient of determination are two measures among many others to quantify the fitness of the simulations to the experimental data. The reasons to choose these measures are twofold: first, they are commonly used in the field and, second, it was relatively straightforward to compute them from the available data since most of the values were missing in the publications. Needless to say that other choices of measures could yield different conclusions. Tables 9 and 10 list the relative error and the R 2 coefficient between the experimental data and the values predicted by the PD simulations, respectively. Likewise, Figs. 23 and 24 provide the same information, but shown with respect to the Young modulus of the material. Figure 23 suggests that the relative error is consistently smaller than 20% for aluminum and steel, which shows that simulations can predict well the experimental results. For concrete materials, the relative error ranges between 10 and 30%. For glassy materials, it is between 3 and 71%. In the case of crushing-brittle ice by a rotating cylinder, the relative error lies in the range between 14 and 74% for three different rotating velocities. Conclusions in the case of series of observables, based on Fig. 24, are as follows. For aluminum and steel, the coefficient R 2 is in general greater than 0.9, which indicates a very good agreement between simulations and experiments. For concrete materials, R 2 ranges between 0.5 and 0.9 while for glassy materials, it is between 0.6 and 0.72. In the case of functionally graded materials (FGM) and for sandstone, the coefficient reaches values as high as 0.99 and 0.93, respectively. We thus observe that PD simulations usually agree well with experiments involving metallic materials, such as steel and aluminum, and sandstone.

Extraction of Additional Attributes from Non-local Simulation Data
In peridynamics, information, such as reconstructed stress or strain, is usually provided at the node level. However, quantities measured in experiments are often given at a larger scale. Here, advanced visualization techniques can be utilized to extract missing information or to post-process the solutions in order to be able to compare predictions and experiments.
We collect in Table 11 some information about applications in which visualization techniques were used to Fig. 21 Sketch of the rectangular silicon substrate specimen with three layers of copper metallization embedded in interlayer dielectric (ILD) and etch stops (ES). In the experiment [94], the silicon substrate was indented by an Berkovich diamond, which induced a V-shaped crack [105]. Note that in the simulation an initial 40 • V-shaped crack was introduced similarly as in the experiment. On the normal to the crack faces a velocity of ± 0.5 m s −1 was applied assess peridynamics simulations of fracture phenomena in solids. However, one should remember that visualization techniques involve models to approximate damage quantities, which should also be validated against experiments. For example, the crack growth velocity of the approximated crack surface was compared by Bußler et al. [17] with the one in the Kalthoff-Winkler experiment [68] in order to validate their model. The main objective of this section is to highlight the techniques that can be used for the comparison of peridynamics simulations with experiments.

Visualization of Fragmentation
Material fragmentation is an example that can significantly benefit from visualization post-processing techniques. Fragmentation can be measured in terms of histograms of fragment size or fragment momentum and can thus be straightforwardly compared with experimental data [107,129]. Littlewood et al. [79] developed an algorithm similar to the concept of fragment identification for classical finite element simulations [112] to extract fragments from PD simulations in terms of the nodes that are connected by unbroken bonds. Other attributes, like the mass or the momentum of each fragment, can also be provided. The algorithm was applied to bond-based peridynamics simulations of a brittle disc after impact by a spherical indenter inducing failure in the expanding ductile ring. For this example, Diehl et al. [29] extended a connectedcomponent labeling algorithm [104,110] to identify fragments with respect to broken bonds obtained from a bond-based peridynamics simulation, as shown in Fig. 25a. In this figure, the nodes are colored with respect to the damage attribute (blue means no damage and red means damage based on broken bonds). Figure 25b shows the corresponding fragments associated with their labels, which are the identifiers of the fragments. In both algorithms, fragments are identified based on broken bonds and the assumption that damage arises at broken bonds. However, the algorithms have yet to be validated against experimental data. In our opinion, it would thus be instructive to compare the results of these methods with actual experiments to gain more insight about the peridynamics simulations, e.g., the influence of horizon values and variations in the initial positions on the fragment attributes.

Visualization of Fracture Progression
Bußler et al. [17] developed a height ridge surface extraction approach to visualize fracture progression in peridynamics. The approach consists in approximating the crack surface and extracting other pertinent attributes, such as the fracture velocity or the crack angle, from the PD simulations. For instance, the crack angle and the crack growth velocity are provided in the Kalthoff-Winkler experiment [68]. Using only the "local" information available at the nodes, these quantities of interest are unfortunately not directly accessible from the simulations. Bußler et al. [17] demonstrated the efficiency of their algorithm on two Abbreviations: crack mouth displacement (CMD), crack mouth opening displacement (CMOD), and load point displacement (LPD) examples, namely, the thin disk impacted by a spherical indenter and the Stanford Bunny. Moreover, they validated the extracted fracture growth velocity with the one obtained in the Kalthoff-Winkler experiment (see Table 12). The fracture growth velocity was predicted as 1140 m s −1 , that is a 14 % relative error when compared to the value obtained in the experiment. A similar study using the smoothed-particle hydrodynamics method (SPH) [101] provided a velocity of 1200 m s −1 or a 20 % relative error. Nevertheless, it would be important to further analyze the influence of the initial node positions and the axis-alignment of the crack surfaces on the predicted velocity.

Physically Based Modeling and Rendering
Another interesting application of peridynamics could be for physically based modeling and rendering of visual effects [78,99]. Peridynamics, thanks to the comparative simplicity of the models, could indeed provide a good balance between artistic rendering and accuracy in describing physical behaviors. For example, Levine et al. [78] used bond-based PD for the animation of brittle fracture. Levine et al. [78] animated the impact of a projectile through a glass plate and the fracturing of the Welsh dragon/vase on the floor. Levine et al. [78] reported that as in most physics-based animation methods, tuning parameters can be a tedious task. Chen et al. [22] used state-based PD for the animation of fractures in elastoplastic solids. Chen et al. [22] visualized the isotropic brittle, anisotropic brittle, isotropic ductile, and anisotropic ductile fracture behavior of a sphere travelling through a wall. Next, Chen et al. [22] visualized the complex crack pattern when shooting a bullet through jello, the fragments obtained in the case of the Stanford bunny 3 falling to the ground, and the crack pattern in a thin sheet when pulled apart. Chen et al. [22] concluded that the horizon needed sometimes to be fine-tuned in order to get well-synchronized results.

Concluding Remarks and Perspectives
Applications for validation of bond-based PD and statebased PD were reviewed and summarized in Tables 1  and 3. Overall, we identified 22 applications for bondbased PD and 9 applications for state-based PD. In the case of applications dealing with initiation and propagation of cracks, we chose to classify the contributions with respect to the type of materials, namely composites, aluminum/steel, glassy materials, concrete, and other materials. The Kalthoff-Winkler experiment was used to assess the crack initiation time, the crack propagation speed, and the crack angle. We thus believe that the Kalthoff-Winkler experiment could be a valuable benchmark for the  Table 11 Overview of applications of bond-based (B) and state-based (S) peridynamics simulations using visualization algorithms to assess fracture in solids Animation of brittle fracture [78] Fracture animation in elastoplastic solids [22] Fracture progression [17] Extraction of fragments [29,79] validation of PD models, since it provides three different types of measurement at once. In other words, it is challenging to simultaneously predict the three observables. We also indicated, to some extent, whether the comparative studies between simulations and experiments were quantitative or qualitative. In most cases, we were able to compute either the relative error, for scalar observables, or the coefficient of determination R 2 , for series of observables, to assess the predictability of the PD models.
The major results of this investigation are collected in Tables 9 and 10 and plotted in Figs. 23 and 24. We note that the best correlations were obtained for edge-on impact (EOI) experiment experiments involving aluminum and steel materials.
It is commonly accepted that the size of the PD horizon should be set to three or four times the mesh size. Nevertheless, some authors have studied the influence of the horizon δ and mesh size on the predictions. For example, Zhang [141] investigated the influence of the horizon on the position of the crack path: the values of R 2 were estimated at 0.97 and 0.99 for δ = 1.2 mm and δ = 0.6 mm, respectively. Gu et al. [52] looked at the effect of Table 12 Comparison of the fracture growth velocity obtained in the Kalthoff-Winkler experiment [68] with that predicted by peridynamics (PD) [17] and smoothed-particle hydrodynamics (SPH) [101] Setting Fracture growth velocity m s Applications of computer graphics have usually two objectives: (1) to gain more realistic visual effects using physics-based modeling and rendering; (2) to reconstruct from the solutions large-scale quantities or features, such as number of fragments and shape of the crack surfaces. In all reconstruction approaches, the approximation of damage in PD is based on the notion that damage arises at broken bonds. Unfortunately, validation of these models against experimental data is still lacking at this time. There is maybe one exception in which extraction of crack surfaces in the Kalthoff-Winkler experiment was used with partial success to measure the fracture growth velocity. We believe that advanced visualization techniques could be useful to gain more insight on the effect of initial positions of the discrete PD nodes on quantities of interest, such as the fragment size or alignment of crack surfaces. They could also help investigate optimal ratios between the horizon and mesh size. Fig. 25 Illustration of fragmentation of a thin disk following impact by a spherical indenter using an extended connected-component labeling algorithm [29] Following this comprehensive study, our opinion is that there remain many issues and challenges in order to fully assess peridynamics models with respect to experiments: -Boundary conditions: One of the major challenges pertains to the treatment of boundary conditions when dealing with non-local models [37,50,85,86]. This is an issue that one also finds in molecular dynamics simulations or in smooth particles hydrodynamics. To overcome this issue, several approaches for coupling PD with the finite element method (FEM) or classical continuum mechanics (CCM), e.g., force-based blending [108,109], Arlequin approach [55], or discrete coupling [47,72,82,84,111,138,139] have been developed. One could use a coupling approach in order to enforce boundary conditions, that is, by using CCM or FEM along the boundaries and PD in the region of interest where cracks and fractures arise. -Discretization parameters: For the discrete case, the choice of the nodal spacing and the horizon influences the various convergence scenarios [38,126]. One issue is defining the proper ratio between the horizon and mesh size as numerical results are certainly sensitive [30] to these two parameters. One approach is to adjust the horizon, so that the peridynamic results produce the same dispersion curves as those measured for a specific material [115]. Another approach is to determine the horizon by Griffith's brittle failure criterion, based on critical failure stress which can be found experimentally [32]. -Computational costs: The computational costs for non-local models, like PD, are significant. Some cited references were not able to simulate the full size of the experimental geometry, due to computational limitations. Several massive parallel peridynamics implementations are available: Peridigm [80,97] and PDLammps [98] based on the Message Passing Interface (MPI), peridynamicHPX [31] based on the C++ standard library for parallelism and concurrency (HPX) [56], acceleration card codes based on OpenCL [92] or on CUDA [27,33]. -Calibration versus validation: Simulations are usually compared to experimental data for calibration/fitting or validation. In the first case, the material parameters are calibrated, such that they fit observable quantities in a given experiment, e.g., reproduce the same dispersion curve or the crack growth velocity. In the second case, the calibrated parameters are used for validation against other experimental scenarios, e.g., different loading values or different positions of the initial crack. In the majority of the cited references, the material parameters were calibrated against one experiment and only a few addressed validation.