Statistics of relative velocity for particles settling under gravity in a turbulent flow

We study the joint probability distributions of separation, R, and radial component of the relative velocity, VR, of particles settling under gravity in a turbulent flow. We also obtain the moments of these distributions and analyze their anisotropy using spherical harmonics. We find that the qualitative nature of the joint distributions remains the same as no gravity case. Distributions of VR for fixed values of R show a power-law dependence on VR for a range of VR, exponent of the power-law depends on the gravity. Effects of gravity are also manifested in the following ways: (a) moments of the distributions are anisotropic; degree of anisotropy depends on particle’s Stokes number, but does not depend on R for small values of R. (b) mean velocity of collision between two particles is decreased for particles having equal Stokes numbers but increased for particles having different Stokes numbers. For the later, collision velocity is set by the difference in their settling velocities.


I. INTRODUCTION
Small-sized heavy particles suspended in a turbulent flow are found in many natural phenomena. Some of the examples are dust storms, water droplets in clouds, and astrophysical dust in the protoplanetary disk around a young star. Collisions of these particles play an important role in many processes. For instance, in clouds, small water droplets collide with each other and may coalesce to form larger droplets. Understanding this process of collision and coalescence is important to understand the process of initiation of rain from warm clouds [1][2][3]. The motion of these small droplets is determined by the following two forces acting on them: (1) a hydrodynamic drag force applied by the surrounding fluid; (2) an external force such as gravity in the case of cloud droplets. If radius a of the droplet is very small compared to the Kolmogorov dissipation scale η of the flow and its density ρ p is much larger than the density ρ of underlying fluid, then the equation of motion of such a droplet can be written as: Here x and v denote the position and velocity of a particle, respectively. τ p = (2ρ p /9ρ)a 2 /ν is the characteristic response time of the particle, here ν is the kinematic viscosity of the fluid. g is the acceleration due to gravity andẑ is the unit vector along vertically downward direction. u(x, t) is the flow velocity at the position of particle. This model assumes that fluid-inertia corrections are small and neglects the hydrodynamic interaction between the particles. This model Eq. (1) also neglects the history term present in the Maxey-Riley equation [4]. Recent studies of Refs. [5,6] shows that history term plays an important role in the dynamics of heavy particles and its effects depend on the density ratio ρ p /ρ. This ratio is 1000 for water droplets in the air, in this case, the effects * akshayphy@gmail.com of history term are found to be very small [6]. The inertia of the particle is measured in terms of Stokes number St = τ /τ η , where τ η is the Kolmogorov time scale of the turbulent flow. To measure the relative importance of turbulence and gravity, we define Froude number Fr = uη τηg , where u η = η/τ η , zero gravity corresponds to Fr = ∞. For cloud droplets, St and Fr depend on the mean rate of energy dissipation per unit volume ε. Value of ε varies a lot in clouds [7] and hence St and Fr can have different values in different clouds. Typically for 10 to 50 micrometer-sized droplets St roughly lies between 0.1 to 2.5 and Fr may vary from 0.05 to 0.3 [3,8]. Note that, St and Fr appear by the nondimensionalization of Eq. (1) by using η, τ η , and u η . Many studies use another dimensionless number called Rouse number Ro = u rms /(τ p g) [9,10], where u rms is the root-meansquare speed of the flow. For a given value of Reynolds number, Ro can be expressed in terms of St and Fr. In the present work, we keep the Reynolds number of the flow fixed and vary St and Fr independently.
It is known that small particles cluster in turbulence due to their inertia. Small scale clustering is measured in terms of the probability P (R) of finding two particles within a separation R. It is found that P (R) ∼ R d2 as R → 0 [11], d 2 is called the spatial correlation dimension. If d 2 is smaller than the spatial dimension d, more particles are found at small separation compared to uniformly distributed particles. Inertial particles can detach from the flow streamlines and nearby particles can have high relative velocity, this is referred as sling effect [12]. This can be understood in terms of singularities in the velocity field of the particles [13,14]. These singularities are called caustics. Recently, it has been shown theoretically for smooth-random flows [15,16] and in the direct numerical simulations (DNSs) of the turbulent flows [17,18] that the probability distribution functions (PDFs) of the radial component of relative velocity V R at small separation, have a power-law tail. The exponent of this powerlaw is related to the phase space correlation dimension D 2 . It is also shown that there exists a parameter z * such that the joint PDFs of separation R and V R are independent of V R for V R ≪ z * R and are independent of R for V R ≫ z * R [18]. Here, R and V R are nondimensionalized by using η and u η , respectively. Parameter z * sets a velocity scale for a given separation R and is often referred as the matching scale [16]. Refs. [15,16] gives the following theoretical expression for the joint PDF P(R, V R ), for R < 1, obtained for smooth random flows in the white noise limit: Form of the joint distribution P(R, |V R |) agrees well with the theoretical prediction Eq. (2) for particles suspended in a turbulent flow [18]. Mean collision velocities of particles in a polydisperse suspension are studied in [19]. It is found that the collision velocities between different sized particles can be much higher compared to collision velocities between equal-sized particles.
Many studies have focused on clustering and relative velocities of the particles in turbulent flows without gravity. Relatively less is known about the effect of gravity on clustering and statistics of the relative velocities of the particles. For instance, it is not known, how gravity affects the form of the joint distribution P(R, |V R |) given in Eq. (2)? How anisotropy of the system changes as a function of R and St? How mean collision velocities of the particles are modified by the gravity? This paper focuses on these questions.
Gravity and turbulence together can give rise to many non-trivial phenomena in the dynamics of the heavy particles. In the presence of turbulence, particles are found to settle with a speed that is higher compared to their terminal-speed in still fluid [20][21][22]. Reduction of settling speed by turbulence is also observed in the experiments [10]. Simulations with non-linear drag force on the particles show a similar reduction in the settling speeds. Present work considers only the linear drag case (for which mostly enhancement is observed in simulations), and does not address the issue of enhancement vs hindering of settling by turbulence.
Refs. [22][23][24] shows that the presence of gravity modifies the small scale clustering of heavy particles. This modification depends on the values of St and Fr. For particles having St < 1 clustering is reduced compared to no-gravity case whereas for particles having St > 1 clustering is significantly increased compared to no-gravity case. In a real experiment, gravity is always present hence, these results can be used to understand the experimental observation of particle clustering in turbulent flows [25][26][27].
It is argued that modification of clustering happens because settling particles sample the flow differently compared to particles advected solely by turbulence.
Mean relative velocity of the particles having equal St is found to be reduced by gravity [22,24]. Refs [24,28] study the PDFs of relative velocity at small R in the presence of gravity. These studies find that the PDFs remains non-Gaussian similar to the case when Fr = ∞ but, fluctuations are reduced due to the presence of the gravity. Mean collision rates of equal St particles are found to be reduced by gravity [8,29,30]. This is due to the significant reduction of the mean relative velocity. Anisotropy of the clustering and mean relative velocity is analyzed in Ref. [24] by calculation the radial distribution function and mean relative velocity as a function of R and angle θ between separation vector and direction of gravity. Using spherical harmonic decomposition of these quantities it is shown that the anisotropy is significant at small separation for large values of St. This analysis of anisotropy is done for separations larger than η. Statistics of relative velocities in bi-disperse turbulent suspensions are studied in Ref [31]. It is found that gravity enhances the relative velocities in vertical and horizontal directions. It is shown that for small values of Fr relative velocity is dominated by differential settling but turbulence still plays an important role.
As radii of droplets are always smaller than the dissipation scale η, the separation between the droplets is also much smaller than η when they collide. Therefore, we study the relative velocities of particles for small values of separations. In this paper, we focus on the joint PDFs P(R, |V R |) and its moments for different values of St and Fr. We use length scale η and velocity scale u η to nondimensionalize separation R and relative velocity V R , respectively. For separation R < 1. We show the following results for the first time: • The joint PDF P(R, |V R |), is qualitatively similar to the case of zero gravity; in the sense that in both of these cases: -There exists a parameter z * such that at small R the joint PDF is independent of V R for V R < z * R, and is independent of R for V R > z * R.
-For a fixed R the PDF as a function of V R shows a power-law range with exponent D 2 − 4. Although the phase-space correlation dimension, D 2 , itself is not the same as the zero-gravity case.
• The spatial clustering, described by the 0-th moment of the joint PDF, is anisotropic, in a very special way. A decomposition of the moment into spherical harmonics show that different order harmonics, for ℓ = 0, 2, and 4, scales with the same exponent d 2 but the amplitudes depend on ℓ.
• Similar behaviour is seen for the first moment, which describes the mean relative velocities, too.
• We define collision velocities as the relative velocity of two particles separated by the sum of their radii. The mean and RMS of collision velocities as a function of St 1 and St 2 changes qualitatively from the zero-gravity case; in particular the contours in the St 1 − St 2 plane becomes parallel to the diagonal as Fr decreases, i.e., the collision velocities are a function of | St 1 − St 2 | alone.
• Furthermore, for particles with equal St, the mean collision velocity decreases as Fr decreases, i.e., gravity increases. This decrease is more pronounced at higher St than at lower ones.
• For St 1 = St 2 the qualitative behaviour is opposite, the mean collision velocity increases as Fr decreases. In this case, mean collision velocity is set by the difference in the settling velocities.

II. DIRECT NUMERICAL SIMULATION
The flow velocity u(x, t) is determined by solving the Navier-Stokes equation Here D Dt ≡ ∂ t + u · ∇ is the Lagrangian derivative, p is the pressure of the fluid, and ρ is its density as mentioned above. The dynamic viscosity is denoted by µ ≡ ρν, and S is the second-rank tensor with components Here ∂ k u j are the elements of the matrix A of fluid-velocity gradients. f denotes the external force. To relate pressure p and density ρ, we use the ideal gas equation of state with a constant speed of sound, c s .
To solve Eqs. (3) we use an open-source code called The Pencil Code [32]. This code has been used earlier for the similar studies of particles-laden turbulent flows in Refs [18] and [33]. It uses a sixth-order finite-difference scheme to compute space derivatives and a third-order Williamson-Runge-Kutta [34] scheme for time evolution. We use periodic boundary conditions in all three directions. A white-in-time, Gaussian, forcing f is used that is concentrated on a shell of wavenumber with radius k f in Fourier space [35]. Forcing term is integrated by using the Euler-Marayuma scheme [36]. A statistically stationary state is reached, where the average energy injection by the external forcing f is balanced by the average energy dissipation by viscous forces. The amplitude of the external forcing is chosen such that the Mach number, Ma ≡ u rms /c s is always less than 0.1, i.e., the flow is weakly compressible, where u rms = u 2 is the rootmean-squared velocity of the flow. This week compressibility has no important effect on our results; please see the discussion in Ref. [18], section II and Appendix A in Ref [18] for further details. The same code with similar setting has been used before to study the scaling and intermittency in fluid and magnetohydrodynamic turbulence [37][38][39]. Our simulations are performed in a threedimensional periodic box with sides L x = L y = 2π and L z = 4π. This box is discretised in N equally spaced grid points in each direction.  Due to gravity and periodic boundary conditions, settling particles can travel box length L z in a time interval which is less compared to the correlation time scale of the flow, T eddy in our case. This can produce errors in the statistics measured for the particles [24,28]. To minimize this, we use a longer box in the z direction which is the direction of gravity. From the values of the parameters in TABLE I, we estimate that having L z = 4π is good enough for the values of St considered in this paper.
We introduce the particles into the simulation after the flow has reached a statistically stationary state. Initially, particles are uniformly distributed in the simulation domain with zero initial velocities. To evolve positions and velocities of particles according to Eq. (1), we use a thirdorder Runge-Kutta scheme. To obtain the flow velocity at the positions of the particles, we use a tri-linear interpolation method.
Parameters of the simulations are given in TABLE I. ε ≡ 2νΩ is the mean rate of energy dissipation per unit volume, where the enstrophy Ω ≡ |ω| 2 , and ω ≡ ∇ × u is the vorticity. Taylor microscale Reynolds number is defined as Re λ ≡ (u rms λ)/ν, where u rms is the root-meansquare velocity of the flow averaged over the whole domain, λ ≡ (5u 2 rms )/(2Ω) is the Taylor microscale of the flow, and ν is the kinematic viscosity. The Kolmogorov length is defined as η ≡ (ν 3 /ε) 1/4 , the characteristic time scale of dissipation is given by τ η = (ν/ε) 1/4 and u η ≡ η/τ η is the characteristic velocity scale at the dissipation length scale. The large eddy turnover time is given by T eddy ≡ 1/(k f u rms ). In rest of this paper, unless otherwise stated, we use η, τ η , and u η to non-dimensionalize length, time, and velocity respectively.

III. RESULTS
We compute the joint PDFs P(R, |V R |) and its p−th moments defined as: Here, R and V R are nondimensionalized by scales η and u η . m 0 (R) measures the number of particle pairs having separation between R and R+dR. If there is no clustering m 0 (R) ∼ R 2 in three dimensions. We define V n as |V R | when two particles are moving towards each other (V R < 0) and separation R between their centre of masses is equal to a 1 + a 2 , where a 1 and a 2 are the radii of two particles nondimensionalized by η. Please note that the actual parameters that characterize the particles are St 1 and St 2 . We fix the density ratio ρ p /ρ = 10 3 (which is the case for the water droplets in the air) to get the radii a 1 and a 2 from St 1 and St 2 , respectively [40]. We call V n the collision velocity and consider it as a proxy for actual collision velocity. Mean of V n is defined as: A. Real space clustering Fig. (1)(a) shows plots of m 0 (R)/R 2 as a function of R for St = 3.02 and three different values of Fr. We observe that m 0 (R) has a power-law dependence on R with exponent equals to d 2 −1. For Fr = ∞ data shows no clustering as m 0 (R)/R 2 is independent of R. We observe that for the same St but Fr = 0.3, and 0.1 there is a significant amount of clustering. In Fig. (1)(b), we plot real space correlation dimension d 2 as a function of St for three different values of Fr. We find that for Fr = 0.1 and 0.3 and St > 1, d 2 is small compared to Fr = ∞. This implies that the settling enhances the clustering for St > 1 for the values of Fr considered here. This is consistent with the results of Refs. [22]. We also observe that for Fr = 0.1 and 0.3, d 2 becomes constant for large values of St. This indicates that the clustering does not depends on St for large values of St, this is consistent with the experimental observation of Refs [25,26]. For the range of St and Fr studied here, d 2 remains less than d, hence the phase space correlation dimension D 2 = d 2 .
B. Joint distribution of R and VR Fig. (2) (a) shows the joint distribution P(R, |V R |)/R 2 for St = 2.11 and Fr = 0.1. The dashed blue line shows |V R | = z * R, where z * ≃ 0.1 is found by fitting a line to the data. We find that the value of z * does not depends much on St, this is consistent with the results of [33]. We also observe that z * has a very week dependence on Fr but, we do not have many values of Fr to comment on this dependence. We observe that for |V R | < z * R, contour lines are vertical which implies that the distribution depends only on the separation R. For |V R | > z * R contour lines are horizontal, implying that it depends only on |V R | in this regime. Similar nature of the distributions P(R, |V R |) is observed for zero gravity (Fr = ∞) case in Refs. [18,33]. In Fig. (2) (b), we plot the distribution P(R, V R ) for R = 0.1, St = 2.4, and for three different values of Fr. We find that same as no-gravity case distributions have a power-law tail with scaling exponent D 2 −4. As D 2 changes with Fr see Fig. (1), values of scaling exponents also changes. Observe that as we decrease Fr, power-law tails of the distributions become steeper. This implies that the mean and rms values of |V R | also decrease as gravity increases.
To understand collisions of particles one would like to know the distribution of V R at R = 2a, where a is the radius of the particle. For a given value of St, the radius depends on the density ratio ρ p /ρ (as described at the beginning of Section III). For example, for St = 2.11 and ρ p /ρ = 10 3 , value of a/η is roughly 0.1, hence data for R < 0.1 in Fig. (2) (a) is irrelevant for this density ratio. For a different value of ρ p /ρ and the same value of St, one would need the distribution at a different value of R, therefor covering a range of R is useful to know the statistics of V R (R = 2a) for a range of ρ p /ρ.

C. Effect of anisotropy
Above calculations do not take in to account the anisotropy due to the presence of gravity. To compute the anisotropic contributions, we can obtain the joint distribution P(R, |V R |, cos(θ)), where θ is the angle between R and direction of gravityẑ. We compute the moments m p (R, cos (θ)) of this distribution and use spherical harmonics decomposition. A function f (R, θ, φ) that depends on spherical co-ordinates R, θ, and φ can be expressed as a linear combination of Spherical harmonics Y m ℓ (cos(θ), φ) as: As there is no dependence on the azimuthal angle φ in our system, only m = 0 modes are present. We can write: m p (R, cos(θ)) = ℓ m ℓ p (R)P ℓ (cos (θ)), here P ℓ (x) is the Legendre polynomial of order ℓ. As m p (R, cos(θ)) is a symmetric function of cos(θ), all coefficients m ℓ p (R) for odd values of ℓ vanish. m 0 (R, cos(θ)) is plotted in Fig. (3)(a) for Fr = 0.1 and four different values of St, for R = 0.5. m 0 0 (R), m 2 0 (R), and m 4 0 (R) are plotted in Fig. (3)(b) for St = 3.02 and Fr = 0.1. We find that m ℓ 0 (R) = C ℓ R ζ ℓ . Exponents ζ ℓ are independent of ℓ and are equal to d 2 whereas the coefficients C ℓ is zero for odd ℓ and decrease with ℓ for even ℓ.
To measure the degree of anisotropy we plot m 2 0 /m 0 0 as a function of R for different values of St in Fig. (3)(c). We observe that this ratio does not depend on R for the range of R shown here. This is because both m 2 0 and m 0 0 scales with the same power as a function of R. We also find that anisotropy increases with increasing St. In Fig. (4), we repeat this analysis for m 1 (R, cos(θ)). Fig. (4) (b) shows plots for m 0 1 and m 2 1 as a function of R for St = 3.03 and Fr = 0.1. Both of these coefficients show a power-law behaviour with roughly the same exponent. This is more clear from the plots of m 2 1 /m 0 1 as a function of R shown in Fig. (4) (c). It can be seen that the curves for higher values of St are not constant as a function of R, but dependence is very weak.

D. Mean collision velocity
The mean rate of collision depends on the relative velocity of particles when the separation between their centre of masses is equal to the sum of their radii. We define this as collision velocity V n (see Eq. (5) and text before it). In Fig. (5) we plot the mean of V n as a function of St 1 and St 2 , for Fr = 0.3, 0.1, and ∞. We observe that the qualitative nature of these plot changes due to gravity. Along the diagonal St 1 = St 2 , values for non zero gravity are smaller compared to no gravity case. Away from diagonal for different sized particles mean of V n becomes a function of |St 1 − St 2 | for the settling particles. We also notice that values for settling particles away from the diagonal are much higher compared to particles with Fr = ∞.
To analyze it more quantitatively, we plot V n for St 1 = St 2 = St in Fig. (6) (a), here denotes the mean. We observe that as Fr decrease V n decrease. In Fig. (6) (b) we plot V n as a function of St 2 for St 1 = 0.1. In this case, we find that V n increases significantly as Fr is decreased. Solid lines in this plot show the V n for particles settling under gravity in a fluid at rest.

IV. CONCLUSIONS
We have studied the effect of gravity on the relative velocities of particles suspended in turbulent flow by using direct numerical simulations. We explored a range of St and Fr that is relevant for small droplets in clouds. Our results show that joint PDFs of R and |V R | has a form that is qualitatively similar to the one obtained without gravity. It is shown that for separations smaller than the Kolmogorov scale η (R < 1 in nondimensional units), there exist a scale z * such that joint PDF is independent of |V R | for |V R | < z * R, and is independent of R for |V R | > z * R. We also show that PDFs of |V R | for a fixed value of R scales as |V R | D2−d−1 for some range of |V R |. The phase space correlation dimension D 2 that measures the clustering of the particles in position-velocity phase-space is modified by the presence of gravity. For St > 1, D 2 has a smaller value for Fr = 0.3 and 0.1 compared to D 2 for Fr = ∞. This implies that the positionand phase-space clustering is increased by the gravity for St > 1.
We compute the moments of these joint PDFs of R and |V R | as a function of R and angle θ between the separation vector R and the direction of gravity. We showed that the moments depend on the cos (θ) for a constant value of R. This indicates anisotropy in the system due to the presence of gravity. To characterize this anisotropy, we use a spherical harmonic decomposition and compute coefficients of order ℓ = 0, 2, and 4. We found that all these coefficient have a power-law dependence on R. Exponents of the power-law is found to be the same for all values of ℓ that can be obtained for 0-th and 1-th order moments. We define the degree of anisotropy as the ratio of ℓ = 2 and ℓ = 0 coefficients. As both the coefficients scale with the same exponent, degree of anisotropy is independent of R, for R < 1. A similar analysis of anisotropy is done in Ref. [24] for R > 1. This study shows that for larger R, the degree of anisotropy depends on R and goes to 0 for very large values of R.
We define the velocity of collision V n as the radial component of the relative velocity of two particles separated by the sum of their radii, when two particles are approaching each other. We calculated the mean of V n as a function of St 1 and St 2 and found that it changes qualitatively from the zero-gravity case. Due to the presence of gravity mean of V n becomes a function of | St 1 − St 2 | alone. Furthermore, for particles with equal St, mean V n decreases as Fr decreases, i.e., gravity increases. This decrease is more pronounced at higher St than at lower ones. For particle having different values of St the qualitative behaviour is opposite, the mean collision velocity increase as Fr decreases. For St 1 = St 2 collision velocity is determined by the difference in the settling velocities of two particles.
Joint PDFs of R and |V R | in the presence of gravity are studied for the first time in this paper. Refs. [24,28] have studied the PDFs of V R for small values of R but power-law nature of the PDFs is not shown. Real space correlation dimension d 2 for settling particles is studied in Ref. [22]. This study does not take into account the anisotropy due to the presence of gravity. Our study shows that the scaling exponents d 2 for anisotropic contributions are the same, but the amplitudes are different. Our results also show that gravity can have a significant effect on the collision velocities of particles having St between 0.5 and 3 and hence should be taken in to account in the studies relevant for cloud droplets in this range of St.