Intricate Relations Among Particle Collision, Relative Motion and Clustering in Turbulent Clouds: Computational Observation and Theory

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 $r \sim d\,\,$ (where $d$ is the particle diameter) to small but finite values, while the mean radial-component of particle relative velocities (MRV) increase sharply in magnitudes. 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 ($S_g$) 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.

S1 Further Details of the Direct Numerical Simulation.
The time step in our DNS t is 0.001 s.The courant number is C = 0.073, (where C = t u x + v y + w z max , u etc. are r.m.s.velocities, ∆x etc. are grid spacings).The normalized maximum wavenumber simulated is k max η = 1.2.The turbulent flow is sustained by randomly forcing the two lowest nonzero shells of wave numbers.The integral length scale of the turbulent flow is estimated to be L = 0.646 dm.
We study the statistics of monomers only (i.e. the particle of the same size (d) that we initially introduce into the system and which we later replenish at a constant rate close to the monomer-monomer collision rate).In this sense, the particle (monomers) are naturally lost from our consideration once they collide and become larger particles.Particles that become much larger (St > 21.6) are removed from the DNS at each time step.
S2 Estimation of Leading Order Terms in the Drift Flux, e.g a (1) ik Using the DNS data, we estimate e.g. the value of t −∞ a (1) Note: the averaging is done over fluid particles (the theory assumed St 1 limit, such that all velocity statistics are tied to the fluid's), the integrand is non-vanishing only for t in the vicinity of t − τ η to t (where the turbulent velocity gradient Γ ij retains correlation), thus this quantity may be approximated as: τ η 2 Γ ik (t)Γ lm (t)Γ ml (t) .As shown in Chun et al. (2005), Γ ik (t)Γ lm (t)Γ ml (t) is by definition zero in fully developed turbulence due to the fact that the small-scale statistics of turbulent flows are almost isotropic Kolmogorov (1941).However, the coagulation constraint dictates that at r = d, such averages must be taken with the condition that only fluid-particle pairs with negative radial velocity (w r < 0) are taken into 1 Supplementary Material account (that the inertial particles' motion being tied to the fluid's does not imply that inertial pairs sample the fluid particle pairs's motion uniformly).Under this condition, the DNS data gives τ η 2 Γ ik (t)Γ lm (t)Γ ml (t) ≈ (−0.171 ×10 −3 dm/s)/d * , (d * = 9.49 ×10 −4 dm); here, it is of value to point out that without such constraint or condition, the result for this quantity from the DNS is two orders of magnitude smaller.Similarly, we found for this quantity, the DNS gives roughly the same values with or without the constraint.

S3
Full Definition of the Function f I (R 0 , µ, t f ) in the Model for Non-local Diffusive Flux.Chun et al. (2005), summarized here (with typo corrected), the diffusive action of the turbulence on the particle-pairs is assumed to consist of a random sequence of uniaxial extensional or compressional flows defined, and:

Derived in
where R 0 ≡ r 0 /r, r 0 is the initial separation distance of a particle pair before a straining event, r is the independent variable of the equation for g(r); f + and f − ≡ 1 − f + are the fractions of those flows that are extensional and compressional, respectively.
Comparing with DNS, Chun et al. (2005) calibrated f + and found f + = 0.188 (a result we use here).I ± is an indicator function such that it takes the value +1 (−1) when a secondary particle leaves (enters) a sphere of radius r centered on the primary particle, and otherwise zero.µ is the cosine of the angle between the axis of symmetry of the straining flow event and the displacement vector between the two particles, t f is the lifetime of the event.To obtain a strain rate correlation function that decays exponentially with a characteristic time scale τ S , Chun et al. (2005) set the probability density function for t f to be: The indicator function is used to count the net loss of particles from within the sphere over the duration of an (extensional or compressional) event and can be expressed as: where H(x) is the Heaviside function (zero for x < 0, unity for x ≥ 0), R f ± is the non-dimensional final position of a particle pair with an initial position of R 0 and can be written as: , for uniaxial extension and compression respectively, where: In this work, we deviate crucially 1 from the CK theory Chun et al. (2005) by introducing an extra factor c st (positive, of order unity or less) in the model of non-local diffusion: To determine what c st is (or should be), we begin from an important finding in Chun et al. (2005) that if P is power-law of r, i.e.P = Cr −c1 , then the non-local diffusion q D r can be cast into a differential form (which is usually only true for local 1 'Crucial' refers to the fact that without cst the theory would be inconsistent with previous experimental results (as this section will show) and it would also produces results far from our DNS results.diffusion): where: This, together with: .We now begin from this experimentally validated result and work backward to derive an expression for c st .We plug the power-law form for P into (S2): Comparing with (S1), we have: which is found in experiments (and theories) to be of order 0 to 1 and a function of particle Stokes number St; in words, this means c st is given by the modulus of the power-law exponent of the RDF that would arise in the collision-less case; in the case with collision and sufficiently small particle (d/η 1 ), such as in this study, c st equals the modulus of the power-law exponent of the RDF the range of d r 20η (note: power-laws RDF are empirically observed for r 20η Saw et al. (2008Saw et al. ( , 2012a))).
Note: we have chosen to define c st using the 'modulus' (instead of the 'negative' of the power-law exponent) since it guarantees that q D r is negative (positive) when g(r) is an increasing (decreasing) function of r, so that we are consistent with the fact that q D r is a diffusion flux.We note that both the CK theory and the current modified version assume St 1.
where we have clarified that A τ in our work is defined differently from "A" in (Chun et al., 2005) (we denote the latter as A ck to avoid confusion), and A τ,r d is our A τ evaluated at the large-r limit.In the current context, c 1 maybe obtained via (S4) or alternatively directly from the power-law exponent of g(r) in the range d r 20η as discussed above.Using values of the relevant parameters in our DNS, we found Aτ τη B nl ≈ 2.4St 2 ×.0925 .0397 = 5.6St 2 , which is 15% smaller than the one found in Chun et al. (2005), i.e.A ck B nl ≈ .61St 2 .0926= 6.6St 2 .However, we have observed in our DNS that the direct method (by fitting power-laws to the RDFs in the suitable r-range) gives c 1 which is 3.2 (1.9) times larger than the one obtained using (S4) for the case of St = 0.054 (0.11).
A plausible interpretation of the discrepancy described just above is that there may be another missing dimensionless factor (of order unity, possibly weakly dependent on Reynolds-number) in the correct definition of q D r .This is beyond the scope of this present study (to avoid confusion, we currently restrict ourselves to the least speculative correction only) and is a good subject for future works.However it may be informative to note that, by inspection, we find that if we further include a factor of ∼ 1/3 to 1/2 in the definition of q D r , then the agreement between the theoretical (the integral version) and DNS produced W r is strikingly better in the r d limit, while in the r ∼ d regime, it is slightly better (the former should not come as a surprise as this is the regime of power-law RDFs and the factor of ∼ 1/3 is exactly designed to reproduce the correct c 1 ).
To demonstrate the point just discussed, we show in Fig. 1 the predictions by the theory for the case of St = 0.11.We see that, in the r ∼ d regime, the prediction by the original theory (dotted line in the main figure) is somewhat below the DNS result, while the prediction by the modified theory (with a factor of 1/2 appended to the definition of q D r ), shown as the solid line, is much closer to DNS.In the r d regime, the modified theory's superiority in terms of accuracy is even more pronounced (see inset of Fig. 1).
S6 Relation Between g(r) and P .
In the main text, we state that g(r) ≡ V P , where V is the spatial volume of the full domain of the problem i.e. (2π) 3 in the DNS.Justification: let g(r) be the ratio of probability of finding a second particle at r from a particle, to the probability of such finding in a perfectly random distributed particle population, thus: g(r) ≡ P δxδyδz (δxδyδz)/V ≡ P V .Further, since system is isotropic, g(r) ≡ g(r) .

S7 Modeling of MRV based on Distribution of Particle Approach Angles P (θ).
We imagine the particles are small i.e. d η and St 1.The latter implies their trajectories are almost like fluid particles', while the former implies that, viewed at the scale of interest r ∼ d, their trajectories are almost rectilinear (since the radii of curvature are proportional to η).Thus in the reference frame of a primary particles, no secondary particle could have a trajectory, being straight-line, that has a history of collision with the volume of the primary (otherwise coagulation would have occurred and the secondary particle in question would cease to exist).In trigonometric terms, let θ be the angle between the secondary particle's velocity and its vector position in the rest frame of the primary particle, then we must have: From the above, we could then compute the MRV, w r * based on fluid particles' statistics.Since collision-coagulation affects positive and negative relative particle velocities differently, we begin by writing w r * as a sum of the positive (i.e.w r > 0) and negative branches (with proper statistical weights p ± to account for possible skewness of the probability distribution of velocity): The negative branch p − w r |w r < 0 * is unaffected by collision-coagulation and we thus express it as a simple linear function of r that follows from the K41-phenomenology (Kolmogorov, 1941), i.e. −p − ξ − r, where ξ ± ∼ ε/(15ν), ε is the (kinetic) energy dissipation rate of the flow.For the positive branch, we further assume that the (fluid particles') joint probability density function (PDF) of magnitude of relative velocity (secondary particle relative to primary particle) |w| and approach-angle θ, P (|w|, θ), is separable (note: w r ≡ |w| cos(θ)), hence: , where in the last line, we have replaced the first two integrals, combined, with its K41 estimate, where ξ ± ∼ ε/(15ν) .
S8 Prediction of the Peak Location of the RDF Using the Differential Form of the Drift-Diffusion Equation.
A finite R * c inhibit us from locating the peak of the RDF using (S5) à la Lu et al. (2010) i.e. without knowing g(r), since g(r) could no longer be factored out when ∂g ∂r = 0.However, we argue that (S5) could still give a reasonably accurate account of the peak location.For the case of St = 0.05, at r = 3d (the approximate peak location), we found the DNS data gives −τ η B nl r 4 ∂g ∂r ≈0 + g(r) r 2 W r − A τ r 3 ≈ −1.05×10 −9 and −R * c ≈ −1.01×10 −9 S9 General Analytical Solution for the Differential Form of the Drift-Diffusion Equation.
The general solution for the first-order non-homogenous ordinary differential equation (see e.g.Arfken and Weber (1999)), with w r * given by the model in the main text, is: with q(r) = R * c /(τ η B nl r 4 ); β(r) = exp p(r)dr and p(r) = [ A τ r − w r * ] /(τ η B nl r 2 ).For the current model described in the main text, the integral in (S6) could not be expressed in terms of simpler canonical functions.Hence, for specific applications, we currently anticipate that some sort of power-law expansion or asymptotic reduction (if not numerical integration) would be needed to produce problem specific analytical approximations.We repeat the DNS case of St = 0.054 and St = 0.54 with the particles subjected to gravity (body force), and compares results with the zero-gravity case.Fig. 2 shows the results for case St = 0.054.There is no discernible difference between the cases with and without gravity. 160 For case St = 0.54, the main RDF and MRV results is shown the main text.Here we show only the compensated-RDFs ( g c (r) ), where each g c (r) is calculated via g(r) divide by a power law (c 0 r −c1 ) that resulted from curve-fitting to the original g(r) in the range 0.6η ≤ r ≤ 3η.Fig. 3 compares g c (r) for cases with and without gravity.The fact that there is no discernible difference implies that the uncompensated g(r) could be model as g c ×g g where g c is function that depends only on the particle collision process while g g depends on other factors e.g.gravity and is independent of particle collision.Circles: Sg = 0 (zero gravity); triangles: Sg = 4.9 (nonzero gravity).Interpretation in the text.
by the Theory for Other Stokes Numbers.As mentioned in the main text of this manuscript, even though the theory assumes that MRVs are St-independent for small St's.Nevertheless, it could produce separate predictions for each St.Here we show the predictions for St = 0.001 and 0.11 in Fig.1.The prediction for St = 0.001 (red dash line) agrees well with the DNS results (symbols), but the prediction for St = 0.11 (gold dotted line) deviates significantly, suggesting that finite St effects not captured by the theory start to become significant and thus diminish the accuracy of the theory.S5 Derivation of c st , its Role and Possibility of Further Corrections to The CK Theory.

3Figure S1 .
Figure S1.MRVs of particles.Triangles: DNS result for St = 0.001.Squares: DNS result for St = 0.11.Red-dashed-line: theory's prediction for St = 0.001.Gold-dotted-line: theory's prediction for St = 0.11.Grey-solid-line: modified-theory's prediction for St = 0.11 (details of modification in Sec.S5).Note: the predictions are based on the Aτ values valid in the r ∼ d regime (d = 9.49×10 −4 dm).Inset) Similar plots in logarithmic axes highlighting the large-r regime.Note: the predictions are based on the Aτ values valid in the r d regime.It is clear that the modified-theory's predictions agree much better with the DNS.