Articles | Volume 22, issue 6
Research article
22 Mar 2022
Research article |  | 22 Mar 2022

Intricate relations among particle collision, relative motion and clustering in turbulent clouds: computational observation and theory

Ewe-Wei Saw and Xiaohui Meng

Considering turbulent clouds containing small inertial particles, we investigate the effect of particle collision, in particular collision–coagulation, on particle clustering and particle relative motion. We perform direct numerical simulation (DNS) of coagulating particles in isotropic turbulent flow in the regime of small Stokes number (St=0.001–0.54) and find that, due to collision–coagulation, the radial distribution functions (RDFs) fall off dramatically at scales rd (where d is the particle diameter) to small but finite values, while the mean radial component of the particle relative velocity (MRV) increases sharply in magnitude. Based on a previously proposed Fokker–Planck (drift-diffusion) framework, we derive a theoretical account of the relationship among particle collision–coagulation rate, RDF and MRV. The theory includes contributions from turbulent fluctuations absent in earlier mean-field theories. We show numerically that the theory accurately accounts for the DNS results (i.e., given an accurate RDF, the theory could produce an accurate MRV). Separately, we also propose a phenomenological model that could directly predict MRV and find that it is accurate when calibrated using fourth moments of the fluid velocities. We use the model to derive a general solution of RDF. We uncover a paradox: the past empirical success of the differential version of the theory is theoretically unjustified. We see a further shape-preserving reduction of the RDF (and MRV) when the gravitational settling parameter (Sg) is of order O(1). Our results demonstrate strong coupling between RDF and MRV and imply that earlier isolated studies on either RDF or MRV have limited relevance for predicting particle collision rate.

1 Introduction

The motion and interactions of small particles in turbulence have fundamental implications for atmospheric clouds; specifically, they are relevant for the timescale of rain formation, particularly in warm clouds (Falkovich et al.2002; Wilkinson et al.2006; Grabowski and Wang2013) (a similar problem also applies to planet formation in astrophysics; Johansen et al.2007). They are also important for engineers who are designing future, greener combustion engines, as this is a scenario they wish to understand and control in order to increase fuel efficiency (Karnik and Shrimpton2012). Cloud particles or droplets, due to their inertia, are known to be ejected from turbulent vortices and thus form clusters – regions of enhanced particle density (Wood et al.2005; Bec et al.2007; Saw et al.2008; Karpińska et al.2019); this together with droplet collision is of direct relevance for the mentioned applications. Due to the technical difficulty of obtaining extensive and systematic experimental or field data on particle/droplet collision in turbulent cloud, many of the recent studies rely on direct numerical simulation (DNS), examples of which can be found in, e.g., Onishi and Seifert (2016) and Wang et al. (2008), and reference therein. Up until now, we have not had definitive answers to basic questions such as how to calculate the particle collision rate from basic turbulence-particle parameters and what the exact relation between collision and particle clustering and/or motions is, for, as we shall see, our work reveals that collision–coagulation causes profound changes in particle relative velocity statistics and particle clustering, questioning earlier understanding of the problem. The difficulty of this problem is in part related to the fact that turbulence is, even by itself, theoretically virtually intractable due to its nonlinear and complex nature.

The quest for a theory of particle collision in turbulence started in 1956 when Saffman and Turner (1956) derived a mean-field formula for the collision rate of finite-size, inertia-less particles. In another landmark work (Sundaram and Collins1997), a general relation among the collision rate (Rc), particle clustering and mean particle relative radial velocity was presented: Rc/(n1n2V)=4πd2g(d)wr(d), where g(r) is the particle radial distribution function (RDF), wr is the radial component of relative velocity between two particles, denotes averaging over particle pairs, wr(d)〉 is the mean radial component of relative particle velocity (MRV), ni denotes the global averages of particle number density, V is the spatial volume of the domain and d is the particle diameter. The remarkable simplicity of this finding inspired a “separation paradigm”, which is the idea that one could study the RDF or MRV separately (which is technically easier), and the independent results from the two may be combined to accurately predicts Rc (an idea that we subsequently challenge). Another work of special interest here is the drift-diffusion model by Chun et al. (2005) (hereafter CK theory) (note that there are other similar theories: Balkovsky et al.2001; Zaichik and Alipchenkov2003). The CK theory, derived for non-colliding particles in the limit of a vanishing particle Stokes number, St (a quantity that reflects the importance of the particle's inertia in dictating its motion in turbulence), correctly predicted the power-law form of the RDF (Reade and Collins2000; Saw et al.2008) and has seen remarkable successes over the years, including the accurate account of the modified RDF of particles interacting electrically (Lu et al.2010) and hydrodynamically (Yavuz et al.2018).

Here, we first present results on RDF and MRV for particles undergoing collision–coagulation1. The data are obtained via direct numerical simulation (DNS), which is the gold-standard computational method in terms of accuracy and completeness for solving the most challenging fluid dynamics problem i.e., turbulent flows. It is worth noting that the focus of our work is on the fundamental relationship between collision, RDF and MRV and on highlighting differences from the case with non-colliding particles (Chun et al.2005). To that end, we have designed the DNS to have an idealized setup similar to what was done in Chun et al. (2005), which would allow us to identify the effects of particle collision–coagulation without doubt. As a result, this limits the direct applicability to real systems (these limitations are detailed in Sect. 4.5).

Analysis of the DNS results is followed by a theoretical account of the relations between collision rate, RDF and MRV, which includes mean-field contributions (Saffman and Turner1956; Sundaram and Collins1997) and contributions from turbulent fluctuations (absent from earlier theories, Saffman and Turner1956; Sundaram and Collins1997). The theory is derived from the Fokker–Planck (drift-diffusion) framework first introduced in the CK theory (Chun et al.2005). We shall see that the main effect of collision–coagulation is the enhanced asymmetry in the particle relative velocity distribution2 and that this leads to nontrivial outcomes.

2 Direct numerical simulation (DNS)

To observe how particle collision–coagulation affects RDF and MRV, we performed direct numerical simulation (DNS) of steady-state isotropic turbulence embedded with particles of finite but sub-Kolmogorov size. We solve the incompressible Navier–Stokes equations (Eq. 1) using the standard pseudo-spectral method (Rogallo1981; Pope2000; Mortensen and Langtangen2016) inside a triply periodic cubic box.


where ρ, p, ν and f are the fluid mass density, pressure, kinematic viscosity and imposed forcing respectively. The velocity field is discretized on a 2563 grid. Aliasing resulting from Fourier transform of truncated series is removed via a 2/3 dealiasing rule (Rogallo1981). A statistically stationary and isotropic turbulent flow is achieved by continuously applying random forcing to the lowest wave numbers until the flow's energy spectrum is in steady state (Eswaran and Pope1988). Second-order Runge–Kutta time stepping was employed. Further details of such a standard turbulence simulator can be found in, e.g., Pope (2000), Rogallo (1981) and Mortensen and Langtangen (2016). The accuracy of DNS for turbulent flows has been experimentally validated for decades (see, e.g., the compilation of results in Pope2000).

Table 1Values of the parameters in the DNS. (Note that dm denotes decimeters.) From the left, we have the Taylor-scale Reynolds number, kinematic viscosity of the fluid, root-mean-square of fluid velocity, kinetic energy dissipation rate, Kolmogorov length scale and timescale, length of the simulation cube and particle diameters considered. We have introduced d to represent the specific value 9.49×10-4dm (more details in the text). We choose the units of the length scale (timescale) in the DNS to be in decimeters (seconds), such that ν is nearly its typical value in the atmosphere.

Download Print Version | Download XLSX

Particles in the simulations are advected via a viscous Stokes drag force (Maxey and Riley1983):


where u and v are the local fluid and particle velocity respectively, and τp is the particle inertia response time, defined as τp=118(ρp/ρ-1)(d2/ν), where ρp is the particle mass density, and d is the particle diameter. As mentioned, this work focuses on the fundamental relationship between collision–coagulation, RDF and MRV, as well as on addressing the validity of the theory (to be described). It is thus, beneficial to keep the DNS setting idealized (and in the regime relevant for the theory) for the sake of clarity when interpreting results. To that end, the DNS does not include inter-particle hydrodynamic interactions (HDIs) and gravitational settling, nor does it consider the effects of temperature and humidity variation and phase transitions. Such practice is not uncommon in studies designed to isolate and address fundamental issues related to particles dynamics in turbulence; examples that are closely related to the current setup and/or problem include Sundaram and Collins (1997), Chun et al. (2005), Bec et al. (2007), Salazar et al. (2008), Wang et al. (2008), Woittiez et al. (2009) and Voßkuhle et al. (2013). However, such an approach certainly limits the direct applicability of our results to some realistic problems in the atmosphere; these limitations will be detailed in Sect. 4.5, where a discussion of the effects of gravity and HDI is also given.

In this context, the particle Stokes number, defined as τp/τη, where τη is the Kolmogorov timescale, can be expressed as St=118(ρp/ρ-1)(d/η)2, where η is the Kolmogorov length scale. Time stepping of the particle motion is done using a second-order modified Runge–Kutta method with an “exponential integrator” that is accurate, even for τp much smaller than the fluid's time step (Ireland et al.2013). The particles introduced into the simulation are spherical and are of the same size, the initial number of particles is 107 and they are randomly distributed in space. Particles collide when their volumes overlap and a new particle is formed conserving volume and momentum (Bec et al.2016). We continuously, randomly, inject new particles so that the system is in a steady state after some time. Statistical analysis is done at steady state on monodisperse particles (i.e., particles with the same St). Experimental validation of the accuracy of such a particle simulation scheme in DNS can be found in Salazar et al. (2008), Saw et al. (2012, 2014) and Dou et al. (2018).

Values of key parameters of the DNS are given in Table 1. Values of other parameters and further details can be found in the Supplement.

3 Elements of the drift-diffusion theory

As described in Chun et al. (2005), in the limit of St≪1, particle motions are closely tied to the fluid velocity and, to leading order, completely specified by the particle position and fluid velocity gradients. We consider the Fokker–Planck equation which is closed and deterministic (see, e.g., Appendix J in Pope2000):

(2) P t + ( W i P ) r i = 0 ,

where PP(ri,t|Γij(t)) is the (per volume) probability density function (PDF) for a secondary particle to be at vector position ri relative to a primary particle at time t, conditioned on a fixed and known history of the velocity gradient tensor along the primary particle's trajectory Γij(t), and Wi is the mean velocity of secondary particles relative to the primary, under the same condition. Note that Wi is a conditional average, while wi denotes a realization of relative velocity between two particles.

From this, one could derive an equation for P:

(3) P t + r i ( W i P + W i P ) = 0 ,

where 〈.〉 implies ensemble averaging over primary particle histories (note that Wr〉≡ unconditional mean of wr, averaged over all particle pairs, i.e., the MRV). This equation, however, is not closed due to the correlation between the fluctuating terms Wi and PP-P. The correlation WiP can be written in terms of a drift flux and diffusive flux (detailed derivation is well described in Chun et al.2005), such that we have

(4) P t + r i ( q i d + q i D ) + ( W i P ) r i = 0 ,

where the drift flux is

(5) q i d = - - t W i ( r , t ) W l r l ( r , t ) P ( r , t ) d t ,

and the diffusive flux is

(6) q i D = - - t W i ( r , t ) W j ( r , t ) P r j ( r , t ) d t ,

where r satisfies a characteristic equation: ri/t=Wi(r,t), with boundary condition ri=ri at t=t.

Finally we note that since particles are allowed to collide–coagulate in our theory, we use the conventional definition of MRV: Wrwr. In some works that consider non-colliding (ghost) particles, the conditional mean wr|wr0 must be used for the purpose of calculating the mean collision rate, since wr=0 in isotropic turbulence (Chun et al.2005).

4 DNS results, theory and discussion

We compute the RDF via g(r)=Npp(r)/[12N(N-1)δVr/V], where Npp(r) is the number of particle pairs found to be separated by distance r, and δVr is the volume of a spherical shell of radius r and infinitesimal thickness δr.

Figure 1RDFs (g(r)) of particles that coagulate upon collision. (a) g(r) for various cases of Stokes numbers and particle diameters (d). : St=0.22, d=d, : St=0.54, d=d, St=0.54, d=2d. All g(r) values drop off exponentially when rd (details in text). Inset: g(r) versus rd for the case. This exemplifies the fact that limrdg(r) is nonzero. (b) RDFs versus rd1 (where d1=0.99d) for the case of St=0.054, d=d. : the raw DNS-produced RDF (gDNS(r)). Red line: power-law fit to gDNS(r) (i.e., the plot) in the large r regime (the fit result is 0.890r−0.0535). It is equivalent to gs(r) in the ansatz ga(r)=g0(r)gs(r); i.e., it is the expected RDF for non-colliding particles under the same conditions. : the compensated RDF, i.e., gDNS(r)/gs(r) (note that gs(r) is the red line described earlier); this essentially gives us g0(r), which may be understood as a “modulation” on the RDF due to collision–coagulation. Cyan line: two-piece power-law fits to the compensated RDF (the plot) in the small and large rd1 regimes respectively (fit results: 4.17(rd1)0.212, 1.00(r-d1)-2×10-4); this is an estimate for g0(r). Dashed black line: g0(r)gs(r), (cyan line × red line); this shows that the ansatz accurately reproduces gDNS(r). Inset: RDFs versus rd. : compensated g(r) for St=0.054, d=d, equivalent to the plot in the panel's main figure; : compensated RDF for case St=0.001, d=d, i.e., finite size, almost zero St particles (in this case, the compensated and raw RDFs are the identical). This inset suggests that g0(r) has negligible St dependence.


Figure 1 shows the RDFs obtained for monodisperse particles of various Stokes numbers and sizes. Two cases (St=0.22 and 0.54) are shown in panel a, and two more (St=0.054 and 0.001) are shown in panel b. In this work, we focus on the smaller values of St since the theory which we shall consider is also only applicable in the St≪1 regime. However, we have included the St=0.54 case to demonstrate that the observations to be described also extends to finite St. In all cases, except one, the particles are of the same size d=d, where d represents the specific value of d=9.49×10-4dm, chosen so that the particle sizes are about O(0.1) times the Kolmogorov scale (η), thus allowing us to still observe a regime (3dr30η) of power-law RDFs. To shows the effect of changing particle size, panel a also includes a case of St=0.54,d=2d for comparison. Looking at panel a, apart from the apparent power-law behavior of the RDFs at intermediate values of r, the most striking feature of these RDFs for colliding–coagulating particles is that they fall off dramatically in the rd regime. This is very different from what was seen in earlier studies of non-colliding particles, where g(r) are simple power laws (Chun et al.2005; Saw et al.2008). We also see that as r approaches d, the steepness of the curve (see, e.g., the blue circles) increases as g(r) drops off; this and the fact that the abscissa is logarithmic implies that g/r is increasing exponentially in the process. As a consequence, it is difficult to discern from these plots if the limit of g(r) at particle contact (rd) is still nonzero. This is an important question as limrd[g(r)]=0 implies that the mean-field formula of Sundaram and Collins (1997) has zero contribution towards Rc; i.e., the collision rate is solely due to turbulent fluctuations. It is only by re-plotting g(r) versus rd (see insets in Fig. 1) and using a remarkable resolution that is 103 finer than d that we see a convincing trend supporting a finite g(rd). Also clear in panel a is the observation that with changing particle size (d), the location of the sharp fall-off merely shifts to where the new value of d is.

The strong effect of particle collision on the RDF (also on MRV as we shall see later) challenges the validity of the separation paradigm. We note that similar fall-off of RDF was previously observed (Sundaram and Collins1997), but a complete analysis and theoretical understanding were lacking. Also, a study on multiple collisions (Voßkuhle et al.2013) had hinted at the potential problem with the separation paradigm.

Another observation is that in the power-law regime (3dr30η), the RDFs appear (as expected) as straight lines with slopes (i.e., power-law exponents) that increase with St and are numerically consistent with those found for non-colliding particles (see, e.g., Saw et al.2012).

4.1 Theoretical account via drift-diffusion theory

To theoretically account for the new findings, we make some derivations that are partially similar to the ones in Chun et al. (2005) but under a new constraint due to coagulations: at contact (r=d), the radial component of the particle relative velocities can not be positive3, while with increasing r the constraint is gradually relaxed. The first consequence of this is that the distribution of the radial component of the relative particle velocity (Wr) is highly asymmetric at rd; i.e., the PDFs of positive Wr values are very small (this constitutes the “enhanced asymmetry” mentioned earlier). Thus for rd, Wr must be negative. In Sect. 3, we showed that in the St≪1 limit, one could derive a master equation (Eq. 4), reproduced here for clarity:


where qid is the drift flux (of probability due to turbulent fluctuation) defined in Eq. (5), and qiD is the diffusive flux defined in Eq. (6).

We then expand Wi, Wl/rl and (consequentially) the fluxes as perturbation series with St as the small parameter (details in the Supplement or Chun et al.2005). The coagulation constraint affects the values of the coefficients of these series. For the drift flux, the leading-order terms (in powers of St) are

(7) q i d = - P ( r ) r k - t a i k ( 1 ) St + a k i ( 2 ) St 2 d t ,

with aik(1)=τηΓik(t)Γlm(t)Γml(t) and aki(2)=τη2Γij(t)Γjk(t)Γlm(t)Γml(t), where Γij is the ijth component of the fluid's velocity gradient tensor at the particle position (the aik values are thus related to two time correlations of moments of velocity gradients; Chun et al. (2005) showed that aik(2)S2-R2, where (S2,R2) denotes the average fluid (strain rate tensor, rotation rate tensor) squared at particle positions). As explained earlier, coagulation constraint causes the PDF of relative particle velocities to become highly asymmetric for rd; thus aik(1) is nonzero at this scale. This is very different from the case of non-colliding particles (Chun et al.2005), where aik(1) is always zero due to statistical isotropy. Under the constraint, DNS gives -taik(1)dt-0.18s−1 and -taki(2)dt2.45s−1 (more in the Supplement). Thus for rd, the drift flux is negative for large St but becomes positive4 when St is below the value of ≈0.07; and in the limit of St→0, it is dominated by the first term in Eq. (7).

The term qiD is a “nonlocal” diffusion caused by fluctuations and can be estimated using a model that assumes the particle relative motions are due to a series of random uniaxial straining flows (Chun et al.2005). Chun et al. (2005) showed that, generally, qiD has an integral form (due to non-locality), and only in the special case where g(r) is a power law may it be cast into a differential form (similar to a local diffusion). In view of the nontrivial g(r) observed here, we must proceed with the integral form:


where R0r0/r, with r0 as the initial separation distance of a particle pair before a straining event. F is the probability density function for the duration of each event, fI is determined by relative prevalence of extensional versus compressional strain events (more details in the Supplement or Chun et al.2005) and Ω is the solid angle for the axis of the straining flow. Note that due to coagulation, the R0 integration starts from d/r. We differ crucially from the CK theory via the introduction of the (positive) factor cst, which can be shown to be equal to |c1|, where c1 is the power-law exponent of the RDF the particles would have assuming they are non-colliding (details in the Supplement).

By definition, g(r)≡αP. Periodic boundaries in our DNS imply that α=V (more in the Supplement). Using this and the fact that the problem only has radial (r) dependence, we rewrite Eq. (4) as

(8) r 2 g ( r , t ) t + r r 2 α ( q i d + q i D ) + r 2 W r g ( r , t ) = 0 ,

where the content inside [.] gives the total flux. For a system in steady state, the first term in Eq. (8) is zero, and upon integrating with limits [d,r], we have


where we have identified the total flux at contact (r=d) as the negative of the (always positive) normalized collision rate Rc*Rc/4π[N(N-1)/2]/V, and comparing with Eq. (7), we see that

(10) A τ St - t a i k ( 1 ) d t + St 2 - t a k i ( 2 ) d t ,

with the specific values of the t integrals already given above. For clarity, we reiterate that on the left side of Eq. (9), we have the diffusive flux (qrD), mean-field flux (r2g(r)〈Wr) and drift flux (qrd), while on the right, the total flux is given in terms of the normalized collision rate (Rc*). We note that this equation embodies the full relationship among RDF, MRV and collision rate.

4.2 Ansatz and accuracy of the theory

Simple analytical solutions to Eq. (9) may be elusive due to its integral nature (a consequence of the non-local diffusive flux). However, one could gain insights into it and test its accuracy via numerical solutions. To that end, we begin with a simple ansatz for g(r), then we curve-fit the ansatz to the DNS-produced RDF (gDNS(r)). This enables us to, firstly, verify that the ansatz could accurately represent gDNS(r) and, secondly, to obtain a “calibrated” ansatz that is a numerically accurate representation of gDNS(r). We then show that Eq. (9), supplied with the calibrated ansatz, could numerically predict Wr〉(r) (i.e., MRV) that agrees well with the DNS-produced MRV. In short, we will show that given a “correct” g(r), Eq. (9) produces the “correct” Wr〉(r).

The ansatz has the form ga(r)=g0(r)gs(r), with gs(r)=c0r-c1, i.e., the RDF form for non-colliding particles (Chun et al.2005) under the same conditions. As a first-order analysis, we let g0, which embodies the effects of collision, take the simplest form that could still capture the main features of the RDFs seen in Fig. 1. Specifically, we let g0(r)=c00(r-d1)c10, where c00(r) and c10(r) are each piecewise constant quantities that switch from their small r to large r values at a crossover-scale rc (of the order of d); i.e., g0 is a two-piece power law of rd1. (Note that our earlier finding of g(rd)>0 implies that d1<d.)

From a given DNS-produced RDF (gDNS(r)), we first obtain a calibrated gs by fitting c0r-c1 to gDNS(r) in the power-law regime dr10η (see the red line in Fig. 1b). Next, we compute the DNS estimate of g0 via g0DNS=gDNS(r)/gs(r), which is essentially a compensated RDF (see the plot in Fig. 1b). To get a calibrated g0, we then fit the general form of g0 given above to g0DNS (see the cyan line in Fig. 1b; note that each time, two pieces of power laws are fitted to one g0DNS, and rc results naturally from the intersection of the two). Figure 1b shows the calibrated ansatz for the case of St=0.054 and verifies its accuracy (the red line is gs(r), the cyan line is g0(r) and the dashed black line (g0gs) accurately reproduces gDNS(r)). The inset in Fig. 1b shows that g0DNS(r) (i.e., gDNS(r)/gs(r)) is roughly St-independent for St≪1.

Next, we numerically evaluate the integral in the first term of Eq. (9). The St≪1 assumption allows us to approximate g(r,St) inside the integral by its zero-St cousin g(r,St0) (Chun et al.2005). In practice, we replace g(r,St) with the ansatz fitted to the DNS result of g(r,St=0.001). Next, we use the DNS data to estimate Aτ and compute Rc* and cst (for this case, DNS gives Rc*=9.69×10-10dm3s-1; cst=|c1| as mentioned earlier). Finally we use Eq. (9) to predict Wr〉(r).

Figure 2Mean radial component of relative velocity (MRV) for particles of specific Stokes numbers and some theoretic numerical predictions. (a) Symbols are DNS results with : St=0.001; : St=0.054; : St=0.11. The lines are the numerical predictions by the theories (Eqs. 9 or 13) using the ansatz (details in text). Orange line: Wrrd,St=0.054theory, i.e., the numerical prediction via the integral version of the theory (Eq. 9) for the small r regime (rd). Black line: Wrrd,St=0.054theory, same as the previous but for the large r regime (rd). Green line: prediction of the differential version of the theory (Eq. 13) for the rd regime. Inset: a repeat of the main figure in log–log axes. Exception: cyan line is the prediction of the differential version of the theory but for the rd regime. (b) MRV compared with predictions via the phenomenological model of particle approach angles (Eqs. 11 and 12). DNS results: : St=0.001; : St=0.054. Dotted lines are model predictions of WrSt=0 using Eqs. (11) and (12) with variance K obtained by matching the model's and DNS's transverse-to-longitudinal ratio of structure functions (TLR) of a certain order (from the top: yellow line – order 2, green line – order 4, cyan line – order 6).


Comparison of the predicted Wr〉(r) with the values obtained directly from the DNS is shown in Fig. 2. The prediction shown was made for the case of St=0.054, to be compared with its DNS counterpart (the symbols). (We also show the DNS result for St=0.001 and 0.11 to highlight an observation that Wr〉(r) is almost St-independent in this small St regime.) We showed earlier that for rd, Aτ is given by Eq. (10). However, as stated earlier, as r increases, the (statistical) asymmetry induced by collision–coagulation gradually becomes subdominant to the isotropy of turbulent fluctuation. Statistical isotropy implies aik(1)=0 (Chun et al.2005), a fact our DNS data confirm. Thus, for rd, Aτ equals the order St2 term in Eq. (10), exactly the same as the results of Chun et al. (2005) for non-colliding particles. For this reason, we show two versions of the prediction: Wrrdtheory and Wrrdtheory, which are respectively obtained by setting Aτ to its small r and large r limits (-2.6×10-3s−1, 7.1×10-3s−1) respectively. The agreement between DNS and the predictions is noteworthy, especially for small r. At r≈2d, the DNS result shows a weak tendency to first follow the upward trend of Wrrdtheory and then drops off significantly at r2.5d. The latter is consistent with the fact that Wrrdtheory is below Wrrdtheory, but the drop is sharper than predicted.

4.3 Phenomenological model of MRV

Alternatively, Eq. (9) may be solved for the correct g(r), if Wr is given. As we are assuming St≪1, particle velocity statistics may be approximated by their fluid counterparts (Chun et al.2005); i.e., we may replace Wr with WrSt=0, the latter being the MRV of fluid particles. Hence, if WrSt=0 is known, it may be used, together with Eq. (9), to predict RDF of any finite but small St. Figure 2a shows that WrSt>0 from the DNS do not change significantly for St[0.001,0.1], supporting this approach5.

Here we provide a simple, first-order model for WrSt=0. We limit ourselves to the regime of small particles (dη) and anticipate that Wr is non-trivial (nonzero) only for rd, a fact observable in Fig. 2a. We also assume that the relative trajectories of particles are rectilinear at such small scales. The coagulation constraint then implies that: in the rest frame of a particle (call it P1), a second particle nearby must move in such a way that the angle (θ) between its relative velocity and relative position (seen by P1) must satisfy sin-1(d/r)θπ, under the convention of sin-1(x)[-π2,π2] (more in the Supplement). We can thus write (by treating negative and positive wr separately, applying the K41 theory (Kolmogorov1941) and the bounds on θ, details in the Supplement), for St≪1, that


where . denotes averaging over particle pairs, p+ (p) is the probability of a realization of wr being positive (negative), and Pθ+ is a conditional PDF such that Pθ+Pθ|wr0Pθ|θ[0,π2], θm is the lower bound of θ described above. For a first-order account, we neglect skewness in the distribution of particle relative velocities and set p±=0.5. Following Kolmogorov (1941), we have set wr|wr<0=ξr, where ξ±=Csε/(15ν), (Cs is a Kolmogorov constant, we found Cs=0.76 by matching ξr to the first-order fluid velocity structure-function from the DNS).

A simple phenomenological model for P(θ) may be constructed using the (statistical) central-limit theorem by assuming that the angle of approach θ at any time is the sum of many random-incremental rotations in the past, thus we write

(12) P ( θ ) = N exp [ K cos ( θ - μ θ ) ] sin ( θ ) ,

where Nexp […] is the circular normal distribution, i.e., analog of Gaussian distribution for angular data; sin (θ) results from integration over azimuthal angles (ϕ). We set μθ=π2 (neglect skewness in fluid's relative velocity PDF) and obtain K by matching the transverse to longitudinal ratio of structure functions (TLR) of the particle relative velocities with the ones via the DNS data; N is determined via normalization of P(θ). Figure 2b shows the wr derived via Eqs. (11) and (12), using K calibrated with TLR of second-, fourth-, and sixth-order structure functions respectively. The results have a correct qualitative trend of vanishing values at large r that increase sharply as r approaches d, with the fourth-order's result giving the best agreement with DNS. Currently we do not have a satisfactory rationale to single out the fourth order. The TLR of different orders giving differing results may imply that our first-order model may be incomplete, possibly due to over-simplification in Eq. (12) or to the inaccuracy of the rectilinear assumption (d/η in the DNS may be insufficiently small).

4.4 Differential version of the theory, its validity and its solution

We now discuss an important but precarious theoretical issue. Chun et al. (2005) clearly showed that the non-local diffusion (qrD) may only be converted from its general integral form into a differential version when the underlying RDF is a simple power law. However, Lu et al. (2010) and Yavuz et al. (2018), working in two very different scenarios, found that their predictions using the differential form of the theory agree well with experiments, even when the RDFs involved were clearly not power laws. We shall attempt to remedy this apparent paradox in future work. To examine how well this albeit unjustified method works here, we recast Eq. (9) into its differential form (Chun et al.2005):

(13) - τ η - 1 B n l r 4 g r + g ( r ) [ r 2 W r - A τ r 3 ] = - R c * ,

where Bnl=0.0397 (this value is computed from our DNS), and Bnl is expected to depend on flow characteristics, e.g., Rλ and τη (more in the Supplement). Using Eq. (13), the same gsg0 ansatz, we make another prediction for WrSt=0.054, which is plotted in Fig. 2a (dashed green line). This prediction is far from the DNS at rd but performs as well as the integral version at rd (the jump in the curve is just an artifact from the kink in the ansatz).

One advantage of Eq. (13) is that it allows for a general solution, which we now give, assuming Wr is given by Eqs. (11) and (12):

(14) g ( r ) = 1 β ( r ) β ( r ) q ( r ) d r + C ,

with q(r)=Rc*τη/(Bnlr4), β(r)=expp(r)dr and p(r)=Aτr-wrτη/(Bnlr2) (more in the Supplement).

4.5 Effects of gravity and other limitations

Thus far, we have not considered the effects of gravity on the particles. Here we provide a glimpse into the role of gravity (a detailed analysis is beyond the scope of the present study). In keeping with the scope of current work, we restrict ourselves to the case of monodisperse particles only. For this, we rerun the DNS cases of St=0.054 and 0.54 with gravity (to be compared with the zero-gravity case). The new particle advection equation is dv/dt=(u-v)/τp+g (all other details of the DNS remain unchanged). We choose to have the particle settling parameter Sgτpg/uη (where uη is the Kolmogorov velocity scale) be in the range O(0.1)−O(1) (this is achieved by letting |g|=10dm s−2). As a result, the range of Sg and St explored here is well aligned with measured values in natural clouds (Siebert et al.2010). For the case of St=0.054 (Sg=0.49), we find no discernible difference for both RDF and MRV between the “with gravity” and zero-gravity results (corresponding figures in the Supplement). For the St=0.54 (Sg=4.9) case, Fig. 3 shows the effects of gravity on the RDF and MRV. We see that the slope (exponent of g(r) in the range dr<20η) of the RDF in the gravitational case is reduced by about 15 % compared to the zero-gravity case (Sg=0). However, the shape of the RDF in the collision regime (rd) is approximately preserved, suggesting that a construct of the form gcollision×ggravity may be a good first-order model for the full RDF (close examination of the compensated RDFs gives substantial support for this idea, details in the Supplement). These observations imply that as Sg increases from O(0.1) to O(1), the effects of gravity on RDF grow from negligible to significant but not dominant; the main effect is the reduction of the exponent while the collision-related “modulation” (gcollision) remains largely intact. The inset of Fig. 3 shows that the MRV is also weakened by gravity, though the statistical noise limits the strength of this conclusion. Lastly, it is worth noting that in the complementary DNS by Woittiez et al. (2009) that included gravity but not actual collisions, a much stronger gravitational effect was found on the statistics of bidisperse particles relative to the monodisperse case.

Figure 3 RDFs of particles (St=0.54) subject to action of turbulence and collision–coagulation with and without gravity. Circles: Sg=0 (zero gravity); triangles: Sg=4.9 (nonzero gravity). The latter shows a reduced slope in the power-law regime, while the shape of the two curves is largely similar in the collision regime (rd). Inset: MRVs of the same cases as in the main figure. Gravity weakens the MRV of the particles.


As mentioned, the fundamental focus of our work precludes the DNS and theory from considering a number of complexities relevant to some applications. As a result, this limits the direct quantitative applicability of our results to some realistic problems (e.g., in clouds). Besides gravity, another neglected factor is the hydrodynamic inter-particle force (HDI). Recent works, e.g., Yavuz et al. (2018); Bragg et al. (2022) found that the HDI also has a strong impact on RDF for rd. For monodisperse particles with small to moderate St, the HDI is expected to be more important than gravity. While we expect that the HDI should not alter the qualitative trend that g(r) should fall towards a small value at rd (the same applies to the observed trend of MRV), it is likely that the HDI and collision would affect RDF and MRV in a coupled manner.

Also neglected is the influence of temperature, humidity and the vapor–liquid phase transition, which are important in the atmospheric clouds. These factors have a substantial impact on the polydispersity of small droplets (see, e.g., Kumar et al.2012, 2014). However, for monodisperse statistics considered here, they are likely to play a minor role (they will be more important when future works consider the full polydisperse problem).

One limitation of the theory stems from the assumption of St≪1 and its corollary that particle velocity statistics in this regime are St-independent (Chun et al.2005), which limits the theory's applicability to real systems. This implies that MRV should be St-independent in this regime. Our DNS results (spanning 2 orders of magnitude in St) shown in Fig. 2 give some support to the latter. However, unlike the theoretical prediction for MRV of case St=0.054 (Fig. 2), we have found that the prediction for St=0.11 is discernibly below the DNS result (figure in the Supplement). This could be due to the finite St effect not captured by the theory or could be for other reasons (details in the Supplement). Hence, a finite St extension of the theory is desirable to improve its applicability to real systems.

5 Conclusions

To conclude, we observed that collision strongly affects the RDF and MRV and imposes strong coupling between them6. This challenges the efficacy of a “separation paradigm” and suggests that results from any studies that preclude particle collision have limited relevance for predicting collision statistics. We have presented a theory for particle collision–coagulation in turbulence (based on a Fokker–Planck framework) that explains the above observations and verified its accuracy by showing that Wr could be accurately predicted using a sufficiently accurate RDF. The theory accounts for the full collision–coagulation rate, which includes contributions from mean field and fluctuations, and as such, our work complements and completes earlier mean-field theories (Saffman and Turner1956; Sundaram and Collins1997). We showed that a simple model of particle approach angles could capture the main features of Wr and used it to derive a general solution for RDF from the differential version of the theory. We uncovered a possible paradox regarding the past empirical successes of the differential drift-diffusion equation (see Sect. 4.4). Further shape-preserving reduction of the RDF and MRV was observed when the gravitational settling parameter (Sg) is of order O(1). Our findings provide new perspectives of particle collision and its relation with clustering and relative motion, which have implications for atmospheric clouds or generally for systems involving colliding particles in unsteady flows.

Code and data availability

The data (model codes) may be obtained from the corresponding author upon request.


The supplement related to this article is available online at:

Author contributions

EWS oversaw the conception and execution of the project. EWS did the theoretical derivations in collaboration with XM. XM and EWS conducted the numerical simulation and data analysis. EWS and XM wrote the article.

Competing interests

The contact author has declared that neither they nor their co-author has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


This work was supported by the National Natural Science Foundation of China (grant 11872382) and by the Thousand Young Talents Program of China. We thank Jialei Song for their help. We thank Wai Chi Cheng, Jianhua Lv, Liubin Pan and Raymond A. Shaw for discussion and suggestions.

Financial support

This research has been supported by the National Natural Science Foundation of China (grant no. 11872382).

Review statement

This paper was edited by Graham Feingold and reviewed by two anonymous referees.


Balkovsky, E., Falkovich, G., and Fouxon, A.: Intermittent Distribution of Inertial Particles in Turbulent Flows, Phys. Rev. Lett., 86, 2790, 2001. a

Bec, J., Biferale, L., Cencini, M., Lanotte, A., Musacchio, S., and Toschi, F.: Heavy Particle Concentration in Turbulence at Dissipative and Inertial Scales, Phys. Rev. Lett., 98, 084502, 2007. a, b

Bec, J., Ray, S. S., Saw, E. W., and Homann, H.: Abrupt growth of large aggregates by correlated coalescences in turbulent flow, Phys. Rev. E, 93, 031102, 2016. a

Bragg, A. D., Hammond, A. L., Dhariwal, R., and Meng, H.: Hydrodynamic interactions and extreme particle clustering in turbulence, J. Fluid Mech., 933, A31,, 2022. a

Chun, J., Koch, D. L., Rani, S. L., Ahluwalia, A., and Collins, L. R.: Clustering of aerosol particles in isotropic turbulence, J. Fluid Mech., 536, 219–251, 2005. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x

Dou, Z., Bragg, A. D., Hammond, A. L., Liang, Z., Collins, L. R., and Meng, H.: Effects of Reynolds number and Stokes number on particle-pair relative velocity in isotropic turbulence: a systematic experimental study, J. Fluid Mech., 839, 271–292, 2018. a

Eswaran, V. and Pope, S. B.: An examination of forcing in direct numerical simulations of turbulence, Comput. Fluids, 16, 257–278, 1988. a

Falkovich, G., Fouxon, A., and Stepanov, M. G.: Acceleration of rain initiation by cloud turbulence, Nature, 419, 151, 2002. a

Grabowski, W. W. and Wang, L.-P.: Growth of cloud droplets in a turbulent environment, Annu. Rev. Fluid Mech., 45, 293–324, 2013. a

Ireland, P. J., Vaithianathan, T., Sukheswalla, P. S., Ray, B., and Collins, L. R.: Highly parallel particle-laden flow solver for turbulence research, Comput. Fluids, 76, 170–177, 2013. a

Johansen, A., Oishi, J. S., Mac Low, M.-M., Klahr, H., Henning, T., and Youdin, A.: Rapid planetesimal formation in turbulent circumstellar disks, Nature, 448, 1022–1025, 2007. a

Karnik, A. U. and Shrimpton, J. S.: Mitigation of preferential concentration of small inertial particles in stationary isotropic turbulence using electrical and gravitational body forces, Phys. Fluids, 24, 073301, 2012. a

Karpińska, K., Bodenschatz, J. F. E., Malinowski, S. P., Nowak, J. L., Risius, S., Schmeissner, T., Shaw, R. A., Siebert, H., Xi, H., Xu, H., and Bodenschatz, E.: Turbulence-induced cloud voids: observation and interpretation, Atmos. Chem. Phys., 19, 4991–5003,, 2019. a

Kolmogorov, A. N.: The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers, Dokl. Akad. Nauk SSSR+, 30, 299–303, 1941. a, b

Kumar, B., Janetzko, F., Schumacher, J., and Shaw, R. A.: Extreme responses of a coupled scalar–particle system during turbulent mixing, New J. Phys., 14, 115020, 2012. a

Kumar, B., Schumacher, J., and Shaw, R. A.: Lagrangian mixing dynamics at the cloudy–clear air interface, J. Atmos. Sci., 71, 2564–2580, 2014. a

Lu, J., Nordsiek, H., Saw, E. W., and Shaw, R. A.: Clustering of Charged Inertial Particles in Turbulence, Phys. Rev. Lett., 104, 184505, 2010. a, b

Maxey, M. R. and Riley, J. J.: Equation of motion for a small rigid sphere in a nonuniform flow, Phys. Fluids, 26, 883–889, 1983. a

Mortensen, M. and Langtangen, H. P.: High performance Python for direct numerical simulations of turbulent flows, Comput. Phys. Commun., 203, 53–65, 2016. a, b

Onishi, R. and Seifert, A.: Reynolds-number dependence of turbulence enhancement on collision growth, Atmos. Chem. Phys., 16, 12441–12455,, 2016. a

Pope, S. B.: Turbulent Flows, Cambridge Univ. Press, Cambridge, UK,, 2000. a, b, c, d, e

Reade, W. C. and Collins, L. R.: Effect of preferential concentration on turbulent collision rates, Phys. Fluids, 12, 2530,, 2000.  a

Rogallo, R. S.: Numerical experiments in homogeneous turbulence, vol. 81315, National Aeronautics and Space Administration, 1981. a, b, c

Saffman, P. and Turner, J.: On the collision of drops in turbulent clouds, J. Fluid Mech., 1, 16–30, 1956. a, b, c, d

Salazar, J. P. L. C., de Jong, J., Cao, L., Woodward, S. H., Meng, H., and Collins, L. R.: Experimental and numerical investigation of inertial particle clustering in isotropic turbulence, J. Fluid Mech., 600, 245–256, 2008. a, b

Saw, E. W., Shaw, R. A., Ayyalasomayajula, S., Chuang, P. Y., and Gylfason, A.: Inertial Clustering of Particles in High-Reynolds-Number Turbulence, Phys. Rev. Lett., 100, 214501, 2008. a, b, c

Saw, E.-W., Shaw, R. A., Salazar, J. P., and Collins, L. R.: Spatial clustering of polydisperse inertial particles in turbulence: II. Comparing simulation with experiment, New J. Phys., 14, 105031, 2012. a, b

Saw, E.-W., Bewley, G. P., Bodenschatz, E., Sankar Ray, S., and Bec, J.: Extreme fluctuations of the relative velocities between droplets in turbulent airflow, Phys. Fluids, 26, 111702, 2014. a

Siebert, H., Gerashchenko, S., Gylfason, A., Lehmann, K., Collins, L., Shaw, R., and Warhaft, Z.: Towards understanding the role of turbulence on droplets in clouds: in situ and laboratory measurements, Atmos. Res., 97, 426–437, 2010. a

Sundaram, S. and Collins, L.: Collision statistics in an isotropic particle-laden turbulent suspension. Part 1. Direct numerical simulations, J. Fluid Mech., 335, 75–109, 1997. a, b, c, d, e, f, g

Voßkuhle, M., Lévêque, E., Wilkinson, M., and Pumir, A.: Multiple collisions in turbulent flows, Phys. Rev. E, 88, 063008, 2013. a, b

Wang, L.-P., Ayala, O., Rosa, B., and Grabowski, W. W.: Turbulent collision efficiency of heavy particles relevant to cloud droplets, New J. Phys., 10, 075013,, 2008. a, b

Wilkinson, M., Mehlig, B., and Bezuglyy, V.: Caustic activation of rain showers, Phys. Rev. Lett., 97, 48501, 2006. a

Woittiez, E. J., Jonker, H. J., and Portela, L. M.: On the combined effects of turbulence and gravity on droplet collisions in clouds: a numerical study, J. Atmos. Sci., 66, 1926–1943, 2009. a, b

Wood, A. M., Hwang, W., and Eaton, J. K.: Preferential concentration of particles in homogeneous and isotropic turbulence, Int. J. Multiphas. Flow, 31, 1220, 2005. a

Yavuz, M., Kunnen, R., Van Heijst, G., and Clercx, H.: Extreme small-scale clustering of droplets in turbulence driven by hydrodynamic interactions, Phys. Rev. Lett., 120, 244504, 2018. a, b, c

Zaichik, L. I. and Alipchenkov, V. M.: Pair dispersion and preferential concentration of particles in isotropic turbulence, Phys. Fluids, 15, 1776, 2003. a


Coagulation is, in a sense, the simplest outcome of collision. In the companion paper we shall argue that the major qualitative conclusions of our work also apply to cases with other collisional outcomes.


In the collision-less case, the asymmetry is much weaker and is related to viscous dissipation of energy in turbulence (Pope2000).


In other words, particles may approach each other (and collide), but they can not be created at contact and then separate.


Here a positive qrd merely reflects a deficit in the inward flux of neighboring particles since we find that qrd+qrD is always negative.


This is true in the relatively idealized system simulated but may not apply to the general problem that includes other effects


This statement also holds for other types of collisional outcomes (not only for collision–coagulation), but the details of the specific outcomes should be different from the current case.

Short summary
Collision–coagulation of small droplets in turbulent clouds leads to the production of rain. Turbulence causes droplet clustering and higher relative droplet velocities, and these should enhance the collision–coagulation rate. We find, surprisingly, that collision–coagulation starkly diminishes clustering and strongly alters relative velocities. We provide a theory that explains this result. Our results call for a new perspective on how we understand particle/droplet collision in clouds.
Final-revised paper