See discussions, stats, and author profiles for this publication at: [Link]

net/publication/216546915

Flow Simulations Using Particles

Conference Paper · May 2009

CITATIONS READS

2 324

3 authors:

Petros Koumoutsakos G.-H Cottet

ETH Zurich Grenoble Alpes University
431 PUBLICATIONS 23,214 CITATIONS 147 PUBLICATIONS 3,887 CITATIONS

SEE PROFILE SEE PROFILE

Philippe Chatelain
Catholic University of Louvain
185 PUBLICATIONS 2,434 CITATIONS

SEE PROFILE

All content following this page was uploaded by Petros Koumoutsakos on 31 May 2014.

The user has requested enhancement of the downloaded file.

 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD
10.1146/[Link].37.061903.175753

Annu. Rev. Fluid Mech. 2005. 37:457–87

doi: 10.1146/[Link].37.061903.175753
Copyright

c 2005 by Annual Reviews. All rights reserved

MULTISCALE FLOW SIMULATIONS USING

PARTICLES
Petros Koumoutsakos
Computational Science and Engineering Laboratory, Swiss Federal Institute of
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

Technology, Zurich, CH-8001, Switzerland; email: petros@[Link]

Key Words vortex methods, smooth particle hydrodynamics, molecular dynamics,

continuum-molecular simulations, domain decomposition
■ Abstract Flow simulations are one of the archetypal multiscale problems. Sim-
ulations of turbulent and unsteady separated flows have to resolve a multitude of
interacting scales, whereas molecular phenomena determine the structure of shocks
and the validity of the no-slip boundary condition. Particle simulations of continuum
and molecular phenomena can be formulated by following the motion of interacting
particles that carry the physical properties of the flow. In this article we review La-
grangian, multiresolution, particle methods such as vortex methods and smooth particle
hydrodynamics for the simulation of continuous flows and molecular dynamics for the
simulation of flows at the atomistic scale. We review hybrid molecular-continuum sim-
ulations with an emphasis on the computational aspects of the problem. We identify the
common computational characteristics of particle methods and discuss their properties
that enable the formulation of a systematic framework for multiscale flow simulations.

1. INTRODUCTION
The simulation of the motion of interacting particles is a deceivingly simple, yet
powerful and natural, method for exploring physical systems as diverse as planetary
dark matter and proteins, unsteady separated flows, and plasmas. Particles can be
viewed as objects carrying a physical property of a system, that is being simulated
through the solution of Ordinary Differential Equations (ODEs) that determine
the trajectories and the evolution of the properties carried by the particles. Particle
methods amount to the solution of a system of ODEs:

dx p
N
= u p (x p , t) = K (x p , xq ; ω p , ω q ) (1)
dt q=1

dω p
N
= F(x p , xq ; ω p , ω q ), (2)
dt q=1

0066-4189/05/0115-0457$14.00 457
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

458 KOUMOUTSAKOS

where x p , u p denote the locations and velocities of the N particles, ω p denote

particle properties (such as density, temperature, velocity, vorticity), and K , F
represent the dynamics of the simulated physical system. In flow simulations par-
ticles are implemented with a Lagrangian formulation of the continuum equations,
as in the vorticity formulation of the Navier-Stokes equations, or with systems that
are discrete by nature, as in molecular flows at the nanoscale. Continuum flows,
such as flows in porous media and unsteady separated and turbulent flows, are in-
herently multiscale due to the range of scales that govern the underlying physical
phenomena. The continuum assumption fails in flow regions containing contact
lines and shocks, and suitable molecular descriptions become necessary. A consis-
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

tent and systematic framework is necessary to couple molecular and macroscale

descriptions because the macroscale flows determine the external conditions that
influence the molecular system, which in turn influences the larger scales by mod-
ifying its boundary conditions.
Particle methods such as Vortex Methods (VMs) and Smooth Particle Hydro-
dynamics (SPH) present an adaptive, efficient, stable, and accurate computational
method for simulating continuum flow phenomena and for capturing interfaces
such as vortex sheets. On the other hand, particle methods encounter difficulties in
the accurate treatment of boundary conditions, while their adaptivity is often asso-
ciated with severe particle distortion that may introduce spurious scales. Ongoing
research efforts attempt to address these issues as outlined in the review.
In molecular and mesoscopic simulations particle methods, such as Molecular
Dynamics (MD) and Dissipative Particle Dynamics (DPD), are the methods of
choice because the discrete representation of the underlying physics is inherently
linked to interacting particles. Particle methods for continuum and discrete systems
present a unifying formulation that can enable systematic and robust multiscale
simulations, as we outline in this review.
A remarkable feature of particle methods is that their computational structure
involves a large number of common abstractions that help in their computational
implementation, while at the same time particle methods are distinguished by the
fact that they are inherently linked to the physics of the systems that they simulate.
In this review we focus on updating the reader in methodological advances in
Lagrangian particle methods since the first related such review by Leonard in 1985
(Leonard 1985), with an emphasis toward describing methodologies that enable
multiresolution simulations. In the simulation of discrete systems, starting from the
review of Koplik & Banavar (1995), we focus on hybrid continuum-molecular flow
simulations. In this article we do not discuss particle methods for the simulation
of kinetic equations, for which we refer the reader to Chen & Doolen’s (1998)
work.
The review is structured as follows: We introduce particle methods for con-
tinuous systems by illustrating unifying concepts such as function and derivative
particle approximations. We discuss the fundamental problem of particle distortion
and remeshing associated with the Lagrangian formulation and introduce multires-
olution particle methods. We briefly outline the key characteristics of molecular
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

FLOW SIMULATIONS USING PARTICLES 459

simulations and discuss recent advances in hybrid continuum-molecular simula-

tions. We conclude by describing efficient tools for large-scale simulations using
particle methods and provide an outlook for future developments in particle meth-
ods in this new era of multiscale modeling and simulation.

2. PARTICLE METHODS FOR CONTINUOUS SYSTEMS

Particle methods for continuum flow simulations include VMs and SPH. The
key common characteristic of these methods involves the approximation of the
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

Lagrangian form of the Navier-Stokes equations by replacing the derivative op-

erators through equivalent integral operators that are in turn discretized on the
particle locations.

2.1. Particle Function Approximations

Point particle approximations were the first to attract attention in solving fluid
mechanics problems because their evolution can be formulated in terms of conser-
vation laws. An approximation of a smooth function f in the sense of measures
(Raviart 1986) can be formulated as:


N
fh = w p δ(x − x p ), (3)
p=1

where w p denotes the weights of the particles. Although the point particle approx-
imation has several interesting properties, particularly when considering exact
formulations of conservation equations, smooth function approximations are of-
ten desirable, allowing recovery of the function between particle locations and a
regularized formulation of the particle motion. Smooth function approximations
can be constructed by using a mollification kernel ζε (x):

f ε (x) = f  ζε = f (y) ζε (x − y) dy, (4)

where ε denotes a characteristic length of the kernel.

The particle approximation of the regularized function is defined as


N
f εh (x) = f h  ζε = w p ζε (x − x p ). (5)
p=1

The error introduced by the quadrature of the mollified approximation f εh for the
function f can be distinguished in two parts as
f − f εh = ( f − f  ζε ) + ( f − f h )  ζε . (6)
The first term in Equation 6 denotes the mollification error that can be controlled
by appropriately selecting the kernel properties. The second term denotes the
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

460 KOUMOUTSAKOS

quadrature error due to the approximation of the integral on the particle locations.
Since the early 1980s, mollifier kernels have been developed in VMs with an em-
phasis on the property of moment conservation to comply with vorticity moments
conserved by the Euler equations. The accuracy of these methods is related to the
moments that are being conserved, and a method is of order r when:

 ζ (x) dx = 1

xi ζ (x) dx = 0 if |i| ≤ r − 1 (7)

 r
|x| |ζ (x)| dx < ∞
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

The overall accuracy of the method is then:


hm
f − f εh 0, p ∼ O(ε ) + O
r
. (8)
εm
For equidistant particle locations at spaces h in a d-dimensional space, the weights
can be chosen as: w p = h d f (x p ) with m = ∞ for certain kernels and for positive
kernels such as the Gaussian, r = 2. Kernel cutoffs of arbitrary order (Beale 1986)
are possible by giving up the positivity of the cutoff. These error estimates reveal
an important, albeit often overlooked, fact for smooth particle approximations:
to obtain accurate approximations smooth particles must overlap. Note that the
moment conditions expressed by the integrals of the mollifier functions are not
often well represented for discrete particle sets. These moment conditions can be
ensured by appropriate normalizations (Cottet & Koumoutsakos 2000).

2.1.1. PARTICLE DERIVATIVE APPROXIMATIONS Although the representations of

functions by particles can be considered a post-processing step, the approximation
of derivatives is a key aspect in the development of particle methods for solving
the governing flow equations.
Particle approximations of the derivative operators can be constructed through
their integral approximations. This can be easily achieved by taking the derivatives
of Equation 4 as convolution and derivative operators commute in unbounded or
periodic domains. These approximations can be cast in a conservative formulation
and are extensively employed in SPH.
An alternative formulation involves the development of integral operators that
are equivalent to differential operators such as the Laplacian. Motivated by the
need to construct high-order viscous algorithms for VMs, in 1987 Mas-Gallic
introduced the method of Particle-Strength Exchange (PSE). The PSE scheme can
be derived starting from a straightforward Taylor expansion of f around x:

∂2 f
f (y) = f (x) + (y − x) · ∇ f (x) + (xi − yi )(x j − y j ) + · · · . (9)
i, j
∂ xi ∂ x j

Convolving this expansion with an even function η, the first-order terms and the
cross-terms involving the second-order derivatives of f drop out in the righthand
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

FLOW SIMULATIONS USING PARTICLES 461

side of Equation 9 and, using the normalization conditions,


x
ηε (x) = ε −d η , xi2 η(x) dx = 2 i = 1, · · · , d (10)
ε
leads to the approximation (Degond & Mas-Gallic 1989a):

ε f (x) = ε −2 ( f (y) − f (x)) ηε (y − x) dy, (11)

where ε f (x) denotes the mollified approximation of the Laplacian operator.

by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.

High-order approximations can be obtained by choosing suitable functions η. The

Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

anisotropic extension of this method is defined as

d 


−2
∇ · [B∇ f ](x)  ε ψiεj (x − y)Mi j (x, y) [ f (y) − f (x)] dy, (12)
i, j=1

where Mi j (x, y) are symmetric functions and ψiεj are cut-off functions, related to
each other and to the matrix B through conditions that are detailed in Degond
& Mas-Gallic (1989b). Starting from the PSE formulation Eldredge et al. (2002)
presented a general deterministic integral representation for derivatives of arbitrary
order. The error analysis of particle derivative approximations strengthens the
requirement for particle overlap.
In particle methods the precise connectivity of the computational elements (as,
for example, in finite difference methods) is not required to discretize the gov-
erning equations, but neighboring elements need to overlap to provide consistent
approximations.

2.2. Vortex Methods

Vortex particle methods have been used since the 1930s (Rosenhead 1930) to
describe the evolution of vortical structures in incompressible flows.
Navier-Stokes equations describe the evolution of the vorticity field in 3D,
incompressible, viscous flows in a velocity-vorticity (u, ω = ∇ × u) formulation
as
∂ω
+ (u · ∇) ω = (ω · ∇) u + νω (13)
∂t
The velocity field u is obtained by solving the Poisson equation

∇ 2 u = −∇ × ω (14)

with suitable boundary conditions (Cottet & Koumoutsakos 2000). Velocity calcu-
lations, satisfying explicitly far-field boundary conditions, are based on the Biot-
Savart law:

u = K(x − y) × ω dy + U 0 (x, t), (15)
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

462 KOUMOUTSAKOS

where U 0 (x, t) is the solution of the homogeneous Equation 14, and K(z) denotes
the Biot-Savart kernel for the Poisson equation.
The Navier-Stokes equations can be expressed in a Lagrangian formulation,
leading to a set of equivalent ODEs as:
dx p
= u(x p , t) (16)
dt
dω p
= ∇u(x p , t) ω p + ν ω(x p ), (17)
dt
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

where x p , ω p denote the locations and the vorticity carried by the fluid elements.
The equations need to be supplied with initial and far-field conditions along with
the no-slip condition in the presence of solid boundaries.
The essence of the VMs we describe herein is based on the work of Krasny
(1986) and it amounts to the regularization of the convecting velocity field and the
systematic removal of spurious vortical structures.
VMs are based on this Lagrangian description and use vorticity-carrying parti-
cles with a finite core size ε so the vortex-blob approximation is given by


ω εh (x) = v p ω p ζε (x − x p ), (18)
p

where x p , v p , and ω p , respectively, represent the locations, volumes, and vorticity

of the particles. In VMs, the field is recovered at every location of the domain
only if one considers the collective behavior of all computational elements. In
addition, when particles overlap, the scales of the physical quantities that are re-
solved are determined by the particle core rather than the interparticle distance.
This observation differentiates particle methods from schemes such as finite differ-
ences. Viscous effects are simulated using the method of PSE (Equation 11). Flows
with solid boundaries are treated using a fractional step algorithm (Chorin 1973),
solving in turn the inviscid and viscous parts of the equations. The enforcement
of kinematic boundary conditions (such as no-through flow at solid boundaries)
can be achieved by boundary integral methods whereas enforcement of viscous
boundary conditions (such as no-slip) is translated into a vorticity flux boundary
condition (Cottet & Poncet 2004, Koumoutsakos et al. 1994) complementing the
viscous part of the equations. This boundary condition is enforced using an integral
formulation, resulting in an explicit modification of the particle weights near the
boundary. It can be formulated by adding a forcing term to the vorticity equation,
amounting to an immersed boundary method (see article by Mittal & Iaccarino in
this volume).
The vortex-blob method can be summarized by the following system of ODEs
for the particle locations and vorticities

dx p
N
= vq K ε (x p − xq ) × ω q + U 0 (x p , t) (19)
dt q=1
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

FLOW SIMULATIONS USING PARTICLES 463

 
dω p
N
= vq ∇ K ε (x p − xq ) × ω q ω p (20)
dt q=1

ν
N
+ vq [ω q − ω p ] ηε (|x p − xq |) + F(x p ), (21)
ε2 q=1

where the term F(x p ) accounts for the generation of vorticity at solid boundaries.
Based on these discretizations, particle methods are unconditionally linearly stable.
Nonlinear stability imposes that particle trajectories do not cross, which results in
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

a time-step constraint of the type

t ≤ C∇u−1
∞ , (22)

where the coefficient C depends on the particular numerical scheme. The sta-
bility properties of VMs make them a suitable candidate for the Heterogeneous
Multiscale Methods (HMM) framework introduced by E & Engquist (2003).
In the past, VMs were used extensively for simulations of engineering appli-
cations, such as unsteady bluff body flows, with reasonable agreement between
computational and experimental results (see Sarpkaya 1989 for a thorough review).
This agreement can be explained by analyzing the error (Leonard 1985) introduced
by smooth VMs, indicating that computations with vortex particles implicitly re-
alize some kind of turbulence modeling (Cottet 1996). Ignoring viscous effects,
smooth VMs amount to solving a mollified form of the Euler equations:
∂ω
+ ∇ · uε ω − (ω · ∇)uε = 0, (23)
∂t
where the overbars denote mollification with a smooth kernel. When compared
to the Euler equation for the fields uε and ω ε , this equation involves a truncation
error with a component proportional to ∇ · ([∇uε ] ∇ω), contributing to enstrophy
transfer between different scales. This term is responsible for the emergence of
microstructures in calculations using VMs, because unlike grid-based methods,
their dynamics is not constrained to any minimal scale beyond the initialization
stage. To remove the backscatter, a natural scheme is to formulate the error term as
an integral operator amounting to anisotropic diffusion and to adjust accordingly
the particle weights to compensate for this term. Alternatively, regularization of
the particle locations can compensate for this error.
In the last two decades we have seen a number of theoretical developments and
benchmark flow simulations using VMs. Direct Numerical Simulations (DNS) of
the flow past an impulsively started cylinder (Figure 1) for a range of Reynolds
numbers (Koumoutsakos & Leonard 1995) have demonstrated that VMs can au-
tomatically adapt computational elements in regions of the flow where increased
resolution is necessary to capture unsteady separation phenomena. The results
obtained by VMs for this flow are in excellent agreement with experimental and
analytical studies. Comparisons in terms of the drag coefficient for this flow for
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

464 KOUMOUTSAKOS
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

Figure 1 Left: Vorticity at T = 6.0 for flow past a two-dimensional impulsively

started cylinder of Re = 9500 using the Vortex Method with a Fast Multipole Method
(courtesy of Koumoutsakos & Leonard 1995). Right: Three-dimensional flow past
a cylinder at Re = 300 that was computed using a vortex-in-cell method. Iso-
surfaces of spanwise and transverse vorticity are shown (courtesy of Cottet et al.
2004).

Re = 1000, with high-order finite difference methods (Anderson & Reider 1996),
show that VMs compare favorably in terms of the number of computational ele-
ments for the same accuracy. Comparisons with spectral element methods (Fischer
1997) for Re = 9500 reveal that the results obtained with VMs can be obtained by
spectral element methods, albeit only when using additional elements in critical
parts of the flow. These parts of the flow are not always known a priori and for grid-
based methods suitable criteria need to be devised to add computational elements
in critical regions depending on the physics of the flow. VMs have the advantage
that computational elements are inherently linked to the physics they represent and
thus no such additional criteria are necessary. Simulations using VMs of flow past
a sphere (Ploumhans et al. 2002) and a 3D cylinder (Figure 1) (Cottet & Poncet
2004) have shown the same advantages for 3D flows.
The use of VMs is particularly advantageous for simulations of controlled flows
involving unsteady boundary motions as the Lagrangian formulation of the convec-
tive transport term enables large time steps. In Eulerian-based methods, because
finer resolution is necessary to capture the vortical structures near the boundaries,
smaller time steps are necessary to obey the transport Courant-Friedrichs-Levy
(CFL) condition as the near wall elements experience large velocities induced by
the motion of the boundary (Poncet 2004).
Simulations of homogeneous turbulence show that energy spectra obtained
using VMs are in excellent agreement with those predicted by spectral element
methods (see Figure 2) (Cottet et al. 2002). These simulations indicate that al-
though VMs are less efficient than spectral methods, their computational cost is
not prohibitive for simulations of homogeneous turbulence.
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

FLOW SIMULATIONS USING PARTICLES 465

by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

Figure 2 Comparison of Vortex and Spectral methods of the evolution of en-

strophy and energy spectrum in simulation of homogenous isotropic turbulence
courtesy of Cottet et al. (2002).

The adaptivity and robustness of VMs has enabled simulations of reacting

flows (Ghoniem & Oppenheim 1984, Knio & Ghoniem 1992) and flows in porous
media (Zimmerman et al. 2001). In the latter case, a comparison with grid-based
methods shows that Lagrangian particle methods perform better than several
finite difference methods in capturing highly anisotropic diffusion phenomena
while accurately transporting scalar fields. Simulations of particle-laden flows (see
Figure 3) (Walther & Koumoutsakos 2001) reveal flow structures and flow insta-
bilities that were predicted experimentally but had not been obtained before by
grid-based methods, possibly due to the absence of adaptivity and the dissipation
induced by the discretization of the nonlinear transport term. The method has been
extended to compressible flows (Eldredge et al. 2001), but its relevant advantages
in this field are a subject of ongoing investigations.
It has been well known since the works of Krasny in the 1980s (Krasny 1986)
that particle methods are well suited for interface capturing. Level sets present
today the standard framework to capture interfaces (Osher & Fedkiw 2001) and
a particle level set formulation has been implemented (Enright et al. 2002) to
remedy some problems involved in the evolution of a level set on a fixed grid.
Recently, a novel particle level set method for capturing interfaces was proposed
(Hieber & Koumoutsakos 2004). In this method, the level set equation is solved in a
Lagrangian frame using particles that carry the level set information. A key aspect
of the method involves a consistent remeshing procedure for the regularization
of the particle locations. This Lagrangian description of the level set method is
inherently adaptive and exact in the case of solid body motions. Comparisons on a
set of benchmark problems with existing level set formulations demonstrates that
the proposed particle-level set method achieves superior results using a reduced
number of computational elements.
A detailed description of particle methods with an emphasis on VMs and some
of their applications can be found in the monograph by Cottet & Koumoutsakos
(2000).
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

466 KOUMOUTSAKOS
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

Figure 3 Simulation of particle-laden flow. Vorticity isosurfaces (red and blue) and
solid particles (white and yellow) for a drop of solid particles falling in a fluid with
zero (left) and nonzero (right) initial vorticity field (Walther & Koumoutsakos 2001).

2.3. Smooth Particle Hydrodynamics

The method of SPH was introduced by Lucy in the late 1970s and was further
developed by Monaghan (see Monaghan 1988 and references therein) for grid-
free astrophysics simulations.
In SPH, function approximations are consistent with Equation 4 and particle
weights are selected as w Sp P H = f (x p ) vq = f (x p )m p /ρ(x p ) by invoking the
continuity equation and implicitly bypassing the requirement for an exact calcu-
lation of the volume associated with each particle. In SPH, the key requirements
for the particle kernel are positivity and local support, and the scheme relies in the
conservative approximation of derivative operators using


N
 
D β f (x p ) = f q − f p vq D β W (x p − xq , h), (24)
q

where W (x p − xq , h) is used instead of the mollifier kernel ζε employed in VMs,

with the interparticle distance h taking the role of the mollifier core size. Using
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

FLOW SIMULATIONS USING PARTICLES 467

Equation 24, the continuity and momentum equation can be expressed in the SPH
formulation as

dx p
= up
dt
dρ p
 
= vq uq − u p · ∇W (x p − xq , h)
dt q

du p 1

= vq τ −τ · ∇W (x p − xq , h) + F, (25)
ρp q
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.

dt q p
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

where τ denotes the stress tensor of the flow and F corresponds to external force
fields experienced by the particles. A closure relationship is necessary to express
the stress tensor as a function of known variables. In the past 20 years many
simulations using SPH have been conducted, extending its application range from
gas dynamics in astrophysics to Newtonian and viscoelastic flows (see Ellero et al.
2002, Monaghan 1985b and references therein).
Several open questions remain regarding the enforcement of boundary condi-
tions and the consistency of the method in situations of highly distorted parti-
cle configurations. The particle distortion leads to errors in the approximation of
derivative operators (Belytschko et al. 1996). Flow simulations using SPH involve
an implicit subgrid-scale modeling (although currently no analysis exists for this),
and suitable corrections are necessary to enhance the accuracy of the method.
Several techniques [such as artificial viscosity (Gingold & Monaghan 1983) and
dynamic conditions (Ellero et al. 2002)] have been proposed to compensate for this
problem. Inspired by techniques in VMs, the introduction of regularization of par-
ticle distortion in SPH via remeshing (Chaniotis et al. 2002) has led to second-order
accuracy, but this detracts from the characterization of the method as grid-free.
Recent work has focused on the relationship between SPH and particle methods
developed for solving boundary value problems. These so-called meshless meth-
ods, in order to be distinguished from schemes like finite differences and finite
elements where node connectivity is important, are Galerkin-type methods. They
compute the approximations of derivative operators by solving systems of equa-
tions to construct conservative, particle-based, discrete mollifier kernels. Works
by Duarte & Oden (1996), and Belytschko et al. (1996) provide a unifying frame-
work for methods such as Moving Least Squares, Reproducing Kernel Particle
Methods, and Element-Free Galerkin, and discuss their relationship with SPH. We
refer to review articles by Belytschko, and more recently by Babuska et al. (2002),
on the developments of these methods and their formulation as Partition of Unity
Methods. Meshless methods have been mostly implemented in Galerkin formu-
lations of solid mechanics problems, but their developments carry a number of
concepts, such as multiscale particle representations and formulation of accurate
boundary conditions, that should be further explored to increase the capabilities
of Lagrangian particle methods for flow simulations.
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

468 KOUMOUTSAKOS

3. GRIDS AND PARTICLES

Particle methods are often defined as grid-free methods, making them an attrac-
tive alternative to mesh-based methods for flows past complex and deforming
boundaries. However, the adaptivity provided by the Lagrangian description can
introduce errors and particle methods have to be conjoined with a grid to pro-
vide consistent, efficient, and accurate simulations. The grid does not detract from
the adaptive character of the method and serves as a tool to restore regularity in
the particle locations via remeshing while it simultaneously enables systematic
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.

multiresolution particle simulations (Bergdorf et al. 2004), allows fast-velocity

Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

evaluations (Harlow 1964), and facilitates hybrid particle mesh methods capable
of handling different numerical methods and different equations in various parts
of the domain (Cottet 1990).

3.1. Remeshing for Particle Distortion

Particle methods, when applied to the Lagrangian formulation of convection-
diffusion equations, enjoy an automatic adaptivity of the computational elements
as dictated by the flow map. This adaptation comes at the expense of the regularity
of the particle distribution because particles adapt to the gradients of the flow field.
The numerical analysis of VMs shows that the truncation error of the method is
amplified exponentially in time, at a rate given by the first-order derivatives of the
flow that are precisely related to the amount of flow strain. In practice, particle dis-
tortion can result in the creation and evolution of spurious vortical structures due to
the inaccurate resolution of areas of high shear and to inaccurate approximations
of the related derivative operators.
To remedy this situation, location processing techniques reinitialize the distorted
particle field onto a regularized set of particles and simultaneously accurately
transport the particle quantities. The resulting problem of extracting information
on a regular grid from a set of scattered points has a long history in the fields
of interpolation (Schoenberg 1946) and statistics (see Cleveland & Loader 1996
and references therein). To facilitate the analysis we restrict our attention to a 1D
equispaced regular grid with unit mesh size onto which we interpolate quantities
(qn ) from scattered particle locations (xn ):


Q(x) = qn W (x − xn ). (26)
n

The properties of the interpolation formulas can be analyzed through their behavior
in the Fourier space (Schoenberg 1946). The characteristic function g(k) of the
interpolating function W(x) is defined as
+∞
g(k) = W (x) e−ikx d x.
−∞
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

FLOW SIMULATIONS USING PARTICLES 469

When W decays fast at infinity, g is a smooth function and the interpolation formula
Equation 26 is of degree m if the following two conditions hold simultaneously:
(a) g(k) − 1 has a zero of order m at k = 0 and (b) g(k) has zeros of order m at all
k = 2π n, (n = 0). These requirements translated back in the physical space are
nothing but the moment properties of the interpolant
 
W (y) dy = 1; y α W (y) dy = 0, if 1 ≤ |α| ≤ m − 1.

This is reminiscent of the conditions for accurate function particle approximations

by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.

using moment conserving kernels. In fact, the interpolation accuracy (Hockney

Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

& Eastwood 1988) can be described by splitting the interpolation error into a
convolution and sampling error reminiscent of the smoothing/quadrature error for
function approximations. Hence, good interpolation schemes are those that are
band-limited in the physical space and are simultaneously close approximations
of the ideal low-pass filter in the transformed space. Monaghan (1985b) presents
a systematic way of increasing the accuracy of interpolating functions, such as
B-splines, while maintaining their smoothness properties using extrapolation. He
constructs interpolation formulas such that, if m = 3 or m = 4, the interpolation
will be exact for quadratic functions, and the interpolation will be third- or fourth-
order accurate. One widely used formula involves the so-called M4 function


 0 if |x| > 2

M4 (x) = 12 (2 − |x|)2 (1 − |x|) if 1 ≤ |x| ≤ 2 (27)


 5x 2 3|x|3
1− 2 + 2 if |x| ≤ 1.

Interpolations in higher dimensions can be achieved by tensorial products of these

formulas. However, these tensorial products require particle remeshing on a regu-
lar grid. For non-grid-conforming boundaries, remeshing introduces particles onto
areas that are outside the flow domain and violates the flow boundary conditions.
Remedies such as one-sided interpolation have been proposed and a working solu-
tion can be obtained (Cottet & Poncet 2004, Ploumhans et al. 2002) by eliminating
particles outside the domain and adjusting accordingly the modification of particle
strengths by re-enforcing the boundary conditions in a fractional step algorithm.
Alternatively, weight processing schemes attempt to explicity (Beale 1986) or im-
plicitly (Strain 1997) modify the particle weights in order to maintain the accuracy
of the calculation, but they result in rather costly calculations.

3.1.1. HYBRID METHODS Hybrid methods involve combinations of mesh-based

schemes and particle methods in an effort to combine computational advantages
of each method. The first such method involves the Particle in Cell algorithm
pioneered by Harlow (1964), in which a particle description replaces the nonlinear
advection terms and mesh-based methods can be used to take advantage of the
efficiency of Eulerian schemes to deal with elliptic or hyperbolic problems.
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

470 KOUMOUTSAKOS

Lagrangian-Eulerian domain decomposition methods use high-order grid meth-

ods and VMs in different parts of the domain (Cottet 1990, Ould-Salihi 2000) and
can even be combined with different formulations of the governing equations. A
finite difference scheme (along with a velocity-pressure formulation) can be im-
plemented near solid boundaries, and VMs (in a velocity-vorticity formulation)
can be implemented in the wake to provide the flow solver with accurate far-field
conditions. In this approach Eulerian methods handle the wall boundary conditions
and can be complemented with immersed boundary methods (Mittal & Iaccarino
2005) to handle complex geometries. A rigorous framework for particle-based im-
mersed boundary methods has been developed based on a unified formulation of
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

the equations for flow-structure interaction (Cottet 2002). Simulations involving

this formulation are a subject of ongoing investigations.

4. MULTIRESOLUTION PARTICLE METHODS

The accuracy of smooth particle methods with overlapping cores is determined by
the core size ε of the mollifier. For computational efficiency this core size needs to
be spatially variable to adequately discretize gradients in different parts of the flow,
such as the boundary layer and the wake of bluff body flows. Because particles
must overlap spatially, varying cores imply a corresponding adaptation for the
spacing of the particles. This can be achieved by remeshing the particle locations
on a spatially varying mesh by
■ remeshing on a regular grid corresponding to variable size particles by using
a global (adaptive or nonadaptive) mapping and by
■ remeshing by combining local mappings in a domain decomposition frame-
work.
In VMs, Hou (1990) first introduced a variable-size VM for the 2D Euler
equations by defining a function ε(x) << 1 for the vortex particles so that


ωεh (x) = v p ω p ζε(x p ) (x − xhp ). (28)
p

This leads to a spatially varying mollified velocity kernel for the advancement of
vortex particles


uh (x) = v p ω p Kε(x p ) (x − xhp ). (29)
p

The convergence of the method was proven for the Euler equations under the
assumption that there is a positive, bounded, smooth function Fsuch that ε(x) =
ε F(x).
The straightforward extension of this method to viscous flows by the modifica-
tion of the kernel in Equation 11 leads to an inconsistent approximation. To avoid
this inconsistency, we need to assume a mapping from the physical coordinates,
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

FLOW SIMULATIONS USING PARTICLES 471

with variable-size blobs, to a coordinate system where blobs have a uniform size,
as presented by Cottet et al. (2000). The algorithm involves the mapping of a do-
main  ˆ with uniform blobs, to a domain , where we wish to use variable blob
sizes through a mapping denoted by F:

x = F(x̂); x̂ = G(x); ω(x) = ω̂(x̂).

If J = det[ai j ] = det[ ∂G i
∂x j
] denotes the Jacobian determinant of the inverse map-
ping G, we can express the Laplacian operator in the variable core domain, in
terms of gradient operators in the uniform domain, as:
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

x ω = J divx̂ [B∇x̂ ω̂], (30)


where B is the matrix with entries b jk = J −1 i aki a ji .
Integral approximations for the differential operator in the mapped coordinates,
in the righthand side of Equation 30, can be derived by using Equation 12 with
ψi j = xi x j θ (|x|), where the spherically symmetric kernel θ is normalized such
that xi4 θ (x) dx = d + 2, where d denotes the dimension of the problem. The PSE
scheme for the Navier-Stokes equation that follows from these formulas is given
as
dω p

= νε −4 J (x p ) v̂q (x̂ ip − x̂qi )(x̂ pj − x̂qj )θ ε (x̂ p − x̂q )
dt q,i, j
 
1
x̂ p + x̂q
× bi j − bii δi j (ω q − ω p ). (31)
d +2 i 2

The scheme is conservative because the volumes of the particles in the physical
and mapped spaces are related through v p J (x p ) = v̂ p .
The PSE scheme in Equation 31 relies on the explicit knowledge of a global
invertible mapping in the computational domain. This is easily accomplished for
flows in geometries such as channels and cylinders where such mappings can be
derived. However, for more complex geometries or geometries involving several
bodies such mappings are not available. Similar to r-adaptive mesh-based methods
(Ceniceros & Hou 2001), it is desirable to require enhanced resolution for vortex
particles in areas of high shear. One way to achieve this is to construct adap-
tive maps. Bergdorf et al. (2004) introduced global adaptive mappings through a
particle approximation of a differential and continuous map


M
x(x̂, t) = F(x̂, t) = µ j (t) ϕ j (x̂) . (32)
j=1

The parameters in the map that are changed in the process of adaptation are the
node values {µ j } Mj=1 . Using a map, as described in Equation 32, makes it im-
possible to leap back and forth from physical to reference space. However, its
differentiability enables casting the governing equations into reference space and
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

472 KOUMOUTSAKOS

solving the problem there without needing the inverse map. The adaptivity of the
map presents us with an extra degree of freedom because complementing the con-
vection of the particles with velocity u in physical space, particles are adapted by
the convection/adaptation of the map with a specified velocity
dµi
= U(x̂, t). (33)
dt
The method was implemented successfully for 1D (Burger’s equation) and 2D
(inviscid axisymmetrization of an elliptical vortex) problems (Figure 5). Using
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.

a suitable map velocity that implies that the core size in the physical space is
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

being deformed in the same way the volume is deformed by the flow, the adaptive
particle overlap maintains enhanced resolution in areas of high gradients where
the particles are compressed. This multiresolution VM amounts to specifying rules
for deformation of the particle shapes proportional to the spatially varying spacing
of the particles.
A more versatile technique to achieve multiresolution in particle methods in-
volves the combination of several local mappings with variable blobs associated to
each mapping linked through domain decomposition techniques (Bergdorf et al.
2004, Cottet et al. 2000). In a domain decomposition algorithm involving particle
solvers, interface conditions must be supplied to compute particle velocities and to
update particle strengths. Given a vorticity field, determining the velocity amounts
to solving a Poisson equation, and the Schwarz alternating method is a natural
way to enforce the right interface conditions. In a vortex code, the method iterates
between the boundary source terms that must be added to the Biot-Savart law in
each subdomain. The fact that variable blobs, corresponding to different mappings,
are implemented does not introduce any difficulty, provided the overlapping zone
has a width exceeding the blob sizes in this area.
Once the velocities are evaluated on each subdomain, particle locations and
circulations must be updated. Because VMs are based on explicit time discretiza-
tion of the vorticity convection-diffusion equation, the vorticity transfer from one
domain to another is simply achieved by interpolating the particle strengths in the
overlapping zone.
The essence of the algorithm is independent of the dimension and the geometry
of the domain. To illustrate the algorithm, let us consider the example sketched
in Figure 4. In this example, 1 and 2 are two domains with nonuniform grid
spacing, overlapping with a domain 3 with uniform spacing. During the remeshing
step, starting from a distorted particle configuration in the buffer zone,
■ particles of type 1 and 2 are remeshed onto particles of type 3, using the
cartesian mapping;
■ similarly, particles of type 3 are remeshed onto particles of type 1 and 2 with
polar mappings;
■ and particles are advanced and exchange vorticity in each sub domain via
PSE to account for diffusion.
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

FLOW SIMULATIONS USING PARTICLES 473

by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

Figure 4 Domain Decomposition for flow past two cir-

cular cylinders.

At the end of these steps, the vorticity has been updated in all three domains.
Provided that the overlapping width exceeds the core radius of the PSE kernel, this
procedure allows a consistent transfer of vorticity through diffusion.
In variable blob methods remeshing must be applied in the mapped coordinates
and at a frequency that prevents particles in low-resolution areas to travel too far
in the high-resolution areas between two remeshing steps. Note that the remeshing
of particle methods provides a flexible method for multiresolution representations
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

474 KOUMOUTSAKOS
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

Figure 5 Inviscid simulation of an elliptical vortex. (Bottom left) Initial condition.

(Bottom right) Vorticity field at T = 1.5. (Top right) Close-up of vorticity distribution
as computed by an Adaptive Vortex Method based on Adaptive Mesh Refinement
(AMR). (Top left) Close-up of vorticity distribution as computed by an r-adaptive
Vortex Method based on Adaptive Global Mappings (AGM). (Courtesy of Bergdorf
et al. 2004.)

because it is compatible with hierarchical mesh refinement (Guskov et al. 2002).

When the mesh nodes are translated into particles this provides a flexible mul-
tiresolution representation with no particular restrictions on the connectivity of
the computational elements.
Besides multiresolution particle methods, as described above, a recent work
by Hou (2004) discusses the application of particle-based methods in multiscale
modeling of incompressible flows using ideas from homogenization along with a
formulation for the Lagrangian transport of small scales. The extension of con-
cepts first developed for grid-based methods, such as homogenization, to particle
methods may offer a promising avenue for constructing adaptive and consistent
multiscale formulations.
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

FLOW SIMULATIONS USING PARTICLES 475

5. PARTICLE METHODS FOR DISCRETE SYSTEMS

Discrete models such as molecular dynamics (MD) and dissipative particle dy-
namics (DPD) can describe flows for which the macroscale description through
the Navier-Stokes equations is not adequate. In addition, they can complement
macroscale descriptions by providing a suitable description of fluctuating hydro-
dynamics and by elucidating phenomena such as the validity of the no-slip con-
dition for which only empirical evidence exists. In these methods there is a direct
correspondence between the computational particles and the structures (molecules,
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.

fluid particles) that model the behavior of the system.

Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

5.1. Molecular Dynamics and Fluids

MD simulations are used to model fluids (gas, liquid) characterized by the time
and length scales of molecular motion. MD amounts to computing the trajectories
of particles interacting through classical simplified force fields.
The simulation of fluids in MD dates back to the inception of the method in the
mid-1950s works of Fermi et al. (1955) and Alder & Wainwright (1957), in which
the phase diagram of a hard sphere system was investigated. In 1971, Rahman &
Stillinger (1971) reported the first simulations of liquid water. Over the years, the
development of suitable force fields propelled the current method to become one
of the key interdisciplinary computational tools for investigating large molecular
systems. Please see the book by Schlick (2002) for a comprehensive description.
Today, MD simulations are increasingly popular in the field of fluid mechanics
and in the last decade several review articles have appeared, starting with the work
of Koplik & Banavar (1995), who discussed the formulation of continuum flow
deductions from atomistic simulations. Micci et al. (2001) reviewed nanoscale
flow phenomena related to atomization and sprays.
Today, with increased computational powers, computations involving flows of
complex fluids such as water are becoming routine. MD simulations of water have
elucidated a number of issues associated with fluid mechanics of wetting and
hydrophobicity (Figure 6) at the nanoscale (Walther et al. 2004), water transport
through carbon nanotubes (Hummer et al. 2001), and, particularly, with the validity
of the no-slip boundary condition for water flows in various nanoscale geometries
(Sokhan et al. 2002). Recent MD studies (Walther et al. 2004) of water flows past
carbon nanotubes reveal that the validity of the no-slip boundary condition at the
nanoscale depends not only on the fluid and surface material properties but also
on the particular geometric configuration. Extensive reviews of nanoscale fluid
mechanics can be found in Maruyama (2001) and Koumoutsakos et al. (2004).

5.1.1. MOLECULAR DYNAMICS: FORCE FIELDS AND POTENTIALS The potential en-
ergy function or force field provides a description of the relative energy or forces
of the ensemble for any geometric arrangement of its constituent atoms. This de-
scription includes energy for bending, stretching, and vibrations of the molecules,
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

476 KOUMOUTSAKOS
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

Figure 6 Hydrophobic hydration of carbon nanotubes

in water (from Walther et al. 2001).

and nonbonded interaction energies between the molecules. Classical force fields
are usually built up by the superposition of simple potential energy expressions.
Mostly pair potentials V(ri j ) are used, but in systems where bonds are determining
the structure, multibody contributions V(ri j , rik ), and V(ri j , rik , ril ) may also enter
the expression, thus




U= V (ri j ) + V (ri j , rik ) + V (ri j , rik , ri,l ), (34)
i, j i, j,k i, j,k,l

where ri j = |ri − rj| is the distance between i-th and j-th atoms. The contribution
to the interaction potential can be ordered in two classes: intramolecular and in-
termolecular contributions. Whereas the former describe interactions that arise in
bonded systems, the latter are usually pair terms between distant atoms modeling
electrostatics and Van der Waals interactions.
The study of nonequilibrium processes or dynamic problems, such as fluid
flows in nanoscale geometries, is usually performed by nonequilibrium molecular
dynamics (NEMD). NEMD is based on the introduction of a flux in thermodynamic
properties of the system (Allen & Tildesley 1987). Cummings & Evans (1992)
review NEMD with regard to the computation of transport coefficients of fluids
from the knowledge of pair interactions between molecules. Ryckaert et al. (1989)
compare the performance of NEMD with Green-Kubo approaches to evaluate the
shear viscosity of simple fluids. Tuckerman et al. (1997) present a modified NEMD
approach to ensure energy conservation.

5.1.2. BOUNDARY CONDITIONS FOR MOLECULAR DYNAMICS In multiscale simula-

tions, the MD part of the flow has to interface mesoscopic and macroscale models
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

FLOW SIMULATIONS USING PARTICLES 477

through suitable boundary conditions. For situations involving the simulation of

a solvent, the small volume of the computational box in which solvent and other
molecules of interest are contained can introduce undesirable boundary effects if
the boundaries are modeled as simple walls. To circumvent this problem, the sys-
tem may be placed in vacuum (Allen & Tildesley 1987) or a periodic system may
be assumed. However, periodic boundary conditions imposed on small systems
may introduce artifacts in systems that are not inherently periodic.
Stochastic boundary conditions enable reduction of the size of the system by
partitioning the system into two zones with different functionality: a reaction zone
and a reservoir zone. The reservoir zone is excluded from MD calculations and
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

is replaced by random forces whose mean corresponds to the temperature and

pressure in the system. The reservoir zone is further subdivided into a reaction
zone and a buffer zone. The stochastic forces are only applied to atoms of the
buffer zone. Please see Brunger et al. (1984) and Berkowitz & McCammon (1982)
for the application of stochastic boundary conditions to a water model.

6. CONTINUUM-MOLECULAR SIMULATIONS
Nanoscale flows are often part of larger-scale systems (for example, when nanoflu-
idic channels are interfacing microfluidic domains), and in simulations we are
confronted with an inherently multiscale problem. The simulation of such flows is
challenging because one needs to suitably couple the nanoscale systems with larger
spatial and time scales. The macroscale flows determine the external conditions
that influence the nanoscale system, which in turn influences the larger scales by
modifying its boundary conditions.
Hybrid computational techniques attempt to overcome these problems using,
for the molecular part of the flow, Direct Simulation Monte Carlo (DSMC) for
dilute gases or MD for liquids, coupled with a relevant continuum description
(Flekkøy et al. 2000, Hadjiconstantinou 1999a, Hadjiconstantinou & Patera 1997,
Li et al. 1999, Sun et al. 2004, Werder et al. 2004). An alternative is coarse-grained,
mesoscopic particle models of complex fluids, such as DPD, which do not involve
any explicit coupling of the continuum and the atomistic description.

6.1. Dissipative Particle Dynamics

The initial formulation of the DPD model was given by Hoogerbrugge & Koelman
(1992) and was intended to provide a mesoscale model that enables the simulation
of complex fluids such as colloidal suspensions, emulsions, polymers, and multi-
phase flows. It is based on the notion of fluid particles representing a collection of
atoms or molecules that constitute the fluid. These fluid particles interact pairwise
through three types of forces,


fi = FC (ri j ) + F D (ri j , ui j ) + F R (ri j ) , (35)
j=i
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

478 KOUMOUTSAKOS

where FC represents a conservative force derived from a soft repulsive potential

and F D is a dissipative force depending on the relative particle velocity ui j to model
friction whereas the stochastic force F R models the effect of the suppressed degrees
of freedom in the form of thermal fluctuations.
DPD simulations of complex fluids offer advantages in two respects when
compared to MD. First, the conservative pairwise forces between the DPD par-
ticles are soft repulsive, which makes it possible to extend simulations to longer
timescales, whereas coarse graining in the particle representation allows studies of
larger systems. Second, a special “DPD thermostat” for the canonical ensemble is
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.

implemented in terms of dissipative as well as random pairwise forces such that the
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

momentum is locally conserved, which results in the emergence of hydrodynamic

flow effects on the macroscopic scale. Today, one can reach simulation times of
the order of 100 ns with molecular dynamics, whereas one can routinely study
phenomena on the microsecond scale with DPD. A drawback of the DPD method
is that its thermodynamic behavior is determined by the conservative forces and
is therefore an output of the model and not an input (Serrano & Espanol 2001).
Espanol & Revenga (2003) recently introduced the smoothed dissipative particle
dynamics method (SDPD) starting from a formulation of SPH. In these simulations
every particle has an associated position, velocity, constant mass, entropy and, in
addition, two extensive variables, a volume and an internal energy. The interpolant
used in the SDPD formulation fulfills the second law of thermodynamics explicitly
and thus enables the consistent introduction of thermal fluctuations through the
use of the dissipation-fluctuation theorem.
DPD is particularly well suited for simulations of polymers and surfactant
systems and the reader is referred to recent reviews on mesoscale simulations of
complex fluids using DPD (Glotzer & Paul 2002, Kremer & Müller-Plathe 2002,
Warren 1998).

6.2. Multiscaling: Linking Macroscopic to Atomistic Scales

The Navier-Stokes equations are based on classical Newtonian mechanics and rely
on the continuum approximation as well as on the assumption of thermodynamic
(quasi-) equilibrium. The continuum approximation relies on the formulation of
local flow properties such as density, velocity, and stress as averages over fluid
elements. Thermodynamic (quasi-) equilibrium postulates that when equilibrium
is achieved within small volumes for certain fluid properties, their gradients vary
linearly between these volumes, hence the stress is linearly related to the strain
and the heat flux is linearly related to the temperature gradient. For dilute gas flows
the conditions under which the above assumptions hold are well characterized by
the degree of rarefaction of the fluid measured in terms of the Knudsen number,
Kn = λ/L, where λ is the mean-free path and L a characteristic macroscopic
length (Bird 1994). A flow with Kn < 0.01 is in the continuum regime and can be
well described by the Navier-Stokes with no-slip boundary conditions. For 0.01 <
Kn < 0.1, the slip-flow regime, the Navier-Stokes equations can still be used
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

FLOW SIMULATIONS USING PARTICLES 479

along with tangential slip-velocity boundary conditions. In the transition regime,

for 0.1 < Kn < 10, the constitutive equation for the stress tensor starts to lose its
validity and higher-order corrections such as the Burnett or Woods equations along
with higher-order slip models at the boundary are needed. At even larger Knudsen
numbers (Kn > 10), the continuum assumption fails completely and atomistic
descriptions such as DSMC of the gas flow are needed (Bird 1994). Note that as
the considered system size L shrinks, the thermodynamic equilibrium assumptions
fails before the continuum approximation does (Karniadakis & Beskok 2002).
For liquid flow, the situation is more complicated because the molecules that
constitute a liquid are essentially always in a collision state and, thus, the con-
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

cepts of a mean-free path and Knudsen number are not useful. There is no well
established molecular-based theory for liquids, such as for dilute gases. Therefore,
one must resort to experiments or to a computational analysis using MD, where
fluids are modelled as what they really are—a collection of strongly interacting
molecules.

6.3. Hybrid Atomistic-Continuum Computations

O’Connell & Thompson (1995) described an early attempt to extend the length
scales accessible in molecular dynamics simulations through the combination with
a continuum descriptions. In their simulations, they applied constrained dynamics
in an overlap (X) between the particle (P) and continuum (C) regions to ensure
stress continuity across the P – C interface. O’Connell & Thompson applied this
algorithm to an impulsively started Couette flow where the P – C interface was
parallel to the walls. This ensured that there was no net mass flux across the
MD-continuum interface. Hadjiconstantinou (Hadjiconstantinou 1999, Hadjicon-
stantinou & Patera 1997) pointed out that this scheme decouples length but not
time scales and therefore suggested an iterative procedure based on the Schwarz
alternating method to alleviate this problem. In this iteration, the continuum so-
lution in C provides boundary conditions for a subsequent atomistic solution in P
and vice versa until the solution converges in the overlap region X. The Schwarz
method is inherently bound to steady-state problems; however, for cases with
the hydrodynamic time scale much larger than the molecular time scale, a se-
ries of quasi-steady Schwarz iterations can be used to treat transient problems
(Hadjiconstantinou & Patera 1997). A hybrid formulation of the moving con-
tact line problem served as a test problem. Flekkøy et al. (2000) presented a hy-
brid model which, in contrast to earlier hybrid schemes (Hadjiconstantinou 1999,
O’Connell & Thompson 1995), is explicitly based on direct flux exchange be-
tween the particle and the continuum region. The main difficulty in the approach
of Flekkøy et al. arises in the imposition of the flux boundary condition from the
continuum region on the particle region. The fluxes exhibit significant oscillations
that do not enable fast convergence of the iterative scheme. The scheme was tested
for a 2D Lennard-Jones fluid coupled to a continuum region described by the com-
pressible Navier-Stokes equations. The first test was a Couette shear flow parallel
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

480 KOUMOUTSAKOS

to the P – C interface, and the second test involved a Poiseuille flow where the
flow direction was perpendicular to the P – C interface. They obtained results on
velocity profiles comparable with those obtained by full-scale MD simulations.
Wagner et al. (2002) extended this work to include the energy equation and applied
the technique to flow in a channel. Garcia et al. (1999) proposed a coupling of a
DSMC solver embedded within an adaptive compressible Navier-Stokes solver
to study gas flows. They successfully tested their scheme on systems such as an
impulsively started piston and flow past a sphere.
Li et al. (1998) introduced a method called thermodynamic field estimator
to extract continous fields from particle data based on the concept of maximum
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

likelihood inference. This so-called thermodynamic field estimator method is sub-

sequently used as detector in a feedback loop implemented to impose boundary
conditions in hybrid schemes of the Schwarz iteration type (Li et al. 1999). The
desired boundary condition is obtained
 by resetting the particle velocities in a
buffer region, such as to minimize i |vi − vi |2 , where vi and vi denote the par-
ticle velocities before and after the transformation. An additional buffer layer is
introduced in between the action region and the overlap region to relax the effect
of the artificial disturbance.

6.3.1. DOMAIN DECOMPOSITION ALGORITHMS Recently, Hadjiconstantinou et al.

(2003) derived an estimate for the number of independent samples needed in
molecular systems to obtain a fractional error E q in a quantity q that is measured
in a domain of volume V . The study quantified the fact that the cost in terms
of sampling time to measure fluxes (such as the momentum flux) is orders of
magnitude larger than the cost for densities. Based on this observation which favors
the use of a density-based scheme over flux-based schemes, Werder et al. (2004)
proposed a hybrid multiscale algorithm for simulating dense fluids using argon
molecules to describe the molecular system. The algorithm uses the alternating
Schwarz method (cf., Hadjiconstantinou 1999), and couples a MD and continuum
fluid dynamics system by interfacing the molecular system with a finite volume
mesh covering the continuum computational domain (Figure 7). In contrast to
Hadjiconstantinou’s work (1999), the scheme is not limited to the use of periodic
systems in the molecular part, but enables more general boundary conditions.
These boundary conditions are implemented by a combination of specular walls,
a potential of mean force (Werder et al. 2004), and a particle insertion algorithm
(Delgado-Buscalioni & Coveney 2003).
Using the flow of argon past a carbon nanotube as a test case, the mehod con-
verges in a few iterations, and the drag coefficient is within 30% of the continuum
Stokes-Oseen flow past a circular cylinder for the associated Reynolds number.
In blending MD and continuum simulations it is necessary that a reasonable
equilibration is achieved in the MD part to lead to a convergent algorithm. This
imposes length and time limitations on the MD computations in order to obtain
reasonable average quantities in the interface with the continuum (Werder et al.
2004). One possibility to overcome this problem is to add a mesoscopic domain
between the molecular and macro domains that can accomodate in a systematic
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

FLOW SIMULATIONS USING PARTICLES 481

by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

Figure 7 (a) Computational domain for the reference solution of the flow of ar-
gon around a carbon nanotube using a purely atomistic description. (b) Hybrid atom-
istic/continuum computational domain. Both computational domains have an extent
of 30 nm × 30 nm. (c) Velocity field for the reference solution averaged over 4 ns.
The white lines are streamlines, and the black lines are contours of the speed |u|.
(d) Velocity field of the hybrid solution after 50 iterations. The black square denotes
the location of the atomistic domain. The solution in the atomistic domain is averaged
over 10 iterations. (Courtesy of Werder et al. 2004.)

manner the fluctuations of the molecular system while simultaneously leading to

a consistent approximation of the macroscale system. Ongoing studies involve
the master equation approach introduced by Breuer & Petruccione (1992) for
fluid dynamics. In this approach the fluid is regarded as a stochastic dynamical
system. The velocity of the fluid is related to a stochastic process governed by
an appropriate master equation (van Kampen 1981) acting in a discrete phase
space. The master equation is constructed in such a way that the average of the
velocity field obeys the underlying Navier-Stokes equations (Breuer & Petruccione
1995). Using the master equation formalism presents several advantages because
boundary conditions are easy to implement, the method is robust with respect to
initial conditions, and it is easily parallelizable. Because the interfacing conditions
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

482 KOUMOUTSAKOS

with the macroscale involve velocity or vorticity boundary conditions, it provides

a natural complement to macroscale simulations using vortex particle methods.

7. FAST PARTICLE METHODS

Particle methods constitute an N-body problem with a computational cost that
scales as O(N 2 ) for N particles. Although short-range forces can be calculated
using cutoffs, the most time-consuming aspect in particle simulations is accurately
evaluating the long-range interactions associated with the field equations. A mesh
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.

can be used to solve the field equations (such as the Poisson equation) that pervade
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

the whole space. There are efficient algorithms to reduce the computational cost,
ranging from simple sorting as first implemented by Verlet (1998) to accurate fast
summation techniques such as Ewald summation (Ewald 1921), the Particle-Mesh
Ewald (PME) method (Darden et al. 1993), and the particle-particle particle-mesh
technique (P3 M) (Hockney et al. 1973) to account for particles in close proximity in
terms of the grid spacing. The nominal cost of Ewald summation requires O(N 1.5 )
operations, the PME and P3 M techniques scale as O(N log N ).
In the last 20 years several mesh-free techniques based on the concept of mul-
tipole expansions have been introduced that circumvent the need for simulating
periodic systems and have minimal numerical dissipation. Examples of such meth-
ods include the Barnes-Hut algorithm (Barnes & Hut 1986), the Fast Multipole
Method (FMM) (Greengard & Rokhlin 1988), and the Poisson Integral Method
(PIM) (Anderson 1992). These methods employ clustering of particles and use
expansions of the potentials around the cluster centers with a limited number of
terms to calculate their far-field influence onto other particles. These techniques
rely on tree data structures to achieve computational efficiency. The tree allows a
spatial grouping of the particles, and the interactions of well-separated particles is
computed using their center of mass or multipole expansions for the Barnes-Hut
and FMM algorithms, respectively. Another advantage of using tree-data structures
is that it allows one to incorporate variable time steps (Mathiowetz et al. 1999) to
integrate the particle trajectories. For a comprehensive review of the treatment of
long-range electrostatics in MD simulations, see Sagui & Darden’s work (1999).

8. SUMMARY AND OUTLOOK

In this review we outline advances in multiresolution particle methods for simulat-
ing continuous systems as well as methodologies for interfacing particle methods
describing molecular systems with the continuum. We try to identify the key char-
acteristics of these methods and, where possible, highlight their advantages and
drawbacks compared to other numerical schemes.
Particle methods for continuum flows offer an accurate, stable, and versatile
method to describe incompressible flow phenomena. In the last decade, a number
of benchmark simulations using VMs have demonstrated that the method is ca-
pable of efficient DNS of turbulent, unsteady separated and interfacial flows. In
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

FLOW SIMULATIONS USING PARTICLES 483

addition, a number of open issues are subjects of ongoing research, including the
development of accurate techniques for handling boundary conditions in the con-
tinuum and atomistic scale as well as in the interface between the two descriptions.
In microscale flows as well as in high Re number flows the use of atomistic models
to describe the near wall part of the flow may be critical in developing consistent
and accurate boundary conditions for the macroscale Navier-Stokes equations that
describe the bulk of the flow.
Multiresolution particle methods for continuum flows complement the inherent
adaptive character of the method, making it a powerful alternative to grid-based
methods. In conjuction with research on efficient techniques for spatially varying
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

particle methods, ongoing research involves the incorporation of adaptive time

integrators and the systematic implementation of particle methods in a rigorous
framework to handle practical flow problems with a very high number of scales that
are not separable. The extension of concepts first developed for grid-based meth-
ods, such as homogenization, to particle methods may offer a promising avenue for
constructing adaptive and consistent multiscale formulations. In the future we may
wish to distinguish multiscale simulations by the number of and the separability of
the scales involved. In these simulations a systematic analysis with rigorous error
control is necessary. Particle methods offer a unique and unifying framework to
formulate multiscale flow phenomena and may serve as a starting point for several
seemingly diverse physical systems in a seamless, interdisciplinary fashion. This
framework would exploit the remarkable and unique features of particle methods,
namely that unlike other numerical methods, such as finite differences and finite
elements, they are fundamentally linked to the physics they aim to reproduce.

ACKNOWLEDGMENTS
I wish to acknowledge many inspiring discussions on particle methods over the
years with Georges-Henri Cottet, Tony Leonard, and Jens Walther. Michael
Bergdorf and Thomas Werder provided invaluable help with the preparation of
this manuscript.

The Annual Review of Fluid Mechanics is online at [Link]

LITERATURE CITED
Alder BJ, Wainwright TE. 1957. Phase transi- explicit method for the computation of flow
tion for a hard sphere system. J. Chem. Phys. about a circular cylinder. J. Comput. Phys.
27(5):1208–9 125:207–24
Allen MP, Tildesley DJ. 1987. Computer Sim- Babuska I, Banerjee U, Osborn J. 2002. Mesh-
ulation of Liquids. Oxford: Clarendon less and generalized finite element methods:
Anderson CR. 1992. An implementation of the a survey of some major results. In Meshfree
fast multipole method without multipoles. Methods for Partial Differential Equations,
SIAM J. Sci. Stat. Comput. 13(4):923–47 M Griebel, ed. pp. 1–20. Berlin: Springer
Anderson CR, Reider MB. 1996. A high order Verlag
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

484 KOUMOUTSAKOS

Barnes J, Hut P. 1986. A hierarchical O(N method for the Navier-Stokes equations. J.
log N) force-calculation algorithm. Nature Comput. Phys. 89:301–18
324(4):446–49 Cottet GH. 1996. Artificial viscosity models
Beale JT. 1986. A convergent 3-D vortex for vortex and particle methods. J. Comput.
method with grid-free stretching. Math. Phys. 127:299–308
Comput. 46:401–24 Cottet GH. 2002. A particle model for fluid-
Belytschko T, Krongauz Y, Organ D, Flem- structure interaction. C.R. Acad. Sci. Paris
ing M, Krysl P. 1996. Meshless methods: an Ser. I 335:833–38
overview and recent developments. Comp. Cottet GH, Koumoutsakos P. 2000. Vortex
Meth. Appl. Mech. & Eng. 139(1–4):3–47 Methods – Theory and Practice . New York:
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.

Bergdorf M, Cottet GH, Koumoutsakos P. Cambridge Univ. Press

Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

2004. Multilevel adaptive particle methods Cottet GH, Koumoutsakos P, Salihi MLO.
for convection-diffusion equations. Multi- 2000. Vortex methods with spatially varying
scale Model. Simul. In press cores. J. Comput. Phys. 162(1):164–85
Berkowitz M, McCammon JA. 1982. Mole- Cottet GH, Michaux B, Ossia S, VanderLin-
cular-dynamics with stochastic boundary- den G. 2002. A comparison of spectral and
conditions. Chem. Phys. Lett. 90(3):215–17 vortex methods in three-dimensional incom-
Bird GA. 1994. Molecular Gas Dynamics and pressible flows. J. Comput. Phys. 175:702–
the Direct Simulation of Gas Flows. Oxford: 12
Clarendon Cottet GH, Poncet P. 2004. Advances in direct
Breuer HP, Petruccione F. 1992. A stochastic numerical simulations of 3D wall-bounded
approach to computational fluid dynamics. flows by Vortex-in-Cell methods. J. Comput.
Continuum Mech. Thermodyn. 4:247–67 Phys. 193(1):136–58
Breuer HP, Petruccione F. 1995. How to build Cummings PT, Evans DJ. 1992. Nonequi-
master equations for complex systems. Con- librium molecular dynamics approaches to
tinuum Mech. Thermodyn. 7:439–73 transport properties and non-Newtonian flu-
Brünger A, Brooks CL, Karplus M. 1984. idy rheology. Ind. Eng. Chem. Res. 31:1237–
Stochastic boundary conditions for molec- 52
ular dynamics simulations of ST2 water. Darden T, York D, Pedersen L. 1993. Parti-
Chem. Phys. Lett. 105(5):495–500 cle mesh Ewald: An N · log N method for
Ceniceros HD, Hou TY. 2001. An effi- Ewald sums in large systems. J. Chem. Phys.
cient dynamically adaptive mesh for poten- 98(12):10089–92
tially singular solutions. J. Comput. Phys. Degond P, Mas-Gallic S. 1989a. The weighted
172(2):609–39 particle method for convection-diffusion
Chaniotis AK, Poulikakos D, Koumoutsakos P. equations. Part 1: The case of an isotropic
2002. Remeshed smooth particle hydrody- viscosity. Math. Comput. 53(188):485–507
namics for the simulation of viscous and heat Degond P, Mas-Gallic S. 1989b. The weighted
conducting flows. J. Comput. Phys. 182:67– particle method for convection-diffusion
90 equations. Part 2: The anisotropic case. Math.
Chorin AJ. 1973. Numerical study of slightly Comput. 53(188):509–25
viscous flow. J. Fluid Mech. 57(Part 4):785– Delgado-Buscalioni R, Coveney PV. 2003.
96 Continuum-particle hybrid coupling for
Cleveland W, Loader C. 1996. Smoothing by mass, momentum, and energy transfers in
local regression: principles and methods. In unsteady flow. Phys. Rev. E 67:046704–1–
Statistical Theory and Computational As- 046704–13
pects of Smoothing, W Haerdle, GE Schimek, Duarte CA, Oden JT. 1996. An h-p adaptive
ed. pp.10–49. Berlin: Springer method using clouds. Comp. Meth. Appl.
Cottet GH. 1990. A particle-grid superposition Mech. Eng. 139(1-4):237–62
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

FLOW SIMULATIONS USING PARTICLES 485

E W, Engquist B. 2003. The heterogeneous uation of potential fields in three dimensions.

multiscale methods. Commun. Math. Sci. Lect. Notes Math. 1360:121–41
1(1):87–132 Guskov I, Khodakovsky A, Schröder P,
Eldredge J, Colonius T, Leonard A. 2001. A vor- Sweldens W. 2002. Hybrid meshes: multi-
tex particle method for compressible flows. resolution using regular and irregular refine-
AIAA Pap. 2641:1–9 ment. In Proc. 18th Ann. Symp. Comput.
Eldredge JD, Leonard A, Colonius T. 2002. Geom., pp. 264–72. Barcelona: ACM
A general deterministic treatment of deriva- Hadjiconstantinou NG. 1999a. Combining
tives in particle methods. J. Comput. Phys. atomistic and continuum simulations of con-
180(2):686–709 tact line motion. Phys. Rev. E 59(2):2475–78
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.

Ellero M, Kröger M, Hess S. 2002. Viscoelastic Hadjiconstantinou NG. 1999. Hybrid atomistic-
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

flows studied by smoothed particle dynam- continuum formulations and the moving
ics. J. Non-Newton. Fluid Mech. 105(1):35– contact-line problem. J. Comput. Phys.
51 154:245–65
Enright D, Fedkiw R, Ferziger J, Ian M. 2002. Hadjiconstantinou NG, Garcia AL, Bazant MZ,
A hybrid particle level set method for im- He G. 2003. Statistical error in particle simlu-
proved interface capturing. J. Comput. Phys. ations of hydrodynamic phenomena. J. Com-
183:83–116 put. Phys. 187:274–97
Espanol P, Revenga M. 2003. Smoothed dis- Hadjiconstantinou NG, Patera AT. 1997. Het-
sipative particle dynamics. Phys. Rev. E erogeneous atomistic-continuum representa-
67(4):026705–1–026705–12 tions for dense fluid systems. Int. J. Mod. Phy.
Ewald PP. 1921. Die Berechnung Optischer C 8(4):967–76
und Elektrostatische Gitterpotentiale. An- Harlow FH. 1964. Particle-in-cell computing
nalen der Physik 64:253–87 method for fluid dynamics. Methods Com-
Fermi E, Pasta J, Ulam S. 1955. Studies in non- put. Phys. 3:319–43
linear problems. Los Alamos Rep. LA-1940 Hieber SE, Koumoutsakos P. 2004. A Lan-
Fischer P. 1997. An overlapping Schwarz grangian particle level set method. J. Com-
method for spectral element solution of the put. Phys. Submitted
incompressible Navier-Stokes equations. J. Hockney RW, Eastwood JW. 1988. Computer
Comput. Phys. 133:84–101 Simulation Using Particles. 2nd ed. Bristol,
Flekkøy EG, Wagner G, Feder J. 2000. Hybrid PA: Inst. Phys. Publ.
model for combined particle and continuum Hockney RW, Goel SP, Eastwood JW. 1973.
dynamics. Europhys. Lett. 52(3):271–76 A 10000 particle molecular dynamics model
Garcia AL, Bell JB, Crutchfield WY, Alder BJ. with long-range forces. Chem. Phys. Lett.
1999. Adaptive mesh and algorithm refine- 21:589–91
ment using direct simulation Monte Carlo. J. Hoogerbrugge PJ, Koelman JMVA. 1992. Sim-
Comput. Phys. 154:134–55 ulating microscopic hydrodynamics phe-
Ghoniem AF, Oppenheim AK. 1984. Numeri- nomena with dissipative particle dynamics.
cal solution for the problem of flame propa- Europhys. Lett. 19(3):155–60
gation by the random element method. AIAA Hou TY. 1990. Convergence of a variable blob
J. 22(10):1429–35 vortex method for the Euler and Navier-
Gingold RA, Monaghan JJ. 1983. Shock simu- Stokes equations. SIAM J. Numer. Anal.
lation by the particle method sph. J. Comput. 27:1387–404
Phys. 52(2):374–89 Hou TY. 2004. Multiscale modeling and com-
Glotzer SC, Paul W. 2002. Molecular and putation of incompressible flow. In Applied
mesoscale simulation methods for polymer Mathematics Entering the 21st Century: In-
materials. Annu. Rev. Mater. Res. 32:401–36 vited Talks from the ICIAM 2003 Congress,
Greengard L, Rokhlin V. 1988. The rapid eval- ed. JM Hill, R Moore. SIAM Publ.
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

486 KOUMOUTSAKOS

Hummer G, Rasaiah JC, Noworyta JP. 2001. Mathiowetz AM, Jain A, Karasawa N, God-
Water conduction through the hydropho- dard WA III. 1994. Protein simulations using
bic channel of a carbon nanotube. Nature techniques for very large systems—the cell
414:188–90 multipole method for nonbond interactions
Karniadakis GE, Beskok A. 2002. Micro Flows. and the Newton-Euler inverse mass opera-
Fundamentals and Simulation. New York: tor method for internal coordinate dynamics.
Springer-Verlag Proteins 20(3):227–47
Knio OM, Ghoniem AF. 1992. The three- Micci MM, Kaltz TL, Long LN. 2001. Molec-
dimensional structure of periodic vorticity ular dynamics simulations of atomization
layers under non-symmetrical conditions. J. and spray phenomena. Atomization Sprays/
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.

Fluid Mech. 243:353–92 11(4):351–63

Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

Koplik J, Banavar JR. 1995. Corner flow Mittal R, Iaccarino G. 2005. Immersed
in the sliding plate problem. Phys. Fluids boundary methods for viscous flow. Annu.
7(12):3118–25 Rev. Fluid Mech. 37. In press
Koumoutsakos P, Leonard A. 1995. High- Monaghan JJ. 1985a. Extrapolating B splines
resolution simulation of the flow around for interpolation. J. Comput. Phys. 60(2):
an impulsively started cylinder using vortex 253–62
methods. J. Fluid Mech. 296:1–38 Monaghan JJ. 1985b. Particle methods for hy-
Koumoutsakos P, Leonard A, Pépin F. 1994. drodynamics. Comp. Phys. Rep./ 3:71–123
Boundary conditions for viscous vortex Monaghan JJ. 1988. An introduction to SPH.
methods. J. Comput. Phys. 113(1):52–61 Comp. Phys. Commun. 48(1):89–96
Koumoutsakos P, Zimmerli U, Werder T, O’Connell ST, Thompson PA. 1995. Molecular
Walther JH. 2004. Nanoscale fluid mechan- dynamics-continuum hybrid computations: a
ics. In The Handbook of Nanotechnology. tool for studying complex fluid flow. Phys.
Nanometer Structures. Theory, Modeling, Rev. E 52(6):R5792–95
and Simulation. pp. 319–93 Osher S, Fedkiw RP. 2001. Level set methods:
Krasny R. 1986. A study of singularity forma- an overview and some recent results. J. Com-
tion in a vortex sheet by the point vortex ap- put. Phys. 169(2):463–502
proximation. [Link] Mech. 167:65–93 Ould-Salihi ML, Cottet GH, El Hamraoui M.
Kremer K, Müller-Plathe F. 2002. Multiscale 2000. Blending finite-difference and vortex
simulation in polymer science. Mol. Sim. 28: methods for incompressible flow computa-
729–50 tions. SIAM J. Sci. Comput. 22(5):1655–74
Leonard A. 1985. Computing three-dimen- Ploumhans P, Winckelmans GS, Salmon JK,
sional incompressible flows with vortex el- Leonard A, Warren MS. 2002. Simulation of
ements. Annu. Rev. Fluid Mech. 17:523–59 three-dimensional bluff body flows: applica-
Li J, Liao D, Yip S. 1998. Coupling continuum tions to the sphere at Re = 300, 500 and 1000.
to molecular-dynamics simulation: reflecting J. Comput. Phys. 178:427–63
particle method and the field estimator. Phys. Poncet P. 2004. Topological aspects of the
Rev. E 57(6):7259–67 three-dimensional wake behind rotary oscil-
Li J, Liao D, Yip S. 1999. Nearly exact solu- lating circular cylinder. J. Fluid Mech. In
tion for coupled continuum/MD fluid simu- press
lation. J. Comput.-Aided Mat. Design 6:95– Rahman A, Stillinger FH. 1971. Molecular dy-
102 namics study of liquid water. J. Chem. Phys.
Maruyama S. 2001. Molecular dynamics meth- 55(7):3336–59
ods for microscale heat transfer. In Advances Raviart PA. 1986. Particle approximation
in Numerical Heat Transfer, Vol. 2, ed. WJ of first-order systems. J. Comput. Math.
Minkowycz, EME Sparrow, pp. 189–226. 4(1):50–61
New York: Taylor and Francis Rosenhead L. 1930. The spread of vorticity in
 23 Nov 2004 1:49 AR [Link] [Link] LaTeX2e(2002/01/18) P1: IBD

FLOW SIMULATIONS USING PARTICLES 487

the wake behind a cylinder. Proc. R. Soc. A S, Klein ML. 1997. Modified nonequilib-
127(A):590–12 rium molecular dynamics for fluid flows
Ryckaert JP, Bellemans A, Ciccotti G, Paolini with energy conservation. J. Chem. Phys.
GV. 1989. Evaluation of transport coeffi- 106(13):5615–21
cients of simple fluids by molecular dynam- van Kampen NG. 1981. Stochastic Processes in
ics: comparison of Green-Kubo and nonequi- Physics and Chemistry . Amsterdam: North-
librium approaches for shear viscosity. Phys. Holland
Rev. A 39(1):259–67 Verlet L. 1968. Computer experiments on clas-
Sagui C, Darden TD. 1999. Molecular dy- sical fluids. ii. Equilibrium correlation func-
namics simulations of biomolecules: long- tions. Phys. Rev. 165(1):201–14
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.

range electrostatic effects. Annu. Rev. Bio- Wagner G, Flekkøy E, Feder J, Jossang T.
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

phys. Biomol. Struc. 28:155–79 2002. Coupling molecular dynamics and

Sarpkaya T. 1989. Computational methods with continuum dynamics. Comp. Phys. Commun.
vortices—the 1988 Freeman scholar lecture. 147:670–73
J. Fluids Eng. 111:5–52 Walther JH, Jaffe RL, Kotsalis EM, Werder T,
Schlick T. 2002. Molecular modeling and sim- Halicioglu T, Koumoutsakos P. 2004. Hy-
ulation, Vol. 21. In Interdisciplinary Applied drophobic hydration of C60 and carbon nan-
Mathematics. New York: Springer-Verlag otubes in water. Carbon 42(5–6)1185–94
Schoenberg IJ. 1946. Contribution to the prob- Walther JH, Koumoutsakos P. 2001. Three-
lem of approximation of equidistant data by dimensional particle methods for particle
analytic functions. Q. Appl. Math. 4:45–99, laden flows with two-way coupling. J. Com-
112–41 put. Phys. 167:39–71
Serrano M, Espanol P. 2001. Thermodynam- Walther JH, Werder T, Jaffe RL, Koumoutsakos
ically consistent mesoscopic fluid particle P. 2004. Hydrodynamic properties of car-
model. Phys. Rev. E 64:046115–1–046115– bon nanotubes. Phys. Rev. E 69:062201–1–
18 062201–4
Sokhan VP, Nicholson D, Quirke N. 2002. Fluid Warren PB. 1998. Dissipative particle dynam-
flow in nanopores: accurate boundary condi- ics. Curr. Opin. Colloid Inter. Sci. 3:620–
tions for carbon nanotubes. J. Chem. Phys. 29
117(18):8531–39 Werder T, Walther JH, Koumoutsakos P. 2004.
Strain J. 1997. Fast adaptive 2D vortex meth- Hybrid atomistic-continuum method for the
ods. J. Comput. Phys. 132:108–22 simulation of dense fluid flow. J. Comput.
Sun QH, Boyd ID, Candler GV. 2004. A hy- Phys. Accepted
brid continuum/particle approach for model- Zimmermann S, Koumoutsakos P, Kinzelbach
ing subsonic, rarefield gas flow. J. Comput. W. 2001. Simulation of pollutant transport
Phys. 194(1):256–77 using a particle method. J. Comput. Phys.
Tuckerman ME, Mundy CJ, Balasubramanian 173:322–47
 P1: JRX
November 24, 2004 12:5 Annual Reviews AR235-FM

Annual Review of Fluid Mechanics

Volume 37, 2005

CONTENTS
ROBERT T. JONES, ONE OF A KIND, Walter G. Vincenti 1
by Institute of Mechanics - Chinese Academy of Sciences on 01/17/07. For personal use only.
Annu. Rev. Fluid Mech. 2005.37:457-487. Downloaded from [Link]

GEORGE GABRIEL STOKES ON WATER WAVE THEORY,

Alex D.D. Craik 23
MICROCIRCULATION AND HEMORHEOLOGY, Aleksander S. Popel
and Paul C. Johnson 43
BLADEROW INTERACTIONS, TRANSITION, AND HIGH-LIFT AEROFOILS
IN LOW-PRESSURE TURBINES, Howard P. Hodson
and Robert J. Howell 71
THE PHYSICS OF TROPICAL CYCLONE MOTION, Johnny C.L. Chan 99
FLUID MECHANICS AND RHEOLOGY OF DENSE SUSPENSIONS, Jonathan
J. Stickel and Robert L. Powell 129
FEEDBACK CONTROL OF COMBUSTION OSCILLATIONS, Ann P. Dowling
and Aimee S. Morgans 151
DISSECTING INSECT FLIGHT, Z. Jane Wang 183
MODELING FLUID FLOW IN OIL RESERVOIRS, Margot G. Gerritsen
and Louis J. Durlofsky 211
IMMERSED BOUNDARY METHODS, Rajat Mittal
and Gianluca Iaccarino 239
STRATOSPHERIC DYNAMICS, Peter Haynes 263
THE DYNAMICAL SYSTEMS APPROACH TO LAGRANGIAN TRANSPORT
IN OCEANIC FLOWS, Stephen Wiggins 295
TURBULENT MIXING, Paul E. Dimotakis 329
GLOBAL INSTABILITIES IN SPATIALLY DEVELOPING FLOWS:
NON-NORMALITY AND NONLINEARITY, Jean-Marc Chomaz 357
GRAVITY-DRIVEN BUBBLY FLOWS, Robert F. Mudde 393
PRINCIPLES OF MICROFLUIDIC ACTUATION BY MODULATION OF
SURFACE STRESSES, Anton A. Darhuber and Sandra M. Troian 425
MULTISCALE FLOW SIMULATIONS USING PARTICLES,
Petros Koumoutsakos 457

vii

View publication stats