A History of constitutive modeling via molecular dynamics : Shock waves in fluids and gases

From its inception in the mid-Fifties, the method of molecular-dynamics (MD) computer simulations has been used to probe the foundations of statistical mechanics, first for equilibrium equation-of-state averages, and then for transport properties from equilibrium fluctuations. Traditional statistical mechanical theoreticians were shocked to see that this new-fangled computational physics approach was feasible, even with incredibly tiny samples (on the order of a hundred atoms). When direct measurement of transport coefficients by non-equilibrium molecular dynamics (NEMD) was proposed in the early Seventies, even greater resistance was encountered from the traditionalists – though evidence for convergence with the equilibrium fluctuation method gradually accumulated. In the late Seventies and early Eighties, shock-wave simulations by NEMD made it possible to test directly the principal continuum constitutive theory for fluids, namely, Navier-Stokes viscous flow and Fourier’s Law of heat conduction. To everyone’s surprise – and the consternation of many – NEMD, once again, demonstrated that continuum theory applies at embarrassingly small (atomistic) time and length scales. We pursue this early line of work into the modern era, showing how NEMD shock-wave simulations can still provide surprising insights and improvements upon our understanding of constitutive modeling.


EPJ Web of Conferences
Metropolis, who was in charge of the Los Alamos computers) [1].As far as most of the theoretical statistical mechanics community was concerned, heavily invested in analytic methods for highly nonlinear problems, this was an amusing, but not necessarily useful development.The MC method (associated with Metropolis for alphabetical authorship reasons) followed the "evolution" of atomic "motions" in a Markovian random walk through a correlated ensemble of snapshots sampled from the canonical (NVT) ensemble: N atoms in a periodic-boundary-condition box of volume V and temperature T .If you had been able to make a movie of those snapshots, your eye would have been deceived into thinking that you were actually seeing real-life action at the atomistic scale.But in fact, the ensemble averages provided highly nonlinear equilibrium properties, namely, energy E(V,T ) and pressure P(V,T ) -the so-called "equation of state" (EOS) at equilibrium -given the interaction potential between atoms.No action, however, and therefore, no transport of mass, momentum, and energy that characterizes non-equilibrium dynamics.
The dynamics came only four years later at Edward Teller's Livermore Lab, when Berni Alder and Tom Wainwright [2] "invented" molecular-dynamics (MD) simulations of hard spheres, for which a great deal of theoretical work had been done.But the dynamics meant that the movies that A&W created at great effort and expense were real-life action movies -though maybe not great box office attractions.Again, MD was a spin-off of the Manhattan Project and the ensuing Cold War, and A&W took full advantage of their trusty computer programmer, Maryann Mansigh, who engineered the STEP program to squeeze out precious cpu seconds during idle moments when the workhorse hydrodynamics codes were taking a break.
The output of MD was at first limited to equilibrium time averages from the microcanonical (NVE) ensemble, appropriate to the time-dependent solution of highly nonlinear Newton's equations of motion that conserve energy E, rather than temperature T .Hard-sphere forces are purely impulsive, however, so that tables of potential collisions must be constructed and updated.As long as no compressive shock waves were simulated, some non-equilibrium problems could be solved, however.A&W carried out a hard-sphere simulation of the famous Boltzmann H-theorem, where a cluster of hard spheres was allowed to expand into a large volume, with a demonstrable increase in entropy as the velocity distribution function relaxes.Boltzmann himself would have carried out this fascinating MD demonstration a half-century earlier, if only he had had an electronic computer.(Then he could have checked out the Boltzmann equation and settled many a furious debate in his day.)Berni's demonstration was the highlight of a 1956 conference organized by Prigogine in Brussels [3].
In the next decade, Vineyard and co-workers at Brookhaven Lab, inspired by Berni's demonstration of the H-theorem, sought to apply MD to continuous potentials, with a continuous, secondorder, time-marching integration method (later credited to Verlet).But because Vineyard was interested in radiation damage in crystals, rather than equilibrium averages of the EOS in fluids, much of his team's work remained hidden to the traditional statistical mechanics community [4].There, the interest was in the Lennard-Jones (LJ) 6-12 potential, a model pairwise interaction studied in gas-phase problems.First, Rahman [5], and then Verlet [6], studied the LJ EOS, in 1964 and 1967, respectively.
It should be pointed out that these atomistic MC and MD simulations of equilibrium EOSs were met with enormous resistance from the traditional stat-mech practitioners, most of whom were trained in fancy mathematics, and who viewed these "computer jockeys" as charlatans, at best.Their objections were rooted in the notion that these equilibrium properties needed an enormous number N of atoms to be believably close to the thermodynamic limit, which they imagine had to be a significant fraction of a mole (N = 10 24 ).Therefore, when these "simulators" used 10-100 atoms, the "rigor mortis" school of mathematicians went ballistic.Moreover, the length of time used to obtain these time averages was unbelievably short -rather than seconds (or milliseconds, at least), they considered only 10-100 mean collision times per atom, or at most a few picoseconds (10 −12 s).Then, there was the question of time needed to establish ergodic behavior, which loosely translated meant that MC and MD ensemble differences might be miles from each other.
The latter question was settled by Bill Wood (and Jack Jacobson) at Los Alamos and A&W at Livermore for the hard-sphere melting line -a classic MC vs. MD showdown in back-to-back papers in the Journal of Chemical Physics [7].Needless to say, the answer was positive: the difference between the canonical and microcanonical ensembles (ensemble averages A vs. time averages Ā)

00002-p.2
New Models and Hydrocodes for Shock Wave Processes in Condensed Matter was controllably small: Thus, a system of 100 hard spheres could reproduce the infinite thermodynamic limit to about 1%, at least for equilibrium averages.There was some grumbling, but gradually, the critics were quieted.
Next, in the latter days of the 1960's, A&W addressed the problem of computing fluid transport coefficients for the hard-sphere fluid, in order to make contact with Navier-Stokes hydrodynamics [8].Green-Kubo (GK) theory suggested that equilibrium fluctuations relax on time-scales accessible to MD simulations, and A&W proved the point by computing transport coefficients (viscosity and thermal conductivity) from GK, where the long-time integral of flux autocorrelation functions give the transport coefficients [9].In the equations below, GK self-diffusion and shear-viscosity coefficients are presented, along with the Irving-Kirkwood atomistic expression for the shear component of the pressure tensor in the autocorrelation function for shear viscosity [10].
dt P xy (0) P xy (t) , shear viscosity With another seminal paper in 1970 [11], A&W showed that the self-diffusion coefficient of the hardsphere fluid away from the melting line is enhanced, over and above the level predicted from that predicted from the exponential decay of the velocity-velocity autocorrelation function (VACF) -the so-called "molecular chaos" assumption of Boltzmann, who threw up his hands at trying to follow the chain of collisions between an atom and its neighbors.(Remember, he had no computer.)Because Livermore had hydrocodes that could solve the problem of a fluid element set into motion, the A&W team (mainly Wainwright) set about examining the fluid-flow field that would result, namely, a double vortex in fluid velocity, whose impetus would serve to push the original element from behind, thus enhancing its memory of the original motion.They then plotted a time average of the velocity field in an MD simulation, to see if the same double vortex appeared.In Fig. 1, we reproduce the upper half of the double vortex from Alder and Wainwright's 1970 paper.This demonstration of hydrodynamic behavior at the atomistic scale -both time and distanceforever forced traditional hydrodynamicists to change their view of what is meant by "the hydrodynamic scale."No longer could it be said that extraordinary long times (milliseconds) and distances (millimeters) were required in order to see the fullness of hydrodynamic behavior.And the theoreticians quickly jumped in to verify -once they saw the MD results -that the double vortex indeed could be modeled with hydrodynamic mode-coupling theory [12], confirming Alder's intuition about the d-dimensional t −d/2 decay of the VACF, and calculating its coefficient quantitatively.
Once this demonstration had been made, to the consternation of traditionalists in both the hydrodynamic and statistical mechanics communities, one might have thought that the next obvious development would have been a slam dunk, namely, non-equilibrium molecular-dynamics (NEMD) simulations of steady-state fluid flows of momentum and heat.The fundamental idea of NEMD was to impose experimental driving, but at the atomistic scale.Shear flow, for example, could be obtained by driving two reservoir regions in opposite directions, with sample fluid in between.Then, like the real-world experimental Couette-flow viscometer, the shear viscosity could be measured in the sample as the steady shear stress divided by the shear rate (relative velocity divided by distance between the two reservoirs).In the early 70's, Bill Hoover and his Ph.D. student, Bill Ashurst, undertook the first systematic study of the transport properties (shear viscosity and thermal conductivity) of the LJ fluid, over a wide range of densities and temperatures [13].By extrapolation of their NEMD shear viscosity data to zero strain rate, they were able to predict the correct GK value for the LJ fluid.Unfortunately, previous GK calculations were in error by 30%, which simply fueled the skepticism of "real" theoreticians about the validity of NEMD, who claimed that the driving gradients were orders of magnitude higher than any experiment, and therefore could not possibly represent reality.Of course, the gradients were high, as they had to be, in order to detect a measurable signal in the midst of thermal noise, but the very fact that Hoover and Ashurst were able to show that Ree-Eyring quadratic behavior of shear viscosity as a function of strain rate [14] fit their NEMD data should have been sufficient proof of principle.Their main difficulty was that GK theory was officially "blessed" by theoreticians, despite all of its statistical difficulties and expense in calculation, as compared to the few NEMD simulations required to extrapolate to zero gradient.It took another decade to settle the long-standing dispute -in NEMD's favor [15].(The remaining 5% discrepancy was due an extrapolation based on mode-coupling theory, which predicted an erroneous square-root of strain-rate behavior of viscosity rather than the analytic Ree-Eyring.That was finally resolved in the next decade by very long NEMD simulations at very low strain rates, showing that Ree-Eyring is the correct asymptotic behavior [16].)"Pure" theoreticians also distrusted new, non-Newtonian equations of motion cooked up by NEMD practitioners in the late 70's to simulate boundary-less driving [17], including for bulk viscosity [18].(Ashurst's rigid-boundary reservoirs introduced system-size dependencies, which necessitated extrapolations to infinite system size, but system-size dependence was also a difficulty for GK equilibrium simulations.) In hindsight, my own enlightenment about how to do NEMD more easily, without boundary reservoir regions, occurred twenty years too late to settle these issues [19].By invoking Navier-Stokes (NS) hydrodynamic theory, one can carry out shearing simulations in a periodic-boundary-condition (pbc) system that splits the computational cell into two regions, the left side driven to a constant average downward velocity, and the right side driven upward to the opposite velocity -on average in the region.In other words, the reservoir regions are the entire computational cell, with no artificial reflecting walls, just an average time-dependent force applied to atoms that find themselves on one side or the other.Temperature can be fixed over the whole system by Ashurst's velocity rescaling at each time step.The resulting velocity distribution is sinusoidal (as predicted by NS), and the shear viscosity can be computed from the time-average force applied to the two sides.The results are now obtainable on laptop computers for 1000 particles for a handful of relatively short simulations, and extrapolated to zero strain rate (the GK limit) by an analytic quadratic (i.e., Ree-Eyring) fit.The application to thermal transport is obvious: One side of the pbc system can be fixed at a cooler average temperature, and the other at a hotter temperature, on average; the resulting temperature profile is predicted by Fourier's

00002-p.4
New Models and Hydrocodes for Shock Wave Processes in Condensed Matter Fig. 2. Schematic of boundary-less periodic-boundary-condition (pbc) shear flow NEMD simulations.The Navier-Stokes solution to steady velocity driving (on average) in the two halves of the pbc system is a sinusoidal velocity distribution.A similar approach can be used for heat flow, where one side is (on average) colder and the other hotter, with a sinusoidal temperature field predicted by Fourier's Law.
Law to be sinusoidal, and the thermal conductivity can be determined from the heat injected into the hot side and extracted from the cold side [19].

Atomistic shockwave simulations: Under fire
The resistance to NEMD for steady-state transport processes in fluids was echoed in the skepticism that greeted early shockwave simulations at the atomistic scale.Again, criticism was leveled at the tiny time and distance scales, though the shock process can be very prompt and gradients at the shock front can be very steep -on the order of a few atomic spacings in strong shocks.The notion that continuum theories could be applied to shock waves was widely held to be absurd, since the steep gradients in the five hydrodynamic fields (density, three components of fluid velocity, and energy) would be far beyond the scope of the linear Navier-Stokes treatment of viscous fluid flow and Fourier's Law of heat flow.By "linear," Navier (in 1822), Stokes (in 1842), and Fourier (in 1822) meant that the fluxes of mass, momentum, and energy should be proportional to the corresponding gradients of velocity (or density) and temperature.Almost two centuries later, most theoreticians believed that shock waves would surely require higher-order gradients and/or products of gradients, and therefore that NS (&F) would be inadequate.
In 1978 in the Soviet Union, Klimenko and Dremin performed MD simulations of two strengths of shock waves (weak and moderate) in the Lennard-Jones fluid [20].A year later, believing that NEMD had already shown that hydrodynamics applies at the shockingly small atomistic scale, Bill Hoover at Livermore boldly ignored conventional wisdom and solved Navier-Stokes and Fourier continuum theory for both of the shock waves studied by the Russians [21].The comparison with MD was almost quantitative, which led us at Los Alamos in 1980 to push the limit further by performing an MD simulation of a very strong shock wave [22].(See Fig. 4 for the shockwave density and temperature trajectories of these three shock waves in the phase diagram of the LJ system.)Fig. 3. Strain-rate dependence of NEMD boundary-less shear viscosity (kinematic); the GK limit is predicted to within the error bars of the two approaches, but at considerably less computational expense for NEMD.The shock profiles of density and temperature as functions of distance for the strong LJ fluid shock are shown in Fig. 5.The NS theory is compared to MD and is shown to be too steep.One can therefore conclude that more dissipation is needed in the transport coefficients (longitudinal viscosity and thermal conductivity), in order for a continuum theory to describe adequately the atomistic shockwave results.

00002-p.6
New Models and Hydrocodes for Shock Wave Processes in Condensed Matter The unique feature of shock waves at the atomistic scale, which was discovered by NEMD simulations, is that the velocity distribution at the shock front is very anisotropic -shaped like a Cuban cigar -in contradistinction to the equilibrium (spherical) Maxwell-Boltzmann distribution.Thus, the kinetic temperature in the x-direction of the shock T xx is higher than the temperature in the transverse directions, as shown in Fig. 6.
Not only is T xx higher than the average temperature T (one-third the trace of the kinetic temperature tensor), but it exhibits a distinct peak in the midst of the shock front -a feature that is completely beyond the scope of Navier-Stokes, which utilizes density and (average) temperature for the dependence of the EOS and transport coefficients, and makes no provision for this completely non-equilibrium feature.Of course, the pressure tensor is not isotropic in NS theory, but no explicit consideration is given to the kinetic, as opposed to potential contribution to the pressure tensor.

Yet "atomistic" surely isn't the continuum, is it?
For NS theory, the pressure tensor P αβ under uniaxial compression (the volumetric strain ε is negative, as is the strain rate dε/dt) is given by: The fluid velocity in the x-direction is u; η S is the shear viscosity; η V is the bulk viscosity; η L is the longitudinal viscosity; P eq is the equilibrium pressure from the EOS; τ is the shear stress.In NS theory, the heat-flux vector Q is given by Fourier's Law (κ is the thermal conductivity): For a planar shock wave driven by an infinitely massive piston, the geometry can be represented in three frames of reference: (1) The piston can move to the left at velocity −u p , with the shock moving faster to the left at −u s into cold material at rest (fluid velocity u = 0, density ρ 0 , pressure P 0 , energy E 0 , temperature T 0 ); the hot material near the piston is at the final shocked state (fluid velocity u = −u p , density ρ 1 , pressure P 1 , energy E 1 , temperature T 1 ).( 2) Adding u p to all velocities in picture #1, the piston is then at rest, while cold material rushes toward it at velocity u p , and the shock wave moves to the left at u p − u s .(3) By adding u s to all velocities in picture #1, the shock front is now stationary in this frame of reference, where cold material rushes to the right at u s , and the piston recedes to the right at velocity u s − u p .(See Fig. 7.)

00002-p.7 EPJ Web of Conferences
In the stationary frame of reference for the planar shock wave, the fluxes of mass, momentum, and energy are constant, since the partial time derivatives of density, pressure, and energy per unit mass are zero and the gradients of the fluxes from the hydrodynamic equations are therefore also zero.
As you can see, the EOS and transport properties are functions of density and (average) temperature; there is no provision whatsoever in NS theory for saying what T xx is.We are therefore forced to postulate a relationship between the anisotropic pressure tensor and the kinetic temperature tensor.Holian, Mareschal, and Ravelo [23] proposed to use the so-called "thermodynamic equation of state" from first-year thermodynamics textbooks to solve the riddle: We suppose that P xx can be substituted for P, and T xx for T , under the assumption that the second derivatives of the free energy -(∂P/∂T ) V and (∂E/∂V) T -are isotropic, with possible small volumetric strain-dependent corrections.This ansatz works exactly for the ideal gas, where there are no potential contributions at all, either to the energy per unit mass, or to the pressure tensor.Hence, we can extend this to cover the general case of a dense fluid, thereby determining the difference between anisotropic (shock-direction) Fig. 8. Shockwave profiles of temperature (shock direction T xx , average T ) for the ideal gas [24]; NEMD, circles [25]; NS, dotted lines; use of T xx instead of T -dependence in transport coefficients, solid lines.The x-axis is scaled by the cold mean free path (to the left of the shock front), which is four times smaller at the hot end (to the right).
temperature and the isotropic average temperature from the pressure-tensor: The result is that there is a definite peak in T xx .In the ideal-gas case, T xx is larger than the final (hot) temperature T 1 by a factor of 4/3, and this occurs at fluid velocity u = u s /2.The ideal-gas temperature profiles are shown in Fig. 8 [24].
While it is clear that NS does a remarkably good job of modeling the NEMD results of Salomons and Mareschal [25], the agreement resembles that of the strong shock in the LJ fluid, namely, there appears to be too little dissipation, so that the NEMD profiles are broader.An indication of this is that by substituting T xx for T in the transport coefficients, thereby increasing them, the agreement is noticeably better, though not quite good enough [24].Last year, Holian and Mareschal revisited this problem and proposed two improvements in the heat-flux vector, over and above Fourier's Law: (1) a relaxation term, motivated by the Cattaneo (Maxwell-relaxation) Equation, and (2) a nonlinear strain-rate enhancement of the thermal conductivity, akin to the most important term in the Burnett higher-order treatment (beyond NS&F theory, which is linear in the velocity and temperature gradients) [26].The Cattaneo relaxation term suggests an anisotropic version of the heat-flux vector, compared to the usual isotropic (Fourier) Q 0 , which depends only on the average T .This brings into play the difference between P xx and P, and therefore, Fig. 9. Kinetic temperature profiles (T xx , solid lines; T , dashed lines) for the ideal gas, in units of mu 2 s /k: NS (blue), NSx (NS, with T -dependence in η L and κ replaced by T xx , green), Holian-Mareschal (HM) theory with relaxation term only (red), HM with both relaxation and Burnett terms (black), NEMD (circles) [26].The x-axis is in units of the mean free path in the cold gas, which becomes four times smaller at the hot end.
the by-now familiar difference between T xx and T : The second term to be considered is the nonlinear enhancement of the thermal conductivity: The strain rate du/dx is multiplied by a relaxation (collision) time τ c that depends on density and temperature: τ c = η L /B s , where B S is the isentropic bulk modulus; δ 1 is a dimensionless, order unity constant (if zero, the usual Fourier conductivity is recovered).The resulting heat-flux vector is (δ 2 is another dimensionless, order unity constant, which if zero, also returns Fourier's Law), according to [26]: In Fig. 9, the ideal-gas kinetic temperature profiles are presented.Following that, the heat-flux vector is shown in Fig. 10.Clear improvement is seen with the Holian-Mareschal treatment of the heat-flux vector over NS (&F).Both the shock thickness and the maximum in the heat flux (negative, since heat flows from the hot r.h.s. to the cold l.h.s.) are better reproduced by HM, compared to NS.Nevertheless, we see that the NEMD (reality) takes a noticeably longer time to relax to the final hot state than the theoretical models (either NS or HM).Finally, we show the HM results for the LJ dense-fluid strong shockwave case.The results are strikingly similar to the ideal-gas case, but the potential contributions to the pressure tensor are no longer trivial, as they are (zero) in the ideal gas.This is strong confirmation of the thermodynamic equation of state as applied to the tensors in the shock wave, and also a non-trivial verification of the HM heat-flux vector.[23].Here, the x-axis is in units of the zero-temperature collision diameter of the LJ dimer, or about five times the final (hot) mean free path.

Conclusion: The battle isn't quite over
What can we conclude from these rather old-fashioned constitutive-modeling results, so much more mundane than the fascinating problems of shock response in solids?First of all, the transfer of heat has been shown to be an important component of continuum modeling.Second, we can characterize the highly anisotropic velocity distribution in shock waves (and therefore, the temperature tensor) in terms 00002-p.12New Models and Hydrocodes for Shock Wave Processes in Condensed Matter of the pressure tensor, which is also highly anisotropic.From these first two conclusions comes a new way to look at the heat flow in shock waves.Third, we can say that the battle between the atomistic and continuum viewpoints, now a half century old (or a full century, if you hearken back to poor old Boltzmann), still rages, but only in isolated pockets of the physics community.Molecular dynamics has attained a level of wide acceptance -witness this plenary presentation in a hydrodynamics conference -that does the founding fathers of MD proud.And for that, for the pioneering work of those who went before me, for my co-workers through the years, and for the honor bestowed on me by the organizers of the conference in Paris, I say in my less-than-perfect French, "Merci beaucoup, mes bonnes amis.Comme toujours, la cuisine franc ¸aise et les vins sont superbes, tout comme nos amitiés durables.Et j'espère que ma présentation était cohérent -au moins en anglais." 1

Fig. 1 .
Fig. 1.Demonstration of the double-vortex field of fluid velocity in an equilibrium MD simulation of the hardsphere fluid [11].The long-time tail arises from hydrodynamic flow at the atomistic scale.

Fig. 5 .
Fig. 5. Density and temperature profiles for strong shock (two-fold compression and final shock temperature almost one-hundred times larger than the ambient melting temperature) in LJ fluid [22].

Fig. 6 .
Fig. 6.Kinetic temperature profiles in a strong shockwave in the LJ fluid [22]; T xx is measured in the direction of propagation, (T yy + T zz )/2 in the transverse direction (the difference is labeled T xx − T yy for convenience).

Fig. 11 .
Fig.11.Kinetic temperature profiles for the LJ dense-fluid strong shock wave, in units of ε/k (ε is the LJ dimer bond energy)[23].Here, the x-axis is in units of the zero-temperature collision diameter of the LJ dimer, or about five times the final (hot) mean free path.