PHYSICAL REVIEW E 80, 061402 共2009兲

Direct numerical simulations for non-Newtonian rheology of concentrated particle dispersions
Takuya Iwashita* and Ryoichi Yamamoto†Department of Chemical Engineering, Kyoto University, Kyoto 615-8510, Japan and CREST, Japan Science and Technology Agency,
Kawaguchi 332-0012, Japan
共Received 1 May 2009; revised manuscript received 15 August 2009; published 11 December 2009兲
The non-Newtonian behavior of a monodisperse concentrated dispersion of spherical particles was investigated using a direct numerical simulation method, which takes into account hydrodynamic interactions andthermal fluctuations accurately. Simulations were performed under steady shear flow with periodic boundary
conditions in the three directions. The apparent shear viscosity of the dispersions was calculated at volume
fractions ranging from 0.31 to 0.56. Shear-thinning behavior was clearly observed at high volume fractions.
The low- and high-limiting viscosities were then estimated from the apparent viscosity by fitting these data intoa semiempirical formula. Furthermore, the short-time motions were examined for Brownian particles fluctuating in concentrated dispersions, for which the fluid inertia plays an important role. The mean square displacement was monitored in the vorticity direction at several different Peclet numbers and volume fractions sothat the particle diffusion coefficient is determined from the long-time behavior of the mean square displacement. Finally, the relationship between the non-Newtonian viscosity of the dispersions and the structural
relaxation of the dispersed Brownian particles is examined.
DOI: 10.1103/PhysRevE.80.061402

PACS number共s兲: 82.70.⫺y, 82.20.Wt, 47.50.⫺d, 83.50.Ax

I. INTRODUCTION

The links between macroscopic rheological properties andmicrostructures in colloidal dispersions have been extensively investigated for many systems, including dispersionsof sterically stabilized and charged-stabilized particles inhost fluids 关1,2兴.
The behavior of a monodisperse dispersion composed ofsolid particles immersed in a Newtonian host fluid strongly
depends on the volume fraction of the dispersed particles ⌽and the shear rate ␥˙ . When the shear rate is zero 共␥˙ → 0兲, theshear viscosity of the dispersion is referred to as the zeroshear viscosity ␩0. In the dilute limit 共⌽ Ⰶ 1兲, the zero-shear
viscosity is well approximated by Einstein’s formula 关3兴,

␩0 = ␩共1 + 2.5⌽兲,共1兲

where ␩ is the shear viscosity of the host fluid. In a concentrated dispersion, theoretical difficulties become rather severe, since the behavior of the dispersed particles is complicated by the interactions between the particles and thermal
fluctuations. In particular, the solvent-mediated many-bodyhydrodynamic interactions 共HI兲 between the particles complicate the dynamical behavior. A number of experiments forconcentrated dispersions have been performed to reveal theorigin of the non-Newtonian behavior of these dispersions,and several semiempirical formulas for ␩0 have been proposed to characterize the experimental results. For example,
the Krieger-Dougherty relationship 关4兴,冉

␩0 = ␩ 1 −冊⌽ −关␩兴⌽m,⌽m

共2兲

where 关␩兴 is the intrinsic viscosity and ⌽m is the packing
volume fraction at which the viscosity diverges, is often used*iwashita@cheme.kyoto-u.ac.jp†

ryoichi@cheme.kyoto-u.ac.jp1539-3755/2009/80共6兲/061402共7兲

for fitting the experimental data for uniform colloidal spheres
suspended in nonaqueous media.
When a dispersion is subjected to shear 共␥˙ ⫽ 0兲, the flow
properties of the dispersion show a variety of non-Newtonian
behaviors, such as shear thinning and shear thickening.
These non-Newtonian behaviors are associated with thechanging microstructures of the dispersion. Several physicalmechanisms for these peculiar behaviors have been proposed; for example, shear-induced order-disorder transitions关5–7兴, formations of dynamic clusters of the particles 关8,9兴.
However, a full understanding of the relationships between
the rheological properties and the microstructure has not yet
been obtained, despite extensive studies.
Computer simulations are very powerful tools in the direct investigation of the dynamics of individual particles inconcentrated dispersions. The Stokesian dynamics 共SD兲
method 关10兴 has been widely used to measure the rheology
of dispersions and provides valuable information regarding
the non-Newtonian behavior of flowing dispersions 关14–17兴.
The SD method, however, is based on the Stokes approximation, which assumes that all relaxation times associated withfluid motions are short as compared with those of the particle, i.e., ␶␯ Ⰶ ␶B, where ␶␯ = ␳ f a2 / ␩ and ␶B = 2␳ pa2 / 9␩. Here,␳ f and ␳ p are the density of the fluid and the particle, respectively, and a is the radius of the particle. With this assumption, the HI is treated as the Ronte-Prager-Yamakawa 共RPY兲
tensor and the lubrication correction. Furthermore, it is assumed that the relaxation time associated with the particle’s
inertia is zero 共␶B → 0兲. These approximations are valid forthe motion of the particles at time scales much greater than
the time scales of the relaxation of both the fluid and theparticle inertia. Therefore, the short-time motion of Brownian particles over the kinematic time scale ␶␯ cannot be described by simulation methods based on the Stokesian approximation. On the kinematic time scale, the dynamic
coupling between the fluid motion and the particle motion
remains strong.

061402-1©2009 The American Physical Society

PHYSICAL REVIEW E 80, 061402 共2009兲

TAKUYA IWASHITA AND RYOICHI YAMAMOTO

Sheared dispersions have a nondimensional parameterthat includes the particle’s inertia, e.g., the particle Reynolds
number Rep = ␳ p␥˙ a2 / ␩ = 9␶B / 2␶␥˙ , where ␶␥˙ 共=1 / ␥˙ 兲 represents
the characteristic time due to shear. For most colloidal dispersions, the particle Reynolds number is very small. The SD
method has a particle Reynolds number of zero and, therefore, cannot be applied to problems with finite particle Reynolds numbers. Typical examples include the motions of dispersions subjected to strong shear or composed of large
particles; in these cases, the particle Reynolds number has arelatively large finite value.
Another problem that cannot be treated by the SD method
is the short-time motion of the particle, in which the fluid
inertia becomes significant even if the particle Reynolds
number is very small. For example, the characteristic timescale is ␶␯ ⬃ 1 ␮s for a neutrally buoyant particle of 1 ␮m
in water, and Rep ⬃ 10−8 for ␥˙ = 0.01. The effects of the fluid
inertia appear as memory effects of the particles. For a complete understanding of these flowing dispersions, full timedependent HIs are required.
In recent years, several numerical methods have been developed in order to accurately simulate dispersions in a variety of situations, including those described above. Thesemethods of dispersion modeling are based on the same approach, which involves resolving the fluid motion simultaneously with the particle motion. We refer to this approach as
the direct numerical simulation 共DNS兲 approach. This approach enables us to accurately treat the full time-dependent
HI. The numerical methods differ mainly in the approach
used to resolve the HIs between the fluid and particle motions 关18–26兴.
In this work, we apply a direct numerical scheme basedon the smoothed profile method 共SPM兲 关26–30兴 to a monodisperse concentrated dispersion of repulsive and neutrally
buoyant Brownian particles in a shear flow. In the SPM, theNavier-Stokes equation for the fluid motion is discretized ona regular grid, and the Newtonian equations for the particle
motion are solved simultaneously with the fluid motion. Thisdeveloped scheme accounts for thermal fluctuations, a shear
flow, and memory effects. The SPM can reproduce the correct short-time behavior of a single Brownian particle in ashear flow on the kinematic time scale, and the numericaltests for a single Brownian particle in a shear flow werereported in 关30兴.
Although some groups have studied flowing dispersions
composed of repulsive spherical particles by using a DNS
method based on the lattice Boltzmann method 关31–33兴, all
the simulations have been applied to a dispersion of nonBrownian particles and ignore the thermal fluctuations of theparticles. Systematic analyses have not yet fully carried out
for concentrated dispersions of Brownian particles in a shear
flow at finite Reynolds numbers.
In the present study, we examine the non-Newtonian rheology of concentrated dispersions in shear flow at finite Reynolds and Peclet numbers using a DNS method that takes
into account hydrodynamic interactions and thermal fluctuations accurately. Note that the simulated situations are different from those of the prevailing experiments for colloidal
dispersions, since even the lowest particle Reynolds numberin these simulations is several hundred times larger thanthose of the experiments. We first present the simulation
method and the manner in which the apparent shear viscosityof the dispersions is calculated. Three-dimensional simulations are performed with periodic boundary conditions, andthe non-Newtonian behavior of the shear viscosity is obtained for several volume fractions. Both the high and low
shear limiting viscosities obtained from the simulations arethen compared with the Krieger-Doughty relationship. Moreover, we investigate the short-time motions of Brownian particles in a sheared concentrated dispersion on the kinematictime scales. The mean square displacement 共MSD兲 of
Brownian particles is monitored at several Peclet numbers
and volume fractions, and the long-time diffusion coefficientis determined from the long-time behavior of the MSD. We
also suggest a simple relationship between the nonNewtonian viscosity of the dispersions and the structural relaxation time of the dispersed Brownian particles.
II. SIMULATION METHOD

A direct numerical scheme that implements both a shear
flow and thermal fluctuations is briefly explained. A moredetailed explanation is given in a previous publication 关30兴.
We consider a monodisperse dispersion of N p repulsive
spherical particles of diameter ␴ in a Newtonian host fluid.
The particles interact via a truncated Lennard-Jones 共LJ兲 potential,
ULJ共r兲 =冉冊 冉冊 册
冋
冦
4⑀

␴ 36
␴ 18 1
−
+
4
r
r

0

共r ⱕ 21/18␴兲,
共r ⬎ 21/18␴兲,冧

共3兲

where r is the distance between two particles and the parameter ⑀ characterizes the interaction strength. The position of
the ith dispersed particle is Ri, the translational velocity is Vi,
and the rotational velocity is ⍀i. The time evolution of theith particle with mass M i and moment of inertia Ii is governed by Newton’s equations of motion,CV
M iV̇i = FH
i + Fi + Gi ,Ṙi = V̇i ,⍀Ii · ⍀̇i = NHi + Gi ,

共4兲
共5兲

H
where FH
i and Ni are the hydrodynamic forces and torques
exerted by the host fluid on the particle. FCi is a repulsive

force arising from the potential of Eq. 共3兲, which prevents
the particles from overlapping. GVi and G⍀
i are random forces
and torques, respectively, due to thermal fluctuations. These
random fluctuations are assumed to be Markovian 共white ortime delta correlated兲 and determine the particle temperature
T. The procedure for determining the temperature is described in 关28–30兴.
In the SP method, the velocity and pressure fields, v共x , t兲and p共x , t兲, are defined on three-dimensional Cartesian grids,
which consist of fluid and particle domains. In order to distinguish the particle and fluid domains on the grids, asmoothed function ␾共x , t兲, which is equal to 1 in the particle
domains and 0 in the fluid domains, is introduced. These
domains are separated by a thin interfacial domain of thick-

061402-2

PHYSICAL REVIEW E 80, 061402 共2009兲

DIRECT NUMERICAL SIMULATIONS FOR NON-…

ness ␰. The system size is 关0 , Lx兴 ⫻ 关−Ly / 2 , Ly / 2兴 ⫻ 关0 , Lz兴.
The time evolution of the velocity field is governed by the
Navier-Stokes equation with the incompressibility condition
⵱ · v = 0,

␳ f 共⳵tv + v · ⵱v兲 = ⵱ · ␴ + ␳ f ␾f p + ␳ f fshear ,共6兲

where the stress tensor ␴ = −pI + ␩兵⵱v + 共⵱v兲T其, andfshear共x , t兲 is an external force field that is introduced to enforce a simple shear flow on the system. The flow is imposedin the x direction and the external force is introduced as aconstraint force so that the velocity field satisfies vx共y兲 =−␥˙ Ly / 4 at y = −Ly / 4 and vx共y兲 = ␥˙ Ly / 4 at y = Ly / 4, where ydenotes the distance in the velocity gradient direction. A
simple shear flow with a shear rate of ␥˙ is then approximately produced over a range from y = −Ly / 4 to y = Ly / 4. ␾f p
represents a body force that ensures the rigidity of the particles and the appropriate nonslip boundary conditions at thefluid/particle interface, which is further elaborated upon in
Refs. 关26,27兴.
The unit of length is taken to be the lattice spacing ⌬, andthe unit of time is ␶0 = ␳ f ⌬2 / ␩. Unless otherwise stated, we
set ⌬ = 1, ␶0 = 1, ␩ = 1, ␳ f = 1, ␳ p = 1, a = 4, ⑀ = 1, and ␴ = 8.
Assuming dispersions of neutrally buoyant particles of radius
1 ␮m in water at room temperature, our unit length ⌬ andtime ␶0 correspond to 0.25 ␮m and 0.0625 ␮s, respectively
The second-order Runge-Kutta algorithm is used to integrate
the Newtonian equations. The Navier-Stokes equation is discretized with a Fourier spectral scheme in space and with asecond-order Runge-Kutta scheme in time. The discretized
time step is h = 0.07145. Although this value is chosen fromthe stability condition of the Navier-Stokes equation, it canbe used safely for the particle’s equations motion because his much smaller than the Lennard-Jones time unit ␶ M= 共M i␴2 / ⑀兲1/2 ⯝ 131. The phase diagram of the present
共36:18兲 LJ system depends weakly on the system temperature T as well as the volume fraction ⌽ of the particles. To
avoid the delicate issue of crystallization, all the simulationsin the present paper were carried out at T and ⌽ with which
the system is in amorphous states.
To measure the rheological properties of the particle dispersion in shear flow, we can calculate the apparent stress␴app of the dispersions in the following manner. The momentum equation for the dispersion is formally written asd
共␳tv兲 = ⵱ · ␴dis + ␳tfshear ,dt

共7兲␳t = 共1 − ␾兲␳ f + ␾␳ p ,共8兲

where ␴dis denotes the stress tensor of the dispersion including the inertia, the pressure and the viscous terms. The fullstress tensor s of the flowing dispersion is then defined byintroducing a convective momentum-flux tensor explicitly ass = ␴dis − ␳tvv ,共9兲

where ␳tvv represents momentum transport by the bulk flowof the dispersion.
Although s cannot be calculated directly, we can obtain
the apparent stress ␴app by using the local stress s:

冕
冕
冕
冕 冋冉
冕

␴app =

1
V

dxs,

共10兲

=

1
V

dx关关⵱ · 共sx兲兴T − x ⵱ · s兴

共11兲

=

1
V

dx关− x ⵱ · s兴

共12兲

=

1
V

=

1
V

dx x ␳tfshear −dxx␳tfshear −

1
V

冊册

共13兲

⳵
共␳tv兲,
⳵t

共14兲⳵
共 ␳ tv 兲
⳵t

冕

dxx

with a volume V = LxLyLz. In the derivation of Eq. 共11兲, we
use an second rank identity, s = 关⵱ · 共sx兲兴T − x ⵱ · s. In thesteady state,冓 冔⳵共 ␳ tv 兲
⳵t

共15兲

= 0,
t

where 具 典t denotes time averaging over the steady state. The
time-averaged apparent shear stress of the dispersion canthen be written as
具␴app典t =1
V冓冕

dxx␳tfshear

冔

.

共16兲

t

Then one can obtain the apparent shear stress under steadyshear flow from the external force fshear imposed in the
Navire-Stokes equation.
III. RESULTS AND DISCUSSION

Simulations were performed in a three-dimensional cubic
box, whose side length is L, with periodic boundary conditions. Most present simulations were done with L = 64, where
the number of particles are N p = 300, 400, 450, 500, and 550
for ⌽ = 0.31, 0.41, 0.46, 0.51, and 0.56, respectively. The
index of axis x, y, and z represent the flow, velocity gradient,and vorticity directions, respectively. For a spherical particle,
a = 4, ␰ = 2. The temperature was determined by equilibrium
calculations before the shear flow was imposed. The temperature is kBT = 7. The initial configuration of the particles is
set to be a random distribution. A large number of simulations were performed for volume fractions of 0.31ⱕ ⌽ⱕ 0.56 and shear rates of 5 ⫻ 10−5 ⱕ ␥˙ ⱕ 0.1. These systems
have particle Reynolds numbers Rep ranging from 8 ⫻ 10−4
to 1.6.
For dispersions composed of Brownian particles in thesteady state, the apparent shear viscosity of the dispersion isdefined as␩app =

具␴app
xy 典t
,
␥˙

共17兲

app
. To examine
where ␴app
xy denotes the xy components of ␴
the system size effects in the present simulations, we calcu-

061402-3

PHYSICAL REVIEW E 80, 061402 共2009兲

TAKUYA IWASHITA AND RYOICHI YAMAMOTO4040

Φ = 0.56
0.51

35

35

0.46
0.41

30

25

ηο, ηοοη app

Einstein formlula
Krieger-Dougherty : Φm = 0.63
Krieger-Dougherty : Φm = 0.73

30

0.31
Solvent (η=1)

25
20

20

15

15

10

10

5

5

0

Low limiting viscosity η0
High limitinig viscosity ηοο

-2

10

-1

0

10

Pe

0

1

10

10

lated the apparent shear viscosity for three different systemsizes for a constant ⌽ = 0.41, i.e., L = 32 with N p = 50, L= 128 with N p = 3200, and L = 256 with N p = 25 600. We foundonly negligibly small differences among the values of theapparent shear viscosity calculated for the three systems
within a statistical error.
Figure 1 shows the dependence of the apparent shear viscosity on the Peclet number for several volume fractions.
Here, the Peclet number is defined as Pe= 6␲␩a3␥˙ / kBT. For
the lowest concentration 共⌽ = 0.31兲, the viscosities remain
nearly constant and exhibit Newtonian behavior. For higher
concentrations 共⌽ ⱖ 0.41兲, the dispersions show nonNewtonian behavior. As Pe increases, the shear-thinning behavior is clearly observed from the higher plateau region for
Pe of order 10−2 to a lower plateau region for Pe of about 10.
From both of the plateau values of the viscosity curve, we
can obtain the low shear limiting viscosity 共identified as ␩0兲and the high shear limiting viscosity ␩⬁ for each volume
fraction.
In order to evaluate ␩0 and ␩⬁, we fit our simulation data
into the following simple empirical function of Pe and ⌽,

冉

0.1

0.2

0.3

0.4

0.5

0.6

Volume fraction

FIG. 1. The Peclet number Pe= 6␲␩a3␥˙ / kBT dependence of theapparent shear viscosity of the dispersion for several volume fractions. The solid lines represents ␩ f = ␩⬁ + 共␩0 − ␩⬁兲 / 关1 + b−1共⌽兲Pe兴
for ⌽ ⱖ 0.46, where ␩0 is the low shear limiting viscosity, ␩⬁ is thehigh shear limiting viscosity, and b共⌽兲 is a fitting parameter, withb = 0.20, 0.50, and 0.63 for ⌽ = 0.56, 0.51, and 0.46, respectively.
The dashed line represents the shear viscosity of the host fluid.

␩ f 共Pe,⌽兲 = ␩⬁ +0冊

␩0 − ␩⬁,1 + b−1共⌽兲Pe

共18兲

where b共⌽兲 is a fitting parameter. This empirical function isplotted in Fig. 1 with solid lines for ⌽ ⱖ 0.46. One finds thatit coincides with the simulation data reasonably well. It isalso seen that the onset of shear thinning appears at smaller
Peclet numbers with increasing volume fraction. We note

FIG. 2. 共Color online兲 Volume fraction dependence of the low
shear limiting viscosity ␩0 and the high shear limiting viscosity ␩⬁of the dispersions. The solid line represents the Krieger-Dougherty
relationship with ⌽m = 0.63 for ␩0 and the dashed line ⌽m = 0.73 for␩⬁. The dotted line represents Einsteinfs formula of Eq. 共1兲.

that the particles are randomly distributed in the dispersions
in the present simulations, i.e., all the viscosity data are taken
in situations at which no shear-induced crystallization occurs.
Figure 2 shows the dependence of the low shear limiting
viscosity ␩0 and the high shear limiting viscosity ␩⬁ on thevolume fraction. Both shear limiting viscosities increase
monotonically with the volume fraction. The simulation results for both of the shear limiting viscosities agree well withthe semiempirical relations of Eq. 共2兲 by Krieger-Dougherty,␩0,⬁ = ␩共1 − ⌽ / ⌽m兲−2.5⌽m where ⌽m = 0.63 for ␩0 and ⌽m= 0.73 for ␩⬁. The dotted line represents Einstein’ s formulaof Eq. 共1兲.
The fluid and particle inertia contribution to the apparent
shear stress can be written as=␴inertia
xy

1
V

冕

dx␳tvxvy ,共19兲

which is similar to the Reynolds stress of a uniform fluid inturbulence. This stress represents the strength of the hydrodynamic instability. The inertia contribution to the viscosity,␩inertia = 具␴inertia
典t / ␥˙ , is about two order of magnitude smallerxy
than the apparent shear viscosity. We conclude that the inertia contribution is negligible in the range 8 ⫻ 10−4 ⱕ Rep
ⱕ 1.6.
We next examine the dynamical motion of Brownian particles in a shear flow. We analyze the mean square displacement in the vorticity direction 共z兲 for the Brownian particles,

061402-4N

1 p
具关⌬R 共t兲兴 典 = 兺 具兩Rzj 共t兲 − Rzj 共0兲兩2典.
N p j=1z

2

共20兲

PHYSICAL REVIEW E 80, 061402 共2009兲

DIRECT NUMERICAL SIMULATIONS FOR NON-…
3

3

10

10

Φ = 0.56
0.51
0.46
0.41
0.31

102

Φ = 0.56
0.51
0.46
0.41
0.31

102

analytical solution

analytical solution

(a)
1

(b)

1

10

<∆Rz2 ><∆Rz2 >

10

10

0

10

10-1

10-1

0

τν
-2

10

100

101

τD
102

τν
103

104

t

(a)

-2

10

0

τD

1

10

10

10

(b)

2

10

t

3

10

4

FIG. 3. The mean square displacement in the vorticity direction for several volume fractions 共a兲 at the lowest Peclet number 共Pe= 0.0086 at ⌽ = 0.56 and Pe= 0.022 for ⌽ ⱕ 0.46兲 and 共b兲 at the highest Peclet number 共Pe= 17.2兲. The solid line represents the analyticalsolution for the generalized Langevin equation of a free Brownian particle 关34兴. The arrows indicate the kinematic time ␶␯ = ␳ f a2 / ␩ = 16 andthe diffusion time ␶D = a2 / D0 ⬃ 172.

1
具关⌬Rz共t兲兴2典
t→⬁ 2tDzsim = lim

ever, confirmed that the diffusion coefficients behave as
Dzsim共⌽,L−1兲 = Dz共⌽兲 − ␰kBT/6␲␩L

1.2

1

0.8

0.6

0.4

Pe = 17.23
13.78
1.723

0.2

0.172
0.043

共21兲

and examined the system size effects in Dzsim. In contrast to
the negligibly small system size effects observed in the apparent shear viscosity ␩app, Dzsim shows notable increase withincreasing system size for all volume fractions. It is, how-

共22兲with ␰ ⬃ 2.83, where Dz共⌽兲 represents the value of diffusion
coefficient extrapolated for L → ⬁. The system size effects
observed in the present simulations are essentially the sameas those reported in earlier papers 关11–13兴. We calculated
Dz共⌽兲 from extrapolation of our simulation data Dz共⌽ , L兲with finite L. Figure 4 shows the volume fraction and theDz /D0Figure 3共a兲 shows the MSD of the dispersion with thelowest Peclet number for each volume fraction. These results
are considered to be the MSD at thermal equilibrium. As thevolume fraction is increased, the dynamical behavior of the
MSD varies greatly, as compared with that of the hydrodynamic analytical solution for a single Brownian particle in ashear flow 关34兴. The analytical solution is derived from thegeneralized Langevin equation with memory effects. With
increasing volume fraction, MSD increases more slowly withtime. For ⌽ = 0.56, we found that a plateau region starts toappear around a time scale of order 102, which is slightly
greater than the kinematic time ␶␯. This is a typical behavior
of colloidal dispersions in glassy states, e.g., colloid glasses.
Figure 3共b兲 shows the MSD of the dispersion at the highest Peclet number for each volume fraction. Since Pe⯝ 17.2 here, the shear force is much stronger than the thermal force. The volume fraction dependence of the MSD inthis figure is opposite to the previous case at the lowest Peclet number shown in Fig. 3共a兲. At short times, the MSD
grows more rapidly in time with increasing volume fraction
because HIs between particles are more enhanced at higher
volume fraction. At later times, the MSDs tend to exhibit
diffusive motions with a volume fraction independent diffusion coefficient.
We then calculated the long-time diffusion coefficientDzsim in the vorticity direction via

0
0.25

0.3

0.35

0.4

0.45

0.5

0.55

0.6

Volume fraction
FIG. 4. 共Color online兲 The long-time diffusion coefficient in thevorticity direction for different Peclet numbers and volume fractions. The diffusion coefficient is scaled by D0 = kBT / 6␲␩a.

061402-5

PHYSICAL REVIEW E 80, 061402 共2009兲

TAKUYA IWASHITA AND RYOICHI YAMAMOTO
1.2

Φ = 0.56
0.51

1(ηapp - ηοο)/ (ηο - ηοο )and A is a fitting parameter. The shear thinning starts ataround ␶ p␥˙ ⯝ 1. This behavior is very similar to the nonNewton rheology of supercooled liquids 关35兴.
We confirmed that the present numerical method of introducing the thermal fluctuation successfully reproduces thefluctuation dissipation theorem for time scales longer than
the so-called the Brownian time 关28,29兴. Although it was
confirmed also that the present method works quite well forthe volume fractions 0 ⱕ ⌽ ⱕ 0.56 considered in the present
study, further careful tests must be needed for highly concentrated dispersions, where the fluctuation of a tagged particle
tends to correlate with motions of surrounding particles.

0.46
0.8

0.6

0.4

0.2

IV. CONCLUSIONS
0

-0.2

10-2

10-1

100 .

τp γ

101

102

FIG. 5. Scaled plot of 共␩app − ␩⬁兲 / 共␩0 − ␩⬁兲 against reduced Peclet number ␶ p␥˙ for ⌽ ⱖ 0.46. Line is 1 / 共1 + A␶ p␥˙ 兲 where A = 0.57.

Peclet number dependencies of Dz共⌽兲 normalized by D0= kBT / 6␲␩a. For the low Peclet numbers, the diffusion coefficient decreases with increasing volume fraction, and thisdependence is similar to that of the dispersion at thermal
equilibrium. On the other hand, for the highest Peclet number, the volume fraction dependence of Dz is almost constant.
We can see that the diffusion coefficients at high volume
fractions are strongly affected by shear. The dynamics of the
Brownian particles differs considerably, depending onwhether the thermal or shear force is dominant. The volume
fraction dependence is in qualitative agreement with the numerical results achieved based on the SD method, which isvalid for long-time scales 关14兴.
Finally, the relationship between the non-Newtonian viscosity of the dispersions and the structural relaxation of thedispersed Brownian particles is examined. The structural relaxation time is defined as␶ p共⌽兲 = a2/D0共⌽兲,

共23兲

where D0共⌽兲 is the diffusion coefficient for each volume
fraction at zero shear. This relaxation time represents a timeneeded for a particle to diffuse away a distance comparable
to its radius. Figure 5 shows 共␩app − ␩⬁兲 / 共␩0 − ␩⬁兲 versus ␶ p␥˙for ⌽ ⱖ 0.46. Here, we used Dz共⌽兲 at the lowest Peclet number as D0共⌽兲, and ␶ p␥˙ can be understood as a reduced Peclet
number Pe共⌽兲 = ␥˙ a2 / D0共⌽兲. One can see that the data lies ona single master curve 1 / 共1 + Ax兲 fairly well, where x = ␶ p␥˙

We investigated the rheological properties of monodisperse concentrated dispersions with repulsive spherical particles by using a DNS method that accounted for a shear flow
and thermal fluctuations. Three-dimensional simulations
were performed at Peclet numbers ranging from 0.043 to17.2 and at particle Reynolds numbers from 8 ⫻ 10−4 to 1.6.
The apparent viscosity for ⌽ = 0.31 is almost constant overthe Peclet number change. For ⌽ ⱖ 0.41, the viscosities decrease from the plateau region at the low Peclet number tothe plateau region at the high Peclet number. From the viscosity versus Peclet number curves, we can obtain both the
low and high-limiting viscosities, and these results are ingood agreement with the Krieger-Doughty relationship. The
inertia contribution to the apparent shear viscosity is verysmall throughout the entire range of Peclet numbers and volume fractions examined in the present study. As the volume
fraction increases, the behavior of the MSD in the vorticity
direction deviates from the analytical solution for a single
Brownian particle in a shear flow. For the lowest Peclet number, the MSD develops more slowly in time with increasing
volume fraction. At ⌽ = 0.56, an onset of glassy dynamics
was observed. On the other hand, for the highest Peclet number, the MSD develops more rapidly at short times with increasing volume fraction. Finally, the volume fraction dependence disappears at long times. The diffusion coefficient wascalculated from the long-time behavior of the MSD in thevorticity direction for different Peclet numbers and volume
fractions. For the lowest Peclet number, the diffusion coefficient decreases with increasing volume fraction, while it isalmost constant over the volume fraction change for thehighest Peclet number. We estimated the structural relaxation
time ␶ p from the diffusion coefficient at the lowest Peclet
number. The present non-Newtonian viscosity data agree
well with a simple scaling function 共␩app − ␩0兲 / 共␩⬁ − ␩0兲= 1 / 共1 + A␶ p␥˙ 兲 for ⌽ ⱖ 0.46, similar to the non-Newtonianviscosity of a model supercooled liquid 关35兴.

061402-6

PHYSICAL REVIEW E 80, 061402 共2009兲

DIRECT NUMERICAL SIMULATIONS FOR NON-…
关1兴 W. B. Russel, D. A. Saville, and W. R. Schowalter, Colloidal
Dispersions 共Cambridge University Press, Cambridge, England, 1989兲.
关2兴 R. G. Larson, The Structure and Rheology of Complex Fluids
共Oxford University Press, Oxford, 1999兲.
关3兴 A. Einstein, Ann. Phys. 17, 459 共1905兲.
关4兴 I. M. Krieger and T. J. Doughty, Trans. Soc. Rheol. 3, 137
共1959兲.
关5兴 R. L. Hoffman, Trans. Soc. Rheol. 16, 155 共1972兲.
关6兴 B. J. Ackerson and P. N. Pusey, Phys. Rev. Lett. 61, 1033
共1988兲.
关7兴 M. D. Haw, W. C. K. Poon, and P. N. Pusey, Phys. Rev. E 57,
6859 共1998兲.
关8兴 H. Watanabe, M.-L. Yao, T. Shikata, H. Niwa, Y. Morishima,
N. P. Balsara, and H. Wang, Rheol. Acta 37, 1 共1998兲.
关9兴 B. J. Maranzano and N. J. Wagner, J. Chem. Phys. 117, 10291
共2002兲.
关10兴 J. F. Brady and G. Bossis, Annu. Rev. Fluid Mech. 20, 111
共1988兲.
关11兴 B. Dünweg and K. Kremer, J. Chem. Phys. 99, 6983 共1993兲.
关12兴 I.-C. Yeh and G. Hummer, J. Phys. Chem. B 108, 15873
共2004兲.
关13兴 D. M. Heyes, J. Phys.: Condens. Matter 19, 376106 共2007兲.
关14兴 D. R. Foss and J. F. Brady, J. Fluid Mech. 407, 167 共2000兲.
关15兴 A. Sierou and J. F. Brady, J. Fluid Mech. 506, 285 共2004兲.
关16兴 D. I. Dratler and W. R. Schowalter, J. Fluid Mech. 325, 53
共1996兲.
关17兴 D. I. Dratler, W. R. Schowalter, and R. L. Hoffman, J. Fluid
Mech. 353, 1 共1997兲.
关18兴 A. J. C. Ladd, Phys. Rev. Lett. 70, 1339 共1993兲.
关19兴 A. Malevanets and R. Kapral, J. Chem. Phys. 112, 7260

共2000兲.
关20兴 P. Ahlrichs and B. Dünweg, J. Chem. Phys. 111, 8225 共1999兲.
关21兴 H. Tanaka and T. Araki, Phys. Rev. Lett. 85, 1338 共2000兲.
关22兴 J. T. Padding and A. A. Louis, Phys. Rev. E 74, 031402
共2006兲.
关23兴 P. J. Atzberger, P. R. Kramer, and C. S. Peskin, J. Comput.
Phys. 224, 1255 共2007兲.
关24兴 N. Sharma and N. A. Patankar, J. Comput. Phys. 201, 466
共2004兲.
关25兴 M. Fujita and Y. Yamaguchi, Phys. Rev. E 77, 026706 共2008兲.
关26兴 Y. Nakayama and R. Yamamoto, Phys. Rev. E 71, 036707
共2005兲.
关27兴 Y. Nakayama, K. Kim, and R. Yamamoto, Eur. Phys. J. E 26,
361 共2008兲.
关28兴 T. Iwashita, Y. Nakayama, and R. Yamamoto, J. Phys. Soc.
Jpn. 77, 074007 共2008兲.
关29兴 T. Iwashita, Y. Nakayama, and R. Yamamoto, Prog. Theor.
Phys. 178, 共Suppl.兲, 86 共2009兲.
关30兴 T. Iwashita and R. Yamamoto, Phys. Rev. E 79, 031401
共2009兲.
关31兴 A. Shakib-Manesh, P. Raiskinmaki, A. Koponen, M. Kataja,
and J. Timonen, J. Stat. Phys. 107, 67 共2002兲.
关32兴 J. Kromkamp, D. T. M. van den Ende, D. Kandhai, R. G. M.
van der Sman, and R. M. Boom, J. Fluid Mech. 529, 253
共2005兲.
关33兴 P. M. Kulkarni and J. F. Morris, Phys. Fluids 20, 040602
共2008兲.
关34兴 H. J. H. Clercx and P. P. J. M. Schram, Phys. Rev. A 46, 1942
共1992兲.
关35兴 R. Yamamoto and A. Onuki, Phys. Rev. E 58, 3515 共1998兲.

061402-7