Simulating Colloids and Self-Propelled Particles
with Fully Resolved Hydrodynamics Using the
Smoothed Profile Method (SPM)
Ryoichi Yamamoto, John J. Molina, and Rei Tatsumi
Department of Chemical Engineering, Kyoto University, Kyoto 615-8510, Japan
E-mail: {ryoichi, john, tatsumi}@cheme.kyoto-u.ac.jp
The smoothed profile method (SPM) provides an efficient numerical scheme for coupling continuum fluid dynamics with moving dispersed particles using a smeared interface between thefluids and the particles. The SPM has been successfully applied to directly simulate several
dynamical problems of colloidal dispersions in incompressible fluids, including those involving sedimentation, diffusion, coagulation, rheology, and tumbling motion in shear flow as wellas electrophoresis in external electric fields. More recently, the SPM was extended to two important problems. The first extension simulates colloidal particles in compressible host fluids,
whereas the second extension simulates self-propelled swimming particles. A comprehensive
summary of SPM is provided in this paper.

1

Introduction

Interparticle interactions in colloidal dispersions mainly consist of thermodynamic potential interactions as well as hydrodynamic interactions. Whereas the former applies to both
static and dynamic situations, the latter only applies to dynamic situations. Although thermodynamic interactions in static situations have been studied extensively and are treatedas effective interactions, the nature of dynamic interactions is poorly understood. Because hydrodynamic interactions are essentially long-range, many-body effects, they are
extremely difficult to study using analytical means alone. Numerical simulations can beused to investigate the role of hydrodynamic interactions in colloidal dynamics.
Several numerical methods have been developed to simulate the dynamics of colloidal
dispersions. Two of the most well-known methods include Stokesian dynamics1 and theEulerian–Lagrangian method2 . The former is the most widely used method because of
its proper treatment of hydrodynamic interactions between spherical particles in a Newtonian fluid at zero Reynolds number. Furthermore, it can be implemented as a O(N )scheme for N particles by utilizing the fast multi-pole method3 . However, it is not easy toaddress dense dispersions and dispersions consisting of non-spherical particles by meansof Stokesian dynamics due to the complicated mathematical structures used in Stokesian
dynamics. In contrast, the Eulerian–Lagrangian method is a very natural and sensible approach for stimulating solid particles. It is possible to apply this method to dispersions
consisting of many particles with different shapes. However, numerical efficiencies arise
from the following concerns: i) the re-construction of irregular meshes according to thetemporal particle position is necessary for every simulation step, and ii) the Navier–Stokes
equation must be solved with boundary conditions imposed on the surfaces of all colloidal
particles. Thus, these computational demands are particularly cumbersome for systems
involving many particles, even if the shapes are all spherical.

11

To overcome problems arising from the particle-fluid interface in the Eulerian–
Lagrangian method, we have developed an efficient direct numerical simulation (DNS)
method for colloidal dispersions. This method was named the ”smoothed profile method(SPM)” because the original sharp interface between the colloids and solvent is replaced
by a smeared out, smoothed interface with a finite thickness4–20 . This simple modification
greatly improves the resulting quality of the of numerical computations in comparison withthe original Eulerian–Lagrangian method for the following reasons:
1. Regular fixed Cartesian coordinates can be used for many particle systems by defining
a particle shape instead of providing boundary-fitted coordinates. The particle-fluid
interface has a finite volume (∝ πad−1 ξ, with a and d as the particle radius andsystem dimension) supported by multiple grid points. Thus, the round particles canbe treated in a fixed Cartesian coordinates without any difficulties. The simulation
scheme is thus free from the mesh re-construction problem that significantly suppresses the computational efficiency of the Eulerian–Lagrangian method. In addition,
the simple Cartesian coordinate enables the use of periodic boundary conditions aswell as fast Fourier transformations (FFT).
2. At the particle-fluid interfaces, the velocity component in the direction normal to theinterface of the host fluid must be equal to that of the particle. In the Eulerian–
Lagrangian method, this non-penetration condition is imposed by the Navier–Stokes
equation as the boundary condition defined for the particle-fluid interface. In the SPM,
however, this condition is automatically satisfied by an incompressibility condition on
the entire domain.
3. The computational demands for this method include sensitivity to the number of grid
points (volume of the total system). Nevertheless, because the method is insensitive
to the number of particles, it is suitable for simulating dense colloidal dispersions.
The SPM has been successfully used to directly simulate various dynamical problemsof colloidal dispersion in incompressible fluids, including sedimentation20 , diffusion9, 12, 13 ,coagulation8, 19 , rheology11, 14, 17 , tumbling motion in shear flow15 , and electrophoresis inexternal electric fields7, 10 . Several simulation methods similar in spirit to the SPM havealso been proposed in recent publications22–25 . A comprehensive summary of SPM is provided in this study including the recent key extensions for stimulating colloidal particles incompressible host fluids18 and also for stimulating the self-propelled swimming of particles21 .

2

Colloids in Incompressible Fluids

2.1

Working Equations

The motion of the host fluid is determined by the Navier-Stokes equation with the following
incompressibility condition:∇ · uf = 0
ρ (∂t + uf · ∇) uf = ∇ · σ

12

(1)
(2)

where ρ is the total mass density of the fluid, uf is the host fluid velocity field, σ is thestress tensor
σ = −pI + σ 0t
σ 0 = η ∇uf + (∇uf )

(3)
(4)

and η is the shear viscosity of the fluid. Consider a mono-disperse system containing N spherical particles with a radius a, mass Mp , and moment of inertia Ip = 2/5Mp a2 I (withI the unit tensor). The evolution of colloids is described by the Newton-Euler equations28 ,Q̇i = Qi skew (Ωi )Ṙi = Vi
Mp V̇i = FiH + FiC + Fiext(5)

Ip · Ω̇i = NiH + N ext

where Ri and Vi denote the centre of mass positions and the velocity of the particle
i, respectively, and Qi is the orientation matrixb . Hence, Ωi the angular velocity andskew (Ωi ) is the skew-symmetric angular velocity matrix:
0 −Ωzi Ωyiskew(Ωi ) =  Ωzi 0 −Ωxi (6)−Ωyi Ωxi 0
The forces on the particles are comprised of hydrodynamic contributions arising from fluidparticle interactions F H , colloid-colloid interactions due to the core particle potential F C(which prevents particle overlap), and a possible external field contribution F ext (such as
gravity). Likewise, the torques on the particles can be divided into a hydrodynamic N H
and an external contribution N ext (for simplicity, the particle-particle interactions are assumed to be described by a radial potential). Subsequently, we consider buoyancy-neutral
particles, for which F ext = N ext = 0. Finally, the conservation of momentum between
the fluid and the particles implies the following hydrodynamic force and torque on the i-th
particle:Z
FiH = dSi · σ(7)Z
NiH = (x − Ri ) × (dS i · σ)
(8)R
where dSi indicates an integral over the surface of the particle. In addition, thermal
fluctuations can be introduced by adding a random stress tensor sin Eq. 3, that satisfies thefluctuation-dissipation relation:26
hsik (x, t)sjl (x0 , t0 )i = 2kB T η(δij δkl + δil δkj )δ(x − x0 )δ(t − t0 ),(9)where kB is the Boltzmann constant, T is the temperature. Alternatively, it is also possible to introduce thermal fluctuations by adding Langevin random forces and torque to
Eq. 59, 11–13 .
b For numerical stability, we use quaternion instead of rotation matrices to represent the rigid body dynamics ofthe particles.

13

2.2

Simulation Procedure for Incompressible Fluids

We now present the computational algorithm used to simulate the motion of sphericalparticles using the SPM. We require that all field variables are defined over the entire computational domain (fluid + particle). The concentration field for the colloids is describedas follows:
φ(x, t) =NX

φi (x, t),

(10)

i=1

where φi ∈ [0, 1] is the smooth profile field of particle i. This field is defined as unity within
the particle domain, as zero in the fluid domain, and as a smooth interpolation between thetwo extremes within the interface region. Several possible mathematical forms exist forφi (x), however, we adopted the following definition of φi :
φi (x) = g(|x − Ri |),h((a + ξ/2) − x),g(x) =
h((a + ξ/2) − x) + h(x − (a − ξ/2))
exp −∆2 /x2 x ≥ 0,h(x) =
0x < 0.

(11)
(12)
(13)

where a, ξ, and ∆ are the radius of the particle, the interfacial thickness, and lattice spacing,respectively. The particle velocity field is defined in a similar fashion:
φup (x, t) =N
X

{Vi (t) + Ωi (t) × ri (t)} φi (x, t)

(14)

i=1

with ri = x − Ri , which allows one to define the total fluid velocity field using thefollowing expression:u(x, t) ≡ (1 − φ)uf + φup(15)

where the incompressibility condition is satisfied over the entire domain ∇ · u = 0. The
evolution equation for u is then derived assuming momentum-conservation between fluid
and particles6, 10ρ (∂t + u · ∇) u = ∇ · σ + ρφfp

(16)

where φfp represents the force density field needed to maintain rigidity constraints on theparticle velocity field.
We use a fractional step approach to update the total velocity field. Let un be the fieldat time tn = nh (h is the time interval).
i) We first solve for advection and hydrodynamic viscous stress terms, and we thenpropagate the particle positions (orientations) using the current particle velocities,

14

which affords the following relation:∗

Z tn +hnu =u +
Rin+1 = Rin +Qn+1= Qni +
itn
Z tn +htn
Z tn +h1
ds ∇ ·(−p∗ I + σ 0 ) − uu
ρ



(17)

ds Vi

(18)

ds Qi skew (Ωi )

(19)

tn

where the pressure term p∗ in Eq. 17 is determined by the incompressibility condition∇ · u∗ = 0. The remaining updating procedure imposes a rigidity constraint on thevelocity field.
ii) The hydrodynamic force and torque exerted by the fluid on the colloids is determinedby assuming momentum conservation. The time integrated hydrodynamic force andtorque over a period h are equal to the momentum exchange over the particle domain
"Z
# Z
tn +h

H
ds Fi = dx ρφn+1u∗ − unp
(20)
i
tn

"Z
tn +h

#
ds NiHZ=dx rin+1 × ρφn+1u∗ − unp
i

(21)

tn

Based on this and other forces acting on the colloids, the particles velocities are updated as follows:
"Z
#
"Z
#
tn +h
tn +h

n+1
n
−1
H
−1
C
ext
ds Fi + Fi
(22)
Vi
= Vi + Mp
ds Fi + Mp
tn

Ωn+1
= Ωni + Ip−1 ·
i

"Z

tn +h

tn

#
ds NiH + Ip−1 ·

tn

"Z
tn +h

#
ds Niext(23)tniii) Finally, the resulting particle velocity field φn+1 un+1
is enforced on the total velocitypfield as follows:
"Z
#
tn +h
n+1
∗
u
=u +
ds φfp
(24)
tn

"Z

tn +h

tn

#

 h
ds φfp = φn+1 un+1
− u∗ − ∇pp
p
ρ

(25)

wherein the pressure is due to the rigidity constraint obtained from the incompressibility condition ∇·un+1 = 0. The total pressure field is thus obtained as p = p∗ +pp .
The above procedure defines the consistent time-propagation, {un ; Rin , Qni , Ωni } →{un+1 ; Rin+1 , Qn+1
, Ωn+1}, to simulate colloidal particles in incompressible fluids.
i
i

15

3

Colloids in Compressible Fluids

3.1

Working Equations

The hydrodynamic equations consist of three conservation laws concerning mass, momentum, and energy. The conservation equations of mass and momentum for incompressible
fluids are described by the following:∂ρ+ ∇ · m = 0,
∂t

(26)

∂m
+ ∇ · (mu) = ∇ · σ + ρφfp ,
∂t

(27)

where m(r, t) = ρ(r, t)u(x, t) is the momentum density field. We consider a compressible Newtonian fluid, and the stress tensor is described as follows:


2
(28)
σ = −pI + η[∇u + (∇u)t ] + ηv − η (∇ · u)I,3
where p(r, t) is the pressure, η is the shear viscosity, and ηv is the bulk viscosity. A bodyforce ρφfp is also added to satisfy the rigidity of the particles. Additionally, we assume abarotropic fluid described by p = p(ρ), with a pressure gradient that is proportional to thedensity:∇p = c2 ∇ρ,(29)

where c is the speed of sound in the fluid. Eqs. 26-29 are closed to variables ρ, m, and p;
therefore, energy conservation does not need to be considered for barotropic fluids.
The motion of the dispersed particles is governed by Newton-Euler equations of motion
Eq. 5. The effect of thermal fluctuations on the particles dynamics is important when theparticle size is on the order of a micrometer or smaller. Fluctuations were introduced usinga random stress tensor s, which is added to the stress tensor Eq. 28. The random stress isa stochastic variable satisfying the fluctuation-dissipation relation26 :hsij (r, t)skl (r 0 , t0 )i = 2kB T ηijkl δ(r 0 − r)δ(t0 − t),

(30)

and



2
ηijkl = η(δik δjl + δil δjk ) + ηv − η δij δkl .
3(31)

Brownian motion of the dispersed particles is induced by the random stresses acting on
the fluid. Thermal fluctuations can be introduced using the Langevin approach, whereinrandom forces are exerted on the particles9, 11–13 . However, this approach does not accurately represent the short-time dynamics of the system because the autocorrelation timeof the hydrodynamic force acting on the particles is neglected. Therefore, the fluctuatinghydrodynamics approach is more appropriate for investigating dynamics at a time scale ofsound propagation.

16

3.2

Simulation Procedure for Compressible Fluids

In this section, the time-discretized evolution of the equations is derived for colloidal dispersions in compressible fluids. The time evolution of the fluid is determined using thefollowing steps:i) The mass and momentum density changes associated with sound propagation are calculated using the following equations:
Z tn +h
n+1
n
ρ
=ρ −
ds ∇ · m,(32)tnm∗ = mn − c2

Z tn +hds ∇ρ.

(33)

tn

When we assume a periodic boundary condition and use the Fourier spectral method,a semi-implicit scheme becomes feasible27 . This situation eliminates restrictions ontime increments with a small compressibility factor ε.
ii) The time evolution of the advection and viscous diffusion terms are calculated using
the following equations:
Z tn +hm∗∗ = m∗ +ds ∇ · (σ 0 − mu),(34)
tn

where σ 0 is the dissipative stress defined in Eq. 3.
iii) In concert with the advection of the particle domain, the position (orientation) of eachdispersed particle evolves according to the following equations:
Z tn +hn+1nRi= Ri +
ds Vi
(35)
tn

Qn+1
= Qni +i

Z tn +hds Qi skew (Ωi ) .

(36)

tn

iv) The hydrodynamic force and torque are derived by considering the conservation of
momentum. The time-integrated hydrodynamic force and torque are computed using
the following equations:
Z tn +hZ
ds FiH = dx φn+1(m∗∗ − ρn+1 unp ),
(37)
i
tn

Z tn +hds NiH =Z

dx [(r − Rin+1 ) × φn+1(m∗∗ − ρn+1 snp )].
i

(38)

tn

With these and other forces acting on the particles, the translational and rotational
velocities of each dispersed particle evolve according to the following equations:
Z tn +hn+1n−1Vi= Vi + Mp
ds (FiH + FiC ),
(39)
tn

Ωn+1
= Ωni + Ip−1 ·i

Z tn +htn

17

ds NiH .

(40)v) The updated velocity of the particle region is imposed on the velocity field as thevolume force ρφfp .
Z tn +h
n+1
∗∗
m
=m +
ds ρφfp ,(41)
tn

Z tn +hds ρφfp = φn+1 (ρn+1 un+1
− m∗∗ ).
p

(42)

tn

In the case of an incompressible fluid, the pressure is spontaneously determined by
the solenoid condition of the velocity field. In contrast, in this case, the pressure andmass density variations are independent of the velocity field.
The above procedure defines the consistent time-propagation, {ρn , mn ; Rin , Qni , Ωni } →{ρn+1 , mn+1 ; Rin+1 , Qn+1
, Ωn+1}, to simulate colloidal particles in compressible fluids.
ii

4

Self-Propelled Particles

4.1

Squirmer Model

We consider a simple model of self-propelled spherical swimmers, originally introduced
by Lighthill29 and later extended by Blake30 , which move due to a self-generated surfacetangential velocity us . This specific mechanism was proposed as a model for an ideal
ciliate particle, in which the synchronized beating of the cilia at the surface gives rise to
net motion in the absence of any external fields. If one assumes that the displacements ofthis cilia envelope are purely tangential, then the effective (time-averaged) slip velocity forthese squirmers is described by the following equation30 :us (r̂) =∞X

2Bn (ê · r̂ r̂ − ê) Pn0 (ê · r̂)
n
(n
+
1)
n=1

(43)

where ê is the squirmer’s fixed swimming axis (i.e., we consider that each squirmer carries
with it a fixed coordinate system that determines its preferred swimming direction at each
instant), r̂ is a unit vector from the particle centre to a point on the surface, Pn0 is the derivative of the n-th order Legendre polynomial, and Bn is the amplitude of the corresponding
mode.
When all squirming modes higher than three are neglected, Bn = 0 (n ≥ 3), thefollowing simple expression for the surface tangential velocity as a function of polar angle
θ = cos−1 (r̂ · ê), is obtained:αus (θ) = B1 sin θ + sin 2θ
(44)2
where α = B2 /B1 determines whether the swimmer is a pusher (α < 0) or a puller(α > 0). A schematic representation of the flow profile generated by these two typesof swimmers is provided in Fig. 1. An example of the former include spermatozoa andmost bacteria, whereas the latter includes unicellular algae Chlamydomonas. Although thesquirmer model we adopt does not include a detailed propulsion mechanism, it is capableof distinguishing between pushers/pullers and provides an adequate approximation for the
far-field flow profile generated by these swimmers.

18

U

U

slip velocity

down

up

D

2

PusherD

2

PullerFigure 1. Schematic representation of the flow profiles generated by a pusher (left) and a puller (right). In bothcases, the particle’s swimming direction is towards the top of the page.

For Newtonian fluids, which is the only case considered here, the swimming speed Uof the squirmer is determined uniquely by the first mode B1 , irrespective of the size ofthe particle, as U = 2/3B1 , while the second mode gives the strength of the stresslet31, 32 .
In the Stokes regime, the velocity field generated by a single such squirmer was solved
analytically by Ishikawa et al.31 , providing the following expression in the laboratory frame(fluid at rest far away from the particle):

 2


 

a2 a 2a2αa2u(r) = B1 2
3
cos
θ
−
1
r̂
+
ê + sin θ θ̂ +
−
1
sin
2θ
θ̂
r r 3
2
r2
r2
(45)
where a is the radius of the particle. Notice that for neutral swimmers (α = 0), the velocity
field decays as r−3 , whereas for pushers/pullers (α 6= 0), the velocity field decays as r−2 .
In contrast, the velocity field for a sedimenting particle (or a particle experiencing a net
body force) decays as r−1 33 . This observation will have important consequences on thehydrodynamic interactions describing suspensions of swimmers.
4.2

Simulation Procedure for Squirmers

We now present the computational algorithm used to simulate the motion of sphericalparticles, with a given surface tangential slip velocity us using the SPM. The evolution
equation for u is then derived by assuming momentum-conservation between the fluid andparticles6, 10ρ (∂t + u · ∇) u = ∇ · σ + ρφfp + ρfsq

(46)

where φfp represents the force density field needed to maintain the rigidity constraint
on the particle velocity field and fsq is the force density field generated by the squirming
motion of the particles. The motion of the dispersed particles is governed by Newton-Euler
equations of motion Eq. 5.
We use the fractional step approach to update the total velocity field.

19i) We first solve for the advection and hydrodynamic viscous stress terms, and we thenpropagate the particle positions (orientations) using the current particle velocities.
This operation yields the following results:
Z tn +h1∗n∗0
u =u +(−p I + σ ) − uuds ∇ ·(47)ρ
tn
Z tn +hRin+1 = Rin +
ds Vi
(48)
tn

Qn+1
= Qni +i

Z tn +hds Qi skew (Ωi )

(49)

tn

where the pressure term p∗ in Eq. 47 is determined by the incompressibility condition∇ · u∗ = 0. The remaining updating procedure applies to the slip condition at theparticle boundary as well as the rigidity constraint on the velocity field.
ii) We now consider the momentum change needed to maintain the slip velocity at thesurface of each of the squirmers, where the slip profile us is imposed with respect tothe particle velocities {Vi0 ; Ω0i }, using the previously updated positions and orientations orientations {Rin+1 ; Qn+1}. We note that at this point we do not yet know theicorrect updated particle velocities {Vin+1 ; Ωn+1}, which are the values that shouldibe used when enforcing the surface slip profile Vi0 = Vin+1 (Ω0i = Ωn+1). Therefore,i
we adopt an iterative solution, and as an initial guess, we use the particle velocities atthe previous time step, i.e., Vi0 = Vin (Ω0i = Ωni ). The updated total velocity field is
now obtained using the following:
"Z
#
tn +h
∗∗
∗
u =u +
ds fsq
(50)
tn

"Z

tn +h

#
ds fsq = u∗ +

tn

N
X

ϕi (Vi0 + Ω0i × ri + usi − u∗ )

i=1

+

N
X
i=1

φi (δVi + δΩi × ri ) −

h
∇psq
ρ

(51)

The second term on the right hand side of Eq. 51 imposes a slip velocity profile us at
the surface of each of the squirmers where ϕi ∝ (1 − φi ) |∇φi | is a smooth surface
profile function that is non-zero only within the interface domain of the squirmer (normalized to have a maximum value of one), and zero everywhere else (the red arrows
in Fig. 2). The third term adds a counter-flow entirely within the particle domain, suchthat local momentum conservation is preserved (the blue arrows in Fig. 2). Assuming
rigid-body motion, with velocities δVi and δΩi , this requiresZZ
dx φi (δVi + δΩi × ri ) = − dx ϕi (Vi0 + Ω0i × ri + usi − u∗ ) (52)ZZ
dx ri × φi (δVi + δΩi × ri ) = − dx ri × ϕi (Vi0 + Ω0i × ri + usi − u∗ )
(53)

20

[ 2

aa

tangential slip

propulsion

Figure 2. Schematic representation of the updating scheme used to enforce the slip boundary condition at thesurface of the squirmers. Each particle is considered to exert a force on the fluid at the surface, in order to maintainthe specified flow profile us (red arrows) for the squirming motion. To ensure local momentum conservation, acounter-flow is added within the particle domain (blue arrows).

from which we can easily obtain the counter-flow terms δVi (δΩi ) from the particle
velocities Vi0 (Ω0i ). A schematic representation of the procedure used to enforce thespecific slip-boundary conditions for our model squirmers is shown in Fig. 2. Finally,
the pressure term due to the squirming motion psq is obtained from the incompressibility condition ∇ · u∗∗ = 0. At this point, the momentum conservation relation issolved for the total velocity field.
iii) The hydrodynamic force and torque exerted by the fluid on the colloids (which includes all contributions to the squirming motion) is again derived by assuming momentum conservation. The time integrated hydrodynamic force and torque for a period h are equal to the momentum exchange over the particle domain:
"Z
# Z
tn +h

sq 
H
ds Fi + Fi
= dx ρφn+1u∗∗ − unp
(54)
i
tn

"Z

tn +h

ds


NiH + Nisq

#

Z
=



dx rin+1 × ρφn+1u∗∗ − unp
i

(55)

tn

From this and any other forces on the colloids, the particles velocities are updated
according to the following equations:
"Z
#
"Z
#
tn +h
tn +h

sq 
n+1
n
−1
H
−1
C
ext
Vi
= Vi + Mp
ds Fi + Fi
+ Mp
ds Fi + Fi
tn

tn

(56)
Ωn+1
= Ωni + Ip−1 ·
i

"Z
tn +h

ds NiH + Nisq

tn

#


+ Ip−1 ·

"Z

tn +h

#
ds Niext(57)

tn

We recall that we have imposed the slip profile us with respect to the primed velocities{Vi0 ; Ω0i }, which need not be equal to the final velocities of the particle at step n + 1.
To maintain consistency, we iterate over Eqs. 50-57 until a convergence of velocities
is achieved.

21iv) Finally, the resulting particle velocity field φn+1 un+1
is enforced over the total vep
locity field using the following relations:"Z#
un+1 = u∗∗ +tn +h

ds φfp

(58)

tn

"Z

tn +h

tn

#

 h
ds φfp = φn+1 un+1
− u∗∗ − ∇pp
p
ρ

(59)

wherein the pressure due to the rigidity constraint obtained from the incompressibility
condition ∇ · un+1 = 0. The total pressure field is then given by p = p∗ + pp + psq .
The above procedure defines the consistent time-propagation, {un ; Rin , Qni , Ωni } →{un+1 ; Rin+1 , Qn+1
, Ωn+1}, to simulate self-propelled squirmers in incompressible fluii
ids.
We are aware of two alternative simulation methods that aim to describe these squirmer
suspensions at the same level of description, the first was developed by Ramachandran etal.35 using a Lattice Boltzmann model, and the second was originally introduced by Downton and Stark36 within a multi-particle collision dynamics framework, and later extended
by Götze and Gompper37 to recover the correct rotational dynamics. For the moment
though, these DNS approaches have not been extensively used to study these types ofswimming systems; the most popular methods, which still account for the hydrodynamic
interactions, have usually been based on Stokesian Dynamics31, 34 , and are thus limited to
Newtonian fluids in the Stokes regime.

5

Concluding Remarks

A new computational method named the SPM has been developed to simulate particle
dispersion in fluids4–21 . Utilizing a smoothed profile for particle-fluid boundaries, hydrodynamic interactions in many particle dispersions can be fully taken into account, yielding
both accurate and efficient results. In principle, the SPM can be easily applied to systems
consisting of many particles with different shapes. The reliability and the performance ofthe method was confirmed to be satisfactory by several critical tests4–22 .
Recently, we extended the SPM to particle dispersions in compressible fluids18 . The
validity of the method was confirmed by calculating the velocity relaxation function ofa single spherical particle in a compressible fluid21 . The effect of compressibility on thevelocity relaxation was also observed, revealing a two-stage relaxation process for lowcompressibility fluids and a backtracking motion for high-compressibility fluids. A simulation of the motion of a single spherical particle in a fluctuating fluid was also performed.
The calculated velocity autocorrelation function of the particle showed good agreementwith the analytical solution of the relaxation function, thereby confirming the validity of
the fluctuation-dissipation theorem without any fitting parameters.
We have also shown that SPM can be extended to systems with self-propelled swimming particles, making it possible to describe the actions of squirmers (active swimmersthat move due to self-generated surface tangential velocities)21 . The validity of the method
was confirmed by comparing the simulation data with the exact results for the case of a22single swimmer, wherein the correct swimming speed is recovered and it is possible to accurately reproduce the fluid flow generated by the squirming motion. The advantage of the
SPM for swimming particles in comparison with Stokesian Dynamics (which have beensuccessfully and extensively used to study these systems)31, 34 is its applicability to particle
dispersions in complex fluids. This is relevant in the case of swimming micro-organisms
as the role of nutrients and the presence of a non-Newtonian host fluid must be consideredwhen making comparisons with experimental data.

Acknowledgements
The authors would like to express their gratitude to Dr. K. Kim, Dr. Y. Nakayama,
Dr. T. Iwashita, and Dr. H. Kobayashi for their contributions to the development of the
SPM. This work was supported by Grant No. KAKENHI 23244087.

References
1. J. F. Brady and G. Bossis, Ann. Rev. Fluid Mech., 20, 111 (1988).
2. H. H. Hu, N. A. Patankar, M. Y. Zhu, J. Comput. Phys., 169, 427 (2001).
3. K. Ichiki, J. Fluid Mech., 452, 231 (2002).
4. http://www-tph.cheme.kyoto-u.ac.jp/kapsel/
5. R. Yamamoto, Phys. Rev. Lett., 87, 075502 (2001).
6. Y. Nakayama and R. Yamamoto, Phys. Rev. E, 71, 036707 (2005).
7. K. Kim, Y. Nakayama, and R. Yamamoto, Phys. Rev. Lett., 96, 208302 (2006).
8. R. Yamamoto, K. Kim, Y. Nakayama, K. Miyazaki, and D. R. Reichman, J. Phys.
Soc. Japan, 77, 084804 (2008).
9. T. Iwashita, Y. Nakayama, and R. Yamamoto, J. Phys. Soc. Japan, 77, 074007 (2008).
10. Y. Nakayama, K. Kim, and R. Yamamoto, Eur. Phys. J. E, 26, 361–368 (2008).
11. T. Iwashita and R. Yamamoto, Phys. Rev. E, 80, 061402 (2009).
12. T. Iwashita and R. Yamamoto, Phys. Rev. E, 79, 031401 (2009).
13. T. Iwashita, Y. Nakayama, and R. Yamamoto, Progress of Theoretical Physics, 11,
86–91 (2009).
14. T. Iwashita, T. Kumagai, and R. Yamamoto, Eur. Phys. J. E, 32, 357–363 (2010).
15. H. Kobayashi and R. Yamamoto, Phys. Rev. E, 81, 041807 (2010).
16. S. Jafari, R. Yamamoto, and M. Rahnama, Phys. Rev. E, 83, 026702 (2011).
17. H. Kobayashi and R. Yamamoto, Phys. Rev. E, 84, 051404 (2011).
18. R. Tatsumi and R. Yamamoto, Phys. Rev. E, 85, 066704 (2012).
19. Y. Matsuoka, T. Fukasawa, K. Higashitani, and R. Yamamoto, Phys. Rev. E, 86,
051403 (2012).
20. A. Hamid and R. Yamamoto, J. Phys. Soc. Japan, in print.
21. J. J. Molina, Y. Nakayama, and R. Yamamoto, arXiv:1212.6133.
22. X. Luo, M. R. Maxey, and G. E. Karniadakis, J. Comp. Phys., 228, 1750 (2009).
23. M. Fujita and Y. Yamaguchi, Phys. Rev. E, 77, 026706 (2008).
24. F. B. Usabiaga, I. Pagonabarraga, and R. Delgado-Buscalioni, J. Comp. Phys., 235,
701–722 (2012).
25. F. B. Usabiaga, R. Delgado-Buscalioni, B. E. Griffith, A. Donev, arXiv:1212.6427.

23

26. L. D. Landau and E. M. Lifshitz, Fluid Mechanics, (Pergamon Press, London, 1959).
27. G. Erlebacher, M. Y. Hussaini, H. O. Kreiss, and S. Sarkar, Theor. Comput. Fluid.
Dyn., 2, 73 (1990).
28. J. V. José and E. J. Saletan, Classical Dynamics: A Contemporary Approach, (Cambridge University Press, New York, 1998).
29. M. J. Lighthill, Ann. Rev. Fluid Mech., 1, 413–446 (1969).
30. J. R. Blake, J. Fluid Mech., 46, 199–208 (1971).
31. T. Ishikawa, M. P. Simmonds, and T. J. Pedley, J. Fluid Mech., 568, 119–160 (2006).
32. L. Zhu, E. Lauga, and L. Brandt, Phys. Fluids, 24, 051902 (2012).
33. W. B. Russel, D. A. Saville, and W. R. Schowalter, Colloidal Dispersions, (Cambridge
University Press, Cambridge, 1992).
34. J. W. Swan, J. F. Brady, R. S. Moore, and ChE 174, Phys. Fluids, 23, 071901 (2011).
35. S. Ramachandran, P. B. Sunil Kumar, and I. Pagonabarraga, Eur. Phys. J. E, 20, 151(2006).
36. M. T. Downton and H. Stark, J. Phys.: Condens. Matt., 21, 204101 (2009).
37. I. Götze and G. Gompper, Phys. Rev. E, 82, 041921 (2010).

24