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 correla tion. 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 pr edictive 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 . 1. Introduction The study of brittle crack growth and arrest under dynamic loading condition s 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 mate rials 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 exi sting heterogeneities make s the dynamic fracture process highly sensitive to variations in the local microstructure [2–7]. Having an understanding o f the interaction of the growing crack with its surround ings 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–15]. Experimental studies on crack – flaw interaction for dynamically propagating cracks dates ba ck to the 1970s [16]. Kobayashi et. al. [16] have studied the role of holes as crack arrestors for dynamically propagating cracks, a nd 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–21] on crack -void in teractions 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 angl es, etc. Such experi mental 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–27], peridynamics [28,29] , and phase - field modeling [30]–[33], [34] . These numerical methods have different approach es to capture the sensiti vity of fracture in the local field near the crack tip, where the influence of material softening and the heat release may 2 be crucial [35,36] . Thus, to validate a model’s suitability for design purposes, it has to be tested against several scenarios, while calibrated onl y 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 Hopkinso n 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 2x106 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 evo lution, 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 w ave 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 gen eral agreement between the analysis and experimental data. F inally , the results and their implications are summarized in S ection 6. 2. Experimental Methodology and Data analysis 2.1 Specimen s’ 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 qua si- static 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 w as 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 int roduce 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 s harp 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 scattere d 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 sp ecimen. 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. a) b) c) d) Figure 2. 4 designs used in the experiments and dimensions (mm). Design parameters of all 4 sets are tabulated in Table 1. 4 Design set # Number of holes Radius of the hole (mm) Offset from the notch (mm) 1 (Fig. 2a) 1 2 0 2 (Fig. 2b) 1 1 0 3 (Fig. 2c) 1 2 1 4 (Fig. 2d) 3 2 0 Table 1. The four configuration s used for the study . 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 . Figure 3. Crack deflection angle definition . 2.2 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 st riker is launched at a predetermined veloc ity, 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 tri ggering a high -speed camera (Kirana, Speciali zed Imagin g) 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 ~ 611 10 secx− . a. b. Figure 4. a) An upper view of the experimental setup. b) A side view, demonstrating the positioning of the specimen . 2.3 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 d isplacement field was extracted using the digital image correlation open -source code NCOR R [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, a n 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 5 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 c rack tip is determined using the strain map in the y -direction as shown in Figure 5. a) b) Figure 5. Crack path extraction from the DIC. a) e yy from DIC. b) Obtained c rack path . 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. Consider the domain ℝ2 with the crack tip as the origin, the analytical displacement fields at a point (x, y) ϵ ℝ2 for a crack moving at a speed c are represented by ( , , )xu x y c and ( , , )yu x y c along the abscissa and the ordinate respectively . The transverse displacement ( , , )yu x y c is considered to evaluate the SIF due to the mode I fracture and the loading direction on the specimen and it is represented by a power series expression with N terms (Eq.1). 0 12( , , ) ( , ) ( , , )N n yy nnKuux y c x y F x y c ==+ (1) Where uy0(x,y) is the rigid displacement field , and Kn is the amplitude of each term in the power series and μ is the shear modulus of the material. The function ( , , )nF x y c represents the contour or the shape of the nth term in the expression in the 2D plane while Kn determines the amplitude. ( , , )nF x y c is a dependent function of the space and the instantaneous crack propagation velocity [43] /2 /2 12 1 1 2 22 ( )( , , ) ( 1) sin sin22nn nnn hnF x y c n B r r    = + − +          (2) Where the parameters are, ()2 2 1, tan , 1, 2i i i i r x y i yx  −= + = = () ()22 12 1 , 1lsc c c c = − = − Elastic and inertial properties of the material attribute to ( , , )nF x y c via cl and cd - the longitudinal and shear wave veloci ties of the material . 1, 1lskcc k +== − Where ρ is the density of the material and the factor k is determined based on whether the system is plane str ain or plane str ess. Plane stress analysis is relevant to this study since the displacement fields are observed from the material surface . 6 (3 ) (1 ) : Plane stress 3 4 : Plane strainvvkv−+=− The rest of the factors in Eq. 2 are, 2 1 2 2 2 22 (1 ) : odd() (1 ) / 2 : evennhn n   +=+ 2 2 22 1 2 21 4 (1 )B  +=−+ Using the displacement formulation ( , , )nF 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 field s 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) fie ld 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 become s unstable and sensit ive 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 d isplacement 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 ana lysis to reduce the processing time while having an accurate outcome . a) b) Figure 6. Effect of the number of terms in the displacement expression . a) SIF variation with the number of terms. b ) Normalized error distribution with the number of terms. The experimental and analytical displacement fields are given in Figure 7. The displacement i n 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 7 analytical displacement field and considering only the relevant frequencies for the comparison. So, th e process implicitly removes the random noise which is dominant in the DIC data. A set of point s 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 a ssessed (Eq. 2) using the least square met hod. 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). 0 12( , )N n y nee y y y nKu e u u u F x y =+= − = −  (3) a) b) Figure 7. U y displacement field near the crack tip. a) Experimental data from DIC, b) Formulated displacement field. The location corresponds to the least error is used to determine the crack tip location at a given frame (Figure 8b) and the corresponding SIF (Kn). a) b) Figure 8. Crack tip evaluation. a) Error distribution along the crack path in displac ement field mapping. b) Crack tip corresponding to the best overlap . 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 i f 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 exist s an uncertainty associated with the crack tip which leads to a minor vagueness in the evaluation of SIF. 8 3. 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 . 3.1 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 observ ed 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 y - direction (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 al ong 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 . a) b) c) Figure 9. Fractured specimens - design 1. a) b) Figure 10. a) SIF and b) Velocity of 3 cases during the fracture - design 1. The SIF and the crack propagation velocity during crack propagation are shown in Figure 10. The observed variations in the critical SIF, K 1D 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 velocit y are plotted for a 95% confidence interval. Error margin in the SIF is evaluated by the shortest distance to the closest neighbor. Specimen Impact velocity (m/s) K1D (MPa m1/2) 𝐊̇𝟏 (MPa/s m1/2) Time inside the hole (μs) Deflection angle (0) 1 19.9 1.14 1.53 x105 20.6±0.7 69±1.1 9 2 20.1 1.75 1.84 x105 25.4±0.8 72±1.1 3 19.7 1.40 2.03 x105 24.0±0.7 82±1.1 Table 2. Properties of fracture - design 1. Critical SIF, K 1D 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 a verage 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 740. 3.2 Design set 2 For design set 2, the hole’s radius was decreased from 2mm in design set 1 to 1mm . The f ractured 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 . a) b) c) Figure 11. Fractured specimens - design 2. a) b) Figure 12. a) SIF and b) Velocity of the crack during propagation - design 2 . 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 sa turated velocity is lower for the 1mm radius holes . When c ompared to design set 1, the K1D value is higher for all the cases which leads to higher initia l crack propagation velocity . Specimen Impact velocity (m/s) K1D (MPa m1/2) 𝐊̇𝟏 (MPa/s m1/2) Time inside the hole (μs) Deflection angle (0) 1 18.7 2.65 1.93 x105 5.3±0.7 53±2.3 2 19.4 1.74 1.42 x105 6.0±0.7 36±2.3 3 19.6 2.15 2.02 x105 10.8±0.7 34±2.3 10 Table 3. Properties of fracture - design 2. Once arrested, the crack spend s less time inside the hole when compared to design set 1 and deflects at a smaller angle. The a verage 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 410 (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 . 3.3 Design set 3 Design set 3 takes after design set 1 with an introduced a symmetry . The 2mm radius hole is offset by 1mm from the horizontal symmetry plane of the specimen ( Figure 13) . a) b) c) Figure 13. Fractured specimens - design 3. a) b) Figure 14. a) SIF and b) Velocity of the crack during propagation - design 3 . 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 th e symmetrical case is the lack of saturation in velocity when approaching the hole and the decrease in SIF which is more apparent here. Specimen Impact velocity (m/s) K1D (MPa m1/2) 𝐊̇𝟏 (MPa/s m1/2) Time inside the hole (μs) Deflection angle (0) 1 20.3 2.63 2.49 x105 11.4±0.7 59±1.1 2 20.0 2.51 2.28 x105 29.4±0.7 65±1.1 3 19.4 2.39 2.71 x105 24.6±0.7 72±1.1 Table 4. Properties of fracture - design 3. 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 hori zontally after exiting the hole and continued to grow at an angle. 11 3.4 Design set 4 This set of specimens consist of 3 symmetric holes with a 2mm radius spaced by a uniform gap of 10mm between the m. As shown in Figure 15, in 3 out of 4 specimens crack diverts after the 2nd hole. a) b) c) Figure 15. Fractured specimens - design 4. a) b) Figure 16. a) SIF and b) Velocity of the crack propagation to 1st & 2nd holes - design 4 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 th an 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 fie ld 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 1st hole. In specimen 2, crack initiates at an angle then regains the horizontal path, an d this causes the dip in halfway of the velocity variation. The i ncrease in the SIF after the 1st crack arrest is caused by the higher K 1D 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 1st hole. 1st crack arrest Specimen Impact velocity (m/s) K1D (MPa m1/2) 𝐊̇𝟏 (MPa/s m1/2) Time inside hole-1 (μs) Deflection angle (0) 1 20.6 0.98 1.49 x105 10.0 ± 0.7 7 ± 1.1 2 19.1 1.62 2.03 x105 13.3 ± 0.7 4 ± 1.1 3 20.0 1.58 2.14 x105 8.7 ± 0.7 -2 ± 1.1 Table 5. Properties of fracture - design 4. Crack to the 1st 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 1st crack arrest. Therefor e, the crack continue s to travel in the horizontal direction after being arrested and the measure d 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 12 crack, at a position where the magnitude of the tensile load is higher . The crack arrest time in the 2nd 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. 2nd crack arrest Specimen K1D (MPa m1/2) Time inside hole-2 (μs) Deflection angle (0) 1 4.87 35.3±0.7 29±1.1 2 4.59 30.6±0.7 59±1.1 3 3.51 28.6±0.7 68±1.1 Table 6. Properties of fracture - design 4 . Crack to the 2nd hole. 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 exist s a stage where the hole only experience s compressive loads for a few microseconds before tensio n is reestablished. In Figure 17a, the compression - tension 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 h ole and regenerate the stress singularity required for propagating the crack further . a) b) Figure 17. Compression -tension transition in the specimen. a ) σyy at t=100μsec. b) σyy variation with time at point P (20 mm ahead of the crack tip ). When the crack reaches the 2nd 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 2nd hole supports the compression -tension transition in the local stress field as seen in design set s 1-3. This process increases the crack arrest time in the 2nd hole compared to the 1st hole by 3 times. The c rack deflection angle is smaller compared to design 1 which might be due to the higher kinetic energy which is the result of the higher K ID. Crack propagation velocity stays the same after the 1st crack arrest. Meanwhile, the SIF keeps dro pping 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 2nd 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 1st crack arrest, hence not triggering any asymmetry in the specimen and keeping the crack trajectory horizontal after the 2nd arrest. 13 Figure 18. 3-hole specimen with straight crack after the second hole. a) b) Figure 19. a) SIF and b) Velocity of the crack in the 3-hole specimen where crack passes through all holes. The characteristics of SIF and crack propagation velocity until the 2nd crack arrest are similar to the rest of the cases in the set (Figure 19). The m agnitude of SIF and crack propagation velocity are very close from the 2nd hole to the 3rd hole of design 4 . K1D upon existing holes 1&2 are very close. Hole no. K1D (MPa m1/2) 𝐊̇𝟏 (MPa/s m1/2) Time inside the hole (μs) Deflection angle (0) 1 3.98 2.19 x105 8.7 ± 0.7 -0.5 ± 1.1 2 7.35 -- 32.0 ± 0.7 3.2 ± 1.1 3 7.16 -- 10.7 ± 0.7 13.5 ± 1.1 Table 7. Crack arrest time in each hole. As shown in Table 7, crack arrest time in hole -1 and 3 are very close and 1/3rd 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 K 1D, SIF , and crack propagation velocity which might be due to the identical geometry and the unaltere d local stress field. The c rack arrest time at the 2nd hole is comparable to the rest of the cases in this design set (Table 6). 4. 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 K 1D and K̇1 values are not, strictly speaking, colinear, they do exhibit a linear trend with log(K̇1) and K 1D (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 notch ed specimens without a pre -crack, the observed offset was expected. The s harpness 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 14 set 2 (symmetric hole location of radius 1mm) and design set 3 (hole with offset 1mm and 2mm radius) exhibit relatively higher values of K 1D and K̇1. Figure 20. Comparison of the experimental data with the literature survey. 4.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 spec imen 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 d ynamic material properties of PMMA [41] are taken as 5.76 , 0.42E GPa  == with E being Young's modulus and  Poisson's . The density was taken to be 31180 / Kg m= . The elastic properties of the Maraging steel incident and striker bars are given as 3190.5 , 0.325, 8100 /E GPa Kg m  = = = . 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 maxim ize 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 re sults 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 of elements No void 30608 Void radius = 1mm 38176 Void radius = 2mm 40482 Void radius = 2mm, offset = 1mm 40472 Void radius = 2mm, 3 voids 46144 Table 8. Number of C3D8 elements used to discretize each of the different designs. 15 Figure 21. Finite element model of the Specimen -Hopkinson bar setup 4.1 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 t he 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 l ead 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 35ts= (with 0t= 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 pre - existing 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 t he 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. a) b) c) d) e) Figure 22. Perturbed s tress wave front s demonstrated using the e Von mises stress field at the crack tip vicinity at t= 35 μsec. a) specimen without any hole; b) Design set 2 - 1mm radius hole; c) Design set 1 - 2mm radius hole, d) Design set 3 - 2mm radius hole located at an offset from the symmetry plane; e) Design set 4 -3 holes of radii 2mm. The evolution of the stress field in the local surrounding s of the crack tip with further wave propagation is shown for 35;43;46;49 sect  = in Figure 23 and the calculated SIF as a function of time is shown for the different cases in Figure 24 . Time after impact 35 μ sec 43 μ sec 46 μ sec 49 μ sec No hole Hole radius = 1mm 16 Hole radius = 2mm Hole radius = 2mm Offset = 1mm 3 holes, radius = 2mm Figure 23. Evolution of the Von -Mises stress field near the crack tip for the different specimens shown for 5 instances in time. Figure 24. The evolution of the mode I stress intensity factor as a function of time from impact. From Figures 23 and 24 it is evident that the perturbation observed to result from the 1mm radius hole is rapidly diminishing and the differ ences 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 pac e, resulting in a 14% decrease in 1K at 46 sect= . Further in time ( 49 sect= ) this difference reduces to ~6% and shortly after the two specimens (pristine and Design 1) exhibit the same value of 1K (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 a t 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 1K is observed for design set 4 (3 holes with a radius of 17 2mm). The FE calculations demonstrate a ~ 5 sec delay in the establishment of the mode I stress field (a non - zero value of 1K ). 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. 5. 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 materi als 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]. 5.1 The d amage 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, / Dl= . The damage variable is denoted [0,1]D , where 0D= represents the undamaged state, and D = 1 corresponds to a fully fractured state . All the mechanical fields are assumed to depend on the spatial variables 2 12( , )xx and time t. Following asymptotic homogenization [49], the macroscale field variables are defined. The macroscale equation of motion is 2 2ij i iu tx= (4) And the homogenized stress -strain constitutive equation as a function of the damage field reads: ( ) ( )ij ijkl klCD= u (5) Here Cijkl are the effective coefficients. They depend on the elastic coefficients of the virgin material through the Damage Mechanics linear approximation law: (1 )ijkl ijklC C D=− (6) The damage evolution law is deduced by homogenization [49,50] from the microscopic Griffith criterion for microcracks in the form: 21()( ) ( )2a R ijkl kl ijG C dD dC D dt dD=+ uu (7) The angle bracket 〈⋅〉 represents the positive part of the expression. A characteristic time parameter τc may be introduced cR C= , where CR is Rayleigh waves velocity , and the damage energy release rate is defined by ()( ) ( )ijkl kl ijdC DYdD =− uu such that upon further simplification of Equation (7) we get: 21a cG dD dt Y=− (8) 18 Following [52], we assume that the fracture energy Ga depend on the microcrack tip speed v=(λ/2)dD/dt through the relation Ga = G a0(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 . 5.2 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 us ing quadrilateral elements (edge length varying from 2x10-3m to 4x10-3m). 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-3m near the specimens’ edges to 1x10-5m 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 an d compared with the corresponding experiments . The model calibration was carried out using the experimental results of design set s 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/m2 α = [0.01, 0.035 ] and λ = [3x10-5, 10-3]m. The initial guess of the three parameters was taken to be: are G a0 = 350J/m2, α = 0.02 5, 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/m2, α = 0.028 and λ = 5x10-5m. Here we note that th e value chosen for Ga0 is lower th an the values usually found in the literature, which lie in the range of [200 -1000] J/m2 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 valu es 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 dom inant 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 . 19 The p reliminary analysis showed that a more complex, nonlinear dependency of the fracture energy on the crack propaga tion 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. 5.3 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 ent ire 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. Design SIF Crack velocity 1 2 3 4 Figure 25. SIF and crack velocity comparison of experiments and simulation As can be seen in Figure 2 5, 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 20 (given the experimental scatter) . For design set 4 the SIF shoots up as the crack approaches the first void. As noted, we observed that t he linear relation of Ga and v underestimates Ga at lower velocit ies and overestimates it at high velocit ies. As the crack approached the first hole in design set 4, it accelerat es 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. Crack path Deflection angle Crack a rrest time (μsec) Experiment Simulation Experiment Simulation Experiment Simulation 690, 720, 820 59.80 20.6, 25.4, 24.0 19 340, 360, 530 29.40 5.3, 6.0, 10.8 3 590, 650, 720 61.70 11.4, 29.4, 24.6 24 290, 590, 680 50.30 35, 30.6, 28.6 18 Table 10. Crack path crack emerging angle , and crack arrest time of experiments and simulations. 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 pat hs 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 small er compared to the experiment al values . We expected that the crack arrest time will be severely underestimated in the simulations, due to the value of Ga, however , it mostly f alls 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 t he 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 21 is smaller , as highlighted in the design 2 when the crack approached the 1mm radius void with minor horizontal deviation . 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. 6. Summary Dynamic fracture experiments were conducted on a series of hole s 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 exist ing 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 exper iments, 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 contr olling 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 reasonabl y 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 cra ck growth velocities. 7. Acknowledgments The financial support provided by the European Union's Horizon2020 Program (Excellent Science,Marie -Sklodowska -Curie Actions) under REA grant agreement 675602 (Project OUTCOME), is gratefully acknowledged. The aut hors would like to thank Prof. C. Dascalu .