7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Anisotropic LennardJones fluids in a nanochannel
R.M. Hartkamp and S. Luding
Multiscale Mechanics, Faculty of Engineering Technology, University of Twente, Enschede, 7500AE, The Netherlands
r.m.hartkamp~ utwente.n1 and s.1uding @utwente.n1
Keywords: Molecular dynamics, LennardJones, anisotropy, layering, virial stress
Abstract
During the past few decades molecular dynamics has been a widely applied tool to simulate fluid confined in mi
cro/nano geometries. What makes interfacial fluids fundamentally different from the bulk fluid is the fact that their
density varies considerably over microscopic distances. A class of such strongly inhomogeneous fluids are fluids con
fined in very narrow channels by solid boundaries. In this work, the goal is to study the density and stress terms across
the channel.
We simulate planar Poiseuille flow of a LennardJones fluid in channels of various widths in the nanoscale regime.
A body force and a local thermostat are applied in order to simulate a steadystate flow. Layering and anisotropy
in stress are obtained near the walls of the channel, which leads to nonNewtonian theology. Understanding and
quantifying the nonNewtonian behavior is a first step towards deriving a constitutive model that describes locally the
behavior of a strongly confined fluid.
Nomenclatu re
7 Natural time step
# Coarse graining function
Miscellaneous
(...PParticle property
(. ..)p4 Interaction property
(..) Kinetic stress property
...)"Potential stress property
Introduction
Molecular dynamics simulations have become an im
portant tool for the study of microscopic fluid proper
ties. A channel geometry is often used to study the in
homogeneous behavior of strongly confined fluids [Ko
plik & Banavar (1988); Bitsanis et al. (1988); Todd
& Evans (1995); Travis et al. (1997); Travis & Gub
bins (2000)]. However, our understanding of these non
Newtonian fluid problems is still very limited, while
gaining a deeper insight is becoming increasingly im
portant with the rise of microfluidic and nanofluidic ap
plications, such as labonachip devices [van den Berg
et al. (2009)]. Similar phenomenology (i.e. layering,
anisotropy) is observed in many particle systems [van
Gunsteren et al. (1984); Ghosh et al. (2007)]. In this
study, liquid argon is confined between two flat walls
with a normal in the zdirection (Fig. 1). When an atom
leaves the system in y or zdirection, it reenters at
Roman symbols
b Bin width
f Reduced body force
F Interaction force
m Particle mass
N Number of particles
p Hydrostatic pressure
Interaction distance
re Cutoffradius
t Time
T Period of oscillation
T* Reduced temperature
a Streaming velocity
U Interaction potential
Particle velocity
Fluctuation velocity
V Volume
H Channel width
Greek symbols
SEnergy scale
egj Deformation tensor
v Volume fraction
nII Stress tensor
p* Reduced density
o Length scale
oO aoo o o
ob o 06
.~o 00000
oo 0 0 0 o go
yO x
Figure 1: Liquid argon (grey, 1536 atoms) confined be
tween solid argon walls (red, 512 atoms). The
width of the channel is defined as shown on
the right.
the opposite side due to periodic boundary conditions.
Both walls consist of two layers of argon atoms (each
layer containing 128 atoms) fixed in a face centered cu
bic (fec) lattice (100 direction). The fluidwall inter
action is equal to the fluidfluid interaction and the re
duced density of the system is p* 0.8 (corresponding
to a volume fraction of v 0.415). Physical quanti
ties are nondimensionalized by using the length, energy
and mass scales of liquid argon, which are a 3.405 x
10in m, t 1.67 x 102 J and a t 6.626 x 10"' kg
respectively. A constant body force f* acts on the fluid
in negative zdirection, causing it to flow. The mutual
interaction of the neutral spherically symmetric argon
atoms is modeled via a twobody LennardJones (LJ)
potential:
U(r) =4r (() (%
with r the distance between two atom centers. From the
interaction potential, the force between two atoms can
be calculated:
aU
FL7(r) =
~r
( ( is
.(2)
The force is truncated at re = 2.5a in order to reduce
calculation time, therefore F(r > re) = 0. Further
more, the force is shifted with F(re) in order to maintain
a continuous force at the location of truncation:
( aiL is \
a
FL./Ec). (3)
The natural time step, 7 = (ma2/t)1/2, is propor
tional to the period of oscillation around the potential
minimum. For liquid argon, it is 7 2.14 x 1012
s. From the positions, velocities and interatomic forces,
other physical scalar or tensorial quantities can be cal
culated (e.g. shear rate, stress and viscosity) [Hartkamp
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
005 10 15
Figure 2: Density profile for three different channels
with width 101= 5.985a (blue), if,
11.115a (red) and W3 16.245a (black).
The vertical lines represent the location of the
right wall for that system.
et al. (2010)]. The position and momentum of every
atom and the forces acting on them are used to calcu
late their position and velocity after an increment (At)
in time via the Velocity Verlet algorithm, with a time step
of I 0.0017. The body force leads to an acceleration
of the fluid, while viscous effects retard the flow until a
steady state is reached. Local thermostats near the walls
keep the energy (and thus the temperature) constant in
time (T* 1.0 for the simulations presented here)
[Ghosh et al. (2007); Bartos & Jainosi (2007)]. Atoms
are initially positioned on an fcc lattice. The lattice melts
during the equilibration, followed by a steady state sim
ulation. The steadystate simulation results are averaged
over a period of time of the order of 10007 50007. Fur
thermore, spatial averaging is employed over the direc
tions parallel to the solid walls, whereas, the zdirection
(perpendicular to the solid walls) is divided into equally
spaced bins of width b = 0.083.
In the following sections we first discuss the density
profile and its relation to the channel width. Next, the
computation of the local stresses in the fluid is discussed.
Some results are presented for the stresses across the
channel. Finally, the presented theory and results are
discussed and certain topics and ideas for future work
are briefly motivated.
Density
Fig. 2 shows a volume fraction profile in the channel for
different channel widths. Each of the profiles shows os
cillatory behavior near the channel walls which decays
exponentially away from the walls. The peak closest
to the wall is due to the wall atoms, the density of the
fluid goes to zero there. The period and magnitude of
the most apparent oscillations in the fluid density seems
to be approximately identical for all three simulations
(apart from statistical noise). The oscillatory density
profile near the wall can be captured by a function of
the following form:
(2i7r~  (x x,)
(4)
The f and v signs correspond to a fit near the left and
the right wall respectively. Physical properties, such as
the bulk volume fraction vo and the period (~ wave
length) L of the density oscillations, can be identified.
The exponential decay away from the wall is quantified
by xo and the amplitude of the density oscillations (at
the wall) are fitted with a. It must be noted that, unlike
the period of oscillation and the bulk volume fraction,
the amplitude and decay of the oscillatory peaks, are
strongly dependent on the coarse graining of the data
(in this work, a Gaussian function is used to coarse
grain the information). A leastsquares fit of the data
is made, where the oscillation closest to the wall is not
taken into account since this data does not solely repre
sent the density of the fluid, but also the wall particles
are contributing. Table 1 shows the parameter values for
a leastsquares fit of the left and right half of the channel
of width W3 = 16.245a and Fig. 3 shows the simula
tion data and the exponential fit. The value for L shows
that the approximate period of the oscillations is smaller
than the length scale of the argon atoms for this simu
lation, but larger than the initial lattice spacing, which
is 0.855a. The period of the oscillations, however, does
depend on the average density and temperature in the
system (data not shown). The bulk volume fraction vo
is (almost) equal to the average volume fraction in the
system, these values could deviate for example in a sim
ulation where the parameters for interaction between the
fluid and the wall are different from the interaction be
tween fluid particles.
Table 1: Fitting parameters.
Parameter Left Right
vo 0.4146 0.4146
a0.2324 0.2322
x, 0.0 16.245
L 0.9326 0.9322
xo 1.1509 1.2010
In the largest system shown in Fig. 2, a bulk fluidlike
(i.e. homogeneous) region can be identified in the cen
ter. Six distinct peaks are observed between each wall
and the bulk fluid region of the channel. Decreasing
the system size, but keeping the density and tempera
ture unchanged, the homogeneous region shrinks in size
and finally disappears when the channel width becomes
smaller than approximately 12a. For narrower channels,
a region forms in the center where the oscillations in
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
0.4
E
S0.2
0
0
0.8
c 0.6
o
S0.2
simulation
 fit
10 12
14 16
Figure 3: (a) Volume fraction in the left half of the chan
nel. (b) Volume fraction in the right half of the
channel.
density, induced by both walls, interfere additively; the
oscillatory behavior in this region becomes directly de
pendent on the system size. It will be shown that by
simply taking the density profile in the left and right in
homogeneous region of the largest system, one can pre
dict the density profile in a system that is either larger
or smaller. Consider the density profile in the left and
right half of a channel that is wider than twice the inho
mogeneous region (a channel of width W3 16.245a
is used here) and shift these towards each other until a
channel of width W, 11.115a is realized. The de
viations from the average density are then just summed
up. In Fig. 4, the constructed profile is compared to the
result of a molecular dynamics simulation.
The figure confirms that the interference of the oscil
lations is only \ignilk~.lill in the region that is within 6a
distance from both walls and that the oscillations can be
summed up approximately. The region where the influ
ence of both walls overlaps is small in this system (W,)
and so is the amplitude of the oscillations. If we apply
the same approach to a more narrow channel (W1)(see
Fig. 5), where the overlapping region dominates the sys
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
width = 11.1150
*data
left wall effect
right wall effect
constructed profile
width = 5.9850
*data
left wall effect
right wall effect
1
0.8
0.6
ti0.6
0
5 0 5
x []
2 1 0
1 2
Figure 5: The comparison of a constructed density pro
file and the data of a molecular dynamics sim
ulation for a channel of width W1 = 5.985a.
Stress
Molecular dynamics simulations evolved over the last
decades to a much used tool for studying many types of
smallscale problems. Yet, how to calculate the macro
scopic stress from microscopic information of individual
particles and interactions is still a much debated subject.
An important similarity between continuum mechanics
and molecular dynamics is that they both have to obey
the conservation of mass and momentum.
The virial stress is often used to calculate the stresses
in atomic systems. The virial stress expression was first
derived by Lutsko (1988), based on a generalization of
the virial theorem presented by Clausius (1870). No
formulation of the virial stress is fully accepted to date.
Zhou (2003) wrongly argues that the virial stress is to
tally irrelevant to the mechanical stress and has no phys
ical significance at all. One of the arguments that Zhou
presents is that momentum can only be transported by
mass, and force can only work on a mass. This belief
is contradictory to the generalization made by Lutsko,
according to Zhou. In response, a number of articles
were published, showing the physical significance of the
virial stress (van Dommelen (2003); Subramaniyan &
Sun (2008)). The formulation used here is in accordance
with Todd et al. (1995), apart from details discussed be
low.
The virial stress is given as the sum of a kinetic part,
which accounts for the momentum transport of the par
ticles, and a potential part, which represents the contri
bution of the interaction between all pairs of particles:
IIi "(x) = II," (x) + II (x).(5
The sum of the kinetic and potential contribution is
sometimes multiplied with a Dirac delta function (Eq.
1.7 in Zhou (2003)), which assigns a location to the
E o.0.4 2 1 ( l
0.38
2 x0] 2
Figure 4: (a) Comparison of a constructed density pro
file and the data of a molecular dynamics sim
ulation for a channel of width We = 11.115a.
(b) Closeup of the center of the channel in
(a).
tem (i.e. the density profile in the whole system is di
rectly influenced by both walls), we see a great corre
spondence to the density profile that was obtained by a
molecular dynamics simulation, in both the magnitude
and period of the oscillations.
The approach that was presented here can be used
to predict the oscillatory behavior in a channel with
out having to simulate the system explicitly. However,
more extensive testing of the method is required, espe
cially when the system width decreases below 5.985a
and for systems with other densities and temperatures.
If a quantitative relation between the density profile and
other physical properties can be found for this inhomo
geneous channel flow, the method presented here could
predict layering and the flow behavior of the fluid with
out the need of simulating every individual system width
explicitly.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
atomic stress contribution. While this approach is sim
ple and mathematically convenient, we believe that the
kinetic and the potential term should be coarsegrained
individually as shown below.
The kinetic stress contribution in a bin yields:
inptbin
where mp is the mass of the particles, #9(x) represents
a coarse graining function and if (vP ag(xP)) is
the fluctuation velocity of a particle p (i.e. the differ
ence between the particle velocity and the local stream
ing velocity), where ag (xP) is the local streaming ve
locity subscriptt i indicates the components in the i
direction), defined here as the average velocity of parti
cles in the bin averaged over time. This approach results
in a streaming velocity which slightly depends on the
spatial discretization. This dependence disappears in the
limit of small bins and large averaging time, however,
a proper smoothing of the data is paramount. The fluc
tuation velocity is generally used for the calculation of
kinetic properties such as temperature or the kinetic part
of stress. Some authors have used the total velocity of
the particle rather than the fluctuation velocity, however,
stress obtained with total velocities is not objective.
The potential energy contribution of stress is obtained
by computing the change in the potential energy density
(due to a virtual deformation) with respect to the com
ponents of the deformation tensor (es,):
BU
In,(x) = (7)
Following the above expression, through some calcula
tions it can be shown that the potential stress tensor can
be obtained from the pairwise interactions as:
Inf (x)=M 1 rf 'F,PY999(x)! (8)
ptbin q p
where compressive stresses are defined to be positive
here. The socalled branch vector component r(' =
r(  rf denotes the ith component of the distance vector
between particles p and q and FF'4 the je" component of
the interaction force between particles p and q as defined
in Eq. (3).
The potential stress contribution contains a coarse
graining function, FR4(x), in which the distribution of
information over the branch vector between the particles
is modeled. The formulation presented here does not
impose restrictions on the coarsegraining functions (for
either part of the stress tensor). A more extensive dis
cussion about these functions is given in Hartkamp et al.
(2010).
Figure 6: The kinetic and potential contributions to the
stress term II, across the channel.
The factor 1/2 in Eq. (8) accounts for the fact that
in this expression each pair interaction is taken into
account twice. Since Newton's third law states that:
F~p FP Eq. (8) can be reduced to:
n"i~i Vbin r4Fqp()
ptbin q>p
The presented expression for II@(x) and Inf,(x) are
instantaneous atomic stress contributions. These instan
taneous quantities have rather strong fluctuations, how
ever, averaging both terms over time results in the con
tinuum Cauchy stress (Subramaniyan & Sun (2008)).
Results
The stresses in the system of 16.2450 width, a reduced
density of p* 0.8 and a reduced body force of f*
0.1 are presented here. Figs. 6, 7, and 8 show the nor
mal stresses in the x, y and zdirection respectively. In
the bulk, the normal stresses are approximately equal in
each direction. Near the wall, oscillations occur, similar
to the oscillations in density shown above. The oscilla
tions are not identical for the different normal directions,
indicating anisotropy in this region. The normal stress
II, in the direction of confinement (i.e. perpendicular
to the walls) shows larger oscillations in comparison to
the directions parallel to the walls. The negative values
for stress correspond to tensile stresses, these can occur
locally in an anisotropic fluid, however, more study is
required in order to fully understand this local stress be
havior and the influence of the coarse graining procedure
on it. The hydrostatic pressure (Fig. 9) follows from the
average of the normal stresses p (I1xx+11,y+11zz)/3.
Fig. 10 shows the shear stress II, across the chan
nel. The kinetic contribution is very small compared to
the potential part of the shear stress. The shear stress
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
......kinetic
. .. potential
 total
 kinetic
 potential
 totalI
5 0
x [5]
5 10
5 0
x [(5]
5 10
Figure 7: The kinetic and potential contributions to the
stress term II,, across the channel.
Figure 9: The kinetic and potential contributions to the
hydrostatic pressure across the channel.
............. k inetic
 potential
 totalI
1
0.5
S0
0.5
1

1.5
10
5 0
5 10
Figure 8: The kinetic and potential contributions to the
stress term II,, across the channel.
shows oscillations superimposed on a linear trend across
the channel. The linear trend corresponds to the shear
stress in a Newtonian fluid, whereas, the oscillations
again represent the deviations from a Newtonian fluid.
The kinetic and potential contribution of the other non
diagonal components of the stress tensor, II,, and II,,
are fluctuating around zero, as may be expected.
Discussion and future work
Molecular dynamics simulations are performed for a no
ble fluid in a nanochannel. It is shown here, that the os
cillations in density, close to the wall, can be character
ized by a certain period L, amplitude <1 and decay length
rn. A likely dependence on average density, temperature
and interaction parameters is not studied here.
The calculation of stress in a nanochannel is dis
cussed. Analyzing the local stresses in a strongly con
fined fluid system is important in order to study prop
erties such as anisotropy. Due to the inhomogeneity in
these systems, the flow behavior deviates from that of a
Newtonian fluid. A better understanding of oscillatory
Figure 10l: The kinetic and potential contributions to the
stress term II,, across the channel.
behavior in the density and stress terms is paramount to
deriving a new constitutive model that captures the non
Newtonian behavior of the fluid. In addition to the total
virial stress, results are shown for the kinetic and poten
tial contribution to the virial stress. The results show
that, while the potential part dominates over the kinetic
part, neither contribution is negligibly small for the nor
mal stresses in the present situation. Only the kinetic
contribution to II., is very small compared to the poten
tial contribution.
The local value of various physical quantities depend
on a the spatial discretization and a proper coarse grain
ing of the information, as mentioned above. The width
of the bins determines the resolution of the averaged
data. Decreasing the bin size, results in fewer atoms
per bin, and thus worse statistics. The quality of the
statistics can be improved by coarse graining of the data
with some type of nonlocal smoothing function. The
physical justification for this action is the fact that par
ticles have a finite size. Since atoms don't have a spec
ified exact radius, a Gaussian distribution is suitable to
distribute the information of an atom over space (e.g.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
van Gunsteren W.F., Berendsen H.J.C., Postma J.P.M.,
Di~ola A. and Haak J.R., Molecular Dynamics with
Coupling to an external bath, Joumnal of Chemical
Physics, Vol. 81, pp. 3684, 1984
Hartkamp R.M., Ghosh A. and Luding S., Anisotropy of
a fluid confined in a nanochannel, In preparation, 2010
Koplik J. and Banavar J.R., Molecular dynamics simu
lation of microscale poiseuille flow and moving contact
lines, Physical Review Letters, Vol. 60, pp. 1282, 1988
Lutsko J.F., Stress and elastic constants in anisotropic
solids: Molecular dynamics techniques, Joumnal of Ap
plied Physics, Vol. 64, pp. 1152, 1988
Subramaniyan A.K. and Sun C.T., Continuum interpre
tation of virial stress in molecular simulations, Inter
national Joumnal of Solids and Structures, Vol. 45, pp.
43404346, 2008
Todd B.D. and Evans D.J., The heat flux vector for
highly inhomogeneous nonequilibrium fluids in very
narrow pores, Joumnal of Chemical Physics, Vol. 103, pp.
9804, 1995
Todd B.D., Evans D.J. and Daivis, P.J., Pressure ten
sor for inhomogeneous fluids, Phys. Rev. E, Vol. 52, pp.
1627, 1995
Travis K.P., Todd B.D. and Evans D.J., Departure from
NavierStokes hydrodynamics in confined liquids, Phys
ical Review E, Vol. 55, pp. 4288, 1997
Travis K.P. and Gubbins K.E., Poiseuille flow of
LennardJones fluids in narrow slit pores, Joumnal of
Chemical Physics, Vol. 112, pp. 1984, 2000
Zhou M., A new look at the atomic level virial stress: on
continuummolecular system equivalence, Proceedings
of the Royal Society of London. Series A: Mathematical,
Physical and Engineering Sciences, Vol. 459, pp. 2347
2392, 2003
for ## in Eq. (6)). Similarly, the force between an in
teracting particle pair does not act solely on the cen
ter of the atom. Thus, also a coarse graining function
is needed to describe the microscopic force field (#pG
in Eq. (9)), since atoms do not have a force concentra
tion on a welldefined contact point, such as interacting
macroscopic particles have. Little is known about how
to distribute the information of an atomic (longrange)
interaction over space. Goldhirsch (2010) developed a
promising method to coarsegrain microscopic informa
tion of interacting particles. A more elaborate analysis
of this method will be given in later work.
The results from the virial stress can be verified by
calculating the force (rate of change in momentum) on
the walls of a simulation system. Dividing this force by
the area of the walls gives the corresponding stress vec
tor and stress components. The obtained values can be
compared to those obtained from the virial stress. This
method is pretty straightforward for a system in which
the particles have no average streaming velocity (i.e. no
body force acts on the particles).
Acknowledgements
This work was financially supported by MicroNed grant
4A2.
References
Bartos I. and J~nosi I.M., Side pressure anomalies in 2D
packing of frictionless spheres, Granular Matter, Vol. 9,
pp. 81, 2007
van den Berg A., Sparreboom W. and Eijkel J.C.T., Prin
ciples and applications of nanofluidic transport, Nature
Nanotechnology, Vol. 4, pp. 713, 2009
Bitsanis I., Vanderlick T.K., Tirrell M. and Davis H.T.,
A tractable molecular theory of flow in strongly inho
mogeneous fluids, Joumnal of Chemical Physics, Vol. 89,
pp. 3152, 1988
Clausius R., On a mechanical theory applicable to heat.,
Philosophical Magazine, Vol. 40, pp. 122127, 1870
van Dommelen L., Physical interpretation of the
virial stress, http://www.eng. fsu. edu/ domme
len/papers/viriallindex.pdf, 2003
Ghosh A., Paredes R. and Luding S., Poiseuille flow in a
nanochannel use of different thermostats, Intemnational
Congress of Particle Technology, CD Proceedings, pp.
14, 2007
Goldhirsch I., Stress, stress asymmetry and couple
stress: from discrete particles to continuous fields, Gran
ular Matter, Vol. 12, in press, 2010
