PHYSICAL REVIEW LETTERSPRL 96, 208302 (2006)

week ending
26 MAY 2006

Direct Numerical Simulations of Electrophoresis of Charged Colloids
Kang Kim, Yasuya Nakayama,* and Ryoichi Yamamoto
Department of Chemical Engineering, Kyoto University, Kyoto 615-8510, Japan,
and PRESTO, Japan Science and Technology Agency, Kawaguchi 332-0012, Japan(Received 27 January 2006; published 26 May 2006; publisher error corrected 8 June 2006)
We propose a numerical method to simulate electrohydrodynamic phenomena in charged colloidal
dispersions. This method enables us to compute the time evolutions of colloidal particles, ions, and host
fluids simultaneously by solving Newton, advection-diffusion, and Navier-Stokes equations so that theelectrohydrodynamic couplings can be fully taken into account. The electrophoretic mobilities of charged
spherical particles are calculated in several situations. The comparisons with approximation theories showquantitative agreements for dilute dispersions without any empirical parameters; however, our simulation
predicts notable deviations in the case of dense dispersions.
DOI: 10.1103/PhysRevLett.96.208302

PACS numbers: 82.70.Dd, 47.65.d, 82.20.Wt, 82.45.h

Electrohydrodynamic phenomena are of great importance in physical, chemical, and biological science, andalso in several engineering fields [1]. In the case of electrophoresis of charged particles, for example, the particles
start to move on the application of external electric fields.
The electric double layer, i.e., the cloud of counterions
around charged particles, tends to be deformed and its
distribution becomes anisotropic because of the appliedexternal field and also of the friction between ions andfluids. The electrophoretic mobility of a single colloidal
particle is then determined by the balance between theelectrostatic driving force and the hydrodynamic frictional
force acting on the particle. In this situation, the time
evolutions of the colloidal particles, the ions, and thehost fluids are described by coupled equations of hydrodynamics (Navier-Stokes) and electrostatics (Poisson) withproper boundary conditions imposed on the surfaces of thecolloidal particles. However, the usual numerical techniques of partial differential equations are hopeless indealing with the dynamical evolution of many-particle
systems since the moving particle-fluid boundary condition
must be treated at every discrete time step.
In recent years, efforts to resolve hydrodynamic interactions in colloidal dispersions have attracted a lot of
attention. Various advanced methods have been proposedsuch as the Stokesian dynamics (SD) [2], the finite element method (FEM) [3], the lattice Boltzmann method(LBM) [4], the stochastic rotation dynamics [5], the fluid
particle dynamics (FPD) [6], and yet another method
which treats solid-fluid interaction efficiently [7].
Pioneering approaches have been proposed also to simulate charged colloidal dispersions without hydrodynamics [8–12]. Although extensions have been done totake into account the hydrodynamics by using SD [13],
FEM [14], LBM [15–19], and FPD [20], quantitatively
reliable simulations have not yet been performed successfully for many-particle dispersions due to the complexityof the system.
Recently, we developed a reliable and efficient numerical method, called the smoothed profile (SP) method
0031-9007=06=96(20)=208302(4)[21,22], to resolve the hydrodynamic interactions acting
on solid particles immersed in Newtonian host fluids. In the
SP method, the original sharp boundaries between colloids
and host fluids are replaced with diffuse interfaces withfinite thickness . This simple modification greatly improves the performance of numerical computations since
it enables us to use the fixed Cartesian grid even for the
problems with moving boundary conditions.
The SP method is not only applicable to the dispersions
in Newtonian fluids, but particularly suitable for the particle dispersions in complex fluids. It has already been
applied successfully to liquid crystal colloidal dispersions[23,24] and charged colloidal dispersions [25]. Fieldparticle hybrid simulations were performed where theaverage direction of the liquid crystal molecules and thedensity of ions were treated as coarse-grained continuum
objects while colloids were treated explicitly as particles.
The interaction between fields and particles was taken
through the diffuse interface. The above methods for thedispersions in complex fluids are, however, not yet appropriate for simulating dynamical phenomena since hydrodynamic effects are completely neglected. The purpose of
the present study is to establish an efficient and reliable
simulation method applicable for electrohydrodynamic
phenomena such as electrophoresis by combining our SP
methods for hydrodynamic [21,22] and electrostatic [25]interactions.
In the present Letter, we briefly outline our numericalmodeling for charged colloidal dispersions and then demonstrate the reliability of the combined SP method by comparing our numerical results with classical approximation
theories [26 –29]. Finally, comparisons are made for theelectrophoretic mobilities of dense dispersions, where thesimulation results show notable deviations from a meanfield–type theory according to the cell model [30,31].
Let us consider N spherical particles with radius a, themass Mp , and the inertia tensor Ip in a host fluid consistingof multicomponent ions of species  with charges Z e,
where e is the unit charge. The local number density of ions is C r~; t at a time t. The total charge on a colloidal

208302-1© 2006 The American Physical Society

PRL 96, 208302 (2006)

particle is Ze, distributed uniformly on its surface. The~ r; t. The position, thevelocity field of the host fluid is v~translational velocity, and the angular velocity of the ith~ i , respectively. We define theparticle are R~i , V~i , and Poverall profile function ~r; t  N
r; t, wherei1 i ~
i 2 0; 1 is the ith particle’s profile field, which is unity
in the particle domain jr~  R~i j < a  =2, zero in the fluid
domain j~r  R~i j > a  =2, and has a continuous diffuse
interface within the thin interface domain a  =2 < jr~ 
R~i j < a  =2 whose thickness is . The mathematicaldefinition of i is given in Ref. [21]. We define the spacial
distribution of the surface charge eq~r  Zejrj=4a2
using ; then the local density of the total charge is~ 
represented smoothly everywhere in the system by e r
P Z eC  eq. The complete dynamics of the system isobtained by solving the following time evolution equations[21,22].
(i) The Navier-Stokes equation:
@t  v~ rv~  rp  r2 v~  e r  ex   f~p ;(1)
with the incompressible condition r v~  0, where  isthe mass density, p is the pressure,  is the shear viscosityof the host fluid, ex  E~ r~ is the external electric~ and f~p
potential due to the constant electric field E,
represents the body force arising from the rigidity of theparticles [22]. The electrostatic potential r~ is to be
determined by solving the Poisson equation r2  e with the dielectric constant  of the host fluid.
(ii) Newton’s and Euler’s equations of motion:_
Hc_
R~ i  V~ i ; Mp V~ i  F~ i  F~ i ;week ending
26 MAY 2006

PHYSICAL REVIEW LETTERS~_ i  N~ H
Ip 
i ;

  kB T lnC  Z e  ex ;

where kB is the Boltzmann constant and T is the temperature. If we set v~  0 in Eq. (3), the equilibrium (t ! 1)
ionic density is given by the Boltzmann equationC  C  expZ e  ex =kB T:~H
where F~ H
i and N i are the hydrodynamic force and torquec[22], and F~ i is the force arising from the excluded volumeof particles which prevents colloids from overlapping.
Hereafter, soft-core potential of the truncated LennardJones potential is adopted for F~ ci . We include the electric
driving force due to E~ in F~ Hi .
(iii) Advection-diffusion equation:~ C r ; (3)@t C  r C v~  f1 r I  n~ n
where f is the ionic friction coefficient, I is the unit
tensor, and n~ is a unit-vector field defined by n~ r=jrj. In our method, the actual density fields ofions are defined as C r~; t  1  r~; tC ~r; t using
the auxiliary density field C r~; t. This definition avoids
penetration of ions into colloids explicitly without using
artificial potentials, which requires smaller time incre~ in Eq. (3) ensures the conments. The operator I  n~ n
servation of C since the no-penetration conditionn~ r  0 is directly assigned
at the diffuse interface.
R
Then the charge neutrality e dr~  0 of the total system
is guaranteed automatically. Based on the density functional theory [32,33], the chemical potential is defined as(5)

The bulk salt concentrationq
C  is related to the DebyeP1 1= 4 B  Z2 C  , where B screening length
e2 =4kB T is the Bjerrum length, which is approximately
0.72 nm in water at 25 C.
Simulations have been performed in a three-dimensional
cubic box with the periodic boundary conditions. The
linear dimension is L=  64, where  is the lattice spacing chosen as a unit length, which is taken related to the
Bjerrum length as   4 B . We use the particle radiusa  5 and the thickness parameter   2 throughout thepresent simulations. The host fluid contains 1:1 electrolytecomposed of monovalent counterions (  ) and coions(  ). The unit of the energy and the electrostatic
potential is kB T and kB T=e, respectively. The later corresponds to 25.7 mV at 25 C. The nondimensional parameter m  2kB Tf =3e2 is set to m  m  0:184,
which corresponds to the value of KCl solution at 25 C.
Our unit of time  2 f =kB T corresponds to 0:44 s.
We first consider a single charged particle moving withthe drift velocity V~  V; 0; 0 in a constant electric fieldE~  E; 0; 0. The electrophoretic mobility V=E is related
to the zeta potential , which is defined as the electrostatic
potential on a slipping plane as
V=E  f =

(2)

(4)

(6)

when is small [1]. The Smoluchowski equation f  1 isvalid in the limiting case a ! 1 [27], while the Hückel
equation f  2=3 is valid in the opposite case a ! 0[28]. Henry derived an expression f  fH  a which isvalid for a general value of a [29]. These equations
indicate that the mobility is proportional to ; however,this relation tends to fail for larger where the relaxation
effect due to deformations of electric double layer becomesnotable. O’Brien and White proposed an approximation
theory which is valid also for larger [26].
We have performed simulations for electrophoresis of asingle particle in linear response regimes E & 0:15 andcompared them with the O’Brien-White theory. A constantuniform electric field E  0:1, which corresponds to
2:85 103 V=cm, was applied. The terminal V was calculated for 50 Z 750 with 1  10, corresponding
to the salt concentration 11 mol=l at 25 C in water. We
chose  =  5, so the Reynolds number Re  aV=
remains small. In both the O’Brien-White theory and oursimulations, the zeta potential is assumed to be the electrostatic potential at the particle surface,  jsurface . In our
simulations, the surface charges were chosen as Z  50,
100, 200, 300, 400, 500, and 750, corresponding to y  0:525, 1.044, 2.035, 2.927, 3.692, 4.332, and208302-2

PRL 96, 208302 (2006)

PHYSICAL REVIEW LETTERS

5.510, respectively. Here the dimensionless zeta potentialy  e =kB T is introduced. A relationship between thesurface charge jZje and the dimensionless zeta potential
y is shown in Fig. 1(a), where our numerical results areplotted together with an analytic solution of the nonlinear
Poisson-Boltzmann (PB) equation [34] and the DebyeHückel linearized solution  jZje=4a2  1  a1 .
We see that our numerical results are consistent with thenonlinear PB theory. In Fig. 1(b), the dimensionless mobility Em  3eV=2kB TE is plotted as a function of thedimensionless zeta potential with a  0:5. It is clearly
demonstrated that our method reproduced the O’BrienWhite theory almost perfectly including the nonlinear
regime y  2 within only a few percent error. We emphasize that such a precise agreement with the theory has neverbeen obtained by any simulation methods so far proposed.
The distributions of charge density due to counterions andcoions are shown in Fig. 1(c) for y  1:044 and (d) for y 
3:692. One can see that the electric double layer is deformed considerably in the nonlinear regime (d), while it isalmost isotropic in the linear regime (c). As mentioned
before, the relaxation effect due to the deformed double
layer causes the nonlinearity in Em .
Our simulation method is easily applicable to dense
dispersions consisting of many particles. We thus examined the effect of the particle concentration on the electrophoretic mobility. The linearized theory for a single

FIG. 1 (color). Relationship between surface charge jZje anddimensionless zeta potential y (a). Our numerical data follow
nicely on the analytic solution of the nonlinear PB equation [34]but deviates notably from the Debye-Hückel linearized theory.
Dimensionless mobility Em of a single particle is plotted in (b) as
a function of dimensionless zeta potential y. For comparison,
results of Smoluchowski, Henry, and O’Brien-White for a 
0:5 are plotted. The colorPcontours in (c) and (d) represent the
total ionic charge density  eZ C around a single particle for(c) Z  100 and (d) Z  500. The electric field is applied in
the horizontal (  x) direction.

week ending
26 MAY 2006

particle Eq. (6) is still valid for dense dispersions when E
is small; however, f is now depending on both a and ’.
Simulations were carried out with Z  100 and E  0:1
for various particle volume fractions ’  4a3 N=3L3 to
calculate f a; ’  V= E. The Debye length 1 istaken to be 5 and 10, which correspond to a  1 and 0.5,
respectively. The corresponding salt concentration is
44 mol=l for 1  5 and 11 mol=l for 1  10,
respectively. Figure 2 shows typical snapshots of the systems with (a) fcc, (b) bcc, and (c) random configurations[35]. The horizontal color map represents the charge density for a  1 on a cross section perpendicular to the zaxis. In the cases of fcc and bcc, E was applied perpendicular to the (1,0,0) and (1,1,1) faces, but we obtained
very small differences within only 1%.
The mobility coefficient f a; ’ for a  1 and 0.5 isplotted as a function of ’ in Figs. 3(a) and 3(b), respectively. We found that f decreases rapidly with increasing’. Furthermore, the overall behavior looks almost independent of the particle configurations. A theoretical model
has been proposed by Levine and Neale for the electrophoretic mobility of dense dispersions by using the cell
model [30]. They assumed the situation that a sphericalparticle with radius a is located at the center of a sphericalcontainer (cell) with radius b and calculated V as a function of a and ’  a=b3 . Ohshima proposed a simpler
expression for the mobility coefficient f according to thecell model [31]. Ohshima’s prediction is shown in
Figs. 3(a) and 3(b) together with our numerical results.
The overall agreement between our simulation andOhshima’s theory is better in (a) with a smaller Debye
length 1  5  a than in (b) with a larger one 1 
10. In both (a) and (b), the simulation results tend to belarger than the theory as ’ increases. We expect that thedeviation is attributable to the occurrence of overlapping of
the electric double layers for larger ’ because such aneffect is totally neglected in the theory. To this end, we
estimated the effective radius a  1 of the ionically
dressed particles and defined the effective volume fraction’eff  4a  1 3 N=3L3  1   a1 3 ’. As isclearly seen in Figs. 3(a) and 3(b), our results agree wellwith Ohshima’s theory for ’eff 1 where the effect of

FIG. 2 (color). Snapshots of the electrophoresis of dense dispersions with (a) fcc, (b) bcc, and (c) random particle configurations.
The color map represents the total ionic charge densityPeZC  in a plane perpendicular to the z axis. The electricfield is applied in the x direction normal to (1,0,0) face for fcc
and bcc. See movies in Ref. [35].

208302-3

PHYSICAL REVIEW LETTERSPRL 96, 208302 (2006)
0

10

10

0

(a)

10

(b)

−1

−1

10

ϕeff =1ffϕeff =1

10

−2

−2

10
fcc
bcc
randomOhshima(κa=1)

−3

10 0

0.2

ϕ

0.4

fcc
bcc
randomOhshima(κa=0.5)

−3

0.6

week ending
26 MAY 2006

10 0

0.2

ϕ

0.4

0.6

FIG. 3 (color online). The mobility coefficient f a; ’ as afunction of the volume fraction ’ in (a) a  1 and(b) a  0:5. The solid lines represent the approximation theory
proposed by Ohshima [31]. The theory is confirmed to beaccurate for ’eff 1, but tends to deviate from our numericalresults for ’eff > 1 where overlapping of the electric double
layers becomes notable.

overlapping is small. However, for ’eff > 1 where theoverlapping of the electric double layers becomes large,
deviations between our simulations and the theory become
notable. We emphasize that the present study is the firstsuccessful simulations which provide quantitative datanecessary to examine the reliability of the Ohshima’s cell
model calculations including their boundary conditions forelectrophoresis in dense colloidal dispersions. Our results
are consistent with recent studies which also to take into
account the effects of double layer overlapping [19,36,37].
In summary, we have developed a unique numericalmethod for simulating electrohydrodynamic phenomena
in colloidal dispersions. The method was first applied tosimulate electrophoresis of a single spherical particle, andwe found that our method can reproduce quantitatively thereliable analytical theory proposed by O’Brien and White.
Simulations were then performed for electrophoresis ofcolloids in dense dispersions, and we compared themwith the theoretical analysis based on the cell model. We
found that the cell model is reliable when overlapping ofelectric double layers is small but less reliable as overlapping becomes larger.

*Present address: Department of Chemical Engineering,
Kyushu University, Fukuoka 819-0395, Japan.
[1] W. B. Russel, D. A. Saville, and W. R. Schowalter,
Colloidal Dispersions (Cambridge University Press,
Cambridge, 1989).
[2] J. F. Brady and G. Bossis, Annu. Rev. Fluid Mech. 20, 111(1988).
[3] H. H. Hu, N. A. Patankar, and M. Y. Zhu, J. Comput. Phys.
169, 427 (2001).
[4] A. J. C. Ladd and R. Verberg, J. Stat. Phys. 104, 1191
(2001).
[5] A. Malevanets and R. Kapral, J. Chem. Phys. 110, 8605
(1999).
[6] H. Tanaka and T. Araki, Phys. Rev. Lett. 85, 1338 (2000).

[7] T. Kajishima, S. Takiguchi, H. Hamasaki, and Y. Miyake,
JSME Int. J., Ser. B 44, 526 (2001).
[8] M. Fushiki, J. Chem. Phys. 97, 6700 (1992).
[9] H. Löwen, P. A. Madden, and J.-P. Hansen, Phys. Rev.
Lett. 68, 1081 (1992).
[10] H. Löwen, J.-P. Hansen, and P. A. Madden, J. Chem. Phys.
98, 3275 (1993).
[11] J. Chakrabarti, J. Dzubiella, and H. Löwen, Phys. Rev. E
70, 012401 (2004).
[12] J. Dobnikar, Y. Chen, R. Rzehak, and H. H. von Grünberg,
J. Chem. Phys. 119, 4971 (2003).
[13] Y. W. Kim and R. R. Netz, Europhys. Lett. 72, 837 (2005).
[14] T. Yamaue, M. Sasaki, and T. Taniguchi, Multi-Phase
Dynamics Program ‘‘Muffin’’ User’s Manual (OCTA
Users Group, Tokyo, 2005), http://octa.jp.
[15] J. Horbach and D. Frenkel, Phys. Rev. E 64, 061507(2001).
[16] R. Z. Wan, H. P. Fang, Z. Lin, and S. Chen, Phys. Rev. E
68, 011401 (2003).
[17] V. Lobaskin and B. Dünweg, J. Phys. Condens. Matter 16,
S4063 (2004).
[18] A. Chatterji and J. Horbach, J. Chem. Phys. 122, 184903(2005).
[19] V. Lobaskin, B. Dünweg, M. Medebach, T. Palberg, andC. Holm, cond-mat/0601588.
[20] H. Kodama, K. Takeshita, T. Araki, and H. Tanaka,
J. Phys. Condens. Matter 16, L115 (2004).
[21] Y. Nakayama and R. Yamamoto, Phys. Rev. E 71, 036707(2005).
[22] Y. Nakayama, K. Kim, and R. Yamamoto, cond-mat/
0601322.
[23] R. Yamamoto, Phys. Rev. Lett. 87, 075502 (2001).
[24] R. Yamamoto, Y. Nakayama, and K. Kim, J. Phys.
Condens. Matter 16, S1945 (2004).
[25] K. Kim and R. Yamamoto, Macromol. Theory Simul. 14,278 (2005).
[26] R. W. O’Brien and L. R. White, J. Chem. Soc., Faraday
Trans. 2 74, 1607 (1978).
[27] M. Smoluchowski, Z. Phys. Chem. 92, 129 (1917).
[28] E. Hückel, Phys. Z. 25, 204 (1924).
[29] D. C. Henry, Proc. R. Soc. A 133, 106 (1931).
[30] S. Levine and G. H. Neale, J. Colloid Interface Sci. 47, 520(1974).
[31] H. Ohshima, J. Colloid Interface Sci. 188, 481 (1997).
[32] J.-P. Hansen and H. Löwen, Annu. Rev. Phys. Chem. 51,
209 (2000).
[33] J.-L. Barrat and J.-P. Hansen, Basic Concepts for Simple
and Complex Liquids (Cambridge University Press,
Cambridge, 2003).
[34] H. Ohshima, T. W. Healy, and L. R. White, J. Colloid
Interface Sci. 90, 17 (1982).
[35] See EPAPS Document No. E-PRLTAO-96-051621 foranimation corresponding to Fig. 2. For more information
on EPAPS, see http://www.aip.org/pubservs/epaps.html.
[36] F. Carrique, F. J. Arroyo, M. L. Jiménez, and Á. V.
Delgado, J. Phys. Chem. B 107, 3199 (2003).
[37] T. Palberg, M. Medebach, N. Garbow, M. Evers, A. B.
Fontecha, H. Reiber, and E. Bartsch, J. Phys. Condens.
Matter 16, S4039 (2004).

208302-4