Constant-stress Hugoniostat method for following the dynamical evolution of shocked matter

We present an alternative equilibrium molecular dynamics method—the uniaxial constant-stress Hugoniostat—for following the dynamical evolution of condensed matter subjected to shock waves. It is a natural extension of the recently developed uniaxial constant-volume Hugoniostat [Maillet et al., Phys. Rev. E 63, 016121 (2001)]. Integral feedback is employed to reach the Hugoniot (final) state of the shock process by controlling both the normal component of the stress tensor and internal energy. The finite strain rate imposed on the system is closely related to that inherent in the front of a shock wave. The method can easily identify phase transitions along the Hugoniot shock states, even those that exhibit multiple wave structures. As an example of the method, we have simulated the Hugoniot of a Lennard-Jones crystal shocked along the $⟨110⟩$ direction. The results agree well with multi-million-atom nonequilibrium molecular-dynamics simulations.

Figure 1

(Color online) Time evolution of normal stress ${\sigma }_{zz}$ and shear stress $\tau$ $\left[2\tau ={\sigma }_{zz}-\left({\sigma }_{xx}+{\sigma }_{yy}\right)∕2\right]$ for a $25200\text{−}\mathrm{atom}$ $\mathrm{LJ}$ crystal shocked along the $⟨110⟩$ direction to a normal pressure $P_{zz}=40$, close to the Hugoniot elastic limit $\left(\mathrm{HEL}\right)$ for this crystallographic direction $\left(P_{zz}^{HEL}=41\right)$. The molecular-dynamics equilibration to this final state has been done for two cases: no damping $\left({\beta }_H=0={\beta }_p\right)$ and with damping $\left({\beta }_H=30,{\beta }_p=5\right)$. In both cases, the rate parameters have been set to ${\nu }_H=30$, ${\nu }_p=1$. The large initial fluctuations, evident in the undamped case, lead to plastic deformation below the actual $\mathrm{HEL}$ with an accompanying reduction in the shear stress. No plastic deformation is observed in the damped case. The inset shows the early time evolution of the normal stress for the damped case, which reaches the target equilibrium value at significantly earlier times than the undamped case.

Figure 2

(Color online) $\mathrm{NPzzHug}$ trajectory of normal stress ${\sigma }_{zz}$ as a function of volume $V$ from a simulation in which the target $P_{zz}=160$. Also shown is the Rayleigh line connecting the initial and final states. The final state is in the overdriven region of the Hugoniot. The hump in the trajectory is caused by initial overshoot of the $\mathrm{HEL}$ followed by plastic deformation.

Figure 3

(Color online) Heat into the $\mathrm{NPzzHug}$ system $\Delta Q$ and work done by the piston $-\Delta W$ as a function of time, as obtained from numerically integrating Eqs. (20, 21), respectively. The sum of the two terms is equal to $\Delta E=\frac{1}{2}\left(P_1+P_0\right)\left(V_0-V_1\right)$. Most of this energy change ($\approx 90\%$) is due to the work done by the piston.

Figure 4

(Color online) A generic Hugoniot. The dashed curve between the Hugoniot elastic limit $\left(\mathrm{HEL}\right)$ and the overdriven $\left(\mathrm{OD}\right)$ points is a region of the Hugoniot inaccessible from the initial state $O$. It is characterized by an elastic wave traveling at $u_s^{HEL}$ ($O\text{−}\mathrm{HEL}$ line) and a plastic wave whose velocity $u_s<u_s^{HEL}$.

Figure 5

(Color online) $\mathrm{NPzzHug}$ and $\mathrm{NEMD}$ shock velocity $u_s$ (normalized by the zero-pressure longitudinal sound speed $c_0=10.3$) vs $P_{zz}$ in $\mathrm{fcc}$ $⟨110⟩$-direction shocks. Label (0) indicates system started from initial zero-pressure state, while $\left(\mathrm{HEL}\right)$ refers to states obtained from re-centering at the $\mathrm{HEL}$. The horizontal dashed line has been drawn to help identification of the $\mathrm{OD}$ point, where the plastic wave speed equals $u_s^{HEL}$.

Figure 6

(Color online) $\mathrm{NPzzHug}$ and $\mathrm{NEMD}$ shock velocities $u_s$ as a function of particle velocity $u_p$ normalized to the longitudinal sound speed at zero pressure $\left(c_0=10.3\right)$. Also shown is a linear $u_s-u_p$ fit using the Hugoniostat data points above the $\mathrm{OD}$: $u_s=7.98+2.03u_p$.

Figure 7

(Color online) Comparison between $\mathrm{NEMD}$ and $\mathrm{NPzzHug}$ shock temperature as a function of normal $P_{zz}$.

Figure 8

(Color online) $\mathrm{NPzzHug}$, $\mathrm{NVHug}$, and $\mathrm{NEMD}$ normal pressure $P_{zz}$ vs volume of the $\mathrm{LJ}$ system shocked along the $⟨110⟩$ direction. $\mathrm{NVHug}$ simulations (triangles) in the unstable region of the Hugoniot give equilibrium temperatures and normal pressures below those at the $\mathrm{HEL}$. Several simulations are required in order to identify the stable region beyond the $\mathrm{HEL}$, where $P_{zz}⩾P_{zz}^{HEL}$. In constrast, $\mathrm{NPzzHug}$ gives only stable phases along the Hugoniot curve. Both methods agree well with $\mathrm{NEMD}$ simulations.

Figure 9

Cross-sectional profiles $\left(25\times 25r_0^2\right)$ of atomic configurations at $26\%$ compression from three types of simulations: (a) $\mathrm{NVHug}$ (Ref. 13), (b) $\mathrm{NPzzHug}$, and (c) $\mathrm{NEMD}$. The structure shown in (a) is not melted, but amorphous. The instantaneous compression in $\mathrm{NVHug}$ produces very large initial shear stress values and results in amorphization on a short time scale (within $2t_0$). In constrast, no such amorphization is seen in $\mathrm{NPzzHug}$ (b) or $\mathrm{NEMD}$ (c) simulations over the entire range of normal pressures reported here (i.e., below shock-induced melting). In addition, the observed structural deformation in the $\mathrm{NPzzHug}$ simulations closely resemble those produced in large-scale $\mathrm{NEMD}$ simulations.