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

J. Bec, Raphael Chétrite

► To cite this version:

J. Bec, Raphael Chétrite. Toward a phenomenological approach to the clustering of heavy particles in turbulent flows. New Journal of Physics, 2007, 9, pp.77. 10.1088/1367-2630/9/3/077. [hal-00388079](https://hal.science/hal-00388079v1)

HAL Id: hal-00388079

https://hal.science/hal-00388079v1

Submitted on 17 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
OPEN ACCESS

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

To cite this article: Jérémie Bec and Raphaël Chétrite 2007 *New J. Phys.* **9** 77

View the [article online](#) for updates and enhancements.

## Related content

- *Diffusion in time-dependent random environments: mass fluctuations and scaling properties*
Giorgio Krstulovic, Rehab Bitane and Jérémie Bec

- *Superdiffusion of massive particles induced by multi-scale velocity fields*
R. M. Baram, P. G. Lind, J. S. Andrade jr. et al.

- *The Eulerian description of dilute collisionless suspension*
G. Boffetta, A. Celani, F. De Lillo et al.

## Recent citations

- *Push and pull: attractors and repellors of a dynamical system can localize inertial particles*
Guy Metcalfe

- *Lacunarity exponents*
Michael Wilkinson et al

- *Scale dependence of multiplier distributions for particle concentration, enstrophy, and dissipation in the inertial range of homogeneous turbulence*
Thomas Hartlep et al
# Toward a phenomenological approach to the clustering of heavy particles in turbulent flows

Jérémie Bec<fnref n="1" />,<fnref n="3" /> and Raphaël Chétrite<fnref n="2" />

New Journal of Physics 9 (2007) 77

Received 17 January 2007

Published 30 March 2007

Online at http://www.njp.org/

doi:10.1088/1367-2630/9/3/077

**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 is a measure of 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 a power-law behaviour 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. Contrary to what has been observed in direct numerical simulations of turbulent flows (Bec *J et al* 2007 *Phys. Rev. Lett.* **98** 084502), 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.

<footnote n="1">¹ Laboratoire Cassiopée, CNRS, OCA; BP4229, 06304 Nice Cedex4, France</footnote>

<footnote n="2">² Laboratoire de Physique, ENS Lyon, 46 Allée d'Italie, 69007 Lyon, France
E-mail: jeremie.bec@obs-nice.fr</footnote>

<footnote n="3">³ Author to whom any correspondence should be addressed.</footnote>
Contents

<table><tr><td>1. Introduction</td><td>2</td></tr><tr><td>2. Ejection of heavy particles from eddies</td><td>4</td></tr><tr><td>3. A simple mass transport model</td><td>5</td></tr><tr><td>4. Distribution of mass</td><td>8</td></tr><tr><td>5. Coarse-grained mass distribution</td><td>12</td></tr><tr><td>6. Conclusions</td><td>15</td></tr><tr><td>Acknowledgment</td><td>15</td></tr><tr><td>References</td><td>15</td></tr></table>

1. 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 modelled 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

$$ \frac{d^2 X}{dt^2} = -\frac{1}{\tau} \left[ \frac{d X}{dt} - u(X, t) \right], \qquad (1) $$

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 $\delta$-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 $\delta$-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]–[15].

Recent numerical studies in fully developed turbulent flows [16] showed that the spatial distribution of particles at length scales within the inertial range is 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. These numerical investigations confirmed a prediction made in [17] regarding the presence of such a power-law behaviour at small masses, which was based on studying the growth rate of mass moments in the limit of small particle response times. To understand at least from a qualitative and phenomenological viewpoint how voids form and how they affect the statistics of the mass distribution, 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 $\omega$ 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 $\omega$ toward the neighbouring cells where the vorticity vanishes. The amount of mass escaping to neighbours is at most a fixed fraction $\gamma$ 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 $\gamma$ and the response time $\tau$ of the heavy particles. Section 3 consists of 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.
## 2. Ejection of heavy particles from eddies

The goal of this section is to give some phenomenological arguments explaining why the model which is briefly 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 $\gamma$ and the response time $\tau$ entering the dynamics of heavy particles. For this we focus on the two-dimensional case and consider a cell of size $\ell$ where the fluid vorticity $\omega$ is constant and the fluid velocity vanishes at the centre 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

$$ \frac{d^2 X}{dt^2} = -\frac{1}{\tau} \frac{dX}{dt} + \frac{\omega}{\tau} \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix} X. \qquad (2) $$

It is easily checked that the four eigenvalues of the evolution matrix are the following complex conjugate

$$ \lambda_{\pm, \pm} = \frac{-1 \pm \sqrt{1 \pm 4i\tau\omega}}{2\tau}. \qquad (3) $$

Only $\lambda_{+, -}$ and $\lambda_{+, +}$ have a positive real part which is equal to

$$ \mu = \frac{-1 + (1/2)\sqrt{2\sqrt{1 + 16\tau^2\omega^2} + 2}}{2\tau}. \qquad (4) $$

This means that the distance of the particles to the center of the cell increases exponentially fast in time with a rate $\mu$. 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\mu$. Namely the mass of particles which are still in the cell at time $T$ is

$$ m(T) = m(0)(1 - \gamma) = m(0) \exp \left[ -\frac{T}{\tau} \left( -1 + \frac{1}{2} \sqrt{2\sqrt{1 + 16\tau^2\omega^2} + 2} \right) \right]. \qquad (5) $$

The rate $\gamma$ at which particles are expelled from the cell depends upon the response time $\tau$ 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 $\omega$. The second time scale is the turnover $\omega^{-1}$ which measures the strength of the eddy. There are hence two dimensionless parameters entering the ejection rate $\gamma$: the Stokes number $St = \tau\omega$ giving a measure of inertia and the Kubo number $Ku = T\omega$ which is the ratio between the correlation time of structures and their eddy turnover time. We hence obtain the following estimate of the ejection rate

$$ \gamma = 1 - \exp \left[ -\frac{Ku}{St} \left( -1 + \frac{1}{2} \sqrt{2\sqrt{1 + 16St^2} + 2} \right) \right]. \qquad (6) $$
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Fraction of the mass of particles that are uniformly distributed and ejected from an eddy of arbitrary size $\ell$ as a function of the Stokes number $St = \tau\omega$. The various curves refer to different values of the Kubo number $Ku$ as labelled.</figcaption></figure>

The graph of the fraction of particles ejected from the cell as a function of the Stokes number is given in figure 1 for three different values of the Kubo number. The function goes to zero as $Ku St$ in the limit $St \to 0$ and as $Ku St^{-1/2}$ in the limit $St \to \infty$. It reaches a maximum which is an indication of a maximal segregation of the particles, for $St \approx 1.03$ independently of the value of $Ku$.

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 $\omega$. There are however two main difficulties. The first is that in three dimensions the eigenvalues of the strain matrix in rotating regions are no longer 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 $\omega$. 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.

### 3. 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 cells 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; $\Omega_j = 1$ with probability $p$ and $\Omega_j = 0$ with probability $1-p$. The
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Sketch of the dynamics in the one-dimensional case: the fluxes of mass are represented as arrows. A cross means no flux.</figcaption></figure>

evolution of mass between times $n$ and $n + 1$ is given by:

$$m_j(n + 1) = \begin{cases} m_j(n) - \frac{\gamma}{2} [2 - \Omega_{j-1} - \Omega_{j+1}] m_j(n) & \text{if } \Omega_j = 1, \\ m_j(n) + \frac{\gamma}{2} [\Omega_{j-1} m_{j-1}(n) + \Omega_{j+1} m_{j+1}(n)] & \text{if } \Omega_j = 0. \end{cases} \quad (7)$$

In other terms, when $\Omega_j = 1$, the $j$th cell loses mass if $\Omega_{j-1} = 0$ or $\Omega_{j+1} = 0$, and when $\Omega_j = 0$, it gains mass if $\Omega_{j-1} = 1$ or $\Omega_{j+1} = 1$. The flux of mass between the $j$th and the $(j+1)$th cell is proportional $\Omega_j - \Omega_{j+1}$ (see figure 2). In particular, if $\Omega_j = \Omega_{j+1}$, no mass is transferred 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 $\sum_j m_j = N$. Spatial homogeneity of the random process $\Omega_j$ implies that $\langle m_j \rangle = 1$ for all later times, where the angular brackets denote average with respect to the realizations of the $\Omega_j$'s.

A noticeable advantage of such a model for mass transportation is that the mass field $m = (m_1, \dots, 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

$$\begin{aligned} p_N(m, n+1) &= \int d^N m' p_N(m', n) P[m' \rightarrow m] \\ &= \int d^N m' p_N(m', n) \int d^N \Omega p(\Omega) P[m' \rightarrow m | \Omega], \end{aligned} \quad (8)$$

where $P[m' \rightarrow m | \Omega]$ denotes the transition probability from the field $m'$ to the field $m$ conditioned on the realization of $\Omega = (\Omega_1, \dots, \Omega_N)$. In our case it takes the form

$$P[m' \rightarrow m | \Omega] = \prod_{j=1}^{N} \delta[m_j - (m'_j + \mu_{j-1}(n) - \mu_j(n))]. \quad (9)$$

The variable $\mu_j$ denotes here the flux of mass between the $j$th and the $(j+1)$th cell. It is a function of $\Omega_j$, $\Omega_{j+1}$, and of the mass contained in the two cells. It can be written as

$$\mu_j(n) = \frac{\gamma}{2} [\Omega_j(n)(1 - \Omega_{j+1}(n))m'_j(n) - \Omega_{j+1}(n)(1 - \Omega_j(n))m'_{j+1}(n)]. \quad (10)$$
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Left panel: sketch of the ejection model in two dimensions; the rotating cells are designated by small eddies; the flux of mass (blue arrows) is from the rotating cells to those without any eddy. Right panel: snapshot of the distribution of mass in the statistically steady regime for a domain of $50^2$ cells with $p = 0.75$; white squares are almost empty cells and the darker regions correspond to cells where the mass is concentrated.</figcaption></figure>

In the particular case we are considering, the joint probability of the $\Omega_j$'s factorizes and we have

$$p(\Omega_j) = p\delta(\Omega_j - 1) + (1 - p)\delta(\Omega_j), \quad (11)$$

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 neighbours (see figure 3 left panel). 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\le i,j\le N}$. The transition probability reads in that case

$$P[M' \to M|\Omega] = \prod_{i=1}^{N} \prod_{j=1}^{N} \delta[m_{i,j} - (m'_{i,j} + \mu_{i-1,j}^{(1)} - \mu_{i,j}^{(1)}) + \mu_{i,j-1}^{(2)} - \mu_{i,j}^{(2)}]. \quad (12)$$

where the fluxes now take the form

$$\mu_{i,j}^{(1)} = \frac{\gamma}{4} [\Omega_{i,j}(1 - \Omega_{i+1,j})m'_{i,j} - \Omega_{i+1,j}(1 - \Omega_{i,j})m'_{i+1,j}], \quad (13)$$

and $\mu_{i,j}^{(2)}$ defined by inverting $i$ and $j$ in the definition of $\mu_{i,j}^{(1)}$.

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 one corresponding to a uniform distribution; strong deviations at small masses can be qualitatively observed (see figure 3 right panel).

The model can be easily generalized to arbitrary space dimensions. 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.
<table><caption>**Table 1.** Enumeration of all possible configurations of the spin vorticity Ω in three neighbouring cells, together with their probabilities and the associated mass fluxes.</caption><thead><tr><th>Ω<sub>j−1</sub></th><th>Ω<sub>j</sub></th><th>Ω<sub>j+1</sub></th><th>Probability</th><th>μ<sub>j−1</sub></th><th>μ<sub>j</sub></th></tr></thead><tbody><tr><td>0</td><td>0</td><td>0</td><td>(1 − p)<sup>3</sup></td><td>0</td><td>0</td></tr><tr><td>0</td><td>0</td><td>1</td><td>p(1 − p)<sup>2</sup></td><td>0</td><td>−γm'<sub>j+1</sub>/2</td></tr><tr><td>0</td><td>1</td><td>0</td><td>p(1 − p)<sup>2</sup></td><td>−γm'<sub>j</sub>/2</td><td>γm'<sub>j</sub>/2</td></tr><tr><td>0</td><td>1</td><td>1</td><td>p<sup>2</sup>(1 − p)</td><td>−γm'<sub>j</sub>/2</td><td>0</td></tr><tr><td>1</td><td>0</td><td>0</td><td>p(1 − p)<sup>2</sup></td><td>γm'<sub>j−1</sub>/2</td><td>0</td></tr><tr><td>1</td><td>0</td><td>1</td><td>p<sup>2</sup>(1 − p)</td><td>γm'<sub>j−1</sub>/2</td><td>−γm'<sub>j+1</sub>/2</td></tr><tr><td>1</td><td>1</td><td>0</td><td>p<sup>2</sup>(1 − p)</td><td>0</td><td>γm'<sub>j</sub>/2</td></tr><tr><td>1</td><td>1</td><td>1</td><td>p<sup>3</sup></td><td>0</td><td>0</td></tr></tbody></table>

## 4. 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$

$$p_1(m_j) = \int dm'_{j-1} dm'_j dm'_{j+1} p_3(m'_{j-1}, m'_j, m'_{j+1}) \int d\Omega_{j-1} d\Omega_j d\Omega_{j+1} \times p(\Omega_{j-1}) p(\Omega_j) p(\Omega_{j+1}) \\ \times \delta[m_j - (m'_j + \mu_{j-1} - \mu_j)]. \quad (14)$$

We then explicit all possible fluxes, together with their probabilities, by considering all possible configurations of the spin vorticity triplet $(\Omega_{j-1}, \Omega_j, \Omega_{j+1})$. The results are summarized in table 1. This leads to rewriting the one-point PDF as

$$p_1(m) = [p^3 + (1-p)^3]p_1(m) + \frac{2p^2(1-p)}{1-\gamma/2}p_1\left(\frac{m}{1-\gamma/2}\right) + \frac{p(1-p)^2}{1-\gamma}p_1\left(\frac{m}{1-\gamma}\right) \\ +2p(1-p)^2 \int_0^{2m/\gamma} dm' p_2(m', m - \frac{\gamma}{2}m') \\ +p^2(1-p) \int_0^{2m/\gamma} dm' \int_0^{2m/\gamma-m'} dm'' p_3(m', m - \frac{\gamma}{2}(m' + m''), m''). \quad (15)$$

The first term on the right-hand side comes from realizations with no flux. The second term is ejection to one neighbour and the third to two neighbouring cells. The fourth term involving an average over the two-cell distribution is related to events when mass is transferred from a single neighbour to the considered cell. Finally, the last term accounts for transfers from the two direct neighbours. Note that, in order to write (15), we made use of the fact that $p_2(x, y) = p_2(y, x)$.

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 $\gamma$ 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 PDFs display an algebraic
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Log-log plot of the one-point PDF of the mass in one dimension for $p = 1/2$ and different values of the parameter $\gamma$ as labelled. The integration was done on a domain of $2^{16} = 65536$ cells and time average we performed during $10^6$ time steps after a statistical steady state is reached.</figcaption></figure>

behaviour 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

$$P^<(m) = \mathrm{Prob}(m_j(n) < m) \geq \mathrm{Prob}(A), \qquad (16)$$

where $A$ is a set of space-time realizations of $\Omega$ 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 one and its two neighbours have zero. The mass at time $n$ is related to the mass at time $n - N$ by

$$m_j(n) = (1 - \gamma)^N m_j(n - N), \quad \text{that is} \quad N = \frac{\log[m_j(n - N)/m_j(n)]}{\log(1 - \gamma)}. \qquad (17)$$

The probability of such a realization is clearly $p^N (1 - p)^{2N}$. Replacing $N$ by the expression obtained in (17), we see that

$$\mathrm{Prob}(A) = \left[ \frac{m_j(n)}{m_j(n - N)} \right]^\beta \quad \text{with} \quad \beta = \frac{\log[p(1 - p)^2]}{\log(1 - \gamma)}. \qquad (18)$$
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Scaling exponent $\alpha$ as a function of the ejection rate $\gamma$ for three different values of $p$ as labelled. The solid lines represents the prediction given by (21); the error bars are estimated from the maximal deviation of the logarithmic derivative from the estimated value. Inset: difference between the numerical estimation and the value predicted by theory.</figcaption></figure>

After averaging with respect to the initial mass $m_j(n - N)$, one finally obtains

$$P^<(m) \geq Am^{\beta}. \qquad (19)$$

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 \to 0$, i.e. that $p_1(m) \simeq Cm^{\alpha}$. To obtain how the exponent $\alpha$ behaves as a function of the parameters $\gamma$ and $p$, this ansatz is injected in the stationary version of the Markov equation (15). One expects that the small-mass behaviour involves only the terms due to ejection from a cell, namely the three first terms in the rhs of (15), and that the terms involving averages of the two-point and three-point PDFs give only sub-dominant contributions. This leads to

$$Cm^{\alpha} \approx C[p^3 + (1-p)^3]m^{\alpha} + C\frac{2p^2(1-p)}{1-\gamma/2}\left[\frac{m}{1-\gamma/2}\right]^{\alpha} + C\frac{p(1-p)^2}{1-\gamma}\left[\frac{m}{1-\gamma}\right]^{\alpha}. \qquad (20)$$

Equating the various constants we finally obtain that the exponent $\alpha$ satisfies

$$\frac{2p}{(1 - \gamma/2)^{\alpha+1}} + \frac{(1 - p)}{(1 - \gamma)^{\alpha+1}} = 3. \qquad (21)$$

Note that the actual exponent $\alpha$ given by this relation is different from the lower-bound $\beta + 1$ obtained above in (18) and (19). However it is easily checked that $\alpha$ approaches the lower bound when $p \to 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 $\gamma$ 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 $\gamma = 0.9$ in figure 5. 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\gamma$. Notice that for all values of *p*, one has $\alpha \le 0$ when $\gamma \ge 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 \approx 1$.

It is much less easy to get from numerics the behaviour 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) \propto -m \log m$ when $m \gg 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 from a large number *M* of neighbouring 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 ($\Omega_{j-1}$, $\Omega_j$, $\Omega_{j+1}$) = (0, 1, 1). After *N* time steps with this configuration, the fraction of mass transfered is $1 - (1 - \gamma/2)^N$. This process is then repeated for moving mass to the second neighbour, and so on. After order *M* iterations, the mass in the *M*th neighbour is

$$m = \frac{1 - [1 - (1 - \gamma/2)^N]^M}{(1 - \gamma/2)^N} \qquad (22)$$

This means that

$$M = M(m, N) = \frac{\log[1 - m(1 - \gamma/2)^N]}{\log[1 - (1 - \gamma/2)^N]} \qquad (23)$$

with the condition that $N > -(\log m)/[\log(1 - \gamma/2)]$. The probability of this whole process of mass transfer is

$$\mathcal{P} = [p^2(1-p)]^{NM} = \exp[\log(p^2(1-p))N M(m, N)]. \qquad (24)$$

All the processes of this type will contribute terms in the right tail of the mass PDF. The dominant behaviour 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. $N \gg - (\log m)/[\log(1 - \gamma/2)]$), then $NM(m, N) \gg -m(\log m)/[\log(1 - \gamma/2)]$. On the other hand when *N* is chosen of the order of $\log m$, then $NM(m, N) \propto m \log m$. This suggests that the minimum is attained for $N^* \propto \log m$. Finally, such estimates lead to the prediction that the right tail of the mass PDF behaves as

$$p_1(m) \propto \exp[-C m \log m], \qquad (25)$$

where *C* is a positive constant that depends upon the parameters *p* and *γ*. As seen in figure 6, such a behaviour 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 $\alpha$
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Lin-log plot of the one-point PDF of the mass *m* in one dimension represented as a function of *m* log *m* for *p* = 1/2 and various values of the parameter *γ* as labelled; the different colours and symbols are the same as those used in figure 4. Inset: behaviour of the constant *C* appearing in (25) as a function of the ejection rate *γ* for three different values of the fraction of space *p* occupied by eddies (blue crosses: *p* = 0.1, black times: *p* = 0.5, red circles: *p* = 0.9).</figcaption></figure>

of the algebraic behaviour at small masses is given as a solution of

$$ \frac{4p^3}{(1 - \gamma/4)^{\alpha+1}} + \frac{6p^2(1-p)}{(1 - \gamma/2)^{\alpha+1}} + \frac{4p(1-p)^2}{(1 - 3\gamma/4)^{\alpha+1}} + \frac{(1-p)^3}{(1 - \gamma)^{\alpha+1}} = 5(1 - p + p^2). \quad (26) $$

By arguments which are similar to the one-dimensional case and which are not detailed here, one obtains also that $\log p_1(m) \propto -m \log m$. Numerical experiments in two dimensions confirm these behaviours of the mass probability distribution. As seen from figure 7 an algebraic behaviour of the left tail of the PDF of *m* is observed and the value of the exponent is in good agreement with (26).

## 5. 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

$$ \bar{m}_L = \frac{\ell}{L} \sum_{j=-K}^{K} m_j, \quad \text{where } K = L/2\ell. \qquad (27) $$

As seen from the numerical results presented on figure 8, the functional form of the PDF $p_L(\bar{m})$ is qualitatively similar to that of the mass in a single cell. In particular for various values of *L*
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Log-log plot of the one-point PDF of the mass in two dimensions for $p = 1/2$ and various values of the parameter $\gamma$; the symbols and colour refer to the same values of $\gamma$ as in figure 4. The integration was done on a domain of $1024^2$ cells and time average were performed during $3 \times 10^4$ time steps after a statistical steady state is reached. Inset: exponent $\alpha$ of the algebraic left tail as a function of the ejection rate $\gamma$ for three different values of $p$ (blue crosses: $p = 0.1$, black times: $p = 0.5$, red circles: $p = 0.9$). The solid lines show the analytic values obtained from (26).</figcaption></figure>

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 behaviour of the exponent.

For this, we consider the cumulative probability $P_L^<(\bar{m})$ to have $\bar{m}_L$ smaller than $\bar{m}$. We first observe that in order to have $\bar{m}_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 neighbours. For this the three central cells must have vorticities $(\Omega_{-1}, \Omega_0, \Omega_1) = (0, 1, 0)$, respectively, during $N$ time steps. After that the second step consists in transferring the mass toward the next neighbours; the best realization is then to have during $N$ time steps $(\Omega_{-2}, \Omega_{-1}, \Omega_0, \Omega_1, \Omega_2) = (0, 1, 1, 1, 0)$. The transfer toward neighbours is repeated recursively. At the $j$th step, the best transfer is given by choosing $(\Omega_{-j-1}, \Omega_{-j}, \Omega_{-j+1}) = (0, 1, 1)$ and $(\Omega_{-j-1}, \Omega_j, \Omega_{j+1}) = (1, 1, 0)$ during a time $N$. One can easily check that for large $N$ and after repeating this procedure $K$ times, the mass which remains in the $2K + 1$ cells forming the coarse-grained cell is $\bar{m}_L \simeq (1 - \gamma/2)^N$. The total probability of this process is $[p^4(1-p)^2]^{KN}$, which leads to estimating the cumulative probability of $\bar{m}_L$ as

$$P_L^<(\bar{m}) \propto \bar{m}^{\alpha_L}, \quad \text{with} \quad \alpha_L \approx \frac{L}{2 \log(1 - \gamma/2)} \log[p^4(1-p)^2]. \qquad (28)$$
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Log-log plot of the PDF of the coarse-grained mass $\bar{m}_L$ in one dimension for $p = 2/3$ and various values of $\gamma$ and $L$ associated to three different ratios $\kappa = L/\log(1 - \gamma/2)$ as labelled.</figcaption></figure>

This approach guarantees that the PDF $p_L(\bar{m}) = dP_L^< /d\bar{m}$ of the coarse-grained mass $\bar{m}$ 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(\bar{m})$, one expects a behaviour similar to that obtained in the case of the one-cell mass distribution, namely $\log p_L(\bar{m}) \propto -\bar{m} \log \bar{m}$ for $\bar{m} \gg 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 behaviours of $p_L(\bar{m})$ as a function of $L$ and of the ejection rate $\gamma$. For instance one can argue whether the limits $L \to \infty$ and $\gamma \to 0$ are equivalent. The estimation (28) suggests that the exponent $\alpha_L$ depends only on the ratio $\kappa = L/\log(1 - \gamma/2)$. Note that the limit of small $\gamma$ 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 mass $\bar{m}_L$ for various values of $L$ and $\gamma$ chosen such that the ratio $\kappa$ is 8, 10 or 14. While the left-hand 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.

<figure><img src="image_10.png" /></figure>
## 6. 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 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 frozen (and occupy the same cell) during their whole lifetime. The model could be extended by adding to the dynamics random hops between cells of these structures. Another extension could consist of 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 Kolmogorov 1941 dimensional arguments to be used 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.

## Acknowledgment

This work benefited from useful discussions with M Cencini, D Mitra and S Musacchio who are warmly acknowledged.

## References

[1] Dimotakis P 2005 Turbulent mixing *Annu. Rev. Fluid Mech.* **37** 329–56

[2] Maxey M R and Riley J 1983 Equation of motion of a small rigid sphere in a nonuniform flow *Phys. Fluids* **26** 883–9

[3] Eaton J K and Fessler J R 1994 Preferential concentration of particles by turbulence *Int. J. Multiphase Flow* **20** 169–209
[4] Falkovich G, Fouxon A and Stepanov M 2002 Acceleration of rain initiation by cloud turbulence *Nature* **419** 151–4

[5] Lewis D and Pedley T 2001 The influence of turbulence on plankton predation strategies *J. Theor. Biol.* **210** 347–65

[6] Seinfeld J 1986 *Atmospheric Chemistry and Physics of Air Pollution* (New York: Wiley)

[7] de Pater I and Lissauer J 2001 *Planetary Science* (Cambridge: Cambridge University Press)

[8] Bec J 2005 Multifractal concentrations of inertial particles in smooth random flows *J. Fluid Mech.* **528** 255–77

[9] Maxey M R 1987 The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields *J. Fluid Mech.* **174** 441–65

[10] Kraichnan R H 1968 Small-scale structure of a scalar field convected by turbulence *Phys. Fluids* **11** 945–53

[11] Piterbarg L I 2002 The top Lyapunov exponent for stochastic flow modeling the upper ocean turbulence *SIAM J. Appl. Math.* **62** 777–800

[12] Mehlig B and Wilkinson M 2004 Coagulation by random velocity fields as a Kramers problem *Phys. Rev. Lett.* **92** 250602

[13] Duncan K, Mehlig B, Östlund S and Wilkinson M 2005 Clustering by mixing flows *Phys. Rev. Lett.* **95** 240602

[14] Derevyanko S, Falkovich G, Turitsyn K and Turitsyn S 2006 Explosive growth of inhomogeneities in the distribution of droplets in a turbulent air *Preprint nlin.CD/0602006*

[15] Bec J, Cencini M and Hillerbrand R 2007 Heavy particles in incompressible flows: the large Stokes number asymptotics *Physica D* **226** 11–22

[16] Bec J, Biferale L, Cencini M, Lanotte A, Musacchio S and Toschi F 2007 Heavy particle concentration in turbulence at dissipative and inertial scales *Phys. Rev. Lett.* **98** 084502

[17] Balkovsky E, Falkovich G and Fouxon A 2001 Intermittent distribution of inertial particles in turbulent flows *Phys. Rev. Lett.* **86** 2790–3