PHYSICS OF FLUIDS 25, 046101 (2013)

Propagation of hydrodynamic interactions between
particles in a compressible fluid
Rei Tatsumia) and Ryoichi Yamamotob)
Department of Chemical Engineering, Kyoto University, Kyoto 615-8510, Japan(Received 21 November 2012; accepted 11 March 2013; published online 23 April 2013)

Hydrodynamic interactions are transmitted by viscous diffusion and sound propagation, and the temporal evolution of hydrodynamic interactions by both mechanisms
is studied using direct numerical simulation in this paper. The hydrodynamic interactions for a system of two particles in a fluid are estimated using the velocity
correlation of the particles. In an incompressible fluid, hydrodynamic interactions
propagate instantaneously at the infinite speed of sound followed by a temporal evolution due to viscous diffusion. Conversely, sound propagates in a compressible fluid
at a finite speed, which affects the temporal evolution of the hydrodynamic interactions through an order-of-magnitude relationship between the time scales of viscous
diffusion and sound propagation. The hydrodynamic interactions are characterized
by introducing the ratio of these time scales as an interactive compressibility factor.
C 2013 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4802038]


I. INTRODUCTION

In particle dispersions, the motion of each particle in a fluid solvent affects the motion of theother particles. Such dynamical interactions are called hydrodynamic interactions; these interactions
produce complex dynamical behavior for dispersions that can be observed in particle aggregationand sedimentation phenomena.
Hydrodynamic interactions correspond to momentum exchange among particles through theambient fluid and are transmitted by two mechanisms: viscous diffusion and sound propagation.
These two mechanisms occur at different time scales. The time scale for viscous diffusion over a
distance equal to the particle size is τ ν = a2 /ν, while the time scale for sound propagation is τ c= a/c, where a is the particle radius, ν is the kinematic viscosity, and c is the speed of sound in thefluid. To study dynamical effects at the scale of particle size, the relative significance of the soundpropagation mechanism in hydrodynamic interactions is assessed from the ratio of the two timescales described above:
ν
τc
.
(1)
=
ε=
τν
ac
This dimensionless quantity is called the compressibility factor. Generally, sound propagation ismuch faster than viscous diffusion, and the compressibility factor is quite small; for example, ina dispersion of water and particles of radius a = 100 nm, the compressibility factor is estimatedto be ε ≈ 7 × 10−3 . In this case, the assumption of incompressibility is fully justified. Thus,
fluids are assumed to be incompressible in many theoretical studies of hydrodynamic interactions,and only viscous diffusion is considered as the temporal evolution mechanism for hydrodynamic
interactions.1, 2 However, the speed of sound in a liquid is around 103 m/s irrespective of the liquidconsidered, and this factor can produce a large compressibility factor for a dispersion of a highly
viscous fluid solvent, e.g., ε ≈ 0.7 for olive oil and ε ≈ 10 for corn syrup.

a) Electronic mail: tatsumi@cheme.kyoto-u.ac.jpb) Electronic mail: ryoichi@cheme.kyoto-u.ac.jp

1070-6631/2013/25(4)/046101/15/$30.00

25, 046101-1

C 2013 AIP Publishing LLC

046101-2

R. Tatsumi and R. Yamamoto

Phys. Fluids 25, 046101 (2013)

In recent years, the temporal evolution of hydrodynamic interactions between two particles has
been directly observed.3–6 In these experimental studies, particles were trapped by optical tweezers,
and the correlations between the positional fluctuations of particles were measured. The authors
reported the temporal evolution of hydrodynamic interactions in the viscous diffusion regime, whichcoincided with analytical predictions that assumed the fluid was incompressible. Evidence forhydrodynamic interactions caused by sound propagation was also observed3, 7, 8 but could not be
captured due to the extremely short time-scale of sound propagation.
In the present study, we investigate the propagation process of hydrodynamic interactions
using direct numerical simulation. Within this approach, the hydrodynamic interactions are directlycomputed by simultaneously solving for the motions of the fluid and the particles with appropriate
boundary conditions. We use the smoothed profile method (SPM),9, 10 which can be applied to acompressible fluid as well as an incompressible fluid.11 By considering the fluid compressibility,
sound propagation can be captured.
We consider a system of two particles in a fluid and investigate the correlated motion of theparticles. In particular, we estimate the velocity relaxation of one particle in response to the exertionof an impulsive force on the other particle. This velocity cross-relaxation function is equivalent tothe velocity cross-correlation function from the fluctuation-dissipation theorem. The velocity crossrelaxation function may be interpreted in terms of the temporal evolution of the flow field around theparticle. We first consider an incompressible fluid to identify the separate characteristics of soundpropagation and viscous diffusion in the hydrodynamic interactions. We then consider a compressible
fluid to investigate the effect of the order-of-magnitude relationship between the sound propagationand viscous diffusion time scales on the temporal evolution of the hydrodynamic interactions. In
addition, we examine the validity of analytical solutions within the Oseen approximation, which areoften compared with experimental results.
II. MODEL
A. Basic equations

We model a dispersion as a system in which rigid spherical particles are dispersed in a Newtonian
fluid. The motion of the particles is governed by Newton’s and Euler’s equations of motion, which
can be written for the i-th particle asMid
V i = F iH + F iC + F iE ,
dt
Ii ·

d
Ri = V i ,dtdi = N iH + N iE ,
dt

(2)
(3)

where Ri , V i , and i are the position, the translational velocity, and the rotational velocity of thei-th particle, respectively. The particle has a mass Mi and a moment of inertia I i . The fluid exerts
a hydrodynamic force F iH and a torque N iH on the particle, while a force F iC is exerted throughdirect interactions among the particles. A force F iE and a torque N iE are externally applied. The
hydrodynamic force and torque are evaluated by considering the momentum conservation between
the particle and the fluid.
The fluid dynamics are governed by the following hydrodynamic equations:∂ρ+ ∇ · (ρv) = 0,
∂t

(4)

∂ρv
+ ∇ · (ρvv) = ∇ · σ + f R ,∂t(5)

where ρ(r, t) and v(r, t) are the mass density and velocity fields, respectively, of the fluid. The stresstensor is given by


2
T
(6)
σ = − p I + η[∇v + (∇v) ] + ηv − η (∇ · v)I,3

046101-3

R. Tatsumi and R. Yamamoto

Phys. Fluids 25, 046101 (2013)

where p(r, t) is the pressure, η is the shear viscosity, and ηv is the bulk viscosity. A body forcef R (r, t) is added to satisfy particle rigidity. We also assume a barotropic fluid described by p = p(ρ)
with a constant speed of sound c such thatdp
= c2 .
dρ(7)Equations (4)–(7) are closed for the variables ρ, v, and p without considering energy conservation.
We use the SPM to perform direct numerical simulations. In this method, the boundaries between
particle and fluid are modeled by a continuous interface. For this purpose, a smoothed profile function
φ(r, t) ∈ [0, 1] is introduced to distinguish between the particle and fluid domains, i.e., φ = 1 in theparticle domain and φ = 0 in the fluid domain. These two domains are smoothly connected througha thin interfacial region of thickness ξ . The body force for the particle rigidity f R is expressed as
ρφ f p . The detailed mathematical expressions of φ and φ f p have been previously given.9
B. Linear formulation

The hydrodynamic equations can be linearized for low Reynolds number flows in a dispersion.
Correspondingly, the equations of motion of the particles are also linear. Most approximate theories
are based on this formulation.
Here, for an Np -particle system, we define the following column vectors, which are 6Np dimensional vectors, to describe the equations concisely:
⎛ E ⎞
⎞
⎛
⎛ H ⎞
F1
V1
F1
⎜
⎟
⎟
⎜
⎜ .. ⎟⎜ .. ⎟⎜ .. ⎟
⎜ . ⎟
⎜
⎟
⎜. ⎟
.
⎜
⎟
⎜
⎟
⎟
⎜
⎜ FH ⎟⎜⎟⎟⎜E⎜ Np ⎟F⎜⎜ V N p⎟
Np ⎟
⎜ H ⎟
⎜
⎟
⎜
U =⎜
⎟, E = ⎜ E ⎟⎟, H = ⎜⎟.
⎜ N1 ⎟⎜ N1 ⎟⎜ 1 ⎟⎜⎟⎜⎟⎟⎜⎜ .. ⎟⎜ .. ⎟⎜ .. ⎟
⎜ . ⎟
⎜
⎜. ⎟
. ⎟
⎝
⎠
⎝
⎠
⎠
⎝
H
N Np
N EN p
N p
Neglecting direct particle interactions, the equations of particle motion, i.e., Eqs. (2) and (3), aresummarized bydU(t) = H(t) + E(t),(8)
dtwhere the matrix M represents the direct sum of the mass tensor and the moment of inertia tensor.
For spherical particles of equal radii and densities, the matrix M reduces toM·M=

M I 3N p002
Ma 2 I 3N p
5

,

(9)

where M and a are the mass and the radius, respectively, of each particle. By performing a Fouriertransform on Eq. (8), the corresponding equation for the Fourier components with time factor e−iωt
is obtained as−iωM · Û(ω) = Ĥ(ω) + Ê(ω).

(10)

The following relation holds true for the linearized hydrodynamic equations:2
Û(ω) = −μ̂(ω) · Ĥ(ω),(11)

where μ̂ is the 6Np × 6Np mobility matrix. The mobility matrix depends on the configuration of allof the particles {Ri }, so we further assume, for simplicity, that the particles do not move significantly
over the time scale considered, i.e., the configuration of particles {Ri } is independent of time t.
Substituting Eq. (11) into Eq. (10) yields an explicit expression for the translational and rotational

046101-4

R. Tatsumi and R. Yamamoto

Phys. Fluids 25, 046101 (2013)

Particle 1

Particle 2

R

a

z
y
x

FIG. 1. Geometry of the present model system. Two spherical particles of equal radius a are situated in a fluid with acenter-to-center distance R = |R|. This system is axisymmetric about the R = R2 − R1 direction such that two directions
are of interest: parallel and perpendicular to the center-to-center axis of the particles. In this figure, the y and z directions aredegenerate perpendicular directions.

particle velocities:
Û(ω) = [I − iωμ̂(ω) · M]−1 · μ̂(ω) · Ê(ω).

(12)

Therefore, for an externally applied force and torque, the dynamics of particles in a fluid can be
described by the mobility matrix.
The components of the mobility matrix are given byμ̂(ω) =

μ̂tt (ω)μ̂tr (ω)μ̂rt (ω) μ̂rr (ω)

,

(13)
αβ

where μ̂αβ is the 3Np × 3Np matrix, which is composed of the mobility tensors μ̂i j for each pairingof particles i and j. The mobility tensors describe the mutual coupling of the translational androtational motions between particles, and the superscript represents the mode of the motion: “t” fortranslation and “r” for rotation. The mobility matrix is symmetric as a consequence of the Lorentz
reciprocal theorem such that the following relations are satisfied:12μ̂tt = (μ̂tt )T , μ̂rr = (μ̂rr )T , μ̂tr = (μ̂rt )T .

(14)

Now, we consider a system of two identical spherical particles in a fluid. The configuration of
the particles is described in Fig. 1. Numbers are assigned to the particles: 1 for the particle on theleft and 2 for the particle on the right. We define a vector R = R2 − R1 to describe the geometryof this system. The particle center-to-center distance is denoted by R = |R|. The unit vector along
the line of centers is described by R̂ = R/R. Due to axisymmetry about the R̂ axis, each mobility
tensor is described by at most two scalar functions:13, 14
tt

(15a)

μ̂irrj (R, ω) = μ̂i j (R, ω) R̂ R̂ + μ̂irr⊥
j (R, ω)(I − R̂ R̂),

rr

(15b)

μ̂itrj (R, ω) = μ̂itr⊥
j (R, ω) R̂ × I.

(15c)

μ̂ittj (R, ω) = μ̂i j (R, ω) R̂ R̂ + μ̂itt⊥
j (R, ω)(I − R̂ R̂),The superscripts  and ⊥ denote the directions of motion parallel and perpendicular to thesymmetry axis, respectively. In the direction parallel to the symmetry axis, the translational androtational motions are decoupled as indicated by Eq. (15c). Interchanging the particle numbers 1
and 2 corresponds to an inversion of the direction of R̂, which causes a sign inversion for thetranslation-rotation cross-mobility tensor μ̂itrj without changing the sign of the other cross-mobility

046101-5

R. Tatsumi and R. Yamamoto

Phys. Fluids 25, 046101 (2013)

tensors, μ̂ittj and μ̂irrj . The mobility tensors Eqs. (15) within the Oseen approximation are describedin the Appendix.
C. Temporal evolution of the flow field

We first consider a particle dispersion with an incompressible fluid solvent. Incompressibility
corresponds to an infinite speed of sound in the fluid such that the solenoidal condition is imposedon the fluid velocity field from Eqs. (4) and (7) as∇ · v = 0.

(16)

The velocity field can be decomposed into two contributions as
v = w − ∇ϕ.

(17)

The vector field w(r, t) is the juxtaposition of the velocity fields in the fluid and particle domains;
however, the fluid-particle impermeability boundary condition, which demands the continuity of thenormal velocity on the boundary, is not satisfied. An irrotational flow field described by the scalarpotential ϕ(r, t) is added so that the velocity field v will satisfy the boundary condition. Due to thesolenoidal condition on the velocity field v, the scalar potential ϕ obeys the Poisson equation∇ 2 ϕ = ∇ · w.

(18)

On the right-hand side, ∇ · w is zero except at the fluid-particle boundary at which singularities
exist. From Eq. (18), the scalar potential is given by15
ϕ(r, t) = −

1
4π

dr 

(r − r  )· w(r  , t).
|r − r  |3

(19)

If we consider an isolated single particle, the scalar potential and the corresponding velocity field
can be simply described due to axisymmetry:ϕ(r, t) =Q(t) · r
,
r3

(20)

Q(t) · (3r̂ r̂ − I)
,
(21)
r3
where r̂ = r/r is a unit vector in the radial direction. The velocity field of Eq. (21) represents adoublet flow, which corresponds to the electric field generated by an electric dipole. The vector Q(t)is parallel to the particle velocity and depends on time to the same degree as the particle motion
relative to the fluid. The analytical form of Q(t) has been derived previously,16 and the strength
Q(t) = | Q(t)| is time-independent for a neutrally buoyant particle. The doublet flow is generated
instantaneously to satisfy the solenoidal condition; thus, the doublet flow is interpreted to expand at
the infinite speed of sound.
The vector field w represents the shear flow due to viscous diffusion. When the Reynolds number
is sufficiently low, the hydrodynamic equations can be linearized as follows:−∇ϕ(r, t) =∂v
= −∇ p + η∇ 2 v + f R ,(22)∂twhere ρ0 is the fluid mass density at equilibrium.
Under the condition of incompressibility, the pressure gradient imposes fluid-particle impermeability on the time-derivative of the velocity field, and the pressure relates to the scalar potential as
ρ0

∂ϕ
.
∂t
Therefore, the temporal evolution of the vector field w is given byp = ρ0ρ0∂w= η(∇ 2 w − ∇∇ · w) + f R .
∂t

(23)

(24)

046101-6

R. Tatsumi and R. Yamamoto

Phys. Fluids 25, 046101 (2013)

Modifying the diffusion term on the right-hand side to be solenoidal results in the doublet flow being
gradually offset by the diffusion of the vector field w.
In short, the exertion of an impulsive force on a particle in an incompressible fluid instantaneously generates a doublet flow. The shear flow then spreads out diffusely, which cancels out the
doublet flow. This picture is validated in the numerical simulation results below.
For a compressible fluid, a second contribution is added to Eq. (17):v = w − ∇ψ − ∇ϕ,(25)

where the scalar potential ψ(r, t) corresponds to the volume compression flow. Here, we assume
that the dynamics of the fluid can be described by the linearized hydrodynamic equations:∂p= −ρ0 c2 ∇ · v,
∂t


∂v
1
ρ0
= η∇ 2 v +η + ηv ∇∇ · v − ∇ p + f R .
∂t
3

(26)
(27)

The time evolution equation of the total scalar potential ϕ  = ϕ + ψ is derived from the equations
above asηl1 ∂ 2ϕ
ηl
2 
2 ∂ϕ
2 ∂
= 1+
(∇ · w),(28)∇∇∇ ϕ − 2 2 +
c ∂t
ρ0 c2
∂t
ρ0 c2 ∂twhere ηl = (4/3)η + ηv denotes the longitudinal viscosity. This is a damped wave equation in which
the source is given on the right-hand side of the equation. Therefore, a potential flow propagates
at a finite speed of sound in a compressible fluid. In the long-time limit, Eq. (28) reduces to the
Poisson equation for an incompressible fluid as given by Eq. (18), and the potential flow reduces toan instantaneous doublet flow.
III. NUMERICAL RESULTS

Numerical simulations are performed for a three-dimensional box with periodic boundary
conditions. The space is divided into meshes of length , which is a unit length. The units of theother physical quantities are defined by combining  with η = 1 and ρ 0 = 1. The system size is Lx
× Ly × Lz = 256 × 256 × 256. The other parameters are set to a = 4, ξ = 2, ρ p = 1, ηv = 0, and h= 0.05 where ρ p is the particle mass density and h is the time increment of a single simulation step.
We consider the system of two identical spherical particles in a fluid whose geometry is describedin Fig. 1. Both particles have the same density and radius, which also means that they have the same
mass and moment of inertia. We investigate the time-dependence of the velocity for particle 1
following the exertion of an impulsive force at the center of particle 2, for various center-to-center
distances between the particles R∗ = R/a. This cross-relaxation function is a manifestation of thetemporal evolution of the hydrodynamic interactions between particles 1 and 2. The impulsive force
is assumed to be sufficiently small such that the Reynolds number of the flow is sufficiently low to
exclude the effects of turbulence. We set the impulsive force to produce an initial particle Reynolds
number of Rep = 10−3 . Additionally, in this condition, the displacement of particles is negligible
because the particle displacement before stopping, scaled by the particle radius a, is comparable to(ρ p /ρ 0 )Rep . Therefore, direct interactions between particles, which include overlap repulsion forces,
are not considered because particle collisions do not occur.
The cross-relaxation tensor is introduced as the normalized change in velocity of particle 1:
V 1 (R1 , t; R2 ) =P2· γ 12 (R, t),M⊥
γ 12 (R, t) = γ12 (R, t) R̂ R̂ + γ12(R, t)(I − R̂ R̂),

(29)
(30)

where P 2 is the impulsive force exerted on particle 2 at t = 0. Due to the axisymmetry, thecross-relaxation tensor is characterized by only two directions of motion, which are parallel andperpendicular to the center-to-center vector R, and the motions in both directions are decoupled.

046101-7

R. Tatsumi and R. Yamamoto

Phys. Fluids 25, 046101 (2013)

From the fluctuation-dissipation theorem, the cross-relaxation tensor is equivalent to the velocity
cross-correlation function of the fluctuating system:
γ 12 (R, t) =MV 1 (0, 0)V 2 (R, t) ,
kB T

(31)

where kB is the Boltzmann constant and T is the thermodynamic temperature. We examine both⊥by adjusting the direction of thethe parallel correlation γ12 and the perpendicular correlation γ12
impulsive force P 2 . As predicted by the formulation in Sec. II B, translation-rotation coupling in
the particle motion is observed for motion perpendicular to the axis. However, we focus entirely ontranslation-translation coupling because the influence of sound propagation on the rotational motion
is expected to be small. We also compare the simulated results with approximate solutions wherein
the analytical solution for the self-mobilities of isolated single particles and the Oseen approximation
for the cross-mobilities are applied (see the Appendix).
A. Incompressible fluid

First, we investigate the temporal evolution of hydrodynamic interactions in an incompressible
fluid. The dynamics of the incompressible fluid is governed by the hydrodynamic equations composed of Eqs. (5), (6), and (16). The center-to-center distance between particles takes the values of R∗= 3, 4, and 7. The simulation results for the velocity cross-relaxation functions are shown in Fig. 2.
⊥(t)While the parallel correlation γ12 (t) is positive at all times, the perpendicular correlation γ12
is initially negative and subsequently becomes positive. Comparable results have been reported inexperimental studies3–5 where the difference between the parallel and perpendicular correlations
was attributed to the flow field around the particle.4, 7 We discuss the time dependence of the velocity
cross-relaxation functions in more detail below.
The temporal evolution of the flow field was described in Sec. II C. Here, the velocity fields
generated by the impulsive force exerted on a single particle are shown in Fig. 3. These simulation
results are obtained for a single particle system with size Lx × Ly × Lz = 128 × 128 × 128. At earlytimes such as t/τ ν = 2.5 × 10−2 , the doublet flow is dominant. The doublet flow is characterized byloop streamlines flared in the direction perpendicular to the particle motion, and backflow is observed
as described by Eq. (21). As discussed in Sec. II C, the doublet flow appears instantaneously andcan be interpreted as the propagation of an infinite-speed sound wave. Conversely, shear flow due toviscous diffusion, whose strength is given by the vorticity ∇ × v, is only observed in the vicinity ofthe particle. However, at a time t/τ ν = 1.2, the shear flow has diffused over a large range following(a) 10-1(b) 10-1
R* = 3R* = 4
R* = 7tailR* = 3R* = 4
R* = 7
tail

-2

10

-2

//
12 (t)⊥
12 (t)

10

-3

10

-3

10

-4

10

-4

10

-5

-2

10

-1

0

10

10

t/

10

1

10

-2

10

-1

0

10

10

10

1

t/

FIG. 2. Velocity cross-relaxation function in an incompressible fluid in (a) parallel and (b) perpendicular directions relative tothe symmetry axis of the particles. The distances between the two particles are R∗ = 3, 4, and 7. Simulation results are shown
as bold solid lines (positive values) and dashed-dotted lines (negative values). Solutions within the Oseen approximation are
represented by thin solid lines (positive values) and dashed-dotted lines (negative values). The broken line represents thelong-time tail with an algebraic power law decay Bt−3/2 as given by Eq. (32).

046101-8

R. Tatsumi and R. Yamamoto

Phys. Fluids 25, 046101 (2013)

FIG. 3. Temporal evolution of the velocity field around an isolated single particle due to an impulsive force exerted on theparticle in an incompressible fluid. Cross-sections parallel to the impulsive force direction including the particle center are
shown. The impulsive force is exerted at time t = 0, and the simulation results are given at (a) t/τ ν = 2.5 × 10−2 and(b) t/τ ν = 1.2. The direction of the impulsive force is to the right in the pictures, and the particle is represented by a black
circle. The vorticity of the velocity field ∇ × v is described by a color scale, which goes from negative (darker) to positive(lighter) vorticity. The vorticity is normalized by τ ν /Rep .

Eq. (24); therefore, there is a spreading fluid region flowing in the same direction as the particle
motion and the loop streamlines get away from the particle.
The velocity cross-relaxation function is related to the temporal evolution of the flow fieldaround particle 2. The doublet flow propagated by sound waves produces an instantaneous velocity
correlation between the particles. The doublet flow direction shown in Fig. 3 produces a positive
parallel correlation and a negative perpendicular correlation. Almost no change in the cross-relaxation
functions occurs in the early stages, which reflects the time-independence of the strength of the
doublet flow due to the neutrally buoyant particle.16
The reduction of the parallel correlation and the sign inversion of the perpendicular correlation
occur at about the same time, which depends on the inter-particle distance R∗ . Because shear flowby viscous diffusion can cause such changes in the cross-relaxation functions, we estimate the timescale on which the shear flow generated by the particle 2 arrives at particle 1 by the viscous diffusion
time scale over the length L: τν∗ = L 2 /ν. In this case, the characteristic length is the distance between
the particle surfaces, L = R − 2a; thus, for R∗ = 3, 4, and 7, the viscous diffusion time scales areτν∗ /τν = (R ∗ − 2)2 = 1, 4, and 25, respectively. The time at which the parallel correlation begins todecrease and the perpendicular correlation has a sign inversion roughly corresponds to τν∗ in eachcase.
In an incompressible fluid, hydrodynamic interactions are instantly propagated at the infinite
speed of sound and are subsequently transmitted by viscous diffusion on the time scale τν∗ . The
temporal evolution of the flow fields around two particles separated by R∗ = 7 is shown in Fig. 4.
At early times in the flow, t/τ ν = 6.25 × 10−1 , the shear flow has not arrived at particle 1, and its
dynamics are governed by the doublet flow. Later, at t/τ ν = 12.5, the shear flow has diffused overtime to arrive at particle 1, and the motion of particle 1 is governed by the shear flow. Here, thestresslet flow field is observed in the vicinity of particle 1. This flow is generated to maintain particle
rigidity against deformation in a shear flow. In the perpendicular correlation, rotational motion ofparticle 1 is also observed.
In the final stages of the flow, the momentum of particle 2 has completely diffused away, andthe particles and the fluid move collectively. These dynamics are manifested in a long-time tail witha power law decay t−3/2 of the cross-relaxation functions shown in Fig. 2, which are given by171 ρp
,⊥γ12 (t) = √9 π ρ0τν 3/2
as t → ∞.
t(32)

046101-9

R. Tatsumi and R. Yamamoto

Phys. Fluids 25, 046101 (2013)-1(a) t/τν = 6.25 x10 (parallel)(c) t/τν = 12.5 (parallel)

x10-2
8

(b) t/τν = 6.25 x10-1 (perpendicular)

x10-2
8

4

4

0

0

-4

-4

-8

-8

x10
4

-4

(d) t/τν = 12.5 (perpendicular)

x10
4

2

2

0

0

-2

-2

-4

-4

-4

FIG. 4. Temporal evolution of the velocity field around the particles. Cross-sections including the center-to-center axis
between particles are shown. The center-to-center distance between two particles is R∗ = 7. The impulsive force is exerted
at time t = 0 on the particle on the right, and the simulation results are presented at (a) and (b) t/τ ν = 6.25 × 10−1 and(c) and (d) t/τ ν = 12.5. The direction of the impulsive force is (a) and (c) to the right or (b) and (d) upwards. The particles
are represented by black circles. The vorticity of the velocity field ∇ × v is described by a color scale, which ranges fromnegative (darker) to positive (lighter) vorticity. The vorticity is normalized by τ ν /Rep .

Discrepancies between the simulation results and the Oseen approximation are observed especially at t/τ ν  3 as shown in Fig. 2. The approximate cross-relaxation functions clearly changewith time in the early stages, and the time at which the effect of the shear flow becomes apparent islater than in the simulation results. Particles are assumed to be points in the Oseen approximation;
therefore, it is implied that the particle separation is much larger than the particle radius as R∗= R/a1, and the time scale of observation is much longer than that of viscous diffusion over
1. Neglecting the particle size results in a
the length of the particle radius as t/τ ν = tν/a2larger characteristic length L = R and a longer diffusion time scale τν∗ /τν = R ∗2 . Consequently, theOseen approximation can only be accurately applied at long times and large particle separations. The
simulation results shown in Fig. 2 confirm that the validity of the Oseen approximation is enhanced
by increasing inter-particle distance R∗ . In previous experimental studies, the measurements wereconducted over the range for which the Oseen approximation is valid, i.e., particle separations R∗  5
and measurement frequencies corresponding to t/τ ν  30. Consequently, the experimental results
showed good agreement with the Oseen approximation.3, 4

B. Compressible fluid

Next, we consider a compressible fluid in which sound propagates at a finite speed. The velocity
cross-relaxation functions for inter-particle distances R∗ = 3 and 7 are shown in Figs. 5 and 6,
respectively. Due to the periodic boundary conditions imposed on the simulation box, a sound pulse
generated by the particle motion can continue to affect the simulation results after the sound pulse
arrives at the edge of the simulation box. The oscillatory structure of the cross-relaxation functions,especially at ε = 0.1, is one of the striking artifacts of the periodic boundary conditions.

046101-10

R. Tatsumi and R. Yamamoto

Phys. Fluids 25, 046101 (2013)(a) 10-1(b) 10-1

R* = 3incompressible
= 0.6
= 1.5

R* = 3incompressible
= 0.1

-2

-2

//
12 (t)

10//
12 (t)

10

-3

-3

10

10

-4

10

-4

-2

10

-1

0

10

10

10

10

1

-2

10

-1

t/

10

1

(d) 10-1

R* = 3incompressible
= 0.6
= 1.5

R* = 3incompressible
= 0.1

-2

-2

10
⊥
12 (t)

10⊥
12 (t)

10t/(c) 10-1

-3

10

-4

-3

10

-4

10

10

-5

10

0

10

-5

-2

10

-1

0

10

10

t/

10

1

10

-2

10

-1

0

10

10

10

1

t/

FIG. 5. Velocity cross-relaxation function in a compressible fluid in the (a) and (b) parallel and (c) and (d) perpendicular
directions to the symmetry axis. The center-to-center distance between two particles is R∗ = 3. The compressibility factor
takes the values of ε = 0.1, 0.6, and 1.5. Simulation results are shown as bold solid lines (positive values) and dashed-dottedlines (negative values). Solutions using the Oseen approximation are represented by thin solid lines (positive values) anddashed-dotted lines (negative values).

The velocity cross-relaxation functions are zero in the early stages and then suddenly change tonon-zero values. This time lag increases with the inter-particle separation and is likely to be due tothe sound propagation from particle 2 to particle 1. As with the diffusion of shear flow, we estimate
the time scale for sound propagation over a characteristic length L = R − 2a as τc∗ = L/c. For
R∗ = 3 and 7, the sound propagation time scales are τc∗ /τν = (R ∗ − 2)ε = ε and 5ε, respectively.
The results shown in Figs. 5 and 6 confirm the correspondence between the time scale of soundpropagation and the peak of the cross-relaxation functions. The peak in the cross-relaxation function
broadens as sound propagation is dissipated by the longitudinal viscosity according to Eq. (28). The
velocity fields around particles separated by R∗ = 7 at t/τ ν = 6.25 × 10−1 are shown in Fig. 7.
Expanding doublet flows with strengths given by the divergence of the velocity field ∇ · v areobserved. For a small compressibility factor of ε = 0.1, the doublet flow has just arrived at particle
1. This situation corresponds to a peak in the cross-relaxation functions shown in Figs. 6(a) and6(c). Conversely, for a larger compressibility factor of ε = 0.6, the doublet flow has not yet arrived
at particle 1, and the cross-relaxation functions are, correspondingly, zero.
When the fluid compressibility is as small as ε = 0.1, the velocity cross-relaxation functions
superimpose onto those for the incompressible fluid after a time τc∗ . This behavior only describes
hydrodynamic interactions by sound propagation before viscous diffusion effects come into play.
Conversely, for a large compressibility of ε = 1.5, the behavior of the cross-relaxation functions is
considerably different except in the hydrodynamic long-time tail; in particular, there is no negative
perpendicular correlation. For a fluid with a medium compressibility such as ε = 0.6, a negative

046101-11

R. Tatsumi and R. Yamamoto

Phys. Fluids 25, 046101 (2013)(a) 10-1(b) 10-1

R* = 7

incompressible
= 0.6
= 1.5

R* = 7

incompressible
= 0.1

-2

-2

//
12 (t)

10//
12 (t)

10

-3

-3

10

10

-4

10

-4

-2

10

-1

0

10

10

10

10

1

-2

10

-1

10

1

(d) 10-1

R* = 7

incompressible
= 0.6
= 1.5

R* = 7

incompressible
= 0.1

-2

-2

10
⊥
12 (t)

10⊥
12 (t)

10t/(c) 10-1

-3

10

-4

-3

10

-4

10

10

-5

10

0

10

t/

-5

-2

10

-1

0

10

10

10

1

10

-2

10

t/

-1

0

10

10

10

1

t/

FIG. 6. Velocity cross-relaxation function in a compressible fluid in the (a) and (b) parallel and (c) and (d) perpendicular
directions to the symmetry axis. The center-to-center distance between two particles is R∗ = 7. The compressibility factor
takes values of ε = 0.1, 0.6, and 1.5. Simulation results are shown as bold solid lines (positive values) and dashed-dottedlines (negative values). Solutions using the Oseen approximation are represented by thin solid lines (positive values) anddashed-dotted lines (negative values).

perpendicular correlation is observed only for the inter-particle separation of R∗ = 7. Because
sound propagation and viscous diffusion generate negative and positive perpendicular correlations,
respectively, a balance of the two time scales τν∗ and τc∗ is expected to characterize the generalbehavior of the cross-relaxation function. Therefore, we define the interactive compressibility factor
using a characteristic length L = R − 2a as
ε∗ =

τc∗
ε
.
= ∗
∗
τν
R −2

(33)

As the inter-particle distance increases, the two time scales separate to reduce the interactive compressibility factor. In the present simulations, the interactive compressibility factors for R∗ = 3 and 7
are ε∗ = ε and 0.2ε, respectively. When the interactive compressibility factor is small, hydrodynamic
interactions propagating at the speed of sound arrive first at the other particle in about a time τc∗ .
Following these interactions are those propagating by viscous diffusion, which arrive at about atime τν∗ . The temporal evolution is qualitatively the same as in an incompressible fluid except for
the time lag due to sound propagation at a finite speed. This situation applies to the results where anegative perpendicular correlation is observed in the early stage. Conversely, the order of arrival ofsound propagation and viscous diffusion should be reversed if the interactive compressibility factor
is sufficiently large. In this situation, the perpendicular correlation is positive from the start due toviscous diffusion, and the effect of sound propagation on the hydrodynamic interactions produces a
local minimum or even temporary negative values in the cross-correlation function.

046101-12

R. Tatsumi and R. Yamamoto(a) ε = 0.1 (parallel)

Phys. Fluids 25, 046101 (2013)

x10-2
2

(b) ε = 0.1 (perpendicular)

x10-2
2

1

1

0

0

-1

-1

-2

-2

(c) ε = 0.6 (parallel)(d) ε = 0.6 (perpendicular)
0.4

0.4

0.2

0.2

0

0

-0.2

-0.2

-0.4

-0.4

FIG. 7. Velocity field around particles in a compressible fluid. The simulation results were obtained at t/τ ν = 6.25 × 10−1for compressibility factors with values (a) and (b) ε = 0.1 and (c) and (d) ε = 0.6. The divergence of the velocity field ∇ · v
is described by a color scale, which goes from negative (darker) to positive (lighter) divergence. The value of divergence isnormalized by τ ν /Rep .

In Fig. 8, the velocity cross-relaxation functions for various center-to-center distances R∗ andcompressibility factors ε are classified according to the sign of the perpendicular correlation in thefirst stage. We can find that the perpendicular correlation is positive from the start when ε∗ ≥ 0.175
is satisfied. In this condition, the hydrodynamic interactions are propagated by viscous diffusionahead of sound propagation.

1.5

1

0.5

0

2

4

6

8

R*
FIG. 8. Classification of calculated cross-relaxation functions for various center-to-center distances of particles and compressibility factors. The cross-relaxation functions are classified according to the sign of the perpendicular correlation in thefirst stage: negative and positive are represented by downward and upward triangles, respectively. The broken line represents
Eq. (33) with ε ∗ = 0.175.

046101-13

R. Tatsumi and R. Yamamoto

Phys. Fluids 25, 046101 (2013)

In a compressible fluid, the time condition for the validity of the Oseen approximation is given byt max (τ ν , τ c ) because there is a finite sound propagation time scale as well as a viscous diffusion
time scale. As with the diffusion of shear flow, sound propagation from particle 2 to particle 1 isdelayed within the Oseen approximation as shown in Figs. 5 and 6. This delay occurs because thesound propagation time scale is increased within this approximation so that τc∗ /τν = R ∗ ε. Thus,
the discrepancies between the Oseen approximation results with the simulation results can be fairly
well-resolved by increasing the particle separation.

IV. CONCLUSION

In the present study, we investigated the temporal evolution of hydrodynamic interactions for asystem of two particles using SPM to perform direct numerical simulations. Fluid compressibility
was considered in examining temporal evolution by sound propagation. Hydrodynamic interactions
were estimated by the velocity cross-relaxation functions, which are equivalent to the velocity
cross-correlation functions in a fluctuating system.
In an incompressible fluid, hydrodynamic interactions were observed to propagate instantaneously at the infinite speed of sound followed by subsequent temporal evolution by viscous
diffusion. Theoretical analysis showed that a doublet flow and a shear flow are generated by soundpropagation and viscous diffusion, respectively. Because the cross-relaxation function reflects thecharacteristics of each flow field, sound propagation and viscous diffusion are associated with negative and positive perpendicular correlations between the particle velocities, respectively. Therefore,
the time at which the behavior of the cross-relaxation function changes is related to the time scaleof viscous diffusion over the particle separation τν∗ , which is the only time scale characterizing the
temporal evolution of hydrodynamic interactions in an incompressible fluid.
In a compressible fluid, sound propagates between particles in a finite time that scales as τc∗ .
The effect of the order-of-magnitude relationship between the two time scales τν∗ and τc∗ on the
cross-relaxation function was observed. An interactive compressibility factor ε∗ was defined as theratio of the two time scales given by Eq. (33). In our simulation results, a reversal of the order ofarrival of sound propagation and viscous diffusion at the other particle was observed at ε∗ ≥ 0.175;in this case, the perpendicular correlation is positive from the start and hydrodynamic interactions
were largely governed by viscous diffusion.
The temporal evolution of hydrodynamic interactions does not qualitatively change from thatfor an incompressible fluid as long as the interactive compressibility factor is small. Only when theinteractive compressibility factor is sufficiently large can differences in the temporal evolution ofhydrodynamic interactions be expected. Such a situation could be realized in a highly viscous fluid,
and experimental validation of the results presented here is desirable.

ACKNOWLEDGMENTS

This work was supported by KAKENHI 23244087 and the JSPS Core-to-Core Program “International research network for non-equilibrium dynamics of soft matter.”

APPENDIX: MOBILITY MATRIX IN THE OSEEN APPROXIMATION

Hydrodynamic interactions among particles in low Reynolds number flows are described by themobility matrix discussed in Sec. II B. However, there is no known closed-form solution even for a
two-particle system. In this section, we introduce an approximate analytical form for the mobility
tensors.
For low Reynolds number flows, the hydrodynamic equations can be linearized as given by
Eqs. (26) and (27), and the corresponding equations for the Fourier components with the time factor
e−iωt are obtained as−iω p̂ω + ρ0 c2 ∇ · v̂ ω = 0,

(A1)

046101-14

R. Tatsumi and R. Yamamoto

Phys. Fluids 25, 046101 (2013)
−iωv̂ ω = η∇ v̂ ω +2

1R
η + ηv ∇∇ · v̂ ω − ∇ p̂ω + f̂ ω .
3(A2)

The fluid velocity field generated by the body force f R is expressed asv̂ ω (r) =R

dr  Ĝ(r − r  , ω) · f̂ ω (r  ),where the Green’s function is given by18, 19
 −αr

e
1
eiμr − e−αr
I + α −2 ∇∇,Ĝ(r, ω) =
4π η
r
r

(A3)

(A4)

with
α = (−iωρ0 /η)1/2 ,andμ = ω/c̃,

(A5)



1/2
iω 4
c̃ = c 1 −
η
+
η
.
v
ρ0 c2 3

(A6)

For an isolated single spherical particle in a fluid with a body force constraint f R to satisfy thecondition of particle rigidity, the self-mobilities can be analytically derived18, 20, 21 as follows:μ̂tt11 (ω) = μt0

92x 2 (1 − i y) − (1 + x)y 2 − x 2 y 2,
2x 2 (1 + x)(9 − 9i y − 2y 2 ) + x 2 (1 − i y)3(1 + x)μ̂rr11 (ω) = μr0,
3 + 3x + x 2

(A7)
(A8)

where x = αa, y = μa, μt0 = (6π ηa)−1 , and μr0 = (8π ηa 3 )−1 . When the speed of sound is assumed
to be infinite, i.e., μ = 0, the solution for an incompressible fluid is obtained. The translation-rotation
coupling does not appear in the self-mobility, i.e., μ̂rt11 = 0.
To calculate the cross-mobilities, we introduce the Oseen approximation in which the particles
are regarded as points. In the Oseen approximation, a particle exerts a Stokeslet (point force);
therefore, according to Eq. (A3), the translation-translation cross-mobility is the Green’s function
itself:μ̂tt12 (R, ω) = Ĝ(R, ω).

(A9)The components of the parallel and perpendicular directions are explicitly given by
3a[(2 − 2iμR − μ2 R 2 )eiμR − (2 + 2α R)e−α R ],(A10a)
2α 2 R 3
3at[(−1 + iμR)eiμR + (1 + α R + α 2 R 2 )e−α R ].
(A10b)μ̂tt⊥
12 (R, ω) = μ02α 2 R 3
The Stokeslet also generates a rotational motion with an angular velocity of (1/2)∇ × v̂ ω ; therefore,
the translation-rotation cross-mobility is given bytt

μ̂12 (R, ω) = μt01
μ̂rt12 (R, ω) = − ∇ × Ĝ(R, ω),2which has a perpendicular-only component

(A11)

a3
(1 + α R)e−α R .
(A12)R2There are no effects from the fluid compressibility. This result can also be derived by considering thevelocity field generated by a rotlet (point torque).22 The rotation-rotation cross-mobility is obtained
from the angular velocity generated by a rotlet:r
μ̂rt⊥
12 (R, ω) = μ0μ̂rr12 (R, ω) =1
∇ × ∇ × Ĝ(R, ω).
4

(A13)

046101-15

R. Tatsumi and R. Yamamoto

Phys. Fluids 25, 046101 (2013)

Then, each component is given byrr

a3(1 + α R)e−α R ,R3(A14a)

a3(1 + α R + α 2 R 2 )e−α R .
2R 3(A14b)μ̂12 (R, ω) = μr0r
μ̂rr⊥
12 (R, ω) = −μ0

Using the analytical solutions for the self-mobilities of an isolated single particle and theOseen approximation for the cross-mobilities, the approximate velocity cross-relaxation functions
are obtained by Eq. (12).
1 W. van Saarloos and P. Mazur, “Many-sphere hydrodynamic interactions: II. Mobilities at finite frequencies,” Physica A

120, 77–102 (1983).
2 H. J. H. Clercx and P. P. J. M. Schram, “Retarded hydrodynamic interactions in suspensions,” Physica A 174, 325–354(1991).
3 S. Henderson, S. Mitchell, and P. Bartlett, “Propagation of hydrodynamic interactions in colloidal suspensions,” Phys. Rev.

Lett. 88, 088302 (2002).
4 M. Atakhorrami, G. H. Koenderink, C. F. Schmidt, and F. C. MacKintosh, “Short-time inertial response of viscoelastic

fluids: Observation of vortex propagation,” Phys. Rev. Lett. 95, 208302 (2005).
5 S. Martin, M. Reichert, H. Stark, and T. Gisler, “Direct observation of hydrodynamic rotation-translation coupling between

two colloidal spheres,” Phys. Rev. Lett. 97, 248301 (2006).
6 Y. von Hansen, A. Mehlich, B. Pelz, M. Rief, and R. R. Netz, “Auto- and cross-power spectral analysis of dual trap opticaltweezer experiments using Bayesian inference,” Rev. Sci. Instrum. 83, 095116 (2012).
7 P. Español, M. A. Rubio, and I. Zúñiga, “Scaling of the time-dependent self-diffusion coefficient and the propagation ofhydrodynamic interactions,” Phys. Rev. E 51, 803–806 (1995).
8 P. Español, “On the propagation of hydrodynamic interactions,” Physica A 214, 185–206 (1995).
9 Y. Nakayama and R. Yamamoto, “Simulation method to resolve hydrodynamic interactions in colloidal dispersions,” Phys.

Rev. E 71, 036707 (2005).
10 Y. Nakayama, K. Kim, and R. Yamamoto, “Simulating (electro) hydrodynamic effects in colloidal dispersions: Smoothed

profile method,” Eur. Phys. J. E 26, 361–368 (2008).
11 R. Tatsumi and R. Yamamoto, “Direct numerical simulation of dispersed particles in a compressible fluid,” Phys. Rev. E

85, 066704 (2012).
12 J. Happel and H. Brenner, Low Reynolds Number Hydrodynamics (Kluwer Academic, Boston, 1983).
13 D. J. Jeffrey and Y. Onishi, “Calculation of the resistance and mobility functions for two unequal rigid spheres in low-

Reynolds-number flow,” J. Fluid Mech. 139, 261–290 (1984).
14 M. Reichert and H. Stark, “Hydrodynamic coupling of two rotating spheres trapped in harmonic potentials,” Phys. Rev. E

69, 031407 (2004).
15 J. D. Jackson, Classical Electrodynamics, 3rd ed. (Wiley, New York, 1999).
16 B. U. Felderhof, “Transient flow caused by a sudden impulse or twist applied to a sphere immersed in a viscous incom-pressible fluid,” Phys. Fluids 19, 073102 (2007).
17 B. Cichocki and B. U. Felderhof, “Long-time tails in the solid-body motion of a sphere immersed in a suspension,” Phys.

Rev. E 62, 5383–5388 (2000).
18 D. Bedeaux and P. Mazur, “A generalization of Faxén’s theorem to nonsteady motion of a sphere through a compressible

fluid in arbitrary flow,” Physica 78, 505–515 (1974).
19 B. U. Felderhof, “Transient dipolar sound wave in a viscous compressible fluid,” Physica A 389, 5602–5610 (2010).
20 B. U. Felderhof, “Effect of fluid compressibility on the flow caused by a sudden impulse applied to a sphere immersed ina viscous fluid,” Phys. Fluids 19, 126101 (2007).
21 B. U. Felderhof, “Force density induced on a sphere in linear hydrodynamics: II. Moving sphere, mixed boundary

conditions,” Physica A 84, 569–576 (1976).
22 S. Kim and S. J. Karrila, Microhydrodynamics (Butterworth-Heinemann, Boston, 1991).