Turbulence modulation in heavy-loaded suspensions of tiny particles

TL;DR

This paper investigates how tiny particles in heavy-loaded suspensions influence turbulence, revealing a new energy spectrum scaling law and demonstrating the impact on energy transfer and particle clustering through a robust regularised model and numerical simulations.

Contribution

It introduces a physically-based regularisation approach to describe small-scale turbulence modulation by tiny particles, capturing the $k^{-4}$ energy spectrum scaling and energy exchange mechanisms.

Findings

Particles modify the energy spectrum to follow a $k^{-4}$ scaling law.

The regularised model accurately captures small-scale turbulence physics.

Particle clustering at small scales is reduced, affecting collision dynamics.

Abstract

The features of turbulence modulation produced by a heavy loaded suspension of small solid particles or liquid droplets are discussed by using a physically-based regularisation of particle-fluid interactions. The approach allows a robust description of the small scale properties of the system exploiting the convergence of the statistics with respect to the regularisation parameter. It is shown that sub-Kolmogorov particles/droplets modify the energy spectrum leading to a scaling law, $E(k)\propto k^{-4}$, that emerges at small scales where the particle forcing balances the viscous dissipation. This regime is confirmed by Direct Numerical Simulation data of a particle-laden statistically steady homogeneous shear flow, demonstrating the ability of the regularised model to capture the relevant small-scale physics. The energy budget in spectral space, extended to account for the inter-phase momentum exchange, highlights how the particle provide an energy sink in the production range that turns into a source at small scales. Overall, the dissipative fluid-particle interaction is found to stall the energy cascade processes typical of Newtonian turbulent flows. In terms of particle statistics, clustering at small scale is depleted, with potential consequences for collision models.

Equations

$$\begin{array}[]{l}\displaystyle\nabla\cdot{\bf u}=0\\ \\ \displaystyle\rho\left(\frac{\partial{\bf u}}{\partial t}+{\bf u}\cdot\nabla{\bf u}\right)=-\nabla p+\mu\nabla^{2}{\bf u}+{\bf F}\end{array}$$

$$\mathbf{F}(\mathbf{x},t)=-\sum _{p=1}^{N_p}{\mathbf{D}}_p(t)\delta \left[\mathbf{x}-{\mathbf{x}}_p(t)\right]$$

$${\mathbf{F}}_R(\mathbf{x},t)=-\sum _{p=1}^{N_p}{\mathbf{D}}_p(t-{\epsilon }_R)g\left[\mathbf{x}-{\mathbf{x}}_p(t-{\epsilon }_R),{\epsilon }_R\right],$$

$$\begin{array}[]{l}\displaystyle\frac{d{\bf x}_{p}}{d\,t}={\bf v}_{p}\\ \\ \displaystyle\frac{d{\bf v}_{p}}{d\,t}=\frac{{\bf D}_{p}}{m_{p}}=\frac{{\bf u}\big{|}_{{\bf x}_{p}}-{\bf v}_{p}}{\tau_{p}}\,,\end{array}$$

$$\frac{\partial E(k)}{\partial t}=T(k)+P(k)-D(k)+\Psi (k)$$

$${\int }_0^{\infty }\Psi (k)dk=-{ϵ}_e<0$$

$$\begin{array}[]{l}\left.\begin{array}[]{l}\displaystyle\nabla\cdot{\bf u}=0\\ \displaystyle\frac{\partial{\bf u}}{\partial t}+{\bf u}\cdot\nabla{\bf u}=-\frac{1}{\rho_{f}}\nabla{\rm p}+\nu\nabla^{2}{\bf u}\end{array}\right\}\qquad{\bf x}\in{\cal D}\backslash\Omega(t)\\ {}\\ \displaystyle{\bf u}\lvert_{\partial\Omega_{p}(t)}={\bf v}_{p}({\bf x})\lvert_{\partial\Omega_{p}(t)}\qquad\qquad p=1,\ldots,N_{p}\\ {}\\ \displaystyle{\bf u}({\bf x},0)={\bf u}_{0}({\bf x})\qquad\qquad{\bf x}\in{\cal D}\backslash\Omega(0)\ .\end{array}$$

$$\begin{array}[]{l}\begin{array}[]{l}\displaystyle\nabla\cdot{\bf w}=0\\ \displaystyle\frac{\partial{\bf w}}{\partial t}+{\bf F}=-\frac{1}{\rho_{f}}\nabla\pi+\nu\nabla^{2}{\bf w}\\ \end{array}\\ \displaystyle{\bf w}({\bf x},0)={{\bf u}}_{0}({\bf x})\,,\end{array}$$

$${\bf F}=\left\{\begin{array}[]{ll}{\bf u}\cdot\nabla{\bf u}&\qquad\mbox{for}\,\,{\bf x}\in{\cal D}\backslash\Omega(t)\\ \\ {\bf v}_{p}\cdot\nabla{\bf v}_{p}&\qquad\mbox{for}\,\,{\bf x}\in\Omega(t)\end{array}\right.$$

$$\begin{array}[]{l}\left.\begin{array}[]{l}\displaystyle\nabla\cdot{\bf v}=0\\ \displaystyle\frac{\partial{\bf v}}{\partial t}=-\frac{1}{\rho_{f}}\nabla{\rm q}+\nu\nabla^{2}{{\bf v}}\end{array}\right\}\qquad{\bf x}\in{\cal D}\backslash\Omega(t)\\ {}\\ \displaystyle{\bf v}\lvert_{\partial\Omega_{p}(t)}={\bf v}_{p}({\bf x})\lvert_{\partial\Omega_{p}(t)}-{\bf w}\lvert_{\partial\Omega_{p}(t)}\qquad p=1,\ldots N_{p}\\ {}\\ \displaystyle{\bf v}({\bf x},0)=0\qquad{\bf x}\in{\cal D}\backslash\Omega(0)\,.\end{array}$$

$$v_i(\mathbf{x},t)=-\sum _p{\int }_0^t{D}_j^p(\tau )G_{ij}(\mathbf{x},{\mathbf{x}}_p,t,\tau )d\tau ,$$

$$\begin{array}[]{l}\displaystyle\frac{\partial{\bf v}}{\partial t}-\nu\nabla^{2}{\bf v}+\frac{1}{\rho_{f}}\nabla{\rm q}=-\frac{1}{\rho_{f}}\sum_{p}{\bf D}_{p}(t)\,\delta\left[{\bf x}-{\bf x}_{p}(t)\right]\\ \\ \displaystyle{\bf v}({\bf x},0)=0\,.\end{array}$$

$$\mbox{\boldmath$\zeta$}({\bf x},t)=\frac{1}{\rho_{f}}\int_{0}^{t}{\bf D}_{p}(\tau)\times\nabla g\left[{\bf x}-{\bf x}_{p}(\tau),t-\tau\right]d\tau\,,$$

$$\mbox{\boldmath$\zeta$}_{R}({\bf x},t)=\frac{1}{\rho_{f}}\int_{0}^{t-\epsilon_{R}}{\bf D}_{p}(\tau)\times\nabla g\left[{\bf x}-{\bf x}_{p}(\tau),t-\tau\right]d\tau$$

$$\displaystyle\frac{\partial\mbox{\boldmath$\zeta$}_{R}}{\partial t}$$

$$\frac{\partial {\mathbf{v}}_R}{\partial t}$$

$$\displaystyle\frac{\partial{\bf v}_{\mbox{\boldmath$\zeta$}_{R}}}{\partial t}$$

$$\frac{\partial \mathbf{u}}{\partial t}$$

$${\mathbf{D}}_p(t)=6\pi \mu a_p\left[\tilde{\mathbf{u}}({\mathbf{x}}_p,t)-{\mathbf{v}}_p(t)\right]$$

$${\mathbf{v}}_R(\mathbf{x},t_{n+1})$$

$$\frac{\partial \mathbf{u}}{\partial \tau }+\mathbf{u}\cdot \nabla \mathbf{u}=-\frac{1}{\rho }\nabla p+\nu {\nabla }^2\mathbf{u}-S{u}_2{\mathbf{e}}_1.$$

Full text

Particle laden turbulent flows are central in many physical and technological contexts. In astrophysics [1, 2] the turbulence is known to influence the aggregation of dust particles in protoplanetary (accretion) disks, see [3, 4, 5] and reference therein. Similarly, in warm clouds, the turbulence controls the growth by condensation of small droplets [6], and ultimately speeds-up rain formation [7, 8]. In the combustion of liquid fuels [9, 10], the turbulence determines the effectiveness of atomization, evaporation and mixing [11]. All these examples show that turbulence strongly interacts with the transported phase. Less understood is the reciprocal effect expected on the basis of the action-reaction principle by which the transported phase alters the turbulence. An extreme example of this reciprocal effect arises in the environmental context, where small active organisms such as plankton [12] or bacteria [13] induce small-scale chaotic flows which affects the chemical and the biological activity. Significant alteration of the turbulent flow is also found in bubbly grid-generated flows, [14]. In general, significant back reaction effects are expected in all the other contexts mentioned above. Concerning in particular technological applications, in a typical diesel engine, see e.g. [15], the mass of fluid injected per cycle per cylinder in the form of small droplets is about $3\;\times 10^{-4}\,{\rm kg}$. Considering a four stroke, $2.5$ litre engine with $4$ cylinders, back of the envelope calculations immediately give a mass loading of about $\phi\simeq 0.4$ and a volume fraction of the order of $\phi_{v}\simeq 6\,\times 10^{-3}$. In modern common-rail injection systems the diameter of the droplets is about $d_{p}\simeq 0.1-10\mu m$ whilst the Kolmogorov scale in a combustion chamber can be estimated on the order of $\eta\simeq 30\mu m$. According to the accepted classification, see  [16, 17], the suspension must then be considered dilute (no direct interaction among droplets) even though the inter-phase momentum coupling is particularly significant.

Among the different regimes of a particle laden flow [18], the present Communication addresses conditions like those mentioned above where i) the dimensions of the single suspended particle are much smaller than the relevant macroscopic scale of the turbulent flow, $d_{p}/\eta\ll 1$; ii) the particles are extremely diluted with negligible direct particle-particle interaction, i.e. the volume fraction is small; iii) the mass loading of the suspension (particle to fluid mass ratio) is significant, $\phi=m_{p}/m_{f}={\cal O}(1)$, implying that a considerable particle-induced force is exerted on the flow. In these conditions, beside turbulence-induced particle clustering already observed at small mass loading [19, 20, 21, 22, 23, 25, 26], new phenomena associated to turbulence modulation are expected, defining a still poorly understood realm of multiphase turbulence. In particular, the standard Kolmogorov-like paradigm [27], which assumes that the turbulence is forced at large scales and eventually dissipated at small scales with a universal direct energy cascade [28] emerging in the inertial range, is expected to fail.

In the new conditions the particle population forces the fluid across the entire range of available scales, posing several new questions concerning the structure and the dynamics of turbulence under significant back-reaction effects. The first class of questions is methodological: how can the effect of many sub-Kolmogorov particles be modelled in a physically consistent manner in Direct Numerical Simulations (DNS)? Is a numerical simulation which truly couples the discrete, point-like phase with a continuum fluid feasible with the present state-of-the-art numerical tools? Can the coupling be made realistic yet affordable from the computational point of view? Are the singularities arising from the coupling amenable of rigorous treatment? As will be shown, answers to these methodological issues can be found in the context of a newly designed inter-phase momentum coupling strategy, the Exact Regularized Point Particle approach (ERPP) [29]. The second family of questions, is more physical: what are the effects of the back-reaction on the turbulence dynamics? How the disperse phase affects the energy cascade processes and, in turns, the energy spectrum? What is the resulting effect of the coupling on the particle population? Can we trust the numerical predictions, particularly at small scales, where most of the particle-fluid interaction is expected to occur?

This Communication provides an answer to all these questions, discussing the results of new simulations based on the ERPP approach that are free of the bias that hampers other available techniques aimed at realising the particle-fluid interaction. Among others, the crucial advancements over the present state-of-the-art concern: a) a physically-based, grid-independent regularisation of the singular response of point-like particles; b) the possibility to take a weak limit for the statistics with the regularisation parameter approaching zero; c) the ability to exactly remove from the field the unphysical self-induction velocity of each single particle in the calculation of the hydrodynamic force; d) the recovery of the exact momentum balance in the force coupling of each particle with the fluid; e) the convergence of the coupling scheme also when a fixed number of particles, independent of grid size, is considered.

In order to address these issues in the cleanest form, the flow should be as simple as possible. Traditionally homogeneous and isotropic turbulence is the elective choice. However it requires an external forcing acting at large scales to provide the energy dissipated by viscosity. Although this is not an issue for classical Newtonian turbulence, the external forcing introduces undesired features in the context of particle laden flows in presence of backreaction. The reason is that, as shown below, the particle forms long clusters spanning the entire range of scales, up to the integral scale. The external forcing interferes with the large scales of the clusters and their backreaction on the fluid, thereby introducing dynamical artefacts. A flow able to self-sustain the turbulence with no artificial external forcing which still retains a substantial simplicity, e.g. statistical spatial homogeneity and stationarity, is the homogeneous shear flow, where a linear average shear is enforced on turbulence fluctuations. This flow, described in detail in the Supplemental Material (SM, [30]), will be exploited below to discuss generic features of particles laden flows under strong loading.

When $d_{p}\ll\eta$, the carrier flow is described by the incompressible Navier-Stokes equations

[TABLE]

where

[TABLE]

is the (singular) field representing the back-reaction of the point-like particles on the flow. In equation (2), $N_{p}$ denotes the number of particles, ${\bf D}_{p}$ the hydrodynamic force acting on the p-th particle and the Dirac delta function localises the force at the particle position ${\bf x}_{p}(t)$. Clearly equations (1-2) need to be regularized to be amenable to numerical treatment.

In the classical Particle in Cell approach, see e.g. [31], the singularity is removed averaging the feedback on the computational cell, giving rice to several drawbacks, see e.g. [32, 18, 33]. Typically, convergence can be achieved only at constant number of particles per computational cell, implying that the number of particles should increase (at constant mass loading) as the grid size is reduced. Additionally, the particles are affected by their own self-induced disturbance, which introduces errors in the hydrodynamics force. This source of error gets more and more pronounced as the number of particles per cell is reduced, as always happens under grid refinement. These drawbacks do not affect the ERPP method where the Dirac delta function is regularized in a physically consistent manner. The disturbance due to each point-like particle is evaluated in a closed analytical form exploiting the exact solution of a local unsteady Stokes problem and the viscosity of the fluid naturally takes care of regularising the fluid response to the particle forcing. In turbulence, when $d_{p}\ll\eta$, is natural to set the regularisation length on the order of the Kolmogorov length-scales $\eta$ or below. The singular forcing (2) is effectively replaced by its (exact) regularized counterpart,

[TABLE]

where the Gaussian function $g$ consistently emerges from the small scale diffusion of the particle disturbance field described by the unsteady Stokes operator [29]. The spatial cut-off scale $\sigma_{R}=\sqrt{2\nu\varepsilon_{R}}$ is directly related to the diffusion time-scale $\varepsilon_{R}$ which represent the typical time needed by the singular vorticity produced by the particle at time $t-\varepsilon_{R}$ to spread over the resolved length scale at time $t$, see [29] and Supplemental Material [30].

The dispersed particles follow Newton’s equations,

[TABLE]

where ${\bf x}_{p}$ and ${\bf v}_{p}$, $p=1,\dots,N_{p}$, are the particle positions and velocities, respectively, $m_{p}$ is the particle mass and, in the conditions considered here, ${\bf D}_{p}$ reduces to the Stokes drag [34, 35] proportional to the fluid-particle relative velocity with {\bf u}\big{|}_{{\bf x}_{p}} the fluid velocity at the particle position. In the jargon of particle laden flows, such relative velocity is sometimes called the slip velocity. Since the particle modifies the fluid velocity, care should be taken not to contaminate {\bf u}\big{|}_{{\bf x}_{p}} with the particle self-disturbance. Otherwise, at decreasing the grid size, the spurious contribution would dominate the overall particle-fluid interaction. The Stokes number, $St_{\eta}=\tau_{p}/\tau_{\eta}$, where $\tau_{\eta}$ is the Kolmogorov time scale of the turbulence and $\tau_{p}=(\rho_{p}/\rho)d_{p}^{2}/(18\nu)$ is the Stokes relaxation time, is a central control parameter which, e.g., determines the intensity of particle clustering, that is the trend to segregate [19, 23, 7, 25, 32, 6, 20] in long, tiny structures.

Figure 1 shows a slice of an instantaneous configuration of particle distribution and feedback force field in a turbulent homogeneous shear flow. The turbulence, at $Re_{\lambda}=\lambda u_{{\rm rms}}/\nu=80$, with $\lambda=u_{{\rm rms}}\sqrt{\nu/\epsilon}$ and $u_{{\rm rms}}=\sqrt{\langle(u-{\bar{u}})^{2}\rangle}$, $St_{\eta}=1$ and $\phi=0.4$, is sustained by a constant mean shear $S=dU_{x}/dy$, see Supplemental Material [30] for details. The energy is extracted from the mean flow by the Reynolds shear stresses $-\langle u\,v\rangle$ which force the turbulent fluctuations at scales larger than the shear scale $L_{S}=\sqrt{\epsilon/S^{3}}$ [36]. Typical of unitary Stokes number flows, the disperse phase forms elongated clusters, apparent in the plot. They are oriented by the mean flow which imprints on them a strong anisotropy. The clusters span a range of scales from their width, of the order of the dissipative scale, up to their length, comparable with the integral scale of the flow [25, 32]. The force feedback ${\bf F}_{R}$ is strongly correlated with the clusters and affect the same range of scales. This kind of distributed, effective field differs substantially from the classical Kolmogorov scenario where the forcing is designed to prevent the flow from dissipating, it is confined to the large scales to avoid contamination of the cascade and is assumed to be statistically independent of the flow.

It is instrumental to look at the flow in spectral space where, adopting index notation, the interphase momentum coupling is described by the Fourier transform $\cal F$ of the correlation $\Psi_{ij}({\bf k})={\cal F}\langle F_{R,i}({\bf x})\,u_{j}({\bf x}+{\bf r})\rangle$ between the back-reaction and the fluid velocity. The quantity $\Psi(k)=\int_{\Omega}\Psi_{ii}({\bf k})\,k^{2}d\Omega$, where the integral is taken over the solid angle $\Omega$ in wavenumber space and $k^{2}={\bf k}\cdot{\bf k}$, forces the equation for the turbulence spectrum $E(k)$ according to

[TABLE]

where $E(k)=\int_{\Omega}E_{ii}({\bf k})\,k^{2}d\Omega$ and $E_{ij}({\bf k})={\cal F}\langle u_{i}({\bf x})\,u_{j}({\bf x}+{\bf r})\rangle$. Equation (5) is the extension to particle laden flows of the classical equation for the spectral balance of turbulent kinetic energy, sometimes called the Kolmogorov-Onsager-von Weizsäcker-Heisenberg equation [37, 28]. In equation (5) the energy transfer term $T(k)$ is defined as $T(k)=\int_{\Omega}\imath k_{j}T_{j}({\bf k})\,k^{2}d\Omega$ where the Fourier transform of the triple correlation function is $T_{j}({\bf k})={\cal F}\langle u_{i}({\bf x})u_{i}({\bf x}+{\bf r})u_{j}({\bf x})-u_{i}({\bf x}+{\bf r})u_{i}({\bf x})u_{j}({\bf x}+{\bf r})\rangle$. The non-linear triadic interactions among different Fourier modes conserves energy, $\int_{0}^{\infty}T(k)\,dk=0$, and ultimately originate the energy cascade. $P(k)=-SE_{uv}(k)$ is the production of turbulent kinetic energy at wavenumber $k$ where $E_{uv}=E_{12}(k)$ is the energy cospectrum and $D(k)=2\nu k^{2}E(k)$ the dissipation spectrum. Note that once integrated across the entire range of wavenumbers the energy cospectrum returns the Reynolds shear stresses $-\langle u\,v\rangle=\int_{0}^{\infty}E_{uv}(k)\,dk$, and the dissipation spectrum gives the viscous dissipation $\epsilon=\int_{0}^{\infty}D(k)\,dk$. In statistically steady conditions the time derivative of the energy spectrum vanishes.

Concerning eq. (5), one of the simulative issues with particle laden flows in the two coupling regime, is the sensitivity of small scale observables to the numerical implementation of the particle feedback. The approach here proposed allows for obtaining a clean asymptotic also for small scale observables. This is achieved in the limit $\sigma_{R}\to 0$, where the limit is to be understood in the weak sense, i.e. first the statistics is acquired as a function of the regularisation parameter and only after the limit is taken on the averages. This process is illustrated in figure 2 where turbulent kinetic energy spectra are shown for the same particle population and two different Reynolds number at decreasing $\sigma_{R}/\eta$. Apparently the data nicely collapse and a well defined energy distribution emerges at decreasing $\sigma_{R}$. This is expected at large scales which soon become independent of the regularisation parameter. A new feature emerges at small scales (large wavenumber) where a well definite scaling range eventually appears at $k\eta\simeq 1$. The right panel shows the compensated plot, $k^{4}E(k)$ vs $k\eta$. About one decade of $k^{-4}$ scaling is detected for the smallest $\sigma_{R}/\eta$ we have considered. The scaling range approximately extend from about $k_{\lambda}\eta\simeq 0.45$, which is order of the Taylor micro-scale where the dissipation spectrum peaks, to the cut-off $k\sigma_{R}\simeq 1$ corresponding to $k\eta\simeq 4$.

We may note that data in absence of particle feedback show a completely different trend, consistently with the behaviour expected in the dissipation range. This result shows that the regularisation procedure we have put forward can be used to obtain physically significant and numerically convergent information on the small scale statistics of the system. Indeed, by reducing $\sigma_{R}$ at given turbulence intensity, we can approach any given small scale in the system. This is important in view of taking into account interactions between particles, such as collisions, lubrication effects, short range attraction or repulsion between particles, e.g. Van der Walls forces, which arise at the inner length scale $d_{p}$ of the particles.

For comparison, the right panel of figure 2 reports the compensated spectra obtained with the PIC approach operated in the same conditions, namely $Re_{\lambda}=55$ and $\sigma_{R}/\eta=0.5$. Mass loading $\phi=0.4$ and Stokes number $St_{\eta}=1$ fix the number of particles $N_{p}=595520$, corresponding to few particles per cell, namely $N_{p}/N_{c}\simeq 0.04$ where $N_{c}$ is the number of computational cells. The PIC approach is reasonably able to describe the behaviour of the compensated spectrum at $k\eta\simeq 1$ where a glimpse of a short plateau seems to appear. However, at smaller scales, the trend reveals a clear departure from the $k^{-4}$ scaling law. The reason is that the high wave number modes are badly behaved due to the non-smooth and grid dependent numerical feedback field, see e.g. [32]. This hampers reaching progressively smaller and smaller scales. The behaviour gets worser and worser when finer grids are used (data not shown).

The spectral budget, eq. (5), is shown in figure 3. The main panel focuses on the range of wave-numbers where the $k^{-4}$-scaling is observed (see the inset for a global view). The production $P(k)$ and the transfer term $T(k)$ vanish where $k\eta\simeq 1$, showing that the dominant balance is between the inter-phase coupling $\Psi(k)$ – the only energy source present at those scales in absence of the energy transfer – and the viscous dissipation $D(k)$. The back-reaction has overwhelmed the inertial transfer and stalled the energy cascade, right panel with the comparison of the energy transfer with, $T(k)$, and without, $T(k)_{\phi=0}$, coupling. The reduced transfer is replaced by the energy injected by the particles which, in turn, drain from the large scales the energy $P(k)$ extracts from the mean flow. As a consequence, the energy feeding the cascade is reduced by the amount drained by the disperse phase. The balance between energy intercepted by the particles at large scales and the energy released at small scales is negative,

[TABLE]

implying a dissipative effect of the particles. Considering the overall budget, including fluid and particles, $-S\langle u\,v\rangle=\epsilon+\epsilon_{e}$, the energy produced by the Reynolds stresses is turned into the sum of viscous dissipation and the extra-dissipation due to the particles, $\epsilon_{e}$. In other words, the disperse phase provides an alternative dissipation channel.

The data just discussed show that the $k^{-4}$ scaling range corresponds to the region where $\Psi(k)\simeq D(k)$. Note that in a periodic box any term in eq. (5), defined as the Fourier transform of the relevant correlation, can be replaced by the average product of the corresponding Fourier coefficients, e.g. $\Psi(k)=\langle{\hat{F}}_{R,i}(k){\hat{u}}_{i}^{*}(k)\rangle$. In order to get a deeper insight into the origin of the new scaling law, it is useful to consider the spectrum of the particle back-reaction field $F(k)=\langle\hat{F}_{R,i}(k)\hat{F}^{*}_{R,i}(k)\rangle$. Figure 4 shows $F(k)$ for the case at $Re_{\lambda}=55$ and $\sigma_{R}/\eta=0.25$, which is the case with the largest separation between Kolmogorov and regularisation scale we have considered. In the range of wavenumbers centred at $k\eta\simeq 1$ which are not yet affected by the regularisation, i.e. $k\sigma_{R}<1$, $F(k)\simeq{\hat{F}}^{2}_{0}$ is roughly constant. This result is somehow expected since the field $F_{R,i}({\bf x},t)$ is the superposition of Gaussians with variances still significantly smaller than the considered scales, see eq. (3). The Fourier transform reads $\displaystyle\hat{F}_{R,i}=-\sum_{p=1}^{N_{p}}D_{p,i}(t-\varepsilon_{R})e^{-\frac{1}{2}k^{2}\sigma_{R}^{2}}\,e^{-\imath k_{j}x_{p,j}(t-\varepsilon_{R})}$ which, apart from the phase, is proportional to $e^{-1/2k^{2}\sigma_{R}^{2}}$, hence almost constant for $k\sigma_{R}<1$. The inter-phase momentum coupling $\Psi(k)$ is also reported in the figure in comparison with the estimate $\sqrt{F(k)\,E(k)}$ (squares). The data show that, where $F(k)\simeq{\hat{F}}^{2}_{0}$, $\Psi(k)$ closely matches the curve $\sqrt{F(k)\,E(k)}$. It follows that $\Psi(k)\sim\sqrt{F(k)\,E(k)}\sim{\hat{F}}_{0}\sqrt{E(k)}$. Then, given the observed $k^{-4}$ scaling for the spectrum, we infer $\Psi(k)\propto k^{-2}$, as confirmed by the collapse of the data represented by circles ($\Psi$), squares ($\sqrt{F(k)\,E(k)}$) and solid line ($k^{-2}$). In other words, at these scales, the Fourier transform of velocity and backreaction are found to be uncorrelated. This suggests that a purely dimensional argument can be put forward: neglecting force-velocity correlations in the Fourier modes at small scales, assuming $\Psi(k)\sim{\hat{F}}_{0}\,\hat{u}$, and introducing the ansatz $\hat{u}\propto k^{\alpha/2}$, the balance of backreaction $\Psi(k)$ and dissipation $D(k)=2\nu k^{2}E(k)\sim k^{2}{\hat{u}}^{\alpha}$ leads to the observed scaling law $E(k)\propto k^{-4}$.

From previous studies in the one-way-coupling regime it is well known that clustering peaks at $St_{\eta}={\cal O}(1)$ [38, 19]. Clustering is also observed in the two way coupling regime. It is however substantially reduced by the back reaction, as measured by the radial distribution function (RDF, see 111The radial distribution function $g_{00}(r)$ is the density of particle pairs in a ball ${\cal B}_{r}$ of radius $r$ normalised with the density pairs $n_{0}=0.5N_{p}(N_{p}-1)/V_{0}$ in the whole domain $V_{0}$, namely $g_{00}(r)=1/(4\pi r^{2}n_{0})dN_{r}/dr$, where $N_{r}$ is the number of pairs in the ball ${\cal B}_{r}$. The small scales divergence of the radial distribution function corresponds to the occurrence of small scale clustering. In fact, whenever a scaling law $g_{00}(r)\propto r^{-\alpha}$ with positive $\alpha$ occurs, the scaling exponent $\alpha$ measures the correlation dimension ${\cal D}_{2}=3-\alpha$ of the multi-fractal measure associated with the particle density [39].) of the particles shown in figure 5. Clustering increases the overall probability that particles could collide. Beside clustering, the collision frequency is determined by the mean relative velocity of close particles - a further crucial small scale property of the system that needs accurate modelling. Technically, the relevant statistical quantity is the average longitudinal velocity difference between two particles $Q_{00}=\langle\delta v_{\parallel}(r)|\delta v_{\parallel}(r)<0\rangle$ where the average is conditioned to negative relative velocity $\delta v_{\parallel}$ 222The collision rate, i.e. the number of collision per unit time and volume is given by $\Gamma=2\pi\sigma^{2}g_{00}(r=\sigma)Q_{00}(r=\sigma)$ where $\sigma=d_{p}$ is the collision radius, $g_{00}(r=\sigma)$ is the RDF evaluated at collision and $Q_{00}(r=\sigma)$ is mean relative velocity of the colliding pair, see e.g. [38]., right panel of figure 5. The collision probability is proportional to the product $g_{00}\times Q_{00}$, [38] evaluated at contact ($r=d_{p})$. This object is reported in the inset of the right panel of the figure as a function of separation. The present data show that, in the relevant range of scales below $\eta$, the two-way coupling may deplete the collision frequency since the decrease of the clustering intensity prevails on the slight increase of the relative velocity.

In conclusion the present Communication highlights new features of turbulence in highly loaded suspensions of tiny, heavy particles. The particles are found to drain energy from the carrier flow at the large scales and release it back at the small scales. It follows that, in this kind of multiphase flows, turbulent fluctuations are unusually forced in the dissipative range. The back-reaction stalls the energy cascade and enforces a newly observed $E(k)\propto k^{-4}$ scaling law for the energy spectrum at scales order of $\eta$, where the particle-injected energy is immediately dissipated by viscosity. Noteworthy, small scale clustering is depleted by the particle-fluid interaction while the relative particle velocity is slightly modified. Consequently, the collision probability turns out to be reduced. In more general terms, it has been shown that the coupling strategy described in the Communication provides a viable technique to robustly evaluate small scale statistics in highly loaded particle laden flows. The approach, relying on a physical regularisation of the singular force feedback, provides convergent result with respect to the regularisation parameter allowing a safe evaluation of central observables for heavy loaded dilute suspensions. The approach can be easily extended to turbulent flow laden with micro-bubbles and to wall bounded flows.

The research received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement no. [339446]. Support from PRACE, projects FP7 RI-283493 and grant no. 2014112647, is acknowledged.

Supplemental Material

0.1 The Exact Regularized Point Particle method

This section summarizes the physical model used to couple carrier fluid and disperse phase for the specific case of periodic boundary conditions considered in the simulation of the homogenous turbulent shear flow. The reader can refer to [29] for additional details and more general conditions.

The carrier fluid fills the domain ${\cal D}\backslash\Omega$ where $\cal D$ is the periodic flow domain and $\Omega(t)=\cup_{p}\Omega_{p}(t)$ denotes the region occupied by the $N_{p}$ rigid particles, with $\Omega_{p}(t)$ the – small but still finite – domain occupied by the $p$th particle. The motion of the carrier fluid is described by the incompressible Navier-Stokes equations with the no-slip condition at the particle boundaries and periodic boundary conditions on $\partial{\cal D}$,

[TABLE]

In equations (6), ${\bf u}_{0}({\bf x})$ is the velocity field at time $t=0$, $\rho_{f}$ denotes the fluid density, $\nu$ is the kinematic viscosity, and ${\bf v}_{p}({\bf x})$ the velocity of the particle boundary. In presence of a number of small particles the idea is to relocate the boundary condition at the particle surface on a properly defined correction flow field for which, in the limit of small particles, an analytical solution can be provided. The carrier fluid velocity ${\bf u}$ is decomposed into two parts, ${\bf u}({\bf x},t)={\bf w}+{\bf v}$ where the periodic (background) field ${\bf w}({\bf x},t)$ is assumed to satisfy the equations

[TABLE]

where ${\bf x}\in{\cal D}$ and

[TABLE]

is a field reproducing the convective term of the Navier-Stokes equation in ${\cal D}\backslash\Omega$. For the time being, in (7) the convective term is considered as prescribed and the no-slip condition at the particle surface has been removed. In fact, the no-slip boundary condition at the particle surface is recovered when considering the particle perturbation field ${\bf v}({\bf x},t)$ which satisfies the linear unsteady Stokes problem (the complete non-linear term has been retained in the equation for ${\bf w}$)

[TABLE]

The boundary integral representation of the solution to the unsteady Stokes equations (9) can be expressed in terms of multipoles. In the limit of small particles the far field reduces to

[TABLE]

where G_{ij}({\bf x},\mbox{\boldmath\xi},t,\tau) is the unsteady Stokeslet, i.e. the fluid velocity ($i$th direction) at position ${\bf x}$ and time $t$ due to the singular forcing \delta({\bf x}-\mbox{\boldmath\xi})\delta(t-\tau) ($j$th direction) applied at point $\xi$ and at time $\tau$ and $D^{p}_{j}(\tau)$ are the Cartesian components of the hydrodynamic force on the particle. The partial differential equation whose solution is given by (10) reads

[TABLE]

A regularized solution of equation (11) can be achieved by reasoning in terms of the vorticity field \mbox{\boldmath\zeta}=\nabla\times{\bf v} which obeys a (vector) diffusion equation. The solution is

[TABLE]

where g({\bf x}-\mbox{\boldmath\xi},t-\tau) is a Gaussian function with time dependent variance $\sigma(t-\tau)=\sqrt{2\nu(t-\tau)}$. The (still) singular field $\zeta$ is regularized using a temporal cut-off $\epsilon_{R}$ leading to its splitting into a regular and a singular component \mbox{\boldmath\zeta}({\bf x},t)=\mbox{\boldmath\zeta}_{R}({\bf x},t;\epsilon_{R})+\mbox{\boldmath\zeta}_{S}({\bf x},t;\epsilon_{R}) where

[TABLE]

is smooth with smallest spatial scale given by $\sigma_{R}=\sigma(\epsilon_{R})=\sqrt{2\nu\epsilon_{R}}$. It obeys a forced diffusion equation where the forcing is applied at the slightly earlier time $t-\epsilon_{R}$,

[TABLE]

with \mbox{\boldmath\zeta}_{R}({\bf x},0)=0. The associated velocity field obeys the forced unsteady Stokes equation

[TABLE]

for the solenoidal field ${\bf v}_{R}$ that can be split in terms of a pseudo-velocity,

[TABLE]

governed be the unsteady diffusion operator plus a gradient correction required to enforce solenoidality. The highly localized singular contribution ${\bf v}_{S}$, which cannot be represented on a discrete grid, is successively reintroduced in the field as soon as it diffuses sufficiently to reach the smallest physically relevant scales. The regularized (solenoidal) fluid velocity in presence of the particles is ${\bf u}_{R}={\bf w}+{\bf v}_{R}$ and is governed by the equation

[TABLE]

0.2 The Removal of Particle Self-interaction in the Evaluation of the

Hydrodynamic Force

The hydrodynamic force acting on a small particle of diameter $d_{p}$ and density $\rho_{p}\gg\rho_{f}$ reduces to the Stokes drag [34, 35],

[TABLE]

The velocity $\tilde{\bf u}({\bf x}_{p},t)$ is the fluid velocity, at the particle position, in absence of the particle self-disturbance, i.e. $\tilde{\bf u}_{p}$ must account for the background turbulent flow altered by the disturbances generated by all the other particles except the $p$th one. In the two-way coupling regime where the particle back-reaction modifies the carrier flow, the calculation of $\tilde{\bf u}_{p}$ needs the removal from the field ${\bf u}({\bf x},t)$ of the particle self-interaction contribution. In the ERPP approach the (regularised) disturbance flow induced at time $t$ and position ${\bf x}$ by a particle located at ${\bf x}_{0}$, i.e. ${\bf v}_{R}({\bf x}-{\bf x}_{0},t)$, is known in closed form. The actual hydrodynamic force on the $p$th particle can be evaluating by subtracting from ${\bf u}({\bf x}_{p},t)$ the value ${\bf v}_{R}[{\bf x}_{p}(t)-{\bf x}_{p}(t-\Delta t),\Delta t]$ induced at time $t$ at the current particle position ${\bf x}_{p}(t)$ by the same particle when it was placed at ${\bf x}_{p}(t-\Delta t)$,

[TABLE]

where ${\bf D}^{n}={\bf D}(t_{n}-\epsilon_{R})$, ${\bf r}={\bf x}-{\bf x}_{p}(t_{n}-\epsilon_{R})$, the hat denotes the unit vector $\hat{{\bf r}}={\bf r}/r$, $\eta=r/\sqrt{2}\sigma$ is the dimensionless distance with $\sigma=\sqrt{2\nu(\epsilon_{R}+\Delta t)}$ and $\displaystyle f(\eta)=\frac{\sqrt{\pi}}{2}\mbox{erf}(\eta)-\eta e^{-\eta^{2}}$, see [29] for the formal derivation of equation (19). This procedure can be straightforwardly extended to the Runge-Kutta algorithm employed in the present Letter to integrate in time the equations of the carrier and of the disperse phase.

0.3 The homogeneous shear flow

The homogeneous shear flow consists in a turbulent flow into a periodic box in which velocity fluctuations are fed by an imposed mean velocity profile, see [40, 41] for more details. The velocity field ${\bf v}$ is decomposed into a mean flow ${\bf U}=Sx_{2}\,{\bf e}_{1}$ and a fluctuation ${\bf u}$ where ${\bf e}_{1}$ is the unit vector in the $x_{1}$ direction (streamwise), $x_{2}$ denotes the coordinate in the direction of the mean constant shear $S$ (transverse direction) and $x_{3}$ is the spanwise direction. The dynamics of the velocity fluctuations in a deforming coordinate system convected by the mean flow according to the transformation of variables $\xi_{1}=x_{1}-Stx_{2};\quad\xi_{2}=x_{2};\quad\xi_{3}=x_{3};\quad\tau=t$ is described by the incompressible Navier-Stokes equations [42]

[TABLE]

In the homogeneous shear flow the Reynolds shear stresses $\langle u_{1}\,u_{2}\rangle$ extract energy from the mean flow and feeds the energy cascade up to viscous dissipation according to the balance $-S\langle u_{1}\,u_{2}\rangle=\epsilon$ where $\epsilon$ is energy dissipation rate. Turbulent fluctuations are spatially homogeneous and statistically stationary in time. Beyond the integral scale $L_{0}=\left(2K\right)^{3/2}/\epsilon$ and Kolmogorov dissipative scale $\eta=\left(\nu^{3}/\epsilon\right)^{1/4}$, the homogeneous shear flow features the so-called shear scale $L_{S}=\sqrt{\epsilon/S^{3}}$, where $K$ is the average turbulent kinetic energy. Due to the shear, turbulent fluctuations are strongly anisotropic at large scales $L_{S}\ll r\ll L_{0}$ (production range) where the production associated with the Reynolds stresses overwhelms the other mechanisms. At small scales, $\eta\ll r\ll L_{S}$ (isotropy recovery range) inertial energy transfer prevails leading to the classical Kolmogorov energy cascade. Two dimensionless parameters characterise the flow, the Corrsin parameter, $S_{c}=\sqrt{S^{2}\nu/\epsilon}=\left(\eta/L_{S}\right)^{2/3}$, and the shear strength, $S^{*}=2KS/\epsilon=\left(L_{0}/L_{S}\right)^{2/3}$, which ratio corresponds to the classical Taylor-Reynolds number, $Re_{\lambda}=2K/\sqrt{\nu\epsilon}=S^{*}/S_{c}$.