86

Progress of Theoretical Physics Supplement No. 178, 2009

Velocity Autocorrelation Function of Fluctuating Particles
in Incompressible Fluids
Toward Direct Numerical Simulation of Particle Dispersions
Takuya Iwashita,1 Yasuya Nakayama2 and Ryoichi Yamamoto1,3
1 Department of Chemical Engineering, Kyoto University, Kyoto 615-8510, Japan
2 Department of Chemical Engineering, Kyushu University,

Fukuoka 819-0395, Japan

3 CREST, Japan Sience and Technology Agency, Kawaguchi 332-0012, Japan

Motions of fluctuating Brownian particles in an incompressible viscous fluid have beenstudied by coupled simulations of Brownian particles and host fluid. We calculated the velocity autocorrelation functions of Brownian particles and compared them with the theoreticalresults. Extensive discussions have been made on the time scales for which our numericalmodel is valid.

§1.

Introduction

The motions of small ﬂuctuating particles in viscous ﬂuids have been studied
for a long time. Although theoretical or numerical analysis based on the coupledmotions of the particles and the host ﬂuid are very complicated, it becomes rather
simple if one considers only the particles’ motions by assuming that the host ﬂuiddegree of freedom can be safely projected out from the entire degree of freedom ofthe dispersions. One of such models is the well-known generalized Langevin equation(GLE) for Brownian particles, i.e.,
 t

dVi
=ds
Γij (t − s)Vj (s) + Gi (t),(1.1)
Midt−∞jGi (t) · Gj (0) = 3kB T Γij (t),(1.2)

where Mi and Vi denote the mass and the translational velocity of the i-th particle,
respectively. Γij (t) is a friction tensor, which represents the eﬀect of hydrodynamic
interactions (HI) between i-th and j-th particles. Gi is the random force acting on
the i-th particle induced by thermal ﬂuctuations of the solvent, kB is Boltzmann
constant, and T is the temperature of Brownian particles.
For a single spherical particle (i = 1) immersed in a inﬁnitely large host ﬂuid,
the analytic form of the time-dependent friction1) is known as(1.3)
dtΓ11 (t) exp(−iωt) = Γ̂11 (−iω) = 6πηa(1 + a −iω/ν − iωa2 /9ν) ,where Γˆ11 (−iω) is the Fourier transform of Γ11 (t) and ω is the angular frequency. The
ﬁrst term corresponds to the normal Stokes friction for a spherical particle of radiusa in a Newtonian ﬂuid whose viscosity is η. The second term represents the memory

Velocity Autocorrelation Functions of Fluctuating Particles

87

eﬀect, which is related to the momentum diﬀusion in a viscous medium. Here thekinematic viscosity is deﬁned as ν = η/ρf with ρf being the density of the ﬂuid. The
third term corresponds to the eﬀect of the acceleration of the host ﬂuid surroundingthe tagged particle when the particle is accelerated through the host ﬂuid. Using
Eq. (1.3), the hydrodynamic GLE can be solved analytically. The translational
velocity autocorrelation function (VACF) Vi (t) · Vi (0)/3 then obtained is known toexhibit the characteristic power-law relaxation for long-time region, which is widely
known as the “hydrodynamic long-time tail”.2)–5)
For dispersions composed of many particles interacting via HI, the situation isstill not straightforward because we do not know the true analytic expression for thehydrodynamic friction tensor Γij (t). Some approximated expressions, such as Oseen
or Rotne-Prager-Yamakawa (RPY) tensor, can be obtained by introducing the Stokes
approximation. However, those expressions completely neglect the memory eﬀect
that corresponds to the second term of Eq. (1.3). This means that the hydrodynamic
long-time tail cannot be reproduced correctly with Oseen or RPY tensor.
In the present study, we developed a numerical method to take into accountthe eﬀects of hydrodynamics directly by simultaneously solving the Navier-Stokes
equation for the host ﬂuid with the Brownian motions of the particles. We ﬁrst
examined the VACF for a single Brownian particle and compared it with the analytical form mentioned above. Secondly, we examined the rotational motions of a single
Brownian particle. We furthermore examined the motions of Brownian particles inharmonic potentials to check the validity of our method.
§2. Simulation method

Here we brieﬂy explain the basic equations of our numerical model since thoseare explained in detail elsewhere.6) A smooth proﬁle function 0 ≤ φ(x, t) ≤ 1 isintroduced to deﬁne ﬂuid (φ = 0) and particle (φ = 1) domains on a regular Cartesiangrid. Those two domains are separated by thin interface regions whose thickness isξ. The position of the i-th particle is Ri , the translational velocity is Vi , and therotational velocity is Ωi . The motion of i-th particle with mass Mi and the momentof inertia Ii is governed by the following Langevin-type equations,dVi= FiH + FiC + Fiex + GVi ,dtdΩi= NiH + GΩIi ·i ,dt
MidRi= Vi ,dt(2.1)(2.2)

where FiH and NiH are the hydrodynamic forces and torques acting on the i-th
particle due to HI, respectively. FiC and Fiex denote the direct particle-particle
interaction and external force. GVi and GΩ
i are the random force and torque due tothermal ﬂuctuations deﬁned stochastically as
GVi  = GΩi  = 0,
GVi (t) · GVj (0) = 3kB T αV δ(t)δij ,(2.3)
(2.4)

ΩΩ
GΩ
i (t) · Gj (0) = 3kB T α δ(t)δij ,(2.5)

88

T. Iwashita, Y. Nakayama and R. Yamamoto

where αV and αΩ are parameters to control the temperature T .
The motions of the host ﬂuid are governed by the Navier-Stokes equation
ρf (∂t v + v · ∇v) = −∇p + η∇2 v + ρf φfp(2.6)

with the incompressible condition ∇ · v = 0, where v and p are the velocity and thepressure ﬁelds of the host ﬂuid, respectively, and φfp is the body force deﬁned sothat the rigidity of the particles is automatically satisﬁed. Note that FiH and NiH
are determined from the body force φfp .7), 8)§3.

Results and discussion

A single spherical particle ﬂuctuating in a Newtonian ﬂuid was simulated in
the absence of external forces Fiex = 0 as depicted in Fig. 1. We take the meshsize Δ and τ = Δ2 ρf /η as the units of space and time. Simulations have been
performed with η = 1, a = 5, and ξ = 2 in a three-dimensional cubic box composedof 64 × 64 × 64 grid points. The particle and ﬂuid densities are identically set to beunity, ρp = ρf = 1.
Figure 2 shows our simulation results () of VACF for a single Brownian particle ﬂuctuating in a Newtonian host ﬂuid at kB T = 0.83. The temperature T wasdetermined by comparing the long-time diﬀusion coeﬃcient Dsim obtained from simulations with DV = kB T /6πηaK(Φ), where K(Φ) takes into account the eﬀects ofﬁnite volume fraction9) and Φ denotes the volume fraction. The volume fraction of asingle particle is Φ = 0.002. One ﬁnds that the VACF approaches asymptotically tothe power-law line with the exponent −3/2, and the long-time behavior of our simulation agrees well with the analytical solution5) of the hydrodynamic GLE rather
than the Markovian VACF which neglects memory eﬀects. This behavior indicates

ZYX-0.005-0.002

0.002

0.005

Fig. 1. A snapshot of a single Brownian particle immersed in a Newtonian fluid. The one eighthof the entire system is graphically displayed. The color map on the horizontal plane shows thevalue of the local fluid velocity in the x-direction.

Velocity Autocorrelation Functions of Fluctuating Particles

89-2< Ω i (t)・Ω i (0)> /310-310<Vi(t)・Vi(0)>/3,10-4

10-5<Vi (t)・Vi (0)>/3-1.5

Bt
translational analytical solution
Markovian VACF
< Ω i (t)・Ω i (0)> /3-6

10-2.5Ct
rotational analytical solution
Markovian RVACF

10-710-1100

t

101

102

Fig. 2. The translational velocity autocorrelation function Vi (t) · Vi (0)/3 (triangle) and the rotational velocity autocorrelation function ˙i (t) · ˙i (0)/3 (circle) for a single Brownian particle
fluctuating in a Newtonian fluid. The simulation data was taken at kB T = 0.83. The solidlines indicate the analytic results for the translational5) and the rotational motions. The dotted lines show power-laws, Bt−3/2 with B = kB T /12ρf (πν)3 for the translational motions andCt−5/2 with C = πkB T /32ρf (πν)5/2 for the rotational motions. The dashed lines indicate
the Markovian VACF and RVACF, which decay exponentially as exp(−t/τB ) and exp(−t/τr ),respectively.

that the memory eﬀects are accurately taken into account. Similar to the translational motions, we have studied the rotational motions of the Brownian particle in
the host ﬂuid. The GLE of the rotational motions for a single spherical particle canbe written as tdsμ(t − s)Ωi (s) + Gi (t),(3.1)
Ii Ω̇i = −−∞Gi  = 0,

Gi (t) · Gi (0) = 3kB T μ(t),(3.2)

wherefriction μ(t) has the form μ̂(−iω) = 8πηa3 [1 − iω/3ν(1 + the time-dependent10)
in Fourier space. The ﬁrst term in μ̂ is the Stokes friction and thea −iω/ν)]second term represents the memory eﬀect due to the kinematic viscosity of the ﬂuid.
The GLE can be solved analytically, and the analytical solution of the rotational
velocity autocorrelation function (RVACF) is obtained in the following form:Ωi (t) · Ωi (0)3kB T ν= −
8πηa5

 ∞0

y 3/2
dyexp(−yt/τν ), (3.3)
3π[1 − ( ττr + 13 )y]2 + y(1 − ττr y)2νν

90

T. Iwashita, Y. Nakayama and R. Yamamoto

where τν = a2 /ν and τr = Ii /8πηa3 . In Fig. 2, simulation results () of RVACF
are also plotted. The RVACF clearly shows the asymptotic approach to the hydrodynamic long-time tail with the exponent −5/2 which agrees well with the analytical solution (3.3) rather than a simple Markovian RVACF. By comparing thepresent simulation results with the corresponding analytical solutions more in detail, one may notice that some discrepancies become notable for t < τB or t < τr ,where τB = Mi /6πηa = 2a2 ρp /9η  5 is the Brownian relaxation time and τr =Ii /8πηa3 = 3τB /10  1.5 is the Brownian rotational relaxation time. For opposite
cases t > τB or t > τr , however, the agreements between the numerical results andthe analytical solutions are excellent. This is because we neglected memory eﬀects
in thermal noises GVi and GΩ
i . We however believe that the long-time behavior of
our numerical model is valid for t > τB since τB is much longer than the memory
times of the thermal noises.
There exist many other characteristic time-scales in particle dispersions. Important ones are the kinematic time-scale τν = a2 ρf /η = 25 which measures the momentum diﬀusion over the particle size and the diﬀusion time-scale τD = a2 /D  3 × 103
which measures the particle diﬀusion over the particle size. As one can see in Fig. 2,
the present model works quite well for the time-scales comparable to τν and τD ,while it becomes inaccurate for t < τB .
In order to test the validity of our method for the long-time behavior of Brownian particles, we next applied the present model to simulate Brownian particles
ﬂuctuating in external harmonic potentials. The potentials are introduced with theformFi ex = −k(Ri − Rieq ) = −kΔRi ,(3.4)

where Rieq is the i-th particle’s equilibrium position and k is the spring constant.
Figure 3 shows the positional autocorrelation function ΔRi (t)·ΔRi (0)/3 of two
Brownian particles in harmonic potentials whose minimum positions are separated
by a ﬁxed distance of 5a. The pair of particles are interacting only hydrodynamically,
and there exists no direct interactions between them. The spring constant is set tok = 10, and the temperature is kB T  0.0066, which was determined by the averagepotential energy kB T = kΔRi2 /3. The simulation results () agree well withthe hydrodynamic analytical solution5) in harmonic potentials which account for
the eﬀects of ﬁnite volume fraction. The analytical solution was derived by solving
the GLE of a single Brownian particle in a harmonic potential which includes themodiﬁed Stokes friction ζ = 6πηaK(Φ). The correlation functions decay much slower
than the Markovian relaxation functions. We also conﬁrmed that the validity of ourmethod is excellent for t ≥ τν .
§4.

Conclusion

We proposed a numerical model to simulate Brownian particles ﬂuctuating in
Newtonian host ﬂuids. To test the validity of the model, the translational velocity
autocorrelation function (VACF), the rotational velocity autocorrelation function

Velocity Autocorrelation Functions of Fluctuating Particles

91

0.007
Hydrodynamic analytical solution
present simulation
Markovian

0.006<ΔRi (t)・ΔRi (0)>

0.005
0.004
0.003
0.002
0.001
0-0.001050

100

150200t

Fig. 3. The positional autocorrelation function ΔRi (t) · ΔRi (0)/3 (circle) for a system composedof two Brownian particles in harmonic potentials −k(Ri − Rieq ). The distance between theirequilibrium positions is |R1eq − R2eq | = 5a. The solid line shows the analytical solution.5) The
dotted line shows the Markovian functions, i.e. the solution of Mi V̇i = −6πηaVi − kRi + Gi .

(RVACF), and the positional autocorrelation function of ﬂuctuating Brownian particles were calculated in some simple situations for which analytical solutions wereobtained. We compared our numerical results with the analytical solutions and foundexcellent agreements between them specially for long-time regions t > τB while somediscrepancies were found for short time regions t < τB . This is because our model isdesigned to simulate correct long-time behaviors of Brownian particles in host ﬂuids.
Applications of the present method for more complicated situations are in progress.
Acknowledgements
This work was partly supported by a grant from Hosokawa powder technology
foundation.
References
1) E. M. Lifshitz and L. D. Landau, Fluid Mechanics (Addison-Wesley, Reading, 1959).
2) B. J. Alder and T. E. Wainwright, Phys. Rev. A 1 (1970), 18.
3) A. Widom, Phys. Rev. A 3 (1971), 1394.
4) E. H. Hauge and A. Martin-Löf, J. Stat. Phys. 7 (1973), 259.
5) H. J. H. Clercx and P. P. J. Schram, Phys. Rev. A 46 (1992), 1942.
6) T. Iwashita, Y. Nakayama and R. Yamamoto, J. Phys. Soc. Jpn. 77 (2008), 074007.
7) Y. Nakayama and R. Yamamoto, Phys. Rev. E 71 (2005), 036707.
8) Y. Nakayama, K. Kim and R. Yamamoto, Eur. Phys. J. E 26 (2008), 361.
9) A. A. Zick and G. M. Homsy, J. Fluid Mech. 115 (1982), 13.
10) S. H. Lamb, Hydrodynamics (Dover Publications, New York, 1932).