Kolmogorov-size particles in homogeneous and isotropic turbulence
Alessandro Chiarini Simone Tandurella Marco Edoardo Rosti

UAV-Mounted Movable Antenna: Joint Optimization of UAV Placement and Antenna Configuration
    Xiao-Wei Tang, Yunmei Shi, Yi Huang, and Qingqing Wu
Xiao-Wei Tang, Yunmei Shi, Yi Huang ({xwtang, ymshi, and huangyi718b}@tongji.edu.cn) are with the Department of Information and Communication Engineering, Tongji University, Shanghai, China.
Qingqing Wu (qingqingwu@sjtu.edu.cn) is with the Department of Electronic Engineering, Shanghai Jiao Tong University, Shanghai 200240, China.

    3 September 2024
======================================================================================================================================================================================================================================================================================================================================================================================================




§ ABSTRACT 

We investigate the fluid-solid interaction of suspensions of Kolmogorov-size spherical particles moving in homogeneous isotropic turbulence at a microscale Reynolds number of Re_λ≈ 140. Two volume fractions are considered, 10^-5 and 10^-3, and the solid-to-fluid density ratio is set to 5 and 100. 
We present a comparison between interface-resolved (PR-DNS) and one-way-coupled point-particle (PP-DNS) direct numerical simulations. 
We find that the modulated energy spectrum shows the classical -5/3 Kolmogorov scaling in the inertial range of scales and a -4 scaling at smaller scales, with the latter resulting from a balance between the energy injected by the particles and the viscous dissipation, in an otherwise smooth flow. An analysis of the small-scale flow topology shows that the particles mainly favour events with axial strain and vortex compression. The dynamics of the particles and their collective motion studied for PR-DNS are used to assess the validity of the PP-DNS. 
We find that the PP-DNS predicts fairly well both the Lagrangian and Eulerian statistics of the particles motion for the low-density case, while some discrepancies are observed for the high-density case. Also, the PP-DNS is found to underpredict the level of clustering of the suspension compared to the PR-DNS, with a larger difference for the high-density case. 








§ INTRODUCTION



Particle-laden turbulent flows have been extensively investigated over the years, because of their relevance from both the fundamental and applicative viewpoints. They are indeed ubiquitous in several natural and engineering scenarios <cit.>, such as in volcanic ash and cloud droplets in atmospheric turbulence, dust particles in protoplanetary disks, sandstorms, ocean microplastics, and fuel droplets in spray combustion.



 §.§ The Maxey-Riley-Gatignol (MRG) Equation


In particle-laden turbulent flows the turbulent scales of the fluid phase are coupled in a non trivial manner with the solid phase. Properly resolving the flow around each particle is thus crucial to capture the fluid-solid interaction and describe the dynamics of the particles. Due to the prohibitive computational cost, however, most of the theoretical and numerical studies rely on approximations and models. 
As an example, in the context of particle clustering and preferential sampling — i.e. the tendency of particles to explore flow regions with specific properties — we mention the recent theoretical works by <cit.>, <cit.>, <cit.> and <cit.>. Most of the models used are based on the seminal works of <cit.> and <cit.>, where the equation for small rigid spheres in a non-uniform flow (hereafter referred to as MRG equation) has been derived by exploiting linear perturbation theory. In these models, each particle is treated as a mathematical point source of mass, momentum and energy. In point-particle models, particles are indeed assumed to be much smaller than any structure of the flow, as the MRG equation holds when the fluid velocity field does not show a turbulent behaviour at the particle scale. In other words, the Reynolds number based on the particle diameter and the particle-fluid relative velocity has to be small, i.e. 𝒪(1).

The models based on the MRG equation do not resolve the flow around the particles, and the influence of the solid phase on the fluid phase has to be modelled; see for example <cit.>, <cit.> and <cit.>. However, the back-reaction of the particles on the carrier flow and the inter-particle collisions are usually negligible in the limit of very dilute regimes with Φ_V = V_s/(V_s + V_f) ≤ 10^-5 (where Φ_V is the volume fraction, and V_s and V_f are the volumes of the solid and fluid phases, respectively), small particles D_p ≤η (where η is the turbulent Kolmogorov scale) and small solid-to-fluid density ratios ρ_p/ρ_f <cit.>. In this case, the influence of the solid phase on the carrier fluid can be neglected and the fluid-particle interaction is often modelled with one-way coupling models, in which the particles move under the action of the flow, but they do not modulate it. Starting from the MRG equation, several corrections have been proposed over the years to account for several effects, and extend its range of validity to a wider range of parameters. For example, <cit.> introduced a lift force which is crucial to properly model the particle dynamics in the presence of a linear shear. This force has been later extended by <cit.> to fit the numerical data of <cit.> at large Reynolds numbers. Other corrections have been introduced to account for finite values of the particle Reynolds number. For example, we refer to the corrections to the drag term reported in <cit.>, and to the different convolutional kernels for the Basset time-history force contribution proposed by <cit.>. 

The point particle approximation coupled with Direct Numerical Simulations of the Navier–Stokes equations (PP-DNS) has been widely used to investigate turbulent particle-laden flow in the one-way coupling regime in several scenarios <cit.>. Despite the large number of studies based on PP-DNS, however, clear understanding of the range of validity of the underlying assumptions of the point-particle model is not yet present, and studies that investigate its limitations are needed. In this respect, a step forward has been done in the last years thanks to the introduction of several numerical methods (PR-DNS) which resolve the flow around each particle by coupling the direct numerical simulations of the Navier–Stokes equations with the immersed boundary method (IBM); see for example <cit.>, <cit.>, <cit.>, <cit.>, <cit.> and <cit.>. Unlike experiments, PR-DNS allows to explain the underlying physical mechanisms and investigate the limitations of the point-particle approximation.

In the context of forced homogeneous isotropic turbulence (HIT) with fixed particles, <cit.> and <cit.> found that, when compared with PR-DNS, the two-way coupled PP-DNS based on the Schiller-Naumann drag correlation captures fairly well the turbulence attenuation and the fraction of the turbulence dissipation rate due to the particles. 
Surprisingly, they found a good agreement also for D_p/η = 1.5 <cit.> and D_p/η≈ 2.2 <cit.>, although the model is expected to largely under-perform for these particle sizes. <cit.> used PR-DNS to assess the validity of PP-DNS in a decaying isotropic turbulent particle-laden flow, focusing on the particle acceleration model. They found that the predictions of the PP-DNS models they considered are in excellent agreement compared with PR-DNS for small Stokes numbers. For large Stokes numbers, however, they found that PP-DNS under-predicts the true particle acceleration and that second moment quantities are not properly captured. They showed that the predictions improve once considering finite Reynolds number corrections to the model. <cit.> tested the one-way point-particle approximation in a turbulence channel flow laden with small inertial particles, with high particle-to-fluid density ratios. They considered a volume fraction of Φ_V ≈ 10^-5 to ensure that the feedback of the particles on the fluid phase was negligible. They found that in the bulk of the channel the model predicts fairly well the statistic of the particles velocity. Close to the wall, however, they observed that the model fails, as it is not able to capture the shear-induced lift force acting on the particles, which instead is well predicted by the lift force model introduced by <cit.>. 

In this work, we further address the limit of the one-way point-particle approximation, and we use PR-DNS to investigate the reliability of one-way coupled PP-DNS in the context of forced homogeneous isotropic turbulence laden with Kolmogorov-size particles.



 §.§ Particles in homogeneous isotropic turbulence


The fluid-solid interaction of suspension of spherical particles moving in homogeneous isotropic turbulence has been widely investigated over the years by means of both simulations and experiments. In the following, we list some of the main contributions, and we refer the interested reader to <cit.> and <cit.> for a comprehensive review. 

The majority of the numerical studies dealing with small particles D_p ≪η are based on the point-particle approximation, as particle-resolved methods are prohibitively expensive in this case. Although the related model error is not known, the available numerical studies have shown that small particles may either amplify or damp the turbulence of the carrier flow, in agreement with previous experimental studies <cit.>. <cit.> used two-way coupled PP-DNS to investigate homogeneous isotropic turbulence laden with small heavy spherical particles at Re_λ = u' λ / ν≈ 38 (where u' is the average velocity fluctuation, λ is the Taylor length scale, and ν is the kinematic fluid viscosity). 
Compared to the particle-free case, they reported a significant attenuation of the fluid kinetic energy and of the dissipation rate. By looking at the energy spectrum, they found that the addition of the particles results into a relative enhancement of the energy at the small scales compared to the energy content at the large scales. They also showed that heavier particles cause a less selective modification of the turbulence properties. Heavy particles are indeed more uniformly dispersed by the turbulence, and cause a more homogeneous modification of the flow properties compared to lighter particles, that instead show a stronger preferential collection in regions of low vorticity and high strain <cit.>. <cit.> used two-way coupled PP-DNS to investigate decaying homogeneous isotropic turbulence laden with small spherical particles. Besides confirming the non-uniform modulation of the energy spectrum, they observed that the energy enhancement at the small scales is accompanied by an increase of the viscous dissipation rate and, thus, by an enhancement of the rate of energy transfer from larger to smaller scales. <cit.> made use of two-way coupled PP-DNS to study the influence of small and heavy particles on forced homogeneous isotropic turbulence at Re_λ = 62. 
They reported that the influence of the particles changes with their inertia, with the small-scale energy content being attenuated/enhanced by large/small particles. By investigating the spectrum of the fluid-particle energy exchange rate, they observed that particles act as a sink of kinetic energy at large scales, while they add kinetic energy to turbulence at the smallest scales. The large scale motions of the fluid drag the particles, while the small-scale fluctuations are driven by the presence of the solid phase.
<cit.> investigated the influence of small and heavy particles in decaying isotropic homogeneous turbulence at the initial Reynolds number of Re_λ = 30 and Re_λ = 50, with a focus on 
particles with very small inertia and small relaxation time. Unlike the previous works, they found that the turbulent kinetic energy and the viscous dissipation rate increase at all times compared to the particle-free case. The presence of the particles, indeed, largely enhances the small-scale energy content, while slightly reduces the large-scale energy content, with positive integral variation. 

This was later confirmed by <cit.>, that investigated by PP-DNS the influence of particles on decaying homogeneous isotropic turbulence with an initial Reynolds number of Re_λ = 75. 













For large particles with D_p>η the numerical studies are based on PR-DNS, as in this case the point-particle approximation does not hold. <cit.> and <cit.> investigated the influence of particles with size D_p ≈λ (where λ is the Taylor length scale) in decaying homogeneous turbulence. They observed that in contrast to what happens when D_p<η, the presence of the particles damps the turbulent kinetic energy of the fluid compared to the particle-free case at all times, and that the two-way coupling rate of change is always positive. <cit.> investigated the influence of particles with a solid-to-fluid density ratio of ρ_p/ρ_f = 1.15 and 1.73 on forced homogeneous isotropic turbulence at Re_λ = u' λ / ν = 61, varying the volume fraction between Φ_V = 0.02 and Φ_V = 0.1. They found that 
the energy spectrum is enhanced for wavenumbers κ > κ_p ≈ 0.75 κ_D, where κ_D = 2 π /D_p, while it is attenuated for κ < κ_p. These results were later confirmed by <cit.> at Re_λ≈ 60, ρ_p/ρ_f = 1.4 and Φ_V = 0.06. <cit.> considered particles with D_p/η≈ 5-8 and ρ_p/ρ_f = 1.5 at the larger Reynolds number of Re_λ≈ 130. They focused on the dynamics of the particles, and observed that finite-size inertial particles exhibit a moderate level of clustering, as later confirmed also by <cit.>. <cit.> considered the effect of particles on homogeneous isotropic turbulence at the larger Reynolds number of Re_λ≈ 400, which ensures a well developed inertial range of scales. They set the volume fraction at Φ_V = 0.079, and investigated the turbulence modulation by particles with size D_p/η=123 and solid-to-fluid density ratio between 1.3 ≤ρ_p/ρ_f ≤ 100. They showed that the solid phase modifies the energy cascade described by Richardson and Kolmogorov; the fluid-solid coupling drives the energy cascade at large scales, while the classical energy cascade is restored at scales smaller than the particle size. <cit.> studied the turbulence modulation due to spherical particles with 7.8 ≤ D_p/η≤ 64 by setting the volume fraction at Φ_V = 8.1 × 10^-3 and the reference Reynolds number at Re_λ≤ 100. They found that the turbulent kinetic energy content monotonically decreases with D_p, due to the increase of the energy dissipation rate in the wake of the particles. More recently, <cit.> and <cit.> investigated by PR-DNS how the flow modulation changes with D_p and ρ_p. They set the Reynolds number to Re_λ≈ 400 and the volume fraction to Φ_V = 0.079, and varied the particles size and the solid-to-fluid density ratio in the 16 ≤ D_p/η≤ 123 and 1.3 ≤ρ_p/ρ_f ≤ 100 range. <cit.> found that, in presence of an inhomogeneous mean shear, particles might modulate the largest scales of the flow towards an anisotropic state, while preserving homogeneity and isotropy at the smaller scales. <cit.> observed that interface-resolved particles enhance flow intermittency favouring events with large localised velocity gradients. For the smallest and heaviest particles, they found that the classical energy cascade is subdominant at all scales, and that the energy transfer is completely driven by the fluid-solid coupling term. <cit.> investigated the effect of the Reynolds number on the flow modulation by finte-size particles in homogeneous isotropic turbulence. Notably, they observed that the modulation of the turbulent kinetic energy has little dependence on Re_λ, and that particles modulate turbulence also at the smallest Reynolds numbers.

While a relatively larger body of literature has investigated the dynamics of particle of size larger and smaller than Kolmogorov size, fewer works have considered Kolmogorov-size particles with D_p ≈η, which are the focus of the present work. From an experimental point of view the D_p ≈η case is complex, as it requires a resolution of sub-Kolmogorov scales when measuring the velocity perturbations near the particles <cit.>. Numerical schemes based on the MRG equation, which are commonly used for D_p ≪η, are generally thought of as not valid when D_p ≈η <cit.>. PR-DNS, on the other hand, becomes prohibitively expensive as D_p decreases when Re is sufficiently large, due to the extra resolution required to properly resolve both the flow perturbations induced by the particles and all the turbulence scales. Among the few works available, we mention <cit.> that experimentally investigated the influence of a dilute dispersion of particles with D_p ≈η in forced homogeneous isotropic turbulence at Re_λ≈ 230. They observed that Kolmogorov-size particles attenuate the turbulent global kinetic energy and the viscous dissipation rate up to 40% and 50% for a mass loading Φ_M = ρ_p V_p/(ρ_f V_f) of Φ_M = 0.3. <cit.> studied by PR-DNS the interaction of decaying isotropic turbulence with finite-size D_p ≈η particles. They set the initial Reynolds number of the flow at Re_λ = 79, and varied the solid-to-fluid density ratio between 40 ≤ρ_p/ρ_f ≤ 5000 and the mass loading between 0.01 ≤Φ_M ≤ 1. They observed that in the vicinity of the particles the viscous dissipation rate of the fluid is amplified due to the large velocity gradients that are generated by the boundary conditions at the surface of the particles <cit.>. Particles also release kinetic energy to the fluid by locally accelerating the surrounding flow, similarly to what seen by <cit.> for larger particles. From a global viewpoint, <cit.> observed that, for large ρ_p/ρ_f, particles with D/η≈ 1 induce local velocity disturbances that significantly modulate the distribution and the decay of the fluid kinetic energy at all scales. 




 §.§ Present study


In this study, we investigate the fluid-solid interaction of suspension of Kolmogorov-size spherical particles moving in homogeneous isotropic turbulence at the relatively large microscale Reynolds number of Re_λ≈ 140 by use of direct numerical simulations. The study is based on both PR-DNS and PP-DNS. The objective of the present study is twofold; we aim (i) to investigate the modulation of forced homogeneous isotropic turbulence by finite Kolmogorov-size particles at a Reynolds number which is large enough to ensure a proper separation of scales, and (ii) to address the limits and the range of validity of the one-way-coupled PP-DNS in the simplest configuration of homogeneous isotropic turbulence. To do this, we consider a portion of the parameter space which is on the edge of the range of validity of the one-way-coupled PP-DNS <cit.>. 


The structure of the work is as follows. After this introduction, the computational set up and the numerical methods are described in <ref>. Then, section <ref> is devoted to the assessment of the flow modulation, and discusses the results of the PR-DNS. The influence of the particles on the energy spectrum, on the scale-by-scale energy budget and on the local structure of the flow are discussed. Sections <ref> and <ref> deal respectively with the dynamics of the particles and with the inhomogeneity of their distribution in the flow. In these sections, we assess the validity of the one-way coupled PP-DNS. Eventually, concluding remarks are provided in <ref>.




§ MATHEMATICAL FORMULATIONS AND NUMERICAL METHOD



We consider a turbulent flow in a triperiodic box of size L=2π laden with N spherical particles; see figure <ref>. The carrier flow is governed by the incompressible Navier–Stokes equations





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


where u=(u,v,w) is the fluid velocity, p is the reduced pressure, and ρ_f and ν are the fluid density and kinematic viscosity. At the right-hand-side of the momentum equation f is the Arnold-Beltrami-Childress (ABC) cellular forcing <cit.> used to inject energy at the largest scales and sustain turbulence, while f^← p is the force the particles exert on the fluid phase.



 §.§ Particle-resolved simulations (PR-DNS)


The motion of a rigid particle can be described using the translational velocity u_p and the rotational velocity ω_p of its centre of mass, that obey the classical Newton-Euler equations for rigid body dynamics,


    m_p du_p/d t = f^← f  + f_p^↔ p,  
    I_p dω_p/d t = L_p^← f,


where m_p = πρ_p D_p^3/6 and I_p = m_p D_p^2/10 are the mass and inertial moment of the particle, with ρ_p being the particle density and D_p the particle diameter. Here f^↔ p is the force due to particle collisions, while f^← f and L_p^← f are the force and momentum due to the fluid-solid interaction, namely


    f^← f   = ∮_∂ V_pσ·ndA, L_p^← f = ∮_∂ V_pr×( σ·n) d A,


where σ= -p ℐ + 2 μ𝒟 is the Cauchy stress tensor, with ℐ being the identity tensor, μ the fluid kinematic viscosity, 𝒟 the strain rate tensor, and n the unit vector normal to the surface of the particle.



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


The interface-resolved simulations are complemented with one-way-coupled point particle simulations (PP-DNS). In this work, we consider the complete governing equation for a point particle as introduced by <cit.> and <cit.>, i.e.


    [t]
        ρ_p V_p du_p/dt   = 3 π D_p ρ_f ν(u - u_p + 1/6( D_p/2)^2 ∇^2 u) _Stokes drag
       + ρ_f V_p/2( 3Du/Dt - du_p/dt + 1/10( D_p/2)^2 d/d t∇^2 u)_Added mass
          + 3/2 D_p^2 ρ_f √(πν)∫_-∞^t K_B(t-τ) ( du/dτ - du_p/dτ+ 
           1/6( D_p/2)^2 d/d t∇^2 u) dτ_Basset force.


Here, V_p = πρ_p D_p^3/6 is the volume of the particle, and D/Dt denotes the material derivative. Note that the Faxén correction <cit.> proportional to ∇^2 u has been included in the Stokes drag, added mass and Basset forces. According to <cit.>, this correction reproduces dominant finite-size effects on velocity and acceleration fluctuations for neutrally buoyant particles with diameter up to D_p/η≈ 4. For the added mass, we have used the form described by <cit.>. The computation of the Basset history force presents some challenges, and its evaluation can become extremely time consuming and memory demanding; indeed, this term requires at each time step the computation of an integral over the complete time history of the particle. Over the years, several attempts have been proposed to approximate this term; see for example <cit.>, <cit.> and <cit.>. In this work, we resort on the second-order and memory-efficient algorithm developed by <cit.>, whose details are briefly reported in appendix <ref> for completeness.



 §.§ Computational details


We consider a single-phase micro-scale Reynolds number of Re_λ = u' λ / ν≈ 140 to ensure a relatively large inertial range of scales; here u' is the root mean square of the velocity fluctuations and λ is the Taylor length scale. The particle diameter is set to D_p/η≈ 0.9, where η is the Kolmogorov length scale for the single-phase case. Two volume fractions are considered, i.e. Φ_V = V_s/(V_s + V_f)= 10^-5 and Φ_V = 10^-3 for a total number of particles of N=742 and N=74208, respectively. For each volume fraction two values of the particle density are considered, ρ_p/ρ_f = 5 and ρ_p/ρ_f=100 to consider both light and heavy particles. This leads to a total of four PR-DNS. In PP-DNS the volume fraction is not a parameter, and only two simulations have been carried out for the different density ratios.

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

For the PR-DNS the governing equations for the particles are dealt with the immersed boundary method introduced by <cit.>. The fluid-solid coupling is achieved in an Eulerian framework, and accounts for the inertia of the fictitious fluid inside the solid phase, so as to properly reproduce the particles' behaviour in both the neutrally-buoyant case and in the presence of density difference between the fluid and solid phases. The soft sphere collision model <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 becomes large.

For the one-way-coupled PP-DNS the governing equation for the particle velocity is advanced in time using the second-order Adams-Bashforth time scheme. At each time step, the fluid velocity is evaluated at the position of the particle using a second-order linear interpolation.

For the PR-DNS the fluid domain is discretised on a cubic domain using N_p = 2048 points along each direction, leading to η / Δ x ≈ D_p /Δ x ≈ 6-7, where Δ x denotes the grid spacing. The accuracy of the results has been verified by running an additional simulation on a coarser grid with N_p=1440 for the case with Φ_V = 10^-3 and ρ_p/ρ_f = 100, resulting in a negligible difference in the scale-by-scale fluid energy spectrum and budget in figures <ref> and <ref>, and in the Lagrangian and Eulerian particles' statistics in figures <ref> and <ref>. For the single-phase case and the PP-DNS the fluid domain is discretised using N_p = 1024 points in the three direction, leading to η/Δ x ≈ 3-4. Excluding the initial transient period, all simulations are advanced in time for approximately 50 τ_f, where τ_f = ℒ/√(2E/3) is the average turnover time of the largest eddies; ℒ =π/(4E/3)∫_0^∞ℰ(κ)/κdκ is the fluid integral scale with ℰ(κ) being the energy spectrum, E(x,t) is the local and instantaneous fluid kinetic energy, and the · and · operators denotes average in space and time, respectively. 

Details of the PR-DNS are reported in table <ref>. Note that, when looking at the bulk quantities, the flow modulation due to the solid phase is rather low.






§ FLOW MODULATION



In this section we discuss the PR-DNS results and focus on the influence of the particles on the carrier flow. First, we show the influence of the solid phase on the energy spectrum, on the structure functions and on the the scale-by-scale energy budget. Next, the influence of the particles on the local structure of the small scales of the flow is addressed.



 §.§ Energy Spectrum





Figure <ref> shows the influence of the solid phase on the energy spectrum ℰ(κ), and highlights how Kolmogorov-size particles modulate the turbulent fluctuations scale-by-scale. The top and bottom panels are for Φ_V = 10^-5 and Φ_V=10^-3, respectively. The black solid line refers to the reference unladen case; note the inertial range of scales where ℰ∼κ^-5/3 extends for more than one decade of wavenumbers, confirming that the present Reynolds number is large enough to ensure a proper separation of scales. For validation purposes, in the bottom panel we also plot with symbols the energy spectrum obtained for the Φ_V=10^-3 and ρ_p/ρ_f=100 case with the coarser grid, showing good agreement with those obtained with the standard grid, thus ensuring the suitability of the chosen grid resolution (see <ref>). 

At large scales, the spectra of the particle-laden cases substantially overlap with the unladen spectrum. We indeed observe only a weak depletion of the energy content at the intermediate scales; see the insets in the two panels. At scales smaller than a certain wavenumber κ_p, the energy spectra of the particle-laden cases deviate above the reference spectrum. Solid particles enhance the energy content of the small scales by means of their wake. Notably, this mechanism is amplified as Φ_V and/or ρ_p/ρ_f increase, as conveniently visualised by the larger values of ℰ for κ>κ_p and by the shift of κ_p towards smaller wavenumbers. However, notice that due to the low values of Φ_V considered, the flow modulation is rather low for all cases, being substantially negligible for Φ_V=10^-5 and/or ρ_p/ρ_f=5.

Figure <ref> shows that the modulated energy spectrum exhibits multiscaling behaviour. The classical ℰ(κ) ∼κ^-5/3 decay in the inertial range of scales is indeed followed by a steeper decay ℰ(κ) ∼κ^-4 for κ > κ_p.
A similar steep decay has been observed in bubbly flows at scales smaller than the bubble diameter, by means of both experiments <cit.> and simulations <cit.>. 
Additionally, a similar multiscaling behaviour has also been observed in a turbulent planar Couette flow laden with small particles <cit.>, and in homogeneous isotropic turbulence laden with slender fibres <cit.>. 

In the context of bubbly flows, the emergence of the κ^-α decay with α≥ 3 has been attributed to the wakes the bubbles generate in a otherwise smooth flow <cit.>; here the velocity fluctuations produced by the bubbles are directly dissipated by viscosity <cit.>. Accordingly, this scaling has been indeed observed only when the bubble Reynolds number is large enough and is in the 10 ≤ Re_bub≤ 1000 range <cit.>. 

To compare, we computed the local particle Reynolds number using the relative velocity between the particle and the surrounding flow, and found that it is in the 0.25 ⪅ Re_p ⪅ 2.5 range (see table <ref>). Further investigation on the link between the fluctuations induced by particles and the ℰ(κ) ∼κ^-4 decay is provided in <ref> by looking at the scale-by-scale energy budget.

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




 §.§ Structure function and intermittency





We extend the analysis done in the spectral domain by computing the longitudinal structure functions defined as S_p(r) = δ u(r)^p where δ u(r) = (u(x+r) - u(x) )·r/r and r = |r|. In particular, figure <ref> plots S_2, S_4 and S_6 as a function of r for (top) Φ_V=10^-5 and (bottom) Φ_V = 10^-3. In the single-phase case, we observe that S_p ∼ r^p/3 in the inertial range of scales is in agreement with the Kolmogorov prediction <cit.>, and that S_p ∼ r^p at the small scales, as a result of the differentiability of the fluid velocity field <cit.>. Recall that although S_2(r) is commonly referred to as scale energy <cit.>, its meaning slightly differs from ℰ(κ); indeed while ℰ(κ) dκ refers to the amount of energy associated with the scale r = 2π/κ, S_2(r) can be interpreted as the amount of energy associated with scales up to r.
In agreement with the modulation of the energy spectrum, figure <ref> shows that particles enhance the energy content at small scales compared to the unladen case. The energy enhancement is more intense for larger Φ_V and ρ_p/ρ_f, and becomes more evident when considering higher order structure functions. At the same time, the presence of the particles decreases the amount of energy stored at the largest scales (as seen by the blue and red curves laying below the black one). By interacting with the vortical structures of the flow, particles drain energy from scales larger than D and reinject it back at smaller scales by means of their wake.





Structure functions are commonly employed to quantify the flow intermittency, i.e. the relevance of extreme events that are localised in space and time and break the Kolmogorov similarity hypothesis <cit.>. In figure <ref>, we use the extended self similarity introduced by <cit.>, and plot S_6/S_2^3 as a function of r. In the limit case where extreme events do not occur, the S_6 ∼ S_2^3 power law holds, i.e. S_6/S_2^3 ∼ constant, and deviations from this behaviour are a measure of the flow intermittency. Accordingly with the intrinsic intermittent nature of turbulent flows, figure <ref> shows that S_6/S_2^3 deviates from the Kolmogorov prediction also in the single-phase case, and this deviation increases in the particle-laden cases, similarly to what is found for suspension of particles with size in the inertial range of scales <cit.>. This is due to the no-slip and no-penetration boundary conditions at the particles surface that give origin to localised and intense velocity gradients. Figure <ref> shows that for ρ_p/ρ_f = 5 the deviation from the single-phase case is rather small for both Φ_V = 10^-5 and Φ_V=10^-3, in agreement with the low level of flow modulation discussed above. For ρ_p/ρ_f = 100, instead, the S_6/S_2^3 curve significantly deviates from the single-phase case at small scales. For the considered parameters, heavy Kolmogorov-size particles increase the intermittency of the small scales of the flow. 



 §.§ Scale-by-scale energy budget



We now detail the influence of the particles on the organisation of the fluctuations, by studying the scale-by-scale energy budget equation. For the present case with three homogeneous directions, the energy balance reads


    P(κ) + Π(κ) + Π_fs(κ) - D_v(κ) = 0,


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



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


where ·̂ denotes the Fourier transfrom operator, and the superscript ·^* denotes complex conjugate. The term Ĝ is the Fourier transform of the non linear term ∇· ( uu ). Note that here we integrate all the terms from κ to ∞. Π(κ) and Π_fs(κ) do not produce nor dissipate energy at any scale, but redistribute it among scales by means of the classical energy cascade and of the fluid-solid interaction. Also, note that since we integrate the viscous term from κ to ∞, D_v(0) = ϵ. For the complete derivation we refer the reader to <cit.>.



Figure <ref> plots the Π, Π_fs and D_v terms of equation <ref> as a function of κ for the four particle-laden cases investigated. For validation purposes, we also plot with symbols the terms obtained with the coarser grid for Φ_V=10^-3 and ρ_p/ρ_f=100. In agreement with the multiscaling behaviour shown in the energy spectra (see figure <ref>), the energy budgets exhibit two distinct behaviours. Energy is injected at the largest scales at a rate that equals the dissipation rate P(0) = ϵ (not shown). In the inertial range of scales κ_L < κ< κ_p, the fluid-solid coupling term is subdominant and


    Π(κ) ≈ D_v(κ) ≈ϵ.


Thus, Π(κ) ≈ϵ is constant with κ at these scales, exhibiting a plateau. In agreement with the Kolmogorov theory, the viscous effects are negligible 2 νκ^2 ℰ(κ) ≈ 0 and energy is transferred from larger to smaller scales at a rate that matches the energy injection rate ϵ. This corresponds to the range of scales where ℰ(κ) ∼κ^-5/3. Similarly to what observed in the energy spectrum, the range of scales where this relation holds shrinks as Φ_V and ρ_p/ρ_f increase. For the small scales with κ > κ_p where the spectrum shows the ℰ(κ) ∼κ^-4 decay, instead, the non linear flux is negligible Π(κ) ≈ 0. In this range of scales, the viscous effects and the fluid-solid coupling term are not negligible, and the energy budget reduces to


    Π_fs(κ) ≈ D_v(κ).


Here the fluctuations that are produced by the fluid-solid interaction are directly dissipated by viscosity. 
In agreement with the energy spectrum, the range of scales where this regime holds widens as Φ_V and/or ρ_p/ρ_f are increased. 



 §.§ The local structure of the flow



As shown in <ref>, Kolmogorov-size particles mainly modify the organisation of the velocity fluctuations at the smallest scales. Particles indeed modulate the energy spectrum and the structure functions for κ⪆κ_p only. Here we characterise the smallest scales of the flow to provide new insights of the influence of the dispersed phase on the structure of the velocity fluctuations. For the sake of brevity, in this section we consider the Φ_V = 10^-3 cases only, 
and we investigate how particles modify the velocity gradient field A_ij = ∂ u_j/∂ x_i. In the neighbourhood of a given point (x_0,t), the velocity field can be approximated as u_i(x,t) = u_i(x_0,t) + A_ij(x_0,t)(x_j-x_0,j) + 𝒪(|x-x_0|^2).

This linear expansion is valid in the region around x_0, where the fluid is sufficiently smooth and the variations of A_ij are small <cit.>; for a turbulent flow the extent of this region is of the order of the Kolmogorov scale η. Based on these arguments, we study the influence of the particles on the smallest scales of the flow, by inspecting their effect on the A_ij tensor.

We decompose A_ij into its symmetric and antisymmetric parts, namely the strain-rate tensor S_ij = (A_ij + A_ji)/2 and the rotation rate tensor W_ij = (A_ij - A_ji)/2. The field of the velocity gradient is completely addressed when knowing: (i) the three principal rates of strain α≥β≥γ, i.e. the three eigenvalues of S_ij, (ii) the magnitude of the vorticity ω^2 = ω·ω, i.e. the enstrophy, and (iii) the orientation of ω relative to the three principal axes of strain, i.e. the eigenvectors of S_ij <cit.>. Note that, due to the incompressibility constraint α + β + γ = 0, meaning that α is always nonnegative, γ is always nonpositive while β can have any sign depending on the local straining state.



We start by characterising α, β and γ in figure <ref>. Following <cit.>, we use s^* which is defined as


    s^* = - 3 √(6)αβγ/(α^2 + β^2 + γ^2)^3/2.


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


In the absence of particles, the distribution peaks at s^*=1, in agreement with the fact that for purely Newtonian turbulence the most likely state of straining is an axisymmetric extension <cit.>. Figure <ref> shows that the addition of Kolmogorov-size particles with Φ_V ≤ 10^-3 leads to a rather small variation of the distribution of s^* for the lighter particles with ρ_p/ρ_f = 5. We observe instead that particles with ρ_p/ρ_f = 100 decrease the probability of events with large positive s^* and increase the probability of events with s^* ≤ 0. This agrees with the observation of <cit.> who found that particles with size in the inertial range favour the occurrence of events with two-dimensional straining states and with axisymmetric compression. These events indeed are associated with the shear layers that separate from the surface of the particles, that strengthen as ρ_p/ρ_f increases. 



Figure <ref> shows the influence of the solid phase on (top) the square of the vorticity magnitude ω^2, i.e. the enstrophy, and on (bottom) the alignment between the vorticity ω and the principal axes of strain. The distribution of ω^2 shows that the tail becomes longer and the probability of large events increases due to the presence of the particles. This agrees with the enhanced flow intermittency discussed in <ref>. The tail of the distribution is longer for the ρ_p/ρ_f = 100 case, as the velocity gradients at the particles surface are more intense because of the larger relative velocity between the particles and the surrounding fluid phase. However, even the lighter particles with ρ_p/ρ_f = 5 are able to produce a non-negligible change of the distribution.

Instead, the bottom panel of figure <ref> shows that the alignment between ω and the principal axes of strain are only slightly influenced by the presence of the particles. The results for the single-phase case perfectly overlaps with that of ρ_p/ρ_f = 5. 
The presence of the heavy particles, instead, slightly reduces the alignment between ω and ê_β, as well as the anti-alignment between ω and ê_γ. Our results suggest that for D_p/η≈ 1 the perturbation field induced by the heavy particles is characterised by events where vorticity is more aligned with the direction of extension (ê_α) and compression (ê_γ), and more anti-aligned with the intermediate eigenvector (ê_β). Interestingly, this differs from what observed by <cit.>, who found the opposite trend for particles with size in the inertial range. We presume that the difference is due to the different particle Reynolds number, which results into a different kind of wake. 

We now move and consider the entire velocity gradient tensor A_ij. Any second-order tensor possesses three invariants P, Q and R, which are directly related to its eigenvalues λ by the characteristic polynomial function


    λ^3 + P λ^2 + Q λ + R = 0.


Following <cit.>, it can be shown that


    P =    α + β + γ, 
    
    Q =     - 1/2( α^2 + β^2 + γ^2 ) + ω^2/4, 
    
    R =     - 1/3( α^3 + β^3 + γ^3 ) - 1/4ω_i ω_j S_ij.


Note that, P=0 due to the incompressibility constraint. The Q and R invariants are commonly used to distinguish between regions of intense vorticity and regions of strong strain. In particular, the discriminant of equation <ref> Δ = 27/4R^2 + Q^3 is used to distinguish between regions where motions are mainly vortical (i.e. regions where Δ<0, meaning that A_ij has one real and two complex conjugate eigenvalues) or characterised by a node-saddle streamline pattern (i.e. regions where Δ >0, and all the eigenvalues are real). When Q is large and negative the strain is intense, while the vorticity is weak; in this case, R ∼ - (α^3 + β^3 + γ^3) = - αβγ <cit.>, and a positive R implies a region of biaxial strain (γ<0, α> β > 0), while a negative R implies a region of axial strain (γ < β < 0 and α>0). When instead Q is large and positive the strain is locally weak and R ∼ - ω_i ω_j S_ij. In this case, a positive R implies vortex compression, while a negative R implies vortex stretching. 






Figure <ref>, plots the joint distribution of Q and R for Φ_V = 10^-3; the black solid lines denote the left- and right-Vieillefosse tails with Δ=0. For completeness the bottom panels report the distributions of Q (left) and R (right). In absence of the particles, the Q-R joint distribution takes a tear-drop pattern, with a clear point at the right-Vieillefosse tail with R>0 and Q<0 <cit.>. The distribution is skewed towards positive Q, but rather evenly distributed among positive and negative values of R (see the bottom panels). The largest probability is observed in the second and fourth quadrants, i.e. Q<0 R>0 and Q>0 and R<0. Thus, in a purely Newtonian turbulent flow there is a strong negative correlation between Q and R, and the two most common states are vortex stretching ω_i ω_j S_ij>0, and biaxial strain αβγ <0 <cit.>. Also, the points in the Q-R map are distributed around the origin, since the mean values of Q and R are zero in homogeneous flow <cit.>. The presence of the particles enlarges the range of possible Q and R values, in agreement with the increase of the probability of events with intense velocity gradients associated with the boundary conditions at the particles' surface. In particular, Kolmogorov-size particles mainly favour events that lay in the first and third quadrant, resulting into a joint distribution that is more symmetric with respect to an inversion of the R axis (see the bottom right panel in figure <ref>). Compared to the unladen case, particles mainly promote events with axial strain αβγ > 0 (R<0 and large negative Q<0) and with vortex compression ω_i ω_j S_ij<0 (R>0 and Q>0). The probability of events with R<0 and Q<0 is particularly enhanced, as visualised in the top panel of figure <ref> by the occurrence of a point at the left-Vieillefosse tail. This is consistent with the increase of the probability of events with s^*=-1 shown in figure <ref>. 
A visual investigation of the Q and R fields around the particles (figure <ref>) helps to explain this effect by highlighting the local contribution of the particles. By comparing the two fields in the surroundings of the particles, we observe that they are both almost axisymmetric along the axis aligned with their travelling direction. However, across the velocity-normal median plane, Q is symmetric, while R is antisymmetric. This relation between the two fields implies contributions across all four quadrants of the Q-R distribution. The effect, however, is particularly apparent in the third quadrant, as the region is otherwise not explored by the single-phase flow. In particular, regions of Q<0 and R<0 are associated with axial strain and appear to be found at one of the two ends of the particles' along their travelling direction. Consistently with the above-mentioned symmetries, at the opposite end of the particles, a region of biaxial strain (Q<0 and R>0) is found.

This scenario only partially agrees with the results of <cit.> for decaying homogeneous isotropic turbulence laden with Kolmogorov-size particles. They also found that particles favour events with axial strain (Q<0 and R<0) and biaxial strain (Q<0 and R>0), as shown by the occurrence of a point at the left-Vieillefosse tail and by the more pronounced right-Vieillefosse tail in figure 20 of their paper. However, they did not report an increase of the probability of events with Q>0 like in the present case. It is worth mentioning, however, that they do report a global increase of the vortex stretching. Compared to larger particles with size in the inertial range of scales, the scenario is completely different: in fact, <cit.> found that when large particles are added both Q and R are reduced. Besides energising the small scales, large particles indeed behave also as obstacles for the large flow structures, largely weakening thus the energy content at the large scales. <cit.> only observed an increase of the probability in the strain-dominated region (that we also observe), which is an indication of the intense dissipation regions that arise around the particles.



Additional insights are provided by looking separately at the invariants of S_ij and W_ij; see figure <ref>. In particular, we consider their second invariants, i.e.


    Q_S = - 1/2( α^2 + β^2 + γ^2 )  andQ_W = 1/4ω^2.


These invariants are related with the fluid dissipation ϵ = -4 ν Q_S and with the fluid enstrophy ω^2 = 4Q_W. Therefore, the Q_S-Q_W joint distribution determines whether the flow is dominated by dissipation (extensional dominated regions with Q_S>Q_W) or by enstrophy (rigid rotation regions with Q_W>Q_S). In shear dominated regions, dissipation and enstrophy balance and -Q_S = Q_W <cit.>. For simplicity, we follow <cit.> and introduce 𝒦 = (-Q_W/Q_S)^1/2; when 𝒦 = 0 the flow is extension dominated, when 𝒦 = ∞ the flow is dominated by rigid rotation events, and when 𝒦 = 1 the rotation and the stretching are equal, as typical of vortex sheets and shear layers. In the unladen case, events with Q_W>-Q_S (𝒦→∞) are more frequent, meaning that for purely Newtonian turbulence the flow is mainly dominated by rigid rotations; see also the distributions of Q_S and Q_W in the bottom panels of figure <ref>. In the presence of the particles the scenario slightly changes. A first observation is that the probability of events with large Q_S (or 𝒦 = 0) increases, indicating that the perturbation field induced by these small particles is extensional dominated. A second observation is that the solid phase favours also events with -Q_S ≈ Q_W (𝒦≈ 1), that agrees with the presence of the shear layers induced by the presence of the particles. 



























§ PARTICLE DYNAMICS



This section is devoted to the dynamics of the particles, and we compare the results of the PR-DNS with those of the PP-DNS. Besides characterising the motion of the particles, indeed, the objective is to address the reliability of the one-way-coupled PP-DNS at the present parameters. 



 §.§ Lagrangian velocity increments






The Lagrangian statistics of the particles motion are of fundamental importance in the understanding of transport and mixing.
In order to investigate this, we study the Lagrangian velocity increments, defined here as δ_τ u_p,i = u_p,i (t + τ) - u_p,i(t), with u_p,i(t) being the instantaneous velocity of a particle along direction i at time t. The symmetries of the present problem make the statistics of δ_τ u_p,i independent on both t and i; for simplicity hereafter we drop the i index. Figures <ref> and <ref> describe the particle dynamics at different time scales, by plotting δ_τ u_p for different values of the time scale τ in the 0.2 τ_η≤τ≤ 30 τ_η range, where τ_η = (ν/ϵ)^1/2 is the Kolmogorov time scale. For small time scales, the velocity increment provides information about the particle acceleration, i.e. δ_τ u_p ∼ a_pτ. A first observation is that the distributions are symmetric, in agreement with the symmetries of the flow. The probability density function of δ_τ u_p continuously deforms from the Gaussian at large time scales (see τ≈ 30 τ_η) to the development of stretched exponential tails at dissipative time scales (see τ≈ 0.2 τ_η), which are the statistical signature of an intermittent Lagrangian dynamics; see <cit.>, <cit.> and <cit.> for small tracers and <cit.> for finite-size neutrally buoyant particles. The wide stretched exponential tails for the smallest τ show that the finite-size particles with D_p ≈η experience very high acceleration events, with a probability which is higher than Gaussian, similarly to what was found for small tracers by <cit.> and for finite-size particles by <cit.>. 

We start looking at the influence of ρ_p/ρ_f and Φ_V on the Lagrangian intermittency of Kolmogorov-size particles. 
The left panels of figure <ref> show that heavier particles (ρ_p/ρ_f=100) are less likely to experience intermediate values of the acceleration compared to lighter particles (ρ_p/ρ_f=5). In contrast, they are more likely to exhibit very low or very intense accelerations. On one side, the larger inertia of these particles opposes to large accelerations and favours small values of a_p. On the other side, heavy particles enhance the flow intermittency (see <ref>), promoting extreme events in the flow that are in turn responsible for rare but large particle accelerations. It is worth noticing that for Φ_V=10^-5 the latter effect is barely visible (see  figure <ref>), in agreement with the weak flow modulation shown in <ref>. For light particles ρ_p/ρ_f=5 figure <ref> shows that the distribution of δ_τ u_p obtained for Φ_V=10^-5 and 10^-3 overlap almost perfectly, in agreement with the low level of backreaction. 




The left panels of figure <ref> show that for ρ_p/ρ_f=5 the distributions obtained by means of PP-DNS and PR-DNS almost perfectly overlap for all time scales τ: for light particles the complete MRG equation properly predicts the Lagrangian intermittency of the particles dynamics. For heavier particles the match between the PP-DNS and the PR-DNS is rather good at large time scales, but differences are observed at small τ, where the PP-DNS does not predict the large tails for τ⪅ 2 τ_η. As discussed above, these extreme events are associated with the flow modulation which is not modelled in our PP-DNS.

The comparison between the PP-DNS and PR-DNS results is further detailed in figure <ref> where the evolution of the δ_τ u_p distribution with τ is quantified by means of the excess kurtosis 𝒦_δ_τ u_p(τ) = δ_τ u_p^4/δ_τ u_p^2 ^2 -3. At large scales 𝒦_δ_τ u_p≈ 0 in agreement with the Gaussian-like shape of the distribution, while it steeply increases at small scales. For ρ_p/ρ_f=5 the good agreement between the PP-DNS and the PR-DNS is again clear, with a small deviation for the Φ_V=10^-3 case, which is due to the non-zero flow modulation. For ρ_p/ρ_f=100 the agreement is good at large scales, while the three curves substantially deviate for small τ, in agreement with the larger tails of δ_τ u_p found in the PR-DNS. 




 §.§ The particle-velocity Structure function


We now consider the statistics of the particle-particle relative velocity δu_p = u_p(x_p,j(t),t) - u_p(x_p,i(t),t), where x_p,i(t) and x_p,j(t) denote the position of any two particles i and j at time t. The distribution of δu_p across all particle couples plays a key role in several theories regarding the tendency of particles to form clusters; see for example <cit.>, <cit.> and <cit.>.



Figure <ref> plots the second-order structure function of the particle velocity, i.e.


    S_2,p(r) = ( δu_p(r) ·r/r)^2 ,


where r is the separation vector between particle i and j, and r = |r|. The top panel is for ρ_p/ρ_f=5, while the bottom panel is for ρ_p/ρ_f=100. 


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

The bottom panel of figure <ref> deals with the ρ_p/ρ_f=100 cases. A first observation is that the results from the PR-DNS with Φ_V=10^-5 and Φ_V=10^-3 do not collapse; this is consistent with the larger flow modulation observed for the larger volume fraction, and agrees with the above discussed results for the single-particle statistics. Notably, for heavy particles S_2,p differs from the fluid structure function S_2 at the small scales. According to both the PP-DNS and the PR-DNS, S_2,p does not exhibit a r^2 scaling at the smallest scales, being substantially flat at small r. The relative motion between couples of heavy particles placed at a small distances r is substantially uncorrelated as well as decoupled from the small scale fluid motion due to their large inertia. Note that, the absence of the S_2,p∼ r^2 scaling indicates that at small scales the Eulerian particle velocity field cannot be described with a Taylor expansion. Notably, figure <ref> shows that S_2,p recovers the S_2-slope at larger r, exhibiting the classical Kolmogorov r^2/3 scaling in the inertial range of scales. This suggests that, despite the large inertia, the relative particle-particle velocity δu_p between two particles is driven by turbulent eddies having size comparable to r, provided that r is large enough.

When comparing the results of the PP-DNS with those of the PR-DNS, we note that the slope of S_2,p matches for small (r/η⪅ 5) and large (r ⪆ 20) scales. For intermediate scales 5 ⪅ r/η⪅ 20, instead, the PR-DNS predict a steeper slope for both Φ_V=10^-5 and Φ_V=10^-3. The finite size of the particles does not influence S_2,p for large and small scales where r/D_p = 𝒪(100) and r/D_p =𝒪(1), but it does for intermediate scales r/D_p = 𝒪(10). 













Figure <ref> sheds further light on the relative particle-particle velocity by plotting the distribution of δu_p ·r/r, i.e. the component of δu_p projected along the vector separating the two particles, for different values of r. When δu_p ·r/r>0, the two particles depart, whereas they get closer when δu_p ·r/r<0. We consider the case with ρ_p/ρ_f = 100 to provide further insights of the distribution of S_2,p, shown in the right panels of figure <ref>. A first observation is that, similarly to what found for larger particles by <cit.>, the distribution of δu_p ·r/r is left skewed, with a slightly positive mode and a long negative tail. The distribution becomes progressively more flat when increasing r, in agreement with a lower level of the correlation between the velocity of the two particles. When comparing the distributions for Φ_V = 10^-5 and Φ_V = 10^-3, figure <ref> shows that for all r the tails are shorter for the larger Φ_V, with the difference decreasing as r increases. This is consistent with the stronger flow modulation that globally leads to a weaker level of the flow fluctuations; see table <ref>. Also, in agreement with the evolution of S_2,p with r (see figure <ref>), for Φ_V=10^-5 the distribution with the PP-DNS collapses nicely with that obtained with the PR-DNS at small scales (see the top panels), with some substantial differences arising for intermediate scales when the finite-size effects are relevant.



§ THE COLLECTIVE MOTION OF THE PARTICLES



In this section we focus on the collective motion of the particles. First, we investigate whether Kolmogorov-size particles agglomerate and form clusters. Then, we relate the presence of the clusters with the tendency of particles to preferentially sample particular regions of the flow.



 §.§ Clustering


Over the years, several tools have been used to characterise the spatial arrangement of the particles in the flow <cit.>. We use the Voronoï tessellation, which has been extensively used by several authors <cit.>. The position of each particle is identified with its centre and the computational domain is divided in subdomains, such that each grid cell is associated with the closest particle. The Voronoï volume V_V of each particle is thus defined as the collective volume of grid cells that are closer to it than to other particles. The inverse of the Voronoï volumes provides a measure of the local concentration: particles placed in void regions possess a large Voronoï volume, while particles that are part of a cluster have a small Voronoï volume. Based on these observations, the intensity of clustering of a suspension can be measured by comparing the distribution of its Voronoï volumes to that of a control consisting of an equivalent, uniformly random, non-overlapping suspension of particles. More intense clustering leads to a variance of the distribution of V_V that is larger than that of the control. 



Figure <ref> presents the clustering intensity for the different values of Φ_V and ρ_p/ρ_f considered. A first observation is that the PP-DNS underestimate the level of clustering for all cases. Our computations show that the discrepancy between the PP-DNS and the PR-DNS increases with ρ_p/ρ_f and/or Φ_V (see <ref> for further details). 

We now move to the effect of the volume fraction Φ_V and of the particle density ρ_p/ρ_f. As expected, the level of clustering increases with Φ_V. When fixing Φ_V, instead, figure <ref> shows that heavier particles with ρ_p/ρ_f = 100 cluster more than lighter particles with ρ_p/ρ_f =5. For light particles, the low level of clustering σ/σ_rand≈ 1 agrees with the previous results of <cit.>, <cit.> and <cit.>, who considered larger particles 5 ≤ D_p/η≤ 123 over a wide range of Reynolds numbers 105 ≤ Re_λ≤ 430. Light particles have small inertia and are less likely to drift from the trajectories of the fluid elements. In contrast, the larger level of clustering observed when increasing the particle density from ρ_p/ρ_f = 5 to ρ_p/ρ_f=100 is not consistent with what found for larger particles. For suspensions of particles with size in the inertial range of scales, indeed, <cit.> found that the level of clustering exhibits a non-monotonous dependence on ρ_p/ρ_f, with the maximum occurring at intermediate densities (see figure 24 of their paper), and the minimum being at the largest density ratio they considered, i.e. ρ_p/ρ_f = 100. However, they found that in the D_p-ρ_p space of parameters the maximum level of clustering moves towards larger ρ_p as D_p decreases, suggesting that the tendency of particles to cluster is driven by the Stokes number of the particles, rather than by their density. Accordingly, their data show that the level of clustering is maximum when St = 𝒪(1-10). This agrees with the early works of <cit.>, <cit.> and <cit.>, and it is consistent with the our results. 



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

A different type of information regarding the spatial arrangement of the particles can be provided by means of the radial distribution function g(r), also referred to as pair correlation function (see figure <ref>). It describes how the particle density varies as a function of the distance away from a reference particle. In other words, it is a measure of the probability of finding particles at a distance r relative to that of a homogeneous distribution. Following <cit.> and <cit.>, the radial distribution function is defined as


    g(r) =  N_s(r)/Δ V(r) / N_pa/V ,


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



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



 §.§ Preferential sampling



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

In the limit where the point-particle approximation holds, <cit.> has shown that particles with large density tend to collect in regions of high strain rate and low vorticity. In the presence of a vortex, indeed, heavy particles cannot  follow the flow streamlines because of their large inertia, and tend to drift from the vortex core. Similarly, in the case of a pure straining flow, heavy particles drift towards the stagnation point at the centre. 

We quantify the tendency of particles to sample regions of high strain by using the second invariant of the deformation tensor (see <ref>), i.e. Q = - S_ijS_ij/2 + ω^2/4; recall that regions where Q is large and positive are regions of high vorticity (Q ∼ω^2/4), while regions where Q is large and negative are regions of high strain (Q ∼ - S_ijS_ij). We investigate the particle preferential sampling by computing the probability density function of Q at the particles position <cit.>. For the sake of brevity, in this section we limit the investigation to the largest volume fraction Φ_V = 10^-3. For PP-DNS the value of Q at each particle position is obtained after linear interpolation. For PR-DNS, instead, the value of Q seen by each particle is estimated as the average value within a spherical shell centred with the particle and having a radius R_sh>R_p, where R_p=D_p/2 is the radius of the particles. It is important to note, however, that due to the particles' backreaction, in the case of PR-DNS, the value of Q seen by each particle is actually the result of three different effects: (i) the larger scale flow properties of the region that the particle is sampling, (ii) the smaller scale influence of the particle on the surrounding flow, and (iii) the effect of nearby particles on the flow. Also, a suitable choice of R_sh should be made: when R_sh is too small, only the influence of the particle on the surrounding flow is considered <cit.>, while when R_sh is too large, spurious contributions that do not affect the particle location are instead captured. In order to obtain a more complete picture, we have tested different values of R_sh.



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

Let us now focus on the differences between the PR-DNS and PP-DNS results. For particles with ρ_p/ρ_f=5, the Q_ℓ distribution obtained with PP-DNS overlaps almost perfectly with that obtained with PR-DNS; the point-particle approximation predicts fairly well the tendency of light particles to sample the Q<0 regions of the flow. Note that this is consistent with the good agreement found in figure <ref> when discussing the distribution of the Voronoï volumes. For ρ_p/ρ_f = 100, instead, the Q_ℓ distribution obtained with the PP-DNS significantly deviates from that obtained with the PR-DNS. For all R_sh, the Q_ℓ distribution obtained with PR-DNS shows a shorter right tail and predicts a higher probability of negative Q: finite-size heavy particles are less/more likely to see positive/negative values of Q. Based on this, one may conclude that, at the present parameters, the PP-DNS underestimates the tendency of the particles to preferentially sample regions of high-strain, and this may explain the larger level of clustering predicted by the PR-DNS for the Φ_V = 10^-3 and ρ_p/ρ_f = 100 case (see figure <ref>). 

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




To provide additional insights regarding the relation between the presence of clusters and the particle preferential sampling, figure <ref> shows the joint probability density function of Q_ℓ and V_V for Φ_V = 10^-3. Based on the above discussion, here we set R_sh/R_p = 3 for the computation of Q_ℓ, since it is large enough to account for the particle preferential location and small enough to avoid spurious contributions. We recall that according to <cit.>, particles with V_V ≤ V_th,l are part of a cluster, while particles with V_V ≥ V_th,r are in void regions of the flow. For all cases, the most negative values of Q_ℓ well correlate with small and intermediate Voronoï volumes, with V_V ⪅ V_th,r. Particles that are in void regions and are not part of a cluster are less likely to see large negative values of Q. This agrees with the above observation that the tendency of Kolmogorov-size particles to form clusters is governed by the centrifuge mechanism. We now focus on ρ_p/ρ_f=100 (see the bottom panels). The joint distribution shows that the larger probability of negative Q_ℓ predicted by the PR-DNS is concentrated at the smallest Voronoï volumes with V_V ⪅ V_th,l. Again, this shows that the higher level of clustering detected in this case with the PR-DNS well correlates with the larger tendency of finite-size particles to sample regions of the flow with intense strain.




































































§ CONCLUSION



We have investigated by direct numerical simulations the fluid-solid interaction of suspensions of Kolmogorov-size spherical particles moving in homogeneous isotropic turbulence. The work is based on both interface-resolved (PR-DNS) and one-way-coupled point-particle (PP-DNS) direct numerical simulations. In PR-DNS the presence of the particles is dealt with the immersed boundary method introduced by <cit.>. In PP-DNS the motion of the particles is described by solving the complete Maxey-Riley-Gatignol equation <cit.>, including the time-history Basset term and the Faxén correction. The objective of the study is twofold. On one side, we aim to shed light on how Kolmogorov-size particles influence the organisation of the velocity fluctuations. Few work have indeed considered particles with D_p/η≈ 1 due to the intrinsic complexity of the problem: in experiments it requires the access to sub-Kolmogorov measurements, and in simulations it requires an extremely fine grid with a resulting prohibitive computational cost. On the other side, we aim to assess the limits of the one-way-coupled PP-DNS that, despite the large number of works present in literature, have not been completely addressed yet. For this reason, we consider a portion of the parameter space that is on the edge of the range of validity of the one-way coupled point-particle models <cit.>. The micro-scale Reynolds number is Re_λ≈ 140, being large enough to ensure a proper separation of scales. The volume fraction of the suspension has been set to the small values of Φ_V=10^-5 and Φ_V = 10^-3 to guarantee that the backreaction of the solid phase on the carrier flow is low. Two solid-to-fluid density ratios are considered, i.e. ρ_p/ρ_f = 5 and ρ_p/ρ_f = 100, to investigate the role of inertia. 

The PR-DNS shows that at the present parameters the modulation of the flow is rather low and mainly involves the smallest scales. The modulated energy spectrum ℰ(κ) shows a multi-scaling behaviour: the classical -5/3 scaling in the inertial range of scales is indeed followed by a steeper -4 scaling, that resembles what has been observed by several authors in the context of bubbly flows <cit.>. Accordingly, the scale-by-scale energy budget shows two different regimes: in the inertial range of scales the fluid-structure interaction term is negligible and the non linear term equals the dissipation rate, i.e. Π(κ) ∼ϵ; at these scales energy is transferred from larger to smaller scales by means of the classical energy cascade described by Richardson and Kolmogorov. At small scales, where ℰ(κ) ∼κ^-4, the non linear term is negligible and the fluid-structure interaction term is balanced by the viscous dissipation, i.e. Π_fs(κ) ∼ D_v(κ); at these scales the energy injected into the flow by the particles is directly dissipated by viscosity. 
The small-scale topology of the flow has been investigated by inspecting the influence of the particles on the invariants of the velocity gradient tensor A_ij = ∂ u_i/∂ x_j <cit.>. The effect of the solid phase on the eigenvalues of the strain-rate tensor shows that the presence of the particles favours axisymmetric compression rather than axisymmetric extension. The joint p.d.f of the second and third invariants of A_ij reveals that, particles mainly enhance events with axial strain and vortex compression. Accordingly, the inspection of the joint p.d.f. of the second invariants of the symmetric and antisymmetric parts of A_ij indicates that the presence of the particles favours dissipation events dominated by extensional events rather than rotational ones.

The limits of the one-way-coupled point-particle models have been addressed by looking at the dynamics of the particles and at their collective motion. We find that the PP-DNS predicts fairly well the Lagrangian and Eulerian statistics of the particles velocity field for the low-density case. For heavy particles, however, some discrepancies are observed, particularly for the larger volume fraction. These differences are due to a combination of the finite-size of the particles and of the flow modulation, that are not accounted for in the PP-DNS. By using the Voronoï tessellation method and the radial distribution function, we find that the PP-DNS underpredicts the level of clustering; the discrepancy with the PR-DNS results increases with the volume fraction and the particle density. In the attempt to have a clearer picture, we have investigated the tendency of the particles to preferentially sample particular regions of the flow. By plotting the distribution of the second invariant Q of the fluid velocity gradient tensor at the particle position, we find that, according to both PP-DNS and PR-DNS, the particles preferentially sample regions of high strain rate. This suggests that the presence of the clusters is driven by the centrifuge mechanism introduced by <cit.>. Accordingly with the larger level of clustering, we find that PR-DNS shows a larger tendency of the particles to sample these regions of the flow compared to PP-DNS. Note, however, that some care is needed when dealing with these results. In PR-DNS, indeed, the value of Q seen by each particle is the result of three different contributions that cannot be easily isolated, i.e. (i) the larger scale flow properties of the region that the particle is sampling, (ii) the smaller scale influence of the particle on the surrounding flow, and (iii) the effect of nearby particles on the flow.




By characterising the fluid-solid interaction of Kolmogorov-size particles in homogeneous isotropic turbulence, the present study aims to serve as a stepping stone for further investigations. A natural extension of this work is to use the present PR-DNS database to assess the limits of two-way-coupled PP-DNS, that account for the backreaction of the solid phase to the carrier flow; see for instance the models introduced by <cit.> and <cit.>. In addition, the present results may be used as a ground truth for studies in the spirit of <cit.>, that investigate the relevance of each term at the right-hand-sides of the MRG equation in predicting the different statistics of the particles. This knowledge will help guide the choice of suitable models for engineering applications. Eventually, it would be of interest to investigate whether Kolmogorov-size particles modulate the energy spectrum also in the inertial range of scales at larger volume fractions, influencing thus the -5/3 scaling range and the classical energy cascade as observed for larger particles <cit.>. Despite the computational challenge that would such a study would provide, the field would benefit from the investigation. Overall, the present results can be exploited for the development of improved point-particle models for the one- and two-way-coupling regimes. 











§ ACKNOWLEDGMENTS
 
The authors acknowledge the computer time provided by the Scientific Computing and Data Analysis section of the Core Facilities at Oist and the computational resources on Fugaku provided by HPCI (project ID: hp230536). A.C. acknowledges Marildo Kola for discussions and suggestions.
  


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


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



§ THE BASSET TIME HISTORY FORCE



In the PP-DNS we resort on the second-order and memory-efficient algorithm developed by <cit.> to deal with the Basset time history force. 

The Basset force is split into two parts, denoted as window and tail. In particular, at time t̃ the first part consists of a numerical integration over the t̃-t_w≤ t ≤t̃ interval, considering thus N_w=t_w/Δ t previous steps. The second integral, instead, considers the -∞≤ t ≤t̃ - t_w interval and is approximated using recursive exponential functions.
The kernel K_B(t) in equation <ref> is thus replaced with an approximated kernel K(t) such that 

    K(t) =  
                      K_B(t)    if  t < t_win
    
                      K_tail(t)    if  t ≥ t_twin,

with

    lim_t → +∞ K_tail(t) = 0.

The Basset force, therefore, reads 

    𝐅_B(t) = c_B ∫_t-t_win^t K_B(t-τ) g(τ) dτ_𝐅_B-win(t) +
    c_B ∫_-∞^t-t_win K_tail(t-τ) g(τ) dτ_𝐅_B-tail(t),

where c_B = 3/2 D_p^2 ρ_f √(πν) and g(t) = d(𝐮-𝐮_p + (1/6) (D_p/2)^2 ∇^2 𝐮) / dt. The window term is integrated in time using the diffusive Basset kernel. The integration exploits a modified trapezoidal method, which allows the kernel's singularity to be taken into account. Thus, following the work of <cit.>, the window contribution reads:

    𝐅_B-win   = 4/3 c_B √(Δ t)𝐠_0 + 4/3 c_B √(Δ t)∑_n = 1^N_w-1[  (n-1)√(n-1) - 2n√(n) + (n+1)√(n+1)] 𝐠_n 
       + c_B √(Δ t)[ 4/3(N_w - 1)√(N_w-1) + (2-4/3N_w)√(N_w)]𝐠_N_w,
    

where g_n = g(t - n Δ t ) with  n = 0,1,..., N_w. Here Δ t = t_win/N_w and N_w is the number interval in which the window is discretised.

As stated above, the tail term is integrated in a recursive manner and, exploiting exponential kernels, it reads:


    𝐅_B-tail(t) = ∑_i=1^m a_i 𝐅_i(t) = ∑_i=1^m  a_i ( 𝐅_i-di(t) + 𝐅_i-re(t) ),


where 𝐅_i-di is computed directly as

    𝐅_i-di(t) =  2 c_B √(e t_i)exp(-t_win/2t_i) {𝐠_N [ 1 - ϕ( Δ t/2t_i) ] +  𝐠_N+1[ϕ(- Δ t/2t_i) - 1 ]  },


and 𝐅_i-re is computed recursively as

    𝐅_i-re(t) = exp( -Δ t/2 t_i) 𝐅_i(t-Δ t ).


Here, ϕ(z) = (exp(z)-1)/z, and for a given value of m the coefficients {a_i,t_i }_i=1^m are chosen to minimize the error. For a detailed explanation, the reader is referred to <cit.> and <cit.>. We choose m=10 and set the values of the a_i and t_i parameters to the ones proposed in the work of <cit.>. As suggested by <cit.> the point-particle equation is written in a semi-implicit manner to guarantee numerical stability when integrating in time.



jfm

101
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

[Alised et al.(2002)Alised, Cartellier, Hainaux &
  Lasheras]aliseda-etal-2002
Alised, A., Cartellier, A., Hainaux, F. & Lasheras,
  J.C. 2002  Effect of preferential concentration on the settling
  velocity of heavy particles in homogeneous isotropic turbulence.  J.
  Fluid Mech.  468,  77–105.

[Alméras et al.(2017)Alméras, Mathai, Lohse &
  Sun]almeras-etal-2017
Alméras, E., Mathai, V., Lohse, D. & Sun, C.
  2017  Experimental investigation of the turbulence induced by a
  bubble swarm rising within incident turbulence.  J. Fluid Mech.
  825,  1091–1112.

[Auton et al.(1988)Auton, Hunt &
  Prud’Homme]auton-hunt-1988
Auton, T. R., Hunt, J. C. R. & Prud’Homme, M. 1988
   The force exerted on a body in inviscid unsteady non-uniform rotational
  flow.  J. Fluid Mech.  197,  241–257.

[Balachandar(2009)]balachandar-2009
Balachandar, S. 2009  A scaling analysis for
  point–particle approaches to turbulent multiphase flows.  Int. J.
  Multiph. Flow  35,  801–810.

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

[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.

[Betchov(1956)]betchov-1956
Betchov, R. 1956  An inequality concerning the production
  of vorticity in isotropic turbulence.  J. Fluid Mech.  1 (5),
  497–504.

[Boivin et al.(1998)Boivin, Simonin &
  Squires]boivin-simonin-squires-1998
Boivin, M., Simonin, O. & Squires, K.D. 1998
  Direct numerical simulation of turbulence modulation by particles in
  isotropic turbulence.  J. Fluid Mech.  375,  235–263.

[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 J. Phys.  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.  Phys. Rev. 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.

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

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

[Cannon et al.(2024)Cannon, Olivieri &
  Rosti]cannon-olivieri-rosti-2024
Cannon, I., Olivieri, S. & Rosti, M.E. 2024
  Spheres and fibers in turbulent flows at various reynolds numbers.
  Phys. Rev. Fluids  9,  064301.

[Casas et al.(2018)Casas, Ferrer &
  Oñate]casas-ferrer-onate-2017
Casas, G., Ferrer, A. & Oñate, E. 2018
  Approximating the basset force by optimizing the method of van hinsberg
  et al.  J. Comput. Phys.  352,  142–171.

[Chevillard et al.(2003)Chevillard, Roux, Levêque, Mordant,
  Pinton & Arneodo]chevillard-etal-2003
Chevillard, L., Roux, S. G., Levêque, E., Mordant,
  N., Pinton, J.-F. & Arneodo, A. 2003  Lagrangian
  Velocity Statistics in Turbulent Flows: Effects of
  Dissipation.  Phys. Rev. Lett.  91 (21),  214502.

[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.

[Chiarini & Rosti(2024)]chiarini-rosti-2024
Chiarini, A. & Rosti, M.E. 2024  Finite-size inertial
  spherical particles in turbulence.  J. Fluid Mech.  988,
  A17.

[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.

[Costa et al.(2020)Costa, Brandt &
  Picano]costa-brandt-picano-2020
Costa, P., Brandt, L. & Picano, F. 2020
  Interface-resolved simulations of small inertial particles in turbulent
  channel flow.  J. Fluid Mech.  883,  A54.

[Cundall & Strack(1979)]cundall-strack-1979
Cundall, P. A. & Strack, O. D.L. 1979  A discrete
  numerical model for granular assemblies.  Geotechnique  29 (1),
  47–65.

[Davidson(2004)]davidson-2004
Davidson, P.A. 2004 Turbulence: An Introduction for
  Scientists and Engineers.  Oxford University Press.

[Davidson & Pearson(2005)]davidson-pearson-2005
Davidson, P. A. & Pearson, B. R. 2005  Identifying
  Turbulent Energy Distributions in Real, Rather than Fourier,
  Space.  Phys. Rev. Lett.  95 (21),  214501, publisher:
  American Physical Society.

[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.

[Dorgan & Loth(2007)]dorgan-loth-2007
Dorgan, A.J. & Loth, E. 2007  Efficient calculation
  of the history force at finite reynolds numbers.  Int. J. Multiphase
  Flow  33,  833–848.

[Druzhinin(2001)]druzhinin-2001
Druzhinin, O. A. 2001  The influence of particle inertia on
  the two-way coupling and modification of isotropic turbulence by
  microparticles.  Phys. Fluids  13 (12),  3738–3755.

[Dung et al.(2023)Dung, Waasdorp, Sun, Lohse &
  Huisman]dung-etal-2022
Dung, O.-Y., Waasdorp, P., Sun, C., Lohse, D. &
  Huisman, S.G. 2023  The emergence of bubble-induced scaling in
  thermal spectra in turbulence.  J. Fluid Mech.  958,  A5.

[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.

[Faxén(1922)]faxen-1922
Faxén, H. 1922  Der widerstand gegen die bewegung einer
  starren kugel in einer zähen flüssigkeit, die zwischen zwei parallelen
  ebenen wänden eingeschlossen ist.  Annalen der Physik
  373 (10),  89–119.

[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.

[Fessler et al.(1994)Fessler, Kulick &
  Eaton]fessler-etal-1994
Fessler, J. R., Kulick, J.D. & Eaton, J.K. 1994
  Preferential concentration of heavy particles in a turbulent channel
  flow.  Phys. Fluids  6 (11),  3742–3749.

[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.

[Gatignol(1983)]gatignol-1983
Gatignol, R. 1983  The faxen formulae for a rigid particle
  in an unsteady non-uniform stokes flow.  J. méc. théor. appl.
  2 (2),  143–160.

[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.

[Gualtieri et al.(2015)Gualtieri, Picano, Sardina &
  Casciola]gualtieri-etal-2015
Gualtieri, P., Picano, F., Sardina, G. & Casciola,
  C. M. 2015  Exact regularized point particle method for
  multiphase flows in the two-way coupling regime.  J. Fluid Mech.
  773,  520–561.

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

[van Hinsberg et al.(2011)van Hinsberg, ten Thije Boonkkamp &
  Clercx]hinsberg-boonkkamp-clercx-2011
van Hinsberg, M.A.T., ten Thije Boonkkamp, J.H.M. & Clercx,
  H.J.H. 2011  An efficient, second order method for the
  approximation of the basset history force.  J. Comput. Phys.
  230,  1465–1478.

[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.

[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.

[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-frolich-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,  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(1941)]kolmogorov-1941
Kolmogorov, A.N. 1941  The Local Structure of
  Turbulence in an Incompressible Viscous Fluid for Very Large
  Reynolds Numbers.  Dokl. Akad. Nauk. SSSR  30,  301–305.

[La Porta et al.(2001)La Porta, Voth, Crawford, Alexander &
  Bodenschatz]laporta-etal-2001
La Porta, A., Voth, Greg A., Crawford, Alice M.,
  Alexander, Jim & Bodenschatz, Eberhard 2001  Fluid
  particle accelerations in fully developed turbulence.  Nature
  409 (6823),  1017–1019.

[Lance & Bataille(1991)]lance-betaille-1991
Lance, M. & Bataille, J. 1991  Turbulence in the
  liquid phase of a uniform bubbly air–water flow.  J. Fluid Mech.
  222,  95–118.

[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.

[Lund & Rogers(1994)]lund-rogers-1994
Lund, Thomas S. & Rogers, Michael M. 1994  An
  improved measure of strain state probability in turbulent flows.  Phys.
  Fluids  6 (5),  1838–1847.

[Martinez Mercado et al.(2010)Martinez Mercado, Chehata Gomez,
  Van Gils, Sun & Lohse]mercado-etal-2010
Martinez Mercado, J., Chehata Gomez, D., Van Gils, D.,
  Sun, C. & Lohse, D. 2010  On bubble clustering and
  energy spectra in pseudo-turbulence.  J. Fluid Mech.  650,
  287–306.

[Matsuda et al.(2024)Matsuda, Yoshimatsu &
  Schneider]matsuda-etal-2024
Matsuda, K., Yoshimatsu, K. & Schneider, K. 2024
  Heavy particle clustering in inertial subrange of high–reynolds number
  turbulence.  Phys. Rev. Lett.  132,  234001.

[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.

[McLaughlin(1991)]mclaughlin-1991
McLaughlin, J.B. 1991  Inertial migration of a small sphere
  in linear shear flows.  J. Fluid Mech.  224,  261–274.

[Mehrabadi et al.(2018)Mehrabadi, Horwitz, Subramaniam &
  Mani]mehrabadi-etal-2018
Mehrabadi, M., Horwitz, J. A. K., Subramaniam, S. &
  Mani, A. 2018  A direct comparison of particle-resolved and
  point-particle methods in decaying turbulence.  J. Fluid Mech.
  850,  336–369.

[Mei(1992)]mei-1992
Mei, R. 1992  An approximate expression for the shear lift
  force on a spherical particle at finite reynolds number.  Int. J.
  Multiph. Flow  18,  145–147.

[Mei & Adrian(1992)]mei-adrian-1992
Mei, R. & Adrian, R. J. 1992  Flow past a sphere with
  an oscillation in the free-stream velocity and unsteady drag at finite
  reynolds number.  J. Fluid Mech.  237,  323–341.

[Meneveau(2011)]meneveau-2011
Meneveau, C. 2011  Lagrangian Dynamics and Models
  of the Velocity Gradient Tensor in Turbulent Flows.  Annu. Rev.
  Fluid Mech.  43 (1),  219–245.

[Michaelides(1992)]michaelides-1992
Michaelides, E.E. 1992  A novel way of computing the basset
  term in unsteady multiphase flow computations.  Phys. Fluids A
  4 (7),  1579–1582.

[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.

[Mordant et al.(2001)Mordant, Metz, Michel &
  Pinton]mordant-etal-2001
Mordant, N., Metz, P., Michel, O. & Pinton, J.-F.
  2001  Measurement of Lagrangian Velocity in Fully Developed
  Turbulence.  Phys. Rev. Lett.  87 (21),  214501.

[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.  Int. J. Multiph. Flow  85,
  369–379.

[Nomura & Post(1998)]nomura-post-1998
Nomura, K.K. & Post, G.K. 1998  The structure and
  dynamics of vorticity and rate of strain in incompressible homogeneous
  turbulence.  J. Fluid Mech.  377,  65–97.

[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.(2022a)Olivieri, Cannon &
  Rosti]olivieri-cannon-rosti-2022
Olivieri, S., Cannon, I. & Rosti, M. E.
  2022a  The effect of particle anisotropy on the
  modulation of turbulent flows.  J. Fluid Mech.  950,  R2.

[Olivieri et al.(2022b)Olivieri, Mazzino
  & Rosti]olivieri-mazzino-rosti-2022
Olivieri, S., Mazzino, A. & Rosti, M. E.
  2022b  On the fully coupled dynamics of
  flexible fibres dispersed in modulated turbulence.  J. Fluid Mech.
  946,  A34.

[Olivieri et al.(2014)Olivieri, Picano, Sardina, Iudicone &
  Brandt]olivieri-etal-2014
Olivieri, S., Picano, F., Sardina, G., Iudicone, D. &
  Brandt, L. 2014  The effect of the basset history force on
  particle clustering in homogeneous and isotropic turbulence.  Phys.
  Fluids  26 (4),  041704.

[Pandey et al.(2022)Pandey, Mitra &
  Perlekar]pandey-mitra-perlekar-2022
Pandey, V., Mitra, D. & Perlekar, P. 2022
  Turbulence modulation in buoyancy-driven bubbly flows.  J. Fluid
  Mech.  932,  A19.

[Pandey et al.(2023)Pandey, Mitra &
  Perlekar]pandey-mitra-perlekar-2023
Pandey, V., Mitra, D. & Perlekar, P. 2023
  Kolmogorov turbulence coexists with pseudo-turbulence in buoyancy-driven
  bubbly flows.  Phys. Rev. Lett.  131,  114002.

[Pandey et al.(2020)Pandey, Ramadugu &
  Perlekar]pandey-ramadugu-perlekar-2020
Pandey, V., Ramadugu, R. & Perlekar, P. 2020
  Liquid velocity fluctuations and energy spectra in three-dimensional
  buoyancy-driven bubbly flows.  J. Fluid Mech.  884,  R6.

[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.

[Prakash et al.(2016)Prakash, Martínez Mercado, van
  Wijngaarden, Mancilla, Tagawa, Lohse & Sun]prakash-etal-2016
Prakash, V. N., Martínez Mercado, J., van Wijngaarden, L.,
  Mancilla, E., Tagawa, Y., Lohse, D. & Sun, C. 2016
   Energy spectra in turbulent bubbly flows.  J. Fluid Mech.
  791,  174–190.

[Prasath et al.(2019)Prasath, Vasan &
  Govindarajan]prasath-etal-2019
Prasath, S. Ganga, Vasan, Vishal & Govindarajan, Rama
  2019  Accurate solution method for the maxey–riley equation, and
  the effects of basset history.  J. Fluid Mech.  868,
  428–460.

[Qureshi et al.(2007)Qureshi, Bourgoin, Baudet, Cartellier &
  Gagne]qureshi-etal-2007
Qureshi, N.M., Bourgoin, M., Baudet, C., Cartellier,
  A. & Gagne, Y. 2007  Turbulent Transport of Material
  Particles: An Experimental Study of Finite Size Effects.
  Phys. Rev. Lett.  99 (18),  184502.

[Riboux et al.(2010)Riboux, Risso &
  Legendre]riboux-risso-legendre-2010
Riboux, G., Risso, F. & Legendre, D. 2010
  Experimental characterization of the agitation generated by bubbles
  rising at high reynolds number.  J. Fluid Mech.  643,
  509–539.

[Risso(2018)]risso-2018
Risso, F. 2018  Agitation, mixing, and transfers induced by
  bubbles.  Ann. Rev. Fluid Mech.  50,  25–48.

[Saffman(1965)]saffman-1965
Saffman, P. G. 1965  The lift on a small sphere in a slow
  shear flow.  J. Fluid Mech.  22 (2),  385–400.

[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.

[Schneiders et al.(2017)Schneiders, Meinke &
  Schröder]schneiders-etal-2017
Schneiders, L., Meinke, M. & Schröder, W. 2017
  Direct particle–fluid simulation of
  Kolmogorov-length-scale size particles in decaying isotropic turbulence.
   J. Fluid Mech.  819,  188–227.

[Schumacher et al.(2007)Schumacher, Sreenivasan &
  Yakhot]schumacher-etal-2007
Schumacher, J., Sreenivasan, K.R. & Yakhot, V. 2007
  Asymptotic exponents from low-Reynolds-number flows.  New J.
  Phys.  9 (4),  89.

[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.

[Soria et al.(1994)Soria, Sondergaard, Cantwell, Chong &
  Perry]soria-etal-1994
Soria, J., Sondergaard, R., Cantwell, B. J., Chong,
  M. S. & Perry, A. E. 1994  A study of the fine‐scale
  motions of incompressible time‐developing mixing layers.  Phys.
  Fluids  6 (2),  871–884.

[Squires & Eaton(1990)]squires-eaton-1990
Squires, K.D. & Eaton, J.K. 1990  Particle response
  and turbulence modification in isotropic turbulence.  Phys. Fluids A
  2 (7),  1191–1203.

[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.

[Truesdell(1954)]truesdell-1954
Truesdell, C. 1954 The Kinematics of Vorticity.
  Indiana University Press.

[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.

[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 & 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.

[Vreman(2016)]vreman-2016
Vreman, A. W. 2016  Particle-resolved direct numerical
  simulation of homogeneous isotropic turbulence modified by small fixed
  spheres.  J. Fluid Mech.  796,  40–85.

[Wang et al.(2023)Wang, Jiang & Sun]wang-etal-2023
Wang, C., Jiang, L. & Sun, C. 2023  Numerical
  study on turbulence modulation of finite-size particles in plane-couette
  flow.  J. Fluid Mech.  970,  A7.

[Wang & Maxey(1993)]wang-maxey-1993
Wang, L.-P. & Maxey, M.R. 1993  Settling velocity and
  concentration distribution of heavy particles in homogeneous isotropic
  turbulence.  J. Fluid Mech.  256,  27–68.

[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.



