Experimental and numerical study of the interaction between dynamically loaded cracks and pre-existing flaws in edge loaded PMMA specimens.

Dynamic fracture tests are carried out for four groups of hole-containing edge loaded specimens. The crack growth velocity, crack path, and dynamic toughness are extracted from the experiments using high-speed photography and digital image correlation. The importance of the interaction between the in-coming stress wave and the pre-existing hole is revealed and analyzed. A micromechanical damage model is calibrated to the experimental data from two of the specimens' designs and evaluated for its predictive capabilities using the other experimental configurations. The studied model is shown to be in reasonable agreement with the experimental data, and its limits are discussed.


Introduction
The study of brittle crack growth and arrest under dynamic loading conditions has posed an interesting challenge to the scientific community since the early days of fracture mechanics. Brittle (or quasi-brittle) materials are frequently used in impact-loaded components. Dynamic crack growth in brittle materials is largely dominated by the interaction of the growing crack with pre-existing flaws, such as voids and micro-cracks, as well as their evolution under the applied load and the interactions among them [1,2]. The interaction of the waves radiating from the propagating crack tip with existing heterogeneities makes the dynamic fracture process highly sensitive to variations in the local microstructure [2][3][4][5][6][7]. Having an understanding of the interaction of the growing crack with its surroundings may assist in developing a better design scheme considering energy dissipation, directional strength, and control of the fracture path. As closed-form analytical solutions are scarce and limited, experimental and numerical approaches are usually used to gain insights into this complex problem [8][9][10][11][12][13][14][15].
Experimental studies on crackflaw interaction for dynamically propagating cracks dates back to the 1970s [16]. Kobayashi et. al. [16] have studied the role of holes as crack arrestors for dynamically propagating cracks, and observed that the crack arrest capability of a hole is strongly correlated with the amount of strain energy released when the crack penetrated the hole. Once the crack has penetrated the hole, the strain and kinetic energy are gradually developing into a new local stress distribution which will trigger a new crack emanating from the crack arrestor. As noted by Milios and Spathis [17], even holes lying away from the predicted crack trajectory may attract a growing crack and momentarily lead to crack arrest. More recent studies, such as the work of Yang et. al. [18][19][20][21] on crack-void interactions in dynamic fracture scenarios have focused on modifying the conditions leading to crack arrest. Using the caustic method, Yang et al. analyzed the stress intensity factor (SIF) and crack propagation velocity during the fracture of a three-point bending specimen with a void and proposed an empirical formulation for the fracture parameters. Yang et. al. have considered multiple scenarios by varying the location and radius of the void and formulated the influence of the void on the fracture process. Crack propagation velocity was observed to increase along with a decrease in the SIF as the crack approaches the void. In [22], the effect of an inclusion's stiffness on its interaction with a rapidly growing crack was studied experimentally and was found to correlate with the degree of mode mixity evolving during the fracture process.
The role of the experimental efforts briefly mentioned above goes beyond just observing trends and measuring crack arrest time deflection angles, etc. Such experimental campaigns, serve as a backbone for modern computational methods, as they allow for complex scenarios to be validated against experiments in a clear, detailed manner [23,24]. Brittle and quasi-brittle material models have been developed to study the dynamic fracture of brittle materials based on the cohesive element method [25][26][27], peridynamics [28,29], and phasefield modeling [30]- [33], [34]. These numerical methods have different approaches to capture the sensitivity of fracture in the local field near the crack tip, where the influence of material softening and the heat release may be crucial [35,36]. Thus, to validate a model's suitability for design purposes, it has to be tested against several scenarios, while calibrated only on a minimal set of experiments [23].
The present study is a detailed analysis of the crack-hole interactions for dynamically loaded PMMA fracture specimens. A Hopkinson bar apparatus is used to facilitate crack growth using reverse tapered double cantilever beam specimens [37]. Taking advantage of the well-predicted crack path in the unperturbed configuration, cylindrical holes are pre-machined into the crack path and the crack-hole interaction is studied using high-speed photography (up to 2x10 6 fps and 924x768 pixels). Before conducting the experiments, the specimens are coated with a random speckle pattern and the Digital Image Correlation technique is used to extract the evolution of stress intensity factor during crack growth along with the crack velocity evolution.
The rest of the paper is organized as follows: In section 2, the experimental apparatus, specimens' design, and data reduction techniques are detailed. Section 3 contains a summary of the experimental results obtained for 4 specimen's designs where the SIF evolution, crack velocity, crack arrest time, and crack deflection angle are summarized for each configuration of pre-existing holes. In section 4, linear-elastic finite element calculations are used to evaluate the effect of the pre-drilled holes on stress wave propagation in the specimens during the loading stage. In Section 5, a damage model is calibrated against one of the specimens designs used for this study, and its prediction capabilities are evaluated using the other sets of experiments, showing a general agreement between the analysis and experimental data. Finally, the results and their implications are summarized in Section 6.

Specimens' design
The material chosen for the presented investigation is PMMA. The specimen geometry selected for this study takes after the Reverse-Tapered Double Cantilever Beam (RT-DCB) specimen originally proposed for quasistatic fracture tests due to its stable crack growth and predicted crack path [38,39]. A modified version of the RT-DCB was later proposed by Chen et. al. [37] for the dynamic fracture of brittle materials ( Figure 1). The main criteria for choosing this geometry was the ease of specimen alignment which will allow us to study the crack-holes interaction in a rather consistent manner by minimizing errors arising from the activation of a mixed-mode scenario. Furthermore, it was shown in [37] that the resulting crack path is very predictable which allows us to place the holes along a known path. In Chen et. al. [37], the specimen was notched with a blunt notch and no pre-crack was introduced, similar published experiments on the dynamic fracture of PMMA have previously used the same approach [40,41]. In the experiments presented here, a sharp pre-crack was introduced to the already notched specimens by inserting a razor blade into the notch (0.7mm wide) and lightly impacting it, resulting in sharp cracks of approximately 0.4 ± 0.1 mm in length. Unfortunately, this method of pre-cracking did not always introduce a straight crack front, and the small deviations in crack length and crack front curvature led to variations in the measured fracture properties.
3 Figure 1. Schematic of the RT-DCB specimen. A sharp pre-crack is introduced at the tip of the machined notch.
As stated in the introduction, the effect of pre-existing defects on the measured crack growth resistance and crack growth initiation is the result of several factors: i. The incoming stress wave is scattered by the existing holes, resulting in a different loading condition per defect geometry, thus it is expected that different holes arrangements will affect the crack growth initiation toughness and the rate at which the K dominant filed is generated; ii. As the crack grows, stress waves are emanating from the crack tip, interacting with the defect lying ahead of it; iii. Once entering the holes, the crack is expected to arrest and only continue to grow once enough strain energy has been accumulated in the specimen. The stress field leading to this scenario is also influenced by the presence of existing defects in the specimen.
Four specimen designs were chosen to illustrate the interaction of a growing crack with pre-existing holes (Table 1, Figure 2). For each design, several specimens were tested and the 3 specimens showing the smallest errors resulting from the data reduction procedure are presented. Design parameters of all 4 sets are tabulated in Table 1. The fracture experiments were analyzed by extracting the initiation toughness along with the crack growth resistance over time using a high-speed camera and digital image correlation (DIC). Similarly, the crack growth velocity, the duration of crack arrest in the holes, and the deflection angles ( Figure 3) upon exiting the holes were extracted.

Experimental apparatus
Following [37], the specimens were loaded using a Hopkinson bar in the 1-bar configuration (19.7 mm diameter, C300 Maraging steel bar). A schematic of the loading apparatus is shown in Figure 4. The specimens are positioned such that the apex is fully in contact with the incident bar. The striker is launched at a predetermined velocity, kept constant (as possible) throughout all the experiments. The strain gauge glued on the incident bar, served to measure the loading pulse to be entering the specimen, while at the same time triggering a high-speed camera (Kirana, Specialized Imaging) which is set to start imaging at a delay corresponding to the time it takes the stress wave to travel from the strain gauge to the specimen. 180 images at a resolution of 924x768 pixels are then taken at a frame rate of ~6

Crack tip identification and data reduction
Prior to the experiments, the PMMA specimens were speckled using acrylic spray paint. The resulting black and white patterns were then used to extract the displacement fields from the images taken by the KIRANA. The displacement field was extracted using the digital image correlation open-source code NCORR [42]. The DIC parameters slightly varied between experiments but nominally kept constant, where the subset radius for image correlation is 10 (pixels) with zero subset spacing. High strain analysis is used for the correlation while enabling the subset truncation. To extract the stress intensity factors (SIF) from the experimental data, an algorithm was developed and implemented in MATLAB to detect the crack tip and calculate the corresponding SIF for each frame in the experiment.
To assist in the crack tip detection a 2 steps procedure was followed. In the first step, the final image, where the crack has transverse the entire frame, is used to restrict the possible locations of the crack tip. Next, the obtained displacement field is transformed into the Fourier space to be further analyzed and compared to the analytical filed.
i. Crack path detection A window limiting the possible locations of the crack tip is determined using the strain map in the y-direction as shown in Figure 5. ii.
Evaluating the crack tip location and SIF. The crack tip is scanned along the crack path using the displacement field from the DIC. The experimental displacement field at an instance is compared with the analytical displacement field developed using the linear elastic material model.   Where ρ is the density of the material and the factor k is determined based on whether the system is plane strain or plane stress. Plane stress analysis is relevant to this study since the displacement fields are observed from the material surface.
Using the displacement formulation ( , , ) n F x y c is evaluated and compared with the experimental displacement to estimate the corresponding Ki for each term. This process is carried out for the selected possible crack tips as origin then the overlap of the analytical and experimental displacement fields are quantified. The best overlap or least error in the overlap corresponds to the crack tip of that instance.
The method chosen for the comparison of the displacement fields to obtain the SIF is a modified form of the procedure proposed by Roux and Francois [44]. Roux and Francois introduced a way to consider the overall displacement (or stress) field as a superposition 8 basis fields including the rigid body motion and multiple fracture modes. Using this approach, the SIF of multiple modes and coefficients in the power series expression could be evaluated.
Since the number of terms in the displacement expression increases the analysis time, the terms are reduced by examining the data from the numerical analysis. Although the accuracy increases with the number of terms, the process becomes unstable and sensitive to the noise. A Finite element analysis of a similar system with the Hopkinson bar is carried out with a striker velocity of 19.6m/s. The displacement of the domain near the fully developed mode I crack tip is taken for the analysis. As shown in Figure 6, the accuracy of SIF and the error has saturated with two terms in the displacement expression. So, the first two terms of the displacement power series expression are used for the numerical analysis to reduce the processing time while having an accurate outcome. The experimental and analytical displacement fields are given in Figure 7. The displacement in the cartesian coordinate system is converted into the Fourier space during the analysis, thus allowing for better noise removal and accuracy. Instead of comparing the individual data points which contain the information about one location, the characteristic frequencies which contain a trend in the domain are examined. Each point in the Fourier space contains a frequency and an amplitude that contributes a trend to the displacement field. The main advantage of implementing the analysis in the Fourier space is efficient noise removal and it is attained by looking at the analytical displacement field and considering only the relevant frequencies for the comparison. So, the process implicitly removes the random noise which is dominant in the DIC data.
A set of points along the crack path ( Figure 5b) is considered with the corresponding analytical displacement field and they are overlapped with the experimental displacement field. The precision of the overlap is quantified with the error in mapping and it is assessed (Eq. 2) using the least square method. The error is normalized by the number of data points in the displacement field. The variation of overlapping error is examined along the crack path to obtain the least error which corresponds to the best overlap ( Figure 7a). This process of identifying the crack tip from the overlap of displacement fields is carried out for each frame during the fracture and the corresponding SIFs are evaluated. Once the crack tips are evaluated, the crack propagation velocity is found using the crack tip locations and the time interval between the frames.
Although this method is efficient and accurate, crack tip detection is challenging if the displacement field is not fully developed or the domain of interest contains irregularities such as holes or bad speckles. In those cases, the analysis is corrected and carried out by comparing the crack path manually with the experiment data. Due to the finite resolution of the image, there exists an uncertainty associated with the crack tip which leads to a minor vagueness in the evaluation of SIF.

Results and observations
The experimental procedure and data reduction described in the previous section were carried out for multiple specimens from each design while maintaining the velocity of the striker bar at 19-20 m/s for all the cases. All of the examined specimens in the same design-set exhibited similar trends and the best three (in terms of DIC and image quality) are presented here.

Design set 1
As tabulated in Table 1, this set contains a hole of a 2mm radius lying on the symmetry plane of the specimen. Three fractured specimens are shown in Figure 9. Despite preserving the symmetry of the specimen and loading, in all the cases the crack is observed to deflected to one direction after the crack arrest in the hole. The impact stress wave consists of an axially compressed wave front in the x-direction, leading to a tensile wave in the ydirection (see [37] for a full discussion). After the crack arrest, the interaction of compressive and tensile waves in the specimen compels the crack to deflect along the lateral direction for a short period before it re-aligns in the horizontal direction. In case 1( Figure 9a) the crack is deflected downwards, due to slight misalignment of the pre-crack and it enters the hole below the hole's center. The SIF and the crack propagation velocity during crack propagation are shown in Figure 10. The observed variations in the critical SIF, K1D are due to the sensitivity of the crack growth initiation to the length and shape of the pre-crack. The crack velocity is evaluated considering 23 neighbor data points (frames from the DIC) with linear regression. Error limits in the crack propagation velocity are plotted for a 95% confidence interval. Error margin in the SIF is evaluated by the shortest distance to the closest neighbor.   Critical SIF, K1D varies between cases, due to the sensitivity to the initial sharp notch and slight variations in the hole which affect the scattering of the incoming wave. The measured K1D and K̇1 shows a monotonic relation with a similar trend to previous measurements in the literature ( Figure 20), even though the values are spread as previously discussed. The average crack arrest time in the hole is 23.3μs. Specimen 1 exhibits the lowest crack arrest time, since the crack entered the hole at an angle, causing a rapid stress accumulation at one point of that hole circumference. Deflection angle in all specimens stays in reasonable margin with an average deflection angle 74 0 .

Design set 2
For design set 2, the hole's radius was decreased from 2mm in design set 1 to 1mm. The fractured specimens of this set are shown in Figure 11. As before, the sensitivity of the fracture process to the initial pre-crack leads to some scatter in the results. Similar to design 1, SIF of all three specimens increase just after the crack initiation, but at a lower rate ( Figure  12a). Unlike in Design set 1, both SIF and crack propagation velocity saturate after ~10 mm of crack propagation from the notch, but the saturated velocity is lower for the 1mm radius holes. When compared to design set 1, the K1D value is higher for all the cases which leads to higher initial crack propagation velocity.  Once arrested, the crack spends less time inside the hole when compared to design set 1 and deflects at a smaller angle. The average crack arrest time is 7.3μs, which almost 1/3 of the time found for design set 1. The average crack deflection angle is 41 0 (this is due to the larger deflection observed for specimen 1). The angle at which the crack enters the hole highly influences the crack arrest time and the crack emerging angle. Since the hole radius is small, the crack entering angle is more sensitive than in design set 1. The K̇1 values, appear to be higher than in design set 1, which we attribute to the smaller disturbance to the stress wave created by the smaller hole.

Design set 3
Design set 3 takes after design set 1 with an introduced asymmetry. The 2mm radius hole is offset by 1mm from the horizontal symmetry plane of the specimen ( Figure 13). Here, the crack growth velocity is accelerating in an almost linear manner up to reaching the hole (Figure 14b). The major difference with respect to the symmetrical case is the lack of saturation in velocity when approaching the hole and the decrease in SIF which is more apparent here. As before, the specimen in which the crack penetrated the hole at an angle exhibit the smallest arrest time (Table  4). Also, for this specimen, the crack did not realign itself to grow horizontally after exiting the hole and continued to grow at an angle.

Design set 4
This set of specimens consist of 3 symmetric holes with a 2mm radius spaced by a uniform gap of 10mm between them. As shown in Figure 15, in 3 out of 4 specimens crack diverts after the 2 nd hole. The evolution of SIF and crack growth velocity are shown in Figure 16. For all specimens 1-3, the initial crack growth velocity is smaller than the one observed in previous designs. For all three cases, the SIF at initiation is smaller than in previous designs when considering the rather large loading rate. It is important to emphasize that the interaction of the stress waves with the holes, which are now closer to the initial crack tip, generates stress field with stronger compressive stress in the x-direction thus affecting the conditions at the crack tip. A significant increase in the crack growth velocity and SIF are visible in all 3 cases, before the crack entering the 1 st hole. In specimen 2, crack initiates at an angle then regains the horizontal path, and this causes the dip in halfway of the velocity variation. The increase in the SIF after the 1 st crack arrest is caused by the higher K1D of the hole compared to the sharp notch. The perturbation to the stress wave front is amplified in this design with 3 holes. Nonetheless, these specimens exhibit similar behavior to design 2 until the crack reaches the 1 st hole. The time duration of the incoming stress wave was taken to be long enough so that the hole is still under mode I during the 1 st crack arrest. Therefore, the crack continues to travel in the horizontal direction after being arrested and the measured deflection angles are rather small ( Table 5). The crack arrest time is observed to be much smaller than in the previous cases, which was expected as the hole is placed much closer to the initial crack, at a position where the magnitude of the tensile load is higher. The crack arrest time in the 2 nd hole, which is located at the same distance from the pre-crack as the hole in designs 1-3, are rather similar and revolve around 30μs. Similarly, the deflection angle obtained is close to the one observed for design set 1. KID and the crack arrest time are significantly higher at the second hole (Table 6) where a transition in stress state is observed. Prior to the crack entering the hole, the overall stress state is tensile in the y-direction, however, due to the interaction with the incoming compressive wave, there exists a stage where the hole only experiences compressive loads for a few microseconds before tension is reestablished. In Figure 17a, the compressiontension transition is shown to occur 100μsec after the incident bar has impacted the specimen. The transverse stress component σyy, which enables the horizontal crack propagation in the specimen, consists of tension and compression segments (Figure 17b). This complex scenario, resulting from wave interactions with the crack the holes and the specimen boundaries, greatly affects the time it takes for the kinetic and strain energy to redistribute themselves in from of the hole and regenerate the stress singularity required for propagating the crack further. When the crack reaches the 2 nd hole, the stress field around the hole is not strong enough to initiate the crack hence the crack must "wait" for the stress accumulation to attain the threshold. The deflection of the emerging crack from the 2 nd hole supports the compression-tension transition in the local stress field as seen in design sets 1-3. This process increases the crack arrest time in the 2 nd hole compared to the 1 st hole by 3 times. The crack deflection angle is smaller compared to design 1 which might be due to the higher kinetic energy which is the result of the higher KID. Crack propagation velocity stays the same after the 1 st crack arrest. Meanwhile, the SIF keeps dropping at this stage, an observation we attribute to the local change in the stress field mentioned above.
In design 4, 1 out of 4 sets showed a different crack path from others while having the same striker velocity (19.17m/s). As shown in Figure 18, crack keeps a straight path after the 2 nd hole. The symmetric nature of the fracture is noticeable from the crack path in the specimen. This specimen has nearly perfect horizontal crack path before and after the 1 st crack arrest, hence not triggering any asymmetry in the specimen and keeping the crack trajectory horizontal after the 2 nd arrest. The characteristics of SIF and crack propagation velocity until the 2 nd crack arrest are similar to the rest of the cases in the set ( Figure 19). The magnitude of SIF and crack propagation velocity are very close from the 2 nd hole to the 3 rd hole of design 4. K1D upon existing holes 1&2 are very close.  As shown in Table 7, crack arrest time in hole-1 and 3 are very close and 1/3 rd of the arrest time of hole-2 due to the compression-tension transition in the impact stress wave. The crack propagation initiating from hole 1 and hole 2 are similar in terms of the K1D, SIF, and crack propagation velocity which might be due to the identical geometry and the unaltered local stress field. The crack arrest time at the 2 nd hole is comparable to the rest of the cases in this design set (Table 6).

Wave front perturbations and the resulting mode I stress field
A comparison of the experimentally observed values of K1D and K̇1 with previously published data is given in Figure 20. while the K1D and K̇1 values are not, strictly speaking, colinear, they do exhibit a linear trend with log(K̇1) and K1D (Figure 20), however, the linear behavior observed in the experiments appears to be at an offset from the previously published values [37,40,41,45]. Given the fact that previously reported data points were obtained for notched specimens without a pre-crack, the observed offset was expected. The sharpness of the notch reduces the critical SIF for the crack initiation as studied by A. Faye et al. [46] where the energy required for the crack initiation was shown to decrease with the crack tip radius. As noted, specimens from the design set 2 (symmetric hole location of radius 1mm) and design set 3 (hole with offset 1mm and 2mm radius) exhibit relatively higher values of K1D and K̇1.

Finite element model
Explicit finite element (FE) analysis was utilized to study the effect of the pre-existing holes on the stress-wave front and the resulting SIF. For that purpose, the 4 different experimental designs were modeled along with a pristine specimen, i.e. a specimen containing no drilled holes. In all cases presented here, a seam was introduced at the notch tip, having a length of 0.5mm to mimic the experimental pre-crack. The simulations were held using Abaqus/Explicit V6.14 [47]. A full 3D model of the experimental setup was constructed (specimen, incident bar, and striker) and all of its components were modeled using linear elasticity. To mimic the experimental boundary conditions as closely as possible, a steel support brick was modeled below the specimen. The dynamic material properties of PMMA [41] The contact between the striker-incident bar and incident bar-specimen was modeled as hard contact with no friction. The FE model was meshed using eight nodded brick elements (C3D8 elements in Abaqus) Figure 21. The mesh density was distributed in the specimen to maximize the accuracy at critical locations (i.e. the stress concentrator), without compromising the overall processing time, resulting in element edge length in the range of 0.2-2mm. Nonetheless, a mesh study was made for both element size and type, and the results reported in this section were found to remain essentially the same. Since the number, size, and location of the holes varied between the four specimen designs, the number of elements on the specimen was not constant and is summarized in Table 8.

Case
Number

Wave front pinning due to pre-existing holes
The vast majority of experimental data available in the literature regarding dynamically loaded structures with pre-existing flaws (e.g. [19,22]) are conducted in a manner in which the load is applied on the side of the specimen which is opposite to the notched/pre-cracked edge. Under such loading conditions, the propagating stress wave will have to interact with the pre-existing flaws before reaching the crack tip. Since the characteristic lengths of such specimens are not necessarily sufficiently large for the dynamic version of the Saint-Venant's principle to be invoked [48], the perturbation to the stress wave front, caused by the pre-existing flaw, will not self-correct by the time it reaches the crack tip, and as such may lead to local variation the established mode I field.
In figure 22, the calculated Von-Mises stress field, extracted from the FE calculations are presented for 35 ts  = (with 0 t = being the moment the stress wave entered the specimen) is presented for the five specimen geometries. When comparing Figure 22a (stress distribution in a pristine specimen) with the wave fronts calculated for the 4 experimental designs, it is evident that stress wave is both delayed and perturbed by the preexisting holes, an effect which seems to become more dominant with the increase in hole size and number. The perturbation and pinning of the stress wave front, affects not only the stress evolution at the crack tip but also the distribution of strain energy in the specimen. This, in turn, may play a role in determining the crack arrest time in the holes, as well as the re-emergence angle. in Figure 23 and the calculated SIF as a function of time is shown for the different cases in Figure 24.   From Figures 23 and 24 it is evident that the perturbation observed to result from the 1mm radius hole is rapidly diminishing and the differences between the SIF evolution for this case (Design set 2) and the pristine specimens are negligible in terms of crack growth initiation. Upon increasing the hole size to 2mm in radius, the mode I stress field is observed to grow at a slightly slower pace, resulting in a 14% decrease in 1 K at 46 sec t  = . Further in time ( 49 sec t  = ) this difference reduces to ~6% and shortly after the two specimens (pristine and Design 1) exhibit the same value of 1 K ( Figure 24). Surprisingly, despite the variance in experimental results between designs 2 and 3 (2mm radius with and without offset), the evolution of SIF seems to be identical, and the asymmetry, introduced by offsetting the hole at 0.5 of its radius away from the symmetry plane does not influence the rate of the build-up of the mode I stress field, the only observable effect is a slight distortion (broken symmetry) in the stress filed shape very close to the crack tip. As anticipated, the strongest perturbation (figure 23e) and as a result the strongest effect on 1 K is observed for design set 4 (3 holes with a radius of 2mm). The FE calculations demonstrate a ~5 sec  delay in the establishment of the mode I stress field (a nonzero value of 1 K ). As can be inferred from figure 24, the pinning of the stress wave front has led to a larger build-up of strain energy further away from the crack, leading to a higher loading rate once the mode I field is finally established.

Fracture simulations with a two-scale dynamic damage model
In [49,50] a damage evolution law was proposed for simulating dynamic crack propagation in brittle materials. The multiscale model in [49,50] is based on a Griffith type criterion for micro-cracks, homogenized to produce a rate-sensitive continuum damage model. Recently, the model was compared with experimental results available in the literature and was demonstrated yield satisfactory agreement. More specifically, loading rate effects were shown to be captured for several materials and specimen geometries [51], and a modified version of the model was shown to be in agreement with the experimentally measured thermal evolution at the vicinity of the crack tip [52].

The damage model
The multiscale approach assesses the macroscopic system variables such as stress, strain, and damage from the evolution of microcracks in the close vicinity of the crack tip (i.e. the process zone). The process zone is assumed to have a uniform periodic microcrack distribution with a spatial interval λ. The length of the microcrack is represented by l and the corresponding local damage is defined as, (1 ) ijkl ijkl The damage evolution law is deduced by homogenization [49,50] from the microscopic Griffith criterion for microcracks in the form: The angle bracket 〈⋅〉 represents the positive part of the expression.
Following [52], we assume that the fracture energy Ga depend on the microcrack tip speed v=(λ/2)dD/dt through the relation Ga = Ga0(1+αv), with the constant value Ga0 and α encapsulates the linear increment in Ga with the crack velocity.
In the simulations of the impact test, the system (4-9) is numerically solved for the elastic and damage fields and the fracture of the specimen will be the consequence of the damage field evolution.

Finite element implementation and model calibration
The continuum damage model described above was implemented in Abaqus/Explicit via VUMAT subroutine. The experimental setup, including the striker, incident bar, and specimens were modeled in 2D to avoid high computational cost. All other dimensions (including the pre-crack and hole diameter and locations) are kept as detailed in previous sections. Similarly, the loading conditions and contact definitions are identical to those described in Section 4. The striker, incident bar, and support brick, were meshed using quadrilateral elements (edge length varying from 2x10 -3 m to 4x10 -3 m). The fracture specimens were meshed using 2D linear triangular elements. The mesh density was increased in the central region of the specimens encapsulating the notch, and pre-drilled holes for accurately describing the crack propagation and crack-hole interactions. The elements' edge length in the specimens' meshes varied from 2x10 -3 m near the specimens' edges to 1x10 -5 m in the denser mesh region. The number of elements used for each specimen design is given in table 9.

Case
Number of elements Void radius = 1mm 2714197 Void radius = 2mm 1950181 Void radius = 2mm, offset = 1mm 2757142 Void radius = 2mm, 3 voids 1571555 Table 9 Number of elements used to discretize each of the different designs.
For remaining consistent with the experimental data analysis methodology described in Section 2, A Python script was utilized to extract the displacement and damage fields from the simulation. The data was then fed into the same Matlab script used for analyzing the experiments and the SIF, crack velocity, crack path, and crack arrest time were calculated and compared with the corresponding experiments.
The model calibration was carried out using the experimental results of design sets 1&4. Three parameters are required to calibrate the chosen continuum damage model, namely: Ga0, α, and λ. These parameters were determined by minimizing the errors between the experimental observations and numerical predictions. The range of parameters considered for the calibration was Ga0 = [20, 600] J/m 2 α = [0.01, 0.035] and λ = [3x10 -5 , 10 -3 ]m. The initial guess of the three parameters was taken to be: are Ga0 = 350J/m 2 , α = 0.025, and λ = 3.5x10 -4 m, following the parameters used in [51]. The final values, which were deemed to provide a reasonable estimation of the experimental data are Ga0 = 100J/m 2 , α = 0.028 and λ = 5x10 -5 m. Here we note that the value chosen for Ga0 is lower than the values usually found in the literature, which lie in the range of [200 -1000] J/m 2 for dynamic fracture experiments on PMMA. Indeed, in our simulations, the initiation toughness is underestimated, however, higher values of Ga0 did not yield satisfactory agreement for the propagation stage.
The following observations were made during the calibration process: • While larger values of Ga0 will result in a better agreement with the toughness, it reduces the crack propagation velocity and increases the crack arrest time inside the holes. • An increase in the value of Ga0 accompanied by a reduction of α also could not balance the increased resistance for crack propagation. In this scenario, λ comes into the picture. Due to the coupling between the effect of the parameters on the damage evolution rate (Equation 9) • Lower values of λ increase the material resistance to crack growth, resulting in the slowing down of crack and increase in the crack arrest time. Similarly, an increase in λ leads to faster crack propagation, and beyond a certain value, crack branching becomes a dominant mechanism.
In fact, during the calibration process, we observed that the damage rate, and hence the fracture toughness, crack propagation velocity and crack arrest, can be easily calibrated for a range of crack propagation velocities, however since our experimental data covers the range of 0-700m/s, one set of variables could not provide an accurate estimation over the entire range.
The preliminary analysis showed that a more complex, nonlinear dependency of the fracture energy on the crack propagation velocity may be necessary to obtain a shape similar to that reported by Scheibert et al. [53] as a better estimation of the entire range. However, this is not incorporated in the present analysis and is currently under investigation.

Damage model validation
As previously noted, the calibrated parameters and the choice of the linear approximation for the fracture energy were found to underestimate the crack growth initiation toughness in comparison with the experimental results on design set 4. Nonetheless, a satisfactory agreement was observed for the entire propagation regime as well as for the crack deflection post crack arrest. We thus chose to proceed with this set of parameters and examine its predictive capabilities by comparing the simulation predictions with the experimental results for design sets 1-3. The obtained predictions are compared with the experimental data in Figure 25 and Table 10. As can be seen in Figure 25, the calibrated damage model provides a reasonably good agreement with the experiments. In most specimens, the SIF did not deviate from the experimental values by a significant margin (given the experimental scatter). For design set 4 the SIF shoots up as the crack approaches the first void. As noted, we observed that the linear relation of Ga and v underestimates Ga at lower velocities and overestimates it at high velocities. As the crack approached the first hole in design set 4, it accelerates almost linearly, resulting in a sharp increase of the SIF in the simulation which is not apparent in the experiments. Larger values of λ were observed to remedy this effect but resulted in extensive branching which was not observed in the experiments. Crack growth velocity was found to be slightly higher than observed experimentally while following the same trends. The pronounced acceleration phase after the fracture initiation, which was not observed in the experiments, is also the consequence of the approximation considered for the fracture energy. The crack emerging angle and crack arrest time from the experiments and simulations are tabulated for comparison in Table 10, along with snapshots of the crack path. The crack paths in all the cases are replicated in the simulations remarkably well. The crack emerging angles measured from the experiments are comparable to the simulation results, with the largest difference observed for design 2 (hole radius of 1mm,) where the crack arrest time and the crack emerging angle are significantly smaller compared to the experimental values. We expected that the crack arrest time will be severely underestimated in the simulations, due to the value of Ga, however, it mostly falls in the range of the experimental results or close to it. The crack deflection from the symmetry plane, upon emerging from the hole, is mostly attributed to the tension-compression transition of the incoming stress waves, as illustrated in figure 17. If the crack initiates too early, the emerging angle would be lower and vice versa. The reduction in the resistance for crack initiation encourages early crack propagation and thus lowers the crack emerging angle. The horizontal crack path before reaching the void is unbent in experiments thanks to the sturdiness provided by the crack front curve and the 3D effects. Whereas these elements are absent in the 2D numerical analysis, hence having a slight deviation just before the crack arrest. This deviation is significant in determining the crack emerging and crack arrest time especially if the void radius is smaller, as highlighted in the design 2 when the crack approached the 1mm radius void with minor horizontal deviation.

Crack path
Although the model could not capture all the nuances in the crack propagation of a 3D experiment by conducting a 2D analysis, it has provided a reasonably good agreement with the experiments.

Summary
Dynamic fracture experiments were conducted on a series of holes containing RT-DCB PMMA specimens. The experimental and numerical results presented in this paper indicate a complex interaction between the stress waves in the specimen arising both from the pinning of the incoming wave due to the pre-existing flaws as well as from the interaction between the propagating crack and the holes in front of it. Variations in the crack growth initiation toughness, crack velocity and crack deflection angle upon existing the holes were attributed to the stress wave interaction and their effect on the accumulation of strain energy.
The effect of the pre-existing holes on the incoming stress wave, which was previously ignored in the literature dealing with similar experiments, provides new and important insights as to the design of crack arrestors for rapidly propagating cracks. Our results indicate, that by proper geometrical design, it is feasible to create crack arrestors for dynamically propagating cracks while controlling their emergence angle. Stress reversal (tension/compression) may serve as a key attribute in the design of such geometries.
The experimental data, collected via high-speed imaging followed by a customized data-reduction plan, was used to calibrate a two-scale continuum damage model, previously proposed by Dascalu et al. [49,50]. The implemented damage model was observed to yield reasonably accurate results, under the relatively simplified hypothesis of linear dependence of the fracture energy on the crack velocity, responsible for some discrepancies between the experiments and simulations. In future work, we intend to study different forms of this behavior, aiming to achieve correspondence between the experimental and simulated results across a wider range of loading conditions and crack growth velocities.