Glass Fracture Upon Ballistic Impact: New Insights From Peridynamics Simulations

: Most glasses are often exposed to impact loading during their service life, which may lead to the failure of the structure. While in situ experimental studies on impact-induced damage are challenging due to the short timescales involved, continuum-based computational studies are complicated by the discontinuity in the displacement field arising from the propagation of cracks. Here, using peridynamics simulations, we investigate the role of the mechanical properties and geometry in determining the overall damage on a glass plate subjected to ballistic impact. In particular, we analyze the role of bullet velocity, bullet material, and elastic modulus, fracture energy, and radius of the plate. Interestingly, we observe an power-law dependence between the total damage and the fracture energy of the glass plate. Through an auto-regressive analysis of the evolution of cracks, we demonstrate that the self-affine growth of cracks leads to this power-law dependence. Overall, the present study illustrates how peridynamics simulations can offer new insights into the fracture mechanics of glasses subjected to ballistic impacts. This improved understanding can pave way to the design and development of glasses with improved impact-resistance for applications ranging from wind-shields and smart-phone screens to ballistics.


Introduction:
Glass, an archetypical brittle material, is ubiquitously used in everyday life from car wind shields to smart phone screens (Mauro, 2014;Mauro et al., 2014;Mauro and Zanotto, 2014;Wang et al., 2017a).One of the major causes of failure in glass is the sudden load, such as impact, to which they are exposed to during the service period.These impact loads could vary from small stone pellets to bullets, resulting in a minor crack to a catastrophic failure of the glass structure, respectively (Bobaru et al., 2012;Hu et al., 2013).Understanding the role of material and geometric properties of the structure on impact-induced damage is crucial to design next generation glassware with improved impact-tolerance.
Extensive studies have been carried out in the past decades using experimental techniques on the impact-induced damage in brittle materials (Ball and McKenzie, 1994;Barefoot et al., 2013;Bless and Chen, 2010;Bourne et al., 1995;Chaudhri and Walley, 1978;Field, 1988;Je et al., 1989;Knight et al., 1977).These studies focused both on the dynamic aspects and post-mortem of cracks in typical glass-based structures subjected to impact.However, the experimental studies were limited by two major aspects.Firstly, the crack propagation and associated dynamic phenomena during impact occur at timescales barely accessible to the experiments, even with the fastest available recording devices (Tang et al., 2014).As such, most of the experimental measurements do not have a direct access to the crack propagation and energy releasing mechanisms.Secondly, a parametric study, analyzing the effect of a single material property while keeping others constant, is next to impossible using experiments.This is because, even in materials where continuous variation of compositions are possible, such as glasses or cement, properties do not vary uniformly and continuously (Bauchy et al., 2017;Krishnan et al., 2016Krishnan et al., , 2017;;Li et al., 2017).This poses a fundamental challenge in analyzing critically the role of material and geometric properties on impact-induced failure through experiments.
On the other hand, in silico "numerical experiments" provide a unique means to control the initial and boundary conditions and material properties with high fidelity.Further, it provides complete access to the spatio-temporal evolution of displacement and stress fields, and energy in the structure.To study the impact-induced damage in materials, numerical simulations across the length-scales-from molecular dynamics (MD) (Anders et al., 2012;Holmström et al., 2011;Saiz and Gamero-Castaño, 2012;Samela and Nordlund, 2010) to finite element (FE) (Doner et al., 2019;Krishnan et al., 2010;Peng et al., 2013;Pyttel et al., 2011;Timmel et al., 2007)-have been carried out in the past.Although MD simulations can provide insight into the atomistic mechanisms of crack propagation (Buehler et al., 2003;Wang et al., 2016), they cannot be used to study the damage occurring at the macroscale, e.g., bullet impact.This is attributed to the fact that the system size, in MD simulations, is restricted to a few hundred nanometers, and the timescale to a few nanoseconds.In order to simulate failure at larger length scales, FE simulations are preferred.While traditional FE methods are notorious for their inability to simulate crack, advanced techniques such as extended FE methods (XFEM) (Belytschko et al., 2003;Zi and Belytschko, 2003), FE with element deletion method (Ismail et al., 2012), and cohesive-zone model (Elices et al., 2002) have been proposed as an alternative to investigate macroscopic failure in brittle materials.Each of these models carry their own deficiencies in modeling realistic crack propagations (Agwai et al., 2011;Nayak et al., 2019).The major reason for the inability of FE-based methods to predict realistic crack propagation is that the mathematical basis for these methods rely on the classical continuum mechanics.The basic assumption in continuum theory that the body remains continuous when it undergoes deformation is violated during crack initiation and propagation.This discontinuity in the displacement field, due to crack formation, renders the governing differential equations meaningless.As such, although these methods can be used to study various aspects of impact-induced failure, a start-to-end simulation of impactinduced crack propagation and failure in glasses is impeded due to the aforementioned reasons.
To address this inherent limitation in continuum mechanics, Silling et al. proposed a non-local formulation of continuum mechanics (Silling, 2000).This method, inspired from the MD-based methods (Seleson et al., 2009), proposes an alternate approach to formulate the continuum equations of motion (see the section Methodology).Thanks to the non-local integral-based formulation of peridynamics, the method has been extremely successful to simulate discontinuous phenomena such as crack propagation (Ha and Bobaru, 2011), failure phase separated glasses (Tang et al., 2018), impact (Bobaru et al., 2012;Parks et al., 2008), blast (Fan and Li, 2017;Wang et al., 2017b), and other forms of extreme failures such as landslides (Lai et al., 2015) and earthquakes.
In the present work, we use peridynamics simulations to analyze the damage in glasses subjected to a projectile impact.We investigate the role of mechanical and geometric properties such as elastic modulus, fracture energy, and plate radius in determining the overall damage in the glass structure.An interesting power-law dependence is observed between the overall damage and the fracture energy of the glass plate.Through a detailed analysis on the initiation and propagation of cracks, we show that the power-law dependence is closely related to the crackbranching instability during dynamic crack propagation.

Peridynamics:
Peridynamics is a reformulation of continuum mechanics into integral equations of motion (Le and Bobaru, 2017;Silling, 2000;Silling et al., 2007;Silling and Lehoucq, 2008).This allows peridynamics to handle cracks and discontinuities inherently, as the differentiability of displacement field is not a necessary condition for the equations of motion to be valid (Agwai et al., 2011;Foster, 2009;Ha and Bobaru, 2011;Silling and Askari, 2004).In peridynamics, the continuum domain is discretized using particles or lattice points which interact with each other based on a pair-wise force.A more generalized version of peridynamics with state-based variables has also been proposed to simulate complex material response including plasticity and viscoelasticity (Mitchell, 2011b(Mitchell, , 2011a;;Silling et al., 2007).In peridynamics, the classical assumption of "local force" is replaced by "force acting at a distance" through the introduction of a characteristic length scale.This characteristic length scale, known as the horizon, gives the subdomain within which all particles interact with the central particle.Based on this, the peridynamics equation of motion is given by where ρ is the mass density, ü is the acceleration vector field, f(.) is the pairwise force function for particles at locations x and x', b is the body force intensity, and Hx defines the subdomain in the body, i.e., the horizon.The form of the pairwise force function depends on the material model used in the peridynamic formulation.Here, we restrict our discussion to bond-based material model, as this model has been shown to reproduce the characteristics of fracture in brittle materials successfully (Bobaru et al., 2012;Ha and Bobaru, 2011;Hu et al., 2013).For a detailed overview on peridynamics and different material models, please refer to Refs.(Mitchell, 2011b(Mitchell, , 2011a;;Silling, 2000;Silling et al., 2007).
For a bond-based material namely, prototype microelastic brittle (PMB) material (Parks et al., 2008), the pairwise force is given by where ξ = x'−x is the relative position of the particle, η = u(x',t) − u(x,t) is its relative displacement, c(ξ) is the micromodulus function, and  = ║ + ║ -║║

║𝛏║
gives the relative bond elongation.Using a condition for the equivalency of elastic strain energy density, the micromodulus can be derived from the bulk modulus K, and horizon  as (Ha and Bobaru, 2011) Equation (3) Bonds between the lattice points in PMB are purely elastic with a critical stretch s0, after which they will break.The broken bond corresponds to a local material damage and allows for the initiation and propagation of cracks in the material.The critical bond stretch s0 is derived from the fracture energy of the material so as to represent a realistic strain-energy release rate as (Ha and Bobaru, 2011) where   is the fracture energy of the material.It should be noted that a material particle is connected to multiple particles within the horizon.To account for this while computing damage of a point, a Boolean function (, , ) is defined which has a value of 0 corresponding to a broken bond, and 1 otherwise.Thus, the damage at a material point x is defined as Thus, a material point with all the bonds intact will have a damage of 0, while the one with all the bonds broken have a damage of 1. Peridynamics simulations of impact damage on a glass plate are carried out using the opensource LAMMPS package (Parks et al., 2008(Parks et al., , 2011;;Plimpton, 1995).To this extent, a wellestablished methodology (Bobaru et al., 2012;Ha and Bobaru, 2011;Hu et al., 2013;Parks et al., 2008Parks et al., , 2011) ) which is validated against several experimental studies is used.Figure 1 shows the geometry of the plate and the bullet setup.The plate has a radius of 37 mm and thickness of 2.5 mm.The domain is discretized using a simple cubic lattice with lattice parameter of 0.5 mm.The horizon is kept constant with a value of 2.0 mm, four times the lattice parameter.The value is chosen based on a convergence study considering three different horizons (see Supplementary Material) and is in agreement with an earlier study (Bobaru and Hu, 2012).The material properties of the plate are taken equivalent to that of a soda-lime silicate glass with an elastic modulus of 94.5 GPa and fracture energy of 4.3 J/m 2 (Wang et al., 2016).Note that the material properties and geometry mentioned here are chosen as the standard simulation.In other words, while investigating the effect of material properties and geometry on impact-induced damage, the constant parameters are set to the values as in the standard simulation.The values taken by the variable parameter is detailed in the relevant sections later.

Simulation details:
A spherical bullet with a radius of 3 mm is used in the present study.Further, to ensure the generality of results, two different bullets are considered.Both the bullets are modelled as linearelastic materials.The first one, named as copper, has notably higher stiffness of 115 GPa and fracture energy of 2665 J/m 2 than the target glass plate.This ensures that the bullet suffers close to zero damage upon impact.The second bullet has the same material properties as that of the target plate, and is named as a glass bullet.It should be noted that the behavior of the first bullet may not correspond to that of a real copper bullet.This is because the effects of plasticity are not included in the model.However, the aim of the study is to ensure that the results presented herein are valid irrespective of the stiffness or fracture energy of the bullet.Both copper and glass bullets are projected with a velocity of 100 m/s in the standard case.Initially, the bullet is separated from the glass plate by a distance two times the horizon.This is to avoid any interaction between the bullet and the glass plate at the initiation of the simulation.Peridynamic equations of motion are integrated using a timestep of 10 -7 s (Parks et al., 2008).

Computation of crack velocity
To compute the velocity with which cracks propagate, we focus on the major radial cracks.To this extent, we first plot histograms of damage with respect to angular position of all the lattice points at every timestep.In the histogram, peaks correspond to the major radial cracks and the angles at which these peaks occur correspond to the angular position of the radial crack.In other words, the cracks propagate radially outward at this angle on the glass plate.Now, around these angles at which cracks propagate, we compute the damage in the lattice points along the radial direction.Note that propagation of cracks manifests as increased damage in the lattice points.Based on previous studies, we assume that a point with damage close to 0.5 forms part of a cracked surface (Agwai et al., 2011;Bobaru and Hu, 2012;Foster, 2009).This is based on the assumption that since the "bonds" of a lattice point are uniformly distributed in three dimensions, a failure of 50% of the bonds correspond to an exposed surface.Thus, at each timestep, we compute the radial position of the cracked lattice point to obtain the crack velocity as where () corresponds to the position to the maximum radial position of a lattice point with damage greater than or equal to 0.5, and ∆ is the timestep.It should be noted that these regions with damage greater 0.5 may not correspond to fully separated crack surfaces.It is possible that the damaged regions may still have bridging bonds left across a crack that is yet to break (Xu et al., 2018).However, we observe that the velocity of the crack is not sensitive to the cutoff 0.5 and is consistent with slightly lower and cutoffs (see Supplementary Material).

Autoregressive modeling of damage evolution
To understand the nature of crack-branching and growth upon bullet impact, we model the damage evolution as a recursive time-series, namely, first order auto-regressive (AR) model (Ljung, 1998).Note that the auto-regressive model inherently suggests an exponential growth, as the value at any time t depends on the value at t-1.Thus, the values are compounded at every instant following the growth rate of the AR model.Here, the damage at any time instant t, D(t) is modeled auto-regressively, that is, evolving with respect to damage the previous time instant, D(t-1), as () = ()( − 1) + ( − 1) Equation ( 7) where,  is the unmodeled dynamics or the noise, which follows a Gaussian distribution with zero mean and non-zero variance, and () is the term governing time evolution, called as the damage growth rate.We use data sets that correspond to 5 different fracture energy levels, namely, 4.3, 15, 50, 100, 150 kJ/m 2 for both glass and copper bullets for performing the modeling.To exclude the effects of boundary conditions and other dynamics, we perform the fitting considering the damage evolution in the initial region.
In order to ensure realistic simulation of impact, we compare our simulation with experiments conducted by Knight et al. (Knight et al., 1977) on the impact of steel balls on soda-lime glass.
Here, all the properties and parameters are chosen to match exactly with the ones in the experiment.As such, a spherical ball of radius 0.5 cm is projected towards a soda-lime glass plate with dimensions 10 cm x 10 cm x 5 cm, with different initial velocities varying from 30 m/s to 300 m/s, namely 30,40,60,80,100,110,150,200, and 300 m/s.The Young's modulus and fracture energy of the steel ball are considered to be 200 GPa and 10000 J/m 2 , respectively, and that of the glass plate are 72 GPa, and 7 J/m 2 , respectively.For the peridynamics simulations, the discretization of the glass plate and ball is carried out with lattice units of 0.1 mm and a constant horizon of 0.4 mm is maintained.The micromodulus and critical stretch are computed using Equations ( 3) and (4).Note that both the steel ball and the glass plate are considered to isotropic.First, we focus on the coefficient of restitution (CR) of the bullet.Note that CR, being a ratio of approaching velocities to separating velocities, essentially captures the energy dissipated in the impact process.Thus, CR ranges from 0 to 1, the former representing a perfectly plastic collision with complete energy dissipated in crack formation and the latter representing a perfectly elastic collision with no energy dissipation.Figure 2(a) shows the variation of CR with respect to initial velocity of the bullet.For low velocities of the bullet (< 100 m/s), CR remains constant and close to 1, suggesting an elastic collision for this velocity range.For higher velocities (> 100 m/s), CR reduces drastically suggesting energy dissipation in the collision, almost following a power-law with respect to the initial velocity.This is consistent with the experiments (Knight et al., 1977) where CR is observed to be close to 1 for velocities lower than 100 m/s with little damage on the glass plate.Figures 2(b) and 2(c) show the cross-sectional view of the glass plate subjected to impact in experiment (Knight et al., 1977) and simulation, respectively, at the same time (4 μs).
We note that the depth of the damage observed in both the simulation and experiment is similar.Further, both the simulation and experiment produce debris upon impact, which is flying away from the surface.However, we observe that the peridynamics simulation over-estimate the pulverization behavior in the glass plate in comparison to the experiment.This could be due to the inherent deficiency of particle-based methods, which exhibits an increased tendency for particle separation.Note that the damage pattern simulated by peridynamics exhibits an asymmetry (see Fig. 2(c)), which could be ascribed to the following reasons.(i) The bullet atom is not perfectly spherical due to the discretization used in peridynamics.(ii) A slight perturbation is provided for the lattice points used in the discretization to avoid the formation of perfectly symmetric crack patterns.As such, the impact on the glass plate may result in an asymmetrical damage pattern.Overall, we observe that the peridynamic simulations exhibit a reasonable match with respect to experiments in predicting impact induced damage.the glass plate subjected to impact at four different times, namely, 0, 15, 17, and 30 0 μs.The coloring scale, representing the local damage, varies from 0% (blue) to 100% (red).Now, we focus on the crack formation behavior when a bullet impact occurs on a brittle glass plate.The glass plate has a radius of 0.037 mm, thickness of 0.0025 m, elastic modulus of 94.5 GPa, and fracture energy of 4.3 J/m 2 and is subjected to the impact of a copper bullet with a velocity of 100 m/s.The overall fraction of cracks present in the material can be quantitatively represented in the peridynamics framework by the average damage per particle.An average damage of 100% corresponds to all the lattice points broken and separated away from each other and 0% correspond to all the points intact.

Crack formation under ballistic impact:
Figure 3(a) shows the evolution of the average damage per particle and total mechanical energy in the glass plate.We observe that the average damage increases quickly upon impact and then eventually becomes constant, close to 48%, representing a saturated state.Note that the saturation of damage occurs very quickly-around 40 s after the bullet comes into contact with the plate.Now, we focus on the evolution of the total mechanical energy gained by the plate upon impact.Note that this energy includes both the elastic energy and kinetic energy of the plate.From Fig. 3(a), we note that the energy initially increases, which is attributed to the elastic potential energy of the plate due to deformation and kinetic energy due to the transverse and longitudinal waves produced.Subsequently, the energy reaches a peak value, and then reduces slightly before reaching saturation.The slight decrease in energy before saturation is attributed to the release of elastic energy by the formation of a network of cracks.Note that the saturation of energy occurs concurrent to the saturation of damage.A visual representation of the bullet impact on the glass plate, followed by the crack formation, is shown in Fig. 3(b).The crack pattern observed is qualitatively similar to the ones observed in earlier experimental studies (Chaudhri and Walley, 1978;Knight et al., 1977), confirming the realistic nature of the simulations.We note that, upon impact, the cracks are first formed in the plate at the point where the bullet strikes (see Fig. 3(b)).Major radial cracks then form which, while spreading outward, branch to form minor cracks.These cracks then spread to the entire glass plate forming a network of cracks in the form of radial and tangential cracks.Overall, we observe that the energy gained by glass plate upon impact is released through the formation of cracks and debris.Note that the results correspond to that of a glass plate with a radius of 37 mm, thickness of 2.5 mm, elastic modulus of 94.5 GPa, and fracture energy of 4.3 J/m 2 .

Effect of bullet material and velocity:
In order to ensure the generality of the crack formation, we consider two different bulletscopper and glass.The copper bullet is significantly stiffer and has a fracture energy which is four orders of magnitude higher than that of the glass plate.Thus, upon impact, the copper bullet remains perfectly intact.On the contrary, the glass bullet shares the same material properties as the glass plate.During the impact, the bullet shatters along with the glass plate.From Fig. 4(a), we note that the evolution of the average damage per particle exhibits similar trend in the case of both copper and glass bullets.Overall, we note that, given the geometry and material properties of the plate remain constant, the crack formation behavior in the glass plate is similar irrespective of the bullet material.Now, we focus on the effect of the bullet velocity on the overall damage on the glass plate.Figure 4(b) shows the maximum percentage damage of the glass plate when impacted by bullets of different velocities.To this extent, six different bullet velocities are considered, viz., 5 m/s, 10 m/s, 15 m/s, 25 m/s, 50 m/s, and 100 m/s.Note that the maximum percentage damage in each case corresponds to the saturated value of the damage after impact.We observe that the overall damage increases with the bullet velocity.While the increase is minor for lower velocities (< 25 m/s), it is notable higher for higher velocities.Realistically, a velocity of 5 m/s would be similar to that of a stone striking the wind shield, while 100 m/s would be closer to that of a bullet impact.Further, we note that the velocity of 100 m/s creates maximum damage, among the selected velocities, due to high energy during impact.Therefore, we fix the bullet velocity as 100 m/s for the following studies.To investigate the effect of plate geometry, we vary the plate radius.We focus on the effect of plate radius by considering four different radii, viz., 20 mm, 37 mm, 50 mm, 100 mm.The projectile is modeled using both glass and copper bullets.Figure 5 shows the final cumulative damage (left axis) and average damage per particle (right axis) with respect to the plate radius.
The cumulative damage initially increases with increasing plate radius up 0.05 m, after which it tends to plateau.Note that a slight decrease in the cumulative damage is observed in the case of both glass and copper bullets for a radius of 0.1 m in comparison to 0.05 m.The corresponding average damage per particle is also shown in Fig. 5.For a plate radius of 20 mm, ~80% final average damage is observed.Further, the average damage per particle decreases monotonically with increasing plate radius, eventually reaching ~8% for a radius of 0.1 m.The results are comparable in the case of both copper and glass bullets, although copper bullets cause slightly increased damage.
Reconciling the cumulative and percentage damage, we note that for plate with smaller radii, the energy gained during impact is not completely dissipated by crack formation due to its smaller size.With increasing plate radius, the overall size increases, ultimately leading to the formation of more cracks, as suggested by the increasing cumulative damage.The saturated state of cumulative damage corresponds to the state wherein the maximum area of cracks that can be formed with the impact energy has formed.In the present case, this occurs for a radius of 0.05 m.Further, with increasing plate radius, increasing kinetic energy is absorbed by the plate for dynamic motion leaving lower energies for crack formation.As such, the cumulative damage decreases slightly after the radius of 0.05 m.

Effect of material properties:
Figure 6.Average damage per particle of the glass plate with respect to the elastic modulus of the plate.
To understand the effect of the material properties, we focus on the elastic modulus and fracture energy of the plate.First, we investigate the effect of elastic modulus by considering eight moduli values spanning over two orders of magnitude, viz., 1 GPa, 5 GPa, 10 GPa, 25 GPa, 50 GPa, 75 GPa, 100 GPa, and 150 GPa. Figure 6 shows the average damage per particle in the glass plate with respect to the elastic modulus when impacted by glass and copper bullets.Note that all other parameters including the plate geometry, fracture energy, and bullet velocity are kept constant during the study.We observe that the average damage per particle remains fairly constant over the entire range of elastic moduli in the case of the copper bullet.In the case of the glass bullet, there is a slight increase in damage of the plate with increasing elastic modulus.However, the percentage damage is restricted within a range of 43% to 49% for an increase in the elastic modulus by two orders of magnitude.It should be noted that the elastic modulus is related to the fracture toughness of the material.An increase in elastic modulus would result in an increase in the fracture toughness of the material as well.This suggests that the stress required for the initiation of damage would increase with increasing elastic modulus.However, since the fracture energy of the plate-that is, the energy released upon the formation of cracks-is constant, the overall damage is not significantly affected by the variation in the elastic modulus.
Figure 7. Final average damage per particle in the glass plate with respect to fracture energy for copper and glass bullets.The lines represent a power-law fit in the least-squared sense.Now, we focus on the effect of the fracture energy of the glass plate.To this end, we consider five different fracture energies spanning over two orders of magnitude, viz., 4.3 J/m 2 , 15 J/m 2 , 50 J/m 2 , 100 J/m 2 , and 150 J/m 2 .As in the case of the elastic modulus, other parameters are kept constant during the study.Figure 7 shows the final average damage per particle in the glass plate with respect to the fracture energy.Interestingly, we observe an power-law dependence of the average damage with respect to the fracture energy of the glass plate for both copper and glass bullets.It should be noted that, in peridynamics, the overall damage is proportional to the new surface area created by the cracks (Silling, 2000;Silling and Askari, 2004).Further, in the case of thin plates subjected to impact, the surface area is proportional to the overall length of the cracks, including all the major and minor crack branches.In other words, the overall damage is proportional to the total length of the cracks formed upon impact.This suggests the existence of a power-law relationship between the total crack length and the fracture energy of the material.
In other words, the total crack length decreases exponentially with increase in fracture energy.
To establish the power law, we apply a logarithmic fit (see Fig. 7), minimizing the least-squared error, of the form,  =   − Equation (1) where,  represents the percentage damage (%),  is a pre-factor,   is the fracture energy (J/m 2 ), and  is the power-law exponent.The value of the power-law exponent is found to be the same, 0.526 in this case, for both glass and copper bullets.This suggests a universal nature for the relationship between total crack length and fracture energy in plates subjected to impact, which is similar to the self-affine nature of cracks in fracture (Borodich, 1997).We now discuss the power-law relationship observed between the fracture energy and crack length in the glass plate.To this extent, we focus on the radial cracks that grow from the center outwards when the bullet impacts the plate (see Fig. 3(b)).Note that the velocity of a crack depends on the energy flowing to the crack tip and the fracture energy required to create surface.However, the maximum velocity with which the cracks can propagate is limited for any given material, which corresponds to the Rayleigh wave speed VR of the material (Sharon and Fineberg, 1999).In brittle materials, a crack tip instability has been observed above a critical velocity vc, that is 0.36VR, which leads to local crack branching (Sharon et al., 1995).Such branching reduces the velocity of the cracks due to increased energy dissipation.However, as soon as one of these branches terminate, the major radial crack again accelerates leading to further branching.Thus, the crack tip velocity exhibits a cyclic behavior resulting in an increased energy dissipation (Sharon et al., 1995).Here, we suggest that two complementary mechanisms are at the origin of the power-law dependence.On the one hand, for low fracture energies, the radial cracks that propagate outwards undergo crack branching due to increased crack velocities leading to crack tip instability (Fineberg et al., 1991;Sharon et al., 1995;Sharon and Fineberg, 1999).This leads to multiplication of cracks resulting in an exponential growth of the damage with a high growth rate.On the other hand, increased fracture energy decreases the crack velocity.This reduces the crack branching and the growth rate resulting in a decreased cumulative damage.

Discussion:
To establish the hypothesis, we first analyze whether the crack tip velocity exhibits a cyclic behavior representative of the crack tip instability or not.To this extent, we quantify the distribution of cracks throughout the glass plate using a polar plot as shown in Fig. 8(a).The rest of the histogram with comparable peaks represent the axisymmetric crushing of the plate, along with some tangential cracks.Now, we compute the velocity of the major radial cracks formed during impact (see Methodology for details).Note that the intention here is to observe the relative variation in the crack tip velocity rather than the exact value itself.Therefore, we normalize the velocity values with the maximum value Vmax obtained during the simulation.Figure 9(a) shows the normalized velocity of propagation of one of the major radial cracks in the glass plate with a fracture energy 4.3 J/m 2 .We observe that the velocity of the crack tip exhibits an oscillatory behavior before finally going to zero, thereby confirming the existence of the microbranching instability.Moreover, the minimum value during the oscillation is fairly close to 0.36Vmax.This suggests that the crack branching reduces the velocity of the crack up to 0.36Vmax, after which it starts accelerating again.Figure 9(b) shows the maximum velocity attained by a major radial crack during the impact with respect to the fracture energy of the glass plate.We observe that the maximum velocity attained by the major radial crack decreases monotonically, thereby further affirming our hypothesis.Now, we model the damage evolution in a glass plate using an AR-based model (see Methodology).Note that an AR-based model represents a self-affine transformation essentially leading to an exponential growth of the damage.The crack branching behavior is quantified using the damage growth rate  (see Methodology). Figure 9(b) shows the damage growth rate with respect to the fracture energy of the glass plate.Inset shows the AR fit compared to the measured evolution of damage for a glass plate with a fracture energy of 4.3 J/m 2 .We observe that the AR exhibits an excellent fit for the evolution of the damage confirming the exponential growth.
Further, the damage growth rate monotonically decreases with increase in fracture energy.This is qualitatively in agreement with the trend exhibited by the maximum crack velocity, which decreases monotonically with increasing fracture energy of the material.This suggests that there is a decreased tendency for crack branching with increasing fracture energy, as evident from the growth rate.Finally, Fig. 10 shows the saturated damage map in the glass plates, with varying fracture energies, subjected to bullet impact.In agreement with the hypothesis, we observe visually that crack branching is increasingly suppressed with increase in fracture energy.Overall, through independent analysis of crack tip velocity, maximum crack velocity, and the AR growth rate, we corroborate that the crack branching and growth behavior is at the origin of the powerlaw dependence.

Conclusions:
In the present study, based on peridynamics simulations, the roles of geometric and material properties on the impact-induced damage in a glass plate are explored.To this extend, approximately hundred high-throughput peridynamics simulations of bullet impact on a glass plate are performed.This methodology, validated against experimental results earlier (Bobaru et al., 2012;Ha and Bobaru, 2011;Hu et al., 2013), provides a unique way to analyze to the effect of each of the properties individually which is otherwise impossible using experiments.We observe that the plate geometry does not play a major role in the overall damage as long as the geometry is large enough to accommodate energy dissipation by the formation of cracks.
Similarly, the elastic modulus of the material does not have a significant effect on the overall damage, although the damage tends to increase slightly with increasing values of modulus.On the contrary, the fracture energy of the glass plate directly controls the overall damage through a power-law relationship.Through a detailed analysis on crack evolution using auto-regressive model and crack-propagation velocity, we demonstrate that the power-law originates from a self-affine growth of the major radial cracks resulting from a crack tip instability.Overall, the present study suggests that peridynamics can offer new insights into the fracture behavior of glasses upon ballistic impact, which opens the door toward the development of highperformance bullet-proof glasses for applications in next generation automobile windshields to armors in ballistic applications.

Figure 1 .
Figure 1.(a) Front and (b) top view of the plate-bullet simulation setup.Note that the glass plate and bullet are colored in blue and red, respectively.

Figure 2 .
Figure 2. (a) Simulated coefficient of restitution (filled triangles) of the impact compared with the previous experiments by Knight et al. (Knight et al., 1977) (empty squares).The dotted line is a guide to eye.(b) Impact of steel sphere on soda-lime glass with an impact velocity of 200 m/s by Knight et al.(Knight et al., 1977).(c) Simulation of the impact experiment by Knight et al.(Knight et al., 1977) using peridynamics.Red color represents damaged particles.

Figure 3 .
Figure 3. (a) Average damage per particle (left axis, red) and total mechanical energy (right axis, blue) of the glass plate, subjected to an impact by a copper bullet, with respect to time.Note that the results correspond to that of a glass plate with a radius of 37 mm, thickness of 0.0025 m, elastic modulus of 94.5 GPa, and fracture energy of 4.3 J/m 2 .(b) Damage map in

Figure 4 .
Figure 4. (a) Average damage per particle of the glass plate with respect to time upon impact with copper and glass bullets having a velocity of 100 m/s.(b) Average damage per particle, after saturation, of the glass plate with respect to bullet velocity for copper and glass bullets.Note that the results correspond to that of a glass plate with a radius of 37 mm, thickness of 2.5 mm, elastic modulus of 94.5 GPa, and fracture energy of 4.3 J/m 2 .

Figure 5 .
Figure 5. Cumulative damage (left axis) and average damage per particle (right axis) of the glass plate with respect to plate radius for both copper (square symbol) and glass (circle symbol) bullets.

Figure 8 .
Figure 8.(a) Distribution of cracks in the glass plate in a polar coordinate system, (b) histogram of the distribution of cracks with respect to the angular position θ°.Note that the peaks indicate the major radial cracks.
Figure8(b)  shows the histogram of the angular distribution of cracks in one of the glass plates subjected to impact.The distribution of radial cracks can be obtained from the peaks of the histogram in the angular distribution.We note that the major radial cracks are represented by the taller peaks , while the minor radial cracks and branches are represented by smaller peaks.

Figure 9 .
Figure 9. (a) Evolution of normalized crack velocity of a major radial crack in a glass plate with a fracture energy 4.3 J/m 2 for a copper bullet.Note that the crack velocity is normalized with respect to the maximum crack velocity during the propagation of the crack.(b) Maximum crack velocity (left axis) and damage growth rate (right axis) with respect to the fracture energy of the glass plate.INSET shows the AR fit compared to the measured evolution of damage in peridynamics simulation of impact on a glass plate with fracture energy 4.3 J/m 2 .

Figure 10 .
Figure 10.Damage map in the glass plates with different fracture energies, namely, (a) 4.3 J/m 2 , (b) 15 J/m 2 , (c) 50 J/m 2 , (d) 100 J/m 2 , and (e) 250 J/m 2 , subjected to impact.(f) The coloring scheme used to denote the local damage in percent.Note that 100 corresponds to aparticle with all the bonds broken, and 0 to a particle with all the bonds intact.