Phenomenology of the Maximum Fragment Mass Dependence Upon Ballistic Impact Parameters

Molecular dynamics simulations of the ballistic Taylor test are used to explore correlation between the largest fragment mass and the impact energy of a projectile as well as a set of selected state variables. Flat-ended, monocrystalline, nanoscale bars collide with a rigid wall with striking velocities ranging from 0.27 km/s to 60 km/s. The investigation emphasis is on two border regions of the emerging nonlinear phenomenological model identified with two transitions: the damage-fragmentation transition and the shattering transition. In between these two nonlinear regions, the maximum fragment mass is largely inversely proportional to the impact energy, and the maximum values of the pressure, temperature, and the square of the effective strain. A reverse-sigmoid phenomenological model is proposed to capture the unifying features of this nonlinear and saturable dependence. A crystallographic orientation dependence of the damage-fragmentation transition parameters is investigated.


INTRODUCTION
The high-velocity ballistic Taylor test (Taylor, 1948) is a time-honored procedure of exploration of the dynamic response of materials.A series of two-dimensional (2D) traditional MD (molecular dynamics) simulations of this classic experiment is performed in this study by using nanoscale projectiles made of the Lennard-Jones 6-12 (LJ) monocrystalline solid, under a tacit assumption that this, admittedly rather simple, model is sufficient to capture some essential features of the investigated phenomenon.An extension and refinement of an earlier analysis (Mastilovic, 2015a) result in collision of nanoscale projectiles with a rough rigid wall with impact velocities (v) varying in a wide range from 0.27 km/s to 60 km/s that reveal a nonlinear and saturable character of the maximum fragment mass dependence upon selected ballistic parameters.The resulting reverse-sigmoid phenomenological model, suggested in the present article, consists of two nonlinear border regions identified with two phase transitions: the damage-fragmentation transition (v = v0) and the shattering transition (v = v1), with largely linear region in between.The accompanying hyper-exponential fragment mass distribution typical of instantaneous fragmentation of the ductile solids was discussed earlier (Mastilovic, 2015a) within a narrower striking velocity range.The ultrafast flat-end collision of the projectile with the rigid target is an extremely intense loading event belonging to the realm of akrology within the study of materials' physics.The extremely steep gradients of state variables (e.g., Mastilovic, 2016a), well documented by shock experiments, may cause phase transitions and lead to sequential fractures that culminate eventually in energetic expulsion of fragment debris.Since a shock wave excitation is inherently ultrafast, the present MD method requires a femtosecond time resolution to observe the collective dynamics of material on picosecond time scales, which renders simulations extremely time-consuming.
The fundamental principles of dynamic fragmentation of solids were investigated extensively, both experimentally and theoretically, and the substantial literature is compiled, among others, by Grady (2006), Elek and Jaramaz (2009), Ramesh et al. (2015).Investigations of small-scale fragmentation induced by hypervelocity impact are relevant for evaluation of the risk of space debris and dust impacts in earth's orbit.Although perhaps miniscule in size, space debris collisions (with striking velocities v > 10 km/s, up to 70 km/s for micrometeoroids) are environmental factors of growing concern since they can deliver sufficient impact energy to "compromise or deplete" functionality of space assets (Lamberson et al., 2012).Since the extreme loading rates in conjunction with the nonlinear constitutive relationship render an analytical treatment of the problem extremely difficult (Alves and Yu, 2005), various computational methods are often employed-in addition to experimental techniques-to gain insight into salient features of the impact and dynamicfragmentation phenomena.The small-scale 2D-MD simulations were used frequently in the last three decades to study these irreversible, nonlinear, nonlocal, and far-from-equilibrium processes.Holian and Grady (1988) were first to use MD to explore the fragmentation phenomena by simulating a homogeneous adiabatic expansion of condensed matter.Similar computation techniques were adopted subsequently to explore the 2D explosive fragmentation (Diehl et al., 2000, Astrom, 2000) and the fragmentation of grooved target under flyer-plate impact (He et al., 2015).Recently, the small-scale 2D-MD simulations were also used by Sator and his collaborators (Sator et al., 2008, Sator andHietala, 2010) to investigate generic behaviors and damage evolution in the instantaneous point fragmentation of the LJ brittle solid colliding with a wall.In addition to MD, examples of computational techniques utilized recently in the dynamic fracture and fragmentation investigations include particle models (Baker and Warner, 2012;Kumar and Ghosh, 2015), discrete element models (Wittel et al., 2008;Iturrioz et al., 2009;Timar et al., 2010Timar et al., , 2012;;Paluszny et al., 2014), finite element methods (Levy and Molinari, 2010;Ugrcic, 2013), and meshfree methods (Wu et al., 2014;Li et al., 2012Li et al., , 2015)).
The investigation of dependence of the maximum fragment mass upon various state variables is extended and considerably refined in the present article in comparison to the preceding investigation which resulted in the piecewise-linear approximation (Mastilovic, 2015a).First, the analysis of the onset of the damage-fragmentation transition for slender projectiles, based on the maximum and Latin American Journal of Solids and Structures 14 (2017) 1529-1546 average fragment mass, is refined to identify more transparently this continuous phase transition (Kun and Herrmann, 1999;Timar et al., 2012).Second, an attempt is made to capture so-called shattering transition (Kun and Herrmann, 1999) -the elusive terminal fragmentation defined by the uniformly monatomic debris (mmax ≡ 1) and predicted by the piecewise-linear model to correspond to the impact velocity of approximately 45 km/s.Most importantly, emerging unifying features of the nonlinear and saturable dependence of the maximum fragment mass upon the set of impact parameters and ballistic state variables are captured by a single reverse-sigmoid phenomenological model.

COMPUTER SIMULATION TECHNIQUE
The present investigation is based on the traditional MD in which the dynamic state of the atomic system is defined by laws of classical mechanics with atomic motions being uniquely determined by an empirical potential (Allen and Tildesley, 1996).The model is described in detail in preceding studies (Mastilovic, 2015a(Mastilovic, , 2016a(Mastilovic, , 2016b)), thus, a succinct summary is deemed sufficient herein.A monatomic system is comprised of atoms of equal masses mi = m0 that form an ideal defect-free triangular lattice (without any quenched disorder) and interact with their nearest neighbors according to the LJ potential to mimic a monocrystalline, flat-nosed projectile.The three LJ model parameters used to match, as close as possible, physical properties of tungsten (74W) are the atomic mass m0 = 3.1×10 -25 kg (183.85 u), the atomic radius 1.4 Å (≡ r0/2 where r0 is the equilibrium interatomic distance), and the strength of attraction ε = 7.5×10 -20 J.The coordination number of bulk atoms in the reference configuration is six and the potential energy per atom is -2.96 ε, which is slightly in excess of the bulk value (-3 ε) due to the surface effects.Since the shock wave excitation is inherently ultrafast, the Cauchy problem is solved numerically by using the Verlet algorithm with the time step of the order of femtoseconds (Mastilovic, 2016a).This extremely small time resolution (required by the ultrahigh power of the simulated event) in conjunction with necessity to approach asymptotically steady states of fragment mass distributions, makes the MD simulations extremely time-consuming (even for the relatively small model size) and effectively limits the maximum achievable striking velocity.
The link between two atoms ruptures when their mutual distance exceeds a predetermined critical value.The cut-off interatomic distance, R ≈ 1.7 r0, is selected herein to be between the first and second nearest neighbors in the reference configuration (the perfect crystal prepared at zero temperature).A fragment is defined as a self-bound cluster of atoms with interatomic distance less than the cut-off distance (rij ≤ R) in a sequential atom-by-atom search for the nearest neighbors (Mastilovic, 2015a(Mastilovic, , 2016a)).
The fragmentation model proposed herein is generic in the sense that it aims to capture the underlying features of the investigated phenomenon.Its simplicity rests primarily on the 2D geometry, LJ potential, and nanoscale projectile dimensions.It is already noticed that the MD fragmentation Latín American Journal of Solids and Structures 14 (2017) 1529-1546 simulation is a slow process, difficult to reach steady state configuration (Holian and Grady, 1988;Astrom et al., 2000;Diehl et al., 2000;Mastilovic, 2015a).Therefore, although the dimensionality of the system is known to influence shock physics and the universality classes of fragmentation phenomena (Timar et al., 2010), the 2D choice is necessitated by extremely laborious MD computations and justified by a qualitative character of the study.1 Furthermore, it has been recently demonstrated by a finite size scaling approach (Sator and Hietala, 2010) that-for a similar MD simulation technique-"the fragmentation features of the system are not sensitive to the number of particles" and that generic behaviors seem to be shared by fragmenting systems regardless of the details of their interaction potentials.Nonetheless, the size of the present model exceeds those used in the recent point-impact studies (Sator and Hietala, 2010;Timar et al, 2012) that similarly utilize 2D-MD simulations to investigate universality and generic behavior in the impact fragmentation.Last but not least, the plasticity in confined dimensions is a fascinating and rapidly developing research area in itself at present (e.g., Kraft et al., 2010;Greer and De Hosson, 2011;Rinaldi, 2011).

OBSERVATIONS AND DISCUSSION
The knowledge of the maximum fragment mass dependence upon the impact energy (the initial kinetic energy) of the projectile (K  v2 ) is of obvious interest for engineering applications.While the mean fragment mass illustrates the average character of the fragmentation process, the maximum fragment mass, not so frequently encountered in the existing literature, is potentially of considerable importance for structural-survival and risk analyses since it provides a lower bound for definition of the secondary-impact design events.The simulation data presented in Figure 1a suggests that, beyond the fragmentation onset velocity (v  vfo), 2 the maximum fragment mass is inversely proportional to the kinetics energy of the projectile.A piecewise-linear approximation for the maximum fragment mass dependence upon the striking velocity is proposed by Mastilovic (2015a) for the hypervelocity impact range up to vL  30 km/s.Based on this approximation, the terminal fragmentation (mmax ≡ 1) is anticipated to be v = v1 = 45 km/s.The investigation of two nonlinear border regions is refined in the present article by a set of additional simulations.A part of this effort is directed toward the upper-end of hypervelocity impact range to explore limits of validity of the linear extrapolation.Notably, the elusive terminal fragmentation is not achieved even for the v = 60 km/s.
In addition to the striking velocity (the impact energy), three macroscopic observables explored henceforth in connection to the maximum fragment mass include: the average normal stress, P = (σx + σy)/2; the instantaneous kinetic temperature T (Wagner et al., 1992;Mastilovic, 2016a); and the effective strain, . The evolutions of these three 2D state variables are recorded at twelve evaluation areas mimicking the measurement gages (details are available in Mastilovic, 2015aMastilovic, , 2015b)).The first mechanical stress invariant, P, is a measure of the force interaction be-Latin American Journal of Solids and Structures 14 (2017) 1529-1546 tween material points inside averaging areas while the instantaneous kinetic temperature, T, is a measure of the intensity of vibratory motion.Importantly, the mechanical stress, defined by interatomic forces and atomic positions (Zhou, 2003), becomes physically ill-defined as a measure for the mean mechanical force between material points when the averaging area, in the course of projectile distortion and fragmentation, becomes incompletely occupied by atoms.It has been verified by the present simulations that the maximum values of the average normal stress (Pmax) reported henceforth were achieved much before this took place for every single evaluation area.On the other hand, it should be noted that the definition of the instantaneous kinetic temperature (Wagner et al., 1992;Mastilovic, 2016a) is not based intrinsically on the space averaging over a certain evaluation area but rather on averaging over all atoms belonging to the evaluation area; which makes it a state variable less sensitive to distortion and fragmentation than the mechanical stress. 3s an example, the scaling relation of the form Latín American Journal of Solids and Structures 14 (2017) 1529-1546 is obtained with a reasonable confidence based on the simulation results presented in Figure 1b for the lower and intermediate part of the hypervelocity impact range (vfo < v < vL).As pointed out by Mastilovic (2015a), the scaling relation ( 1) is in agreement with the prediction of the Grady's classic model of dynamic fragmentation due to shear banding in the shock-compressed ductile materials.Similarly to the functional dependence of the maximum fragment mass upon the striking velocity (Figure 1a), these simulation data suggest the divergence from linearity at vL  30 km/s.
The logarithmic plots of the simulation data presented in Figures 1c and 1d also indicate that the maximum fragment mass is inversely proportional to the maximum values of the pressure and the instantaneous kinetic temperature within the lower and intermediate part of the hypervelocity impact range, which implies linearity between pressure and temperature in the ejected plasmatic debris in agreement with the classic ideal gas law.(The maximum values of the selected state variables represent the arithmetic mean of the values evaluated at the evaluation points positioned at the projectile centerline.) The following scaling relation captures the elucidated linearity observations within the lower and intermediate part of the hypervelocity impact range, with ξ ≈ 1 for the generic state variable Ξ  {K, Pmax, Tmax, 2 max  }.Accord- ing to the simulation results presented in Figure 1, the upper bound for validity of the scaling relation (2) is vL  30 km/s.
The dependence of the maximum fragment mass upon the aforementioned state variable Ξ for the entire (non-negligible, Ξ  Ξ0) fragmentation range explored in the present study is schematically depicted in the logarithmic space by the following expression (3) and illustrated by Figure 2. The nonlinear and saturable empirical formula (3) disregards the minuscule fragmentation below the onset of the damage-fragmentation transition, Ξ < Ξ0 (Figure 1a).In (3), the three uppercase alphabetic letters (A,B,C) designate the fitting parameters (for example, in Figure 1a, A=2.2, B=7, C=4).The parameters denoted by subscript "0" correspond to the fragmentation threshold (Ξ0, mmax0) as illustrated in Figures 1 and 2. The damage-fragmentation transition is a continuous phase transition discussed recently in literature for the bulky impactors (e.g., Kun and Herrmann, 1999;Timar et al., 2012).For a slender projectile of any given aspect ratio, the onset of the damage-fragmentation transition is dependent upon its cross sectional dimension.The impact fragmentation in confined spatial dimensions is expected to result in increased fragmentation thresholds due to the small-scale hardening of the material as indicated, for example, by Rinaldi (2011).Timar et al. (2012) and Mastilovic (2016b) proposed recently scaling forms of the critical velocity v0 in terms of the system size and determined critical velocities of the infinite system.
It is obvious from expression (3) and Figure 2 that the three fitting parameters A-C are not mutually independent but need to satisfy the following condition from the approximately linear domain The inflection point ( 42) is uniquely determined from the condition that it corresponds to the second derivative of (3) being equal to zero.
Finally, an alternative form of the nonlinear phenomenological model for the maximum fragment mass can be derived straightforwardly by combining ( 2) and ( 4).
The empirical expressions (3) and ( 5) cannot formally capture the terminal fragmentation, mmax ≡ 1, except as the limit case: defining the horizontal asymptote illustrated in Figure 2. As previously mentioned, the lower bound of the shattering-transition threshold is extrapolated to v1 ≈ 45 km/s based on the piecewise-linear approximation (Mastilovic, 2015a).It has been emphasized in the same article that only a miniscule part of the impact energy of the hypervelocity impact is spent on the fracture process, and noted that it is not obvious that any impact velocity in the present simulation setup would actually result in the terminal fragmentation (all fragments are the smallest unit size fragments; mmax = 1 ≡ mmax1).The additional MD simulations performed herein at three ultrahigh striking velocities v = (45, 50, 60) km/s reach the fragmentation steady state with mmax = (3, 2, 2) > mmax1 ≡ 1, respectively.Thus, the elusive terminal fragmentation of 15×110-nm 74W projectile is not achieved even after 900,000 time steps (90 ps).The nonlinear data fit illustrated in Figure 1a points roughly to v1 > 150 km/s, but such increase of the striking velocity would require the reduction of the time step that precludes the effort to clarify the issue at present.Katsuragi et al. (2004) recognized, in their experimental study of 2D flat-impact fragmentation, that although this terminal fragmentation state is extremely difficult to achieve it can exist as an ideal limit case that clearly defines the saturation state.
Furthermore, the simulation-data fitting by ( 5), illustrated by solid lines in Figure 1, reveals not only that the fitting parameters B and C have the same values for all four impact parameters and state variables Ξ  {K, Pmax, Tmax, Thus, (5) can be rewritten in a form (Note that based on simulation results presented in Figure 1, D is estimated to be 0.20 for Ξ and 0.45 for v K  .)Consequently, the nonlinear phenomenological model for the maximum fragment mass (8) is captured by the same curve when the argument, ln Ξ, is shifted by A as illustrated by Figure 3.A along the abscissa in the hypervelocity impact range (where A  {-35.9, 20.3, 5.0} in respective SI units).Structures 14 (2017) 1529-1546 3.1 A geometric estimate of the orientation dependence of m max0 in perfect monocrystals

Latin American Journal of Solids and
The maximum fragment mass corresponding to the fragmentation threshold (mmax0) is the model parameter that determines the upper horizontal asymptote defining the baseline fragmentation response.Any specific application of empirical relations (3, 5, 8) hinges on the determination of the critical point (Ξ0, mmax0), considering the saturable character of the proposed phenomenological dependence and the uniquely-defined: (i) the lower horizontal asymptote (ln mmax1 = 0) and (ii) the linear intermediate hypervelocity range with slope ξ  1 (Figure 3).
In the present 2D-MD framework, it is intuitively clear that mmax0 (being the fragmentation outcome of the minimum impact energy that is sufficient to brake off a pair of fragments) should be highly dependent on the geometry of the problem; most notably: (i) the global symmetry and (ii) the crystal orientation of the slender monocrystalline impactor.Figure 4 illustrates schematically one typical instance in the time sequence of stress fields generated in the Taylor projectile by interaction of the reflected compressive waves (c) and the lateral release tensile waves (t) (adopted from Grady and Kipp, 1993;also Mayers, 1994).An outcome of this complex stress field is the lateral mass transfer in the process zone inherently related to the damage (d) nucleation and evolution.Consequently, the minimum energy sufficient to disconnect a pair of fragments is assumed herein to result in the fragment area (the volume of unit thickness) shaded yellow in Figure 4, inspired by the damage-evolution boundaries implied by the snapshots such as those presented in Figure 5.These snapshots of the highly-distorted projectiles reveal typical features of the ductile fracture with voids that nucleate, grow and coalesce by plastic deformation (Grady and Kipp, 1989;Woodward et al., 1992;Meyers , 1994, p. 490).(The relative thickness of contour lines reflect qualitatively the field intensities.).Structures 14 (2017) 1529-1546 The perfect triangular lattice used in the present investigation, due to its 60° rotational symmetry, offers two distinct symmetric directions of the highest linear atomic density on the given closed-packed crystallographic directions that are convenient for the following rough geometric estimate of mmax0 and its orientation dependence.These two crystal orientations are defined by ±60° and ±30° angle of the closed-packed direction with respect to the impact direction (denoted by λ in Figure 4 and outlined in Figure 5).

Latín American Journal of Solids and
Thus, based on the projectile geometry and the stress field illustrated in Figure 4, the maximum fragment area (volume of unit thickness) at the brake-off moment is  10) provide a rough estimate of the maximum fragment mass immediately after the brake-off; the shock induced process of thermal attrition reduces somewhat the fragment until the stable configuration is asymptotically reached.Due to the 60° rotational symmetry of the perfect triangular lattice, there are two symmetric impact configurations corresponding to λ = (30°, 60°).The diameter of the slender projectile used in the present investigation is approximately 53 0  D which, based on Eq. ( 101), results in mmax0 = (2100, 700) for to λ = (30°, 60°).The comparison of the maximum fragment masses corresponding to the various striking velocities immediately after the damage-fragmentation transition for these two crystallographic orientations of the flat-end monocrystalline projectile is presented in Table 1.The maximum fragment mass, mmax0, for seven striking velocities (v v0).The 2D-MD results are obtained for two symmetric crystallographic orientations of the flat-end monocrystalline projectile, λ = (60°, 30°).Note that striking velocities 0.357 m/s and 0.470 m/s correspond to the damage-fragmentation transition, respectively (i.e., N/A designates the values below the critical striking velocity corresponding to the fragmentation threshold, v < v0).
The MD simulation results presented in Table 1 agree rather well with the rough estimates of Eq. ( 10).All values following immediately the critical striking velocity-v0 = (0.357 km/s, 0.470 km/s) corresponding respectively to λ = (60°, 30°)-belong to the region close to the upper horizontal asymptote defining the baseline fragmentation response.With reference to the data presented in Table 1, it cannot be overemphasized that mmax0 is, due to the shock-induced thermal attrition, expected to be somewhat below the rough geometrical estimate of Eq. ( 10).Nonetheless, based on the simulation data, this fragment mass reduction due to the thermal equilibration is well within the inherent randomness of the initial fragmentation evident from Table 1 (compare, as an example, the striking velocities 0.357 km/s and 0.360 km/s.).Stochasticity of the maximum fragment masses in the phase transition region is expected based on the well-known aleatory variability inherent to low-energy fracture events (e.g., Mastilovic, 2011).The ratio of mmax0 values for the two orientations is-under the circumstances-reasonably close to the value Ctg 60°/Ctg 30° = 3 predicted by Eq. ( 101).Also, the second largest fragment masses are typically close to the mmax0 values in agreement with the symmetric-failure assumption of the fragmentation scenario of Figure 4.
Furthermore, the ratio of the impact energies corresponding to the damage-fragmentation transition for the two single-crystal orientations scales, fortuitously or not, with the ratio of the inclined failure lengths indicating that the failure energy necessary to brake-off the maximum fragment is dominated by the inclined-separation term.
The extreme orientation dependence of the perfect monocrystalline projectiles from the standpoint of both the critical velocity (v0) and the maximum fragment mass (mmax0), is a consequence of the increased ductility of single crystals in the case of more favorably oriented close-packed directions.
Finally, an extension of the present concept to an estimate of mmax0 for 3D Taylor impact would involve more complex analysis since the projectile fracture under multiaxial dynamic loading is a challenging problem.Following severe distortion, the impact surface of the Taylor projectile (made of ductile material) eventually petals (e.g., Meyers, 1994).The details of this sunflower-like petaling define mmax0 (under highly stochastic circumstances of the critical-impact-energy neighborhood).Based on the present 2D-MD simulation observations (Mastilovic, 2015a(Mastilovic, , 2016b)), it is justi-Latín American Journal of Solids and Structures 14 (2017) 1529-1546 fied to assume that the onset of the damage-fragmentation transition (Ξ0, mmax0) would correspond to the minimum impact energy sufficient to cause the non-negligible fragmentation.Thus, the extension is straightforward but the complexities involved would probably be most efficiently addressed by experimental determination of the critical point coordinates.

Validation of the MD Simulation Model
The 2D simulation model used in the present investigation is based on the time-honored MD techniques and utilizes the most-frequently used empirical potential (LJ 6-12) to explore some general trends of the small-scale impact fragmentation of a slender projectile at ultrahigh striking velocities.The motivation for the choice of the simulation technique is to study inhomogeneous fragmentation with no assumptions made about underlying processes and mechanisms.This reflects favorably on the model validation process since the first requirement-the validation of model assumptions-is satisfied by definition.It cannot be overstated that this robust MD model requires only three material input parameters (the atomic mass and diameter, and the strength of attraction); obviously, only a limited quantitative agreement with experimental data could be expected under such circumstances.The main limitations of the present model have been already discussed in the preceding article (Mastilovic, 2016b).
One of the most advertised advantages of the computational models of discontinua in general is the ability to push exploration of physical phenomena beyond current experimental limits.Unfortunately, this advantage inherently implies scarcity of data available for the model validation, which is frequently the reason why attempts to validate such models leave necessarily something to be desired (Sargent, 2011).Notably, the present 2D-MD model is developed and used intensively over the last tree years and it has been verified that the conceptual ideas have been correctly implemented to the extent that the program modules perform as expected.Thus, from the verification standpoint the simulation output is deemed acceptable for the investigation purpose.
Unfortunately, the model validation is adversely affected by lack of the Taylor test data consistent with the simulation results reported throughout this study.Namely, the ballistic Taylor test is the most commonly used direct impact experiment originally developed as a method of estimating the dynamic compressive strength of ductile materials.More recently, it has been used to verify material constitutive relations by comparing numerical predictions with experimental data (Field, 2004).Notably, the reported Taylor-test experimental data do not include either fragment distributions or pressure or temperature fields which preclude the direct comparison with the simulation results.The following model validation is, therefore, out of necessity, performed-to a certain extent-in a roundabout way.
First, the abovementioned LJ model parameters (m0 = 3.1×10 -25 kg, r0 = 2.8 Å, ε = 7.5×10 -20 J) are selected to match as close as possible elastic constants of tungsten (74W).Consequently, the value of the modulus of elasticity is captured in the MD simulation with excellent accuracy (within 1%).On the other hand, the velocity of longitudinal wave propagation, estimated by the present MD model to be C0MD = 4.35 km/s (Mastilovic, 2016a), is 8% higher than the reported experimental value for tungsten, C0exp = 4.03 km/s (e.g., Hixson and Fritz, 1992;Mayers, 1994) Mayers, 1994).Note that the accuracy of the modulus of elasticity estimate can be indirectly verified by rearranging the well-known expression for the velocity of longitudinal wave propagation, E MD /E exp ρ MD /ρ exp C0MD /C0exp 2 0.99.Further, the experimentally-verifiable simulation results include the fragmentation onset velocity (vfo), whose knowledge is of great practical importance in ballistics for the design of both the target and the projectile.Within the framework of the present simulation technique, vfo is determined to be roughly in the striking velocity range 2 km/s < v < 3 km/s (Mastilovic, 2015a), which agrees well with observations reported by Livingstone et al. (2001).Note that the above MD estimate could be further refined with, practically, arbitrary precision by performing iterative computations at additional striking velocities, which was not necessary for the objective of the study.
Next, the present MD model applied to low-velocity rigid-anvil simulations by Mastilovic and Krajcinovic (1999) reproduced the Taylor's experimental observations (1948) showing that the relative shortening of the slender projectile (L1/L0) is independent of the slenderness ratio (L0/D0).Also, these MD simulation results are in agreement with the classic analysis, originated by Taylor (and refined by Wilkins and Guinan; e.g., Meyers, 1994), which suggested the scaling relation between the relative projectile shortening and the impact energy , L1/L0 ∝ exp (K).
Due to the scarcity of quantitative Taylor test data available for direct comparison, it is necessary to rely on semi-quantitative and qualitative observations for the purpose of model validation.For example, the 2D-MD simulation results (Mastilovic, 2015a(Mastilovic, , 2016a(Mastilovic, , 2016b) ) confirm that fragment mass distributions are generally of the Poisson hyper-exponential type as experimentally verified for ductile materials (e.g., Mayers, 1994).As far as the strain rate effect on the fragment mass is concerned, the MD simulation result presented by Eq. ( 1), as already mentioned, agree well with the classic theoretical estimates.Succinctly, the Grady's classic model (1982) of dynamic fragmentation of ductile materials during shock compression determines the shear band spacing to be inversely proportional to the strain rate.Since the mean fragment mass, for 2D fragmentation, is proportional to the square of the shear band spacing, it follows that the mean fragment mass is expected to be, approximately, inversely proportional to the square of strain rate, which is a scaling relation of the same form as (1) (Mastilovic, 2015a).Moreover, the deformed projectile configurations (such as those of Figure 5) are in qualitative agreement with experimental observations of large number of void nucleation, growth and aggregation near the impact interface of the ductile-material projectile (e.g., Grady and Kipp, 1989;Woodward et al. 1992).
Finally, MD simulation results for various striking velocities are used to perform a detailed program of statistical-hypothesis testing to verify that the shock-induced vibrational velocities, used in the temperature calculation, belong to the Maxwell-Boltzmann distribution, which defines the local thermal equilibrium (Mastilovic, 2015b).The null hypothesis that the vibrational velocity is distributed according to the Maxwell-Boltzmann distribution is rejected at the 5 percent level (based on the Cramér-von Mises test) if the parameter presented in areas selected along the centerline of the projectile (for more details refer to Mastilovic, 2015aMastilovic, , 2015b)).Admittedly, the abovementioned examples of experimentally verifiable simulation results fall short of the verification objective since they do not include quantitative estimates of thermodynamic state parameters.Thus, an attempt is made to validate the conversion of simulation data generated at the nanoscale level (atomic positions and velocities, and interatomic forces) to macroscopic observables (stress and temperature) reported herein by devising an additional MD simulation model, which allows direct comparison with experimental observations for tungsten (Hixson and Fritz, 1992).This model mimics a plate impact by simulating a laterally-confined Taylor test (CTT) with the entire length of the slender projectile loosely fitted into a rigid hole at the moment of collision.
In the CTT configuration the projectile diameter matches that of the hole (i.e., it is ideallyloosely fitted with minimal clearance) and the projectile motion is friction-free along the lateral boundary mimicking a perfect lubrication.(It is illustrative to dubbed it "the worm in the hole" configuration, inspired by the de Gennes' memorable "the ant in the labyrinth" of statistical mechanics.)Thus, six additional CTT simulations are performed for verification purpose and the mean stress (the first stress invariant) values are compared with the corresponding shock pressures in Table 3.The experimental and computational peak pressure values are in a very good agreement, especially keeping in mind that the rigid-anvil simulations are two-dimensional.Unfortunately, the velocity range of the experimental study is relatively narrow, up ∈ (2.07, 3.91) km/s.(Note that the more recent laser-based experiments resulting in terapascal pressure levels (e.g., Crowhurst et al, 2011) are not available for tungsten.)Hixon and Fritz (1992) for tungsten with a two-stage light-gas gun with the CTT simulation results.(up is the initial particle velocity.)Note that lateral motion of the lateral boundary atoms is constrained by rigid walls in the CTT configuration to mimic the plate impact.
The verified pressure agreement is important not only in itself but also since it provides a means necessary to verify (alas, roughly and indirectly) the temperature simulation output.Namely, while the pressure is commonly reported in shock experiments the corresponding temperature data is not available.Thus, it should be noted that the linearity between pressure and temperature (implied by Figures 1 and 3) is in agreement with the classic ideal gas law, which is indicative keeping in mind the LJ-solid phase transition accompanying ultrahigh-velocity impact (Mastilovic, 2016a).

CONCLUSIONS
The simple 2D-MD model of the Taylor experiment aims to provide a phenomenological picture depicting influence of various ballistic state variables on the maximum fragment mass resulting from the inhomogeneous fragmentation.The variables include the impact energy and the maximum values of pressure, temperature and (square of) effective strain.The nonlinear and saturable character of the revealed functional dependences is described by a single reverse-sigmoid empirical relation that covers the complete fragmentation range from the damage-fragmentation transition (v = v0) to the shattering transition (v = v1).The simulation results suggest that, within the low and intermediate hypervelocity range (vfo < v < vL), the rigid-anvil impact is characterized by the maximum fragment mass approximately inversely proportional to the aforementioned state variables.For the present simulation setup, this linearity breaks down at the striking velocity vL  30 km/s.The unified phenomenological model offered in this article clearly outlines the limit of validity of the previously formulated piecewise-linear approximation by capturing the divergence from linearity at the high-energy end, corresponding to the asymptotic approach to the terminal (shattering) fragmentation.An estimate of the maximum fragment mass corresponding to the damage-fragmentation transition, mmax0, is offered based on geometrical considerations.Since the simulation data suggest a common slope (ξ ≈ 1) for all four state variables identified in this article (K, Pmax, Tmax, 2 max  ), once mmax0 is determined, this nonlinear model requires only one fitting constant per variable to capture functional dependence of the maximum fragment mass within the entire fragmentation range.The critical-point location (e.g., (v0, mmax0)) is, obviously, expected to be sensitive to both the projectile geometry and the textural features on the spatial scale that dominates the dynamic response of the particular ductile material (e.g., the grain morphology, various inhomogeneities, cavities...).On the Latín American Journal of Solids and Structures 14 (2017) 1529-1546 other hand, the extremely high strain rates, inherent to the hypervelocity impact, render the material response increasingly insensitive to the subtle microstructural features and the phenomenological observations obtained from the nanoscale MD simulations may become increasingly indicative of the macroscale response.

Figure 1 :
Figure 1: Logarithmic plot of the maximum fragment mass dependence upon the striking velocity of the projectile in the entire impact range (a) and the maximum values of the three selected state variables: the (square of the) effective strain rate (b), the pressure (c), and the instantaneous kinetic temperature (d) in the hypervelocity impact range.(The data point depicted by the solid square corresponds to v = 3 km/s, which is slightly above the hypervelocity impact threshold, vfo; while the data point at the extreme right corresponds to v = 60 km/s.Fortuitously or not, v× ≈ 1.2 km/s(Mastilovic, 2015a) corresponds to the balance between the initial impact energy and the potential energy per atom.).

Figure 2 :
Figure 2: Schematic plot of the maximum fragment mass vs. the generic state variable Ξ  {K, Pmax, Tmax, 2 max  } in the logarithmic space where the scaling parameter in the approximately linear intermediate hypervelocity range is ξ  1.This model disregards the minuscule fragmentation below the onset of the damage-fragmentation transition, Ξ < Ξ0 (Fig. 1a).

Figure 3 :
Figure 3: Logarithmic plot of the maximum fragment mass vs. the state variable Ξ  {K, Pmax, Tmax} shifted by

Figure 4 :
Figure 4: Schematic representation of the proximal region of the symmetric impact of a slender flat-nose bar on the rigid boundary.The bar diameter is designated by D0 and the crystallographic orientation of the monocrystalline projectile by λ (the angle of the symmetric directions of the highest linear atomic density with respect to the impact direction).The blue dotted contours depict the advancing compression wave (c) enclosing the trailing region of the compressed material, the red solid contours depict the tension (t) field generated by lateral release waves bouncing off each other at the symmetry line, and the dashed black contours depict the resulting damage (d) field.

Figure 5 :
Figure 5: Snapshots of two impact configurations corresponding to: (a) λ = ±60° (v = 357 m/s), (b) λ = ±30° (v = 480 m/s).The void alignments in the designated directions imply the typical ductile distortion and fracture propagation along inclined separation lines (l60°, l30°) roughly enveloping tension/damage fields.Another feature typical of the ductile dynamic fracture is the evident independent void formation at many sites.
. The discrepancy between C0 values is related to the inherent inability of the triangular MD lattice to model completely accurately any crystalline plane of the BCC lattice.Consequently, the geometry Latin American Journal of Solids and Structures 14 (2017) 1529-1546 of the MD hexagonal unit cell results in mass density experimental value ρexp = 19250 kg/m 3 (

Figure 6 :
Figure 6: Histogram and corresponding probability density function of the Maxwell-Boltzmann distribution obtained by the simulation data fit of atomic vibratory velocities in the averaging area M3 for v = 7 km/s at t = 17 ps.

American Journal of Solids and Structures 14 (2017) 1529-1546
Table 2 is less than 0.05.A statistical testing results indicate that the Maxwell-Boltzmann distribution is established well before the steady state is reached.

Table 2 :
Results of the Cramér-von Mises statistical hypothesis testing with null hypothesis being that the vibrational velocity is distributed according to the Maxwell-Boltzmann distribution.The time instances presented correspond to v = 7 km/s.The designators M1, M2 and M3 mark three averaging (measurement)

Table 3 :
Comparison of the experimental shock pressures obtained by