PHYSICAL REVIEW E 85, 066704 (2012)

Direct numerical simulation of dispersed particles in a compressible fluid
Rei Tatsumi* and Ryoichi Yamamoto†Department of Chemical Engineering, Kyoto University, Kyoto 615-8510, Japan(Received 26 April 2012; published 12 June 2012)
We present a direct numerical simulation method for investigating the dynamics of dispersed particles in acompressible solvent ﬂuid. The validity of the simulation is examined by calculating the velocity relaxation of
an impulsively forced spherical particle with a known analytical solution. The simulation also gives information
about the ﬂuid motion, which provides some insight into the particle motion. Fluctuations are also introducedby random stress, and the validity of this case is examined by comparing the calculation results with theﬂuctuation-dissipation theorem.
DOI: 10.1103/PhysRevE.85.066704

PACS number(s): 47.11.−j, 82.70.Dd

I. INTRODUCTION

In the case of a dispersed particle of radius a = 100 nm inwater, for instance, the compressibility factor is evaluatedfrom ν = 1.0 × 10−6 m2 /s and c = 1.5 × 103 m/s to beε = 6.7 × 10−3 . In many cases, researchers are interested
in phenomena progressing over the time scale of viscous
diffusion or even longer time scales, such as those relating to shear properties, electrophoresis, and sedimentation.
Therefore, the assumption of an incompressible host ﬂuid isvalid, and most direct numerical simulation methods, including
SPM, have been developed on the premise of incompressible
ﬂuids. However, when we investigate phenomena associatedwith sound propagation, such as sonic agglomeration, acoustic
spectroscopy, and electroacoustic measurements, the consideration of compressibility is required.
In the present study, we extend the SPM to compressible
ﬂuids. Some aspects of the dynamics of a single particle
in a compressible ﬂuid have been theoretically analyzed[3–7], and we compare the simulation results obtained hereinwith analytical solutions to determine the accuracy of thesimulation. In particular, we consider the velocity relaxationof a spherical particle after an impulsive force is added. The
numerical simulation also gives information regarding ﬂuid
motion for which the analytical solution is unknown, and the
dynamics of the particle can be investigated from the viewpointof the ﬂuid dynamics. In addition, we also consider a systemwith thermal ﬂuctuations by introducing random stress, and thevelocity autocorrelation function is compared with the analytical solution according to the ﬂuctuation-dissipation theorem.

Particle dispersions have various unique properties, and anunderstanding of these properties is important in many ﬁeldsof science and engineering. These properties originate fromthe dynamics of particles, which are extremely complicatedbecause of the hydrodynamic interactions among particles
mediated by the motion of the surrounding ﬂuid. Therefore,
several numerical approaches have been formulated to investigate the dynamics of such dispersions. As one of theseapproaches, a direct numerical simulation has been developed,
wherein the hydrodynamic interactions are directly computedby simultaneously solving for the motion of the ﬂuid and the
motion of the particle. In recent years, we have developed an
efﬁcient direct numerical simulation scheme for dispersions,which is called the smoothed proﬁle method (SPM) [1,2].
In this scheme, the sharp interface between the ﬂuid andthe particles is replaced by a smoothed interface using acontinuous proﬁle function. The Navier-Stokes equations are
solved for the ﬂuid motion on a ﬁxed square grid, and Newton’sand Euler’s equations of motion for the particles are solved
simultaneously while considering the momentum exchange
between the ﬂuid and the particles.
The hydrodynamic interactions are transmitted in two ways:
via viscous momentum diffusion and via sound propagation.
The time scale of viscous diffusion over the particle size isτν = a 2 /ν, and that of sound propagation is τc = a/c, wherea is the particle radius, ν is the kinematic viscosity, and cis the speed of sound in the ﬂuid. Here, we will deﬁne thecompressibility factor as the ratio of the two time scales:
ν
τc
=
ε=
.
(1)
τν
ac
This factor represents the degree of inﬂuence of compressibility in the dynamics of the dispersed particles. According to
Eq. (1), the compressibility becomes increasingly importantas the particle size decreases. Indeed, the compressibilityof a system has been considered an important factor inmolecular scale dynamics of monatomic liquids studied usinghydrodynamic theories [3,4]. On the other hand, in studies ofparticle dispersions, compressibility has rarely been considered, and an incompressible host ﬂuid is generally assumed.

II. SIMULATION METHOD
A. Equations

In the SPM, the particle-ﬂuid boundary is replaced by acontinuous interface. For this purpose, the smoothed proﬁlefunction φ(r,t) ∈ [1,0] is introduced. This function represents
the boundary between the ﬂuid and particle regions: φ = 1
at the particle domain, and φ = 0 at the ﬂuid domain. The
two regions are smoothly connected through thin interfacial
regions with a thickness characterized by ξ . The mathematicalexpression of φ can be found in a previous paper [1]. The total
velocity ﬁeld is deﬁned asv = (1 − φ)v f + φv p ,*†tatsumi@cheme.kyoto-u.ac.jp
ryoichi@cheme.kyoto-u.ac.jp1539-3755/2012/85(6)/066704(6)(2)

where v f (r,t) is the ﬂuid velocity ﬁeld and v p (r,t) is theparticle velocity ﬁeld constructed from rigid motions of the066704-1©2012 American Physical Society

REI TATSUMI AND RYOICHI YAMAMOTO

PHYSICAL REVIEW E 85, 066704 (2012)

particles [1,2]. When a compressible host ﬂuid is considered,
the ﬂuid mass density is altered, and we deﬁne the mass density
ﬁeld asρf = (1 − φ)ρ.

(3)

The auxiliary mass density ﬁeld ρ(r,t) is deﬁned over theentire domain. However, the physical ﬂuid mass densityρf (r,t) must be zero within the particle domain, and thisrequirement is satisﬁed by multiplying by (1 − φ).
The equations governing the dynamics of the dispersion
system are given as hydrodynamic equations with the additionof a body force. The hydrodynamic equations consist of threeconservation laws concerning mass, momentum, and energy.
The conservation equations of mass and momentum aregiven as∂ρ
+ ∇ · m = 0,
(4)
∂t
∂m
(5)
+ ∇ · (mv) = ∇ · σ + ρφ f p ,∂twhere m(r,t) = ρ(r,t)v(r,t) is the momentum density ﬁeld.
We consider a compressible Newtonian ﬂuid, and the stresstensor is given by
σ = −p I + η[∇v + (∇v)T ] + ηv − 23 η (∇ · v)I,(6)
where p(r,t) is the pressure, η is the shear viscosity, and ηv
is the bulk viscosity. A body force ρφ f p is added so that therigidity of the particles is satisﬁed. Additionally, we assume
a barotropic ﬂuid described by p = p(ρ), and the pressure
gradient is proportional to that of the mass density:∇p = c2 ∇ρ.

(7)Equations (4)–(7) are closed for the variables ρ, m, and p;
therefore, energy conservation does not need to be considered.
The motion of the dispersed particles is governed by
Newton’s and Euler’s equations of motion:ddCRi = V i ,(8)
Mi V i = f H
i + Fi ,
dt
dt
d
(9)
I i · i = N Hi ,dtwhere Ri , Vi , and i are the position, translational velocity,and rotational velocity of the ith particle, respectively. The
particle has a mass Mi and a moment of inertia I i . The
H
hydrodynamic force F H
i and torque N i are exerted on
the particle by the ﬂuid, and the force F Ci is exerted throughdirect interactions among the particles.
The effect of thermal ﬂuctuations on the dynamics ofparticles is important when the particle size is on the order of a
micrometer or smaller. Fluctuations are introduced througha random stress tensor s, which is added in Eq. (6). The
random stress is a stochastic variable satisfying the ﬂuctuationdissipation relation [8]:sij (r,t)skl (r  ,t  ) = 2kB T ηij kl δ(r  − r)δ(t  − t),(10)where kB is the Boltzmann constant, T is the temperature, and(11)ηij kl = η(δik δj l + δil δj k ) + ηv − 23 η δij δkl .

introduced using the Langevin approach, where random forces
are exerted on the particles [9]. However, this approach does
not accurately represent the short-time dynamics of the systembecause the time autocorrelation of the hydrodynamic force
acting on the particles is neglected. Therefore, the ﬂuctuatinghydrodynamics approach is more appropriate for investigating
dynamics over the time scale of sound propagation.
B. Simulation procedure

Here, the time-discretized evolution of the equations isderived. The time step tn represents the nth discretized time,
and the time step change from tn to tn+1 = tn + h will beconsidered. The time evolution of the system is determinedthrough the following steps.
(i) The mass and momentum density changes associatedwith sound propagation are calculated as tn +hρ n+1 = ρ n −dt∇ · m,(12)tn∗

m = m −cn2

 tn +h
dt∇ρ.

(13)

tn

When we assume a periodic boundary condition and use theFourier spectral method, a semi-implicit scheme becomesfeasible [10]. This situation eases the restriction on the timeincrement for a small compressibility factor ε.
(ii) The time evolution of the advection and viscous
diffusion terms is calculated as tn +hm∗∗ = m∗ +dt∇ · (τ − mv),
(14)tn

where τ is the dissipative stress given by σ = −p I + τ .
(iii) In concert with the advection of the particle domain,
the position of each dispersed particle evolves as= Rni +
Rn+1i

 tn +hdt V i .

(15)tn(iv) The hydrodynamic force and torque are derived by considering the conservation of momentum. The time-integratedhydrodynamic force and torque are computed as tn +h(16)
dt F H=d rφin+1 m∗∗ − ρ n+1 v np ,i
 tn +htntn

dt N H
i =





× φin+1 m∗∗ − ρ n+1 v np .
d r r − Rn+1i
(17)With these and other forces acting on the particles, thetranslational and rotational velocity of each dispersed particle
evolve as tn +h


n+1
−1
n
C
(18)
dt F HV i = V i + Mi
i + Fi ,Brownian motion of the dispersed particles is induced byrandom stress acting on the ﬂuid. Thermal ﬂuctuations can be
066704-2tn

n+1= ni + I −1ii ·

 tn +htn

dt N Hi .

(19)

DIRECT NUMERICAL SIMULATION OF DISPERSED . . .

PHYSICAL REVIEW E 85, 066704 (2012)(v) The updated velocity of the particle region is imposedon the velocity ﬁeld as the body force ρφ f p .
 tn +hn+1∗∗=m +dtρφ f p ,
(20)
m

tn

tn



dtρφ f p = φ n+1 ρ n+1 v n+1− m∗∗ .
p

10

10

0

-1

(21)

In the case of an incompressible ﬂuid, the pressure isspontaneously determined by the solenoidal condition of thevelocity ﬁeld. On the other hand, in the present case, thepressure or mass density varies independently of the velocity
ﬁeld.

(t)

 tn +h

(a)

10

10

-2

= 0.1
= 0.6
analytical solution (incompressible)
analytical solution ( = 0.1)
analytical solution ( = 0.6)

-3

III. SIMULATION RESULTS

A. Relaxation

Let us consider a single spherical particle in a ﬂuid at rest.
We will investigate the relaxation of the particle velocity after
exerting an impulsive force at the center of the particle. The
impulsive force is assumed to be small; in other words, wewill consider the motion of the particle and the ﬂuid for low
Reynolds and Mach numbers. We set the impulsive force toachieve an initial particle Reynolds number of Rep = 0.01.
Here, we introduce the velocity relaxation function as thenormalized velocity change of the particle asPγ (t),(22)
M
where P is the impulsive force added at t = 0. The analyticalform of the relaxation function is obtained from a Stokes
approximation (see the Appendix).
V (t) =

10

(b)

10

10

| (t)|

Numerical simulations are performed for a threedimensional box with periodic boundary conditions. The spaceis divided by meshes of length , which is the unit length.
The units of the other physical quantities are deﬁned by
combining η = 1 and ρ0 = 1 with , where ρ0 is the ﬂuid massdensity at equilibrium. The system size is Lx × Ly × Lz =128 × 128 × 128. The other parameters are set as a = 4,
ξ = 1, ρp = 1, ηv = 0, and h = 0.01, where ρp is the particle
mass density. We performed simulations of the dynamics ofan isolated spherical particle in a ﬂuid in two situations.
First, we investigated the relaxation response of a particle
with an impulsive force. Second, we consider the velocity
autocorrelation function of a particle with thermal ﬂuctuations.
Because the input particle radius a = 4 is not necessarilyequal to the effective hydrodynamic radius a ∗ of the particle
represented by the smoothed proﬁle function φ with a fuzzy
interface of thickness ξ , we calculated a ∗ from the drag force
acting on the spherical particle moving at velocity V , which isanalytically given by FD = 6π ηaK(φ)V , where K(φ) represents the effect of the periodic boundary condition dependingon the volume fraction of the particle φ [11]. The effective
radius was evaluated as a ∗ = 3.87 in the present simulations,and therefore, the momentum diffusion time is estimated tobe τν = a ∗2 /ν in this case. When we compare the present
simulation results with the analytical solutions, the correctedparticle density ρp∗ = (a ∗ /a)3 ρp and the compressibility factorε∗ = (a/a ∗ )ε of the analytical solutions are employed.

10

10

10

-3

10

-2

10

-1

t/

10

0

10

1

0

-1

-2

= 1.0
= 1.5
analytical solution (incompressible)
analytical solution ( = 1.0)
analytical solution ( = 1.5)

-3

-4

10

-3

10

-2

10

-1

t/

10

0

10

1

FIG. 1. (Color online) The velocity relaxation function for various compressibility factors: (a) ε = 0.1, 0.6, (b) 1.0, and 1.5. The
bold lines illustrate the analytic solutions for various compressibility
factors, and the bold dashed double-dotted line shows the analytic
solution for an incompressible ﬂuid.

The simulations are performed with compressibility factorsof ε = 0.1, 0.6, 1.0, and 1.5. The condition ε > 1 implies thatthe sound propagation occurs more slowly than the viscous
diffusion, according to the deﬁnition given in Eq. (1). The
simulation results agree quite well with the analytical solutions
for all of the compressibility factors ε, as shown in Fig. 1. The
oscillation of t/τν  2 at ε = 0.1 arises from the periodic
boundary conditions. In this case, the sound pulse arrives atthe end of the system at t/τν = 1.6; afterward, the sound pulse
returns and affects the particle motion.
The velocity relaxation results highlight some remarkableproperties of dynamics in a compressible ﬂuid. In the caseof an incompressible ﬂuid, part of the particle momentum isinstantly carried away by the propagation of the inﬁnite-speed
sound wave, and the particle moves as if its mass wereM ∗ = M + Mf /2, where Mf = 4π a 3 ρ0 /3 is the mass of

066704-3

REI TATSUMI AND RYOICHI YAMAMOTO(a) ε = 0.1

t /τν = 0.13(b) ε = 1.5

t /τν = 0.13

PHYSICAL REVIEW E 85, 066704 (2012)

x10-4
2
1.5
1
0.5
0
-0.5
-1
-1.5
-2

x10-4

t /τν = 0.27

0.6
0.4
0.2
0
-0.2
-0.4
-0.6

x10-3
1.5
1

x10-4

t /τν = 0.53

x10-3

t /τν = 0.27

2
1.5
1
0.5
0
-0.5
-1
-1.5
-2

0.5
0
-0.5
-1
-1.5

0.2
0.15
0.1
0.05
0
-0.05
-0.1
-0.15
-0.2

x10-3

t /τν = 2.67

1.5
1
0.5
0
-0.5
-1
-1.5

FIG. 2. (Color online) Time evolution of the ﬂuid density deviation around the particle. The compressibility factors are (a) ε = 0.1 and(b) ε = 1.5. The color scale (gray scale) represents negative (darker) to positive (lighter) density deviation. The black circle represents theparticle. The direction of the initial particle velocity is right in these pictures.

the displaced ﬂuid [12]. Therefore, the velocity relaxation
function for an incompressible ﬂuid at the initial time isγ (+0) = M/M ∗ ; afterward, the particle velocity gradually
decreases because of momentum diffusion in the ﬂuid causedby its viscosity. Eventually, this decay obeys the power law
t −3/2 in the long-time region. On the other hand, in the case of
a compressible ﬂuid, the velocity relaxation function indicatesthat the relaxation due to sound propagation occurs in a ﬁnite
time interval. For the compressibility factor ε < 1, the tworelaxation processes are almost separate, and the relaxation
function coincides with that of an incompressible ﬂuid after
the relaxation due to the sound propagation occurs. With anincrease in the compressibility factor, nonmonotonic behavior
is observed in the relaxation function, and ﬁnally, inversion of
the particle velocity is observed, as shown in the investigationof the analytical solution [7].
The density deviation around the particle is shown in Fig. 2,
for which an analytical form has not been obtained. For asmall compressibility ε = 0.1, the sound wave pulse expands
from the particle very quickly, which corresponds to a rapid
decrease in the particle momentum due to sound propagation.
On the other hand, for a large compressibility ε = 1.5, thesound pulse does not spread. The pulse remains within thevicinity of the particle and gradually decays through viscous
diffusion. The continuation of high ﬂuid density in front of theparticle is expected to cause backtracking motion.
The velocity ﬁeld around the particle is also affected by
the compressibility. The velocity ﬁeld can be formally decomposed into two components: an incompressible (solenoidal)
component v I and a compressible component v C . The former
is a solenoidal vector ﬁeld when ∇ · v I = 0, and the latter
is an irrotational vector ﬁeld when ∇ × v C = 0. Therefore,
the divergence of the velocity ﬁeld ∇ · v gives information
regarding the compressible component, and the rotation of

the velocity ﬁeld ∇ × v gives information pertaining to theincompressible component. The rotation and divergence ofthe velocity ﬁeld are shown in Figs. 3 and 4, respectively.
The velocity ﬁeld is also depicted in these ﬁgures. An
obvious vortex ring is observed around the particle for a small
compressibility ε = 0.1, while the vortex ring is not clear fora large compressibility ε = 1.5. The rotation corresponds tothe intensity of the vorticity, and the map in Fig. 3 represents
a pair of vortex rings that is inherent to an incompressible
ﬂuid [13]. The rotation remains nearly constant, regardless ofthe compressibility, which indicates that the incompressible
component is not affected by the compressibility. The vectorﬁelds v I and v C inﬂuence each other only through the nonlinear
terms in Eqs. (4) and (5); however, in the present simulations,a low Reynolds number ﬂow is assumed, and the additivity(a) ε = 0.1(b) ε = 1.5
0.4
0.3
0.2
0.1
0
-0.1
-0.2
-0.3
-0.4

0.4
0.3
0.2
0.1
0
-0.1
-0.2
-0.3
-0.4

FIG. 3. (Color online) The ﬂuid velocity ﬁeld and its rotation,∇ × v, at t/τν = 0.27. The values are normalized by the factor
Ma ∗ /| P|. The component normal to the ﬁgure plane is depicted. The
compressibility factors are (a) ε = 0.1 and (b) ε = 1.5. The color
scale (gray scale) represents negative (darker) to positive (lighter)rotation. The black circle represents the particle. The direction of theinitial particle velocity is right in these pictures.

066704-4

DIRECT NUMERICAL SIMULATION OF DISPERSED . . .

(b) ε = 1.5
0.08
0.06
0.04
0.02
0
-0.02
-0.04
-0.06
-0.08

10

0

0.15
0.1
0.05

M V·V(t) /3kBT

0
-0.05

›

(a) ε = 0.1

PHYSICAL REVIEW E 85, 066704 (2012)

-0.1

-1

‹

-0.15

10

FIG. 4. (Color online) The ﬂuid velocity ﬁeld and its divergence,∇ · v, at t/τν = 0.27. The values are normalized by the factor
Ma ∗ /| P|. The compressibility factors are (a) ε = 0.1 and (b) ε = 1.5.
The color scale (gray scale) represents negative (darker) to positive(lighter) divergence. The black circle represents the particle. The
direction of the initial particle velocity is right in these pictures.

of the incompressible and compressible components is almost
valid. Therefore, the effect of compressibility is observed only
from the divergence of the velocity ﬁeld depicted in Fig. 4. The
region of positive divergence represents the source of the ﬂow,
and the negative region represents the sink of the ﬂow. Fromthe equation of continuity Eq. (4), the divergence is equal to
the reverse sign of the time change rate of the mass density:∇·v = −

D
ln ρ.
Dt

(23)

According to this relation, the source corresponds to thedensity decrease, and the sink corresponds to the density
increase. Therefore, the source and the sink move according
to the sound propagation. The propagation speed decreases
with increasing compressibility, which corresponds to thedifference in Fig. 4 related to the compressibility. The patternof the total velocity ﬁeld is described as the superposition of
the vortex convection of the incompressible component v I andthe source-sink ﬂow of the compressible component v C . The
compressibility factor governs the relative time evolution ofeach component to produce various ﬂow patterns.
B. Fluctuation

Thermal ﬂuctuations are introduced through random stressin the host ﬂuid. Computationally, a random stress term
satisfying Eq. (10) is added to the stress tensor given by
Eq. (6). We consider a single particle in the ﬂuid, whichmoves randomly as a result of thermal ﬂuctuations. From theﬂuctuation-dissipation theorem, the velocity autocorrelation
function of the particle is related to the relaxation function as
γ (t) =MV (0) · V (t).
3kB T

(24)

The accuracy of the ﬂuctuating system can be conﬁrmed by
the validity of this relation. In the numerical procedure forsolving a stochastic differential equation, the numerical error
can be larger than that in an ordinary differential equation due
to the truncation error in the time integration of the random
noise term [14]. This error is decreased with the decrease in
the time increment, and we set the smaller time increment
than that in the relaxation case as h = 0.01. The system sizeis therefore scaled down to Lx × Ly × Lz = 64 × 64 × 64 to

10

= 0.1
analytical solution (incompressible)
analytical solution ( = 0.1)

-2

10

-3

10

-2

10

-1

t/

10

0

10

1

FIG. 5. (Color online) The velocity autocorrelation function at
ε = 0.1 and kB T = 10−4 . The bold solid line represents the analytic
solution of the velocity relaxation function. The bold dashed doubledotted line shows the analytic solution for an incompressible ﬂuid.

compensate the increased computational demand due to thesmall h explained above. The simulation results for ε = 0.1
are shown in Fig. 5, which shows good agreement with theanalytical solution of the relaxation function. The consistency
between the input and calculated temperatures is also tested.
We can evaluate the temperature from the average kinetic
energy of a ﬂuid or a particle, i.e., kB Tf = 3 ρv 2 /3 orkB Tp = MV 2 /3. We evaluated these two temperatures as
the ratio to the input temperature T . From the ﬂuid motion,
the temperature was evaluated as Tf /T = 1.06, while thatof the particle motion was Tp /T = 0.93. The overestimationof the ﬂuid temperature Tf is simply due to the truncation error
in the time integration of the stochastic differential equation.
On the other hand, the particle temperature Tp is slightly below
the input temperature T . This discrepancy is considered to be
due to the small numerical inaccuracy introduced in treating
the momentum transfer through the particle-ﬂuid interface
using the SPM. Further improvements on the treatments of
ﬂuctuations in the particle-ﬂuid interfacial region will be
discussed in the future.
IV. CONCLUSION

We extended the SPM to particle dispersions in compressible ﬂuids. The validity of the method was conﬁrmed
by calculating the velocity relaxation function of a single
spherical particle in a compressible ﬂuid. The effect ofcompressibility on the velocity relaxation was also observed,
showing two-stage relaxation in a low-compressibility ﬂuid
and backtracking motion in a high-compressibility ﬂuid.
These particle motions were considered by investigating theﬂuid density deviation. The propagation of the sound pulse
around the particle is governed by the compressibility, and theinﬂuence of the sound disappears in a low-compressibility ﬂuid
but is maintained in a high-compressibility ﬂuid. The effect ofcompressibility on the ﬂuid velocity ﬁeld was also observed,
which was essentially understood to arise from changes in thetime evolution of the source-sink ﬂow component caused by
the compressibility.

066704-5

REI TATSUMI AND RYOICHI YAMAMOTO

PHYSICAL REVIEW E 85, 066704 (2012)

A simulation of the motion of a single spherical particle ina ﬂuctuating ﬂuid was also performed. The calculated velocity
autocorrelation function of the particle showed good agreement with the analytical solution of the relaxation function,and the validity of the ﬂuctuation-dissipation theorem was
conﬁrmed.

surface of the sphere, the memory kernel is obtained from thelinearized hydrodynamic equations (Stokes approximation)as [5,6]
 ∞ζ (t)eiωt dt,(A2)ζ̂ (ω) =−∞(1 + x)(9 − 9iy − 2y 2 ) + x 2 (1 − iy)
4π
ηax 2
,
=
3
2x 2 (1 − iy) − (1 + x)y 2 − x 2 y 2(A3)

ACKNOWLEDGMENTS

The authors would like to express their gratitude to
Dr. S. Yasuda, Dr. H. Kobayashi, and Dr. Y. Nakayama for theiruseful comments and discussions. This work was supported by
Grant No. KAKENHI 23244087 and the JSPS Core-to-Core
Program “International research network for non-equilibrium
dynamics of soft matter.”withx = a(−iωρ0 /η)1/2 , y = aω/c̃,
and




iω 4
c̃ = c −η + ηv
ρ0 3

1/2

2

APPENDIX: VELOCITY RELAXATION OF A PARTICLE
IN A COMPRESSIBLE FLUID

The equation of motion for a spherical particle under anexternal force E(t) in a compressible ﬂuid can be expressedby the linear response theory as
 tddsζ (t − s)V (s) + E(t),(A1)
M V =−
dt
−∞

(A4)

.

(A5)

According to Eq. (A1), the particle velocity is linearly
dependent on the external force in frequency representation:ˆV̂ (ω) = (ω)Ê(ω),(A6)

ˆ
where the admittance (ω)
is given byˆ(ω)= [−iωM + ζ̂ (ω)]−1 .

(A7)where ζ (t) is the memory kernel of the friction force.
The memory kernel is analytically expressed in frequency
representation. Assuming a stick boundary condition on the

In the case that an impulsive force is exerted on the sphere, theexternal force in frequency representation is given as a constantvector Ê(ω) = P, and the velocity relaxation function is givenby γ (t) = M(t) according to Eq. (22).

[1] Y. Nakayama and R. Yamamoto, Phys. Rev. E 71, 036707(2005).
[2] Y. Nakayama, K. Kim, and R. Yamamoto, Eur. Phys. J. E 26,
361 (2008).
[3] R. Zwanzig and M. Bixon, Phys. Rev. A 2, 2005 (1970).
[4] H. Metiu, D. W. Oxtoby, and K. F. Freed, Phys. Rev. A 15, 361(1977).
[5] D. Bedeaux and P. Mazur, Physica 78, 505 (1974).
[6] B. U. Felderhof, Phys. Fluids 19, 126101 (2007).
[7] B. U. Felderhof, J. Chem. Phys. 123, 044902 (2005).

[8] L. D. Landau and E. M. Lifshitz, Fluid Mechanics (Pergamon,
London, 1959).
[9] T. Iwashita, Y. Nakayama, and R. Yamamoto, J. Phys. Soc. Jpn.
77, 074007 (2008).
[10] G. Erlebacher, M. Y. Hussaini, H. O. Kreiss, and S. Sarkar,
Theor. Comput. Fluid Dyn. 2, 73 (1990).
[11] H. Hashimoto, J. Fluid Mech. 5, 317 (1959).
[12] R. Zwanzig and M. Bixon, J. Fluid Mech. 69, 21 (1975).
[13] B. U. Felderhof, Phys. Fluids 19, 073102 (2007).
[14] R. L. Honeycutt, Phys. Rev. A 45, 600 (1992).

066704-6