PHYSICAL REVIEW E 71, 036707 共2005兲

Simulation method to resolve hydrodynamic interactions in colloidal dispersions
1

Yasuya Nakayama1,* and Ryoichi Yamamoto1,2
PRESTO, Japan Science and Technology Agency, 4-1-8 Honcho Kawaguchi, Saitama 332-0012, Japan2
Department of Chemical Engineering, Kyoto University, Kyoto 615-8510, Japan
共Received 2 March 2004; published 25 March 2005兲

A computational method is presented to resolve hydrodynamic interactions acting on solid particles immersed in incompressible host fluids. In this method, boundaries between solid particles and host fluids are
replaced with a continuous interface by assuming a smoothed profile. This enabled us to calculate hydrodynamic interactions both efficiently and accurately, without neglecting many-body interactions. The validity ofthe method was tested by calculating the drag force acting on a single cylindrical rod moving in an incompressible Newtonian fluid. This method was then applied in order to simulate sedimentation process of colloidal dispersions.
DOI: 10.1103/PhysRevE.71.036707

PACS number共s兲: 83.85.Pt, 82.70.⫺y, 83.10.Pp, 82.20.Wt

I. INTRODUCTION

There are a number of useful systems consisting of small
solid particles dispersed in host fluids. Among them, colloidal dispersions are most common to our daily life and are ofgreat importance, particularly in the fields of engineering andbiology 关1,2兴. Colloidal dispersions have been reported to
exhibit several unusual phenomena, such as long-range correlations in sedimenting particles 关3兴, long-range anisotropicinteractions in liquid-crystal colloidal dispersions 关4兴, transient gel formations during phase separations of colloidal
suspensions 关5兴, and electrorheological effects in particle
suspensions of nonconductive fluids 关6兴.
Since the dynamics of colloidal dispersions is very complicated, it is extremely difficult to investigate their dynamic
properties by means of analytical methods alone. Computational approaches are necessary in order to elucidate the truemechanisms of dynamic phenomena in a variety of situations. Colloidal dispersions, however, have a typical multiscale problem. The molecules comprising host fluids aremuch smaller and move much faster than colloidal particles.
From a computational point of view, performing fully microscopic molecular simulations for this kind of multiscale system is extremely inefficient. An alternative, which is generally considered much better than microscopic simulations, isto treat host fluids as coarse-grained continuum media.
Several numerical methods have been developed in aneffort to simulate colloidal dispersions. Two of the mostwell-known methods are the Stokesian dynamics 关7兴 and theEulerian-Lagrangian method. The former is thought to be themost efficient method, capable of treating hydrodynamic interactions properly. Furthermore, it can be implemented as an
O共N p兲 scheme for N p particles by utilizing the fast multipole
method 关8兴. However, it is extremely difficult to deal with
dense dispersions and dispersions consisting of nonsphericalparticles by means of Stokesian dynamics due to the complicated mathematical structures used in Stokesian dynamics.
On the other hand, the Eulerian-Lagrangian method is a very

*Electronic address: nakayama@cheme.kyoto-u.ac.jp1539-3755/2005/71共3兲/036707共7兲/$23.00

natural and sensible approach to simulate solid particles with
arbitrary shapes. A number of kinds of tailor-made mesh,including unstructured mesh, overset mesh, and boundaryfitted coordinates, have been applied to specific problems, sothat the shapes of the particles are properly expressed in thediscrete mesh space. Thus, in principle, it is possible to apply
this method to dispersions consisting of many particles withany shape. However, a numerical inefficiency arises from thefollowing: 共i兲 reconstructions of the irregular mesh are necessary at every simulation step according to the temporal
particle position, and 共ii兲 the Navier-Stokes equation must be
solved with boundary conditions imposed on the surfaces of
all colloidal particles. The computational demands thus areenormous for systems involving many particles, even if the
shapes are all spherical.
Thus, our goal is to develop an efficient simulation
method that can be applied to particle dispersions in complex
fluids. Since host fluids are considered incompressible insuch systems, an efficient simulation must address how toefficiently and accurately evaluate hydrodynamic interactions. As a first step towards this goal, we attempted to develop a method to simulate colloidal dispersions in simple
Newtonian fluids. The reliability of this method was tested
by calculating the drag force acting on a cylindrical object ina flow. Its performance was subsequently demonstrated bysimulating the sedimentation processes of colloidal particles
in a Newtonian fluid within a small Reynolds number regime.
II. SIMULATION METHOD

In order to overcome the problems arising at the solidfluid interface in the Eulerian-Lagrangian method, rather
than the original discontinuous rectangle profile 共interfacial
thickness, ␰ = 0兲 schematically depicted in Fig. 1, a smoothed
profile was introduced to the interface 共␰ ⬎ 0兲. This simple
modification greatly benefits the performance of numericalcomputations, compared to the original Eulerian-Lagrangian
method for the following reasons.
共i兲 Regular Cartesian coordinates can be used for many
particle systems with any particle shape, rather than

036707-1©2005 The American Physical Society

PHYSICAL REVIEW E 71, 036707 共2005兲

Y. NAKAYAMA AND R. YAMAMOTO

FIG. 1. An example of the smoothed profile 共dashed line兲. The
original rectangular profile is also shown for comparison 共solidline兲.

boundary-fitted coordinates. The solid-fluid interface has afinite volume 共⬀␲ad−1␰, with a and d as the particle radiusand system dimension兲 supported by multiple grid points.
Thus, the round particle shape can be treated in fixed Cartesian coordinates without difficulty. The simulation scheme isthus free from the mesh reconstruction problem that significantly suppresses the computational efficiency of theEulerian-Lagrangian method. In addition, the simple Cartesian coordinate enables use of periodic boundary conditions
as well as the fast Fourier transformation 共FFT兲.
共ii兲 At the interfaces, the velocity component in the direction normal to the interface of the host fluid must be equal tothat of the particle. This kinetic condition is imposed in the
Navier-Stokes equation, as the boundary value condition defined for the solid-fluid interface in typical methods. In our
method, however, this condition is automatically satisfied byan incompressibility condition on the entire domain, which
will be subsequently explained in detail.
共iii兲 The computational demands for this method includesensitivity to the number of grid points 共volume of the total
system兲; however, it is insensitive to the number of particles.
Thus, our method is thought to be suitable for simulating
dense colloidal dispersions.
The nonzero interfacial thickness ␰ is the only approximation used in the present method. Thus, interparticle hydrodynamic interactions can be fully resolved within the approximation of the nonzero thickness in the present method.
There have been two similar methods developed by different authors 关9,10兴. The basic ideas of these methods is touse a fixed grid and represent the particles not as boundary
conditions to the fluid, but by a body force or Lagrange
multipliers in the Navier-Stokes equation. The essential difference in our approach to these two methods is to introduce
an explicit diffuse interface in a smoothed profile. As long asa fixed grid is used, the moving boundary is inevitably represented as diffuse as grid spacing. The introduction of a
explicit smoothed profile allows us to present a clear formulation of a numerical algorithm and is advantagous when it isapplied to complex fluids 关11,12兴
The lattice Boltzmann 共LB兲 method 关13兴 has attracted
much attention in recent years to simulate colloidal disper-

sions with hydrodynamic interactions 关14兴. The LB equation
was proved to offer a faithful discretization of the NavierStokes equation, and colloidal dispersions are simulated in
the Eulerian-Lagrangian manner. In practical viewpoints, the
LB method is forumlated on a fixed Cartesian lattice and iswell adapted to parallel computation. The LB approach, although the formulation is not intuitive and its treatment of
the moving solid-fluid boundary is somewhat complicated,
has several similar merits to the present method.
The “fluid particle dynamics” 共FPD兲 method was proposed earlier and is similar to our method in spirit 关15兴. In
this method, although a similar smoothed profile wasadopted, there are several important differences between
FPD and the present method. The most significant differenceis that particles are modeled as a highly viscous fluid withviscosity ␩c, much greater than the fluid viscosity ␩s in FPD.
This enables the rigidity of the particles to be sustained approximately by artificial diffusivity ⌬␩␾共x , t兲 共⌬␩ ⬅ ␩c − ␩sⰇ ␩s兲 within the particle domain. While this model is physically correct, a practical problem remains in that a larger
viscosity requires smaller time increments. In contrast, thepresent method treats colloidal particles as undeformable solids 共i.e., ⌬␩ → ⬁兲, and thus no additional constraint arises in
the numerical implementations.
A. Basic working equations

Colloidal dispersions are considered in a simple Newtonian liquid. The motion of the host fluid is governed by the
Navier-Stokes equation with the incompressibility condition
共⳵t + u f · ⵱ 兲u f =1
⵱ · ␴f ,␳⵱ · u f = 0,

共1兲
共2兲

where u f is the fluid velocity and ␳ is the fluid density. The
stress tensor is represented by␴ f = − pI + ␩兵⵱u f + 共⵱u f 兲T其,共3兲

where p is the pressure and ␩ is the fluid viscosity.
The colloidal particles are assumed to be rigid and spherical, and their positions Ri are tracked in a Lagrangian reference frame,Ṙi = Vi ,共4兲with the translational momentum equation
PPg
M iV̇i = FH
i + Fi + Fi

共5兲

and the angler momentum equation
Ii · ⍀̇i = NH
i ,

共6兲

where Ri, Vi, ⍀i, M i, and Ii are the position, the translational
velocity, the angular velocity, the mass, and the inertia tensor
of the ith particle, respectively. The hydrodynamic force FHiand torque NH
i acting on a particle can be obtained by integrating the stress tensor over the surface as

036707-2

PHYSICAL REVIEW E 71, 036707 共2005兲

SIMULATION METHOD TO RESOLVE HYDRODYNAMIC…

FH
i =

冕

␴ f · dS,

⵱ · up = ⵱ ·

共7兲

Si

NH
i =

冕

=
r ⫻ 共␴ f · dS兲,

共8兲

Si

where r is the relative position vector from the center ofrotation to the colloid surface. Furthermore, FiPP is the force
due to direct particle-particle interactions, and Fgi = M i共1− ␳−1* 兲g is the buoyant force where ␳* = ␳c / ␳ is the mass density ratio of the particles to the host fluid and g is the gravitational acceleration. Relevant dimensionless parameters in
the above equations include the Reynolds number Re= UL / ␯, the Froude number Fr= U / 冑gL, the mass density
ratio ␳*, and the volume fraction ␣. Here U and L represent
the typical velocity and length scales specific to the systems
under consideration, respectively. The kinematic viscosity ␯= ␩ / ␳ and the mass density of colloidal particles, ␳c, areassumed to be constant.
In typical methods, the above set of equations should be
solved using proper boundary conditions defined at the solidfluid interface. In the present method, however, the solidfluid boundary condition is replaced with a body force andan incompressibility condition on a total velocity defined on
the entire domain.

⵱N␾共x,t兲 ⬅ 兺 ␾i共x,t兲,

共9兲

⵱ · u = 共⵱ ␾兲 · 共u p − u f 兲.

共10兲
The total 共fluid+ particle兲 velocity field is then given by
共11兲

Since the particle velocity field u p is constructed from therigid motions of particles, ⵱ · u p = 0 is verified as

共12兲

共14兲

The gradient of ␾ is proportional to the surface-normal vector and has a support on the interfacial domains. Therefore,
the incompressibility condition on the total velocity ⵱ · u = 0
means the solid-fluid impermeability condition at the solidfluid interface.
We are to derive the evolution of the total velocity u. To
make the points clearer, we first consider the problem assuming that the motions of particles 兵Ri共t兲 , Vi共t兲 , ⍀i共t兲其 are
given. In Eq. 共11兲, only the fluid velocity u f is to be solved.
The evolution equation of the total velocity is split as
共⳵t + u · ⵱ 兲u =1
⵱ · ␴,
␳⳵ tu = ␾ f p ,

共15兲
共16兲

where the stress tensor is␴ = − p*I + ␩兵⵱u + 共⵱u兲T其.

共17兲

By integrating Eqs. 共15兲 with 共17兲, the total velocity is predicted as u = u*. The pressure p* in Eq. 共17兲 is determined tofulfill the incompressibility condition ⵱ · u* = 0. Then, thebody force ␾f p is added to enforce Eq. 共11兲 and the solidfluid inpermeability condition. Therefore, the time-integrated
body force ␾f p is determined as冕t+hds␾f p = ␾共u p − u*兲 −t1
⵱ p ph,␳

共18兲

where the pressure p p is determined to fulfill ⵱ · u= 共⵱␾兲 · 共u p − u*兲 = 0. By solving Eq. 共16兲 with the body force,
Eq. 共18兲, we finally arrive at Eq. 共11兲 where the fluid velocityis共1 − ␾兲u f = 共1 − ␾兲u* −i=1

u共x,t兲 ⬅ 共1 − ␾兲u f + ␾u p = u f + ␾共u p − u f 兲.

␾兺i 兵Vi + ⍀i ⫻ 共x − Ri兲其 · ⵱ ␾i ,Assuming the incompressibility of the fluid velocity u f , thedivergence of the total velocity isN␾u p共x,t兲 ⬅ 兺 兵Vi共t兲 + ⍀i共t兲 ⫻ 关x − Ri共t兲兴其␾i共x,t兲.

冎
共13兲

i=1

using the ith particle’s profile function ␾i共x兲 centered at Ri.
Several possible mathematical forms for ␾i共x兲 exist, however; some typical functions are listed in the Appendix.
The continuum velocity field u p is defined for the solid
particles using 兵V1 , . . . , VN其, 兵⍀1 , . . . , ⍀N其, and ␾i as␾

␾i 共⵱ ␾i兲␾ − ␾i共⵱ ␾兲 共⵱ ␾i兲␾i − ␾i共⵱ ␾i兲
=
=
= 0.
␾
␾2
␾2

B. Modified working equations

Assuming a smoothed profile with a finite thickness ␰ tothe solid-fluid interface, we here derive the body force which
accurately takes the interactions between solids and fluids
due to the motions of colloids in an incompressible fluid into
consideration. The present study considers a monodisperse
system consisting of N spherical particles with radius a. The
positions of the particles 兵R1 , . . . , RN其 are first transformed toa continuous field,

再兺i Vi + ⍀i ⫻ 共x − Ri兲 ␾i1
⵱ p ph.
␳

共19兲

We note that the nonslip condition at the solid-fluid interface
is fulfilled in this time evolution of the total velocity. Since
the viscous stress 共17兲 acts on the entire domain including
the interfacial domain, the tangential velocity difference between u f and u p is reduced. In other words, the nonslip or
slip condition can be imposed in the definition of the stress ␴used in Eq. 共15兲. When compared to FPD 关15兴, their bodyforce is ␾f p = ⌬␩ ⵱ · ␾兵⵱u + 共⵱u兲T其 with artificial diffusivity⌬␩. In the limit of ⌬␩ → ⬁, the particle becomes rigid; how-

036707-3

PHYSICAL REVIEW E 71, 036707 共2005兲

Y. NAKAYAMA AND R. YAMAMOTO

ever, this limit cannot be achieved by the numerical scheme
used in FPD. In contrast to FPD, the body force ␾f p guarantees the rigidity of the solid particle without additional large
artificial diffusivity.
We complete the time evolution by deriving the hydrodynamic force FH
i acting on the particles. The hydrodynamic
force is defined as the momentum flux between the fluid andsolid. Thus, the hydrodynamic force is simply the counteraction from the fluid:
FH
i =−

冕

␳␾if pdx.

共20兲

C. Simulation procedure

共i兲 For a given particle configuration 兵Rni 其, velocity 兵Vni 其,and angular momentum 兵⍀ni 其 共i = 1 , . . . , N兲, where the superscript n denotes the time step and h is the time increment, thefluid velocity at a time t = nh is predicted as冕

tn−1+hds ⵱ ·

tn−1

冉

冊

1
␴ − uu ,␳

共21兲

under the incompressibility condition ⵱ · u* = 0 which determines the intermediate pressure p*.
共ii兲 At step n, the total velocity un should be equal to u p
within the particle domain and the surface-normal velocity
components of particles and fluid should be match in theinterfacial domain. Thus, it must be corrected by the bodyforce f p defined by␾f np = 兵␾n共unp − u*兲其/h −1
⵱ pp .
␳

共22兲

The correcting pressure p p is determined to make a resultant
total velocity incompressible. This leads to the Poisson equation of p p:
ⵜ2 p p = ␳ ⵱ · 兵␾共unp − u*兲其/h.

共23兲

Finally, we have
un = u* + ␾f nph,

共24兲

pn = p* + p p .

共25兲

共iii兲 The hydrodynamic force and torque acting on eachcolloidal particle are now computed using the volume integrals,

NH
i =−

冕

冕

␳␾if npdx,

共26兲

␳共x − Rni 兲 ⫻ ␾if npdx.

共27兲

FH
i =−

1
Mi

冕

tn+h

PP
g
ds共FH
i + Fi + Fi 兲,

共28兲

tn

= ⍀ni + I−1
⍀n+1
i
i ·

冕

tn+h

dsNH
i ,

共29兲

tn

= Rni +
Rn+1
i

冕

tn+h

dsVi .

共30兲

tn

In contrast to the force in Eq. 共7兲 which is expressed as thesurface integral, the above force is given as the volume integral. The volume integral is much advantageous in meshbased discretization compared to the surface integral since nogeneration of body-fitted mesh is needed.

u* = un−1 +Vn+1= Vni +iand the velocity, angular velocity, and position of each colloidal particles at the current step n are updated to 共n + 1兲 as

Since the same ␾f p is used for both the host fluid 共24兲 andthe colloidal particles 共26兲 and 共27兲 through the interface, no
excess or shorts for solid-fluid interactions exist. The same
type of solid-fluid interaction has been previously proposed
in Ref. 关16兴, though the treatment used in the present method
is more general. It is important to note that the timing of the
updates for fluid 共21兲 and particles 共28兲–共30兲 has beenshifted; the particles always go one step ahead of the fluid.
This shift is primarily due to a technical issue related to thesolid-fluid boundary condition. In general, the boundary
value condition, which is replaced with the force due tosolid-fluid interactions in the present method, is necessary to
update the fluid. Otherwise, in order to update both the fluid
and particles simultaneously, the implicit scheme for particles must be used. Consider the problem to integrate the
n−1
n−1
n−1
, to the共n − 1兲th step variables, 兵Rn−1
i , Vi , ⍀i 其 with unnth step. The particle velocity of the next step, u p, in Eq. 共22兲is needed to update the total velocity un−1 to un before
n−1
n−1
兵Rn−1
i , Vi , ⍀i 其 is updated. This situation requires the implicit treatment. The use of the implicit scheme complicatesthe algorithm and reduces the efficiency. The timing shift istherefore necessary to realize the full explicit scheme described above. The initial condition 兵Rni , Vni , ⍀ni 其 with un−1
can be generically constructed. From the total velocity un−1
satisfying the initial boundary condition given by
n−1
n−1
兵Rn−1
i , Vi , ⍀i 其, the hydrodynamic force and torque arecomputed as surface integrals and the particle trajectory is
integrated from the 共n − 1兲th to nth step. Together with theconstruction of the initial condition, the timing-shift algorithm is realized without loss of generality.
In the above formulation, although the solid-fluid interaction force 共22兲 was computed on the basis of the Euler
scheme for simplicity of presentation, implementations ofhigher-order schemes are straightforward for Eqs. 共28兲–共30兲.
Furthermore, in order to update the host fluid in Eq. 共21兲, norestrictions exist for the time discretization and any conventional scheme can be used for incompressible fluids.
III. NUMERICAL RESULTS

The present method has been applied to two specific problems: One is the calculation of the drag force acting on an
infinitely long cylindrical rod moving in a Newtonian fluid inorder to check the validity. The method was also applied tosimulations of many sedimenting particles in a twodimensional fluid in order to demonstrate the performance.
In the present simulations, the Navier-Stokes equation was
discretized with a de-aliased Fourier spectral scheme in

036707-4

PHYSICAL REVIEW E 71, 036707 共2005兲

SIMULATION METHOD TO RESOLVE HYDRODYNAMIC…

FIG. 3. The relative error in the drag coefficient CD as a function of the interfacial thickness ␰ / ⌬.
FIG. 2. 共Color online兲 Schematic representation of the crosssection of geometry around the rod. ⌬ is the lattice spacing, a is therod radius, and ␰ is the interfacial thickness. The rod surface now
has a finite volume ⬃2␲a␰ supported by several grid points on the
fixed Cartesian coordinate.

space and a second-order Runge-Kutta scheme 共the Heun
scheme兲 in time. For the colloidal particles, the velocity andangular velocity were integrated with the Heun scheme andthe position was integrated with the Crank-Nicolson scheme.
The external boundary condition on the edge of the systems
was imposed in the same manner as the fluid-solid boundary
condition on the particle surface. The simulation code is remarkably simple due to such unified treatment for all boundary conditions.

The velocity at the external boundary r = L was set
tou共r = L , ␪兲 = U / 关1 − 2 ln共a / L兲兴关兵1 − 共a / L兲2 − 2 ln共a / L兲其ex− 2兵1 − 共a / L兲2其cos ␪er兴 where ex and er are the unit vectors inthe x and r directions, respectively, and tan ␪ = y / x. An analytical solution for the Stokes equation is known for thisboundary condition, and the drag force is given by FD= 8␲␩U / 关1 − 2 ln共a / L兲兴. The computed CD using the present
method as a function of Re is shown in Fig. 4 and is in good
agreement with the theoretical Stokes law in Re艋 1 within
5%.
The accuracy of the present method using the finite interfacial thickness ␰ / ⌬ = 1 was determined to be acceptable forsimulating colloidal dispersions for NRe 艋 1 based on the numerical results above.
B. Sedimentation

A. Drag force on a cylindrical rod

The drag force acting on an infinitely long cylindrical rod
with radius a was computed by solving the Navier-Stokes
equation around the rod in order to check the accuracy of thepresent method. Figure 2 shows a cross section of the geometry around the rod with finite thickness ␰ at the interface.
First, the effects of the finite thickness on the drag force
are examined in the square domain of L2. An uniform stream
U in the x direction was assigned to the edge of the domain
as the boundary condition. Here the Reynolds number was
defined by Re= 2aU / ␯. The drag coefficient was calculated
as CD = FD / ␳U2a, where the drag force FD was computedfrom Eq. 共26兲 for various values of ⌬, a, L, and U. Figure 3
shows the relative error 兵CD共Re, ␰ / ⌬兲 − CD共Re, ␰ / ⌬= 0兲其 / CD共Re, ␰ / ⌬ = 0兲 as a function of the interfacial thickness ␰ / ⌬, where CD共Re, ␰ / ⌬ = 0兲 was estimated by extrapolating the measured curve of CD共Re, ␰ / ⌬兲 to ␰ / ⌬ → 0. The
relative error in CD was observed to increase with increasing␰ / ⌬; however, it tended to converge within 5% for several
values of a / ⌬ for ␰ / ⌬ = 1 and 0 ⬍ Re⬍ 20. Thus, ␰ / ⌬ = 1 was
set for further simulations.
Next, the drag coefficient CD was calculated. The rod was
fixed at the origin in the circular domain with radius L.

The performance of the present method was examined bysimulating sedimentation processes of monodisperse particles in a two-dimensional Newtonian fluid in a rectangular

FIG. 4. Comparisons of the drag coefficient CD 共plus兲 from our
method with the theoretical curve of the Stokes law 共solid line兲.

036707-5

PHYSICAL REVIEW E 71, 036707 共2005兲

Y. NAKAYAMA AND R. YAMAMOTO

FIG. 5. The normalized lubrication force acts on two approaching infinitely long cylindrical rods as a function of the nearest distance between the two surfaces. Different symbols denote different
approaching velocities U, ranging from 2.5⫻ 10−6 to 9.6⫻ 10−2 by
U = 2.5⫻ 10−6 ⫻ 2.61n, n = 0 , . . . , 11, which almost collapsed. The
observed scaling behavior F = ␩Uf共h兲 with scaling function f共¯兲 ischaracteristic of Stokes flow, due to Reynolds numbers 2aU / ␯ ⬍ 1.

box surrounded by nonslip walls with ␳* = 1.1 and ␣ = 0.143.
The dimensionless parameters were taken to be Re= 0.0916
and Fr= 0.0512, where the settling velocity and the diameterof particle were taken as the characteristic velocity andlength. Other computational parameters were chosen as ⌬= 1, ␰ / ⌬ = 1, a / ⌬ = 10, Lx / ⌬ = 512, and Ly / ⌬ = 1024, where
the y axis is in the direction of gravity. In order to prevent theparticles from overlapping within the core radius ⯝a, theforce was added, FiPP = −⳵E PP / ⳵Ri, due to direct particleparticle interaction using the repulsive part of the LennardN−1 N
兺 j=i+1关共2a / Rij兲12
Jonespotential
E PP = 0.4兺i=1
6
7/6− 共2a / Rij兲 兴␪共2 a − Rij兲, where ␪共¯兲 is the step function
and Rij ⬅ 兩Ri − R j兩. The direct interaction FiPP is not very important when the particles are moving around because theparticles never overlap due to the lubrication effect, even
without FiPP. Figure 5 shows the lubrication force acting on
two approaching rods computed using the present method.
The lubrication force is always repulsive in this case andthus prevents the rods from approaching each other. The
strength of the repulsion increases with increasing velocity
U. On the other hand, when the particles are stacked on thebottom wall during the later stage of sedimentation, FiPP isrequired to sustain the stacking against gravity. In fact, therepulsion vanishes for immobile pairs of rods.
At the initial configuration, all the particles were placednear the upper wall and both fluid and particle velocities
were set to zero, as depicted in Fig. 6共a兲. A typical snapshot
during sedimentation is shown in Fig. 6共b兲. Regions withswirled particles were observed, in which the particle velocities were highly correlated as a result of long-range interparticle hydrodynamic interactions. A simulation with periodic
boundary conditions in the horizontal 共x兲 direction was also
conducted. In this simulation, swirls were still developed;
however, they were smaller than those observed with nonslip
walls. The effect of confinement in the nonslip walls there-

FIG. 6. 共Color兲 Snapshots of 240 colloidal disks sedimenting ina two-dimensional Newtonian fluid obtained using the present
method. The magnitude of the host fluid velocity is indicated incolor; the change of color from blue to red corresponds to a changeof the fluid velocity from small to large.

fore enhances the velocity correlation. The computational demand required for the present simulation is less than one dayof processing on a normal PC.
IV. CONCLUDING REMARKS

A computational method has been developed to simulate
particle dispersion in fluids. Utilizing a smoothed profile forsolid-fluid boundaries, hydrodynamic interactions in many
particle dispersions can be taken fully into account, both accurately and efficiently. In principle, the present method canbe easily applied to systems consisting of many particles
with any shape. The reliability of the method was examinedby calculating the drag force acting on a cylindrical object ina flow. The performance of the method was demonstrated tobe satisfactory by simulating sedimentations of particles in a
Newtonian fluid.
Another primary benefit of using the smoothed profile
arose when the method was extended to colloidal dispersions
in complex fluids with an internal degree of freedom, such as
the molecular orientation or ionic density. In complex fluids,
interparticle interactions can be mediated by the internal degree of freedom of the fluid. In such cases, the fluid-particle
interactions at the colloid surface could be more efficiently
handled by utilizing a smoothed profile. Previous studies onparticle dispersions in liquid-crystal solvents demonstrate astriking example of this efficiency 关11,12兴. Although the hydrodynamic effects were neglected in these simulations, extensions to implement the hydrodynamic effects by incorporating the present method are currently underway.
APPENDIX: SELECTION OF SMOOTHED PROFILES

The specific form of the smoothed profile should be selected according to the convenience of the physical modeling

036707-6

PHYSICAL REVIEW E 71, 036707 共2005兲

SIMULATION METHOD TO RESOLVE HYDRODYNAMIC…

of systems under consideration. In the present study, an infinitely differentiable function with compact support was used.
We adopted ␾ defined as␾i共x兲 = g共兩x − Ri兩兲,g共x兲 =

共A1兲

h关共a + ␰/2兲 − x兴,h关共a + ␰/2兲 − x兴 + h关x − 共a − ␰/2兲兴

h共x兲 =

再exp共− ⌬2/x2兲 x 艌 0,0x ⬍ 0.

冎

共A2兲

共A3兲

␾i共x兲 =

冉冊a − 兩x − Ri兩
1
tanh
+1 .
2
␰

共A4兲

This choice was used in Refs. 关11,12,15兴. This ␾ is also
infinitely differentiable as well as analytically easy to handle.
However, the support is not compact and the separation of
the three domains is ambiguous. Furthermore, the ambiguityof domain separation tends to be more enhanced for higherorder derivatives. For practical implementation, due to exponential decay of the hyperbolic function, a proper cutoff radius is adopted for the calculation of the integrals in Eqs.
共26兲 and 共27兲.
The third possible choice is given by␾i共x兲 = s共a − 兩x − Ri兩兲,

共A5兲

where Ri, a, ␰, and ⌬ were the position of the particle, theradius of the particle, the interfacial thickness, and the latticespacing, respectively. This choice is shown in Fig. 2. While
this ␾ may appear somewhat complicated compared to othermore simple choices, this ␾ has the following benefits: 共i兲three domains—solid, fluid, and interface—are explicitly
separated; namely; ␾ = 1 is the solid domain 共兩x − Ri兩 ⬍ a− ␰ / 2兲, ␾ = 0 is the fluid domain 共a + ␰ / 2 ⬍ 兩x − Ri兩兲, and 0⬍ ␾ ⬍ 1 is the interfacial domain 共a − ␰ / 2 ⬍ 兩x − Ri兩 ⬍ a+ ␰ / 2兲, 共ii兲 high-order derivatives of ␾i with respect to x canbe analytically calculated, and 共iii兲 due to its support compactness, the integrals in Eqs. 共26兲 and 共27兲 remain local,
which contributes greatly to the efficiency of the computation.
The second possible choice iswhich has the property of exact separation of the three domains; however, the second derivative of ␾ is discontinuousat the fluid-interface boundary. Therefore, it is not recommended for computational models requiring derivatives
higher than the second order of ␾.
The detailed choice of ␾ does not affect the results of thepresent simulations because only the first-order derivative of
␾ is required in the present case. However, care must taken if
higher-order derivatives are required.

关1兴 W. B. Russel, D. A. Saville, and W. R. Schowalter, Colloidal
Dispersions 共Cambridge University Press, Cambridge, England, 1989兲.
关2兴 R. G. Larson, The Structure and Rheology of Complex Fluids
共Oxford University Press, New York, 1999兲.
关3兴 P. N. Segré, F. Liu, P. Umbanhowar, and D. A. Weitz, Nature
共London兲 409, 594 共2001兲.
关4兴 H. Stark, Phys. Rep. 351, 387 共2001兲.
关5兴 H. Tanaka, Phys. Rev. E 59, 6842 共1999兲.
关6兴 W. M. Winslow, J. Appl. Phys. 20, 1137 共1949兲.
关7兴 J. F. Brady and G. Bossis, Annu. Rev. Fluid Mech. 20, 111
共1988兲.
关8兴 K. Ichiki, J. Fluid Mech. 452, 231 共2002兲.

关9兴 R. Glowinski, T. W. Pan, T. I. Hesla, D. D. Joseph, and J.
Périaux, J. Comput. Phys. 192, 363 共2001兲.
关10兴 K. Höfler and S. Schwarzer, Phys. Rev. E 61, 7146 共2000兲.
关11兴 R. Yamamoto, Phys. Rev. Lett. 87, 075502 共2001兲.
关12兴 R. Yamamoto, Y. Nakayama, and K. Kim, J. Phys.: Condens.
Matter 16, S1945 共2004兲.
关13兴 S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics
and Beyond 共Clarendon Press, Oxford, 2001兲.
关14兴 A.J. C.Ladd and R. Verberg, J. Stat. Phys. 104, 1191 共2001兲.
关15兴 H. Tanaka and T. Araki, Phys. Rev. Lett. 85, 1338 共2000兲.
关16兴 T. Kajishima, S. Takiguchi, H. Hamasaki, and Y. Miyake,
JSME Int. J., Ser. B 44, 526 共2001兲.

冦冉0,s共x兲 =1␲xsin+ 1 , 兩x兩 ⬍ ␰/2,2
␰

1,036707-7冊

x ⬍ − ␰/2,x ⬎ ␰/2,

冧

共A6兲