New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Nanohydrodynamics simulations: An atomistic view of the Rayleigh–Taylor instability

Contributed by Berni J. Alder, February 23, 2004
Abstract
Nanohydrodynamics simulations, hydrodynamics on the nanometer and nanosecond scale by molecular dynamics simulations for up to 100 million particles, are performed on the latest generation of supercomputers. Such simulations exhibit Rayleigh–Taylor instability, the mixing of a heavy fluid on top of a light in the presence of a gravitational field, initiated by thermal fluctuations at the interface, leading to the chaotic regime in the longtime evolution of the mixing process. The earlytime behavior is in general agreement with linear analysis of continuum theory (Navier–Stokes), and the latetime behavior agrees quantitatively with experimental observations. Nanohydrodynamics provides insights into the turbulent mixing process that are inaccessible to either continuum calculations or to experiment.
In this article it is shown that one can make quantitative investigations of complex chaotic hydrodynamic flows, starting at the most fundamental level of atomic motion by using the latest generation of supercomputers. To do so requires solving the simultaneous Newtonian equations of motion for an enormous number of classical particles, largescale molecular dynamics (MD), which in effect solves the Liouville equation. Such hydrodynamical simulations by MD are limited to less than a billion atoms for a nanosecond, which means systems in three dimensions are, at most, 1,000 particles wide (hundreds of nanometers); hence, one can characterize hydrodynamic flows at these time and distance scales as “nanohydrodynamics.” Recent history in the development of massively parallel computers suggests that MD will not reach the scale of micrometers and microseconds in the foreseeable future (1).
Rayleigh–Taylor (RT) Instability
The RT instability occurs when a heavy fluid lies on top of a light fluid in the attendance of a gravitational field g; the fluids will subsequently mix in a moreorless turbulent process. The RT process is the classical example of turbulence and was first investigated by Lord Rayleigh in 1883 (2) and later theoretically investigated by Taylor (3). Its relevance ranges from astrophysical supernova explosions to geophysical formations like salt domes and volcanic islands (4, 5), all the way down to inertial confinement fusion (6) and to the general turbulent mixing of fluids.
For these purposes MD simulations have been performed on the well studied RT instability (7) as an example and validation of nanohydrodynamics. RT simulations require a large number of particles in an MD calculation (Fig. 1) to resolve fluid structures at late times and the imposition of an enormous gravitational field to reduce the most unstable wavelength to the nanoscale and to mix the fluids on the nanosecond timescale. This most unstable mode of wavelength λ (wavenumber, k = 2π/λ) must be several times smaller than the width of the simulation cell to develop several bubbles and spikes. From linear stability analysis, the most unstable wavelength decreases with increasing g (8), and thus a value of g ≈ 1 × 10^{10} g _{Earth} is required. Whether such a large g value distorts the instability process depends on whether g is a scalable variable up to that magnitude, which can be validated by comparison with experimental and continuum descriptions of the hydrodynamical problem as described by the Navier–Stokes (NS) equations, which assume scalability (4, 6, 8–14). Furthermore, supernovae explosions with g » g _{Earth} suggest scalability, and although most simulations were performed with that large value of g, an order of magnitude smaller value did not change the results within statistical fluctuations. The large influence of fluctuations, which have been minimized in the initial conditions by using an atomically flat interface, might lead to different results when compared with continuum descriptions where they are ignored. Here, the results of these nanohydrodynamical calculations compare quantitatively with macroscopic experiments (5, 15), and we compare them with theoretical continuum approaches.
Atomistic Simulations
From the fundamental Liouville equation (16) that governs the evolution of an ensemble of a conservative Hamiltonian system, one can derive the usual hydrodynamic NS equations by neglecting physical nonlinearities and thermal fluctuations. Solving the Liouville equation by MD has a further advantage over solving the continuum NS equations, because one can impose much more realistic boundary conditions, as in the onset of the RT instability (2, 3), from a surface roughness generated by natural thermal fluctuations. Moreover, in MD simulations, one can control the interactions between atoms, for example, making fluids perfectly miscible or immiscible at the flick of a switch. Numerically, the advantages of MD are also significant, because MD is a gridless method that is unconditionally stable, with no continuum interfacetracking algorithms and their attendant approximations, for example. The disadvantage of MD is that it is computerintensive, but that disadvantage is mitigated because the number of continuum grid points required to spatially resolve complex hydrodynamic flows is similar to the number of atoms required to express the structural features of the flow.
MD has played an exceptionally important role in the history of our understanding of fundamental issues in statistical mechanics, both equilibrium and nonequilibrium. First, in the early 1950s, MD helped establish the fact that to reach the thermodynamic limit was not at all the “astronomical number” (10^{24} particles) that had been assumed for more than a century; indeed, no more than 100 particles was necessary to evaluate equilibrium averages (17). Second, in the late 1960s, MD simulations demonstrated that fluctuations in an equilibrium fluid decay quantitatively by hydrodynamic modes, at a scale of 10 nm and a few picoseconds, resulting in a surprising hydrodynamical enhancement of the memory of a particle's velocity that was not foreseen by Boltzmann in his assumption of molecular chaos (18). Third, in the early 1970s, nonequilibrium MD simulations of simple steadystate flows of mass, momentum, and energy were evaluated as responses to external driving forces, similar to those carried out in real laboratory experiments (19), except that the forces were 6 orders of magnitude larger for the measured signal to overcome the noise caused by atomistic thermal fluctuations. Nevertheless, extrapolation of nonequilibrium MDshearing experiments to zero shear rate, even in the shearthinning regime, gave quantitative agreement with the equilibrium fluctuation–dissipation MD results for shear viscosity (20). Thus, the danger facing nonequilibrium MD simulations is that, for driving forces that are too large, the physics can change, such as in shear thinning at extremely high rates. One must assess the characteristics of the flow in response to the large gradients imposed to be sure that unwanted transitions have not occurred. However, nonequilibrium MD simulations of strong shock waves in fluids do not depart far from the NS description (21) even though such strong shocks exhibit significant departures from local thermodynamic equilibrium caused by steep gradients in the hydrodynamic fields. All these findings encourage an effort to extend MD from simple steadystate flows to violent chaotic flows. Previous nanohydrodynamics simulations for Rayleigh–Benard and Taylor–Couette flow (22, 23) have been restricted to linear stability analysis that is rigorously valid at short times, but they have not pursued the longtime behavior into the fully chaotic regime.
Methods
Atomic Interactions and Reduced Units. The atoms within each fluid interact by means of the intermediaterange cubic spline Lennard–Jones (LJ) potential (24): φ^{LJ}(r) = ε[(r _{0}/r)^{12} – 2(r _{0}/r)^{6}], for an r less than the inflection point; for a larger r, the potential is a cubic spline with a finite range. The interactions between the two fluids was either the same potential in the miscible case or purely repulsive in the immiscible case: φ^{rep} = φ^{LJ}(r) + ε, for r < r _{0}, and 0 otherwise. We use reduced LJ units: bond energy, ε; bond length, r _{0}; and the mass of the light fluid, m _{1} (except for the largest density ratio where we set m _{1} = 0.25). The mass of the heavy fluid was adjusted to give the desired A. If the light fluid is methanol (m _{1} = 32 atomic mass unit, ε/k _{B} = 507 K, r _{0} = 0.402 nm, k _{B} is Boltzmann's constant), the units for the time, t _{0} = r _{0}(m/ε)^{1/2}, and acceleration, , become 1.11 ps and 3.26 × 10^{14} m/s^{2} = 3.3 × 10^{13} g _{Earth}, respectively.
Simulation Setup. As the initial condition we set up a mechanically stable compressed fluid at temperature T = 1.3 ε/k _{B} under the influence of a gravitational field g. Starting with a given pressure P _{0} on top of the sample, we integrate the equation dP = –ρ(h)gdh numerically from the top to the bottom along the height h, filling up the space with a uniaxially compressed facecentered cubic (fcc) lattice. The atomic volume V is determined from an approximation of the isothermal equation of state, P(V,T = 1.3 ε/k _{B}). The quality of the approximation of the equation of state is not important for the results presented here. To speed up the melting of the fcc lattice, the positions were then randomized such that the closest possible atomic distance is no closer than 99% of the LJ potential minimum. Gravity ranged from to and gave, for a typical system height of 612.5 r _{0}, pressures on top from to and at the bottom from to . A further increase of g would reduce the wavelength of the most unstable mode and therefore the computational costs; however, this increase leads also to a solidification of the liquid near the bottom due to high pressures.
Linear Stability Analysis. A linear stability analysis of the hydrodynamic equations of a heavy fluid on top of a light fluid in the presence of a gravitational field was first presented by Harrison (25) in 1908, and later derived in a more general context by Chandrasekhar (ref. 8, see p.442, equation 113). To compare this linear stability analysis to atomistic simulations we used the surface tension and the viscosities (26) for the densities at the interfaces and solved the implicit expression derived by Harrison and Chandrasekhar numerically. [We calculated the surface tension (27) γ as a function of the pressurevolume PV at , 1.2 ε < PV < 1.7 ε.] This calculation gives the growth rate n as a function of k (Fig. 2), from which one can determine the mode of maximum growth.
Reynolds Number. The Reynolds number , where h is the height of the mixing zone and is the massaveraged shear viscosity of the two fluids (A. W. Cook, private communication) does not exceed 1,500 in the simulations, whereas experiments range from 3,700 to 10^{5} when spikes reach the bottom. Assuming h ∝ gt ^{2} provides and shows that to increase Re, either g or the system height H (which limits the possible fall height h) has to be increased. However, increasing g or h increases the pressure near the bottom such that solidification can occur. To increase the range of Re in atomistic simulations, g has to be decreased and the sample has to be widened (since the most unstable wavelength increases). Thus, the sample can be higher because the increase of pressure along the height is reduced.
Modeling. The scalable parallel shortrange MD code (refs. 28 and 29; http://bifrost.lanl.gov/MD/MD.html) with pairwise additive LJ interatomic interactions was used in the simulations. The particles in each of the two fluids interact with the same potential, except that one is more massive than the other. Particles in the two different fluids either interact with the same potential for the miscible case, where the surface tension is zero, or with only the purely repulsive part of the LJ potential in the immiscible case. The heavy fluid with density ρ_{2} and kinematic viscosity ν_{2} is initially placed on top of the light fluid of density ρ_{1} and kinematic viscosity ν_{1} under the influence of the gravitational field g. The density ratio is characterized by the Atwood number, A = (ρ_{2} – ρ_{1})/(ρ_{2} + ρ_{1}). The initial conditions were set up such that the system was as close as possible, given the equation of state, to a mechanically stable compressed fluid under the influence of the gravitational field. The interface roughness was then due only to atomicscale thermal fluctuations. Most runs used a thinslab geometry (1,400 r _{0} width, 612.5 r _{0} height, 5.25 r _{0} thickness, with the bond length r _{0}) containing 3 million particles, small enough and thus short enough runs (≈12 h on 200 processors of a Linux cluster) to investigate a wide range of parameter space. However, a few runs were also performed in a fully 3D geometry; an example containing 100,589,840 atoms (544.25 r _{0} × 525 r _{0} × 544.25 r _{0}) is given in Figs. 1 and 5, which is published as supporting information on the PNAS web site.
Computing Resources. Using 1,600 central processing units (CPUs) of the ASCI Q computer system (Compaq) at Los Alamos, we simulated 100,589,840 atoms for 250,000 integration time steps (δt = 0.01t _{0}) (400,000 CPU hours). Smaller simulations were performed on the ASCI QSC computer system (1,024 CPUs) and the 256CPU Linux cluster Grendels (Intel Xeon CPUs) at Los Alamos. The simulations presented here needed an estimated total of 600,000 CPU hours. Accounting for different processor speed on the used systems this is more than a week on all 4,096 CPUs of the ASCI Q computer system.
Results and Discussion
EarlyTime Mixing. To demonstrate the quantitative validity of the MD approach, the earlytime growth of an imposed single sinusoidal mode perturbation of the interface with wavelength λ was compared with analytical results of linear stability analysis based on continuum hydrodynamics (8, 25). The amplitude of the disturbance grows as e^{nt} (Fig. 6, which is published as supporting information on the PNAS web site), where the growth rate n depends on g, the surface tension, and the viscosities of the two fluids and is a function of the imposed wavenumber. This dependence on k reaches a maximum k _{max}, which is the fastest growing or most unstable mode, as shown in Fig. 2. It is important for quantitative purposes to solve the linear stability equation numerically rather than to depend on inaccurate analytic approximations (8). The MD simulations require an initial mode with an amplitude of 0.057λ to override the fluctuations always present in MD. The growth of the interface has been fit from t/δt = 120,000 to t/δt = 300,000 time steps (δt = 0.01 t _{0}), and the amplitude grows as much as 0.15λ at k _{max}. At earlier times the interface exhibits rather strong fluctuations, as it equilibrates from the artificial initial condition. The comparison between linear stability analysis of continuum theory and atomistic simulations in Fig. 2 shows good agreement with the value of the most unstable mode and deviations in the growth rate, which can be qualitatively explained by the nonlinearity introduced by the large initial amplitude. A numerical solution of the continuum equations for the immiscible but inviscid case (hence, no maximum) shows that an initial amplitude of 0.10(0.15)λ lowers the growth at k _{max} by 11(27)% (A. W. Cook, private communication). This figure also shows that fluctuations become more pronounced the smaller the structures are (i.e., the larger the k becomes). Here, thermal fluctuations have a larger influence on the growth than for the larger structures.
The growth from an interface perturbed only by fluctuations was found to be dominated quickly by the wavenumber corresponding to k _{max}, which then determines the initial number of bubbles and spikes. For example, the fully 3D system should, according to k _{max}, have 3.5 spikes along the edge of the simulation cell, whereas four were actually observed at early times (see Fig. 1). In some cases the deviations from the theoretically most unstable mode predictions were larger because the presence of multiple modes and the finite edge length.
LateTime Mixing. At later times the spikes and bubbles develop mushrooms at their tips, which eventually interact and merge in a process that can be described as turbulent (Movie 1, which is published as supporting information on the PNAS web site). The penetration distance h of the spikes and bubbles is characterized at long times by the only relevant distance scale applicable, namely, gt ^{2} multiplied by the coefficient αA, so that h = αAgt ^{2}, where h is determined at the position where the fluid concentration has changed by 95% (see Fig. 3). The fact that rise and fall of bubbles and spikes in the RT scenario has a constant acceleration rather than a constant velocity v (which is expected for a fall of an isolated droplet) can be explained by assuming the Newton friction force (–Cv ^{2}) rather than Stokes friction (–Cv), in addition to a constant buoyancy and inertia in the equation of motion for a bubble or spike. The friction coefficient C reduces with the fall, h, as C ∝ 1/h, since the structures become thinner during the fall (30). The longtime growth of the bubbles and spikes, as characterized by α, is thus determined from the slopes of the evolution of the boundaries shown in Fig. 3. Somewhat surprisingly, no statistically significant difference occurred between the average value of α, as determined from the large number of simulations done with the slab geometry at A = 0.867, and a single fully 3D system.
The growth of the bubbles and spikes, as given by the value of α determined by the simulations, were found to agree with experiment (5, 15) and Youngs' theoretical model (30) (see Fig. 4). This remarkable quantitative agreement with experiment at late times and in three dimensions in the turbulent regime is a significant validation of nanohydrodynamics. At the longest time of penetration of the spikes in the simulations, the depth reached 30% of the experimental value (5), measured in terms of the most unstable wavelength. However, the Reynolds number in the simulation reached only 1,500, compared with experimental values ranging from 3,700 to 10^{5}; NS calculations reached Reynolds numbers of 5,500 (6). Different continuum hydrodynamic calculations show a range of α values for bubbles between 0.03 and 0.08, compared with an experimental range of 0.05 to 0.08 (13, 14). Recently, it was suggested that the relatively large NS variation in α is due to the different initial perturbations along the interface (31, 32). Also, one should emphasize that most NS calculations use lower Atwood numbers than presented here. For spikes, α is not independent of A, as it is for bubbles, but tends toward the free fall limit at A = 1 of 0.5 (30). In this limit, the spikes are not hindered by any friction, whereas for bubbles, the dynamics is determined by a balance between friction and buoyancy.
Dependence on Initial Conditions. With MD simulation it is possible to quantitatively determine how α depends on initial conditions, such as the surface tension and imposed disturbances. It was found that the α values of the spikes in the miscible case were reduced compared with the immiscible case at a high Atwood number. In the miscible case, the interface gets blurred by diffusion. This result can be interpreted as a reduction of the effective Atwood number, leading to a reduction in the α of the spikes but not bubbles, since their α value is nearly independent of A. It was also found that an initial singlemode perturbation of the interface enhances the mixing process at short times, an effect also observed experimentally (5). However, at distances large when compared with the initial perturbing mode, the effect is washed out, and α is similar to the value from a fluctuationperturbed interface (Fig. 3). This result suggests a memory loss about the initial perturbations in time in the atomistic simulations due to the chaotic nature of the underlying equations of motions. However, continuum descriptions seem to depend on the initial perturbations on the α values (31, 32).
In any case, it is questionable whether α alone characterizes the longtime behavior (33). For example, from experiments (5) and the simulations presented in Fig. 3, oscillations in the bubble growth about the limiting slope are observed. These oscillations might have to do with the compressibility of the fluid, sometimes ignored in continuum calculations, coupled with viscous effects that have a major influence on bubble dynamics. Even more substantive are the observations of spikes in a slab that were five times higher than reported so far (3062.5 r _{0}). The purpose of that run (Movie 1) was to reach longer times, since in previous runs the spikes reached the bottom of the simulation cell. This simulation showed αAgt ^{2} behavior in the first half of the run but, after that, behaved dramatically differently: the mushrooms at the tips of the spikes broke off and formed droplets, an artifact that could either be due to an insufficient number of particles in the simulation or to the large value of g.∥ But it could also be real, when the thickness of the structures reach the nanoscale, since thinning of the spikes as they advance was observed experimentally (30). Furthermore, thinning and breakup of nanojets has been observed theoretically by including fluctuations, either by a full atomistic approach or by an additional stochastic term in the continuum equations (34).
Conclusions
The largescale atomistic simulations of the RT mixing process presented here have tested and quantitatively validated longstanding efforts in complex hydrodynamical flows by a completely different approach. Our simulations compare favorably with experimental observations and can help to interpret experimental results. Furthermore, these atomistic simulations pave the way to a detailed fundamental insight into turbulent flows in general, because they can capture the smallest physical scale, the atomistic scale. Initial perturbations arising from thermal fluctuations, which are unattainable by either NS calculations or experiments, can be studied by the nanohydrodynamics approach we have reported here. A detailed comparison with continuum calculations should be undertaken to see if hydrodynamic averaging over atomistic behavior introduces any significant differences into the description of the turbulent mixing process.
Acknowledgments
We thank Andy Cook (Lawrence Livermore National Laboratory) for discussions about hydrodynamic calculations and Leo Kadanoff for critical discussions, and we acknowledge QSC access through the Institutional Computing Program. This work was supported by Department of Energy Grant LDRD20020053DR and by the Advanced Simulation and Computing Program. Los Alamos is operated by the University of California under U.S. Department of Energy Contract W7405Eng36.
Footnotes

↵ † To whom correspondence should be addressed. Email: kkadau{at}lanl.gov.

Abbreviations: RT, Rayleigh–Taylor; MD, molecular dynamics; LJ, Lennard–Jones; NS, Navier–Stokes; CPU, central processing unit.

↵ ∥ This effect might be triggered by the fact that the density difference between top and bottom is huge for this tall sample, so that much less heavy fluid is in the simulation cell than light fluid. After some time, not enough heavy fluid exists to “feed” the spikes, which break up into droplets. To circumvent this problem (and to see the whether this effect can be observed in experiments), gravity must be reduced. Small gravity consumes much more computer time, since mixing takes longer and the most unstable wavelength increases.
 Copyright © 2004, The National Academy of Sciences
References

↵
Kadau, K., Germann, T. C. & Lomdahl, P. S. (2004) Int. J. Mod. Phys. C 15 , in press.

↵
Rayleigh, Lord (1883) Proc. London Math. Soc. 14 , 170–177.

↵
Taylor, G. I. (1950) Proc. R. Soc. London Ser. A 201 , 192–196.

↵
Di Prima, R. C. & Swinney, H. L. (1981) Hydrodynamic Instabilities and the Transition to Turbulence, eds. Swinney, H. L. & Gollup, J. P. (Springer, Berlin).
 ↵

↵
Cook, A. W. & Zhou, Y. (2002) Phys. Rev. E 66 , 026312112.

↵
Dzwinel, W., Alda, W., Pogoda, M. & Yuen, D. A. (2000) Physica D 137 , 157–171.

↵
Chandrasekhar, S. (1961) Hydrodynamic and Hydromagnetic Stability (Oxford Univ. Press, Oxford).

Harlow, F. H. & Welch, J. E. (1966) Phys. Fluids 9 , 842–851.

Youngs, D. L. (1990) Phys. Fluids A 3 , 1312–1320.
 ↵

↵
George, E., Glimm, J., Marchese, A., Li, X. L. & Xu, Z. L. (2002) Proc. Natl. Acad. Sci. USA 99 , 2587–2592. pmid:11854452
 ↵

↵
Kirkwood, J. G. (1947) J. Chem. Phys. 15 , 72–76.
 ↵

↵
Alder, B. J. & Wainwright, T. E. (1970) Phys. Rev. A 1 , 18–21.

↵
Hoover, W. G. & Ashurst, W. T. (1975) Adv. Theor. Chem. 1 , 1–51.
 ↵
 ↵
 ↵
 ↵
 ↵

↵
Harrison, W. J. (1908) Proc. London Math. Soc. 6 , 396–405.

↵
Ashurst, W. T. (1974) Ph.D. thesis (Sandia National Laboratories, Livermore, CA).

↵
Thompson, P. A., Brinckerhoff, W. B. & Robbins, M. O. (1993) J. Adhesion Sci. Technol. 7 , 535–554.

↵
Lomdahl, P. S., Tamayo, P. GronbechJensen, N. & Beazley, D. M. (1993) in Proceedings of Supercomputing 93, ed. Ansell, G. S. (IEEE Computer Society Press, Los Alamitos, CA), pp. 520–527.

↵
Beazley, D. M. & Lomdahl, P. S. (1997) Comput. Phys. 11 , 230–238.
 ↵

↵
Dimonte, G., Youngs, D. L., Dimits, A., Weber, S., Marinak, M., Wunsch, S., Garasi, C., Robinson, A., Andrews, M. J., Ramaprabhu, P., et al. (2004) Phys. Fluids 16 , in press.

↵
Dimonte, G. (2004) Phys. Rev. E 69 , in press.
 ↵

↵
Moseler, M. & Landman, U. (2000) Science 289 , 1165–1169. pmid:10947978
Citation Manager Formats
More Articles of This Classification
Physical Sciences
Related Content
 No related articles found.
Cited by...
 New directions for RayleighTaylor mixing
 Numerical simulations of twofluid turbulent mixing at large density ratios and applications to the RayleighTaylor instability
 Atomistic methods in fluid simulation
 RichtmyerMeshkov instability: theory of linear and nonlinear evolution
 The importance of fluctuations in fluid mixing