EPJ E
Soft Matter andBiological Physics

EPJ .org
your physics journal

Eur. Phys. J. E 032, 357–363 (2010)

DOI: 10.1140/epje/i2010-10638-7

A direct numerical simulation method for complex modulus of particle dispersions
T. Iwashita, T. Kumagai and R. Yamamoto

Eur. Phys. J. E 32, 357–363 (2010)
DOI 10.1140/epje/i2010-10638-7

THE EUROPEAN
PHYSICAL JOURNAL E

Regular Article

A direct numerical simulation method for complex modulus ofparticle dispersions
T. Iwashita1,2,a , T. Kumagai1 , and R. Yamamoto1,2,b
1
2Department of Chemical Engineering, Kyoto University, Kyoto 615-8510, Japan
CREST, Japan Science and Technology Agency, Kawaguchi 332-0012, Japan
Received 31 August 2009 and Received in ﬁnal form 30 November 2009c EDP Sciences / Società Italiana di Fisica / Springer-Verlag 2010
Published online: 11 August 2010 – 
Abstract. We report an extension of the smoothed proﬁle method (SPM) (Y. Nakayama, K. Kim, andR. Yamamoto, Eur. Phys. J. E 26, 361 (2008)), a direct numerical simulation method for calculating thecomplex modulus of the dispersion of particles, in which we introduce a temporally oscillatory external
force into the system. The validity of the method was examined by evaluating the storage G′ (ω) andloss G′′ (ω) moduli of a system composed of identical spherical particles dispersed in an incompressible
Newtonian host ﬂuid at volume fractions of Φ = 0, 0.41, 0.46, and 0.51. The moduli were evaluated atseveral frequencies of shear ﬂow; the shear ﬂow used here has a zigzag proﬁle, as is consistent with theusual periodic boundary conditions. The simulation results were compared with several experiments forcolloidal dispersions of spherical particles.

1 Introduction
The viscoelastic properties of the dispersion of solid particles are of particular importance in several scientiﬁc,engineering, and industrial ﬁelds. These properties arestrongly aﬀected not only by direct interactions between
particles, but also by thermal ﬂuctuations of the systemand the hydrodynamic interactions (HI) acting on dispersed particles, mediated by the surrounding ﬂuid [1, 2].
The viscoelasticity of materials are most commonly characterized by the complex modulus, which represents themechanical response of materials to temporally oscillating small shear deformations of frequency ω. The complex
modulus consists of elastic and viscous components (i.e.,the storage modulus G′ (ω) and the loss modulus G′′ (ω)).
In general, materials tend to lose elastic energy by dissipation when the frequency of external deformation is low,
but can store elastic energy when the frequency is high.
Understanding viscoelastic behavior has been the subject of both fundamental and technological interest. Especially, the linear and non-linear responses of colloidal systems, including colloidal glasses and aggregated colloidal
suspensions and gels, have attracted much attention inrecent years [3–5]. While extensive experimental studies
have been carried out to determine the viscoelastic properties of particle dispersions over a wide range of volume
fractions, from dilute (Φ ≃ 0) to very dense (Φ ≥ 0.5),
theoretical studies encounter a fundamental problem whenabe-mail: tiwashit@utk.edu
e-mail: ryoichi@cheme.kyoto-u.ac.jp

the density of dispersion is high, due to many-body eﬀects
and the long-range nature of the HI among dispersed particles. Developing analytical theories becomes even morediﬃcult when the host ﬂuids are complex, such as electrolytes or viscoelastic media. Computer simulations canthus be very powerful tools in theoretical investigations ofdense dispersions in general.
The widely used simulation technique is the Stokesian
dynamics (SD) method [6], which is based on the Stokes
approximation (Re → 0), involves near-ﬁeld lubrication
forces and far-ﬁeld many-body HIs. The SD method has
succeeded in simulating the motions of particle dispersions at steady shear. However the viscoelasticity of particle dispersions at oscillatory shear has not been examined.
Furthermore, the SD method is valid only for simulating
the particles in a Newtonian ﬂuid and is not applicable tothe motions of particles in complex ﬂuids, such as charged
particles in electrolytes or viscoelastic media.
Various numerical methods have recently been proposed for simulating dense dispersions, including somecases where the host ﬂuids are complex. In order to accurately track the motions of host ﬂuids as well as the motions of dispersed particles, a variety of numerical methods have been developed. One of these methods, which
we call the direct numerical simulation (DNS) method
for particle dispersions, involves solving the Navier-Stokes(NS) equation for host ﬂuids in a manner consistent withboundary conditions deﬁned according to the particle positions and motions. Aside from several successful implementations of DNS methods for particle dispersions [7–15],
other methods have been developed that do not rely on

358

The European Physical Journal E

the Navier-Stokes equation to resolve ﬂuid motions. The
coupling methods of particle dynamics with the latticeBoltzmann (LB) method [16] or stochastic rotational dynamics (SRD) [17, 18] are the most popular alternatives
to DNS, because LB and SRD are believed to be morecomputationally eﬃcient.
Problems studied here have a very small Reynold number and the full NS equations for the host ﬂuid are notnecessarily solved to simulate the motions of particles in a
Newtonian ﬂuid. It is suﬃcient to use the methods basedon the Stokes approximation. However there are several
advantages for solving the NS equations for the host ﬂuid.
Those techniques mentioned above enable us to simulate
the short-time motions of particle dispersions, where thecoupling between the particle motion and the ﬂuid motion
remains strong, even if the Reynolds number is very small.
One of the typical behavior is the power law decay in thevelocity correlations of particles, which is known as “longtime tail”. In addition, those methods have succeeded insimulating several complex phenomena, such as the structure formation in colloidal gels where the pressure ﬁled of
the ﬂuid plays an important role [19] and electrophoresisof charged colloidal dispersions [20, 21].
Although extensions of those methods have alreadybeen developed for steady shear ﬂow with DNS [22],
LB [23–25], and SRD [26], there exists no successful attempt for calculating the complex moduli of the dispersions with an imposed oscillatory shear ﬂow. The purpose of the present paper is to report our successful reproductions of experimentally observed behaviors of thestorage G′ (ω) and loss G′′ (ω) moduli in temporally oscillating shear ﬂow.
To check the validity of the method, we ﬁrst apply themethod to a Newtonian host ﬂuid of viscosity η, whichshould exhibit a purely viscous response, G′ (ω) = 0 andG′′ (ω) = ηω. We then apply the method to dense dispersions composed of identical spherical particles in a Newtonian host ﬂuid. Since several experimental measurementsof the frequency dependence of G′ (ω) and G′′ (ω) for dense
colloidal dispersions have already been reported [5, 27–29],
we compare our numerical results with those experiments,
to examine the validity of our method for dense particle
dispersions.

2 Simulation method
Let us consider a dispersion composed of identical spherical particles of radius a in a Newtonian host ﬂuid subjected to oscillatory shear. The host ﬂuid is described by
the velocity ﬁeld v(r, t) and the pressure ﬁeld p(r, t). The
i-th dispersed particle is described by {Ri , Vi , Ωi }, where
Ri is the position of the particle, Vi is the translational
velocity, and Ωi is the rotational velocity. The coupling
scheme between the ﬂuid motions and the particle motions is based on the smoothed proﬁle method (SPM),
which introduces a particle density ﬁeld φ(r, t) ∈ [0 : 1]on the entire ﬁeld, φ = 1 for the particle domains and
φ = 0 for the ﬂuid domains. These domains are separated
by a thin interfacial domain of thickness ξ. The SPM is aneﬃcient method to resolve the HIs between the ﬂuid andparticle motions; its details are given in [10, 11, 13].
The time evolution of the host ﬂuid is governed by the
Navier-Stokes equation:
ρf (∂t v + v · ∇v) = ∇ · σ + ρf φfp + ρf f shear ,(1)with the incompressibility condition ∇ · v = 0, where ρfis the density of the ﬂuid, η is the shear viscosity, thestress tensor σ = −pI + η{∇v + (∇v)T }, and f shear (r, t)
is an external force ﬁeld that is introduced to enforce anoscillatory shear ﬂow on the entire system. The body force
φfp is introdued to ensure the rigidity of particles andthe appropriate boundary condition at the ﬂuid/particle
interface [10, 11, 13].
The time evolution of the i-th particle with mass Mi
and moment of inertia Ii is governed by Newton’s equations of motion:
Mi V̇i = FiH + FiC + GVi ,Ii · Ω̇i = NiH + GΩ
i ,Ṙi = Vi ,

(2)
(3)

where FiH and NiH are the hydrodynamic forces andtorques, respectively, exerted by the host ﬂuid on theparticle. FiC is the repulsive force that prevents the particles from overlapping, and a truncated Lennard-Jones
potential, V (rij ) = 4[(σ/rij )36 − (σ/rij )18 + 1/4] forrij < 21/18 σ or V (rij ) = 0, is adopted in this work.
Here, σ = 2a and rij = |Ri − Rj |. GVi and GΩ
i are
the random force and torque, respectively, due to thermal ﬂuctuations. These random ﬂuctuations are assumed
VV
to be GVi  = GΩi  = 0, Gi (t) · Gi (0) = 3αV δ(t) andΩΩGi (t) · Gi (0) = 3αΩ δ(t), where  denotes time averaging, and αV and αΩ are numerical parameters to control
the particle temperature T . The procedure for determining the temperature is described in [30].
The apparent stress σ app of the dispersion is written
as


1
d
1
(4)
dxxρt f shear −dxx (ρt v),
σ app =
VV
dt
where ρt = (1 − φ)ρf + φρp , ρp is the density of the particles, and V is the volume of the system. The derivationof eq. (4) was reported in [22].
The apparent stress σ app consists of two terms: theﬁrst term is a stress tensor including the external force
and the second term is a stress tensor arising from theacceleration of the dispersion. In experimental viscoelastic measurements of dispersions, the acceleration term canbe ignored, because the relaxation time scales related tothe acceleration are considerably smaller than the experimental time scales. On the other hand, in the simulations based on the DNS approach, the acceleration termstrongly aﬀects the apparent stress of the dispersion.
There are two key points for calculating the apparent
stress σ app of the dispersion: i) how to calculate the acceleration term in σ app , and ii) how to introduce the external
force f shear . The ﬁrst point is straightforward. For a simulation time step h, the acceleration term can be simply

T. Iwashita et al.: A direct numerical simulation method for complex modulus of particle dispersions


calculated as 1/V dxx(ρt v(x, t + h) − ρt v(x, t))/h. Next
the external force f shear is introduced as a body force, toenforce the following oscillatory velocity ﬁeld over the entire system:⎧⎪⎨γ̇(t)(−y − Ly /2), (−Ly /2 < y ≤ −Ly /4),0vx (y) = γ̇(t)y,(−Ly /4 < y ≤ Ly /4), (5)⎪⎩γ̇(t)(−y + Ly /2), (Ly /4 < y ≤ Ly /2),

f shear = (ρt (vx0 (r) − vx (r))/h, 0, 0),r ∈ V.

(7)

If the external force is introduced as a boundary force,then the development of the velocity from the boundary to bulk (i.e., the propagation modes) is observed
for short time scales, such as the kinematic time scaleτν (= ρf a2 /η). In the simulation, the propagation modes
inﬂuence the overall viscoelastic properties of the dispersion. Conversely, in experimental measurements the propagation modes are ignored. Using the external forces mentioned above, we can eliminate the propagation modes numerically.
To measure the storage modulus G′ (ω) and loss modulus G′′ (ω), we monitor the xy component of the apparent
app
shear stress, σxy, and shear rate γ̇ as a function of time.
In general, the xy component of the stress is written asapp
σxy= σ0 cos(ωt − δ),(8)where σ0 is the amplitude of the stress and δ is the phaseappand γ̇. In our model, the shear
diﬀerence between σxy
rate is an externally controlled parameter. By using theobtained σ0 and δ, we can determine the following dynamic viscoelastic moduli:
σ0 sin δω,γ̇0G′′ (ω) =σ0 cos δ
ω.
γ̇0

(9)

3 Simulation results
Three-dimensional simulations were performed at several
frequencies ranging from 0.0005 to 0.2. The amplitude of π/2ω
strain, γ0 = 0γ̇(s)ds = γ̇0 /ω, was set to be 0.2. The
spatial and temporal units are expressed in terms of thelattice spacing ∆ and ρf ∆2 /η, respectively. The system
lengths are Lx = Ly = Lz = 64. The parameters of thesimulation are a = 4, ξ = 2, ρp = 1, ρf = 1, η = 1, andkB T = 7. The dispersed particles are initially randomly
distributed.

(a)

0.01

shear stressshear rate

0

-0.01

-0.02
500

550

600

650

700

750

800

t
101

(6)

G’’
ηω
100

3

where y denotes distance in the velocity gradient direction, and γ̇(t) is the oscillatory shear rate, with amplitude γ̇0 and frequency ω. Here, the ﬂow is imposed in
the x-direction and Ly is the length of the system in they-direction. This zigzag velocity proﬁle was ﬁrst used tosimulate dispersions in a steady shear [22, 30]. Note thatf shear is deﬁned to be a body force that constrains thevelocity ﬁeld of the dispersions, and its explicit form is
written as

G′ (ω) =0.02

G’’a/kBT

γ̇(t) = γ̇0 cos(ωt),

359

(b)

10-1

10-2

10-3
10-1

100

101

102

2

ωa/D0
Fig. 1. The complex muduli of the host ﬂuid at Φ = 0:app(a) Time evolutions of the shear stress σxy
and the shear rateγ̇ at ω = 0.1. (b) Frequency dependence of the loss modulus
G′′ of the host ﬂuid. The modulus and frequency are scaled by
a3 /kB T and a2 /D0 , respectively, and D0 = kB T /6πηa. The
solid line represents the loss modulus of the host ﬂuid, ηω.

If we consider a particle of radius 0.4 µm in water atroom temperature, the units of space and time correspondto be 0.1 µm and 0.0112 µs, respectively. In this case thesimulated range of the frequency, f = ω/2π, is between
7.1 and 2842 kHz.

3.1 Test of the simulation method
In order to test the validity of our method, we applied themethod to a Newtonian host ﬂuid that does not containdispersed particles. Figure 1(a) shows the time evolutions
appof the shear rate γ̇ and shear stress σxy
at ω = 0.1. It
can be seen that both curves develop in time with thesame phase, i.e., δ = 0. This behavior represents the typical features of Newtonian ﬂuids. The loss modulus G′′
was then calculated using eq. (9) for diﬀerent frequencies. Figure 1(b) displays the frequency dependence of theloss modulus G′′ for the host ﬂuid. The modulus and frequency are non-dimensionalized by a3 /kB T and a2 /D0 ,respectively, where D0 is the self-diﬀusion coeﬃcient of a
Brownian particle at inﬁnite dilution, D0 = kB T /6πηa.
The loss modulus increases linearly with ω, and its slope
is equal to the viscosity of the host ﬂuid. The solid line

360

The European Physical Journal E
0.15

(a)

0.1

increases linearly with ω. As ω increases, G′ grows monotonically until it reaches a plateau region, while G′′ develops up to a linear region, and its slope is smaller than
the slope at low frequencies. The slope at high frequencies′, which is given
represents the high-frequency viscosity η∞byshear stressshear rate

0.05
0
-0.05
-0.1
-0.15
500

550

600

650

700

750

800

′
η∞
= lim G′′ (ω)/ω.
ω→∞

t

(10)

102

From ﬁg. 2 the high-frequency viscosity was roughly esti′≃ 5.2. Next we focus on the high-frequency
mated as η∞
elastic shear modulus G′∞ , which is deﬁned as

G’

101

∼ η’∞ω

(b)
G ’∞

3

G’a/kBT and G’’a/kBT

G’’

10

-1

1

ω

3

10

-2

10

2

ω

ωp
10-3
10-2

lim G′ (ω) = G′∞ .

ω→∞

0

10-1

100

101

102

2

ωa/D0
Fig. 2. The complex moduli of the dispersion at Φ = 0.41:
app(a) Time evolutions of the shear stress σxy
and the shear rateγ̇ at ω = 0.1. (b) Frequency dependence of the storage modulus G′ () and loss modulus G′′ () for the dispersion. The
modulus and frequency are scaled by a3 /kB T and a2 /D0 , respectively. G′∞ is the high-frequency elastic modulus, which
was obtained by the numerical integration of a theoretical ex′pression for G′∞ (dotted line) [31], and η∞is the high-frequency
viscosity. The down arrow indicates the inverse of the particle
relaxation time, ωp = τp−1 = 6D0 (Φ)/a2 .

indicates the loss modulus of the host ﬂuid, ηω. These results show that the viscoelastic properties of the host ﬂuid
are correctly reproduced by the simulation.
We next investigated a concentrated dispersion composed of spherical particles ﬂuctuating in the host ﬂuid.
The volume fraction of the dispersed particles was Φ =0.41. Figure 2(a) shows the time evolutions of the shear
app
rate γ̇ and shear stress σxy
at ω = 0.1 and Φ = 0.41.
The behavior of the shear stress of the dispersion diﬀers
from that of the host ﬂuid shown in ﬁg. 1. We can seethat there are phase diﬀerences between the shear rateand shear stress (i.e., δ = 0), and the amplitude of theshear stress becomes greater than that of the shear rate.
From the obtained δ, σ0 and γ̇0 , we can calculate thestorage modulus G′ and the loss modulus G′′ at diﬀerent frequencies. Figure 2(b) shows the frequency dependence of the storage modulus and the loss modulus forthe dispersion at Φ = 0.41. The modulus and frequency
are scaled by a3 /kB T and a2 /D0 , respectively. For low
frequencies, G′ (ω) increases linearly with ω 2 and G′′ (ω)(11)

We compared our results with a theoretical expression forG′∞ of hard spheres, which has been derived by Lionberger
and Russel [31]. To evaluate the theoretical value, we per′formed a numerical integration by using the obtained η∞and Percus-Yevick distribution function for Φ = 0.41. The
calculated high-frequency modulus, G′∞ , is shown in ﬁg. 2.
We can see that the values of G′ at high frequencies approach the theoretical value.
Furthermore, the loss modulus is larger than the storage modulus, and the dispersion behaves like a viscous
ﬂuid. Here, we deﬁne the particle relaxation time to beτp = a2 /6D0 (Φ), where D0 (Φ) is the diﬀusion coeﬃcientof a spherical particle in the dispersion at thermal equilibrium. From the equilibrium calculations (γ̇ = 0), we canestimate the diﬀusion coeﬃcient at Φ = 0.41, resulting inωp ≡ τp−1 ∼ 1.5 × 10−2 . The characteristic frequency ωp is
also indicated by an arrow in ﬁg. 2. At ω ∼ ωp (ωτp ∼ 1),
the onset of elasticity is clearly observed. These behaviors accurately represent the typical viscoelastic featuresof concentrated dispersions [5].
Finally, we examined the dispersion behavior of particles at higher concentrations, Φ = 0.46 and 0.51. The
frequency dependence of G′ and G′′ for diﬀerent volume
fractions is shown in ﬁg. 3; the data at Φ = 0 and 0.41 arealso plotted for comparison. For Φ = 0.46, the structuresof the particles are still randomly distributed under shear,
wheres for Φ = 0.51 the system forms crystalline phases
completely. For Φ = 0.46 and 0.51, each modulus shifts
to higher values than those at Φ = 0.41. At low frequencies, the Φ-dependence of both moduli is remarkable: thestorage modulus G′ rises rapidly with increasing Φ, causing the storage modulus to become larger than the loss
modulus at low frequencies and Φ = 0.51. The dispersion
then becomes elastic. Furthermore, the loss modulus G′′
at Φ = 0.51 is a concave curve, with a minimum at low
frequencies.
We found that the viscoelastic properties of the concentrated dispersions depend strongly on the volume fraction. In addition, at Φ = 0.51, the storage modulus increases very slowly with increasing ω, and at ω = 0.02,
a crossover from elastic- to viscous-dominant regions isobserved.

T. Iwashita et al.: A direct numerical simulation method for complex modulus of particle dispersions
102

101

G’a/kBT and G’’a/kBT

(a)
0

101

3

10

10-1

G’’
G’
G ’∞
G’’ 0.37 Shikata
G’ 0.37 Shikata

100

10-1

3

3

G’a/kBT

361

10-2

-2

G’ : Φ = 0. 51
G’ :
0. 46
G’ :
0. 41
10-3
10-2

10-1

100

101

10

102

3

101

102

103

Fig. 4. The frequency-dependent shear moduli of the simulations at Φ = 0.41 (open symbol) and the experimental resultsat Φ = 0.37 for 53 nm silica particles in ethylene glycol/glycol(ﬁlled symbol), which were measured by Shikata and Pearson [5], and the high frequency elastic modulus from the simulation (dotted line).

(b)

100

G’’a/kBT

100

ωa/D0

ωa/D0

10-1

G’’ : Φ = 0. 51
0. 46
0. 41
0. 0

10-2

10-3
10-2

10-1

2

2

101

10-3
10-2

10-1

100

101

102

2

ωa/D0
Fig. 3. Frequency dependence of (a) the storage modulus G′
and (b) the loss modulus G′′ at diﬀerent frequencies: Φ = 0.51
(black symbols), Φ = 0.46 (gray symbols), Φ = 0.41 (white
symbols), and Φ = 0.0 (∗). The modulus and frequency are
scaled by a3 /kB T and a2 /D0 , respectively. The amplitude ofstrain is set to 0.2.

3.2 Comparison of simulations and experiments
The simulation results obtained in this work were compared with experimental results for the dynamic viscoelasticity of colloidal dispersions, which were measured by several groups [5, 27–29]. These experiments use colloidal dispersions of spherical particles dispersed in solvent. All experimental data were also plotted as a function of dimensionless variables, G′ a3 /kB T , G′′ a3 /kB T , and ωa2 /D0 . In
ﬁg. 4, the simulation results for elastic and loss modulus
at Φ = 0.41 were plotted together with those measured
experimentally by Shikata and Pearson [5]. No ﬁtting parameters are used. The viscoelastic responses of the simulations agree well with those of the experiment, although
the volume fraction between them is not exactly the same.

This clearly conﬁrms our simulation can provide the typical viscoelastic behavior of colloidal dispersions in ﬂuid
states.
Comparison of simulation results and several experiments, including concentrated dispersions with higher Φ,
would give us comprehensive information for the understanding of the dynamical behavior of colloidal dispersions. Figure 5 shows the experimental data on the elastic
and loss modulus of concentrated dispersions over a wide
range of volume fractions from 0.37 to 0.56 [5, 27–29], andour simulation results were also plotted. These experimental data were measured over a range of ωa2 /D0 from 10−5
to 103 , and our simulation results lie within the range.
The experimental curve shifts at lower values of ωa2 /D0
as the volume fraction increases. This reﬂects the characteristic time scales of the system becoming longer, and the
dynamics of dispersed particles slows down.
The experimental data for high volume fraction vary
widely depending on each experiment, and it is considered that the viscoelastic response becomes sensitive tothe details of direct interactions between particles. The
repulsive potential used in this study can produce the features of dispersions with relatively small volume fraction(see ﬁg. 4). The simulation at Φ = 0.51 has a crystalline
structure over the whole frequency range, where a colloidal crystal is formed. Thus the simulation results at
Φ = 0.51 are completely diﬀerent from the experiments,
for which the structures are in ﬂuid and glasslike states.
Recently Crassous et al. [27] have examined the eﬀect of
crystallization on the viscoelasticity of the dispersions in
the low-frequency regime in the vicinity of the glass transition, and they found the elastic modulus is greater thanthe loss modulus at low frequencies even in ﬂuidlike states.

362

The European Physical Journal E
102

100

3

G’a/kBT101

0.51 Simulation
0.46 Simulation
0.41 Simulation
0.37 Shikata
0.48 Shikata
0.51 Mason
0.52 Mason
0.53 Mason
0.54 Crassous
0.56 Crassous
0.40 Watanabe
0.53 Watanabe

10-1

-2

10

10-3
10-4

10-3

10-2

10-1

100

101

102

103

2

ωa/D0
102

100

3

G’’a/kBT

101

10-1

0.51 Simulation
0.46 Simulation
0.41 Simulation
0.37 Shikata
0.51 Mason
0.52 Mason
0.53 Mason
0.54 Crassous
0.56 Crassous
0.40 Watanabe
0.53 Watanabe

-2

10

10-3
10-5

10-4

10-3

10-2

10-1

100

101

102

under several non-equilibrium conditions, such as electrophoresis of charged spherical particles under external
electric ﬁelds [20]. The method was then modiﬁed to introduce thermal ﬂuctuations into the dispersions [12], so that
one can also simulate situations where the thermal ﬂuctuations and hydrodynamic interactions acting among dispersed particles are both important. We have carried out
systematic simulations for dispersions composed of spherical particles with and without steady shear ﬂow, in orderto analyze the diﬀusion process of dispersed particles indetail [30] and to investigate the nonlinear viscosity of thesystem [22].
In the present paper, we report an important extensionof the SPM for analyzing the viscoelastic properties of particle dispersions immersed in host ﬂuids, by introducing atemporally oscillatory external force into the system. To
be consistent with the usual periodic boundary conditions,a shear ﬂow with zigzag proﬁle was employed. The validity of the method was examined by evaluating the storage G′ (ω) and loss G′′ (ω) moduli of a system comprisingidentical spherical particles dispersed in an incompressible Newtonian host ﬂuid at volume fractions of Φ = 0,0.41, 0.46, and 0.51, for ﬂow rates spanning a range of frequencies, 0.005 ≤ ω ≤ 0.2. We conﬁrmed that the method
could successfully reproduce 1) purely viscous responses
for Φ = 0 and 2) typical viscoelastic responses for Φ = 0.41
and 0.46, in excellent agreement with experimental dataobtained for colloidal dispersions.
To our knowledge, the present study is the ﬁrst successful attempt to calculate the complex modulus of particle dispersions using DNS-type methods. Further applications of our DNS method to more complex systems, suchas dispersions of non-spherical particles, particles in polymer matrices, or dispersions of aggregating particles, arepromising.

103

2

ωa/D0
Fig. 5. Data for the elastic and loss modulus of colloidal
dispersions in an oscillatory shear ﬂow. The experimental
data were collected from several papers (Shikata and Pearson [5]; 53 nm silica sphere in ethylene glycol/glycol, Mason [28]; 210 nm silica sphere in ethylene glycol, Crassous [27];thermosensitive particles in water, and Watanabe [29]; 46 nm
silica sphere in ethylene glycol/glycol). The numbers indicate
the volume fraction of particles.

This behavior is qualitatively similar to the simulation results at Φ = 0.51, and it represents solidlike responses. In
the present work the dispersions in glassy states were notstudied.

4 Conclusion
We have developed a DNS method for simulating dynamics of solid particles dispersed in simple and complex
ﬂuids [10, 11]. This method, called SPM, has also beensuccessfully applied to simulate properties of dispersions

References
1. W.B. Russel, D.A. Saville, W.R. Schowalter, Colloidal Dispersions (Cambridge University Press, Cambridge, UK,1989).
2. R.G. Larson, The Structure and Rheology of Complex Fluids (Oxford University Press, Oxford, 1999).
3. H.M. Wyss, K. Miyazaki, J. Mattsson, Z. Hu, D.R. Reichman, D.A. Weitz, Phys. Rev. Lett. 98, 238303 (2007).
4. V. Carrier, G. Petekidis, J. Rheol. 53, 245 (2009).
5. T. Shikata, D.S. Pearson, J. Rheol. 38, 601 (1994).
6. J.F. Brady, G. Bossis, Annu. Rev. Fluid Mech. 20, 111(1988).
7. H.H. Hu, N.A. Patankar, M.Y. Zhu, J. Comput. Phys. 192,
427 (2001).
8. H. Tanaka, T. Araki, Phys. Rev. Lett. 85, 1338 (2000).
9. T. Kajishima, S. Takiguchi, H. Hamasaki, Y. Miyake,
JSME Int. J., Ser. B 44, 526 (2001).
10. Y. Nakayama, R. Yamamoto, Phys. Rev. E 71, 036707(2005).
11. Y. Nakayama, K. Kim, R. Yamamoto, Eur. Phys. J. E 26,
361 (2008).
12. T. Iwashita, Y. Nakayama, R. Yamamoto, J. Phys. Soc.
Jpn. 77, 074007 (2008).

T. Iwashita et al.: A direct numerical simulation method for complex modulus of particle dispersions
13. X. Lio, M.R. Maxey, G.E. Karniadakis, J. Comput. Phys.
228, 1750 (2009).
14. P.J. Atzberger, P.R. Kramer, C.S. Peskin, J. Comput.
Phys. 224, 1255 (2007).
15. M. Fujita, Y. Yamaguchi, Phys. Rev. E 77, 026706 (2008).
16. A.J.C. Ladd, Phys. Rev. Lett. 70, 1339 (1993).
17. A. Malevanets, R. Kapral, J. Chem. Phys. 112, 7260
(2000).
18. J.T. Padding, A.A. Louis, Phys. Rev. E 74, 031402 (2006).
19. R. Yamamoto, K. Kim, Y. Nakayama, K. Miyazaki, D.R.
Reichman, J. Phys. Soc. Jpn. 77, 084804 (2008).
20. K. Kim, Y. Nakayama, R. Yamamoto, Phys. Rev. Lett.
96, 208302 (2006).
21. T. Araki, H. Tanaka, EPL 82, 18004 (2008).
22. T. Iwashita, R. Yamamoto, Phys. Rev. E 80, 061402(2009).

363

23. A. Shakib-Manesh, P. Raiskinmaki, A. Koponen, M.
Kataja, J. Timonen, J. Stat. Phys. 107, 67 (2002).
24. J. Kromkamp, D.T.M. van den Ende, D. Kandhai. R.G.M.
van der Sman, R.M. Boom, J. Fluid Mech. 529, 253 (2005).
25. P.M. Kulkarmi, J.F. Morris, Phys. Fluids 20, 040602(2008).
26. M. Hecht, J. Harting, H.J. Herrmann, Phys. Rev. E 74,
021403 (2006).
27. J.J. Crassous, M. Siebenbürger, M. Ballauﬀ, M. Drechsler,
D. Hajnal, O. Henrich, M. Fuchs, J. Chem. Phys. 128,204902 (2008).
28. T.G. Mason, D.A. Weitz, Phys. Rev. Lett. 75, 2770 (1995).
29. H. Watanabe, M.L. Yao, T. Shikata, H. Niwa, Y. Morishima, Rheol. Acta 38, 2 (1999).
30. T. Iwashita, R. Yamamoto, Phys. Rev. E 79, 031401(2009).
31. R.A. Lionberger, W.B. Russel, J. Rheol. 38, 1885 (1994).