Predicting shock hugoniot and equation of state of single-crystal molybdenum by molecular dynamics simulation

In the present paper, we employ non-equilibrium molecular dynamics approach to simulate impact-induced shock propagation in single-crystal Molybdenum. Shock hugoniot, generated by simulation of impact with varied strength shows excellent agreement with experimental data in the strong shock regime. The resulting hugoniot parameters obtained by linear fitting of shock velocity vs. particle velocity data are then used to estimate ambient pressure Grüneisen coefficient. Finally Mie-Grüneisen equation of state (EOS) with hugoniot as reference state, is employed to express pressure-volume-energy relationship for [001] single-crystal Mo. The influence of different analytical forms for volume dependence of Grüneisen parameter on EOS is investigated.


Introduction
Investigation of material properties at extreme thermodynamic conditions of pressure and temperature is of interest to many fields dealing with high energy density, such as astrophysics, geophysics, inertial confinement fusion, explosive and impact driven phenomena, etc. Theoretical study of these dynamic systems require detail knowledge of Equation of State (EOS) at wide range of pressures and energies.
Molybdenum is a refractory metal which is known to be stable in a body-centred cubic structure up to high pressures.Mo, its alloys and other body-centered cubic (BCC) refractory metals find wide range of applications in high pressure physics, e.g.planetary systems, and hypervelocity impacts of rocky planets, aerospace and armor/anti-armor applications, because of their excellent strength, mechanical stability, creep resistance and ductile nature.Most of the existing theoretical and experimental study on Mo are aligned towards examining structural stability, phase transition and melting under static pressures [1].Current work is motivated by recent interest in understanding response of single-crystal Mo under shock compression -a relatively unexplored domain.While most of the shock compression experiments carried out in laboratory are limited to strain rates ∼ 10 6 /s, ultrashort pulsed laser extends this limit to a great extent.However, the same for still higher strain rates, as encountered in many high energy density systems remains unexplored largely due to lack of experimental facilities.For this purpose large-scale nonequilibrium molecular dynamics (NEMD) study is preferred which facilitates extracting information about shock evolution and associated microstructural changes at such exotic environment.
Due to recent advances in computational resources classical MD simulations are being increasingly used to study the shock compression behaviour of crystalline solids.In this regard, recently we have examined the elastic-plastic behavior of Ta, another BCC metal, at ultra-high strain rates using NEMD simulation [2].Added advantage is that MD simulation provides atomistic insight of underlying physical phenomena.
Aim of the present work is to gain microscopic insights of shock propagation, construction of shock hugoniot of single-crystal Mo and utilize it to generate off-hugoniot states in the compressed region through Mie-Grüneisen EOS (MGEOS) (P − V − E relation).

Simulation details
All the results presented in this work are generated using open source classical molecular dynamics simulation package, LAMMPS [3].Atomistic visualization of shock propagation is facilitated by OVITO software [4].The simulation box contains 32 × 32 × 320 unit cells which corresponds to a dimension of 100 nm along z-direction and about 10 nm along x and y directions.As required for NEMD simulation, free boundary condition is applied along shock propagation direction (z), while periodic boundary conditions are used for the transverse directions.The interaction between atoms is treated by N-body semi-empirical potential parameterized by Ackland et al [5].Before shock simulations, the sample is initially thermalized at 300 K and zero pressure for 30 ps under NPT ensemble.The time step of integration is 1 fs.
Shock is imparted to the sample by setting a pre-decided constant velocity to atoms in first four unit cells of the simulation box on one end.This atomic layer behaves like a rigid-piston that compresses the sample.Thus generated shock is allowed to propagate along positive z-direction which coincides with the (001) crystallographic orientation of Mo crystal.Determination of resultant thermodynamic and macroscopic flow variables are described in Ref. [2] .

Shock Hugoniot
Shock is a propagating discontinuity traveling at certain high velocity (more than the sound velocity of the medium), resulting in an abrupt and irreversible change in state variables, such as density (ρ = 1/V ), pressure (P H ) and temperature (T H ).
The thermodynamic variables on either side of the wave are related through conservation relations with the parlance Rankine-Hugoniot jump conditions, mathematically expressed as: Hugoniot is the locus of all the points attained in single shock experiments performed at various impact velocities.In the present work, piston velocity is tuned to 1.8 -3.0 km/s so as to capture only pure plastic wave for the chosen crystal orientation.For illustration, we present visualization plot in Fig. 1(a) demonstrating the generation of shock and its propagation at different instances within the simulation box due to impact at 2.2 km/s.Shock front travels from left to right which can be identified from the abrupt change in kinetic energy of the atoms.Corresponding spatio-temporal evolution of impact generated pressure in the target (bin-averaged values) is displayed in 2D contour plot of Fig. 1(b).
Shock velocity (U s ) vs. particle velocity (U p ) hugoniot is displayed in Fig. 2(a).In case of uniaxial compression it is appropriate to use the component of pressure along shock direction, i.e., P zz as hugoniot pressure.Pressure (P H )-particle velocity hugoniot is shown in Fig. 2(b).Both the hugoniot obtained in NEMD simulation show good agreement with experimental data [6,7] pertaining to polycrystalline sample.Available single-crystal data for [001] direction [8] is also shown.This establishes the validity of the chosen interatomic potential for the present requirement in high-pressure applications.The linear-fitting of U s − U p curve, U s = C + S U p gives the hugoniot parameters C, i.e. bulk sound speed, and slope, S. Using our simulations carried out on Mo crystal, we get C = 5.24 km/s and S = 1.27.Using linear (U s − U p ) relation in Eq. 1 we can write hugoniot pressure as (2) 1300 ( 2024

Mie-Grüneisen EOS
Previously, we had constructed thermodynamically complete global EOS model [9] that utilizes additive property of helmholtz free energy.EOS data in tabular form was generated for range of materials including Mo.Thermal EOS which relates thermodynamic variables for materials can offer insightful details about their phase diagrams, melting, etc. Mie-Grüneisen EOS is one such analytical form of global EOS that utilizes shock hugoniot or cold isotherm as reference state.This specially finds its practicality in hydrodynamic simulations of shock-induced phenomena [10].
With hugoniot as reference, MGEOS can be written as In Eq. 3, γ(V ) is the Grüneisen parameter and is a function of specific volume V = 1/ρ.The analytical form of Grüneisen parameter is the well known expression, The constant 't' takes the value of 0, 1 and 2 for Slater, Dugdale-MacDonald and Vashchenko-Zubarev models respectively [11].P c (V ) is the pressure along zero temperature isotherm.We employ D-M model (t = 1) to determine ambient pressure Grüneisen coefficient (γ 0 ) given by wherein we have assumed hugoniot curve approaches cold isotherm at V → V 0 .Taking first and second derivative of Eq. 2 and putting V = V 0 we get, Using the above calculated value of S, we obtain γ 0 = 1.54.Incidentally, this value is in close agreement with experimental value (1.58) of Hixson [7].Numerous empirical relations have been proposed to describe the volume dependence of γ(V ).In the present work, we have used three models of Ref. [11,12,13] to obtain MGEOS data given as: Following Ref. [11], we take g 1 = 0.95, g 2 = 0.19 and α = 6.8 for Mo.The volume dependence of γ as predicted by Eqs. 7, 8 and 9 is shown in Fig. 3(a).It can be noted that γ 0 estimated using Ref. [11] [12] predicts faster decay with volume.Above three models of γ(V ) are then used to generate the P-V-E data (Eq.3) in the compressed regime as displayed in Fig. 3(b), (c) and (d).
Table 1 shows quantitative comparison of pressure values estimated by three models of γ for representative energy values of 5.0 and 15.0 MJ/kg.There is no appreciable difference among the pressure values estimated by three models.

Conclusion
The shock hugoniot of Mo single-crystal along [001] orientation has been obtained by NEMD simulation which shows good agreement with experimental results.The hugoniot parameters are utilized to estimate the Grüneisen coefficient at ambient pressure.Finally, MGEOS has been constructed using different models for volume dependent Grüneisen parameter.

Figure 1 :
Figure 1: (a) Atomistic visualization of shock propagation due to impact at U p = 2.2 km/s for three different times post-impact.Atoms are colored according to their kinetic energies with red being maximum and blue minimum (unshocked).(b)Spatio-temporal evolution of pressure inside Mo target for same impact velocity.

4 Figure 2 :
Figure 2: Hugoniot curves of single-crystal Mo obtained by NEMD simulation: (a) U s − U p ; (b) P H − U p .Experimental data for single and polycrystalline Mo are also shown.

Figure 3 :
Figure 3: (a) Variation of Grüneisen parameter with compression (V 0 /V ) for three different models.(b, c, d) MGEOS data for three different models of γ(V ) as given in Eqs. 7, 8 and 9.

Table 1 :
is higher than other two models, whereas γ when Comparison of P −V −E data for three different models of Grüneisen parameter