PHYSICAL REVIEW E 87, 022310 (2013)

Direct numerical simulations of anisotropic diffusion of spherical particles in sedimentation
Adnan Hamid* and Ryoichi Yamamoto†Department of Chemical Engineering, Kyoto University, Kyoto 615-8510, Japan(Received 12 December 2012; published 27 February 2013)
We investigated the scaling of the hydrodynamic velocity fluctuations and self-diffusion in the sedimentationof monodispersed spherical particles via direct numerical simulations using the smoothed profile method over amoderate range of volume fractions (0.01  φ  0.12). Hydrodynamic velocity fluctuations are visible at large
Peclet numbers (Pe), and they scale as (φL/a)1/2 at low volume fractions (φ  0.04). Their characteristics become
independent of volume fraction at moderate volume fractions (0.06  φ  0.12). Both vertical and horizontal
self-diffusion coefficients scale as (L/a)3/2 φ 1/2 at low volume fractions. At moderate volume fractions, the vertical
diffusion scales as (L/a)3/2 φ −1/2 ; in contrast, the horizontal diffusion is saturated with respect to volume fraction.
The diffusion anisotropy increases with increasing Pe and saturates at high Pe values. The saturated value remains
unchanged at low volume fractions, whereas further increase in the volume fraction decreases this anisotropy.
The reduction of this anisotropy is attributed to the φ −1/2 scaling of the vertical relaxation time at moderate
volume fractions; however, the horizontal relaxation time is independent of the volume fraction at this regime.
DOI: 10.1103/PhysRevE.87.022310

PACS number(s): 82.70.Dd, 31.15.xv, 46.15.−x, 82.20.Wt

I. INTRODUCTION

The settling of particles in a viscous fluid is a fundamental
phenomenon of nonequilibrium suspension dynamics, but the
topic is extremely challenging to examine because accurately
characterizing the long-range many-body interactions between
the particles proves difficult. These interparticle interactions,caused by fluid flow, are referred to as hydrodynamic interactions (HIs). These HIs may decay as slowly as 1/r,
strongly affecting the dynamic behavior of the suspension.
Although many experimental [1–6], theoretical [7–9], andnumerical [10–25] investigations have been performed toprovide a basic understanding of sedimentation, knowledge ofthe HIs remains incomplete. Significant progress was madeby Batchelor [26], who explained the hindered settling of
the sedimenting particles and obtained an expression of theφ=0
φ=0form Vs = Vs (1 − kφ) in which Vs= 29 ga 2 (ρp − ρf )/η
represents the Stokes velocity [27] of a single particle, a is theparticle radius, ρp is the particle density, ρf is the fluid density,
g is the gravitational acceleration, η is the fluid viscosity, φ isthe particle volume fraction, and k is a constant that shows theeffect of fluid backflow [28].
Due to the constantly changing configuration of the suspension microstructure, the velocity of the individual particles
fluctuates around their mean settling velocity. Inconsistencies
exist between the experiments [1,2], and the simulations [13]and theory [7,8]; the velocity fluctuations predicted by thetheory and numerical simulations depend on the system size,
but the experiments exhibit no such dependency. Subsequentstudies have shown that if the simulations exhibit dependence
on container size, some screening or cutoff length scalemust exist in the suspension microstructure that depends on
the particle volume fraction [29]. This inconsistency wassolved by Segre et al. [4], who introduced the concept ofa characteristic swirl size or correlation length: Velocity
fluctuations are explained as a function of cell size only when*†hamid@cheme.kyoto-u.ac.jp
ryoichi@cheme.kyoto-u.ac.jp1539-3755/2013/87(2)/022310(11)the cell size is less than the swirl size. These results werelater confirmed by Nguyen and Ladd [14], using the lattice
Boltzmann method (LBM). These velocity fluctuations arestrongly anisotropic, with the vertical velocity fluctuations
larger than the horizontal.
The disagreement between the experiments and the early
simulations stimulated further work to address the natureof the screening mechanism for the velocity fluctuations(i.e., the manner in which correlations in the fluctuations
decay over time and space). Shaqfeh and Koch [8] proposeda Debye-like screening, and Ramaswamy [30] proposed astochastic convection-diffusion model; however, these ideas
were not confirmed by the simulations of Ladd et al. [10–14].
Hinch [31] suggested another mechanism in which the bottom
wall acts as a sink for fluctuations. This mechanism wasconfirmed using simulations [13,14]. Some authors [6,14]
posited that stratification and polydispersity also play key rolesin screening. More details can be found in Guazzelli’s review
paper [32].
Self-diffusion is one of the most important parameters of
sedimentation due to its role in mixing and other chemical processes. This diffusion is referred to as hydrodynamic diffusion
because it originates in the HIs between the particles. As withthe velocity fluctuations, hydrodynamic diffusion obtained
from the simulations is affected not only by the size of thesimulation box but also by its shape. Diffusion is also strongly
anisotropic, with vertical diffusion greater than horizontal.
This anisotropic behavior arises from the anisotropic nature of
the velocity fluctuations and their relaxation times [33]. Early
self-diffusion simulations [10] yielded a large value for thisanisotropy. A more realistic value of the anisotropy was later
obtained by Ladd [12] and Cunha [19]; the former increased
the aspect ratio of the periodic box, and the latter used dynamic
simulations with an impenetrable lower boundary. Despite theimportance of hydrodynamic diffusion in industrial chemical
research, few simulations [10,11,18,19] have been performed
to address its behavior.
Under sedimentation, the microstructure of the suspension
is determined by the long-range HIs, which are characterized
by large time and length scales, features that render the022310-1©2013 American Physical Society

ADNAN HAMID AND RYOICHI YAMAMOTO

PHYSICAL REVIEW E 87, 022310 (2013)

nonequilibrium properties of the colloidal particles challenging to simulate. Different methods [10,18,19,34–36] areavailable to simulate the nonequilibrium colloidal phenomena,e.g., Stokesian dynamics [34], LBM [10–15], dissipative
particle dynamics [35], stochastic rotation dynamics (SRD)[16], and direct numerical simulations (DNS) [19]. Moredetails on these methods are presented by Padding et al. [17]and Ladd et al. [15]. In this study, we used DNS with thesmoothed profile (SP) method, which replaces the original
sharp boundaries between the particle and the host fluidwith a diffuse interface of finite thickness. This approachenabled us to use a fixed Cartesian grid, thereby significantly
improving the speed of the numerical computations. A similarsmoothed profile [37] was adopted in the previously proposed
fluid-particle dynamics (FPD) method in which particles aremodelled by a highly viscous fluid. We treated the particles asnondeformable solids such that no additional constraints arose.
Unlike most previous studies [3,5,13,19], which were
performed in a non-Brownian regime (Pe → ∞), we focused
our work in the finite Peclet number (Pe) regime with different
volume fractions and system sizes. Sedimentation at finite Pe
values describes the relative effects of thermal fluctuations(caused by collisions between colloids and solvent particles)
and hydrodynamic fluctuations (caused by hydrodynamic
flows of the solvent). We investigated the imprints of these
relative effects on the anisotropic behavior of the velocity
fluctuations and diffusion. Padding et al. [16–18] were the firstto simulate sedimentation at a finite Pe using coarse-grained
SRD simulations. The simplified dynamics allowed accurate
calculation of the transport properties; however, due to its
stochastic nature, SRD is less effective than pre-averaged LBM
for large Pe. In SRD, the highest achievable Pe is limitedby 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. Padding [16–18] focused his work in the low
volume fraction regime. The present work not only validates
the work of Padding at low volume fraction but also explores
the effects of moderate volume fraction on velocity fluctuations
and anisotropic diffusion.
The main objective of the present study is to investigate
the relative effects of thermal and hydrodynamic forceson velocity fluctuations and particle diffusion for different
volume fractions. In the previous studies [38,39], we haveexamined the finite Pe effects at a single volume fraction(φ = 0.02), whereas in the present study, we investigated thescaling relations of hydrodynamic velocity fluctuations andself-diffusion at different system sizes and volume fractions.
Working along the lines of Caflish and Luke [7], Hinch [31]previously suggested the scaling of velocity fluctuations andself-diffusion, which is valid for a low volume fraction.
However, the high volume fraction regime requires additional
attention. This paper provides insights into the anisotropic
behavior of diffusion at finite Pe. We also suggest valid scaling
relations of the velocity fluctuations, their relaxation times, andself-diffusion for a moderate volume fraction regime. In thispaper, Sec. II explains the simulation method, and Sec. III
describes the selection of the working parameters. The results
are provided in Sec. IV, and we present conclusions based on
our work in Sec. V.

II. SIMULATION METHOD

In the SP method, the colloid surface is not treated as a
sharp interface lacking thickness, but an intermediate region is
introduced at the surface. This intermediate region or interface
has a width comparable to or larger than the grid spacing,
and the colloid’s density profile is defined such that it changes
smoothly within that interface. Quantities such as the velocity
and pressure are defined over the entire computational domain,which includes the colloid and the solvent. A smoothed
profile function [36] 0  φ(x,t)  1 is used to distinguishbetween the fluid and particle domains in which x denotes
the particle spacial position. Here, φ = 0 stands for fluid,
0 < φ < 1 describes the interface and φ = 1 represents theparticle domain. In this section, we briefly explain the salientfeatures of our method. The detailed formulas, algorithm, andapplicability of the SP method can be found in the literature[36,40–44].
The motion of the ith colloidal particle is obtained bysolving Newton’s equations of motion:Vcext
Mi V̇ i = F H
i + Fi + Fi + Gi ,Ṙi = V i

˙ i = N H + N ext + G  ,Ii · 
i
i
i

(1)
(2)

in which Ri is the position of the particle, and V i and iare the translational and rotational velocity of the particle,
respectively. The mass and moment of inertia are denotedby Mi and I i , respectively. The hydrodynamic torque andforce exerted by the solvent on the particle are representedH
extby N Hand N ext
are thei and F i , respectively. The F iiexternal force and torque, respectively. The G VandGii are therandom force and torque due to thermal fluctuations, which
can be described as G ni (t)G nj (0) = α n Iδ(t)δij in whichG ni (t) = 0 and α n (n ∈ V,) is the parameter that controls
the temperature of the system T . The actual value of the particle
temperature is determined using the long-time diffusion coefficient of the equilibrium system. When simulating a Brownian
particle using HIs, the diffusion coefficient is affected by
the finite size effects (i.e., artifacts arise when the systemsize is smaller than the correlation length). These effects canbe accounted for by D0 φ=0 = D0 K(φ), in which D0 φ=0 is thethermal diffusion coefficient of a Brownian particle at infinite
dilution, D0 is the thermal diffusion coefficient of Brownian
particles obtained for φ = 0, and K(φ) is the coefficientthat represents the effects of the finite volume fraction ofdispersed particles under periodic boundary conditions [45].
Finally, the temperature of the system can be determinedusing the Stokes-Einstein equation, kB T = 6π ηaD0 φ=0 . The
detailed implementation of Brownian motion in some test
cases is presented by Iwashita et al. [41]. Direct interparticle
interactions are denoted by F ci , and we represent theseinteractions as a truncated Lennard-Jones (LJ) potential withthe large powers of 24:12, which can be defined as follows:  24  12 
1
4 rσij+ rij  2 12 σ ,− rσijULJ (rij ) =
(3)


1
0
rij > 2 12 σ ,in which rij = |Ri − Rj |. The parameters σ and denote thelength and energy units of the LJ potential, respectively, withσ = 2a representing the particle diameter.

022310-2

DIRECT NUMERICAL SIMULATIONS OF ANISOTROPIC . . .

PHYSICAL REVIEW E 87, 022310 (2013)

In the SP method, the solvent motion is obtained by solving
the modified Navier-Stokes equation:
ρf (∂t v + v · ∇v) = −∇p + η∇ 2 v + ρf φ f p

3.5
3

(4)

Pe=0
Pe=5
Pe=9
Pe=24
Pe=47
Pe=94

2.5

FIG. 1. (Color online) Radial distribution function g(r) for
different Pe values. Simulations are performed in a cubic periodic
box of length L/a = 32 with a particle size of a = 4, keeping φ and
kB T constant at ≈0.08 and 0.3, respectively.

III. SIMULATION PARAMETERS

The Peclet number, which is the ratio of convection todiffusion, can be defined as follows:

A number of simulations were performed under periodic
boundary conditions with a volume fraction ranging from0.01 to 0.12. In the present study, we divide the range ofvolume fractions into a low volume fraction regime defined as
φ  0.04 and a moderate volume fraction regime defined as
0.06  φ  0.12. This division helps to explain the scaling of
the hydrodynamic velocity fluctuations and the self-diffusion
in an effective way, which will be explained later. Threedifferent cubic periodic boxes of dimensions L = 64 , 128 ,and 256 were used and the particle size was held constantat 4 for all simulations in which the grid spacingistaken as the unit of length. For given values of η and ρf ,the remaining units of mass, time, pressure, and energy aredefined as ρf 3 , ρf 2 /η, η2 /ρf 2 , and η2 /ρf , respectively.
The thermal fluctuations kB T were maintained at ∼0.3, and theof the LJ potential was set at 2.5. Gravity was introduced in
the z direction, and the Pe was increased by increasing the gravitational force. We held the particle-to-fluid-density ratio at 5 tohave a large range of Pe. The center-of-mass of the whole fluid
is constant to avoid an indefinite acceleration of the system.
In this study, the range of Pe is achieved by increasing
the value of gravity, which is introduced in the z direction.
The increase in gravity also increases the Reynolds number,
φ=0
Re = ρf aVs /η from 10−3 to 0.2, which is low enoughto ignore the inertial effects and remains within the Stokes
regime. In contrast to the simulations studies [18,19], mostof the experiments [1,2,4,32] have been performed at low Re
values, usually of the order of 10−5 or less. Because a relative
deviation from the Stokes regime scales with the square of Re,
we can probe relatively high Re values without straying fromthe relevant experimental conditions. In addition, maintaining
such a low value of Re is computationally more expensive.
For a given initial configuration, particles require a certain
time to acquire their steady-state velocities and positions.
We monitored the data, and only that corresponding to thesteady-state velocities is used. The simulations are run up to
500ts for Pe = 0.7 and to 1200ts for Pe > 50; sufficient toyield statistically meaningful data for the analysis in which ts
denotes the Stokes time.

g(r)for the total velocity field v and the pressure field p under theincompressibility condition (∇ · v = 0). The rigidity of theparticle and the fluid/particle nonslip boundary conditions areincorporated through the body force φ f p [36,40]. Moreover,an estimation of the error [46], a comparison of the SP
method with the force-coupling method (FCM) [47], and ahigh resolution DNS [48] based on high-order spectral/hp
element discretization on hybrid grids revealed that, comparedwith FCM and aforementioned DNS scheme, the SP method iscomputationally less expensive and advantageous for simulating moving particles because it avoids complex discretizations
around the particles. This study simulated many test cases andfound that the SP method accurately resolves the near-field
and far-field flows.

2
1.5
1
0.5

(φ=0.08, L/a=32)

0
1.5

2

2.5

3

3.5

r/a

φ=0

Pe =

aVs
.
D0

(5)

One can expect HIs to play a dominant role at higher Pe values,
thereby leading to changes in the microstructure. The HIs are
increased with an increase in 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 the radial distribution function (RDF),
2L3 g(r) = 2
δ(r − r ij ) ,
(6)
N
i<j
in which N is the total number of particles, r = |r|,r ij = Ri − Rj , · · · denotes an ensemble average, and thesummation i<j is taken over all particle pairs. The definitionof g(r) is such that 4πρg(r)r 2 r represents the mean numberof particles in a shell of radius r and thickness r thatsurrounds a particle at the origin [49]. Figure 1 shows the g(r)for different Pe values for a system of size L/a = 32 and φ =0.08. The peak of the function increases with increasing Pe,
demonstrating the formation of aggregates or particle clusters
in close contact. This clustering is induced by the dominanceof hydrodynamic forces with increasing Pe, which, in turn,progressively reduces the effects of random motion. A similar
phenomena has also been observed by Brady et al. [50] in theirinvestigation of the relative effects of thermal and shear forces.
They reduced the effects of thermal fluctuations by increasing
the shear rate and observed a similar cluster formation. Apartfrom the initial peak, the results are indistinguishable from theequilibrium results for Pe  9, demonstrating the strong effectof the thermal fluctuations. For Pe > 9, these results are differentiable, providing evidence that the HIs exert a progressively
dominating effect. A relatively large change in the RDF at Pe 
24 indicates that the HIs overpower the sedimentation phenomena. Moreover, the peak of g(r) at Pe = 0 is higher than thepeak value provided by the well-known exact solution [51,52]of the Percus-Yevick (PY) integral equation for the RDF.
In our results, this discrepancy arises from an overlapping

022310-3

ADNAN HAMID AND RYOICHI YAMAMOTO

PHYSICAL REVIEW E 87, 022310 (2013)
1

Pe=6
Pe=12
Pe=29
Pe=58
Pe=87
Pe=116Cx(t)/Cx(0)IV. RESULTS AND DISCUSSION
A. Hydrodynamic velocity fluctuations and their scaling

One of the main objectives of the current study is to developa scaling for the hydrodynamic velocity fluctuations and theirrelaxation times, which are two of the key parameters for
diffusion in sedimentation. To establish a scaling relationshipthat accounts for the effects of system size, volume fraction,and thermal fluctuations on the hydrodynamic velocity fluctuations, we calculate the temporal autocorrelation functions of
the velocity fluctuations as follows:

10-2
-3

10

0

5

10

15

0.2
0
0

10

20

30

40

50

t/ts
1

(φ=0.04, L/a=32)
0
10

(7)
0.8

(b)

)
Cz(t)/(VsCz(t)/Cz(0)(8)

in which δViz = Viz − Vsed , Viz and Vix are the temporalvelocities of the ith particle in the z and x direction,
respectively, and Vsed = Viz  is the mean settling velocityof the particles.
Figure 2 shows the time decay of the correlation functions
for the z and x components of the velocity fluctuations at
φ = 0.04 and L/a = 32. The main plots in Fig. 2 use the initial
fluctuations Cz (0) and Cx (0) as normalization constants tovisualize the relaxation of the velocity fluctuations, whereas
the figure insets use the square of the Stokes velocity tonormalize the correlations, highlighting the effects of thenonequilibrium hydrodynamic fluctuations. Time is normalφ=0ized by the Stokes time ts such that ts = a/Vs .
The main plots in Figs. 2(a) and 2(b) indicate that thelarge difference in relaxation times for the vertical andhorizontal correlations, whereas the experiments [1] show asmall difference in these time scales. The theory proposed by
Koch [33] suggests that this difference between the relaxation
times originates from the vertical periodic boundary conditions
because less time is required for a particle to sample all vertical
positions than is required to sample the horizontal positions.
This difference in time scales can be reduced by increasing
the aspect ratio of the periodic box or by increasing thevolume fraction of the particles. The simulations performedby Ladd [10] and Padding [18] observed the same phenomenausing vertical periodic boundary conditions.
At low Pe, the strong effects of the thermal fluctuations
cause a rapid decay in the correlations, whereas at higherPe, an exponential decay is evident as shown in Fig. 2.
In an experimental study, Nicolai et al. [1] observed an
exponential relaxation of the temporal correlations of velocityββfluctuations of the form Cβ (t) = ( VH )2 exp(−t/τH ) in whichββ
VH (β ∈ x,z) and τH denote the amplitude and relaxation
time of the hydrodynamic velocity fluctuations, respectively.
The present DNS results support the following forms for the
autocorrelations of the velocity fluctuations as the simulations
are performed at finite Pe:
2



(9)
Cx (t) = C0 (t) + VHx exp −t τHx ,




2
(10)
Cz (t) = C0 (t) + VHz exp −t τHz ,

)

0.4

and
Cz (t) = δViz (t)δViz (0),Eq. (9)φ=0 2

Cx (t) = Vix (t)Vix (0)

0.6

φ=0 2

0.8

(a)

-1

10

Cx(t)/(Vsof the smoothed profiles of two particles in close contact.
Similar discrepancies are also reported in the literature [18].

0.6
0.4

Eq. (10)
-1

10

-2

10
0.2

0

20

40

60

0
0

40

80

120

160

200

t/ts
FIG. 2. (Color online) Temporal autocorrelation functions of thevertical [Cz (t)] and horizontal [Cx (t)] velocity fluctuations. The
main plots show the correlation functions normalized by the initial
fluctuations, whereas the temporal correlation functions shown in
the insets are normalized by the square of the Stokes velocity
on a semilog scale. The insets also indicate that at large valuesof Pe, these correlation functions relax exponentially according toββCβ (t) = ( VH )2 exp(−t/τH ) as shown by the dotted lines in whichββ
VH (β ∈ x,z) and τH denote the amplitude and the relaxation
time of the hydrodynamic velocity fluctuations, respectively. Timeis normalized by the Stokes time as ts = a/Vsφ=0 .

in which C0 (t) represents the velocity autocorrelation functionin the presence of thermal fluctuations but without gravity,
which becomes negligible at higher Pe. These equations represent the summation of the pure thermal and hydrodynamic
forces; the former is dominant at low Pe, and the latter plays
a key role at high Pe. In addition to these two Pe regimes,
a transition regime is also expected to exist in which neitherof these two forces is dominant. Hence, the scaling relations
for the diffusion coefficients based on Eqs. (9) and (10) should
deviate from simulation results in this transition regime, which
will be tested later.
At finite Pe, the velocity fluctuations consist of a thermal
and a hydrodynamic component. The dominant role of thehydrodynamic fluctuations at higher Pe causes the correlation
functions to fall on the same curve when scaled with the squareof the Stokes velocity as suggested by Nicolai (see Fig. 2
insets). Our data demonstrate good agreement with the fitted
curves.
ββ
We obtained the values of VH and τH at different
volume fractions and system sizes from the exponential

022310-4

DIRECT NUMERICAL SIMULATIONS OF ANISOTROPIC . . .

1

0

10

τHx/ts

ΔVHx/Vsφ=0

10

-1

10

-2

10

10-2

φ

10-1

101

L/a=16
L/a=32
L/a=64

τHz/ts

ΔVHz/Vsφ=0

102

0

10

-1

10

-2

10

-1

φ

10

FIG. 3. (Color online) Scaling of the hydrodynamic velocity
fluctuations (left scale) and their relaxation times (right scale) as
a function of volume fraction, using three different system sizes:
βφ=0
L/a = 16, 32,√ and 64. The solid lines indicate a scaling to VH /Vsβof the form A1 Lφ/a,whereas the dashed lines indicate a scaling toββ√τH /ts of the form A2 L/(φa). The simulation results are represented
by points. The symbols and lines are color coded with respect to thesystem size.

fits. These values are plotted in Fig. 3 for three different
system sizes. Hinch [31] proposed scaling relationships for theβhydrodynamic velocity fluctuations ( VH ) and their relaxationβtimes (τH ). Cunha et al. [19] explained this scaling by
hypothetically dividing the simulation box into two equal parts.
The imbalance in the particle weight caused by variations in
the particle number in both parts is balanced by the Stokes
βdrag.
They predicted the scaling relations as VH /Vs φ=0 =
√
√
β
β
β
β
Aβ1 Lφ/a and τH /ts = A2 L/(φa) in which A1 and A2
are constants that depend on the system parameters (e.g., thecontainer shape, swirl size, particle shape, and polydispersity).
This simple scaling is expected to work well for low volume
fractions.
ββThe prefactors A1 and A2 are obtained by fitting thescaling relationships to the simulation data, as shown in Fig. 3,
which demonstrates the scaling of the hydrodynamic velocity
fluctuations and their relaxation times as a function of thevolume fraction for three different system sizes both in the x
and z direction. Our simulation results are in good agreementwith the scaling relations in the low volume fraction regime,
but deviations are evident in the moderate volume fraction for
all system sizes. Figure 3 also depicts that, for a given systemβsize, VH and τHx reach a saturation value at moderate volume

PHYSICAL REVIEW E 87, 022310 (2013)

fractions. The simulations of Padding et al. [18], Cunha et al.
[19], and the rationalized results [9] (considering the side wall
effects) of Segre et al. [4] also demonstrated similar scaling
at low volume fractions. Cunha showed the transition fromφ 1/2 , but unfortunately the scaling was unclear at high volume
fraction. In contrast, Brenner [9] theoretically proposed thetransition from φ 1/2 to φ 1/3 , which is not confirmed by anyother simulation or experimental study. In contrast to τHx , τHz
follows the φ −1/2 scaling for the full range of volume fraction.
Kuusela et al. [24] and Padding et al. [18] found a similar
scaling relation in their studies. This different decay patternin the vertical direction is attributed to enhanced memory
effects in the direction of flow. In contrast, experiments of
Nicolai et al. [1] did not find any clear decay pattern. It
is well known [9] that experiments are affected by many
nonideal conditions, such as secondary flow, instability, sidewall effects, and polydispersity. These nonideal conditions,especially polydispersity and side walls, greatly affect therelaxation time [9,14].
The velocity fluctuations and their relaxation times in both
directions scale as (L/a)1/2 for the full range of volume fractions studied, indicating that our simulations are in a spatiallycorrelated regime. The relevant scaling effects induced by thesize dependence of this finite system are discussed in detail in
Sec. IV C. The summary of the scaling is presented in Table I.
These findings at moderate volume fraction lead us to
split the scaling relationships into two: one for low volume
fractions (φ  0.04), similar to that observed in other studies[9,18,19], and a second for moderate volume fractions (0.06 φ  0.12) as shown in Eqs. (11)–(13):⎧⎨ Aβ Lφ/a,(φ  0.04)
1
β
φ=0
VH Vs
= 
(11)
⎩
V
β
A1 Lφsatβ /a, (0.06  φ  0.12)⎧ √⎨Ax2 L/(φa),(φ  0.04)
 , (12)
τHx ts =τ⎩Ax2 L/ φsatx a (0.06  φ  0.12)
and 
τHz ts = Az2 L/(φa), (φ  0.12),

(13)

V

in which φsatβ provides the saturation value of the hydrodynamic velocity fluctuations with respect to the volume fraction
in both the x and z directions at moderate volume fractions.
TABLE I. Summary of the scaling obtained for the hydrodynamic
velocity fluctuations ( VHβ /Vsφ=0 ) and their relaxation times (τHβ /ts )in both the vertical and horizontal directions. In addition to the aboveββscaling, both VH /Vsφ=0 and τH /ts are scaled as (L/a)1/2 with systemsize in the low to moderate volume fraction regimes.
Scaling of

VHx /Vsφ=0
VHz /Vsφ=0
τHx /ts
τHz /ts

022310-5

Low φ(φ  0.04)∼φ 1/2∼φ 1/2∼φ −1/2∼φ −1/2

ββ

VH /Vsφ=0 and τH /ts
Moderate φ(0.06  φ  0.12)
∼φ 0
∼φ 0
∼φ 0
∼φ −1/2

ADNAN HAMID AND RYOICHI YAMAMOTO

PHYSICAL REVIEW E 87, 022310 (2013)
TABLE II. A quantitative comparison of the fitting coefficientsobtained by fitting Eqs. (11)–(13) onto the simulation data. The
difference between the x- and z-direction coefficients indicates theanisotropic nature of the phenomena.

70
60
x

τH /tsEq. (12)z
τH /ts40βτH /ts50

30

Eq. (13)

20(L/a=32)

10
0
0

0.02

0.04

0.06

φ

0.08

0.1

0.12

FIG. 4. (Color online) Hydrodynamic velocity fluctuation relaxation times for both the vertical and horizontal directions as a functionof volume fraction. The dashed line indicates the scaling of τHx /ts andthe dot-dashed line indicates the scaling of τHz /ts . The simulation data
are represented by points. This plot shows that the difference in thevertical and horizontal relaxation times decreases with increasing
volume fraction at moderate volume fraction regimes. Our data show
good agreement with the scaling Eqs. (12) and (13).
τx
Furthermore, φsat
is the saturation value of the relaxation
time in the x direction only, as vertical relaxation times
obey the corresponding scaling for the full range of volume
fractions as evident in Eqs. (12) and (13). Similar phenomena
have been observed by Ladd [10] in a simulation study inwhich no change was found in the horizontal relaxation timeat φ = 0.05 and φ = 0.25, but the vertical relaxation timedecreased at these volume fractions. In the same study, Laddalso reported a saturation of the vertical relaxation times athigh volume fractions (φ  0.45). In an experimental study,
Bernard et al. [53] explored the dependency of the relaxation
time on the volume fraction and system size and determinedthat the relaxation time scaled as (L/a)0.65±0.1 , in agreementwith the (L/a)0.5 scaling observed in this study. Unlike thepresent study, Bernard found no systematic dependency of therelaxation time on the volume fraction.
The different scaling of the vertical and horizontal relaxation times at moderate volume fractions leads to a decrease in
the difference between the vertical and horizontal relaxation
times with increasing volume fraction as shown in Fig. 4.
This finding is important because it aids in characterizing theanisotropic nature of diffusion. Further discussion on the effectof this decrease is presented in Sec. IV B. Our data show goodconsistency with the scaling Eqs. (11)–(13).
ββThe prefactors A1 and A2 are determined from the fits tozz
be A1 ≈ 0.151, A2 ≈ 0.9857, Ax1 ≈ 0.0253, and Ax2 ≈ 0.239

Dz /D0 =Scaling coefficients

Values

Ax1
Ax2
Az1
Az2
Vx
φsat
V
φsatz
τx
φsat

0.02530.239
0.151
0.98570.05140.037
0.0526at low volume fractions. For moderate volume fractions, weVτxby fitting the scaling
extracted the values of φsatβ and φsat
Vz
Vx
relationships and obtained φsat = 0.037, φsat
= 0.05135, andτxφsat = 0.0526. The values of these fitting coefficients are also
summarized in Table II. The differences in the vertical andhorizontal prefactors indicate the anisotropic behavior of thevelocity fluctuations, with vertical fluctuations larger than thehorizontal ones. This anisotropy in the velocity fluctuations
varies from 2.5 to 4 in the simulations [10,12,19,21]. The
experiments [1,4] have provided a value of ≈2.5. Notably,
the ratio of the vertical to horizontal hydrodynamic velocity
fluctuations in this study is VH z / VH x ≈ 2.44 at low
volume fractions, and 2.07 at moderate volume fractions,which are in good agreement with the experimental [1,4]and simulation [21] results. This anisotropic behavior in thevelocity fluctuations is believed to be due to the asymmetry of
the system induced by gravity.
B. Scaling and anisotropic behavior of diffusion

Diffusion occurs due to the fluctuating motion of particles.
Individual particles lose the memory of their velocity after
experiencing hydrodynamic interactions with surrounding particles and thereafter follow a random-walk diffusion process.
To find the scaling of a long-time steady-state self-diffusion
coefficient, we consider the total diffusion coefficient D as
the sum of the thermal contribution D0 and the hydrodynamic
contribution DH , where the latter is estimated asββ 2 βDH ≈ VH τH .
(14)
βBased on the aforementioned scaling relationships for VHβand τH , for both low and moderate volume fraction regimes,
we can also find the scaling for the diffusion coefficients in thez and x directions:
1 + Az1 Az2 Pe (L/a)3/2 φ 1/2 = 1 + 0.149Pe (L/a)3/2 φ 1/2 ,V
1 + Az1 Az2 φsatz Pe (L/a)3/2 φ −1/2 = 1 + 0.0056Pe (L/a)3/2 φ −1/2 ,(φ  0.04)(0.06  φ  0.12)


(φ  0.04)
1 + Ax1 Ax2 Pe (L/a)3/2 φ 1/2 = 1 + 0.006Pe (L/a)3/2 φ 1/2 ,Dx /D0 =
.


−1/2
Vx
τx
φsat
Pe (L/a)3/2 = 1 + 0.0013Pe (L/a)3/2 , (0.06  φ  0.12)
1 + Ax1 Ax2 φsat

022310-6

,

(15)

(16)

DIRECT NUMERICAL SIMULATIONS OF ANISOTROPIC . . .
TABLE III. Summary of the scaling obtained for the selfdiffusion coefficients (Dβ /D0 ) in both the vertical and horizontal
directions. Additionally, diffusion coefficients in both directions scaleas (L/a)3/2 with system size and linearly with Pe.

PHYSICAL REVIEW E 87, 022310 (2013)

directions as follows:

Dz (t) =Scaling of Dβ /D0
Low φ(φ  0.04)∼φ 1/2∼φ 1/2

Dx /D0Dz /D0Moderate φ(0.06  φ  0.12)∼φ 0∼φ −1/2

1
(Rix (t) − Rix (0))2 ,
2t

(17)

1
(Riz (t) − Riz (0) − Vsed t)2 .
2t

(18)

Dx (t) =The equilibrium self-diffusion coefficient can be obtainedusing the Einstein relationship:
D0 = lim

1t→∞ 6t

These scaling relations indicate that the vertical andhorizontal diffusion coefficients increase linearly with Pe, butwith a smaller prefactor in the horizontal direction. In addition,
they scale as (L/a)3/2 with the system size for the full range ofvolume fractions that we have considered. Furthermore, boththe vertical and horizontal diffusion coefficients scale with φ 1/2
at low volume fraction, but at moderate volume fraction the
diffusion coefficient in the vertical direction scales as φ −1/2 ,and the horizontal diffusion coefficient becomes independentof volume fraction. This additional dependency on φ −1/2
in the vertical direction causes a decrease in the diffusion
anisotropy with increasing volume fraction. The scaling
relations for the diffusion coefficients are summarized inTable III.
The long-time steady-state self-diffusion coefficient [Dβ =limt→∞ Dβ (t)] is obtained from the linear growth of the meansquare displacement (MSD) in the horizontal and vertical(19)

Figure 5 shows a comparison of the long-time steadystate self-diffusion coefficients with the scaling relations.
Figures 5(a) and 5(b) show the diffusion coefficient as afunction of Pe in the horizontal direction for low and moderate
volume fractions, respectively. Figures 5(c) and 5(d) show
the diffusion coefficient in the vertical direction as a functionof Pe for low and moderate volume fractions, respectively.
Figures 5(a) and 5(c) indicate that the diffusion coefficientincreases with increasing volume fraction both in the horizontal and vertical directions at low volume fractions. In
contrast with the low volume fraction regime, the vertical
diffusion coefficient decreases with increasing volume fraction
at moderate volume fractions as shown in Fig. 5(d). Figure 5(b)demonstrates no change in the horizontal diffusion coefficient.
This decrease in the vertical diffusion coefficient is attributed
to the φ −1/2 scaling of τHz . However, due to the saturation of τHx
and VHx at moderate volume fraction, no change is observedin the horizontal diffusion coefficient.

2

2

10

10

(L/a=32)

(a)

Dx/Do

Dx/Do

(Ri (t) − Ri (0))2 .

1

10

100 -1
10

100

101

(b)

1

10

100 -1
10

102

100

Pe

101

102

Pe

3

3

10

10

(c)

Dz/Do

102

Dz/Do

102

(d)

1

10

1

10

φ=0.01
φ=0.02φ=0.040

10

φ=0.06
φ=0.08φ=0.10φ=0.12

0

-1

10

0

1

10

10

2

10

Pe

10

-1

10

0

1

10

10

2

10

Pe

FIG. 5. (Color online) A comparison of the long-time steady-state self-diffusion coefficients in both the vertical and horizontal directions
with those predicted by the scaling Eqs. (15) and (16) as a function of Pe. The self-diffusion coefficients are normalized by the equilibrium
diffusion coefficient D0 . (a) and (b) show the scaling of Dx /D0 for low and moderate volume fractions, respectively, whereas (c) and (d) show
the scaling of Dz /D0 for low and moderate volume fractions, respectively. The simulation results are represented by points. The lines indicate
the scaling relations. The symbols and lines are color coded with respect to the volume fraction.
022310-7

ADNAN HAMID AND RYOICHI YAMAMOTO

PHYSICAL REVIEW E 87, 022310 (2013)
2

10

25
(L/a=32)

15

Dz/Dx

Dz/Dx

20

10
φ=0.01
φ=0.02φ=0.045

φ=0.06
φ=0.08φ=0.101

10Scaling
Present Study
Nicolai [1]Padding [18]φ=0.120

Cunha [19]Ladd [12]Kuusela [24]

0

0

20

40

60

80

100

10

120

Pe
FIG. 6. (Color online) The anisotropy of the vertical (Dz ) andhorizontal (Dx ) self-diffusion coefficients as a function of Pe fordifferent volume fractions. The solid lines indicate the scaling of
Dz /Dx at low volume fractions, whereas the dashed lines indicate thescaling at moderate volume fractions. The scaling lines are obtainedby dividing Eq. (15) with Eq. (16) for the respective volume fraction
regime. Lines and points of the same color indicate the results for thesame volume fraction.

Our simulation results indicate that the diffusion coefficients increase linearly with increasing Pe, both parallel andperpendicular to gravity. A deviation from the scaling relations
is evident at low Pe. Because scaling relations are derived byadding pure thermal and hydrodynamic forces, these relations
are expected to show good agreement with the simulation
results when either thermal or hydrodynamic fluctuations
dominate the phenomena. Hence, this deviation indicates thetransition regime in which neither force is dominant. Our
simulations show good agreement with Eqs. (15) and (16)at high Pe, demonstrating the dominance of the hydrodynamic
interactions.
The majority of studies [1,4,10,19,23,53] on particle diffusion in sedimentation have focused on the non-Brownian
regime. Particle diffusion and its anisotropic nature at finite
Pe values have not yet been explored. We have attempted toinvestigate this anisotropic behavior and found that, for a given
volume fraction, the anisotropy increases with increasing Pe,
becoming saturated at higher Pe as shown in Fig. 6. Thisobservation suggests that the effect of thermal fluctuations is
significant at low Pe before the HIs begin to dominate thesedimentation phenomenon. We can predict the anisotropic
behavior of the diffusion at finite Pe with scaling relationships.
Our data show good agreement with the predicted diffusion
anisotropy as indicated in Fig. 6. Furthermore, the scaling
relationships suggest that the saturated or steady-state value ofthe anisotropy remains unchanged at low volume fractions anddecreases with increasing volume fraction at moderate volume
fractions as shown in Fig. 6. This decrease in anisotropy is
attributed to the decrease in the difference between the verticaland horizontal relaxation times. This difference is decreasedby the φ −1/2 scaling of the vertical relaxation time at moderate
volume fractions. The horizontal relaxation time is independent of volume fraction in this regime as shown in Table I.
Our simulation results are also compared with the previous
experimental and simulation results in Fig. 7, which shows
the diffusion anisotropy as a function of the volume fraction.

0

0.02

0.04

0.06

φ

0.08

0.1

0.12

FIG. 7. (Color online) Comparison of the diffusion anisotropy
with the scaling and previously published experimental and simulated
data. Nicolai [1] obtained the low anisotropy value in an experimental
study in a spatially screened regime, whereas Cunha [19] and Padding[18] achieved these results from the DNS and SRD using a periodic
box with an aspect ratio of 3. A relatively large value of anisotropy isachieved due to cubic periodic box and spatially unscreened regime.
Our data show good agreement with the simulation study of Kuusela
et al. [24] and the scaling.

Nicolai [1] obtained these results in a spatially screened
regime, whereas Cunha [19] and Padding [18] performed
DNS and SRD simulations in an elongated simulation box
with an aspect ratio of 3. We have obtained a relatively highvalue of the anisotropy at a low volume fraction. This higher
value originates from the difference in vertical and horizontal
relaxation times. Because diffusivity is the product of VH2
and τH , a larger difference in relaxation time leads to a higher
diffusion anisotropy. The diffusion anisotropy can be reduced
by increasing the aspect ratio of the simulation box [33].
An increase in the aspect ratio reduces this difference andthus the diffusion anisotropy, as exhibited by Cunha [19] andPadding [18]. In addition, increasing the aspect ratio reduces
the prefactors of Eqs. (15) and (16), especially in the verticaldirection [33], which in turn reduces the diffusion anisotropy.
Hence, diffusion anisotropy is not only dependent on volume
fraction, it also depends on the aspect ratio of the cell. We
obtained an anisotropy of ≈24, whereas Cunha and Padding
found this anisotropy to be ≈10 and 7, respectively, at low
volume fractions. In another simulations study [24], Kuusela
et al. obtained large diffusion coefficients in both directions
using a 2D periodic square simulation box, but the ratio ofdiffusion coefficients shows good agreement with our results.
In contrast with the present study, Nicolai [1] obtained ananisotropy of ≈7 and found no decrease in the anisotropy atmoderate volume fractions, primarily because they obtained
their results in a large container, which is many times greaterthan the correlation length. In simulations, the use of such alarge system is impossible due to the large computational cost.
Another possible explanation for the discrepancy between thesimulations and the experimental results is the presence of aside wall in the experiments, which creates microstructural
inhomogeneities over time [9,54]. In addition, the presence ofpolydispersity, even to a small extent, can temper the diffusion.
Ladd [12] obtained a high anisotropy ≈77 at φ = 0.05. This

022310-8

DIRECT NUMERICAL SIMULATIONS OF ANISOTROPIC . . .

large value is attributed to the small system size and the use offull periodic boundary conditions. Similar to the present study,
Ladd [12] demonstrated that diffusion anisotropy decreases
with increasing volume fraction. The anisotropy of diffusion,
which was ≈77 at φ ≈ 0.05, decreased to ≈39 at φ ≈ 0.107.
We have many experimental and simulation studies to comparewith at low volume fractions, but unfortunately, we could not
find a study which had addressed the scaling of diffusion
anisotropy with respect to volume fraction at moderate ortransition regime.
C. Finite size effects

Theoretical arguments [7] and simulations [12] have longexhibited the strong dependency of velocity fluctuations on
the system size for a random suspension of particles. In
contrast, experiments [1–3,5] show no such divergence. Thisdisagreement was solved by Segre et al. [4] who suggested
that the velocity fluctuations increase with the system size only
when the system size is smaller than the correlation length, andabove this correlation length, the simulations and experiments
should be in good agreement. Experimentation has the freedom
to use large system sizes, whereas in simulations, large systems
require enormous resources. Thus, most simulation studies[10,18,19,21] are affected by the artifacts induced by thefinite system size. To observe the system size dependency
in our results, we defined the spatial correlation function of
the velocity fluctuations, which defines the flow pattern of theparticles as follows:2L3 Cz (r) = 2
δViz δVj z δ(r − r ij ) .
(20)
N
i<j
We define Cz (z) and Cz (x), with respect to the distancevector r, in either the vertical r = zδ z or horizontal r = xδ x
direction.
We plotted the spatial correlation function of the z component of the velocity as a function of the distance perpendicular
to gravity in Fig. 8 using three different cubic boxes L/a =16,
32, and 64, at φ = 0.10. Figure 8 shows that the velocity
1

Cz(x)/Cz(0)0.8

L/a=16
L/a=32
L/a=640.6
0.4(φ=0.10)

0.2
0
-0.2
0

5

10

15

20

25

x/a
FIG. 8. (Color online) Spatial correlation function of the z
component of velocity as a function of the distance perpendicular
to gravity, demonstrating the effects of system size on velocity
fluctuations. All simulations were performed at a volume fractionof φ ≈ 0.10 and Pe ≈ 80.

PHYSICAL REVIEW E 87, 022310 (2013)fluctuations become long ranged with an increase in the systemsize, which is in agreement with the tendencies observedin similar simulation studies [18,24]. Our simulation results
support the experimental findings of Segre [4] when the systemsize is smaller than the correlation length. Therefore, theparticle motion should be affected by the finite system inour simulations. Figure 8 shows that a large system size isneeded to reproduce the saturation of the velocity fluctuations,but this requires enormous computational time and resources.
Therefore, most of the computational studies [10,11,14,18,19]
are similarly limited. Systematic experimental [2,4] andsimulation [11,13] studies, however, have revealed that theeffects of finite system size on velocity fluctuations can beeffectively explained using the concept of finite-size scaling[18]. We believe that there is value in performing critical tests
on the scaling concept. In addition, the experiments [2,4] andsimulations [11,13] addressing the effects of finite system sizeon velocity fluctuations have also established that finite systemsize may lead to smaller velocity fluctuations [18].
V. CONCLUSION

We studied steady-state sedimentation at a finite Peclet
number using the direct numerical simulations with thesmooth-profile method in a spatially unscreened regime.
The present study focused on the scaling of hydrodynamic
velocity fluctuations and self-diffusivities with respect tovolume fraction and system size. The study also examined
the relative effects of thermal and hydrodynamic fluctuations.
We observed a clear transition from a Brownian-motiondominant regime to a hydrodynamic-fluctuations-dominant
regime. Moreover, the hydrodynamic velocity fluctuations increased with increasing volume fraction at low volume fraction(φ  0.04) in accordance with theoretical predictions [7,9]and previous studies [18,19], whereas at moderate volume
fraction (0.06  φ  0.12), their behavior was independentof φ, similar to Nicolai et al. [1] and Climent et al. [25].
We concluded that the amplitudes of the velocity fluctuation
correlations scale with the square of the Stokes velocity at largePe, with vertical hydrodynamic velocity fluctuations that are≈2.5 times larger than the horizontal hydrodynamic velocity
fluctuations at low volume fractions, in good agreement withother studies [1,18,19,25]. At moderate volume fractions, thisvalue drops to ≈2. In addition, we tested the scaling relations
for the hydrodynamic velocity fluctuations suggested by Hinch[31] for different system sizes and volume fractions. We foundthat this scaling worked well for low volume fractions anddeviated at moderate volume fraction regimes.
Consistent with the theoretical predictions [7,9], we foundthat both the vertical and horizontal velocity fluctuations scaleas (φL/a)1/2 at low volume fractions, saturating with respectto φ at moderate volume fractions. Similarly, the horizontal
velocity fluctuation relaxation time scales as (L/aφ)1/2 at low
volume fractions, saturating with respect to φ at moderate
volume fractions. This result is in contrast with that ofthe vertical relaxation time, which scales as (L/aφ)1/2 forboth regimes. Padding et al. [18] and Kuusela et al. [24]observed a similar decay behavior for vertical relaxation time.
Furthermore, we found that the difference in the vertical andhorizontal hydrodynamic velocity fluctuation relaxation times

022310-9

ADNAN HAMID AND RYOICHI YAMAMOTO

PHYSICAL REVIEW E 87, 022310 (2013)

decreased with increasing volume fraction at moderate volume
fractions.
Based on the scaling of the hydrodynamic velocity fluctuations and their relaxation times, we inferred the scaling
relations of the long-time steady-state self-diffusion coefficient for both low and moderate volume fraction regimes.
Theoretical arguments [7,9,19] suggest that both the verticaland horizontal self-diffusion coefficients scale as (L/a)3/2 φ 1/2
at low volume fractions. Our results have shown the similar
behavior. In contrast to the low volume fraction regimes, thevertical self-diffusion coefficient scales as (L/a)3/2 φ −1/2 atmoderate volume fractions, due to a decrease in the verticalrelaxation time with volume fraction, whereas the horizontal
diffusion coefficient is saturated with respect to φ in thisregime. Our simulation results found good consistency withthe scaling relations. These relations in hydrodynamic velocity
fluctuations, their relaxation times, and the self-diffusion
coefficients are summarized in Tables I and III. The scaling of
the diffusion coefficients allows us to predict the anisotropyof the vertical and horizontal diffusion. We found that the
diffusion anisotropy increases with increasing Pe, saturating
at high Pe, in accordance with Ref. [18]. This saturated value
remains unchanged at low volume fractions and decreases withincreasing volume fraction at moderate volume fractions. The

decrease in the anisotropy with volume fraction is induced by
the decrease in the difference in magnitudes of the verticaland horizontal hydrodynamic relaxation times. We obtained a
relatively high anisotropy at low volume fraction. This higher
value is attributed to the use of a cubic periodic box, as
an increase in aspect ratio reduces the diffusion anisotropy[18,19,33].
The present simulations were performed maintaining
Re < 0.2, so that the inertial effects are negligible. Extensive
simulations are underway to study the effects of inertial
forces at moderate Re values. Additionally, we are attempting
to simulate even higher volume fractions for non-Brownian
particles to study the behavior of velocity fluctuations andparticle diffusion.

[1] H. Nicolai, B. Herzhaft, E. J. Hinch, L. Oger, and E. Guazzelli,Phys. Fluids 7, 12 (1995).
[2] H. Nicolai and E. Guazzelli, Phys. Fluids 7, 3 (1995).
[3] H. Nicolai, Y. Peysson, and E. Guazzelli, Phys. Fluids 8, 855(1996).
[4] P. N. Segre, E. Helbolzheimer, and P. M. Chaikin, Phys. Rev.
Lett. 79, 2574 (1997).
[5] P. N. Segrè, F. Liu, P. Umbanhowar, and D. A. Weitz, Nature(London) 409, 594 (2001).
[6] P. N. Segre, Phys. Rev. Lett. 89, 254503 (2002).
[7] R. E. Caflish and J. H. C. Luke, Phys. Fluids 28, 759 (1985).
[8] D. L. Koch and E. S. G. Shaqfeh, J. Fluid Mech. 224, 275 (1991).
[9] M. P. Brenner, Phys. Fluids 11, 754 (1999).
[10] A. J. C. Ladd, Phys. Fluids A 5, 299 (1993).
[11] A. J. C. Ladd, Phys. Rev. Lett. 76, 1392 (1996).
[12] A. J. C. Ladd, Phys. Fluids 9, 491 (1997).
[13] A. J. C. Ladd, Phys. Rev. Lett. 88, 048301 (2002).
[14] N. Q. Nguyen and A. J. C. Ladd, J. Fluid Mech. 525, 73 (2005).
[15] A. J. C. Ladd and R. Verberg, J. Stat. Phys. 104, 1191 (2001).
[16] J. T. Padding and A. A. Louis, Phys. Rev. Lett. 93, 220601(2004).
[17] J. T. Padding and A. A. Louis, Phys. Rev. E 74, 031402 (2006).
[18] J. T. Padding and A. A. Louis, Phys. Rev. E 77, 011402 (2008).
[19] F. R. Cunha, G. C. Abade, A. J. Sousa, and E. J. Hinch, J. Fluids
Eng. 124, 957 (2002).
[20] S. Y. Tee, P. J. Mucha, M. P. Brenner, and D. A. Weitz, Phys.
Fluids 19, 113304 (2007).
[21] M. C. Miguel and R. Pastor-Satorras, Europhys. Lett. 54, 45(2001).
[22] Y. Xiaolong and D. L. Koch, Phys. Fluids 20, 043305 (2008).
[23] E. Kuusela, J. M. Lahtinen, and T. Ala-Nissila, Phys. Rev. E 69,
066310 (2004).

[24] E. Kuusela and T. Ala-Nissila, Phys. Rev. E 63, 061505 (2001).
[25] E. Climent and M. R. Maxey, Int. J. Multiphase Flow 29, 579(2003).
[26] G. K. Batchelor, J. Fluid Mech. 52, 245 (1972).
[27] G. G. Stokes, Mathematical and Physical Papers (Johnson
Reprint Corporation, New York, 1966).
[28] J. M. Ham and G. M. Homsy, Int. J. Multiphase Flow 14, 533(1988).
[29] R. H. Davis, J. Fluid Mech. 310, 325 (1996).
[30] A. Levine, S. Ramaswamy, E. Frey, and R. Bruinsma, Phys.
Rev. Lett. 81, 5944 (1998).
[31] E. J. Hinch, in Disorder and Mixing, edited by E. Guyon,
Y. Pomeau, and J. P. Nadal (Kluwer Academic, Dordrecht,1988).
[32] E. Guazzelli and J. Hinch, Annu. Rev. Fluid Mech. 43, 97 (2011).
[33] D. L. Koch, Phys. Fluids 6, 2894 (1994).
[34] J. W. Swan, J. F. Brady, V. N. Michailidou, and G. Petekidis,
Phys. Rev. Lett. 109, 068302 (2009).
[35] P. J. Hoogerbrugge and J. M. V. A. Koelman, Europhys. Lett.
19, 155 (1992).
[36] Y. Nakayama and R. Yamamoto, Phys. Rev. E 71, 036707 (2005).
[37] H. Tanaka and T. Araki, Phys. Rev. Lett. 85, 1338 (2000).
[38] A. Hamid and R. Yamamoto, J. Phys. Soc. Jpn. 82, 024004(2013).
[39] A. Hamid and R. Yamamoto, AIP Conf. Proc. 1518, 444 (2013).
[40] Y. Nakayama, K. Kim, and R. Yamamoto, Eur. Phys. J. E 26,
361 (2008).
[41] T. Iwashita, Y. Nakayama, and R. Yamamoto, J. Phys. Soc. Jpn.
77, 074007 (2008).
[42] T. Iwashita and R. Yamamoto, Phys. Rev. E 79, 031401 (2009).
[43] Y. Nakayama, K. Kim, and R. Yamamoto, Adv. Powder Technol.
21, 206 (2010).

ACKNOWLEDGMENTS

This work was supported by KAKENHI 23244087 and
the Global COE Program, “International Center for Integrated
Research and Advanced Education in Materials Science” from
the Ministry of Education, Culture, Sports, Science and Technology of Japan. The author would like to extend his gratitude
to Dr. Takahiro Murashima, Dr. Hideyuki Mizuno, and Dr.
John Molina for their useful discussions and suggestions.

022310-10

DIRECT NUMERICAL SIMULATIONS OF ANISOTROPIC . . .
[44] KAPSEL website, http://www-tph.cheme.kyoto-u.ac.jp/kapsel/.
[45] A. A. Zick and G. M. Homsy, J. Fluid Mech. 115, 13(1982).
[46] X. Luo, R. Martin, M. R. Maxey, and G. E. Karniadakis, J.
Comp. Phys. 228, 1750 (2009).
[47] M. R. Maxey and B. K. Patel, Int. J. Multiphase Flow 27, 1603
(2001).
[48] R. M. Kirby, T. C. Warburton, I. Lomtev, and G. E. Karniadakis,
Appl. Numer. Math. 33, 393 (2000).

PHYSICAL REVIEW E 87, 022310 (2013)[49] D. C. Rapaport, in The Art of Molecular Dynamics Simulation,
2nd ed. (Cambridge University Press, Edinburgh, 2004),
pp. 90–91.
[50] G. Bossis and J. F. Brady, J. Chem. Phys. 87, 5437 (1987).
[51] M. S. Wertheim, Phys. Rev. Lett. 10, 321 (1963).
[52] E. Thiele, J. Chem. Phys. 39, 474 (1963).
[53] G. Bernard-Michel, A. Monavon, D. Lhuillier, D. Abdo, and
H. Simon, Phys. Fluids 14, 2339 (2002).
[54] Y. Peysson and E. Guazzelli, Phys. Fluids 10, 44 (1998).

022310-11