KAPSEL: Kyoto Advanced Particle Simulator for
ELectrohydrodynamics
─Toward Direct Numerical Simulations of Colloidal Dispersions─†R. Yamamoto1, K. Kim2, and Y. Nakayama3
Department of Chemical Engineering, Kyoto University＊

Abstract
We have released a colloid simulator named KAPSEL implemented the “Smoothed Profile (SP)method” which has been developed by ourselves for direct numerical simulations of particulate flow
providing a way to couple continuum fluid dynamics with rigid-body dynamics through smoothed
profile of colloidal particle. KAPSEL enebles us to simulate multi-component fluids, such systems as
charged colloids in electrolyte solutions. Dynamics of colloidal dispersions is solved as much computational cost as required for solving non-particulate flows. KAPSEL computes the fluid velocity andthe electrostatics potential by solving both Navier–Stokes and Poisson equations directly. The time
evolutions of the colloidal particles and the density of counter ions are then determined by solving
Newton’s equation of motion and advection-diffusion equation, respectively, in a consistent mannerso that the electro-hydrodynamic coupling can be fully taken into account. The electrophoretic mobility of spherical colloidal particles is calculated in several situations including those in concentrated
dispersions. The comparisons with theories show excellent quantitative agreements.
Keywords: Colloidal dispersion, Simulation, Rheology, Electrophoresis, Hydrodynamics, Sedimentation

1. Introduction
Interparticle interactions in colloidal dispersions
mainly consist of thermodynamic potential interaction and hydrodynamic interactions 2). Whereas theformer works in both static and dynamic situations,the latter works solely in dynamic situation. Although
the thermodynamic interactions in static situations
have been studied extensively and summarized asa concept of effective interactions, the nature of dynamic interactions are poorly understood. Since the
†
＊
１

２

３

Accepted: June 23, 2006
Kyoto 615-8510, Japan
Corresponding author,
TEL&FAX: +81-75-383-2682
E-mail: ryoichi@cheme.kyoto-u.ac.jp
Present address: Institute for Molecular Science, Okazaki 444-8585, Japan
Present address: Department of Chemical Engineering,
Kyushu University, Fukuoka 819-0395, Japan

KONA No.24 (2006)

hydrodynamic interaction is essentially long-ranged,
many-body, it is extremely difficult to study its roleby means of analytical way alone. Numerical simulations can aid to investigate fundamental role of thehydrodynamic interaction in colloidal dynamics.
Between November 2002 and March 2006, the“Development of Hybrid Molecular Dynamics Simu”4
lation Project, sponsored by the Japan Science andTechnology Agency,5 was carried out initially at the
Department of Physics and Astrophysics,6 Graduate
School of Science, Kyoto University and later at the
Depar tment of Chemical Engineering, 7 Graduate
School of Engineering, Kyoto University. As the ma-

４

５
６
７

http://www-tph.cheme.kyoto-u.ac.jp/hymd/index-e.
html
http://www.jst.go.jp/EN/
http://www.scphys.kyoto-u.ac.jp/index-e.html
http://www.cheme.kyoto-u.ac.jp/english/

167

jor objective of this project, we developed a reliable
and ef ficient numerical method, called smoothed
profile (SP) method 5) 6) 10) and 11), to resolve thehydrodynamic interactions acting on solid particles
immersed in complex host fluids. In the SP method,the original sharp boundaries between colloids andhost fluids are replaced with diffuse interfaces withfinite thickness ξ . This simple modification greatly
improve the performance of numerical computations
since it enables us to use the fixed Cartesian grideven for problems with moving boundary conditions.
In 2006, we have released a program package
named KAPSEL, the Kyoto Advanced Particle Simulator for Electrohydrodynamics, which is based on
the SP method and applicable to several electrohydrodynamic phenomena such as electrophoresisof charged colloidal particles. The entire KAPSEL
source code is viewable on the web 8 and anyone
agreeing to the “OCTA Free Software License Agreement”9 can download it freely for use. This document
outlines the basic principles of the SP method and tutorials describing how to install and execute KAPSEL
using a few sample test cases. The latest information
will be appearing online, so the reader is recommended to refer to the web if the instructions described inthis paper return unsuccessful results.
2. Installation
2.1 Installation of OCTAKAPSEL is designed to be used as an add-on engine of “OCTA10,” an integrated simulation system forsoft materials developed by the joint project of industr y and academia funded by Ministr y of Economy,
Trade and Industry (METI), Japan. To use KAPSEL,
one needs a graphical user interface called “Gourmet”
distributed within the “OCTA” program package,
whose newset version 11 can be obtained from theweb. One must install “OCTA” properly by following
the installation procedure of “OCAT” on each machine.
Because KAPSEL needs library called “libplatform,”
whose source codes are attached to “Gourmet” toimplement data input/output through UDF files, onemust follow the instructions below to build “libplatform” on each machine in advance of compilation of
KAPSEL.

% cd src% make% make install
Windows: For Windows machines, first install thefree software Cygwin12 that emulates UNIX on the
Windows. In doing so, be sure to include “gcc,” “gccmingw,” and “make.” Subsequently, start up the Cygwin shell and build “libplatform” as follows.
% ln -s /cygdrive/c/OCTA2006/ /usr/local/.
% cd /cygdrive/c/OCTA2006/GOURMET_2006% ./GOURMET2006_SRC.EXE% cd src% make -f Makefile.cygmingw% make -f Makefile.cygmingw install% cd ../lib/cygwin% ln -s ../win32/libplatform.a .
Mac: For MacOS X machines, follow the instruction below to build “libplatform.” Here we assume“OCTA2006” is installed in “/Application/OCTA/”
and two envir onment variables, “PATH” and“PF_FILES” , are set properly.
% cd /Application/OCTA/GOURMET_2006/src% make% make install
2.2 Installation of KAPSELThe complete source codes of the latest version of
KAPSEL can be obtained from the web as a “gzippedtar” archive file under the name of “kapsel.tar.gz.”
The file can be decompressed into the “./kapsel/” using the following command.
% tar xvfpz kapsel.tar.gz
Linux and Windows: The compilation of KAPSEL
can be made as follows using the supplied “Makefile”
without any changes.
% cd kapsel% makeMac: For MacOS X machines, put the following in
the “Makefile,” and then perform “make.”
ARCH=macosx
CC=gcc
CXX=g++
８

Linux: For Linux machines, follow the instruction
below to build “libplatform.”% cd /usr/local/OCTA2006/GOURMET_2006% tar xvfpz gourmet*_src.tar.gz

168

９
10
11

12

http://www-tph.cheme.kyoto-u.ac.jp/kapsel/
http://octa01.t.u-tokyo.ac.jp/OCTA/license_eng.html
http://octa.jp/
At this moment, the newest version is “OCTA2006.” We
further assume “OCTA2006” is properly installed.
http://cygwin.com/

KONA No.24 (2006)

CCOPT=-DMACOSX -03 -fno-inline
LINKS=-lm -lplatform
2.3 User support
The BBS of the OCTA home page http://octa.jp
contains the bulletin board of KAPSEL. Enter questions or bug reports about KAPSEL in the bulletinboard.

3.1 Basic Equations
3.1.1 Equation of Motion of Colloidal Particles
Let us consider the motion of N spherical particles
in a Newtonian fluid. The position, velocity, and angular velocity {Ri, Vi, Ωi} of the ith colloidal particle
having mass Mp evolve over time according to theequations of motion:
Ri=Vi,(2)
MpVi=F Hi +F other
,
(3)
i
H
Ip･Ωi=N i ,(4)
Ip=(2/5) Mpa2I denotes the moment of inertia of thecolloidal particle, with I denoting the unit tensor. In
the above equations, FHi and NHi respectively denotethe force and torque received by the colloidal particle
from the fluid, expressing the momentum balance.
For the implementation of this calculation, see Refs.
5, 6). Fotherdenotes direct interaction between colloii
dal particles and buoyancy.
44

4

3. Rheology of Colloidal Dispersions
Although the rheology of suspensions is importantin terms of practical application, its theoretical generalization is too complex a process. (For example,
refer to Ref. 2) In the low concentration limit of colloidal particles, the theoretical result shows that the apparent viscosity of a suspension ηeff increases as given
in equation (1) (see 2), for example.),ηeff
5=1+ φ,(1)η
2
where η denotes the viscosity of the solvent and φ
the volume fraction of the colloidal particles. This is
Einstein’s well-known viscosity equation, although it
only holds for the range with the volume fraction φ <<0.1. As the particle concentration increases, the effectof particle interactions becomes visible with changes
in particle structures resulting from the interaction,
vitrification and crystallization occurring. This makes
it difficult to predict the rheology of a suspension theoretically on the basis of the basic laws. Although affected by the aforementioned difficulties, the importance of the subject means a number of experimental
studies have been performed, with several empiricalor semi-empirical formulas for the apparent viscosity
having been proposed to summarize such studies.
(See the Doughherty-Krieger equation 3), Refs. 4),
and references therein.)
As the par ticle concentration increases, the involvement of diverse effects other than the colloid
volume fraction in the phenomena necessitates
analysis of an individual phenomenon rather than theformation of a cohesive theory. The Smoothed Profile(SP) method provides a simulation method allowing
solvent-induced many-body interactions between colloidal particles to be calculated directly; it also allows
the rheology of various kinds of suspensions to beanalyzed. KAPSEL so far allows simulation and rheology measurement in a simple shear flow condition tobe performed for a monodisperse system in a Newtonian fluid.

KONA No.24 (2006)

3.1.2 Solvent Motion: Navier-Stokes Equationin Shear Flow
Solvent motion is in the form of incompressible
flow (∇･υ=0). The velocity field υ of the entire system of the solvent and the colloidal particles follows
the Navier-Stokes equation:
1(∂t+υ･∇)υ=- ∇p+ν∇2υ+φfp+f shear,(5)ρρ, η=ρν, and p denote the density, viscosity, and pressure of the solvent, respectively. The body force φfp
takes into account the solid/non-slip boundary condition on surfaces of colloidal particles. For the detailed
implementation, see Ref. 5, 6).
3.2 Rheology Measurement
The external force f shear controls υx so that it is kept
at the profile of a simple shear flow given by
Ly
υx=S -y .
(6)2
In equation (6), S denotes the shear rate. This driving method is the product of the application of theSchumacher and Eckhardt method 9), devised foruniform flow, to the Smoothed Profile (SP) method13.
In this driving method, periodic boundary conditions
are applied in the flow direction (x-) and span direc-

( )

13

To solve a simple shear flow problem, we can use amethod involving giving the velocity/force at the boundar y wall as the boundar y conditions, Rogallo’s moving grid method 8, 9), analogous to the Lees-Edwards
boundary condition in molecular dynamics simulation,
and another method; however, we used the aforementioned calculation method, in which grids do not distortand no boundary problems occur.

169

tion (z-), while a free-slip boundar y condition(∂yυx,z( y=0, L y)=0) is imposed on the upper and lower
walls. The upper and lower walls are non-transparent,
that is, condition υy( y=0, L y)=0 holds. In the SP method, f shear is applied as feedback during each step ofcalculation after the advection of colloidal particles
and fluid is completed.
The apparent viscosity of a suspension is definedby
σ eff
ηeff = xy ,(7)
Swhere S (denoting the shear rate) is a control parameter, and the macroscopic shear stress σ exyf f isevaluated on the basis of the calculation data. Let us
derive a stress calculation method for a suspension.
To facilitate the explanation, the case of a uniform fluid is initially considered. Subsequently, the case of acolloidal dispersion (SP method) is described.
3.2.1 Uniform fluids
Let us write the equations of motion for the shear
flow of uniform fluid as follows:∇･υ=0,(8)ρυ=∇･σ+ρf shear.
(9)
In equation (9), f shear denotes the external force
field per unit mass that exerts shear on the system.
This gives the definition of macroscopic stress in thefollowing form:
1
σeff=
dxσ.
(7)
V
This can be calculated through the volume integralof local stresses σ.
Let us consider evaluating σef f by using external
force ρf shear. Using an identity of second orderσ = [ ∇･(σx)] T - x∇･σ,(11)
the following equations are obtained:4ʃʃ dxσ= ʃ dx[ [∇･(σx)] - x∇･σ ]
= dx[ - x∇･σ] ʃ= dx x ρf - ρυ 
ʃ [(
)]
d
= dxxρf dxρυ.
ʃ
dt ʃ
T

shear

shear

44

(12)
(13)
(14)
(15)

Here, we drop the surface term and use the incompressible condition (∇・υ=0). In the steady state,since the change in the total momentum is zero, we
ultimately obtain1
σ eff=
dxxρf shear.
(16)
V

ʃ

170

3.2.2 SP Method for Particle Dispersion
Let us formally write the equation of the momentum of suspension in shear flow as follows:d(ρtυ)=∇･σ suspension +ρt f shear,
(17)dt
ρt=(1-φ)ρf +φρp.
(18)
In equation (17), σ suspension denotes the stress of thesuspension. Since this σ suspension cannot be calculated
directly, the formula to be obtained through the integration of stress (10) cannot be implemented.
On the other hand, the SP-method-based equation
for the velocity field of a suspension under a shear
flow condition is given by1
υ= ∇･σ +φfp+f shear.
(19)ρ
Comparing this with equation (17), we obtain
1
1∇･σ suspension = ∇･σ +φfp,
(20)
ρt
ρ
From this, we obtain1
σ eff=
dxxρt f shear
(21)
V
as an expression of the stress of the suspension.
Although this contains the position vector explicitly,
the macroscopic shear stress to be measured is freefrom how the origin of the system of coordinates isdefined. This is compatible with the fact that thetranslational motion of the entire system is not driven
by f shear. This fact will be confirmed in the following.
When the origin of coordinates is shifted to a position
expressed by x=x’+x0, we obtain4ʃʃ dx(x’+x )ρ f =ʃ dxx’ρ f + x ʃ dxρ f . (22)
0

t

shear

t

shear

0

t

shear

Since the external force f shear is controlled such as
to produce a shear flow, rather than a translationalmotion of the entire system, ʃ dxρt f shear=0. The SP
method uses equation (21) to evaluate macroscopic
stress.
3.3 UDF Description
3.3.1 Input Parameters
Selecting “Shear_Navier_Stokes” for “constitutive_eq”
makes it possible to simulate the motion of colloidal
particles in a Newtonian fluid under shear flow14. In
the following, the variables to be specified in an input

14

Simulation can be performed similarly when other alternatives for “constitutive_eq” are selected: for example,“Navier_Stokes” for the motion of colloidal particles ina Newtonian fluid, such as gravitational sedimentationof particles, and “Electrolyte” for the motion of charged
colloids, electrophoresis in an external electric field, andsimilar.

KONA No.24 (2006)

UDF file are described.
The section “constitutive_eq.Shear_Navier_Stokes”
specifies the characteristics of the solvent:●constitutive eq.Shear Navier Stokes.DX...Grid width
Δ as the unit quantity for length.
●constitutive eq.Shear Navier Stokes.RHO... Densityof the solvent.
●constitutive eq.Shear Navier Stokes.ETA... Viscosityof the solvent.
●constitutive eq.Shear Navier Stokes.Shear rate...
Shear rate.
The section “particle_spec[],” specifies the characteristics concerning the particle:●particle spec[].Particle number... Number of particles.
●particle spec[].MASS RATIO... Ratio of the densityof the particles to that of the fluid.
The pieces of information on the particle radiusand boundary thickness are entered as follows:●A_XI... Boundary thickness ξ. When “Shear_Navier_Stokes”
is selected, at least ξ ≥ 1 must be selected.
●A… Particle radius.
The section “gravity,” enter the pieces of information on gravity:●gravity.G... Gravitational acceleration.
●gravity.direction... Select the direction in which to
apply gravity.
Specify the type of Lennard-Jones potential for interparticle force.
●EPSILON... the energy unit for the Lennard-Jones
potential.
●LJ powers... the power exponent of the LennardJones potential.
The section “mesh,” specifies the simulation size
(The current version is only compatible with a periodic boundary condition in a rectangular parallelepiped cell.):●mesh.NPX... The size in the x-direction is Lx=2NPX.
●mesh.NPY... The size in the y-direction is L y=2NPY.
●mesh.NPZ... The size in the z-direction is Lz=2NPZ.
The section “time_increment,” specifies the timeincrement:●time_increment... When you select “auto,” specify
the time increment upper limit Tdump=ρ/ηk2max, where
kmax is the wave number determined by the grid
width DX. When you select “manual,” specify thevalue manually.
●time_increment.auto.factor... The time increment is
specified as Δt=Tdump×factor.
●time_increment.manual.delta_t… Enter Δt manually.
The section “switch,” specifies the simulation con-

KONA No.24 (2006)

dition:●switch.ROTATION... Select “ON” to solve an equation of motion for particle rotation.
●s w i t c h . H Y D R O _ i n t . . . S e l e c t “ C o r r e c t ” t osolve a Navier-Stokes equation for a fluid; for“Shear_Navier_Stokes,” however, the “Correct”
switch is the only one implemented.
●s witch.Stokes... Select “with advection” to takethe inertia term into consideration when solving a
Navier-Stokes equation. Select “w/o advection” toperform Stokes approximation not taking the inertia
term into consideration; for “Shear_Navier_Stokes,”
however, the option “w/o advection” is not implemented yet.
●swithch.LJ truncate... When handling force between
particles arising from the Lennard-Jones potential,
select “OFF” for a normal form containing the attractive force term, “ON” for repulsive force, notcontaining the attractive force term, or “NONE”
when no force is applied at all.
●s witch.INIT_distribution... With this switch, select the initial configuration of par ticles from“unifor m_random” (random),“random_walk”(shifted from a square lattice randomly), “FCC”(on an FCC lattice), “BCC” (on a BCC lattice), and“user_specify” (the coordinates and velocity areuser-specified.). When you select “user_specify,”
enter the initial particle position and velocity into
“user_specify.Par ticles[].R” and “user_specify.
Particles[].v,” respectively. If the number of lists
to be entered is smaller than that specified in“Par ticle_number,” increase “user_specify.Particles[]” by implementing the “Edit → Add an array
Element” operation on “Gourmet” or edit the UDF
file directly.
●switch.FIX_CELL... Set the direct current component of the full velocity at zero, υk=0=0. To introduce
gravitational acceleration to measure the settling
velocity, set this to “ON.”The section “output,” specifies the pieces of information on the data output:●output.GTS... Number of data output interval steps.
●o utput.Num_snap...Frequency of data output; inother words, the total number of steps is specified
by GTS×Num_snap.
●output.AVS... Select “ON” to carry out data output
in the AVS format.
●output.AVS.ON.Out_dir... Insert the director y for
outputting AVS format data. For example, when“data” is specified, it is necessar y to create the“./data” and “./data/avs/” directories in advance.
The AVS field file is outputted to the “./data/” direc-

171

tory under the filename of “data.fld”; the data file isoutputted to the “./data/avs/” directory under the
filename of “data_*.dat.” The number of steps is entered in the asterisk *.
●output.AVS.ON.FileType... Select between “Binary”
and “ASCII” for formatting the AVS data file.
●output.UDF... Select “ON” for output UDF.
As many pieces of record data of the coordinatesand velocities of particles as specified in “output.
Num_snap” are stored after “Particles[]” of the output UDF.
In the section “resume,” select from the following
options to resume an interrupted calculation for the
resume UDF:
●resume.caclucation... Select “NEW” to start a new
calculation. Select “CONTINUE” to restart calculation after reading the pieces of information concerning the number of calculation steps, particle coordinates, velocity, angular velocity, solvent velocity
field, and ion density field obtained at the end of theprevious calculation. The data is saved in a location
below “resume.CONTINUE.Saved_Data” of the restart UDF file. In this case, calculation is resumed
with the restart UDF created at the end of the previous calculation being used as the input UDF. How
to use this input UDF will be described in detail in
Section 5.5.
3.3.2 Units of Space and TimeThe grid width Δ is used as the unit of length.
The following points should be noted concerning
the unit of time:●A system of units is adopted that makes the relation
ρ=η=Δ=1 hold in the Navier-Stokes equation.
●When the entries RHO=A, ETA=B, and DX=C are
made in the input udf file, the maximum wave number kmax is given using C, with the upper limit of thetime increment given as Tdump = (A/B)/k2max on thebasis of the diffusion time of the momentum. The
time increment Δt is adjusted by Tdump×factor. However, the unit for time itself is the same as in the caseof ρ=η=Δ=1.
4. Electrophoresis of Charged Colloids
When colloidal particles are dispersed in a solventwith very large permittivity such as water, ions are
emitted from leaving groups on colloid surfaces, withthe particle surfaces becoming charged. The emitted
ions are electrostatically attracted to particle surfaces, and at the same time, are diffused by thermal fluctuation and form a cloud-like ionic atmosphere called172

an electric double layer around the colloidal particles.
The properties of a colloidal disperse system in equilibrium are described by a Poisson-Boltzmann equation. The Debye-Hückel approximation, obtained bylinearizing the Poisson-Boltzmann equation, allows
an unobstructed understanding to be obtained. In
contrast, in so-called electrokinetic phenomena, suchas electrophoresis, the behavior of particles and ionic
distribution is determined by competition betweenhydrodynamic and electrostatic interactions. As aresult of the competition, the ionic distribution becomes incapable of following the motion of particles,with electric double layers distorted out of sphericalsymmetry, resulting in a non-equilibrium state. The
only means to analyze such a complex state is tointroduce a theoretically simple approximation. In addition, even computer-based simulation has seldom
yielded a correct reproduction.
To overcome this situation, this project involved
developing a new simulation method called the
Smoothed Profile (SP) method 5, 10, 6, and 11).
In the SP method, colloids are handled as particle
pictures. In contrast, the solvent and ions are presented as a coarse-grained density field in the stateof a continuum. In a continuum picture, the interface
between a colloidal particle and a solvent, usually expressed through a step function, is expressed safely
on a fixed orthogonal grid, using an interface function to overcome difficulties in numerical calculation.
Accordingly, the SP method uses both particle andcontinuum pictures to describe their degrees of freedom and the interaction between them with physicalcorrectness. This is an outstanding feature of the SP
method that allows multi-scaled complex phenomena
to be analyzed. In addition, the SP method solveshydrodynamic interactions by carr ying out direct
numerical calculations of a Navier-Stokes equation 5,
6). To implement these calculations, we constructed
a formalism to consistently solve the three degrees
of freedom, namely, colloidal particles (the Newton’s
equation of motion), ion distribution (an advectiondiffusion equation), and a solvent velocity field (theNavier-Stokes equation), based on the assumptionof a smooth interface function. Using this, we firstly
succeeded in simulating electrophoresis quantitatively 11). In the following section, the basic equations
taken into consideration in KAPSEL are described,
with a brief mention of the underlying theory.
4.1 Basic Equations
As described in the Refs. 6, 11), let us take a general view of the basic electrohydrodynamics equa-

KONA No.24 (2006)

tions necessary for simulation. N colloidal particles
with a radius of a, dispersed in an electrolytic solvent,are considered. The permittivity of the solvent � isassumed to be spatially uniform. The SP method expresses the interface between a particle and solventin the form of a smooth function φ(r)∈[0,1] of finite
thickness ξ. On a fixed orthogonal grid, the particle
domain is expressed by φ=1, the solvent domain by
φ=0, and the interface domain by 0 <φ< 1. Other
methods using φ(r) on an orthogonal grid are proposed by Tanaka-Araki 12) and Kajishima et al. 13).
Introducing an interface function allows remarkably
improved calculation efficiency in comparison withthe finite element method using nonstructural lattices. The surface of a colloidal particle is assumed to
be charged uniformly and the quantity of charge perparticle to be Ze. In an ordinary continuum picture,
the particle surface charge distribution is expressedusing a delta function. Consequently, a suitableboundary-compatible lattice is used in the finite element method; which significantly encumbers thecalculation efficiency. In contrast, the SP method alsouses a smooth distribution for this particle surface
charge distribution eq(r). For example, when we use
the differential of the first order of the interface functionφ, eq(r) is expressed as
Ze|∇φ(r)|eq(r)=.
(23)
4πa2
Just as the interface function φ(r) is structured
such as to be reduced to a step function with ξ→0,
q (r) is structured such as to be reduced to a delta
function with ξ→0.
4.1.1 Advection-diffusion Equation
Let us define the density distribution Cα of α typeions with a valence of Zα in the entire calculation domain as shown below:
Cα(r, t)=(1−φ(r, t))C＊（24）
α(r, t).
The domain where ions exist is expressed by (1φ). C＊
α (r, t) is the auxiliary variable for calculation,
being defined as smooth throughout the entire calculation domain. C＊
α in the colloid domain (φ=1) has
no physical meaning. The total charge distribution,
including the particle surface charge distribution, isgiven byΣαρe(r)=e

＊
＊
∂tC＊
μα).
(26)
α =-∇･Cα υ+Γα∇･(Cα ∇This equation consists of two terms, one expressing the advection, due to the solvent velocity field υ,and the other expressing the diffusion due to the gradient of the chemical potential μα. Since C＊
α follows＊
the advection-diffusion equation, ʃ drCα is conservative.
The absence of ion permeation into colloidal particles is considered to be the condition for the component normal to the interface of the ion diffusion flowrate to become zero. In other words, it is expressedby n ･ ∇μα =0 6, 11). In this equation, n denotes the(outward) normal vector on the particle surface. Using the interface function, n can be expressed in theform of n=- ∇ φ/ | ∇ φ | . Γα is the Onsager transportcoefficient of an α-type ion, being related to the friction and diffusion coefficients of the ion as fα=1/Γα,
Dα=kBT Γα. The chemical potential of the ion is assumed to be expressed in the form 16):μα=kBT lnC＊(27)
α +Zαe(Ψ-E･r).
In the equation, E denotes the external electricfield, and the electrostatic potential Ψ(r) is obtainedby solving the Poisson equation2�∇ Ψ=-ρe.
(28)This chemical potential assume a Poisson-Boltzmann
distribution in equilibrium.

4.1.2 Navier-Stokes Equation
Solvent flow is incompressible ( ∇ ･ υ=0). The velocity field υ of a solvent follows the Navier-Stokes
equationρ(∂t+υ･∇)υ=-∇p+η∇2υ-ρe(∇Ψ-E)+φfp. (29)In the equation, ρ, η, and p denote the density,
the viscosity of a solvent, and pressure, respectively.
Note here that electrostatic force -ρe(∇Ψ-E) acts on
the fluid. The external force term φfp ensures the rigidity of a particle. In other words, φfp takes the nonslip boundary conditions of the particle surface into
consideration. The details are defined in Refs. 5, 6).
4.1.3 Equation of Motion of Colloidal Particles
The position, velocity, and angular velocity of theith colloidal particle having mass Mp evolve over time
according to the equations of motion
Ri=Vi,
(30)
H
other
H
MpVi=Fi +Fi , Ip･Ωi=N i .
(31)
In equation (31), FHi denotes the force actingon theparticle induced by the fluid, expressing the momentum balance between the solid and fluid 6); F otherdeinotes the force due to the potential between particles,
such as the Lennard-Jones potential. The rotationalmotion of the particles is also taken into consider44ZαCα(r)+eq(r).

(25)

In the equation, the initial distribution is constructed such as to meet the electrical neutrality conditionʃ ρedr=0.
The time evolution of the auxiliary ion density C＊
α isassumed to follow the advection-diffusion equation

KONA No.24 (2006)4173

ation, though not described here 6). The above are
the equations necessary for the simulation.
4.2 Electric Double Layer
4.2.1 Poisson-Boltzmann Equation
The basis for quantitative treatment of the structure of an electric double layer is the PoissonBoltzmann equation. Let us determine the equilibrium distribution of ions under equation (27) for thecase where the external electric field does not exist(E=0). When the chemical potential has become uniform, that is, μα=constant, we obtain_ZαeΨ(r)C*α(r)=C αexp (32)kBT
as the equilibrium ion distribution. This is a
Boltzmann distribution under the electrostatic potential Ψ. Combined this with equation (28), called the
Poisson-Boltzmann equation.

()

4.2.2 Debye-Hückel Approximation and Debye
Screening Length
Let us consider a spherical colloidal particle ina z:z-symmetr y electrolytic solvent. The PoissonBoltzmann equation
for this condition is given by 2)_2zeCzeΨ(r)∇2Ψ(r)=
sinh
.
(33)
�
kBT
The boundary conditions
at a point at infinity are_Ψ|r=∞ =0 and C ＊ |r=∞ =C . The boundary condition on
the surface of the particle is given byσe∇Ψ|surface=- ,
(34)

(

)

�

where the sur face charge density is given by
σe=Ze/4πa2. In other words, a boundar y conditionof constant surface charge is considered. A methodof approximation, in which ze Ψ /kBT ≪1 is assumed
and equation (33) is linearized, is called the DebyeHückel approximation.
It gives_
2z2e2C2∇ Ψ(r)=
Ψ=κ2Ψ.
(35)
kBT�
In connection with this equation, the constant having a length dimension
1 _
κ-1=
(36)
√￣￣￣̅
8πλBz2C
is called the Debye screening length. In equation(36), λB=e2/4πkBTε is called the Bjerrum length. In an
ordinary electrolyte,
1
_
κ-1=
(37)
√￣￣￣￣￣̅
4πλBΣαZ α2 C αholds. Since the system is spherically symmetric, it
is sufficient to consider the radius r=|r| only, then we
obtain174d2Ψ 2 dΨ
+
=κ2Ψ.
(38)
r dr
dr2
The general solution of this equation is given by a
Yukawa-type potentialaΨ(r)=Ψ0 exp[-κ(r -a)].
(39)r
The electrostatic force arising from colloidal
charge is screened to the extent of κ-1. the Debye
screening length κ-1 can be considered the distanceat which a balance is struck between the effect of the
Coulomb force, via which a surface charge attracts
an ion of the opposite sign, and the thermal diffusion
tending to scratch out the local presence of an ion.
For this reason, κ-1 can be regarded as the thicknessof an electric double layer. When the temperature
is high and hence the thermal energy kBT is large,
κ-1 becomeslarge. In addition, as the ionic strength_2(ΣαZ αC α/2) increases, the screening effect increases
and κ-1 decreases. For example,
let us consider the_bulk salt concentration C throughout a z:z-symmetricelectrolytic solvent. With the medium assumed to bewater with a temperature of 25℃, the Bjerrum length
is given by λB=0.72nm; substituting the actual value
into equation (36), we obtain
0.3
κ-1= _ (nm).
(40)
z √̅
C
_
-1
When
_ z=1 and C =0.1M, κ =1nm; and when z=1-1
and C = 0.001M, κ =10nm.
The surface potential Ψ(r=a)=Ψ0 is determined on
the basis of the boundar y condition dΨ/dr(r=a)=-σe/�. Ψ and the surface charge density σe are related byσe=�κΨ0(1+(κa)-1).

(41)

In other words, as the surface charge increases,
the surface potential increases linearly. The DebyeHückel approximation cannot be applied to locations
where the surface potential is high. In reality, as thesurface charge increases, the increase in the surface
potential becomes more difficult. The relation for thissituation can be obtained by solving the nonlinear
Poisson-Boltzmann equation. Approximate solutions
for a 1:1 electrolyte solvent are proposed by LoebOverbeaak-Wiersema 2) and Oshima-Healy-White
15).
4.3 Backgrounds of Electrophoresis
When an electric field E is applied from outside,a colloidal particle charged to Ze is subjected to theelectrostatic force ZeE, and starts moving. Although
the particle is accelerated by the electrostatic force,
it is subject to viscous resistance from the fluid;

KONA No.24 (2006)

and reaches the steady velocity V when a balance isstruck between the electrostatic force and the viscous resistance. Assuming the viscous resistance to
be Stokes resistance, 6πηaV, this allows the balance
between forces to be written in the form of
ZeE=6πηaV.

(42)

From this, the electrophoresis mobility
Ze
V
=
(43)
E 6πηais derived. However, this is unsuitable, and theactual electrophoresis mobility is smaller than thisvalue. This is due to the lack of consideration of the
motion of the ion atmosphere and electrostatic force,
despite the electrostatic force acting on the ion atmosphere around the colloid. Accompanied by the
motion of the ion atmosphere, colloid migration exhibits a much slower migration velocity. In addition,
if the electric double layer distorts out of sphericalsymmetr y, the resulting force must be taken into
consideration. This will make theoretical consideration extremely complicated. For this reason, several
analyses of electrophoresis have been made on simplified models 2).
4.3.1 Smoluchowski Theor y
In the case of κa≫1, where the radius of a particle
a considerably exceeds the thickness of an electric
double layer κ-1, Smoluchowski’s equation is applicable. Since the thickness of an electric double layer
is infinitely small in this case, the particle surface
can be regarded as a plane with the curvature of theparticle surface neglected. An external electric fieldEx, parallel to the plane (assumed to be in the x-direction), is assumed to be applicable. When the motionof the fluid is observed from above the particle, in acoordinate system fixed to the latter, the fluid, apartfrom the particle, should have a velocity of−V. Considering the balance between the viscous resistanceand the electrostatic force, we obtain∂2υxη=eCαE x=0.
(44)
∂y2
α
Since the ion distribution is related to the second
order dif ferentiation of the electrostatic potential
through the Poisson equation, the following relation
holds:
∂2υx
∂ 2Ψ
η
E x.
(45)
2 =�
∂y
∂y2
Performing integration under the boundary conditions whereby the velocity gradient, potential gradient, and potential are all equal to zero for y→∞, thisyieldsη[υx(y)+V]-�EΨ(y)=0.

Subsequently, we obtain Smoluchowski’s equation
V �ζ
＝
(47)
E η
because the velocity is zero at the particle surface(y = 0). In equation (47), the zeta potential ζ is originally defined as the potential at the slip plane, but it
is replaced by the potential at the particle surface Ψ(0).
4.3.2 Hückel Theor y
Let us consider a limit opposite to the condition
represented by Smoluchowski, κa ≪1, at which theparticle is extremely small compared with the thickness of the electric double layer. In other words,
when Ze is a limit of a point charge, Hückel’s equation is applicable. When the potential on the particle
surface is replaced by the Coulomb potential
Ze
ζ=
,
(48)
4π� a
Hückel’s equation
V 2 �ζ
=(49)
E 3 η
is obtained.
4.3.3 Henr y and O’Brien-White Theories
For an ordinary κa value, Henry’s equation�ζV=f(κa)(50)ηE
is known as a coupler of equations (47) and (49).
f(κa) is called Henry’s coefficient, being defined asfollows:

1

Σ

KONA No.24 (2006)(46)‘henry.dat’

0.95

f (κa)

0.9
0.85
0.8
0.75
0.7
0.1

1

10

100

1000

κa
Fig. 1 Henry’s coefficient f(κa).

175

f(κa)=1-5exp(κa)E7(κa)+2exp(κa)E5(κa)
2 (κa)2 5(κa)3 (κa)4 (κa)5
= +
+
24
72
144
144
3
4
6
(κa)
(κa)
+
exp(κa)E1(κa).
12
144

[

]

(51)

(52)

In equation (51), En (κa) is an exponential-integral
function of the nth order. Smoluchowski’s equation
corresponds to f=1(κa→∞) and Hückel’s equation tof = 2/3 (κa→0). (Fig. 1)
Henr y’s equation (50) is a linear relation to thezeta potential, applicable only when the zeta potential
is low. As the zeta potential increases, considerationof the deformation effect of the electric double layer(called the relaxation effect) becomes necessary. Using numerical analysis, O’Brien and White have proposed, for ordinary κa and ζ values, the relationshipbetween the electrophoretic mobility and the zeta potential 16). Subsequently, Ohshima, Healy, and White
have proposed an analytical equation for the range ofκa ≥ 10 17).
4.4 UDF Description
4.4.1 Input Parameters
Selecting “Electrolyte” as “constitutive_eq” allows
the simulation of motion of charges colloids and anelectrophoretic phenomenon under an external electric field15. In the following, the variables to be specified in an input UDF file are described.
Enter pieces of information on the solvent and ion
distribution in the section “constitutive_eq.Electrolyte”:●constitutive_eq.Electrolyte.DX... Grid width Δ as
the unit quantity for length●constitutive_eq.Electrolyte.RHO... Solvent density.
●constitutive_eq.Electrolyte.ETA... Solvent viscosity.
●constitutive_eq.Electrolyte.kBT... Thermal energy
at the temperature T.
●constitutive_eq.Electrolyte.Dielectric cst... Solvent
permittivity.
●constitutive_eq.Electrolyte.Init_profile... Choose between “Uniform” and “Poisson_Boltzmann.” When
you choose “Boltzmann,” note that considerable
15

Simulations can be performed also when alternatives tothe above parameters are selected for “constitutive_eq”: for example, “Navier_Stokes” for the motion of colloidal particles in a Newtonian fluid, such as gravitational
sedimentation of particles and “Shear_Navier_Stokes”
for the motion of colloidal particles in a Newtonian fluid
under a shear flow condition, rheology of a dispersion
system, and similar.

176

time may be needed in induction periods of simulations to handle a multi-particle problem, in particular.
●constitutive_eq.Electrolyte.Add_salt... When you
use “saltfree,” consider a single component of ion
having the sign opposite to that of the particle surface. If you use “salt,” consider two types of ions,
positive and negative.
●c onstitutive_eq.Electrolyte.Add_salt.saltfree.
Valency_counterion... Valency of a counterion.
●constitutive_eq.Electrolyte.Add_salt.saltfree.Onsag
er_coeff_counterion... Onsager transport coefficientof a counterion.
●constitutive_eq.Electrolyte.Add_salt.salt.Valency_
positive_ion... Valency of a positive ion.
●constitutive_eq.Electrolyte.Add_salt.salt.Valency_
negative ion... Valency of a negative ion.
●constitutive_eq.Electrolyte.Add_salt.salt.Onsager_c
oeff_positive_ion... Onsager transport coefficient ofa positive ion.
●constitutive_eq.Electrolyte.Add_salt.salt.Onsager_c
oeff_negative_ion... Onsager transport coefficient ofa negative ion.
●c onstitutive_eq.Electrolyte.Add_salt.salt.Debye_
length... Specifying the Debye screening length
causes the corresponding salt concentration to be
specified.
●constitutive_eq.Electrolyte.Electric_field... Choose
between “ON” and “OFF” of the external electricfield.
●c onstitutive_eq.Electrolyte.Electric_field.ON...
Choose between DC (direct current electric field)
and AC (alternating current electric field).
●constitutive_eq.Electrolyte.Electric_field.ON.DC.
Ex... Intensity of the x-direction electric field; thesame applies in the case of an alternating current
electric field.
●constitutive_eq.Electrolyte.Electric_field.ON.DC.
Ey... Intensity of the y-direction electric field; thesame applies in the case of an alternating current
electric field.
●constitutive_eq.Electrolyte.Electric_field.ON.DC.
Ez... Intensity of the z-direction electric field; thesame applies in the case of an alternating current
electric field.
●constitutive_eq.Electrolyte.Electric_field.ON.AC.
Frequency... Frequency of an alternating current
electric field.
The section “particle_spec[],” specifies pieces of
information on the particle type:●particle_spec[].Particle_number... Number of particles.

KONA No.24 (2006)

●particle_spec[].MASS_RATIO... Specific weight of aparticle relative to the fluid.
●particle_spec[].Surface_charge... Total amount ofparticle surface charges.
Enter the pieces of information on the par ticle
radius and interface thickness into the following locations:●A_XI... Interface thickness ξ. When “Electrolyte”
is selected, the differential of the first order of theinterface function ∇φ is used; consequently, thethickness of the smallest interface is ξ=2.
●A... Particle radius.
The section “gravity,” specifies the pieces of information on gravity:●gravity.G... Gravitational acceleration.
●g ravity.direction... Select the direction in which
gravity is applied.
Determine the type of the Lennard-Jones potential
for interparticle force.
●“EPSILON”... Specify the unit for the energy of the
Lennard-Jones potential.
●“LJ_powers”... Specify the power exponent of the
Lennard-Jones potential.
The section “mesh,” specifies the simulation size
(The current version is compatible only with a periodic boundary condition in a rectangular parallelepiped cell.):●mesh.NPX... The size in the x-direction is Lx=2NPX.
●mesh.NPY... The size in the y-direction is Ly=2NPY.
●mesh.NPZ... The size in the z-direction is Lz=2NPZ.
The section “time_increment,” specifies the timeincrement:●time_increment... When you select “auto,” specify
the time increment upper limit Tdump=min(ρ/ηk2max,
1/kBT Γα k2max), where kmax is the maximum wave
number determined by the grid width DX. Whenyou select “manual,” specify the value manually.
●time_increment.auto.factor… The time incrementis specified as Δt=Tdump ×factor using “factor,” thevalue to be entered.
●t ime_increment.manual.delta_t… Enter the Δtvalue manually.
The section “switch,” specifies the simulation conditions:●switch.ROTATION... Select “ON” to solve an equation of motion for particle rotation.
●s witch.HYDRO int... Select “Correct” to solve a
Navier-Stokes equation for a fluid; for “Electrolyte,”
however, the “Correct” switch is the only one imple-

KONA No.24 (2006)

mented.
●switch.Stokes... Select “with advection” to take theinertia term into consideration when solving a Navier-Stokes equation. Select “w/o advection” to perform a Stokes approximation not taking the inertia
term into consideration; for “Shear_Navier_Stokes,”
however, the option “w/o advection” is not implementedyet.
●s withch.LJ_truncate... When handling force between par ticles arising from the Lennard-Jones
potential, select “OFF” for a normal form containing
the attractive force term, “ON” for repulsive force
not containing the attractive force term, or “NONE”
when no force is applied at all.
●s witch.INIT_distribution... With this switch, select the initial configuration of par ticles from“uniform_random” (random), “random_walk”(shifted from a square lattice randomly), “FCC”(on an FCC lattice), “BCC” (on a BCC lattice), and“user_specify” (the coordinates and velocity areuser-specified.). When you select “user_specify,”
enter the initial particle position and velocity into
“user_specify.Par ticles[].R” and “user_specify.
Particles[].v,” respectively. If the number of lists
to be entered is smaller than that specified in“Par ticle_number,” increase “user_specify.Particles[]” by carrying out the “Edit → Add an array
Element” operation on “Gourmet” or edit the UDF
file directly.
●switch.FIX_CELL... Set the direct current component of the full velocity at zero, that is, υk=0=0. To
introduce gravitational acceleration to measure thesettling velocity, set this to “ON.”The section “output,” specifies the pieces of information on the data output:●output.GTS... Number of data output interval steps●output.Num snap... Number of times of data output;
in other words, the number of total steps is specified by “GTS×Num_snap.”●output.AVS... Select “ON” to carry out data output
in the AVS format.
●output.AVS.ON.Out_dir... Insert the director y for
outputting AVS format data. For example, when“data” is specified, the “./data” and “./data/avs/” directories must be created in advance. The AVS field
file is outputted to the “./data/” directory under the
filename of “data.fld”; and the data file is outputtedto the “./data/avs/” under the filename of “data_*.dat.” The number of steps is entered in the asterisk*.
●output.AVS.ON.FileType... Select between “Binary”

177and “ASCII” for formatting the AVS data file.
●output.UDF... Select “ON” for output UDF.
The number of pieces of record data of coordinates and velocities of particles specified in “output.
Num_snap” is stored after “Particle[]” of the output
UDF.
In the section “resume,” select from the following
options to resume an interrupted calculation for the
resume UDF:
●resume.caclucation... Select “NEW” to start a new
calculation. Select “CONTINUE” to restart calculation after reading the pieces of information on thenumber of calculation steps, particle coordinates,
velocity, angular velocity, solvent velocity field, andion density field obtained at the end of the previous
calculation. The data is saved in a location below“resume.CONTINUE.Saved_Data” of the restart
UDF file. In this case, calculation is resumed withthe restart UDF created at the end of the previous
calculation being used as the input UDF. How to use
this input UDF will be described in detail in Section
5.5.
4.4.2 Units of Space and TimeThe grid width Δ is used as the unit of length.
The following points should be noted for the unit oftime:●A system of units is adopted that makes the relation
ρ=η=Δ=1 hold in the Navier-Stokes equation.
●If the entries RHO=A, ETA=B, and DX=C are madein the input UDF file, the maximum wave number
kmax is given using C, with the upper limit of the timeincrement given as Tdump=(A/B)/k2max on the basis of
the diffusion time of the momentum. The time increment Δt is adjusted by Tdump × factor. However,the unit for time itself is the same as in the case of
ρ=η=Δ=1.
●In the case of “Electrolyte,” use either (A/B)/k2max or(1/kBTΓα)/k2max as Tdump, whichever is smaller. Use aunit for time so that ρ=η=Δ=1 holds.
5 Tutorial
5.1 How to Run KAPSELKAPSEL is executed as follows:% kapsel -Iinput.udf -Ooutput.udf -Ddefine.udf-Rrestart.udf
KAPSEL is executed by using -I to designate theinput UDF, -O to designate the output UDF, -D to

178

Fig. 2

Flow of a suspension under shear; arrows indicate the velocity
field on a plane normal to the flow.

designate the definition UDF, and -D to designate therestart UDF. Note that no option can be omitted. The
restart UDF will be described in Section 5.5.
5.2 Sample Simulations
5.2.1 Particle Dispersion under Shear“shear.udf” can be used as a sample input UDF file;
to be more specific, initially ensure that “./data” and“./data/avs/” have been created, and then executethe following:% kapsel -Ishear.udf -Ooutput.udf -Ddefine.udf-Rrestart.udf
In this example (shown in Fig. 2), calculation iscarried out on the assumption that the mesh size is
64×64×64, the shear rate S=0.01, the particle radius
a=4, the interface thickness ξ=1, and the number ofparticles is 288. Under this condition, the particle volume fraction φ is 0.29.
With “Shear_Navier_Stokes” selected as“constitutive_eq,” the data as shown below is outputted in the standard error output:#1:dev_sigma_yx 2:shear_rate_obs 3:shear_rate ...
0.00647023 -0.00449194 0.01 ...
.
.
.
The data in the second and the eighth columns
show the temporal shear rate and the temporal shear
stress, respectively. Since the data in the 17th column
represents the shear strain multiplied by L y, one may
draw the stress-strain curve by plotting “-$8” versus“-$17/L y,” where “$n” denotes the data in the nthcolumn. The data in each column is accumulated over

KONA No.24 (2006)

GraphSheet[]0“plot.dat” using 2:6

–0.02

vx0

–0.04
–0.06
–0.08
–0.1
–0.12
–0.14

0

2000

4000

6000

8000

10000

time step
Fig. 4 Evolution of the migration velocity of a colloidal particle overtime.
Fig. 3A colloidal particle (colored blue) in electrophoretic motion under an external electric field applied in the +x direction; in the
figure, light and shade represent the electric charge density distribution and the arrows the solvent velocity fields.

simulation time excluding the initial one fourth of thetotal duration, and the averaged value is printed in
the last line.
5.2.2 Electrophoresis: a single particle
To obtain an example of calculation for the simulation of single-particle electrophoresis, use “sample1.
udf ” as a sample input UDF file; to be more specific,
first make sure that “./data” and “./data/avs/” havebeen created, and then execute the following:% kapsel -Isample1.udf -Ooutput.udf -Ddefine.udf-Rrestart.udf
In this example, a 1:1 electrolytic solvent is considered in a mesh size of 64×64×64, and the calculation is carried out on the assumption that the radius
a=5, the interface thickness ξ=2, the colloidal charge
Z=-100, the electric field E x=0.1, and the Debye
screening length κ-1=10. Fig. 4 shows the evolutionof the migration velocity over time. The figure showsthat a steady velocity is realized at the point of equilibrium between the electrostatic force and the fluid
resistance.

Fig. 5

A colloidal particle (colored blue) in electrophoretic motion under an external electric field applied in the +x direction.

Fig. 6

Positively charged colloidal particles (colored red) and negatively charged colloidal particles (colored blue) in electrophoretic motion under an external electric field applied in the +x
direction; positive-signed colloids migrate in the +x direction
and negative-signed ones in the -x direction.

5.2.3 Electrophoresis: many particles
The simulation of electrophoresis with the numberof particles assumed to be 32 can be calculated with“sample2.udf.” Other parameters are identical tothose in the previons case where the number of particles is one.
5.2.4 Electrophoresis: +/- binar y particles
The electrophoresis of a mixed system, in which
64 positively charged colloidal particles and as many

KONA No.24 (2006)

179

negatively charged particles exist, can be calculatedwith “sample3.udf.” Other parameters are identical tothose in the previous cases.
5.3 Visualization
5.3.1 Using AVS/Express
With “output.AVS” set at “ON,” an AVS-format
file is outputted. Using AVS/Express16 to read the“avs_charge.v” file and designating the field file “data.
fld” as the “Read_field” module allows the particles,
solvent velocity field, and charge distribution to bevisualized respectively. The particles are displayed incolors according to the surface charge value Ze. The
visualization shown in Figs. 3, 5, and 6 is obtainedby using AVS/Express with “avs_charge.v.”
5.3.2 Using “Gourmet”
When one starts up “Gourmet” to read “output.
udf,” loads a Python script “show_field.py” suppliedwith KAPSEL on the Python panel of “Gourmet,” andexecutes “Run,” the graphic window opens to visualize the particles, solvent velocity field, and charge
density distribution17. And clicking on the playback
button at the bottom of the graphic window causes
animation to start. Since particles are visualized on
the basis of the particle coordinate data outputtedto the output UDF as record data, setting “output.
UDF” at “ON” is necessary. Since the solvent velocity
field and the charge density distribution have read
“Binary” files outputted in the AVS format, it is necessar y to set “output.AVS” at “ON” and “output.AVS.
ON.FileType” at “Binary.”
5.4 Plotting Data: gnuplot
Using “gnuplot18,” it is possible to plot the time series of the position and velocity of a particle in graphform on the basis of the record data outputted to
UDF 1).
When one starts up “Gourmet” to read “output.
udf,” loads a Python script “plot.py” supplied withKAPSEL on the Python panel of “Gourmet,” andexecutes “Run,” one can create a graph sheet on thebasis of the record data storing the positions andvelocities of particles. To plot those data, change the
16
17

18

http://www.avs.com/
Versions earlier than “Gourmet_2003” have some
glitches and it has been shown that this script does notwork normally. In such cases, use the supplied Python
script “particleshow.py” for visualization. The solvent
velocity field and the charge density distribution are notdisplayed; only particles are displayed.
http://www.gnuplot.info/

180

check mark placed on “Tree” in the “View” box to“Table,” select “GraphSheet” in the list of variables,and browse the “GraphSheet” data. Subsequently,
executing the “Make” operation on the “Plot” panelof “Gourmet” and executing the “Plot” operation allows the graph shown in Fig. 4 to be drawn. If anexcessive number of line types hinders an easy viewof a graph, one may delete unnecessar y por tions
from the “plot” command created in the editor of the“Plot” panel. How to call “gnuplot” from “Gourmet” isdescribed in detail in Chapter 3 of the document 1).
5.5 Restart UDF
Select “NEW” in “resume.Calculation” in the input
UDF and end the calculation in the sequence% kapsel -Iinput.udf -Ooutput.udf -Ddefine.udf-Rrestart.udf
and the pieces of information concerning the number of steps, particles, solvent velocity field, and ion
distribution at the end of the calculation will be outputted to “restart.udf.” To restart calculation using
these pieces of data, follow the steps below:
1. Copy “restart.udf” to “input.udf.”
2. Open “input.udf” on “Gourmet” and select “CONTINUE” in “resume.Calculation.”
3. Specify and store the number of steps to be calculated as a single unit after increasing “output.
Num_step.”
4. If the sequence% kapsel -Iinput.udf -Ooutput.udf -Ddefine.udf-Rrestart.udf
is followed again, the data is read from “input.udf”
to restart calculation.
5. For both the AVS data and “output.udf,” the output
data of the resumed calculation can be added withthe previously calculated data as it is.
6. The data at the completion of calculation is stored
in “restart.udf” again; therefore, repeating the steps
above allows calculation to be continued.
5.6 Known Problems
1. When “Electrolyte” is selected for “constitutive_eq,”
a segmentation error or abortion may occur at thetime of ending a program run (at the time of manipulating the restart UDF file) depending on themachine used. However, the result of the execution
itself seems to pose no problems.
2. “show_field.py” seems not to work normally on“Gourmet_2003.”
3. KAPSEL that has been built in the Windows Cygwin + MinGW environment seems to end due to errors during execution. When “MinGW” is not used,

KONA No.24 (2006)

no problems occur.
References
１) Doi, M. and Y. Nishio: “Graphical User Inter face
GOURMET Operation’s Manual”, OCTA Users Group,Tokyo, Japan (2005), http://octa.jp.
２) Russel, W.B., D.A. Saville, and W.R. Schowalter: “Colloidal Dispersions”, Cambridge University Press,
Cambridge, England, (1989).
３) Krieger, I.M.:. Rheology of monodisperse latices, Adv.
Colloid Interface Sci., 3. 111 (1972).
４) Hsueh, C.H. and P.F. Becher: Effective viscosity ofsuspensions of spheres, J. Am. Ceram. Soc., 88. 1046
(2005).
５) Nakayama Y. and R. Yamamoto: Simulation method
to resolve hydrodynamic interactions in colloidal dispersions, Phys. Rev. E, 71. 036707 (2005).
６) Nakayama, Y., K. Kim, and R. Yamamoto: Simulating(electro)hydrodynamic effects in colloidal dispersions:smoothed profile method, http://arxiv.org/abs/condmat/0601322 (2006).
７) Schumacher J. and B. Eckhardt: On statistically stationar y homogeneous shear turbulence, Europhys.
Lett., 52. 627 (2000).
８) Rogallo, R.S.: Numerical experiments in homogeneous
turbulence: NASA Tech. Memo., 81315 (1981).
９) Onuki, A: A new computer method of solving dynamic
equations under externally applied deformations, J.
Phys. Soc. Jpn., 66. 1836 (1997).

KONA No.24 (2006)

10) Kim, K. and R. Yamamoto: Ef fcient simulations ofcharged colloidal dispersions: A density functional appraoch, Macromol. Theory Simul., 14. 278 (2005).
11) Kim, K., Y. Nakayama, and R. Yamamoto: Direct Numerical Simulations of Electrophoresis of Charged
Colloids, Phys. Rev. Lett., 96. 208302 (2006).
12) Tanaka H. and T. Araki: Simulation method of colloidal
suspensions with hydrodynamic interactions: Fluid
particle dynamics, Phys. Rev. Lett., 85. 1338 (2000).
13) Kajishima, T., S. Takiguchi, H. Hamasaki, and Y. Miyake: Turbulence structure of particle-laden flow in a
vertical plane channel due to vortex shedding, JSMEInt. J. Ser. B, 44. 526 (2001).
14) Barrat, J.-L. and J.-P. Hansen: “Basic Concepts for
Simple and Complex Liquids”, Cambridge University
Press, Cambridge, England (2003).
15) Ohshima, H., T. W. Healy, and L. R. White: Accurate
analytic expressions for the surface charge density/
surface potential relationship and double-layer potential distribution for a spherical colloidal particle, J. Colloid Interface Sci., 90. 17 (1982).
16) O’Brien, R. W. and L. R. White: Electrophoretic mobility of a single colloidal particle, J. Chem. Soc. Faraday
Trans. 2, 74. 1607 (1978).
17) Ohshima, H., T. W. Healy, and L. R. White: Approximation analytic expressions for the electrophoretic mobility of spherical colloidal particles and the conductivity of their dilute suspensions, J. Chem. Soc. Faraday
Trans. 2, 79. 1613 (1983).

181

Author’s short biography
Kang Kim
Kang Kim obtained his B.Sc., M.Sc. and Ph.D degrees from Kyoto University in
1998, in 2000 and in 2003, respectively. He then joined a research project, “Hybrid
Molecular Dynamics Simulation for Soft Matters,” as a post doctoral fellow under
the direction of Prof. R. Yamamoto until June, 2006. He is currently a Research Associate at Institute for Molecular Science. His research interests include computer
simulations for analyzing dynamics in simple and complex liquids.

Yasuya Nakayama
Yasuya Nakayama received his PhD in Informatics from Kyoto University in 2003,
MSc from Kyoto University in 2000 and BSc from Kyushu University in 1998. He
joined a research project, “Hybrid Molecular Dynamics Simulation for Soft Matters,” as a post doctoral fellow under the direction of Prof. R. Yamamoto till October, 2005. Currently he is a Research Associate at Kyushu University. His current
activities are focused on physics and simulations of dispersions and polymer processing.
Ryoichi Yamamoto
Ryoichi Yamamoto received his PhD in Engineering from Kyoto University in 1996,
MEn and BEn from Kobe University in 1992 and 1988. He became a Research Associate at Kobe University (1994) and moved to Department of Physics, Kyoto University where he became a Research Associate (1996) and a Lecturer (2000). He isan Associate Professor at Department of Chemical Engineering, Kyoto University
since 2004. He concerns mainly to understand the nature of unique properties ofsoft matters and complex fluids by using computer simulations. He is interested
also in developments of new simulation methods applicable to complex materials.

182

KONA No.24 (2006)