Molecular Dynamics Simulation on Stability of Converging Shocks

Last 40 years MD method was successfully applied to many problems in different fields such as statistical and chemical physics. With increasing of the computer power, simulations of many millions of atoms in 3D space are now feasible. We believe that nowadays is time to apply MD method to fluid hydrodynamic problems.

MD approach is based on tracking of the atom motions and hence it can provide us whole information of a system behavior including phase transition in the matter, diffusion, thermoconductivity and other nonequilibrium processes far from thermodynamic steady state. Due to direct computation of the atom movements and the atom-atom interactions MD method has an fundamental advantages over numerical hydrodynamic methods which assumes shocks as a structureless discontinuity and requires an equation of state. Moreover, the different hydrodynamic codes show significant deviations in the simulation flow and a higher mesh resolution can not improve the situation [2], but the various MD programs give the same simulation results.

Fluid simulations of the Richtmyer-Meshkov instability by using numerical hydrodynamic methods have so far failed to provide quantitative agreement with experiments due to grid induced numerical instability. For this reason, it is necessary to apply an alternate method, namely MD approach as a very promising method of hydrodynamic simulation. It is the purpose of this work to carry out such numerical experiments. Simulations were performed by own full vectorized MD code -Molecular Dynamics Solver ( 11 seconds per one time step for 2 21 LJ atoms, and 92 ns per a pair of atoms on the NEC SX-5 supercomputer ).

In our MD simulation the atoms of a liquid interact via modified Lennard-Jones short-range pair potential with cut-off and smoothing at r cut = 2.5σ [3]:

where r 0 = 2 1/6 σ , and σ , ǫ are the usual LJ parameter, and a 2 = -3.5289 × 10 -3 , a 3 = 5.75868 × 10 -4 . We use these parameters and the atomic mass m a /48 as reduced molecular dynamic units. Here and after all quantities are in MD units (mdu). For argon atoms σ = 3.405 Ȧ, ǫ/k B = 119.8 K, time mdu = σ m a /48ǫ = 3.113 × 10 -13 s, velocity mdu = 1094 m/s, and the unit of pressure is 0.0419 GPa . At the beginning atoms are placed into rectangular MD cell inside cylindrical domain surrounded by a piston as a cylinder wall. The piston is simulated by an external potential ∼ [r -R(φ, t)] 2 and its position depends on angles in case of a perturbed boundary as well as on time to generate shock waves. The total number of atoms in simulation is 2097148 . The computational MD cell has dimensions L x × L y × L z , where L x = L y = 549.7 and the thickness of MD cell is L z = 21.8 , and periodical boundary conditions are imposed on the system along the cylinder axis (z-axis).

At the preparation stage the cylinder radius is fixed R 0 = 197.25 and all atoms have the same atomic weight m a = 48 . The Langevin thermostat is applied to prepare initial cold LJ liquid at the thermodynamic equilibrium with given temperature T = 0.72 and number density n = 0.79 . In case of a perturbed cylinder boundary the initial position of the piston is defined as R(φ) = R 0 [1 + δ sin(mφ)] where m is the mode number, and δ = 0.05 is the relative perturbation, ∆ = 2δR 0 is the initial perturbation amplitude.

In the end of the preparation stage an equilibrium state with uniform mass density is reached. At the beginning of the simulation stage in order to simulate two different materials the atoms outside of the interface change instantly its mass to heavy one m b = 16m a and its velocities to v b = v a /4 to hold uniform temperature and pressure. Interaction between light-heavy atoms Typical snapshots of MD system in the case of the perturbed converging shock are shown in Figure 1. It is found that for enough large initial perturbation the curved shock front generates Mach stems at very early stage of converging. Then the first generation of Mach stems entails the secondary Mach stems and so on. In fact polygons converge to the center instead of smooth curved shock front. There are m -side polygons for mode m or 2m -sides at phase inversion stage. Also the Mach stems generate intricate structures of reflected shocks in downstream flow which result in supersonic jets running away from the target boundary. We suppose these jets can amplify instability of ablation front in the real ICF experiments. Figure 2 shows amplitudes of the measured shock front perturbation as a function of a mean shock radius and the perturbation growth rate for different mode n

…(Full text truncated)…