Phenomenology of the Maximum Fragment Mass Dependence Upon Ballistic Impact Parameters

Abstract

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.

Keywords:
Taylor test; impact fragmentation; maximum fragment mass; hypervelocity impact

1 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 = v ₀) and the shattering transition (v = v ₁), 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 dynamic-fragmentation 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 and Hietala, 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., 2010, 2012; Paluszny et al., 2014), finite element methods (Levy and Molinari, 2010; Ugrcic, 2013), and meshfree methods (Wu et al., 2014; Li et al., 2012, 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 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 (m _max ≡ 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.

2 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, 2016a, 2016b), thus, a succinct summary is deemed sufficient herein. A monatomic system is comprised of atoms of equal masses _{^{mi =}} m ₀ 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 (₇₄W) are the atomic mass m ₀ = 3.1×10^-25 kg (183.85 u), the atomic radius 1.4 Å (≡ r ₀/2 where r ₀ 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 slender LJ projectile (D ₀×L ₀ = 15×110 nm, m _p = 7.0×10^-21 kg), prepared at zero temperature, impacts a rough rigid target represented by a set of immovable atoms. The conversion of simulation data generated at the nanoscale level (atomic positions and velocities, and interatomic forces) to macroscopic observables (temperature, stress and strain) is firmly established nowadays (e.g., Wagner et al., 1992; Allen and Tildesley, 1996; Zhou, 2003; Buehler, 2008; Mastilovic, 2016a).

The link between two atoms ruptures when their mutual distance exceeds a predetermined critical value. The cut-off interatomic distance, R ≈ 1.7 r ₀, 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, 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 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.¹ 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).

3 OBSERVATIONS AND DISCUSSION

The knowledge of the maximum fragment mass dependence upon the impact energy (the initial kinetic energy) of the projectile (K ∝v ²) 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 ( v _fo),² 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 v _L ( 30 km/s. Based on this approximation, the terminal fragmentation (m _max ≡ 1) is anticipated to be v = v ₁ = 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.

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, v _fo; 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.).

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, ( ² _eff = (( ² _x = ( ² _y ) / 2. The evolutions of these three 2D state variables are recorded at twelve evaluation areas mimicking the measurement gages (details are available in Mastilovic, 2015a, 2015b). The first mechanical stress invariant, P, is a measure of the force interaction between 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 (P_max) 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.³

As an example, the scaling relation of the form

$m_{\max }\propto {\left({\dot{\epsilon }}_{\max }^2\right)}^{-1}$ (1)

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 (v _fo < v < v _L). 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 v _L ≈ 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

$m_{\max }\propto {\Xi }^{-\xi }$ (2)

captures the elucidated linearity observations within the lower and intermediate part of the hypervelocity impact range, with ξ ≈ 1 for the generic state variable Ξ ∊{K, P_max, T_max, ${\dot{\epsilon }}_{\max }^2$}. According to the simulation results presented in Figure 1, the upper bound for validity of the scaling relation (2) is v _L ≈ 30 km/s.

The dependence of the maximum fragment mass upon the aforementioned state variable Ξ for the entire (non-negligible, Ξ ≥ Ξ₀) fragmentation range explored in the present study is schematically depicted in the logarithmic space by the following expression

$\ln m_{\max }=\ln m_{\max 0}\cdot \exp \left[{\left(\frac{\ln {\Xi }_0-A}{B}\right)}^C-{\left(\frac{\ln \Xi -A}{B}\right)}^C\right],\left(\Xi \ge {\Xi }_0\right)$ (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, Ξ < Ξ₀ (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 (Ξ₀, m _max0) 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 v ₀ in terms of the system size and determined critical velocities of the infinite system.

Figure 2:
Schematic plot of the maximum fragment mass vs. the generic state variable Ξ ∊ {K, P_max, T_max, ${\dot{\epsilon }}_{\max }^2$} 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, Ξ < Ξ₀ (Fig. 1a).

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

$\xi =\ln m_{\max 0}\cdot \frac{C}{B}\cdot \exp \left[-1+\frac{1}{C}+{\left(\frac{\ln {\Xi }_0-A}{B}\right)}^C\right],\ln {\Xi }_{\xi }=A+B\cdot {\left(1-\frac{1}{C}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$C$}\right.}$ (4)

The inflection point (4₂) 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

$\ln m_{\max }=M_0\cdot \exp \left[-{\left(\frac{\ln \Xi -A}{B}\right)}^C\right],M_0=\frac{B}{C}\cdot \exp \left(1-C^{-1}\right)\cdot \xi$ (5)

can be derived straightforwardly by combining (2) and (4).

The empirical expressions (3) and (5) cannot formally capture the terminal fragmentation, m _max ≡ 1, except as the limit case:

$\underset{\Xi \to \infty }{\lim }\left(m_{\max }\right)=1$ (6)

defining the horizontal asymptote illustrated in Figure 2. As previously mentioned, the lower bound of the shattering-transition threshold is extrapolated to v ₁ ≈ 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; m _max = 1 ≡ m _max1).

The additional MD simulations performed herein at three ultrahigh striking velocities v = (45, 50, 60) km/s reach the fragmentation steady state with m _max = (3, 2, 2) > m _max1 ≡ 1, respectively. Thus, the elusive terminal fragmentation of 15×110-nm ₇₄W projectile is not achieved even after 900,000 time steps (90 ps). The nonlinear data fit illustrated in Figure 1a points roughly to v ₁ > 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, P_max, T_max, ${\dot{\epsilon }}_{\max }^2$} but also that

$B\approx \ln m_{\max 0},C\approx 2\xi$ (7)

Thus, (5) can be rewritten in a form

$\ln m_{\max }=\ln m_{\max 0}\cdot \exp \left[D-{\left(\frac{\ln \Xi -A}{\ln m_{\max 0}}\right)}^{2\xi }\right],D=1+\ln \frac{1}{2}-\frac{1}{2\xi }$ (8)

(Note that based on simulation results presented in Figure 1, D is estimated to be 0.20 for Ξ and 0.45 for v ∝ $\sqrt{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.

Figure 3:
Logarithmic plot of the maximum fragment mass vs. the state variable Ξ ∊{K, P_max, T_max} shifted by A along the abscissa in the hypervelocity impact range (where A ∊{-35.9, 20.3, 5.0} in respective SI units).

3.1 A geometric estimate of the orientation dependence of m _max0 in perfect monocrystals

The maximum fragment mass corresponding to the fragmentation threshold (m _max0) 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 (Ξ₀, m _max0), considering the saturable character of the proposed phenomenological dependence and the uniquely-defined: (i) the lower horizontal asymptote (ln m _max1 = 0) and (ii) the linear intermediate hypervelocity range with slope ξ ≈ 1 (Figure 3).

In the present 2D-MD framework, it is intuitively clear that m _max0 (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).

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 D ₀ 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.(The relative thickness of contour lines reflect qualitatively the field intensities.).

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 (l _60°, l _30°) roughly enveloping tension/damage fields. Another feature typical of the ductile dynamic fracture is the evident independent void formation at many sites.

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 m _max0 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).

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

$V_{\max 0}=\frac{3}{8}{D}_0^{2}Ctg\lambda$ (9)

Since the area per atom in the perfect triangular lattice is $\sqrt{3}{r}_0^2/2$, the maximum fragment mass of the fragmentation threshold is

$m_{\max 0}=\frac{\sqrt{3}}{4}{\overline{D}}_0^{2}Ctg\lambda ,{\overline{m}}_{\max 0}=\frac{3Ctg\lambda }{8\left(L_0/D_0\right)}$ (10)

where ${\overline{D}}_0=D_0/r_0$, m _max0 is given in number of atoms (m ₀) constituting the fragment, ${\overline{m}}_{\max 0}$ is the maximum fragment mass normalized by the initial projectile mass, and L ₀/D ₀ stands for the projectile slenderness ratio. It is important to recognize that Eqs. (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 ${\overline{D}}_0=53$ which, based on Eq. (10₁), results in m _max0 = (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.

Table1:
The maximum fragment mass, m _max0, for seven striking velocities (v ≥ v ₀). 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 < v ₀).

λ	v [km/s]
	0.357	0.360	0.374	0.408	0.442	0.470	0.490
60°	717	682	620	729	628	623	815
30°	N/A	N/A	N/A	N/A	N/A	1798	2060

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-v ₀ = (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 m _max0 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 m _max0 values for the two orientations is-under the circumstances-reasonably close to the value Ctg 60°/Ctg 30° = 3 predicted by Eq. (10₁). Also, the second largest fragment masses are typically close to the m _max0 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

${\left(\frac{{\left(v_0\right)}_{{30}^{\circ }}}{{\left(v_0\right)}_{{60}^{\circ }}}\right)}^2\propto \left(\frac{l_{{30}^{\circ }}}{l_{{60}^{\circ }}}\right)=\sqrt{3}$ (11)

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 (v ₀) and the maximum fragment mass (m _max0), 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 m _max0 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 m _max0 (under highly stochastic circumstances of the critical-impact-energy neighborhood). Based on the present 2D-MD simulation observations (Mastilovic, 2015a, 2016b), it is justified to assume that the onset of the damage-fragmentation transition (Ξ₀, m _max0) 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.

3.2 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 (m ₀ = 3.1×10^-25 kg, r ₀ = 2.8 Å, ε = 7.5×10^-20 J) are selected to match as close as possible elastic constants of tungsten (₇₄W). 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 C _{0 MD} = 4.35 km/s (Mastilovic, 2016a), is 8% higher than the reported experimental value for tungsten, C _{0 exp} = 4.03 km/s (e.g., Hixson and Fritz, 1992; Mayers, 1994). The discrepancy between C ₀ 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 of the MD hexagonal unit cell results in mass density ${\rho }_{MD}=\left(2\sqrt{3}/3\right)\left(m_0/r_0^3\right)=16300kg/m^3$ that by 15% underestimates the experimental value ϱ_exp = 19250 ^_kg/m3 (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, ^_EMD /^_Eexp = (^_ρMD /^_ρexp )( C _{0 MD} /C _{0 exp} )² = 0.99.

Further, the experimentally-verifiable simulation results include the fragmentation onset velocity (v _fo), 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, v _fo 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, 2016a, 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 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.

Table2:
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 M₁, M₂ and M₃ mark three averaging (measurement) areas selected along the centerline of the projectile (for more details refer to Mastilovic, 2015a, 2015b).

t [ps]	M1	M2	M3
6	0.000	0.116	0.001
11	0.165	0.805	0.119
13	0.660	0.394	0.210
17	0.377	0.198	0.955

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 M₃ for v = 7 km/s at t = 17 ps.

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 ideally-loosely 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.)

Table3:
Comparison of the experimental shock pressures obtained by Hixon and Fritz (1992) for tungsten with a two-stage light-gas gun with the CTT simulation results. (u _p 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.

up [km/s]	Pmax [GPa]		Relative Difference [%]
	simulation	experiment
2.07	258.	263.88	2.2
2.34	316.	315.87	0.0
2.82	422.	408.33	3.3
3.19	525.	499.50	5.1
3.60	645.	603.44	6.9
3.91	742.	675.93	9.7

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).

4 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 = v ₀) to the shattering transition (v = v ₁). The simulation results suggest that, within the low and intermediate hypervelocity range (v _fo < v < v _L), 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 v _L ≈ 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, m _max0, 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, P_max, T_max, ${\dot{\epsilon }}_{\max }^2$), once m _max0 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., (v ₀, m _max0)) 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 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.

References

• Allen, M.P., Tildesley, D.J. (1996). Computer Simulation of Liquids. Oxford University Press.
• Alves, M.,Yu, J.L., (2005). Material influence on the response of impacted beams. Latin American Journal of Solids and Structures 2: 167-193.
• Astrom, J.A., Holian, B.L., Timonen, J., (2000). Universality in fragmentation. Physical Review Letters 84(14): 3061-3064.
• Baker, K.L., Warner, D.H., (2012). Simulating dynamic fragmentation processes with particles and elements. Engineering Fracture Mechanics 84: 96-110.
• Buehler, M.J. (2008). Atomistic Modeling of Materials Failure. Springer.
• Crowhurst, J.C., Armstrong, M.R., Knight, K.B., Zaug, J.M., Behymer, E.M. (2011). Invariance of the Dissipative Action at Ultrahigh Strain Rates Above the Strong Shock Threshold. Physical Review Letters, 107: 144302.
• Diehl, A., Carmona, H.A., Araripe, L.E., Andrade Jr, J.S., Farias, G.A., (2000). Scaling behavior in explosive fragmentation. Physical Review E 62(4): 4742-4746.
• Dos Santos, F.P.M., Barbosa, V.C., Donangelo, R., Souza, S.R., (2011). Experimental analysis of lateral impact on planar brittle material. Physical Review E 84: 026115.
• Elek, P., Jaramaz, S., (2009). Fragment Mass Distribution of Naturally Fragmenting Warheads. FME Transactions 37: 129.
• Field, J.E., Walley, S.M., Proud, W.G., Goldrein, H.T., Siviour, C.R., (2004). Review of experimental techniques for high rate deformation and shock Studies. International Journal of Impact Engineering 30, 725-775.
• Grady, D.E., Kipp, M.E., (1989). Fragmentation of solids under dynamic loading, in: “Structural Failure”, eds. Wierzbicki, T. and Jones, N., pp. 1-40. John Wiley & Sons, New York.
• Grady, D.E. (2006). Fragmentation of Rings Shells. Springer.
• Grady, D.E., (1982). Local inertial effects in dynamic fragmentation. Journal of Applied Physics 53(1): 322-325.
• Grady, D.E., Kipp, M.E., (1993). Dynamic Fracture and Fragmentation, in: “High-Pressure Shock Compression of Solids”, eds. Asay, J.R. and Shahinpoor, M., pp. 265-322, Springer, Berlin.
• Greer, J.R., De Hosson, J. Th. M., (2011). Plasticity in small-sized metallic systems: Intrinsic versus extrinsic size effect. Progress in Material Science, 56: 654-724.
• He, A-M., Wang, P., Shao, J-L., (2015). Molecular dynamics simulations of ejecta size distributions for shock-loaded Cu with a wedged surface groove. Computational Material Science 98: 271-277.
• Hixson, R.S., Fritz, J.N., (1992). Shock compression of tungsten and molydbenum. Journal of Applied Physics 71(4): 1721-1728.
• Holian, B.L., Grady, D.E., (1988). Fragmentation by molecular dynamics: the microscopic “big bang”. Physical Review Letters 60: 1355-1358.
• Iturrioz, I., Fleck, L., Miguel, F., Riera, J.D., (2009). Dynamic fracture analysis of concrete or rock plates by means of the Discrete Element Method. Latin American Journal of Solids and Structures 6: 229 - 245
• Kadono, T., (1997). Fragment Mass Distribution of Platelike Objects. Physical Review Letters 78(8): 1444.
• Kadono, T., Arakawa, M., (2002). Crack propagation in thin glass plates caused by high velocity impact. Physical Review E 65: 035107(R).
• Katsuragi, H., Sugino, D., Honjo, H., (2004). Crossover of weighted mean fragment mass scaling in two-dimensional brittle fragmentation. Physical Review E 70: 065103(R).
• Kraft, O., Gruber, P.A., Monig, R., Weygand, D., (2010). Plasticity in confined dimensions. Annual Review of Material Research, 40: 293-317.
• Kumar, V., Ghosh, A., (2015). Non-linear dynamic fragmentation using Cracking Particles Method. Computational Material Science 98:117-122.
• Kun, F., Herrmann, H.J., (1999). Transition from damage to fragmentation in collision of solids. Physical Review E 59 (3): 2623-2632.
• Lamberson, L., Eliasson, V., Rosakis, A.J., (2012). In Situ Optical Investigations of Hypervelocity Impact Induced Dynamic Fracture. Experimental Mechanics 52 (2): 161-170.
• Levy, S. and Molinari, J.F., 2010. Dynamic fragmentation of ceramics, signature of defect s and scaling of fragment sizes. Journal Mechanics and Physics of Solids 58: 12 -26.
• Li, B., Kidane, A., Ravichandran, G., Ortiz, M., (2012). Verification and validation of the Optimal Transportation Meshfree (OTM) simulation of terminal ballistics. International Journal of Impact Engineering 42: 25-36.
• Li, B., Pandolfi, A., Ravichandran, G., Ortiz, M., (2015). Material-point erosion simulation of dynamic fragmentation of metals. Mechanics of Materials 80: 288-297.
• Livingstone, I.H.G., Verolme, K., Hayhurst, C.J., (2001). Predicting the fragmentation onset velocity for different metallic projectiles using numerical simulations. International Journal of Impact Engineering 26: 453-464.
• Mastilovic, S., (2011). Some observations regarding stochasticity of dynamic response of 2D disordered brittle lattices. International Journal of Damage Mechanics 20: 267-277.
• Mastilovic, S., (2015a). Impact fragmentation of nanoscale projectiles at ultrahigh striking velocities. Meccanica, 50: 2353-2367.
• Mastilovic, S., (2015b). Ballistic impact fragmentation at high striking velocities: Part II: Maximum fragment mass. Proceedings of the 5^th International Congress of Serbian Society of Mechanics (Eds. Spasic, D.T., Lazarevic, M., Grahovac, N., Zigic, M.) 978-86-7892-715-7, Arandjelovac, Serbia.
• Mastilovic, S., (2016a). Molecular-dynamics simulations of the nanoscale Taylor test under extreme loading conditions. Mathematics and Mechanics of Solids 21(3): 326-338.
• Mastilovic, S., (2016b). Damage-fragmentation transition: Size effect and scaling behavior for impact fragmentation of slender projectiles. International Journal of Damage Mechanics, DOI: 10.1177/1056789516671775.
  » https://doi.org/10.1177/1056789516671775
• Mastilovic, S., Krajcinovic, D., (1999). Statistical models of brittle deformation: Part II: computer simulations, International Journal of Plasticity 15, 427-456.
• Meyers, M.A., (1994). Dynamic Behavior of Materials. John Wiley & Sons.
• Paluszny, A., Tang, X.H., Zimmerman, R.W., (2014). Fracture and impulse based finite-discrete element modeling of fragmentation. Computational Mechanics 52(5): 1071-1084.
• Ramesh, K.T., Hogan, J.D., Kimberley, J., Stickle, A., (2015). A review of mechanisms and models for dynamic failure, strength, and fragmentation. Planetary and Space Science 107: 10-23.
• Rinaldi, A., (2011). Effects of dislocation density and sample size on plastic yielding at the nanoscale: a Weibull-like framework. Nanoscale 3(1): 4817-4823.
• Sargent, R.G., (2011). Verification and validation of simulation models. Proceedings of the 2011 Winter Simulation Conference (Eds. Jain, S., Creasey, R.R., Himmelspach, J., White, K.P., and Fu, M.C.) 978-1-4577-2109-0/11 Phoenix, AZ, USA.
• Sator, N., Hietala, H., (2010). Damage in impact fragmentation. International Journal of Fracture 163: 101.
• Sator, N., Mechkov, S., Sausset, F. (2008) Generic behaviors in impact fragmentation. Europhysics Letters 81: 4402.
• Taylor, G.E., (1948). The use of flat-ended projectiles for determining dynamic yield stress I. Theoretical considerations. Proceedings of the Royal Society of London A 194 (1038): 289-299.
• Timár, G., Blömer, J., Kun, F., Herrman, H.J., (2010). New Universality Class for the Fragmentation of Plastic Materials. Physical Review Letters 104: 095502.
• Timár, G., Kun, F., Carmona, H.A., Herrman, H.J., (2012). Scaling laws for impact fragmentation of spherical solids. Physical Review E 86 (4): 016113.
• Ugrcic, M., (2013). Numerical Simulation of the Fragmentation Process of High Explosive Projectiles. Scientific Technical Review 63(2): 47-57.
• Wagner, N.J., Holian, B.L., Voter, A.F., (1992). Molecular-dynamics simulations od two-dimensional materials at high strain rates. Physical Review A 45(12): 8457-8470.
• Wittel, F.K., Carmona, H.A., Kun, F., Herrmann, H.J., (2008). Mechanisms in impact fragmentation. International Journal of Fracture 154:105-117.
• Woodward, R.L., O’Donnell, R.G., Flockhart, C.J., (1992). Failure mecanisms in impacting penetrators. Journal of Materials Science 27(23): 6411-6414.
• Wu, Y., Wang, D., Wu, C-T., (2014). Three dimensional fragmentation simulation of concrete structures with a nodally regularized meshfree method. Theoretical and Applied Fracture Mechanics 72: 89-99.
• Zhou, M., (2003). A new look at the atomic level virial stress: on continuum-molecular system equivalence. Proceedings of the Royal Society of London A 459: 2347-2392.

• 1
  Note that although the impact fragmentation is commonly a 3D phenomenon, experiments were performed with 2D experimental setup as well (Kadono, 1997; Kadono and Arakawa, 2002; Dos Santos et al., 2011).

• 2
  The fragmentation onset velocity, vfo, defined as the threshold velocity just sufficient to fully fragment the projectile, is determined within the present simulation framework to be roughly between 2 and 3 km/s (Mastilovic, 2015a).

• 3
  Additionally, it has been verified that the stress and temperature simulation results are objective with respect to the evaluation area size.

Publication Dates

• Publication in this collection
  16 Aug 2017

History

• Received
  05 May 2016
• Accepted
  14 June 2017