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

HAL
open science

# Turbulence, Particles and Magnetic fields

Rémi Zamansky

► To cite this version:

Rémi Zamansky. Turbulence, Particles and Magnetic fields. Fluid Dynamics [physics.flu-dyn]. Institut National Polytechnique de Toulouse, 2023. tel-04488981

HAL Id: tel-04488981

https://hal.science/tel-04488981

Submitted on 4 Mar 2024

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Université de Toulouse

Mémoire présenté en vue de l'obtention de

l'Habilitation à Diriger des Recherches

par

Rémi Zamansky

Maître de Conférences à l'Institut National Polytechnique de Toulouse
Institut de Mécanique des Fluides de Toulouse (IMFT)

Turbulence, particles and magnetic fields

Soutenance : 7 novembre 2023

Composition du jury :

<table><tr><td>Alain Pumir</td><td>Directeur de recherche (CNRS - ENS Lyon)</td><td>Rapporteur</td></tr><tr><td>Olivier Simonin</td><td>Professeur des Universités (INPT - IMFT)</td><td>Rapporteur</td></tr><tr><td>John Christos Vassilicos</td><td>Directeur de recherche (CNRS - LMFL)</td><td>Rapporteur</td></tr><tr><td>Aurore Naso</td><td>Directrice de recherche (CNRS - LMFA)</td><td>Examinatrice</td></tr><tr><td>Joachim Peinke</td><td>Professor (University Oldenburg)</td><td>Examinateur</td></tr><tr><td>Nicolas Rimbert</td><td>Professeur des Universités (ENSEM - LEMTA)</td><td>Examinateur</td></tr><tr><td>Dominique Legendre</td><td>Professeur des Universités (INPT - IMFT)</td><td>Correspondant</td></tr></table>
<figure><img src="image_3.png" /></figure>
# Contents

<table><tr><td>1</td><td>About this manuscript</td><td>5</td></tr><tr><td>2</td><td>Turbulence</td><td>9</td></tr><tr><td>2.1</td><td>Universality, cascade and intermittency</td><td>9</td></tr><tr><td>2.2</td><td>Lagrangian dynamics</td><td>17</td></tr><tr><td>2.3</td><td>Scaling laws of the acceleration</td><td>20</td></tr><tr><td>2.4</td><td>Stochastic modeling of the fluid-particle dynamics</td><td>30</td></tr><tr><td>2.5</td><td>Discussion and perspectives</td><td>44</td></tr><tr><td>3</td><td>Particles and turbulence</td><td>47</td></tr><tr><td>3.1</td><td>The Euler-Lagrange approach</td><td>48</td></tr><tr><td>3.2</td><td>Force and acceleration statistics for particles in turbulence</td><td>49</td></tr><tr><td>3.3</td><td>Stochastic modeling of the dynamics of a particle and application to LES</td><td>66</td></tr><tr><td>3.4</td><td>Turbulence induced by a swarm of rising bubbles</td><td>77</td></tr><tr><td>3.5</td><td>Perspectives</td><td>95</td></tr><tr><td>4</td><td>Magnetic fields, particles and turbulence</td><td>97</td></tr><tr><td>4.1</td><td>Convection of a liquid metal subject to alternating magnetic field</td><td>97</td></tr><tr><td>4.2</td><td>Perturbation of eddy-currents by inclusions in liquid metal</td><td>106</td></tr><tr><td>4.3</td><td>Perspectives</td><td>113</td></tr><tr><td>5</td><td>Bibliography</td><td>115</td></tr></table>
<figure><img src="image_4.png" /></figure>
# 1 About this manuscript

This manuscript describes my recent work related to the Lagrangian modeling of turbulent flows, to the modeling of particle motion in turbulent flows and the generation of turbulence induced by a dispersed phase, and finally on flows involving a liquid metal with non-conducting particles subjected to an alternating magnetic field. It is clear that the content of the manuscript is not as general as the title might suggest, and this title rather gives the keywords associated with the work presented in this dissertation.

In the first chapter, after an introduction on energy cascade of turbulent flow and intermittency, I present recent development on the statistical analysis of the material derivative of the velocity in isotropic turbulence and on the stochastic modeling of the fluid particles dynamics. Obtaining numerical solutions of the Navier-Stokes equation with supercomputers, I highlight a relationship between the force on a fluid particle (the acceleration), the kinetic energy and the dissipated power. This relation connects the idea of the energy cascade, and its intrinsic multi-scale character to the dynamics of the tracer in turbulent flows. By relying on this relation, I propose a vectorial stochastic model that effectively describes the interaction of fluid particles with all other fluid particles of the flow. This model presents some of the remarkable features of the fluid particle dynamics, in particular the occurrence of extreme events, the temporal asymmetry and the emergence of anomalous scaling law.

The second chapter starts with a discussion on the Euler-Lagrange approach used to model dispersed phase flows and on the coarse graining implied by this framework. Then, statistical descriptions of the acceleration and hydrodynamic forces of particles that are heavy or light, small or large is presented. Subsequently, we consider stochastic modeling of the high frequency part of the dynamics of these particles, as well as its coupling with the large eddy simulation approaches. Finally, the last part of the chapter focus on the numerical modeling of the turbulent agitation induced by the rising of a swarm of bubbles, and to the analysis of the specific properties of this flow resulting from the interactions between wakes.

In the last chapter, we discuss our work on the magnetohydrodynamic of liquid metal subjected to an alternating magnetic field. In the first part, I present the experimental setup and the numerical simulations developed to study thermo-magneto-convection flows induced by the Joule heating and the Laplace force. The second part focus on a modeling approach, confronted to an experimental approach, of the perturbation of an alternating magnetic field caused by the presence of bubbles in a liquid metal.

The common thread between these three chapters is the focus on the modeling of the physical phenomena. I find that the 11 successive states of lithography "Le Taureau", realized by Pablo Picasso in 1945, reproduced in the figure 1.1, illustrates quite well the process of proposing an adequate model, as simple and smart as possible.<fnref n="1" /> How many

<footnote n="1">¹Frederic Risso indicated to me that Christophe Clanet had used this same artwork to discuss the role of modeling during a plenary Lecture at the international conference on multiphase flows and Wladimir</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 1.1: Pablo Picasso, The 11 successive states of lithography "Le Taureau", 1945.</figcaption></figure>

pencil strokes are needed to instantly recognize a bull or how many degrees of freedom are needed to describe a system? The answer obviously depends on the level of detail one want to reproduce: Is it an animal? Is it a cow or a bull? What is the color of the animal's coat? Is it wounded?... I reproduce an excerpt from Fernand Mourlot's text, "Gravés dans ma mémoire", Ed. Robert Laffont, 1979, which details Picasso's creative process for his engraving of the Bull, and which underlies the necessity to well grasp the physical behavior of the system under study to propose an adequate model: "L'opération a duré quinze jours. Le 5 décembre 1945, un mois après son arrivée rue de Chabrol, Picasso a dessiné au lavis un taureau. Un taureau magnifique, très bien fait, gentil même. Et puis on lui a donné l'épreuve; nous en avons tiré à peine deux ou trois, ce qui fait que ce taureau est extrêmement rare. Une semaine après, il revient et il demande une nouvelle pierre; il reprend son taureau au lavis et à la plume; il recommence le 18. Troisième état, le taureau est repris au grattage à plat, puis à la plume en accentuant fortement les volumes; le taureau est devenu un animal terrible avec des cornes et des yeux effroyables. Bon, ça n'allait pas, Picasso éxécute un quatrième état, le 22 décembre, et un cinquième, le 24; à chaque fois il simplifie le dessin qui devient de plus en plus géomètrique avec des aplats noirs."

Sixième et septième états, les 26 et 28 décembre, puis, après le retour de Picasso, quatre autres états, onze en tout, les 5, 10 et 17 janvier. Le taureau est réduit à sa plus simple expression; quelques traits d'une maîtrise exceptionnelle qui symbolisent comme un jeu de signes ce malheureux taureau avec sa petite tête d'épingle et ses cornes ridicules en forme d'antenne. Les ouvriers se désolaient d'avoir vu un taureau aussi magnifique transformé.

<footnote>Bergez also reported to me that Lydéric Bocquet had also shown it to illustrate the reduction of complexity during his inaugural lesson at the Collège de France. Although it is not as original as I thought, I left it because I think it has its place here.</footnote>
en une espèce de fourmi...

... C'est Célestin qui a trouvé le mot de la fin : "Picasso, il a fini par là où, normalement, il aurait dû commencer." C'est vrai, seulement pour arriver à son taureau d'une seule ligne, il a fallu qu'il passe par tous les taureaux précédents. Et quand on voit son onzième taureau on ne peut s'imaginer le travail qu'il lui a demandé".
<figure><img src="image_6.png" /></figure>
# 2 Turbulence

## Contents

<table><tr><td>2.1 Universality, cascade and intermittency</td><td>9</td></tr><tr><td>2.2 Lagrangian dynamics</td><td>17</td></tr><tr><td>2.3 Scaling laws of the acceleration</td><td>20</td></tr><tr><td>2.3.1 Methodology</td><td>20</td></tr><tr><td>2.3.2 Conditional statistics given the dissipation and the kinetic energy</td><td>20</td></tr><tr><td>2.3.3 Similarity of the conditional statistics given the dissipation</td><td>23</td></tr><tr><td>2.3.4 Reynolds number dependence of the unconditional acceleration variance</td><td>27</td></tr><tr><td>2.3.5 Multiplicative cascade for the acceleration</td><td>28</td></tr><tr><td>2.4 Stochastic modeling of the fluid-particle dynamics</td><td>30</td></tr><tr><td>2.4.1 Model formulation</td><td>30</td></tr><tr><td>2.4.2 Parameters and numerical approach</td><td>35</td></tr><tr><td>2.4.3 Results</td><td>35</td></tr><tr><td>2.5 Discussion and perspectives</td><td>44</td></tr></table>

## 2.1 Universality, cascade and intermittency

Turbulence is fascinating with all its eddies and vortices forming, deforming and dancing together. As beautiful is the theoretical analysis of turbulent flows, looking for order and coherence in this apparently erratic motion. Everything indicates that the spatio-temporal evolution of turbulent flows can be described by the Navier-Stokes equation, and in particular the numerical simulations which, in this context, are called “direct numerical simulation” (or DNS). Nevertheless, bringing into play nonlinear and nonlocal interactions, it still resists complete analytical method, despite having been proposed 200 years ago [190]. So understanding all the details of turbulent flow remains, to say the least, a formidable problem. The situation bears similarity with critical phenomena, in which the interactions among the large numbers of degrees of freedom induce collective organization at macroscopic scales that cannot be figured out from the behavior of individual particles. The connection between the small-scale velocity fluctuations in turbulence and the large-scale flow properties implies scaling relations. The cornerstone of the statistical description of turbulent flows was advanced in the 40’s by Kolmogorov proposing a universal picture of turbulence based on the idea of a cascade of kinetic energy and its dissipation into heat. Completing the Kolmogorov statistical description approach and in particular its bridging with the Navier-Stokes equations, remains one of the great challenges in the study of turbulence.
## Cascade

Following the idea of the energy cascade due to Richardson [221], Taylor [250, 251] and Kolmogorov [126, 125, 124], energy is injected at some length scale $L$ where it causes vortices of size $L$. As these large vortices do not dissipate energy, structures on smaller scales must form. It is convenient to imagine these new length scales appearing as the large vortices break up into smaller ones, giving rise to eddies of all possible sizes, at least in an interval between $L$ and $\eta$, the dissipative scale. The vortices at a given scale $\ell$ are supplied with energy by the larger ones, at the same time they lose energy by disintegrating themselves, after a time $\tau_\ell$ into smaller vortices. The key assumption here is that the viscosity is not relevant for this process, and so if no dissipation occurs the energy transfer rate (per unit of mass) $\epsilon_\ell$ is constant across scales in a steady regime. Assuming that the eddies of size $\ell$ are characterized be a velocity $u_\ell$, we find scale invariance in the energy cascade:

$$ \frac{u_{\ell}^{2}}{\tau_{\ell}} = \frac{u_{\ell}^{3}}{\ell} = \epsilon_{\ell} = \text{cst} = \langle \epsilon \rangle . \qquad (2.1) $$

This inertial cascade takes place in a range of scales in which the viscosity has negligible effect, that is to say as long as $Re_\ell \gg 1$, where $Re_\ell = u_\ell \ell / \nu$ the Reynolds number of the flow at scale $\ell$. From this Reynolds number, one defines the dissipative length scale $\eta = \nu^{3/4} \langle \epsilon \rangle^{-1/4}$ which correspond to $Re_\eta = 1$, corresponding to the end of the inviscid cascade. And by conservation, in the original Kolmogorov theory, the characteristic energy transfer rate is taken to be the average rate of dissipation of kinetic energy into heat $\langle \epsilon \rangle$. This cascade scaling finds a remarkable confirmation in the “4/5 law”, one of the few exact results in turbulence, which relates the third moment of the velocity increments $(\delta_r u = u(x+r) - u(x))$ to the average dissipation rate [125]: $\langle (\delta_r u)^3 \rangle = -\frac{4}{5} \langle \epsilon \rangle r + \dots$. For the second moment of the velocity increments, Kolmogorov also proposed an empirical relation derived from similarity hypotheses [126]. From dimensional analysis we have

$$ \langle (\delta_r u)^2 \rangle = \text{fct}(r, \nu, \langle \epsilon \rangle, L, \dots) = (\langle \epsilon \rangle r)^{2/3} \Phi(r/L, r/\eta) . \qquad (2.2) $$

The first hypothesis proposed that for high enough Reynolds numbers, turbulence is universal. For $r \ll L$, it is assumed that the velocity increment statistics are independent of the specific forcing mechanism which can be expressed mathematically as $\lim_{r/L \to 0} \Phi(r/L, r/\eta) = \phi(r/\eta)$. And in the second hypothesis, it is proposed that for $r \gg \eta$ velocity increment statistics become independent of the viscosity ($\lim_{r/\eta \to \infty} \phi(r/\eta) = \text{cst}$). With these two hypotheses, which correspond to a complete similarity hypothesis using the Barenblatt terms [14], the variance of the velocity increments scales like

$$ \langle (\delta_r u)^2 \rangle \sim (\langle \epsilon \rangle r)^{2/3} \qquad (2.3) $$

for $\eta \ll r \ll L$. Although this relation is not derived from the Navier-Stokes equation, it received many experimental confirmations.

## Intermittency

Actually with similar arguments, one could generalize the previous scaling relation (2.3) to higher moments of the velocity increment and obtain: $\langle |\delta_r u|^p \rangle \sim (\langle \epsilon \rangle r)^{p/3} \sim (\langle (\delta_r u)^2 \rangle)^{p/2}$. Such a scaling relation would imply that $\delta_r u$ presents a Gaussian distribution with variance $\sigma^2 \sim (\langle \epsilon \rangle r)^{2/3}$. The problem is that although for very large separation $r$, $\delta_r u$ presents a
Gaussian distribution, experiments and numerical simulations have reported systematic deviations from the Gaussian distribution when $r$ is reduced [7, 49, 111]. This is known as intermittency, because on small scales turbulent flows present weak fluctuations with intermittent bursts of much larger fluctuations. And indicates that the large scales of the flow keep influencing small-scale motion of the flow, or equivalently, that the structure of the flow on small scales presents a Reynolds number dependence.

This observation is directly related to the Landau's remark [138], which pointed out that the 1941 picture of a universal inertial range depending only on the mean energy dissipation rate is flawed, because the dissipation rate presents very large, intermittent, spatiotemporal fluctuations. Therefore, the global average of $\epsilon$ is not the relevant scale, since locally it can take values hundreds of times smaller or larger than $\langle \epsilon \rangle$.

These remarks led Kolmogorov and Oboukhov to the refined similarity hypothesis by introducing a local scale for the energy transfer rate based on a locally averaged dissipation rate on a sphere of radius $r$ [128, 192]

$$ \epsilon_r(\mathbf{x}) = \frac{3}{4\pi r^3} \int_{|h|<r} \epsilon(\mathbf{x} + \mathbf{h})d\mathbf{h}. \qquad (2.4) $$

Keeping the idea of the central role of the transfer rate, a local velocity scale is derived $v_r = (\epsilon_r r)^{1/3}$ as well as a local Reynolds number $Re_r = v_r r/\nu$.

Using these local scales, they proposed local similarity hypothesis, that states that $\delta_r u/v_r$ presents universal distribution, only dependent on $Re_r$. Further, it was proposed that for $Re_r \gg 1$ the distribution of $\delta_r u/v_r$ does not dependent on $Re_r$. This translates in the following scaling relation for the structure function when $Re_r \gg 1$:

$$ \langle (\delta_r u)^p \rangle \sim r^{p/3} \langle \epsilon_r^{p/3} \rangle . \qquad (2.5) $$

Further assuming a log-normal distribution of $\epsilon_r$, with parameters that are scale-dependent, for $\eta \ll r \ll L$ the structure function can be expressed as:

$$ \langle (\delta_r u)^p \rangle = c_p (r \langle \epsilon \rangle)^{p/3} (L/r)^{\alpha_p} \qquad (2.6) $$

where $c_p$ is a prefactor that depends on $p$ and possibly on the large-scale features of the flow, and $\alpha_p$ is called intermittency exponent or anomalous exponent by analogy with so-called critical points in phase transitions [90, 73, 16]. The later characterizes the absence of scale separation that manifests as an influence of the large scales on the small ones through the intermittency. This refined scaling of the velocity increments and the intermittency exponent has been studied in much detail, in particular within the framework of multifractal analysis [129, 79, 188, 49, 83, 262, 243].

With these refined hypotheses, the structure of the velocity field loses its universal behavior as it depends on the large scales through the long-range correlation of $\epsilon$ which is set by the specific forcing mechanism. However, in this view, what remains universal is the cascade mechanism itself, which, as fragmentation processes, presents scale-to-scale transformations with local similarities [51, 50, 48, 83, 37, 93, 235].

To summarize, the small-scale quantities, typically the velocity derivatives, present a Reynolds number dependence<fnref n="1" />. They primarily depend on the local dissipation rate whose statistics depend on the Reynolds number through the intermittency of the cascade (the

<footnote n="1">¹Taking $r = \eta$ in (2.6) we have $\langle (\partial_x u_x)^p \rangle \approx \langle (\delta_r u)^p / \eta^p \rangle \approx c_p (\langle \epsilon \rangle / \nu)^{p/2} Re^{3\alpha_p / 4}$ with $L/\eta = Re^{3/4}$.</footnote>
depth of the cascade varies like ln Re$_{λ}$, as discussed in the next subsection). This scaling theory has been proposed in an Eulerian framework as it focuses on the length scale and the question of the dynamics (time evolution) is not a central aspect of this description. Conversely, the Lagrangian description of the turbulence naturally focuses on the dynamics and could be a bridge between Navier-Stokes theory and cascade scaling. More specifically, the acceleration of fluid particles (i.e. the material derivative of the velocity) is classically considered as a small-scale quantity [178, §21.5], meaning that it is determined by the dissipative scales of the turbulence. Indeed, it is intimately linked to the dissipation rate of turbulence. However, various experimental and numerical works have also reported its dependence on the kinetic energy (a quantity carried by the large scales). This is generally attributed to long-range correlations, or intermittency, of the dissipation rate, as discussed above. Below, in sections 2.2-2.3, we show that in addition to this intermittency effect, there is a direct effect of the large scales of the flow structure on the acceleration. This questions the pertinence of considering the acceleration as a small-scale quantity only, and open the possibility of modeling the dynamics of a tracer in a turbulent flow as presented in section 2.4.

Before presenting this aspect of the Lagrangian dynamics, we briefly present the multiplicative cascade model for the energy transfers rate, and the stochastic modeling of the dissipation rate in the Lagrangian framework, which has been proposed consistently with the cascade picture. Note that the reminder of this chapter is largely based on the paper [274] published last year, with just minor additions and corrections<fnref n="2" />.

## Dissipation as multiplicative cascade process

The image of the energy cascade is naturally associated with multiplicative processes [191, 270, 127, 178, 157, 25, 80]. Such a model proposes to express the locally averaged dissipation $\epsilon_n$ over a volume of size $\ell = L\lambda^n$, with L the large-scale of the flow and $\lambda < 1$, as:

$$ \epsilon_n = \epsilon_0 \frac{\epsilon_1}{\epsilon_0} \dots \frac{\epsilon_n}{\epsilon_{n-1}} = \epsilon_0 \prod_{i=1}^{n} \xi_i . \qquad (2.7) $$

Assuming that $\xi_i = \epsilon_i/\epsilon_{i-1}$ are independent positive random numbers with identical distribution across scales we write:

$$ \ln \frac{\epsilon_n}{\epsilon_0} = \sum_{i=1}^{n} \ln \xi_i . \qquad (2.8) $$

Therefore according to the central limit theorem the term on the right tends to a normal distribution with parameters $\mu = n\mu_\xi$ and $\sigma^2 = n\sigma_\xi^2$. The parameters $\mu_\xi$ and $\sigma_\xi^2$ appear as fundamental unknowns<fnref n="3" />. Setting $\ell = \eta$ (i.e. $n = \ln(\eta/L)/\ln \lambda \sim \ln Re_\lambda$) we obtain a model for the local dissipation rate. The log-normal distribution for $\epsilon$ has been confirmed for example by DNS of [271] or experiments [188] and is shown in figure 2.1.

Moreover the variance of the logarithm of the local dissipation rate is then given by $\sigma^2 = \frac{\sigma_\xi^2}{\ln \lambda} \ln \eta/L = A + B \ln Re_\lambda$ as predicted by Kolmogorov and Oboukhov [128, 192].

<footnote n="2">²Actually it is the preparation of this chapter of the memoir that initiated the work published in Ref. [274].</footnote>

<footnote n="3">³Nevertheless, they are related by this relation $\mu_\xi = -\sigma_\xi^2/2$ obtained from the moments of a log-normal variable in order to guarantee that the average energy flux is conserved throughout the cascade.</footnote>
<figure><img src="image_7.png" /><figcaption>Figure 2.1: PDF of pseudodissipation (solid line) (left) and log of pseudodissipation (solid line) (right) centered and normalized by $\mu = \log(\varepsilon/\langle\varepsilon\rangle) = -\sigma^2/2$ and $\sigma^2 = \langle\log^2\varepsilon\rangle - \langle\log\varepsilon\rangle^2$ for $Re_\lambda = 50, 90, 150, 230$ and 380 from orange to black. Comparison with the actual dissipation for $Re_\lambda = 380$ in black dashed line and with $\exp(-x^2/2)/\sqrt{2\pi}$ in gray dashed lines on the right panel.</figcaption></figure>

Such evolution for $\sigma^2$ has been also confirmed by the DNS of Ref. [271] showing that $\sigma^2 \approx 3/8 \ln Re_\lambda/10$. Note that alternative Reynolds number dependence have been proposed in Ref. [49].

<figure><img src="image_8.png" /><figcaption>Figure 2.2: Schematic representation of multiplicative cascades (image from geidav.wordpress.com)</figcaption></figure>

Implicitly the discrete model presented above assumes that the scale ratio $\lambda$ is large enough, so that successive values of the scaling factor $\xi_i$ can be considered as statistically independent. A continuous description of the evolution of energy transfer rate $\varepsilon_r$ across scales, i.e. for infinitesimal $\lambda$, by a Markovian stochastic process has been proposed in Ref. [189], or similarly for the evolution of the velocity increments through scale $r$ in Ref. [81].

Such a multiplicative process, that assumed that the coupling is local in scale, also implies a logarithmic evolution of the spatial correlation of the dissipation rate as explained by Mandelbrot [156]. To show that we consider the transfer rate at two points A and B, $\varepsilon_n^A$ and $\varepsilon_n^B$, both defined on the same scale $n$. The points A and B are separated by a distance $r$ (with $r$ in the inertial range $L > r > \eta$) from each other and we note $k = \ln(r/L)/\ln \lambda$, then $0 < k < n$. Clearly, as illustrated in the figure 2.2, the greater the distance between
the two points, the larger the scale of their common root in the cascade:

$$ \varepsilon_n^A = \varepsilon_0^{AB} \frac{\varepsilon_1^{AB}}{\varepsilon_0^{AB}} \cdots \frac{\varepsilon_k^{AB}}{\varepsilon_{k-1}^{AB}} \frac{\varepsilon_{k+1}^A}{\varepsilon_k^A} \cdots \frac{\varepsilon_n^A}{\varepsilon_{n-1}^A} \quad (2.9) $$

$$ \varepsilon_n^B = \varepsilon_0^{AB} \frac{\varepsilon_1^{AB}}{\varepsilon_0^{AB}} \cdots \frac{\varepsilon_k^{AB}}{\varepsilon_{k-1}^{AB}} \frac{\varepsilon_{k+1}^B}{\varepsilon_k^B} \cdots \frac{\varepsilon_n^B}{\varepsilon_{n-1}^B}. \quad (2.10) $$

In the two previous equations, we have distinguished by the exponents A and B the variables which are specific to points A and B and by AB those which are common. This can be expressed as:

$$ \ln \frac{\varepsilon_n^A}{\varepsilon_0} = \sum_{i=1}^{k} \ln \xi_i^{AB} + \sum_{i=k+1}^{n} \ln \xi_i^A \quad (2.11) $$

$$ \ln \frac{\varepsilon_n^B}{\varepsilon_0} = \sum_{i=1}^{k} \ln \xi_i^{AB} + \sum_{i=k+1}^{n} \ln \xi_i^B. \quad (2.12) $$

The correlation between $\ln \varepsilon_n^A$ and $\ln \varepsilon_n^B$ is defined as

$$ R_{\ln \varepsilon}(r) = \langle (\ln \frac{\varepsilon_n^A}{\varepsilon_0} - \mu)(\ln \frac{\varepsilon_n^B}{\varepsilon_0} - \mu) \rangle = \langle \ln \frac{\varepsilon_n^A}{\varepsilon_*} \ln \frac{\varepsilon_n^B}{\varepsilon_*} \rangle, \quad (2.13) $$

where we noted $\varepsilon_* = \varepsilon_0 e^{\mu}$. Introducing similarly $\xi_* = e^{\mu\chi}$ and $\xi' = \xi/\xi_*$ we express the correlation as:

$$ \begin{aligned} R_{\ln \varepsilon} &= \langle \sum_{i=1}^{n} (\ln \xi_i^A - \mu_\xi) \sum_{j=1}^{n} (\ln \xi_j^B - \mu_\xi) \rangle = \langle \sum_{i=1}^{n} \ln \xi_i'^A \sum_{j=1}^{n} \ln \xi_j'^B \rangle \\ &= \langle \left( \sum_{i=1}^{k} \ln \xi_i'^{AB} + \sum_{i=k+1}^{n} \ln \xi_i'^A \right) \left( \sum_{j=1}^{k} \ln \xi_j'^{AB} + \sum_{j=k+1}^{n} \ln \xi_j'^B \right) \rangle \\ &= \sum_{i=1}^{k} \sum_{j=1}^{k} \delta_{ij} \sigma_\xi^2 = k \sigma_\xi^2. \end{aligned} \quad (2.14) $$

To obtain this results we used the hypothesis that within the same branch, the events at a given scale are independent of those at another scale, $\langle \ln \xi_i'^{AB} \ln \xi_j'^{AB} \rangle = \delta_{ij} \sigma_\xi^2$, as well as the absence of correlation between branches A and B: $\langle \ln \xi_i'^A \ln \xi_j'^B \rangle = 0$. This gives a logarithmic evolution of the correlation coefficient $\rho_{\ln \varepsilon} = R_{\ln \varepsilon}/\sigma^2$, in the range $\eta < r < L$:

$$ \rho_{\ln \varepsilon} = \frac{\langle \ln \frac{\varepsilon_n^A}{\varepsilon_*} \ln \frac{\varepsilon_n^B}{\varepsilon_*} \rangle}{\langle \ln^2 \frac{\varepsilon_n}{\varepsilon_*} \rangle} = \frac{k}{n} = \frac{\ln L/r}{\ln L/\eta} = 1 - \frac{\ln r/\eta}{\ln L/\eta}. \quad (2.15) $$

## Stochastic modeling of the dissipation

Although not trivial, this result can be transposed for the temporal correlation along a particle path [240, 107, 169, 35]. The logarithmic behavior of the correlation is confirmed by DNS, as it can be seen in Ref. [146] or in figure 2.3, where the evolution of the
<figure><img src="image_9.png" /><figcaption>Figure 2.3: Autocorrelation of the log of the dissipation rate along the path of fluid particles from the database of Refs. [21, 140] at $Re_λ = 420$, comparison with $ρ_{lnε} = 1 - ln(τ/τ_η)/ln(τ_ε/τ_η)$.</figcaption></figure>

Lagrangian correlation of the logarithm of the dissipation is plotted from the DNS data of [21, 140].

It has been proposed to model the dissipation rate as a stochastic multiplicative process. Such a process can be generically expressed as :

$$dε = εΠdt + εΣdW \qquad (2.16)$$

with $dW$ the increment of the Wiener process ($\langle dW \rangle = 0$ and $\langle dW^2 \rangle = dt$) and where $Π$ and $Σ$ are to be determined.

Considering that $ε$ follows a log-normal distribution with parameter $σ^2$ and $μ = -σ^2/2$, we define the standard normal variable $χ$ (Gaussian random variable with zero mean and unit variance) as:

$$\frac{ε}{\langle ε \rangle} = \exp \left( \sigmaχ - \frac{σ^2}{2} \right) . \qquad (2.17)$$

A stochastic process for $χ$ has to be given in order to obtain the stochastic process for $ε$, via the Ito transformation.

Pope and Chen [203] proposed to obtain $χ$ from an Orstein-Uhnlebbeck process with a characteristic time $τ_ε$:

$$dχ = -\frac{χ}{τ_ε}dt + \sqrt{\frac{2}{τ_ε}}dW . \qquad (2.18)$$

According to the Ito formula, this gives for $Π$ and $Σ$:

$$Π = - \left( \ln \frac{ε}{\langle ε \rangle} - \frac{σ^2}{2} \right) / τ_ε, \qquad (2.19)$$

$$Σ = \sqrt{2σ^2/τ_ε} . \qquad (2.20)$$

This process gives, as expected, a log-normal distribution for $ε$ (normal distribution for $χ$) as well as an exponential decrease of the correlation of $lnε$ with a characteristic time $τ_ε$. This exponential behavior is not consistent with the multiplicative cascade model as
discussed above. It rather corresponds to a direct energy transfer from large to small-scales.

To ensure a logarithmic decorrelation of the dissipation, Chevillard [54] proposed to adapt the Gaussian multiplicative chaos introduced by Mandelbrot [156]. This leads to a multi-fractal model in which the increment of the Wiener process in Eq. (2.18) is replaced by a fractional Gaussian noise:

$$d\chi = -\frac{\chi}{\tau_{\varepsilon}}dt + \frac{1}{\Lambda}dW_{\tau_c}^0, \qquad (2.21)$$

where, $dW_{\tau_c}^0$ is formally a fractional Gaussian noise with a 0 Hurst exponent, regularized at a timescale $\tau_c$, and $\Lambda$ is a normalization factor ensuring unit variance for $\chi$. The value of $\Lambda$ is dependent on the specific regularization of $dW_{\tau_c}^0$. As explained in Ref. [54], this process can be reexpressed as

$$d\chi(t) = \left(-\frac{\chi}{\tau_{\varepsilon}} + \frac{\Gamma}{\Lambda}\right) dt + \frac{1}{\sqrt{\Lambda^2 \tau_c}} dW \qquad (2.22)$$

with $dW$ the increments of a standard Wiener process and $\Gamma$ corresponds to a convolution of the Wiener increments:

$$\Gamma = -\frac{1}{2} \int_{-\infty}^{t} (t - s + \tau_c)^{-3/2} dW(s) \qquad (2.23)$$

where $dW(s)$ are the increments of the same realization of the Wiener process as in Eq. (2.22). In Eq. (2.23), the regularization time $\tau_c$ prevents the divergence of the kernel when $s \to t$. The normalization factor $\Lambda$ is estimated as $\Lambda = \langle X^2 \rangle$ where $X$ obey the stochastic Eq. (2.22) in which $\Lambda$ has been set to 1.

The stochastic process (2.22) gives a logarithmic correlation for $\chi$: $\langle\chi(t)\chi(t-s)\rangle \sim \ln \frac{\tau_{\varepsilon}}{s}$ for $\tau_c \ll s \ll \tau_{\varepsilon}$, as illustrated in the Fig. 2.4.

<figure><img src="image_10.png" /><figcaption>Figure 2.4: Correlation of $\chi$ for various values of $\tau_{\varepsilon}/\tau_c$.</figcaption></figure>

With the Ito transformation, we obtain the process for $\varepsilon$ from (2.22). This gives for $\Pi$ and $\Sigma$:

$$\Pi = \left( -\ln \frac{\varepsilon}{\langle \varepsilon \rangle} + \frac{\sigma^2}{2\Lambda^2} \left( \frac{\tau_{\varepsilon}}{\tau_c} - \Lambda^2 \right) + \frac{\sigma}{\Lambda} \Gamma \tau_{\varepsilon} \right) / \tau_{\varepsilon}, \qquad (2.24)$$

$$\Sigma = \sqrt{\frac{\sigma^2}{\Lambda^2 \tau_c}}. \qquad (2.25)$$
## 2.2 Lagrangian dynamics

With the advances of experimental techniques and the increase in computing power of the last decades, remarkable features of the dynamics of fluid particles in turbulent flows have been discovered. Among them, the measurement of the probability distribution of the acceleration of these tracers has been shown to be clearly non Gaussian with a high frequency of observing very intense events [134, 264, 181, 180]. Even for moderate Reynolds numbers, it is relatively common to observe accelerations more than 100 times greater than their standard deviation. In addition, the components of acceleration and its norm present very different correlation times, the ratio of these characteristic times increasing with the Reynolds number [199, 182, 183] showing that the dynamics of the tracers is influenced by the full spectrum of turbulence scales. On the one hand, the short-time correlation of the acceleration component is connected to the centripetal forces in intense vorticity filaments [30, 183]. On the other hand, the acceleration norm has been shown to be directly correlated with the local dissipation rate of turbulence [217, 271, 106, 31], in accordance with Kolmogorov's hypotheses. Nevertheless, various experimental and numerical works have also reported its dependence on the local kinetic energy [180, 241, 29, 62, 36, 217, 8], which is generally attributed to nonindependence of the small-scale dynamics from those of the large-scales. In this view, the Lagrangian acceleration is essentially given by the local gradients of the velocity, but the latter present correlation with the kinetic energy caused by direct energy transfers between large- and small-scales in the energy cascade [258]. The absence of a proper scale separation explains that the Lagrangian correlation functions present power-laws with anomalous exponents which can be described by the multifractal formalism [55, 9, 28, 239, 139] as the signature of intermittency and persistence of viscous effects. To end this list, we mention the asymmetry of the fluctuations of the mechanical power received or given up by a fluid particle reflecting the temporal irreversibility of its dynamics [198, 75, 268, 208, 68, 66].

Such a complex phenomenology must be attributed to the collective and dissipative effects. Indeed according to the Navier-Stokes equation, the acceleration of a fluid particle is essentially given by the pressure field which is determined by the motion of all the other particles [65, 259]. Moreover, although the Laplacian term in the Navier-Stokes equation is of order $Re^{-1}$ smaller than the pressure gradient term, the viscosity cannot be neglected. Indeed, as a small force integrated over a long period could be significant, the viscosity insidiously affects the fluid tracer velocity. Which in turn influences the particle acceleration through modifications of the pressure gradient, and local interactions are intrinsically inseparable from the nonlocal ones. This is manifested by the persistence of the effect of the Reynolds number on the acceleration statistics, even for very large Reynolds numbers. Such a scenario is supported by Refs. [61, 196, 56] who showed that adding some noise to an inviscid Lagrangian flow causes the irreversibility of the dynamics.

Following the Kolmogorov first hypothesis [128, 126] stating that locally homogenous turbulent flows are universal, it should be possible, in principle, to propose a stochastic model that reproduces the dynamics of a single fluid particle by effectively accounting for the interactions with all the other fluid particles. Let us note that the Kolmogorov first hypothesis received some support from recent studies [249, 143, 37]. To propose such a stochastic model, our main assumption is to write the increments of the acceleration vector of a fluid particle as $da_i = M_i dt + D_{ij} dW_j$. Both M and D depend on particle acceleration $a$ and velocity $u$. The latter is simply given by the kinematic relation of
a fluid particle $u_i = \int a_i dt$. It is indeed a necessary condition that $a$ depends on $u$ to present a restoring effect that can counteract the diffusion in the velocity space and get a statistically stationary dynamics of the fluid particle. We will propose closed expressions for $M$ and $D$ from basic considerations using, as a starting point the acceleration statistics conditioned on both the local values of the dissipation rate and the kinetic energy observed from direct numerical simulations (DNS), which are presented below. It will be shown that introducing a “maximal winding hypothesis” associated with a nondiagonal diffusion tensor, this simple stochastic model reproduces all the statistical features of the Lagrangian dynamics presented above, without any adjustable parameter.

Let us first review some previous works on the stochastic modeling of the Lagrangian dynamics (see also Ref. [8]). Among the pioneering works, Sawford [237] proposed a scalar Gaussian model for the acceleration presenting a feedback term proportional to the velocity. Pope and Chen [203] devised a Langevin-like equation for the velocity coupled with a log-normal model for the dissipation through the introduction of conditional statistics. Similarly, Refs. [216, 215, 23] proposed an extension of the Sawford model leading to a non-Gaussian scalar model for the acceleration. This work was further refined by Ref. [136] who also advanced a non-Gaussian scalar model for the dynamics by prescribing an ad hoc shape of the conditional acceleration statistics with the dissipation along with a linear dependence on the velocity. The model introduced in Ref. [218] describes the increments of the derivative of the acceleration in a so-called third-order model to better account for the Reynolds number dependence on the acceleration statistics. Recently, Ref. [261] proposed generalization to an infinite number of layers leading to smooth 1D trajectory along with a multifractal correction to account for intermittency, as introduced in Refs. [10, 157, 118]. An acceleration vector model has been proposed in Ref. [214] by imposing an empirical correlation between velocity and acceleration, with additive noise leading to Gaussian statistics for the acceleration. Likewise, Ref. [202] presented a 3D Gaussian model, with a linear dependence on the velocity as well as an extension to non-homogenous flows. To account for intermittency effect, in Refs. [92, 230, 279, 228, 17, 91] the 3D acceleration vector is given by the product of two independent stochastic processes, one for the acceleration norm and the other for its orientation. In these models, the velocity feedback on the dynamics is realized by a coupling with a large eddy simulation framework. To summarize, to our knowledge, a 3D vectorial model for the tracer dynamics that is autonomous and reproducing the essential features of Lagrangian turbulence (irreversibility, non-Gausianty, multifractality) has not yet been proposed.

The essential building block of previously cited models is the conditional acceleration statistics. Previous studies have focused on conditional statistics with either the velocity or the dissipation rate separately. From the extensive analysis of Ref. [271], one can conclude that the acceleration variance conditioned on the dissipation rate $\varepsilon$ presents a power-law behavior at large values of $\varepsilon$ with a Reynolds-number-dependent exponent reflecting that the small-scale dynamics are not independent of large scales.

Regarding the links between the fluid particle acceleration and their velocity, Biferale et al. [29] argued that according to the Heisenberg-Yaglom scaling for the acceleration $\langle a^2 \rangle \sim a_\eta^2 = \langle \varepsilon \rangle^{3/2} \nu^{-1/2} = \langle K \rangle^{9/4} L^{-3/2} \nu^{-1/2}$, with $\nu$ the kinematic viscosity, $K = 1/2 u_i u_i$ the kinetic energy and $L$ the characteristic size of the large structures, one should expect that the variance of the velocity-conditioned acceleration behaves like $\langle a^2 | K \rangle \propto K^{9/4}$. Then on the basis of the multifractal formalism, they proposed a very close scaling law, $\langle a^2 | K \rangle \sim K^{2.3}$. The proposed relation was observed to be in agreement with DNS for
large velocity, typically $|u| > 3\sigma_u$ with $\sigma_u = \sqrt{2/3\langle K \rangle}$. These events remain very rare since the PDF of the fluid velocity is Gaussian so the range of validity of the power-law is, at best, very limited. Alternatively, Sawford et al. [241] propose that $\langle a_x^2|u_x \rangle \sim u_x^6$ based on a mechanism involving vorticity tubes. This scaling law which seems compatible with the first measurements of the acceleration conditioned on velocity in Ref. [180], is confirmed neither by the DNS of Ref. [29] nor in a second experimental paper by Crawford et al. [62] which gave more credit to the $K^{9/4}$ law. As mentioned above, it has been proposed that the dependence of the acceleration on the velocity arises through the dependence of the dissipation rate on the kinetic energy due to intermittency effect [29]. Additionally, Ref. [258] proposed that the dependence on velocity is a consequence of direct and bidirectional coupling of large and small scales caused by kinematic relations related to nonlocal interactions.

Here, we study the acceleration statistics conditional on both the kinetic energy and the dissipation rate. To our knowledge such doubly conditional statistics of the acceleration have never been presented. It will be shown that the variance can be expressed as $\langle a^2|\epsilon, K \rangle \sim \exp(\alpha K/\langle K \rangle + \gamma \ln \epsilon/\langle \epsilon \rangle)$. This result is clearly in contrast with the previously proposed power-law dependence on velocity. It shows that the influence of the large-scales through the intermittent distribution of the dissipation rate, which manifests through the Reynolds number dependence of the coefficient $\gamma$, is supplemented by an explicit dependence on the local kinetic energy. This direct dependence on the large-scale characteristics is of a kinematic nature as it appears independent of the Reynolds number. The behavior of the doubly conditional acceleration can be interpreted as a consequence of a scaling symmetry for the fluid-particle acceleration incorporating both the intermittency and the kinematic effects of the flow structure. We also propose to apply the incomplete similarity framework introduced by Barenblatt to explain the dependence of the statistics of the acceleration conditioned to the dissipation rate on the Reynolds number and to account for the intermittency effect. That enables us to provide as well new scaling relations for the unconditional variance in good agreement with the DNS. Eventually the doubly conditional statistics of the acceleration which gives a relation between the force, the energy, and the power will serve as a corner stone to build the stochastic model for the dynamics of a fluid particle mentioned above. Although such a model could be of interest for practical applications, its construction is relevant to study the specificities of the Lagrangian description of turbulence by linking the cascade picture to the fluid particles dynamics on the basis of the behavior of the conditional statistics obtained by the DNS of the Navier-Stokes equations.

In Sec. 2.3 we present the statistics of the acceleration conditioned on the local values of the dissipation rate and kinetic energy obtained from DNS. Then we show that the Reynolds number dependence on the acceleration conditioned on the dissipation rate can be described using the Barenblatt incomplete similarity. We deduce a new relation for the unconditional acceleration variances. To end this section, we show that these new results can be interpreted as a multiplicative cascade for the acceleration with scale dependent kinematic effects. Then in Sec. 2.4 we give the derivation of the stochastic model for a single fluid-particle dynamics taking as an initial step the doubly conditional acceleration variance, and present the outcome of the model for the Reynolds number up to $Re_\lambda = 9000$, along with comparison with DNS results when available.
## 2.3 Scaling laws of the acceleration

### 2.3.1 Methodology

We present in this section results concerning the statistics of the acceleration of a fluid particle. These results have been obtained from five direct numerical simulations (DNS) of isotropic turbulence in a periodic box with Taylor-scale Reynolds numbers of $Re_λ = 50, 90, 150, 230$ and $380$. We used pseudospectral code as detailed in Refs. [275, 279, 144]. The DNS was carried with resolutions of $128^3, 256^3, 512^3, 1024^3$ and $2048^3$ with the large scale forcing proposed by Ref. [131]. For each simulation we have $η/Δx = 1$ with $η = ⟨ε⟩^{-1/4}ν^{3/4}$ the Kolmogorov length scale and $Δx$ the grid size. The statistics are computed from 40 3D fields sampled at roughly each large-eddy-turnover time.

We will show statistics of the acceleration conditioned by the dissipation rate and the kinetic energy. Note that we consider here the pseudodissipation $\tilde{\epsilon} = ν(\partial_j u_i)^2$, which is the second invariant of the velocity gradient tensor multiply by the viscosity rather than the dissipation $\epsilon = \frac{1}{2}ν(\partial_j u_i + \partial_i u_j)^2$. We prefer to show here the statistics of the pseudodissipation to be consistent with the next section, in which we will use the log-normal distribution hypothesis for the dissipation. Indeed, this property is very well verified for the pseudodissipation whereas it is only approximate for the dissipation [271]. Nevertheless, the statistics presented below have also been computed considering the dissipation, $\epsilon$, and no significant differences were observed. To lighten the writing, in the sequel, we will drop the tilde in the notation of the pseudodissipation, as well, in the text, we will write dissipation instead of pseudodissipation.

### 2.3.2 Conditional statistics given the dissipation and the kinetic energy

To illustrate the relationships between acceleration, energy dissipation and kinetic energy, we show in Fig. 2.5 visualizations of these quantities at the same instant obtained from our DNS. We notice that $ln a^2/⟨a^2⟩$ and $ln ε/⟨ε⟩$ show a fairly marked correlation although the acceleration appears more diffuse than dissipation. We also notice that to some extent the kinetic energy and the dissipation rate appears correlated. In addition, it seems that some areas of the flow where the kinetic energy is high also correspond to regions of high acceleration magnitude.

Figure 2.6 presents the variance of the acceleration of a fluid particle conditioned to the local values of the kinetic energy and the dissipation rate: $⟨a^2|ε, K⟩$. In Fig. 2.6 (top), the levels of the logarithm of the conditional variance are shown as a function of $K$ and of $ε$. We see that the conditional variance of the acceleration depends on these two quantities and that the dependence on $K$ seems somewhat similar to that of $ln ε$. In a more quantitative way, we show in Fig. 2.6 (left) the variance of the acceleration as a function of $ε$ for different values of $K$. We can see that the shape of the curves remains globally unchanged when $K$ varies and also presents the same shape as the variance conditioned by $ε$ only as also presented in this figure. Essentially, it is observed that the conditional variance presents power-law behavior for $ε \gg ⟨ε⟩$ with an exponent close to 3/2 and a prefactor depending on $K$. As discussed in more detail below, we observe a slight deviation of the scaling law compared to the acceleration conditioned only by the dissipation.

Figure 2.6 (right) shows the variance of the acceleration as a function of $K$ for different values of $ε$. As expected, we find that the variance of the acceleration increases with $K$.
<figure><img src="image_11.png" /><figcaption>Figure 2.5: Visualization of the instantaneous fields of the square of the acceleration, of the dissipation and of the kinetic energy in a cut $y-z$ of the flow by DNS at $Re_λ = 380$. Left: $\ln(a^2/\langle a^2 \rangle)$; middle: $\ln(\varepsilon/\langle \varepsilon \rangle)$ and right: $K/\langle K \rangle$.</figcaption></figure>

We clearly notice an exponential growth of the variance over the whole range of $K$ with a growth rate $\alpha$ which appears independent of $\varepsilon$:

$$ \langle a^2 | \varepsilon, K \rangle = c_{\varepsilon} a_{\eta}^2 \exp(\alpha K / \langle K \rangle) \quad (2.26) $$

with $a_{\eta}^2 = \langle \varepsilon \rangle^{3/2} \nu^{-1/2} = \langle \varepsilon \rangle / \tau_{\eta}$ the so-called Kolmogorov acceleration and the prefactor $c_{\varepsilon}$ depending on $\varepsilon$. From our DNS it appears that $\alpha \approx 1/3$ for all the Reynolds numbers considered here. We also find the same value of $\alpha$ from the database of [21, 140] obtained for $Re_λ = 400$ suggesting that the value of $\alpha$ is independent of the Reynolds number. This exponential behavior contrasts with the references mentioned in the introduction in which power-laws behavior for the variance conditioned on $K$ solely had been proposed. Nevertheless, we can notice that exponential growth does not seem to disagree with the data presented in these references. Interestingly, this relationship only depends on a characteristic velocity, (not a time and a length scale separately). The absence of characteristic time is attributed to the scale separation between large structures and small ones (the large structures of the flows appear as quasi stationary and infinite to the smallest ones such that only their relative velocity matters). The independence of the coefficient $\alpha$ on the Reynolds number tends to confirm that the velocity scale used for the nondimensionalization of the argument of the exponential is appropriate.

In appendix of [274] we propose to estimate the factor $c_{\varepsilon}$ as:

$$ c_{\varepsilon} \approx A \langle a^2 | \varepsilon \rangle / a_{\eta}^2 \quad (2.27) $$

where $A = \left(1 - \frac{2}{3}\alpha\right)^{3/2}$, which is equal to $A = 7\sqrt{7}/27 \approx 0.686$, for $\alpha = 1/3$, neglecting a small logarithmic dependence on $\varepsilon/\langle \varepsilon \rangle$.

Consequently, for large Reynolds numbers, the doubly-conditioned variance of the fluid-particle acceleration is expressed as

$$ \langle a^2 | \varepsilon, K \rangle = A \langle a^2 | \varepsilon \rangle \exp(\alpha K / \langle K \rangle). \quad (2.28) $$
<figure><img src="image_12.png" /><figcaption>Figure 2.6: Variance of acceleration conditioned on the local dissipation rate and kinetic energy obtained from DNS at $Re_λ = 380$. Top: Map of $\ln\langle a^2|\epsilon, K\rangle/a_\eta^2$ versus $\ln\epsilon/\langle\epsilon\rangle$ and $K/\langle K\rangle$. Left: Plot, in logarithmic scales, of $\langle a^2|\epsilon, K\rangle/a_\eta^2$ against $\epsilon/\langle\epsilon\rangle$ for $K/\langle K\rangle = 0.025, 0.1, 0.5, 1, 2, 3, 5, 6.5 \pm 30\%$ from orange to black. Comparison with $\langle a^2|\epsilon\rangle/a_\eta^2$ in gray dashed line and with the power-law $(\epsilon/\langle\epsilon\rangle)^{3/2}$ in gray dotted line. Right: Plot, in semilogarithmic scales, of $\langle a^2|\epsilon, K\rangle/a_\eta^2$ against $K/\langle K\rangle$ for $\epsilon/\langle\epsilon\rangle = 0.01, 0.05, 0.1, 0.3, 0.5, 1, 3, 5, 10, 50 \pm 30\%$ from orange to black. Comparison with $\langle a^2|K\rangle/a_\eta^2$ in gray dashed line and with $\exp(\alpha K/\langle K\rangle)$ with $\alpha = 1/3$ in gray dotted line.</figcaption></figure>

This relation is confirmed in Fig. 2.7 which presents the conditional variance of the acceleration normalized by $Aa_\eta^2 \exp(\alpha K/\langle K\rangle)$ as a function of $\epsilon$ for different values of $K$ as well as normalized by $A\langle a^2|\epsilon\rangle = \langle a^2|\epsilon, K=0\rangle$ as a function of $K$ for different values of $\epsilon$. It can be seen that a fairly good overlap of the various curves is obtained, confirming the self-preserving character of the acceleration conditioned on both the kinetic energy and the dissipation rate. We see in this relation an explicit effect of the local kinetic energy on the acceleration. Since the argument of the exponential depends on $K/\langle K\rangle$ not on a local Reynolds number, it suggests pure kinematic effects for the acceleration which is likely associated to the divergence free constrain and the nonlocality of the pressure gradient. There is also an indirect effect through the dependence of the dissipation rate on the large-scale structures. The later is manifested as Reynolds number dependence of the conditional acceleration on the dissipation rate solely. This intermittency effect is analyzed further in the next section. We postpone to Sec. 2.3.5 further comments on the behavior of the doubly-conditioned variance.
<figure><img src="image_13.png" /><figcaption>Figure 2.7: Normalized variance of acceleration conditioned on the local dissipation rate and kinetic energy obtained from DNS at $Re_λ = 380$. (Left) Plot of $\langle a^2|\epsilon, K \rangle / Aa_\eta^2 \exp(\alpha K/\langle K \rangle)$ against $\epsilon/\langle\epsilon\rangle$ for various values of K. Comparison with $\langle a^2|\epsilon\rangle / a_\eta^2$ in gray dashed line and with the power-law $(\epsilon/\langle\epsilon\rangle)^{3/2}$ in gray dotted line. (Right) Plot of $\langle a^2|\epsilon, K \rangle / \langle a^2|\epsilon, K = 0 \rangle$ against $K/\langle K \rangle$ for various values of $\epsilon/\langle\epsilon\rangle$. Comparison with $\langle a^2|K \rangle / a_\eta^2$ in gray dashed line and with $\exp(\alpha K/\langle K \rangle)$ with $\alpha = 1/3$ in gray dotted line. The ranges of the fixed values of K and $\epsilon$ for both plots are the same as in Fig. 2.6.</figcaption></figure>

### 2.3.3 Similarity of the conditional statistics given the dissipation

We propose now to focus with more details on the scaling law of the acceleration variance conditioned on the dissipation rate only, $\langle a^2|\epsilon \rangle$. For that we consider the DNS data from Yeung et al. [271], along with our DNS data. Figure 2.8 (Left) presents the conditional acceleration variance for Reynolds numbers in the range $Re_λ = 40$ to 680. We first notice that for weak values of the dissipation rate ($\epsilon \ll \langle\epsilon\rangle$) the value of the conditional acceleration variance tends towards an asymptotic value, which depends on the Reynolds number. The saturation of the conditional acceleration shows that the local acceleration is not only determined by the microstructure of the flow, and that it presents somehow effects of the large structures of the flow which dominates in low dissipative regions. We denote by $a_0^2$ the asymptotic value of the conditional variance:

$$a_0^2 = \lim_{\epsilon \to 0} \langle a^2 | \epsilon \rangle. \qquad (2.29)$$

Assuming that the acceleration of fluid particles in low dissipative regions is mainly influenced by large scales, we can estimate $a_0^2$ as $a_0^2 \sim \langle K \rangle / \tau_L^2$ with $\tau_L$ the integral timescale of the flow. This leads to the following estimate:

$$a_0^2/a_\eta^2 \sim \tau_\eta/\tau_L \sim Re_\lambda^{-1}. \qquad (2.30)$$

We test this scaling law for $a_0$ in Fig. 2.8(right) by presenting $a_\eta^2/a_0^2$ as a function of $Re_λ$ from the different DNS datasets. We observe a linear growth rate of $a_0^2/a_\eta^2$ with $1/Re_λ$. For large values of $\epsilon$, we notice in Fig. 2.8(Left), as already reported in [271], that the conditional variance presents a power-law behavior with $\epsilon$. The exponent of this scaling law is seen to evolve continuously with the Reynolds number, and seems to tend asymptotically
towards $\epsilon^{3/2}$. From dimensional analysis we define $f$ as:

$$ \frac{\langle a^2 | \epsilon \rangle}{\epsilon^{3/2} \nu^{-1/2}} = f(\epsilon / \langle \epsilon \rangle, \text{Re}_\lambda). \qquad (2.31) $$

In the inset of Fig. 2.8(Left), it is seen that $f$ seems to admit an asymptotic constant value for $\epsilon \gg \langle \epsilon \rangle$ only in the limit of a very large Reynolds number. For intermediate Reynolds numbers, $f$ presents power-law behavior with $\epsilon$ for $\epsilon \gg \langle \epsilon \rangle$ but with a Reynolds number dependent exponent. This implies an absence of similarity of the flow when the Reynolds number is changed and the persistence of the Reynolds number effect, even for large Reynolds numbers, which highlights an absence of proper scale separation suggesting direct coupling between large- and small-scales. This is reminiscent of the incomplete similarity framework proposed by Barenblatt [15, 14, 13]. Following Barenblatt, we assume that $f$ presents an incomplete similarity in $\epsilon/\langle \epsilon \rangle$ and absence of similarity in $\text{Re}_\lambda$. Accordingly we write

$$ f(\epsilon/\langle\epsilon\rangle, \text{Re}_\lambda) = B(\epsilon/\langle\epsilon\rangle)^\beta \qquad (2.32) $$

where the anomalous exponent $\beta$, and the prefactor $B$ are both functions of $\text{Re}_\lambda$. Arguing for a vanishing viscosity principle, it can be assumed that the critical exponent becomes independent of the Reynolds number in the limit of asymptotically large Reynolds number. Finally arguing that the dependence of $B$ and $\beta$ on the Reynolds number is small, Barenblatt further proposed that they present inverse logarithmic dependence on $\text{Re}_\lambda$, which is also in agreement with the log-similarity proposed by Refs. [51, 82]. Expanding $\beta$ and $B$ in power of $1/\ln(\text{Re}_\lambda)$ yields, keeping only the leading-order term in $\text{Re}_\lambda$:

$$ \beta = \beta_0 + \beta_1 / \ln \text{Re}_\lambda \qquad (2.33) $$

$$ B = B_0 + B_1 / \ln \text{Re}_\lambda . \qquad (2.34) $$

To have a finite limit, consistently with the vanishing viscosity principle, we need $\beta_0 = 0$. The remaining constants $B_0$, $B_1$ and $\beta_1$ are then determined by comparison with the DNS data. From the inset of Fig. 2.8(Left) we see that both $\beta$ and $B$ are increasing functions of $\text{Re}_\lambda$ implying that both $B_1$ and $\beta_1$ are negative. In Fig. 2.9(Left) we assess the relations (2.32)-(2.34) by plotting

$$ \chi = \frac{1}{\gamma} \ln(1/B \langle a^2 | \epsilon \rangle / a_{\eta}^2), \qquad (2.35) $$

with

$$ \gamma = 3/2 + \beta, \qquad (2.36) $$

against $\ln(\epsilon/\langle\epsilon\rangle)$ for various Reynolds numbers. It is observed that with $B_0 = 17.1$, $B_1 = -54.7$ and $\beta_1 = -1$, all the DNS data collapse on the line $\chi = \ln(\epsilon/\langle\epsilon\rangle)$ (the bisectrix of the graph) for $\epsilon \gg \langle\epsilon\rangle$, validating the scaling relation.

We can go a step further by using the low dissipative limit of the conditional acceleration.

For that we introduce $\chi_0 = \lim_{\epsilon \to 0} \chi = \frac{1}{\gamma} \ln(1/B a_0^2/a_\eta^2)$. With this definition, $\chi - \chi_0 = \ln[(\langle a^2 | \epsilon \rangle / a_0^2)^{1/\gamma}]$ tends to 0 in the low dissipative regions ($\epsilon \ll \langle \epsilon \rangle$). On the other hand, for $\epsilon \gg \langle \epsilon \rangle$, $\chi - \chi_0$ should be equal to $\chi - \chi_0 = \ln(\epsilon/\langle \epsilon \rangle) - \chi_0 = \ln[\epsilon/\langle \epsilon \rangle (Ba_\eta^2/a_0^2)^{1/\gamma}]$. This is seen in Fig. 2.9 (right) that presents the evolution of $(\langle a^2 | \epsilon \rangle / a_0^2)^{1/\gamma}$ against

$$ \zeta = \epsilon/\langle\epsilon\rangle (Ba_{\eta}^{2}/a_{0}^{2})^{1/\gamma} \qquad (2.37) $$
for the various Reynolds numbers considered here.

It is interesting to note that the curves are all overlapping even for intermediate values of $\epsilon$, suggesting that the conditional acceleration variance can be cast in a self-similar form:

$$ \langle a^2 | \epsilon \rangle = a_0^2 (\phi(\zeta))^\gamma \qquad (2.38) $$

with $\phi$ a universal function of only one argument $\phi = \phi(\zeta)$ with the asymptotics $\phi(\zeta) = 1$ for $\zeta \ll 1$ and $\phi(\zeta) = \zeta$ for $\zeta \gg 1$. Making a Taylor expansion of $\phi$ around $\zeta = 0$ and using a matching asymptotic argument, simply yields

$$ \phi(\zeta) = 1 + \zeta. \qquad (2.39) $$

It is seen in the inset of the Fig. 2.9 (right) that the proposed expression for $\phi$ gives a good approximation of the data over the whole range of $\epsilon$ and $Re_\lambda$. We can indeed observe more than 5 decades of quasilinear growth of $(\langle a^2 | \epsilon \rangle / a_0^2)^{1/\gamma} - 1 = \exp(\chi - \chi_0) - 1$ with $\zeta$. The nondimensional function $f$ introduced in Eq. (2.31) can, in consequence, be expressed as:

$$ f(\epsilon/\langle\epsilon\rangle, Re_{\lambda}) = B (\epsilon/\langle\epsilon\rangle)^{\beta} \left(1 + \frac{1}{\zeta}\right)^{3/2+\beta}, \qquad (2.40) $$

where the term within the brackets is interpreted as a correction factor for small dissipative regions. Accordingly, we obtain the following expression for the conditional acceleration variance:

$$ \langle a^2 | \epsilon \rangle = B a_{\eta}^2 \left( \left( \frac{1}{B} \frac{a_0^2}{a_{\eta}^2} \right)^{1/\gamma} + \frac{\epsilon}{\langle \epsilon \rangle} \right)^{\gamma}. \qquad (2.41) $$

As $\gamma = 3/2 + \beta$ and $B$ evolve slowly with $Re_\lambda$, their expressions remain speculative and would require a much larger range of Reynolds numbers to be validated.
<figure><img src="image_14.png" /></figure>

Figure 2.8: (Left) Acceleration variance conditioned on the local dissipation rate normalized by the Kolmogorov acceleration $\langle a^2 | \epsilon \rangle / a_{\eta}^2$. The continuous line are for our DNS for $Re_{\lambda} = 50, 90, 150, 230$ and 380 from orange to black; the dashed lines correspond to the DNS of Yeung et al. [271] for $Re_{\lambda} = 40, 139, 238, 385, 680$, from orange to black. Comparison with the power-law $(\epsilon/\langle\epsilon\rangle)^{3/2}$. Inset conditional acceleration normalized by $\epsilon^{3/2}\nu^{-1/2}$. (Right) $a_{\eta}^2/a_0^2$ as a function of $Re_{\lambda}$ with $a_0^2 = \lim_{\epsilon \to 0} \langle a_i^2 | \epsilon \rangle$ the acceleration variance in low dissipative regions. Gray dots for the DNS of Yeung et al., black crosses for our DNS. Comparison with the line $0.0028Re_{\lambda} + 1.16$ in dashed lines.

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

Figure 2.9: (Left) Evolution of $\chi = \frac{1}{\gamma} \ln(1/B \langle a^2 | \epsilon \rangle / a_{\eta}^2)$, with $\gamma = 3/2 + \beta$ versus $\ln \epsilon / \langle \epsilon \rangle$ with $B = B_0 + B_1 / \ln(Re_{\lambda})$ and $\beta = \beta_1 / \ln(Re_{\lambda})$ and $B_0 = 17.1$, $B_1 = -54.7$ and $\beta_1 = -1$, for various $Re_{\lambda}$: continuous line for our DNS at $Re_{\lambda} = 50, 90, 150, 230$ and 380 from orange to black; and dashed lines correspond to the DNS of Yeung et al. [271] for $Re_{\lambda} = 40, 139, 238, 385, 680$, from orange to black. Comparison with the line $\chi = \ln \epsilon / \langle \epsilon \rangle$ in gray dashed line. (Right) Evolution of $(\langle a^2 | \epsilon \rangle / a_0^2)^{1/\gamma} = \exp(\chi - \chi_0)$ against $\zeta = \epsilon / \langle \epsilon \rangle (Ba_{\eta}^2/a_0^2)^{1/\gamma}$ for the various Reynolds numbers. Inset: plot of $(\langle a^2 | \epsilon \rangle / a_0^2)^{1/\gamma} - 1$ against $\zeta = \epsilon / \langle \epsilon \rangle (Ba_{\eta}^2/a_0^2)^{1/\gamma}$.
### 2.3.4 Reynolds number dependence of the unconditional acceleration variance

Assuming the distribution of the dissipation, we can integrate relation (2.41) to obtain the unconditional variance

$$ \langle a^2 \rangle = \int \langle a^2 | \varepsilon \rangle P(\varepsilon) d\varepsilon, \qquad (2.42) $$

and thus propose an alternative formula to the empirical relations proposed in Refs. [271, 241, 103]. We consider that $\varepsilon/\langle\varepsilon\rangle$ presents a log-normal distribution with parameter $\sigma^2 \approx 3/8 \ln Re_\lambda/R_c$ with $R_c = 10$ as shown by Ref. [271] from DNS data, consistently with the proposition of Kolmogorov and Oboukhov [128, 192]. Notice nevertheless that other expressions for $\sigma^2$ have been proposed in the literature reflecting the vanishing viscosity limits [49]. Taking for $\langle a^2 | \varepsilon \rangle$ the expression (2.41), we perform the integration numerically with the expression (2.33) and (2.34) for $\beta$ and $B$ with the values of $B_0$, $B_1$ and $\beta_1$ and the expression of $a_0^2/a_\eta^2$ proposed above. The resulting evolution of the acceleration variance with the Reynolds number is presented in Fig. 2.10. It is seen that the predicted acceleration variance is in very good agreement with the DNS of Ref. [271], our DNS, as well as the DNS of [42]$^4$ and is also very close to the relation proposed by Ref. [241] for $Re_\lambda < 1000$.

<figure><img src="image_16.png" /><figcaption>Figure 2.10: Evolution of the acceleration variance normalized by the Kolmogorov acceleration with the Reynolds number. From our DNS data (black crosses), the DNS data from [271] (gray dots), the DNS data from [42] (blue triangles) and from the stochastic model (red plus) and comparison with the relation $3 \times 1.9 Re_\lambda^{0.135}(1 + 85 Re_\lambda^{-1.35})$ from [241] (gray dash dot line), with the numerical integration of $\langle a^2 \rangle = \int \langle a^2 | \varepsilon \rangle P(\varepsilon) d\varepsilon$ with $\langle a^2 | \varepsilon \rangle$ given by (2.41) and $P(\varepsilon)$ log-normal (continuous red line), with its approximation of (2.45) (red dotted lines), with the large-Reynolds number limit relation (2.44) (red dashed line) and with the asymptotic power-law (2.46) (gray dotted lines).</figcaption></figure>

The first term within the brackets in Eq. (2.41) is the footprint of the large-scale structures whose effects are vanishing if the local dissipation rate is larger than $\varepsilon/\langle\varepsilon\rangle \gg (B a_\eta^2/a_0^2)^{-1/\gamma}$ and therefore can be neglected when the Reynolds number is large since $a_\eta^2/a_0^2 \sim Re_\lambda$.

<footnote>$^4$It worth noting that the Ref. [42] was published after the presentation of results in Ref. [274]. We can thus see a stronger validation in the agreement between the model and the DNS of Ref. [42] over a larger range of Reynolds number than the one initially available.</footnote>
Hence, for large Reynolds numbers, Eq. (2.41) reduces to:

$$ \langle a^2 | \epsilon \rangle / a_{\eta}^2 = B \left( \frac{\epsilon}{\langle \epsilon \rangle} \right)^{\gamma}. \qquad (2.43) $$

With this expression, the acceleration variance is simply estimated from the moments of the log-normal distribution as:

$$ \frac{\langle a^2 \rangle}{a_{\eta}^2} = B (Re_{\lambda}/R_c)^{9/64+3\beta(1+\beta/2)/8}. \qquad (2.44) $$

This expression, also presented in Fig. 2.10, is shown to converge to the previous estimate as the Reynolds number increases.

In appendix of [274], we show that the integral (2.42) can be expressed from the generalized binomial series expansion. We further obtain the following estimation for the acceleration variance keeping only the first two terms:

$$ \frac{\langle a^2 \rangle}{a_{\eta}^2} = B \left( \frac{Re_{\lambda}}{R_c} \right)^{3/16\gamma(\gamma-1)} \left( 1 + \gamma \left( \frac{1}{B} \frac{a_0^2}{a_{\eta}^2} \right)^{1/\gamma} \left( \frac{Re_{\lambda}}{R_c} \right)^{-3/8(\gamma-1)} \right). \qquad (2.45) $$

In Fig. 2.10, it is seen that this relation almost overlaps with the direct numerical calculation of the variance through Eq. (2.42). The term within the brackets enables us to measure the contribution from small Reynolds number effects. At $Re_{\lambda} \approx 100$, the two estimates (2.44) and (2.45) for the variance differ by about 20%, while there is about 8% in difference at $Re_{\lambda} \approx 500$. That confirms that the term containing $a_0^2$ is indeed vanishing at large Reynolds numbers.

The previous estimations of the acceleration variance tend asymptotically to the following power-law:

$$ \langle a^2 \rangle / a_{\eta}^2 = 7.62 Re_{\lambda}^{9/64}, \qquad (2.46) $$

where we have used Eq. (2.33) to obtain the value of the prefactor, $R_c^{-9/64} B_0 \exp(3\beta_1/8) \approx 7.62$. This expression is presented as well in Fig. 2.10, confirming that the convergence toward the power-law is very slow, and that Eq. (2.44) should be considered as an intermediate asymptotic expression for the acceleration variance.

### 2.3.5 Multiplicative cascade for the acceleration

Substituting Eq. (2.43) into Eq. (2.28), we can eventually estimate the doubly conditional acceleration variance for large Reynolds numbers as:

$$ \langle a^2 | \epsilon, K \rangle / a_{\eta}^2 = C \exp(\alpha K / \langle K \rangle) \left( \frac{\epsilon}{\langle \epsilon \rangle} \right)^{\gamma}, \qquad (2.47) $$

where $C = AB$, with expressions of $A$ and $B$ determined above.

As mentioned above and apparent in the previous formula, the acceleration depends on the local value of the kinetic energy, along with the local dissipation rate. The acceleration, being mainly due to the pressure gradient, it presents a nonlocal behavior. The fact that the acceleration depends on the local kinetic energy but not on a local Reynolds number reflects that its nonlocality is a purely kinematic effect. Further, the exponential dependence on the kinetic energy suggests that the acceleration can respond to the structures of all sizes.
To illustrate this point, we discuss a multiplicative cascade model for the acceleration that incorporates the effect of the full spectrum of the flow structure. Fluctuations of the locally-space-averaged dissipation rate can be modeled by multiplicative cascades [270, 127, 178, 157, 25, 80]. Such model proposes to express the local dissipation over a volume of size $\ell = L\lambda^n$, with $\lambda < 1$ and $L$ being the large-scale of the flow, as the product of $n$ random numbers $\xi_i$:

$$ \varepsilon_{\ell} = \langle \varepsilon \rangle \prod_{i=1}^{n} \xi_{i}. \qquad (2.48) $$

Typically, for large $n$, this yields log-normal distribution of $\varepsilon_{\ell}$ assuming the $\xi_i$ are independent and identically distributed (and have as well finite variance).

We propose likewise to write the squared acceleration, coarse-grained at scale $\ell$, as:

$$ a_{\ell}^{2} = a_{0}^{2} \prod_{i=1}^{n} \theta_{i}. \qquad (2.49) $$

The scale-to-scale factor $\theta_i$ is given by:

$$ \theta_i = \exp \left( \frac{\alpha}{\langle K \rangle} \frac{1}{2} v_i^2 \right) (\xi_i)^{\gamma} = \exp \left( \frac{\alpha}{\langle K \rangle} \frac{1}{2} v_i^2 + \gamma \ln \xi_i \right), \qquad (2.50) $$

where $v_i$ is here the velocity of eddies of size $\ell_i = L\lambda^i$, which is also a fluctuating quantity. The exponential modulation is then interpreted as an entrainment acceleration due to these structures.

With this expression we obtain:

$$ a_{\ell}^{2} = a_{0}^{2} \exp \left( \frac{\alpha}{\langle K \rangle} \sum_{i=1}^{n} \frac{1}{2} v_{i}^{2} + \gamma \sum_{i=1}^{n} \ln \xi_{i} \right). \qquad (2.51) $$

Setting $n = \ln(\eta/L)/\ln(\lambda) \sim \ln Re_{\lambda}$, $\eta$ being the Kolmogorov length scale, we have $K = \sum_{i=1}^{n} \frac{1}{2} v_i^2$ due to the additive nature of the kinetic energy. Thereby using Eq. (2.48), we obtain back Eq. (2.47) by taking the conditional average of Eq. (2.51). The order of magnitude of the eddy velocities can be estimated from the Kolmogorov relation, $(\varepsilon_{\ell}\ell)^{1/3}$, showing that the sum is a priori dominated by the large-scales but, on the other hand, because of the intermittent behavior of $\varepsilon_{\ell}$, it may well happen that the inertial-scale structures can be dynamically important.

Note that in this multiplicative model, we have transposed the statistical relation (2.47) to an instantaneous version. Such idealization, find support in the invariance of the conditional PDF, which is shown in appendix of [274]. Another important point to mention is that although we assume that the local acceleration depends both on $K$ and $\varepsilon$ it is not assumed that those two variables are independent.

The dissipation presents large fluctuations leading to very important accelerations and, even if the acceleration orientation is changing rapidly, it can cause a local increase of the velocity. When the kinetic energy becomes significantly larger than its averaged value, then the modulation of the acceleration by the exponential term becomes preponderant, thus offering a feedback mechanism allowing obtaining the normal fluctuations of the velocity.

This dynamic scenario appears consistent with the recent DNS analysis of Ref. [197] showing that the fluid-particles can undergo energy gains in intense dissipative regions and is developed in the next section.
## 2.4 Stochastic modeling of the fluid-particle dynamics

### 2.4.1 Model formulation

The foregoing multiplicative model suggests that the acceleration norm can be determined from the local kinetic energy and dissipation rate. The relation (2.47) is pointwise and so it applies equally well to both Lagrangian and Eulerian descriptions. However, in the Lagrangian framework, the kinematic relation between velocity and acceleration allows proposing a model for the acceleration depending only on the local dissipation rate. The evolution of the later along the particle path is to be obtained from a stochastic process. For the derivation of this model, we will rely on the relation (2.47), in which the contribution from low dissipative events is neglected:

$$a^2 = f(K, \varepsilon) = a_{\eta}^2 C \left( \frac{\varepsilon}{\langle \varepsilon \rangle} \right)^{\gamma} \exp \left( \alpha \frac{K}{\langle K \rangle} \right). \quad (2.52)$$

We express the increments of $a^2$ as a second order Taylor expansion in $K$ and $\varepsilon$,

$$da^2 = a^2 \left( \alpha \frac{dK}{\langle K \rangle} + \gamma \frac{d\varepsilon}{\varepsilon} + \frac{\alpha^2}{2} \frac{dK^2}{\langle K \rangle^2} + \frac{\gamma(\gamma-1)}{2} \frac{d\varepsilon^2}{\varepsilon^2} + \gamma \alpha \frac{dK}{\langle K \rangle} \frac{d\varepsilon}{\varepsilon} \right). \quad (2.53)$$

We consider $\varepsilon$ as stochastic variable reflecting the very large number of degrees of freedom that control them. In a fairly general way, we consider that the dissipation $\varepsilon$ follows a multiplicative stochastic process:

$$d\varepsilon = \varepsilon\Pi dt + \varepsilon\Sigma dW, \quad (2.54)$$

where $dW$ are the increments of the Wiener process ($\langle dW \rangle = 0$; $\langle dW^2 \rangle = dt$). We specify the terms $\Pi$ and $\Sigma$ below. Substituting Eq. (2.54) into Eq. (2.53) one obtains, following the Ito calculus, at first order in $dt$:

$$da^2 = a^2 \left[ \frac{\alpha}{\langle K \rangle} P + \gamma \Pi + \frac{\gamma(\gamma-1)}{2} \Sigma^2 \right] dt + \gamma a^2 \Sigma dW. \quad (2.55)$$

We used the identity $dK = u_i du_i = u_i a_i dt = P dt$, where $P$ is the mechanical power per unit of mass exchanged by the fluid particle. Even if $\Pi$ and $\Sigma$ are given, eq. (2.55) is not closed, as it remains to estimate $P = a_i u_i$, which requires the knowledge of $a_i$ and $u_i$.

As mentioned in the introduction, we introduce a vectorial stochastic model for the dynamics of a fluid particle. We are looking for a stochastic process of the form:

$$du_i = a_i dt, \quad (2.56)$$

$$da_i = M_i dt + D_{ij} dW_j, \quad (2.57)$$

where $dW_j$ are the increments of the jth component of a tridimensional Wiener process ($\langle dW_j \rangle = 0$; $\langle dW_i dW_j \rangle = dt \delta_{ij}$). a priori, the vector $M$ and the tensor $D$ depend on the vectors $a$ and $u = \int adt$. Indeed, $M$ must depend on $u$ to allow the particle velocity to reach a statistically steady state.

Now, we propose expressions for $M_i$ and $D_{ij}$. For this, we want to impose, on the one hand, that the model is isotropic ($\langle a_i a_j \rangle = 0$ for $i \neq j$) and, on the other hand, that
its norm $a^2 = a_i a_i$ is compatible with the Eq. (2.55). We therefore write the stochastic equation for $a_{ij}^2 = a_i a_j$ derived from Eq. (2.57), thanks to the Ito formula:

$$ 
\begin{aligned}
da_{ij}^2 &= a_j da_i + a_i da_j + da_i da_j \\
&= (M_i a_j + M_j a_i + D_{ik} D_{jk}) dt + (a_j D_{ik} + a_i D_{jk}) dW_k.
\end{aligned} 
\quad (2.58) $$

For the square of the norm $a^2 = a_i a_i$, we have:

$$ da^2 = (2a_i M_i + D_{ij} D_{ij}) dt + 2a_i D_{ij} dW_j. \quad (2.59) $$

We then proceed by identification between Eqs. (2.59) and (2.55), in a similar way as Refs. [89, 195], by identifying first the square of the diffusion term and then the drift term.

## Identification of the diffusion term and the maximum winding hypothesis

Quite generally, we can decompose the diffusion tensor into:

$$ D_{ij} = c_1 \delta_{ij} + S_{ij} + \Omega_{ij}, \quad (2.60) $$

where $S_{ij}$ is a zero-trace symmetric tensor and $\Omega_{ij}$ is an antisymmetric tensor. The latter can be written as $\Omega_{ij} = \epsilon_{ijk}\omega_k$ with $\epsilon_{ijk}$ the Levi-Civita permutation symbol and $\omega_k$ a pseudovector. $S_{ij}$ must be zero in order to guarantee the statistical isotropy of the acceleration. But $\Omega_{ij}$ can be different from 0. Indeed, the experimental results of Ref. [181] and numerical results of Ref. [199] have shown that the acceleration presents a scale separation between the evolution of the components and its norm, and that this separation can be modeled using processes for the acceleration norm and its orientation vector [276, 92, 228]. A stochastic model for orientation can be formulated as a diffusion process with a rotational part in the diffusion tensor [92, 266]. Since the model for the dynamics (2.56) - (2.57) involves only two vectors, $a$ and $u$, we propose to form the pseudovector $\omega$ from these two vectors in order to get a closed model: $\omega_k = c_2 \epsilon_{klm} a_l u_m$. The model remains statistically isotropic and the chirality of the flow is not broken either since the odd moments of $dW_j$ are zero (Gaussian with zero mean). In other words, the sign of $c_2$ does not matters. We then have:

$$ D_{ij} = c_1 \delta_{ij} + c_2 (a_i u_j - a_j u_i). \quad (2.61) $$

It is to be noted that $c_1$ and $c_2$ are not constant.

By identifying the square of the diffusion term between (2.55) and (2.59) we find

$$ \gamma^2(a^2)^2 \Sigma^2 = 4a_i a_j D_{ik} D_{jk}. \quad (2.62) $$

Expanding it by using expression (2.61), we find

$$ \gamma^2(a^2)^2 \Sigma^2 = 4a^2(c_1^2 + c_2^2(2a^2K - P^2)), \quad (2.63) $$

which gives for $c_1$:

$$ c_1^2 = \frac{\gamma^2}{4} a^2 \Sigma^2 - c_2^2 (2a^2 K - P^2) = a^2 \left( \frac{\gamma^2}{4} \Sigma^2 - 2c_2^2 K \left( 1 - \frac{a_T^2}{a^2} \right) \right), \quad (2.64) $$
where we have introduced the tangential acceleration $a_T$, as the projection of the acceleration vector in the direction of the velocity vector: $a_T = a_i u_i / \sqrt{u^2} = P / \sqrt{2K}$. Equation (2.64) imposes a constraint on $c_2$ to guarantee the positivity of $c_1^2$:

$$c_2^2 2K \leq \frac{\gamma^2}{4} \Sigma^2, \qquad (2.65)$$

since $0 \leq 1 - \frac{a_T^2}{a^2} \leq 1$. So, to guarantee the positivity of $c_1^2$ whatever $K$, $c_2^2$ must be proportional to $1/K$. Introducing a parameter $c_R$ as $c_2^2 = \frac{\gamma^2}{4} \Sigma^2 \frac{c_R^2}{2K}$, with the constraint $c_R^2 \leq 1$, we obtain

$$c_1^2 = \frac{\gamma^2}{4} \Sigma^2 (a^2(1 - c_R^2) + c_R^2 a_T^2). \qquad (2.66)$$

Subsequently, we only consider the limit $c_R = 1$ that corresponds to the maximum rotational part of the diffusion tensor. We will discuss this choice in more detail below in Sec. 2.4.3, when presenting the results.

Finally, from Eq. (2.61) and the expressions of $c_1$ and $c_2$, we write the components of the diffusion tensor as

$$D_{ij} = \sqrt{\frac{\gamma^2}{4} \Sigma^2} \left[ \sqrt{a_T^2} \delta_{ij} + \sqrt{a_N^2} \epsilon_{ijk} b_k \right], \qquad (2.67)$$

where we introduced the normal component $a_N$ of the acceleration $a_N^2 = a^2 - a_T^2$, and the bi-normal unit vector<fnref n="5" /> $b_k = \epsilon_{klm} u_l a_m / |\epsilon_{hij} u_i a_j|$.

Note that $b_k$, $a_T$ and $a_N$ are not well defined when $K = 0$. However, $c_2$ must vanish when $u = 0$ and we can therefore consider that $c_R = 0$ in that case.

## Determination of the drift term

Identifying the drift term between Eqs. (2.59) and (2.55), we get:

$$2a_i M_i + D_{ij} D_{ij} = a^2 \left( \frac{\alpha}{\langle K \rangle} P + \gamma \Pi + \frac{\gamma(\gamma-1)}{2} \Sigma^2 \right). \qquad (2.68)$$

From Eq. (2.67) the term $D_{ij}D_{ij}$ is computed as <fnref n="6" />

$$D_{ij} D_{ij} = \frac{\gamma^2}{4} \Sigma^2 (2a^2 + a_T^2). \qquad (2.69)$$

We then have

$$a_i M_i = a^2 \left( \frac{\alpha}{2\langle K \rangle} P + \frac{\gamma}{2} \Pi - \frac{\gamma}{4} \Sigma^2 \right) - a_T^2 \frac{\gamma^2}{8} \Sigma^2. \qquad (2.70)$$

<footnote n="5">⁵To obtain this relation we notice that $a_i u_j - a_j u_i = \epsilon_{ijk} \epsilon_{klm} a_l u_m$ and that the vector $b_k$ is the unit vector collinear to $\epsilon_{klm} a_l u_m$: $b_k = \epsilon_{klm} a_l u_m / |\epsilon_{hij} a_i u_j|$. By expanding the norm, we have: $(\epsilon_{hij} a_i u_j)^2 = 2a^2 K - P^2$. We therefore write: $\epsilon_{klm} a_l u_m = b_k \sqrt{2a^2 K - P^2} = b_k \sqrt{2K} \sqrt{a^2 - a_T^2} = b_k \sqrt{2K} \sqrt{a_N^2}$.</footnote>

<footnote n="6">⁶$D_{ij}D_{ij} = \frac{\gamma^2}{4} \left( a_T^2 \underbrace{\delta_{ij} \delta_{ij}}_{3} + a_N^2 \underbrace{\epsilon_{ijk} b_k \epsilon_{ijl} b_l}_{2\delta_{kl} b_k b_l} \right)$, $\delta_{kl} b_k b_l = 1$ since $b$ is a unit vector and with $a^2 = a_T^2 + a_N^2$, we obtain the result.</footnote>
To go further we must now specify the terms $\Pi$ and $\Sigma$ used for the stochastic process for $\epsilon$. Various models for the dissipation have been proposed. Pope and Chen [203] proposed a simple model based on the exponential of an Orstein-Uhlenbeck process (see (2.18)). Here, we rely on the model proposed in Refs. [54, 195]. This non-Markovian log-normal model presents a logarithmic decrease in the correlation of $\epsilon$ which is consistent with the idea of a turbulent cascade and a multiplicative process (see (2.22)), unlike the Pope and Chen model which gives an exponential decrease, see also the discussion in Ref. [146]. As presented in (2.24)-(2.25), the drift and diffusion terms are written, respectively, as:

$$ \Pi = \frac{1}{\tau_{\epsilon}} \left( -\ln \frac{\epsilon}{\langle \epsilon \rangle} + \frac{\sigma^2}{2\Lambda^2} \left( \frac{\tau_{\epsilon}}{\tau_c} - \Lambda^2 \right) + \frac{\sigma}{\Lambda} \hat{\Gamma} \tau_{\epsilon} \right), \quad (2.71) $$

and

$$ \Sigma = \sqrt{\frac{\sigma^2}{\Lambda^2 \tau_c}}, \quad (2.72) $$

with $\sigma^2$ the variance of the logarithm of $\epsilon$, $\tau_{\epsilon}$ the correlation time of $\epsilon$, $\tau_c$ the regularization timescale of the process (taken equal to the Kolmogorov dissipative time $\tau_{\eta}$), $\Lambda^2$ a normalization factor, and $\hat{\Gamma}$ the convolution of the Wiener increments with a temporal kernel, ensuring the non-markovian property of the process. In the process for $\epsilon$ proposed by Refs. [54, 195], the latter corresponds to a fractional Gaussian noise with 0 Hurst exponent [159] regularized at scale $\tau_c$. The expression of the convolution kernel proposed by Refs. [54, 195], (see (2.23)), applies to a scalar noise since the dissipation rate is a scalar, whereas the acceleration model involves a vectorial noise. Therefore, the kernel in $\hat{\Gamma}$ includes a projection in order to apply to the vectorial Wiener increments:

$$ \hat{\Gamma} = -\frac{1}{2} \int_{-\infty}^{t} \frac{1}{(t - s + \tau_c)^{3/2}} \mathcal{P}_j dW_j(s). \quad (2.73) $$

By proceeding in a similar way as Ref. [195], the projection operator is obtained by identification between the diffusion terms of Eqs. (2.55) and (2.59):

$$ \mathcal{P}_j = \frac{2a_i D_{ij}}{\gamma a^2 \Sigma} = \frac{\sqrt{a_T^2 - a_T}}{a^2} a_j + e_j \quad (2.74) $$

where we have used the relation recall in footnote 5 and where $e_j$ is the unit vector tangent to the trajectory, $e_i = u_i/\sqrt{2K}$. It is interesting to remark that the rotational part of the diffusion tensor induces an asymmetry of the projector between positive and negative power exchange (recall that $P = \sqrt{2K}a_T$). Indeed for $P \ge 0$, $\mathcal{P}_j = e_j$ while for $P < 0$ one has $\mathcal{P}_j = (1 - 2p^2)e_j - 2p\sqrt{1 - p^2}b_j$ with $p = P/\sqrt{2K}a^2$. In both cases, as it can be readily checked, $\mathcal{P}$ is a unit vector.

Substituting the expressions (2.71) and (2.72) for $\Pi$ and $\Sigma$ in Eq. (2.70) we have

$$ a_i M_i = a^2 \left( \frac{\alpha}{2\langle K \rangle} P - \frac{\gamma}{2\tau_{\epsilon}} \left( \ln \frac{\epsilon}{\langle \epsilon \rangle} + \frac{1}{2}\sigma^2 - \frac{\sigma}{\Lambda} \hat{\Gamma} \tau_{\epsilon} \right) - a_T^2 \frac{\gamma^2}{8} \frac{\sigma^2}{\Lambda^2 \tau_c} \right). \quad (2.75) $$

According to Eq. (2.52), we can write:

$$ \ln \left( \frac{\epsilon}{\langle \epsilon \rangle} \right) = \frac{1}{\gamma} \left( \ln \frac{a^2}{a_{\eta}^2} - \ln C - \alpha \frac{K}{\langle K \rangle} \right), \quad (2.76) $$
which gives, once substituted into Eq. (2.75),

$$a_i M_i = a^2 \left[ \frac{\alpha}{2\langle K \rangle} \left( P + \frac{K}{\tau_{\varepsilon}} \right) - \frac{1}{2\tau_{\varepsilon}} \left( \ln \left( \frac{a^2}{a_{\eta}^2} \right) - \underbrace{\ln C + \frac{\gamma}{2}\sigma^2 - \gamma \frac{\sigma}{\Lambda} \hat{\Gamma} \tau_{\varepsilon}}_{-\hat{\Gamma}_*} \right) \right] - \frac{a_T^2}{\tau_c} \underbrace{\frac{\gamma^2}{8} \frac{\sigma^2}{\Lambda^2}}_{\sigma_*^2} \quad (2.77)$$

To simplify the notations, we have introduced $\hat{\Gamma}_* = \gamma \frac{\sigma}{\Lambda} \hat{\Gamma} \tau_{\varepsilon} + \ln C - \frac{\gamma}{2} \sigma^2$ and $\sigma_*^2 = \frac{\gamma^2}{8} \frac{\sigma^2}{\Lambda^2}$.

It is interesting to notice that in Eq. (2.77) the terms $P + \frac{K}{\tau_{\varepsilon}} = \frac{dK}{dt} + \frac{K}{\tau_{\varepsilon}}$ acts as a penalization leading the correlation of the kinetic energy to decay exponentially.

We then propose for $M_i$ an expression compatible with Eq. (2.77). Proceeding by identification, we have the following relation:

$$M_i = \frac{\alpha}{2\langle K \rangle} \left( \lambda a_i P + (1-\lambda)a^2 u_i + a_i \frac{K}{\tau_{\varepsilon}} \right) \\ -a_i \left( \ln \left( \frac{a^2}{a_{\eta}^2} \right) - \hat{\Gamma}_* \right) \frac{1}{2\tau_{\varepsilon}} \\ - \frac{\sigma_*^2 a_T^2}{\tau_c a^2} a_i + B_i, \quad (2.78)$$

where we have introduced the vector $B_i$, such that $B_i a_i = 0$ as well as the factor $\lambda$ that both account for the indeterminacy inherent to the inverse projection. By assuming again that there are only two vectors at our disposal, we can take $B_i = \frac{\alpha}{2\langle K \rangle} \lambda' (Pa_i - a^2 u_i)$ by introducing the factor $\lambda'$. Note that from the point of view of the projection, the factors $\lambda$ and $\lambda'$ are arbitrary in the sense that the scalar product of $a_i$ and Eq. (2.78) gives Eq. (2.77) whatever their values. We can nevertheless notice that the terms involving $\lambda$ and $\lambda'$ can be combined, and, by noting $c_u = \lambda + \lambda'$, we get:

$$M_i = \frac{\alpha}{2\langle K \rangle} \left( a_i(c_u P + \frac{K}{\tau_{\varepsilon}}) - (c_u - 1)a^2 u_i \right) \\ -a_i \left( \ln \left( \frac{a^2}{a_{\eta}^2} \right) - \hat{\Gamma}_* \right) \frac{1}{2\tau_{\varepsilon}} \\ - \frac{\sigma_*^2 a_T^2}{\tau_c a^2} a_i. \quad (2.79)$$

We can notice that the terms of the first line correspond to the coupling with the velocity, those of the second take into account the log-normal and non-Markovian character of the dissipation and the last term is due to the rotational part of the diffusion tensor. The diffusion term (2.67) becomes, by using expression (2.72),

$$D_{ij} = \sqrt{\frac{2\sigma_*^2}{\tau_c}} \left[ \sqrt{a_T^2} \delta_{ij} + \sqrt{a_N^2} \epsilon_{ijk} b_k \right]. \quad (2.80)$$

We have thus specified our stochastic model for the dynamics of a fluid particle. It is given by Eqs. (2.56), (2.57), (2.79) and (2.80).
## 2.4.2 Parameters and numerical approach

From a dimensional point of view, to determine the physical parameters of the stochastic model, one must specify time and velocity scales as well as a Reynolds number. This amounts for example to imposing the average kinetic energy $\langle K \rangle$, the average dissipation rate $\langle \varepsilon \rangle$ and the viscosity $\nu$. From these physical parameters, we calculate $a_{\eta}^2 = \langle \varepsilon \rangle^{3/2} \nu^{-1/2}$, $\tau_{\eta} = \langle \varepsilon \rangle^{-1/2} \nu^{1/2}$. We can also get the Reynolds number based on the Taylor scale $Re_{\lambda} = u'\lambda/\nu = 2\sqrt{15}/3 \langle K \rangle / \sqrt{\langle \varepsilon \rangle} \nu$ with $u' = \sqrt{2\langle K \rangle}/3$ and $\lambda^2 = 15\nu u'^2 / \langle \varepsilon \rangle$. We then deduce the Lagrangian integral times scale $\tau_L$ as $\tau_L = 0.08 Re_{\lambda} \tau_{\eta}$ from the DNS results reported by Refs. [101, 238].

The parameter $\sigma^2$ is estimated using the relation given by Ref. [271]: $\sigma^2 \approx 3/8 \ln Re_{\lambda}/R_c$ with $R_c \approx 10$ compatible with the prediction of Kolmogorov and Obhoukov [128, 192]. As mentioned in Refs. [128] and [178], the specific value of $R_c$ is depending on the large-scales. Since the influence of the large-scales is neglected in our modeling (see Sec. 2.3.4), we choose in the following simply $\sigma^2 \approx 3/8 \ln Re_{\lambda}$. We set as well $\alpha = 1/3$ and $\gamma = 3/2 + \beta$ with $\beta = -1/\ln Re_{\lambda}$ in accordance with the results of the DNS presented above. The prefactor $C$ is computed as $C = c_0 AB$ where $A = \left(1 - \frac{2}{3}\alpha\right)^{3/2} \approx 0.686$, $B = 17.1 - 54.7/\ln Re_{\lambda}$ as determined by DNS. The term $c_0$ is introduced so that the predicted acceleration variance follows Eq. (2.44), as one would expect from the construction of the stochastic model, despite the fact that we take $\sigma^2 = 3/8 \ln Re_{\lambda}$ instead of $\sigma^2 = 3/8 \ln Re_{\lambda}/R_c$. Consequently, we have $c_0 = (1/Re_c)^{9/64+3\beta/8(1+\beta/2)}$.

For simplicity we have used $\tau_c = \tau_{\eta}$ and $\tau_{\varepsilon} = \tau_L$. From $\tau_{\varepsilon}$ and $\tau_c$ we calculate the value of the normalizing constant $\Lambda$ as explained after (2.23). Finally, for the parameter $c_u$, which is the only free parameter of the model, we have determined numerically that with $c_u = 5.22$ the ratio $K/\langle K \rangle$ is 1 on average for all values of the Reynolds number.

The sample paths of this model are obtained by numerical integration of the stochastic differential equation. Numerical integration is made with an explicit Euler scheme by taking a time-step $dt = \tau_{\eta,min}/100$ with $\tau_{\eta,min} = \sqrt{\nu/\varepsilon_{max}}$, an estimation of the minimum dissipative timescale likely to happen during the simulation. This is estimated from the log-normal distribution of the dissipation: $\tau_{\eta,min} = \tau_{\eta} \exp(-x\sigma^2/2 + \sigma^2/4)$, with $x = 6$ by considering that the probability that a random number following the normal distribution reaches a value of 6 standard deviation is sufficiently low (see Eq. (2.17)).

For the calculation of the convolution term $\hat{\Gamma}$ appearing in Eq. (2.79), we propose in appendix of [274] an efficient algorithm. Note also that a simple Python script presenting the algorithm used to integrate the proposed stochastic model is available in supplemental material of Ref. [274].

## 2.4.3 Results

We show in Fig. 2.11 a realization of this process for $Re_{\lambda} = 1100$. We see the temporal evolution of the components of acceleration and velocity. There is a very intermittent acceleration with an alternation of periods in which the acceleration of the fluid particle is almost zero with phases of very intense activity. This results in fluid-particle trajectories, obtained by integration of the velocity $x_i = \int u_i(t)dt$, in long quasiballistic periods with typical length of the order of the integral scale ($L \approx \langle K \rangle^{3/2}/\langle \varepsilon \rangle$) and short-term disruptions during which the trajectory rolls up on itself.
<figure><img src="image_17.png" /><figcaption>Figure 2.11: A realization of the stochastic process for $Re_λ = 1100$. Top left: evolution of the acceleration with time, $a_x$: red, $a_y$: green, $a_z$: blue, $|a|$: black. Bottom left: evolution of the velocity with time, $u_x$: red, $u_y$: green, $u_z$: blue, $|u|$: black. Right: 3D trajectory of a fluid particle for a duration of 100 $τ_L$.</figcaption></figure>

We have simulated the stochastic model for 15 different Reynolds numbers between $Re_λ = 70$ and 9000. In each case, we have computed 26,000 realizations. The simulations are carried out over a period of $120τ_L$, over which we exclude an initial transitional regime of $20τ_L$ for the calculation of the statistics. In all cases, the initial value of the components of acceleration and velocity are sampled from the normal distribution having a standard deviation of $10^{-9}a_η$ for the acceleration and $10^{-9}\sqrt{2\langle K \rangle}/3$ for the velocity. We can indeed notice from Eqs. (2.79) and (2.80) that, if the acceleration is exactly zero, the stochastic model predicts that the acceleration would remain so. However, it should be noted that this event has a zero probability, and that for arbitrarily small, but nonzero, accelerations, the model presents evolution towards a nontrivial stationary state. This is illustrated in Fig. 2.12, which presents the temporal evolution of the variance of the velocity and of the acceleration for $Re_λ = 1100$, calculated from all the realizations.

Figure 2.13 shows the evolution with the Reynolds number of the mean kinetic energy in the stationary state. In this figure, we see that the average kinetic energy is equal, within the statistical convergence, to the value prescribed to the model. We note that the value of the average kinetic energy is directly related to the value of the parameter $c_u$ in Eq. (2.79) as mentioned above.

Regarding the variance of the acceleration, we expect, by construction of the stochastic model, that the predicted value follows the log-normal relation (2.44). We observe in Fig. 2.10, that it is indeed the case, with only slight deviations for the largest Reynolds numbers which are attributed to numerical errors. We recall that the underestimation of the acceleration variance at small Reynolds numbers compared to the DNS or Eq. (2.45) stems from the fact that the model is based on the relation (2.47) in which the effect of low dissipative and large-scale structures are neglected (see the discussion in Sec. 2.3.3). This simplification enables us to obtain the analytical formulation of the model proposed here.
<figure><img src="image_18.png" /></figure>

Figure 2.12: Evolution of the acceleration and velocity variance with time for the stochastic model at $Re_λ = 1100$ starting from an initial condition for the acceleration and velocity with very small magnitude. Comparison in dashed gray line with the expected values: Eq. (2.44) for the acceleration and the prescribed value of $u' = \sqrt{2\langle K \rangle}/3$ for the velocity.

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

Figure 2.13: Evolution with the Reynolds number, in the stationary regime, of the kinetic energy obtained by the stochastic model normalized by the prescribed kinetic energy $\langle K \rangle$.

Figure 2.14 compares the autocorrelation of the components of the acceleration and of its norm calculated from the stochastic model for $Re_λ = 400$ with the calculations from the DNS of [140, 21]. It can be seen that the characteristic times of these two quantities are very different and that it is in good agreement with the DNS. It should be mentioned that the scale separation between the components and the norm results from the rotational part of the diffusion tensor. Indeed, no scale separation is found when $c_R$ is set to zero in Eq. (2.66) (corresponding then to a diagonal diffusion tensor). We see in Eq. (2.67) that considering this rotational part, leads to the decomposition of the acceleration into its normal and tangential component. The former is associated with the intense rotation that rapidly changes the acceleration direction, whereas the second is associated with the variation of the kinetic energy of the particle.

Figure 2.14 also presents the evolution of the autocorrelation coefficient of the velocity components and of the power. It can be seen here also that the agreement with the DNS is relatively good. In Fig. 2.14, we also show the evolution of the characteristic correlation times for these four quantities with Reynolds numbers in the range $Re_λ = 70 - 9000$ as
predicted by the stochastic model. The characteristic correlation time for the velocity, the acceleration norm, the acceleration components and the power are $\tau_u = \int \rho_{u_i}(\tau)d\tau$, $\tau_{|a|} = \int \rho_{|a|}(\tau)d\tau$, $\tau_{a_i} = \int |\rho_{a_i}|(\tau)d\tau$ and $\tau_P = \int |\rho_P|(\tau)d\tau$. It can be seen in Fig. 2.14 that the scales for the norm of the acceleration and for the velocity normalized by $\tau_L$ remains quasiconstant with the Reynolds number and that the ratio between the correlation scale for the velocity and $\tau_L$ is of order 1. Note that the characteristic time entering the model formulation is $\tau_\epsilon$ (the correlation time of the dissipation rate following the path of a fluid particle). For the calculation of the model, we simply set $\tau_\epsilon = \tau_L$ arguing that the two quantities should be closed. It is therefore interesting to remark that the integral time of the velocity is very close to the prescribed one $\tau_u \approx \tau_L$. Regarding the correlation scales for a component of acceleration and for the power normalized by $\tau_L$, they both present a variation close to $1/\text{Re}_\lambda$, as expected.

We also show in Fig. 2.14 autocorrelation coefficient of $a_i$, $(a^2)^{1/2}$, $P = a_i u_i$ and $u_i$ for Reynolds numbers in the range $\text{Re}_\lambda = 70-9000$ obtained from the model. It is seen that, when the time shift is normalized by the corresponding integral timescale, the correlation coefficients of the power and of the acceleration component remains nearly unchanged with the Reynolds number. We observe as well that the shape of the autocorrelation obtained from DNS is well reproduced, although the decay predicted by the model is too fast at very short time lag. This is attributed to the fact that the dissipative region is only taken into account in the model via the cutoff $\tau_c = \tau_\eta$ of the kernel $\hat{\Gamma}$. We observe that the correlation for the acceleration norm presents a logarithmic decrease, reflecting the absence of characteristic time for its evolution. As expected, the correlation norm exhibits a lower slope as the Reynolds number increases. This is directly attributed to the use of the non-Markovian process of Ref. [54] for the dissipation rate, which proposes a logarithmic evolution of the autocorrelation in agreement with the underlying model of the turbulent energy cascade as discussed in section 2.1 (see figure 2.3).

The shape of the velocity correlation from the model is overall close to the DNS. At small $\tau$, it presents some dependence on the Reynolds number, while at large time shift (i.e. $\tau$ of the order of $\tau_L$) the correlation decreases exponentially, as it can be seen in the inset of Fig. 2.14(f), in agreement with DNS and experiments. The exponential relaxation results from the terms $P + K/\tau_\epsilon = dK/dt + K/\tau_\epsilon$ appearing in the drift part of the stochastic model (2.78). It is interesting to remark that the presence of this term in the model is a direct consequence of the exponential dependence of the conditional acceleration variance on the kinetic energy (2.47). This term leads to the Reynolds number dependence on the velocity correlation observed at small $\tau$, which is connected to the logarithmic decorrelation of the acceleration norm, to vanish at large $\tau$ at which it relaxes exponentially. This suggests therefore anomalous scaling at intermediate time lag.

We show in Fig. 2.15 the velocity spectrum for $\text{Re}_\lambda$ between 400 and 9000, which we compare with the DNS of Ref. [21] for $\text{Re}_\lambda = 400$. We see a good agreement between the DNS and the stochastic model. For higher Reynolds numbers, we clearly see that a power-law behavior develops at intermediate scales. We see that the slope of the power-law deviates from the Hinze spectra [256] predicted by dimensional arguments similar to those presented by Kolmogorov, with spectra less stiff than $\omega^{-2}$. This shows that the proposed stochastic model leads to an anomalous scaling that reflects the persistent influence of the Reynolds number in the inertial-scales. We further notice that the slope that develops at intermediate scales are close to $-2+0.14$, where 0.14 is the exponent of the asymptotic power-law of the acceleration variance with the Reynolds number determined in Eq. (2.46)
(see also Fig. 2.10). We see here a confirmation of the relation between the acceleration scaling and the anomalous scaling of the velocity spectra proposed by Ref. [75].

We present in Fig. 2.16 the PDFs of the velocity and of the acceleration for $Re_λ = 400 \sim 9000$, as well as the comparison with the DNS of Ref. [21]. First, we find that the velocity distribution is very close to a Gaussian distribution for all Reynolds numbers, while the acceleration presents a much more stretched distribution. For $Re_λ = 400$ the acceleration PDF is in very good agreement with the DNS, and, the model predicts an increase of the stretching of the tails with increasing the Reynolds number. We also show in this figure the PDF of the velocity increments for different time shifts $\delta_τ u_i = u_i(t + \tau) - u_i(t)$ at $Re_λ = 400$. We observe that the distribution gradually returns to a Gaussian distribution as the time shift increases, and that at each time shift the agreement with the DNS of Ref. [21] is very good. This is confirmed by the presentation of the flatness of the velocity increments for $Re_λ = 400 \sim 9000$, which reflects the strongly non-Gaussian behavior on small-scales which decreases to 3 for the larger-scales. Here also we notice a good agreement with the DNS of [21] for $Re_λ = 400$. We also show in the inset, a quasi-linear increase of the flatness of the acceleration with the Reynolds number.

Finally, in Fig. 2.17 we show the second and third moments of the power $P = a_i u_i$. It is observed that the increases of both moments with the Reynolds number are in close agreement with the power-law supported by the DNS results of Ref. [268], $\langle P^2 \rangle / \langle \epsilon \rangle^2 \sim Re_\lambda^{4/3}$ and $-(\langle P^3 \rangle / \langle \epsilon \rangle^3) \sim Re_\lambda^2$. Clearly, the third order moment is negative, meaning that the time irreversibility of the dynamics of a fluid particle in a turbulent flow is correctly reproduced by the proposed stochastic model. The skewness of the power, $S = \langle P^3 \rangle / \langle P^2 \rangle^{3/2}$, seems to converge to -0.5 as the Reynolds number increases, as reported in Ref. [268].
<figure><img src="image_20.png" /><figcaption>Figure 2.14: (a) Evolution of the autocorrelation of $a_i$ (black), $(a^2)^{1/2}$ (red), $u_i$ (blue) and $P = a_i u_i$ (green) from the stochastic model for $Re_\lambda = 400$ and comparison with the DNS data from Ref. [21] in dashed lines. (b) Evolution of the integral timescale of $a_i$ (black), $(a^2)^{1/2}$ (red), $u_i$ (blue) and $P = a_i u_i$ (green) normalized by $\tau_L$ with the Reynolds number. (c,d,e,f) Evolution of the autocorrelation of $a_i$, $(a^2)^{1/2}$, $P = a_i u_i$ and $u_i$ respectively, for $Re_\lambda = 400, 567, 800, 1130, 1600, 2263, 3200, 4526, 6400$ and $9051$ from orange to black and comparison with the DNS data from Ref. [21] in dashed lines. In these plots the time lag is normalized by the corresponding integral timescale. For panel f, inset: logarithmic scaling of the y-axis and comparison with $\exp(-\tau/\tau_u)$ in dotted line.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 2.15: Velocity spectra from the stochastic model for $Re_λ = 400$ to $Re_λ = 9000$ from orange to black and comparison with the DNS data from [21] at $Re_λ = 400$ (gray dashed line), with the Hinze spectra $\omega^{-2}$ (gray dot-dashed line), and with the power-law with anomalous exponent $\omega^{-2+9/64}$ (gray dotted lines).</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure 2.16: PDF of $a_i$ (top left) and comparison with the DNS data from Ref. [21], and PDF of $u_i$ (top right) and comparison with the normal distribution, for $Re_λ$ = 400 to 9000 from orange to black. (Bottom left) PDF of the velocity increments for various times shift: $τ/τ_η$ = 610, 200, 70, 22, 7.3, 2.4, 0.8 and 0.25 each time lag is shifted upward by one decade for $Re_λ$ = 400 comparison with the PDF of the acceleration from the model (red) and from the DNS of Ref. [21] (gray) and the PDF velocity (blue). (Bottom right) Evolution of the Flatness of the velocity increments versus the time shift for $Re_λ$ = 70 to 9000 from orange to black and evolution of the flatness of the acceleration with the Reynolds number and comparison with the linear law in the inset.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure 2.17: Evolution of $\langle P^2 \rangle / \langle \varepsilon \rangle^2$ (continuous black line) and $-\langle P^3 \rangle / \langle \varepsilon \rangle^3$ (dashed black line) with the Reynolds number, and comparison with the power-laws $Re_\lambda^{4/3}$ and $Re_\lambda^2$.</figcaption></figure>
## 2.5 Discussion and perspectives

We have analyzed the behavior of the acceleration statistics conditioned on both local dissipation rate and local kinetic energy, which to our knowledge have not been considered before. We have reported that the doubly conditional variance is proportional to the acceleration variance conditional on the dissipation rate solely, with the proportionality factor depending exponentially on the kinetic energy: $\langle a^2|\epsilon, K \rangle = A \exp(\alpha K/\langle K \rangle) \langle a^2|\epsilon \rangle$. For large enough Reynolds number, we show that $A = (1 - 2\alpha/3)^{3/2}$ and we proposed that the $\alpha$ coefficient is independent of the Reynolds number and its value $\alpha = 1/3$ was obtained from the DNS.

This expression shows a direct effect of the kinetic energy, a large-scale quantity, on the Lagrangian acceleration. Furthermore since the argument of the exponential depends on $K/\langle K \rangle$, not on a local Reynolds number, it suggests a kinematic effect for the acceleration which may be due to the nonlocality of the pressure. More specifically, these effects can be due to the interaction of vorticity and strain [65]. In case of persistent large-scale strain, intense vorticity tube would be generated and align with the principal direction of the strain [175]. It was shown that such vortical structure can produce significant acceleration in the direction of the vorticity [17, 145]. Anyway, although the proper physical mechanism leading to the exponential dependence of the acceleration on the kinetic energy deserves further studies, it is an additional effect to the influence of the large-scales on the acceleration through the intermittency of the dissipation rate. To study this later effect, we subsequently have proposed to account for the Reynolds number dependence of the acceleration variance conditional on the dissipation rate within the intermediate asymptotic framework [13] leading to: $\langle a^2|\epsilon \rangle = a_{\eta}^2 B(\epsilon/\langle\epsilon\rangle)^{3/2+\beta}$ for $\epsilon \gg \langle\epsilon\rangle$ with $B$ and $\beta$ depending logarithmically on the Reynolds number as the signature of the intermittency and the persistence of viscous effects. Further, we advance an expression for the conditional acceleration variance valid for the whole range of fluctuations of $\epsilon$ by accounting for the dominant effect of the large-scale structures in low dissipative regions (2.41). From this finding we determine the evolution of the unconditional acceleration variance with the Reynolds number (2.45) and show that it is in good agreement with DNS, which gives another empirical validation of the incomplete similarities assumption used to obtain these results.

Eventually, for large Reynolds numbers, we propose to express the doubly conditional variance as $\langle a^2|\epsilon, K \rangle = C a_{\eta}^2 \exp(\alpha K/\langle K \rangle + \gamma \ln \epsilon \langle\epsilon\rangle)$, $\gamma = 3/2 + \beta$, which can be viewed as the results of a multiplicative process for the acceleration. Such a process can be interpreted as a momentum fluctuation cascade that includes kinematic effects by eddies all along the turbulence spectrum.

Based on these results we propose a 3D stochastic model for the dynamics of a fluid particle that reproduce the essential features of the Lagrangian dynamics observed from DNS and experiments. To obtain such a model, (i) we have assumed, inline with the Kolmogorov universality hypothesis, that the dynamics can be described as a set of stochastic differential equation $da_i = M_i dt + D_{ij} dW_j$; $du_i = a_i dt$, with $M_i$ and $D_{ij}$ depending on the velocity and acceleration along with Reynolds number dependent parameters. (ii) We used the doubly conditional acceleration variance presented first, to model the instantaneous relation of the dynamics between acceleration (or force), kinetic energy, and energy dissipation. This amounts to consider that the remaining degree of freedom can be discarded in procedure similar to an adiabatic elimination [84] as discussed by Ref. [49]. (iii) We
introduce a nondiagonal diffusion tensor along with a maximum winding hypothesis to ensure its physical realizability. (iv) We consider that the dissipation rate along the trajectory is given by the non-Markovian log-normal process proposed recently by Ref. [54], giving logarithmic correlation consistently with the turbulent cascade picture. The model is closed by using the relation $DK/Dt = P = a_i u_i$. For the model, it implies dependence of $K$ on $\epsilon$ through the dependence of $a^2$ on $\epsilon$. This can be interpreted as feedback of the small scales on the large scales. On the other hand, the influence of the large-scales on the small-scales is accounted for in the model through the intermittent cascade model for the dissipation rate. With these four hypotheses, we obtain the model given by Eqs. (2.56), (2.57), (2.79), and (2.80) which reads:

$$da_i = \left[ \frac{\alpha}{2\langle K \rangle} \left( a_i(c_u P + \frac{K}{\tau_{\epsilon}}) - (c_u - 1)a^2 u_i \right) - a_i \left( \ln\left(\frac{a^2}{a_{\eta}^2}\right) + \hat{\Gamma}_* \right) \frac{1}{2\tau_{\epsilon}} - \frac{\sigma_*^2 a_T^2}{\tau_c a^2} a_i \right] dt \\ + \sqrt{\frac{\sigma_*^2}{\tau_c}} \left[ \sqrt{a_T^2} \delta_{ij} + \sqrt{a_N^2} \epsilon_{ijk} b_k \right] dW_j. \qquad (2.81)$$

We show that the proposed model predicts Lagrangian dynamics presenting non-Gausssianity, long-range correlations, anomalous scaling and time irreversibility. Moreover statistics obtained from the stochastic model are in good agreement with the DNS.

In Eq. (2.81) the term proportional to $\alpha$, which follows directly from the exponential dependence of the conditional acceleration on the kinetic energy, involves the coupling between velocity and acceleration and leads to the exponential relaxation of the velocity correlation for large time lag along with one-time Gaussian distribution for the velocity. Introducing a rotational part in the diffusion tensor naturally leads to decomposition of the acceleration vector into its tangential part and its normal components. Since the normal part is associated with the curvature of the trajectory, the rotational part of the diffusion leads to the emergence of a time-scale separation between the correlation of the norm and the components of the acceleration. The term associated with the non-Markovianity of the dissipation along with the rotational part produce irreversible dynamics, as seen by the skewness of the exchanged power and ensures a scale separation between velocity and acceleration. These three points can be easily checked, by taking $\alpha = 0$ or $c_R = 0$ in Eq. (2.66) or by using for $\Pi$ the Markovian log-normal dissipation model proposed by Ref. [203] rather than the non-Markovian one of Ref. [54].

It is worth noting that from the conditional acceleration statistics obtained from DNS of the Navier-Stokes equation, it is possible to establish, in a fairly natural way, that is to say without using any other hypothesis than the cascade picture, a link between the refined Kolmogorov assumption and the dynamics of fluid particles. It would be interesting to analyze further the stochastic equation to demonstrate the irreversibility of the dynamics, the emergence of anomalous scaling or to study the geometry of the particle trajectory e.g. its curvature and torsion, as well as to further test the conditional statistics between the acceleration and the velocity. Also interesting could be the improvement of the modeling of the high frequency part of the spectrum. Indeed, the dissipative part of the spectrum is not well reproduced by the model of Ref. [54] which intends to model the dissipation rate in the inertial range.

To simplify the construction of the model, we have not taken into account the nonlocal effects of the largest structures of the flow, arguing that their effect vanish as the Reynolds number increases (term with $a_0^2$ in eq. (2.41)). Based on the relation (2.41) it is possible
to account for the large-scale in the stochastic modeling. However, since this term is de-
pendent on the Reynolds number, it is likely that it also depends on the flow configuration
and boundary conditions. The proposed stochastic model could be further generalized
to address shear flows [17] and improve Reynolds-averaged simulations [200, 109]. This
model could be used among other things to improve the calculation of the dynamics of
a dispersed phase with the large eddy simulation (LES) approach [279, 92, 280]. Finally,
let us mention that an interesting extension could be the coupling of the proposed model
with stochastic model for the evolution of the velocity gradients as proposed in Refs.
[89, 168, 114, 195].
# 3 Particles and turbulence

## Contents

<table><tr><td>3.1</td><td>The Euler-Lagrange approach</td><td>48</td></tr><tr><td>3.2</td><td>Force and acceleration statistics for particles in turbulence</td><td>49</td></tr><tr><td>3.2.1</td><td>Small particles</td><td>50</td></tr><tr><td>3.2.2</td><td>Large particles</td><td>63</td></tr><tr><td>3.3</td><td>Stochastic modeling of the dynamics of a particle and application to LES</td><td>66</td></tr><tr><td>3.3.1</td><td>The large eddy simulation approach</td><td>66</td></tr><tr><td>3.3.2</td><td>Large eddy simulations for two-phase flows</td><td>68</td></tr><tr><td>3.4</td><td>Turbulence induced by a swarm of rising bubbles</td><td>77</td></tr><tr><td>3.4.1</td><td>Simulation of the bubble swarm</td><td>78</td></tr><tr><td>3.4.2</td><td>Analysis of the bubble swarm</td><td>84</td></tr><tr><td>3.4.3</td><td>Summary</td><td>94</td></tr><tr><td>3.5</td><td>Perspectives</td><td>95</td></tr></table>

In this section, I present some works on the modeling of turbulent flows with dispersed phases. Several problems are addressed: the a priori estimation of the statistics of the forces experienced by a particle in turbulent flows, the stochastic modeling of the instantaneous value of the forces and its application to the large eddy simulations (LES), and the turbulence induced by the rising of a cloud of bubbles.

One speaks of a dispersed phase when the evolution of each of these components, typically a bubble, a drop or a solid particle, can be obtained from the Newton's equation, by applying the sum of the forces to the center of mass of the object

$$ \frac{dx_p(t)}{dt} = v_p ; \quad m_p a_p = m_p \frac{dv_p(t)}{dt} = F_p . \qquad (3.1) $$

We can decompose this resulting force by separating the hydrodynamic forces $F_{\text{fluid}}$ resulting from the exchange of momentum between the continuous phase (the carrier phase) and the particle, from the external forces $F_{\text{ext}}$, typically the weight: $F_p = F_{\text{fluid}} + F_{\text{ext}}^1$. From a formal point of view, the hydrodynamic force on a particle is determined by integrating the fluid stresses $\Sigma$ (normal stress and shear) on the surface of the particle $F_{\text{fluid}} = \int_{S_p} \Sigma \cdot n_p ds$.

<footnote n="1">¹There is also the force $F_I$ which translates the direct interactions between objects of the dispersed phase, typically the collisions, but it can also be an electromagnetic force if the particles are charged or if the medium conducts electricity...</footnote>
## 3.1 The Euler-Lagrange approach

In the Euler-Lagrange approach, the fluid is described in an Eulerian way and each object of the dispersed phase in a Lagrangian way, based directly on Newton's equation. One of the limitations of this approach is that it requires the resolution of a differential equation for each object of the dispersed phase, which is impossible if they are too many. The alternative, logically called Euler-Euler approaches, is to give a continuous description of the dispersed phase, by proposing statistical closures [33, 142, 245, 260, 77, 71, 78, 170, 211]. However, the main difficulty remains, that of estimating the resulting force on each particle. Indeed, the calculation of the stress that the fluid exerts on a particle requires to know precisely (i.e. at a scale smaller than the particle) the velocity and pressure field. When the particles are small, or when the boundary layers developing around them are small, in comparison of the other characteristic scales of the flow, it becomes impractical. The Euler-Lagrange approach is based on an effective description of the fluid at a scale $l$ larger than that of the particle. That is to say that the details of the velocity field are removed: $u_f \to \bar{u}_f$, and with respect to this filtered velocity field, the particle is pointwise. The resultant of the hydrodynamic forces is then estimated by a model which depends on the filtered velocity field, without requiring to have the full details of the velocity field. Formally, we can introduce a functional $\mathcal{F}$ giving the force on a particle from the (filtered) velocity field of the fluid, the velocity of the particles $u_p(t)$, their positions $x_p(t)$, as well as their orientation $\theta_p(t)$ and rotation speed $\Omega_p(t)$ ...

$$F_{\text{fluid}} = \mathcal{F}[\bar{u}_f(x, t); u_p(t), x_p(t), \Omega_p(t), \theta_p(t) \dots] . \quad (3.2)$$

In general, the expression of such a functional is unknown. As an example, in the case of a fixed object in a uniform and stationary flow, we can make $l$ go to infinity, and obtain the force (or at least its average value in time) from the value of the velocity far upstream and its orientation relative to the object by relying on relations (most often empirical) for the drag and lift coefficients.

For a turbulent flow which is obviously neither stationary nor uniform, the use of the upstream velocity does not make sense. The method generally adopted is rather to use the value of the velocity field at the position of the particle. Indeed, the field $\bar{u}_f$ is defined everywhere, including at the particle position. This is fully justified if we restrict ourselves to particles with a diameter $d$ much smaller than the smallest scale of the turbulence (the dissipative scale of the turbulence $\eta$) $d \ll \eta$. In this limit, one can consider that at an intermediate scale $l$, $\eta \gg l \gg d$, the velocity field $\bar{u}_f(x, t)$ is quasi-uniform<fnref n="2" /> at scale $l$. With this scale separation, the expression of the fluid force must depend only on the local value of the fluid velocity field and its derivatives described at scale $l$.

In this way, we can write that the hydrodynamic force must be a functional of the value of the velocity and its derivatives at the position of the particle (and not of the whole field as in (3.2)):

$$F_{\text{fluid}} = f[\bar{u}_f(x=x_p(t),t) - u_p(t), \nabla \bar{u}_f|_{x_p}, D_t \bar{u}_f|_{x_p}, d_t u_p(t)] . \quad (3.3)$$

Since the force must be Galilean invariant, only the relative velocity of the fluid and the particle can be used to determine the force. For the same reason, it is the material

<footnote n="2">²This uniformity assumption also implies that we consider an isolated particle, or in practice sufficiently diluted suspension.</footnote>
derivative of the velocity of the fluid that must intervene. Note that for simplicity, we
have omitted the terms involving the orientation and rotation of the particle.

The functional (3.3), and thus the expression of the hydrodynamic force, is only well
known for spherical objects moving in a quasi-uniform flow with a very small relative
velocity, so that their Reynolds number is small: $Re_p = |\bar{u}_f - u_p|d/\nu \ll 1$. In this
case we can decompose the hydrodynamic force into drag force, history force (unsteady
effect of the drag)<fnref n="3" />, added-mass force and inertia force [86, 164, 152, 74]. When $Re_p$
is larger, it is generally assumed that this decomposition remains valid if one corrects
certain force expressions (e.g. inertia effect on the drag force) and adds other forces such
as lift forces (note that several effects are grouped under the name of lift force). At high
Reynolds numbers $Re_p$, the validity of the previous equation is more empirical than really
rigorously shown, see for example [171] which considers the case of a fixed particle in a
turbulent flow.

In practice the scale $\ell$ can be taken as the resolution scale of the numerical simulation
$\ell \approx \Delta x$ and $\bar{u}_f$ corresponds then to the velocity resolved by the simulation. If, in addition,
the mesh is fine enough to resolve the turbulent fluctuations, we have $\eta \gg \ell \gg d$. This
setting is sometimes described as quasi-DNS although the dynamics of the small particles
are obtained from a model for the hydrodynamic force. In Sec. 3.2, we discuss the
statistics of acceleration and particle forces in a turbulent flow obtained assuming such
scale separation.

In the case where $\ell \gg \eta$, the mesh cannot capture all the turbulent scales of the flow, we
then speak about Large Eddy Simulations. In this case the model for the hydrodynamic
force must take into account the fluctuations that are not resolved by the mesh. In Sec.
3.3 we discuss the stochastic model formulation for such an approach for small particles
$\ell \gg \eta \gg d$ as well as for large particles $\ell \gg d \gg \eta$.

Finally, we consider the case where the particle size is not negligible compared to the size
of the mesh (the scale separation $\ell \gg d$ is not verified anymore). In this case, due to the
presence of the particle, the velocity of the fluid $\bar{u}_f$ close to the particle can be altered
considerably [231, 86, 164]. Therefore the fluid velocity $\bar{u}_f$ at the particle position must
be corrected to reliably calculate the hydrodynamic force. This case can be encountered
when studying the effect of the particles on the carrier phase (this is called “Two-Way
coupling”), at least when one is interested in the effect of each particle more than in a
global effect resulting from a large mass or volume loading. We will discuss in Sec. 3.4 the
simulation of the turbulence induced by the rising of a cloud of bubbles and the correction
that we have proposed for the calculation of the force for bubbles generating a significant
wake.

## 3.2 Force and acceleration statistics for particles in turbulence

We present models for the variance of the acceleration of spherical particles in a homo-
geneous and isotropic turbulent flow. In this section, we consider that the particles are
isolated (negligible volume fraction) and that we can neglect the modifications they cause
on the flow of the carrier phase. We are interested in particles that can be heavy or

<footnote n="3">³Note that if history effects are neglected, the hydrodynamic force is given by the instantaneous value of the velocity and its derivatives at the position of the particle so (3.3) can be understood as a function (and not a functional).</footnote>
light with respect to the carrier phase and smaller or larger than the Kolmogorov scale. Three dimensionless numbers are involved in the description of this problem: the ratio of the particle size to the dissipative scale of the flow $d/\eta$, the density ratio $\rho_p/\rho_f$ and the Reynolds number of the flow, for example the Taylor scale Reynolds number $Re_\lambda$.

### 3.2.1 Small particles

We first consider the case of small particles in the sense that both $d < \eta$ and $Re_p < 1$. In this limit, it is reasonable to neglect lift and history force [153, 86, 164]. We can thus express the acceleration on these small particles as:

$$ \frac{du_p}{dt} = -\frac{u_p - u_f}{\tau_p} + \beta \frac{Du_f}{Dt}. \qquad (3.4) $$

The coefficient $\beta$ depends on the density ratio of the particle:

$$ \beta = \frac{1 + C_M}{\rho_p / \rho_f + C_M} \qquad (3.5) $$

with $C_M$ the added-mass coefficient which is $C_M = 1/2$ for a sphere in an infinite medium [153, 19]. For a particle of high density $\rho_p/\rho_f \gg 1$ we have $\beta \to 0$ and for an extremely light particle $\rho_p/\rho_f = 0$ (typically the case of a bubble) we have $\beta = 3$. The relaxation time of the particle $\tau_p$ is:

$$ \tau_p = (\rho_p / \rho_f + C_M) \frac{d^2}{18\nu}. \qquad (3.6) $$

The expression is given here for the case where there is no slip of the fluid at the interface of the object. In the case where there is free shear at the interface, it is appropriate to replace the factor 18 by 12 as in Ref. [280]. This last condition corresponds in some cases to a bubble with an uncontaminated interface, but we draw attention to the fact that the interfacial dynamic is very rich is that a bubble is much more complicated than a light particle with a slip condition [149, 207].

To describe the dynamics of such a small particle, we prefer to use the following dimensionless numbers:

$$ St = \tau_p / \tau_\eta; \quad \beta \qquad (3.7) $$

with $\tau_\eta = \eta^2/\nu$ the dissipative timescale of the flow structure. Note that in Eq. (3.4) buoyancy has been neglected, which is valid for sufficiently intense turbulence and low gravity, i.e. if $Fr/St \gg 1$ [162], where $Fr = a_\eta/g$ [110, 22].

### Small heavy particles

We consider first the case of small and heavy particles so that in the limit $\beta \to 0$ only drag force is important. In this case, the equation of motion of the particle simply becomes:

$$ \frac{du_p}{dt} = -\frac{u_p - u_f}{\tau_p}. \qquad (3.8) $$

This regime has attracted considerable attention for several decades, most certainly because of its interest in the understanding of many problems both industrial and natural, but probably also because of the simplicity of the model.
Among the remarkable aspects that have been highlighted is the formation of areas of very high particle density (clusters) [12, 45, 246, 69]. Due to their inertia, the particles tend to be ejected from areas of high rotation and accumulate in regions of high shear [163]. This cluster formation is maximum for St numbers around unity. For larger Stokes numbers, particles are expected to respond to structures with longer lifetimes, which are less intense [94].

It is important to note that Eq. (3.8) shows that the inertia of the particle acts as a first-order filter: high-frequency fluctuations in the fluid velocity along the particle trajectory are filtered out.

As shown in Ref. [20] for weakly inertial particles ($St < 1$), the inertia does not have a very large impact on the statistics of the particle acceleration which remain very close to the rate of variation of the fluid velocity along the trajectory. However the acceleration is sufficiently different from the acceleration of a fluid particle that would occupy the same place, to deviate the trajectory of the inertial particle from that of fluid particles and leads to cluster formation.

For larger inertia ($St > 1$) the fluctuations of the fluid along the trajectory of a particle are effectively filtered out and the particle only responds to low frequencies. This results in a much reduced variance of the particle acceleration compared to that of a fluid particle occupying the same position, and the trajectories become much straighter, thus reducing the phenomenon of preferential concentration. Particles with high inertia can pass through the high rotation zones without being significantly affected.

As proposed in Ref. [92] this results in the following relations for the variance of the acceleration of a particle. For $St \ll 1$, the variance of the acceleration remains very close to that of a fluid particle, so neglecting the preferential concentration, and the effect of intermittency (which was discussed at length in the first chapter, but which we leave aside for now) we have:

$$ \langle \mathbf{a}_p^2 \rangle / a_\eta^2 = O(1), \qquad (3.9) $$

where $a_\eta^2 = \langle \epsilon \rangle / \tau_\eta$. For very large inertia, $St \gg \tau_L / \tau_\eta$ where $\tau_L$ is the Lagrangian integral scale of the flow, we can consider that the velocity fluctuations of the particle appear decorrelated from the fluid ones, leading to a “thermalization” of the particle with the flow. And because of the particle inertia the variance of the relative velocity will be equal to that of the fluid: $\langle (\mathbf{u}_p - \mathbf{u}_f)^2 \rangle \approx \langle \mathbf{u}_f^2 \rangle$. This gives the following estimate for the variance of the particle acceleration:

$$ \langle \mathbf{a}_p^2 \rangle / a_\eta^2 \sim Re_\lambda St^{-2}. \qquad (3.10) $$

To obtain this relation we used that $\tau_L \sim \langle \mathbf{u}_f^2 \rangle / \langle \epsilon \rangle$ and that $\tau_L / \tau_\eta \sim Re_\lambda$. For intermediate inertia $1 \ll St \ll \tau_L / \tau_\eta$ one can make the assumption, similar to Kolmogorov, that the variance of the particle acceleration depends only on the average dissipation rate $\langle \epsilon \rangle$ and on the response time of the particle and obtain:

$$ \langle \mathbf{a}_p^2 \rangle / a_\eta^2 \sim St^{-1}. \qquad (3.11) $$

A comparison is made with results computed from the DNS database [21, 141] obtained with the inertial particle model (3.8). We observe that the agreement is only partial and that the three scaling laws seem to give rather the envelope of the evolution of the variance with $St$. This comes from the fact that the Lagrangian inertial zone is very small even at large Reynolds (at least large for the current capabilities of DNS). Indeed from DNS we can see that $\tau_L / \tau_\eta = 0.08 Re_\lambda$ [238, 101]. Thus for $Re_\lambda = 400$ we have only
$\tau_L/\tau_\eta = 32$ which strongly limits the inertial zone where $1 \ll St \ll \tau_L/\tau_\eta$. Thus, without invalidating the proposed scaling laws, it shows that the way of matching these laws has a great importance to obtain an accurate estimate of the variance of the acceleration.

<figure><img src="image_24.png" /><figcaption>Figure 3.1: Variance of the acceleration of an inertial particle normalized by the variance of the acceleration of a fluid particle as a function of the Stokes number, obtained from the DNS database of Refs. [21, 140] at $Re_\lambda = 400$. Comparison with (3.9), (3.10) and (3.11) relations in dotted line, as well as with (3.14) in dashed lines.</figcaption></figure>

A more systematic approach taking advantage of the linearity of (3.8) proposes to calculate the associated transfer function [253, 104, 273, 279]. By taking the Fourier transform of (3.8) we obtain a relation between the spectrum of the particle velocity and that of the fluid velocity along the trajectory:

$$E_p(\omega) = \frac{1}{1 + \omega^2 \tau_p^2} E_f(\omega) \qquad (3.12)$$

with $E_p(\omega) = \hat{u}_p \hat{u}_p^*$ and $E_f(\omega) = \hat{u}_f \hat{u}_f^*$ the spectral density of velocity of the particle and of the fluid. The evolution of the spectrum of the fluid velocity along the trajectory of the particle with $St$ from DNS is shown in figure 3.2. Without going into details, we can see that, on the one hand, for very weakly inertial particles, the spectrum of the velocity seen by the particle is relatively well described by the Hinze spectrum [256, 184]:

$$E_f(\omega) \approx \frac{k_0 \langle \varepsilon \rangle \tau_L^2}{1 + \tau_L^2 \omega^2}, \qquad (3.13)$$

at least for $\omega \ll \omega_\eta = 2\pi/\tau_\eta$. Note that this spectrum simply corresponds to an exponential decay of the Lagrangian autocorrelation of the velocity.

Moreover, we note that at frequencies lower than $\omega < 2\pi/\tau_p$ we can consider that the spectra of the fluid velocity are invariant with $St$. Specifically, the value of $k_0$ seems constant.
Using the normalization condition $\langle u_f^2 \rangle = 2 \int E_f(\omega) d\omega$ we find that $k_0 = \frac{1}{\pi} \frac{\langle u_f^2 \rangle / \langle \epsilon \rangle}{\tau_L}$ where $\langle u_f^2 \rangle$ represents the variance of the velocity of the particle along the trajectory. We measured from DNS $k_0 = 1.04$. The dependence on $St$ at higher frequency reflects, on the one hand, the effect of the preferential concentration (moderate $St$) and, on the other hand, the fact that the trajectory becomes more ballistic and thus cuts the turbulent structures developing a more energetic spectrum at high frequency in agreement with the sweeping mechanism proposed by Ref. [255].

<figure><img src="image_25.png" /><figcaption>Figure 3.2: Spectrum of the fluid velocity along the trajectory of an inertial particle for $St = 0.24, 0.9, 1.5, 3., 4.5, 7. 5, 15, 30, 45, 60, 75$ (from black to orange), obtained from the DNS database of Refs. [21, 140] at $Re_\lambda = 420$. Comparison with the relation $\langle \epsilon \rangle \omega^{-2}$ and the relation (3.13).</figcaption></figure>

The spectrum of the particle acceleration is simply obtained by multiplying that of its velocity by $\omega^2$. We obtain an estimate of the variance of the acceleration by integration:

$$ \langle \boldsymbol{a}_p^2 \rangle = 2 \int_0^\infty \omega^2 E_p(\omega) d\omega \\ \approx \frac{c_0 a_\eta^2}{1 - St_L^2} \left[ \frac{\tan^{-1}(c_1 St)}{c_1 St} - \frac{\tan^{-1}(c_1 \tau_L / \tau_\eta)}{c_1 \tau_L / \tau_\eta} \right] . \quad (3.14) $$

The constant $c_1 = 2.8$ appears, because the previous integral is truncated for $\omega > c_1/\tau_\eta$ and we introduce $c_0 = 2c_1k_0$ to simplify the notations. We also introduced $St_L = \tau_p/\tau_L = St/0.08Re_\lambda$. The relation (3.14) is compared with the DNS in the figure 3.1. In this figure the variance of the acceleration is normalized by the variance of the acceleration of a fluid particle, which is estimated from (3.14) with $St = 0$. We can see a fairly good agreement with the DNS, although around $St \approx 1$ the variance of the acceleration is underestimated, because the effect of the preferential concentration on the spectrum of the fluid velocity has been neglected. One can check that the expression (3.14) is in agreement with the
relations (3.9), (3.10) and (3.11). For example, in the case of large Reynolds numbers, $\tau_L/\tau_\eta = Re_\lambda/0.08 \gg 1$ and $\tau_L/\tau_\eta \gg St$, the relation (3.14) simplifies into

$$ \langle \mathbf{a}_p^2 \rangle \approx c_0 a_\eta^2 \frac{\tan^{-1}(c_1 St)}{c_1 St} . \qquad (3.15) $$

This relation behaves indeed as $St^{-1}$ for $St \gg 1$ and tends to a constant for $St \ll 1$. In the limit where the Stokes number is large $St \gg 1$ (and also $\tau_L/\tau_\eta \gg 1$) we find:

$$ \langle \mathbf{a}_p^2 \rangle \approx a_\eta^2 \frac{k_0 \pi / 2}{St(1 + St_L)} . \qquad (3.16) $$

This corresponds to (3.10) and (3.11). We can notice that in the previous relation the constant $c_1$ is eliminated meaning that for large Stokes numbers the dissipative zone has no more direct effect, because it appears anyway filtered by the inertia of the particle.

Moreover, it should be noted that $k_0$ presents dependence with the Reynolds number in order to account for the intermittency of the flow (as discussed in the first chapter, see the relation (2.43)) reflecting that for very high Reynolds numbers even at high Stokes numbers, one should expect that the large fluctuations at small scales might influence the dynamics of the particle.

## Small, heavy or light particles

We are now interested in particles with arbitrary density, small or large, but whose size remains small ($d/\eta \ll 1$ and $Re_p \ll 1$). We mention that bubbles can be assimilated to very light particles $^4$. It is therefore a priori necessary to take into account the effect of the inertia of the fluid, which includes the added-mass force and the Tchen force. The dynamics of the particle is therefore given by (3.4). We note

$$ F_I = \beta D_t \mathbf{u}_f, \qquad (3.17) $$

the inertia force per unit of displaced mass, and

$$ F_D = -(\mathbf{u}_p - \mathbf{u}_f)/\tau_p . \qquad (3.18) $$

The drag force also per unit of displaced mass, with the relaxation time always defined by (3.6).

The first question concerns the importance of the inertia force of the fluid compared to the drag force. In the case of very small inertia $St \to 0$, the particle behaves like a tracer. Therefore in this limit, its acceleration will be identical to the acceleration of a fluid particle $\mathbf{a}_p = \mathbf{a}_f = D_t \mathbf{u}_f$. Then from (3.4), and using (3.17), we have $F_D = (1 - \beta)\mathbf{a}_f$. It implies that for $\beta > 0$

$$ F_D = \frac{1 - \beta}{\beta} F_I . \qquad (3.19) $$

For example, for very light particles ($\beta = 3$), we expect that the two forces are of the same order of magnitude, but of opposite directions.

<footnote>$^4$There are however some precautions to take, such as contamination of the interface by surfactants which change the dynamics condition on the interface (free slip or no-slip), the deformation of the bubble, characterized for example by the Weber number ...</footnote>
As the inertia of the particle increases, the drag force (per unit mass) is expected to decrease, as seen earlier. This suggests that for inertial particles with $St \gg 1$, the fluid inertia force can be important even if the density ratio is not very low. Moreover, for light particles we expect the variance of the particle acceleration to increase with $St$, since $F_I$ and $F_D$ are opposite.

These behaviors are confirmed in figure 3.3 showing the evolution of the variance of the particle acceleration for different values of $\beta$. The reported values were obtained by Zhen-tong Zhang from DNS [279] for $\beta = 3$, for $\beta = 2.5$ and 0.14 by those of Ref. [46], and for $\beta = 0$ by Ref. [21]. We see that for $\beta > 0$, the variance of the acceleration tends to a constant value when $St$ becomes large, and that the asymptotic value depends on $\beta$. This reflects the fact that the drag force becomes negligible at large $St$. For particles with density less than the fluid ($\beta > 1$) we also see that the variance becomes larger than for a fluid particle. In fact, we will show below that this asymptotic value is $\beta^2$.

A remark is necessary: the model (3.4) is based on the hypothesis that $d/\eta < 1$ and $Re_p < 1$ but this imposes limits for the value of $St$ and $\beta$ as we can see with the relation:

$$ \left(\frac{d}{\eta}\right)^2 = \frac{18\beta St}{(1 + C_M)}. \qquad (3.20) $$

For $\beta = 3$, the Stokes number must be less than 5 to guarantee that the diameter remains small enough for neglecting the finite size of the particles (corresponding roughly to $d/\eta < 10$). Indeed on the figure 3.3, we see that for $\beta$ close to 3, the experimental results show a decrease of the variance when $St$ increases, this is attributed to finite size effects of the particle. These effects are discussed in paragraph 3.2.2. We come back to the condition $Re_p \ll 1$ a little later.

The approach of Tchen can be generalized [104, 166, 3, 279] to account for the fluid inertia force in order to propose an estimate of the variance of the acceleration of a small, heavy or light particle as a function of $St$, $\beta$, and $\tau_L/\tau_\eta \approx 0.08 Re_\lambda$:

$$ \frac{\langle \mathbf{a}_p^2 \rangle}{a_\eta^2} \approx c_0 \left[ \beta^2 + \frac{1 - \beta^2}{1 - St_L^2} \frac{\tan^{-1}(c_1 St)}{c_1 St} - \frac{1 - \beta^2 St_L^2}{1 - St_L^2} \frac{\tan^{-1}(c_1 \tau_L / \tau_\eta)}{c_1 \tau_L / \tau_\eta} \right] = \Gamma_a(St, \beta, \tau_L / \tau_\eta). \qquad (3.21) $$

Again with $St_L = \tau_p/\tau_L = St\tau_\eta/\tau_L$ and $c_1 = 2.8$. This estimate is compared to the results of DNS in figure 3.3. We see that the agreement is relatively good, at least as long as the finite size effects are negligible. We can notice that for $\beta = 0$, the relation (3.21) is indeed identical to (3.14). For $Re_\lambda \gg 1$ and $St_L \ll 1$ we can simplify the relation:

$$ \frac{\langle \mathbf{a}_p^2 \rangle}{a_\eta^2} \approx c_0 \left[ \beta^2 + (1 - \beta^2) \frac{\tan^{-1}(c_1 St)}{c_1 St} \right]. \qquad (3.22) $$

We notice that for $St \ll 1$ and $Re_\lambda \gg 1$ the variance of the acceleration tends to a constant value $\langle \mathbf{a}_p^2 \rangle \approx c_0 a_\eta^2 \beta^2$.

The equation (3.21) is obtained from the following transfer function

$$ E_p(\omega) = H_u^2(\omega) E_f(\omega) \quad ; \quad H_u^2(\omega) = \frac{1 + \beta^2 \omega^2 \tau_p^2}{1 + \omega^2 \tau_p^2}, \qquad (3.23) $$

and by estimating, as before, that the spectrum of the fluid velocity at the particle position is given by (3.13). This transfer function is obtained at the cost of an approximation on
<figure><img src="image_26.png" /><figcaption>Figure 3.3: Variance of acceleration normalized by the acceleration of a fluid particle from the DNS of [279] for a bubble $\beta = 3$, $Re_\lambda = 215$ (• and continuous line), the DNS of [21, 140] for dense particles $\beta = 0$, $Re_\lambda = 400$ (▲ and continuous line), the DNS of [46] for $\beta = 2.5$ and $Re_\lambda = 180$, with or without finite size correction (+ with dash-dotted and dotted lines respectively), and $\beta = 0.14$ (+ and continuous line), experimental measurements of [205] for $Re_\lambda = 145 - 230$ and $\beta \approx 3$ (× and dashed lines). Comparison with the equation (3.21): $\langle \alpha_b^2 \rangle_{(St,\beta,\tau_L/\tau_\eta)} / \langle \alpha_b^2 \rangle_{(0,1,\tau_L/\tau_\eta)}$ with $c_1 = 2.8$ and $Re_\lambda = 200$ for $\beta = 0, 0.14, 0.5, 1, 1.5, 2, 2.5$ and 3 from black to orange, in dash line, and with $c_1 = 2.8$ and $Re_\lambda = 400$, $\beta = 0$ in black dotted line. Inset: the transfer function (3.23) as a function of $\omega\tau_b$ for $\beta = 0, 0.5, 1, 2$ and 3 from black to orange.</figcaption></figure>

the total acceleration of the fluid at the particle position. We have indeed replaced the total acceleration of the fluid at the position of the particle $Du_f/Dt$ by the Lagrangian derivative along the trajectory of the particle $du_f/dt = \partial_t u_f + u_p \cdot \nabla u_f$. This allows the use of the linear response approach, since it is equivalent to considering that the velocity of the particle responds only to the fluctuations of the fluid velocity along the trajectory and its time derivative. This assumption can be evaluated by writing the material derivative of the fluid in the following form

$$D_t u_f = \partial_t u_f + u_f \nabla u_f = d_t u_f - (u_p - u_f) \cdot \nabla u_f. \quad (3.24)$$

We thus see that our approximation is a priori valid for Stokes numbers small enough for the relative velocity to be small (as discussed in the previous paragraph about the drag force). Figure 3.4 compares the variance of $D_t u_f$ and $d_t u_f$, as well as the variance of their difference, for the extreme cases $\beta = 0$ and $\beta = 3$. We notice that $\langle (D_t u_f)^2 \rangle$ evolves only slightly with $St$, although it presents a slightly lower value around $St = 1$ for both $\beta = 0$ and $\beta = 3$. On the other hand, the Lagrangian derivative of the velocity along the trajectory shows a clear increase with $St$. However, it should be noted that our estimate is not so bad, because we have also neglected the phenomenon of preferential concentration
(in fact the effect of $St$ in the velocity spectrum seen by the particle $E_f$), and thus the variance of the inertial term predicted by the approach does not depend on $St$ which is finally in relatively good agreement with the DNS. This is the result of two errors which, by chance, compensate each other.

<figure><img src="image_27.png" /><figcaption>Figure 3.4: Evolution with $St$ of the variance of the material derivative of fluid velocity at the particle position $D_t u_f$ (black) and of its Lagrangian derivative along the particle trajectory $d_t u_f$ (blue/red), for inertial particles $\beta = 0$ at $Re_\lambda = 400$ (left) and light particles $\beta = 3$ at $Re_\lambda = 200$ (right). The quantities are normalized by the variance of the acceleration of a fluid particle computed on the whole domain. Inset: variance of $D_t u_f - d_t u_f$.</figcaption></figure>

We just mention that in Ref. [279] we have also studied by DNS, the PDFs of the particle acceleration for $St = 3$. We observe that for $St = O(1)$, the flatness of the PDFs is much larger than that for fluid particles. This is because the two forces $F_D$ and $F_I$ are of the same order of magnitude and mostly anti-aligned, as the drag force responds quickly to the solicitation imposed by the inertial force, leading to a relatively low acceleration. However, when the two forces are aligned, this leads to very high accelerations. Using the response function approach, we can estimate the variance of the drag force $F_D$ as well as of the fluid inertia force $F_I$:

$$ \frac{\langle F_D^2 \rangle}{a_\eta^2} \approx c_0 \frac{(1-\beta)^2}{1-St_L^2} \left[ \frac{\tan^{-1}(c_1 St)}{c_1 St} - \frac{\tan^{-1}(c_1 \tau_L / \tau_\eta)}{c_1 \tau_L / \tau_\eta} \right] = \Gamma_D(St, \beta, \tau_L / \tau_\eta), \quad (3.25) $$

$$ \frac{\langle F_I^2 \rangle}{a_\eta^2} \approx c_0 \beta^2 \left[ 1 - \frac{\tan^{-1}(c_1 \tau_L / \tau_\eta)}{c_1 \tau_L / \tau_\eta} \right] = \Gamma_I(St, \beta, \tau_L / \tau_\eta). \quad (3.26) $$

We see that indeed the estimate of the predicted inertia force is independent of $St$, as mentioned above.

The previous relations were obtained by assuming a Stokes regime for the flow around the particle. This results in a Reynolds number of the particles that must necessarily be very small in front of 1. Using the relation of the variance of the drag force (3.25), one can estimate the evolution of the Reynolds number or more exactly of its second order moment with $St$:
$$ \langle Re_p^2 \rangle = \langle \boldsymbol{F}_d^2 \rangle \tau_p^2 d^2 / \nu^2 = \frac{18\beta}{1+C_M} St^3 \Gamma_D (St, \beta, \tau_L / \tau_\eta) . \quad (3.27) $$

## Effect of the finite Reynolds number of the particles

The relative velocity between the particle and the fluid increases with the inertia of the particle. The inertial effects of the flow around the particle are characterized by the Reynolds number of the particle $Re_p$. These are generally considered to be present from $Re_p = 0(1)$, although the drag force remains dominated by viscous effects until $Re_p = 0(100)$. Such inertial effects are not taken into account in (3.4) which assumes a Stokes regime for the flow around the particle giving a drag force proportional to the relative velocity. One can introduce a correction $\phi(Re)$ to the drag force to take into account finite Reynolds number effects:

$$ \boldsymbol{F}_D = -\phi(Re_p) \frac{\boldsymbol{u}_p - \boldsymbol{u}_f}{\tau_p} . \quad (3.28) $$

We choose<fnref n="5" /> by definition that for a solid particle with $Re_p \ll 1$ we have

$$ \phi(Re) = 1. \quad (3.29) $$

For a solid particle, we obtain the correction from the expression of the drag coefficient proposed by Ref. [242], valid for $Re_p < 800$:

$$ \phi(Re_p) = 1 + 0.15 Re_p^{0.687} . \quad (3.30) $$

This empirical relation is based on experimental data for the mean drag and tends towards the Stokes solution in the limit $Re_p \to 0$. The $\phi$ correction can also integrate other effects such as the mobility of the interface or its contamination, the viscosity of the fluid in the case of drops or even the shape of the particle. For example, we can also consider the case of clean spherical bubbles (spheres with a free shear condition at the interface), with the Mei drag law valid for all $Re_p$:

$$ \phi(Re_p) = \frac{2}{3} + \frac{2}{3} \left( \frac{8}{Re_p} + \frac{1}{2} \left( 1 + \frac{3.315}{Re_p^{1/2}} \right) \right)^{-1} . \quad (3.31) $$

This relation is based on DNS and has been constructed so as to find in the limit $Re_p \to 0$ the viscous solution and the Oseen development [102, 227, 252]. At high Reynolds numbers the relation (3.31) gives $\phi(Re_p) = 2$ and $\phi(Re_p) = 2(1 - 2.211Re_p^{-1/2})$ corresponding respectively to the Levich viscous potential solution [147] and the Moore correction [179] taking into account the development of a boundary layer. The $\phi$ correction coefficient is related to the drag coefficient by

$$ \phi(Re) = \frac{C_D}{C_{D,0}}, \quad (3.32) $$

with $C_{D,0} = 48/Re_p$ the drag coefficient of a solid particle in the Stokes regime as a reference. The evolution of $\phi$ with $Re_p$ corresponding to (3.29), (3.30) and (3.31) is shown

<footnote n="5">⁵This choice is different in [280] where we had chosen as reference the case of a clean spherical bubble, i.e. with a free slip condition at the interface.</footnote>
<figure><img src="image_28.png" /><figcaption>Figure 3.5: On the left: $\phi(Re_p)$ as a function of $Re_p$. The case of a solid particle: (gray dotted lines) Stokes solution $\phi = 1$ (relation (3.29)); (thin gray dotted lines) the Oseen solution $\phi = 1 + 3Re/8$; (black dotted lines) the Schiller and Naumann relation (3.30). For a spherical bubble: (black line) relation (3.31); (black dotted lines) Stokes regime $\phi = 2/3$; (thin gray dotted lines) the Taylor and Acrivos solution $\phi = 2/3(1 + Re/8)$; (thin gray dotted lines) the Moore relation $\phi(Re_p) = 2[1 - 2.211Re_p^{-1/2}]$; (gray dotted lines) Levich's solution $\phi(Re_p) = 2$.</figcaption></figure>

Right: Evolution of $\langle Re_p \rangle \phi'((\langle Re_p \rangle))/\phi((\langle Re_p \rangle))$ as a function of $\langle Re_p \rangle$. The symbols correspond to the results of the DNS: (○) for the law (3.31) and (×) for (3.30). The lines are the analytical results.

in Figure 3.5. The expressions for $\phi$ are based on the implicit assumption, generally accepted, that the expression for the mean drag force is applicable to the instantaneous drag force and that the decomposition of the force balance on the particle remains unchanged. Such an assumption is certainly valid if the fluctuations of the Reynolds number are small compared to its mean value and the transient effects (related to the history force) are very brief.

In Ref. [280] we considered simulations for the three drag laws (3.29), (3.30) and (3.31), with $\beta = 3$. We found that the instantaneous Reynolds number of the particles in these three cases presents a distribution close to a log-normal with a mean of the same order of magnitude as its standard deviation:

$$ \langle Re_p \rangle \approx \sqrt{\langle Re_p^2 \rangle / 2}. \qquad (3.33) $$

To propose an estimate of the variance of the acceleration and forces, taking into account the inertial effects in the response of the particle, we can decompose the Reynolds number of the particle into an average part and a fluctuating part $Re_p = \langle Re_p \rangle + Re'_p$ and approach the drag force (3.28) thanks to a Taylor development of the function $\phi$ around $\phi(\langle Re_p \rangle)$:

$$ F_D = -\phi(\langle Re_p \rangle) \frac{u_p - u_f}{\tau_p} - Re'_p \phi'((\langle Re_p \rangle)) \frac{u_p - u_f}{\tau_p} + \dots \qquad (3.34) $$

with $\phi'$ the derivative of $\phi$ with respect to $Re_p$. Although the Reynolds number of the particle presents important fluctuations, we show in figure 3.5 by plotting the ratio $\langle Re_p \rangle \phi'((\langle Re_p \rangle))/\phi((\langle Re_p \rangle))$, that we can retain only the first term of the Taylor expansion,
at least as long as we are interested in basic statistics (typically mean or variance). When the Stokes number becomes very large, the second order is no longer negligible, but the first term remains dominant.

The relation (3.34) allows introducing an effective relaxation time $\tau_p^*$, as proposed for example in Refs. [77, 26, 280], which allows to take into account the effect of the finite Reynolds number of the particles:

$$ \tau_p^* = \tau_p / \phi(\langle Re_p \rangle), \qquad (3.35) $$

and an effective Stokes number

$$ St_* = St / \phi(\langle Re_p \rangle). \qquad (3.36) $$

To illustrate the relevance of introducing $St_*$, we present on figure 3.6 the evolution of the variance of the drag force and of the acceleration of the particles as a function of $St$ and $St_*$ for the DNS performed with the drag laws (3.29), (3.30) and (3.31) and $\beta = 3$. It can be seen that the evolution with the Stokes number differs quite markedly when the drag law is modified. When these evolutions are plotted as a function of $St_*$ we notice that the different curves almost merge. This means that the effects of finite $Re_p$ can be taken into account, at first order, thanks to the effective Stokes number.

Moreover, we have shown in the previous paragraph that for $\phi(Re_p) = 1$ we can estimate the evolution of the variance of the forces and the acceleration as a function of $St$ and $\beta$ thanks to the relations (3.21), (3.25), (3.26). Thus, the (approximate) self-similarity of the variances as a function of $St_*$ in figure (3.6) indicates that we can estimate the variance of the particle forces and acceleration taking into account inertial effects ($Re_p > 1$) with the formulas (3.21), (3.25), (3.26) evaluated with $St_*$ instead of $St$:

$$ \frac{\langle a_p^2 \rangle}{a_\eta^2} \approx \Gamma_a(St_*, \beta, \tau_L/\tau_\eta) \qquad (3.37) $$

$$ \frac{\langle F_D^2 \rangle}{a_\eta^2} \approx \Gamma_D(St_*, \beta, \tau_L/\tau_\eta) \qquad (3.38) $$

$$ \frac{\langle F_I^2 \rangle}{a_\eta^2} \approx \Gamma_I(St_*, \beta, \tau_L/\tau_\eta) \qquad (3.39) $$

To use these relations one needs to estimate $St_*$, which requires to know the mean value of the particle Reynolds number. By using (3.33), (3.34), (3.38) and the definition of $Re_p$, we can obtain an implicit relation allowing estimating, to the first order, the average Reynolds number

$$ \langle Re_p \rangle \approx \frac{St}{\phi(\langle Re_p \rangle) \eta} \left( \frac{1}{2} \Gamma_D(St/\phi(\langle Re_p \rangle), \beta, \tau_L/\tau_\eta) \right)^{1/2} \qquad (3.40) $$

This relation can be solved iteratively by taking for example as initial value $\langle Re_p \rangle = 0$. A validation of (3.40) can be found in Ref. [280] based on comparisons with the DNS for different drag laws.

Fluid inertia forces have a dominant role in the dynamics of light particles ($\rho_p/\rho_f \ll 1$), in contrast, it is generally considered that for very dense particles ($\rho_p/\rho_f \gg 1$), the fluid inertial force is irrelevant. This suggests that one can conclude about the relevance of the
<figure><img src="image_29.png" /><figcaption>Figure 3.6: Left: Variance of the drag force normalized by the variance of the particle acceleration. Right: Variance of the particle acceleration normalized by the variance of the fluid particles. Top: Evolution as a function of $St$; bottom: Evolution as a function of $St_* = St/\phi(\langle Re_p \rangle)$. With $\phi(Re_p)$ given by (3.29) (○), (3.30) (*) and (3.31) (Δ). The dotted lines correspond to estimates by (3.38) and (3.37).</figcaption></figure>
inertial force based only on the density. Using the estimates (3.38) and (3.39) for the drag force and inertia force, we show in what follows that this is more subtle and also depends on the particle size.

Figure 3.7 shows the evolution of the ratio $\langle F_I^2 \rangle / \langle F_D^2 \rangle$ estimated from (3.38) and (3.39) for a range of density and particle size. For this figure we chose the $Re_p$-correction (3.30) corresponding to the drag force of a solid sphere given by the Schiller and Naumann relation. The figure 3.7(a) presents 3 levels of the force ratio $\langle F_I^2 \rangle / \langle F_D^2 \rangle = 1\%, 10\%$ and 100% as a function of $\rho_p/\rho_f$ and $St$. The shaded area indicates the region of parameter space in which the fluid inertia force is not negligible (arbitrarily set to $\langle F_I^2 \rangle / \langle F_D^2 \rangle = 10\%$). As expected, we observe that for light particles ($\rho_p/\rho_f < 1$), the inertia force of the fluid is dominant, whereas for very heavy particles, $\rho_p/\rho_f \gg 1$, we observe that the inertia force of the fluid can be neglected. However, we observe that when the Stokes number of the particle increases, the density threshold from which the inertia force is negligible increases. Typically for $St < 1$, this force is negligible if $\rho_p/\rho_f > 10$, while for $St = 100$ it requires $\rho_p/\rho_f > 100$. The non-negligible role of the fluid inertia force for very heavy particles and of inertia (or size, both are related, see (3.20)) is simply explained by the observation that at first order the fluid inertia force (per unit of displaced mass) is independent of particle size, as long as finite size effects can be ignored, while the magnitude of the drag force (again per unit of displaced mass) decreases with particle size, as shown in figures 3.1 and 3.6.

Moreover, considering the evolution of the force ratio as a function of $\rho_p/\rho_f$ and $d/\eta$ given in figure 3.7(b), we can notice that for particles of size $d/\eta \approx 3$ the inertial force remains important even for very large density ratios. This observation depends on the Reynolds number of the flow. Indeed, for $Re_\lambda = 400$, the inertia force of the fluid should not be neglected for particles larger than $d/\eta \approx 7$. Since it has been proposed by Ref. [46] that the finite size effect can be ignored for particles smaller than $d/\eta \approx 10$, the results presented in this section indicate that the added-mass force may turn out to be of the order of magnitude of the drag force, when the inertia of the particles becomes large, even for high density ratio. For higher Reynolds numbers, one has to keep in mind that intermittency is not taken into account, and since the fluctuations of the fluid acceleration can be much larger than its standard deviation, this contributes a priori to reinforce the effect of the fluid inertia force.
<figure><img src="image_30.png" /><figcaption>Figure 3.7: Contours level $\langle F_I^2 \rangle / \langle F_D^2 \rangle = 0.01, 0.1$ and 1, estimated from (3.38) and (3.39) for $Re_\lambda = 100$ (continuous black lines) and $Re_\lambda = 400$ (dashed lines) as a function of $\rho_p/\rho_f$ and $St$ (a) and of $\rho_p/\rho_f$ and $d/\eta$ (b). The red lines give the iso values of the Reynolds number obtained from (3.40) for $Re_\lambda = 100$. The green lines give the iso values of $d/\eta$ (a) and $St$ (b). The shaded region corresponds to $\langle F_I^2 \rangle / \langle F_D^2 \rangle > 0.1$.</figcaption></figure>

### 3.2.2 Large particles

We are considering the statistics of the acceleration of finite size particles, i.e. whose diameter is not negligible in comparison with the Kolmogorov scale. Experimental works of Refs. [264, 210, 209, 263, 38, 40, 267] have shown that for particles with $d > \eta$ and having a density ratio close to 1 the variance of their acceleration varies as $\langle a_p^2 \rangle \sim \langle \epsilon^{4/3} \rangle d^{-2/3}$. The experiments of Refs. [209, 38] have studied in detail the effect of size and density for size ratios between $d/\eta = 10$ and 26 and density between 1 and 100. They showed a complicated evolution of the variance of the acceleration with size and density. It also appears that the tails of the PDFs of particle acceleration tend to shrink with increasing size [263]. It should be noted, however, that this behavior is not observed in the results of [209] which has PDFs that appear to be nearly invariant with size, this may be due to the resolution limitations of the measurements.

On the other hand, in numerical simulations of Refs. [105, 57], at moderate Reynolds numbers, it was observed that the variance varies as $d^{-4/3}$. This would result from the fact that the particle dynamics is dominated by large flow structures leading to $\langle a_p^2 \rangle \sim \epsilon^{2/3} d^{-4/3} \langle (u_f - u_p)^2 \rangle$. It should be noted that studies on finite size particles are still mainly experimental and that DNS at high Reynolds numbers are not yet available.

To model these different aspects, we can assume first that the decomposition of the hydrodynamic force (drag, mass added inertia force ...) remains relevant for particles much larger than the Kolmogorov scale. This remains, of course, hypothetical, because there is no rigorous derivation of the hydrodynamic force for such particles (except for the added-mass force). Many works are currently looking for empirical expressions of such forces [150]. The main difficulty in this path being to define the fluid velocity seen by the particle that appears in the expression of the forces. We return to this point in the paragraph 3.4. In this hypothesis, we consider that the hydrodynamic force on the particle can generally be given by the drag force and the inertia force, as in the previous part.

We discuss here, how the expressions of these two forces can be adapted to take into
account the finite size effect.

To determine the drag force it is necessary to take into account that the flow is non-homogeneous at the scale of the particle, thus leading to exchanges of momentum between the fluid and the particle, and which can be related to an additional drag. In Ref. [92] it has been proposed to model this aspect by introducing an effective viscosity to account for inertial effects in the drag calculation. The effective viscosity takes into account the momentum flux induced by fluctuations at the particle scale. Using the Prandtl mixing length model, it is estimated as:

$$ \nu_{eff} = \nu + u'd = \nu + c\epsilon_d^{1/3}d^{4/3} \qquad (3.41) $$

where $u'$ is the characteristic velocity at the scale of the particle that we estimated from the Kolmogorov scaling law. Thus, $\epsilon_d$ must be understood as the dissipation rate averaged at the particle scale. Putting aside for the moment the aspects related to $\nu_{eff}$ fluctuations and intermittency, we write that $\nu_{eff}/\nu = 1+(d/d_c)^{4/3}$ where we have introduced $d_c \approx 10\eta$ the diameter of the particle from which the effects of finite size are felt. From this turbulent viscosity, we can construct a turbulent response time for the particle:

$$ \tau_p^* \rightarrow \tau_{p,eff} = \tau_p^* \frac{\nu}{\nu_{eff}}. \qquad (3.42) $$

We have used here $\tau_p^* = \tau_p/\phi(\langle Re_p \rangle)$ introduced in the previous paragraph (3.35) in order to take into account also the finite Reynolds number effects (which are also inertial effects). The idea is then to estimate the drag force, per unit mass, from the velocity field filtered at the scale of the particle $\eta \ll d \ll L$, as $(\mathbf{u}_f - \mathbf{u}_p)/\tau_{p,eff}$.

The inertia forces of the fluid, which is directly proportional to the acceleration of the fluid, can be seen as a multiplicative cascade (as discussed in Part 1). Thus this term will be directly influenced by the filtering effect operated by the particle size. This is equivalent to considering that from the point of view of the particle the turbulent cascade stops at a scale $d > \eta$. Using a dimensional argument, we can write that the acceleration of the fluid at scale $d$ thus varies as $a_d \sim \epsilon_d^{2/3}d^{-1/3}$. This view is consistent with the Faxen formulation used in [45] also see [135, 254, 148]. Thus the proposal is to replace the material derivative of the fluid entering the expression of the fluid inertia term by a filtered acceleration:

$$ D_t \mathbf{u}_f \rightarrow \overline{D_t \mathbf{u}_f} \sim D_t \mathbf{u}_f (d/d_c)^{-1/3}. \qquad (3.43) $$

In this relation we used again $d_c \approx 10\eta$ as introduced previously.

With the two previous ingredients, we propose to model the motion of large particles by the following equation:

$$ d_t \mathbf{u}_p = - \frac{\mathbf{u}_p - \mathbf{u}_f}{\tau_{p,eff}} + \beta_{eff} D_t \mathbf{u}_f \qquad (3.44) $$

with $\beta_{eff} = \beta(d/d_c)^{-1/3}$.

This amounts to introducing two dimensionless numbers $St_{eff} = St_* \frac{\nu}{\nu_{eff}}$ and $\beta_{eff}$ instead of $St_*$ and $\beta$ which are the relevant dimensionless numbers for the case of small particles. One can notice that for $d/\eta \gg 1$, $St_{eff} \sim (d/\eta)^{2/3}$. Thus when the inertia force of the fluid is negligible we expect to find (from what we saw with the equation (3.11)) that the variance of the acceleration behaves as $\langle a_p^2 \rangle \approx a_\eta^2/St_{eff}^{-1} \sim a_\eta^2(d/\eta)^{-2/3}$. On the other hand,
when the inertia force dominates we saw with (3.22) that the variance of the acceleration was given by $\langle \mathbf{a}_p^2 \rangle \approx a_\eta^2 \beta_{eff}^2$, which again gives the scaling law in $(d/\eta)^{-2/3}$.

In fact, one can argue that for large particles $\eta \ll d \ll L$ the various forces are dominated by inertia and thus must exhibit a scaling law as $\epsilon_d^{2/3} d^{-1/3}$. This can be formalized by considering the momentum exchanged between the particle and the fluid per unit time $dP/dt$. Placing ourselves in the reference frame of the particle, we write that $dP$ is the mass of fluid swept by the particle $(\rho_f u' \pi d^2 dt/4)$ multiplied by $u'$. Thus, with $u' \sim (\epsilon_d d)^{1/3}$, we find that the acceleration of the particle, i.e. the momentum exchanged per unit of time and displaced mass $((\rho_p + C_M \rho_f)\pi d^3/6)$ varies as:

$$a_p^2 \sim (\rho_p / \rho_f + C_M)^{-2} \epsilon^{4/3} d^{-2/3}. \qquad (3.45)$$

Finally, assuming that the relation (3.21) giving the variance as a function of $St$ and $\beta$ proposed previously remains valid, providing we substitute $\beta$ by $\beta_{eff}$ and $St$ by $St_{eff}$, one can thus obtain an estimate of the variance of the acceleration of large particles, given here in the limit of large Reynolds numbers:

$$\begin{aligned} \frac{\langle \mathbf{a}_p^2 \rangle}{a_\eta^2} &= \Gamma_a(St_{eff}, \beta_{eff}, \tau_L/\tau_d) \\ &\approx \beta_{eff}^2 + (1 - \beta_{eff}^2) \frac{\tan^{-1}(St_{eff})}{St_{eff}}. \end{aligned} \qquad (3.46)$$

The evolution of the estimate for the variance of the acceleration as a function of $d/\eta$ and $\rho_p/\rho_f$ is shown in figure 3.8. This estimate is compared with the experimental results of Refs. [205, 263, 209]. It can be seen that the agreement is reasonable over the range of diameters and density considered experimentally. For $d/\eta \gg 1$, we find that (3.46) behaves like the asymptotic expression (3.45). The expression (3.46) is a priori valid only for large enough $d/\eta$ as $\beta_{eff}$ diverges for $d/\eta \to 0$. However, we can note that (3.46) behaves as (3.22) in the limit $d/\eta \to 0$. This can be shown with the series expansion of $\tan^{-1}(x)/x$ around 0. This can also be observed in figure 3.8. Nevertheless, the behavior in the transition zone between small particles ($d < d_c$) and large ones ($d > d_c$) is undoubtedly much more complex as can be seen in figure 5 of Ref. [209]. It is quite likely that in addition to preferential concentration effects, there is also the effect of lift forces.
<figure><img src="image_31.png" /></figure>

Figure 3.8: Variance of the acceleration of finite-size particles. Comparison with the experiments of Refs. [205, 263, 209] and the model (3.46) for $ρ_p/ρ_f = 0, 0.1, 1, 2, 5, 10, 20, 50$ and $100$ (from orange to black). The gray dashed lines correspond to the power law $d/η^{-2/3}$, and the gray dotted line to the power law $d/η^{-0.81}$.

## 3.3 Stochastic modeling of the dynamics of a particle and application to LES

In this section, we are interested in taking into account the fluctuations at subgrid scale in the dynamics of particles.

### 3.3.1 The large eddy simulation approach

To compute a turbulent flow by the direct numerical simulation (DNS) method from the Naviers-Stokes equations, all scales of the flow must be solved. The smallest scale of the turbulence, the Kolmogorov scale $\eta$, imposes the resolution of the mesh, while the size of the integration domain must include the largest structures. It is then usually considered that the necessary number of meshes points $N$ varies as $N \sim (L/\eta)^3 \sim Re^{9/4}$. The strong dependence of the number of operations to be performed with the Reynolds number, makes that the computational power remains the limiting factor to simulate flows at high Reynolds numbers, despite an exponential increase in the computational power. This is illustrated in figure 3.9 showing the time evolution of the largest number of mesh points and the maximum Reynolds number, based on the Taylor scale, $Re_\lambda$ for DNS of incompressible and isotropic turbulent flow. We can see in this figure the counterpart of Moore's law for fluid mechanics: since the 70s, the number of mesh points has doubled every 2 years and the maximum Reynolds number $Re_\lambda$ has doubled about every 8 years$^6$. Nevertheless, we can note that lately, although the mesh number follows an exponential growth, this does not seem to be the case for the Reynolds number reached by these simulations. Indeed Refs. [68, 272] draw attention to the fact that the number of mesh points required for a uniform mesh is certainly underestimated by the previous relation, because of the occurrence of extreme events in the dissipation leading to fluctuations at scale well below $\eta$. Thus, at this rate, even if the power of computers continues to grow exponentially, the DNS at Reynolds numbers relevant to the applications will remain

<footnote>$^6$The Reynolds number based on the large scales $Re \sim Re_\lambda^2$ doubles every 4 years</footnote>
<figure><img src="image_32.png" /><figcaption>Figure 3.9: Evolution of the mesh size and Reynolds number of the largest DNS over the years. From [193, 226, 64, 262, 234, 243, 95, 119, 113, 41]. Update from Ref. [52]. Comparison with exponential growth at a rate ln(2)/8 for Reλ vs. years.</figcaption></figure>

inaccessible for a long time.

The alternative proposed by the large eddy simulation (LES) approach is to decouple the required mesh size from the Reynolds number of the flow, by not resolving the smallest structures. This approach is justified because the largest structures (those that are resolved) contain most of the kinetic energy, are the most efficient in the transport phenomena while being very dependent on the flow configuration (its forcing and boundary conditions). The small (unresolved) structures are considered universal, at least when the turbulence can be assumed to be locally homogeneous and isotropic (e.g. far from the walls). Nevertheless, a model must be used for the effective action of small scales on large-scale dynamics.

The equation for the filtered velocity field $\bar{u}$ is obtained in two steps, similar to those used to study critical phenomena with the renormalization group technique [117]. The first step of “decimation” consists in eliminating the small scales: $u_i(\boldsymbol{x}, t) \rightarrow \bar{u}_i(\boldsymbol{x}, t) = \int H_\Delta(\boldsymbol{x}')u_i(\boldsymbol{x} - \boldsymbol{x}', t)$. where $H_\Delta$ represents the kernel of the filter at the $\Delta$ scale, typically the mesh scale. The second step of “renormalization” consists of modifying the interaction between the remaining degrees of freedom to leave unchanged the kinetic energy of the large scales and its rate of change (transfer rate). Classically, this is obtained by introducing an effective viscosity. We thus obtain the LES equation:

$$ \frac{\partial \bar{u}_i}{\partial t} + \bar{u}_j \frac{\partial \bar{u}_i}{\partial x_j} = - \frac{1}{\rho} \frac{\partial \bar{p}}{\partial x_i} + \frac{\partial}{\partial x_j} 2\nu_\Delta \bar{S}_{ij}, \quad (3.47) $$

with $S_{ij} = 1/2(\partial_j u_i + \partial_i u_j)$ the rate of strain of the flow. The effective viscosity $\nu_\Delta$ models the effect of unresolved scales on the resolved scales. It can be obtained from the mixing length model 7: $\nu_\Delta = \nu + (c_s\Delta)^2|S_{ij}|$ with $c_s \approx 0.2$ for HIT flows [87, 58].

When $\Delta \gg \eta$ the computation cost is thus much less than for DNS. But this has a price! At high Reynolds number, the dissipation is very intermittent. These very intense fluctuations of the velocity gradients over very short distances (and whose statistics depend on the Reynolds number) are not reproduced by the LES. The LES has a much lower effective Reynolds number $Re_{eff} = Re\nu/\nu_\Delta$. Is this problematic? Yes, if one is interested

<footnote>7This is the Smagorinski model, there are many other models for inhomogeneous flows [201, 232]</footnote>
in the coupling of the flow with other phenomena that depends on the small scales of the flow.

### 3.3.2 Large eddy simulations for two-phase flows

The interactions between particles and turbulence depend a priori on all scales of the flow. In order to calculate the hydrodynamic forces on the particles, it is necessary to have access to the turbulent fluctuations of the velocity field at the particle scale. But with a LES simulation the velocity field is filtered at the scale of the mesh $\Delta$, so that the fluctuations of the field between $\Delta$ and $d$ are not available. The naive approach then consists of simply ignoring the problem and calculating the evolution of the particles directly from $\bar{u}$. For example, in the case of a small particle ($d < \eta$ and $Re_p < 1$) discussed earlier:

$$ \begin{aligned} \boldsymbol{a}_p = \frac{d\boldsymbol{u}_p}{dt} &= \frac{\overline{d\boldsymbol{u}_p}}{dt} \\ &= -\frac{\boldsymbol{u}_p - \bar{\boldsymbol{u}}_f}{\tau_p} + \beta \frac{\overline{D\boldsymbol{u}_f}}{Dt} \end{aligned} \qquad (3.48) $$

With such an approach one expects that the highest frequencies are eliminated. In order to illustrate the problems of this approach, we compare in figures 3.10 and 3.11 the variance and the PDF of the particle acceleration obtained by LES with the integration of (3.48) with the one obtained with a finely resolved velocity field (DNS). This test is performed for small, heavy $\beta = 0$ and light $\beta = 3$ particles, for different Stokes numbers and different LES resolutions. We find that the acceleration of the particles is largely underestimated with the LES approach and that the LES resolution has a strong effect on the acceleration statistics. Moreover, we see that the PDF of the acceleration tends to a Gaussian shape when the small scales are not taken into account in the calculation of the forces. This reflects that the dynamics of the particles is clearly influenced by the small scales of the flow, and shows that it is necessary to introduce the effect of the small scales not resolved by the LES.

Different approaches have been considered to account for unresolved scales in the calculation of the drag force on small inertial particles. Many of them stochastically reconstruct the subgrid scales of the fluid velocity carrying the particles [27, 204, 39, 76, 172, 115, 237, 236, 194]. Most of these models are energetic in nature and aimed primarily at reconstructing a realistic spectrum of the velocity seen by the particles but do not exhibit an explicit dependence on the Reynolds number. Therefore, the effect of the intermittency of the flow in the dynamics of the particles cannot be reproduced. Moreover, none of the cited approaches has addressed the issue of the fluid inertia force, which is essential in the dynamics of light particles.

On the other hand, approaches have been proposed to reconstruct an Eulerian field which will then be interpolated to move the dispersed phase [43, 44, 123, 88]. Among these approaches, Ref. [229] introduced the decomposition of the instantaneous fluid acceleration field into a filtered contribution and a random contribution to account for intermittency at small unresolved scales (see also Refs. [230, 276, 228]).

In Refs. [92, 279] it was proposed to consider a similar decomposition for the instantaneous
<figure><img src="image_33.png" /><figcaption>Figure 3.10: Variance of the acceleration of small particles as a function of the Stokes number. Left: for heavy particles ($\beta = 0$) in THI for $Re_\lambda = 420$, dashed lines: DNS (with $N = 2048^3$ mesh), circles: LES without model, and cross: LES with stochastic model, for $N = 64^3$, $96^3$ and $128^3$ from black to gray. Right: for light particles ($\beta = 3$) in THI for $Re_\lambda = 200$, blue stars: DNS (with $N = 1024^3$ meshes), solid symbols: LES without model, open symbols: LES with stochastic model, with $N = 64^3$ (diamonds), $48^3$ (triangles), $32^3$ (circles)</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Figure 3.11: Normalized PDFs of small particle acceleration for DNS, LES without model and LES with model. Left: for heavy particles ($\beta = 0$) in THI for $Re_\lambda = 420$ for $St = 0.016, 0.16, 0.6, 1, 2, 3, 5, 10, 20$ (respectively, shifted up by one decade relative to each other for clarity), black dashed lines: DNS (with $N = 2048^3$ meshes), continuous black lines: LES without model, and continuous gray lines: LES with stochastic model, for $N = 64^3$. Right: for light particles ($\beta = 3$) in THI for $Re_\lambda = 200$ for $St = 0.02, 0.074, 0.20, 0.45, 1.01, 1.55$ and $2.07$. Continuous black line: DNS (with $N = 1024^3$ meshes), blue dotted lines: LES without model, continuous red lines: LES with stochastic model, for $N = 64^3$.</figcaption></figure>
acceleration of particles:

$$ \boldsymbol{a}_p = \underbrace{\frac{d\boldsymbol{u}_p}{dt}}_{\text{Filtered acceleration}} + \underbrace{\boldsymbol{a}_*}_{\text{Stochastic acceleration}} . \qquad (3.49) $$

The first term represents the response of the particles at large scales governed by the fluid velocity field resolved by the LES. The second, the residual acceleration, reproduces the effect of fluctuations at unresolved scales in the particle momentum balance. According to the Kolmogorov 62 hypothesis, fluctuations at unresolved scales are mainly attributed to fluctuations in the rate of energy transfer to smaller scales.

Based on the separation of timescales for the evolution of the norm of the acceleration and its orientation [199, 184], we consider here that the residual acceleration of the particles is given by the product of two independent stochastic processes, one for the norm, the other for the orientation.

$$ \boldsymbol{a}_* = a_* e_* \qquad (3.50) $$

The norm of the residual acceleration is obtained from the local value of the energy transfer rate “seen” by the particle along its trajectory. This transfer rate at the particle scale is not accessible from the velocity field computed by LES and will be obtained by a stochastic log-normal model, from the energy transfer rate resolved by the LES. The stochastic process for the norm should present a correlation over a time of the order of the integral scale. The process for the orientation is obtained from a random walk on a sphere and presents a fast decorrelation on the order of the dissipative time.

## Model for the norm

More specifically for the stochastic process for the norm we consider, according to the local similarity hypothesis (K62), that the main source of fluctuations comes from the fluctuations of the dissipation rate [128, 49] and that the statistics of the acceleration of the particle conditioned to the value of the dissipation seen by the particles are universal and depend only on the characteristics of the particle (size, density)

$$ \langle \boldsymbol{a}_p^2 \rangle = \int \langle \boldsymbol{a}_p^2 | \varepsilon_p \rangle P(\varepsilon_p) d\varepsilon_p . \qquad (3.51) $$

In this relation $\varepsilon_p$ represents the dissipation rate filtered on the particle timescale. It has been shown that for particles and bubbles [92, 279] the acceleration conditioned by $\varepsilon_p$ evolves as a power law for $\varepsilon_p/\langle\varepsilon\rangle \gg 1$:

$$ \langle \boldsymbol{a}_p^2 | \varepsilon_p \rangle / \langle \boldsymbol{a}_p^2 \rangle \approx (\varepsilon_p / \langle \varepsilon \rangle)^{3/2} . \qquad (3.52) $$

This relation is consistent with the case of fluid particles at large Reynolds number shown in the first part (see figure 2.8). We then propose to write the norm of the subgrid acceleration as

$$ \begin{aligned} a_* &= \langle \boldsymbol{a}_p^2 | \varepsilon_* \rangle^{1/2} \\ &= a_\eta \Gamma_a^{1/2} (St, \beta, Re_\Delta)^{1/2} \frac{\varepsilon_*^{3/4}}{\langle \varepsilon \rangle^{3/4}} . \end{aligned} \qquad (3.53) $$
To obtain the second line in the previous equation we used the estimate of the acceleration variance $\Gamma_a(St, \beta, Re_\Delta)$ introduced previously (3.21) (or (3.37) to take into account the finite Reynolds number effects). In order to take into account only the non-resolved part of the fluctuations, the function $\Gamma_a$ is here evaluated with the Reynolds number of the subgrid scale of the flow ($\tau_L/\tau_\eta \to \tau_\Delta/\tau_\eta = (\Delta/\eta)^{2/3} = Re_\Delta^{1/2}$, where we introduced the characteristic time of the largest unresolved scales $\tau_\Delta = \Delta^2/\nu_\Delta$.

The dissipation rate along the trajectory of the particle $\epsilon_*$ is considered as a stochastic variable. Its evolution is obtained here with the stochastic log-normal model of Pope [203] (also discussed in the first part (2.18)). From the Ito transformation, we obtain the stochastic process for $\epsilon_*^{3/4}$ which depends on the local value of the transfer rate $\epsilon_\Delta(t) = \nu_\Delta(\nabla \bar{u})^2$ computed at the scale of the LES mesh and interpolated to the particle position.

$$ \frac{d\epsilon_*^{3/4}}{\epsilon_*^{3/4}} = \frac{d\epsilon_\Delta^{3/4}}{\epsilon_\Delta^{3/4}} - \left( \ln\left(\frac{\epsilon_*^{3/4}}{\epsilon_\Delta^{3/4}}\right) - \frac{3}{16}\sigma^2 \right) \frac{dt}{\tau_\Delta} + \sqrt{\frac{9\sigma^2}{8\tau_\Delta}} dW. \quad (3.54) $$

In this equation $dW$ is the increment of the Wiener process and $d\epsilon_\Delta^{3/4}$ is the increment of $\epsilon_*^{3/4}$ along the particle trajectory. This stochastic process ensures that $\langle \epsilon_* \rangle = \langle \epsilon_\Delta \rangle$ (see [92]). The parameter $\sigma^2 = 1/2 \ln \frac{\tau_\Delta}{\tau_p+\tau_\eta}$ gives the depth of the cascade. The value of the coefficient in front of the logarithm is set to reproduce the Reynolds number dependence reported by [271]. The timescale $\tau_\Delta$ imposes the temporal correlation of $\epsilon_*$. As discussed previously this model does not reproduce satisfactorily all the constraints of the energy cascade, but it is used here for simplicity and also because the work presented in the first part, which is an improvement of this process, was done afterwards.

To illustrate the process (3.54) we compare in figure 3.12 a realization of the dissipation rate $\epsilon$ along the trajectory of a solid particle at $St = 1$ and $Re_\lambda = 420$ obtained the DNS of [140] with the evolution of $\epsilon_\Delta$ obtained by LES under the same conditions with $N = 64^3$ and the realization of $\epsilon_*$ obtained from stochastic process (3.54) and $\epsilon_\Delta$. We observe that $\epsilon_\Delta$ (LES) shows much smaller fluctuations than $\epsilon$ (DNS), but that the stochastic model is able to reproduce well the intermittency of the dissipation rate for $\epsilon_*$ while ensuring some correlation with $\epsilon_\Delta$.

## Model for the orientation

With the stochastic model for orientation, we first want to impose the temporal correlation of the components of the acceleration vector. In the case of weakly inertial particles, this one persists over a time of the order of the Kolmogorov scale, much shorter than the correlation scale for the norm. Thus the correlation of the components of the acceleration will be given by that of its orientation. We thus propose to model the orientation of the subgrid scale force by a random walk on a sphere with a unit radius:

$$ de_* = e_* \times \alpha dt \quad (3.55) $$

In this model the position $e_*$ is modified by a random angular displacement $\alpha dt$. In the simplest case, when the inertia force is neglected ($\beta = 0$), the angular velocity $\alpha$ follows a diffusive process including restoring, damping and diffusion terms:

$$ d\alpha = -e_* \times \frac{d\mathbf{u}_p/dt}{\Delta} - \alpha \frac{dt}{\tau_e} + \sqrt{\frac{\sigma_e^2}{\tau_e}} d\mathbf{W} \quad (3.56) $$
<figure><img src="image_35.png" /><figcaption>Figure 3.12: (Left) A realization of the evolution of $\epsilon$ along a particle trajectory at $St = 1$ and $Re_\lambda = 420$, from the DNS of [141]. (right) In red, a realization of the evolution of $\epsilon_\Delta$ along a particle trajectory at $St = 1$ and $Re_\lambda = 420$, from a LES with a mesh size of $N = 64^3$; in black, a realization of the stochastic process for $\epsilon_*$.</figcaption></figure>

The diffusion term, in which $dW$ is a 3D Wiener process, ensures the tendency to return to a statistically isotropic orientation by allowing the vector $e_*$ to take all possible orientations (i.e. to visit the whole sphere). The time to return to isotropy corresponds to the characteristic time for the “walker” to invade the whole sphere and is therefore related to its characteristic (angular) velocity $\sigma_e$. The damping term is necessary to control the rate of return to isotropy at short times and thus guarantee that the orientation remains correlated over a time of the order of $\tau_e$. We chose $\tau_e = \frac{1}{2}(\tau_p + \tau_\eta)) = 1/\sigma_e$ so that in the case of fluid particles ($\tau_p = 0$) the orientation correlation time is of the order of the Kolmogorov time, and that for very inertial particles ($\tau_p \gg \tau_\eta$) the orientation of the particle remains correlated over a time of the order of $\tau_p$ reflecting the more and more ballistic character of its trajectory. Finally, the restoring term tends to align $e_*$ with the orientation of the resolved contribution $\overline{d_t u_p}/|\overline{d_t u_p}|$, and thus counteracts the diffusive tendency. Since the return to small-scale isotropy results from the cascade, the alignment between the resolved acceleration and the subgrid acceleration depends on the ratio between the resolved scale and the Kolmogorov scale. We assume that the strength of the restoring force is estimated as $|\overline{d_t u_p}|/\Delta$. For $\tau_e \ll (|\overline{d_t u_p}|/\Delta)^{1/2}$ the orientation becomes independent of the large-scale orientation, which is consistent with the assumption of local isotropy on small scales at large Reynolds number. On the other hand, for $\tau_e \gg (|\overline{d_t u_p}|/\Delta)^{1/2}$ the orientation is imposed by the resolved scales.

An illustration of this diffusive process on the sphere is shown in figure 3.13.

In the case where the inertial force is not neglected ($\beta > 0$), we wish to take into account the correlation between the inertia and drag forces. Indeed, in Ref. [279], we showed that the strong increase in the flatness of the particle acceleration around $St = 1$ for light particles ($\beta = 3$) should be attributed to the correlation between these forces and their relative orientation rather than to preferential concentration effects. For $St \approx 1$ the two forces $F_D$ and $F_I$ are typically of the same order of magnitude and generally anti-aligned, leading to a relatively small acceleration, whereas when these forces tend to align a high acceleration will result. To reproduce this geometric effect we introduced the orientation
<figure><img src="image_36.png" /><figcaption>Figure 3.13: A realization of the stochastic process (3.55) and (3.56) on the sphere.</figcaption></figure>

of these two forces, respectively $e_I^*$ and $e_D^*$:

$$a_* = |F_D|e_D^* + |F_I|e_I^* \qquad (3.57)$$

where $|F_D|$ and $|F_I|$ represent the norm of the subgrid contribution of the drag and inertia forces. We then write:

$$a_* = a_\eta \frac{\epsilon_*^{3/4}}{\langle \epsilon \rangle^{3/4}} (\Gamma_D^{1/2} e_D^* + \Gamma_I^{1/2} e_I^*) \qquad (3.58)$$

where the intensity of the drag and inertia forces are estimated by (3.25), and (3.26) and considering that their conditional variance both vary as $\epsilon_*^{3/4}$. By identifying (3.58) with (3.50) and (3.53) we can write the vector $e^*$ as $e^* = (\Gamma_D^{1/2}e_D^* + \Gamma_I^{1/2}e_I^*)/\Gamma_a^{1/2}$. It should be noted, however, that in this case the vector $e^*$ is no longer an orientation vector because it does not necessarily have a unit norm. It depends on the orientation between the two forces and their correlation:

$$e_*^2 = (\Gamma_D + \Gamma_I + 2\Gamma_D^{1/2}\Gamma_I^{1/2}e_D^*e_I^*) / \Gamma_a = 1 + \frac{2\Gamma_D^{1/2}\Gamma_I^{1/2}}{\Gamma_a} (e_D^*e_I^* - \langle e_I.e_D \rangle) \qquad (3.59)$$

To obtain the second equality, we assumed the separation of the timescales of orientations and norms which allows us to deduce that the average orientation between $e_I$ and $e_D$ is given by:

$$\langle e_I.e_D \rangle = \frac{\langle F_D.F_I \rangle}{(\langle F_D^2 \rangle\langle F_I^2 \rangle)^{1/2}} = \frac{1}{2} \frac{\Gamma_a - \Gamma_D - \Gamma_I}{(\Gamma_D\Gamma_I)^{1/2}}. \qquad (3.60)$$

This assumption was verified in Ref. [279].

To respect the constraint (3.60) in the modeling of $e_*$, we have introduced two, coupled, stochastic processes for $e_I^*$ and $e_D^*$. The evolution of the two orientation vectors $e_I^*$ and $e_D^*$ is given by a model analogous to (3.55)-(3.56), but the restoring term of $e_D^*$ imposes an equilibrium configuration where $e_I^*$ and $e_D^*$ are anti-aligned. The details are given in Ref. [279]. This model recovers the very strong increase in flatness observed by DNS around $St = 1$. It is interesting to note that such a mechanism has already been proposed concerning the antialignment between $\partial_t u$ and $u \nabla u$ for a fluid particle [259]. Note that in the limit $\beta \to 0$ this model degenerates to the orientation model detailed above.
## Contribution of the resolved scales

The first term of the equation (3.49) represents the contribution of large scales to the acceleration. It corresponds to the acceleration of the particles calculated from the LES field filtered spatially and temporally: In the case where the particles are small we write:

$$ \frac{d\overline{u}_p}{dt} = G_{\Delta, \tau_{\Delta}} * \frac{d\overline{u}_p}{dt} = \frac{\overline{u}_f - u_p}{\max(\tau_{\Delta}, \tau_p)} + \beta \frac{D\overline{u}_f}{Dt} \quad (3.61) $$

where $G_{\Delta, \tau_{\Delta}}$ represents the convolution by a space-time filter eliminating the wave numbers higher than $\Delta$ and the frequencies higher than $\tau_{\Delta} = \Delta^2/\nu_{\Delta}$. In the context of the LES, we can consider that the spatial filtering is implicit because of the mesh resolution with scale $\Delta$. On the other hand, the time step of the simulation, which is fixed by the CFL criterion, can be much smaller than $\tau_{\Delta}$ ($\tau_{\Delta}/dt \sim Re_{eff}^{1/4}$). Thus to avoid that the resolved term develops high frequencies in response to the stochastic forcing $a_*$, we explicitly filter the high frequencies of this term. This is simply done by replacing the response time of the drag term by the characteristic time of the smallest resolved scales $\tau_{\Delta}$, if $\tau_p < \tau_{\Delta}$. Note that in the absence of stochastic forcing this time filter does not modify the dynamics of the particle since the particle inertia effectively filters the flow fluctuations.

## Results for small particles

Figures 3.10 and 3.11 present the variance and PDFs of particle acceleration obtained by LES with the stochastic model compared to those obtained by DNS and by LES without the subgrid model. We consider in these two figures the case of heavy ($\beta = 0$) and light ($\beta = 3$) particles and a range of Stokes numbers. We can see that with the stochastic model the LES reproduces these statistics in very good agreement with the DNS, contrary to the LES without submesh model, as already discussed at the beginning of the paragraph. Moreover, we find that the agreement with the DNS remains good even when the resolution of the LES is strongly reduced. In the case of inertial particles, we can see on figures 3.10 and 3.11 (left) that the LES without model approaches the DNS when $St$ increases, meaning that particles with high inertia filter out high frequencies of the flow.

Other interesting results have been obtained. It is shown in Ref. [279] that the variance of the particle velocity increments (equivalent to the Lagrangian velocity spectra) obtained by the LES with the stochastic model becomes very close to the DNS while without the model the LES gives a strong underestimation.

We have studied the variance of the particle acceleration conditioned on the instantaneous value of the dissipation rate $\epsilon_*$ estimated from the stochastic model. It was observed in Ref. [279] that, similarly to the DNS results, the particle acceleration from the LES with the model increases as $\epsilon_*^{3/2}$ for large values of the dissipation rate, and is independent of $\epsilon_*$ for small values of the latter. This behavior shows that the acceleration of particles in weakly dissipative regions is given by the resolved contribution $\overline{d_t u_p}$, while larger fluctuations in the acceleration are correctly reproduced with the stochastic model.

We are also interested in the power exchanged between the particles and the fluid $P = a_p \cdot u_p$. This quantity presents very strong fluctuations, much larger than $\langle \epsilon \rangle$. With the stochastic model the LES is able to reproduce these large fluctuations, both the variance and the shape of the PDFs as can be seen on figure 3.14, whereas the LES clearly underestimates this quantity when the Stokes number is small. Nevertheless, we can notice
<figure><img src="image_37.png" /><figcaption>Figure 3.14: (Left) Variance of the power exchange with the fluid per unit mass of a particle normalized by the mean dissipation rate in logarithm scale. Crosses: LES with the proposed model; circles: LES without the model; for the three meshes: 64³ (black), 96³ (dark gray) and 128³ (light gray). Comparison with the DNS of [21] in dashed lines. Inset zoom for small Stokes number in linear scale. (Right) PDF of the power exchange with the fluid normalized by its variance for $St = 0.016, 0.16, 0.6, 1, 2, 3, 5, 10, 20$ (respectively shifted upward by one decade from each other for clarity). gray: LES at 64³ with the proposed model; black: LES at 64³ without the model; black dashed line: DNS of [21] and Gaussian distribution in gray dashed line.</figcaption></figure>

that the LES with the stochastic model does not succeed in reproducing the asymmetry of the PDF observed for small Stokes numbers, observed by DNS. This asymmetry is linked to the temporal irreversibility of the flow (the power $P$ being the variation of the kinetic energy of the particles per unit of time). We have shown in the first part that for fluid particles this power asymmetry can be reproduced by a stochastic model and that it is linked to both the winding of the particles and the non-Markovian character of the dissipative model. This observation opens up new possibilities of improvement for the subgrid models by adapting the developments proposed in the first part to the case of inertial particles. On the other hand, it is interesting to note that when the Stokes number increases, the asymmetry of the PDF becomes much more marked and that the LES is able to reproduce it even without the addition of a stochastic model.

## Results for large particles

For particles that are large compared to the Kolmogorov scale, we propose to start from (3.44) based on the effective relaxation time of the particle taking into account its finite size. In the case of the LES we write, assuming that the size of the particle remains much