Molecular Dynamics Simulations of Shock Wave Propagation through

Subscriber access provided by CORNELL UNIVERSITY LIBRARY

Article

Molecular Dynamics Simulations of Shock Wave Propagation Through the Crystal-Melt Interface of (100)-Oriented Nitromethane Shan Jiang, Thomas D. Sewell, and Donald L Thompson J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.6b07002 • Publication Date (Web): 08 Sep 2016 Downloaded from http://pubs.acs.org on September 11, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry C is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Molecular Dynamics Simulations of Shock Wave Propagation through the CrystalMelt Interface of (100)-oriented Nitromethane

Shan Jiang, Thomas D. Sewell,1 and Donald L. Thompson2,* Department of Chemistry, University of Missouri-Columbia, Columbia, MO 65211-7600

1 2,

Electronic mail: [email protected]; Tel: 573-882-7725 Electronic mail: [email protected]; Tel: 573-882-0051

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 37

ABSTRACT: Molecular dynamics (MD) simulations were used to study shock wave passage with normal incidence through the equilibrium interface between (100)-oriented nitromethane and the melt. The simulations were performed using the fully flexible, non-reactive SRT force field (Sorescu, D. C.; Rice, B. M.; Thompson, D. L. The Journal of Physical Chemistry B 2000, 104, 8406-8419). The local kinetic energies (intermolecular, intramolecular, and total) and stress states differ significantly in the liquid and crystal regions, and depend on whether the shock is initiated in the crystal or liquid. The number and spatial distribution of shock-induced molecular disorientations in the crystal for shocks initiated in the crystal are similar to these obtained for analogous simulations for a completely crystalline sample; however, substantial differences in the extent and distribution of shock-induced molecular disorientations in the crystal region were observed when the shock was initiated in the liquid. All three measures of kinetic temperature in the crystal region are higher when the shock is initiated in the crystal than when it is initiated in the liquid. Kinetic temperature profiles exhibit features in the vicinity of the interface considerably different from those in either bulk phase. The shock-induced local mechanical states (von Mises stress) indicate that the crystal is less able to support shear stresses when the shock is initiated in the crystal than when it is initiated in the liquid. There is a strong reflection back into the liquid when the shock wave passes through the liquid and encounters the interface with the crystal. This causes a large increase in the potential energy of the liquid, limits the amount of energy transmitted into the crystal, which also limits the molecular disorientations in the crystal. Thus, a shock from liquid to crystal yields less inelastic deformation in the crystal.


Page 3 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


INTRODUCTION We are interested in understanding the fundamental dynamics of organic materials in response

to thermal and mechanical perturbations.1-12 We have often used the simplest nitroalkane, nitromethane, as a prototype for studies of both homo-phase and hetero-phase properties and processes in molecularly complex condensed phase systems.2,4-7,10-13 We have reported molecular dynamics (MD) studies of the solid14 and liquid15 phases of nitromethane in which we developed the flexible-molecule Sorescu-Rice-Thompson force field (SRT FF). We have used the SRT FF in studies of melting initiated at the crystal-melt interface13,16 and the crystal-vacuum interface.1,17 Simulations by Agrawal et al.16 of coexisting liquid and crystal nitromethane using the SRT FF predicted a pressure-dependent melting temperature in good agreement with experiment. Using that same force field Siavosh-Haghighi and Thompson1 investigated melting initiated at solid-vacuum interfaces and Siavosh-Haghighi et al.13 simulated the crystallization of nitromethane from the melt. Those studies were focused on thermodynamic processes. Here we report the results of a study of the effects of a mechanical perturbation, that is, shocks. Nitromethane is a liquid under standard ambient conditions. Typically, for liquid explosives shock wave detonation initiation occurs via a homogeneous mechanism wherein material behind the shock front is initially uniformly compressed and heated. Ignition of chemistry is essentially thermal and initiation occurs when local reaction fronts, which propagate supersonically compared to the uncompressed liquid, coalesce and eventually reach and reinforce the lead shock front resulting in a detonation wave. This process is clearly described in the early studies of Campbell et al.18 They reported that chemistry begins in the vicinity of the liquid-solid interface between nitromethane and the metallic or polymeric shock attenuator that serves to transmit shock energy from a donor explosive into the nitromethane. The presence of bubbles, suspended particulates, or



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 37

surface roughness of the shock attenuator can decrease the threshold shock strength for initiation. These effects have been studied in some detail using modern diagnostics by Dattelbaum et al.,19 who studied the detonation initiation behavior of liquid nitromethane (gelled using a low concentration of guar gum) in the presence of solid silica beads and hollow glass microballoons. Their studies included the effects of particle size, impurity concentration, inter-balloon spacing, and projectile velocity. Some MD studies of shock waves and shock-induced chemistry in organic energetic materials containing interfaces have been reported. Shi and Brenner20 used MD to study shocks traversing interfaces with various degrees of nanoscale faceting between an inert binder (represented by a single-component Lennard-Jones model) and crystalline nitrogen cubane (described by a reactive empirical bond-order model). Their results demonstrate that a nanometer-scale faceted-structure interface between a binder and energetic crystal can influence sensitivity. Long et al.21 developed ab initio based pair potentials for HMX-TATB and TATB-graphite interfaces (HMX: octahydro1,3,5,7-tetranitro-1,3,5,7-tetrazocine; TATB: 1,3,5-triamino-2,4,6-trinitrobenzene), and used them in MD shock simulations of TATB-coated HMX and graphite-coated HMX that showed that the hotspots tend to form at the interface and are heated by plastic work.22 An et al.23,24 reported MD simulations using the ReaxFF25,26 reactive force field of a plastic-bonded explosive (PBX) consisting of RDX and HTPB23,24 and another consisting of PETN/Si-PETN and HTPB27 (RDX: hexahydro-1,3,5-trinitro-1,3,5-s-triazine; PETN: pentaerythritol tetranitrate; HTPB: hydroxylterminate polybutadiene). They studied shocks propagating from opposite directions through sawtooth-shaped interfaces and identified the regions where hotspots formed for a given shock direction. An et al.23 proposed an approach to de-sensitize the RDX-HTPB composite by applying a thin coating of low-density (0.48 g/cm3) HTPB on the RDX, yielding a model material that they


Page 5 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


called interface-matched PBX (IM-PBX). Their simulations showed that the IM-PBX design concept dramatically suppresses hotspot formation. These studies were of interfaces between two solids, and showed that the interfaces significantly affect the shock physics. We are interested in how a solid-liquid interface affects the physics. Several MD simulations of physical and chemical processes that focused on shocked nitromethane have been reported. Seminario et al.28 performed MD simulations of liquid nitromethane to investigate structural changes and vibrational frequency shifts that occur upon compression from 1 bar, 300 K to 143 kbar, 600 K. Reed et al.29 predicted a transient semi-metallic layer in detonating nitromethane from 16-molecule quantum MD simulations. Siavosh-Haghighi et al.4,5 studied the deposition of energy into molecular normal modes upon shocking of crystalline nitromethane. He et al.6,7 studied plane-specific effects in the shock wave propagation and postshock energy redistribution for several different crystal orientations. Rivera-Rivera et al.10 carried out MD simulations of shocked (100)-oriented nitromethane to determine the time and space scales for the approach to local thermal equilibrium behind a shock wave. These simulations have dealt with nitromethane initially in a single phase; here we report a study of the effects of shocking across the nitromethane crystal-melt interface. We used MD simulations to study shock passage with normal incidence through the interface between a (100)-oriented crystal in equilibrium with the melt from both the crystal and liquid ends of the system. The aim is to better characterize and understand the fundamental effects of shocks at interfaces. Toward this end, we have computed as functions of time and location in the sample a locally averaged orientational order parameter and two-dimensional (2D) mean square displacement (MSD), inter-, intra-, and total-molecular kinetic temperatures, potential energy, and mean and von Mises stress. Understanding shock wave propagation through the nitromethane crystal-liquid



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 37

interface may serve as a reference point for studies of shocks across interfaces between nitromethane and other substances, such as a solvent, chemical or physical sensitizing agent, or a containment vessel.30 COMPUTATIONAL METHODS The simulations were performed using the LAMMPS31 code with the flexible molecule SRT FF,14,16 which has been shown to accurately predict the static and dynamic properties of nitromethane in the crystalline and liquid states as well as transitions between them.2,4,6,7,13-16 The intramolecular potential is described using the Morse function for covalent bonds, harmonic functions for the bond angles, and cosine terms for dihedral and improper dihedral angles.14,16 Non-bonded interactions, which only apply between atoms in different molecules, are described by the Buckingham exponential-6 potential plus the Coulomb potential for fixed partial charges.14 The functional forms and parameters can be found in Agrawal et al.16 The repulsion and dispersion interactions and the real part of the Coulomb interaction were truncated at 11.0 Å. Long-range Coulomb forces were evaluated using the particle-particle particlemesh (PPPM) algorithm32 with the relative error in forces set to 10-6. Nitromethane crystallizes in the orthorhombic space group P212121 in a unit cell containing four molecules.33-35 The experimental36,37 and theoretical16 values of the normal melting point are 244 K and 266.5±7.7 K, respectively. The theoretical value was obtained from MD simulations of void-nucleated melting of the three-dimensional periodic crystal. Two-phase MD simulations yielded a crystal-liquid coexistence temperature of 255.5 K at 1 atm.16 The starting point for the shock wave simulations was an equilibrated interface between the (100) crystal face and the melt. We used a procedure for constructing a co-existing crystal and melt similar to one described previously.16 Starting from the experimental lattice parameters reported by Trevino et al.35 (a = 5.244 Å, b = 6.320 Å, and c = 8.730 Å at 224 K and 1 atm), a 3D periodic orthorhombic simulation cell containing 69,120 molecules was constructed by replicating


Page 7 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


, and crystal directions, respectively. The , the unit cell by factors of 40, 24, and 18 along the resulting supercell was equilibrated for 50 ps at 256 K and 1 atm in the isobaric-isothermal (NPT) ensemble using the Nosé-Hoover thermostat and barostat38,39 with relaxation constants 10 fs and 100 fs, respectively. The same values of the relaxation constants were used in the other NPT and isochoric , and were , isothermal (NVT) simulations described below. The pressure components along controlled independently by the barostat. The simulation box was constrained to be orthorhombic during the NPT simulation. Starting from the configuration at the end of the NPT simulation for the crystal, a liquid simulation cell was generated by re-setting the atomic velocities to values selected from the Maxwell distribution for T = 400 K and then equilibrating the system for 20 ps in the NPT ensemble at 1 atm and 400 K. Next, new atomic velocities were selected from the Maxwell distribution for 256 K and the system was simulated in the NPT ensemble for 4.0 ps during which additional velocity reselections from the Maxwell distribution for 256 K were performed at 1.0 ps intervals (including the final time, t = 4.0 ps). Starting from phase space point at the end of that simulation 10.0 ps of additional NPT equilibration at 256 K was performed without reselection of velocities. Finally, the transverse simulation cell lengths for the liquid were adjusted to match those for the crystal supercell by varying them linearly over the course of a 20.0 ps NVT simulation during which the longitudinal cell length was adjusted to preserve the system volume as a function of time. The crystal and liquid supercells with identical transverse edge lengths were combined to yield a 3D periodic supercell containing 138,240 molecules with two crystal-melt interfaces oriented perpendicular to the x-axis. Prior to combining the crystal and liquid supercells the molecules were unwrapped along the x-axis, and a small (≈5 Å) vacant space was introduced between the liquid and



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 37

crystal to avoid atomic overlaps. The resulting system was relaxed in three steps, each with a 0.2 fs time step: (1) NVT equilibration at 256 K for 24 ps, with velocity reselection at t = 0 and at 1 ps intervals for the next 4 ps; (2) NPT equilibration at 256 K and 1 atm for 20 ps, which led to a small difference in temperature for the liquid and crystal regions; and (3) starting from the final phase point of step (2), velocity reselection from the 256 K Maxwell distribution to eliminate the small temperature difference observed in (2) followed by system equilibration for 20 ps in the NVT ensemble. After the third step both ends of the two-phase sample were exposed to vacuum, and then additional NVT equilibration was performed for 25 ps with a 0.5 fs time step. Velocity reselection was performed every 1 ps during the first 5 ps to eliminate breathing motions along the shock direction that result from the sudden introduction of the free surfaces. Figure 1 illustrates the simulation supercell. The liquid is on the right and the crystal on the left with vacant spaces at either end. Frame (a) illustrates the configuration for the liquid-to-crystal shock and (b) for the crystal-to-liquid shock. The shock wave propagation directions are shown by solid red arrows. The crystal is oriented so the shock waves pass with normal incidence through the (100) face. Shocks propagating parallel to the x-axis (i.e., parallel or anti-parallel to the [100] crystal direction) were generated using a reverse-ballistic configuration in which the sample is impacted with speed |up| = 1.0 km/s onto a stationary piston of rigid molecules. The piston was created by fixing the velocities and forces at zero on all atoms in molecules with center-of-mass positions within 15 Å of the end of the sample from which the shock was initiated. The impact velocity vector up = (±1.0, 0.0, 0.0) km/s (blue arrows in Figure 1) was added to the instantaneous thermal velocities of the atoms in the flexible slab, after which the system was propagated in the microcanonical (NVE) ensemble with time step 0.1 fs.


Page 9 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Fixed-volume bins of width ≈4 Å parallel to the shock direction were defined for the analysis of some properties.4,6 The bins are defined in the lab frame and are Eulerian in the sense that material flows through them as the simulation proceeds. The temporal-spatial distributions of several properties were obtained by averaging over all the molecules in a given bin at a given time. The center-of-mass positions of the molecules were used to determine to which bin they belonged. Prior to shock compression, each bin in the crystal region contained ≈1,400 molecules and each bin in the liquid region contained ≈1,200 molecules (the crystal and liquid densities differ). Similar to procedures used in a previous study of shock waves in single crystal nitromethane,6 the instantaneous shock front location and the shock speed were determined by monitoring the intermolecular, intramolecular, and total kinetic energies and the stress tensor in each fixed-volume bin as functions of time. Although the simulations are for supported shocks, due to the multiple waves generated when the shock crosses the interface it is not feasible to improve the statistical precision of the results by transforming the raw simulation data from the lab frame to a translating reference frame centered on the shock front as was done in earlier studies in which the shocks were essentially steady.6 Thus, except where specifically noted below, such averaging was not performed. Usually in a simulation of an explicit shock wave (as is studied here) the simulation is stopped when the shock front reaches a free surface, because after that a rarefaction wave will propagate back through the sample. For the system sizes studied here, this limits the time of the simulations to ≈11 ps. To extend the time for which the shocked material could be studied we applied the shock-

front absorbing boundary condition (SFABC)3,40 in some cases. In the SFABC method the shock wave is propagated until it completely traverses the sample (i.e., until the instant of maximum compression), at which time a second piston is created by setting to zero and holding at zero the velocities and forces of all molecules with centers of mass within a prescribed distance (here, 15 Å)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 37

of the end of the sample. This captures the shock-compressed material in the volume and energy state at the instant of maximum compression. Details of the SFABC method can be found in Bolesta et al.40 and Cawkwell et al.3 In the present case, for which multiple waves are generated due to scattering of the incident shock wave at the crystal-liquid interface, application of the SFABC results in complicated wave structures due to reflections from both pistons combined with numerous wave-interface and wave-wave interactions within the sample. However, for our purpose the ability to sustain the material in the NVE state initially prepared by shock passage outweighs the effects of artifacts introduced by the SFABC. We used a locally averaged orientational order parameter to determine the structural environment (i.e., locally crystalline or amorphous) of each molecule in the system as a function of time. A similar approach was used by Benet et al.41 to characterize a water crystal-melt interface. The molecular orientational order parameter for the ith molecule is identical to that used by Siavosh-Haghighi and Thompson,1 namely, the absolute value of the cosine of the angle between the direction of the C–N bond and the c lattice vector: = | ∙ |,

(1)

where is the unit vector parallel to the C-N bond at a given time and is the unit vector parallel to c. The locally averaged orientational order parameter for each molecule j is defined as

= ∑ ,

(2)

where N(j) is the number of neighboring molecules with center-of-mass positions within a cutoff radius (7 Å) from the center of mass of molecule j (including the jth molecule) at a given instant. For the (100)-oriented nitromethane crystal is close to 1.0 because the angles between the C–N bonds and lattice vector c are close to 0° and 180°. The value of decreases to ≈0.5 for molecules in an amorphous environment. Note that because only provides local structure and 10 ACS Paragon Plus Environment

Page 11 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


not dynamics, it does not trivially distinguish a local liquid environment from a local glassy state. Intermolecular, intramolecular, and total kinetic temperatures Tinter, Tintra, and Ttotal, respectively, were calculated for the molecules in fixed-volume bins at fixed time intervals. The Tinter, Tintra, and Ttotal are defined as6 = . = 5.6 =

∑()* ! "〈!〉%& ' , +,- 4

∑()* ∑( 0)* /0 120 "! 3

(3) '

+ 4 " ,-

44

∑( 0)* /0 120 "〈!〉%& 3 +,- 44

,

(4)

'

;

(5)

where 7 and 8 are the mass and center-of-mass velocity of molecule i; 9 and 2 are the mass and velocity of atom j in molecule i; 〈!〉:/ is the center-of-mass velocity of all the molecules in the fixed-volume bin where the center of mass of molecule i is located; ; is the number of molecules in the bin; ;′ = 7, the number of atoms in a nitromethane molecule; ;′′ = 7N; and => is the Boltzmann constant. We emphasize that eqs (3)-(5) are kinetic temperatures; that is, kinetic energies scaled to temperature units. As noted above, the locally averaged orientational order parameter does not distinguish between liquid- and glass-like phases; thus to distinguish between them the transverse 2D MSD was calculated in each fixed-volume bin in the initial crystal region as a function of the post-shock time using MSDB = 〈|C B − C 0 |F 〉,

(6)

where C 0 and C B are the center-of-mass location of a given molecule at the time origin (defined as the instant at which the shock front reaches the bin containing the given molecule) and subsequent times, respectively. The quantity C is the 2D projection of the 3D position vector into 11 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 37

the shock front plane (i.e., the y-z plane in the lab frame). Because only the crystal region is considered and because the speed of the lead shock is essentially constant in the crystal, it is valid for this property to perform the analysis in a translating reference frame for which the zero of time in a given fixed-volume bin is defined as the time the lead shock front reaches the boundary of that bin. With this choice the time t in eq 6 is the time molecules in a given bin have been in a shocked state and the ensemble average is over the set of molecules that have been in a shocked state for that time t. Stress states were calculated as functions of location along the shock direction by averaging over the atomic stress tensor for atoms belonging to the molecules in a given fixed-volume bin at a given time. We considered two types of stress quantities, namely the volumetric stress GH (also known as the hydrostatic stress or mean stress):

GH = + G + GFF + G++ ,

(7)

where G , GFF , and G++ are the diagonal components of the stress tensor; and the von Mises stress GHJK (also known as the equivalent stress ), given by

GHJK = L3NO = PF QG − GFF F +G − G++ F +GFF − G++ F + 6GF F + G+ F + GF+ F S,

(8)

where NO is the second deviatoric stress invariant and the G = 1,2,3; = 1,2,3 are the components of the stress tensor.

RESULTS AND DISCUSSION Figure 2 shows Tinter, Tintra, and Ttotal (solid, dashed, and dotted curves, respectively) as functions of position along the length of the sample (see Figure 1) at various times during the simulations. Panels (a) and (b) are for liquid-to-crystal and crystal-to-liquid shock propagations, respectively.


Page 13 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


For both shock directions Tinter rises rapidly as the shock front passes and then decreases toward approximately constant values behind the shock on the 10 ps time scale of the simulations. The initial increase in Tinter is ≈90% larger when the shock is initiated in the liquid (≈370 K) than when it is initiated in the crystal (≈195 K), while Tintra and Ttotal increase more slowly for both the crystal and the liquid regions. The curves for Tintra and Ttotal in both the crystal and the liquid increase approximately monotonically to long-time values. This behavior is consistent with the initial shock excitation energy being deposited mainly in low-frequency motions that are dominated by intermolecular motions, followed by equipartitioning of energy among all the degrees of freedom.6,42 Generally, all three kinds of kinetic energy in the initial crystal region reach higher values when the shock is from crystal-to-liquid (see Figure 2b) than when from liquid-to-crystal (see Figure 2a). For example, focusing on t = 10 ps in both cases (green curves) for the liquid-to-crystal shocks (see Figure 2a) when the shock has nearly passed through the sample Ttotal in the crystalline region approaches ≈310 K; whereas Ttotal in the same region for the crystal-to-liquid case (see Figure 2b) is ≈370 K. This dependence of thermal (and mechanical) properties along the shock direction at the same impact speed on whether the shock is initiated in the crystal or the liquid is revisited below (see Figures 4, 7, and 8) when we discuss the results using time-position (t-x) diagrams. Figure 3 shows instantaneous profiles of Tinter for several times that bracket passage of the shock wave through the interface. Profiles in the vicinity of the interface exhibit features considerably different from those in either bulk phase. As can be seen in Figure 3a, which is for the liquid-tocrystal shock, Tinter decreases as the shock moves through the interface, but when the shock is crystal-to-liquid (Figure 3b) Tinter increases as the shock passes through the interface. This can be explained by the difference in the densities of the crystal and the melt and the associated different



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 37

amounts of P-V heating. However, as will be shown below, there are interesting effects on the microscopic states of the shocked crystal region and the associated mechanical properties behind the shock depending on whether it a liquid-to-crystal or crystal-to-liquid shock. Time-position (t-x) diagrams illustrate the evolution of local properties of the system as simultaneous functions of time and space. They are commonly used for the interpretation of shock dynamics in 1D systems (which is essentially the case here), and are used in Figures 4, 7, and 8, where the abscissa denotes position in the lab frame parallel to the shock direction and the ordinate is the elapsed time since initiation of the shock; the color at a given (t, x) point corresponds to the value of the property at time t averaged over the volume of material in a given fixed-volume bin centered at position x in the sample. Figures 4, 7, and 8 contain, respectively, results for the intermolecular temperature (Tinter), average potential energy per atom (PE), and volumetric and von Mises stresses (GH and GHJK ) for both liquid-tocrystal shocks (top panels) and crystal-to-liquid shocks (bottom panels). The x = 0 is defined to coincide with the interface plane at the start of the shock simulations, and the crystalline and liquid regions are located to the left and right, respectively, so that liquid-to-crystal shocks propagate from right to left and crystal-to-liquid shocks from left to right. Dark blue regions at the left and right ends of the t-x diagrams correspond either to empty space or to pistons. Results for Tinter for open boundaries along the shock direction, for which the system undergoes shock compression followed by decompression, are shown in the left-hand column, panels (a) and (b), of Figure 4. Corresponding results obtained when the SFABC is imposed starting at the time of maximum compression are shown in the right-hand column, panels (c) and (d). Thus, horizontally adjacent results in panels (a) and (c) and similarly in panels (b) and (d) are identical up to the times at which the SFABC was applied. For the liquid-to-crystal shock where the SFABC was not used (see Figure 4a) Tinter in the liquid is much higher than in the crystal. This difference persists until significant release of the


Page 15 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


compressive state occurs at approximately 20 ps into the simulation. When for the same shock direction the SFABC is applied at 11 ps (see Figure 4c) there is no release and Tinter for the crystal and liquid do not approach the same value within the 27 ps duration of the simulation. Note the partial reflection of the lead shock when it encounters the interface. And, also note that when the SFABC is applied additional and enhanced-amplitude reflected waves are generated due to imperfect capturing of the material in the non-uniform state between the two pistons. When the shock is crystal-to-liquid (see Figures 4b and 4d), the maximum Tinter (red and orange regions from 0 ps to ≈20 ps) is still much lower than that reached when the shock is liquid-to-crystal (pink and white regions from ≈5 ps to ≈15 ps). Similar to Figure 4a, the temperature difference across the interface finally disappears due to the expansion and relaxation in the sample, for the open boundary crystal-to-liquid case (green regions from ≈20 ps to the end; see Figure 4b). However, this temperature difference between the two regions persists to 27 ps, the end of the simulation, 17 ps after applying the SFABC (see Figure 4d). This effect, which is seen for both shock directions, can be explained by the fact that the second piston prevents the tensile strain of the sample and more of the work of the shock loading is transferred into heat rather than strain energy as would occur under tensile expansion. Particularly apparent in Figure 4d is the generation of several waves that originate both from partial reflection at the interface and from imperfect SFABC capture of the lead shock when it reaches the end of the liquid region. There are substantial differences in the spatial distributions of shock-induced molecular disorientations in the crystal for the two different directions of shock propagation. The instantaneous molecular center-of-mass locations are shown in Figure 5 for (a) a liquid-to-crystal shock, (b) a crystalto-liquid shock, and (c) a shock through a perfect crystal. Molecules are colored based on the values of the locally averaged orientational order parameter (see eq 2): Red if the parameter is greater than 0.95 (an empirically determined threshold criterion for belonging to the crystal structure), yellow if less than 0.95



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 37

but larger than 0.8, green if less than 0.8 but larger than 0.5, and blue if 0.5 or smaller. All three systems contain the same number of molecules and the cross-sectional areas are all identical. The initial length of the perfect crystal sample is slightly less than twice those containing an interface because the density of the crystal is higher than that of the liquid. The SFABC was applied in all three cases; at ≈11 ps in panel (a) and at ≈10 ps in panels (b) and (c). Molecules initially in the crystal that become disoriented are regarded as crystal defects. As shown in Figure 5a, which is for the liquid-to-crystal shock, the defects in the crystal nucleate at and close to the crystal-melt interface. Note the “roughness” at the interface, which increases with increasing time as the defects grow into the crystal. Also notice how defects eventually propagate all the way across the crystalline region. The defect pattern consists mainly of disoriented molecules located on planes at approximately 45° relative to the shock propagation direction. Although extensive molecular disorientations occur (green), much of the crystal structure remains (red). The crystalline region of the sample clearly does not amorphize completely. Corresponding results for the crystal-to-liquid shock are shown in Figure 5b, the center column. In this case the defects (green) nucleate at or close to the piston then propagate into the crystal resulting in a severe inelastic deformation state; the inelastic deformation in the crystal structure is much greater in this case than for the liquid-to-crystal shock. Note that the locally crystalline regions that exist at the time of application of the SFABC tend to persist for the remainder of the simulation. It is interesting to compare the results in Figure 5b with those for a similar shock simulation of an entirely crystalline sample; the results for which are shown in Figure 5c. Focusing on the left half of the pure crystal sample that corresponds to the initially crystalline region (denoted using inscribed rectangles in the top and bottom frames) in the crystal-melt sample, it can be seen that for both cases for which the shock is initiated in the crystal the patterns of the shock-induced molecular disorientations are quite similar. Although the left-


Page 17 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and right-hand ends of the initially purely crystalline sample in Figure 5c exhibit noticeable differences at longer times (20 ps – 30 ps), the results in Figures 5b and 5c suggest that on the short time- and spacescales studied, the presence of the crystal-melt interface does not strongly affect short-time defect evolution in the crystal when the shock is crystal-to-liquid. This conclusion is further supported by the results shown in Figure 6. For the analysis in Figure 6, all molecules initially located in the crystalline regions (the inscribed rectangles in Figure 5) were used to calculate the fractions of disoriented molecules and 2D MSD as functions of time in any shocked state. For this analysis the zero of time for each molecule is defined as the time the lead shock reached a given bin in which that molecule is located, and the dependent-variable values shown in Figure 6 at a given post-shock time were obtained by averaging over the subset of shocked molecules that had been in any shock state for that amount of time. It can be seen from Figure 6a that more molecular disorientations are generated in the crystal for the crystal-to-liquid (blue curve with light-blue circles) shock than for the liquid-to-crystal (red curve with orange squares) shock. However, for both cases for which the shock is initiated in the crystal (blue and green curves) the fraction of the disoriented molecules eventually reaches approximately the same value at approximately 6 ps after shock passage. These results are consistent with the molecular disorientations shown in Figure 5, and suggest that less severe inelastic deformation will occur in a crystal if the shock wave is transmitted liquid-tocrystal. This implies that the more compliant liquid acts as a “shock absorber” that protects the crystal and thus limits the generation of defects in it. The MSD results in Figure 6b, calculated as described in connection with eq 6 and which reflect molecular translational mobility, are broadly consistent with those for the fraction of disoriented molecules shown in Figure 6a. Figure 7 shows t-x diagrams for the potential energy (see panels (a) and (b)) and the potential energy increase (see panels (c) and (d)) for the liquid-to-crystal shock and the crystal-to-liquid shock. The top



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 37

panels are for the liquid-to-crystal shock and the bottom panels for crystal-to-liquid. The increase in potential energy for each atom (denoted as ∆PE(i)) is obtained from the difference in the potential energy at time t and t = 0, and then the increase in potential energy for each bin is calculated by averaging ∆PE(i) over the numbers of atoms in that bin. Similar to the results for temperature shown in Figure 4, the profile of the potential energy across the length of the sample depends on the shock direction. The potential energy increase of the crystal for the liquid-to-crystal shock (note the green region in Figure 7c between ≈6 ps and ≈9 ps and the corresponding change from dark blue to dark green in Figure 7a) is substantially less than that for the crystal-to-liquid shock (note the purple region in Figure 7d between ≈1 ps and ≈4 ps and the corresponding change from dark blue to light green and yellow in Figure 7b). Correspondingly, the magnitude of the potential energy increase of the liquid (note the light purple region in Figure 7c between ≈6 ps and ≈9 ps and the corresponding change from light green to dark orange and red in Figure 7a) for the liquid-to-crystal shock is higher than that for the shock from the other direction (note the dark purple and red region in Figure 7d between ≈6 ps and ≈9 ps and the corresponding change from light blue to light orange in Figure 7b). By construction, the total external energy input to the system is nearly the same for both shock directions. The results in Figures 7a,c show that there is a strong reflection back into the liquid when the shock wave passes through the liquid and encounters the interface with the crystal, leading to a larger potential energy increase in the liquid; which in turn limits the amount of energy transmitted into the crystal compared to when the lead shock travels in the opposite direction. The increase of the potential energy reflects the strain energy gained in the material. Thus, the larger increase of the potential energy in the liquid demonstrates that the liquid absorbs more energy when the shock wave is transmitted from the liquid into crystal, which explains why it can act as a “shock absorber” and thereby results in fewer disoriented molecules in the crystalline region (see Figure 5a and Figure 6a) than when the shock is transmitted from the crystal into the liquid.


Page 19 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Stress states were evaluated for both shock directions; they are shown in Figure 8. Figures 8a,b show the time evolution of the volumetric stress GH profiles calculated using eq 7 for systems with SFABC applied. Because GH , which is the negative of the mean pressure, does not include the deviatoric stress components, it only reflects changes in the volume of the stressed body. One can see that after the shock wave passes through the entire sample (from ≈10 ps to ≈15 ps) the absolute volumetric stress in the liquid is much higher for the liquid-to-crystal shock (dark purple region with value ≈ -4.5 GPa) than for the crystal-to-liquid shock (light purple and red regions with value around -4 GPa). Consistent with the result shown in Figure 7a, this shows that the liquid absorbs more energy, resulting in a higher pressure state for the liquid-to-crystal shock direction than for the opposite case. These results are also consistent with the conclusions drawn from the results in Figure 6. The von Mises stress GHJK is shown in Figures 8c,d for, respectively, the liquid-to-crystal and crystal-to-liquid shocks with SFABC. It is commonly used to predict the yield behavior of materials and satisfies the underlying property that two stress states with equal distortion energy have equal von Mises stresses.43 Thus, von Mises stress is also called equivilent stress. In solid mechanics the von Mises criterion suggests that yielding of the crystal begins when the elastic energy of distortion reaches a critical value, irrespective of detailed stress states.44 Because GHJK is directly determined by the second invariant of the deviatoric stress tensor, which is generally close to zero for liquids, it can be used to discriminate the liquid regions from locally crystalline or glassy enviroments. The profiles of GHJK in Figures 8c,d illustrate that there is a large difference between the crystal and the liquid. Because the crystal can resist shear deformation, it exhibits much larger GHJK (especially immediately after the shock front passes) than the liquid does, for both shock directions. Comparing panels (c) and (d), one can see that GHJK in the crystal remains large for a much longer time when the shock is initiated in the liquid (see the green region from ≈5 ps to ≈25 ps panel (c)) than when it is initiated in the crystal (see the green



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 37

region from 0 ps to ≈10 ps in panel (d)). This is because for the liquid-to-crystal shock more of the initially crystalline region in a given fixed-volume bin remains crystalline than for the opposite shock direction, and therefore can support larger shear stresses. This is consistent with the more extensive inelastic deformation that develops in the initially crystalline region when the shock is initiated in the crystal than when it is initiated in the liquid. It also can be found in Figure 8d that when the shock is initiated from the crystalline side, the relatively large magitude of GHJK (> 2.5 GPa) in the crytal only lasts for about 1 ps after the shock front passes, and then decreases to smaller values (

Molecular Dynamics Simulations of Shock Wave Propagation through

Recommend Documents