2.1. Particle-resolved simulations (PR-DNS)

The motion of a rigid particle can be described using the translational velocity $\boldsymbol {u}_p$ and the rotational velocity $\boldsymbol {\omega }_p$ of its centre of mass, which obey the classical Newton–Euler equations for rigid-body dynamics:

(2.2) \begin{equation} m_p \frac {\text {d} \boldsymbol {u}_p}{\text {d} t} = \boldsymbol {f}^{\leftarrow f} + \boldsymbol {f}_p^{\leftrightarrow p}, \quad I_p \frac {\text {d} \boldsymbol {\omega }_p}{\text {d} t} = \boldsymbol {L}_p^{\leftarrow f}, \end{equation}

where $m_p = \pi \rho _p D_p^3/6$ and $I_p = m_p D_p^2/10$ are the mass and inertial moment of the particle, with $\rho _p$ being the particle density and $D_p$ the particle diameter. Here $\boldsymbol {f}^{\leftrightarrow p}$ is the force due to particle collisions, while $\boldsymbol {f}^{\leftarrow f}$ and $\boldsymbol {L}_p^{\leftarrow f}$ are the force and momentum due to the fluid–solid interaction, namely

(2.3) \begin{equation} \boldsymbol {f}^{\leftarrow f} = \oint _{\partial V_p} \boldsymbol {\sigma } \cdot \boldsymbol {n} \text {d}A \quad \text {and} \quad \boldsymbol {L}_p^{\leftarrow f} = \oint _{\partial V_p} \boldsymbol {r} \times \left ( \boldsymbol {\sigma } \cdot \boldsymbol {n} \right ) \text {d} A, \end{equation}

where $\boldsymbol {\sigma }= -p \mathcal {I} + 2 \mu \mathcal {D}$ is the Cauchy stress tensor, with $\mathcal {I}$ being the identity tensor, $\mu$ the fluid kinematic viscosity, $\mathcal {D}$ the strain rate tensor and $\boldsymbol {n}$ the unit vector normal to the surface of the particle.

2.2. One-way-coupled point-particle simulations (PP-DNS)

The interface-resolved simulations are complemented with PP-DNS. Here the particles move under the action of the fluid flow, but they do not modulate it; the particle–particle interactions are also neglected. In this work, we consider the complete governing equation for a point particle as introduced by Maxey & Riley (Reference Maxey and Riley1983) and Gatignol (Reference Gatignol1983), i.e.

(2.4) \begin{equation} \begin{aligned}[t] \rho _p V_p \frac {\text {d}\boldsymbol {u}_p}{\text {d}t} &= \underbrace {3 \pi D_p \rho _f \nu \left (\boldsymbol {u} - \boldsymbol {u_p} + \frac {1}{6} \left ( \frac {D_p}{2} \right )^2 \boldsymbol {\nabla }^2 \boldsymbol {u} \right ) }_{\textit{Stokes drag}} \\ &\quad+ \underbrace {\frac {\rho _f V_p}{2} \left ( 3\frac {D\boldsymbol {u}}{Dt} - \frac {\text {d}\boldsymbol {u}_p}{\text {d}t} + \frac {1}{10} \left ( \frac {D_p}{2} \right )^2 \frac {\text {d} }{\text {d} t} \boldsymbol {\nabla }^2 \boldsymbol {u} \right )}_{\textit{Added mass}} \\ &\quad + \underbrace {\frac {3}{2} D_p^2 \rho _f \sqrt {\pi \nu } \int _{-\infty }^t K_B(t-\tau ) \left ( \frac {\text {d}\boldsymbol {u}}{\text {d}\tau } - \frac {\text {d}\boldsymbol {u_p}}{\text {d}\tau }+ \frac {1}{6} \left ( \frac {D_p}{2} \right )^2 \frac {\text {d}}{\text {d} t} \boldsymbol {\nabla }^2 \boldsymbol {u} \right ) \text {d} \tau }_{\textit{Basset force}}. \end{aligned} \end{equation}

Here, $V_p = \pi \rho _p D_p^3/6$ is the volume of the particle and $D/Dt$ denotes the material derivative. Note that the Faxén correction (Faxén Reference Faxén1922) proportional to $\boldsymbol {\nabla }^2 \boldsymbol {u}$ has been included in the Stokes drag, added mass and Basset forces. According to Homann & Bec (Reference Homann and Bec2010), this correction reproduces dominant finite-size effects on velocity and acceleration fluctuations for neutrally buoyant particles with diameter up to $D_p/\eta \approx 4$ . For the added mass, we have used the form described by Auton et al. (Reference Auton, Hunt and Prud’Homme1988). The computation of the Basset history force presents some challenges, and its evaluation can become extremely time-consuming and memory-demanding; indeed, this term requires at each time step the computation of an integral over the complete time history of the particle. Over the years, several attempts have been proposed to approximate this term (see e.g. Michaelides Reference Michaelides1992; Dorgan & Loth Reference Dorgan and Loth2007; Prasath et al. Reference Prasath, Vasan and Govindarajan2019). In this work, we resort to the second-order and memory-efficient algorithm developed by van Hinsberg et al. (Reference van Hinsberg, ten Thije Boonkkamp and Clercx2011), the details of which are briefly reported in Appendix A for completeness.

2.3. Computational details

We consider a single-phase microscale Reynolds number of $Re_\lambda = u' \lambda / \nu \approx 140$ to ensure a relatively large inertial range of scales; here $u'$ is the root mean square of the velocity fluctuations and $\lambda$ is the Taylor length scale. The particle diameter is set to $D_p/\eta \approx 0.9$ , where $\eta$ is the Kolmogorov length scale for the single-phase case. Two volume fractions are considered, i.e. $\varPhi _V = V_s/(V_s + V_f)= 10^{-5}$ and $\varPhi _V = 10^{-3}$ for a total number of particles of $N=742$ and $N=74208$ , respectively. For each volume fraction two values of the particle density are considered, $\rho _p/\rho _f = 5$ and $\rho _p/\rho _f=100$ , to consider both light and heavy particles. This leads to a total of four PR-DNS. In PP-DNS the particles do not modulate the flow and do not interact; therefore, the volume fraction is not a parameter. Because of this, only two PP-DNS simulations have been carried out for the different density ratios.

The governing equations are numerically integrated in time using the in-house solver Fujin (https://groups.oist.jp/cffu/code). It solves the Navier–Stokes equations using an incremental pressure-correction scheme. The governing equations are written in primitive variables on a staggered grid, and second-order finite differences are used in all the directions. The Adams–Bashforth time scheme is used for advancing the momentum equation in time. The Poisson equation for the pressure enforcing the incompressibility constraint is solved using a fast and efficient approach based on the fast Fourier transform.

For the PR-DNS the governing equations for the particles are dealt with by the immersed boundary method introduced by Hori et al. (Reference Hori, Rosti and Takagi2022). The fluid–solid coupling is achieved in an Eulerian framework, and accounts for the inertia of the fictitious fluid inside the solid phase, so as to properly reproduce the particles’ behaviour both in the neutrally buoyant case and in the presence of density difference between the fluid and solid phases. The soft sphere collision model (Cundall & Strack Reference Cundall and Strack1979; Tsuji et al. Reference Tsuji, Kawaguchi and Tanaka1993) is used to prevent the interpenetration between particles and compute the $\boldsymbol {f}^{\leftrightarrow p}$ term in (2.2). In this model the particles are allowed to slightly penetrate. The collision is modelled as a spring and dashpot dynamical system, with the collision force being proportional to the penetration depth between the particles and to their relative normal velocity. A fixed-radius near-neighbours algorithm (see Monti et al. (Reference Monti, Rathee, Shen and Rosti2021), and references therein) is used for the particle interaction to avoid an otherwise prohibitive increase of the computational cost when the number of particles becomes large. The PR-DNS are performed without any additional lubrication correction, i.e. we consider only the lubrication that naturally arises from the method. For the one-way-coupled PP-DNS the governing equation for the particle velocity is advanced in time using the second-order Adams–Bashforth time scheme. At each time step, the fluid velocity is evaluated at the position of the particle using a second-order linear interpolation.

For the PR-DNS the fluid domain is discretised on a cubic domain using $N_p = 2048$ points along each direction, leading to $\eta / \Delta x \approx D_p /\Delta x \approx 6{-}7$ , where $\Delta x$ denotes the grid spacing. The accuracy of the results has been verified by running an additional simulation on a coarser grid with $N_p=1440$ for the case with $\varPhi _V = 10^{-3}$ and $\rho _p/\rho _f = 100$ , resulting in a negligible difference in the scale-by-scale fluid energy spectrum and budget in figures 2 and 5, and in the Lagrangian and Eulerian particles’ statistics in figures 14 and 15. For the single-phase case and the PP-DNS the fluid domain is discretised using $N_p = 1024$ points in the three directions, leading to $\eta /\Delta x \approx 3{-}4$ . In this case the flow around the particles does not need to be solved. Excluding the initial transient period, all simulations are advanced in time for approximately $50 \tau _f$ , where $\tau _f = \mathcal {L}/\sqrt {2 \! \left \langle {E} \right \rangle \!/3}$ is the average turnover time of the largest eddies; $\mathcal {L} =\pi /(4 \! \left \langle {E} \right \rangle \!/3)\int _0^\infty \mathcal {E}(\kappa )/\kappa \text {d} \kappa$ is the fluid integral scale with $\mathcal {E}(\kappa )$ being the energy spectrum; $E(\boldsymbol {x},t)$ is the local and instantaneous fluid kinetic energy; and the $ \! \left \langle {\cdot } \right \rangle \!$ operator denotes averages in space and time.

Details of the PR-DNS are reported in table 1. Note that, when looking at the bulk quantities, the flow modulation due to the solid phase is rather low. For comparison, we mention that at a reference Reynolds number of $Re_\lambda \approx 240$ , Hwang & Eaton (Reference Hwang and Eaton2006) experimentally report an energy attenuation of $ \! \left \langle {E} \right \rangle \!/ \! \left \langle {E_0} \right \rangle \! = 0.78$ for a suspension of spherical particles with $D_p/\eta = 0.9$ , $\rho _p/\rho _f = 2048$ and $\varPhi _V \approx 10^{-4}$ , which corresponds to a mass loading of $\varPhi _M = \rho _p V_p/(\rho _f V_f) = 0.1$ (the subscript ‘0’ refers to the single-phase case). Although those parameters differ from ours, the value of the mass loading matches that of the $\rho _p/\rho _f = 100$ and $\varPhi _V = 10^{-3}$ case ( $\varPhi _M \approx 0.098$ ), for which we, however, found a lower energy attenuation, i.e. $ \! \left \langle {E} \right \rangle \!/ \! \left \langle {E_0} \right \rangle \! = 0.9 \pm 0.26$ . This suggests that the value of the mass loading alone is not sufficient to predict the level of the turbulence attenuation by Kolmogorov-size particles. For completeness, we mention that our predictions are close to those obtained by two-way-coupled PP-DNS for similar parameters, although at lower Reynolds numbers. For example, Boivin et al. (Reference Boivin, Simonin and Squires1998) predict an attenuation of $ \! \left \langle {E} \right \rangle \!/ \! \left \langle {E_0} \right \rangle \! \approx 0.88$ for $D_p/\eta = 0.11$ , $Re_\lambda = 62$ and $\varPhi _M = 0.2$ . At $Re_\lambda = 38$ , $D_p/\eta \lt 1$ , $\rho _p/\rho _f \gg 1$ and $\varPhi _M = 0.1$ , instead, Squires & Eaton (Reference Squires and Eaton1990) report an attenuation of $ \! \left \langle {E} \right \rangle \!/ \! \left \langle {E_0} \right \rangle \! = 0.93$ .

Figure 2. Energy spectrum for (a) $\varPhi _V = 10^{-5}$ and (b) $\varPhi _V = 10^{-3}$ . The black line refers to the single-phase case, the red/blue light line is for $\rho _p/\rho _f = 5$ and the red/blue dark line is for $\rho _p/\rho _f = 100$ . The symbols in (b) are from the simulation carried out with the coarser grid. The filled circle on the $\kappa$ axis denotes the wavenumber associated with the particle size. The solid line is for the Kolmogorov $\kappa ^{-5/3}$ scaling. The dashed line is for $\kappa ^{-4}$ . Here $\kappa _L = 2 \pi /L$ .

Table 1. Details of the PR-DNS considered in the present parametric study. Here $\epsilon$ is the dissipation rate; $St$ is the Stokes number defined as $St=\tau _p/\tau _f$ , where $\tau _p=(\rho _p/\rho _f)D_p^2/(18 \nu )$ is the relaxation time of the particle and $\tau _f=\mathcal {L}/\sqrt {2 \! \langle {E} \rangle \!/3}$ is the turnover time of the largest eddies; $Re_p$ is the particle Reynolds number defined as $Re_p = |\Delta \boldsymbol {u}| D /\nu$ , where $\Delta \boldsymbol {u} = \boldsymbol {u}_p - \boldsymbol {u}_f$ is the fluid–particle relative velocity. Here, $\boldsymbol {u}_f$ is the fluid velocity seen by the particle evaluated as the average of the fluid velocity within a shell centred with the particle and with radius $R_{sh} = 3(D_p/2)$ (see Uhlmann & Chouippe Reference Uhlmann and Chouippe2017; Chiarini & Rosti Reference Chiarini and Rosti2024).

3.1. Energy spectrum

Figure 2 shows the influence of the solid phase on the energy spectrum $\mathcal {E}(\kappa )$ , and highlights how Kolmogorov-size particles modulate the turbulent fluctuations scale by scale. Figures 2(a) and 2(b) are for $\varPhi _V = 10^{-5}$ and $\varPhi _V=10^{-3}$ , respectively. The black solid line refers to the reference unladen case; note the inertial range of scales where $\mathcal {E}(\kappa ) \sim \kappa ^{-5/3}$ extends for more than one decade of wavenumbers, confirming that the present Reynolds number is large enough to ensure a proper separation of scales. For validation purposes, in figure 2(b) we also plot with symbols the energy spectrum obtained for the $\varPhi _V=10^{-3}$ and $\rho _p/\rho _f=100$ case with the coarser grid, showing good agreement with that obtained with the standard grid, thus ensuring the suitability of the chosen grid resolution (see § 2). We compute the energy spectra (and the structure functions in the next section) using all the grid points of the computational domain including those that are inside the particles. However, as mentioned in Chiarini & Rosti (Reference Chiarini and Rosti2024), with the present immersed boundary method the results do not change when the points within the particles are neglected.

At large scales, the spectra of the particle-laden cases substantially overlap with the unladen spectrum. We indeed observe only a weak depletion of the energy content at the intermediate scales (see the insets in figure 2). At scales smaller than a certain wavenumber $\kappa _p$ , the energy spectra of the particle-laden cases deviate above the reference spectrum. Solid particles enhance the energy content of the small scales by locally deforming the fluid flow around them. Notably, this mechanism is amplified as $\varPhi _V$ and/or $\rho _p/\rho _f$ increase, as conveniently visualised by the larger values of $\mathcal {E}(\kappa )$ for $\kappa \gt \kappa _p$ and by the shift of $\kappa _p$ towards smaller wavenumbers. However, notice that due to the low values of $\varPhi _V$ considered, the flow modulation is rather low for all cases, being substantially negligible for $\varPhi _V=10^{-5}$ and/or $\rho _p/\rho _f=5$ .

Figure 2 shows that the modulated energy spectrum exhibits multiscaling behaviour. The classical $\mathcal {E}(\kappa ) \sim \kappa ^{-5/3}$ decay in the inertial range of scales is indeed followed by a steeper decay $\mathcal {E}(\kappa ) \sim \kappa ^{-4}$ starting at wavenumbers close to $\kappa _p$ . A similar steep decay has been observed in bubbly flows at scales smaller than the bubble diameter, by means of both experiments (Martinez Mercado et al. Reference Martinez Mercado, Chehata Gomez, Van Gils, Sun and Lohse2010; Riboux et al. Reference Riboux, Risso and Legendre2010; Prakash et al. Reference Prakash, Martínez Mercado, van Wijngaarden, Mancilla, Tagawa, Lohse and Sun2016; Dung et al. Reference Dung, Waasdorp, Sun, Lohse and Huisman2023) and simulations (Pandey et al. Reference Pandey, Ramadugu and Perlekar2020, Reference Pandey, Mitra and Perlekar2022). Additionally, a similar multiscaling behaviour has also been observed in a turbulent planar Couette flow laden with small particles (Wang et al. 2023), and in homogeneous isotropic turbulence laden with slender fibres (Olivieri et al. Reference Olivieri, Cannon and Rosti2022a ,Reference Olivieri, Mazzino and Rosti b ). In the context of bubbly flows, the emergence of the $\kappa ^{-\alpha }$ decay with $\alpha \geqslant 3$ has been attributed to the wakes that the bubbles generate in an otherwise smooth flow (see e.g. Alméras et al. Reference Alméras, Mathai, Lohse and Sun2017; Pandey et al. Reference Pandey, Ramadugu and Perlekar2020); here the velocity fluctuations produced by the bubbles are directly dissipated by viscosity (Lance & Bataille Reference Lance and Bataille1991; Risso Reference Risso2018). Accordingly, this scaling has been indeed observed only when the bubble Reynolds number is large enough and is in the $10 \leqslant Re_{bub} \leqslant 1000$ range (Pandey et al. Reference Pandey, Ramadugu and Perlekar2020). To compare, we computed the local particle Reynolds number using the relative velocity between the particle and the surrounding flow, and found that it is in the $0.25 \lessapprox Re_p \lessapprox 2.5$ range (see table 1). More recently, Zamansky et al. (Reference Zamansky, De Bonneville and Risso2024) observed that in bubbly flows the $\kappa ^{-3}$ range constitutes a transition between the production-dominated (anisotropic) and the inertial/dissipative (isotropic) ranges of scales, and proposed that the specific $\alpha =3$ scaling results from the mean shear rate imposed by the bubbles which controls the rate of return to isotropy. In view of this, we mention that unlike what is commonly found in bubbly-induced agitation where the $\kappa ^{-3}$ subrange is between the energetic and the inertial range of scales (see also figure 5 of Risso (Reference Risso2018)), in our case the $\kappa ^{-4}$ scaling is observed at the small dissipative scales. It is also worth mentioning that in a recent work Ramirez et al. (Reference Ramirez, Burlot, Zamansky, Bois and Risso2024) investigated the effect of singularities (e.g. induced by the particles) of various orders on the energy spectrum. They showed that singularities may affect the spectrum beyond the imposition of simple oscillations (Lucci et al. Reference Lucci, Ferrante and Elghobashi2010), and actually cause power-law scalings at wavenumbers $\kappa \geqslant \kappa _p$ , with the exact slope being dependent on the order of the singularity. However, in most of the cases investigated by us, the start of the $\kappa ^{-4}$ scaling region appears at wavenumbers sensibly lower than $\kappa _p$ . Further investigation on the link between the fluctuations induced by particles and the $\mathcal {E}(\kappa ) \sim \kappa ^{-4}$ decay is provided in § 3.3 by looking at the scale-by-scale energy budget.

Overall, figure 2 shows that for the present parameters, particles almost do not modulate the inertial range of scales, where the classical energy cascade described by Richardson and Kolmogorov is preserved, but mainly affect the (otherwise smooth) smallest scales of the flow.

3.2. Structure function and intermittency

Figure 3. Structure functions for (a) $\varPhi _V = 10^{-5}$ and (b) $\varPhi _V = 10^{-3}$ . For each panel, from bottom to top the plots are for $S_2$ , $S_4$ and $S_6$ . The dashed lines represent $r^p$ , while the dash-dotted ones $r^{p/3}$ . The insets show a zoom of $S_2$ for $10 \leqslant r/\eta \leqslant 40$ .

We extend the analysis done in the spectral domain by computing the longitudinal structure functions defined as $S_p(r) = \! \left \langle {\delta u(r)^p} \right \rangle \!$ , where $\delta u(r) = (\boldsymbol {u}(\boldsymbol {x}+\boldsymbol {r}) - \boldsymbol {u}(\boldsymbol {x}) )\cdot \boldsymbol {r}/r$ and $r = |\boldsymbol {r}|$ . In particular, figure 3 plots $S_2$ , $S_4$ and $S_6$ as a function of $r$ for $\varPhi _V=10^{-5}$ and $\varPhi _V = 10^{-3}$ . In the single-phase case, we observe that $S_p \sim r^{p/3}$ in the inertial range of scales is in agreement with the Kolmogorov prediction (Kolmogorov Reference Kolmogorov1941), and that $S_p \sim r^p$ at the small scales, as a result of the differentiability of the fluid velocity field (Schumacher et al. Reference Schumacher, Sreenivasan and Yakhot2007). Recall that although $S_2(r)$ is commonly referred to as scale energy (Davidson Reference Davidson2004; Davidson & Pearson Reference Davidson and Pearson2005), its meaning slightly differs from that of $\mathcal {E}(\kappa )$ . Indeed while $\mathcal {E}(\kappa ) \text {d}\kappa$ refers to the amount of energy associated with the scale $r = 2\pi /\kappa$ , $S_2(r)$ can be interpreted as the amount of energy associated with scales up to $r$ . In agreement with the modulation of the energy spectrum, figure 3 shows that particles enhance the energy content at small scales compared with the unladen case. The energy enhancement is more intense for larger $\varPhi _V$ and $\rho _p/\rho _f$ , and becomes more evident when considering higher-order structure functions. At the same time, for $\rho _p/\rho _f=100$ the presence of the particles decreases the amount of energy stored at the larger scales (as seen by the dark blue and dark red curves lying below the black one in the insets of figure 3). By interacting with the vortical structures of the flow, the (heavy) particles drain energy from scales larger than $D$ and reinject it back at smaller scales.

Figure 4. Extended self-similarity for (a) $\varPhi _V = 10^{-5}$ and (b) $\varPhi _V = 10^{-3}$ . Here $S_6/S_2^3$ is plotted against $r$ .

Structure functions are commonly employed to quantify the flow intermittency, i.e. the relevance of extreme events that are localised in space and time and break the Kolmogorov similarity hypothesis (Frisch Reference Frisch1995; Pope Reference Pope2000). In figure 4, we use the extended self-similarity introduced by Benzi et al. (Reference Benzi, Ciliberto, Tripiccione, Baudet, Massaioli and Succi1993), and plot $S_6/S_2^3$ as a function of $r$ . In the limit case where extreme events do not occur, the $S_6 \sim S_2^3$ power law holds, i.e. $S_6/S_2^3 \sim \text{const.}$ , and deviations from this behaviour are a measure of the flow intermittency. Accordingly with the intrinsic intermittent nature of turbulent flows, figure 4 shows that $S_6/S_2^3$ deviates from the Kolmogorov prediction also in the single-phase case, and this deviation increases in the particle-laden cases, similarly to what is found for suspension of particles with size in the inertial range of scales (Chiarini & Rosti Reference Chiarini and Rosti2024). This is due to the no-slip and no-penetration boundary conditions at the particle surface that give origin to localised and intense velocity gradients. Figure 4 shows that for $\rho _p/\rho _f = 5$ the deviation from the single-phase case is rather small for both $\varPhi _V = 10^{-5}$ and $\varPhi _V=10^{-3}$ , in agreement with the low level of flow modulation discussed above. For $\rho _p/\rho _f = 100$ , instead, the $S_6/S_2^3$ curve significantly deviates from the single-phase case at small scales. For the considered parameters, heavy Kolmogorov-size particles increase the intermittency of the small scales of the flow.

3.3. Scale-by-scale energy budget

We now detail the influence of the particles on the organisation of the fluctuations, by studying the scale-by-scale energy budget equation and characterising thus the dominant energetic mechanisms in the different regimes identified with the energy spectra (see figure 2). For the present case with three homogeneous directions, the energy balance reads

(3.1) \begin{equation} P(\kappa ) + \varPi (\kappa ) + \varPi _{fs}(\kappa ) - D_v(\kappa ) = 0, \end{equation}

where $P(\kappa )$ is the scale-by-scale turbulent energy production due to the external forcing, $\varPi (\kappa )$ is the energy flux associated with the nonlinear convective term, $\varPi _{fs}(\kappa )$ is the fluid–solid coupling term and $D_v(\kappa )$ is the scale-by-scale viscous dissipation. Specifically, these terms are defined as

(3.2) \begin{align} P( \kappa ) = & \int _\kappa ^\infty \frac {1}{2} \left ( \,\hat {\boldsymbol {\!f}} \cdot \hat {\boldsymbol {u}}^* + \,\hat {\boldsymbol {\!f}}^* \cdot \hat {\boldsymbol {u}} \right ) \text {d} k, \end{align}

(3.3) \begin{align} \varPi (\kappa ) = & \int _\kappa ^\infty -\frac {1}{2} \left ( \hat {\boldsymbol {G}} \cdot \hat {\boldsymbol {u}}^* + \hat {\boldsymbol {G}}^* \cdot \hat {\boldsymbol {u}} \right ) \text {d} k, \end{align}

(3.4) \begin{align} \varPi _{fs}(\kappa ) = & \int _\kappa ^\infty \frac {1}{2} \left ( \,\hat {\boldsymbol {\!f}}^{\leftrightarrow p} \cdot \hat {\boldsymbol {u}}^* + \,\hat {\boldsymbol {\!f}}^{\leftrightarrow p,*} \cdot \hat {\boldsymbol {u}} \right ) \text {d} k, \end{align}

(3.5) \begin{align} D_v(\kappa ) = & \int _\kappa ^\infty \left ( 2 \nu k^2 \mathcal {E} \right ) \text {d} k, \end{align}

where a hat denotes the Fourier transform operator and the superscript asterisk denotes complex conjugate. The term $\hat {\boldsymbol {G}}$ is the Fourier transform of the nonlinear term $\boldsymbol {\nabla } \cdot ( \boldsymbol {u} \boldsymbol {u} )$ . Note that here we integrate all the terms from $\kappa$ to $\infty$ . Terms $\varPi (\kappa )$ and $\varPi _{fs}(\kappa )$ do not produce nor dissipate energy at any scale, but redistribute it among scales by means of the classical energy cascade and of the fluid–solid interaction. Also, note that since we integrate the viscous term from $\kappa$ to $\infty$ , $D_v(0) = \epsilon$ . For the complete derivation we refer the reader to Pope (Reference Pope2000).

Figure 5. Scale-by-scale energy budget for (a,b) $\varPhi _V = 10^{-5}$ and (c,d) $\varPhi _V = 10^{-3}$ . Plots for (a,c) $\rho _p/\rho _f = 5$ and (b,d) $\rho /\rho _f = 100$ . The production term $P$ acts at the largest scales $\kappa /\kappa _L \leqslant 1$ only, being $P=0$ for $\kappa /\kappa _L \geqslant 2$ (not visible as the $y$ axes are in log scale). The filled symbols in (d) are from the simulation carried out with the coarser grid; circles are for $\varPi /\epsilon$ , triangles for $\varPi _{fs}/\epsilon$ and squares for $D_v/\epsilon$ .

Figure 5 plots the terms of equation (3.1) as a function of $\kappa$ for the four particle-laden cases investigated. For validation purposes, we also plot with filled symbols the terms obtained with the coarser grid for $\varPhi _V=10^{-3}$ and $\rho _p/\rho _f=100$ . In agreement with the multiscaling behaviour shown in the energy spectra (see figure 2), the energy budgets exhibit two distinct behaviours. Energy is injected at the largest scales at a rate that equals the dissipation rate $P(0) = \epsilon$ . In the inertial range of scales $\kappa _L \lt \kappa \lt \kappa _p$ , the fluid–solid coupling term is subdominant and

(3.6) \begin{equation} \varPi (\kappa ) \approx D_v(\kappa ) \approx \epsilon . \end{equation}

Thus, $\varPi (\kappa ) \approx \epsilon$ is constant with $\kappa$ at these scales, exhibiting a plateau. In agreement with the Kolmogorov theory, the viscous effects are negligible $2 \nu \kappa ^2 \mathcal {E}(\kappa ) \approx 0$ and energy is transferred from larger to smaller scales at a rate that matches the energy injection rate $\epsilon$ . This corresponds to the range of scales where $\mathcal {E}(\kappa ) \sim \kappa ^{-5/3}$ . Similarly to what is observed in the energy spectrum, the range of scales where this relation holds shrinks as $\varPhi _V$ and $\rho _p/\rho _f$ increase. For the small scales with $\kappa \gt \kappa _p$ where the spectrum shows the $\mathcal {E}(\kappa ) \sim \kappa ^{-4}$ decay, instead, the nonlinear flux is negligible $\varPi (\kappa ) \approx 0$ . In this range of scales, the viscous effects and the fluid–solid coupling term are not negligible, and the energy budget reduces to

(3.7) \begin{equation} \varPi _{fs}(\kappa ) \approx D_v(\kappa ). \end{equation}

Here the fluctuations that are produced by the fluid–solid interaction are (on average) directly dissipated by viscosity and are not transferred among scales. In agreement with the energy spectrum, the range of scales where this regime holds widens as $\varPhi _V$ and/or $\rho _p/\rho _f$ are increased.

3.4. The local structure of the flow

As shown in § 3.1, Kolmogorov-size particles mainly modify the organisation of the velocity fluctuations at the smallest scales. Particles indeed modulate the energy spectrum and the structure functions for $\kappa \gtrapprox \kappa _p$ only. Here we characterise the smallest scales of the flow to provide new insights of the influence of the dispersed phase on the structure of the velocity fluctuations. For the sake of brevity, in this section we consider the $\varPhi _V = 10^{-3}$ cases only, and we investigate how particles modify the velocity gradient field $A_{ij} = \partial u_j/\partial x_i$ . In the neighbourhood of a given point $(\boldsymbol {x}_0,t)$ , indeed, the velocity field can be approximated as $u_i(\boldsymbol {x},t) = u_i(\boldsymbol {x}_0,t) + A_{ij}(\boldsymbol {x}_0,t)(x_j-x_{0,j}) + \mathcal {O}(|\boldsymbol {x}-\boldsymbol {x}_0|^2)$ . This linear expansion is valid in the region around $\boldsymbol {x}_0$ , where the fluid is sufficiently smooth and the variations of $A_{ij}$ are small (Meneveau Reference Meneveau2011); for a turbulent flow the extent of this region is of the order of the Kolmogorov scale $\eta$ . Based on these arguments, we study the influence of the particles on the smallest scales of the flow, by inspecting their effect on the $A_{ij}$ tensor.

We decompose $A_{ij}$ into its symmetric and antisymmetric parts, namely the strain-rate tensor $S_{ij} = (A_{ij} + A_{ji})/2$ and the rotation-rate tensor $W_{ij} = (A_{ij} - A_{ji})/2$ . The field of the velocity gradient is completely addressed when knowing: (i) the three principal rates of strain $\alpha \geqslant \beta \geqslant \gamma$ , i.e. the three eigenvalues of $S_{ij}$ , (ii) the magnitude of the vorticity $\omega ^2 = \boldsymbol {\omega } \cdot \boldsymbol {\omega }$ , i.e. the enstrophy, and (iii) the orientation of $\boldsymbol {\omega }$ relative to the three principal axes of strain, i.e. the eigenvectors of $S_{ij}$ (Davidson Reference Davidson2004). Note that, due to the incompressibility constraint $\alpha + \beta + \gamma = 0$ , meaning that $\alpha$ is always non-negative, $\gamma$ is always non-positive while $\beta$ can have any sign depending on the local straining state.

Figure 6. Probability density function (PDF) of $s^*$ for $\varPhi _V = 10^{-3}$ , with $\rho _p/\rho _f=5$ and $\rho _p/\rho _f=100$ .

We start by characterising $\alpha$ , $\beta$ and $\gamma$ in figure 6. Following Lund & Rogers (Reference Lund and Rogers1994), we use $s^*$ which is defined as

(3.8) \begin{equation} s^* = - \frac {3 \sqrt {6} \alpha \beta \gamma } {(\alpha ^2 + \beta ^2 + \gamma ^2)^{3/2}}. \end{equation}

For a random velocity gradient field with no preferred structure, the distribution of $s^*$ is uniform. When $s^*=1$ , $\alpha = \beta = - \gamma /2$ , meaning that the state of straining is an axisymmetric extension, in which a small spherical fluid element moving in the flow extends symmetrically in two directions and contracts in the third one, forming thus a disk-like structure. When $s^*=-1$ , instead, $\gamma = \beta = - \alpha /2 \lt 0$ , and the state of straining is an axisymmetric compression, in which a small fluid element contracts in two directions and extends in the third one, forming thus a vortex tube. Finally, when $s^*=0$ we have $\beta =0$ , meaning that the straining state is two-dimensional, as typical for shear dominated regions.

In the absence of particles, the distribution peaks at $s^*=1$ , in agreement with the fact that for purely Newtonian turbulence the most likely state of straining is an axisymmetric extension (Davidson Reference Davidson2004; Meneveau Reference Meneveau2011). Figure 6 shows that the addition of Kolmogorov-size particles with $\varPhi _V \leqslant 10^{-3}$ leads to a rather small variation of the distribution of $s^*$ for the lighter particles with $\rho _p/\rho _f = 5$ . We observe instead that particles with $\rho _p/\rho _f = 100$ decrease the probability of events with large positive $s^*$ and increase the probability of events with $s^* \leqslant 0$ . This agrees with the observation of Cannon et al. (Reference Cannon, Olivieri and Rosti2024) who found that particles with size in the inertial range favour the occurrence of events with two-dimensional straining states and with axisymmetric compression. These events indeed are associated with the shear layers that separate from the surface of the particles, that strengthen as $\rho _p/\rho _f$ increases.

Figure 7. (a) Probability density function (PDF) of the enstrophy. (b) Alignment of the vorticity vector $\boldsymbol {\omega }$ with the principal axes of strain. Here, $\hat {\boldsymbol {e}}_\alpha$ is aligned with the direction of maximum elongation of the flow, $\hat {\boldsymbol {e}}_\gamma$ is aligned with the direction of compression of the flow and $\hat {\boldsymbol {e}}_\beta$ is orthogonal to the previous two directions. The data shown are for $\varPhi _V = 10^{-3}$ .

Figure 7 shows the influence of the solid phase on ( figure 7 a) the square of the vorticity magnitude $\omega ^2$ , i.e. the enstrophy, and on (figure 7 b) the alignment between the vorticity $\boldsymbol {\omega }$ and the principal axes of strain. The distribution of $\omega ^2$ shows that the tail becomes longer and the probability of large events increases due to the presence of the particles. This agrees with the enhanced flow intermittency discussed in § 3.2. The tail of the distribution is longer for the $\rho _p/\rho _f = 100$ case, as the velocity gradients at the particle surface are more intense because of the larger relative velocity between the particles and the surrounding fluid phase. However, even the lighter particles with $\rho _p/\rho _f = 5$ are able to produce a non-negligible change of the distribution. Instead, figure 7(b) shows that the alignment between $\boldsymbol {\omega }$ and the principal axes of strain is only slightly influenced by the presence of the particles. The results for the single-phase case perfectly overlap with those of $\rho _p/\rho _f = 5$ . The presence of the heavy particles, instead, slightly reduces the alignment between $\boldsymbol {\omega }$ and $\hat {\boldsymbol {e}}_\beta$ , as well as the anti-alignment between $\boldsymbol {\omega }$ and $\hat {\boldsymbol {e}}_\gamma$ . Our results suggest that for $D_p/\eta \approx 1$ the perturbation field induced by the heavy particles is characterised by events where vorticity is more aligned with the direction of extension ( $\hat {\boldsymbol {e}}_\alpha$ ) and compression ( $\hat {\boldsymbol {e}}_\gamma$ ), and more anti-aligned with the intermediate eigenvector ( $\hat {\boldsymbol {e}}_\beta$ ). Interestingly, this differs from that observed by Cannon et al. (Reference Cannon, Olivieri and Rosti2024), who found the opposite trend for particles with size in the inertial range. We presume that the difference is due to the different particle Reynolds number, which results in a different kind of deformation of the local fluid flow.

We now move and consider the entire velocity gradient tensor $A_{ij}$ . Any second-order tensor possesses three invariants, $P$ , $Q$ and $R$ , which are directly related to its eigenvalues $\lambda$ by the characteristic polynomial function

(3.9) \begin{equation} \lambda ^3 + P \lambda ^2 + Q \lambda + R = 0. \end{equation}

Following Davidson (Reference Davidson2004) and Meneveau (Reference Meneveau2011), it can be shown that

(3.10) \begin{equation} \begin{aligned} P &= \alpha + \beta + \gamma , \\ Q &= - \frac {1}{2} \left ( \alpha ^2 + \beta ^2 + \gamma ^2 \right ) + \frac {\omega ^2}{4}, \\ R &= - \frac {1}{3} \left ( \alpha ^3 + \beta ^3 + \gamma ^3 \right ) - \frac {1}{4} \omega _i \omega _j S_{ij}. \end{aligned} \end{equation}

Note that $P=0$ due to the incompressibility constraint. The $Q$ and $R$ invariants are commonly used to distinguish between regions of intense vorticity and regions of strong strain. In particular, the discriminant of equation (3.9) $\varDelta = 27/4R^2 + Q^3$ is used to distinguish between regions where motions are mainly vortical (i.e. regions where $\varDelta \lt 0$ , meaning that $A_{ij}$ has one real and two complex conjugate eigenvalues) or characterised by a node-saddle streamline pattern (i.e. regions where $\varDelta \gt 0$ , and all the eigenvalues are real). When $Q$ is large and negative the strain is intense, while the vorticity is weak; in this case, $R \sim - (\alpha ^3 + \beta ^3 + \gamma ^3) = - \alpha \beta \gamma$ (Davidson Reference Davidson2004), and a positive $R$ implies a region of biaxial strain ( $\gamma \lt 0$ , $\alpha \gt \beta \gt 0$ ), while a negative $R$ implies a region of axial strain ( $\gamma \lt \beta \lt 0$ and $\alpha \gt 0$ ). When instead $Q$ is large and positive the strain is locally weak and $R \sim - \omega _i \omega _j S_{ij}$ . In this case, a positive $R$ implies vortex compression, while a negative $R$ implies vortex stretching.

Figure 8. The $Q$ – $R$ map for the (a) single phase and (b) $\varPhi _V=10^{-3}$ and $\rho _p/\rho _f = 100$ . For both panels the $11$ isolines have a logarithmic spacing between $5\times 10^{-5}$ and $10^2$ . Distribution of $Q$ (c) and $R$ (d) for $\varPhi _V=10^{-3}$ .

Figure 9. Volumetric rendering of the (a) $Q$ and (b) $R$ fields for $\varPhi _V = 10^{-3}$ and $\rho _p/\rho _f=100$ . Orange and magenta regions are associated with large, positive values of $Q$ and $R$ , i.e. $0.2072 \lessapprox Q \tau _\eta ^2 \lessapprox 2.072$ and $0.2109 \lessapprox R \tau _\eta ^3 \lessapprox 0.4218$ . Indigo and green regions indicate large, negative values of $Q$ and $R$ , i.e. $-2.072 \lessapprox Q \tau _\eta ^2 \lessapprox -0.2072$ and $-0.4218 \lessapprox R \tau _\eta ^3 \lessapprox -0.2109$ . The particle travelling direction relative to the local fluid velocity in a shell of radius $R_{sh}=5$ is indicated for each particle by a pointer. See figure 10 for a schematic representation of the flow.

Figure 8 plots the joint distribution of $Q$ and $R$ for the $\varPhi _V = 10^{-3}$ and $\rho _p/\rho _f = 100$ case (figure 8 b) and for the single phase (figure 8 a); the black, thin, solid (curved) lines denote the left- and right-Vieillefosse tails with $\varDelta =0$ . For completeness the figure shows the distributions of $Q$ (figure 8 c) and $R$ (figure 8 d). In the absence of the particles, the $Q$ – $R$ joint distribution takes a tear-drop pattern, with a clear point at the right-Vieillefosse tail with $R\gt 0$ and $Q\lt 0$ (Meneveau Reference Meneveau2011). The distribution is skewed towards positive $Q$ , but rather evenly distributed among positive and negative values of $R$ (figure 8 c,d). The largest probability is observed in the second and fourth quadrants, i.e. $Q\lt 0$ and $R\gt 0$ , and $Q\gt 0$ and $R\lt 0$ . Thus, in a purely Newtonian turbulent flow there is a strong negative correlation between $Q$ and $R$ , and the two most common states are vortex stretching $\omega _i \omega _j S_{ij}\gt 0$ and biaxial strain $\alpha \beta \gamma \lt 0$ (Betchov Reference Betchov1956; Davidson Reference Davidson2004). Also, the points in the $Q$ – $R$ map are distributed around the origin, since the mean values of $Q$ and $R$ are zero in homogeneous flow (Nomura & Post Reference Nomura and Post1998). The presence of the particles enlarges the range of possible $Q$ and $R$ values, in agreement with the increase of the probability of events with intense velocity gradients associated with the boundary conditions at the particles’ surface. In particular, Kolmogorov-size particles mainly favour events that lie in the first and third quadrants, resulting in a joint distribution that is more symmetric with respect to an inversion of the $R$ axis (see figure 8 d). Compared with the unladen case, particles mainly promote events with axial strain $\alpha \beta \gamma \gt 0$ ( $R\lt 0$ and large negative $Q\lt 0$ ) and with vortex compression $\omega _i \omega _j S_{ij}\lt 0$ ( $R\gt 0$ and $Q\gt 0$ ). The probability of events with $R\lt 0$ and $Q\lt 0$ is particularly enhanced, as visualised in figure 8 b by the occurrence of a point at the left-Vieillefosse tail. This is consistent with the increase of the probability of events with $s^*=-1$ shown in figure 6. A visual investigation of the $Q$ and $R$ fields around the particles (figure 9) helps to explain this effect by highlighting the local contribution of the particles. By comparing the two fields in the surroundings of the particles, we observe that they are both almost axisymmetric along the axis aligned with their travelling direction. However, across the velocity-normal median plane, $Q$ is symmetric while $R$ is antisymmetric. This relation between the two fields implies contributions across all four quadrants of the $Q$ – $R$ distribution (figure 8). The effect, however, is particularly apparent in the third quadrant, as the region is otherwise not explored by the single-phase flow. In particular, regions of $Q\lt 0$ and $R\lt 0$ are associated with axial strain and appear to be found at the downstream end of the particles along their travelling direction. Consistently with the above-mentioned symmetries, a region of biaxial strain ( $Q\lt 0$ and $R\gt 0$ ) is found at the upstream end of the particles. This is conveniently visualised in the schematic of figure 10.

Figure 10. Schematic representation of the distributions of $Q$ and $R$ around the particles based on figure 9. The flow, represented by continuous streamlines, is incoming from the left in the particle’s reference frame.

This scenario only partially agrees with the results of Schneiders et al. (Reference Schneiders, Meinke and Schröder2017) for decaying homogeneous isotropic turbulence laden with Kolmogorov-size particles. They also found that particles favour events with axial strain ( $Q\lt 0$ and $R\lt 0$ ) and biaxial strain ( $Q\lt 0$ and $R\gt 0$ ), as shown by the occurrence of a point at the left-Vieillefosse tail and by the more pronounced right-Vieillefosse tail in figure 20 of their paper. However, they did not report an increase of the probability of events with $Q\gt 0$ like in the present case. It is worth mentioning, however, that they do report a global increase of the vortex stretching. A possible explanation of this difference is the lower Reynolds number they considered ( $Re_\lambda \approx 79$ at the initial time) that does not ensure a proper separation of scales, in particular at larger times when turbulence decays. Compared with larger particles with size in the inertial range of scales, the scenario is completely different: in fact, Cannon et al. (2024) found that when large particles are added both $Q$ and $R$ are reduced. Besides energising the small scales, large particles indeed behave also as obstacles for the large flow structures, largely weakening thus the energy content at the large scales. Cannon et al. (Reference Cannon, Olivieri and Rosti2024) only observed an increase of the probability in the strain-dominated region (that we also observe), which is an indication of the intense dissipation regions that arise around the particles.

Figure 11. The $-Q_S$ – $Q_W$ map for the single-phase case (a) and $\varPhi _V=10^{-3}$ and $\rho _p/\rho _f=100$ (b). For both panels the $20$ isolines have a logarithmic spacing between $10^{-7}$ and $10^1$ .

Additional insights are provided by looking separately at the invariants of $S_{ij}$ and $W_{ij}$ (see figure 11). In particular, we consider their second invariants, i.e.

(3.11) \begin{equation} Q_S = - \frac {1}{2} \left ( \alpha ^2 + \beta ^2 + \gamma ^2 \right ) \ \text {and} \ Q_W = \frac {1}{4} \omega ^2. \end{equation}

These invariants are related to the fluid dissipation $\epsilon = -4 \nu Q_S$ and with the fluid enstrophy $\omega ^2 = 4Q_W$ . Therefore, the $Q_S$ – $Q_W$ joint distribution determines whether the flow is dominated by dissipation (extensional-dominated regions with $Q_S\gt Q_W$ ) or by enstrophy (rigid rotation regions with $Q_W\gt Q_S$ ). In shear-dominated regions, dissipation and enstrophy balance and $-Q_S = Q_W$ (Soria et al. Reference Soria, Sondergaard, Cantwell, Chong and Perry1994). For simplicity, we follow Truesdell (Reference Truesdell1954) and introduce $\mathcal {K} = (-Q_W/Q_S)^{1/2}$ ; when $\mathcal {K} = 0$ the flow is extension-dominated, when $\mathcal {K} = \infty$ the flow is dominated by rigid rotation events and when $\mathcal {K} = 1$ the rotation and the stretching are equal, as typical of vortex sheets and shear layers. In the unladen case, events with $Q_W\gt -Q_S$ ( $\mathcal {K} \rightarrow \infty$ ) are more frequent, meaning that for purely Newtonian turbulence the flow is mainly dominated by rigid rotations. In the presence of the particles the scenario slightly changes. A first observation is that the probability of events with large $Q_S$ (or $\mathcal {K} = 0$ ) increases, indicating that the perturbation field induced by these small particles is extensional-dominated. A second observation is that the solid phase favours also events with $-Q_S \approx Q_W$ ( $\mathcal {K} \approx 1$ ), which agrees with the presence of the shear layers induced by the presence of the particles.

4.1. Lagrangian velocity increments

The Lagrangian statistics of the particles’ motion are of fundamental importance in the understanding of transport and mixing. In order to investigate this, we study the Lagrangian velocity increments, defined here as $\delta _\tau u_{p,i} = u_{p,i} (t + \tau ) - u_{p,i}(t)$ , with $u_{p,i}(t)$ being the instantaneous velocity of a particle along direction $i$ at time $t$ . The symmetries of the present problem make the statistics of $\delta _\tau u_{p,i}$ independent of both $t$ and $i$ ; for simplicity hereafter we drop the $i$ index. Figures 12 and 13 describe the particle dynamics at different time scales, by plotting $\delta _\tau u_p$ for different values of the time scale $\tau$ in the $0.2 \tau _\eta \leqslant \tau \leqslant 30 \tau _\eta$ range, where $\tau _\eta = (\nu /\epsilon )^{1/2}$ is the Kolmogorov time scale. For small time scales, the velocity increment provides information about the particle acceleration, i.e. $\delta _\tau u_p \sim a_p\tau$ . A first observation is that the distributions are symmetric, in agreement with the symmetries of the flow. The probability density function of $\delta _\tau u_p$ continuously deforms from the Gaussian at large time scales (see $\tau \approx 30 \tau _\eta$ ) to the development of stretched exponential tails at dissipative time scales (see $\tau \approx 0.2 \tau _\eta$ ), which are the statistical signature of an intermittent Lagrangian dynamics; see Mordant et al. (Reference Mordant, Metz, Michel and Pinton2001), La Porta et al. (Reference La Porta, Voth, Crawford, Alexander and Bodenschatz2001) and Chevillard et al. (Reference Chevillard, Roux, Levêque, Mordant, Pinton and Arneodo2003) for small tracers and Qureshi et al. (Reference Qureshi, Bourgoin, Baudet, Cartellier and Gagne2007) for finite-size neutrally buoyant particles. The wide stretched exponential tails for the smallest $\tau$ show that the finite-size particles with $D_p \approx \eta$ experience very high acceleration events, with a probability which is higher than Gaussian, similarly to what was found for small tracers by La Porta et al. (Reference La Porta, Voth, Crawford, Alexander and Bodenschatz2001) and for finite-size particles by Qureshi et al. (Reference Qureshi, Bourgoin, Baudet, Cartellier and Gagne2007).

Figure 12. Probability density function (PDF) of the Lagrangian velocity increments of the particles $\delta _\tau u_p$ for (a) $\tau = 0.2 \tau _\eta$ and (b) $\tau = 2 \tau _\eta$ . The data are for the case with $\varPhi _V=10^{-3}$ and $\rho _p/\rho _f=5$ and $\rho _p/\rho _f=100$ .

Figure 13. Probability density function (PDF) of the Lagrangian velocity increments of the particles $\delta _\tau u_p$ for (a,b) $\tau = 0.2 \tau _\eta$ , (c,d) $\tau = 2 \tau _\eta$ and (e,f) $\tau = 30 \tau _\eta$ . Plots for (a,c,e) $\rho _p/\rho _f = 5$ and (b,d,f) $\rho _p/\rho _f = 100$ . The black solid lines are the results of the PP-DNS.

We start looking at the influence of $\rho _p/\rho _f$ and $\varPhi _V$ on the Lagrangian intermittency of Kolmogorov-size particles. Figure 12 a shows that heavier particles ( $\rho _p/\rho _f=100$ ) are less likely to experience intermediate values of the acceleration compared to lighter particles ( $\rho _p/\rho _f=5$ ). In contrast, they are more likely to exhibit very low or very intense accelerations. On the one side, the larger inertia of these particles opposes large accelerations and favours small values of $a_p$ . On the other side, heavy particles enhance the flow intermittency (see § 3), promoting extreme events in the flow that are in turn responsible for rare but large particle accelerations. It is worth noticing that for $\varPhi _V=10^{-5}$ the latter effect is barely visible (see figure 13), in agreement with the weak flow modulation shown in § 3. For light particles $\rho _p/\rho _f=5$ figure 13 shows that the distributions of $\delta _\tau u_p$ obtained for $\varPhi _V=10^{-5}$ and $10^{-3}$ overlap almost perfectly, in agreement with the low level of backreaction.

Figure 14. Evolution of the excess kurtosis factor $\mathcal {K}_{\delta _\tau u_p}(\tau )$ for the distributions of the time increments of the particle velocity: (a) $\rho _p/\rho _f=5$ ; (b) $\rho _p/\rho _f=100$ . The open circles in (b) are from the simulation with the coarser grid at $\varPhi _V = 10^{-3}$ and $\rho _p/\rho _f = 100$ , and are shown for validation purpose.

Figure 13(a,c,e) shows that for $\rho _p/\rho _f=5$ the distributions obtained by means of PP-DNS and PR-DNS almost perfectly overlap for all time scales $\tau$ : for light particles the complete MRG equation properly predicts the Lagrangian intermittency of the particles dynamics. For heavier particles the match between the PP-DNS and the PR-DNS is rather good at large time scales, but differences are observed at small $\tau$ , where the PP-DNS does not predict the large tails for $\tau \lessapprox 2 \tau _\eta$ . As discussed above, these extreme events are associated with the flow modulation which is not modelled in our PP-DNS. The comparison between the PP-DNS and PR-DNS results is further detailed in figure 14 where the evolution of the $\delta _\tau u_p$ distribution with $\tau$ is quantified by means of the excess kurtosis $\mathcal {K}_{\delta _\tau u_p}(\tau ) = \langle { \delta _\tau u_p^4} \rangle \!/ \langle { \delta _\tau u_p^2 } \rangle \!^2 -3$ . At large scales $\mathcal {K}_{\delta _\tau u_p} \approx 0$ in agreement with the Gaussian-like shape of the distribution, while it steeply increases at small scales. For $\rho _p/\rho _f=5$ the good agreement between the PP-DNS and the PR-DNS is again clear, with a small deviation for the $\varPhi _V=10^{-3}$ case, which is due to the non-zero flow modulation. For $\rho _p/\rho _f=100$ , instead, the agreement is good at large scales, while the three curves substantially deviate for small $\tau$ , accordingly with the larger tails of $\delta _\tau u_p$ found in the PR-DNS.

4.2. The particle-velocity structure function

We now consider the statistics of the particle–particle relative velocity $\delta \boldsymbol {u}_p = \boldsymbol {u}_p(\boldsymbol {x}_{p,j}(t),t) - \boldsymbol {u}_p(\boldsymbol {x}_{p,i}(t),t)$ , where $\boldsymbol {x}_{p,i}(t)$ and $\boldsymbol {x}_{p,j}(t)$ denote the position of any two particles $i$ and $j$ at time $t$ . The distribution of $\delta \boldsymbol {u}_p$ across all particle couples plays a key role in several theories regarding the tendency of particles to form clusters (see e.g. Gustavsson & Mehlig Reference Gustavsson and Mehlig2011; Bragg & Collins Reference Bragg and Collins2014; Bragg et al. Reference Bragg, Ireland and Collins2015).

Figure 15. Second-order structure function based on the particle velocity field $S_{2,p}$ : (a) $\rho _p/\rho _f=5$ ; (b) $\rho _p/\rho _f = 100$ . The open circles in (b) are from the simulation with the coarser grid at $\varPhi _V = 10^{-3}$ and $\rho _p/\rho _f = 100$ , and are shown for validation purpose. The inset in (b) shows the compensated $S_{2,p}r^{-2/3}$ in the $10 \leqslant r/\eta \leqslant 50$ range.

Figure 15 plots the second-order structure function of the particle velocity, i.e.

(4.1) \begin{equation} S_{2,p}(r) = \! \left \langle { \left ( \delta \boldsymbol {u}_p(\boldsymbol {r}) \cdot \frac {\boldsymbol {r}}{r} \right )^2 } \right \rangle \!, \end{equation}

where $\boldsymbol {r}$ is the separation vector between particle $i$ and $j$ , and $r = |\boldsymbol {r}|$ . Figure 15(a) is for $\rho _p/\rho _f=5$ , while figure 15(b) is for $\rho _p/\rho _f=100$ .

For $\rho _p/\rho _f=5$ , $S_{2,p}$ resembles the fluid second-order structure function $S_2$ (see figure 3 a): light particles have small inertia and closely follow the fluid motion. Function $S_{2,p}$ exhibits the $r^2$ scaling at small scales and the $r^{2/3}$ scaling predicted by the Kolmogorov theory in the inertial range of scales. In agreement with the negligible flow modulation, the results of the PP-DNS match almost perfectly those of the PR-DNS for these parameters.

Figure 15(b) deals with the $\rho _p/\rho _f=100$ cases. A first observation is that the results from the PR-DNS with $\varPhi _V=10^{-5}$ and $\varPhi _V=10^{-3}$ do not collapse; this is consistent with the larger flow modulation observed for the larger volume fraction, and agrees with the above discussed results for the single-particle statistics. Notably, for heavy particles $S_{2,p}$ differs from the fluid structure function $S_2$ at small scales. According to both the PP-DNS and the PR-DNS, $S_{2,p}$ does not exhibit an $r^2$ scaling at the smallest scales, being substantially flat at small $r$ . The relative motion between couples of heavy particles placed at a small distances $r$ is substantially uncorrelated as well as decoupled from the small-scale fluid motion due to their large inertia. Note that the absence of the $S_{2,p} \sim r^2$ scaling indicates that at small scales the Eulerian particle velocity field cannot be described with a Taylor expansion. Notably, figure 15 shows that $S_{2,p}$ recovers the $S_2$ slope at larger $r$ , exhibiting the classical Kolmogorov $r^{2/3}$ scaling in the inertial range of scales. This suggests that, despite the large inertia, the relative particle–particle velocity $\delta \boldsymbol {u}_p$ between two particles is driven by turbulent eddies having size comparable to $r$ , provided that $r$ is large enough. When comparing the results of the PP-DNS with those of the PR-DNS, we note that the slope of $S_{2,p}$ matches for small ( $ r/\eta \lessapprox 5$ ) and large ( $r \gtrapprox 30$ ) scales. For intermediate scales $10 \lessapprox r/\eta \lessapprox 30$ , instead, the PR-DNS predicts a steeper slope for both $\Phi _V=10^{-5}$ and $\Phi _V=10^{-3}$ (see the inset in figure 15 b). The finite size of the particles does not influence $S_{2,p}$ for large and small scales where $r/D_p = \mathcal {O}(100)$ and $r/D_p =\mathcal {O}(1)$ , but it does for intermediate scales $r/D_p = \mathcal {O}(10)$ .

Figure 16. Probability density function (PDF) of the radial particle–particle relative velocity $\delta \boldsymbol {u}_p \cdot \boldsymbol {r}/r$ for $\rho _p/\rho _f=100$ : (a) $r \approx 3.5D_p = 0.0627$ , (b) $r \approx 7D_p =0.1254$ , (c) $r \approx 21D_p = 0.3761$ and (d) $r \approx 67D_p = 1.2$ .

Figure 16 sheds further light on the relative particle–particle velocity by plotting the distribution of $\delta \boldsymbol {u}_p \cdot \boldsymbol {r}/r$ , i.e. the component of $\delta \boldsymbol {u}_p$ projected along the vector separating the two particles, for different values of $r$ . When $\delta \boldsymbol {u}_p \cdot \boldsymbol {r}/r\gt 0$ , the two particles depart, whereas they get closer when $\delta \boldsymbol {u}_p \cdot \boldsymbol {r}/r\lt 0$ . We consider the case with $\rho _p/\rho _f = 100$ to provide further insights of the distribution of $S_{2,p}$ , shown in figure 15(b). A first observation is that, similarly to that found for larger particles by Chiarini & Rosti (Reference Chiarini and Rosti2024), the distribution of $\delta \boldsymbol {u}_p \cdot \boldsymbol {r}/r$ is left skewed, with a slightly positive mode and a long negative tail. The distribution becomes progressively more flat when increasing $r$ , in agreement with a lower level of the correlation between the velocity of the two particles. When comparing the distributions for $\varPhi _V = 10^{-5}$ and $\varPhi _V = 10^{-3}$ , figure 16 shows that for all $r$ the tails are shorter for the larger $\varPhi _V$ , with the difference decreasing as $r$ increases. This is consistent with the stronger flow modulation that globally leads to a weaker level of flow fluctuations (see table 1). Also, in agreement with the evolution of $S_{2,p}$ with $r$ (see figure 15), for $\varPhi _V=10^{-5}$ the distribution obtained with the PP-DNS collapses nicely with that obtained with the PR-DNS at small scales (see figure 16 a,b), with some substantial differences arising for intermediate scales when the finite-size effects are relevant.

5.1. Clustering

Over the years, several tools have been used to characterise the spatial arrangement of particles in a flow (for an overview, see Brandt & Coletti (Reference Brandt and Coletti2022)). We use the Voronoï tessellation, which has been extensively used by several authors (see e.g. Monchaux et al. Reference Monchaux, Bourgoin and Cartellier2010, Reference Monchaux, Bourgoin and Cartellier2012). The position of each particle is identified with its centre and the computational domain is divided in subdomains, such that each grid cell is associated with the closest particle. The Voronoï volume $V_V$ of each particle is thus defined as the collective volume of grid cells that are closer to it than to other particles. The inverse of the Voronoï volume provides a measure of the local concentration: particles placed in void regions possess a large Voronoï volume, while particles that are part of a cluster have a small Voronoï volume. Based on these observations, the intensity of the clustering of a suspension can be measured by comparing the distribution of its Voronoï volumes with that of a control consisting of an equivalent, uniformly random suspension of particles. More intense clustering leads to a variance of the distribution of $V_V$ that is larger than that of the control. In the context of PR-DNS, the overlap between particles is not allowed also in the reference random arrangement.

Figure 17. Comparison of the variance of the Voronoï volumes obtained with PR-DNS and PP-DNS: (a) $\varPhi _V=10^{-5}$ ; (b) $\varPhi _V = 10^{-3}$ .

Figure 17 presents the clustering intensity for the different values of $\varPhi _V$ and $\rho _p/\rho _f$ considered. A first observation is that the PP-DNS underestimate the level of clustering for all cases. Our computations show that the discrepancy between the PP-DNS and the PR-DNS increases with $\rho _p/\rho _f$ and/or $\varPhi _V$ (see § 5.2 for further details).

We now move to the effect of the volume fraction $\varPhi _V$ and of the particle density $\rho _p/\rho _f$ . As expected, the level of clustering increases with $\varPhi _V$ . When fixing $\varPhi _V$ , instead, figure 17 shows that heavier particles with $\rho _p/\rho _f = 100$ cluster more than lighter particles with $\rho _p/\rho _f =5$ . For light particles, the low level of clustering $\sigma /\sigma _{rand} \approx 1$ agrees with the previous results of Fiabane et al. (Reference Fiabane, Zimmermann, Volk, Pinton and Bourgoin2012),Uhlmann & Chouippe (Reference Uhlmann and Chouippe2017) and Chiarini & Rosti (Reference Chiarini and Rosti2024), who considered larger particles $ 5 \leqslant D_p/\eta \leqslant 123$ over a wide range of Reynolds numbers $105 \leqslant Re_\lambda \leqslant 430$ . Light particles have small inertia and are less likely to drift from the trajectories of the fluid elements. In contrast, the larger level of clustering observed when increasing the particle density from $\rho _p/\rho _f = 5$ to $\rho _p/\rho _f=100$ is not consistent with that found for larger particles. For suspensions of particles with size in the inertial range of scales, indeed, Chiarini & Rosti (Reference Chiarini and Rosti2024) found that the level of clustering exhibits a non-monotonic dependence on $\rho _p/\rho _f$ , with the maximum occurring at intermediate densities (see figure 24 of their paper), and the minimum being at the largest density ratio they considered, i.e. $\rho _p/\rho _f = 100$ . However, they found that in the $D_p{-}\rho _p$ space of parameters the maximum level of clustering moves towards larger $\rho _p$ as $D_p$ decreases, suggesting that the tendency of particles to cluster is driven by the Stokes number of the particles, rather than by their density. Accordingly, their data show that the level of clustering is maximum when $St = \mathcal {O}(1{-}10)$ . This agrees with the early works of Wang & Maxey (Reference Wang and Maxey1993),Fessler et al. (Reference Fessler, Kulick and Eaton1994) and Aliseda et al. (Reference Aliseda, Cartellier, Hainaux and Lasheras2002), and it is consistent with our results (see table 1 for the particles’ Stokes number).

Figure 18. Probability density function (PDF) of the Voronoï volumes for (a) $\varPhi _V = 10^{-5}$ and (b) $\varPhi _V = 10^{-3}$ for all $\rho _p/\rho _f$ . In (b), the dashed lines are used for the PP-DNS. The black lines refer to the Voronoï volumes for a random distribution of particles. Note that particles are not overlapping for the PR-DNS also in the random distribution.

The complete distributions of the Voronoï volumes are provided in figure 18. Compared with the corresponding random arrangement of particles, the tails of the $V_V$ distribution are longer, and grow ever more so with increasing $\varPhi _V$ and/or $\rho _p/\rho _f$ . This is in agreement with the above discussion, since stronger clustering corresponds to a larger number of small and large Voronoï volumes. Note that the PP-DNS underestimation of the level of clustering is visualised in figure 18 with the shorter tails. Figure 18 can be used to determine which particles are part of clusters and which are part of void regions (Monchaux et al. Reference Monchaux, Bourgoin and Cartellier2010). This information is used in § 5.2, when discussing particle preferential sampling. In the presence of clusters, two cross-over points arise between the $V_V$ distributions of the actual suspension and that of the corresponding random arrangement of particles. Particles with a Voronoï volume smaller than the left cross-over point $V_{th,l}$ are part of a cluster, while those with a Voronoï volume larger than the right cross-over point $V_{th,r}$ are part of void regions. Particles that are part of a cluster and have Voronoï volumes that share at least one vertex are part of the same cluster. Note that, as the level of clustering increases, the threshold of the Voronoï volume that delimits the particles entrapped in clusters decreases.

A different type of information regarding the spatial arrangement of the particles can be provided by means of the radial distribution function $g(r)$ , also referred to as pair correlation function (see figure 19). It describes how the particles’ density varies as a function of the distance away from a reference particle. In other words, it is a measure of the probability of finding particles at a distance $r$ relative to that of a homogeneous distribution. Following Salazar et al. (Reference Salazar, Jong, Cao, Woodward, Meng and Collins2008) and Saw et al. (Reference Saw, Shaw, Ayyalasomayajula, Chuang and Gylfason2008), the radial distribution function is defined as

(5.1) \begin{equation} g(r) = \frac { N_s(r)/\Delta V(r) }{ N_{pa}/V }, \end{equation}

where $N_s(r)$ is the number of particle pairs separated by a distance between $r-\Delta r$ and $r+\Delta r$ , $\Delta V(r)$ is the volume of the spherical shell of inner and outer radius $r-\Delta r$ and $r+\Delta r$ , respectively, $N_{pa}$ is the total number of particle pairs present in the system, $N_{pa} = N(N-1)/2$ , and $V$ is the volume of the computational domain. In a uniform distribution where overlapping between particles is allowed, $g(r)=1$ for all $r$ .

Figure 19. Radial distribution function for (a) $\varPhi _V = 10^{-5}$ and (b) $\varPhi _V = 10^{-3}$ . Solid lines are for PR-DNS, while dashed lines are for PP-DNS.

The radial distribution function (see figure 19) shows that for all cases the accumulation is maximum at the smallest distances. Note that the maximum of $g(r)$ occurs at $r \approx D_p$ for PR-DNS, as the overlap between particles is not allowed. In agreement with the above discussion, the heavy particles with $\rho _p/\rho _f = 100$ show a larger level of clustering compared with the lighter ones with $\rho _p/\rho _f = 5$ . The level of accumulation is also slightly larger for $\varPhi _V = 10^{-3}$ at all $r$ . Similarly to what observed with the Voronoï tessellation, figure 19 shows that the PP-DNS underpredicts the level of particle accumulation. As clearly visible for $\varPhi _V = 10^{-3}$ and $\rho _p/\rho _f = 100$ , the discrepancy between the PP-DNS and PR-DNS is maximum at the smallest separations.

5.2. Preferential sampling

In the previous section we have shown that the solid phase is not homogeneously distributed in space, and that the particles exhibit a mild level of clustering. In this section we relate the presence of the clusters to the tendency of the particles to preferentially sample certain regions of the flow. In doing this, we also provide a possible explanation of the different level of clustering predicted by the PR-DNS and PP-DNS for the $\varPhi _V = 10^{-3}$ and $\rho _p/\rho _f = 100$ case. Over the years several mechanisms have been proposed as governing the particles’ preferential sampling of the flow, most of them justified using the MRG equation and thus, strictly speaking, valid only in the context of sub-Kolmogorov particles. In the following we use the centrifuge mechanism (Maxey Reference Maxey1987) to explain the tendency of Kolmogorov-size particles to form clusters.

In the limit where the point-particle approximation holds, Maxey (Reference Maxey1987) has shown that particles with large $\rho _p/\rho _f$ tend to collect in regions of high strain rate and low vorticity. In the presence of a vortex, indeed, heavy particles cannot follow the flow streamlines because of their large inertia, and tend to drift from the vortex core. Similarly, in the case of a pure straining flow, heavy particles drift towards the stagnation point at the centre. We quantify the tendency of particles to sample regions of high strain by using the second invariant of the deformation tensor (see § 3.4), i.e. $Q = - S_{ij}S_{ij}/2 + \omega ^2/4$ ; recall that regions where $Q$ is large and positive are regions of high vorticity ( $Q \sim \omega ^2/4$ ), while regions where $Q$ is large and negative are regions of high strain ( $ Q \sim - S_{ij}S_{ij}$ ). We investigate the particle preferential sampling by computing the probability density function of $Q$ at the particles position (Squires & Eaton Reference Squires and Eaton1990). For the sake of brevity, in this section we limit the investigation to the largest volume fraction $\varPhi _V = 10^{-3}$ . For PP-DNS the value of $Q$ at each particle position is obtained after linear interpolation. For PR-DNS, instead, the value of $Q$ seen by each particle is estimated as the average value within a spherical shell centred with the particle and having a radius $R_{sh}\gt R_p$ , where $R_p=D_p/2$ is the radius of the particles. It is important to note, however, that due to the particles’ backreaction, in the case of PR-DNS, the value of $Q$ seen by each particle is actually the result of three different effects: (i) the larger-scale flow properties of the region that the particle is sampling, (ii) the smaller-scale influence of the particle on the surrounding flow and (iii) the effect of nearby particles on the flow. Also, a suitable choice of $R_{sh}$ should be made: when $R_{sh}$ is too small, only the influence of the particle on the surrounding flow is considered (see e.g. Kidanemariam et al. Reference Kidanemariam, Chan-Braun, Doychev and Uhlmann2013; Oka & Goto Reference Oka and Goto2022), while when $R_{sh}$ is too large, spurious contributions that do not affect the particle location are instead captured. In order to obtain a complete picture, we have tested different values of $R_{sh}$ .

Figure 20. Probability density function (PDF) of the second invariant $Q$ of the velocity gradient tensor evaluated at the particle position: (a) $\varPhi _V=10^{-3}$ and $\rho _p/\rho _f=5$ ; (b) $\varPhi _V=10^{-3}$ and $\rho _p/\rho _f = 100$ . Here $R_{sh}$ indicates the radius of the spherical shell used to estimate the value of $Q$ seen by the particles in the PR-DNS. The black line indicates the distribution according to the PP-DNS.

Figure 20 shows the probability density function of $Q_\ell$ , i.e. $Q$ computed at the particle position. We start by looking at the dependence of the PR-DNS results on the radius $R_{sh}$ of the shell. For $\rho _p/\rho _f = 5$ , the curves computed for values of $R_{sh}$ in the range $2 \leqslant R_{sh}/R_p \leqslant 7$ show an almost perfect overlap. This is consistent with the low level of modulation discussed in § 3, and indicates that for these light particles the main contribution to $Q_\ell$ comes from the larger-scale properties of the flow region sampled by the particles (note that due to the low volume fraction the influence of the particle–particle interaction is negligible). For heavy particles with $\rho _p/\rho _f = 100$ , instead, the distribution of $Q_\ell$ largely varies with $R_{sh}$ : as $R_{sh}$ decreases, the left tail of the distribution becomes longer, meaning that particles are more likely to see large negative values of $Q$ . These large negative values of $Q$ are, at least partially, the result of the influence of the particles on the neighbouring flow; see the $Q{-}R$ map in figure 8. Note that, for $R_{sh}/R_p \gtrapprox 5$ , the distribution of $Q_\ell$ shows only marginal variations, as for these $R_{sh}$ the large-scale flow contribution dominates. This suggests that for $\rho _p/\rho _f = 100$ the influence of particles on the surrounding flow extends for less than $5R_p$ . Overall, for both light and heavy particles the distribution of $Q_\ell$ is left-skewed and shows an almost null probability of positive $Q_\ell$ . For the present parameters, both PP-DNS and PR-DNS give evidence that Kolmogorov-size particles preferentially sample regions of high strain rate. This is also visible in the instantaneous snapshot shown in figure 1, with particles sampling regions with low $\omega ^2$ . In other words, in the context of Kolmogorov-size particles the formation of clusters is, at least partially, governed by the centrifuge mechanism.

Let us now focus on the differences between the PR-DNS and PP-DNS results. For particles with $\rho _p/\rho _f=5$ , the $Q_\ell$ distribution obtained with PP-DNS overlaps almost perfectly with that obtained with PR-DNS; the point-particle approximation predicts fairly well the tendency of light particles to sample the $Q\lt 0$ regions of the flow. Note that this is consistent with the good agreement found in figure 18 when discussing the distribution of the Voronoï volumes. For $\rho _p/\rho _f = 100$ , instead, the $Q_\ell$ distribution obtained with the PP-DNS significantly deviates from that obtained with the PR-DNS. For all $R_{sh}$ , the $Q_\ell$ distribution obtained with PR-DNS shows a shorter right tail and predicts a higher probability of negative $Q$ : finite-size heavy particles are less (more) likely to see positive (negative) values of $Q$ . Based on this, one may conclude that, for the present parameters, the PP-DNS underestimates the tendency of the particles to preferentially sample regions of high strain, and this may explain the larger level of clustering predicted by the PR-DNS for the $\varPhi _V = 10^{-3}$ and $\rho _p/\rho _f = 100$ case (see figure 17). We speculate that the discrepancy between the PR-DNS and the PP-DNS observed for the $\varPhi _V=10^{-3}$ and $\rho _p/\rho _f=100$ case is due to the flow modulation, rather than to the particles’ finite-size effects. Unlike for the light particles ( $\rho _p/\rho _f=5$ ), indeed, for the present parameters the influence of the heavy particles ( $\rho _p/\rho _f = 100$ ) on the fluid phase is not negligible (see § 3).

A last comment regards the influence of $\rho _p/\rho _f$ on the distribution of $Q_\ell$ . According to both PP-DNS and PR-DNS, heavier particles exhibit a larger tendency to sample regions with more negative $Q$ , as shown by the left tail being longer for the $\rho _p/\rho _f=100$ case. Due to their larger inertia, indeed, heavier particles enhance the centrifuge mechanism, being more likely to drift from the high-vorticity regions of the flow.

Figure 21. Joint probability density function of $Q$ and $V_V$ for (a,b) $\varPhi _V = 10^{-3}$ and for (c,d) $\rho _p/\rho _f=5$ and $\rho _p/\rho _f=100$ . (a,c) The PR-DNS and (b,d) the PP-DNS. White/black colour denotes minimum/maximum probability. The dashed lines represent $V_{th,l}$ , while the dash-dotted lines $V_{th,r}$ .

To provide additional insights regarding the relation between the presence of clusters and the particle preferential sampling, figure 21 shows the joint probability density function of $Q_\ell$ and $V_V$ for $\varPhi _V = 10^{-3}$ . Based on the above discussion, here we set $R_{sh}/R_p = 3$ for the computation of $Q_\ell$ , since it is large enough to account for the particle preferential location and small enough to avoid spurious contributions. We recall that according to Monchaux et al. (Reference Monchaux, Bourgoin and Cartellier2010), particles with $V_V \leqslant V_{th,l}$ are part of a cluster, while particles with $V_V \geqslant V_{th,r}$ are in void regions of the flow. For all cases, the most negative values of $Q_\ell$ well correlate with small and intermediate Voronoï volumes, with $V_V \lessapprox V_{th,r}$ . Particles that are in void regions and are not part of a cluster are less likely to see large negative values of $Q$ . This agrees with the above observation that the tendency of Kolmogorov-size particles to form clusters is governed by the centrifuge mechanism. We now focus on $\rho _p/\rho _f=100$ (see figure 21 c,d). The joint distribution shows that the larger probability of negative $Q_\ell$ predicted by the PR-DNS is concentrated at the smallest Voronoï volumes with $V_V \lessapprox V_{th,l}$ . Again, this shows that the higher level of clustering detected in this case with the PR-DNS well correlates with the greater tendency of finite-size particles to sample regions of the flow with intense strain.