2 Formulation of the problem

Many studies have considered the dynamics of heavy inertial particles, much denser than the carrying fluid. When the particles are very heavy and at the same time very small (point particles) the motion is simply determined by Stokes’ drag. The dynamics of light particles by contrast is also affected by pressure gradients of the unperturbed fluid, and added-mass effects. Neglecting the effect of gravitational settling (and thus buoyancy) the equation of motion reads:

(2.1) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}l@{}}\dot{\boldsymbol{r}}_{t}=\boldsymbol{v}_{t},\\ \dot{\boldsymbol{v}}_{t}={\it\beta}\text{D}_{t}\boldsymbol{u}(\boldsymbol{r}_{t},t)+(\boldsymbol{u}(\boldsymbol{r}_{t},t)-\boldsymbol{v}_{t})/{\it\tau}_{s},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{r}_{t}$ is the particle position at time $t$ , $\boldsymbol{v}_{t}$ , is the particle velocity, $\boldsymbol{u}(\boldsymbol{r}_{t},t)$ is the velocity field of the undisturbed fluid and $\text{D}_{t}\boldsymbol{u}=\partial _{t}\boldsymbol{u}+(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{u}$ is the Lagrangian derivative. The dimensionless constant ${\it\beta}=3{\it\rho}_{f}/({\it\rho}_{f}+2{\it\rho}_{p})$ accounts for the contrast between particle density ${\it\rho}_{p}$ and fluid density ${\it\rho}_{f}$ , while the Stokes number is defined as $St={\it\tau}_{s}/{\it\tau}_{K}$ where the particle response time is ${\it\tau}_{s}=R^{2}/(3{\it\nu}{\it\beta})$ , $R$ is the particle radius, ${\it\nu}$ the kinematic viscosity of the flow and ${\it\tau}_{K}=\sqrt{{\it\nu}/{\it\epsilon}}$ the Kolmogorov time defined in terms of the fluid energy dissipation, ${\it\epsilon}$ . Many studies have employed this model (Babiano et al. Reference Babiano, Cartwright, Piro and Provenzale2000; Bec et al. Reference Bec, Celani, Cencini and Musacchio2005; Calzavarini et al. Reference Calzavarini, Cencini, Lohse and Toschi2008a , Reference Calzavarini, Volk, Bourgoin, Leveque, Pinton and Toschi2009; Volk et al. Reference Volk, Calzavarini, Verhille, Lohse, Mordant, Pinton and Toschi2008). The model takes into account added-mass effects but neglects buoyancy forces. It is an open question under which circumstances this is a quantitative model for the acceleration of small particles in turbulence.

The following analysis is based on (2.1). It is important to first understand this case before addressing more realistic situations (inclusion of buoyancy, finite particle size, collisions and feedback on the flow). The sources of the difficulties are twofold. First, even in the much simpler case of a non-turbulent Eulerian velocity field $\boldsymbol{u}(\boldsymbol{r},t)$ the particle dynamics is still complicated and often chaotic, simply because equation (2.1) are nonlinear. Second, turbulence makes the problem even harder due to the existence of substantial spatial and temporal fluctuations. This results in chaotic Lagrangian dynamics of fluid elements. Note that even though Lagrangian fluid elements sample space uniformly, their instantaneous motion is in general correlated with structures in the underlying flow. This implies that multi-time Lagrangian flow correlation functions evaluated along tracer trajectories do not coincide with the underlying Eulerian correlation functions which are evaluated at fixed positions in space. Particles that are heavier or lighter than the fluid may detach from the flow if they have inertia. This leads to the inertial preferential sampling mentioned in the introduction. As a result, the flow statistics experienced by an inertial particle differs from that of a Lagrangian fluid element.

Figure 1. Parameter space showing available DNS data sets with $Re_{{\it\lambda}}=185$ and $Re_{{\it\lambda}}=400$ . For values of $St$ smaller than $0.05$ the $x$ axis is linear, while it is logarithmic for values of $St$ larger than $0.05$ . For each data set with $Re_{{\it\lambda}}=185$ we have analysed a total of 130 000 trajectories of duration $6T_{L}$ and for each data set with $Re_{{\it\lambda}}=400$ we have analysed a total of 200 000 trajectories of duration $2.5T_{L}$ . Level curves $St(1-{\it\beta})={\it\epsilon}$ for constants ${\it\epsilon}=\{-1,-0.1,-0.01,0,0.01,0.1,1\}$ are plotted as black lines. Parameter families: with $Re_{{\it\lambda}}=185$ : ${\it\beta}=0$ (red, ○), ${\it\beta}=0.25$ (magenta, ▵), ${\it\beta}=0.5$ (cyan, ▿), ${\it\beta}=0.75$ (green, ▹), ${\it\beta}=1$ (black, ▫), ${\it\beta}=1.25$ (brown, ♢), ${\it\beta}=1.5$ (purple, ▿), ${\it\beta}=2$ (orange, ▹), ${\it\beta}=2.5$ (dark green, ◃), ${\it\beta}=3$ (blue, ▵). With $Re_{{\it\lambda}}=400$ : ${\it\beta}=0$ (red,  

). Additional data from Calzavarini et al. (Reference Calzavarini, Volk, Bourgoin, Leveque, Pinton and Toschi2009), Prakash et al. (Reference Prakash, Tagawa, Calzavarini, Mercado, Toschi, Lohse and Sun2012):  ${\it\beta}=3$ (blue,  

).

3 Lagrangian closure

Let us notice that the dynamics determined by equation (2.1) tends to the evolution of a tracer in both limits $St\rightarrow 0$ and ${\it\beta}\rightarrow 1$ . Indeed, in the limit ${\it\tau}_{s}\rightarrow 0$ , imposing a finite Stokes drag leads to $\boldsymbol{v}_{t}=\boldsymbol{u}+O({\it\tau}_{s})$ . For ${\it\beta}=1$ , by contrast we may approximate the material derivative along fluid-trace trajectories with the derivative along the trajectory of a particle, $\text{D}_{t}\boldsymbol{u}\sim \dot{\boldsymbol{u}}$ , and consistently check that the evolution of (2.1) leads to an exponential relaxation of the evolution of the particle to the trajectory of the tracer. Hence, the idea is to approximate the effects of inertial forces on the particle trajectory starting from the evolution of tracers. This approach cannot be exact. It is for example known that vortices are preferentially sampled by inertial particles. Nevertheless, it is important to understand and quantify how big the difference is as a function of the distance in the parameter phase space, $(St,{\it\beta})$ , from the two lines $St=0$ or ${\it\beta}=1$ where the closure must be exact (see figure 1). The starting point is to evaluate (2.1) along tracer trajectories:

(3.1) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}l@{}}\dot{\boldsymbol{r}}_{t}^{(L)}=\boldsymbol{u}(\boldsymbol{r}_{t}^{(L)},t),\\ \dot{\boldsymbol{v}_{t}}={\it\beta}\text{D}_{t}\boldsymbol{u}(\boldsymbol{r}_{t}^{(L)},t)+(\boldsymbol{u}(\boldsymbol{r}_{t}^{(L)},t)-\boldsymbol{v}_{t})/{\it\tau}_{s},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{r}_{t}^{(L)}$ denotes the Lagrangian trajectory of a tracer particle. This approximation is a closure in the sense that the first equation in (3.1) is independent of the second one, and the second equation can be solved in terms of the Lagrangian velocity statistics $\boldsymbol{u}(\boldsymbol{r}_{t}^{(L)},t)$ of the underlying flow. We solve (3.1) for $\boldsymbol{v}_{t}$ (disregarding initial conditions here and below because these do not matter for the steady-state statistics) to obtain

(3.2) $$\begin{eqnarray}\displaystyle \boldsymbol{v}_{t}={\it\beta}\boldsymbol{u}(\boldsymbol{r}_{t}^{(L)},t)+\frac{1-{\it\beta}}{{\it\tau}_{s}}\int _{0}^{t}\text{d}t_{1}\,\text{e}^{(t_{1}-t)/{\it\tau}_{s}}\boldsymbol{u}(\boldsymbol{r}_{t_{1}}^{(L)},t_{1}). & & \displaystyle\end{eqnarray}$$

Using this expression, the particle acceleration follows from the second of equation (3.1)

(3.3) $$\begin{eqnarray}\displaystyle \boldsymbol{a}_{t}={\it\beta}\frac{\text{D}\boldsymbol{u}}{\text{D}t}(\boldsymbol{r}_{t}^{(L)},t)+\frac{1-{\it\beta}}{{\it\tau}_{s}}\int _{0}^{t}\text{d}t_{1}\,\text{e}^{(t_{1}-t)/{\it\tau}_{s}}\frac{\text{D}\boldsymbol{u}}{\text{D}t}(\boldsymbol{r}_{t_{1}}^{(L)},t_{1}). & & \displaystyle\end{eqnarray}$$

Using equation (3.3) we express two-time acceleration statistics in terms of the two-point Lagrangian acceleration correlation function

(3.4) $$\begin{eqnarray}\displaystyle C_{L}(t)\equiv \langle \text{D}_{t}\boldsymbol{u}(\boldsymbol{r}_{t}^{(L)},t)\boldsymbol{\cdot }\text{D}_{t}\boldsymbol{u}(\boldsymbol{r}_{0}^{(L)},0)\rangle . & & \displaystyle\end{eqnarray}$$

In the steady-state limit we find

(3.5) $$\begin{eqnarray}\displaystyle \langle \boldsymbol{a}_{t}\boldsymbol{\cdot }\boldsymbol{a}_{0}\rangle & = & \displaystyle {\it\beta}^{2}C_{L}(t)+\frac{1-{\it\beta}^{2}}{{\it\tau}_{s}}\left[\cosh \left[\frac{t}{{\it\tau}_{s}}\right]\int _{t}^{\infty }\text{d}t_{1}\,\text{e}^{-t_{1}/{\it\tau}_{s}}C_{L}(t_{1})\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,\text{e}^{-t/{\it\tau}_{s}}\int _{0}^{t}\text{d}t_{1}\,\cosh \left[\frac{t_{1}}{{\it\tau}_{s}}\right]C_{L}(t_{1})\right].\end{eqnarray}$$

The acceleration variance is obtained by letting $t\rightarrow 0$ in (3.5)

(3.6) $$\begin{eqnarray}\displaystyle \langle \boldsymbol{a}^{2}\rangle ={\it\beta}^{2}C_{L}(0)+\frac{1-{\it\beta}^{2}}{{\it\tau}_{s}}\int _{0}^{\infty }\text{d}t_{1}\,\text{e}^{-t_{1}/{\it\tau}_{s}}C_{L}(t_{1}). & & \displaystyle\end{eqnarray}$$

Similarly, the fourth moment of the particle acceleration is obtained as:

(3.7) $$\begin{eqnarray}\displaystyle \langle |\boldsymbol{a}|^{4}\rangle & = & \displaystyle {\it\beta}^{4}C_{L}(0,0,0)-4\frac{({\it\beta}-1){\it\beta}^{3}}{{\it\tau}_{s}}\int _{0}^{\infty }\text{d}t_{1}\,\text{e}^{-t_{1}/{\it\tau}_{s}}C_{L}(-t_{1},0,0)\nonumber\\ \displaystyle & & \displaystyle +\,6\frac{({\it\beta}-1)^{2}{\it\beta}^{2}}{{\it\tau}_{s}^{2}}\int _{0}^{\infty }\text{d}t_{1}\,\int _{0}^{\infty }\text{d}t_{2}\,\text{e}^{-(t_{1}+t_{2})/{\it\tau}_{s}}C_{L}(-t_{1},-t_{2},0)\nonumber\\ \displaystyle & & \displaystyle -\,\frac{({\it\beta}-1)^{3}(3{\it\beta}+1)}{{\it\tau}_{s}^{3}}\int _{0}^{\infty }\text{d}t_{1}\,\int _{0}^{\infty }\text{d}t_{2}\,\int _{0}^{\infty }\text{d}t_{3}\,\text{e}^{-(t_{1}+t_{2}+t_{3})/{\it\tau}_{s}}C_{L}(-t_{1},-t_{2},-t_{3}).\quad \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

Here isotropy of the acceleration components $\langle a_{i}a_{j}a_{k}a_{l}\rangle =[{\it\delta}_{ij}{\it\delta}_{kl}+{\it\delta}_{ik}{\it\delta}_{jl}+{\it\delta}_{il}{\it\delta}_{jk}]\langle a_{1}^{4}\rangle /3$ was used to express $\langle |\boldsymbol{a}|^{4}\rangle$ in terms of the four-point Lagrangian correlation function

(3.8) $$\begin{eqnarray}\displaystyle C_{L}(t_{1},t_{2},t_{3})\equiv \frac{d(d+2)}{3}\langle \text{D}_{t}u_{1}(\boldsymbol{r}_{t_{1}}^{(L)},t_{1})\text{D}_{t}u_{1}(\boldsymbol{r}_{t_{2}}^{(L)},t_{2})\text{D}_{t}u_{1}(\boldsymbol{r}_{t_{3}}^{(L)},t_{3})\text{D}_{t}u_{1}(\boldsymbol{r}_{0}^{(L)},0)\rangle , & & \displaystyle\end{eqnarray}$$

where $\text{D}_{t}u_{1}$ is a component of the fluid acceleration and $d$ is the spatial dimension. Equations (3.5)–(3.7) express the fluctuations of inertial-particle accelerations in terms of Lagrangian correlation functions of the underlying flow, and predict how the inertial-particle accelerations depend on ${\it\beta}$ and $St$ . The integrals in (3.6) and (3.7) are hard to evaluate numerically for very small and for very large values of  $St$ . When $St$ is small, the exponential factors in the right-hand side of (3.6) and (3.7) become singular and one needs a very high sampling frequency of the fluid acceleration along the particle trajectory to evaluate the integrals reliably. When $St$ is large, on the other hand, the integrals sum up large-time contributions with large fluctuations due to the finite length of experimental and numerical trajectories.

Before we proceed to a quantitative assessment of the model let us make a few general remarks about the range of applicability and accuracy of the approximation made. First, in (3.6) and (3.7) we need to evaluate the acceleration correlation function of the tracers for general values of $t$ . This might be seen as a problem because the trajectories of tracers and of inertial particles must depart on a time scale of the order of the Lyapunov time. On the other hand, the acceleration correlation function of the tracer is known to decay on a time of the order of the Kolmogorov time, ${\it\tau}_{K}$ , which is ten times smaller than the typical Lyapunov time of heavy particles (Bec et al. Reference Bec, Biferale, Boffetta, Cencini, Musacchio and Toschi2006b ) and stays close to zero in the inertial range of scales (Falkovich et al. Reference Falkovich, Xu, Pumir, Bodenschatz, Biferale, Boffetta, Lanotte and Toschi2012). Therefore we do not see any problem in the evaluation of the integrals in the above equations. Second, we are only interested in stationary properties of the inertial particle statistics and we always assume that the initial conditions for all particles are chosen from their stationary distribution functions. As a result no influence of the initial condition should appear in the closure. The Lagrangian closure adopted here can also be applied to correlation functions of other observables of the particle, or of the flow, provided that the underlying Lagrangian correlation functions decay quickly enough. Here, we focus on the acceleration statistics because of their highly non-Gaussian properties and high sensitivity to the parameters ${\it\beta}$ and $St$ . Finally, let us stress that our closure is fully based on Lagrangian properties and does require the use of any Eulerian correlation function, as opposed to closures based on the fast Eulerian approach by Ferry & Balachandar (Reference Ferry and Balachandar2001) or the Lagrangian–Eulerian closure to predict two-particle distributions (Derevich Reference Derevich2000; Zaichik & Alipchenkov Reference Zaichik and Alipchenkov2003; Alipchenkov & Zaichik Reference Alipchenkov and Zaichik2007; Pan & Padoan Reference Pan and Padoan2010). To make further progress we need to determine the Lagrangian correlation functions. In § 4 we present the most important new results, we determine the Lagrangian correlation functions by DNS of inertial particle dynamics in turbulence, substitute into (3.5)–(3.7) and assess the accuracy of these equations in predicting particle acceleration fluctuations and correlations.

4 Direct numerical simulations

4.1 Simulation method

We present here the analysis of data obtained from DNS of a homogeneous and isotropic turbulent flow seeded with point-like particles with different inertia and particle–fluid density ratios (see figure 1 for a summary of the $(St,{\it\beta})$ values available). The data set was previously obtained by Bec et al. (Reference Bec, Biferale, Lanotte, Scagliarini and Toschi2010b ). The flow obeys the Navier–Stokes equations for an incompressible velocity field $\boldsymbol{u}(\boldsymbol{x},t)$ :

(4.1) $$\begin{eqnarray}\displaystyle \partial _{t}\boldsymbol{u}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}=-\boldsymbol{{\rm\nabla}}p+{\it\nu}{\rm\nabla}^{2}\boldsymbol{u}+\boldsymbol{f},\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0. & & \displaystyle\end{eqnarray}$$

The external forcing $\boldsymbol{f}$ is statistically homogeneous, stationary and isotropic, injecting energy in the first low wavenumber shells by keeping their spectral content constant (Chen et al. Reference Chen, Doolen, Kraichnan and She1993). The viscosity ${\it\nu}$ is set such that the Kolmogorov length scale ${\it\eta}\approx {\it\delta}x$ , where ${\it\delta}x$ is the grid spacing. The numerical domain is $2{\rm\pi}$ -periodic in the three directions. We use a fully dealiased pseudospectral algorithm with second-order Adams–Bashforth time stepping. For details see Bec et al. (Reference Bec, Biferale, Boffetta, Celani, Cencini, Lanotte, Musacchio and Toschi2006a ), Cencini et al. (Reference Cencini, Bec, Biferale, Boffetta, Celani, Lanotte, Musacchio and Toschi2006). Two series of DNS are analysed: Run I, with a numerical resolution of $512^{3}$ grid points and the Reynolds number at the Taylor scale $Re_{{\it\lambda}}\approx 200$ ; Run II, with $2048^{3}$ resolution and $Re_{{\it\lambda}}\approx 400$ . Details can be found in table 1.

Figure 2. Acceleration variance (a) and flatness (b) for DNS data at changing ${\it\beta}$ , $St$ and $Re$ . (a) Points correspond to the DNS data (same symbols of figure 1). Solid lines (labelled with their corresponding value of ${\it\beta}$ ) show the closure scheme prediction for the acceleration variance of light and heavy particles, (3.6), normalized with the fluid variance, for all data from Run I; dashed line corresponds to the closure for Run II. (b) The flatness measured on the DNS data and the one predicted by the closure scheme, (3.7), for heavy and light particles. Thin dashed black line shows the limit of normal distributed acceleration components. Additional data for ${\it\beta}=3$ (blue,  

) are omitted in (b) because the flatness was not evaluated in Calzavarini et al. (Reference Calzavarini, Volk, Bourgoin, Leveque, Pinton and Toschi2009) and Prakash et al. (Reference Prakash, Tagawa, Calzavarini, Mercado, Toschi, Lohse and Sun2012).

Table 1. Eulerian parameters for the two runs analysed in this article: Run I and Run II in the text. $N$ is the number of grid points in each spatial direction; $Re_{{\it\lambda}}$ is the Taylor-scale Reynolds number; ${\it\eta}$ is the Kolmogorov dissipative scale; ${\it\delta}x=\mathscr{L}/N$ is the grid spacing, with $\mathscr{L}=2{\rm\pi}$ denoting the physical size of the numerical domain; ${\it\tau}_{K}=\left({\it\nu}/{\it\varepsilon}\right)^{1/2}$ is the Kolmogorov dissipative time scale; ${\it\varepsilon}$ is the average rate of energy injection; ${\it\nu}$ is the kinematic viscosity; $t_{\mathit{dump}}$ is the time interval between two successive dumps along particle trajectories; ${\it\delta}t$ is the time step; $T_{L}=L/U_{0}$ is the eddy turnover time at the integral scale $L={\rm\pi}$ and $U_{0}$ is the typical large-scale velocity.

4.2 Comparison for acceleration variance and flatness against DNS

In figure 2(a) we show the comparison of the acceleration variance, $\langle a^{2}\rangle$ , to DNS data as a function of $St$ for different values of ${\it\beta}$ . Consider first heavy particles ( ${\it\beta}<1$ ). It is clear that the closure (3.1) captures the general trend and it becomes better and better for larger Stokes numbers. Similarly, this approximation must become exact as $St\rightarrow 0$ , but the DNS data set does not contain values of $St$ small enough to reach this limit. At intermediate Stokes numbers the Lagrangian closure described in § 3 does not match the DNS results. This mismatch for Stokes numbers between $0.1$ and $1$ is certainly due to preferential sampling, as already remarked by Bec et al. (Reference Bec, Biferale, Boffetta, Celani, Cencini, Lanotte, Musacchio and Toschi2006a ). Yet, the closure predicts a small Reynolds-number dependency (compare solid and dashed red lines for ${\it\beta}=0$ ) inherited by the dependency of the fluid tracers. Such a small variation is not detectable within the accuracy of our numerical data. Now consider light particles ( ${\it\beta}>1$ ). Figure 2(a) shows that also in this case the agreement between the DNS data and the closure scheme is good, except around $St\sim 10$ where the closure underestimates the DNS data. In this case inertial preferential sampling must be important. The reason is that light particles are drawn into vortex filaments where they experience high accelerations. Nevertheless, inertial preferential sampling must become irrelevant in the limit $St\rightarrow \infty$ as shown by the trend for very large $St$ in the same figure. It is interesting to remark that the closure scheme (solid line) becomes better and better the closer ${\it\beta}$ is to unity, and/or the smaller the Stokes number is, suggesting the possibility to develop a systematic perturbative expansion in the small parameter ${\it\epsilon}=St(1-{\it\beta})$ . Furthermore it is important to note that while the acceleration variance increases monotonically as both ${\it\beta}$ and $St$ increase, the flatness of the acceleration, defined as

(4.2) $$\begin{eqnarray}\displaystyle F_{St,{\it\beta}}=\langle |\boldsymbol{a}|^{4}\rangle /\langle \boldsymbol{a}^{2}\rangle ^{2}, & & \displaystyle\end{eqnarray}$$

has a non-monotonic dependency on $St$ . In figure 2(b) we show that light particles have a maximum in their flatness at $St\sim 0.5$ and heavy particles have minimum flatness for $St>1$ (except for the case of very heavy particles with ${\it\beta}=0$ ). Importantly, our closure approximation predicts these extrema qualitatively, indicating that the non-Gaussian tails observed numerically and experimentally (Volk et al. Reference Volk, Calzavarini, Verhille, Lohse, Mordant, Pinton and Toschi2008) in the acceleration probability distribution function of bubbles are not only due to inertial preferential sampling, which is neglected by the closure. Let us also note that the non-monotonicity shown by the model for $St\sim 0.1$ and ${\it\beta}=0$ , $0.25$ and $3$ are artefacts due to insufficient accuracy in the numerical evaluation of the integral in (3.7). In conclusion, remarkably, the closure approximation is in reasonable agreement with the DNS data, with the exception of those values of Stokes number where preferential sampling is important. This is an important result, supporting the idea that many key properties, including deviations from Gaussian statistics, of the acceleration statistics of inertial particles are due to the kinematic structure of the equation of motion together with the non-trivial Lagrangian properties of the fluid tracers evolution, and not only due to inertial preferential sampling.

Figure 3. Relative error ${\rm\Delta}a^{2}$ (4.3) between the DNS data and the closure scheme prediction for heavy particles (a) and light particles (b). Markers according to figure 1. Lines are drawn as a guide to the eye.

In order to quantify the accuracy of the closure described in § 3 we plot in figure 3 the relative error in the prediction of the acceleration variance for both light and heavy particles:

(4.3) $$\begin{eqnarray}\displaystyle {\rm\Delta}a^{2}=\left|1-\frac{\langle \boldsymbol{a}^{2}\rangle }{\langle \boldsymbol{a}^{2}\rangle _{DNS}}\right|, & & \displaystyle\end{eqnarray}$$

where with $\langle \boldsymbol{a}^{2}\rangle$ we refer to the expression (3.6) and with $\langle \boldsymbol{a}^{2}\rangle _{DNS}$ to the numerical value obtained from the DNS results. As one can see the approximation is never very bad, with a maximum discrepancy of the order of 30–40 % at those values where inertial preferential sampling is important, i.e. for $St\sim O(1)$ for heavy particles and for $St\sim O(10)$ for light particles.

5 Statistical Eulerian velocity model

Mathematical analysis and numerical studies of the particle dynamics become easier when the turbulent fluctuations of $\boldsymbol{u}(\boldsymbol{r},t)$ are approximated by a stochastic process. Following Gustavsson & Mehlig (Reference Gustavsson and Mehlig2016) we use a smooth, homogeneous and isotropic Gaussian random velocity field with root-mean-squared speed $u_{0}$ and typical length and time scales ${\it\eta}$ and ${\it\tau}$ . The model is characterized by a dimensionless number, the Kubo number, $Ku={\it\tau}/({\it\eta}/u_{0})$ , that measures the degree of persistence of flow structures in time. Very small Kubo numbers correspond to a rapidly fluctuating fluid velocity field. In this limit, the closure approximation described in § 3 is exact, inertial preferential sampling is negligible and the Lagrangian correlation functions of tracer particles are well approximated by the Eulerian correlation functions. In this limit it is also possible to perform a systematic perturbative expansion (Gustavsson & Mehlig Reference Gustavsson and Mehlig2016). In this paper we are interested in comparing the validity of the Lagrangian closure for the statistical model at $Ku\sim O(1)$ , where no analytical results can be obtained, and to further compare them with the DNS results shown in § 4. The motivation is the following. The statistical model has no ‘internal intermittency’, i.e. there is no Reynolds-number dependency on the acceleration statistics (there is not even the meaning of a Reynolds number). Nevertheless, once the Gaussian Eulerian velocity field is prescribed, we can calculate the acceleration probability density function of the fluid tracers. It turns out that this is not Gaussian and that it depends on the Kubo number, due to the effect of the quadratic advection term, $\boldsymbol{a}_{f}=\partial _{t}\boldsymbol{u}+Ku[\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\text{T}}]\boldsymbol{u}$ . As a result, we expect that many of the properties shown by the acceleration distribution of inertial particles evolved in real turbulent flows are shared by particles evolved in a Gaussian random flow. Finally, a comparison between DNS data and the statistical model will allow us to assess further the importance of internal intermittency.

5.1 Construction of the random-velocity field

For simplicity we discuss only the two-dimensional case. Generalization to three dimensions is straightforward. The velocity field is given in terms of the streamfunction: $\boldsymbol{u}(\boldsymbol{r},t)=\boldsymbol{{\rm\nabla}}{\it\psi}(\boldsymbol{r},t)\wedge \hat{\boldsymbol{e}}_{3}\,$ , which is defined as a superposition of Fourier modes with a Gaussian cutoff,

(5.1) $$\begin{eqnarray}\displaystyle {\it\psi}(\boldsymbol{x},t)=\frac{{\it\eta}^{2}u_{0}}{\sqrt{{\rm\pi}}L}\mathop{\sum }_{\boldsymbol{k}}a_{\boldsymbol{k}}(t)\text{e}^{\text{i}\boldsymbol{x}\boldsymbol{\cdot }\boldsymbol{k}-k^{2}{\it\eta}^{2}/4}. & & \displaystyle\end{eqnarray}$$

Here the system size $L$ is $10{\it\eta}$ , $k_{i}=2{\rm\pi}n_{i}/L$ and $n_{i}$ are integers with an upper cutoff $|n_{i}|\leqslant 2L/{\it\eta}$ because higher-order Fourier modes are negligible. The resulting spatial correlation function of ${\it\psi}$ is Gaussian, if $L\gg {\it\eta}$ we have: $\langle {\it\psi}(\boldsymbol{x},0){\it\psi}(\mathbf{0},0)\rangle =(u_{0}^{2}{\it\eta}^{2}/2)\,\text{e}^{-x^{2}/(2{\it\eta}^{2})}$ . The random coefficients $a_{\boldsymbol{k}}(t)$ in ${\it\psi}$ are drawn from random Gaussian distributions with zero means, smoothly correlated in time. To do that we used an Ornstein–Uhlenbeck process convolved with a Gaussian kernel of the form $w(t)\equiv \exp [-t^{2}/(2t_{0}^{2})]/(t_{0}\sqrt{2{\rm\pi}})\,$ to have a smooth correlation function also for the acceleration. The parameter $t_{0}$ must be small, $t_{0}\ll {\it\tau}$ , in order for the flow field to decorrelate at long times in a similar fashion as in fully developed turbulence. The Eulerian autocorrelation function of $\boldsymbol{u}$ is:

(5.2) $$\begin{eqnarray}\displaystyle \langle \boldsymbol{u}(\boldsymbol{x}_{0},t)\boldsymbol{\cdot }\boldsymbol{u}(\boldsymbol{x}_{0},0)\rangle =\frac{u_{0}^{2}}{2\sqrt{{\rm\pi}}}\text{e}^{-t^{2}/(4t_{0}^{2})}\left(\mathscr{F}\left[\frac{t_{0}}{{\it\tau}}+\frac{t}{2t_{0}}\right]\right.+\left.\mathscr{F}\left[\frac{t_{0}}{{\it\tau}}-\frac{t}{2t_{0}}\right]\right), & & \displaystyle\end{eqnarray}$$

where $\mathscr{F}(x)=\sqrt{{\rm\pi}}\exp (x^{2})\text{erfc}(x)$ . Note that the flow field is homogeneous in space and also in time and $\langle \boldsymbol{u}(\boldsymbol{x}_{0},t)\boldsymbol{\cdot }\boldsymbol{u}(\boldsymbol{x}_{0},0)\rangle$ is a function of $|t|$ only. For $|t|\ll t_{0}$ this correlation function is Gaussian and for $|t|\gg t_{0}$ the correlation function is exponential: $\langle \boldsymbol{u}(\boldsymbol{x}_{0},t)\boldsymbol{\cdot }\boldsymbol{u}(\boldsymbol{x}_{0},0)\rangle \sim \text{e}^{-|t|/{\it\tau}}$ .

5.2 Results for the random-velocity model

We consider first the acceleration variance. Simulation results for the statistical model are compared with the Lagrangian closure (3.6)–(3.7) in figure 4 by changing both $St,{\it\beta}$ for $Ku=5$ (panel a). We observe a good agreement between the Lagrangian closure and the numerical simulations, comparable to what was observed for the DNS data. In figure 4(b) we show the results for the flatness. It is important to stress two facts. Also here, both heavy and light particles depart from the corresponding fluid value with a qualitative trend similar to that observed for the DNS case in the previous section. Also, for the random velocity field, the Lagrangian closure works qualitatively well. The departure from the numerical data is a signature of the corresponding importance of preferential sampling at those Stokes numbers. Let us notice nevertheless an important difference with respect to the DNS data. Here the absolute values of the flatness are much smaller, due to the absence of internal intermittency. In the stochastic signal the acceleration of the fluid is non-Gaussian only because of kinematic effects. In real flows, the acceleration is more intense and more fluctuating because of the vortex stretching mechanism and of the turbulent energy cascade.

Figure 4. Same as in figure 2 but for the statistical model given by (5.1) at $Ku=5$ .