Relative velocities in bidisperse turbulent suspensions

We investigate the distribution of relative velocities between small heavy particles of different sizes in turbulence by analysing a statistical model for bidisperse turbulent suspensions, containing particles with two different Stokes numbers. This number, ${\rm St}$, is a measure of particle inertia which in turn depends on particle size. When the Stokes numbers are similar, the distribution exhibits power-law tails, just as in the case of equal ${\rm St}$. The power-law exponent is a non-analytic function of the mean Stokes number $\overline{\rm St}$, so that the exponent cannot be calculated in perturbation theory around the advective limit. When the Stokes-number difference is larger, the power law disappears, but the tails of the distribution still dominate the relative-velocity moments, if $\overline{\rm St}$ is large enough.

The dynamics of small heavy particles in turbulence plays a crucial role in many scientific problems and technological applications. Any model of the particle dynamics must refer to the turbulence the particles experience. This is a challenge for realistic modelling of such systems, for their direct numerical simulation (DNS), and for experiments. Novel particle-tracking techniques and improved DNS algorithms have made it possible to uncover striking phenomena in turbulent aerosols. For example, heavy particles tend to avoid the vortices of the turbulent fluid, and they form small-scale fractal patterns. Nearby particles can have very high relative velocities, an effect caused by 'caustic' singularities in the particle dynamics. The analysis of statistical models has led to substantial progress in explaining these phenomena, reviewed in Ref. [1]. This analysis has offered fundamental insights about how caustics shape the distribution of relative velocities of nearby particles [2][3][4][5][6][7][8][9][10].
These results apply only to 'monodisperse' suspensions of identical particles. The question is therefore how particles of different sizes cluster and move relative to each other. Spatial clustering of such 'bidisperse' suspensions was analysed in Refs. [11,12]. It was found that the particles cluster onto two distinct attractors, and that the fractal distribution of separations between differentlysized particles is cut off at a small spatial scale, r c , that can be much larger than the particle size.
In this Letter we analyse the distribution of relative velocities of particles with different sizes. Our results are important for the physics of turbulent aerosols, because the distribution of relative velocities determines the speeds of colliding particles, their collision rate, and collision outcomes [13][14][15].
In a dilute suspension of small, heavy, spherical particles the dynamics of a single particle is approximately given by Stokes law with constant γ ≡ 9ρ f ν/(2a 2 ρ p ): Here x and v are particle position and velocity, a is the particle size, ν is the fluid viscosity, and ρ f and ρ p are fluid and particle densities. We model the turbulent fluid Dotted lines show cut-off scales rc and vc. b Dependence of rc on θ (green, ), and of vc on θ (red, •), for St = 1. c (vr, r)/r d−1 evaluated at r < rc, for St = 1 and θ = 0.001 (green line), θ = 0.01 (solid black line), and θ = 0.1 (magenta line). Cross-over scales vc (arrows). d Power-law exponent µc versus St for small θ. Exponent for power-law tails in vr for fixed r with θ = 0.001 (green, ) and θ = 0.01 (red,•). Exponent for power-law tails in r for fixed vr with θ = 0.01 (blue,♦). Numerical data for min(D2, d + 1) where D2 is the phase-space correlation dimension (solid blue line).
velocities by a random Gaussian velocity field u(x, t) with zero mean, correlation time τ , correlation length η, and typical speed u 0 , representing the universal small spatial scales of turbulence [16], neglecting intermittent, non-Gaussian features [1]. Eq. (1) assumes that the particle and shear Reynolds numbers are small. Gravitational settling [17][18][19][20][21] is disregarded, this is valid when the turbulence is intense enough. The dimensionless parameters of the model are the Stokes number St ≡ 1/(γτ ), a measure of particle inertia, and the Kubo number, Ku ≡ u 0 τ /η, measuring the persistence of the flow.
Statistics of relative velocities and separations. Fig. 1a shows results of statistical-model simulations for Ku = 1 in two spatial dimensions. Panel a shows the distribution (v r , r) of relative radial velocity v r and separation r between particles with slightly different Stokes numbers. For fixed small v r , the distribution exhibits a power law as a function of r for large values of r, and becomes uniform for small values of r. This cross over defines the cut-off scale r c , observed and discussed in Refs. [11,12]. Fig. 1a shows that there is a second cut-off scale, v c , that distinguishes between a power law as a function of v r for large v r , and a plateau for small v r . Our numerical results demonstrate that this new scale v c does not depend on r (Fig. 1a). This is in marked difference to the monodisperse case, where the distribution has no plateau in the limit of r → 0. Panel b shows that v c depends linearly on the parameter θ ≡ |St 1 − St 2 |/(St 1 + St 2 ) measuring the difference between the Stokes numbers. Panel c shows distributions for r r c as functions of v r , for different values of θ. The dashed line is the power law |v r | µc−d−1 , where d = 2 is the spatial dimension. Panel d demonstrates that the exponent µ c is approximately equal to min(D 2 , d + 1) where D 2 (St) is the phase-space correlation dimension of a suspension of identical particles with Stokes number St ≡ (γτ ) −1 and γ = 1 2 (γ 1 +γ 2 ). The shape of the distribution shown in Fig. 1c has a strong effect on the moments of v r . For bidisperse suspensions, the uniform part of (v r , r) (r < r c and |v r | < v c ) yields an r d−1 -contribution to the moments |v r | p , just as the caustic-contribution for identical particles [22], caused by uniformly distributed, uncorrelated particles with large relative velocities originating from caustic singularities [3,4]. For small mean Stokes number, St 1, when caustics are rare, the dominant contribution to the moments of v r at small r comes from the uniform part of the distribution. But for St ∼ 1, caustics are abundant and for p ≥ 1 the moments |v r | p are dominated by the tails of the distribution.
The uniform part of the bidisperse distribution below r c and v c affects the tails indirectly. Comparing the normalised mono-and bidisperse distributions we see that the tails for θ = 0 must lie above those for θ = 0. This means that the moments |v r | p must have a minimum for θ = 0 (at fixed St and fixed small separation). This is consistent with the findings of Refs. [23][24][25][26][27].
Analysis of the white-noise limit. We now show that all of the above observations can be explained qualitatively by analysing a one-dimensional white-noise model. We is neglected here, this disregards preferential sampling which is not important in the white-noise limit [1]. The gradient A(t) has zero mean and correlation function . Neglecting the spatial dependence we write B(t) ≡ 2u(x, t). The two noise terms A(t) and B(t) are uncorrelated, and Dropping the tildes we find: where ∆v ≡ v 2 − v 1 is the relative velocity. The whitenoise limit is taken by letting St → ∞ and Ku → 0 such that ε 2 ≡ 3Ku 2 St stays finite. In this limit we obtain . The white-noise problem has three coupled dynamical variables, ∆x and ∆v as in the monodisperse case, and the mean velocity v ≡ 1 Cross-over scale v c . For small values of θ, the mean velocity v is approximately governed by an Ornstein-Uhlenbeck process driven by B. In the equation for ∆v we identify three competing terms: the two noise terms V 1 ≡ A∆x and V 2 ≡ θ(B − 2v), and the damping −∆v.
For small relative velocities (|∆v| → 0), the damping term is negligible. In this limit, the nature of the dynamics depends on whether |∆x| θ or |∆x| θ. In the first case V 2 dominates over V 1 , in the second case V 1 dominates. We find the spatial cut off r c by estimating the value of ∆x where V 1 and V 2 are of the same order. This yields r c ∝ Φ V2 /Φ A ∝ θ, where Φ X is the integrated correlation function of X. Now consider small separations (|∆x| r c ) and finite ∆v. In this case V 1 is negligible. Now we must consider whether |∆v| θ or |∆v| θ. In the first case the noise term V 2 dominates, so that the dynamics is diffusive. When |∆v| θ, on the other hand, the dynamics is deterministic, and the noise plays no role. Between these two regimes we find the cutoff scale v c ∝ Φ V2 ∝ θ. These considerations explain the linear θ-dependencies of the cutoff-scales r c and v c , Figs. 1c and Fig. 2a. Finally, consider the case where V 2 is negligible (|∆x| r c or |∆v| v c ). In this limit the dynamics is essentially that of monodisperse suspensions, but with mean Stokes number St. For |∆x| > r c or |∆v| > v c the distribution is expected to show the same power-law tails as the monodisperse case (Fig. 1c).
Power-law tails. We now discuss how the power-law tails observed in bidisperse suspensions for small θ are related to the power laws observed in the monodisperse case. For θ 1 we decompose the joint distribution P (∆v, ∆x,v) = 0 (∆v, ∆x)p(v) + δP (∆v, ∆x, v), where 0 (∆v, ∆x) equals the distribution for monodisperse particles (θ = 0). It follows from the analysis above that the relative dynamics is diffusive for |∆x| r c and |∆v| v c . In this case δP (∆v, ∆x,v) cannot be neglected. For |∆x| r c or |∆v| v c , by contrast, the dynamics is that of the monodisperse suspension and 0 (∆v, ∆x) domi-nates over δP (∆v, ∆x,v). To study the tails it thus suffices to consider 0 (∆v, ∆x). We follow Ref. [7] and make the separation ansatz: 0 (∆x, z) = µ A µ g µ (∆x)Z µ (z) with z ≡ ∆v/∆x . Inserting this ansatz into the Fokker-Planck equation corresponding to Eq. (2) we obtain g µ (∆x) = |∆x| µ−1 and d dz (z+z 2 +ε 2 d dz )Z µ = µzZ µ (3) to lowest order in θ. This is the equation for a suspension of identical particles [7] with Stokes number St. In a slightly different form, Eq. (3) was used to model the effect of wall collisions in turbulent suspensions [28].
In the bidisperse case, this divergence is regularised since diffusion dominates for |∆x| r c and |∆v| v c . More precisely, δP (∆v, ∆x,v) regularises P (∆v, ∆x, v). Thus a power-law steady-state distribution is obtained, even for ε < ε c . This is analogous to the effect of smallscale diffusion that regularises a power-law steady state distribution of separations [30].
The full solution of (3) consists of a sum of two terms, A 0 ∆x −1 Z 0 (z) + A µc |∆x| µc−1 Z µc (z). As was discussed in [7], the separation ansatz for (∆x, z) is strictly valid only for |∆x| 1. The A 0 -term is therefore subdominant for µ c < 0. For µ c > 0 we expect A 0 ∼ θ A µc since A 0 must vanish in order for the distribution to be normalisable over ∆x = 0 as θ → 0, consistent with the numerical small-θ results in Fig. 2a. Therefore we disregard the µ = 0-term for all values of ε. Fig. 2b shows that for θ 1 the exponents of both the ∆x-and ∆v-distributions (symbols) follow the numerically calculated µ c (ε) very precisely (solid line). The asymptote Z µc (z) ∼ |z| µc−2 for large |z| gives rise to the power-law tails in the relative-velocity distribution of the form |∆v| µc−2 . In d spatial dimensions a similar argument yields tails on the form v µc−d−1 r . This explains the power laws observed in Figs. 1c and 2a.
In summary the white-noise analysis qualitatively explains not only the existence of the plateau in the relativevelocity distribution and the corresponding cutoff scales r c and v c . It also yields the linear dependence of these scales upon θ, and explains the power-law tails. Thus, the one-dimensional white-noise analysis qualitatively explains all the observed features in Fig. 1.
Failure of perturbation theory in ε. Fig. 2b demonstrates that the exponent µ c depends very sensitively on ε, for small values of ε. Here we show that this dependence is non-analytic. This is important, because it means that perturbation theory in ε fails. This is, in turn, a likely reason [31] why Borel resummation of perturbation theory does not yield accurate results for the correlation dimension D 2 (Fig. 1 in Ref. [32]). There is to date no theory explaining these observations.
The exponent µ c appears as a generalised eigenvalue in Eq. (3), but exact solutions to Eq. (3) are known only for µ = 0 [29], µ = −1 [33], and in the large ε-limit [34]. The physical solution at small values of ε corresponds to the eigenvalue µ c (ε). The solution for µ = −1 is unphysical since it does not obey particle-exchange symmetry.
Numerical analysis reveals that µ c (ε) → −1 as ε → 0. We write µ c = −1 + δµ, where δµ ≥ 0 is an 'eigenvalue splitting' that vanishes as ε → 0, and attempt to find the physical solution by perturbation theory in δµ: Substituting (4) into (3) we obtain a hierarchy of differential equations for Z (n) (Supplemental Material [35]). The lowest-order solution is the unphysical one, Z (0) = Z −1 . Its general form is known [33]: Here U (z) ≡ ( 1 3 z 3 + 1 2 z 2 )/ε 2 , and C 1 and C 2 are constants. The first-order solution reads: The integrand involves Z −1 , and we must take C 1 = 0 in Eq. (5) for Z −1 since the first term in (5) is not normalisable. We further choose C 2 = −ε 2 because it normalises the large-z tails of Z −1 to Z −1 ∼ −z −3 . To evaluate (6) we use the fact that δµ → 0 as ε → 0 and apply a WKB approximation [36] of Z −1 that becomes exact in the limit ε → 0. The equation for Z −1 has two 'turning points' at z = −1, 0 where the WKB approximation breaks down. We find locally exact solutions and use these to match the WKB solutions across the turning points. Details of this asymptotic-matching method are given in the Supplemental Material [35]. In this way we find Z (1) ∼ |z| −3 log |z| for large negative values of z, and Z (1) ∼ |z| −3 [2πe 1/(6ε 2 ) − log |z|] for large positive z. Imposing particle-exchange symmetry yields This result is shown as a dashed line in Fig. 2b. It is in good agreement with the numerical data for small ε.
We emphasise that the non-analytic dependence (7) cannot be obtained by the ε-perturbation theory described in Refs. [1,37,38]. That approach starts by rescaling z → εz, and solves the resulting equation d dz (z + εz 2 + d dz )Z µ = εµzZ µ perturbatively in ε. The solution remains essentially Gaussian. Even for very small ε, however, the z-dynamics can escape from z ≈ 0 to −∞ by forming caustics [1]. This gives rise to power-law tails that are not described by the local theory. Higher orders in ε improve the accuracy only locally, inside an ε-sized 'boundary layer' [39] around z = 0, but not in the tails. Consequently, perturbation theory fails to give the non-analytical dependence (7).
Our theory, by contrast, provides a uniform approximation valid for small ε and arbitrary z. It allows to compute the exponentially small eigenvalue splitting δµ because it captures the formation of caustics. For small ε, the rate J of caustic formation is J ∼ (2π) −1 e −1/6ε 2 [1]. It follows that µ c ∼ 2J − 1 at small values of ε.
We now show that caustics cause similar problems in small-St expansions [1,11,[40][41][42][43]. In one dimension, the equation for z reads in dimensional units The Stokes number used in DNS of heavy particles in turbulence [44,45] is κ ≡ u 0 /(ηγ) [1]. Therefore we consider the small-κ limit of (8), taking A to be constant (in this limit the local fluid-velocity gradients are roughly constant during the time η/u 0 [1]). As , localised around zero. As in the white-noise case, this local approximation fails when caustics allow z to escape to −∞. Kramers-escape theory for coloured noise [46,47] reveals that caustics occur in the persistent limit when A < −γ/4, with probability p = −γ/4 −∞ dA P (A). Evaluating this integral for κ 1, one obtains a non-analytic dependence, p ∝ exp[−1/(96κ 2 )], consistent with Ref. [48]. In DNS a similar non-analytic dependence is found [49], albeit of a slightly different form because the turbulent velocitygradients are non-Gaussian. We conclude that the perturbation theories in the white-noise limit and in the persistent limit fail for the same reason. Both expansions are valid only near z = 0 and fail in the presence of caustics.
Conclusions. We analysed the distribution of relative velocities in turbulence between small, heavy particles with different Stokes numbers St. We demonstrated that the difference in St causes diffusive relative motion at small separations, giving rise to a plateau in the distribution for relative velocities smaller than a cut off v c . This is in qualitative agreement with DNS [50]. Fig. 2

in
Ref. [50] shows the v r -distribution for a bidisperse suspension. The cut off v c depends linearly on θ, the parameter characterising the difference in Stokes numbers.
At small values of θ the distribution exhibits algebraic tails |v r | µc−d−1 as in the monodisperse case. The exponent µ c is determined by the phase-space correlation dimension of a monodisperse system with mean Stokes number St. When St is O(1) or larger, the moments of relative velocities at small separations are dominated by the tails of the distribution. At small θ the tails are power laws, but as θ becomes larger, our statistical-model simulations show that the tails gradually disappear, consistent with v c ∝ θ. The DNS of Ref. [50] are for quite large θ. It would be of interest to perform DNS for small values of θ to find the power-law tails we predict here.
Our analysis shows how sensitive the distribution is to polydispersity. This is important for experiments tracking the dynamics of µm-size particles in turbulence [9], where strict monodispersity is difficult to achieve.
We explained these observations using a one-dimensional statistical model. The difference in Stokes numbers regularises the distribution, so that the bidisperse model has a power-law steady-state at small ε. This makes it possible to compute how the power-law exponent µ c depends on the inertia parameter ε. We demonstrated that the dependence of µ c upon ε is not analytic.
Our theory explains why small-ε expansions fail to give non-analytic contributions to physical quantitites in the white-noise limit. It is likely that non-analytic terms in white-noise expansions for the correlation dimension [31] and Lyapunov exponents [37,38] in two and three dimensions have similar origins, and we speculate that Eq. (7) contains only the first two terms in an infinite series of the form [51] ∞ k=0 e −k/(6ε 2 ) ∞ m=0 b (k) m ε 2m . There could also be logarithmic terms, log n ε 2 .
More generally we explained that small-St expansions for heavy particles in turbulence [1,11,[40][41][42][43] suffer from similar problems when caustics occur. This indicates that matched asymptotic expansions (or similar methods) are required to explain the characteristic minimum [45] of the correlation dimension as a function of St.
Our predictions can be directly tested by experiments or by DNS of turbulent bidisperse turbulent suspensions. We note, however, that our analysis pertains to the dynamics in the dissipative range. At higher Stokes numbers, when separations between particle pairs explore the inertial range [52,53], we expect corrections to the powerlaw exponents derived here.