Toward a phenomenological approach to the clustering of heavy particles in turbulent flows

A simple model accounting for the ejection of heavy particles from the vortical structures of a turbulent flow is introduced. This model involves a space and time discretization of the dynamics and depends on only two parameters: the fraction of space-time occupied by rotating structures of the carrier flow and the rate at which particles are ejected from them. The latter can be heuristically related to the response time of the particles and hence measure their inertia. It is shown that such a model reproduces qualitatively most aspects of the spatial distribution of heavy particles transported by realistic flows. In particular the probability density function of the mass $m$ in a cell displays an power-law behavior at small values and decreases faster than exponentially at large values. The dependence of the exponent of the first tail upon the parameters of the dynamics is explicitly derived for the model. The right tail is shown to decrease as $\exp (-C m \log m)$. Finally, the distribution of mass averaged over several cells is shown to obey rescaling properties as a function of the coarse-grain size and of the ejection rate of the particles. Contrarily to what has been observed in direct numerical simulations of turbulent flows (Bec et al., http://arxiv.org/nlin.CD/0608045), such rescaling properties are only due in the model to the mass dynamics of the particles and do not involve any scaling properties in the spatial structure of the carrier flow.


Raphaël Chétrite
Laboratoire de Physique, ENS Lyon, 11 Allée d'Italie, 69007 Lyon, France Abstract. A simple model accounting for the ejection of heavy particles from the vortical structures of a turbulent flow is introduced. This model involves a space and time discretization of the dynamics and depends on only two parameters: the fraction of space-time occupied by rotating structures of the carrier flow and the rate at which particles are ejected from them. The latter can be heuristically related to the response time of the particles and hence measure their inertia. It is shown that such a model reproduces qualitatively most aspects of the spatial distribution of heavy particles transported by realistic flows. In particular the probability density function of the mass m in a cell displays an power-law behavior at small values and decreases faster than exponentially at large values. The dependence of the exponent of the first tail upon the parameters of the dynamics is explicitly derived for the model. The right tail is shown to decrease as exp(−Cm log m). Finally, the distribution of mass averaged over several cells is shown to obey rescaling properties as a function of the coarse-grain size and of the ejection rate of the particles. Contrarily to what has been observed in direct numerical simulations of turbulent flows (Bec et al., nlin.CD/0608045), such rescaling properties are only due in the model to the mass dynamics of the particles and do not involve any scaling properties in the spatial structure of the carrier flow.

Introduction
Understanding the dynamics of small-size tracer particles or of a passive field transported by an incompressible turbulent flow plays an important role in the description of several natural and industrial phenomena. For instance it is well known that turbulence has the property to induce an efficient mixing over the whole range of length and time scales spanned by the turbulent cascade of kinetic energy (see e.g. [1]). Describing quantitatively such a mixing has consequences in the design of engines, in the prediction of pollutant dispersion or in the development of climate models accounting for transport of salinity and temperature by large-scale ocean streams.
However, in some settings, the suspended particles have a finite size and a mass density very different from that of the fluid. Thus they can hardly be modeled by tracers because they have inertia. In order to fully describe the dynamics of such inertial particles, one has to consider many forces that are exerted by the fluid even in the simple approximation where the particle is a hard sphere much smaller than the smallest active scale of the fluid flow [2]. Nevertheless the dynamics drastically simplifies in the asymptotics of particles much heavier than the carrier fluid. In that case, and when buoyancy is neglected, they interact with the flow only through a viscous drag, so that their trajectories are solutions to the Newton equation : where u denotes the underlying fluid velocity field and τ is the response time of the particles. Even if the carrier fluid is incompressible, the dynamics of such heavy particles lags behind that of the fluid and is not volume-preserving. At large times particles concentrate on singular sets evolving with the fluid motion, leading to the appearance of strong spatial inhomogeneities dubbed preferential concentrations. At the experimental level such inhomogeneities have been known for a long time (see [3] for a review). At present the statistical description of particle concentration is a largely open question with many applications. We mention the formation of rain droplets in warm clouds [4], the coexistence of plankton species [5], the dispersion in the atmosphere of spores, dust, pollen, and chemicals [6], and the formation of planets by accretion of dust in gaseous circumstellar disks [7]. The dynamics of inertial particles in turbulent flows involves a competition between two effects: on the one hand particles have a dissipative dynamics, leading to the convergence of their trajectories onto a dynamical attractor [8], and on the other hand, the particles are ejected from the coherent vortical structures of the flow by centrifugal inertial forces [9]. The simultaneous presence of these two mechanisms has so far led to the failure of all attempts made to obtain analytically the dynamical properties or the mass distribution of inertial particles. In order to circumvent such difficulties a simple idea is to tackle independently the two aspects by studying toy models, either for the fluid velocity field, or for the particle dynamics that are hopefully relevant in some asymptotics (small or large response times, large observation scales, etc.). Recently an important effort has been made in the understanding of the dynamics of particles suspended in flows that are δ-correlated in time, as in the case of the well-known Kraichnan model for passive tracers [10]. Such settings, which describe well the limit of large response time of the particles, allows one to obtain closed equations for density correlations by Markov techniques. The δ-correlation in time, of course, rules out the presence of any persistent structure in the flow; hence any observed concentrations can only stem from the dissipative dynamics. Most studies in such simplified flows dealt with the study of the separation between two particles [11,12,13,14,15].
Recent numerical studies in fully developped turbulent flows [16] showed that the spatial distribution of particles at lengthscales within the inertial range are strongly influenced by the presence of voids at all active scales spanned by the turbulent cascade of kinetic energy. The presence of these voids has a noticeable statistical signature on the probability density function (PDF) of the coarse-grained mass of particles which displays an algebraic tail at small values. To understand at least from a qualitative and phenomenological viewpoint such phenomena, it is clearly important to consider flows with persistent vortical structures which are ejecting heavy particles. For this purpose, we introduce in this paper a toy model where the vorticity field of the carrier flow is assumed piecewise constant in both time and space and takes either a finite fixed value ω or vanishes. In addition to this crude simplification of the spatial structure of the fluid velocity field we assume that the particle mass dynamics obeys the following rule: during each time step there is a mass transfer between the cells having vorticity ω toward the neighboring cells where the vorticity vanishes. The amount of mass escaping to neighbors is at most a fixed fraction γ of the mass initially contained in the ejecting cell. We show that such a simplified dynamics reproduces many aspects of the mass distribution of heavy particles in incompressible flow. In particular, we show that the PDF of the mass of inertial particles has an algebraic tail at small values and decreases as exp(−A m log m) when m is large. Analytical predictions are confirmed by numerical experiments in one and two dimensions.
In section 2 we give some heuristic motivations for considering such a model and a qualitative relation between the ejection rate γ and the response time τ of the heavy particles. Section 3 consists in a precise definition of the model in one dimension and in its extension to higher dimensions. Section 4 is devoted to the study in the statistical steady state of the PDF of the mass in a single cell. In section 5 we study the mass distribution averaged over several cells to gain some insight on the scaling properties in the mass distribution induced by the model. Section 6 encompasses concluding remarks and discussions on the extensions and improvements of the model that are required to get a more quantitative insight on the preferential concentration of heavy particles in turbulent flows.

Ejection of heavy particles from eddies
The goal of this section is to give some phenomenological arguments explaining why the model which is shortly described above, might be of relevance to the dynamics of heavy particles suspended in incompressible flows. In particular we explain why a fraction of the mass of particles exits a rotating region and give a qualitative relation between the ejection rate γ and the response time τ entering the dynamics of heavy particles. For this we focus on the two-dimensional case and consider a cell of size ℓ where the fluid vorticity ω is constant and the fluid velocity vanishes at the center of the cell. This amounts to considering that the fluid velocity is linear in the cell with a profile given to leading order by the strain matrix. Having a constant vorticity in a two-dimensional incompressible flow means that we focus on cases where the two eigenvalues of the strain matrix are purely imaginary complex conjugate. The particle dynamics reduces to the second-order two-dimensional linear system It is easily checked that the four eigenvalues of the evolution matrix are the following complex conjugate Only λ +,− and λ +,+ have a positive real part which is equal to This means that the distance of the particles to the center of the cell increases exponentially fast in time with a rate µ. If we now consider that the particles are initially uniformly distributed inside the cell, we obtain that the mass of particles remaining in it decreases exponentially fast in time with a rate equal to −2µ. Namely the mass of particles which are still in the cell at time T is The rate γ at which particles are expelled from the cell depends upon the response time τ of the particles and upon two characteristic times associated to the fluid velocity. The first is the time length T of the ejection process which is given by the typical life time of the structure with vorticity ω. The second time scale is the turnover ω −1 which measures the strength of the eddy. There are hence two dimensionless parameters entering the ejection rate γ: the Stokes number St = τ ω giving a measure of inertia and the Kubo number Ku = T ω which is the ratio between the correlation time of structures and their eddy turnover time. One hence obtain the following estimate of the ejection rate The graph of the fraction of particles ejected from the cell as a function of the Stokes number is represented in figure 1  In three dimensions, one can extend the previous approach to obtain an ejection rate for cells with a uniform rotation, i.e. a constant vorticity ω. There are however two main difficulties. The first is that in three dimensions the eigenvalues of the strain matrix in rotating regions are not anymore purely imaginary but have a real part given by the opposite of the rate in the stretching direction. Such a vortex stretching has to be considered to match observation in real flows. The second difficulty stems from the fact that the vorticity is now a vector and has a direction, so that ejection from the cell can be done only in the directions perpendicular to the direction of ω. These two difficulties imply that the spectrum of possible ejection rates is much broader than in the two-dimensional case. However the rough qualitative picture is not changed.

A simple mass transport model
We here describe with details the model in one dimension and mention at the end of the section how to generalize it to two and higher dimensions. Let us consider a discrete partition of an interval in N small cells. Each of these cell is associated to a mass which is a continuous variable. We denote by m j (n) the mass in the j-th cell at time t = n. At each integer time we choose randomly N independent variables; Ω j = 1 with probability p and Ω j = 0 with probability 1 − p. The evolution of mass between times n and n + 1 is given by: In other terms, when Ω j = 1, the j-th cell looses mass if Ω j−1 = 0 or Ω j+1 = 0, and when Ω j = 0, it gains mass if Ω j−1 = 1 or Ω j+1 = 1. The flux of mass between the j-th and the (j + 1)-th cell is proportional Ω j − Ω j+1 (see figure 2). In particular, if Figure 2. Sketch of the dynamics in the one-dimensional case: the fluxes of mass are represented as arrows. A cross means no flux.
Ω j = Ω j+1 , no mass is transfered between cells. When the system is supplemented by periodic boundary conditions between the cells N and 1, it is clear that the total mass is conserved. Hereafter we assume that the mass is initially m j = 1 in all cells., so that the total mass is j m j = N. Spatial homogeneity of the random process Ω j implies that m j = 1 for all later times, where the angular brackets denote average with respect to the realizations of the Ω j 's.
A noticeable advantage of such a model for mass transportation is that the mass field m = (m 1 , . . . , m N ) defines a Markov process. Its probability distribution p N (m, n + 1) at time n + 1, which is the joint PDF of the masses in all cells, is related to that at time n by a Markov equation, which under its general form can be written as where P [m ′ → m|Ω] denotes the transition probability from the field m ′ to the field m conditioned on the realization of Ω = (Ω 1 , . . . , Ω N ). In our case it takes the form The variable µ j denotes here the flux of mass between the j-th and the (j + 1)-th cell.
It is a function of Ω j , Ω j+1 , and of the mass contained in the two cells. It can be written as In the particular case we are considering, the joint probability of the Ω j 's factorizes and we have so that the Markov equation (8) can be written in a rather explicit and simple manner.
The extension of the model to two dimensions is straightforward. The mass transfer out from a rotating cell can occur to one, two, three or four of its direct nearest neighbors (see figure 3 left). One can similarly derive a Markov equation which is similar to (8) for the joint PDF p N,N (M, n) at time n of the mass configuration M = {m i,j } 1≤i,j≤N . The transition probability reads in that case where the fluxes now take the form and µ (2) i,j defined by inverting i and j in the definition of µ i,j . After a large number of time steps, a statistically steady state is reached. The stationary distribution is obtained assuming that p N,N (M, n) = p N,N (M) is independent of n in the Markov equation (8). In this stationary state the mass fluctuates around its mean value 1 corresponding to a uniform distribution; strong deviations at small masses can be qualitatively observed (see figure 3 right). The model can be easily generalized to arbitrary space dimension. However, as we have seen in previous section, besides its interest from a purely theoretical point of view, the straightforward extension to the three-dimensional case might not be relevant to describe concentrations of inertial particles in turbulent flows.

Distribution of mass
Let us consider first the one-dimensional case in the statistically stationary regime. After integrating (8), one can express the single-point mass PDF p 1 in terms of the three-point mass distribution p 3 at time n We then explicit all possible fluxes, together with their probabilities, by considering all possible configurations of the spin vorticity triplet (Ω j−1 , Ω j , Ω j+1 ). The results are summarized in table 1. This leads to rewrite the one-point PDF as The first term on the right-hand side comes from realizations with no flux. The second term is ejection to one neighbor and the third to two neighboring cells. The fourth term involving an average over the two-cell distribution is related to events when mass is transfered from a single neighbor to the considered cell. Finally, the last term accounts for transfers from the two direct neighbors. Note that, in order to write (15), we made use of the fact that p 2 (x, y) = p 2 (y, x). Table 1. Enumeration of all possible configurations of the spin vorticity Ω in three neighboring cells, together with their probabilities and the associated mass fluxes.
Numerical simulations of this one-dimensional mass transport model are useful to grab qualitative information on p 1 . Figure 4 represents the functional form of p 1 in the stationary regime for various values of the ejection rate γ and for p = 1/2. The curves are surprisingly similar to measurements of the spatial distribution of heavy particles suspended in homogeneous isotropic turbulent flows [3,16]. This gives strong evidence that, on a qualitative level, the model we consider reproduces rather well the main mechanisms of preferential concentration. More specifically, a first observation is that in both settings the probability density functions display an algebraic behavior for small masses. This implies that the ejection from cells with vorticity one has a strong  statistical signature. A second observation is that at large masses, the PDF decays faster than exponentially, as also observed in realistic flows. As we will now see these two tails can be understood analytically for the model under consideration.
We here first present an argument explaining why an algebraic tail is present at small masses. For this we exhibit a lower bound of the probability P < (m) that the mass in the given cell is less than m. Namely, we have where A is a set of space-time realizations of Ω such that the mass in the j-th cell at time n is smaller than m. For instance we can choose the set of realizations which are ejecting mass in the most efficient way: during a time N before n, the j-th cell has spin vorticity 1 and its two neighbors have 0. The mass at time n is related to the mass at time n − N by The probability of such a realization is clearly p N (1 − p) 2N . Replacing N by the expression obtained in (17), we see that After averaging with respect to the initial mass m j (n − N), one finally obtains Hence the cumulative probability of mass cannot have a tail faster than a power law at small arguments. It is thus reasonable to make the ansatz that p 1 (m) have an algebraic tail at m → 0, i.e. that p 1 (m) ≃ Cm α . To obtain how the exponent α behaves as a function of the parameters γ and p, this ansatz is injected in the stationary version of the Markov equation (15). One expects that the small-mass behavior involves only the terms due to ejection from a cell, namely the three first terms in the r.h.s. of (15), and that the terms involving averages of the two-point and three-point PDFs give only sub-dominant contributions. This leads to (20) Equating the various constants we finally obtain that the exponent α satisfies Note that the actual exponent α given by this relation is different from the lower-bound β + 1 obtained above in (18) and (19). However it is easily checked that α approach the lower bound when p → 0. As seen from figure 5, formula (21) is in good agreement with numerics. Note that the large error bars obtained for p small and γ large are due to the presence of logarithmic oscillations in the left tail of the PDF of mass. This log periodicity is slightly visible for γ = 0.9 in figure 4. It occurs when the spreading of the distribution close to the mean value m = 1 is much smaller than the rate at which mass is ejected. This results in the presence of bumps in the PDF at values of m which are powers of 1 − 2γ. Notice that for all values of p, one has α ≤ 0 when γ ≥ 2/3. However, according to the estimate (6), values of the ejection rate larger than 2/3 can be attained only for large enough Kubo numbers. This is consistent with the fact that power-law tails with a negative exponent were not observed in the direct numerical simulations of turbulent fluid flows [16] where Ku ≈ 1.
It is much less easy to get from numerics the behavior of the right tail of the mass PDF p 1 (m). As seen from figure 4, there was no events recorded where the mass is larger than roughly ten times its average. We however present now an argument suggesting that the tail is faster than exponential, and more particularly that log p 1 (m) ∝ −m log m when m ≫ 1. We first observe that in order to have a large mass in a given cell, one needs to transfer to it the mass coming form a large number M of neighboring cells. Estimating the probability of having a large mass is equivalent to understand the probability of such a transfer. For moving mass from the j-th cell to the (j − 1)-th cell, the best configuration is clearly (Ω j−1 , Ω j , Ω j+1 ) = (0, 1, 1). After N time steps with this configuration, the fraction of mass transfered is 1 − (1 − γ/2) N . This process is then repeated for moving mass to the second neighbor, and so on. After order M iterations, the mass in the M-th neighbor is This means that All the processes of this type will contribute terms in the right tail of the mass PDF. The dominant behavior is given by choosing N = N ⋆ such that N ⋆ M(m, N ⋆ ) is minimal. Such a minimum cannot be written explicitly. One however notices that, on the one hand, if N is much larger than its lower bound (i.e. This suggests that the minimum is attained for N ⋆ ∝ log m. Finally, such estimates lead to predict that the right tail of the mass probability density function behaves as where C is a positive constant that depends upon the parameters p and γ. As seen in figure 6, such a behavior is confirmed by numerical experiments. The estimations of the left and right tails of the distribution of mass in a given cell can be extended to the two-dimensional case. The results do not qualitatively change. The exponent α of the algebraic behavior at small masses is given as a solution of By arguments which are similar to the one-dimensional case and which are not detailed here, one obtains also that log p 1 (m) ∝ −m log m. Numerical experiments in two dimensions confirm these behaviors of the mass probability distribution. As seen from figure 7 an algebraic behavior of the left tail of the PDF of m is observed and the value of the exponent is in good agreement with (26).

Coarse-grained mass distribution
We investigate in this section the probability distribution of the mass coarse-grained on a scale L much larger than the box size ℓ, which is defined as As seen from the numerical results presented on figure 8, the functional form of the PDF p L (m) is qualitatively similar to that of the mass in a single cell. In particular for various values of L it also displays an algebraic tail at small arguments with an exponent which depends both on L and on the parameters of the model. We here present some heuristic arguments for the behavior of the exponent. For this, we consider the cumulative probability P < L (m) to havem L smaller than m. We first observe that in order to havem L small, the mass has to be transfered from the bulk of the coarse-grained cell to its boundaries. Assume we start with a mass order unity in each of the 2K + 1 sub-cells. The best realization to transfer mass is to start with ejecting an order-unity fraction of the mass contained in the central cell with index j = 0 to its two neighbors. For this the three central cells must have vorticities (Ω −1 , Ω 0 , Ω 1 ) = (0, 1, 0), respectively, during N time steps. After that the second step consists in transferring the mass toward the next neighbors; the best realization is then to have during N time steps (Ω −2 , Ω −1 , Ω 0 , Ω 1 , Ω 2 ) = (0, 1, 1, 1, 0). The transfer toward neighbors is repeated recursively. At the j-th step, the best transfer is given by choosing (Ω −j−1 , Ω −j , Ω −j+1 ) = (0, 1, 1) and (Ω j−1 , Ω j , Ω j+1 ) = (1, 1, 0) during a time N. One can easily check that for large N and after repeating K times this procedure, the mass which remains in the 2K + 1 cells forming the coarse-grained cell ism L ≃ (1 − γ/2) N . The total probability of this process is [p 4 (1 − p) 2 ] KN ,which leads to estimate the cumulative probability ofm L as This approach guarantees that the probability density function p L (m) = dP < L /dm of the coarse-grained massm behaves as a power law at small values. Note that only the contribution from realizations with an optimal mass transfer is here evaluated and the actual value of the exponent should take into account realizations of the vorticity which may lead to a lesser mass transfer. However we expect the estimation given by (28) to hold for L sufficiently large, because the contribution from realizations with a sub-dominant mass transfer become negligible in this limit.
As to the right tail of p L (m), one expects a behavior similar to that obtained in the case of the one-cell mass distribution, namely log p L (m) ∝ −m logm form ≫ 1. Indeed, the probability of having a large mass in a coarse-grained cell should clearly be of the same order as the probability of having a large mass in a single cell. This, together with the estimates (28) for the exponent of the left tail, gives a motivation for looking, at least in some asymptotics, for possible rescaling behaviors of p L (m) as a function of L and of the ejection rate γ. For instance one can argue whether the limits L → ∞ and γ → 0 are equivalent. The estimation (28) suggests that the exponent α L depends only on the ratio κ = L/ log(1 − γ/2). Note that the limit of small γ should mimic that of small response time of the heavy particles. Rescaling of the distribution of the coarse-grained mass was observed in direct numerical simulations of heavy particles in turbulent homogeneous isotropic flows [16]. Such a rescaling is confirmed by numerical simulations. Figure 8 represents the PDF of the coarse-grained massm L for various values of L and γ chosen such that the ratio κ is 8, 10 or 14. While the left-tail and the core of the distribution are clearly collapsing, the rescaling is much less evident for the right-hand tail. Obtaining better evidence would require much longer statistics in order to resolve the distribution at large masses.

Conclusions
We introduced here a simple model for the dynamics of heavy inertial particles in turbulent flow which solely accounts for their ejection from rotating structures of the fluid velocity field. We have shown that this model is able to reproduce qualitatively most features of the particle mass distribution which are observed in real turbulent flows. Namely the probability density of the mass in cells is shown to behave as a power-law at small arguments and to decrease faster than exponentially at large values. Moreover, we studied how this distribution depends on the parameters of the model, namely the ejection rate of particles from eddies and the fraction of space occupied by them. Such dependence reproduce again qualitatively observations from numerical simulations in homogeneous turbulent flows. Finally, we have seen that coarse-graining masses on scales larger than the cell size is asymptotically equivalent to decrease the ejection rate related to particle inertia. This gives evidence that there exists a scaling in the limit of large observation scale and small response time of the particles, even if the flow has no scale invariance.
There are several extensions that need to be investigated in order to gain from the study of such models a more quantitative information on the distribution of particles in real flows. The most significant improvement is to give a spatial structure to the fluid velocity. This can be done by introducing a spatial correlation between the vorticities of cells. Preliminary investigations suggest that such a modified model could be approached by taking its continuum limit. Another effect that may be worth taking into consideration is random sweeping of structures by the fluid flow. We assumed that the eddies are frozzen (and occupy the same cell) during their whole life time. The model could be extended by adding to the dynamics random hops between cells of these structures. Another extension could consist in investigating in a systematic manner three-dimensional versions of the model. As stated above, many statistical quantities may depend on how ejection from rotating regions is implemented in three dimensions. Finally, it is worth mentioning again that the main advantage of such models is to give a heuristic understanding of the relations between the properties of the fluid velocity field and the mass distribution of particles. This first step is necessary in the development of a phenomenological framework for describing the spatial distribution of heavy particles in turbulent flows. This would allow for using Kolmogorov 1941 dimensional arguments to understand how the particle dynamical properties depend on scale. Moreover, such a framework could be used to obtain refined predictions accounting for the effect of the fluid flow intermittency and to describe the dependence upon the Reynolds number of the spatial distribution of particles.