Finite-size inertial spherical particles in turbulence
Alessandro Chiarini Marco Edoardo Rosti

A Novel Channel Coding Scheme 
 for Digital Multiple Access Computing




    Xiaojing Yan, Saeed Razavikia, Carlo Fischione

	School of Electrical Engineering and Computer Science, KTH Royal Institute of Technology, Stockholm, Sweden

	Email: {xiay, sraz, carlofi}@kth.se

    Accepted XXX, Received YYY; in original form ZZZ
=======================================================================================================================================================================================================





We investigate by direct numerical simulations the fluid-solid interaction of non-dilute suspensions of spherical particles moving in triperiodic turbulence, at the relatively large Reynolds number of Re_λ≈ 400. The solid-to-fluid density ratio is varied between 1.3 and 100, the particle diameter D ranges between 16 ≤ D/η≤ 123 (η is the Kolmogorov scale), and the volume fraction of the suspension is 0.079. Turbulence is sustained using the Arnold-Beltrami-Childress cellular-flow forcing.  
The influence of the solid phase on the largest and energetic scales of the flow changes with the size and density of the particles. Light and large particles modulate all scales in a isotropic way, while heavier and smaller particles modulate the largest scales of the flow towards an anisotropic state. Smaller scales are isotropic and homogeneous for all cases. The mechanism driving the energy transfer across scales changes with the size and the density of the particles. For large and light particles the energy transfer is only marginally influenced by the fluid-solid interaction. For small and heavy particles, instead, the classical energy cascade is subdominant at all scales, and the energy transfer is essentially driven by the fluid-solid coupling. The influence of the solid phase on the flow intermittency is also discussed. Besides, the collective motion of the particles and their preferential location in relation with properties of the carrier flow are analysed. The solid phase exhibits moderate clustering; for large particles the level of clustering decreases with their density, while for small particles it is maximum for intermediate values. 







§ INTRODUCTION



Particle-laden turbulent flows have attracted the attention of many scholars over the last decades because of their relevance that goes beyond a fundamental interest, and encompasses several applications in natural, industrial and geophysical fields
<cit.>.

The inclusion of solid particles in the flow introduces additional momentum in the suspension that may result in modulation of the flow structures. When the suspension is dilute enough, the carrier flow is almost not altered by the presence of the particles. When the suspension is non-dilute, instead, the carrier flow undergoes macroscopic changes, and the particle-fluid interaction can not be neglected <cit.>.

The rheological properties of suspensions have been first studied in the limit of viscous Stokesian regime, where the inertial effects are negligible. The presence of the particles affects the deformation of the surrounding fluid, and the effective viscosity μ_e of the particle-fluid mixture depends on the dynamics of the dispersed phase. For a dilute suspension with negligible inter-particle interactions, <cit.> derived an analytical linear relation for the effective viscosity, i.e. μ_e = μ(1+2.5Φ_V), with Φ_V = V_p/(V_p + V_f) being the volume fraction of the dispersed phase; V_p and V_f are the volumes of the solid and fluid phases. Later, a quadratic correction was proposed by <cit.> and <cit.> to account for mutual particle interaction for slightly higher volume fractions. For higher concentrations, where the inter-particle interactions play a crucial role, only semi-empirical formulas exist <cit.>.

In the turbulent regime the particle size has been found to play a crucial role in determining the effect of the dispersed phase on the carrier flow. <cit.> considered together data from particulate turbulent pipe flows and jets, and found that the turbulence modulation depends on the ratio of the particle diameter D and the flow integral scale ℒ. When D/ℒ>0.1 the turbulent intensity of the carrier phase increases with respect to the single phase case, while when D/ℒ<0.1 it decreases. Large particles produce fluctuations in their wake and enhance the turbulent activity, while small ones drain energy from the large-scale turbulent eddies of the flow. Indeed, <cit.> experimentally considered an air-solid two-phase flow in a horizontal pipe with diameter of 30 mm, and investigated the influence of small and large particles with diameter of 0.2 mm and 3.4 mm on the carrier flow, varying the section average air velocity between 6 and 20 m/s. They observed that large particles markedly increase the turbulent activity, while small ones reduce it. These results were later confirmed by the same authors, with experiments of an air-solid two-phase flow in a vertical pipe <cit.>. More recently, <cit.> studied the turbulence modulation of a nearly isotropic flow field due to the presence of single glass particles with varying diameter, and observed that the critical ratio of the particle size to integral scale of D/ℒ = 0.41 separates the regimes of attenuation and enhancement of turbulent intensity.

Large steps towards a complete description of the interaction between rigid particles of different size and density and the fluid phase have been performed over the last years thanks to the increase in the supercomputers' power, which enables the study of such problem via direct numerical simulations (DNS) <cit.>. One of the most common approaches used today to model many particle-laden flows is based on the point-particle approximation. In general, each particle is treated as a mathematical point source of mass, momentum and energy, with the particles assumed to be much smaller than any structure of the flow. The point-particle approximation is valid in the limit in which the volume of each particle is negligible, i.e. the case with vanishing disperse-phase volume fraction. Indeed, the point-particle assumption requires that the fluid velocity field does not display a turbulent behaviour at the scale of the particle, meaning that D ≤η (η being the Kolmogorov scale).  
According to this method, the fluid-solid interaction is completely described by a forcing term that should be accurately modelled; see for example <cit.>. The theoretical developments of the point-particle approximation, however, are limited to a relatively narrow range of parameters. For example, investigating the motion of a single particle in homogeneous isotropic turbulence, <cit.> found that the particle velocity variance reflects what predicted by the point-particle approximation for D/ η≤ 4 only. Within the point-particle approximation, <cit.> and <cit.> showed the importance of the particle Stokes number in determining the turbulence modulation. <cit.> tested the point-particle approximation in a turbulent channel flow laden with small inertial particle, with high particle-to-fluid density ratio. They considered a volume fraction of Φ_V ≈ 10^-5 to ensure that the particle feedback on the flow is negligible, and discussed the validity of models that approximate the particle dynamics considering only the inertial and nonlinear drag forces. They observed that in the bulk of the channel these models predict pretty well the statistics of the particle velocity. Close to the wall, however, the models fail as they are not able to capture the shear-induced lift force acting on the particles in this region, which instead is well predicted by the lift force model introduced by <cit.>.

A complete understanding of the fluid-solid interaction problem over a wider range of parameters <cit.> requires to properly resolve the flow around each particle, without the need of relying on models. In this respect, over the last years, several numerical methods based on the coupling between direct numerical simulations and the immerse boundary method (IBM) have appeared, and have been used to study particulate turbulence; see for example the numerical methods described in <cit.>. Based on these developments, several works have carried out DNS of particulate turbulence, but mainly considering small Reynolds numbers and/or large particles, due to the extremely large computational cost; see for example <cit.> for homogeneous isotropic turbulence, <cit.> for channel flow, <cit.> for duct flow and <cit.> for pipe flow.

In this work, we consider non-dilute suspensions of finite-size spherical particles in periodic turbulence, and investigate the fluid-solid interaction in the two-dimensional parameter space of particle size and density. This flow configuration has been investigated by several authors over the years. <cit.> investigated the influence of suspensions of finite-size particles with a solid-to-fluid density ratio of ρ_p/ρ_f ≈ 1.5 and Φ_V ≈ 0.02-0.1 on periodic turbulence at a microscale Reynolds number of Re_λ = u' λ / ν = 61 (u' is the average velocity fluctuation, λ is the Taylor length-scale and ν is the fluid kinematic viscosity). They found that the energy spectrum is enhanced for wavenumbers κ>κ_p≈ 0.72 κ_d, where κ_d = 2 π /D is the wavenumber corresponding to the particle diameter, and that it is attenuated for κ<κ_p. The particles drain energy from the large scales of the flow, and inject it at smaller scales by means of their wake. <cit.> experimentally studied the influence of heavy particles with a diameter similar to the Kolmogorov scale on homogeneous isotropic turbulence, setting the Reynolds number at Re_λ =230. They observed that, the fluid turbulent kinetic energy and the viscous dissipation decrease when increasing the mass loading. <cit.> numerically confirmed the results by <cit.> at Re_λ≈ 60 with a particle-to-fluid density ratio of ρ_p/ρ_f = 1.4 and a volume fraction of Φ_V = 0.06. <cit.> studied the influence of particles of Taylor length-scale on decaying isotropic turbulence. They found that, in contrast to what happens when using particles with diameter smaller than the Kolmogorov scale, the turbulent kinetic energy is always smaller than that of the single-phase flow, and that the two-way coupling rate of change is always positive. Later, using the same framework, <cit.> found that the Stokes number should not be used as an indicator to estimate at which extent particles of the Taylor length-scale modulate the turbulent carrier flow, as particles with same response time but different diameter or density may have different effect on the flow, unlike what observed for sub-Kolmogorov particles <cit.>. <cit.> investigated the influence of spherical particles with diameter of approximately 5 and 8 times the Kolmogorov length and particle-to-fluid density ratio of ρ_p/ρ_f=1.5 on forced homogeneous isotropic turbulence at Re_λ≈ 130. They observed that the disperse phase exhibits clustering with moderate intensity, and that this clustering decreases with the particle diameter. They also suggest that small and light finite-size particles follow the expansive directions of fluid acceleration field, as happens also for sub-Kolmogorov particles <cit.>. 
More recently, <cit.> considered homogeneous isotropic turbulence with a well developed inertial range of scales (they set the Reynolds number at Re_λ≈ 400), and investigated the turbulence modulation by non-dilute suspensions (Φ_V=0.079) of spherical particles with size that lies in the inertial range (D/η=123), varying the particle-to-fluid density ratio between 1.3 ≤ρ_p/ρ_f ≤ 100. Besides the energy depletion at large scales and the energy enhancement at the small ones, they characterise how particles modify the classical energy cascade described by Richardson and Kolmogorov. They observed that the energy transfer is governed by the fluid-solid coupling at scales larger than the particle diameter, and that the classical energy cascade is recovered only at smaller scales.
<cit.> performed a parametric DNS study of particulate homogeneous isotropic turbulence, and investigated the turbulence modulation due to spherical solid particles with different Stokes number at a volume fraction of Φ_V = 8.1 × 10^-3 and reference Reynolds number of Re_λ≤ 100. They varied the particle diameter in the 7.8 ≤ D/η≤ 64 range, and found that the turbulent kinetic energy decreases when reducing D, due to the additional energy dissipation rate in the wake. 
Similarly, <cit.> used a multiple-relaxation-time lattice Boltzman model, and studied the influence of the particle-to-fluid density ratio and particle diameter on the turbulence modulation at Re_λ≈ 75, varying the parameters between 8 ≤ D/η≤ 16 and 5 ≤ρ_p/ρ_f ≤ 20. Besides confirming the turbulent kinetic attenuation for smaller and heavier particles, they observed that the region affected by the particles depends on their diameter, and that the dissipation close to their surface increases with the particle density due to the larger slip velocity and particle Reynolds number.



Despite the large number of studies in the literature, further investigations are necessary for a complete understanding of particle-laden turbulent flows, even in the simple framework of a triperiodic box. Most of the numerical works, indeed, have considered relatively low Reynolds numbers, where the inertial range of turbulence is not completely developed, and relatively low volume fractions, that lead to a rather weak flow modulation (see figure <ref>). As reported in <cit.>, reliable studies at larger Reynolds numbers and at higher concentrations are needed, to bridge the gap between the studies of small, heavy particles and large, weakly buoyant particles, and to promote the development of models that go beyond the point-particle approximation. In this work, we do a step in this direction. By means of a massive parametric study based on DNS and IBM, we investigate the fluid-solid interaction in a triperiodic turbulent flow laden with spherical particles with size that lies within the inertial range. We consider a Reynolds number of Re_λ≈ 400, which ensures an extensive inertial range of scales, and a volume fraction of Φ_V = 0.079, which is large enough for the suspension to be non dilute, but small enough for the particle-particle interactions to be subdominant; see figure <ref>. The fluid-solid interaction is studied in the two-dimensional parameter space of particle size D and particle-to-fluid density ratio ρ_p/ρ_f. The particle size is varied between 16 ≤ D/η≤ 123, while the particle density is varied between 1.3 ≤ρ_p/ρ_f ≤ 100. The specific goals of the paper are: (i) to provide an extensive characterisation of the carrier flow modulation in the D-ρ_p/ρ_f space, extending the previous works by <cit.> and <cit.> to a wider range of parameters and to a larger Reynolds number; (ii) to describe the near-particle flow modulation and provide insights on how the fluid-solid interaction changes with D and ρ_p/ρ_f, which is relevant for modelling purposes; (iii) to investigate the collective motion and preferential location of the particles, extending the work by <cit.> to larger and heavier particles and to a larger Reynolds number.

The paper is structured as follows. In section <ref>, the physical model and the numerical method are briefly presented. Section <ref> describes the influence of the solid phase on the carrier flow. In section <ref>, the effect of the particle size and density on the near-particle flow modulation is addressed. The dynamics of the particles and their collective motion are then discussed in section <ref>. Eventually, in section <ref> concluding remarks are provided.




§ NUMERICAL METHODS





The fluid-solid interaction of non-dilute suspension of solid spherical particles is investigated by means of a massive parametric study based on direct numerical simulations. We consider an ensemble of N finite-size spherical particles with size D and density ρ_p, suspended in a periodic turbulent flow; see figure <ref>. Particles are released in a cubic domain of size L=2π, having periodic boundary conditions in all directions, where turbulence is generated and sustained using the Arnold-Beltrami-Childress (ABC) cellular-flow forcing <cit.>. 

The carrier flow is governed by the incompressible Navier–Stokes equations for a Newtonian fluid, i.e.


    ∂u/∂ t + ∇·uu = -1/ρ_f∇ p + ν∇^2 u  + f + f^← p
    ∇·u = 0.


Here, u=(u,v,w) is the velocity vector field, p is the reduced pressure, and ρ_f and ν are the fluid density and kinematic viscosity. In the momentum equation, the term f^← p is the force due to the solid phase, and f is the external ABC body force used to sustain turbulence

    f=(2π/L)^2 F_o
    [ A sin(2π/Lz) + C cos(2π/Ly); B sin(2π/Lx) + A cos(2π/Lz); C sin(2π/Ly) + B cos(2π/Lx);                             ],


where A=B=C=1, and F_o is a constant parameter. 

The dynamics of the spherical particles is governed by the Euler-Newton equations,


    m du_p/d t = f^← f + f^↔ p
    
    I dω_p/d t = L_p^← f,


where u_p and ω_p are the velocity and angular velocity of the particles, and m=πρ_p D^3/6 and I = m D^2/10 are the mass and inertial moment of the particles. Also, f^← f is the force due to the fluid and f^↔ p is the force due to the interaction between particles. In this study gravitational forces are not considered.

We set ρ_f=1, ν=1/400 and F_o=5 to achieve in the single-phase case a micro-scale Reynolds number of Re_λ = u' λ / ν≈ 435, (u' is the root mean square of the velocity fluctuations and λ is the Taylor length scale), which ensures an extensive inertial range of scales. 
The particle diameter D is varied between 16 ≤ D/η≤ 123 (η is the Kolmogorov length scale for the single-phase case). The number N of particles is increased as the diameter decreases, to maintain the volume fraction Φ_V =  V_p/(V_p + V_f) constant to a value of Φ_V = 0.0792, as in <cit.>. The particle density ρ_p is varied between 1.29 ≤ρ_p/ρ_f ≤ 104.7 to consider both light and heavy particles, corresponding to a variation of the mass fraction M = ρ_p V_p / (ρ_f V_f + ρ_p V_p) between 0.1 ≤ M ≤ 0.9. The smallest value of ρ_p/ρ_f=1.29 and the largest value of ρ_p/ρ_f = 104.7 have been chosen based on the results of <cit.>. As shown in figures 2 and 3 of their paper, indeed, the flow modulation only marginally changes when further increasing the density ratio ρ_p/ρ_f. 

The governing equations are numerically integrated in time using the in-house solver Fujin (<https://groups.oist.jp/cffu/code>), which is based on an incremental pressure-correction scheme. It considers the Navier–Stokes equations written in primitive variables on a staggered grid, and uses second-order finite differences 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. The solver is parallelised and is based on a Cartesian block decomposition of the domain in the x and y directions, and uses the Message Passing Interface (MPI) library for maximum portability. The governing equations for the particles are dealt with the immersed boundary method introduced by <cit.>. The fluid-solid coupling is achieved in an Eulerian framework <cit.>, and accounts for the inertia of the fictitious fluid inside the solid phase, to properly reproduce their behaviour in both the neutrally-bouyant case and in presence of density difference between the fluid and solid phases. The soft sphere collision model, first proposed by <cit.>, is used to prevent the interpenetration between particles. A fixed-radius near neighbours algorithm <cit.> is used for the particle interaction to avoid an otherwise prohibitive increase of the computational cost when the number of particles increases. 

The fluid domain is discretised using N_point=1024 points in the three directions, to ensure that all the scales down to the smallest dissipative (Kolmogorov) ones are solved, leading to η/Δ x= O(1), where Δ x denotes the grid spacing. At the initial time the particles are randomly distributed within the domain. Excluding the initial transient period needed to reach the statistically steady state, all simulations are advanced for approximately 15T, where T=L/u̅'̅ is the turnover time of the largest eddies. For the particle-laden cases with D/η=32, the simulations have been advanced for a longer period, i.e. approximately 35T, to ensure convergence of the temporal averages discussed in section <ref>. Note that the particle-laden cases with D/η=123 are the same considered in <cit.>. Details of the numerical simulations are provided in table <ref>. The adequacy of the grid resolution is assessed in appendix <ref>.



To demonstrate the independence of the results on the external forcing, additional simulations have been carried out using the forcing introduced by <cit.> to sustain turbulence; see appendix <ref>. Unlike the ABC forcing, indeed, the Eswaran & Pope's forcing has a random component, and does not generate an inhomogeneous shear at the largest scales. As shown in section <ref> and in appendix <ref>, the effect of the solid phase on the largest and energetic scales of the flow changes with the kind of forcing considered. At smaller scales, however, the results do not depend on the external forcing. At small scales the flow modulation does not depend on how energy is injected in the system, and the effect of the solid phase is substantially the same for the two considered forcings.





§ THE CARRIER FLOW





 §.§ Integral quantities



In this section we investigate the influence of the particles on the fluid phase. 
Figure <ref> shows the dependence of the fluid average kinetic energy E (x,t) on M and D. We define the fluid kinetic energy as


    E(t) = 1/2 | u(x,t) - U(x)|^2,


while · and · indicate average in time and along the homogeneous directions respectively; U(x) ≡u( x,t) is the three-dimensional fluid velocity field averaged in time. 



As discussed by <cit.>, due to the presence of the inhomogeneous mean shear induced by the external ABC forcing, the dispersed phase modulates the largest scales of the carrier flow in a way that substantially changes with the size and density of the particles. Similarly to what found by previous authors <cit.>, for large and/or light particles, i.e. D/η > 64 and/or M < 0.45, the kinetic energy E of the carrier flow monotonically decreases when the mass fraction increases and/or the particle size decreases; we call this regime A in figure <ref>. When fixing D, indeed, an increase of the density of the particles leads to an increase of the inertia of the fluid-particle system, and the same large-scale external forcing generates weaker fluctuations. When fixing M, instead, smaller particles are more effective in attenuating the fluid velocity fluctuations. A decrease of D corresponds to an increase of the total solid surface area and, therefore, to an increase of the high dissipation rate regions that form around the particles (see section <ref>). On the other hand, when particles are smaller and heavier, i.e. D/η≤ 64 and M ⪆ 0.45, the scenario changes: the kinetic energy E sharply increases, and does not have a monotonous dependence on D and M. This is what we call regime B in figure <ref>. The value of M that delimits regimes A and B increases when the particle size decreases, suggesting that the occurrence of regime B is mainly driven by the inertia of the single particles. Note that for D/η=32 and 0.45 ≤ M ≤ 0.6 and D/η=64 and M=0.45 the presence of the solid phase enhances the total fluid energy compared to the single-phase case. 

In this work, regimes A (circles in figure <ref>) and B (squares in figure <ref>) are only briefly described; we refer the interested reader to <cit.> for more details. To characterise the two regimes, we use the temporal average operator and isolate the influence of the dispersed phase on the largest and smaller scales of the flow. We focus on D/η=32, being the case that shows the strongest flow modulation, and for which E undergoes the maximum enhancement. We decompose the complete velocity field u(x,t) into its temporal mean field U(x) and the fluctuating field u'(x,t) = u(x,t) - U(x), and plot the variances of their three components in figure <ref>.



Light particles (regime A) influence the flow without introducing a preferential direction, modulating both the mean and fluctuating fields in a isotropic way: here UU≈VV≈WW and u'u'≈v'v'≈w'w'. Heavy particles (regime B), instead, modulate differently the largest and smaller scales of the flow. They attenuate the fluctuating field (i.e. the smaller scales of the flow) in a isotropic way, but modulate the mean flow field (i.e. the largest scales of the flow) towards a more energetic and anisotropic state. In this case, the presence of the dispersed phase enhances two components of the mean flow (U and V) and attenuates the third one (W). Here UU≈VV≫WW and u'u'≈v'v'≈w'w'. The increase of the total fluid energy E observed in figure <ref> for 0.45 ≤ M ≤ 0.6 is entirely due to the enhanced mean-flow contribution. In regime B the scenario is the following. Due to their large inertia, particles are not able to follow the inhomogeneous mean shear generated by the ABC forcing, and deviate towards almost straight trajectories that lay in x-y plane (see section <ref>), showing anomalous transport. In doing this, they modulate the flow towards an anisotropic and almost two-dimensional state. Here the mean-flow velocity components that lay in the plane of the trajectories of the particles (U and V) are enhanced, while the out-of-plane component (W) is attenuated. As detailed in <cit.>, this anisotropic modulation of the largest scales is due to the interaction of the particle with the inhomogeneous mean shear induced by the ABC forcing. When the inhomogeneous mean shear is not present the anisotropic modulation of the largest scales is not observed; see appendix <ref>. 

Note that when presenting the results, we deliberately choose a reference system such that for all cases the z axis is aligned with the mostly attenuated mean-flow velocity component, similarly to what done in <cit.>. In the reference system of the simulation, however, the direction aligned with the attenuated mean-flow velocity component changes with the initial condition.




 §.§ Mean flow modulation



Figure <ref> plots in the x=L/2 plane the three components of the mean-flow velocity field, i.e. (left) U, (middle) V, and (right) W, for (top) the single-phase case, and for the D/η=32 particulate cases with (middle) M=0.3 and (bottom) M=0.6.



In the single-phase case the mean flow resembles the laminar ABC profile, i.e. (U,V,W) ≈ (A sin(z) + C cos(y), B sin(x) + A cos(z), C sin(y) + B cos(x)) with A=B=C=1, similarly to what observed for the Kolmogorov flow by other authors <cit.>; note for example the cellular pattern in the map of U. In the particle-laden case with M=0.3 (regime A), the mean flow retains a similar pattern, in agreement with the isotropic and marginal mean-flow modulation observed in figure <ref>. In regime B, instead, the structure of the mean flow changes; for M=0.6 figure <ref> shows that the U and V velocity components are strongly enhanced, while the W component is attenuated. Moreover, the mean flow does not show the ABC cellular pattern anymore (see the U velocity component): the dependence on x and y is almost lost, while the sinusoidal dependence on z is retained, i.e. (U,V,W) ≈ (sin(z),cos(z),0) <cit.>.

For a quantitative description of the mean-flow modulation, figure <ref> shows the histogram of the spatial distribution of the three velocity components U(x), V(x) and W(x) for the D/η=32 particle-laden cases. Recall that U(x) ≡u(x,t) is the space-dependent fluid velocity field averaged in time.



In regime A (M < 0.45) the modulation of U(x) is isotropic and the distribution of the three velocity components almost overlap, showing an unimodal distribution centred in Û = V̂ = Ŵ = 0. In regime B (M ≥ 0.45), the W component retains a similar unimodal distribution, but the distribution becomes narrower in agreement with the progressive attenuation of its spatial fluctuations. In agreement with the spatial distribution shown in figure <ref>, instead, the U and V velocity components show a symmetric bimodal distribution. The modes ±Û = ±V̂ change with M, and can be used as an estimate of the mean-flow anisotropy; see the bottom panel of figure <ref>.



 §.§ Turbulence modulation




  §.§.§ Energy spectra









To shed light on the mechanism with which the particles modulate the flow at smaller scales, we investigate the influence of the solid phase on the scale-by-scale energy content of the carrier flow, i.e. the energy spectrum.

It is worth recalling that, although for some M and D the solid phase modulates the largest scales (κ/κ_L=1, where κ is the wave number and κ_L = 2 π/L) towards an anisotropic state (regime B), smaller scales with κ/κ_L>1 are isotropic and homogeneous for all cases (see figure <ref>). The flow modulation for κ/κ_L>1, indeed, does not depend on how energy is injected in the system, and on how particles modify the structure of the largest scales. This is clearly shown in appendix <ref>, where we present results from additional simulations carried out using the forcing introduced by <cit.> that, unlike the ABC forcing, does not generate a coherent and inhomogeneous shear at the largest scales κ/κ_L=1.

The presence of the solid phase modifies the energy spectrum of the carrier flow in a way that largely depends on the size and density of the particles. As first observed by <cit.>, particles drain energy from the fluid phase at scales larger than D, and inject it back in the fluid at smaller scales by means of their wake. This results into an energy depletion at low wavenumbers, i.e. large scales, and energy enhancement at large wavenumbers, i.e. small scales. Figures <ref> and <ref> consider separately the effect of the particle size and density on the energy spectrum. Note in passing that the single phase energy spectrum (black line) shows that the Reynolds number considered in this work leads to a proper separation of scales, with an inertial range of scales that extends to almost two decades of wavenumbers. Also, the oscillations at large wavenumbers appear as we compute the spectra using the fluid velocity field at all the mesh points of the computational domain, including those that are inside the particles <cit.>. However, when dealing with structure functions in the real space (see section <ref>), we have verified that the results do not change when the points within the particles are neglected.

We start investigating the dependence of the energy spectrum on the particle size (figure <ref>). For large particles, the energy depletion is limited at the largest scales κ < κ_p,1⪅κ_d = 2 π/D, the energy spectrum recovers the κ^-5/3 decaying predicted by the Kolmogorov theory for intermediate κ, and energy is enhanced at the smallest scales κ > κ_p,2. For example, for D/η=123 the energy depletion is observed for κ/κ_L ≤κ_p,1/κ_L ≈ 5, while energy enhancement occurs for κ/κ_L ≥κ_p,2/κ_L ≈ 230. Smaller particles amplify the overall mechanism. They drain energy from a wider range of scales, and, therefore, inject back a larger amount of energy at a wider range of small scales. For smaller particles, indeed, κ_p,1 increases, while κ_p,2 decreases. For M=0.6, for example, we measure that κ_p,1/κ_L ≈ 49, 25, 17 and 8 and κ_p,2/κ_L ≈ 189, 199, 219 and 230 for D/η=16, 32, 64 and 123; κ_p,1 well correlates with the particle diameter wavenumber, being κ_d/κ_L ≈ 96, 45, 25 and 12.5 respectively. Interestingly, for D/η≤ 32 the κ^-5/3 decay is not observed, indicating that in non-dilute suspensions of small particles the inertial energy cascade is substantially modified (see section <ref>). 

We now consider the effect of the mass fraction on the energy spectrum (figure <ref>). When fixing the size of the particles, the cutoff wavenumbers κ_p,1 and κ_p,2 almost do not vary with the mass fraction: the range of scales where particles drain or release energy does not change.

However, heavier particles interact more effectively with the fluid phase, leading to a stronger large-scale energy depletion — that explains the larger attenuation of the fluctuating energy discussed in section <ref> —, and to a larger energy enhancement at the smallest scales <cit.>. Heavier particles, indeed, result into an increase of the fluid-particle system inertia <cit.>. 

Note that, for D/η=32 and 0.45 ≤ M ≤ 0.6, the energy content at κ/κ_L = 1 is larger compared to the single-phase case, in agreement with the enhancement of the mean flow (largest scales) previously discussed in section <ref>.



  §.§.§ Scale-by-scale energy transfer









For a more detailed insight of the influence of the dispersed phase on the energy distribution mechanism, we look at the scale-by-scale energy transfer balance. It is obtained after some manipulations of the Fourier-transformed form of the Navier–Stokes equations <ref>. The energy balance can be compactly written as


    P(κ) + Π(κ) + Π_fs(κ) + D_v(κ) = ϵ,


where P(κ) is the scale-by-scale turbulent energy production due to the external forcing, Π(κ) and Π_fs(κ) are the scale-by-scale energy fluxes associated with the non linear convective term and with the fluid-solid coupling term, and D_v(κ) is the scale-by-scale viscous dissipation. They are defined as


    P(κ)    = ∫_κ^∞1/2( f̂·û^* + f̂^* ·û) dk, 
    Π(κ)    = ∫_κ^∞ -1/2( Ĝ·û^* + Ĝ^* ·û) dk, 
    Π_fs(κ)    = ∫_κ^∞1/2( f̂^↔ p·û^* + f̂^↔ p,*·û) dk, 
    
      D_v(κ)    = ∫_0^κ( -2 ν k^2 ℰ) d k.


Here ·̂ denotes the Fourier transform operator, and the superscript ·^* denotes complex conjugate. Ĝ is the Fourier transform of the non linear term ∇· (uu). For a detailed derivation of the budget equation we refer the reader to <cit.>. The production term and the fluxes are obtained integrating from κ to ∞, while the dissipation term is integrated from 0 to κ, to obtain a positive quantity that matches ϵ. Π(κ) and Π_fs(κ) do not produce nor dissipate energy at any scale, but redistribute it among scales by means of the classical energy cascade and the fluid-solid interaction. As such, when integrating from κ = 0 to κ = ∞, both Π(0) and Π_fs(0) have to be null, meaning that, in a statistical sense, the amount of energy the fluxes drain from the flow at certain scales has to match the amount of energy they release back to the flow at other scales (Π_fs has to be null for κ = ∞, as in the present case the particle-particle interaction is subdominant).
Figures <ref>, <ref> and <ref> show the dependence of the scale-by-scale energy budget on the particle size for M=0.3, M=0.6 and M=0.9.


In the single-phase case, the energy transfer mechanism is qualitatively and quantitatively well described by the Kolmogorov theory <cit.>. Energy is injected at the largest scales of the flow by the external forcing P(κ), and is transferred by means of the classical direct energy cascade — identified in equation <ref> by the non linear convection term Π(κ) — to the smallest scales, where it is dissipated by viscosity D_v(κ); see the dotted lines in the top left panels of figures <ref>, <ref> and <ref>. The presence of the solid phase introduces an additional energy flux, associated with the fluid-solid interaction Π_fs(κ), and the relevance of Π(κ) and Π_fs(κ) at the different scales changes with D and M, in agreement with the different modulation of the energy spectrum.

For light particles, M ≤ 0.3, Π_fs is small for all κ and the fluid-solid interaction only marginally influences the overall scale energy transfer (see figure <ref>). For large particles with D/η≥ 32, indeed, the fluid-solid coupling term is subdominant at all scales. For small particles with D/η=16, instead, Π_fs slightly overcomes Π for κ/κ_L ⪆ 34.

When heavier particles are considered (M ≥ 0.6), the relevance of Π_fs increases and the scale energy transfer mechanism due to the fluid-solid interaction progressively gains importance for all particle sizes (see figures <ref> and <ref>).

For large and heavy particles with D/η≥ 64 and M ≥ 0.6, the fluid-solid coupling contribution Π_fs(κ) dominates at small wavenumbers κ≤κ_p, while the non linear term Π(κ) dominates at larger wavenumbers κ>κ_p, with the cross-over κ_p increasing with the mass fraction and the particle wavenumber. For M=0.6 and M=0.9, we measure κ_p/κ_L ≈ 18 and 24 for D/η=64, and κ_p/κ_L ≈ 6 and 8 for D/η=123. Both the Π(κ) and Π_fs(κ) fluxes present a plateau, meaning that at the corresponding wavenumbers (in average) they do not drain nor release energy from the flow, but simply transfer it from larger to smaller scales. The scenario is the following: particles drain energy from the flow at scales larger than D, where -dΠ_fs/dκ<0, and released it by means of their wake at smaller scales with κ⪆κ_p ≈κ_d, where -dΠ_fs/dκ>0. The non linear energy transfer Π(κ), instead, drives the energy transfer at smaller scales κ > κ_p, where the classical energy cascade is recovered; see in figures <ref> and <ref> that here the spectrum follows the Kolmogorov κ^-5/3 decay. Π(κ) drains energy from the flow (-dΠ/dκ<0) for κ⪅κ_d, and transfers it to the smallest and dissipative scales (-dΠ/dκ>0). 

For smaller particles with D/η≤ 32 and M ≥ 0.6 the relevance of the fluid-solid coupling term increases. In fact, the plateau of Π_fs progressively increases, and the range of scales dominated by the fluid-particle interaction widens. For small and heavy particles (D/η≤ 32 and M>0.6), the fluid-solid coupling term Π_fs dominates at all scales, while the non linear term Π is largely attenuated. In this case the classical energy cascade is almost annihilated, and the solid phase generates a direct link between the largest and energetic scales of the flow and the smallest and dissipative ones. The energy injected at the largest scales is drained by the particles, and by means of their wake it is directly transferred to the smallest scales where it is dissipated. This is consistent with the modification of the energy spectrum observed in figure <ref>.

A last comment regards the influence of the size and density of the particles on the dissipative range of scales. We denote with κ_diff the wavenumber above which the dissipative term D_v(κ) dominates. Overall, our data show that κ_diff is only marginally influenced by D and M. In fact, for large particles with D/η≥ 32, κ_diff/κ_L ≈ 60 for all D and M. For small particles with D/η=16, instead, κ_diff increases with M up to κ_diff/κ_L ≈ 94: for heavier particles the dissipative range of scales shrinks, consistently with the enlargement of the plateau of Π_fs.



  §.§.§ Structure functions and intermittency





To extend the analysis done in the spectral domain, we compute the longitudinal structure functions, defined as S_p(r) = ( δ u(r) )^p, where p is the order of the structure function and δ u(r) = u(x+r)-u(x) is the velocity increment across a length scale r; see figure <ref>. 

For the single-phase case, S_p(r) ∼ r^p at small scales and S_p(r) ∼ r^p/3 in the inertial range, as predicted by the Kolmogorov theory <cit.>, with some deviations due to the flow intermittency <cit.>.

With the dispersed phase, the structure functions deviate from the single phase behaviour and the deviation progressively increases as the mass fraction increases and the particles size decreases, accordingly with the modulation of the energy spectrum. The even order structure functions flatten for r ≥ r_d ≈ D, progressively approaching a r^0 power-law, denoting that velocity fluctuations with scale larger than the particle size decorrelate and lose their coherency.

For the single-phase, the deviation of the high-order moment S_p(r) from the Kolmogorov theory is commonly associated with the intermittency of the flow <cit.>, i.e. extreme events which are localised in space and time that break the Kolmogorov similarity hypothesis. In different words, as stated by <cit.>, the flow intermittency implies that the probability distribution function of the velocity fluctuations deviates from a Gaussian distribution and that this deviation increases by decreasing the scale. These extreme events correspond to tails in the velocity increment distributions and make vast contribution to the high-order moments. The larger deviation from the Kolmogorov behaviour observed for the particle-laden cases, suggests that the dispersed phase influences the flow intermittency as well. To estimate this, we use the extended self-similarity form <cit.>, and plot one structure function against the other. In figures <ref> and <ref>, we plot S_q against S_2 for q=4 and 6, and investigate the effect of M and D on the flow intermittency. 





In the limit case where extreme events do not occur, the S_q ∼ S_2^q/2 power-law holds; the deviation from this behaviour is a measure of the flow intermittency. Before investigating the influence of the dispersed phase, it is worth noting that, as expected, S_4 and S_6 deviate from the theoretical prediction already in the single-phase case. However, at least for these values of q, the corrections proposed by <cit.> and <cit.> provide a good prediction of the data. 

Moving to the effect of the dispersed phase, figure <ref> shows that an increase of the mass fraction generally leads to a larger deviation from the S_q ∼ S_2^q/2 scaling. For D/η=123 the deviation from the single-phase is monotonous: heavier particles lead to a stronger flow intermittency. For smaller D/η, instead, the deviation is not monotonic. The deviation from the single phase increases when the mass fraction is increased up to M=0.6, but decreases for further heavier particles. This happens because the increase of the intermittency is due to the no-slip and no-penetration boundary conditions at the surface of the particles, which lead to strong velocity gradients events and widen the tails of the velocity increment distribution. The level of intermittency is, therefore, expected to increase with the relative velocity between the fluid and the solid phase |Δu|. As shown in the following, however, |Δu| does not show a monotonic dependence on M, but for the largest particles (see figure <ref>). For D/η=123, it increases monotonically with M, while for D/η=32 it increases for mass fractions up to M=0.6, and decreases for larger M, in agreement with the trends shown in figure <ref>.

Figure <ref> shows the dependence of the flow intermittency on the particle size using M=0.9 as an example. When the particle size decreases the deviation from the S_q ∼ S_2^q/2 scaling progressively increases, and the flow becomes more intermittent. Indeed, when the particle size decreases the number of particles increases (recall that Φ_V is constant), together with the total surface area of the solid phase. Therefore, the number of events associated with the particles boundary conditions increases, and the tails of the velocity increments distribution become more relevant. 





§ NEAR-PARTICLE FLOW



In this section we focus on the influence of a single particle on the surrounding flow, since the characterisation of the fluid modulation around isolated (and non isolated) particles plays a relevant role for the development of accurate particle-fluid interaction models. The classical point-particle models, indeed, are based on the simplifying assumption of separation between all turbulence scales and the particle size, the absence of hydrodynamic interactions, and the existence of a simple parametrisation of the fluid-particle interaction that does not change with the particle size and density <cit.>. These models properly perform for very dilute suspensions of small particles (Φ_V ⪅ 10^-4 and D/η⪅ 1), but are unable to describe the complex fluid-solid interaction of non-dilute suspensions of finite-size particles <cit.>. In these conditions, particles modify the surrounding flow in a way that changes with their size (D/η) and density (ρ_p/ρ_f), and with the volume fraction of the suspension <cit.>. For small and light particles, for example, the particle-fluid relative velocity is weak, and the flow does not separate from their surface. In this case, the effect of the particles on the surrounding flow is limited in space in the vicinity of the surface of the particles, and large part of the energy drained by the particles from the fluid is dissipated within the boundary layer that develops over their surface. For larger and heavier particles, instead, the fluid-particle relative velocity increases, and the flow separates from the surface of the particles with vortices being possibly shed downstream in the wake. In this case, part of the energy drained by the particles is injected back in their wake, and is then dissipated farther away from the particle. In case of non-dilute suspensions, this increased kinetic energy can also interact with downstream particles, affecting their wake and modulating their shedding. All these effects are not captured by the classical point-particle models. Also, the flow around a particle moving in a turbulent flow shows a more complex dynamics than the flow past a sphere in free stream <cit.>. Indeed, different mechanisms are expected to be at play, depending on the inertia of the particles and on the intensity of the flow fluctuations <cit.>.

We investigate the near-particle flow modulation by means of a conditional average of the flow. For each particle, we define a local Cartesian reference system (ξ,η,ζ), which is, at each time step, centred with the particle and has the ξ axis aligned with the relative fluid-particle velocity Δu. We use spherical coordinates (r,ϕ,θ), where ξ= r cos(ϕ) cos(θ), η = r cos(ϕ) sin(θ) and ζ = r sin(ϕ)  (r ∈ [0, ∞), θ∈ [0, 2π), ϕ∈ [-π/2, π/2]), and define the conditional average of the generic quantity a, as

    a_cp( r, θ, ϕ) = lim_T →∞1/T∫_t^t+T( 1/N_p∑_ip=1^N a_ip( r, ϕ,θ,τ) ) dτ.
 

Radial profiles away from the particle surface are then obtained by averaging over the surface of spherical shells centred with the particle. For the generic quantity a, we define the radial profile a_cp,r as


    a_cp,r(r) = 1/ 4 π r^2∫_0^2 π∫_-π/2^π/2a_cp(r,θ,ϕ) r^2 sin(ϕ) dϕdθ.


In this section we present the results in terms of ρ_p/ρ_f, instead of M. Note however, that the two quantities are interchangeable as Φ_V is maintained constant.



 §.§ Relative velocity and particle Reynolds number






The modulation of the flow around the particles depends on the relative motion between the fluid and the solid phase, which, indeed, is a key quantity of interest for modelling the motion of the particles. We define the velocity difference between the i_th particle and the fluid as Δu^i = u_p^i - u_f^i, where u_p^i is the velocity of the i_th particle, and u_f^i is the fluid velocity seen by the particle. Different approaches have been proposed to estimate u_f^i, see for example <cit.>. We follow the works by <cit.>, and evaluate u_f^i as the average of the fluid velocity within a shell centred with the particle at x_p^i, and delimited by x_p^i + R and x_p^i + R_sh. The external R_sh radius should be large enough to avoid the average fluid velocity to equal the particle velocity as the no-slip boundary is approached, and small enough to avoid considering portions of the fluid that are not relevant for the particle motion. Here we chose R_sh = 3 R as in <cit.>, but we have verified that the results qualitatively do not change when considering 2 R ≤ R_sh≤ 3 R.

Figure <ref> shows the dependence of |Δu|_p on ρ_p/ρ_f and D; here ·_p indicates average over particles and in time. For all densities, the fluid-solid relative velocity increases with the particle size, in agreement with the increase of the particle inertia, and the fact that particles become less able to follow the fluid. Generally, |Δu|_p increases also with the density of the particle when fixing D. Note, however, that for D/η≤ 64, |Δu|_p decreases for ρ_p/ρ_f ⪆ 17 (M ⪆ 0.6), consistently with the results shown in section <ref>.



We define the particle Reynolds number as Re_p = |Δu|_p D / ν. Figure <ref> shows the probability density functions of Re_p. Similarly to what shown by <cit.>, for small ρ_p/ρ_f the probability density function of Re_p is skewed since the components of |Δu|_p have large tails. When ρ_p/ρ_f increases, however, the probability density function becomes less skewed for all particle sizes. Heavier particles, indeed, are less able to follow the flow and thus the fluid-particle relative velocity increases, and the relevance of the tails of |Δu|_p decreases. Overall, the range of Re_p largely varies with the particle size and density, suggesting that in the considered range of parameters, the flow regime changes. Note, moreover, that for all cases considered, large values of Re_p, that are sufficient for the formation of vortical structures in the vicinity of the particles, are instantaneously attained (see the following discussion). This suggests that a quasi-steady linear drag assumption, commonly adopted in point-particle models, is not sufficient for modelling suspensions at the present conditions.



 §.§ Radial profiles


We start investigating the average radial profiles of the kinetic energy and dissipation rate. The aim is to quantify the attenuation of the energy fluctuations, and the enhancement of the dissipation rate at the particle surface for different values of D and ρ_p/ρ_f.

Figure <ref> plots the average radial flow energy distribution away from the particle surface e_k_cp,r.



The turbulent kinetic energy is strongly attenuated in the neighbourhood of the particles. The region of strong energy reduction extends out to less than 0.5 the particle radius, with a size that widens as D decreases, being approximately 0.2R, 0.25R, 0.4R and 0.5R for D/η =123, D/η=64, D/η=32 and D/η=16; see the regions where e_k_cp,r/E<1, which can be regarded as the flow region mostly influenced by the particles. This local attenuation of the kinetic energy is due to the particle-fluid density difference, and scales with the average relative fluid-particle velocity |Δu|_p. Particles with larger inertia respond only to the turbulent eddies with large time- and length-scales, and act as obstacles for fluctuations with smaller scales, attenuating therefore more effectively the fluid fluctuations. The local attenuation compared to the box average value, indeed, increases with ρ_p/ρ_f and D, up to approximately 80% for ρ_p/ρ_f = 104.7 and D/η=123. Note, however, that for D/η = 32 the near-particle energy attenuation does not have a monotonic dependence on ρ_p/ρ_f (M), as it decreases when the flow moves to the anisotropic and more energetic state (regime B) discussed in section <ref>. The energy attenuation increases with the density of the particles when 1.29 ≤ρ_p/ρ_f ≤ 4.98 (0.1 ≤ M ≤ 0.3) and 9.51 ≤ρ_p/ρ_f ≤ 104.7 (0.45 ≤ M ≤ 0.9), but it decreases when increasing ρ_p/ρ_f from ρ_p/ρ_f =4.98 (M=0.3) to ρ_p/ρ_f = 9.51 (M=0.45), i.e. when moving from regime A to regime B. It is worth stressing that, while single particles with larger D are more effective in reducing the turbulent fluctuations of the surrounding flow, when fixing the volume fraction of the suspension, smaller particles are more effective in reducing the global turbulent energy (see figure <ref>), due to the increase of the number of particles.

Figure <ref> shows the radial dissipation rate profile away from the particles surface.



Due to the large velocity gradients that develop within the boundary layer, a region of high dissipation rate arises in the neighbourhood of each particle. Part of the turbulent kinetic energy driving the relative motion between the fluid and the particle is dissipated here, giving rise to an unique cross-scale energy transfer from scales larger than D to the particle-surface boundary scales. The intensity of the near-particle dissipation rate changes with the particle size and density, similarly to what observed by <cit.>. ϵ_cp,r/ϵ increases with ρ_p/ρ_f and D up to approximately 11% for D/η=123 and ρ_p/ρ_f = 104.7. This is consistent with the increase of the the fluid-particle relative velocity with ρ_p/ρ_f and D, that leads to more intense velocity gradients at the particle surface. The width of the ϵ_cp,r/ϵ>1 region is a measure of the particle boundary layer thickness <cit.>, and its dependence on D and ρ_p/ρ_f is relevant for modelling implications. The ϵ_cp,r/ϵ>1 regions enlarges when the particle size and particle Reynolds number increase, while it only marginally changes with the density of the particles. This trend recalls what observed for the cutoff scale 2 π/κ_p,2 in section <ref>, suggesting a relation between these two quantities. Note, however, the wavenumber corresponding to the ϵ_cp,r/ϵ>1 region decreases with the particle size, while κ_p,2 increases.



 §.§ Near-particle flow


For a more detailed characterisation of the near-particle flow modulation, we now account for the direction of the relative particle-fluid velocity, and investigate the conditionally averaged flow field without averaging along θ and ϕ. When choosing the local reference system of the i_th particle, here we evaluate the fluid velocity seen by the particle u_f^i as the average fluid velocity within a shell centred with the particle and delimited by R ≤ r ≤ R_sh, where R_sh corresponds to the distance from the particle centre at which ϵ_cp,r/ϵ = 1 (see figure <ref>). We show the D/η=32 cases as an example; the results for the other particle sizes are similar and are omitted for the sake of conciseness.

Figure <ref> shows the average fluid velocity seen by the particle in the local reference system, i.e. u_cp,ℓ = u - u_f_cp; the ·_ℓ subscript refers to quantities evaluated in the (moving) local reference system of the particle. Due to the symmetries of the flow we consider only the (ξ,ζ) plane with η=0 (we show the maps for ζ≥ 0, and use the ζ≤ 0 half plane to double the statistical sample). The left and right panels are for u_cp,ℓ and w_cp,ℓ respectively (here u and w denote the velocity components aligned with the ξ and ζ axes) while from top to bottom the density of the particle is increased from ρ_p/ρ_f = 1.29 to ρ_p/ρ_f = 104.7. Figure <ref> shows instead the variance of the local velocity components, i.e. (left) u'u'_cp,ℓ and (right) w'w'_cp,ℓ, to characterise the relative flow unsteadiness; here the apex refers to fluctuations around the conditional average value. 





When approaching the particle, the averaged flow progressively decelerates, having null velocity at the (ξ,ζ)=(-R,0) stagnation point. Then, it accelerates along the upstream side of the particle, due to the favourable pressure gradient, and decelerates along the downstream side where, depending on the flow regime, it may separate giving rise to a recirculating region. The flow regime substantially changes with ρ_p/ρ_f, in agreement with the large variation of the particle Reynolds number shown in figure <ref>. For small densities, ρ_p/ρ_f=1.29, the mean flow separates and a small recirculating region arises, see e.g. the regions with u_cp,ℓ<0 and w_cp,ℓ>0 just downstream the particle. However, no vortex shedding is detected in this case. In fact, both u'u'_cp,ℓ and w'w'_cp,ℓ are relatively small, and do not show the typical localised peaks in the wake region, as shown for example in figure 13 of <cit.> and in figure 8 of <cit.>. For ρ_p/ρ_f=4.98, instead, the mean flow recirculation is not detected downstream the particle, but a localised peak of u'u'_cp,ℓ is observed along the particle side at (ξ,ζ) ≈ (0,1.5R). Although the absence of a localised peak of w'w'_cp,ℓ indicates that the classical vortex shedding does not occur, this peak of u'u'_cp,ℓ suggests that some unsteadiness has arisen. Interestingly, the maps of the fluctuations are compatible with the flow pattern shown in figure 9 of <cit.>, that discusses the flow past a sphere with an unsteadiness driven by the incoming flow fluctuations, rather than by a global instability of the wake. For ρ_p/ρ_f ≥ 17.45, instead, the maps of the mean flow and the variances are compatible with an unstable wake, and with an alternate shedding of vortices. Indeed, along the downstream side of the particle a mean recirculating region is observed, together with localised peaks of u'u'_cp,ℓ and w'w'_cp,ℓ. When the particle inertia increases, the mean recirculating region enlarges and the intensity of the fluctuations increases, revealing a stronger vortex shedding. In passing, note that, unlike in the flow past a sphere in free-stream <cit.>, here the peak of u'u'_cp,ℓ is much larger than that of w'w'_cp,ℓ: in this case the vortex shedding is substantially different, being modulated by the particle motion and the fluid velocity fluctuations.



We now move to the distribution of the fluid kinetic energy; see the left panels in figure <ref>. 
In agreement with figure <ref>, the flow energy is reduced at the particle surface, due to the no-slip and no-penetration boundary conditions, while it becomes larger than the box-average value when moving away. For small densities, ρ_p/ρ_f=1.29, the energy distribution is almost symmetric, in agreement with the above discussed lack of vortex shedding in the wake. When ρ_p/ρ_f increases, instead, the low energy region with e_k_cp/E<1 becomes asymmetric and extends mainly downstream the particle, in agreement with the visualisations of previous authors <cit.>. This is consistent with the distribution of the mean flow and of the velocity fluctuations shown in figure <ref>. For large ρ_p/ρ_f a region with e_k_cp/E>1 arises over the lateral sides of the particle, (ξ,ζ) ≈ (0, 1.5D), being the trace of the larger mean flow acceleration and larger fluctuating energy discussed above. 

The right panels of figure <ref> consider the fluid dissipation rate. In agreement with figure <ref>, a region of high dissipation rate is observed in the immediate neighbourhood of the entire particle surface, being associated with the large velocity gradients that develop in the particle boundary layer. The intensity of ϵ_cp/ϵ in this region increases with ρ_p/ρ_f, consistently with the increase of the particle Reynolds number and with the progressive decorrelation of the particle velocity and the local flow. For small ρ_p/ρ_f, ϵ_cp/ϵ≈ 1 outside this region, meaning that the influence of light particles on the flow dissipation rate is relatively limited in space. For larger ρ_p/ρ_f, instead, ϵ_cp/ϵ < 1 outside the particle boundary layer, indicating a wider influence of the particle. Aside from this region, two further regions of relatively high dissipation rate are observed. The first is observed for all ρ_p/ρ_f and is placed upstream the particle; here ϵ_cp/ϵ is maximum. In fact, the largest velocity gradients are attained close to the front stagnation point, where the flow deceleration is strong. Note that this region of large dissipation has been also observed by <cit.>, when considering particles with 1 ≤ρ_p/ρ_f ≤ 4 in homogeneous isotropic turbulence at a lower Reynolds number Re_λ≈ 70. The second region is observed only for large ρ_p/ρ_f, and is associated with the separated flow in the wake region. Note that it is barely visible for ρ_p/ρ_f=17.45, while clearly visible for ρ_p/ρ_f=104.7. By means of the vortex shedding, indeed, heavy particles modulate the flow dissipation rate over a wider portion of the wake region.





§ PARTICLE DYNAMICS





 §.§ Particles velocity and trajectories







In this section we briefly characterise the dynamics of the particles in the A and B flow regimes mentioned in section <ref>; we refer the reader to <cit.> for a detailed discussion. Figure <ref> shows typical trajectories for the D/η=32 particle-laden cases with M=0.3 (left) and M=0.6 (right), which are representative of the motion of the particles in the two regimes. Figure <ref>, instead, quantitatively characterises the velocity of the particles for D/η=32 and 0.1 ≤ M ≤ 0.9.

In regime A (D/η > 64 and M ≤ 0.3) particles follow, at least partially, the cellular motion imposed by the external ABC forcing, and exhibit complex trajectories in agreement with the chaotic Lagrangian structure of the ABC flow <cit.>; see the left panel of figure <ref>. In this case, particles do not select a preferential direction, and progressively span the complete computational domain. Accordingly, the probability density functions of the three velocity components almost collapse, and show a symmetric unimodal distribution with the modes being û_p = v̂_p = ŵ_p = 0. When M and/or D are increased, the three distributions narrow (i.e. the variances decrease) due to the larger inertia, and the particle velocity experiences smaller fluctuations.

In regime B (D/η≤ 64 and M ≥ 0.45) particles show a completely different behaviour. Due to the large inertia, they deviate from the cellular flow path imposed by the ABC forcing, and move along almost straight trajectories showing anomalous transport <cit.>, see the right panel of figure <ref>. Recall that, in regime B, we denote with z the direction orthogonal to the plane where the trajectories of the particles lay, which corresponds to the direction aligned with the attenuated mean flow velocity component (see section <ref>). Compared to regime A, here the fluctuations of the in-plane u_p and v_p particle velocity components are strongly enhanced, while the fluctuations of the out-of-plane w_p component are attenuated, see figure <ref>. The direction of the trajectories in the x-y plane changes with z, as the mean flow and particles velocity retain the sinusoidal dependence on z inherited by the ABC forcing (see also section <ref>). Accordingly, the probability density function of w_p exhibits a very narrow symmetric unimodal distribution centred in ŵ_p=0 (see figure <ref>), while the distributions of u_p and v_p show a symmetric bimodal distribution with the modes ±û_p = ±v̂_p that change with the particle size and density.



 §.§ Clustering



After the characterisation of the dynamics of a single particle, in this section we focus on their collective motion and investigate how it changes with the particle size and density. We look at the local concentration of the suspensions to inspect the presence of clusters. Different metrics have been proposed over the years for detecting clusters <cit.>. Voronoï tesselletion is a computationally efficient tool for the analysis of the spatial arrangement of the dispersed phase in the carrier flow, and has been used by several authors <cit.>. The position of each particle is determined by its centre, and the computational box is partitioned such that each grid point is associated with the closest particle. The Voronoï cell of each particle, therefore, consists of the ensemble of grid cells that are closer to it. The inverse of the volume of a Voronoï cell provides a measure of the local concentration: a smaller Voronoï cell corresponds to a denser particle arrangement, while a larger Voronoï cell to a more sparse distribution. <cit.> have shown that when point particles are randomly drawn in a box, the probability density function of the corresponding Voronoï cell volumes exhibits a gamma distribution with parameters that can be evaluated analytically. However, for finite-size particles this is not the case, as they can not overlap. This means that the parameters describing a random distribution of finite-size particles need to be computed separately for different particle sizes, volume fractions and box sizes. For particle suspensions that exhibit clustering, the variance of the Voronoï cell volumes is larger than that computed for the corresponding random distribution <cit.>.



Figure <ref> shows the time average of the variance of the volume of the Voronoï cells for different particles sizes and mass fractions. It is clearly visible that for all mass fractions smaller particles exhibit stronger clustering. Due to their lower inertia, indeed, smaller particles are more able to follow the fluid motion, and this may promote their clustering due to the particle preferential sampling (see section <ref>). In contrast, the influence of the mass fraction on the clustering is not monotonous, and changes with the particle size. For large particles, D/η=123, heavier particles exhibit weaker clustering; for M=0.9 we compute σ/σ_rand≈ 1 indicating the almost absence of clusters. For D/η≤ 64, instead, the level of clustering is maximum for intermediate mass fractions; σ/σ_rand is maximum at M=0.3 for D/η=32 and D/η=64, and at M=0.6 for D/η=16. As shown in the following (see figure <ref> and relative discussion), this non monotonous dependence of the level of clustering on M is related to the different trajectories of the particles in the A and B regimes discussed in section <ref>.

The low level of clustering observed for M=0.1 recalls the results of <cit.> that report no clustering for neutrally buoyant particles with size D/η≈ 20. The present results, however, give evidence that clusters are present in non-dilute suspension of relatively small and heavy finite-size particles. 

<cit.> used the Voronoï tessellation to provide an objective definition of a cluster.



As shown in figure <ref> for D/η=32, the fact that the variance of the actual simulation is larger than the corresponding random arrangement of particles leads to two cross-overs between the two associated probability density functions of the Voronoï cells volume. These intersection points can be used to identify clusters and void regions. <cit.> proposed that all particles associated with a Voronoï cell with volume smaller than the lower cross-over point are part of a cluster, while all particles associated with a Voronoï cell with volume larger than the larger cross-over point are part of void regions. Moreover, particles with Voronoï cells smaller than the lower cross-over point that share at least one vertex, are part of the same cluster.



Following this approach, figure <ref> provides a qualitative view of the clustering for D/η=32; note that qualitatively the same results have been found also for the other particle sizes. The figure shows the clusters found in an instantaneous snapshot in the x=L/2 plane for 0.1 ≤ M ≤ 0.9. Different colours indicate Voronoï cells associated with different particles. Accordingly with the previous discussion, the level of clustering is strong for M=0.3, while rather weak for M=0.1 and M=0.9. The spatial arrangement of the clusters changes with the trajectories of the particles (see section <ref>), explaining the non monotonous dependence of σ/σ_rand on M for D/η≤ 64. For M=0.3, indeed, the cluster arrangement recalls the cellular shape inherited by the ABC forcing, and does not reveal a preferential direction. For M=0.6, instead, clusters are stretched and aligned in the y direction, accordingly with the almost two-dimensional motion of heavy particles previously discussed. The low level of clustering for M=0.9 shows that heavy particles move along almost straight lines in a isolated manner.

An alternative way for characterising clustering is the radial distribution function g(r) <cit.> defined as


    g(r) = N_s(r)/Δ V_i(r)/N_p/V.


Here, N_s(r) is the number of particle pairs separated by a distance within r-Δ r/2 and r + Δ r/2, Δ V_i is the volume of the discrete shell located at r, N_p=N(N-1)/2 is the total number of particle pairs in the sample, and V is the total volume of the system. This quantity describes the probability of having a pair of particles at a given mutual distance, and its magnitude characterises the strength of clustering at scale r. In a perfectly uniform distribution of point particles (where the overlap between particles is possible) g(r)=1 for all r. 



Figure <ref> shows the dependence of the radial distribution function g(r) on the mass fraction. For large particles, D/η≥ 64, it turns out that the accumulation at small length scales r ≈ D is maximum for light particles, while progressively decreases with M. At large length scales, instead, the accumulation is maximum for M=0.3, with the minimum found for M=0.9. The overall scenario agrees with the analysis of the Voronoï tessellation, where the maximum level of clustering was found for M=0.3. Note that besides the dominant peak at r = D, a second less evident peak is observed for r ≈ 2D. They are respectively the statistical trace of the first and second coordination shells around the test particle. For smaller particles D/η=32 and D/η=16, the accumulation is largest at all length scales for M=0.3 and M=0.6, respectively. This agrees with the large increase of the variance of the Voronoï volumes seen in figure <ref>. Overall, the radial distribution function confirms that for D/η≤ 64 the level of clustering increases with M for small mass fraction, but decreases for larger M, when due to their larger inertia particles move along almost two-dimensional and straight trajectories.



Figure <ref> details the effect of the particles size on the radial distribution function. It is clearly visible that for all mass fractions the level of accumulation increases at all length scales when the particle size decreases. The same trend has been found by <cit.> when considering finite-size particles with D/η=5 and D/η=11 and ρ_p/ρ_f = 1.5, while the opposite trend, i.e. an increase of the accumulation with the particle size, has been reported by <cit.> for particles smaller than the Kolmogorov length scale D < η. 



We have also observed that for small mass fractions, M = 0.1, the decay of the radial distribution function with the distance r from the test particle approximately follows a power law, which indicates a self-similar spatial distribution <cit.>. Despite theoretical arguments support this formulation for dissipative separations r/η<1 and for small St only, several authors provided numerical and experimental evidence that the power-law form continues to hold also for larger r <cit.>. For point particles, <cit.> observed that the clustering mechanisms operating in the inertial range are analogous to those operating in the dissipation range <cit.>. They argued that when St_r ≪ 1 (where St_r = τ_p/τ_r is the Stokes number based on
τ_r, i.e. the eddy turnover time at scale r defined as τ_r = ϵ^-1/3r^2/3), preferential sampling of the coarse-grained fluid velocity gradient tensor at scale r generates the inward drift and clustering. When, instead, St_r > O(1) the non local, path-history mechanism contributes to the clustering, breaking thus the self-similarity. In figure <ref>, we consider M = 0.1 and 16 ≤ D/η≤ 123, for which the Stokes numbers are between 0.1 ⪅ St ⪅ 10. Although the results for D/η = 16 suggest that self-similarity is broken, for larger D/η we measure a power law of approximately g(r) ∼  (r/η)^-1 up to r/η≈ O(10^2).




It is instructive to investigate the distribution of the particle-particle relative velocity δu_p. In the context of point-particles, indeed, several theories regarding the clustering mechanisms rely on the exact equation for the radial distribution function g(r), in which the distribution of δu_p plays an important role; see for example <cit.>. Figure <ref> shows the distribution of δu_p ·δx_p/|δx_p|, i.e. the component of the relative velocity between two particles δu_p projected along their separation vector δx_p. When δu_p ·δx_p/|δx_p|>0 the two particles depart, while when δu_p ·δx_p/|δx_p|<0 they get closer. For conciseness, figure <ref> considers only the D/η=32 case, but the following discussion holds also for the other particle sizes. Note that for all |δx_p|/D, the average value of the distribution is zero, due to the statistically steady state condition considered in this work. However, the distribution of δu_p ·δx_p/|δx_p| is not symmetric, in agreement with the above discussed tendency of particles to form clusters. For all |δx_p|/D, indeed, the distributions are left skewed: the mode is slightly positive, but the negative tails are longer. When increasing |δx_p|/D the distribution becomes progressively more flat, in agreement with a less correlated motion of the particles. Figure <ref> shows that the dependence of the distribution on M agrees with what found for g(r). For M=0.9, indeed, the positive and negative tails of the distribution almost overlap, consistently with the low value of clustering detected in figure <ref>. A similar effect is found when varying the particle size (not shown): when fixing M, an increase of D leads to a more symmetric distribution, in agreement with the above discussed decrease of the level of clustering.






 §.§ Preferential particle location



In the previous section we have shown that the solid phase is not randomly distributed in space and that, depending on the size and density of the particles, the suspension features a mild level of clustering. In this section we investigate the probability density functions of the acceleration and vorticity vectors in the region surrounding the particles, to speculate whether the two principal mechanisms that are known to govern the particle preferential location for sub-Kolmogorov particles may play a role also in the case of particles with size that lays in the inertial range. We consider the centrifuge <cit.> and the sweep-stick <cit.> mechanisms. The centrifuge mechanism is based on the hypothesis that (i) the particle inertia is sufficiently large, (ii) the particle paths deviate from the fluid path, and (iii) the excess inertia does not place the particles in the ballistic regime. In this case, particles are centrifuged out from regions of high vorticity (vortex cores), preferring regions of large strain.
The sweep-stick mechanism, instead, links the particle locations with the fluid acceleration a_f. In the framework of the one-way coupling point-particle model, following the work by <cit.>, <cit.> showed that the particles accumulate in points that satisfy the following criterion e_1 ·a_f = 0 and λ_1=0, where λ_1 is the largest eigenvalue of the symmetric part of the acceleration gradient tensor, and e_1 is the associated eigenvector. This results from the observation that the slip velocity between a point-particle and the fluid is proportional to the fluid acceleration, i.e. u_p - u_f ≈ - τ_p a_f. Later, <cit.> found that point particles in homogeneous and isotropic turbulence accumulate preferentially in the vicinity of zero-acceleration points, having observed that this simpler criterion is more strongly correlated with the location of their point-particles. 

It is worth mentioning that, by analysing the governing equation for the radial distribution function, <cit.> suggested that the clustering mechanisms in the inertial and dissipative ranges are analogous. For any r which is less than the integral length scale of the flow, they propose that the clustering mechanism for small St is associated with centrifuging out effect of eddies of that scale. For St > O(1), instead, they found that non local mechanisms contribute to the clustering, generating a statistical asymmetry of the path history of approaching and separating particle pairs. However, despite the universality of the clustering mechanism across the range of scales disagrees with the presence of different mechanisms, they observe that in the inertial range for St ≪ 1 the mechanism they propose is essentially equivalent to the sweep-stick mechanism.

We investigate whether the centrifuge and sweep-stick mechanisms play a role in determining the preferential location of particles with size in the inertial range in non-dilute suspensions. This analysis follows the work by <cit.>, that investigated the preferential sampling of light and small particles (D/η = 5-11 and ρ_p/ρ_f=1.5) in more dilute suspensions (Φ_V=0.005), and at smaller Reynolds numbers (Re_λ≈ 100). The idea is to compute the probability density functions of the modulus of the acceleration and vorticity vectors in shells around the particles, and to compare them with the probability density functions of the same quantities evaluated for the complete fluid phase. The difference between the probability density functions may reveal whether there is a link with these mechanisms or not. In case the sweep-stick mechanism plays a role in the particle locations, we expect an increase of the probability of small |a_f| events in the region surrounding the particles. Similarly, an increase of the probability of small |ω_f| events in the neighbourhood of the particles may suggest a link with the centrifuge mechanism. Before going on, it is worth stressing that the interpretation of these results in our case is not straightforward. In fact, unlike what considered by the point-particle approach for sub-Kolmogorov particles in dilute suspensions, finite-size particle substantially influence the surrounding flow (see section <ref>), and can potentially collide. Therefore, the differences between the global and local probability density functions are due to a combination of three effects: (i) the preferential particle location, (ii) the influence of the particle on the surrounding flow and (iii) the influence of the particle-particle interaction on the surrounding fluid phase. A further consideration regards the radius of the shell R_sh around the particles that is used for the local probability density functions. If R_sh is too small, only the local perturbations of the particles on the surrounded flow are considered. If R_sh is too large, spurious contributions of the fluid phase that do not influence the particle location would be taken into account. Following section <ref>, we choose R_sh = 1.5 R_p, but we have verified that small variations of R_sh qualitatively do not influence the results. Clearly, for larger R_sh the differences between the local and global probability density functions progressively decrease. 

Figure <ref> plots the local and global probability density functions of the modulus of the fluid acceleration |a_f| for D/η=64 and D/η=123 as representative cases.



We start considering small mass fractions M=0.1 (ρ_p/ρ_f = 1.29). In this case, for all D considered, the probability of large acceleration in the shell decreases compared to the rest of the fluid phase, while the probability of small values increases suggesting a link with the sweep-stick mechanism. Note, however, that for this M the difference between the local and global probability density functions is relatively small, in agreement with the low level of clustering above discussed. This result is in line with the results of <cit.>, that for particles with D/η=5-11 and ρ_p/ρ_f = 1.5 found a non negligible correlation between the particle position and the sticky points. Therefore, despite the theory developed by <cit.> rigorously holds only in the limit of the point-particle approach, our results indicate that a link between the preferential location of light particles with size within the inertial range and low fluid acceleration areas may exist.

When larger M are considered the scenario changes and the correlation between the particle location and regions with low fluid acceleration is less clear. The peak of the global probability density functions moves towards smaller |a_f|, in agreement with a decrease of the mean fluid acceleration. In the shell surrounding the particles, instead, the probability density function becomes progressively more flat as M increases. Comparing the two probability density functions, the probability of both low and large |a_f| values increases in the neighbourhood of the particles, while the probability of intermediate values decreases. The increased probability of large values of |a_f| with M in the region surrounding the particles is consistent with the increase of the particle inertia and of the relative velocity between the fluid phase and the particle, that results into an increase of the momentum exchange between the fluid and solid phases.

To investigate whether the sweep-stick mechanism actually plays a role for light and/or heavy particles, we compute the probability g_ps(r) for each particle to have at a certain distance r a point where the fluid acceleration is smaller than a threshold a_th. When g_ps(r)> 1 the probability is larger compared to a random distribution, whereas when g_ps(r)<1 it is lower. In figure <ref>, we show the dependence of g_ps(r) on M for D/η=64. Following the work by <cit.> we have set a_th=0.05 |a|_rms, and we have verified that the results do not change significantly by slightly decreasing or increasing this value. 



For all mass fractions there is a large probability of having fluid points with low acceleration immediately close to the particles. This is a direct effect of inertial particles on the surrounding flow, and is associated with the no-slip and no-penetration condition on the surface. Due to their inertia, particles face a lower acceleration with respect to the fluid phase and, therefore, the same holds for the boundary layer surrounding them. Note that the increase of g_ps with the mass fraction for small r agrees with figure <ref>, that shows that the probability density function of |a|_f in the shell surrounding the particles progressively flattens when M increases. Moving away from the particle boundary layer, i.e. at larger r, the distribution of g_ps(r) changes with the mass fraction. For M=0.1 (ρ_p/ρ_f=1.29), g_ps remains larger than 1, with a value of approximately g ≈ 1.2, up to r ≈ 3R, while for M>0.3 (ρ_p/ρ_f > 4.98) g_ps<1 for all r. This clarifies the above discussion: for small mass fractions a link between the particle location and regions of low fluid acceleration is possible also for particles with size within the inertial range. This can be explained by the fact that a particle of size D behaves to leading order like a point particle with respect to eddies of scale ℓ≫ D. Hence, from a qualitative point of view, particles of size D can exhibit sweep-stick like phenomena when the scale of the flow that drives their motion is of size ℓ≥ O(D). For larger M, instead, this is not the case. Indeed, the increase of the probability of small |a| values seen in figure <ref> is only due to the increase of the probability of negative values in the boundary layer attached to the particle surface, while g_ps is less than one for 0.14 ⪅ r/R_p-1 ⪅ 1.4, indicating a lower probability of having points with |a|<a_th compared to a random distribution.

We now move to the centrifuge mechanism.



For all mass fractions and particle sizes, the probability of small vorticity values in the region surrounding the particles is smaller than in the complete box, while the probability of large values increases (see figure <ref>). For larger mass fraction and larger particles we observe that the difference between the two probability density functions progressively increases, in agreement with the increase of the particle Reynolds number and with the occurrence of vortex shedding in the particle wake (see section <ref>). 
Based on this observation, one may be tempted to exclude that the centrifuge mechanism plays a role in determining the preferential location of finite-size particles. However, it is possible that a particle of size D is centrifuged out by a vortex of size ℓ≫ D. While doing this, a strong vorticity may generate on the surface of the particle due to the no-slip and no-penetration boundary conditions, potentially leading to the above effect on the local probability density function.








§ CONCLUSION



In this work we have investigated the fluid-solid interaction of suspensions of solid spherical particles in triperiodic turbulence, at the relatively large micro-scale Reynolds number of Re_λ≈ 400. The study is based on direct numerical simulations coupled with an immersed boundary method, to properly resolve the flow around each particle. The fluid-solid interaction is studied in the two-dimensional parameter space of the particle size and density. The particle size is varied between 16 ≤ D/η≤ 123, while the solid-to-fluid density ratio is varied between 1.29 ≤ρ_p/ρ_f ≤ 104.7, corresponding to a variation of the mass fraction between 0.1 ≤ M ≤ 0.9. Non-dilute suspensions with a volume fraction of Φ_V = 0.079 are considered. Turbulence is sustained with the Arnold-Beltrami-Childress (ABC) cellular forcing <cit.>, that generates a three-dimensional and inhomogeneous shear at the largest scales of the flow.

In the considered portion of the parameter space, the solid phase modulates the largest scales of the flow in a way that changes with the size and the density of the particles. Depending on the ratio between the particle size and the length-scale of the inhomogeneous mean shear induced by the forcing, the solid phase may modulate the largest scales of the flow towards an anisotropic and almost two-dimensional state <cit.>. The smaller scales, instead, are homogeneous and isotropic for all considered cases, and their modulation does not depend on the external forcing. For these scales, the independence of the results on the external forcing has been assessed using the forcing introduced by <cit.>. The influence of the solid phase on the energy spectrum indicates that the mechanism driving the energy transfer across scales changes with the particle size and density. By means of the scale-by-scale energy budget, we have shown that the two mechanisms that drive the energy transfer, i.e. the classical energy cascade described by the Kolmogorov theory and the energy transfer associated with the fluid-solid coupling, play a different role depending on the size and density of the particles. For large and light particles (D/η≥ 64 and M ≤ 0.3), the flow modulation is weak and the classical energy cascade dominates, being only marginally altered by the solid phase. For large and heavier particles (Dη≥ 64 and M ≥ 0.6), the two mechanisms coexist: the fluid-solid coupling term drives the energy transfer at scales larger than the particle size, while the classical energy cascade is recovered at smaller scales. For small and heavy particles (D/η≤ 32 and M>0.6), instead, the classical energy cascade is subdominant, and the overall energy transfer across scales is driven by the fluid-solid coupling term, which links directly the largest scales where energy is injected, and the smallest scales where energy is dissipated. By means of the extended similarity <cit.>, we have shown that the solid phase increases the flow intermittency in a way that depends on the relative velocity between the particles and the surrounding flow.

By means of a conditional average of the flow in a shell surrounding the particles, we have shown that the near-particle flow pattern largely changes with D and ρ_p/ρ_f, with interesting implications for modelling. Light particles (ρ_p/ρ_f < 5) partially follow the flow, and their influence on the near-flow region is limited within the boundary layer that develops over their surface. In this case, vortex shedding does not occur, and the energy drained by the particle-fluid relative velocity is mainly dissipated at the particle surface. When the particle density increases (ρ_p/ρ_f ≥ 5), instead, the influence of the particles on the surrounding flow becomes more intense, and the particle wake becomes unsteady. Here, part of the energy drained by the particle is dissipated away from the particle, and this requires a more sophisticated approach for modelling. For intermediate particle sizes, our results suggest that an actual vortex shedding occurs only for heavy particles with ρ_p/ρ_f ≥ 17, and that for lighter particles with ρ_p/ρ_f ≈ 5 the wake unsteadiness is essentially due to the relative motion between the particles and the fluid phase <cit.>.

We have also quantified the particle clustering with the aid of the Voronoï tessellation <cit.> and the radial distribution function <cit.>. Both methods show that, for all mass fractions, the level of clustering increases when the particle size decreases. The dependence of the clustering on the particle density, instead, is not monotonic and changes with D. For large particles (D/η=123) the level of clustering decreases with the particle density, being almost negligible for the largest mass fractions. For smaller particles (D/η≤ 64), instead, the level of clustering is not monotonic, being maximum for intermediate ρ_p/ρ_f. Following the work by <cit.>, we have then explored the possibility that at the considered parameters the preferential location of the particles might be determined by similar effects as for point-particles. We have investigated whether the centrifuge <cit.> and the sweep-stick <cit.> mechanisms play a role in determining the preferential location of light and heavy particles, with size that lays in the inertial range. Overall, our results indicate that a link between the preferential location of the particles and areas of very small fluid acceleration is possible, but only for light particles with ρ_p/ρ_f = 1.3; however, particular care is needed when looking at these results because, unlike in the limit of point particles, here the particles are altering the flow in their surroundings, thus different effects are mixed in the analysis. 



§ ACKNOWLEDGMENTS
 
The authors acknowledge the computer time provided by the Scientific Computing & Data Analysis section of the Core Facilities at OIST and the computational resources of the supercomputer Fugaku provided by RIKEN through the HPCI System Research Project (Project IDs: hp210229 and hp210269).
  


§ FUNDING
 
The research was supported by the Okinawa Institute of Science and Technology Graduate University (OIST) with subsidy funding from the Cabinet Office, Government of Japan.
  


§ DECLARATION OF INTERESTS
 
The authors report no conflict of interest.



§ AUTHOR ORCIDS


Alessandro Chiarini, https://orcid.org/0000-0001-7746-2850; 

Marco Edoardo Rosti,  https://orcid.org/0000-0002-9004-2292.



§ DEPENDENCE OF THE RESULTS ON THE EXTERNAL FORCING



In this section the dependence of the results on the external forcing is assessed. We show that the flow modulation for κ/κ_L>1 does not depend on the inhomogeneous mean shear induced by the ABC forcing at κ/κ_L=1. In this respect, we have carried out additional simulations using the forcing introduced by <cit.>, that does not generate an inhomoegeneous mean shear at the largest scales, and looked at the flow modulation induced by the addition of particles. The particle size has been fixed at a value of D/L = 0.0207, the volume fraction has been kept constant at Φ_V = 0.079, and the mass fraction has been varied between 0.1 ≤ M ≤ 0.9. For these additional simulations the Reynolds number has been set to a smaller value of Re_λ≈ 260 for the single phase case, meaning that D/η≈ 18. The grid resolution is kept constant with N_point=1024 points along each direction.



Due to the absence of the inhomogeneous mean shear, when using the <cit.>'s forcing, particles modulate all scales in an isotropic way for all M, and the energy enhancement observed in figures <ref> and <ref> for regime B does not occur. This is shown in the inset of figure <ref> where E/E_0 is plotted as a function of M.

We now look at the scale-by-scale energy spectrum. Figure <ref> plots ℰ for these additional simulations, to be compared with figure <ref>. For κ/κ_L>1 the way the particles modulate the carrier flow does not change with the external forcing. Indeed, figure <ref> confirms the conclusions drawn with the ABC forcing. Particles drain energy at scales larger than D (see the energy depletion for κ < κ_p,1) and release it at smaller scales by means of their wake (see the energy enhancement for κ > κ_p,2), with heavier particles being more effective. When comparing the energy spectra obtained with the different forcings, the only qualitative difference regards the largest scale. Indeed, when using the ABC forcing, for 0.45 ≤ M ≤ 0.75 ℰ(κ/κ_L=1) becomes larger than the single-phase value, due to the above described mean flow enhancement.



In addition, in figure <ref> we report the scale-by-scale energy transfer balance computed with these additional simulations, to be compared with figures <ref>, <ref> and <ref>. Again, the results are in agreement with what found using the ABC forcing. For small mass fractions the non linear term Π dominates at all scales, and the fluid-solid coupling term only marginally influences the scale energy transfer. For large M, instead, Π_fs dominates at all scales and the classical energy cascade is substantially annihilated.




§ SENSITIVITY TO THE GRID RESOLUTION





In this section the adequacy of the grid resolution is investigated. An additional simulation has been carried out for D/η=16, i.e. the smallest particle considered, on a coarser grid. The mass fraction has been set to M=0.3, and turbulence is sustained by means of the ABC forcing. In this coarser grid the fluid domain is discretised using N_point=512 points in the three directions, leading to η/Δ x ≈ 2. In doing this, the number of points across each particle decreases from 16 (in the standard grid) to 8. Figure <ref> shows that the energy spectra obtained with the standard and coarse grids collapse pretty well, indicating that the resolution of the standard grid is adequate.



jfm
94
natexlab#1#1#1#1 #1#1 #1#1#1#1#1#1
  #1#1#1#1 #1#1 #1#1
  #1#1#1#1#1#1#1#1

[Bagchi & Balachandar(2003)]bagchi-balachandar-2003
Bagchi, P. & Balachandar, S. 2003  Effect of
  turbulence on the drag and lift of a particle.  Phys. Fluids
  15 (11),  3496–3513.

[Balachandar & Eaton(2010)]balachandar-eaton-2010
Balachandar, S. & Eaton, John K. 2010  Turbulent
  Dispersed Multiphase Flow.  Annu. Rev. Fluid Mech.
  42 (1),  111–133.

[Batchelor(1970)]batchelor-1970
Batchelor, G. K. 1970  The stress system in a suspension of
  force-free particles.  J. Fluid Mech.  41 (3),  545–570.

[Batchelor & Green(1972)]batchelor-green-1972
Batchelor, G. K. & Green, J. T. 1972  The
  determination of the bulk stress in a suspension of spherical particles to
  order c2.  J. Fluid Mech.  56 (3),  401–427.

[Benzi et al.(1993)Benzi, Ciliberto, Tripiccione, Baudet,
  Massaioli & Succi]benzi-1993
Benzi, R., Ciliberto, S., Tripiccione, R., Baudet, C.,
  Massaioli, F. & Succi, S. 1993  Extended self-similarity
  in turbulent flows.  Phys. Rev. E  48 (1),  R29–R32.

[Biferale(2003)]biferale-2003
Biferale, L. 2003  Shell Models of Energy Cascade
  in Turbulence.  Annu. Rev. Fluid Mech.  35 (1),
  441–468.

[Borue & Orszag(1996)]borue-orszag-1996
Borue, V. & Orszag, S. A. 1996  Numerical study of
  three-dimensional Kolmogorov flow at high Reynolds numbers.  J.
  Fluid Mech.  306,  293–323.

[Botto & Prosperetti(2012)]botto-prosperetti-2012
Botto, L. & Prosperetti, A. 2012  A fully resolved
  numerical simulation of turbulent flow past one or several spherical
  particles.  Phys. Fluids  24 (1),  013303.

[Bragg & Collins(2014)]bragg-collins-2014
Bragg, A. D. & Collins, L. R. 2014  New insights from
  comparing statistical theories for inertial particles in turbulence: I.
  Spatial distribution of particles.  New Journal of Physics
  16 (5),  055013, publisher: IOP Publishing.

[Bragg et al.(2015)Bragg, Ireland &
  Collins]bragg-ireland-collins-2015
Bragg, A. D., Ireland, P. J. & Collins, L. R. 2015
  Mechanisms for the clustering of inertial particles in the inertial range
  of isotropic turbulence.  Physical Review E  92 (2),
  023029, publisher: American Physical Society.

[Brandt & Coletti(2022)]brandt-coletti-2022
Brandt, L. & Coletti, F. 2022  Particle-Laden
  Turbulence: Progress and Perspectives.  Annu. Rev. Fluid
  Mech.  54 (1),  159–189.

[Breard et al.(2016)Breard, Lube, Jones, Dufek, Cronin,
  Valentine & Moebis]breard-etal-2016
Breard, E. C. P., Lube, G., Jones, J. R., Dufek, J.,
  Cronin, S. J., Valentine, G. A. & Moebis, A. 2016
  Coupling of turbulent and non-turbulent flow regimes within pyroclastic
  density currents.  Nat. Geosci.  9 (10),  767–771.

[Breugem(2012)]breugem_2012a
Breugem, W. P. 2012  A second-order accurate immersed
  boundary method for fully resolved simulations of particle-laden flows.
  J. Comput. Phys.  231 (13),  4469–4498.

[Burton & Eaton(2005)]burton-etal-2005
Burton, T. M. & Eaton, J. K. 2005  Fully resolved
  simulations of particle-turbulence interaction.  J. Fluid Mech.
  545,  67–111.

[Chen et al.(2006)Chen, Goto &
  Vassilicos]chen-goto-vassilicos-2006
Chen, L., Goto, S. & Vassilicos, J. C. 2006
  Turbulent clustering of stagnation points and inertial particles.
  J. Fluid Mech.  553,  143–154.

[Chiarini et al.(2024)Chiarini, Cannon &
  Rosti]chiarini-etal-2023
Chiarini, A., Cannon, I. & Rosti, M. E. 2024
  Anisotropic mean flow enhancement and anomalous transport of finite-size
  spherical particles in turbulent flows.  Phys. Rev. Lett.  132,
  054005.

[Cisse et al.(2015)Cisse, Saw, Gibert, Bodenschatz &
  Bec]cisse-etal-2015
Cisse, M., Saw, E.-W., Gibert, M., Bodenschatz, E. &
  Bec, J. 2015  Turbulence attenuation by large neutrally
  buoyant particles.  Phys. Fluids  27 (6),  061702.

[Coleman & Vassilicos(2009)]coleman-vassilicos-2009
Coleman, S. W. & Vassilicos, J. C. 2009  A unified
  sweep-stick mechanism to explain particle clustering in two- and
  three-dimensional homogeneous, isotropic turbulence.  Phys. Fluids
  21 (11),  113301.

[Constantinescu & Squires(2004)]constantinescu-squires-2004
Constantinescu, G. & Squires, K. 2004  Numerical
  investigations of flow over a sphere in the subcritical and supercritical
  regimes.  Phys. Fluids  16 (5),  1449–1466.

[Costa et al.(2020)Costa, Brandt &
  Picano]costa-brandt-picano-2020
Costa, Pedro, Brandt, Luca & Picano, Francesco 2020
  Interface-resolved simulations of small inertial particles in turbulent
  channel flow.  Journal of Fluid Mechanics  883,  A54.

[Costa et al.(2021)Costa, Brandt &
  Picano]costa_brandt_picano_2021a
Costa, P., Brandt, L. & Picano, F. 2021
  Near-wall turbulence modulation by small inertial particles.  J.
  Fluid Mech.  922,  A9.

[Crowe et al.(1996)Crowe, Troutt & Chung]crowe-etal-1996
Crowe, C. T., Troutt, T. R. & Chung, J. N. 1996
  Numerical Models for Two-Phase Turbulent Flows.  Annu. Rev.
  Fluid Mech.  28 (1),  11–43.

[De Lillo et al.(2014)De Lillo, Cencini, Durham, Barry,
  Stocker, Climent & Boffetta]delillo-etal-2014
De Lillo, F., Cencini, M., Durham, W. M., Barry, M.,
  Stocker, R., Climent, E. & Boffetta, G. 2014
  Turbulent Fluid Acceleration Generates Clusters of Gyrotactic
  Microorganisms.  Phys. Rev. Lett.  112 (4),  044502.

[Dombre et al.(1986)Dombre, Frisch, Greene, Hénon, Mehr &
  Soward]dombre-etal-1986
Dombre, T., Frisch, U., Greene, J. M., Hénon, M.,
  Mehr, A. & Soward, A. M. 1986  Chaotic streamlines in
  the ABC flows.  J. Fluid Mech.  167,  353–391.

[Einstein(1906)]einstein-1906
Einstein, A. 1906  Eine neue Bestimmung der
  Moleküldimensionen.  Annalen der Physik  324 (2),
  289–306.

[Elghobashi(1994)]elghobashi-1994
Elghobashi, S. 1994  On predicting particle-laden turbulent
  flows.  Appl. Sci. Res.  52 (4),  309–329.

[Elghobashi & Truesdell(1993)]elghobashi-truesdell-1993
Elghobashi, S. & Truesdell, G. C. 1993  On the
  two-way interaction between homogeneous turbulence and dispersed solid
  particles. I: Turbulence modification.  Phys. Fluids A: Fluid
  Dyn.  5 (7),  1790–1801.

[Eswaran & Pope(1988)]eswaran-pope-1988
Eswaran, V. & Pope, S. B. 1988  An examination of
  forcing in direct numerical simulations of turbulence.  Comput. Fluids.
   16 (3),  257–278.

[Falkinhoff et al.(2020)Falkinhoff, Obligado, Bourgoin &
  Mininni]falkinhoff-etal-2020
Falkinhoff, F., Obligado, M., Bourgoin, M. & Mininni,
  P. D. 2020  Preferential Concentration of Free-Falling Heavy
  Particles in Turbulence.  Phys. Rev. Lett.  125 (6),
  064504.

[Ferenc & Néda(2007)]ferenc-neda-2007
Ferenc, J.-S. & Néda, Z. 2007  On the size
  distribution of Poisson Voronoi cells.  Phys. A: Stat. Mech.
  385 (2),  518–526.

[Ferrante & Elghobashi(2003)]ferrante-elghobashi-2003
Ferrante, A. & Elghobashi, S. 2003  On the physical
  mechanisms of two-way coupling in particle-laden isotropic turbulence.
  Phys. Fluids  15 (2),  315–329.

[Fiabane et al.(2012)Fiabane, Zimmermann, Volk, Pinton &
  Bourgoin]fiabane-etal-2012
Fiabane, L., Zimmermann, R., Volk, R., Pinton, J.-F.
  & Bourgoin, M. 2012  Clustering of finite-size particles in
  turbulence.  Phys. Rev. E  86 (3),  035301.

[Frisch(1995)]frisch-1995
Frisch, U. 1995 Turbulence: The Legacy of A.
  N. Kolmogorov.  Cambridge University Press.

[Gao et al.(2023)Gao, Samtaney & Richter]gao-etal-2023
Gao, W., Samtaney, R. & Richter, D. H. 2023
  Direct numerical simulation of particle-laden flow in an open channel
  at.  J. Fluid Mech.  957,  A3.

[Gore & Crowe(1989)]gore-crowe-1989
Gore, R. A. & Crowe, C. T. 1989  Effect of particle
  size on modulating turbulent intensity.  Int. J. Multiph. Flow.
  15 (2),  279–285.

[Goto & Vassilicos(2008)]goto-vassilicos-2008
Goto, S. & Vassilicos, J. C. 2008  Sweep-Stick
  Mechanism of Heavy Particle Clustering in Fluid Turbulence.
  Phys. Rev. Lett.  100 (5),  054503.

[Gustavsson & Mehlig(2011)]gustavsson-mehlig-2011
Gustavsson, K. & Mehlig, B. 2011  Distribution of
  relative velocities in turbulent aerosols.  Physical Review E
  84 (4),  045304, publisher: American Physical Society.

[Homann & Bec(2010)]homann-bec-2010
Homann, H. & Bec, J. 2010  Finite-size effects in the
  dynamics of neutrally buoyant particles in turbulent flow.  J. Fluid
  Mech.  651,  81–91.

[Hoque et al.(2016)Hoque, Mitra, Sathe, Joshi &
  Evans]hoque-etal-2016
Hoque, M. M., Mitra, S., Sathe, M. J., Joshi, J. B. &
  Evans, G. M. 2016  Experimental investigation on modulation of
  homogeneous and isotropic turbulence in the presence of single particle using
  time-resolved PIV.  Chem. Eng. Sci.  153,  308–329.

[Hori et al.(2022)Hori, Rosti &
  Takagi]hori-rosti-takagi-2022
Hori, N., Rosti, M. E. & Takagi, S. 2022  An
  Eulerian-based immersed boundary method for particle suspensions with
  implicit lubrication model.  Comput. Fluids  236,  105278.

[Huang et al.(2007)Huang, Shin & Sung]huang-etal-2007
Huang, W.-X., Shin, S. J. & Sung, H. J. 2007
  Simulation of flexible filaments in a uniform flow by the immersed
  boundary method.  J. Comput. Phys.  226 (2),  2206–2228.

[Hwang & Eaton(2006)]hwang-eaton-2006
Hwang, W. & Eaton, J. K. 2006  Homogeneous and
  isotropic turbulence modulation by small heavy particles.  J. Fluid
  Mech.  564,  361–393.

[Johnson & Patel(1999)]johnson-patel-1999
Johnson, T. A. & Patel, V. C. 1999  Flow past a
  sphere up to a Reynolds number of 300.  J. Fluid Mech.  378,
   19–70.

[Kajishima et al.(2001)Kajishima, Takiguchi, Hamasaki &
  Miyake]kajishima-etal-2001
Kajishima, T., Takiguchi, S., Hamasaki, H. & Miyake,
  Y. 2001  Turbulence Structure of Particle-Laden Flow in a
  Vertical Plane Channel Due to Vortex Shedding.  JSME Int. J.
  B-Fluid T.  44 (4),  526–535.

[Kempe & Fröhlich(2012)]kempe-etal-2012
Kempe, T. & Fröhlich, J. 2012  An improved
  immersed boundary method with direct forcing for the simulation of particle
  laden flows.  J. Comput. Phys.  231 (9),  3663–3684.

[Kidanemariam et al.(2013)Kidanemariam, Chan-Braun, Doychev
  & Uhlmann]kidanemariam-etal-2013
Kidanemariam, A. G., Chan-Braun, C., Doychev, T. &
  Uhlmann, M. 2013  Direct numerical simulation of horizontal
  open channel flow with finite-size, heavy particles at low solid volume
  fraction.  New J. Phys.  15 (2),  025031.

[Kolmogorov(1941a)]kolmogorov-dissipation-1941
Kolmogorov, A.N. 1941a Dissipation of
  energy in locally isotropic turbulence.  In Dokl. Akad. Nauk
  SSSR, ,  vol. 32,  pp. 16–18.

[Kolmogorov(1941b)]kolmogorov-1941
Kolmogorov, A.N. 1941b  The Local
  Structure of Turbulence in an Incompressible Viscous Fluid for
  Very Large Reynolds Numbers.  Dokl. Akad. Nauk. SSSR  30,
  301–305.

[Kolmogorov(1962)]kolmogorov-1962
Kolmogorov, A.N. 1962  A Refinement of Previous
  Hypothesis Concerning the Local Structure of Turbulence in a
  Viscous Incompressible Fluid at Large Reynolds Number.  J. Fluid
  Mech.  13,  82–85.

[Krieger & Dougherty(1959)]krieger-dougherty-1959
Krieger, I. M. & Dougherty, T. J. 1959  A
  Mechanism for Non-Newtonian Flow in Suspensions of Rigid
  Spheres.  Trans. Soc. Rheol.  3 (1),  137–152.

[Lin et al.(2017)Lin, Yu, Shao & Wang]lin-etal-2017
Lin, Z., Yu, Z., Shao, X. & Wang, L.-P. 2017
  Effects of finite-size neutrally buoyant particles on the turbulent flows
  in a square duct.  Phys. Fluids  29 (10),  103304.

[Lucci et al.(2010)Lucci, Ferrante &
  Elghobashi]lucci-etal-2010
Lucci, F., Ferrante, A. & Elghobashi, S. 2010
  Modulation of isotropic turbulence by particles of Taylor
  length-scale size.  J. Fluid Mech.  650,  5–55.

[Lucci et al.(2011)Lucci, Ferrante &
  Elghobashi]lucci-etal-2011
Lucci, F., Ferrante, A. & Elghobashi, S. 2011
  Is Stokes number an appropriate indicator for turbulence modulation
  by particles of Taylor-length-scale size?  Phys. Fluids
  23 (2),  025101.

[Maxey(2017)]maxey-2017
Maxey, M. 2017  Simulation Methods for Particulate
  Flows and Concentrated Suspensions.  Annu. Rev. Fluid Mech.
  49 (1),  171–193.

[Maxey(1987)]maxey-1987
Maxey, M. R. 1987  The gravitational settling of aerosol
  particles in homogeneous turbulence and random flow fields.  J. Fluid
  Mech.  174,  441–465, publisher: Cambridge University.

[Maxey & Riley(1983)]maxey-riley-1983
Maxey, M. R. & Riley, J. J. 1983  Equation of motion
  for a small rigid sphere in a nonuniform flow.  Phys. Fluids
  26 (4),  883–889.

[Mittal(2000)]mittal-2000
Mittal, R. 2000  Response of the Sphere Wake to
  Freestream Fluctuations.  Theoret. Comput. Fluid Dynamics
  13 (6),  397–419.

[Monchaux et al.(2010)Monchaux, Bourgoin &
  Cartellier]monchaux-etal-2010
Monchaux, R., Bourgoin, M. & Cartellier, A. 2010
  Preferential concentration of heavy particles: A Voronoï analysis.
   Phys. Fluids  22 (10),  103304.

[Monchaux et al.(2012)Monchaux, Bourgoin &
  Cartellier]monchaux-etal-2012
Monchaux, R., Bourgoin, M. & Cartellier, A. 2012
  Analyzing preferential concentration and clustering of inertial particles
  in turbulence.  Int. J. Multiph. Flow  40,  1–18.

[Monti et al.(2021)Monti, Rathee, Shen &
  Rosti]monti-etal-2021
Monti, A., Rathee, V., Shen, A. Q. & Rosti, M. E.
  2021  A fast and efficient tool to study the rheology of dense
  suspensions.  Phys. Fluids  33 (10),  103314.

[Brändle de Motta et al.(2016)Brändle de Motta, Estivalezes,
  Climent & Vincent]brandle-etal-2016
Brändle de Motta, J. C., Estivalezes, J. L., Climent, E.
  & Vincent, S. 2016  Local dissipation properties and
  collision dynamics in a sustained homogeneous turbulent suspension composed
  of finite size particles.  International Journal of Multiphase Flow
  85,  369–379.

[Oka & Goto(2022)]oka-goto-2022
Oka, S. & Goto, S. 2022  Attenuation of turbulence in
  a periodic cube by finite-size spherical solid particles.  J. Fluid
  Mech.  949,  A45.

[Olivieri et al.(2022)Olivieri, Cannon &
  Rosti]olivieri-cannon-rosti-2022
Olivieri, S., Cannon, I. & Rosti, M. E. 2022
  The effect of particle anisotropy on the modulation of turbulent flows.
  J. Fluid Mech.  950,  R2.

[Peng et al.(2019)Peng, Ayala & Wang]peng-etal-2019
Peng, C., Ayala, O. M. & Wang, L.-P. 2019  A
  direct numerical investigation of two-way interactions in a particle-laden
  turbulent channel flow.  J. Fluid Mech.  875,  1096–1144.

[Peng et al.(2023)Peng, Sun & Wang]peng-etal-2023
Peng, C., Sun, Q. & Wang, L.-P. 2023
  Parameterization of turbulence modulation by finite-size solid particles
  in forced homogeneous isotropic turbulence.  J. Fluid Mech.
  963,  A6.

[Petersen et al.(2019)Petersen, Baker &
  Coletti]petersen-etal-2019
Petersen, A. J., Baker, L. & Coletti, F. 2019
  Experimental study of inertial particles clustering and settling in
  homogeneous turbulence.  J. Fluid Mech.  864,  925–970.

[Podvigina & Pouquet(1994)]podvigina-pouquet-1994
Podvigina, O. & Pouquet, A. 1994  On the non-linear
  stability of the 1:1:1 ABC flow.  Phys. D: Nonlinear Phenom.
  75 (4),  471–508.

[Pope(2000)]pope-2000
Pope, S.B. 2000 Turbulent Flows.  Cambridge
  University Press, Cambridge.

[Rodriguez et al.(2011)Rodriguez, Borell, Lehmkuhl, Segarra &
  Oliva]rodriguez-etal-2011
Rodriguez, I., Borell, R., Lehmkuhl, O., Segarra,
  C. D. & Oliva, A. 2011  Direct numerical simulation of the
  flow over a sphere at Re = 3700.  J. Fluid Mech.  679,
  263–287.

[Rosti & Brandt(2020)]rosti_brandt_2020a
Rosti, M. E. & Brandt, L. 2020  Increase of turbulent
  drag by polymers in particle suspensions.  Phys. Rev. Fluids
  5 (4),  041301.

[Saffman(1965)]saffman-1964
Saffman, P. G. 1965  The lift on a small sphere in a slow
  shear flow.  Journal of Fluid Mechanics  22 (2),  385–400,
  publisher: Cambridge University Press.

[Salazar et al.(2008)Salazar, Jong, Cao, Woodward, Meng &
  Collins]salazar-etal-2008
Salazar, J. P. L. C., Jong, J. D., Cao, L., Woodward,
  S. H., Meng, H. & Collins, L. R. 2008  Experimental and
  numerical investigation of inertial particle clustering in isotropic
  turbulence.  J. Fluid Mech.  600,  245–256.

[Saw et al.(2008)Saw, Shaw, Ayyalasomayajula, Chuang &
  Gylfason]saw-etal-2008
Saw, E. W., Shaw, R. A., Ayyalasomayajula, S., Chuang,
  P. Y. & Gylfason, Á. 2008  Inertial Clustering of
  Particles in High-Reynolds-Number Turbulence.  Phys. Rev. Lett.
   100 (21),  214501.

[Sengupta et al.(2017)Sengupta, Carrara &
  Stocker]sengupta-etal-2017
Sengupta, A., Carrara, F. & Stocker, R. 2017
  Phytoplankton can actively diversify their migration strategy in response
  to turbulent cues.  Nature  543 (7646),  555–558.

[Shao et al.(2012)Shao, Wu & Yu]shao-etal-2012
Shao, X., Wu, T. & Yu, Z. 2012  Fully resolved
  numerical simulation of particle-laden turbulent flow in a horizontal channel
  at a low Reynolds number.  J. Fluid Mech.  693,
  319–344.

[She & Leveque(1994)]she-leveque-1994
She, Z.-S. & Leveque, E. 1994  Universal scaling laws
  in fully developed turbulence.  Phys. Rev. Lett.  72 (3),
  336–339.

[Shen et al.(2022)Shen, Peng, Wu, Chong, Lu &
  Wang]shen-etal-2022
Shen, J., Peng, C., Wu, J., Chong, K. L., Lu, Z.
  & Wang, L.-P. 2022  Turbulence modulation by finite-size
  particles of different diameters and particle–fluid density ratios
  in homogeneous isotropic turbulence.  J. Turb.  23 (8),
  433–453.

[Tanaka & Eaton(2010)]tanaka-eaton-2010
Tanaka, T. & Eaton, J. K. 2010  Sub-Kolmogorov
  resolution partical image velocimetry measurements of particle-laden forced
  turbulence.  J. Fluid Mech.  643,  177–206.

[Ten Cate et al.(2004)Ten Cate, Derksen, Portela & van
  Den Akker]tenCate-etal-2004
Ten Cate, A., Derksen, J. J., Portela, L. M. & van
  Den Akker, H. E. a. 2004  Fully resolved simulations of colliding
  monodisperse spheres in forced isotropic turbulence.  J. Fluid Mech.
  519,  233–271.

[Tsuji et al.(1993)Tsuji, Kawaguchi & Tanaka]tsuji-etal-1993
Tsuji, Y., Kawaguchi, T. & Tanaka, T. 1993
  Discrete particle simulation of two-dimensional fluidized bed.
  Powder Technol.  77 (1),  79–87.

[Tsuji & Morikawa(1982)]tsuji-morikawa-1982
Tsuji, Y. & Morikawa, Y. 1982  LDV measurements
  of an air–-solid two-phase flow in a horizontal pipe.  J.
  Fluid Mech.  120,  385–409.

[Tsuji et al.(1984)Tsuji, Morikawa & Shiomi]tsuji-etal-1984
Tsuji, Y., Morikawa, Y. & Shiomi, H. 1984
  LDV measurements of an air-solid two-phase flow in a vertical pipe.
  J. Fluid Mech.  139,  417–434.

[Uhlmann(2005)]uhlmann-2005
Uhlmann, M. 2005  An immersed boundary method with direct
  forcing for the simulation of particulate flows.  J. Comput. Phys.
  209 (2),  448–476.

[Uhlmann(2008)]uhlmann-2008
Uhlmann, M. 2008  Interface-resolved direct numerical
  simulation of vertical particulate channel flow in the turbulent regime.
  Phys. Fluids  20 (5),  053305.

[Uhlmann & Chouippe(2017)]uhlmann-chouippe-2017
Uhlmann, M. & Chouippe, A. 2017  Clustering and
  preferential concentration of finite-size particles in forced
  homogeneous-isotropic turbulence.  J. Fluid Mech.  812,
  991–1023.

[Uhlmann & Doychev(2014)]uhlmann-doychev-2014
Uhlmann, M. & Doychev, T. 2014  Sedimentation of a
  dilute suspension of rigid spheres at intermediate Galileo numbers: The
  effect of clustering upon the particle motion.  J. Fluid Mech.
  752,  310–348.

[Wang et al.(2018)Wang, Abbas &
  Climent]wang-abbas-climent-2018
Wang, G., Abbas, M. & Climent, E. 2018
  Modulation of the regeneration cycle by neutrally buoyant finite-size
  particles.  Journal of Fluid Mechanics  852,  257–282,
  publisher: Cambridge University Press.

[Wang et al.(2015)Wang, Peng, Guo & Yu]wang-etal-2015
Wang, L.-P., Peng, C., Guo, Z. & Yu, Z. 2015
  Flow Modulation by Finite-Size Neutrally Buoyant Particles in a
  Turbulent Channel Flow.  J. Fluids Eng.  138 (4).

[Yang & Shy(2005)]yang-shy-2005
Yang, T. S. & Shy, S. S. 2005  Two-way interaction
  between solid particles and homogeneous air turbulence: Particle settling
  rate and turbulence modification measurements.  J. Fluids Eng.
  526,  171–216.

[Yeo et al.(2010)Yeo, Dong, Climent & Maxey]yeo-etal-2010
Yeo, K., Dong, S., Climent, E. & Maxey, M. R.
  2010  Modulation of homogeneous turbulence seeded with finite size
  bubbles or particles.  Int. J. Multiph. Flow  36 (3),
  221–233.

[Yousefi et al.(2020)Yousefi, Ardekani &
  Brandt]yousefi-etal-2020
Yousefi, A., Ardekani, M. N. & Brandt, L. 2020
  Modulation of turbulence by finite-size particles in statistically
  steady-state homogeneous shear turbulence.  J. Fluid Mech.  899,
   A19.

[Yun et al.(2006)Yun, Kim & Choi]yun-etal-2006
Yun, G., Kim, D. & Choi, H. 2006  Vortical
  structures behind a sphere at subcritical Reynolds numbers.  Phys.
  Fluids  18 (1),  015102.

[Zahtila et al.(2023)Zahtila, Chan, Ooi &
  Philip]zahtila-etal-2023
Zahtila, T., Chan, L., Ooi, A. & Philip, J.
  2023  Particle transport in a turbulent pipe flow: Direct numerical
  simulations, phenomenological modelling and physical mechanisms.  J.
  Fluid Mech.  957,  A1.

[Zeng et al.(2009)Zeng, Najjar, Balachandar &
  Fischer]zeng-etal-2009
Zeng, L., Najjar, F., Balachandar, S. & Fischer, P.
  2009  Forces on a finite-sized particle located close to a wall in a
  linear shear flow.  Phys. Fluids  21 (3),  033302.



