Large-Scale Molecular Dynamics Simulations of Normal Shock Waves


Large-scale computational resources now enable the direct simulation of atomic systems, using realistic interatomic potentials. Molecular Dynamics (MD) simulations do not require any a priori model for the fluid, such as its equation of state, the transport coefficients, or a collision model. The interaction between the atoms is derived by some potential energy function, which is the only modeling effort required. Rapid advancements in the field of computational chemistry combined with continued advances in high performance computing may enable MD to be a very powerful technique to precisely investigate the highly complex nonequilibrium physics, found in high-speed flows, at the molecular level. Such MD simulations could then be used to inform DSMC collision models (a bottom-up approach), as opposed to developing them to match, in the continuum limit, the empirical laws used in CFD computations (a top-down approach). In general, if MD simulations are used to build DSMC collision models, for example, by providing the numerical value for certain parameters, or to determine a functional relation, then DSMC should reproduce the MD results very accurately, and more efficiently. 

Results and Publications

Molecular Dynamics Simulations of Noble Gas Mixtures:

    Valentini, P. and Schwartzentruber, T.E., “Large-scale Molecular Dynamics Simulations of

    Normal Shock Waves in Dilute Argon”, Physics of Fluids, 21, 066101, pp. 1-9, 2009.

    Valentini, P. and Schwartzentruber, T.E., “A Combined Event-Driven/Time-Driven Molecular

    Dynamics Algorithm for the Simulation of Shock Waves in Rarefied Gases”, Journal of

    Computational Physics, 228 (2009) 8766-8778.

    Valentini, P., Tump, P. A., Zhang, C., and Schwartzentruber, T.E., “Molecular dynamics s

    simulations of shock waves in mixtures of noble gases”, Journal of Thermophysics and Heat

    Transfer, Vol. 27, No. 2 (2013), pp. 226-234.

Starting in 2009, Dr. Valentini has pioneered the simulation of nonequilibrium dilute gases using Molecular Dynamics (MD). Research began with shock waves in a simple argon gas (monatomic) and was then extended to monatomic gas mixtures (argon-helium and xenon-helium). A number of experimental data sets exist for flow properties within such shock waves and thus MD simulations could be carefully validated by experimental data. This research demonstrated that pure MD simulations of such shock waves are computationally feasible (~100 cores for a few days), and are able to precisely resolve velocity distribution functions even for mixtures involving low concentrations of certain species. In addition, we developed a novel accelerated MD technique for dilute gases, called the “Event-Driven/Time-Driven” (EDTD) MD algorithm. The EDTD MD algorithm essentially advances molecules directly to their next impending collision (skipping the lengthy free-flight portion of their trajectory) while still integrating each collision with standard Time-Driven MD. The technique also simulates multi-body collisions. As a result, EDTD MD simulations of these shock waves could be performed on a single processor in a few days (50-100 times faster than pure MD). The speedup is directly proportional to the ratio of the mean-collision-time to the timestep required for pure MD simulation (~1 femtosecond).

A New DSMC and CFD Model for Rotational Relaxation (from Molecular Dynamics):

    Valentini, P., Zhang, C., and Schwartzentruber, T.E., “Molecular dynamics simulation of

    rotational relaxation in nitrogen: implications for rotational collision number models”, Physics of

    Fluids, Vol. 24 (2012), pp. 106101-1 – 106101-23.

    Zhang, C., Valentini, P., and Schwartzentruber, T.E., “Nonequilibrium-Direction-Dependent    

    Rotational Energy Model for use in Continuum and Stochastic Molecular Simulation”, AIAA

    Journal, Vol. 52, No. 3 (2014), pp. 604-617.

We have extended our molecular dynamics simulations to diatomic nitrogen for which experimental data on shock wave profiles also exists. At the low temperatures under which the experiments were performed there is no vibrational excitation of the gas, and thus the gas is undergoing translational-rotational relaxation only. The experimental data set of Robben and Talbot includes both density and rotational temperature profiles through the shocks as well as rotational energy distribution functions throughout the shock. Sample results comparing MD solutions to this experimental data are shown in the figures below.

After validating our interatomic potential and MD simulation technique with experimental data, we proceeded to investigate the dominant factors that influence the rotational energy relaxation rate in nitrogen. While prior models parameterized the relaxation collision number as a function of translational temperature only, Zrot(Ttr), our MD simulations clearly showed the collision number to be a strong function of both translational and rotational temperatures, Zrot(Ttr,Trot). More specifically, the relaxation rate was dependent on the direction to the equilibrium state (and magnitude of the deviation from equilibrium). For example, compressing flows (where Ttr > Trot) exhibit fast rotational excitation, whereas expanding flows (where Trot > Ttr) have slower de-excitation, for the same equilibrium temperature. This can be seen in the figures below.

Furthermore, in the limit of thermal equilibrium (Ttr = Trot = T), the collision number predicted by MD simulation, Z(T) extracted from the line in figure above-left, has a much weaker temperature dependence than the widely used Parker model. A variety of results, in this limit of thermal equilibrium, are plotted in the figure below.

Figure (above) depicts species bulk velocities in the flow direction within an argon-helium shock wave. Symbols are experimental data (argon=triangles, helium=circles), and solid lines are results from pure MD simulations.

Figure (above) depicts translational temperatures (Tx-triangles, Ty=Tz-circles) for argon within an argon-helium shock wave. Again, symbols are experimental data (here, argon=circles and helium=triangles), and solid lines are pure MD solutions.

Based on these MD results, we developed a new DSMC model, along with a consistent CFD model, called the Nonequilibrium-Direction-Dependent (NDD) model. Both the probability expression for use in DSMC and the collision number expression required for CFD are shown below. The model is easy to implement, computationally efficient, and most importantly is consistent with available experimental data for shock waves and matches pure MD solutions for a range of compressing and expanding flows. Full details can be found in the above papers.

Rotation-Vibration Coupling and Dissociation:

    Valentini, P., Norman, P.E., Zhang, C., and Schwartzentruber, T.E, “Analysis of rovibrational

    relaxation in nitrogen via direct atomic simulation”, AIAA Paper 2014-1079, Jan. 2014, presented at

    the 52nd Aerospace Sciences Meeting, AIAA SciTech, National Harbor, Maryland.

Most recently, we have extended our MD simulations to nitrogen including rotation, vibration, and dissociation reactions. This work is on-going, however, the general findings so far are: (1) vibrational relaxation does not exhibit the direction-dependent behavior observed in rotation; likely because vibrational energy levels become more-closely spaced as energy is increased, (2) at moderate temperatures (below 10,000 K) the vibrational relaxation rate predicted by MD (using the Ling-Rigby potential) is in good agreement with the experimental data of Millikan & White and rotational and vibrational relaxation processes are effectively decoupled, (3) at higher temperatures, we find that rotational and vibrational relaxation processes are coupled and that an additional relaxation time constant is required to reproduce pure MD results. This is evident from the isothermal relaxation calculations shown in the figure below (left), where direct rotational-vibrational energy transfer is evident as the translational temperature remains constant. The existing Landau-Teller (L-T) and Jeans equations used in CFD calculations cannot account for such energy transfer. The figure below (right) shows the resulting time constants predicted by MD simulations, where pink is the rotation time constant, orange is vibration time constant, and blue is coupled rotation-vibration time constant. These are compared to the Millikan & White models (black lines).

Finally, by using a Morse-type potential (allowing for bond-breaking and dissociation), we have simulated shock waves including the full dissociation of nitrogen using pure Molecular Dynamics. An example MD solution is shown in the figure (image below). We are continuing such simulations with a variety of interatomic potentials to investigate rotation-vibration coupling and dissociation for high-temperature nonequilibrium air flows.