Reprinted from

FULL PAPER
Anisotropic Velocity Fluctuations and Particle Diﬀusion in Sedimentation
Adnan HAMID and Ryoichi YAMAMOTO
J. Phys. Soc. Jpn. 82 (2013) 024004# 2013 The Physical Society of Japan

Person-to-person distribution by the author only. Not permitted for publication for institutional repositories or on personal Web sites.

FULL P APERS

Journal of the Physical Society of Japan 82 (2013) 024004
http://dx.doi.org/10.7566/JPSJ.82.024004

Anisotropic Velocity Fluctuations and Particle Diﬀusion in Sedimentation
Adnan HAMID and Ryoichi YAMAMOTOy
Department of Chemical Engineering, Kyoto University, Kyoto 615-8510, Japan(Received October 6, 2012; accepted December 3, 2012; published online January 17, 2013)

Direct numerical simulations were performed using a smooth proﬁle method to investigate the steady-state
sedimentation of monodisperse spherical particles in an incompressible ﬂuid at ﬁnite Peclet numbers 0  Pe  115.
Hydrodynamic interactions caused strong ﬂuctuations in the instantaneous velocity of the particles around the mean
settling velocity. We found that the amplitude of these velocity ﬂuctuations increases in direct proportion with thesquare of the Stokes velocity at higher Peclet numbers, where sedimentation is dominated by non-equilibriumhydrodynamic ﬂuctuations. The diﬀusive behaviour of the particles was observed to be in a steady state over long timescales, and the steady-state self-diﬀusion coeﬃcient was found to increase linearly with the Peclet number. Our results
provide new insights into the anisotropy of vertical and horizontal diﬀusion. This anisotropy increases with anincreasing Peclet number, plateauing at a high Peclet number.
KEYWORDS: direct numerical simulation, smooth profile method, sedimentation, diffusion anisotropy, Peclet number

1.

Introduction

The sedimentation of monodisperse colloidal particles
is relevant to modern chemical and petroleum industries.
Although many experimental,1–6) theoretical,7–9) and numerical10–24) investigations have been performed to provide a
basic understanding of sedimentation, knowledge of thelong-range interactions between particles remains incomplete. These inter-particle interactions, caused by ﬂuid ﬂow,
are referred to as hydrodynamic interactions (HIs). The ﬁrst
signiﬁcant progress in sedimentation research was made by
Stokes,25) who calculated the sedimentation velocity of asingle sphere as, Vs¼0 ¼ ð2=9Þga2 ðp  f Þ=, where a isthe particle radius, p is the particle density, f is the ﬂuid
density, g is the gravitational acceleration, and  is the ﬂuid
viscosity. The eﬀect of the volume fraction on the mean
settling velocity, in the form Vs ¼ Vs¼0 ð1  kÞ, was later
established by Batchelor26) for the dilute limit   1, where is the particle volume fraction. This decrease in velocity is
attributed to the ﬂuid back ﬂow. The value of k was found to
be 6.55 for a uniform particle distribution, which is suitable
for most of the experimental conﬁgurations.27) Each particle
produces a velocity ﬁeld during sedimentation under acreeping ﬂow condition that aﬀects the motion of theother particles. This velocity ﬁeld causes a change in the
microstructure or position of the particles during sedimentation, thereby resulting in long-range HIs between theparticles. These long-range HIs cause ﬂuctuations in the
instantaneous particle velocity around the mean settling
velocity and largely depend upon the system size, volume
fraction and polydispersity.
Theoretical results7–9) and early simulations13) suggest thatvelocity ﬂuctuations converge to a ﬁnite value depending on
the size of the container, whereas experiments1,2) show nosuch dependency. Furthermore, ﬂuctuation magnitudes arestrongly anisotropic with vertical velocities being greaterthan horizontal velocities. These velocity ﬂuctuations arefound to be strongly dependent on volume fraction and aspectratio (i.e., the ratio of the vertical to horizontal lengths) in
simulations,4,11,12) whereas experiments1) didn’t observe suchdependency. This contradiction between experiment andsimulation results was solved by Segre et al.4) by introducing

the concept of a characteristic swirl size or correlation length;
i.e., velocity ﬂuctuations are explained as a function ofcell size, only when the cell size is less than the swirl size.
These results were later conﬁrmed by Ladd14) using lattice
Boltzmann simulations. Although the quantitative andqualitative behaviour of these velocity ﬂuctuations is wellunderstood, the screening mechanism for the velocity
ﬂuctuations remains unclear (i.e., the manner in whichcorrelations in the ﬂuctuations decay in time and space).
Shaqfeh and Koch8) proposed a Debye-like screening,whereas Ramaswamy28) proposed a stochastic convectiondiﬀusion model; however, these ideas were not conﬁrmedby the computer simulations of Ladd et al.10–14) Hinch9)
suggested another mechanism, in which the bottom wallacts as a sink for ﬂuctuations, that was conﬁrmed using
simulations.13,14) Some authors6,14) posit that stratiﬁcationand polydispersity also play key roles in screening. Additional details can be found in a review paper by Guazzelli.29)
Self-diﬀusion is one of the most important parameters of
sedimentation due to its role in mixing and other chemical
processes. Although the discovery and recognition of theanomalous behaviour of particle diﬀusion in sedimentationis relatively new, the subject has already captured theattention of a growing number of researchers in diverse
ﬁelds, including the study of suspension mechanics,
ﬂuidization, materials processing and granular ﬂows. This
diﬀusion is referred to as hydrodynamic diﬀusion because
its origin lies in the HIs between the particles.
Similar to velocity ﬂuctuations, diﬀusion is also strongly
anisotropic, and the vertical diﬀusion is greater than thehorizontal diﬀusion. This anisotropic behaviour arises fromthe anisotropic nature of the velocity ﬂuctuations and theirrelaxation times.30) Early simulations10) on self-diﬀusion
observed a large value for the anisotropy of vertical andhorizontal diﬀusion. This large anisotropy was later
corrected by increasing the volume fraction and aspectratio of the periodic box.11,20) Despite the importance ofhydrodynamic diﬀusion in chemical industrial research,
only a limited number of simulations10,11,19,20) have been
performed to address its behaviour.
The non-equilibrium properties of colloidal particles are
diﬃcult to simulate because they involve long-range, many-

024004-1#2013 The Physical Society of Japan

Person-to-person distribution by the author only. Not permitted for publication for institutional repositories or on personal Web sites.

FULL PAPERS

J. Phys. Soc. Jpn. 82 (2013) 024004

body HIs that require enormous temporal and spatial scales. A
number of diﬀerent methods10,19,20) are available to simulate
non-equilibrium colloidal phenomena; excellent reviews ofthese methods are presented by Padding et al.18) and Ladd
et al.15) We performed a direct numerical simulation (DNS)
using the smooth proﬁle (SP) method, which replaces theoriginal sharp boundaries between the particle and the host
ﬂuid with a diﬀuse interface of ﬁnite thickness. This approachenables us to use a ﬁxed Cartesian grid, thereby signiﬁcantly
improving the speed of numerical computations. The mainfeatures of our method are explained in Sect. 2.
In contrast to the majority of previous studies,3,5,13,20)
which were performed in a non-Brownian regime(Pe ! 1), we focused our work in the ﬁnite Peclet number(Pe) regime ranging from 0 to 115, which is a signiﬁcantly
higher range than that studied by Padding et al.18,19)(Pe  12). Sedimentation at ﬁnite Pe regime describes therelative eﬀects of thermal ﬂuctuations (caused by collisions
between colloids and solvent particles) and hydrodynamic
ﬂuctuations (caused by hydrodynamic ﬂows of solvent). In
the present study, we investigated the imprints of these
relative eﬀects on the anisotropic behaviour of velocity
ﬂuctuations and diﬀusion. Padding et al.17–19) were the ﬁrst
to study sedimentation in the ﬁnite Pe using coarse-grained
stochastic rotational dynamics (SRD) simulations. In SRD,
the highest achievable Pe is limited by the constraints on the
Mach and Reynolds numbers (Re), whereas in the SP
method, the highest Pe is limited only by Re. Consequently,
a wide range of Pe is achievable, using our method. The
present work validates the work of Padding17–19) and extends
it to Pe  115. This paper also provides new insights into
the anisotropic behaviour of diﬀusion at ﬁnite Pe. Section 2
explains the simulation method, and Sect. 3 describes
selection of the working parameters. The results areprovided in Sect. 4, and we present conclusions based on
our work in Sect. 5.
2.

Simulation Method

In this section, we brieﬂy explain the salient features ofour method. The detailed formulas, algorithm and applicability of the SP method can be found in previous
publications.31–38) In the SP method, the colloid surface is
not treated as a sharp interface lacking thickness; rather,
an intermediate region is introduced at the surface. This
intermediate region or interface has a width comparable to
the grid spacing, and the colloid’s density proﬁle is deﬁned
such that it changes smoothly within that interface.
Quantities such as the velocity and pressure are deﬁned
over the entire computational domain, which includes thecolloid as well as the solvent. A smooth proﬁle function
0  ðx; tÞ  1 is used to distinguish between the ﬂuid andparticle domains. Here,  ¼ 0 stands for ﬂuid domain, and ¼ 1 indicates the particle domain.
For an arbitrary particle i, the motion of the colloidal
particles is obtained by solving Newton’s equations of
motion:_ i ¼ Vi;Mi V_ i ¼ FH þ Fc þ Fext þ GV ; R
ð1Þ
i

i

i

i

I i  _ i ¼ N H þ N ext þ G ;
i

i

i

ð2Þ

where Ri is the position of the particle, V i and i are thetranslational and rotational velocity of the particle, respec-

A. HAMID and R. YAMAMOTO

tively. Mi and I i are the mass and moment of inertia,H
respectively. NH
i and Fi are the hydrodynamic torque andforce, respectively, exerted by the solvent on the particle.
ext
Fext
i and N i are the external force and torque, respectively.
V
Gi and Gi are the random force and torque due to thermal
ﬂuctuations, which can be described as, hGni ðtÞGnj ð0Þi ¼
n IðtÞij , where hGni ðtÞi ¼ 0 and n (n 2 V; ) is theparameter that controls the temperature of the system T .
The actual value of the particle temperature is determined byusing the long-time diﬀusion coeﬃcient of the equilibrium
system. When simulating a Brownian particle with HIs, thediﬀusion coeﬃcient is aﬀected by ﬁnite size eﬀects (i.e.,
the artefacts arises when system size is smaller than thecorrelation length). These eﬀects can be accounted for via
D¼0¼ D0 KðÞ, where D¼0
is the thermal diﬀusion00
coeﬃcient of a Brownian particle at inﬁnite dilution, D0
is the thermal diﬀusion coeﬃcient of Brownian particles
obtained for  6¼ 0 and KðÞ is the coeﬃcient that represents
the eﬀects of the ﬁnite volume fraction of dispersed particles
under periodic boundary conditions.39) Finally, the temperature of the system can be determined using the Stokes–
Einstein equation, kB T ¼ 6aD¼0
0 , where kB denotes
the Boltzmann constant. The detailed implementation ofBrownian motion with some test cases is presented by
Iwashita et al.33–35)
The direct inter-particle interaction Fci is presented as atruncated Lennard-Jones (LJ) potential with the large powersof 24 : 12, which
can be deﬁned as:
8 " 24  12 #><4þ  ðrij  21=12 Þ;
rij
rij
ULJ ðrij Þ ¼
ð3Þ
>
:
1=12
0
ðrij > 2Þ;where rij ¼ jRi  R j j. The parameters and  denote thelength and energy units of the LJ potential, respectively,
where ¼ 2a represents the particle diameter.
In the SP method, the solvent motion is obtained bysolving the modiﬁed Navier–Stokes equation for the total
velocity ﬁeld v and pressure ﬁeld p under the condition ofincompressibility (r  v ¼ 0) as follows:f ð@t v þ v  rvÞ ¼ rp þ r2 v þ f f p :ð4ÞThe rigidity of the particle and ﬂuid/particle non-slip
boundary conditions are incorporated via the body force
f p .31,32)
A similar smooth proﬁle40) was adopted in the previously
proposed ﬂuid particle dynamics (FPD) method, whereparticles are modelled as a highly viscous ﬂuid. In contrast,
we treated particles as non-deformable solids such thatno additional constraint arises. Moreover, Luo et al.41) has
performed error estimation and comparison of SP methodwith force coupling method (FCM)42) and a well-resolved
DNS43) for several problems. They found that the SP method
is resolving the near-ﬁeld ﬂows as well as the far-ﬁeld ﬂows
accurately. Furthermore, they established after the rigorous
analysis that compared to FCM and DNS, SP method iscomputationally faster and advantageous for simulating
moving particles because it avoids the complex discretizations around the particles.
The Peclet number, which is the ratio of convection to
diﬀusion, can be deﬁned as

024004-2#2013 The Physical Society of Japan

Person-to-person distribution by the author only. Not permitted for publication for institutional repositories or on personal Web sites.

FULL PAPERS

J. Phys. Soc. Jpn. 82 (2013) 024004

Pe ¼

aVs¼0
:
D0

2.5

ð5Þ
2

g(r)

One can expect HIs to play a dominant role at higher Pe,
thereby leading to changes in the microstructure.
3.

A. HAMID and R. YAMAMOTO

Simulation Parameters

Pe=0
Pe=6
Pe=29
Pe=57
Pe=115

1.5
1

A number of simulations were performed under periodic
boundary conditions with Pe ranging from 0 to 115 atvolume fraction  ¼ 0:02. A cubic periodic box ofdimension L ¼ 128 and spherical particles of size 4
were used in all of the simulations. This particular particle
size has already been used by Nakayama et al.38) for
diﬀerent test cases and it accurately produced the near andfar-ﬁeld ﬂows, pressure distribution, drag force and lubrication interactions. We deﬁned  as the lattice spacing, whichis also taken as the unit of length. For given values of  andf , the remaining units of mass, time, pressure and energy
are respectively deﬁned as f 3 , f 2 =, 2 =f 2 , and2 =f . The thermal ﬂuctuations kB T were kept constant at
0:3, and the  of the LJ potential was set at 2.5. Gravity
was introduced in the z-direction, and Pe was increased byincreasing the gravity. We constrained the center of mass of
ﬂuid to avoid indeﬁnite accelerating of the system. The
parameters in the Navier–Stokes equation were set such thatparticle Reynolds number, Re ¼ f aVs¼0 =, remained lessthan 0.2 (i.e., creeping ﬂow, Stokes regime). Whenincreasing gravity, it is important that the Re remain within
this creeping ﬂow condition. The majority of the experimental studies1,2,4,29) have been performed at very low Re,
usually of the order of 105 or less. Because a relative
deviation from the Stokes regime will scale with the squareof Re, we can probe relatively high Re values without
straying from the relevant experimental conditions. We kept
the particle to ﬂuid density ratio 5 to have a large range ofPe, keeping Re  0:2. Particles require a certain period oftime after the simulation begins to acquire their steady-state
velocities for a given initial conﬁguration. We monitored
the data, and only that corresponding to the steady-state
velocities was used. Simulations were performed for asuﬃciently long period to yield statistically meaningful data
for this analysis.
The HIs are increased with the increase of Pe. The
additional hindrance caused by these HIs changes the
microstructure of the system. A quantitative measure of
the microstructure at the particle scale is provided by theradial distribution function (RDF)
*
+
2L3 X
gðrÞ ¼ 2
ðr  rij Þ ;
ð6Þ
N
i<j
where N is the total number of particles, r ¼ jrj, rij ¼
Ri  R j , h P
i denotes an ensemble average, and thesummationi<j is taken over all particle pairs. Thisdeﬁnition of gðrÞ implies that 4gðrÞr 2 r represents themean number of particles in a shell of radius r and thickness
r surrounding a particle, set at the origin.44) Figure 1
shows the gðrÞ for diﬀerent Pe. The peak of the function
increases with increasing Pe, demonstrating the formation of
aggregates or clusters. This clustering was expected as
random motion progressively disappears with increasing Pe
which in turn induces the cluster formation. The similar

0.5
0
1

2

3

4

5

r/a
Fig. 1. (Color online) Radial distribution function gðrÞ for diﬀerent Pe.
Simulations are performed in a cubic periodic box of length L ¼ 128
with a particle size of a ¼ 4, keeping  and kB T constant respectively at0:02 and 0.3.

phenomena has also been observed by Brady et al.,16) while
investigating the relative eﬀect of thermal and shear force.
They reduced the eﬀect of thermal ﬂuctuations by increasing
the shear rate and observed the similar cluster formation.
Figure 1 also shows that the results are not distinguishable
from the equilibrium results (Pe ¼ 0) for Pe ¼ 6, showing
the strong eﬀect of thermal ﬂuctuations, whereas for Pe > 6,
these results are diﬀerentiable, evincing a progressively
dominating eﬀect of HIs. A large change in the RDF
between Pe  29 and Pe  57 indicates that HIs overpower
the phenomena.
4.

Results and Discussion

4.1 Spatial correlations of velocity ﬂuctuations
To study the eﬀect of Pe on the ﬂow pattern, we deﬁne
the following spatial correlation functions for the velocity
ﬂuctuations:
*
+
2L3 X
Cz ðrÞ ¼ 2
Viz V jz ðr  rij Þ ;
ð7Þ
N
i<j
where Viz ¼ Viz  Vsed , Viz is the temporal velocity of theith particle in the z-direction, and Vsed ¼ hViz i is the mean
settling velocity of the particles. We deﬁne Cz ðzÞ and Cz ðxÞ,with respect to the distance vector r, in either the verticalr ¼ zz or horizontal r ¼ xx direction.
Figure 2 shows the decay of the spatial correlation
functions of the z-component of velocity as a function of
the vertical (z) and horizontal (x) distances, where thedistance is normalized by the particle radius a. In the mainplots of Figs. 2(a) and 2(b), the correlation functions arenormalized by the initial value of the ﬂuctuations, whereas
the correlation functions in the ﬁgure insets are normalized
by the square of the Stokes velocity. Figure 2(b) shows apersistent positive correlation in the direction of gravity,
whereas Fig. 2(a) shows a rapid decay of all correlations in
the direction perpendicular to gravity followed by a distinct
region of anti-correlation similar to that observed inexperiments4) and simulations.19) Notably, for a suﬃciently
large system, the horizontal correlations should decay
exponentially and become negative before settling to zero.
The minima of these anti-correlations is known as thecorrelation length4)  20a1=3 , i.e.,
is twenty-fold

024004-3#2013 The Physical Society of Japan

Person-to-person distribution by the author only. Not permitted for publication for institutional repositories or on personal Web sites.

FULL PAPERS

J. Phys. Soc. Jpn. 82 (2013) 024004

1

0

0.2
0
0

3

6

9

12

0.4

)

φ=0 2

Cz(z)/(Vs

(b)

0.6

0

10

0

20

1

5

10

15

30

40

50

1
0.5
0.8
0

Pe=6
Pe=11
0.2 Pe=29
Pe=57
Pe=115
0
3

10

t/ts

0.4

0

Eq. (13)-20

15

Cz(t)/Cz(0)Cz(z)/Cz(0)

0.8

-1

10

10-3

0.2

x/a
1

Cx(t)/(Vsφ=0)2

(a)

0.5

6

9

12

10

10

Eq. (14)
-1

-2

0.4

10

0

20

40

60

0.2
0

15

z/a

0

Fig. 2. Spatial correlation function of velocity ﬂuctuations of the
z-component of the velocity as a function of the distance: (a) shows thespatial correlation function perpendicular to gravity and (b) shows thespatial correlation function parallel to gravity. In the main plots, thecorrelation functions are normalized by the initial position ﬂuctuations,whereas the correlation functions in the insets are normalized by the squareof the Stokes velocity. The distance is normalized by the particle radius a.

greater than the inter-particle distance. Although anticorrelations are also clearly visible in our results, a larger
system size is needed to verify the correlation length
obtained experimentally. This determination is, however,beyond the scope of the present study.
The particle motion is aﬀected by the ﬁnite system size inour simulations. A large system size is needed to reproduce
the saturation of the velocity ﬂuctuations, but it requires
enormous computational time and resource. Therefore, mostof the computational studies10,11,14,19,20) are in the same
situation as ours. The systematic experimental2,4) andsimulation11,13) studies, however, addressed that the eﬀects
of ﬁnite system size on velocity ﬂuctuations can beeﬀectively explained using a concept of ﬁnite size scaling.19)
We believe that performing critical tests for the scaling
concept possesses an important value.
4.2 Temporal correlations of velocity ﬂuctuations
To diﬀerentiate between the eﬀects of the hydrodynamic
and thermal ﬂuctuations, we deﬁne the temporal autocorrelation functions of the velocity ﬂuctuations as follows:Cx ðtÞ ¼ hVix ðtÞVix ð0Þi;Cz ðtÞ ¼ hViz ðtÞViz ð0Þi:

0.6

(b)

0

)

0.4

(a)

φ=0 2

0.6

Pe=6
Pe=11
Pe=29
Pe=57
0.6 Pe=115
0.8

Cz(t)/(Vs0.8

1

Cx(t)/Cx(0)Cz(x)/(Vsφ=0)2

1

Cz(x)/Cz(0)A. HAMID and R. YAMAMOTO

ð8Þ
ð9Þ

Figure 3 shows the time decay of the correlation functions
for the z and x components of the velocity ﬂuctuations. Two
diﬀerent normalizations were used: the main plots of Fig. 3

40

80

120

160

200

t/ts
Fig. 3. Temporal autocorrelation functions of the velocity ﬂuctuationsof the vertical ½Cz ðtÞ and horizontal ½Cx ðtÞ components of the velocity.
The main plots show the correlation functions normalized by the initial
ﬂuctuations, whereas the temporal correlation functions shown in the insets
are normalized by the square of the Stokes velocity on a semi-log scale.
Insets also show that at large Pe, these correlation functions relax
exponentially of the form C ðtÞ ¼ ðVH Þ2 expðt= H Þ, where VH( 2 x; z) and H denote the amplitude and the relaxation time of thehydrodynamic velocity ﬂuctuations, respectively. Time is normalized by the
Stokes time ts , where ts ¼ a=Vs¼0 .

show the temporal correlation functions normalized by the
initial values Cz ð0Þ and Cx ð0Þ, whereas the ﬁgure insets
show the temporal correlation functions normalized by thesquare of the Stokes velocity in a semi-log plot. Time wasnormalized by the Stokes time ts , where ts ¼ a=Vs¼0 .
The main plots of Figs. 3(a) and 3(b) show that therelaxation time is diﬀerent for the vertical and horizontal
correlations, whereas the experiments of Nicolai1) did notshow a large diﬀerence at these time scales. The theory
proposed by Koch30) suggests that this diﬀerence in timescales come from the vertical periodic boundary conditions.
Therefore, the time required by a particle to sample all thevertical position is much smaller than that required to sample
horizontal positions. This diﬀerence in time scales can be
reduced by increasing the system size. The computer
simulations of Ladd10) and padding19) also observed thesame phenomena while using vertical periodic boundary
conditions. In contrast to simulations, experiments1,4) use
much larger system which reduce the diﬀerence in verticaland horizontal relaxation time. At low Pe, the strong eﬀects
of thermal ﬂuctuations cause a rapid decay of correlation
functions, whereas at higher Pe, exponential decay isevident, as shown in Fig. 3.

024004-4#2013 The Physical Society of Japan

Person-to-person distribution by the author only. Not permitted for publication for institutional repositories or on personal Web sites.

FULL PAPERS

J. Phys. Soc. Jpn. 82 (2013) 024004

In an experimental study, Nicolai et al.1) observed an
exponential relaxation of the form!
t
2
C ðtÞ ¼ ðVH Þ exp 
;
ð10Þ
H

where VH ( 2 x; z) and H denote the amplitude and therelaxation time of the hydrodynamic velocity ﬂuctuations,
respectively. Hinch9) predicted the scaling relations
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
ð11ÞVH =Vs¼0 ¼ A1 L=a;pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
ð12Þ
H =ts ¼ A2 L=ðaÞ;where A1 and A2 are constants that depend on the speciﬁc
conﬁguration of the system (e.g., the container shape, swirl
size, particle shape, and polydispersity). Cunha et al.20)
explained the scaling relations predicted by Hinch9) by
hypothetically dividing the simulation box into two equalparts. The imbalance in the particle weight caused byvariations in the particle number in both parts is balancedby the Stokes drag. Velocity ﬂuctuations consist of thermaland hydrodynamic component. At high Pe, hydrodynamic
ﬂuctuations dominate the phenomena, whereas thermal
ﬂuctuations play a key role at low Pe. This dominant roleof hydrodynamic ﬂuctuations at higher Pe causes thecorrelation functions to fall on the same curve when scaledwith the square of the Stokes velocity, as shown in the insets
of Fig. 3, which is in good agreement with the simulations.19,20) This normalization highlights the eﬀect of the
non-equilibrium hydrodynamic ﬂuctuations.
The present DNS results support the following forms forthe velocity autocorrelations:tCx ðtÞ ¼ C0 ðtÞ þ ðVHx Þ2 exp  x ;
ð13Þ
H


t
Cz ðtÞ ¼ C0 ðtÞ þ ðVHz Þ2 exp  z ;ð14Þ
H

where C0 ðtÞ represents the velocity autocorrelation functionin the presence of thermal ﬂuctuations but without gravity,
which quantitatively becomes negligible at higher Pe. We
continue our analysis by assuming the scaling arguments20)
as explained in Eqs. (11) and (12). The perfectors A1 andA2 are determined based on exponential ﬁts both on verticaland horizontal correlation functions as Az1  0:156, Az2 
1:08, Ax1  0:028, and Ax2  0:22. The diﬀerence in verticaland horizontal prefactors indicates the anisotropic behaviourof velocity ﬂuctuation with the vertical ﬂuctuations larger
than the horizontal. This anisotropy of velocity ﬂuctuations
varies from 2.5 to 4 in simulations,10,12,20,22) whereas
experiments1,4) have found relatively low value of 2:5.
Notably, the ratio of the vertical to horizontal hydrodynamic velocity ﬂuctuations in this study, VHz =VHx  2:36,
is in good agreement with the experimental1,4) andsimulations22) studies. From the asymmetry of the system
induced by the gravity, velocity ﬂuctuations are expected to
be anisotropic.
4.3 Anisotropic behavior of diﬀusion
Diﬀusion refers to the ﬂuctuating motion of particles.
Individual particles lose the memory of their velocity after
prolonged HIs, and thereafter follow a random-walk
diﬀusion process. This eﬀect is evident in the behaviourof the temporal autocorrelation functions of the velocity

A. HAMID and R. YAMAMOTO

ﬂuctuations and the linear growth of the mean-squaredisplacement (MSD).
The equilibrium self-diﬀusion coeﬃcient can be obtainedusing the Green–Kubo equation
Z1
C0 ð Þ d
ð15Þ
D0 ¼
0

or the Einstein relationshipD0 ¼ lim

1

t!1 6thðRi ðtÞ  Ri ð0ÞÞ2 i:

ð16Þ

To examine the anisotropic particle diﬀusion in sedimentation, we deﬁne the temporal diﬀusion constant both in thehorizontal
1hðRix ðtÞ  Rix ð0ÞÞ2 i
2t

ð17Þ

1
hðRiz ðtÞ  Riz ð0Þ  Vsed tÞ2 i:
2t

ð18Þ

Dx ðtÞ ¼
and vertical directions
Dz ðtÞ ¼

Figure 4 shows the time-dependent diﬀusion coeﬃcientdeﬁned by Eqs. (17) and (18) for various Pe. The diﬀusion
coeﬃcients are normalized by the equilibrium self-diﬀusion
coeﬃcient D0 , which can be found via Eq. (15) or (16) at
Pe ¼ 0, whereas time is normalized by the Stokes time.
One can see in Fig. 4 that self-diﬀusion becomes highly
anisotropic at higher Pe. We determined the long-time selfdiﬀusion coeﬃcients D ¼ limt!1 D ðtÞ by ﬁtting thesimulation data with   t
D ðtÞ ¼ D 1  exp 
;
ð19Þ
r

where D , r , and are the ﬁtting parameters. In Figs. 4 and5, values of Dz =D0 much larger than unity suggest that therole of HIs cannot be overlooked, even at low Pe. A similar
phenomena was observed by Padding et al.19)
To understand the diﬀusivity, we consider the total
diﬀusion coeﬃcient D as the sum of thermal D0 andhydrodynamic diﬀusivity DH , where latter can be estimated
via the aforementioned scaling arguments,
DH  ðVH Þ2 H ¼ A1 A2 Vs¼0 a1=2 ðL=aÞ3=2¼ A1 A2 D0 Pe

1=2

ðL=aÞ

3=2

:

ð20Þ
ð21Þ

The following simple scaling relationships are predicted forthe present simulation results:
Dx =D0 ¼ 1 þ Ax1 Ax2 Pe 1=2 ðL=aÞ3=2 ¼ 1 þ 0:158Pe; ð22Þ
Dz =D0 ¼ 1 þ Az1 Az2 Pe 1=2 ðL=aÞ3=2 ¼ 1 þ 4:313Pe: ð23ÞThese scaling arguments indicate that the vertical andhorizontal diﬀusion coeﬃcients increase linearly with Pe,
but with a smaller prefactor in the horizontal direction.
Figure 5 shows a comparison of the steady-state selfdiﬀusion coeﬃcients with those predicted by the scalingargument. The simulation results show that diﬀusion
coeﬃcients increase linearly with increase of hydrodynamic
ﬂuctuations both parallel and perpendicular to gravity.
A deviation from the predicted scaling equations isevident at Pe < 6, which shows the transition from athermal-ﬂuctuations-dominated regime to a hydrodynamicﬂuctuations-dominated regime. Since, above scaling rela-

024004-5#2013 The Physical Society of Japan

Person-to-person distribution by the author only. Not permitted for publication for institutional repositories or on personal Web sites.

FULL PAPERS

J. Phys. Soc. Jpn. 82 (2013) 024004

A. HAMID and R. YAMAMOTO

3

20

10

16

Eq. (22)
Simulation of Dx /Do
Eq. (23)
Simulation of Dz /Do(a)Dβ /D0Dx /D0

2

10

12
8

1

10
4

0

0
0

60

120

180

240

10

300

10-1

100

101

t/ts

Pe

Dz /D0

500
Pe=6
Pe=11
Pe=29
400
Pe=57
Pe=115300 Eq. (19)102

Fig. 5. A comparison of steady state self-diﬀusion coeﬃcientsD ð 2 x; zÞ both in vertical and horizontal direction with those predictedby scaling Eqs. (22) and (23) as a function of Pe. The self-diﬀusion
coeﬃcients are normalized by the equilibrium diﬀusion coeﬃcient D0 .

(b)

30
200
24

Dz /Dx

100
0
0

60

120

180

240

300

t/ts

18
12
Simulation
Scaling
6

Fig. 4. The horizontal and vertical direction self-diﬀusion coeﬃcientsnormalized by the equilibrium diﬀusion coeﬃcient D0 as a function of time
at various Pe. Both (a) and (b) show the horizontal and vertical timedependent self-diﬀusion coeﬃcients, respectively. Time is normalized by
the Stokes time. The volume fraction is set at 0.02. Dashed lines show the
exponential ﬁt on the simulation data.

tions are derived by adding pure thermal and hydrodynamic
forces. Therefore, scaling relations are expected to work
well at either very low Pe (thermal force dominates) or highPe (hydrodynamic force dominates). Hence, deviation fromthe scaling relations indicates the interplay between these
two forces. Our simulations show good agreement with thepredictive Eqs. (22) and (23), determined from scalingarguments. These scaling arguments also indicate that
anisotropy in diﬀusion increases with the increase of Pe,
reaching a steady value at high Pe. From a simple
dimensional analysis of Eq. (23), one can also predict thatDz =aVs¼0 should collapse to a single value at higher Pe. In
our simulation, this value 4 and is in good agreement withexperimental results.1,2)
The majority of diﬀusion research is focused in the nonBrownian regime, and anisotropic behaviour of diﬀusion at
ﬁnite Pe has not yet been explored. We have attempted toinvestigate this anisotropic behaviour at ﬁnite Pe and foundthat it increases with increasing Pe and settles to a steadyvalue at higher Pe, as shown in Fig. 6. This observationsuggests that the eﬀect of thermal ﬂuctuations is signiﬁcantat Pe  29 before HIs begin to dominate the phenomena.
We can predict the anisotropic behaviour of the diﬀusion at
ﬁnite Pe with scaling arguments found in Eqs. (22) and (23).
Our data shows good agreement with the predicted diﬀusion
anisotropy, as shown in Fig. 6.

0
0

25

50

75

100

125

Pe
Fig. 6. The anisotropy of the vertical (Dz ) and horizontal (Dx ) selfdiﬀusion coeﬃcient as a function of Pe at a volume fraction of 0.02. The
scaling line is obtained by dividing the Eq. (23) with Eq. (22).

Because the diﬀusivity is the product of VH2 and H ,a higher diﬀusion anisotropy value was expected. Thishigher value originates from the diﬀerence in vertical andhorizontal relaxation times. The diﬀusion anisotropy can be
reduced by increasing the aspect ratio of the simulation
box.30) An increase in the aspect ratio reduces the diﬀerencein vertical and horizontal relaxation times and hence thediﬀusion anisotropy. Nicolai1) achieved diﬀusion anisotropy
in a non-Brownian regime of Dz =Dx ’ 7 at  ¼ 0:05,
whereas Padding et al.19) reported the same value insimulations using a periodic box of aspect ratio of 3 at
 ¼ 0:04. Ladd12) observed a much higher value of 77 in a
cubic simulation box for  ¼ 0:05. This higher value wasattributed to the small system size and the full periodic
boundary conditions that were used. Cunha et al.20) observed
the anisotropy to be 10 using an aspect ratio of 3 at
 ¼ 0:03. To quantitatively compare the anisotropy values
using an experimental system, one should consider that oursimulations depend on the system size, whereas experiments
do not. Another possible explanation for the discrepancy
between simulation and experimental results is the presenceof a side wall, which creates microstructure inhomogeneity
over time. The presence of even a small extent of

024004-6#2013 The Physical Society of Japan

Person-to-person distribution by the author only. Not permitted for publication for institutional repositories or on personal Web sites.

FULL PAPERS

J. Phys. Soc. Jpn. 82 (2013) 024004

Table I. Comparison of the observed diﬀusion anisotropy with previously
published experimental and simulated data. Nicolai performed theseexperiments with a ratio of vessel width to particle radius ranging from51 to 212. Ladd performed these simulations in a cubic periodic box at
Re ¼ 0:17 with 32 spheres. Padding and Cunha performed simulations withan aspect ratio of 3.


Dz =DxPresent Study

0.02

25

Nicolai2)

0.05

5–7.4

Ladd12)0.05

77

Padding19)

0.04

7

Cunha20)0.03

10Murashima, Dr. Hideyuki Mizuno, and Dr. John Molina for
their useful discussions and suggestions.

hamid@cheme.kyoto-u.ac.jp
ryoichi@chem.kyoto-u.ac.jp
1) H. Nicolai, B. Herzhaft, E. J. Hinch, L. Oger, and E. Guazzelli: Phys.
Fluids 7 (1995) 12.
2) H. Nicolai and E. Guazzelli: Phys. Fluids 7 (1995) 3.
3) H. Nicolai, Y. Peysson, and E. Guazzelli: Phys. Fluids 8 (1996) 855.
4) P. N. Segrè, E. Helbolzheimer, and P. M. Chaikin: Phys. Rev. Lett. 79(1997) 2574.
5) P. N. Segrè, F. Liu, P. Umbanhowar, and D. A. Weitz: Nature 409(2001) 594.
6) P. N. Segrè: Phys. Rev. Lett. 89 (2002) 254503.
7) R. E. Caﬂisch and J. H. C. Luke: Phys. Fluids 28 (1985) 759.
8) D. L. Koch and E. S. G. Shaqfeh: J. Fluid Mech. 224 (1991) 275.
9) E. J. Hinch: in Disorder and Mixing, ed. E. Guyon, Y. Pomeau, andJ. P. Nadal (Kluwer Academic, Dordrecht, 1988).
10) A. J. C. Ladd: Phys. Fluids A 5 (1993) 299.
11) A. J. C. Ladd: Phys. Rev. Lett. 76 (1996) 1392.
12) A. J. C. Ladd: Phys. Fluids 9 (1997) 491.
13) A. J. C. Ladd: Phys. Rev. Lett. 88 (2002) 048301.
14) N. Q. Nguyen and A. J. C. Ladd: J. Fluid Mech. 525 (2005) 73.
15) A. J. C. Ladd and R. Verberg: J. Stat. Phys. 104 (2001) 1191.
16) G. Bossis and J. F. Brady: J. Chem. Phys. 87 (1987) 5437.
17) J. T. Padding and A. A. Louis: Phys. Rev. Lett. 93 (2004) 220601.
18) J. T. Padding and A. A. Louis: Phys. Rev. E 74 (2006) 031402.
19) J. T. Padding and A. A. Louis: Phys. Rev. E 77 (2008) 011402.
20) F. R. Cunha, G. C. Abade, A. J. Sousa, and E. J. Hinch: J. Fluids Eng.
124 (2002) 957.
21) S. Y. Tee, P. J. Mucha, M. P. Brenner, and D. A. Weitz: Phys. Fluids
19 (2007) 113304.
22) M. C. Miguel and R. Pastor-Satorras: Europhys. Lett. 54 (2001) 45.
23) X. Yin and D. L. Koch: Phys. Fluids 20 (2008) 043305.
24) E. Kuusela, J. M. Lahtinen, and T. Ala-Nissila: Phys. Rev. E 69 (2004)
066310.
25) G. G. Stokes: Mathematical and Physical Papers (Johnson Reprint
Corporation, New York, 1966).
26) G. K. Batchelor: J. Fluid Mech. 52 (1972) 245.
27) J. M. Ham and G. M. Homsy: Int. J. Multiphase Flow 14 (1988) 533.
28) A. Levine, S. Ramaswamy, E. Frey, and R. Bruinsma: Phys. Rev. Lett.
81 (1998) 5944.
29) E. Guazzelli and J. Hinch: Annu. Rev. Fluid Mech. 43 (2011) 97.
30) D. L. Koch: Phys. Fluids 6 (1994) 2894.
31) Y. Nakayama and R. Yamamoto: Phys. Rev. E 71 (2005) 036707.
32) Y. Nakayama, K. Kim, and R. Yamamoto: Eur. Phys. J. E 26 (2008)
361.
33) T. Iwashita, Y. Nakayama, and R. Yamamoto: J. Phys. Soc. Jpn. 77(2008) 074007.
34) T. Iwashita and R. Yamamoto: Phys. Rev. E 79 (2009) 031401.
35) T. Iwashita, Y. Nakayama, and R. Yamamoto: Prog. Theor. Phys.
Suppl. 178 (2009) 86.
36) Y. Nakayama, K. Kim, and R. Yamamoto: Adv. Powder Technol. 21(2010) 206.
37) KAPSEL website [http://www-tph.cheme.kyoto-u.ac.jp/kapsel/].
38) Y. Nakayama, K. Kim, and R. Yamamoto: Eur. Phys. J. E 26 (2008)
361.
39) A. A. Zick and G. M. Homsy: J. Fluid Mech. 115 (1982) 13.
40) H. Tanaka and T. Araki: Phys. Rev. Lett. 85 (2000) 1338.
41) X. Luo, M. R. Maxey, and G. E. Karniadakis: J. Comput. Phys. 228(2009) 1750.
42) M. R. Maxey and B. K. Patel: Int. J. Multiphase Flow 27 (2001) 1603.
43) R. M. Kirby, T. C. Warburton, S. J. Sherwin, A. Beskok, and G. E.
Karniadakis: Proc. 2nd Int. Conf. Computational Technologies for
Fluid/Thermal/Chemical Systems with Industrial Applications, 1999.
44) D. C. Rapaport: The Art of Molecular Dynamics Simulation(Cambridge University Press, Cambridge, U.K., 2004) 2nd ed., p. 90.

ypolydispersity can also temper the diﬀusion. Experimental
and simulated diﬀusion anisotropy results are summarized inTable I.
5.

A. HAMID and R. YAMAMOTO

Concluding Remarks

In conclusion, we studied the steady-state sedimentation
at ﬁnite Pe by DNS using the SP method in a spatiallyunscreened regime. The present study focused on the relativeeﬀects of thermal and hydrodynamic ﬂuctuations on anisotropic behaviour of velocity ﬂuctuation and self-diﬀusion.
We observed that the thermal ﬂuctuations are noticeable at
Pe  29, before HIs signiﬁcantly dominate the phenomena.
We found that the amplitude of the velocity correlation
scales with the square of the Stokes velocity at large Pe, withvertical hydrodynamic velocity ﬂuctuations 2.36 times larger
than the horizontal hydrodynamic velocity ﬂuctuations. In
addition, we successfully tested the scaling for hydrodynamic velocity ﬂuctuations, predicted by Hinch9) on oursystem and used it to evaluate the anisotropic behaviour ofdiﬀusion. We inferred that the long time steady-state selfdiﬀusion coeﬃcient increased with the increase of Pe bothin the direction of gravity and perpendicular to gravity and it
became highly anisotropic at high Pe, with a vertical selfdiﬀusion coeﬃcient higher than the horizontal, in good
agreement with the scaling. This anisotropy of the verticaland horizontal diﬀusion increased as a function of Pe andlevelled oﬀ at higher Pe. Our results endorsed the ﬁndingsof Koch,30) who suggested that the change in the velocity
ﬂuctuations with the box shape is relatively modest, however
diﬀusion coeﬃcients changes rather signiﬁcantly with the
aspect ratio. Similar to Koch ﬁndings, velocity ﬂuctuations
in our simulations are also in good agreement with the othersimulation19,20) and experimental results,1) whereas diﬀusion
anisotropy is higher due to large vertical relaxation time. So
far, researchers have paid little attention to the anisotropic
behaviour of diﬀusion at ﬁnite Pe. We will extend oursimulation studies to investigate the salient features ofthis diﬀusion anisotropy at large volume fractions, various
system sizes, and large Re. Extensive simulations are under
way to study such situations specially to verify the scalingarguments at diﬀerent volume fractions and system sizes.
Acknowledgements
This work was supported by JSPS KAKENHI 23244087.
The author would like to extend his gratitude to Dr. Takahiro

024004-7#2013 The Physical Society of Japan