Heavy particle concentration in turbulence at dissipative and inertial scales

Spatial distributions of heavy particles suspended in an incompressible isotropic and homogeneous turbulent flow are investigated by means of high resolution direct numerical simulations. In the dissipative range, it is shown that particles form fractal clusters with properties independent of the Reynolds number. Clustering is there optimal when the particle response time is of the order of the Kolmogorov time scale $\tau_\eta$. In the inertial range, the particle distribution is no longer scale-invariant. It is however shown that deviations from uniformity depend on a rescaled contraction rate, which is different from the local Stokes number given by dimensional analysis. Particle distribution is characterized by voids spanning all scales of the turbulent flow; their signature in the coarse-grained mass probability distribution is an algebraic behavior at small densities.

Spatial distributions of heavy particles suspended in an incompressible isotropic and homogeneous turbulent flow are investigated by means of high resolution direct numerical simulations. In the dissipative range, it is shown that particles form fractal clusters with properties independent of the Reynolds number. Clustering is there optimal when the particle response time is of the order of the Kolmogorov time scale τη. In the inertial range, the particle distribution is no longer scale-invariant. It is however shown that deviations from uniformity depend on a rescaled contraction rate, which is different from the local Stokes number given by dimensional analysis. Particle distribution is characterized by voids spanning all scales of the turbulent flow; their signature in the coarse-grained mass probability distribution is an algebraic behavior at small densities. Understanding the spatial distribution of finite-size massive impurities, such as droplets, dust or bubbles suspended in incompressible flows is a crucial issue in engineering [1], planetology [2] and cloud physics [3]. Such particles possess inertia, and generally distribute in a strongly inhomogeneous manner. The common understanding of this long known but remarkable phenomenon of preferential concentrations relies on the idea that, in a turbulent flow, vortexes act as centrifuges ejecting particles heavier than the fluid and entrapping lighter ones [4]. This picture was successfully used (see, e.g. [4,5]) to describe the small-scale particle distribution and to show that it depends only on the Stokes number S η = τ /τ η , which is obtained by non-dimensionalizing the particle response time τ with the characteristic time τ η of the small turbulent eddies.
In this Letter, we confirm that this mechanism is relevant to describe the particle distribution at length-scales which are smaller than the dissipative scale of turbulence, η. In particular, maximal clustering is found for Stokes numbers of the order of unity. However, we show that stationary particle concentration experiences also very strong fluctuations in the inertial range of turbulence. In analogy with small-scale clustering, it is expected that for r ≫ η the relevant parameter is the local Stokes number S r = τ /τ r , where τ r is the characteristic time of eddies of size r [6]. Surprisingly, we present evidences that such a dimensional argument does not apply to describe how particles organize in the inertial range of turbulence. We show that the way they distribute depends on a scaledependent rate at which volumes are contracted.
In very dilute suspensions, the trajectory X(t) of small spherical particles much heavier than the fluid evolve ac-cording to the Newton equation [7] τẌ = u(X, t) −Ẋ . (1) The response time τ is proportional to the square of the particles size and to their density contrast with the fluid.
Here, we neglect buoyancy and particle Brownian diffusion. As stressed in [8,9], diffusion may affect concentration of particles. However we assume that the typical lengthscale below which particle diffusion becomes dominant is much smaller than the Kolmogorov scale η of the fluid flow, and than any observation scale considered here. In most situations, particles are so massive that such an approximation is fully justified. The fluid velocity u satisfies the incompressible Navier-Stokes equation p being the pressure and ν the viscosity; F is an external homogeneous and isotropic force that injects energy at large scales L with a rate ε = F · u . We performed direct numerical simulations of (2) by means of a pseudo-spectral code on a triply periodic box of size L = 2π with 128 3 , 256 3 and 512 3 collocation points corresponding to Reynolds numbers (at the Taylor microscale) R λ ≈ 65, 105 and 185, respectively. Once the flow u is statistically stationary, particles with 15 different values of the Stokes number in the range S η ∈ [0. 16, 3.5] (N = 7.5 millions of particles per S η ), are seeded homogeneously in space with velocity equal to that of the fluid, and are evolved according to (1) for about two largescale eddy turnover times. After this time, particle mass distribution reaches a statistical steady state too. Mea- The correlation dimension D2 vs Sη for the three different R λ . Also shown the probability P to find particles in non-hyperbolic (rotating) regions of the flow, for Re λ = 185 (multiplied by an arbitrary factor for plotting purposes). D2 has been estimated taking into account also subleading terms, as described in [12].
surements are then performed. Details on the numerics and on the transient are reported in [10].
Below the Kolmogorov length-scale η where the velocity field is differentiable, the motion of inertial particles is governed by the fluid strain and the dissipative dynamics leads their trajectories to converge to a dynamically evolving attractor. For any given response time of the particles, their mass distribution is singular and generically scale-invariant with fractal properties at small scales [8,11]. In order to characterize particle clusters at these scales, we measured the correlation dimension D 2 , which is estimated through the small-scale algebraic behavior of the probability to find two particles at a distance less than a given r: P 2 (r) ∼ r D2 . The dependence of D 2 on S η and Re λ is shown in Fig. 1. Notice that D 2 depends very weakly on Re λ in the range of Reynolds numbers here explored. A similar observation was done in [5], where particle clustering was equivalently characterized in terms of the radial distribution function. This indicates that τ η , which varies by more than a factor 2 between the smallest and the largest Reynolds numbers considered here, is the relevant time scale to characterize clustering below the Kolmogorov scale η. For all values of Re λ , a maximum of clustering (minimum of D 2 ) is observed for S η ≈ 0.6. For values of S η larger than those investigated here, D 2 is expected to saturate to the space dimension [12]. For small values of S η , particle positions strongly correlate with the local structure of the fluid velocity field. A quantitative measure is also given in Fig. 1 by the probability P to find particles in non-hyperbolic regions of the flow, i.e. at those points where the strain matrix has two complex conjugate eigenvalues. Note that P attains its minimum for values of S η close to the minimum of D 2 . This supports the traditional view relating particle clustering to vortex ejection. As qualitatively evidenced from Fig. 2, particle positions correlate with low values of the fluid acceleration. This phenomenon was already evidenced in [10] where a statistical analysis of the fluid acceleration at particle positions was performed, and also in [13] for 2D turbulent flows.
From Fig. 2, it is clear that fluctuations in the particle spatial distribution extend to scales far inside the inertial range; this confirms the experimental measurements of [14]. Moreover, for Stokes numbers of the order of unity (Figs. 2c and 2d), we note that the sizes of voids span all spatial scales of the fluid turbulent flow, similarly to what observed in 2D inverse cascade turbulence [13,15]. To gain a quantitative insight, we consider the Probability Density Function (PDF) P r,τ (ρ) of the particle density coarse-grained on a scale r inside the inertial range, that is the probability distribution of the fraction of particles in a cube of size r, normalized by the volume r 3 of the cube. For tracers, which are uniformly distributed, and for an infinite number of particles N → ∞, this PDF is a delta function on ρ = 1. For finite N , the probability to have a number n of particles in a box of size r is given by the binomial distribution and leads to a closed form for the PDF of the coarse-grained density ρ = nL 3 /(N r 3 ). In our settings, the typical number of particles in cells inside the inertial range is several thousands. Thus, finite number effects are not expected to affect the core of the mass distribution ρ = O(1), but they may be particularly severe when ρ ≪ 1 because of the presence of voids. To reduce this spurious effect, we computed the quasi-Lagrangian (QL) mass density PDF P (QL) r,τ (ρ), which corresponds to a Lagrangian average with respect to the natural measure [16], and is obtained by weighting each cell with the mass it contains. For statistically homogeneous distribution, QL statistics can be related to the Eulerian ones by noticing that ρ p r QL = ρ p+1 r (see, e.g., [16]).  r,τ (ρ), for various response times τ and at a given scale (r = L/16) within the inertial range. Deviations from a uniform distribution can be clearly observed, they become stronger and stronger as τ increases. This means that concentration fluctuations are important not only at dissipative scales but also in the inertial range. A noticeable observation is that the low-density tail of the PDF (which is related to voids) displays an algebraic behavior P (QL) r,τ (ρ) ∼ ρ α(r,τ ) . This means that the Eulerian PDF of the coarse-grained mass density has also a power-law tail for ρ ≪ 1, with exponent α − 1. The dependence of this exponent for fixed r and varying the Stokes number S η is shown in the inset. For low inertia (S η → 0), it tends to infinity in order to recover the non-algebraic behavior of tracers. At the largest available Stokes numbers, we observe α → 1, indicating a non-zero probability for totally empty areas. Note however that α is expected to go to infinity when S η → ∞ with r fixed, since a uniform distribution is expected for infinite inertia. Algebraic tails are relevant to various physical/chemical processes involving heavy particles. The particle distribution in low density regions is an important effect to be accounted for, e.g., modelling the growth of liquid droplets by condensation of vapor onto aerosol particles [8] and of aerosol scavenging [3].
Fixing the response time τ and increasing the observation scale r reproduces the same qualitative picture as fixing r and decreasing τ . A uniform distribution is recovered in both limits r → ∞ or τ → 0. These two limits are actually equivalent. At length-scales r ≫ η within the inertial range, the fluid velocity field is not smooth: ac-cording to Kolmogorov 1941 theory (K41), velocity increments behave as δ r u ∼ (εr) 1/3 . K41 theory suggests that fluctuations at scale r are associated to time scales of the order of the turnover time τ r = r/δ r u ∼ ε −1/3 r 2/3 . This implies that for any finite particle response time τ and at sufficiently large scales r, the local inertia measured by S r = τ /τ r becomes so small that particles should behave as tracers and distribute uniformly in space [6]. Deviations from uniformity for finite S r are expected not to be scale-invariant [8]. In particular, observations from random δ-correlated in time flows [17] suggest that particle distribution should depend only on the local Stokes number S r . However, as explained below, this argument does not seem to apply to realistic flows.
The presence of inhomogeneities in the spatial distribution of particles is due to a dynamical competition between stretching, folding, and contraction due to their dissipative dynamics. At small scales where the flow is spatially smooth, these mechanisms are described by the Lyapunov exponents and their finite-time fluctuations, associated to the growth and contraction rates of infinitesimal volumes. At larger scales, competition between particle spreading and concentration is also at equilibrium. The spatial distribution of particle at a given scale r ≫ η should only depend on the time scale given by the inverse of the contraction rate γ r,τ of a sizer blob B r of particles with response time τ . However in the inertial range, the flow is not differentiable and an approach based on Lyapunov exponents cannot be used. In particular, the rate γ r,τ should depend on r. To estimate this contraction rate in the limit of small τ , we make use of Maxey's approximation [21] of the particle dynamics [7] meaning that particle trajectories can be approximated by those of tracers in the synthetic compressible velocity field v. The inertial correction is proportional to the fluid acceleration that, in turbulent flows, is itself dominated by pressure gradient. The contraction rate of the blob B r with volume r 3 is γ r,τ = (1/r 3 ) Br d 3 x ∇ · v(x).
As ∇ · v ≃ τ ∇ 2 p, we have γ r,τ ∼ (τ /r 2 ) δ r p, where δ r p denotes the typical pressure increment at scale r. This means that the scale dependence of the contraction rate is directly related to the scaling properties of the pressure field. Dimensional arguments suggest that δ r p ∼ (εr) 2/3 , so that the contraction rate scales as γ r,τ ∼ τ ε 2/3 /r 4/3 . However, as stressed in [19], K41 scaling for pressure is observable only when the Reynolds number is tremendously large (typically for Re λ > ∼ 600). In our simulations, where Re λ < 200, pressure scaling is actually dominated by random sweeping [20], and we observe δ r p ∼ U (εr) 1/3 (as in other simulations at comparable Reynolds, see [19]). This implies that the contraction rate is γ r,τ ∼ τ ε 1/3 U/r 5/3 .  r,τ (ρ) for three choices of γ r,τ obtained from various sets of values of r and τ . The collapse of the different curves strongly supports that the distribution of the coarse-grained mass density depends only upon the scale-dependent contraction rate γ r,τ . In particular, as represented in the inset of Fig. 4, the deviations from unity of the first moment of the distribution collapse for all S η investigated and all scales inside the inertial range of our simulation. This quantity is the same as the Eulerian 2nd-order moment, giving the probability to have two particles within a distance r. We note that particle distribution recovers uniformity at large scales very slowly, much slower than if particle were distributed as Poisson point-like clusters for which ρ 2 −1 ∝ r −3 ∝ γ 9/5 r,τ (also shown in the inset for comparison).
We presented numerical evidence that at moderate Reynolds number the distribution of heavy particles at lengthscales within the inertial range is fully described in terms of a scale-dependent volume contraction rate γ r,τ ∝ τ /r 5/3 . However we expect that γ r,τ ∝ τ /r 4/3 at sufficiently large Reynolds numbers (Re λ > 600). Stateof-the-art numerical simulations can not currently attain these turbulent regimes, so that only experimental measurements of particle distribution can help to confirm such a prediction. We have seen that particle dynamics in the inertial range can be directly related to the structure of the pressure field (and thus of acceleration). Characterizing the distribution of acceleration is thus crucial to understand particle clusters, and conversely, particle distribution could be used as an experimental tool to probe the spatial structure of turbulent flows.
We acknowledge useful discussions with G. Boffetta, A. Celani and A. Fouxon. This work has been partially supported by the EU under contract HPRN-CT-2002-00300 and the Galileo programme on Trasport and dispersion of impurities suspended in turbulent flows. Simulations were performed at CINECA (Italy) and IDRIS (France) under the HPC-Europa programme (RII3- CT-2003-506079). Unprocessed data of this study are freely available from http://cfd.cineca.it.