Quantification of the Kinetic Energy Conversion to Temperature Increase in Metal-on-Metal Impacts up to Hypervelocity Conditions by Molecular Dynamics Simulation

The dynamic impact loading of metals goes along with energy conversion from kinetic energy to internal energy and, ultimately, temperature increase. The fraction of the kinetic energy partitioned into heating is strongly dependent on the impact velocity. Limiting cases are already well characterized, both experimentally and numerically. At low velocities, plastic work is the main source of internal energy increase and contributes to approximately 100% to material heating. Toward high velocities, approaching a hydrodynamic-like condition but still below the threshold for material melting or vaporization, about 50% of the kinetic energy is converted to internal energy. The current work addresses the intermediate regime of mixed phenomenology, where analytical descriptions are hardly feasible and typical simulation methods of impact engineering, namely hydrocodes, fail to make reliable numerical predictions. For this purpose, we here alternatively apply molecular dynamics simulations at the nanometer scale, taking iron as exemplary test case. The results complement early findings by extending them to a broader range of validity.


Quantities
Material density V Volume c v Specific heat capacity at constant volume v Projectile velocity D Projectile diameter L Projectile length P Projectile depth of penetration into the target Q (q) Heat (per volume) E (e) Internal energy (per volume) S (s) Entropy (per volume) W pl (w pl ) Plastic work (per volume) K (k,k) Kinetic energy (per volume, per atom) k B Boltzmann's constant N

Number of atoms a Fe
Atomic mass of iron

Introduction
The topic of thermal effects upon high-strain rate deformation, i.e. for strain-rates ̇≳ 10 2 1/s, is crucial in several fields of study and technological applications covering a wide spectrum of dimensional scales. Among others, we mention examples for crashworthiness of structures [1], penetration mechanics for defense and aerospace applications [2][3][4], manufacturing processes [5], or planetary science [6]. From the perspective of numerical modeling, a correct computation of temperatures is of fundamental importance for reproducing the temperature-dependent behavior of materials, which defines part of their constitutive response.
While more or less sophisticated constitutive models for metals account for thermal scaling and evolution of the flow stresses [7,8], the partition mechanisms of mechanical work into temperature increase in impact events as well as their interplay at different velocities have yet to be fully understood due to limitations in both experimental measurements [9] and numerical computations [10] of such fast-transient events. In particular, the presence of some numerical limitations in hydrocode simulation, such as mass loss associated to element erosion [10,11], affects the energy computation and related partition measurements, compromising a full insight in the above mechanisms within the transition regime between high-velocity and hypervelocity impact, namely ranging from about 2 to 9 km/s depending to the material combination [10,12].
In the context above, molecular dynamics (MD) simulation represents an interesting option by eliminating artifacts of mass loss. Moreover, any influence on the results of the constitutive model and the equation of state would be bypassed by the choice of a proper interatomic potential governing the motion of atoms. This allows a direct correlation of the thermodynamics response to fundamental physical mechanisms at smaller scales. Recently, some works [13,14] have exploited MD to estimate the thermal conversion of plastic work in the dynamic regime at strain-rates associated to low-velocity impact conditions [10], namely in the domain still defined by the Taylor-Quinney coefficient [15,16].
The analysis of the energy conversion over a broader impact velocity regime via MD simulation connecting the already well characterized limiting cases [10] is the subject of the present work.

Basics of Energy Partitioning in Impact Events
In the lower part of the dynamic loading regime ( ̇≲ 10 3 1/s), the material heating is dominated by plastic work. Typically, only a fraction of the plastic work is converted to internal energy associated to temperature increase, reflecting that the typical process is adiabatic ( q = 0 ) but non-reversible: This inequality between work done and internal energy increase is traditionally quantified by the Taylor-Quinney (TQ) coefficient [15,16]: whose upper bound is 1, unless phenomena other than plastic deformation are involved, such as release of latent heat Δw pl , associated to phase transformation [17]. The complement 1 − relates to the increase of the configurational energy of the system, for example in terms of dislocation and twins [18,19]. At very high loading rates toward the high-velocity and hypervelocity impact regimes ( ̇≅ 10 4 1/s to 10 7 1/s), shockwave generation dominates and the volume-dependent part of the internal energy, still negligible for low-rate loading, also contributes to energy partitioning and must not be omitted: As discussed in a recent work by the authors [10] the transition in the energy partitioning mechanisms for metals across different impact velocity regimes is, therefore, generally composed of the plastic work contribution (deviatoric deformation) and that originating from the pressure and volume changes (equation of state) above mentioned. The velocity-dependent temperature increase can be synthetically defined as: where ′ describes the partition due to shock heating [12]; k 0 ∼ p v 2 0 is the initial projectile specific kinetic energy. Equation (1) is widely adopted under several-but equivalent-formulations in different hydrocodes [20,21]. The relationship is valid until incipient melting conditions for involved materials are met [12]. The value of the TQ coefficient is widely accepted in the literature as close to unity at impact velocities below 1 km/s from consolidated experimental measurements and numerical computations [13,14,[17][18][19][22][23][24][25] and this is reflected in its assumption as constant generally between 0.9 to 1 in hydrocode simulations [20,26]. The shock partition factor ′ , accounting for irreversible contribution [12], assumes values of about 0.5 to 0.6 within the transition regime between high-velocity and hypervelocity impact as confirmed by experimental data [4], analytical modeling thereof [12], and hydrocode simulations [10]. Note that the variation in thermal response at different velocities originates from the progressive (relative) decrease of the plastic work contribution Δw pl and parallel increase of the volumetric contribution ( ′ k 0 ). Therefore, the global kinetic energy partition factor ΔE T ∕K 0 , where ΔE T ≡ ∫ Vc v dT is the amount of the internal energy associated to temperature increase, cf. Eq. (3), is expected to move from about 1.0 (low-velocity impact) to 0.5 toward hypervelocity. The behavior in between these well-characterized limit cases is still to be confirmed by either experimental and numerical measurements.

Setup of the Molecular Dynamics Simulation Models
We here reproduce macroscopic penetration experiments at the nano-level by molecular dynamics simulation, attempting to eliminate in this way the possible errors in energy computation originated by mass loss from the system. In MD simulation [27], the interaction between particles (atoms) is given by an appropriately chosen potential. This bond relationship intrinsically incorporates the physical features of the material behavior that in a hydrocode [11] would be modeled by constitutive equations, a failure model, and an equation of state.
Simulations of a (bcc) iron projectile against an iron thick target were performed with the software package LAMMPS [28], replicating at the nanoscale the macroscopic semi-infinite penetration experiments of steel longrod projectiles ( L∕D = 10, cm-size) against steel targets of Ref. [29] and the corresponding hydrocode simulations of Ref. [10]. Such configuration maximizes the conversion of the impact kinetic energy into internal energy by complete deceleration of the projectile and avoiding the generation of target ejecta that would result from full perforation experiments. For this purpose, in this study pure iron is considered to be a good approximation for low-alloy steels typically used in ballistic applications [7,10], such as in the experiments of Ref. [29].
The modeled system consists of a 38.8 nm-thick target, whose boundary conditions are periodic in the in-plane directions. The lateral sizes of the cell are 38.8 × 38.8 nm 2 , cf. Figure S1 in the Supplementary Information, resulting in a total number of atoms in the target of about 5 millions ( N t ). The length of the cylindrical projectile is of 30.3 nm, with L∕D ≅ 10, resulting in N p = 18,851 atoms. Analogous simulations were performed employing a spherical projectile of nearly equivalent mass ( N p = 18,883, D = 7.5 nm) to assess possible shape effects in the energy partitioning. According to previous results in the literature [30], the system is sufficiently large to avoid finite-size effects and correctly predict the dimensional scaling of crater volume up to macroscale. The adopted pair-potential is the embedded atom method relationship by Mendeleev et al. [31,32], which has been extensively proven to accurately describe the behavior of iron at short interatomic distances, such as those originating from extreme pressures, and successfully exploited in shock compression and impact simulations [33][34][35].
The atoms of the system are assigned with random initial velocities consistent with a Gaussian distribution corresponding to a temperature T = 300 K. Thermal equilibration is performed in the isothermal-isobaric (NPT) ensemble toward a target system conditions of T eq = 300 K and zero pressure. The total equilibration time is 2 ps, with a timestep of 1 fs. In a second stage, an initial translational velocity v 0 is assigned to the projectile atoms and the impact simulation is performed in the microcanonical NVE ensamble. The analyzed impact velocities are in the range 0.5 to 3.5 km/s. The total simulation time is set to 50 ps in order to achieve the termination of the impact event, namely the full deceleration of the projectile with a pseudo-equilibration of the system (damping of residual translational kinetic energy). The timestep adopted for the impact stage is of 0.1 fs.

Estimation of the Energy Partition Factor
Snapshots of the projectile penetration at the end of simulation are reported in Fig. 1 for a selection of the adopted impact velocities. It must be noted that the normalized depth of penetration at this scale is much lower than the values expected for impact of similar materials at the macroscale, where, generally, P∕L ≅ √ p ∕ t for (long) rods and P∕L ≅ 2 √ p ∕ t for spheres at v 0 > 2 km/s [37]. However, this difference between MD and macroscopic results [29] is to be expected due to the preponderance of surface over volume forces at this size scale and results are consistent with observation in similar MD studies [30,35]. This is also reflected by the fact that the normalized penetration of the spherical projectile is lower than the rod, as 0.5 km/s 1.5 km/s 2.5 km/s 3.5 km/s 10 nm Fig. 1 Penetration of long-rod and spherical projectiles, top and bottom row respectively, at the end of the simulation ( t = 50 ps) for selected impact velocities. Particles belonging either to projectile or target are identified with different colors. Results are visualized using the software OVITO [36] opposite to macroscale, cf. Fig. 1. Moreover, for a given impact velocity the corresponding strain-rate, ̇∼ v 0 ∕L , is sensibly different with respect to Ref. [29] due to the several orders of magnitude of difference in the characteristic size scale of the systems. According to the equipartition theorem, the computation of the average temperature for a group of N atoms (in our case N = N t + N p ) is defined by their kinetic energies as: where k i are net of the center-of-mass (COM) translational component, and d = 3 is the dimensionality of the simulation. The number of atoms contributing to the temperature computation is constant during the duration of the run, as no mass is lost.
According to Eq. (5), the computation of the current value of the projectile kinetic energy partition into thermal effects can be expressed directly as [38]: where K 0 = 1 2 N p a Fe v 2 0 is the projectile initial impact kinetic energy. Equation (6) expresses the fact that the estimate of the internal energy associated to temperature increase can be reduced to a simple computation of the current values of atomic kinetic energies and the corresponding ones at the reference status. The resulting curves showing the evolution of system temperature and the energy partitioning for different impact velocities are depicted in Fig. 2.
By considering an appropriate time window t 1 < t < t 2 so that almost all of the projectile kinetic energy has been converted (see plots of the displacement and velocity of the Long rod Sphere Fig. 2 Evolution of the system temperature (Eq. (5), top graphs) and of the kinetic energy partition measure (Eq. (6), bottom graphs) at different impact velocities for the two penetrators. System temperature measurements are also reported in detail in Figure S2 in the Supplementary Information. Note that for the purpose of computation of energy partition only the region corresponding to about t ≳ 30 ps is meaningful, cf. the evolutions of the displacement and velocity of the projectile COM in Figures S3 and S4, respectively, in the Supplementary Information projectile COM in Figures S3 and S4, respectively, in the Supplementary Information) the partition factor is then calculated as time average: where t 2 = 50 ps, i.e. the simulation termination time, while t 1 was set so that the difference between the real temperature and the temperature accounting for the translational component of the kinetic energy [38] becomes lower than 10 -2 K, namely the impact kinetic energy of the projectile K 0 can be considered as fully converted and the penetration process terminated.
The results of the evaluation of the partition measure as for Eq. (7) are reported in Fig. 3. As important results, the performed MD simulations confirm that: (i) toward low-velocity impact ⟨ΔE T ∕K 0 ⟩ → 1 and the TQ coefficient framework describes the temperature increase; (ii) toward hypervelocity impact ⟨ΔE T ∕K 0 ⟩ → 0.5 and the partition measure is described by the relationships of the shock problem; (iii) for intermediate velocity ranges with a monotonous decrease of the partition factor with increasing impact velocity.
For the low velocity regime we note that the partition factor is approaching the TQ values computed for the dynamic regime either by several experimental works [17][18][19][22][23][24][25] and MD simulation studies [13,14], namely above about 0.9. Although generally ⟨ΔE T ∕K 0 ⟩ > in this regime, as mechanisms other than plastic work may contribute to temperature increase, the shock contribution is progressively becoming marginal toward the low-velocity regime. Therefore, it can be stated that toward lowvelocity ⟨ΔE T ∕K 0 ⟩ ≅ , namely that most of the impact kinetik energy is being converted to plastic work. The increasing variance of ⟨ΔE T ∕K 0 ⟩ for lower impact velocities, cf. Figs. 2 and 3, is merely a result of the employed method to compute the temperature increase. Indeed, the fact that N p ≪ N t , since the dimensions L and D of the projectiles must be sufficiently smaller with respect to the target size to avoid edge effects in the penetration, yield an increment of the temperature T N that at the lower velocities tends to be comparable to the temperature fluctuations of the system. Toward hypervelocity the energy partition computation is consistent with the experimental measurements by Schneider and Stilp [4] and related hydrocode simulation [10], as well as with our theoretical estimation of the factor ′ from shock-driven penetration mechanism assumptions [12]. Among the various modeling uncertainties, the slight difference with respect to the theoretically predicted value (0.5 vs. 0.6 at 3.5 km/s) could be explained by the fact that the theoretical derivation in Ref. [12] assumes the Hugoniot curve in place of the insentrope for computing the shock unloading path, thus yielding a slight overestimation of the energy retained by the material and converted into heat with respect to the actual behavior. It must also be noted that the progressive reduction of ΔE T ∕K 0 with velocity has to be attributed solely due to a transition in deformation mechanisms and not due to kinetic energy dissipation by phase transition (melting), which was not observed in our MD simulation for v 0 ≤ 3.5 km/s, in agreement with the expected melting thresholds [2,12].

Conclusions
From the results it clearly emerges how the assumption of either TQ approximation or the shock partition calculation for the whole transition regime may yield to a significant discrepancy in the temperature computation-up to a factor of 2 in excess or defect-with respect to the actual temperature increase. Interestingly, the consistency of the results between nano-and macro-scale [4,10] would suggest a multiscale and self-similar nature of such transition.
To conclude, we believe that this work gives a corroboration of the early results on energy partitioning at low-velocity and hypervelocity limit conditions as well as successfully quantifies, for the first time, the transition in the partitioning mechanism in the range of velocities in between these regimes. This estimate is an important ingredient for a more accurate thermo-mechanical understanding of impact events  7) for different impact velocities and projectile geometries for a broad range of applications [39][40][41] and for the derivation of related equation-of-state models.
Funding Open Access funding enabled and organized by Projekt DEAL. Research funding provided by Bundesministerium der Verteidigung (BMVg), Germany.

Conflict of interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work done in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.