Statistical models for predicting pair dispersion and particle clustering in isotropic turbulence and their applications

The purpose of this paper is twofold: (i) to advance and extend the statistical two-point models of pair dispersion and particle clustering in isotropic turbulence that were previously proposed by Zaichik and Alipchenkov (2003 Phys. Fluids15 1776–87; 2007 Phys. Fluids 19, 113308) and (ii) to present some applications of these models. The models developed are based on a kinetic equation for the two-point probability density function of the relative velocity distribution of two particles. These models predict the pair relative velocity statistics and the preferential accumulation of heavy particles in stationary and decaying homogeneous isotropic turbulent flows. Moreover, the models are applied to predict the effect of particle clustering on turbulent collisions, sedimentation and intensity of microwave radiation as well as to calculate the mean filtered subgrid stress of the particulate phase. Model predictions are compared with direct numerical simulations and experimental measurements.


Introduction
Transport of small heavy particles in turbulent fluid flows occurs in a variety of natural and technological systems. Due to inertia, the particle velocity field is compressible even though the turbulent fluid may be considered as incompressible [1,2]. The compressibility of the particle velocity field is of fundamental importance because this causes fluctuations of particle concentration and particle preferable accumulation. The formation of clusters, that is, compact regions of preferential concentration surrounded by low-concentration zones, is one of the most challenging and complicated phenomena caused by the interaction of particles with fluid turbulent eddies [3]- [33]. Considerable progress has been made in computational studies of the clustering of inertial particles in homogeneous isotropic turbulence, but experimental investigations are sparse [13,25,32,33]. One should distinguish between two classes of flows in which clusters may appear, namely, inhomogeneous and homogeneous turbulent flows. The preferential concentration of particles suspended in inhomogeneous turbulent flow arises due to the so-called turbophoresis induced by the gradient of the fluid turbulence intensity [34,35]. Owing to this, heavy particles migrate from regions of high turbulence energy and accumulate in regions of low turbulence energy (for example, in the viscous sub-layer near the wall of a channel). In other words, the particles are pushed by turbophoresis towards the wall and trapped in the viscous wall region. However, the effect of preferential accumulation occurs also in homogeneous turbulence, when the mean gradient of fluid velocity fluctuations and the conventional turbophoresis are absent. Even in homogeneous and isotropic turbulence, inertial particles are not distributed homogeneously, but cluster. In this case, the accumulation 3 effect manifests itself as a tendency of heavy particles to avoid regions of high vorticity and preferentially concentrate in regions of high strain rate. The local accumulation of heavy particles is observed in low-vorticity regions under the action of centrifugal force and is mainly governed by small-scale turbulent structures. Therefore, the clustering of particles is most remarkable when the particle response time is comparable to the Kolmogorov timescale of fluid turbulence and hence the Stokes number is about unity. Fine zero-inertia particles completely follow the fluid velocity fluctuations and are uniformly distributed in space. For large Stokes numbers, the particle concentration distributions become defocused because highinertia particles are weakly responsive to the small-scale vorticity of the fluid. Taking account of the phenomenon of preferential particle concentration requires using two-point modelling approaches. In [36,37], the clustering phenomenon in homogeneous turbulence is treated as a result of a particle migration drift in the separation direction due to the gradient of the radial relative fluctuating velocity. This drift leads to the collection of particle pairs at small separations and hence acts in perfect analogy to an additional attractive force. Thus, in spite of the stochastic nature of turbulence, the spatial distribution of inertial particles even in statistically homogeneous turbulent flow exhibits some regular trends.
The preferential concentration of particles in homogeneous turbulent flow can be of importance in many natural and industrial processes. This effect is capable of making substantial increases in both particle settling velocities [4,13,38,39] and collision-coagulation rates [7,8], [40]- [48]. Preferential concentration of aerosol particles and droplets may play a particularly important role in atmospheric phenomena at large Reynolds numbers. Apparently, this effect is substantially responsible for a rapid growth of droplets in rain clouds [49,50] as well as a deviation from the Beer-Lambert exponential law of radiation extinction [51]. Because of this, the preferential accumulation of particles suspended in homogeneous turbulence attracts meticulous attention and induces a large number of works devoted to its investigation. It is notable that a significant role in increasing the collision-coagulation rate of inertial particles may play a mechanism that differs from spatial clustering. This mechanism is caused by the formation of fold caustics with a multi-valued particle velocity field [50], [52]- [54] and can give an extra contribution to the collision rate due to the so-called sling effect [55].
Various criteria may be employed for specifying preferential particle concentration. For example, because heavy particles are collected in regions of low vorticity (high strain rate), the second invariant of the fluid velocity gradient tensor is suitable for this purpose [3]. A method extensively employed for measuring the accumulation effect (e.g. see [3,5,10,11,16]) is based on a deviation of the probability density of the number of particles located in fixed space cells from the Poisson law, which corresponds to the statistically independent random spatial distribution of particles. However, the most convenient way of quantifying the accumulation effect is to use the radial distribution function (RDF) that is defined as the probability of observing a particle pair normalized by the corresponding value in a uniform suspension, because this quantity directly enters into the equations for the collision rate, the coagulation kernel, the settling velocity, the scattered radiation intensity and so forth.
In this paper, we reconsider and advance the statistical two-point models of pair dispersion and particle clustering proposed in [36,56]. The improvements of the models as compared to those presented in [36,52] consist in correcting the approximation of the Lagrangian two-point structure function of the fluid viewed by particles as well as in extending the applicability from zero-size to finite-size particles. These models are based on a kinetic equation for the two-point probability density function (PDF) of the relative velocity distribution of two particles. The fluid velocity field is assumed to be isotropic, homogeneous and incompressible. The particle fraction is kept small enough so that the effects of turbulence modulation can be neglected. The models advanced in the paper rely substantially on the theory of local similarity that disregards the intermittency effect. It is well known that the feasibility of obtaining universal solutions to the problem of relative dispersion of particle pairs is associated with the phenomenon of intermittency caused by fluctuations of the turbulence dissipation rate [57]- [59]. The classical Kolmogorov similarity hypotheses for small-scale turbulence are justified when disregarding the effect of intermittency. These hypotheses establish the universality of small-scale turbulence, i.e. presume that the characteristics of turbulence in the viscous and inertial ranges at high Reynolds numbers are independent of large-scale vortex structures. Intermittency causes the features of small-scale turbulence to depend on the Reynolds number. Therefore, it is possible to derive universal (self-similar with regard to the Reynolds number) solutions only for nottoo-high-inertia particles whose response times belong to the viscous or inertial interval, with the part played by intermittency being insignificant. The presence of such universality for the RDF is confirmed in direct numerical simulaiopn (DNS) by Collins and Keswani [20], at least for particles with the response time belonging to the viscous interval. Lundgren [60] demonstrated the validity of the classical Kolmogorov theory of local similarity in the limit of high Reynolds numbers. However, in [60], it was found that the intermittency index, which measures anomalous scaling, tends to zero with the Reynolds number very slowly (by a logarithmic law). Because the sensitivity to intermittency increases with the order of the statistical moment of the fluid fluctuating velocity field, our models are legitimate only for predicting the low-order moments of the PDF. Therefore, in view of the fact that intermittency is ignored, the analysis performed is limited to the behaviour of the second-order velocity moments.
The statistical two-point models are applied to predict the effect of particle clustering on turbulent collisions, sedimentation and intensity of microwave radiation. Moreover, as was shown in [61], the statistical approach that stems from a kinetic equation for the two-point PDF of particle-pair relative velocity adequately captures the main findings drawn from the mesoscopic Eulerian formalism for stationary isotropic turbulence. This formalism is based on the partitioning of particle turbulent velocity field into the mesoscopic-correlated and randomuncorrelated components [17,24]. In the present paper, we elucidate a connection between two modelling approaches for predicting particle statistics in evolving isotropic turbulence as well. Lastly, we use the two-point model predictions to analyse the subgrid-scale particulate stress that is required to advance the Eulerian two-fluid approach for calculating two-fluid turbulent flows in the framework of large eddy simulation. This paper addresses homogeneous isotropic turbulent flows laden with small heavy particles. The particle size is assumed to be less than the Kolmogorov length scale and the particle mass fraction is assumed to be low in order for the back-effect of particles on the fluid turbulence to be negligible. Moreover, we consider monodisperse particles. Although it is possible to extend the models to the case of bidisperse particles (as was done in [45]), this is beyond the scope of the present paper.
The paper is organized as follows. In the next section, after recalling the equations of motion, we introduce a kinetic equation for the PDF of the relative velocity distribution of two particles as well as equations for the two-point statistical moments of the PDF. Sections 3 and 4 consider relative velocity statistics and preferential concentration of particle pairs in stationary and decaying isotropic turbulent flow fields. The remainder of the paper is devoted 5 to applications of the two-point statistical models of particle-pair dispersion and accumulation. In section 5, we discuss the effect of particle inertia on the turbulent collision rate. Section 6 and 7 demonstrate the effect of clustering, respectively, on particle settling velocity and radiation scattering. Section 8 describes the connection between the subgrid-scale particulate stress and the particle second-order structure function. A summary of the work is given in section 9.

Theoretical background
The equations governing the relative motion of a pair of identical heavy particles suspended in a homogeneous turbulent flow field are given by where r p and w p are the separation distance and the relative velocity between two particles, u(r p , t) is the increment in fluid velocities at two points where the particles are located, R pα and v pα are the particle position and velocity vectors (α = 1, 2) u(R pα , t) is the local instantaneous velocity of the carrier fluid at the point where a particle is located and τ p is the particle response time. We note that equations (1) hold for particles whose density is much more than that of the carrier fluid and whose size is small as compared to the Kolmogorov length scale. Furthermore, we assume that the turbulent advection dominates the motion of the particles as well as the effects of gravitational settling and interparticle hydrodynamic interaction are not important.

Kinetic equation for the PDF
To proceed from stochastic equations (1) to a statistical description of the relative motion of two particles, the PDF of the relative velocity distribution is introduced where δ(x) symbolizes the Dirac delta function. The particle-pair PDF, P(r, w, t), is defined as the probability of finding a pair of particles separated by a distance r, with a relative velocity w, at time t. The angular brackets in (2) stand for the ensemble average over samples of a random fluid velocity field. Differentiating (2) with respect to time, accounting for (1), and decomposing the increment in fluid velocities at two points into the averaged and fluctuating components ( u(r, t) = U(r, t) + u (r, t)), we derive the following transport equation for the pair PDF: Equation (3) is unclosed because it contains an unknown correlation u i p that quantifies eddy-particle interactions, and hence we are faced with a closure problem. To determine u i p , the fluid relative velocity field is modelled by a Gaussian random process with known two-point correlation moments. Modelling the fluid turbulence by the Gaussian process is the key assumption that allows us to express the correlation between the fluid velocity increment and the particle-pair probability density, u i p , in a closed form. Note that, in actual fact, 6 the PDF of the fluid relative velocity is not the normal distribution [57]- [59]. As is known, the deviation of the fluid relative velocity PDF from the Gaussian distribution is due to intermittency. However, a noticeable departure from Gaussian manifests itself mainly in the PDF tails when the velocity difference magnitudes are quite large, whereas the normal distribution gives a reasonable approximation of the fluid relative velocity PDF for small velocity differences (e.g. see [62]). Since the contribution of the PDF tails to the low-order statistical moments is not very significant, the Gaussian random process can be successfully used for predicting the second-order relative velocity statistics [63]. Nevertheless, even though intermittency modifies only the PDF tails, this phenomenon could be linked to the presence of long-lived coherent structures in the small scales of turbulent flow. These structures can affect the relative motion of particle pairs. Therefore, the issue of the effect of intermittency on pair dispersion and particle clustering calls for further investigation, but this is beyond the scope of the current paper.
Assuming the Gaussian fluid velocity distribution and using the functional formalism, we can represent the correlation between the fluid velocity increment and the particle-pair probability density in the form [56] In order to calculate the integrals in (4), it is necessary to determine the Lagrangian twopoint structure function of fluid velocities viewed along particle trajectories In [56], the Lagrangian two-point structure function is given by Here i j (r, t) and i j (r, t) are the Eulerian two-point strain and rotation tensors, respectively, and σ (τ ) and ω (τ ) are, respectively, the autocorrelation functions of the strain and rotation tensors that describe the time lag of the velocity increments of two fluid elements separated by the distance r ≡ |r|. Obviously, the relation holds with S i j (r, t) being the second-order Eulerian structure function that quantifies the fluid fluctuating velocity increment at two points. In the viscous space range (r η, where η ≡ (ν 3 /ε) 1/4 designates the Kolmogorov length scale, ν is the fluid kinematic viscosity, and ε is the turbulence energy dissipation rate), the strain and rotation tensors are determined as The only dissimilarity from the model presented in [56] consists in determining the 'transport terms' D p i j (r, t)/Dt and D p i j (r, t)/Dt that account for the effects of turbulence unsteadiness as well as of the convective and diffusive transfer of fluid velocity fluctuations along a particle-pair trajectory when specifying the Lagrangian two-point structure function (5). Whereas in [56], it was assumed that the fluid velocity fluctuations are transferred with a particle relative velocity, in this paper we will define the 'transport terms' using the relative fluid velocity of two fluid elements moving along the trajectories of inertial particles (viewed by particles) where U p i denotes the mean relative fluid velocity viewed by inertial particles when they move along their trajectories. The quantities U pi , i jk and i jk will be given below. Note that the refinement made in determining the Lagrangian two-point structure function hardly influences the relative velocity statistics of particle pairs, but this slightly affects the accumulation effect resulting in a moderate enhancement of the RDF. Substituting (5) into (4) and restricting to terms of the first-order spatial derivatives, we obtain As is clear from (9), the quantities F i j and G i j may be treated as the diffusivity tensors that measure the transport of the pair PDF in phase space (r, w). The coefficients f ς , g ς , f ς 1 and l ς quantify a response of a pair of particles to velocity fluctuations of the turbulent fluid, that is, they measure coupling between the particulate and fluid phases. Then by substituting (9) 8 into (3), we obtain the closed kinetic equation for the two-point PDF of the particle-pair relative velocity distribution in homogeneous turbulent flow ∂ P ∂t The left-hand side of (10) describes the evolution in time and the convection in phase space, whereas the right-hand side of (10) quantifies the diffusive transport of the pair PDF due to the interaction between particles and turbulent eddies. Modelling the fluid turbulence by means of a Gaussian process with the Lagrangian autocorrelations of finite timescales enables the eddy-particle interaction to be expressed in the form of the second-order differential operator. The kinetic equation (10) completely controls the relative velocity statistics of the particulate phase. Integrating (10) with respect to velocity subspace w, we can derive a set of governing equations for the two-point statistical moments of the relative velocity PDF.

Equations for the two-point statistical moments of the PDF
The first three moments that follow from (10) are governed by the following equations: Here is the RDF, W i is the mean relative velocity of two particles, S p i j is the secondorder particle velocity structure function (VSF), D r p i j is the relative diffusivity tensor of two particles and D r p i j is the relative diffusivity tensor of two fluid elements moving along the trajectories of inertial particles. The VSF and the RDF quantify, respectively, the relative velocity statistics and inhomogeneous spatial distribution of two particles. To close the equation set stemming from (10) at the second-moment closure level, we invoke a gradient algebraic approximation for the triple particle fluctuating velocity moments. This approximation follows from the corresponding differential equation for the third moments by neglecting time evolution, convection and generation due to mean velocity gradients as well as by using a quasi-Gaussian approximation for the fourth-rank moments [36,56] In view of equation (9), the mean relative fluid velocity seen by the particles is determined as The third-order correlations appearing in (8) are approximated in the same manner as (14) i jk = − In isotropic turbulence with no gravity, spherical symmetry takes place. This implies that the relative velocities and density distributions of particle pairs are independent of the orientation of the separation vector r. Then the mean fluid and particle relative velocity vectors can be expressed through their radial (longitudinal) components Moreover, it should be borne in mind that in isotropic homogeneous turbulence any second-rank tensor may be represented as where M ll and M nn designate the longitudinal and transverse components relative to the direction of separation. By this means, in isotropic turbulence, the set of equations (11)-(16), accounting for (17) and (18) as well as keeping in mind that ll = 0, can be simplified to ∂ ∂t

Particles suspended in stationary isotropic turbulence
First we use the model presented in the previous section for predicting two-point velocity statistics and preferential concentration in a steady-state suspension of particles immersed in homogeneous, isotropic and stationary turbulence. In this case, the pair relative velocities and density distributions may be only dependent on the separation distance, r . Moreover, the mean convective transport in the fluid is supposed to be absent ( U r = 0), and the total number of particles is not changed in time. The latter infers that the balance between the net radial relative inward and outward fluxes takes place and, therefore, the mean relative velocity, W r , is equal to zero.

Governing equations
By this means, the set of equations (19)- (22) reduces to the following system: In (23)-(25), the overbar stands for normalization with the Kolmogorov length scale η or the Kolmogorov velocity scale u k ≡ (νε) 1/4 and St ≡ τ p /τ k symbolizes the Stokes number, where τ k ≡ (ν/ε) 1/2 is the Kolmogorov timescale. The Stokes number measures a response of particles to the perturbation produced by the underlying fluid turbulence. Equation (23) expresses the balance between the migration and diffusion forces in the separation direction between two particles. The migration force is specified by the two first terms of (23), and the diffusion force is identified with the third term of (23). The migration drift may be treated as an 'attractive turbophoretic force' induced by eddy-particle interactions. This mechanism has much in common with particle transport and dispersion in inhomogeneous turbulence due to turbophoresis (see [35]); the structure function for velocity increment being analogous to the spatial variation of fluid velocity and the RDF being analogous to the particle concentration. Equations (24) and (25) represent, respectively, the balance for the longitudinal, S p ll , and transverse, S p nn , components of the VSF, S p i j . Thus, in stationary isotropic turbulence with no gravity, the model under consideration amounts to solving three nonlinear ordinary differential equations involving the RDF and the second-order longitudinal and transverse VSFs. Boundary conditions for equations (23)-(25) are given by whered ≡ d/η is the particle diameter normalized with the Kolmogorov length scale. Conditions (26) appear to be appropriate for displaying the hard-sphere elastic collision model. Note that in the previous papers (e.g. [36,41,56]) we used conditions (26) withr = 0. This case corresponds to elastic collisions of the so-called ghost (zero-size) particles, which are free to occupy any space in the suspension. When using conditions (26) withr =d, we can take into consideration a finite size of particles. Relations (27) express the fact that, at large separations, the longitudinal and transverse VSFs are homogeneous and the particles are uniformly distributed in space. Under the assumption that the relative velocity is distributed normally, the mean relative velocity magnitude can be expressed through the longitudinal VSF as As shown by Wang et al [8], the relative velocity PDF is dependent on particle inertia and becomes Gaussian only at large Stokes numbers. Nevertheless, even for zero-inertia particles, the actual value of |w r | / w 2 r 1/2 simulated by Wang et al [63] is equal to 0.77, which is quite close to the value of √ 2/π = 0.798 predicted by the normal distribution. Henceforth, we will use relation (28) to express |w r | in terms of S p ll .
In order to determine the response coefficients f ς , g ς and f ς1 (ς = σ, ω) in (23)- (25), it is necessary to specify the autocorrelation functions ς (τ ). When using, as in [41,56], the two-scale bi-exponential approximation like that proposed by Sawford [64] for the Lagrangian velocity autocorrelation, the response coefficients take the form where ς ≡ τ p /T σ denotes the particle inertia parameter and z ≡ τ T /T L measures the ratio between the Taylor differential and Lagrangian integral timescales of fluid turbulence. These timescales are given by where a 0 is the acceleration variance, Re λ is the Reynolds number based on the Taylor length scale and u is the rms fluctuating velocity of the fluid. The approximations of a 0 and T L were obtained in [41] by fitting available results of numerical simulations and experiments. Moreover, the approximation of the acceleration variance relied on the Kolmogorov similarity hypotheses in the limit of high Reynolds numbers, and the Lagrangian integral timescale normalized with the Kolmogorov time microscale was assumed to be a linear function of the Reynolds number. When Re λ increases (z → 0), the response coefficients (29) reduce to those that correspond to the single-scale exponential function ς (τ ) = exp (−τ/T ς ), In the viscous space range, according to (7), the longitudinal and transverse components of the strain and rotation tensors as well as the fluid VSFs are represented as ll =r 2

15
,¯ nn =r ,S p nn = 2r 2 15 (30) and the timescales of the strain and rotation rate correlations must be proportional to the Kolmogorov microscale τ k with A σ and A ω being the constants. These are equal, respectively, to 2.3 and 7.2 [65,66].

13
In the external spatial range when the separation r exceeds the integral length scale of fluid turbulence L, T σ and T ω go over into the conventional Lagrangian integral timescale T L . Further let us assume that, in the inertial range (η r L), T σ and T ω coincide and are equal to the Lagrangian two-point timescale of the velocity increment T Lr , which is given in [41] as T Lr = 0.3 ε −1/3 r 2/3 . For continuously describing T σ and T ω in the entire range of r , we employ the approximation [56] The longitudinal component of the spatial strain rate tensor coincides with the Eulerian longitudinal VSF of the fluid and is given by the approximation [67] which obeys Kolmogorov's similarity hypotheses. In the viscous range (r η), the transverse components nn and nn are determined by relations (30). When T σ = T ω , nn and nn appear in the governing equations only as the sum that is equal to the transverse VSF ( nn + nn = S nn ).
Since T σ = T ω in the inertial and external ranges, it would suffice to determine just the transverse VSF in these ranges. This is expressed in terms of the longitudinal VSF as

Low-inertia particles
Let us consider first a solution to equations (23)-(25) for fine low-inertia particles (d 1, St < 1). It should be mentioned that of great interest is the derivation of analytical solutions to the problem of preferential concentration of low-inertia particles at large Reynolds numbers as these can be achieved in atmospheric conditions, because to date it is not possible to reach the atmospheric Reynolds number in DNS. In addition, the values of the Stokes number and the particle size related to the Kolmogorov length scale for typical atmospheric aerosols and rain droplets are less than those realizable in both direct simulations and laboratory experiments. We will solve equations (23)- (25) for ghost particles in the vicinity of the origin r = 0 without regard for boundary conditions (27). This solution is found in the form Substituting (35) and (36) into (23)- (25) and allowing for (29)-(31), we obtain the following set of equations for determining the coefficients α l and α n and the exponent χ :  Figure 1 shows the dependences of α l , α n and χ on St, which were obtained as a result of solution to equations (37) for Re λ = 100. The most important quantity from the standpoint of estimating the possibility of clustering of low-inertia particles is the exponent χ in (36) because this defines the behaviour of the RDF at small separations. As is clear from figure 1, χ is negative, and hence is singular (i.e. tends to infinity) with r → 0. For low-inertia particles, the power-law scaling of onr with a negative exponent was found by both numerical [7,20,21,32,37] and analytical [9,36,37] means. The unrestricted growth of the RDF with the separation tending to zero may be interpreted as clustering event and adopted as the criterion of clustering caused by the 'attractive turbophoretic force'. In figure 1, one can see the nonuniqueness of solution to equations (37) in some range of St. The second solution corresponds to the higher magnitude of χ and hence to the higher singularity of RDF. However, this second solution is hardly realizable physically. It is significant that real solutions to equations (37) and, consequently, the clustering effect in the meaning given above exist only for St < St cr rather than in the entire range of the Stokes number. Thus, the results obtained indicate that the solution to equations (23)- (25), which can be represented at small separations in the form of (35) and (36) (35) and (36) in the vicinity of r = 0, the solution to equations (23)- (25) has the formS where the values ofS p (0) and (0) are dependent on St and Re λ . In the case of St 1, the solution of (37) may be found by means of expansion in terms of the small parameter St. Keeping the terms of the first order of St, one can obtain the following relations: The negative square power dependence of χ on St at St 1 (i.e. χ ∝ −St 2 ) was first derived by Balkovsky et al [9] and later rediscovered in [36,37]. It is obvious that, in accordance with equations (37), χ is a function of St and Re λ . The effect of Re λ on χ is rather weak because we do not take into consideration the intermittency of fluid turbulence. At the same time, χ is strongly influenced by the values of the constants A σ and A ω . If we take A σ = A ω , then, according to (39), χ = −8St 2 /3 at St 1 for not-too-small Reynolds numbers when z 2 St. However, as is demonstrated in the DNS of Falkovich and Pumir [21] and Chun et al [37], |χ|/St 2 has much higher values at St 1 than 8/3. Such high values of |χ|/St 2 may only be obtained if A ω significantly exceeds A σ . When taking A σ = 2.3 and A ω = 7.2, we obtain |χ|/St 2 = 6.22, which is comparable to the values given in [21,37] at St 1. By this means, in order to reach agreement with DNS data for the power-law exponent χ , it is necessary that the timescale of the rotation rate correlation should exceed the respective scale of the strain rate correlation.
The VSF coefficients α l and α n as well as the RDF exponent χ are functions of two parameters: St and Re λ . However, as was already mentioned, the effect of Re λ predicted by the solution of equations (37) is rather weak. Therefore, if St < 0.6 and Re λ > 30, the longitudinal VSF coefficient and the RDF exponent obtained from this solution with A σ = 2.3 and A ω = 7.2 can be approximated by   [37] for Re λ = 47 and [55] for Re λ = 83 as well as with experimental and numerical data [32] obtained in the range of Re λ from 108 to 149. As is seen, χ predicted by equations (37) is rather weakly dependent on Re λ and may be well approximated by (41). The predicted values of χ are in excellent agreement with the DNS [37], but they markedly overestimate those obtained by Salazar et al [32] and underestimates those gained in [55] for relatively large Stokes numbers. In figure 2, the theoretical dependence derived by Chun et al [37] χ = −6.56St 2 (42) and the approximation obtained by Derevyanko et al [68] using the DNS data of Falkovich and Pumir [21,55] are shown as well. As can be observed, dependences (42) and (43)

Arbitrary-inertia particles
Let us now find an asymptotic solution to the problem (23)- (27) at high Reynolds numbers (Re λ → ∞) in the inertial spatial range (η r L) for particles whose response times belong to the inertial temporal interval (τ k τ p T L ). In this limit, the problem (23)-(27) when allowing for (29) and (32)- (34), transforms to the following one: The problem (44) is self-similar with respect to the Stokes number. Figure 3 represents the particle and fluid VSFs as well as the particle RDF obtained by solving (44). It is apparent that, at large separation (ρ 1), the particle longitudinal and transverse VSFs are smaller than those in the fluid and they can be approximately expressed using the local-equilibrium relations: σ p ll = f r σ ll and σ p nn = f r σ nn . However, at small separation, the particle VSFs exceed those of the fluid. This effect is caused by the diffusion transport of velocity fluctuations and takes place only for fairly inertial particles. As is also seen from figure 3, with decreasing ρ, the values of σ p ll and σ p nn approach each other, whereas the RDF tends monotonically to a certain limit. At the origin of coordinates we have the following values: In the high-inertia limit when the particle response time is much longer than the Lagrangian integral timescale (τ p T L ), a simple asymptotic solution to problem (23)- (27) can be obtained.
Owing to intensive transport of velocity fluctuations, the particle longitudinal and transverse VSFs coincide and become near-homogeneous, whereas the RDF deviates hardly from unity The solution to problem (23)- (27) along with (29) and (32)-(34) was obtained by numerical quadrature using the tridiagonal-matrix routine along with an iteration procedure. Figure 4 plots the structure and RDFs calculated for zero-size (d = 0) and finite-size (d = 0.1) particles at various Stokes numbers versus the separation distance between two particles. As can be observed, with increasing St, the particle VSFs deviate more and more from those in fluid (curves 1) and approach the homogeneous asymptotic distributions (46) for high-inertia particles. Although the fluid VSF is equal to zero atr = 0, the VSFs of fairly inertial ghost particles, by virtue of diffusion transport, have nonzero values at the origin. The values ofS p i j at finite-size particles exceed those of zero-size particles, but this difference takes place only at small separations for not-too-high-inertia particles and vanishes at large St. As is also seen from figure 4(c), the RDF of low-inertia ghost particles has a singularity atr = 0. It means that the accumulation (or clustering) effect of fine particles takes place. The accumulation phenomenon is treated as a result of a particle migration drift caused by the gradient of the radial relative fluctuating velocity. This drift leads to the collection of particle pairs at small separations and hence acts analogous to an additional attractive force. At relatively large Stokes numbers, the singularity of the RDF disappears and, in accordance with (38), takes a constant value at r = 0. We note that the value of St cr calculated using the full solution to problem (23)-(27) is slightly more as compared to that predicted from the solution of equations (37). In the limit of high-inertia particles, tends to unity. The RDF of finite-size particles deviates down from that of zero-size particles in the vicinity of the particle diameter. Thus, the preferential concentration decreases with increasing particle size. However, the influence of particle size on the RDF, like that on the VSF, vanishes with the Stokes number.
Relation (47) does not allow for the effect of Re λ because this is not very important for lowinertia particles. It is worth noting that (47) recalls the relation = 1 + 18St 2 , which was obtained in [8] by fitting DNS data for low-inertia particles atr = 1.
Let us now examine the performance of the model (23)- (27) for predicting the velocity statistics of inertial particles dispersed in homogeneous isotropic turbulence. For this purpose, we will compare predictions with results of numerical simulations by Wang et al [8] and Fede and Simonin [69]. Wang et al [8] performed DNS in a frozen turbulence after a stationary state was reached, and the particles ofd = 1 were then introduced in the flow. Fede and Simonin [69] investigated the motion of particles withd = 0.92 in a steady-state forcing turbulence. Predictions were carried out for both finite-size and zero-size particles. Figure 5 illustrates the influence of particle inertia on |w r | over a wide range of Stokes numbers. The mean relative velocity magnitude at particle contact is defined in terms of the longitudinal structure function using relation (28). As is clear, both the predictions and the numerical simulations follow the same trend and indicate a pronounced maximum of |w r | as a function of St. The initial rise in |w r | is attributable to a decrease in the correlation of particle velocities with St. The subsequent decay of |w r | beyond the maximum results from a decrease in the particle fluctuating velocities, since the particles become more sluggish and less responsive to the fluid turbulence. From figure 5, it is also seen that the predictions with d = 1 overestimate the DNS data at small Stokes numbers. The predictions that correspond to ghost particles lead to slightly smaller values of mean relative velocities at small St. However, as it has been already mentioned, a difference in the mean relative velocities of zero-size and finite-size particles takes place only for relatively small Stokes numbers and vanishes at large St. Figure 6 demonstrates comparisons of the predictions of the particle kinetic energy related to the fluid one, k p /k = 15 /2S p ll (∞)/2Re λ , and the radial relative velocity variance scaled by the particle kinetic energy, w 2 r (r ) /k p = 4S p ll (r )/3S p ll (∞), with DNS [69] for Re λ = 34. As is clear from figure 6(a), the model well predicts the fall of the particle kinetic energy with increasing particle inertia without any influence of particle size. Figure 6(b) shows a reasonable agreement between the predictions and the simulations of radial relative velocities. In figure 6(b), the line 4 describes the dimensionless radial relative velocity variance that corresponds to the longitudinal VSF of the fluid in the viscous space range (30) Note that a somewhat better correlation of the predicted relative velocities with the DNS of Fede and Simonin [69] as compared to the DNS of Wang et al [8] may be apparently explained by a distinction between the forcing and frozen turbulent fluid flow fields, in which the particles were immersed.
In what follows we verify the model (23)- (27) comparing the predictions of the RDF with numerical and experimental data given in [8,32]. Figure 7 demonstrates a quite good agreement of the RDF predicted for low-inertia finite-size particles with measurements and simulations by Salazar et al [32]. Note that, when the particles are sufficiently small, the Stokes number and the dimensionless particle diameter are connected by the relation which corresponds to the Stokes law of the flow around a particle with ρ p and ρ f being the particle and fluid densities. Figure 8 shows the influence of particle inertia on the RDF for the separation and the Reynolds numbers, which correspond to the DNS data of Wang et al [8] obtained ford = 1. As expected, in the limiting cases of zero-inertia and high-inertia particles, the concentration field is statistically uniform, and therefore the RDF is equal to unity. In accord with the computations, the predicted RDF goes through a peak as the particle inertia time increases. Thus, there exists a critical particle response time that results in maximum preferential concentration. The value of this critical response time is of the same order as the Kolmogorov timescale, that is, the clustering is more considerable when tuning the particle parameters to the flow Kolmogorov scales. Figure 8 demonstrates a qualitative agreement between the predicted and simulated RDF at particle contact, although the predicted effect of Re λ is not as strong as it was in simulations and maxima of are slightly shifted towards particles with larger inertia. On the basis of comparisons with numerical and experimental results, it is concluded that the model (23)- (27) is able to capture the main features of particle relative velocity statistics and preferential concentration in stationary isotropic turbulence.
At the end of this section, we recall an intimate connection between the two-point statistical approach being considered and the mesoscopic Eulerian formalism that is rooted in the decomposition of the instantaneous particle velocity field into the mesoscopic spatially correlated and random spatially uncorrelated components [17,24]. The former velocity component allows for the spatial correlation of particles due to their interaction with turbulent eddies and the latter represents the statistically independent motion of separate particles. According to this formalism, the total kinetic energy of the particulate phase is represented as the sum of relevant contributions wherek p and δk p are the mesoscopic-correlated and random-uncorrelated energy components. These two parts of the particle kinetic energy are expressed in terms of the second-order VSF of the particulate phase as (51) Figure 9 shows the effect of particle inertia on the correlated and uncorrelated energy components normalized with the total particle kinetic energy. As is clear, the model predictions gained by means of relationships (51) are in good agreement with the results of numerical simulations. For low-inertia particles, the main contribution to the kinetic energy makes the spatially correlated motion and the part of the random-uncorrelated velocity field is negligible. With increasing inertia, the motion of particles becomes less and less correlated, and hence the fraction of the kinetic energy residing in the spatially correlated motion decreases, whereas the fraction residing in the random-uncorrelated component increases.

Particles suspended in non-stationary isotropic turbulence
In this section, we analysse the performance of a statistical second-moment model for predicting particle-pair dispersion in a non-stationary isotropic fluid flow field. For simplicity, the RDF is assumed to be quasi-stationary and only the VSF is treated as truly non-stationary. This simplification relies on the fact that the VSF is not very sensitive to any difference in the RDF, although the RDF is strongly affected by any change in the VSF. The assumption made allows us to amount the set of equations (19)-(22) to the following system: i.e.S 0 p ll (r ) =S ll (r ) andS 0 p nn (r ) =S nn (r ). The Stokes number based on the initial Kolmogorov timescale is equal to 2.2. Figure 10 shows time evolution of the profiles of the longitudinal structure and RDFsdistribution functions. Figure 11 presents the correlated and uncorrelated components of the particle kinetic energy predicted using relations (51). The fluid kinetic energy decreases in time due to viscous dissipation. As is clear from figure 10(a), the particle kinetic energy follows a similar tendency. Figure 11 indicates that there are two distinctive intervals of time evolution of the correlated and uncorrelated energy components. The boundary between these two intervals corresponds to the extreme of the energy components. The first (initial) time interval features an evolution of the VSF and RDF from the appropriate profiles of the fluid to the quasi-stationary distributions of S p i j and . At the second (quasi-stationary) interval, the VSF and RDF are described by the quasi-stationary distributions of S p i j and , which are appropriate to local values of St and Re λ . The particle inertia that is measured by the Stokes number decreases in time, because the Kolmogorov timescale increases. The most remarkable non-homogeneity of particle distribution is achieved at St ≈ 1 (curve 4 in figure 10b). Figure 11 exhibits that, due to particle inertia, the particle velocities become partially uncorrelated in space and δk p /k p begins to increase. However, after the attainment of a maximum, the uncorrelated  Figure 11. Temporal evolution of the mesoscopic-correlated (1, 3) and randomuncorrelated (2, 4) particle energy components normalized with the particle kinetic energy: 1, 2, predictions; 3, 4, DNS [66]. energy component falls due to a decrease in St; this fact is in agreement with the results presented in figure 9 for statistically stationary turbulence. The values ofk p /k p and δk p /k p predicted using the model (52)-(54) reasonably correlate with DNS [70].
The remainder of the paper is concerned with the applications of the two-point statistical models of particle dispersion and accumulation.

Particle collision rate
Let us consider the effect of particle inertia on the collision rate of non-settling particles without hydrodynamic interactions. In this case, the collision rate can be derived from the relationship [7,8] It is clear from (56) that the turbulence-induced collision rate is governed by the mean radial relative velocity as well as by the RDF. Consequently, the interaction of particles with turbulent eddies causes two statistical mechanisms, that contribute to the collision rate, namely, the relative velocity between neighbouring particles (the turbulent transport effect) and the non-uniform particle spatial distribution (the accumulation effect). Assuming the probability distribution of the relative velocity to be Gaussian and expressing the mean radial relative velocity magnitude in terms of the longitudinal VSF (28), we can rewrite (56) as Insertion of (35) and (36) into (57) produces the turbulent collision rate of low-inertia particles The collision rate (58) is valid for not-too-high-inertia particles when (35) and (36) are true, i.e. when St < St cr . In the limiting case of zero-inertia particles, equation (58) reduces to the well-known collision rate of Saffman and Turner [71] β ST = 8π 15 Figure 12 represents the behaviour of the collision rate predicted by (58) along with (40), (41) and (47) as a function of the Stokes number. It is seen that a monotonous rise in β starting from the value given by (59) at St = 0. The model predictions compare favourably with the DNS results obtained by Zhou et al [72] for low-inertia particles. Figure 12 also shows the collision rate of Derevyanko et al [68] which allows for the accumulation and sling effects. In (60), the RDF exponent is specified by (43). As can be observed the collision rates given by (60) are in reasonable agreement with those predicted by (58) and simulated by Zhou et al [72]. As is clear, equation (58) along with (40), (41) and (47), as well as equation (60) along with (43), ignores the effect of fluid turbulence intermittency. This approach is consistent with a weak sensitivity of the RDF to changes in the Reynolds number, which was observed at small Stokes numbers and separations in DNS by Collins and Keswani [20]. A saturation of results with increasing Reynolds number [20] is in harmony with the classical Kolmogorov similarity hypotheses. However, these observations contradict a quite strong growth of the exponent χ with Re λ at small values of St, which was obtained in DNS by Falkovich and Pumir [21]. Thus, the issue of the sensitivity of both the RDF and the collision rate of low-inertia particles to the Reynolds number is yet to be clarified. In the intermediate case, when the particle response time belongs to the inertial interval (1 St T L ) and in the limit of large Reynolds numbers, one can derive where the constant a is determined using (45). In the high-inertia limit (St T L ), the particles move statistically independent and are uniformly distributed in space. Hence, the structure and RDFs are independent of the separation distance between particles and they are equal to Inserting (62) into (57) yields the collision rate of Abrahamson [73] β A = 4 2π k p 3 where f u is the coefficient that measures a response of particles to fluid turbulence.
In the limiting cases of zero-inertia and high-inertia particles, the collision rate is exclusively governed by the turbulent transport mechanism. According to equation (61), the accumulation effect plays a part in particle collisions when the particle response time falls within the inertial interval. However, the influence of preferential concentration on the collision rate and thereby on the coagulation rate of particles is most pronounced, when the particle response time is close to the Kolmogorov turbulent timescale and the RDF reaches its maximum. Figure 13 represents the collision rate found using the solution to problem (23)- (27) in a wide range of Stokes number. The predictions of β are normalized with the Saffman-Turner collision rate (59) and compared with DNS by Wang et al [8] performed ford = 1. As is seen, the model properly reproduces the crucial trends of the DNS results, although the predicted maxima of β are slightly more than those simulated. Figure 13 shows that the influence of Reynolds numbers vanishes for low-inertia particles and the collision rate is well generalized using the Kolmogorov microscales. By contrast, the collision rate of high-inertia particles is governed by the macroscales of turbulence and, with increasing Stokes number, β approaches the Abrahamson collision rate (63). In figure 13, the approximation of Mehlig et al [54] is also depicted which takes into consideration the formation of fold caustics but does not involve the spatial clustering of particles. As is clear, approximation (64) markedly underpredicts β when St ∼ 1.
Note that Wang et al [8], fitting DNS data, derived models for both the radial relative velocity and the RDF. Then, combining both models and using (56), they obtained the collision rate that is in excellent agreement with their numerical simulations at moderate flow Reynolds numbers.

Effect of clustering on particle settling velocity
Maxey [1] first revealed an increase in the mean settling velocity of heavy particles immersed in homogeneous turbulence as compared to their terminal velocity in still fluid, V t . The enhancement of settling velocity due to turbulence was also found numerically by Wang and Maxey [4] and Yang and Lei [74] as well as experimentally by Yang and Shy [75,76]. This effect was explained by 'particle preferential sweeping', that is, by the sweeping of particles preferably towards the regions of downward fluid flow when interacting with turbulent eddies. All of these investigations were limited to dilute suspensions when the influence of particle fraction is not of significance. Aliseda et al [13] measured an increase in the particle settling velocity in nearly homogeneous isotropic turbulence at relatively large values of particle volume fractions. They found considerable higher settling velocities, V s , than those computed and measured in dilute suspensions, V s0 . In this case, the additional enhancement of V s with respect to V s0 was approximately proportional to the particle volume fraction, , and the maximum of V s was reached when the Stokes number, St, was close to unity. Aliseda et al [13] assumed that V s ≡ V s − V s 0 is due to the effect of particle clustering and proposed a simple phenomenological model for its prediction. Numerical simulations by Bosse et al [38] and Wang et al [39] confirmed a linear dependence of V s on for not-too-large particle volume fractions, and thereby they supported the conclusion of Aliseda et al [13] about a decisive part of clusters in increasing the settling velocity of particles in the turbulent fluid. Below we present a clustering-induced statistical model for predicting V s . This model relies on the theory of Batchelor [77] proposed for the gravitational sedimentation of interacting spheres.

Clustering-induced settling velocity
According to Aliseda et al [13], the settling velocity of particles is written as the sum of two terms Here the first term is the settling velocity of isolated particles in the turbulent fluid when → 0. This term equals the terminal velocity plus the enhancement caused by the effect of preferential sweeping by turbulence. The second term describes the additional enhancement of settling velocity due to the effect of preferential concentration.
For small values of the particle volume fraction, V s can be represented as where S is the settling coefficient. To determine S, we involve the theory of Batchelor [77] according to which with ρ ≡ r/a being the distance between the centres of two particles scaled by the particle radius, a ≡ d/2. Although relation (67) was originally deduced for particles suspended in an initially quiescent fluid, this apparently is valid in a turbulent flow if allowing for the effect of turbulence on the RDF. In (67), A 11 , A 12 , B 11 and B 12 denote the mobility functions of a particle pair, which quantify the hydrodynamic interaction of two particles. Detailed expansions of the mobility functions valid for small and large values of ρ were given in [78]. For large separations (ρ → ∞), the following asymptotic expansion takes place: Moreover, according to [79] In terms of relations (68) and (69), it is convenient to rewrite the settling coefficient in (66) in the form In line with (71), S 1 immediately takes account of the hydrodynamic interaction of two particles, which is described by the mobility functions. If the particle size is much less than the Kolmogorov length scale (d η), it is clear from (68) that the hydrodynamic interaction is short ranging because its effect can be of importance only over distances of the order of the particle diameter. By contrast, S 2 is directly independent of hydrodynamic interaction, and its role may be realized only indirectly through the impact on . When neglecting the hydrodynamic interaction, the RDF measures the particle-turbulence interaction resulting in preferential concentration. The particle-turbulence interaction is long ranging since, as is seen from figures 4(c) and 7, the size of radial distribution of non-homogeneity has the order of 10η. Moreover, because figures 4(c) and 7 suggest that > 1, the following inequalities are true: S 1 < 0 and S 2 > 0.
The last term in (70) is caused by an upward flux of the fluid, which compensates the downward current incorporating the fluid volumes replaced and carried away by particles. This back flow hinders the mean settling velocity of the suspension and leads to the appearance of the negative term, −5 V t , in the gravitational sedimentation velocity [79]. It will be shown later that in a turbulent medium, except for the case of very small Stokes numbers, S 2 |S 1 − 5|. This means that the main contribution to the settling coefficient (70) gives S 2 and the part of S 1 that is directly dependent on the hydrodynamic interaction of two particles is not of significance.
Consequently, to predict the settling coefficient, we can use the RDF obtained when neglecting the impact of hydrodynamic interaction. Neglect of this impact introduces a serious error into only at separations ofr = O(d). Since S is determined by integration over the full separationr , the error, which is introduced when neglecting the impact of hydrodynamic interaction into , can be important only for small values of St.

Low-inertia particles
Before analysing the effect of particle inertia on the clustering-induced settling velocity, let us consider the limit of zero-inertia particles (St → 0). In this case, the part of hydrodynamic interaction and thereby the error caused by the neglect of its impact on the RDF have their maxima.
Zero-inertia particles are fully entrained into the motion of fluid turbulence, and hence their spatial distribution when neglecting the hydrodynamic interaction is like that of a passive scalar with = 1. Then, from (66) and (69)-(71), it follows the Batchelor formula [79] derived for non-inertial particles dispersed in an initially quiescent fluid According to (72), the settling velocity of the non-inertial particle suspension, which is uniformly distributed in space, is less than the terminal velocity of a single particle. This collective effect is caused by the back flow of the fluid, which hinders the gravitational settling velocity (the so-called hindering effect). However, in a turbulent flow field, the hydrodynamic interaction results in the appearance of a drift velocity of non-inertial particles one to another (see [66]). Due to this drift, the rise in takes place as ρ decreases with a singularity when ρ → 2. In homogeneous isotropic turbulence, is described by exactly the same radial distribution as in a laminar linear fluid flow field [80] with A and B being monotonically decreasing functions of ρ [81]. The RDF given by (73) diminishes quickly with ρ, and, for ρ → 2, the asymptotic dependence takes place [66] = 0.21 Substituting the RDF (73) into (70) and (71) produces and thus Comparison of (74) with (72) shows that taking into account the influence of hydrodynamic interaction on the RDF leads to an increase in the settling velocity of the suspension, just the same, this is less than the terminal velocity of a single particle. Hence, the effect of clustering of zero-inertia particles caused by their hydrodynamic interaction in isotropic turbulence is not of significance for the hindering effect of the back fluid flow to be completely compensated.
To elucidate the influence of inertia on the settling velocity of fine particles, we use the RDF given by relations (36), (41) and (47). Although these relations are correct atr 1, they will be of use at largerr , bearing in mind that the results obtained by this means have only qualitative interest. For this purpose, we define the radius of the cluster, r cl , as a value of r , at which the RDF determined by relations (36), (41) and (47) is equal to unity. The radius of the cluster is thenr Substitution of (36), (41) and (47) into (71) yields When deducing (77), it was taken into consideration thatr cl d . Figure 14 presents the dependences of S 1 and S 2 on St, which were obtained using equations (76) and (77) along with (49) and (75) for various values of the particle-to-fluid density ratio ρ p /ρ f . It is clear that, for heavy particles (ρ p /ρ f 10 3 ), S 1 makes a contribution to the settling coefficient S only when the Stokes number is very small, and its part is already negligible for St = 0.01. Thus, in the range of 0.01 < St < 0.6 when ρ p /ρ f 10 3 , the clustering-induced settling velocity can be evaluated using the relation It is worth nothing that formula (78) coincides with that predicted by the phenomenological model of Aliseda et al [13], except for the value of the coefficient b. According to (78), b is a function of St, whereas in [13] b is given as 0.3-0.5. The model of Aliseda et al [13] is based on a balance between the gravity and drag forces as applied to a cluster quasi-particle. Because the gravity force is proportional to r 3 cl and the drag one is proportional to r cl , this model leads to V s ∼ r 2 cl . In the present paper, the same dependence is obtained using the Batchelor theory [77].

Arbitrary-inertia particles
Now we present the settling coefficient predicted for fairly inertial particles (St > 0.1). In this case, the major contribution to (70) makes S 2 and the influence of particle hydrodynamic interaction on the RDF may be ignored. By this means, we obtain from (70) and (71) where is found from the solution to the problem (23)−(27) for finite-size particles, and St andd are connected by relation (49). Figures 15 and 16 show, respectively, the settling coefficient (79) as functions of the Stokes number and the dimensionless particle diameter. It is clear that S increases with both the Reynolds number and the particle-to-fluid density ratio. The maximum S takes place when St is close to unity, which is due to the maximum effect of preferential accumulation when fitting the particle response time and the Kolmogorov timescale. In line with (49), as ρ p /ρ f increases, the maximum S shifts to a less value ofd. In figure 15, the settling coefficient given by (77) for small Stokes numbers is presented as well. One can observe that a reasonable agreement between the analytical dependence (77) and the numerical results obtained by integrating (79) occurs only for small St. Figure 17 compares the clustering-induced settling velocity predicted using (66) and (79) with the measurements obtained by Aliseda et al [13] for water droplets settling in an air isotropic turbulent flow. It is seen that, in spite of quantitative agreement, the predictions overestimate V s /V t as compared to the experiments. This discrepancy can be apparently attributed to neglecting the effect of gravitational sedimentation on the RDF. As was noted in [44,55,82], the sedimentation caused by gravity can significantly alter the RDF and,  according to (79), this can result in modifying the clustering-induced settling velocity. However, the inclusion of the effect of gravity on the RDF is beyond the scope of the present paper.

Effect of particle clustering on scattering of microwave radiation
Experimental data reported by Baker et al [83] and Knight and Miller [84] for the reflectivity of developing small cumulus clouds and radar measurements of a smoke plume produced by an intense industrial fire by Rogers and Brown [85] indicated a considerable discrepancy in the spectral reflectivity dependence from the theoretical predictions based on two well-known mechanisms: the incoherent Rayleigh scattering by single particles and the Bragg scattering due to air density fluctuations. According to Erkelens et al [86], this discrepancy may be attributed to a significant contribution of coherent microwave scattering by particle clusters formed by turbulence. In the case of cumulus clouds, the coherent scattering of water clusters can be a predominant mechanism of scattering, especially in the centimetre spectral range [86]. In this section, we estimate the effect of particle clustering on the coherent scattering of microwave radiation for conditions which are typical for atmospheric clouds containing small water droplets.
The problem includes the interference of electromagnetic waves scattered by arbitrary positioned particles and is based on the so-called independent scattering hypothesis when each of particles is assumed to scatter the radiation in exactly the same manner as if all the other particles do not exist. Moreover, the following conditions are assumed to be satisfied: π d/λ 1 and πd|n|/λ 1, where λ is the wavelength and n is the complex index of reflection. Then assuming that the observation point is located in the Fraunhofer zone, i.e. far enough from the scattering volume, one can derive the following expression for the average intensity of scattered radiation: where I is the intensity of scattered radiation, I 0 is the intensity of scattered radiation calculated by neglecting the coherent scattering (i.e. without taking into account the interference of radiation scattered by particles with correlated spatial coordinates), θ is the angle of scattering and K is the coefficient of clustering-induced scattering. The value of K does not depend on radiative properties of particles and is precisely the same as the theory of x-ray scattering by a liquid [87,88] with N being the number of density particles. It is easily seen that the coefficient K can be treated as the Fourier sinus-transform of 4π N [ (r ) − 1]r/s. Taking the Fourier inversion formula, we obtain Equation (82) can be used to identify the RDF from the experimental data for the coefficient of clustering-induced scattering.
In what follows we estimate K for the conditions of clustering of water droplets in atmospheric clouds. For these conditions, the dimensionless particle size,d, is not greater than 0.01 and the Stokes number, St, is less than 0.1 (e.g. see [89]). Thus, inasmuch as this is done to estimate the clustering-induced settling velocity of low-inertia particles in section 6.2, one can use the RDF given by relations (36), (41), and (47). Substituting these relations into (81), we obtain where the cluster radius,r cl , is given by (75).
At small values ofs, equation (83) can be written in the form of a convergent series .
In the case of St 1, equation (84) is reduced to the series As is clear from (85), the coefficient of clustering-induced scattering is proportional to St 2 in the limit of small Stokes numbers. The values of K /N calculated using equations (83) and (85) are shown in figure 18. One can see that (85) yields a good estimate of the clusteringinduced coefficient whens < 0.5 and St 0.05, i.e. in the realistic ranges of these parameters for microwave remote sensing of cumulus clouds. Note that, for cumulus clouds, we haveN ∼ 1 since N ∼ 1 mm −3 and η ∼ 1 mm, so that the values of K /N in figure 18 can be approximately treated as the absolute values of K . The above estimate indicates that the effect of turbulent clustering of small water droplets in atmospheric clouds may lead to a considerable increase in microwave radiation scattering. It means that microwave reflectivity data obtained from radar measurements should be considered as a potential source of important information on turbulent clustering of water droplets in the clouds.
Concluding this section, it should be noted that the statistical model of clustering-induced scattering is true not only in the limit of low-inertia particles but can also be used at arbitrary Stokes numbers in other applications.

Subgrid-scale stress of the particulate phase
Finally, we use the second-order particle VSF to examine the behaviour of the averaged subgridscale particulate stress as functions of particle inertia and filter width. Knowledge of the subgrid stress is required to advance an Eulerian continuum approach for modelling the transport of inertial particles in the framework of large eddy simulation [90]− [94]. For this purpose, we involve the relation between the statistically averaged subgrid stress and the second-order VSF that was suggested for the fluid turbulence by Germano [95]. Obviously, this relation is valid for the particulate phase as well. This is given by where σ r i j is the subgrid particulate stress and G is the filter kernel. Following Germano [95], we perform an analysis in the case of the Gaussian filter where denotes the filter width. Passing in (86) to the coordinates integrating over s as well as over the angle coordinates of r, and accounting for (87), we obtain the following expression for the averaged subgrid particulate stress in isotropic homogeneous turbulence: 2 ) S p ll (r ) + 2S p nn (r ) r 2 dr.
The trace of (88) gives the averaged subgrid kinetic energy of the particulate phase where the explicit expression of the Gaussian filter (87) is used. Figure 19 shows the effect of filter width on the averaged subgrid particle kinetic energy normalized with the total fluid turbulence energy, k. Calculations have been performed by means of formula (89) for the distributions of the VSFs of zero-size (d = 0), which are depicted in figure 4. It is clear that the averaged subgrid fluid energy k r , which corresponds to zeroinertia particles (St = 0), increases monotonically with starting from zero at = 0 and tends to the total fluid turbulence energy. The averaged subgrid particle energy grows with as well. However, with increasing St, the dependence of k r p on flattens due to the nearly homogeneous structure functions at high Stokes numbers. The asymptotic value of k r p /k falls for → ∞ as St increases, because high-inertia particles are weakly responsive to fluid velocity fluctuations. Figure 20 demonstrates the effect of filter width on the averaged subgrid particle kinetic energy normalized with the total particle turbulence energy, k p . Predictions are compared with the filtered particle kinetic energy computed using discrete particle simulation coupled with fluid DNS of homogeneous isotropic decaying turbulence [94]. As is seen, k r p /k p increases monotonically with both /η and St. Moreover, it is apparent that the predictions and simulations presented in figure 20 are in good agreement.

Summary
The main goal of this paper was to present statistical models for predicting relative velocity statistics and preferential concentration of inertial particles suspended in a fluid turbulent flow field. These models start from a kinetic equation for the two-point PDF of the relative velocity distribution of particle pairs. The kinetic equation is closed assuming that the fluid velocity flow field can be considered as a Gaussian process with correlations expressed in terms of the twopoint structure functions and time scales. The validity of the models presented is explored by comparing with numerical simulations and experiments performed in stationary and decaying isotropic turbulent flows. Also, we sought to use these models in various physical applications, in which the effect of particle preferential concentration can be important. The main findings and conclusions are given as follows: 1. The statistical models based on the kinetic theory reproduce by and large all the effects found in the behaviour of both the second-order relative velocity statistics and the nonuniformity distributions of particle pairs, which have been revealed in numerical and experimental investigations. These models are valid over the entire range of particle inertia and are adequate in predicting the most preferential accumulation when the particle response time is comparable to the Kolmogorov time microscale.
2. The effect of preferential concentration is a maximum for zero-size particles. In this case, the unrestricted rise in the RDF with the separation distance tending to zero may be adopted as the criterion of clustering. The singularity of the RDF of zero-size particles at the origin can occur only in the case of lower-than-critical values of the Stokes number. 3. The relative velocity variance grows and the RDF decreases with increasing particle size.
However, the influence of particle size is of importance only at small separations for lowinertia particles. 4. The collision model accounts for two mechanisms influencing inter-particle collisions: the particle relative motion induced by turbulence and the preferential concentration that leads to an additional enhancement of the collision rate. This model is capable of predicting the collision rate over the entire range of Stokes numbers (from zero-inertia to high-inertia particles). 5. Preferential concentration can be crucial in determining particle settling velocity in homogeneous isotropic turbulence. According to experiments, predictions provide a linear rise in settling velocity with particle fraction and a maximum value of settling velocity as the Stokes number is about unity. 6. Particle clustering by fluid turbulence results in an increase in scattering of microwave radiation. The RDF can be determined from measurements solving the inverse problem of coherent Rayleigh scattering. 7. The two-point relative velocity statistics of particle pairs provides the main features of effects of filter width and particle inertia on the subgrid particulate stress. The kinetic theory can be useful in advancing the Eulerian approach for large eddy simulation of particle transport in turbulent flows.