Multiple collisions in turbulent flows
Michel Voßkuhle1 , Emmanuel Lévêque1,2 , Michael Wilkinson3 and Alain Pumir1
1
arXiv:1310.1521v1 [physics.flu-dyn] 5 Oct 2013
2
Laboratoire de Physique, ENS de Lyon and CNRS,
UMR5672, 46, allée d’Italie, F-69007 Lyon, France
Laboratoire de Mécanique des Fluides et d’Acoustique, Ecole Centrale de Lyon and CNRS,
UMR5509, 36 avenue Guy de Collonge F-69134, Ecully, France
3
Department of Mathematics and Statistics, The Open University,
Walton Hall, Milton Keynes, MK6 7AA, England
In turbulent suspensions, collision rates determine how rapidly particles coalesce or react with
each other. To determine the collision rate, many numerical studies rely on the ‘Ghost Collision
Approximation’ (GCA), which simply records how often pairs of point particles come within a
threshold distance. In many applications, the suspended particles stick (or in the case of liquid
droplets, coalesce) upon collision, and it is the frequency of first contact which is of interest. If
a pair of ‘ghost’ particles undergoes multiple collisions, the GCA may overestimate the true collision rate. Here, using fully resolved Direct Numerical Simulations of turbulent flows at moderate
Reynolds number (Rλ = 130), we investigate the prevalence and properties of multiple collisions.
We demonstrate that the GCA leads to a systematic overestimate of the collision rate, which is of
the order of ∼ 15 % when the particle inertia is small, and slowly decreases when inertia increases.
We investigate the probability P (Nc ) for a given pair of ghost particles colliding Nc times. We find
P (Nc ) = βαNc for Nc > 1, where α and β are coefficients which depend upon the particle inertia.
This result is used to explain the discrepancy between the GCA and the true collision rates. We
also investigate the statistics of the times that ghost particles remain in contact. We show that the
probability density function of the contact time is different for the first collision. The difference is
explained by the effect of caustics in the phase space of the suspended particles.
PACS numbers: 47.27.-i, 05.40.-a, 45.50.Tn, 92.60.Mt
I.
INTRODUCTION
Collisions between particles transported by a turbulent flow play a crucial role in several important phenomena. Saffman and Turner [1] suggested that turbulence
in clouds can lead to a very significant enhancement of
the rate of collision between small droplets. This mechanism has been proposed to provide an explanation for
the very fast rate of coalescence reported in warm cumulus clouds [2–4]. Also, collisions between dust grains in a
turbulent circumstellar accretion disc play an important
role in some theories for planet formation [5, 6].
There is a substantial literature on the collision rate of
particles in turbulent flows. Most of the studies devoted
to collisions in turbulent suspensions explicitly deal with
‘geometric collisions’, that is, merely detect when the centers of two particles are separated by a distance d less
than the sum of their radii, a1 + a2 . The time-dependent
separation d(t) may cross a1 + a2 repeatedly, and every
crossing from above is considered as a new collision. In
many applications, however, the colliding particles are
assumed to stick, coalesce, or react on first contact. In
these cases the multiple collisions are spurious, and the
physically relevant collision rate should only count the
frequency with which d(t) decreases below a1 + a2 for
the first time. The ‘ghost collision approximation’ (GCA)
consists in ignoring this aspect, and in following all particles in the flow even after they underwent collisions.
While this approximation is appealing from a numerical
point of view, it has been noticed to lead to questionable
estimates of the collision rate [7, 8], because it includes
spurious multiple collisions [9].
The first aim of this paper is to fully characterize
the shortcomings of the GCA approximation which arise
from the existence of multiple collisions. To this end, we
write the collision rate Γ as a sum of two terms:
ΓGCA = Γ1 + Γm
(1)
where ΓGCA is the collision rate obtained with the GCA.
The term Γ1 corresponds to the collision rate for particles undergoing a first contact, and Γm is the rate for
multiple collisions between a pair of particles that has already collided once. Thus Γm represents the unphysical
contribution to the GCA collision rate due to multiple
encounters between a given pair of particles. One way to
evaluate the true collision rate Γ1 is to count the rate of
multiple collisions, and subtract it from ΓGCA .
What specifically happens to two particles after they
have come into contact depends on the precise physical
problem [10, 11]. We are interested here in problems
where particles have a strong probability of reacting after their first contact. For the sake of the present work,
we assume here that Γ1 , defined by considering only the
first collisions occurring between a pair of particles, is
an appropriate definition of the collision rate. We have
also considered two other processes, which lead to alternative definitions. The first process consists in assuming
explicitly that the two particles do not participate any
longer to the reactions as soon as they have come into
contact, as if they had annihilated. Technically, in order to deal with a steady-state system, we introduce two
2
new particles to compensate for the loss. The implementation of this algorithm requires some care in the case
of inertial particles, as the velocity depends on history.
Another process consists in removing and replacing only
one of the two colliding particles. We demonstrate here
that the collision rates describing these two processes are
equal to Γ1 , in the limit of very dilute suspensions, which
justifies its relevance.
The other objective of our study is to characterize the
statistics of the multiple collisions. We investigate the
statistics of the number of contacts, Nc . We find that
there is a very simple distribution of the number of times
a pair of particles collides: the probability for observing
Nc collisions after one initial collision is found to be well
approximated by
P (Nc |Nc ≥ 1) = β αNc
(2)
where the coefficients α and β depend upon parameters
as discussed below. We use these results to explain the
spurious collision rate, Γm .
It has been observed that tracer particles in a turbulent flow can remain in proximity for a long time [12–14].
This phenomenon may be related to the multiple collisions which we are investigating here, and could be of interest to model classes of reaction occurring with a finite
probability when particles are close enough [11]. With
this motivation, we consider the statistics of the time
that a pair of ghost particles is in contact and relate this
to the distribution of relative velocities of collisions. The
probability distribution function (PDF) of the contact
time exhibits a striking structure. For the first collision,
it follows a power-law at intermediate values, and an exponential decay at longer contact times. The PDF of the
contact time for the first collision are markedly different from those of subsequent collisions (which lack the
power-law behavior and which appear to be independent
of the number Nc of collisions, for Nc ≥ 2). We explain
this discrepancy by appealing to a model for the collision
process introduced in [2] and [4], according to which the
collision rate Γ can be expressed as the sum of two terms,
resulting from two different types of physical processes:
Γ = Γadv + Γcaust .
(3)
The term Γadv represents the rate of collisions of particles
which are advected into contact by shearing motion due
to turbulence. This advective process allows for the possibility of multiple collisions, because the local velocity
gradient can fluctuate as a function of time, so that particles may be brought back into contact after separating
for a while. The term Γcaust results from the possibility
that particles with significant inertia can move relative
to the fluid. If caustics (fold lines) form in the phasespace of the suspended particles, particles occupying the
same position may have different velocities. This phenomenon can also be understood as resulting from particles being centrifuged out of vortices, and has also been
referred to as the ‘sling’ effect [2, 15]. Experiments have
recently demonstrated the existence of the ‘sling’ effect
in a well-controlled laboratory flow [16]. Particles colliding at higher relative velocity are unlikely to undergo
multiple collisions. The higher-velocity collisions generated by caustics contribute to the statistics of the first
collision, but not to the multiple collisions. The use of
the decomposition (3) is supported by work demonstrating that it is valid for a random flow model [17, 18]. A
companion to this paper [19] shows that (3) is a good
model for DNS studies of turbulence, and discusses the
relative importance of the two terms as a function of the
parameters of the turbulent suspension.
We will show that the power-law PDF of contact times
is entirely due to the first collision between particles, and
not to the subsequent ones. The subsequent collisions
appear to be self-similar, in the sense that the distributions of times when particles are within a given distance
remain unchanged, no matter how many times the particles have come close together. We also consider possible
explanations for this observation.
The results presented here very significantly extend a
previous study, using a simpler model of turbulent flow,
namely the Kinematic Simulation approach [20]. While
the results concerning the relative errors made by using
this simplified flow are qualitatively similar to the results
presented in this work, we find that Kinematic Simulations lead to a very significant underestimation of the
collision rates, by almost an order of magnitude.
Our work emphasizes models where the suspended particles coalesce upon contact, so that multiple collisions
are a source of error in the collision rate. These models
are just a limiting case of a larger class of models, which
may also be of physical interest. As an example, the case
of fully elastic collisions was recently considered by [10],
who showed that elastic collisions lead to a finite probability of particles making an infinite number of collisions.
Our investigations of the contact times in the ghost collision approximation involve evaluating the distribution
of the time that one particle in a turbulent flow spends
within a sphere surrounding another particle. We remark
that Jørgensen et al [21] have investigated this distribution experimentally in a different context, where the radius of the sphere is much larger than the Kolmogorov
length of the flow.
This article is organized as follows. In Section II, we
present the numerical schemes used to simulate the carrying turbulent flow and account for the dynamics and
collisions of the suspended particles. We describe the
algorithms used to study the physical problem of particles coagulating during their first contact, and to establish that the reaction rate reduces, in the very limit, to
Γ1 . The nature of the correction due to a finite particle density are discussed in depth in an Appendix. The
quantitative estimates of the errors made by using the
GCA are discussed in Section III, together with the data
justifying the empirical law (2) for a pair of particles to
undergo multiple collisions. Section IV is devoted to the
statistics of the multiple collisions. Section V discusses
3
the explanation for some of the observations in terms of
statistics of the relative velocities upon collision. Finally,
we summarize our results in Section VI.
II.
A.
NUMERICAL METHODS
Direct Numerical Simulation of Navier-Stokes
Turbulence
The work rests on simulating the (incompressible)
Navier–Stokes equations:
∂t u(x, t) + (u(x, t) · ∇)u(x, t) = −∇p(x, t) + ν∇2 u(x, t)
+f (x, t)
(4)
∇ · u(x, t) = 0
(5)
where u(x, t) denotes the Eulerian velocity field, ν is the
viscosity, and f (x, t) is a forcing term; the mass density
is arbitrarily set to unity. These equations are solved in
a cubic box of size 2π with periodic boundary conditions
in the three directions by a pseudo-spectral method. The
pressure p(x, t) is eliminated by taking the divergence of
(4) and by solving the resulting Poisson equation in the
spectral domain [22]. The forcing term acts on Fourier
modes of low wavenumbers, |k| ≤ Kf . It is adjusted in
such a way that the injection rate of energy, ε, remains
constant [23]:
fk = ε P
uk
|k|≤Kf
|uk |2
if
|k| ≤ Kf .
(6)
The simulations discussed in this study have been done
with the following parameters, in code units: ε = 10−3 ,
and ν = 4 × 10−4 . The forcing operates at wavenumbers |k| ≤ Kf = 1.5. With these values, the Kolmogorov scale, ηK = (ν 3 /ε)1/4 , is comparable to the effective spatial resolution ∆x = 2π/256: ∆x/ηK ≈ 1.5 or
kmax η ≈ 2.1, which fulfills standard requirements for the
direct numerical simulation (DNS) of Navier-Stokes turbulence and the integration of particle trajectories [24].
Let us notice that, in practice, the number of grid points
is 384 in each direction according to the two-thirds rule
[25] (used to avoid aliasing errors). However, padded
high-wavenumber modes are not excited and, therefore,
do not contribute to improving the spatial resolution of
the solution. With these values of the parameters, the
value of the Reynolds number based on the Taylor microscale is Rλ ≈ 130 in the statistically stationary longtime limit state.
After spectral truncation, (4) reduces to a set of ordinary differential equations (in time) for the Fourier
modes, which have been integrated using the secondorder Adams-Bashforth scheme. The time step δt has
been chosen so that the Courant number Co = urms ·
kmax δt <
∼ 0.1. With this choice, the relaxation time τp
of the particle dynamics remains of the order of 102 · δt,
which ensures that particle trajectories can be safely integrated from the time-evolving Eulerian velocity field.
B.
Dynamics of particles
We simulated in our flow the motion of small particles,
all of which are assumed to have the same size and mass,
whose motion obeys the following set of equations:
dv
u(x, t) − v
=
.
dt
τp
dx
= v,
dt
(7)
This set of equations is a simplified version of the original set derived by [26, 27], and is appropriate for small
spherical particles with radius a which is much smaller
than the Kolmogorov scale ηK , and whose density is much
larger than the fluid density: ρp /ρf ≫ 1. In order to isolate the role of turbulence, we have explicitly neglected
gravity in (7), despite the fact that it plays an important
role in cloud microphysics, sandstorms and other terrestrial phenomena. Gravitational effects on collision rates
are unimportant in applications to planet formation. The
relaxation time in (7) is determined by the Stokes drag:
τp =
2 ρp a 2
.
9 ρf ν
(8)
This relaxation time is made dimensionless by using the
Kolmogorov time scale, τK = (ν/ε)1/2 , and the Stokes
number
τp
(9)
St =
τK
is a dimensionless measure of the importance of inertial effects in determining the trajectories of the particles. In order to explore parameter ranges which are relevant to cloud microphysics, we used a ratio of densities
ρp /ρf = 103 throughout. Lengths and times are given in
dimensionless units and can be readily scaled to realistic
situations.
The determination of the particle velocity v requires,
according to (7), the evaluation of the fluid velocity u
at the location of the particle. This is done by resorting
to tri-cubic interpolation. The particle trajectories have
been integrated by using the second-order Verlet velocity
algorithm [28].
As we are interested in determining the collision rates
in turbulent flows, we simulated a large number of particles. The number of particles Np used to monitor the
collision rates was chosen in such a way that the volume fraction Φ = 4Np πa3 /(3L3 ), where L is the size
of the system (L = 2π), is either Φ = 4.5 × 10−6 or
Φ = 4.5 × 10−5 . In both cases, collisions involving three
or more particles can be neglected.
In a flow having reached a statistically steady state,
we inserted at a time T = 0 a total number Np′ of particles, initially distributed uniformly in the flow. We then
integrated the equation of motions (4),(7), for a time of
the order of 10 eddy-turnover times, TL , defined by:
Z
L
3π
TL = p
, with L =
k −1 E(k) dk . (10)
2hu2 i
hu2 /3i
4
Past this time, we integrated the equations of motion
for a time ttot larger than 15TL (except for St = 0.2).
All the particle trajectories were saved, with a sampling
time of ∆t ≈ 0.055τK , and processed afterward. Table I
summarizes our runs.
For the case of St = 0, we integrated the trajectories
of Lagrangian tracer particles. Those are point particles
with no extent, but to determine the collision rate it is
necessary to assume they have a finite size. We chose the
radius to be the same as for particles with St = 0.1.
C.
Detection of the collisions
In a monodisperse solution of volume V , containing Np
particles, the average number of collisions Nc per unit of
time and per unit volume, is proportional to the square
of the average number of particles per unit volume, n ≡
Np /V :
Nc = Γc
n2
.
2
(11)
Eq. (11) defines the collision kernel Γc , which depends
both on the fluid motion, and on the physical properties
of the particles.
The particle trajectories determined numerically were
stored and post-processed separately, to determine the
number of collisions and other statistics. Detecting collisions requires checking the mutual distance between Np
particles, which typically requires of the order of Np2 operations. In simulations such as ours involving a large
number of particles, this can lead to a prohibitively large
computational time. We used an algorithm based on a
cell-linked list to speed-up the processing time, as done
in [29].
Each data set was processed in two different ways:
1. Multiple collisions of ghost particles
First we determined the collision rate using the
ghost collision approximation. In this case only
Np of the total Np′ simulated particles were used.
The collision kernel ΓGCA was simply determined
by counting the total number of detected collisions
Nc . In addition, we also determined the rates at
which pairs of particles come into contact for the
Nc th time, described by the collision kernel ΓNc .
To this end, at each collision between two particles, we examine the trajectories leading up to the
collision and determine the number of previous collisions between the same pair. If the pair has undergone Nc − 1 previous collisions with each other, the
collision event contributes to the kernel ΓNc . Because these collision kernels describe an exhaustive
and mutually exclusive decomposition of ΓGCA , we
have
ΓGCA =
∞
X
Nc =1
ΓN c
(12)
The collision kernel of multiple collisions Γm is
simply defined by summing the kernels ΓNc , for
Nc ≥ 2:
Γm =
∞
X
ΓN c .
(13)
Nc =2
2. Collision detection with particle replacement.
We have assumed that, in the case of particles coalescing or reacting on their first contact, Γ1 is the
correct measure of the collision rate. To establish
this result, we compared Γ1 with a collision rate algorithm where particles are removed from the flow
after collision. We determined the collision rate
again with Np particles, but we systematically replaced either one or two of the colliding particles after each collision. The substitutions are carried out
by simply picking one of the (Np′ −Np ) non-colliding
particles, making sure that the newly introduced
particle is not colliding with any other particle at
the moment of its insertion. In so doing, the collision rate is determined at a fixed density. A similar
method has been used in Refs. [7, 30]. The rates determined by this procedure are denoted ΓRe1 and
ΓRe2 , depending upon whether one or two particles
are replaced. We stress that the simplest method,
consisting in generating new particles at random
positions, cannot work in the case of finite inertia
(St 6= 0), as the velocity depends on history. Simulating Np′ > Np particles allows us to deal only with
particles which are already in equilibrium with the
flow. We refer in the following to these algorithms
as the one- or two-particle ‘substitution schemes’.
It is not immediately clear that ΓRe1 = ΓRe2 = Γ1 .
In fact, our numerical results reveal measurable deviations from these identities at the highest particle
volume fraction studied here. However, as we explain in detail in the Appendix, these deviations
decrease linearly when reducing the particle density, thus establishing that Γ1 = ΓRei in the limit
of very dilute suspensions.
In all cases the number of collisions Nc (∆τ ) is measured in consecutive intervals of length ∆τ ∼ TL . For
each interval the collision kernel can be determined as
hΓi∆τ =
2V Nc (∆τ )
,
∆τ Np2
(14)
which is simply the collision rate divided by n2 /2. These
‘instantaneous’ collision kernels vary in time. Determining the level of fluctuations of the number of collisions
recorded over a limited time interval led to an estimate
of the uncertainty of the collision rate. The resulting uncertainty in the numerial value of ΓGCA is less than 2 %,
except for St = 0.2 and St = 0.5, where it is no larger
than 4 %.
5
TABLE I. (Color online) Summary of parameters from our different DNS runs. We tabulate the Stokes number St, the volume
fraction occupied by the particles Φ, the total integration time ttot used to determine the collision rate, expressed in terms of the
large eddy turnover time TL , and the total number of collisions NGCA detected, when using the ghost collision approximation.
St
Φ × 106
ttot /TL
NGCA /104
0.0
4.5
15.5
0.6
0.10
4.5
15.5
1.3
0.20
45
4.2
29
0.30
4.5
31.4
3.3
0.51
45
15.9
200
0.76
45
10.4
180
1.01
45
47.6
400
1.27
45
13.0
81
III. QUANTIFYING THE ERROR INDUCED
BY THE GHOST COLLISION APPROXIMATION
A.
Error of the Ghost-Collision Approximation
In this section, we discuss the error induced by using the ghost collision approximation. The left panel of
Fig. 1 shows the collision rate computed by using the
ghost collision approximation (curve with square symbols), and the collision rate Γ1 (curve with circle symbols). Our estimates of the collision rates compare quantitatively very well with previous numerical work [31, 32].
The right panel shows the relative difference Γm /Γ1 =
(ΓGCA − Γ1 )/Γ1 .
Our results show that among all the collisions recorded,
as many as 15% of them involve pairs of particles that
have already come into contact before. We note that our
results show that the collision rate determined by using the GCA is correctly approximated by the Saffman–
Turner formula when St → 0, consistent with previous
results [7]. However, in the limit of small inertia, the
property that pairs of particles collide more than one
time is very significant. This probability decreases when
the Stokes number increases. This result can be qualitatively understood by using the known fact that when the
Stokes number increases, the particle trajectories differ
more from the fluid trajectories, which allows collisions
between particles with an increasingly large velocity difference. Colliding particles are therefore expected to separate quicker when the Stokes number is large, thus making multiple collisions less likely.
In physical situations where particle pairs react upon
first contact, the collision rate Γ1 is expected to provide
the correct estimate of the reaction rate. To this end, we
have measured the rates ΓRe1 and ΓRe2 , defined by taking
one or two of the colliding particles out of the system
after collisions. The values of ΓRe1 and ΓRe2 are found
to be actually smaller than Γ1 . In fact, as we explain in
the Appendix, the difference between ΓRei (i = 1, 2) and
Γ1 is due to close encounters between three particles, an
effect whose relative importance becomes weaker when
the total density of particles goes to zero. The analysis
presented in the Appendix establishes that in the limit of
a very dilute system, the proper collision rate is indeed
Γ1 .
1.52
45
52.4
250
2.03
45
52.1
150
B.
2.53
45
52.5
99
3.04
45
52.3
72
4.05
45
53.1
42
5.07
45
42.1
22
Statistics of multiple collision
In view of the importance of multiple collisions between a pair of particles, we characterize here the statistical properties of the number of collisions a given pair of
particles undergoes before separating. Fig. 2 shows the
probability distribution function for a particle pair to collide an Nc th time after at least one initial collision. The
clear result from Fig. 2 is that for Nc ≥ 2, the probability
distribution is very well approximated by an exponential
law of the form of equation (2). This remarkably simple functional form leads to the interpretation that once
a pair of particles has undergone more than one collision, it has a probability α to collide once more before
it separates, the quantity α being independent of Nc .
This suggests a Markovian process of multiple collisions,
amenable to a simple modeling.
The exponential law, equation (2), can be used to sum
the series in equation (12) and hence to determine the
error of the ΓGCA estimate. We have ΓNc = Γ1 βαNc for
Nc > 1, so that
Γm = Γ1 β
∞
X
Nc =2
αNc = Γ1
α2 β
.
1−α
(15)
The numerical values obtained from expression (15), with
the values of α and β extracted from the probability of
Nc (see Fig. 2(a)) agrees quantitatively very well with
the value of Γm determined directly (see Fig. 2(b)). (the
two values differ by less than 10 %).
IV.
STATISTICS OF THE CONTACT TIME
BETWEEN PARTICLE PAIRS
Our numerical observation that a given pair of particles may collide many times in a turbulent flow can be
related to some surprising properties of particle trajectories, which have been partly documented before [12–14]
mostly in the case of tracers (St = 0). In this section, we
fully characterize several properties concerning the time
particle pairs spend together, which are relevant to the
subject of the present article.
6
(a)
(b)
100
0.2
0.15
60
Γm /Γ1
ΓτK /(2a)3
80
40
0
0.05
ΓGCA
Γ1
20
0
1
2
3
4
0.1
0
5
0
1
2
St
3
4
5
St
FIG. 1. (Color online) Comparison between the collision kernels ΓGCA and Γ1 . The collision kernel ΓGCA in panel (a) was
obtained according to the ghost collision approximation, taking into account all collisions. The collision kernel Γ1 is restricted
to only first collisions of a same pair. The difference, Γm = ΓGCA − Γ1 , which corresponds to the relative error introduced by
the GCA, is shown in the left panel. This error starts at a value close to ∼ 15% when St → 0 and decreases for higher Stokes
numbers.
(a)
(b)
10−1
10−2
10−3
10−4
10−5
10−6
10−7
α
P (Nc |Nc ≥ 1)
100
2
Γ1 βα2 /(1 − α)
Γm (measured)
5
0 1 2 3 4 5
St
0.5
1.0
2.0
4.0
0
6
0.2
0.15
0.1
0.05
ΓτK /(2a)3
101
4
3
2
1
4
6
8
10
Nc
0
0
1
2
3
4
5
St
FIG. 2. (Color online) Panel (a): The probability that a pair of particles collides an Nc th time conditioned on the fact, that it
collides at least once. The probability of observing Nc collisions goes as βα(St)Nc . The results are shown at different Stokes
numbers. The values of α(St) are shown in the inset of the figure as a function of St. Panel (b): The kernel for multiple
collisions Γm determined from Eq. (15) with the fitting parameters deduced from panel (a) square symbols), and measured
directly in our simulations (circles). The quantitative agreement between both results confirms the consistency of our reasoning
and the quality of our fits.
A.
Particle trajectories can stay close together for
a long time
Fig. 3 shows the distance between two pairs of particles
over a long time. Panel (a) shows the distance in units of
the size of the periodic box during one entire simulation.
As expected, the particles are far from each other for
most of the time. However, we notice that the distance
can reach a very small value at some particular moments.
Whereas one of the pairs of particles (dashed line) is close
to each other only for a short while, the other pair spends
a much longer time close to each other. Panel (b) blows-
up the result in the range of time where the distance is
small. (Here the distance is shown in terms of the collision radius 2a.) The pair whose distance is plotted with
a dashed line reaches the value of 2a, i.e. collides, but
bounces apart immediately. In comparison, the continuous curve shows that the two particles stay together for a
time larger than the large eddy-turnover time. Over this
period of time, the distance fluctuates close to the value
of (2a), which causes multiple collisions between the two
particles. The value of the Stokes number of the particles
shown in Fig. 3 is St = 1.0; the phenomenon shown here
is qualitatively similar at different Stokes numbers.
7
(a)
(b)
0.8
10
8
d(t)/2a
d(t)/2π
0.6
0.4
0.2
0
6
4
2
0
5
0
10 15 20 25 30 35 40 45
t/TL
10
12
14
16
18
20
t/TL
FIG. 3. (Color online) Distance between two different pairs of particles as a function of time. Panel (a) shows the distance
over several eddy turnover times, and shows that particles can be close to each other for a short time (dashed curve), or for
a much longer time (full line), corresponding to several large eddy turnover times TL . Panel (b) magnifies the graph shown
in panel (a) in the region where the distance between the particles is very small. The full line reveals that multiple collisions
occur. The dotted line in panel (a) is the average distance for particle pairs homogeneously distributed in a periodic cube. In
panel (b), the dotted line indicates the collision radius 2a.
ts,1 − te,1
Distribution of contact times
The phenomenon reported here is very reminiscent of
the observation discussed in several earlier works [12–14],
namely that particles can stay close for a very long time.
While this property has been documented mostly in the
case of tracers, our observations suggest that inertial particles can also remain very close for a long time. As it has
been noticed in a slightly different context [11], it is of
general interest to characterize the statistical properties
of the time that particles spend together.
In this subsection, we determine quantitatively the statistical properties of the time particle trajectories remain
close to one another. We illustrate in the two following
subsections the statistics at a fixed value St = 1.5 of the
Stokes number.
Fig. 4 introduces our notation. We consider a pair
of particles which become closer than a threshold distance dc at an instant of time. We denote by te,1 the
first instant for which the distance between the particles
ts,2 − te,2
···
···
2a
te,1
B.
te,3 − te,2
te,2 − te,1
distance
We observe that the timescale for these multiple collision processes is comparable to the turnover time of the
largest eddies in our simulation TL . In the following investigations, we have expressed our results in terms of
a dimensionless time t/TL . However the ratio between
TL and the Kolmogorov time scale, τK in our simulations is not very large: TL /τK ≈ 10.5. At present it is
unclear whether t/TL is the most natural dimensionless
timescale. We remark that a plausible alternative is t λ1 ,
where λ1 is the leading Lyapunov exponent. Because
λ1 ≈ 0.15/τK [33], the variables t/TL and t λ1 are quite
similar in magnitude.
ts,1 te,2
ts,2 te,3 ts,3
time
FIG. 4. (Color online) Illustration of the definition of the
different times: time of first encounter te,1 , time of first separation ts,1 , time of second encounter te,2 and of second separation ts,2 .
becomes less than dc , (that is, the instant of their first
collision) and then, ts,1 the first instant ts,1 > te,1 when
they separate. If particles approach to within a distance
dc again at some later time, we denote by te,2 the first instant for which the distance becomes again smaller than
dc , and ts,2 the time at which the particles separate again.
This notation can be easily generalized when the trajectories become more than two times closer than dc . The
time particles spend together during their first encounter
is denoted ∆T1 ≡ (ts,1 − te,1 ), which again, can be easily
generalized to the time particles spend together during
their nth encounter.
Fig. 5 shows the distribution of the time particles
8
(a)
(b)
102
TL P (ts,1 − te,1 )
100
10−1
10−2
10−3
10−4
10−5
10−6
10−7
0
1
2
3 4 5 6 7
(ts,1 − te,1 )/TL
8
TL P (ts,1 − te,1 )
dc = 4a
dc = 2a
dc = a
dc = a/2
101
9 10
102
101
100
10−1
10−2
10−3
10−4
10−5
10−6
10−7
10−3
dc
dc
dc
dc
= 4a
= 2a
=a
= a/2
10−2
10−1
100
(ts,1 − te,1 )/TL
101
FIG. 5. (Color online) Lin-log (panel a) and log-log (panel b) representations of the distribution of time spent together by pairs
of particles during their first encounter. Several values of the critical distance dc are shown; the similarity between the curves
shows that the dependence on dc is very weak (note that dc ≪ η for all the values of dc shown on the figure). The Stokes
number used here is St = 1.5.
spend together during their first encounter ∆T1 . The
data corresponds to particles with a Stokes number St =
1.5. Panel (a) shows the PDF of ∆T1 in lin-log scale. In
this figure, ∆T1 is expressed in units of the large eddy
turnover time, TL , the correlation time of the flow. Fig. 5
shows an exponential decay of the PDF, with a characteristic time of the order of the large eddy turnover time.
Thus, particles can spend together a time which is comparable to the large eddy turnover time (at least for the
value of Rλ used in our simulations). Fig. 5(a) also indicates that the PDF has a very sharp maximum around
∆T1 ≈ 0. In fact, Panel (b) shows the PDF in log-log
units, and suggests a power law distribution of ∆T1 . The
exponent Rof the power law measured here is of the order
−1.5. As ǫ x−1.5 dx diverges when ǫ → 0, the PDF necessarily saturates for time separations <
∼ ∆t, where ∆t is
the time at which we saved trajectories. Fig. 5 shows the
PDF of ∆T1 , determined with several values of the critical distance dc , very small compared to the Kolmogorov
scale η (we have a/η ≈ 1/12). The PDFs are remarkably independent of dc , at least provided the ratio dc /η
is small, where η is the Kolmogorov length.
Fig. 6 shows the probability distribution function for
∆Tn , with n = 2, 3, 4 (n > 1). The PDFs still exhibit for
very long values of ∆Tn exponentially decaying tails, see
panel (a), with the same decay rate as obtained for ∆T1 .
However, the short time behavior of the PDF does not
exhibit the very large peak seen in Fig. 5. In fact, the
probability distribution does not exhibit any power law
at short values of ∆Tn , as seen in panel (b).
The difference between the statistical properties of
∆T1 and ∆Tn for n ≥ 2 is therefore restricted to the
short time behavior. In physical terms, the probability
that the two particles do not spend much time together
is much larger during the first encounter than during the
following one. Particles can spend a short time when
they are impacting each other with a large velocity difference. One may surmise that in such a case, particles
will separate very fast, and not get into close contact afterward. In other words, the events leading to several
contacts are unlikely to have initially a small value of
∆T1 . To actually check this, Fig. 6 also shows the PDF
of ∆T1 , conditioned on the fact that the two trajectories
will come into contact more than one time, see the curve
with the square symbols. As expected, conditioning the
probability of ∆T1 on the fact that there are more than
one encounter between the trajectories significantly reduces the probability for the particles of separating very
fast. In this spirit, Fig. 7 shows the PDF of ∆T1 conditioned on having several successive encounters between
the trajectories (N > 1), or simply one (N = 1), together with the total PDF of ∆T1 . As the probability
of having multiple encounters between the two particles
remains relatively small, the PDF of ∆T1 is extremely
close to the PDF of ∆T1 conditioned on having N = 1.
Aside from the difference between the PDF of ∆T1 and
∆Tn at short times, it appears that the process leading
to subsequent encounters between particle trajectories is
largely self-similar i.e. does not depend much on the index n > 1. This effect is strengthened by studying the
difference between the time it takes for two trajectories
to come into contact again. To this end, Fig. 8 shows the
probability distribution function of ∆Tne = (te,n+1 −te,n ).
The PDF is found to be remarkably similar (independent
of n). As it was the case for ∆Tn , the PDF has an exponential tail at large values of ∆Tne . The distribution
however peaks at a finite value of ∆Tne ≈ 0.6TL .
9
(a)
(b)
101
101
TL P (ts,1 − te,1 |Nc > 1)
TL P (ts,2 − te,2 )
TL P (ts,3 − te,3 )
TL P (ts,4 − te,4 )
100
10−1
10−2
100
10−1
10−2
10−3
10−3
10−4
10−4
10−5
10−5
10−3
0
1
2
3
4
5 6 7
(ts,i − te,i )/TL
8
9 10
TL P (ts,1 − te,1 |Nc > 1)
TL P (ts,2 − te,2 )
TL P (ts,3 − te,3 )
TL P (ts,4 − te,4 )
10−2
10−1
100
(ts,i − te,i )/TL
101
FIG. 6. (Color online) Lin-log (panel a) and log-log (panel b) representations of the distribution of time spent together by
pairs of particles during second (circle), third (upward pointing triangle), fourth (diamond) encounters. These distributions
are remarkably similar. For comparison, the distribution of time that particles spend together during their first encounter,
conditioned on the fact that they will meet again, is shown (square symbols). The cuspy distribution, observed for ts,1 − te,1 ,
is not seen for these distributions. The Stokes number used here is St = 1.5.
(a)
(b)
102
TL P (ts,1 − te,1 |Nc = 1)
TL P (ts,1 − te,1 |Nc > 1)
TL P (ts,1 − te,1 )
101
100
10−1
10−2
10−3
10−4
10−5
10−6
10−7
0
1
2
3
4
5 6 7
(ts,1 − te,1 )/TL
8
9 10
102
101
100
10−1
10−2
10−3
10−4
10−5
10−6
10−7
10−3
TL P (ts,1 − te,1 |Nc = 1)
TL P (ts,1 − te,1 |Nc > 1)
TL P (ts,1 − te,1 )
10−2
10−1
100
(ts,1 − te,1 )/TL
101
FIG. 7. (Color online) The distribution of the time spent together by two particles during their first encounter, p(ts,1 − te,1 )
(upward pointing triangles), and the same distribution conditioned on the fact that trajectories will separate after their first
encounter and not meet again (Nc = 1; square symbols). Finally the same PDF is shown again, this time conditioned on
the fact that the particles will meet again (Nc ≥ 2; circles). The statistics is dominated by pairs that collide only once. The
large probability, hence the power-law distribution, that two particles spend a short time together comes from trajectories with
Nc = 1. Both conditional probabilities show the exponential tail for large times. The Stokes number used here is St = 1.5.
C.
Stokes number dependence of the contact-time
statistics
The general picture, shown in subsections IV A and
IV B for particles with a Stokes number St = 1.5,
has been found to be qualitatively unchanged when the
Stokes number is varied. However, the details differ quantitatively. This can be seen by representing the distribution of the time difference ∆T1 by an asymptotic fit of
the form:
P (∆T1 ) ≈ N (∆T1 /TL )−ξ exp(−κ∆T1 /TL )
(16)
The coefficient N in (16) is simply adjusted to enforce
that the PDF is properly normalized. The coefficients ξ
and κ are determined by fitting the PDFs. The quality
of the fit is very good, as shown in Fig. 9 (a), at least
for 0.3 <
∼ St. The fitting parameters are found to depend
very significantly on St, see Fig. 9 (b). In the very small
Stokes number limit, the exponent ξ diminishes, suggesting that the distribution P (∆T1 ) is becoming closer to
a purely exponential distribution. In the opposite limit,
the coefficient κ decreases, whereas the power ξ seems to
increase.
10
(a)
(b)
101
i=1
i=2
i=3
100
10−1
TL P (te,i+1 − te,i )
TL P (te,i+1 − te,i )
101
10−2
10−3
10−4
10−5
0
2
4
6
8
(te,i+1 − te,i )/TL
10
i=1
i=2
i=3
100
10−1
10−2
10−3
10−4
10−5
10−6
10−2
12
10−1
100
101
(te,i+1 − te,i )/TL
102
FIG. 8. (Color online) Statistics for successive encounters between two particle trajectories. The PDF of te,i+1 − te,i in lin-log
scales, te,i being defined in Fig. 4. Panel (a) shows the data in lin-log scaling, while panel (b) shows the same data in log-log
scaling. The PDFs do not depend on i, suggesting a self similar process. The Stokes number used here is St = 1.5.
(b)
(a)
TL P (∆T1 /TL )
101
100
10−1
10−2
10−3
10−4
10−5
10−6
10−7
10−8
10−9
Fit coefficients: P (∆T1 ) ∼ e−κ∆T1 /TL
St = 0.5
St = 1.0
St = 2.0
St = 4.0
2.5
∆T1
TL
−ξ
2
1.5
κ
ξ
1
0.5
0
2
4
6
8
∆T1 /TL
10
12
14
0
0
1
2
3
4
5
St
FIG. 9. (Color online) Fit of the PDF of ∆T by (16) at different values of the Stokes number. Panel (a) shows the quality
of the fit, for four values of the Stokes number, indicated in the legend. The quality of the fit degrades at small values of St,
suggesting that a more complicated functional form may be needed. Panel (b) shows the variation of the exponent ξ and of
the coefficient κ determined by fitting the PDF, as a function of St.
V.
VELOCITY DIFFERENCE BETWEEN
COLLIDING PARTICLES
The time that two ghost particles spend in ‘contact’
(that is the length of time over which the separation of
their centers satisfies d ≤ 2a) is determined by two factors: the impact parameter, and the relative velocity of
the collision. The data in section IV show two striking
aspects of the contact-time statistics. Firstly, the statistics of the first contact are different from those of all
the others. Secondly, there is evidence for a power-law
regime in the distribution of the first contact time. In
this section we discuss how both of these observations
can be explained in terms of properties of the distribution of collision velocities.
A.
Multiple collisions happen at slow relative
velocities
The probability of the radial relative velocity, defined
as wr = δv · δr/|δr|, is shown in Fig. 10 for colliding
particles and at St = 1. We remind that the PDF for
colliding particles is different from the one for all particles in contact [34]. In panel (a) this PDF is shown
for two different situations. In the first one, all colliding particles are taken into account, as is the case in
the GCA. In the other one, only pairs that collide for
the first time are taken into account. In both cases the
bulk of the PDF is located at small values of |wr |, and
exhibits an exponential tail at very large collision velocities (see inset). A close comparison of the two probabilities shows that the contribution of small relative veloc-
11
ities corresponding to first collisions is smaller than for
all particle pairs detected with the GCA. This suggests
that the error in the GCA stems mainly from collisions
with small relative velocities. The right panel of Fig. 10
extends this observation to different Stokes numbers by
comparing the mean radial relative collision velocities obtained with both schemes. For Stokes numbers <
∼ 2 the
collisions corresponding to particles in contact for the
first time show on average slightly larger radial relative
velocities. However, the averaged relative velocities for
particle pairs colliding more than one time are significantly smaller than both, the velocities of first collisions
and the velocities in the GCA. This can be seen in the
inset of the right panel. The observations summarized
in Fig. 10 therefore demonstrate that, multiple collisions
occur at small relative velocity.
It has been argued that the collision rate between particles in a turbulent suspension can be resolved into two
components, as represented by equation (3) [2, 4]. In
the decomposition (3), the term Γadv represents collisions
due to shearing motion, which occur with relative velocities of order a/τK , whereas the term Γcaust represents
collisions between particles which are moving relative to
the flow, and on different branches of a phase-space manifold, separated by a caustic. The relative velocity of the
collisions which contribute to Γcaust is much higher, of
order (η/τK )f (St), where f (St) is an increasing function,
which is of order unity at St = 1. The ideas underlying this decomposition also explain why the statistics of
the first contact time are different from those of all the
subsequent contacts. According to this picture, multiple
collisions are almost exclusively due to the advective collision mechanism, and are very unlikely for caustic-induced
collisions because the high relative velocity rapidly moves
the particles out of proximity. In the advective process
multiple collisions arise because of temporal fluctuations
of the shear rate tensor [9]. This model explains why the
first collision has different contact time statistics, and
why the multiple collisions have small relative velocities.
B.
Power-law distribution of contact times
The evidence for a power-law in the distribution of
the contact time for the first collision is one of the conspicuous results of Sec. IV. In order to understand the
origin of such power-laws, we consider the following simplified model. We computed directly the distribution of
the time ∆T that two particles in a gas of particles, with
a Maxwellian distribution of velocity spend within a distance 2a from each other. The root mean square of one
velocity component of these particles is taken to be σ.
This model corresponds to the very large St limit of inertial particles in a turbulent flow. In that case, the gas
of particles reduces to particles each moving with its own
velocity, distributed according to the Maxwell distribution [35]. A simple calculation leads to the following PDF
of the time that two particles spend together:
σ 1
1
P (∆T ) = 2
1 − exp(−ζ 2 ) 1 + ζ 2 + ζ 4
(17)
aζ
2
2a
(18)
with
ζ=
σ ∆T
In the limit of long times, ∆T → ∞, ζ → 0, the PDF in
Eq. (17) reduces to:
P (∆T ) ∝
1
∆T 5
(19)
This result was checked directly by generating a gas of
Maxwellian particles and determining the time particles
spend together (Fig. 11). The distribution of ∆T is found
in excellent agreement with (17). We furthermore note,
that Equation (17) displays another power law for short
contact times. But this behavior is only apparent for
∆T < 2a/σ ≪ τK . In our DNS of the turbulent transport
of particles we do not resolve these time scales.
Let us consider the implications of the model leading
to equation (19) for the distribution of first contact times.
The important point is to observe that a power-law in the
distribution of small relative velocities leads to a powerlaw distribution of long contact times. The data in figure 5 show a power-law distribution in the contact times
at short times. We should, however, remember that the
first contacts are a combination of caustic-mediated collisions (with high relative velocity) and advective collisions
(with low relative velocity). We propose that the powerlaw distribution of the contact time results from the lowvelocity tail of the caustic-mediated collisions having a
PDF which can be approximated by a power-law at small
velocities.
VI.
CONCLUSION
We have studied the collision rate in turbulent suspensions, with the aim of evaluating the systematic errors
made while using the GCA for particles which aggregate upon collision. To this end, we have compared the
results using the GCA, and a more realistic algorithm,
which consists in replacing one of the particles that underwent collision by a particle from a ‘reservoir’. The
error we find by using the GCA is as large as ∼ 15 % at
very small Stokes number, and tends to decrease when
St increases.
We investigated the statistics of the multiple collisions
which are the source of the spurious collisions in the
GCA. We found that, after the first collision, multiple
collisions involving the same pair of particles occur with
a probability that decays exponentially with the number
of collisions, Nc . We also studied the time that particle trajectories spend close to one another. We find the
PDF of the time spent by two particles within a distance
smaller than a critical value dc exhibits an exponential
tail at very long times. At shorter times the contact time
PDF of the first collision obeys a power law.
12
(a)
3.5
uK Pc,i (|wr |)
3
2.5
100
10−2
10−4
10−6
2
1.5
1
0
0.5
0
h|wr |ic,1 /h|wr |ic,GCA
i = GCA
i=1
0
0.2
0.4
0.6
|wr |/uK
10
20
0.8
1
1.14
1.12
1.1
1.08
1.06
1.04
1.02
1
0.98
h|wr |ic,i /uK
(b)
0
1
5
4
3
2
1
0
i=1
i=m
0
1
2
2
3
3
St
4
4
5
5
St
(2a/σ) P (ζ)
FIG. 10. (Color online) Probability of the radial relative collision velocity. In panel (a) the full line corresponds to the PDF
as obtained when taking into account all collisions, as in the GCA. The dashed line gives the PDF of solely first collisions.
The figure has been obtained with a value of the Stokes number St = 1.0. The velocities are expressed in terms of the velocity
uK = (ηK /τK ) at the smallest scale. Panel (b) shows the ratio of the two mean radial relative collision velocities h|wr |ic,1 and
h|wr |ic,GCA for different Stokes numbers. The former takes into account only first collisions, the latter incorporates all collisions
detected within the framework of the GCA. Furthermore the inset in panel (b) shows h|wr |ic,1 and the radial relative velocity
of multiple collisions h|wr |ic,m in terms of uK .
101
100
10−1
10−2
10−3
10−4
10−5
10−6
10−7
ACKNOWLEDGMENTS
We thank Bernhard Mehlig and Kristian Gustavsson
for discussions, and Matthäus Bäbler for pointing out
[21]. AP and EL acknowledge the financial support of
ANR (contract TEC 2). The numerical calculations have
been performed at the PSMN computer center of the
Ecole Normale Supérieure in Lyon. MW and AP were
supported by the EU COST action MP0806 ‘Particles in
Turbulence’.
10−1
100
1/ζ
101
FIG. 11. (Color online) Comparison between the calculation
of the time ∆T spent by two particles together, by simulating
directly a gas of Maxwell particles (symbols), and the result
of the analytic formula (17) (continuous line). The different
symbols correspond to simulations with larger (circles) and
smaller (triangles) particles. The time ∆T is expressed in
terms of ζ, defined by (18). The behavior of the distribution
decays for large values of ∆T as ∼ ∆T −5 . This limit, as well
as the short time behavior ∼ ∆T are shown as dotted lines.
[1] P. G. Saffman and J. S. Turner, J. Fluid Mech. 1, 16
(1956).
[2] G. Falkovich, A. Fouxon, and M. G. Stepanov, Nature
419, 151 (2002).
[3] R. A. Shaw, Annu. Rev. Fluid Mech. 35, 183 (2003).
[4] M. Wilkinson, B. Mehlig, and V. Bezuglyy, Phys. Rev.
Lett. 97, 048501 (2006).
[5] V. S. Safranov, Evolution of protoplanetary cloud and formation of earth and planets (Israel program for scientific
translations, Keter Publishing House, Jerusalem, 1972).
[6] M. Wilkinson, B. Mehlig, and V. Uski, Astrophys. J.
Suppl. 176, 484 (2008).
13
[7] L.-P. Wang, A. S. Wexler, and Y. Zhou, Phys. Fluids
10, 266 (1998).
[8] J. Chun and D. L. Koch, Phys. Fluids 17, 027102 (2005).
[9] K. Gustavsson, B. Mehlig, and M. Wilkinson, New J.
Phys. 10, 075014 (2008).
[10] J. Bec, S. Musacchio, and S. S. Ray, Phys. Rev. E 87,
063013 (2013).
[11] G. Krstulovic, M. Cencini, and J. Bec, J. Stat. Phys. ??,
?? (2013).
[12] M. C. Jullien, J. Paret, and P. Tabeling, Phys. Rev.
Lett. 88, 2872 (1999).
[13] M. P. Rast and Pinton J. F., Phys. Rev. Lett. 107, 214501
(2011).
[14] R. Scatamacchia, L. Biferale, and F. Toschi, Phys. Rev.
Lett. 109, 144501 (2013).
[15] G. Falkovich and A. Pumir, J. Atmos. Sci. 64, 4497
(2007).
[16] G. P. Bewley, E.-W. Saw, and E. Bodenschatz, New J.
Phys. 15, 083051 (2013).
[17] K. Gustavsson and B. Mehlig, Phys. Rev. E 84,
045304(R) (2011).
[18] K. Gustavsson and B. Mehlig, (2013), arXiv:1309.3834
[physics.flu-dyn].
[19] M. Voßkuhle, E. Pumir, A. Lévêque, and M. Wilkinson,
submitted to Physical Review Letters, arXiv:1307.6853
(2013).
[20] M. Voßkuhle, A. Pumir, and E. Lévêque, J. Phys.: Conf.
Ser. 318, 052024 (2011).
[21] J. B. Jørgensen, J. Mann, S. Ott, H. L. Pécseli, and
J. Trulsen, Phys. Fluids 17, 035111 (2005).
[22] U. Frisch, Turbulence: The Legacy of A. N. Kolmogorov
(Cambridge University Press, Cambridge, 1995).
[23] A. G. Lamorgese, D. A. Caughey, and S. B. Pope, Phys.
Fluids 17, 05106 (2005).
[24] E. Calzavarini, R. Volk, M. Bourgoin, E. Lévêque, J.-F.
Pinton, and F. Toschi, J. Fluid Mech. 630, 179 (2009).
[25] S. A. Orszag, J. Atmos. Sci. 28, 1074 (1971).
[26] M. R. Maxey and J. J. Riley, Phys. Fluids 26, 883 (1983).
[27] R. Gatignol, J. Mc. Thor. Appl. 2, 143 (1983).
[28] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and
B. P. Flannery, Numerical Recipes: The Art of Scientific
Computing, 3rd ed. (Cambridge University Press, 2007).
[29] S. Sundaram and L. R. Collins, J. Comput. Phys. 124,
337 (1996).
[30] Y. Zhou, A. S. Wexler, and L.-P. Wang, Phys. Fluids
10, 1206 (1998).
[31] S. Sundaram and L. R. Collins, J. Fluid Mech. 335, 75
(1997).
[32] B. Rosa, H. Parishani, Ayala Orlando, G. Wojciech W,
and Wang Lian Ping, New J. Phys. 15, 045032 (2013).
[33] J. Bec, L. Biferale, G. Boffetta, M. Cencini, S. Musacchio,
and F. Toschi, Phys. Fluids 18, 091702 (2006).
[34] L.-P. Wang, A. S. Wexler, and Y. Zhou, J. Fluid Mech.
415, 117 (2000).
[35] J. Abrahamson, Chem. Eng. Sci. 30, 1371 (1975).
Appendix A
In this appendix we consider the equivalence of two
different approaches to determining the true rate of first
contact collisions, Γ1 , and therefore provide evidence that
Γ1 is indeed the physically appropriate definition of the
collision rate, for a class of systems where particles react
only at first contact. One approach is to use ghost particles and count the rate of multiple collisions, so that
Γ1 = ΓGCA − Γm . The alternative approach is to use one
of the two ‘substitution schemes’, described in subsection
II C, leading to rates ΓRe1 and ΓRe2 . In this appendix we
show that these estimates differ as a result of finite particle density effects. We provide evidence however that
they become equal in the limit of very dilute systems
(that is, as the particle density approaches zero).
Fig. 12 shows the collision rates produced by these algorithms, as well as Γ1 , obtained at several values of the
volume fraction Φ: Φ0 , Φ0 /2, and Φ0 /4. The numerical
results show that
(1) ΓRe2 (Φ) < ΓRe1 (Φ) < Γ1 = ΓGCA − Γm
(2) ΓRe2 (Φ/2) = ΓRe1 (Φ)
A more careful analysis of the data reveals that when
Φ → 0, both ΓRe1 and ΓRe2 tend to ΓGCA − Γm , with
corrections which are linear in Φ. Thus, our numerical
results demonstrate that the results from the three different algorithms agree in the dilute limit, Φ → 0.
The numerical results can be understood in the following way. We introduce the coefficients R = Γn, which
determine the rate of collision of a test particle moving
in a medium containing other particles with number density n. The ergodic assumption implies that the actual
collision rate can be inferred from the long-time behavior
of a single test trajectory.
The quantity nΓGCA is just the total rate of collision
along the test trajectory, and nΓm is the rate for collisions in which the test particle encounters the same target particle more than once. For systems where there is
coalescence on contact, only the rate of first collision is of
interest. The preferential concentration effect enhances
the rate of subsequent collisions with further particles,
once a pair of particles has collided. However, in the limit
n → 0, the time between collision events approaches infinity (∝ 1/n). For this reason, in the limit as n → 0 the
first contact collision rate is precisely nΓRei .
In terms of the long-time average over trajectories, the
algorithm which yields ΓRe2 changes the position of the
trajectory, therefore destroying the possibility of further
collisions with the surrounding environment. In addition,
a different realization of the surrounding background particles is chosen. In contrast, for the algorithm which
yields ΓRe1 , either the position of the test trajectory, or
the background of surrounding particles is changed, each
with a probability 1/2.
Changing the position of the test particles, or the surrounding background results first of all in preventing any
further collision between the pair of particles that came
into contact. The simplest approximation for the probability per unit time of a collision between the test particle
and a third particle, after a contact has been detected,
is proportional to R. This suggests that the approximations ΓRe1 or ΓRe2 miss a correction of order O(Rτ ),
14
(a)
(b)
1
ΓτK /(2a)3
80
60
40
20
0
0.95
ΓGCA (Φ0 ) − Γm (Φ0 )
ΓRe1 (Φ0 )
ΓRe2 (Φ0 )
ΓGCA (Φ0 /2) − Γm (Φ0 /2)
ΓRe1 (Φ0 /2)
ΓRe2 (Φ0 /2)
0
1
2
3
4
0.9
0.85
5
St
ΓRe1,Re2 (Φ0 )/Γ1
ΓRe1,Re2 (Φ0 /2)/Γ1
ΓRe1,Re2 (Φ0 /4)/Γ1
0
1
2
3
4
5
St
FIG. 12. (Color online) Comparison between Γ1 = ΓGCA − Γm , and the collision rates obtained by using the substitution
algorithms (ΓRe1 and ΓRe2 ). The left panel shows the raw data. The right panel shows ΓRei /Γ1 as a function of the Stokes
number, for those runs listed in Table I for which Φ0 = 4.5 10−5 . The open (filled) symbols correspond to the algorithm Re2,
(Re1). The full lines are deduced from each other by multiplication by a factor 2, and suggest that the dependence of the
difference Γ1 − ΓRei behaves linearly with the density of particles in the system.
where τ is a characteristic time over which the surrounding particles rearrange. Moreover, the algorithm yielding
ΓRe2 misses twice as many collisions with a third particle as ΓRe1 , which is also consistent with our numerical
observations. These arguments suggest that
ΓRei = Γ1 [1 − O(iRτ )] .
(A1)
As R ∝ n, this expression justifies the numerical observation that ΓRe1,2 are smaller than Γ1 , and differ from it
by a quantity proportional to n.