## Abstract

Atomistic methods, such as molecular dynamics and direct simulation Monte Carlo, constitute a powerful and growing set of techniques for fluid-dynamics simulation. The more fundamental nature of such methods, which exhibit nonlinear transport effects and small-scale fluctuations, extends their modelling accuracy to a significantly wider range of scales and regimes than the more traditional Navier–Stokes-based continuum fluid-simulation techniques. In this paper, we describe the current state of the art in atomistic fluid simulation, from both a theoretical and a computational standpoint, and outline the advantages and limitations of such methods. In addition, we present an overview of some recent atomistic-simulation results on fluid instabilities and on the physical scaling of atomistic techniques. Finally, we suggest possible avenues of future research in the field.

## 1. Introduction

Traditionally, the dynamics of fluids are described by a set of partial differential equations defining the fluid as a continuum (Drazin 2002). These equations—for example the Navier–Stokes equations—are solved numerically on a grid that defines the resolution of the solution. A pioneering example of continuum simulation for a wave problem was performed by Harlow *et al.* (1965) in the 1960s at Los Alamos National Laboratory. Since then, the range of applications has increased enormously to include complex flows in industrial applications and fluid instabilities such as the Kelvin–Helmholz (Kelvin 1871), Richtmyer–Meshkov (Richtmyer 1960; Meshkov 1969) and Rayleigh–Taylor (RT) (Rayleigh 1882; Taylor 1950) instabilities.

Molecular dynamics (MD; Alder & Wainwright 1957) and direct simulation Monte Carlo (DSMC; Bird 1994) describe fluids on the most fundamental level, i.e. that of atoms or molecules. The dynamics of the atoms or molecules are explicitly calculated, and fluid properties, such as the flow velocity or density field, can be calculated as averages over the trajectories of the particles. These methods are computationally expensive, and the computational effort increases linearly with the number of particles and the physical time scale simulated. Stationary hydrodynamic problems, such as Rayleigh–Benard convection or Taylor–Couette flow, have been investigated by the MD approach and successfully compared with continuum-based simulations, even with a relatively small number of particles (Puhl *et al.* 1989; Rapaport 1989, 2006; Hirshfeld & Rapaport 1998, 2000). For more challenging, non-stationary flows, such as the RT (Kadau *et al.* 2004) and Kelvin–Helmholtz (Glosli *et al.* 2007) instabilities, a larger number of particles is required to describe the full complexity of the mixing process. In addition, MD has been applied to a number of other problems in microscopic fluid dynamics, e.g. shock-driven turbulence (Root *et al.* 2003) and swimming by micro-organisms (Rapaport 2008).

What is the motivation for describing these problems by atomistic methods, given that the largest available supercomputer architectures are often required? Recently, it has been shown that there is a significant difference between the continuum and MD treatment of nanojets. In particular, fluctuations were found to be important in the breakup and stability of such jets (Moseler & Landman 2000). Fluctuations, such as those that stem from molecular disorder and collisions, arise in a natural manner in atomistic methods, but are usually lacking in continuum descriptions. For an instability problem, small fluctuations can be amplified by many orders of magnitude, and dynamic fluctuations—in particular at the interface—may be of importance. Recent large-scale MD and DSMC simulations of the RT instability show similarities and disagreements with continuum descriptions, and compare favourably to experiments. Furthermore, fluctuations were found to be important for the RT problem (Kadau *et al.* 2007).

As computational power is more than doubling every 18 months, larger and larger atomistic simulations can be performed. We have used the high-performance large-scale Scalable Parallel Short-range Molecular (SPaSM) dynamics molecular dynamics code, which recently performed the first trillion atom simulation on the BlueGene/L system (Germann & Kadau 2008), to simulate the RT instability with system sizes up to seven billion atoms (Kadau *et al.* 2007). In this review, we outline the results of these simulations, and present some of the theory behind atomistic simulation.

## 2. Theoretical background

Atomistic methods such as MD represent a framework for describing fluids at a particular level of detail. This framework fits into a hierarchy of theoretical and computational descriptions of fluids, each requiring its own set of approximations and having its own range of applicability. It is instructive to review the various levels on which the dynamics of a fluid may be characterized, and the relationships between them.

Consider a fluid of *N* interacting atoms or molecules (here referred to generically as ‘particles’) having positions and velocities *x*_{j} and *v*_{j}, respectively, where 1≤*j*≤*N*. The most fundamental level at which such a fluid (or any other material) can be described is governed by the time-dependent Schrödinger equation (TDSE),
2.1
Here, *ψ*(*x*_{1}…*x*_{N},*t*) is the wave function for the *N* particles, *V* (*x*_{1}…*x*_{N}) is the interaction potential, *m*_{j} is the mass of the *j*th particle and is Planck’s constant. Although equation (2.1) is, in principle, the most generally applicable description of a fluid, the direct numerical integration of equation (2.1) for more than a few particles is computationally prohibitive. For this reason, application of equation (2.1) in this context has historically been limited to the simulation of very small collections of atoms or molecules via some approximate or reduced version of 2.1 (Öhrn & Deumens 1999).

Fortunately, the full power and generality of equation (2.1) is not often required. In many cases, the size of the quantum-mechanical wave packet representing each particle is reasonably small compared with the average distance between particles. This condition can be characterized by a dimensionless number equal to the ratio of the de Broglie thermal wavelength of a particle to the average distance between particles (Hansen & McDonald 1986),
2.2
Here, *n* is the number density of particles, *T* is the temperature and *k*_{B} is Boltzmann’s constant. As examples, for liquid water at room temperature and liquid hydrogen (H_{2}), *B*≈0.07 and 0.97, respectively.

When *B* is reasonably small—as in the case of water, although not for hydrogen—we can neglect quantum effects in describing the atomic nuclei and pass to the classical limit (Shankar 1994), in which the variables {*x*_{j},*v*_{j}}—nominally random quantities under the TDSE—are replaced by their mean values, to yield the familiar Newtonian equations of motion for MD,
2.3
Since this is a system of ordinary differential equations, rather than a many-dimensional partial differential equation as in equation (2.1), equation (2.3) is comparatively straightforward to integrate numerically. It is for this reason that MD has become such a useful and popular material simulation method in many different contexts.

It is often the case that a fluid spanning macroscopic length scales must be considered. Hybrid methods (Garcia *et al.* 1999), in which the essential parts of the problem are modelled by MD and the remainder by continuum methods, may be capable of capturing such scales in certain contexts. However, such a scenario is currently beyond the reach of MD, due to the large number of particles present in a macroscopic fluid parcel. For example, a cubic metre of water at room temperature contains approximately 3×10^{28} particles, which is 16 orders of magnitude greater than in even the largest MD simulation to date (Germann & Kadau 2008). For this reason, it is often preferable to approximate a fluid as a continuum rather than as a collection of particles. Such a description can be derived directly from equation (2.3). Consider the probability density function *f*({*x*_{j}},{*v*_{j}},*t*) of the positions and velocities of the *N* particles. We define the continuum mass density and velocity fields of the fluid as
2.4and
2.5
respectively. Note that these fields are related to averages of the microscopic densities of various quantities. For this reason, the continuum approximation may just as accurately be referred to as a moment or mean-field approximation.

Via a sequence of mathematical steps known as the Irving & Kirkwood (1950) procedure, it can be shown that **u**(*x*,*t*) obeys
2.6
Here, *σ*=*σ*(*x*,*t*) is the stress tensor at *x*, for which an expression exists that is too complex to reproduce here. Although equation (2.6) exactly describes the evolution of the moment **u**(*x*,*t*) whenever MD is valid, it is not always an accurate description of the fluid in question. The use of equation (2.6) in this context requires that the actual behaviour of the fluid be the same as the dynamics of its mean. In order for equation (2.6) to be accurate at a particular length scale Δ*x*, the statistical fluctuations in the mean velocity field at that scale must be small compared to the average variation in **u**(*x*,*t*) across a distance Δ*x*. The ratio of these two quantities is given by the dimensionless number
2.7
where *S* is the magnitude of the shear rate in the flow under consideration and *ρ*_{0} is its average mass density. If *C* is not small, statistical fluctuations cannot be neglected, and the continuum mean-field approximation is invalid. For a flow of liquid water at room temperature with a shear rate of *S*≈1 *s*^{−1} (a typical macroscopic value in everyday experience), *C*≈0.01 when Δ*x*≈0.13 mm. A cubic volume of this size contains approximately 7×10^{16} water molecules. As an important aside, it is not necessarily always the case that a continuum description such as equation (2.6) can be substituted for MD when *C* is small. As noted in the introduction, it has recently been suggested that microscopic fluctuations may play an important role in fluid instabilities, even on the macroscopic scale (Kadau *et al.* 2007), and the possibility that thermal fluctuations may affect the development of turbulence has been known for decades (Betchov 1961; Hosokawa 1976).

Although equation (2.6)—along with its counterparts describing the evolution of *ρ*(*x*,*t*) and the energy density *e*(*x*,*t*)—is the most general continuum description of a fluid, it is of little use as it is written, since there exists no explicit closure for *σ* in terms of *ρ*, **u**, *e* and the pressure, *p*. In Chapman & Cowling (1960), it is shown that, in cases where the dimensionless *Knudsen number*
2.8
(where *λ* is the mean free path of particles in the fluid) is small, the general expression for the stress tensor *σ* reduces to
2.9
Here, *μ*_{S} and *μ*_{B} are the shear and bulk viscosities, respectively. For liquid water at room temperature, *Kn*≈0.01 when *L*≈10 nm. A cubic volume of this size contains approximately 30 000 water molecules. This expression for the stress tensor characterizes a Newtonian fluid, and results in the familiar compressible Navier–Stokes equation,
2.10
This equation is the most popular framework for the simulation of fluids in general practice.

A schematic Venn diagram is shown in figure 1, detailing the relationships in parameter space between the various fluid descriptions outlined here. Note that the validity of each approximate method depends on the accuracy of the method from which it is derived, and that there exists a dimensionless number (*B*, *C* or *Kn*) characterizing when each approximation is allowable (although the precise ranges of these numbers required for each description will vary widely depending on the problem). Note also that this diagram takes into account only theoretical applicability, and ignores practical concerns, such as the actual length scales or numbers of particles describable via each method, with today’s computational capacity.

The list of fluid descriptions outlined in this section is by no means exhaustive, as there exist numerous hybrid methods that do not fall neatly into such a schema. These include, for example, DSMC (Bird 1994), which is a method similar to MD, except that the interparticle interactions are modelled stochastically via random collisions in which momentum and kinetic energy are exchanged. Due to the simplified nature of the interactions involved, DSMC is generally 10–100 times faster than MD, depending on the fluid density. The main drawback of DSMC is the fact that it is restricted to the description of fluids with an ideal-gas equation of state. DSMC has traditionally been used primarily for the description of rarefied gases in aerospace and industrial applications. More recently, however, DSMC has been applied to more diverse problems in fluid dynamics, such as the Rayleigh–Taylor instability simulations reviewed in this work (Kadau *et al.* 2007), as well as to the Richtmyer–Meshkov instability (Barber *et al.* 2007), vortex shedding (Usami *et al.* 2008) and astrophysical problems (Zhang *et al.* 2003). Since DSMC is nominally derived from the Boltzmann equation, its region of validity straddles the boundary between MD and continuum methods. In addition to DSMC, lattice methods, which are particle-based schemes for the solution of the Boltzmann or Navier–Stokes equations (McNamara & Zanetti 1988), have also been used extensively in fluid-dynamics simulations.

## 3. Rayleigh–Taylor instability

The RT instability occurs when a heavy fluid with density *ρ*_{h} lies on top of a light fluid with density *ρ*_{l} in the presence of a gravitational field *g*. The two fluids subsequently interpenetrate in an archetypical example of turbulent fluid mixing. The incursions of the light fluid into the heavy fluid are generally referred to as bubbles, and those of the heavy fluid into the light fluid are known as spikes. By linearizing the Navier–Stokes equations, one can derive the exponential growth rate *n*(*k*) for modes on the interface as a function of their wave number *k* (Chandrasekhar 1961). This initial exponential growth is well described by both continuum and atomistic methods (Kadau *et al.* 2004; Barber *et al.* 2008; see also figure 2). Once the bubbles and spikes become larger and begin to interact and merge, the penetration depths *h*_{S} and *h*_{B} of the spikes and bubbles, respectively, are usually modelled as growing according to a free-fall type law, *h*_{S/B}(*t*)=*α*_{S/B}*Agt*^{2}, where the *α*_{S/B} are dimensionless coefficients and *A*=(*ρ*_{h}−*ρ*_{l})/(*ρ*_{h}+*ρ*_{l}) is the Atwood number, which characterizes the density difference between the fluids (Sharp 1984).

In a continuum simulation of the RT instability, the otherwise perfectly flat initial interface must be perturbed in order to become unstable, and various studies have examined the dependence of the mixing process on the initial spectrum of perturbations (Cook & Dimotakis 2001; Cabot & Cook 2006). In an atomistic simulation, on the other hand, one has the option to start with an almost perfectly flat interface, perturbed only on the length scale of an interatomic distance (figure 2). Furthermore, thermal fluctuations will excite all possible interfacial modes in a random fashion, and they will subsequently grow (on average) according to the continuum linear-stability prediction (Barber *et al.* 2008). After a short period of time, modes near the mode of maximum instability dominate (figure 2). Subsequently, the flow becomes more complex, and interaction and merging processes between neighbouring spikes/bubbles occur. It is in this stage that the *α*_{S/B} values described above can be measured. MD and DSMC simulations have shown that these values compare favourably with experiments, even though the macroscopic length scales of experiments and the nano- to micrometre scales of the simulations differ by many orders of magnitude (Kadau *et al.* 2004). However, an analysis of atomistic simulations at later stages of the mixing process reveals that the initial orderly merging process becomes increasingly unsynchronized, resulting in spikes that grow thinner and thinner. This leads to a process in which the spikes break up and detach from the main body of the heavy fluid to form droplets. These disconnection events, which have also been observed experimentally (Kadau *et al.* 2007), seem to dominate the front progression at later stages (figures 3 and 4). Most of the atomistic simulations were performed in a quasi two-dimensional thin-slab geometry. The thin out-of-plane dimension allows for three-dimensional transport properties in what is an effectively two-dimensional hydrodynamic flow. A few fully three-dimensional simulations have been performed that did not reveal any significant statistical differences in integrated quantities such as the *α*_{S/B} values or the observed break-up mechanisms. Future studies with larger fully three-dimensional geometries are needed to investigate the difference between two- and three-dimensional hydrodynamic flow.

The flow Reynolds number for an RT instability can be defined as , with *h*=*h*_{S}+*h*_{B} and *ν* the average kinematic viscosity of the two fluids. Fully three-dimensional (3072×3072×(256−3072) grid points, with the vertical grid expanding in time) direct numerical simulation (DNS) of the incompressible Navier–Stokes equation reached *Re*=32 000 on the BlueGene/L supercomputer platform (two weeks on 65536 CPUs; Cabot & Cook 2006). With the compressible atomistic method DSMC, we reached *Re*=12 000 for a quasi two-dimensional thin-slab geometry (1×32768×9856 cells) simulation performed on 65536 CPUs of BlueGene/L for two weeks (figure 3). This suggests about a factor of 50 more effort for DSMC when compared with DNS. The larger effort in DSMC simulations stems from the about 20 particles per cell and the smaller integration time step needed to capture fluctuations.

## 4. Scaling

As alluded to in the introduction, these simulations of the Rayleigh–Taylor instability and other atomistic simulations are confined by the current limits of computational capacity to the description of length scales of the order of microns. In order to ensure that the insights gained from such simulations are relevant on the centimetre or metre length scales of experiments or everyday experience, it is important to investigate the degree to which the results of atomistic computations physically scale. To illustrate what is meant by this notion, consider the Navier–Stokes equation with gravity,
4.1
Here, *ρ* is the mass density, **u**=**u**(**x**,*t*; *g*) is the velocity field of the fluid, *p*=*p*(**x**,*t*; *g*) is the pressure field, *μ* is the shear viscosity and *g* is the acceleration due to gravity. This is a canonical example of an equation with solutions that display various types of scaling. Of particular relevance is the fact that any solution {**u**(**x**,*t*;*g*),*p*(**x**,*t*;*g*)} to equation (4.1) remains a solution under the transformation , , , , , where *a*>0 is any constant. For example, the results of a fluid simulation accomplished via the numerical solution of equation (4.1) will just as accurately represent a system twice as large, provided that the material properties of the fluid (i.e. *ρ* and *μ*) remain the same, the time scale is increased by a factor of 4, gravity is decreased by a factor of 8, the velocity is decreased by a factor of 2 and the pressure is decreased by a factor of 4.

This type of physical scaling—as opposed to the very different notion of algorithmic or computational scaling—ensures that arbitrarily large scales are at least nominally amenable to description by equation (4.1). In the system of ordinary differential equations constituting MD, or in the set of rules describing DSMC, however, no such physical scaling is immediately apparent or rigorously derivable. It is therefore of great utility and interest to empirically check whether the dynamics of a fluid simulated via MD or DSMC displays the same scaling as that described above for the Navier–Stokes equation.

We have performed such an empirical check in the context of DSMC simulation of the RT instability (Kadau *et al.* 2008). The results from such a simulation can be rendered dimensionless through the use of two quantities: , the wavelength of maximum linear instability for perturbations on the interface and *τ* the exponential growth time of that most-unstable mode. We are aware of no closed-form expressions that exist for and *τ* in the general case. However, we can assume that these quantities depend solely upon *g*, *ρ*_{h} and *ρ*_{l}, as well as upon the shear viscosities *μ*_{h} and *μ*_{l} of the heavy and light fluids, respectively. (For a justification of this assumption, see Chandrasekhar (1961)). Under this assumption, it can be proven via dimensional analysis (Kadau *et al.* 2008) that
4.2and
4.3
where is the averaged shear viscosity and *f*_{1} and *f*_{2} are unknown dimensionless functions of the Atwood number *A* and the viscosity ratio *V* =*μ*_{h}/*μ*_{l}. The nature of this dependence on *g* is such that any two systems obeying the above physical scaling should become equivalent when non-dimensionalized in this way. Non-dimensionalized images from several DSMC simulations with different gravities are shown in figure 5 for each of several dimensionless times. Note that intensive quantities—such as the temperature, densities and viscosities—were held fixed between these various simulations, as were algorithmic constants such as the time step. In addition to the gravity *g*, only the size of the domain and the amount of time simulated were allowed to vary. Due to the stochastic nature of DSMC, we should not expect the images from these different simulations to be identical, even in the case of perfect physical scaling. Instead, the equivalence of the dimensionless results is manifested in the fact that the shape of the interface has the same qualitative character at the same dimensionless time in each simulation. This notion can be made more quantitative by considering the fractal dimension of the interface in the last frame shown for each gravity, which—within fluctuations—is approximately the same for all.

The various simulations shown in figure 5 span nearly an order of magnitude in length scale (and more than two orders of magnitude in gravity). However, another even larger simulation (not shown here) is still only 45 *μm* in width (figure 2). To establish whether this physical scaling extends to the macroscopic scale, we have compared the evolution of the spike penetration depth *h*_{S}(*t*) in a magnetic levitation RT experiment (Carlès *et al.* 2006) to that in the corresponding DSMC simulation. The dimensionless results are shown in figure 6. Despite the fact that the two systems differ by factors of approximately 42 million, 70 thousand and 10 billion in length scale, time scale and gravity, respectively, the two curves lie nearly on top of one another. This result demonstrates a strong scaling relationship for atomistic simulations of the RT instability across many orders of magnitude in size. Furthermore, since the RT instability is in no way exceptional among examples of turbulent fluid mixing, we expect this result to hold in a more general context, implying a profound connection between the results of microscopic MD/DSMC simulations and fluid behaviour on the scales with which we are most familiar.

## 5. Conclusion

We have shown that, with today’s computational capacity, it is possible to describe complex non-stationary hydrodynamic flows, such as the RT instability, by atomistic methods achieving Reynolds numbers up to 12 000. Remarkably, we have found that these simulations can be scaled up to describe larger systems, and can even be compared with macroscopic experiments many orders of magnitude larger than the atomistic simulations themselves. The importance of fluctuations in the description of instabilities has been explained (Kadau *et al.* 2007). An alternative approach to the use of atomistic methods is the inclusion of fluctuations into continuum descriptions (Landau & Lifshitz 1959). Preliminary work along these lines has been promising (Moseler & Landman 2000; Bell *et al.* 2007), and such techniques would allow macroscale simulations to include fluctuations, if necessary. Furthermore, as computational capacity continues to increase, greater and greater challenges in fluid dynamics will be accessible by atomistic simulations. These include high-Reynolds-number three-dimensional turbulence and complex mixing processes that include combustion. For example, atomistic methods including chemistry would be promising for the description of inertial confinement fusion (ICF; Rosen 1996; Lindl *et al.* 2004), an important potential source of clean energy production. This problem is not easily describable by traditional methods, as ICF involves a complex coupling between fluid instabilities and combustion processes at very small length and times scales.

## Acknowledgements

This work was carried out under the auspices of the National Nuclear Security Administration of the US Department of Energy at Los Alamos National Laboratory under contract no. DE-AC52-06NA25396 with funding from the ASC project.

## Footnotes

One contribution of 13 to a Theme Issue ‘Turbulent mixing and beyond’.

- © 2010 The Royal Society