Scaling of atomistic fluid dynamics simulations
Kai Kadau,
1,
*
John L. Barber,
1,†
Timothy C. Germann,
1,‡
and Berni J. Alder
2,§
1
Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA
2
Department of Applied Physics, Lawrence Livermore National Laboratory, Livermore, California 94550, USA
Received 11 June 2008; published 27 October 2008
We have performed a series of large-scale atomistic simulations of the Rayleigh-Taylor instability including
up to 5.7 10
9
particles and spanning time and length scales of up to 170 ns and 45 m, respectively. The
results suggest that atomistic fluid dynamics simulations exhibit the same scaling as solutions of the continuum
Navier-Stokes equations. Furthermore, a comparison with macroscopic Rayleigh-Taylor experiments suggests
that the results of such atomistic simulations can, in fact, be scaled up to macroscopic dimensions, even for
complex, nonstationary flows.
DOI: 10.1103/PhysRevE.78.045301 PACS numbers: 47.27.-i, 47.11.Mn, 47.50.Gj
Fluid flows are usually described by continuum models
such as the Navier-Stokes equations. These equations can be
motivated as expressions of the fundamental conservation
laws, supplemented by linear constitutive relations between,
for example, stress and strain. However, continuum models
can also be shown to be strictly equivalent to the ensemble
mean of the underlying molecular dynamics 1. In other
words, solving for the dynamics of the molecules in a fluid is
similar to solving the continuum equations with random fluc-
tuations 2,3. In the case of particulate Brownian motion
several methods for including the random fluctuations in the
surrounding medium have been demonstrated 4 –6. It has
also been observed that fluctuations are important in the con-
text of phase separation in binary fluids 7. For stationary
flows, such as Couette flow or Rayleigh-Bénard convection,
it has previously been shown that the time average of the
particle dynamics is in agreement with the corresponding
continuum prediction, even for very small systems consisting
of only several thousand particles 8–10. This is not entirely
unexpected, as it is known, for example, that the long-time
tail of the velocity autocorrelation function in molecular sys-
tems can be explained by a continuum-based vortex mecha-
nism 11, indicating that the continuum approximation is
quantitatively applicable to such small systems. Further evi-
dence for the extension of the continuum description from
the macroscale 1m to the microscale 10
-6
m is pre-
sented in the work of Libchaber and Maurer 12, in which
Rayleigh-Bénard experiments at an intermediate scale
10
-3
m yield results in agreement with the predictions of
continuum hydrodynamics. However, for some nonstationary
flows it has recently been suggested that atomistic fluctua-
tions may lead to qualitatively different results from those
predicted by continuum models 13. In the present work, we
explore the question of whether atomistic simulations de-
scribing nonstationary flows can be scaled up to describe
macroscopic systems or whether the magnitude of the fluc-
tuations present at microscales alters the flow in such com-
plex problems.
As an example of a nonstationary, complex flow we will
consider the mixing of a heavy fluid density
1
, kinematic
viscosity
1
on top of a light fluid
2
,
2
in the presence of
a gravitational field g. In this Rayleigh-Taylor RT geometry
the density difference is usually characterized by the Atwood
number A =
1
-
2
/
1
+
2
. We assume that there is no in-
terfacial tension between the two fluids. Dimensional analy-
sis shows that the interfacial wavelength of maximum insta-
bility scales with g as
max
g
-1/3
and that the associated
exponential growth time scales as g
-2/3
. Linear stability
analysis yields the actual values as a function of viscosities,
gravity, and Atwood number 14,15. This is consistent with
the spatiotemporal scaling of the solutions of the incom-
pressible Navier-Stokes equations with gravity for the veloc-
ity field ux , t ; ghere, x, t, and P are the position, time, and
pressure, respectively, which remains a solution under the
transformation x → ax, t → a
2
t, g → g / a
3
, and u → u / a, P
→ P / a
2
+ P
0
, where the viscosity is assumed constant 13.
In order to investigate whether this scaling extends all the
way down to the atomistic level, we performed a series of
direct simulation Monte Carlo DSMC16 particle dynam-
ics simulations of RT mixing in a quasi-two-dimensional thin
slab geometry 13. Modern high-performance supercomput-
ing architectures, such as Livermore’s BlueGene/L system
17, enabled simulations containing up to 5.7 10
9
particles
Fig. 1. For most of the runs the Atwood number was 0.67
and the gravity ranged from 10
-3
g
0
to 7.813 10
-6
g
0
, which
corresponds to approximately 2 10
10
g
Earth
and 1.5
10
8
g
Earth
, respectively. By assuming that the light fluid
represents methanol, the units of length r
0
, time t
0
, and gravi-
tation g
0
become 0.402 nm, 1.44 ps, and 1.98 10
13
g
Earth
,
respectively 13,18. Higher gravities produce too much
vertical compression in the system, and lower gravities de-
mand even more computational capacity than is available
even on BlueGene/L, currently one of the world’s fastest
computer systems. The pressure at the interface was kept the
same for all simulations, in order to maintain constant trans-
port coefficients. The initial pressure and density profiles for
these compressible simulations were obtained from the static
equilibrium condition dP / dz =-g. Note, however, that it is
not in general possible to obtain a pressure profile that obeys
the above scaling exactly without altering the values of the
transport coefficients. As computational capacity increases
*
kkadau@lanl.gov
†
jlbarber@lanl.gov
‡
tcg@lanl.gov
§
alder1@llnl.gov
PHYSICAL REVIEW E 78, 045301R2008
RAPID COMMUNICATIONS
1539-3755/2008/784/0453014 ©2008 The American Physical Society 045301-1