Eur. Phys. J. E (2008)
DOI 10.1140/epje/i2007-10332-y

THE EUROPEAN
PHYSICAL JOURNAL E

Simulating (electro)hydrodynamic eﬀects in colloidal dispersions:
Smoothed proﬁle method
Y. Nakayama1,a , K. Kim2 , and R. Yamamoto3,4
1

Department of Chemical Engineering, Kyushu University, Fukuoka 819-0395, Japan
Department of Theoretical and Computational Molecular Science, Institute for Molecular Science, Myodaiji, Okazaki, Aichi
444-8585, Japan3
Department of Chemical Engineering, Kyoto University, Kyoto 615-8510, Japan
4
CREST, Japan Science and Technology Agency - 4-1-8 Honcho Kawaguchi, Saitama 332-0012, Japan
2

Received 18 September 2007 and Received in ﬁnal form 7 May 2008
c EDP Sciences / Società Italiana di Fisica / Springer-Verlag 2008
Published online: 10 June 2008 – 
Abstract. Previously, we have proposed a direct simulation scheme for colloidal dispersions in a Newtonian
solvent (Phys. Rev. E 71, 036707 (2005)). An improved formulation called the “Smoothed Proﬁle (SP)method” is presented here in which simultaneous time-marching is used for the host ﬂuid and colloids.
The SP method is a direct numerical simulation of particulate ﬂows and provides a coupling scheme
between the continuum ﬂuid dynamics and rigid-body dynamics through utilization of a smoothed proﬁlefor the colloidal particles. Moreover, the improved formulation includes an extension to incorporate multicomponent ﬂuids, allowing systems such as charged colloids in electrolyte solutions to be studied. The
dynamics of the colloidal dispersions are solved with the same computational cost as required for solving
non-particulate ﬂows. Numerical results which assess the hydrodynamic interactions of colloidal dispersions
are presented to validate the SP method. The SP method is not restricted to particular constitutive modelsof the host ﬂuids and can hence be applied to colloidal dispersions in complex ﬂuids.
PACS. 47.11.-j Computational methods in ﬂuid dynamics – 82.70.-y Disperse systems; complex ﬂuids –82.20.Wt Computational modeling; simulation

1 Introduction
Interparticle interactions in colloidal dispersions mainly
consist of thermodynamic potential interactions and hydrodynamic interactions [1–3]. Whereas the former interactions occur in both static and dynamic situations,the latter occur solely in dynamic situations. Although
thermodynamic interactions have been studied extensively
and summarized as a concept of the eﬀective interaction [4], the nature of dynamic interactions is poorly understood. Since the hydrodynamic interaction is essentially a long-range, many-body eﬀect, it is extremely diﬃcult to study its role using analytical methods alone. Numerical simulations can aid the investigation of the fundamental role of hydrodynamic interactions in colloidal
dynamics.
In recent decades, various simulations for particulate ﬂow have been developed for most simple situations
in which the host ﬂuid is Newtonian [5–13]. However,these schemes cannot necessarily be applied to problems
with non-Newtonian host ﬂuids or solvents with intera

e-mail: ynakayama@chem-eng.kyushu-u.ac.jp

nal microstructures, which are practically more important
cases. Hydrodynamic simulations of ions and/or charged
colloids have been proposed by making use of some of
the above schemes [14–19]. Nonetheless, the tractability
and/or physical validity of their modeling remain controversial. In reference [14], in solving the motion of the ionicsolutes, the solvent hydrodynamic interaction was incorporated by Rotne-Pragar-Yamakawa–type mobility tensor, which accounted for a long-range part of pair interaction. This scheme misses many-body and/or near-ﬁeld
hydrodynamic interactions between the solutes and surfaces, which can have an eﬀect on the behavior of ionsnear surface. In reference [15], the computational mesh
was arranged to express round shape of a colloidal particle and the boundary value problem of the solvent continuum equations was solved based on the irregular mesh. For
applying this scheme to many-body system, inaccessible
computational resources is inevitable. In references [16–19], they adopted the lattice Boltzmann (LB) method for
explicitly solving the solvent hydrodynamics. Although
the LB method has many advantages for solving large systems, compared to the standard discretization schemes,
the applicability to various constitutive equations for the2

The European Physical Journal E

complex ﬂuids other from the Newtonian ﬂuids is unveiled.
Concerning the interaction between colloidal particles andthe solvent, authors in references [16,17], arranged ﬁnite
interaction points on the surface of the colloids on which africtional coupling was assumed. In this coupling scheme,
friction parameters and the number of the interaction
points were determined phenomenologically and a theoretical basis to determine the parameters in general systems
was not clearly stated.
In this article, we propose a direct simulation scheme
for colloidal dispersions which is applicable to most constitutive models of the host ﬂuids. We call it the “Smoothed
Proﬁle (SP) method” since the original sharp interface
between the colloids and solvent is replaced with an eﬀective smoothed interface with a ﬁnite thickness [7,20]. We
formulate a computational method to couple the particle
dynamics and hydrodynamics of the solvent. A ﬁxed gridis used in both the solvent and the particle domains. Introduction of a smoothed proﬁle makes it possible to realize
stable and eﬃcient implementation of our scheme. The numerical implementation for Newtonian solvents and electrolyte solutions as speciﬁc examples of multi-component
ﬂuids is outlined. Various test cases which verify the SP
method and assess the hydrodynamic interactions are presented.

2 Dynamics of a multi-component solventand colloids

Navier-Stokes equation of incompressible ﬂow with thesource term from solutes:∇ · v = 0,
ρ (∂t +v · ∇) v = −∇p+η∇ v−
2



(4)
Cα ∇μα + ∇ · s, (5)αwhere ρ is the total mass density of the ﬂuid, p is the pressure, η is the shear viscosity of the ﬂuid, and s is a randomstress satisfying the ﬂuctuation-dissipation relation [21]:sik (x, t)sjl (x , t ) =2kB T η(δij δkl + δil δkj )δ(x − x )δ(t − t ).

(6)

The above set of equations is closed when a set of chemical potentials {μα } is given, and describes the dynamics
of a multi-component ﬂuid. For the speciﬁc application ofan electrolyte solution, we consider the Poisson-NernstPlanck equation for the chemical potential
μα ({C1 , . . . , CN }) = kB T log Cα +Zα e(Φ−E ext · x), (7)∇2 Φ = −ρe ,(8)This equation describes the Poisson-Boltzmann distribution for ions at equilibrium, where Zα is the valence of the
α-th ion, e is the elementary charge, Φ is the electrostatic
potential, E ext is the external ﬁeld, is the dielectric constant of the ﬂuid, and ρe is the charge density ﬁeld. Thisset of equations corresponds to the electrokinetic equations which appear in standard textbooks [3].

2.1 Hydrodynamics of multi-component ﬂuids
We ﬁrst give a brief description of multi-component ﬂuid
equations, looking at electrolyte solutions as a speciﬁc application. Consider N (possibly ionic) solute species thatsatisfy the law of conservation for every concentration, Cα ,of the α-th species:∂t Cα + ∇ · Cα v α + ∇ · g α = 0,(1)

where v α is the velocity of the α-th solute and g α is arandom current. Since the inertial time scales of the solute molecules are extremely small, the velocity of the α-thsolute can be decomposed into the velocity of the solventv and the diﬀusive current arising from the chemical potential gradient ∇μα asv α = v − Γα ∇μα ,(2)where Γα kB T is the diﬀusivity of the α-th ion, kB is
the Boltzmann constant and T is the temperature. The
random current should satisfy the following ﬂuctuationdissipation relation [21]:gα,i (x, t)gβ,j (x , t ) =2(kB T )2 Γα δαβ δij δ(x − x )δ(t − t ).

(3)Here, for the sake of simplicity, we do not deal with thecross diﬀusion of diﬀerent solutes. The conservation of momentum implies that the velocity of solvent follows the

2.2 Colloids in electrolyte solutions
The colloid dynamics are maintained by the force exerted
by the solvent. Consider monodisperse spherical colloids
with a radius a, a mass Mp , and a moment of inertia
I p . Momentum conservation between the ﬂuid and thei-th colloid implies the following hydrodynamic force andtorque:HH
F i = (dS i · σ), N i = (x − Ri ) × (dS i · σ), (9)where Ri is the center of mass, dS i (. . . ) indicates thesurface integral over the i-th colloid, and σ is the stresstensor of the ﬂuid. In terms of the electrokinetic equations,the stress reads as
σ = −pI + σ  + σ st + s,

(10)

in which σ  = η(∇v +(∇v)T ) is the dissipative stress, andσ st = {EE − (|E|2 /2)I} is the Maxwell stress for theelectric ﬁeld E = −∇Φ + E ext , where I is the unit tensor.
The evolution of the colloids follows Newton’s equations:Ṙi = V i ,cext
Mp V̇ i = F H
i + Fi + Fi ,ext
I p · Ω̇ i = N H
i + Ni ,

(11)
(12)
(13)

Y. Nakayama et al.: Simulating (electro)hydrodynamic eﬀects in colloidal dispersions: Smoothed proﬁle method

where F ext
and N ext
are the external force and torque,iirespectively, and F ci is the force arising from the core potential of the particles, which prevents the colloids fromoverlapping. Hereinafter, the soft-core potential of thetruncated Lennard-Jones potential is adopted for F ci . For
the charged-colloid system speciﬁcally, the buoyancy is included in the term F ext
i , and the external ﬁeld on thecolloids is accounted for by F Hi .
Having colloidal particles with a ﬁnite volume provides the relevant boundary conditions to the hydrodynamic equations. For the solvent velocity, a no-slip condition is assigned, such that v = V i + Ω i × r i , withr i = x−Ri for the i-th colloid. For the concentration ﬁeld,
a no-penetration condition is assumed, giving n·∇μα = 0,
where n represent the unit normal to the surface of thecolloids. Coupling the hydrodynamics of the solvent withthe dynamics of the colloids deﬁnes the moving-boundarycondition problem above. The usual numerical techniquesof partial diﬀerential equations are hopeless in dealing
with the dynamical evolution of many colloids, since the
sharp interface at the surface of the colloids moves andhenceforth the mesh points at which the boundary condition is assigned vary with each discrete time step. The
moving boundary condition leads to huge computationalcosts.
In contrast, the SP method formulates an eﬃcient
scheme for this kind of moving-boundary-condition problem, and incurs the same level of computational cost as
is required for solving a uniform ﬂuid. The typical computational costs for the ﬂuid and particles are assumed toscale to their degrees of freedom. In the SP method, regular mesh, but not body-ﬁtted mesh, can be used, whichindicates that the inclusion of the dispersed particle phase
does not induce the increase of the grids. In d-dimensional
system, we assume that discretized space contains N dgrids and the variables of the ﬂuid phase are linked to
the grids. Np particles are dispersed in the system andtheir volume scales to Np ad . The number of the particles
are at most Np ad < (N Δx)d , where Δx is the lattice spacing, which results in Np /N d < (Δx/a)d . This inequality
indicates that the computational costs for the dispersed
phase tracking is at largest (Δx/a)d times for the ﬂuids.
In the typical case of a/Δx = 5, the fraction of the computational costs for the particles roughly estimated less thanseveral percent, thus most of the computation should be
for solving the uniform ﬂuid.

3 Computational algorithm
In the SP method, quantities are deﬁned over the entire
domain, which consists of the ﬂuid domain and the particle
domain. To designate the particle domain, we introduce
a concentration ﬁeld for the colloids, given as φ(x, t) =Np
i=1 φi (x, t), where φi ∈ [0, 1] is the proﬁle ﬁeld of thei-th particle, which is unity at the particle domain, zeroat the ﬂuid domain, and which has a continuous diﬀuse
interface of ﬁnite thickness ξ at the interface domain. Withthe ﬁeld φ, the total velocity ﬁeld and concentration ﬁelds

3of the solutes are deﬁned asv = (1 − φ)v f + φv p ,Cα = (1 − φ)Cα∗ ,

(14)
(15)

where (1 − φ)v f represents the velocity ﬁeld of the ﬂuid,
Np
φi (x, t)[V i (t) + Ω i (t) × r i (t)] is theand φv p (x, t) = i=1
velocity ﬁeld of the colloids. The auxiliary concentrationﬁeld Cα∗ is introduced, which can have a ﬁnite value in theparticle domain, whereas Cα , the physical concentrationﬁeld, is forced to be zero through multiplication by (1−φ).
The advection of φ is solved via Ṙi = V i and by mapping {R1 , . . . , RNp } to φ. Henceforth, the volume of theﬂuid and/or the solid is strictly conserved and no numerical diﬀusion of φ occurs. In the SP method, the fundamental ﬁeld variables to be solved are taken as the total
velocity v rather than v f , and Cα∗ rather than Cα . This
choice of variables yields great beneﬁts in terms of allowing
eﬃcient and stable time evolution. The evolution equationof v is derived based on momentum conservation between
the ﬂuid and the particles, and the rigidity of the particle velocity ﬁeld v p . For the solute concentration, Cα andCα∗ diﬀer in terms of whether they exhibit abrupt variation at the interface of the colloids or not: especially forξ → 0, Cα has discontinuity at the solid-ﬂuid interface,
but Cα∗ should not. Since Cα exhibits abrupt variation in
the spatial scale of ξ, in order to solve its advection, it isnecessary to stabilize the evolution over time and to prevent numerical diﬀusion and penetration of Cα into theparticle domain. Compared with Cα , auxiliary concentration Cα∗ is independent of the abrupt variation arising from(1 − φ), and a ﬁnite value of Cα∗ in the particle domain isallowed. Therefore, it is numerically much easier to solve
the advection equation of Cα∗ than Cα .
3.1 Discretization in time
The time-discretized evolution equations are derived asfollows. To simplify the presentation, we neglect random
currents. As an initial condition at the n-th discretized
time step, the position, velocity, and angular velocity of
the colloids, {Rni , V ni , Ω ni } (i = 1, . . . , Np ), are mappedto φn and φn v np and the following conditions are set:v n = (1 − φn )v f + φn v np , satisfying the incompressibility condition on the entire domain, ∇ · v n = 0, andCα∗,n
, satisfying the charge neutrality condition, dx(1 −φn ) α Zα eCα∗,n + dx |∇φ| eσe = 0, where |∇φi |eσe represents the surface charge distribution of the colloids. The
current for the auxiliary concentration ﬁeld is deﬁned as
Cα∗ v α = Cα∗ v + (I − nn) · Cα∗ (−Γα ∇μα ) ,(16)

where n(x, t) is the unit surface-normal vector ﬁeld which
is deﬁned on the interface domain with a ﬁnite thickness ξ.
In this deﬁnition of the current, the no-penetration condition is directly assigned. The auxiliary concentrations areadvected by this current as tn +h
∗,n+1
∗,n
Cα
= Cα −
ds∇ · Cα∗ v α ,
(17)
tn

4

The European Physical Journal E

where h is a time increment, and tn = nh is the n-th
discretized time. The total velocity ﬁeld is updated using
a fractional step approach. First, the advection and thehydrodynamic viscous stress are solved,v∗ = vn + tn +h

Rn+1= Rni +
i


ds∇ ·

tn
 tn +h1(−pI + σ  ) − vv , (18)ρis directly
The resulting particle velocity ﬁeld φn+1 v n+1
p
enforced on the total velocity ﬁeld as
v n+1 = v ∗∗ +



tn +h

dsφf p ,tn +hdsφf p = φn+1 v n+1− v ∗∗ −ptn

dsV i ,

with the incompressibility condition, ∇ · v ∗ = 0. Alongwith the advection of the total velocity, the particle position is updated using the particle velocity. The electrostatic potential for the updated particle conﬁguration isdetermined by solving the following Poisson equation,/ ,∇2 Φn+1 = −ρn+1
e

(20)

=
(1 −
with 
the charge density ﬁeld, ρn+1e
φn+1 ) α Zα eCα∗,n+1 + |∇φn+1 |eσe . The momentum
change as a result of the electrostatic ﬁeld is solved as∇Φn+1 .
v ∗∗ = v ∗ − hρn+1e(21)

At this point, the momentum conservation is entirely
solved for the total velocity ﬁeld. The rest of the updatingprocedure applies to the rigidity constraint on the particle
velocity ﬁeld.
The hydrodynamic force and torque on the colloids
exerted by the ﬂuid are derived by assuming momentum
conservation between the colloids and the ﬂuid. The timeintegrated hydrodynamic force and torque over a period
h are equal to the momentum change over the particle
domain:
tn +hv ∗∗ − v np ,
dsF H
(s)
=
dxρφn+1
i
i

(22)

tn



tn +h
H
v ∗∗ −v np
dsN i (s) = dx r n+1
×ρφn+1ii. (23)

tn

With this and other forces on the colloids, the particle
velocity and angular velocity are updated as+ Mp−1 tn +htn

Ω n+1= Ω ni + I −1
p ·
i

tn +h

dsF H
i

tn

,
ds F ci + F ext
i


tn +h

(24)

dsN H
i

tn

+ I −1p · tn +hdsN ext
i .
tn

h
∇pp ,
ρ

(27)

(19)

tn

V n+1
= V ni + Mp−1
i

(26)

tn

(25)

where φf p represents the force density ﬁeld which imposes the rigidity constraint on the total velocity ﬁeld.
The pressure due to the rigidity of the particle is determined by the incompressibility condition, ∇ · v n+1 = 0,
which leads to the following Poisson equation for pp , viz.,− v ∗∗ )]. We note again that, on∇2 pp = hρ ∇ · [φn+1 (v n+1p
the l.h.s. of equations (22, 23), and (27), the integrands
H
FH
i , N i , and φf p are not explicitly calculated but theirtime integrals are solved. In other words, the solid-ﬂuid
interactions are treated in the form of the momentum
change, namely, momentum impulse.

3.2 Restriction on a time increment
To enable spatial discretization of the hydrodynamic equations, any standard scheme, such as the ﬁnite-diﬀerence
method, ﬁnite-volume method, ﬁnite-element method,
spectral method, lattice Boltzmann discretization and soforth, can be used. The SP method basically deﬁnes a coupling scheme between the hydrodynamic equations for thesolvent and the equations for the discrete colloids. Since
the treatment of the rigidity constraint of the particle velocity does not introduce an additional time scale, the restriction to a time increment h is the same as that in uniform ﬂuid cases. This is advantageous as compared to themethods adopted in references [7,22], where a large viscosity or elasticity is used for the velocity over the particle
domain. For comparison, in the Fluid Particle Dynamics (FPD) method [7], the large viscosity is introduced
for the ﬂuid particle ηc ( η) to realize the rigidity constraint. This means that the required time increment for
FPD should be very small, i.e., η/ηc ( 1), as comparedwith that used in the SP method.
A similar discussion on the restriction of the nopenetration condition in the advection-diﬀusion equationof the solutes can be outlined. One of the simplest treatments of the no-penetration condition in the particle domain is the penalty method adopted in references [23,24], in which an artiﬁcially large potential barrier is introduced for the solutes in the particle domain. The artiﬁcial
potential should at least be larger than the other chemical potentials in order to realize no-penetration of thesolutes. This restriction means that the artiﬁcial potential requires a smaller h. Although this strategy is physically consistent, it is numerically ineﬃcient. In contrast,in the SP method the advection-diﬀusion requires no additional time scale for the inclusion of colloids since the nopenetration condition for the solutes is directly assigned tothe solute current in the ﬁnite interface domain. From theY. Nakayama et al.: Simulating (electro)hydrodynamic eﬀects in colloidal dispersions: Smoothed proﬁle method

Fig. 1. Snapshot of a sedimenting colloid of radius a = 4Δx.
The arrows indicate the velocity ﬁeld.
simulation
Numerical (Zick & Homsy 1982)
low-ϕ asymptotics (Hasimoto 1959)

Drag coefficient: Q=-6πηa V/FD

0.9

FD = −

0.8

6πηaV
,
Q(ϕ)

(28)

where ϕ = (4/3)π(a/L)3 is the volume fraction in a cubicbox of volume L3 . The eﬀect of the boundary condition
extends very far in creep-ﬂow, as the higher the volume
fraction, the larger the drag. Comparison of the drag coeﬃcient given by the SP method with the analytic result
from the Stokes equation by Zick and Homsy [25] veriﬁes
the validity of the SP method over the entire range of ϕ.

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

Fig. 3. (Color online) Snapshot of two approaching colloidsof radius a = 4Δx. The arrows indicate the velocity ﬁeld. The
color map around the colloids represents the pressure distribution: a change in color from red to blue corresponds to achange from high to low pressure.

deﬁned as

1.0

5

0

0.1

0.2

0.3

0.4

0.5

0.6

volume fraction

Fig. 2. Drag coeﬃcient of a periodic array of spheres in steady
Stokes ﬂow as a function of volume fraction ϕ solved using the
SP method as compared with the analytic result [25] (solidline) and low-ϕ asymptotics [26] (dashed line). In this case,ξ is set to unity in the lattice unit. The range of the volume
fraction ϕ was obtained by changing the radius a from 4 to
15.5 for a lattice size Δx with L = 32Δx ﬁxed.

above discussion, we can see that the SP method provides
us with much higher numerical eﬃciency than other methods proposed for direct numerical simulation of colloidal
dispersions.

4 Results and discussions
4.1 Stokes drag on a periodic array of spheres in a
Newtonian ﬂuid
To validate the SP method quantitatively, we measured
the steady-state drag force on a periodic array of spheresin a Newtonian solvent. The velocity distribution aroundthe colloid is depicted in Figure 1. In general, ﬂow arounda colloid occurs as creep-ﬂow with a Reynolds number of
Re = aV /ν  1. Figure 2 shows the drag coeﬃcient Q(ϕ),

4.2 Lubrication interaction in a ﬁnite system
One of the most important eﬀects of solvent ﬂow is the
lubrication interaction between nearby particles with relative motions. The exact solution of the Stokes equation for
two isolated spheres has been found [27] and has provided
much insight into the basic physics of colloidal suspensions. However, its application to many-particle systems
through the method of pairwise addition requires care.
For quantitative prediction of the rheology of concentrated
suspensions, numerical results have identiﬁed many diﬀerences, whether the shear mode of the lubrication interaction is included or not [28,29]. There exists a fundamental
ambiguity in the application of the analytic expression of
two isolated spheres to a many-particle system in the ﬁnite
domain using pairwise addition.
We computed the squeeze lubrication interaction between two approaching spheres in a ﬁnite system. The
velocity and pressure distributions are depicted in Figure 3. Figure 4(a) shows the normalized approach velocity of a pair of particles versus the gap h between two
equal spheres solved by the SP method as compared withother theories. The two asymptotic solutions at h  aand h  a are from the exact solution of the isolated
pair [27] and Rotne-Prager-Yamakawa tensor [30], respectively. The Stokesian Dynamics (SD) solution [5] is basedon an interpolation of these two asymptotic solutions.
The simulation result nicely reproduces not only the two

6

The European Physical Journal E
1.6

3.5
(a)

1.4

3.0
2.5

1.0

eζ/kBT6 π η a |V1-V2| / F

1.2

0.8
0.6

2.0
1.5
1.0

0.4
0.2
0.0
0.0

simulation
RPY tensor(Ewald summed)
2-body exact resistance-tensor R2B
Stokesian Dynamics
0.5

1.0

1.5

2.0

2.5

0.5
3.0

h/a
100
(b)

Fig. 5. Non-dimensionalized zeta potential as a function of
the inverse Debye length κa and the valence of the colloid Z
computed using the semi-analytic formula of Ohshima-HealyWhite [37]. The dotted lines are guides to the eye.

γ1
γ2
γ1:isolated single sphere

γ 1 , γ2
1

0

0.1

0.2
0.3
volume fraction

0.0
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

κa

10

0.1

Z=1
10
50
100
500
1000

0.4

Fig. 4. (a) Relative velocity between two approaching spheres
versus the gap between the sphere surfaces (symbols). The
slight oscillations in the SP results are a result of the ﬁnite
lattice spacing. In this case, ξ is set to unity in the latticeunit. Theoretical curves are shown for far-ﬁeld asymptotics using the Rotne-Prager-Yamakawa tensor (dotted line) [30], anear-ﬁeld exact solution (dashed line) [27], and Stokesian Dynamics (solid line) [5]. (b) The two coeﬃcients representing thesqueeze interaction (scaled to those at inﬁnite dilution) versus
volume fraction. γ1 (square) represents the one-body drag andγ2 (circle) represents the two-body squeeze interaction due torelative motion.

asymptotic regimes but also the crossover which occurs at
h/a ∼ 0.7. It is found that SD underestimates the mobility
in this intermediate regime. This underestimation is dueto the approximation adapted in SD in which the components are just asymptotic two-body solutions. The result
conﬁrms the relevance of our simulation, thus demonstrating the importance of the hydrodynamic interaction in aﬁnite system.
We discuss the dependence of the lubrication interaction on the volume fraction. By assuming the functional
form |V1 −V2 |/F = 1/(Γ1 (ϕ)+2Γ2 (ϕ)/h) from lubrication
theory, where Γ1 (ϕ) = 6πηaγ1 (ϕ) is the one-body dragcoeﬃcient and Γ2 (ϕ) = (3πηa2 /2)γ2 (ϕ) is the squeeze coeﬃcient, the eﬀect of the volume fraction is represented
by γ1 (ϕ) and γ2 (ϕ). These friction coeﬃcients can be extracted by ﬁtting the curve of Figure 4(a) for each ϕ. The
reduced coeﬃcients γ1 and γ2 are plotted in Figure 4(b)as a function of the volume fraction. The solid line in
Figure 4(b) is from Figure 2. Although a periodic arrayof spheres exhibits diﬀerent ﬂow geometry from the caseof two approaching spheres, the volume fraction depen-

dence of γ1 of the two cases is comparable. Moreover, thesqueeze coeﬃcient γ2 is found to be a decreasing functionof the volume fraction. In other words, the squeeze mode
is most enhanced at inﬁnite dilution. In the literature [2],it is pointed out that the squeeze coeﬃcient is at least
smaller than that of the exact solution for isolated pairs.
These results further validate the SP method.
Although the SP method itself is an eﬃcient scheme asa direct simulation of utility in constructing a more coarsegrained model for suspensions, such as in dissipative particle dynamics, constitutive modeling, etc., the calculation
by direct simulation also gives fundamental information
about the hydrodynamic interactions.

4.3 Sedimenting charged colloids
As an example of the speciﬁc application of the method toa multi-component ﬂuid, we compute the hydrodynamic
drag for sedimenting charged colloids in the absence ofan external electric ﬁeld. In this case, the sedimenting
charged colloids induce a ﬂow that determines a charge
distribution which diﬀers from the equilibrium case [3,31].
The skewed ion distribution gives rise to polarization in
the electric double layer. Moreover, double-layer deformation induces an electro-osmotic ﬂow, which makes the ﬂow
around a charged colloid diﬀerent from that of a neutral
colloid.
In this simulation, a periodic array of charged colloids in a 1:1 electrolyte solution under gravity was computed. We specify the valence of the colloids ﬁrst, and thecounterion in the host ﬂuid assures the charge neutralityof the whole system. The bulk concentrations of the 1:1
electrolyte C̄± were determined by specifying the Debye
length κ−1 = {4πλB (C̄+ + C̄− )}−1/2 . The Bjerrum length
4πλB = e2 / kB T was set to unity and ξ was chosen as
twice the lattice unit to resolve the surface charge distribution, which is represented by |∇φ|. For simplicity, thecounterion and coion were set to the Schmidt number 0.5,
whose choice did not aﬀect the qualitative aspect of thefollowing results but only control the transient time to the

Y. Nakayama et al.: Simulating (electro)hydrodynamic eﬀects in colloidal dispersions: Smoothed proﬁle method

7

1.0

V(φ,κa)/V(φ)

0.9

0.8

0.7
Z=1
10
50
100
500
1000

0.6

0.5

0

0.5

1

1.5

2

2.5

3

3.5

4

κa

Fig. 6. Sedimentation velocity of a periodic array of colloidsof diﬀerent valences Z in electrolyte solutions as a functionof the inverse Debye length κa. The ordinate is reduced by
the sedimentation velocity of neutral colloids at the same volume fraction of 0.008 (a = 4Δx, L = 32Δx). The eﬀective
volume fraction including the electric double-layer deﬁned by(4/3)π{(a + κ−1 )/L}3 was chosen to be less than unity. Dottedlines are guides to the eye.

steady state. The Debye length was chosen so that the effective radius of the double layer a + κ−1 was smaller thanhalf of the system size L/2. In order to clarify the situation in terms of a typical colloid science parameter, thenon-dimensional zeta potential eζ/kB T is shown in Figure 5, which should be determined by specifying κa andthe valence of the colloid Z. The corresponding dimension of the zeta potential ζ for this simulation at 20 ◦ C isof the order of 10 mV. We computed the linear response
regime for gravitational driving to observe the eﬀect ofelectro-osmosis on the hydrodynamic drag.
We plot the sedimentation velocity scaled by that of aneutral colloid as a function of the inverse Debye lengthscaled by the radius of colloid in Figure 6. As the valenceof the colloid or the zeta-potential increases, the sedimentation velocity reduces to that of the neutral colloid. The
hydrodynamic drag of the electrolyte solution was mostenhanced when the Debye length was comparable to thesize of the colloids. These facts qualitatively agree withthe analytic result at inﬁnite dilution [32,33].
The charge density and the velocity distributions whenκa = 1 and Z = 1000 are depicted in Figure 7. In this case,the charge distribution was largely uniform and thus thecounteracting electrostatic eﬀect of the double-layer polarization does not have much eﬀect. Therefore, the enhanced
electrohydrodynamic drag can mainly be attributed to thefriction between the solvent and the ions. Because of thiselectrohydrodynamic coupling of the transfer in momentum, the solvent-ﬂow pattern around the colloid was modiﬁed from the case of a neutral colloid. As a result, the viscous drag on the colloid was enhanced. This mechanismof enhanced viscous drag by electro-osmotic ﬂow generally
exists in colloids in electrolyte solutions.
It has been known that in an inﬁnitely dilute system(or the thin double-layer limit) the velocity decays as r−3
in the region of κr  1 in contrast to the r−1 decay of
an inﬁnitely dilute neutral system, where r is the radial

Fig. 7. (Color online) Snapshot of a sedimenting charged colloid when κa = 1 and Z = 1000. The arrows indicate thevelocity ﬁeld. The color map around the particle represents
the charge density ﬁeld: a change in color from red to blue
corresponds to a change from high to low charge density.

distance from the center of the colloid [34,35]. However,in the system size adapted in our simulation this asymptotic regime was not reached. For κr  1, we have thescreened hydrodynamics regime, where the velocity decays as v ∝ e−κr /r [35]. These factors account for why theﬂow patterns in Figures 1 and 7 resemble one another. In
other words, the electrohydrodynamic interaction should
be a pronounced eﬀect in small ﬁnite systems as is seen inneutral systems.
We note that, at a ﬁnite volume fraction and in therange of the Debye length, the charged colloid dynamics
could be eﬀectively assessed through the SP method. Further application of the SP method to electrophoresis ofconcentrated suspensions is reported elsewhere [36].

5 Conclusions
We have presented a new simulation scheme for colloidal
dispersions in a solvent of a multi-component ﬂuid, which
we call the “Smoothed Proﬁle (SP) method”. The SP
method improves and extends a previous method proposedby the authors [38,20,39].
The description of the colloidal systems is based on
the Navier-Stokes equation for the momentum evolutionof the solvent ﬂow, the advection-diﬀusion equation forthe solute distribution, the rigid-body description of colloidal particles, with dynamical coupling of all these elements. Based on momentum conservation between thecontinuum solvent ﬂuid and the discrete rigid colloids,
the time-integrated hydrodynamic force and torque arederived. This expression of the mechanical coupling between a ﬂuid and particles is well suited for numerical simulations in which diﬀerential equations are discretized intime. With this formulation relating the solid-ﬂuid interaction, standard discretization schemes for uniform ﬂuids

8

The European Physical Journal E

are utilized as is; no special care is needed for the solidﬂuid boundary mesh. Since the hydrodynamic interaction
is solved through direct simulation of the solvent ﬂuid, themany-body eﬀect can be fully resolved.
The utility of the SP method was assessed in various test problems, not only using simple ﬂuids but also
simulating charged colloids in electrolyte solutions. The
results conﬁrmed that the SP method is eﬀective forstudying the dynamical behavior of colloids. Although wehave focused on systems of spherical colloids in simple ﬂuids and Poisson-Nernst-Planck electrolytes (i.e. PoissonBoltzmann level description of electrolytes), application
to other types of macromolecules of other shapes, such as
disks [19], rods, and others, and other constitutive models
for solvents is straightforward.
Note added in proofs: After this article was written, incorporation of thermal ﬂuctuation with the SP method byusing a generalized Langevin equation has been discussed
in [40] to simulate the Brownian motion of colloids.

References
1. J. Happel, H. Brenner, Low Reynolds Number Hydrodynamics: With Special Applications to Particulate Media,
2nd edition (Martinus Nijhoﬀ, Dordrecht, 1983).
2. S. Kim, S.J. Karrila, Microhydrodynamics: Principles andSelected Applications (Butterworth-Heinemann, London,1991).
3. W.B. Russel, D.A. Saville, W.R. Schowalter, Colloidal Dispersions (Cambridge University Press, Cambridge, UK,1989).
4. C.N. Likos, Phys. Rep. 348, 267 (2001).
5. J.F. Brady, G. Bossis, Annu. Rev. Fluid Mech. 20, 111(1988).
6. A. Malevanets, R. Kapral, J. Chem. Phys. 110, 8605
(1999).
7. H. Tanaka, T. Araki, Phys. Rev. Lett. 85, 1338 (2000).
8. T. Kajishima, S. Takiguchi, H. Hamasaki, Y. Miyake,
JSME Int. J., Ser. B 44, 526 (2001).
9. H.H. Hu, N.A. Patankar, M.Y. Zhu, J. Comput. Phys. 192,
427 (2001).
10. R. Glowinski, T.W. Pan, T.I. Hesla, D.D. Joseph, J.
Périaux, J. Comput. Phys. 192, 363 (2001).
11. A.J.C. Ladd, R. Verberg, J. Stat. Phys. 104, 1191 (2001).

12. J.T. Padding, A.A. Louis, Phys. Rev. Lett. 93, 220601(2004).
13. M.E. Cates, K. Stratford, R. Adhikari, P. Stansell, J.C. Desplat, I. Pagonabarraga, A.J. Wagner, J. Phys.: Condens.
Matter 16, S3903 (2004).
14. Y.W. Kim, R.R. Netz, Europhys. Lett. 72, 837 (2005).
15. T. Yamaue, M. Sasaki, T. Taniguchi, Multi-Phase Dynamics Program “Muﬃn” User’s Manual, http://octa.jp(2005).
16. V. Lobaskin, B. Dünweg, C. Holm, J. Phys.: Condens.
Matter 16, S4063 (2004).
17. A. Chatterji, J. Horbach, J. Chem. Phys. 122, 184903(2005).
18. F. Capuani, I. Pagonabarraga, D. Frenkel, J. Chem. Phys.
121, 973 (2004).
19. F. Capuani, I. Pagonabarraga, D. Frenkel, J. Chem. Phys.
124, 124903 (2006).
20. Y. Nakayama, R. Yamamoto, Phys. Rev. E 71, 036707(2005).
21. L.D. Landau, E.M. Lifshitz, Fluid Mechanics (Pergamon
Press, London, 1959).
22. C.S. Peskin, D.M. McQueen, J. Comput. Phys. 81, 372(1989).
23. J. Dzubiella, H. Löwen, C.N. Likos, Phys. Rev. Lett. 91,
248301 (2003).
24. H. Kodama, K. Takeshita, T. Araki, H. Tanaka, J. Phys.:
Condens. Matter 16, L115 (2004).
25. A.A. Zick, G.M. Homsy, J. Fluid Mech. 115, 13 (1982).
26. H. Hasimoto, J. Fluid Mech. 5, 317 (1959).
27. D.J. Jeﬀrey, Y. Onishi, J. Fluid Mech. 139, 261 (1984).
28. R.C. Ball, J.R. Melrose, Physica A 247, 444 (1997).
29. J.R. Melrose, R.C. Ball, J. Rheol. 48, 937 (2004).
30. C.W. Beenakker, J. Chem. Phys. 85, 1581 (1986).
31. R.F. Probstein, Physicochemical Hydrodynamics: An Introduction, 2nd edition (John Wiley & Sons, New York,
2003).
32. F. Booth, J. Chem. Phys. 22, 1956 (1954).
33. H. Ohshima, T.W. Healy, L.R. White, R.W. O’Brien, J.
Chem. Soc. Faraday Trans. 2 80, 1299 (1984).
34. J.L. Anderson, Annu. Rev. Fluid Mech. 21, 61 (1989).
35. D. Long, A. Ajdari, Eur. Phys. J. E 4, 29 (2001).
36. K. Kim, Y. Nakayama, R. Yamamoto, Phys. Rev. Lett.
96, 208302 (2006).
37. H. Ohshima, T.W. Healy, L.R. White, J. Colloid Interface
Sci. 90, 17 (1982).
38. R. Yamamoto, Phys. Rev. Lett. 87, 075502 (2001).
39. K. Kim, R. Yamamoto, Macromol. Theory Simul. 14, 278(2005).
40. T. Iwashita, Y. Nakayama, R. Yamamoto, J. Phys. Soc.
Jpn. 77 (2008) in press.