Experimental investigations of dynamic crack propagation in PMMA induced by impact show single cracks running at around 300-400 m/s. Existing numerical models for simulating dynamic fracture in PMMA consistently produce crack propagation speeds significantly higher than those measured experimentally. Here we uncover the reason for this puzzle by showing that localized softening in the fracture process zone (caused by heating due to high strain rates in front of the crack tip), leads to crack propagation speeds that match the observed ones. We introduce a new constitutive model in our peridynamic formulation for PMMA to account for material softening in the crack tip region. With the new model, the computed crack speed and crack length evolution match very closely those found experimentally. Keywords: Dynamic fracture, PMMA, peridynamics, crack propagation speed, nonlinear elasticity 1. Introduction PMMA (polymethyl methacrylate) is a thermoplastic polymer used as an alternative to glass because of its transparency and low density. PMMA’s dynamic fracture behavior make it widely used in many industries: automotive, aerospace, medical, etc. Dynamic fracture of monolithic and bi-layered PMMA samples using Digital Gradient Sensing (DGS) technique was recently studied in [ 1]. Crack propagation speeds in the monolithic sample were found to be in the 300-400 m/s range. Computational modeling of dynamic crack growth in PMMA has consistently shown much higher crack propagation speeds than those measured experimentally. The comparative simulation study in [ 2] showed faster crack speeds found by XFEM and element deletion technique than those measured. It was noted that when making the PMMA fracture energy five-six times higher than the actual one, the computed crack speed is reduced. Dynamic crack propagation in PMMA was also studied in [ 3] where the crack speed obtained using several computational methods (XFEM, peridynamics, cohesive zone model) were compared. Again, all of the methods utilized in [ 3] produced significantly higher crack propagation speeds than in experiments. In another recent paper [ 4], the dynamic crack speed was computed by an approach called bond-particle method (BPM). The BPM also leads to a computed crack speed from dynamic fracture in PMMA that is much higher speed than observed in experiments. 2 In this paper, we try to understand the reasons behind the failure of computational models in capturing the measured crack speed for dynamic fracture in PMMA. We first confirm that a peridynamic model utilizing a linear-elastic with brittle-failure leads to the faster-than-measured crack propagation speed in the impact and dynamic fracture in a monolithic PMMA sample tested experimentally in [ 1]. While the results, in terms of crack path and propagation speed converge in the limit of the horizon going to zero, the crack pattern is different from the one observed in experiments, and the discrepancy in the crack propagation speed is a factor of two. To remove some possible causes for this discrepancy and better understand the dynamic fracture in PMMA, we vary the impact loading for a possibly better match with the actual one. While this effort leads to a single straight crack path, as observed experimentally, it does not resolve the crack velocity problem as the computed crack still runs at almost twice the speed measured. We hypothesize that PMMA’s behavior at high strain rates (like those involved in the impact experiments in [ 1]) is different from the simple model used. Indeed, although similar to glass in some respects, PMMA has a very different microstructure compared to glass, featuring long chemical chains, cross-linked. Before fracture occurs, deformations of the polymer network can cause a measurable temperature rise. The localized heating ahead of the crack tip is caused by the fast straining of macromolecules in the fracture process zone [ 5]. PMMA’s behavior under high strain rate impact by using a split Hopkinson pressure bar was studied in [ 6]. The SEM photomicrographs of fracture surfaces showed some local melting. For this to occur, the local temperature had to reach higher than the softening temperature of PMMA (see Fig. 1). The estimated maximum temperature generated at the crack tip was around 530 [ 6]. The temperature rise near the crack tip in PMMA was studied using ℃ temperature sensitive liquid crystal film and infrared detector in [ 7]. It was indicated that the rise of temperatures (around 500 K) at the crack tip in PMMA did not change much as the crack tip velocity varied between 200 m/s and 600 m/s (see Fig. 2). Fig. 1 SEM photomicrograph of PMMA fracture surface suggesting localized melting (from [ 6]). 3 Color T High Low Fig. 2 Temperature rise along the crack path (left to right) measured by liquid crystal film (blue=higher temperature) (from [ 7]). Computational models utilizing local softening in simulating dynamic fracture in PMMA do not exist. Based on these observations, we introduce a new peridynamic constitutive model that accounts for softening in the fracture process zone. This approach is meant to reflect the actual softening produced by the local heating that is triggered under the high strain rates present in the fracture process zone. We test the new model in terms of crack propagation speed to determine the appropriate size of the softened region. The peridynamic computations reveal a shift in the time scale from the moment of impact until the first data point reported in the experiment. When this time-shift is implemented, we find that the new model is able to explain the observed behavior. This peridynamic model is the first to validate experimental results on dynamic fracture in PMMA. The paper is organized as follows: section 2 presents a brief review of peridynamics; in section 3 we describe briefly the PMMA impact experiment from [ 1], and experimental results that will be used for comparison with those from the peridynamic model; in section 4 we discuss the loading conditions in detail with the help of a full 3D LS-Dyna simulation and show the discrepancies between the regular linear-elastic-brittle PD model and the experimental results; in section 5 we introduce the new constitutive model that accounts for softening (induced by local heating) in the fracture process-zone and determine the appropriate size for the softened elastic zone; we present conclusions in section 6. 2. A brief review of peridynamics Peridynamics was introduced as a nonlocal form of continuum mechanics by Silling in 2000 [ 8] for modeling damage and fracture. Since then, it has been extended to a variety of other problems in which domain changes/discontinuities are part of the problem [ 9-14]. In this theory, each material point is connected through peridynamic bonds to other points within a certain neighborhood region called “the horizon region” [ 15]. The peridynamic bonds transfer forces/heat/mass between points and their failure is used to defines damage at a point. In peridynamics, one replaces the equation of motion by an integro-differential equation in which spatial derivatives are eliminated. This allows it to avoid the mathematical difficulties and inconsistencies present in the classical theory when cracks, for example, develop in the domain. The PD equations for dynamic problems are: 4𝜌𝒖(𝒙,𝑡)=∫ 𝐻𝒇(𝒖 (𝒙,𝑡)‒𝒖 (𝒙,𝑡),𝒙‒𝒙)𝑑𝐴𝒙+𝒃(𝒙,𝑡)(1) where is the acceleration vector field, u is the displacement vector field as a function of point and 𝒖 𝒙 time , b is the body force density, and ρ is the mass density. Also, f is the peridynamic pairwise force 𝑡 density function that describes the interaction between material points. The horizon region is the internal sub-region H (see Fig. 3), defined as: 𝐻={𝒙𝜖ℜ:||𝒙‒𝒙||<𝛿} (2) Fig. 3. The deformation of a peridynamic bond Let be the relative displacement and be the relative position in the 𝜼=𝒖 (𝒙,𝑡)‒𝒖 (𝒙,𝑡) 𝝃=𝒙‒𝒙 reference configuration between two material points of and . From Eq. 2 we have 𝒙 𝒙 ||𝝃|| >𝛿⇒𝒇(𝜼,𝜻) . When the pairwise force derives from a micro-elastic potential , a micro-elastic material is =𝟎 𝑤 defined by the following interaction force: 𝒇(𝜼,𝝃)=∂𝑤(𝜼,𝝃) ∂𝜼(3) A linear micro-elastic material is obtained if we consider: 𝑤(𝜼,𝝃)=𝑐(𝝃)𝑠2||𝝃|| 2(4) where is called the bond micromodulus function and s is the relative elongation of the bond, or the 𝑐(𝝃) bond strain: 𝑠=‖𝜼+𝝃‖‒‖𝝃‖ ‖𝝃‖(5) The pairwise force derived from Eq. 3 and Eq. 4 is: 𝒇(𝜼,𝝃)={𝜼+𝝃 ‖𝜼+𝝃‖𝑐(𝝃)𝑠, |𝝃|≤𝛿 0, |𝝃|>𝛿(6) In this study, we use plane stress conditions, with the conical micromodulus choice (see[ 16]): 𝑐(𝝃)=24𝐸 𝜋𝛿3(1‒𝜈)(1‒||𝝃|| 𝛿)(7) where E is the elastic Young’s modulus, and =1/3 is the Poisson ratio in 2D plane stress conditions. 𝜈 5 In peridynamics, damage is modeled using the critical bond strain concept, allowing a bond to break and no longer sustain a force [ 8, 17] once its strain goes beyond a critical value. In this study, once a peridynamic bond breaks, it remains broken [ 9, 10]. The critical bond strain parameter is obtained from the measured fracture energy, [ 17]. In 2D plane stress, the connection between and is [ 16]: 𝐺0 𝑠0 𝐺0 𝐺0=2∫𝛿 0∫𝛿 𝑧∫𝑐𝑜𝑠‒1(𝑧𝝃) 0[𝑐(𝝃)𝑠02(𝝃)/2]𝝃𝑑𝜃𝑑 𝝃𝑑𝑧(8) For the conical micromodulus given in Eq. 7, one obtains the critical strain as (see [ 18]): 𝑠0=5𝜋𝐺0 9𝐸𝛿(9) The region defining the material is discretized with the one-point Gaussian quadrature, using uniform grid with spacing Δx (nodes are centered in the middle of each cell). The discretized (in space) version of Eq. 1, written at a time-step is: 𝑡𝑛 𝜌𝒖𝑛=∑ 𝑗∈𝐻𝑖𝒇(𝒙𝑗,𝒙𝑖,𝑡𝑛)𝑉𝑗+𝒃𝑖,𝑛(10) in which , and are node numbers, is the volume of node , and n is the time step number. For the 𝑖 𝑗 𝑉𝑗 𝑗 temporal discretization we use the velocity-Verlet explicit scheme [ 19, 20]: 𝒖 𝑛+1 2=𝒖𝑛+∆𝑡 2𝒖𝑛(11a) 𝒖𝑛+1=𝒖𝑛+∆𝑡𝒖 𝑛+1 2(11b) 𝒖𝑛+1=𝒖 𝑛+1 2+∆𝑡 2𝒖𝑛+1(11c) where , and, denote displacement, velocity, and acceleration vectors, respectively. 𝒖𝒖 𝒖 Given that the uniform square discretization grid of the domain does not conform to the disk-shaped horizon region, several algorithms have been introduced to reduce the quadrature error caused by this [21]. The simple HHB algorithm (named so in [ 22] but introduced in [ 21]), provides monotonous convergence to the classical result, when the horizon factor increases. Here we use this algorithm, which corrects the bond length for nodes that are only partially covered by the horizon region. 𝝃 We define the nodal damage index as the ratio of the number of broken bonds and the total number of bonds at that node. The value for the damage index varies between zero (intact material) and one (fully damaged). As we shall see below, in the new model it will be necessary to determine the position of the crack tip as the crack grows. We use a general algorithm that detects the crack tip based on the strain energy density value (see section 4.4). We use a PD code implemented to run on a single CPU with an NVIDIA Tesla P100 GPU [ 23]. The PD simulations using the smallest horizon (and densest grids) over 10k time steps take about 10 minutes to complete on an Intel Xeon E5-2670 2.60-GHz processor with the Tesla P100 GPU. 63. Review of Experimental Test Results The dynamic fracture of monolithic and bi-layered PMMA samples using Digital Gradient Sensing (DGS) technique was studied in [ 1]. Here we only focus on the monolithic case. Material properties for PMMA are shown in Table I. A rod was used as the striker, hitting a long bar which then transferred the loading into the PMMA via two inclined faces of a V-notch. The notch, at its end, has a 2mm-long thinly cut pre-crack (by sharp razor blade). The crack path was recorded by 32 images which allowed for an estimation of the average crack propagation speeds. The dimensions of the PMMA specimen are illustrated in Fig. 4.a. At the impact speeds used in [ 1], a single crack propagates through the monolithic sample in a close to straight-line fashion. Deviations from a straight-line crack trajectory (see Fig. 4.b) could be due to small geometrical imperfections (leading to slight variations in wave reflections and interactions with the growing crack) or slightly asymmetric contact loading at the moment of impact. Histories of the crack tip velocity and the crack length are shown in Fig. 5. (a ) (b) Fig. 4 (a) Geometry of the monolithic PMMA specimen used in [ 1]. Impact happens along the notch surfaces , and (b):Photograph of the fractured specimen showing crack path in a monolithic PMMA sample (from [ 1]). Arrowhead indicates crack growth direction. Table I Material properties for PMMA (from [ 1]). PMMA Young’s Modulus, GPa 5.0 Poisson’s ratio 0.34 Fracture Toughness, MPa.m0.51.31 Density, kg/m31010 thickness, mm 8.6 7 Fig. 5 Crack growth history from experiments (re-drawn from [ 1]) 4. Peridynamic modeling of the dynamic impact 4.1. Loading conditions Since the pressure pulse generated by the impact event was not reported in the experiments in [ 1], we develop a 3-D finite element model in LS-DYNA to estimate this pressure pulse generated on the PMMA V-notch (see Fig. 6) before the fracture event. Although fracture in LS-DYNA can be simulated with the element deletion technique, for example, this approach is less accurate when the model geometry and inertia are important. Therefore, we will not model the fracture event in LS-DYNA. The PMMA has a pre-crack, and we define it by removing a tiny area with the length of 2 mm from the PMMA cross-section, before meshing. Calculated in this way, the pressure pulse will only be valid before fracture initiates in the sample [ 19]. Since the fracture initiation was not clearly defined in the experiment, we use the entire pressure pulse for our first PD simulation. Once we determine the fracture initiation time, , (measured from the time of impact, the start of our PD simulation) from this 𝑡𝑓 PD simulation, we use the computed pressure pulse only up to , and modify it after (see section 𝑡𝑓 𝑡𝑓 4.2). The finite element (FE) model shown in Fig. 6 includes three parts, modeled as linear elastic materials: a striker, a long bar, and a monolithic PMMA, with the same dimension as in [ 1]. A single-surface contact is defined (coefficient of friction of 0.35) between the bar and the PMMA notch surfaces [ 24]. We apply the initial velocity of 15 m/s, used in the experiment, on the striker part and we solve the FE model with the dynamic explicit solver [ 25]. 8 Striker Long bar PMMA Long barPMMA Fig. 6 (a) Geometrical setup in LS-DYNA FE model, (b) the interaction between the impacting long bar and the PMMA sample with the V-notch. The resulting stress wave and its profile along the long bar are shown in Fig. 7.a and Fig. 7.b. The pressure pulse rises within 15 μs, remains constant for about 100 μs, and drops to zero after that. The average peak pressure is 106 MPa and matches the analytical calculation from , where is 𝜎=0.5𝜌𝑐𝑣 𝑐 the longitudinal wave speed, is the material density, and is the part velocity [ 26]. The compressive 𝜌 𝑣 pressure on the V-notch is shown in Fig. 7.c. The pressure rises to 60 MPa within 60 μs and drops after that. Since no fracture is modeled in LS-DYNA FE model, this load history is valid up to failure initiation in PMMA, which happens at . 𝑡𝑓 Fig. 7 (a) Compressive wave in the long bar, (b) the propagation of the compressive stress wave in the long bar, (c) pressure measured on the V-notch surface. We consider a linearized pressure profile as shown in Fig. 7.c, then we convert it to the body force components. We apply this body force on the PD nodes located on parts of the V-notch surfaces (see (b)(a) (c)(a)(b) 9Fig. 8 and [ 27]). It is important to notice that the impactor does not reach to the tip of the notch. We solve the 2-D bond-based PD model with free conditions on the remaining parts of the boundary. 4.2. Convergence study for dynamic fracture PMMA We perform a δ-convergence study (horizon size decreases while keeping the number of nodes covered by a node’s horizon roughly the same, see [ 28]) on the monolithic PMMA to see if the fracture patterns and crack speed converge. We use horizon sizes of 3.2, 1.6, and 0.8 mm with the horizon factor (the ratio of horizon size over the discretization size) of four. Fig. 8 The impact pressure from the bar is applied on PD nodes on parts of the V-notch surfaces as a body force. Using the impact load shown in Fig. 7.c, we calculate histories of crack tip velocities and damage maps with different horizon sizes shown in Fig. 9. Note that the fracture initiates at around =20 μs from the 𝑡𝑓 time of impact (the start of the PD simulation) and the crack tip velocities found by the peridynamic model are approximately twice as fast as those measured experimentally. Interestingly, the crack in the monolithic PMMA branches, while the one in the experiment did not. In order to explain these discrepancies between the experimental observations and the results computed with the PD model, we offer two possible explanations: (1) Since the PMMA fractures at (in the PD simulation), the pressure 𝑡𝑓 load (obtained by the FEM analysis) used after that is not valid and needs to be modified; and (2) dynamic fracture in PMMA could be influenced by factors like temperature rise at the crack tip that triggers softening in that region. 100 50 100 150 Time [s]0100200300400500600700800900Total Crack Tip Velocity [mm/ms]H3.2mm H1.6mm H0.8mm(a) (b) δ=3.2mm δ=1.6mm δ=0.8mm Fig. 9 Convergence study for the peridynamic model for the monolithic PMMA using horizon sizes of δ=3.2, 1.6, and 0.8 mm, respectively. (a) History of crack velocity, (b) deformed configuration at the end of the simulation. The colors show the damage distribution, and the small white disks show the relative size of each horizon region used. Using the pressure profile obtained from the LS-DYNA simulations before , we test different possible 𝑡𝑓 load profiles (see Fig. 10). Keeping the load relatively constant for a duration of about 60 μs after crack initiation seems a good choice (see the behavior reported in a similar experiment in [ 29]). We investigate two possible unloading paths (see Fig. 10), since wave interactions can affect the crack propagation speed [ 30]: Case I: Applied pressure gradually drops to zero; Case II: Applied pressure suddenly drops to zero. With these loading conditions, and a horizon size of 0.8 mm, the crack velocity ( Fig. 10) oscillates before reaching values twice those measured experimentally in the monolithic sample (see Fig. 5). Note that, with the input loading shown in Fig. 10 for Case I, the crack path no longer branches and a single crack, fully splitting the sample, is obtained, as observed experimentally (see Video 1). For Case II loading, the crack arrests at 140 μs, and the sample is not fully split. The significantly faster computed crack propagation speed issue, however, is not resolved by only modifying the loading. Therefore, we test the second hypothesis, that the possible softening of the PMMA around the crack tip, caused by the local heating induced by high strain rates in the crack tip region (see section 1), is responsible for the discrepancy between the computed and measured crack speeds. 110 50 100 150 200 250 300 350 400 Time, [s]02468101214161820Applied Pressure (MPa)Case I Case II(a) (b) (c) Fig. 10 Investigating load conditions with δ=0.8mm. (a) two loading scenarios, (b) crack tip velocities for two loading cases, and (c) damage map and the crack trajectory for case I loading at 300 μm after impact (see also Video 1). 5. A new PD model accounting for softening due to heating near the crack tip The other possible reason for which the PD results show crack speeds double those measured, is the presence of important temperature changes at the tip of the propagating crack and the subsequent softening [ 7]. Temperature rise observed in dynamic fracture experiments on PMMA has been shown to core a thickness of around 3.0 μm [ 7, 31]. In this section, we test whether a PD model that accounts for thermal softening near the tip leads to crack propagation speeds that match the observations. We call the area in which we implement softening (by reducing the micromodulus values of bonds in this region) the “Heat Affected Zone” (HAZ). 5.1. Locating the fracture process zone To apply the HAZ at the crack tip, we need to first determine the crack tip location ( ) at each 𝒙tip solution time-step (k). We use the strain energy density values ( ) at each node calculated as follows: 𝑤𝑛 𝑤𝑛=∫ 𝐻𝑤(𝝃,𝜼) (12) We calculate the value of using for unbroken bonds from the previous solution step (k-1) 𝑤𝑛 𝑤(𝝃,𝜼) to locate , defined here as the node with the highest value (see Fig. 11). After locating , we 𝒙tip 𝑤𝑛 𝒙tip check if a bond is entirely inside the HAZ (if both end-nodes are inside the HAZ). For such bonds we 𝝃 use a smaller micromodulus value, when their strain is a tenth of the critical bond-strain, resulting in a bi-linear force-strain elastic response. For all other bonds, the original linear model is used (see in Eq. 6). Assuming, for simplicity, a circular shape for the HAZ with a radius of (see Fig. 12), we define a 𝑟HAZ parameter that decides whether a bond is inside the HAZ or not, based on the distances and 𝜆HAZ 𝑑𝒙 𝑑𝒙 of nodes and to node as follows: 𝒙 𝒙 𝒙tip 𝜆HAZ(𝒙 , 𝒙)={1, 𝑑𝒙≤𝑟HAZ and 𝑑𝒙≤𝑟HAZ 0, otherwise(13) 12 (a) (b)(c) Fig. 11 Strain energy density (a) and damage index (b) in PMMA near the pre-crack tip. In (c), the superposed damage index and strain energy density maps, showing the location of the highest strain energy density in front of the crack tip. 𝑟𝐻𝐴𝑍 =2𝛿 Crack Tip 𝒙 𝒙𝝃𝑑𝑥 𝑑𝑥 Fig. 12 Finding if a bond is inside the Heat Affected Zone. Here a radius of 2 horizon sizes is used. 5.2. The new peridynamic constitutive model As mentioned in the introduction, the temperature at the tip of a running crack in PMMA can be as high as 530 °C [ 6]. Experimental measurements ([ 24, 32]) report that Young’s modulus drops from 5.02 GPa at room temperature to 0.56 GPa at 659 K (see Table II). The fracture energy of PMMA appears to be independent of the temperature change (see [ 33]). Therefore, we assume fracture energy that does not depend on temperature directly. However, since the critical strain depends on Young’s modulus value (Eq. 8), and the modulus is significantly smaller in the 𝑠0 HAZ, depends on location: the value in the HAZ is different from that in the rest of the sample. The new 𝑠0 constitutive model, that takes into account the presence of the HAZ is set as follows (see Fig. 13), for material inside the horizon of a given point: 𝒇(𝜼,𝝃)={𝜼+𝝃 ‖𝜼+𝝃‖𝑐(𝝃)[0.1s +0.09 𝑠0], if s >0.1𝑠0 and 𝜆𝐻𝐴𝑍 =1 𝜼+𝝃 ‖𝜼+𝝃‖𝑐(𝝃)𝑠, otherwise (14) 13 Note that when . We first assume the PD micromodulus is unchanged when the 𝒇(𝜼,𝝃)=𝟎 |𝝃|>𝛿 bond strain is smaller than (see Fig. 13). However, if the bond is in the HAZ, as the overstress 0.1𝑠0 macromolecules cause local heating, we decrease the micromodulus to , according to the measured 0.1𝑐 behavior in [ 5]. As observed experimentally, the crack runs fast from the beginning of its propagation. Therefore, we activate the HAZ from the onset of the crack growth and, keep it until the end of the simulation (see Fig. 5.a). 0 4 Bond Strain01Bond Force 0.1s0s0sbi-0BilinearLinear Fig. 13. Linear and nonlinear (bi-linear) bond force-strain relation for bonds outside and inside of the HAZ, respectively. To maintain the fracture energy constant, independet whether the bonds are in the HAZ or outside 𝐺0 of it, the following condition needs to be satisfied (areas under the bond force-strain relationships need to match): 1 2𝑐𝑠02=1 2𝑐(0.1𝑠0)2+𝑐(0.1𝑠0)(𝑠bi‒0‒0.1𝑠0)+1 2(0.1𝑐) (𝑠bi‒0‒0.1𝑠0)2(15) We find the critical strain for the bilinear function needs to equal 2.40 . 𝑠bi_0 𝑠0 Table II Young’s modulus of PMMA at different temperatures (from [ 24, 32]). T (K) 200 250 300 350 400 450 500 550 600 650 E (GPa) 5.02 4.29 3.81 3.52 3.11 2.71 1.25 1.25 0.79 0.56 5.3. Determining the extent of the Heat-Affected-Zone We now test the effect the HAZ extent has on our PD simulation results. Since the measured HAZ is much smaller than the horizon size used in this study, we calibrate the HAZ size by comparing the crack velocity with experimental data. One could, of course, use a model in which the horizon size is smaller than the measured heating zone, but that would lead to a much larger computational cost. Ways to reduce the cost of such computations could involve the coupling PD with FEM (see, e.g. [ 34]), but we do not pursue this here. 𝑥𝑡𝑖𝑝 14 Because the starting time for the experimental data on crack length and velocity is not clearly identified in [ 1], we determine its relation to our computational starting time (the moment of impact) by placing the experimental data crack length values on the same graph with our computed values. We notice that a shift of about 90μs in the time axis is required and we use the same time-shift for the velocity data (see Fig. 14). Recall that with the original PD model (see section 2), the crack tip reached velocities about twice the values measured in experiments (see Fig. 10.b). With the new constitutive model, considering the HAZ near the crack tip, the crack tip velocity and the evolution of the crack length are shown in Fig. 14, together with the data from the experiment and the damage map at 400 μs after impact. Using a horizon size of 0.8 mm, three sizes for the radius of the HAZ are investigated with the case I loading (see Fig. 10.a): δ, 2δ, and 3δ. As shown in Fig. 14.a, the crack velocity is in the same range of values as the experimental ones only when the HAZ radius is larger than δ. Increasing the size of the HAZ from 2δ to 3δ does not change the crack speed by much. We can conclude that a HAZ radius of 2δ is a good choice for considering the effect of local heating in PMMA. Video 2 shows the dynamics of the crack growth. Notice that the computed values for the crack propagation speed appear to be more oscillatory than those provided in [ 1], but, on average, the PD results match closely the measured speed values and this is the reason why the evolution of the crack length matches very well with that provided in [ 1] (see Fig. 14). Crack speed fluctuations in the dynamic fracture of PMMA are not surprising: for example, in [ 35- 37], dynamic instabilities are mentioned for crack propagation in PMMA when cracks reach speeds of 330 ±20 m/s. Some of the oscillations seen in Fig. 14 might also be induced by the differences between the real loading profile and the one we used. Sensitivity of crack speed to loading is well-known ([ 30]). For example, we observe (see Appendix A) that by increasing the loading magnitude to reach a peak, at the same =20 μs, that is 30% higher than case I loading, we find that the crack branches (the new PD 𝑡𝑓 model with HAZ is used) and the speed profile changes too (see Fig. 15). 0 50 100 150 200 250 300 Time [s]01020304050607080Crack Length [mm]Experiment PD, no HAZ PD, HAZ=2 0 50 100 150 200 250 Time [s]0100200300400500600700Crack Tip Velocity [m/s]Experiment HAZ= HAZ=2 HAZ=3 (a) (b) (c) Fig. 14 (a) History of the crack velocity of monolithic PMMA with two different HAZ, (b) Comparison of crack Length in monolithic PMMA between experiment and PD simulation with and without HAZ=2δ, (c) Damage map for HAZ=2δ, at 400 μs after impact (see also Video 2). Experimental data taken from [ 1]. All results obtained using δ=0.8 mm. 156. Conclusions In this paper, we found an explanation for a long-unanswered puzzle: the computed crack propagation speed in dynamic fracture of PMMA is significantly higher than that measured in experiments. We introduced a new peridynamic model and validated it for dynamic fracture in PMMA. It was observed that by using a regular linear-elastic-brittle-failure model for simulating dynamic fracture in PMMA, the computed crack speed is almost twice as fast as that measured experimentally. To explain the discrepancy, we hypothesized that material softening near the crack tip might be influencing the fracture behavior. Based on existing experimental evidence for temperature rise around the crack tip in dynamic fracture in PMMA, we considered a model in which the effective Young’s modulus is decreased significantly in the fracture process zone. We considered a Heat Affected Zone (HAZ) around the crack tip in which a nonlinear (bilinear) elastic model is used before peridynamic bonds break. The model with softening around the crack tip reproduced the experimentally measured crack speed and evolution of the crack length with high fidelity. The new peridynamic model allowed us to understand the relevant mechanisms in dynamic fracture of PMMA. With this model, we demonstrated the importance of the softening behavior in front of the crack tip in dynamic fracture of PMMA materials. The model will be used in the future to study the influence of strong and weak interfaces in bi-layered PMMA on crack behavior. 7. Acknowledgments This work has been supported by AFOSR MURI Center for Materials Failure Prediction through Grant Number FA9550-14-1-0073 (program managers Drs. Jaimie Tiley, Ali Sayir, David Stargel, and Fariba Fahroo), and ONR Award #N00014-16-1-2173 (program manager William Nickerson). This work was completed utilizing the Holland Computing Center of the University of Nebraska, which receives support from the Nebraska Research Initiative. Appendix A In this appendix, we investigate the effect of load magnitude on the crack speed and fracture pattern/mechanism. A recent experiment on PMMA showed how sensitive the crack path is to the loading condition [ 38]. We increase by 30% the amplitude of Case I loading (see Fig. 10) and use PD with the HAZ (extent of the HAZ is 2δ) and the new constitutive model. The new impact load leads to a changed fracture pattern and a different crack speed profile (see Fig. 15). Interestingly, in this case, the crack branches early, whereas it did not branch with the lower impact loading. 160 50 100 150 200 250 Time [s]0100200300400500600700Crack Tip Velocity [m/s] 05101520253035404550 Crack Length [mm]Crack tip velocity Crack length (a) (b) Fig. 15 A higher impact magnitude leads to crack branching near the pre-crack tip in the new PD model for PMMA that accounts for the softening due to heating. (a) Histories of crack tip velocity and crack length, (b) Damage map. Results computed using δ=0.8 mm.