Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                

Multiple collisions in turbulent flows

2013, Physical Review E

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.