For the first time, molecular dynamics simulations (MD) are applied to study the evolution of shock-driven Richtmyer-Meshkov instability of the perturbations at a solid surface and at the interface separating two solids of different densities. We study the dependence of the instability evolution on the strength and orientation of the shock, and analyze the structure of the material in a vicinity of the interface perturbations (i.e. the spikes and bubbles)

Abarzhi, S. I.

Bringa, E. M.

Remington, B.

Zhakhovskii, V. V.

Zybin, S. V.

Caltech Authors - Main

 
 
 
MOLECULAR DYNAMICS SIMULATIONS OF THE RICHTMYER-
MESHKOV INSTABILITY IN SHOCK LOADED SOLIDS 
 
S. V. Zybin1, V. V. Zhakhovskii2, E. M. Bringa3, S. I. Abarzhi4, and B. Remington3 
 
1Materials and Process Simulation Center, California Institute of Technology, Pasadena, CA  91125, USA 
2 Institute of Laser Engineering, Osaka University, Osaka, 565-0871, Japan  
3 Lawrence Livermore National Laboratory, Livermore CA  94550, USA 
4 Center for Turbulence Research, Stanford, CA 94305, FLASH Center, The University of Chicago, IL 60637 
 
Abstract. For the first time, molecular dynamics simulations (MD) are applied to study the evolution 
of shock-driven Richtmyer-Meshkov instability of the perturbations at a solid surface and at the 
interface separating two solids of different densities. We study the dependence of the instability 
evolution on the strength and orientation of the shock, and analyze the structure of the material in a 
vicinity of the interface perturbations (i.e. the spikes and bubbles). 
 
Keywords: molecular dynamics, elastic-plastic shock, instability, Richtmyer-Meshkov, shear strength 
PACS:  02.70.Ns, 47.20.Bp, 52.57.Fg, 62.50.+p, 68.35.Ja 
 
 
INTRODUCTION 
 
The dynamics of a solid under deformations of 
very high strain rate and the development of 
instabilities at the solid interfaces are of great 
interest in inertial confinement fusion (ICF), 
detonations in heterogeneous (e.g. polymer-
bonded) explosive materials, ejecta formation, and 
many other applications. It is also one of key issues 
for theories and models, which study material 
strength and phase transitions under a high strain 
rates. The instabilities developed include the 
Rayleigh-Taylor (RTI) and Richtmyer-Meshkov 
(RMI) instabilities, and they occur at interfaces 
between two solids subjected to a sustained and 
impulsive shock-driven acceleration respectively. 
In solids RTI and RMI develop if the value of the 
acceleration (initial shock) exceeds a certain 
threshold, which is set by the material properties 
and can be used to measure these properties under 
high strain rate deformation. Such phenomena have 
been studied both experimentally [1-4] and 
theoretically [5-7], and recently by the molecular 
dynamics (MD) approach [8-9]. Here, we report 
some preliminary results of MD simulations of the 
evolution of the shock-driven RMI in solids. 
 
SIMULATION APPROACH 
  
We conducted MD simulations of solids, when 
a planar shock wave refracts through the perturbed 
interface between 1) two different solids, and 2) a 
solid and a vacuum. For the solid-solid interface, 
we consider a system composed of inert-gas 
crystals with different densities, light (ρ1) and 
heavy (ρ2). The density ratio is ρ1/ρ2 = 1/4. We 
considered the two cases of the shock propagation: 
from the light to heavy solid and from the heavy to 
light solid. The simulations in a system with solid-
vacuum interface correspond to experiments in 
which a metal with corrugations at the free surface 
is shocked from the opposite side. 
For solid-solid case, the atomic interaction is 
modeled via a modified Lennard-Jones (LJ) pair 
potential with an additional spline term to provide 
smooth second derivatives [10]. The potential 
parameters are used to define reduced MD units 
(mdu). For argon, length mdu σ = 3.405Å, energy 
437
Downloaded 21 Aug 2006 to 131.215.240.9. Redistribution subject to AIP license or copyright, see http://proceedings.aip.org/proceedings/cpcr.jsp
pi
st
on
shock
λ η
perturbation
at interface
up us
ρ2ρ1
pi
st
on
pi
st
on
pi
st
on
FIGURE 1. Simulation of a piston-driven planar shock
propagation through the perturbed interface separating
1) two solids of different densities (ρ1<ρ2), and 2) solid
and vacuum (ρ2=0). 
Up=0.7Up=1.0up us1
us2
reflected 
shock
ρ2
FIGURE 2. Shock wave passes through the light-to-
heavy solid interface (ρ2 = 4ρ1, η=λ/3=28 nm) generating
a slower second shock in the denser solid and a reflected
shock running backward into the light material.
mdu ε/kB=119.8K, time mdu= σ((ma/48)/ε))1/2 
=0.311ps, velocity mdu = 1.094 km/s, and pressure 
mdu is 41.9 MPa. The system is built from two LJ 
solids that differ only by the atomic mass: light 
(m1=ma) and heavy (m2=4ma). About 20 million 
atoms were used to build the solid-solid system 
with dimensions 108x15x487 nm. For MD 
simulations of the RM instability of metal-to-
vacuum surface we chose Mishin’s potential for 
copper [11] which has been previously used to 
study shocks in single crystal and polycrystalline 
copper [12-13]. The copper supercell has about 6 
million atoms and sizes of 74x5.1x218nm. Periodic 
boundaries were imposed in both transverse 
directions. 
Figure 1 shows a characteristic setup for 
modeling RMI in planar geometry. A piston wall 
represented by the repulsive potential moves at 
constant velocity up and generates a planar shock 
wave moving at the velocity us in the first solid. 
Then the shock waves passes into the second solid 
(or vacuum) through the interface with a single-
wave perturbation of the amplitude η and 
wavelength λ set by the largest transverse 
dimension (roughly, λ ~ 100 nm). The ratio η/λ in 
most simulations was chosen about 0.25 – 0.33. 
 
RESULTS AND DISCUSSION 
 
In the first series of runs, we have studied the 
planar shock moving from the light to heavy solid. 
Both solids have the [110] orientation along the 
shock direction. Figure 2 shows a snapshot when 
the first shock wave us1 hits the interface and splits 
into two non-planar waves: transmitted shock us2 
and reflected shock moving back into the light 
solid. After reaching the piston this wave reflects 
back and re-shocks the interface. Remarkably, the 
interface perturbation is shrinking when the shock 
passes through it, but then it bounces back and may 
start to grow. The growth rate strongly depends on 
the initial shock strength (see Fig. 3). For instance, 
there is almost no growth at the piston velocity up = 
0.7 mdu (approximately 0.54c1 , where c1 is the 
sound velocity in the light argon solid) until the 
transmitted shock emerges from an open surface of 
the heavy solid and the rarefaction wave moves 
back to the interface (see Fig. 3a). At higher 
velocity of the piston up = 1.0 mdu, the growth 
starts as soon as the shock passes through the 
interface, and then, when the rarefaction wave from 
the open right side reaches the interface, its rate 
increases. Such threshold behavior in solids is 
essentially different from that for RMI in liquids 
and indicates strong dependence of the growth 
conditions on the material properties (e.g. shear 
strength), generation of dislocations, defects, and 
plastic flow. There is a possibility of local melting 
in vicinity of the interface even though a chosen 
piston velocity does not reach the critical value at 
which the light solid could start to melt. This issue 
is currently under investigation and will be 
reported in a subsequent paper. 
In the second series of runs, we have studied the 
evolution of the interface for a shock wave moving 
from heavy to light solid. In this case, instead of 
reflected shock, a non-planar rarefaction wave 
moves backward into the first (heavy) material as 
shown in Fig. 4. The transmitted shock wave us1 
also becomes non-planar, and the front oscillates, 
but the front fluctuations decay soon. The 
438
Downloaded 21 Aug 2006 to 131.215.240.9. Redistribution subject to AIP license or copyright, see http://proceedings.aip.org/proceedings/cpcr.jsp
 
Up=0.5
Up=0.7
Up=0.6
Piston velocity
0 200 400 600 800 1000 1200
Time (mdu)
-80
-40
0
40
80
120
160
200
A
m
pl
itu
de
 (m
du
)
 
FIGURE 5. Evolution of the perturbation amplitude
during heavy-to-light shock passage for three piston
velocities: 0.5, 0.6, and 0.7 mdu; η=λ/3=36 nm. The
amplitude sign change indicates a phase inversion of the
perturbation 
 
Up=1.0
Up=0.7
(a)
0 200 400 600 800 1000 1200
Time (mdu)
80
120
160
200
240
A
m
pl
itu
de
 (m
du
) Piston velocity
Up=0.7
Up=1.0
(b)
 
FIGURE 3. The results for piston velocities 0.7 and 1.0
mdu, η=λ/3=36 nm: (a) snapshots of the perturbation at
the final stage (when passed [110] shock wave emerges
from a surface of the second solid and the rarefaction
wave runs back to the interface); (b) evolution of the
perturbation amplitude. 
us1
rarefaction
up=0.7
us1
ur2
ur2
local melting ? 
shock
phase 
inversion
us2 ρ1
us1
ρ2
rarefaction
ur2
 
FIGURE 4. Snapshots from the simulation of phase
change and growth of the perturbation during shock
propagation in heavy-to-light solid at the piston velocity
up=0.7 mdu; η=λ/3=36 nm. The spike arising from the
valley continues to grow until a rarefaction wave from
the open right side comes down to the piston.  
transmitted shock continues to run forward until it 
reaches the open surface at the right side. After the 
shock passage, the perturbation changes its phase, 
and the spike starts to grow from a hollow. 
However, for a shock wave with piston velocity 
up=0.5 we did not observe the spike growth after 
the hollow collapsed (see the dashed line 
approaching zero in Fig. 5). This indicates a 
threshold-type dynamics of the generation of the 
material flow around the perturbation. In contrast 
to the light-to-heavy system, in this case the 
growth-rate of the spike amplitude is slowing down 
to zero when the rarefaction wave from the open 
surface gets back to the piston. Multiple 
intersections of the transmitted non-planar shock 
waves (e.g. as seen in the second snapshot in Fig. 
4) create local areas with an increased stress and 
temperature around the growing spike (“hot spots”) 
where the material might become either molten or 
amorphous. 
In modeling the shock-induced instability of the 
perturbed copper surface we generate a shock wave 
along the [110] direction. A characteristic 
perturbation in the largest-size simulation (~6x106 
atoms) has a sinusoidal-like shape with the 
wavelength λ = 74 nm and the amplitude η = λ/4. 
The piston velocity is chosen in a range of 1 – 2.8 
km/s in order to keep it below the melting limit of 
439
Downloaded 21 Aug 2006 to 131.215.240.9. Redistribution subject to AIP license or copyright, see http://proceedings.aip.org/proceedings/cpcr.jsp
dislocations
phase 
inversion
(a) (b)
molten nuclei
incipient spall
FIGURE 6.  Two simulations of the RMI of the copper-vacuum interface with different perturbations parameters and piston
velocities: (a) [100] shock, λ=η=36 nm, up=1 km/s, us ≈5.5 km/s, Ps ≈50 GPa; (b) [110] shock, λ=4η=74 nm, up=2.8 km/s, us
≈ 8.2 km/s, Ps ≈ 200 GPa. In both cases, the phase inversion has been observed as a growth of the spike from the hollow.  
~3 km/s for the shocked bulk copper [12]. Figure 6 
shows two different examples of the instability 
growth when the shock wave reaches the open 
surface and the rarefaction wave moves backward 
into the bulk. In the first case (see Fig. 6a), we set 
the ratio λ/η = 1, which corresponds to a high 
degree of non-linearity even at the very early 
stages of instability. As a result, we were able to 
see the spike growth at the velocity of the piston 
up=1 km/s. Remarkably, the cup of the spike 
displays a high density of defects in comparison to 
dislocations in the bulk crystal, indicating possible 
amorphization of the material inside the spike and 
beneath which can facilitate the plastic flow around 
it and can sustain its growth. 
 In case of smaller perturbation amplitude η = 
λ/4 we observed the spike growth for two different 
piston velocities of 2.0 and 2.8 km/s. In the last 
case, with the bulk pressure close to the calculated 
shock melting limit of ~220 GPa [12], we noticed 
an appearance of local spots that might be nuclei of 
molten material. In addition, the void formation 
and nucleation with incipient spall have been 
detected underneath the hills of the initial 
perturbation during its phase inversion and spike 
growth at the unloading surface. Such dynamics is 
not observed in the simulations of shock unloading 
at a planar surface, and we relate it to the 
intersections of the non-planar rarefaction waves 
moving down from the perturbed interface. 
In summary, we have performed the first MD 
simulations of the Richtmyer-Meshkov instability 
in solids for different types of the interface and 
have studied various scenarios of the perturbation 
growth as well as its dependence on the shock 
orientation and strength. Such simulations provides 
unmatched insight into the nanoscale structure of 
the shocked material, and can determine the 
mechanism governing the instability growth in 
solids as well as lead to understanding of the 
material properties at a very high strain rates. 
 
ACKNOWLEDGEMENTS 
 
SVZ is supported by DOE ASC, ONR and ARO 
MURI. Work at LLNL was performed under the 
auspices of the U.S. Department of Energy by 
University of California, Lawrence Livermore 
National Laboratory under contract of No.W-7405-
Eng-48, LDRD 04-ERD-021. SIA is supported by 
NRL and by DOE ASC, contract #B523820. 
 
REFERENCES 
1. Barnes J.F., Blewett P.J., McQuenn R.G., Meyer 
K.A., and Venable D., J. Appl. Phys. 45, 727 (1974) 
2. Barnes J.F., Janney D.H., London R.K., Meyer 
K.A., and Sharp D.H., J. Appl. Phys. 51, 4678 
(1980). 
3. Colvin J.D., Legrand M., Remington B.A., and 
Weber S.V., J. Appl. Phys. 93, 5287-5301 (2003). 
4. Remington B.A. et al, Metal. Mater. Trans. A 35A, 
2587-2607 (2004). 
5. Plohr B.J., Sharp D.H. ZAMP, 49, 786 (1998). 
6. Wouchuk J.G. and Nishihara K., Phys. Plasmas. 4, 
1028 (1997). 
7. Terrones G., Phys. Rev. E, 71, 036306 (2005). 
8. Germann, T.C., Hammerberg, J.E., and Holian, B.L.  
in Shock Compression of Condensed Matter – 2003, 
edited by M. D. Furnish, Y. M. Gupta, and J. W. 
Forbes, AIP Conf. Proc. 706, 285-288 (2004). 
9. Kadau K., et al, PNAS, 101 (16), 5851-5855 (2004). 
10. Zhakhovskii V.V., Zybin S.V., Nishihara K., and 
Anisimov S.I., Phys. Rev. Lett. 83, 1175 (1999). 
11. Mishin et al,  Phys.Rev.B  63, 224106 (2003). 
440
Downloaded 21 Aug 2006 to 131.215.240.9. Redistribution subject to AIP license or copyright, see http://proceedings.aip.org/proceedings/cpcr.jsp
12. Bringa E.M., Cazamias J.U., Erhart P., Stolken J., 
Tanushev N., Wirth B.D., Rudd R.E., and Caturla 
M.J., J. Appl. Phys. 96, 3793 (2004). 
13. Erhart P., Bringa E.M., Kumar M., Albe K. Phys. 
Rev. B 72, 052104 (2005) 
441
Downloaded 21 Aug 2006 to 131.215.240.9. Redistribution subject to AIP license or copyright, see http://proceedings.aip.org/proceedings/cpcr.jsp


Molecular Dynamics Simulations of the Richtmyer-Meshkov Instability in Shock Loaded Solids

https://authors.library.caltech.edu/5044/1/ZYBaipcp06.pdf

Molecular Dynamics Simulations of the Richtmyer-Meshkov Instability in Shock Loaded Solids

Abstract

Similar works

Full text

Available Versions

Caltech Authors - Main