1. Introduction
Particle-turbulence interaction is a topic that has received continued attention from the scientific community for several decades (Maxey & Riley Reference Maxey and Riley1983; Maxey Reference Maxey1987; Squires & Eaton Reference Squires and Eaton1990, Reference Squires and Eaton1991; Elghobashi Reference Elghobashi1994; Crowe, Troutt & Chung Reference Crowe, Troutt and Chung1996; Poelma & Ooms Reference Poelma and Ooms2006; Balachandar & Eaton Reference Balachandar and Eaton2010). The coupling between continuous and dispersed phase, the number of governing parameters, and the wide range of scales, make for a formidable challenge for experimentalists and modellers alike. Driven by the practical relevance of the problem, particle-laden turbulence has attracted renewed attention in the last decade, also thanks to the tremendous progress in our ability to carry out novel measurements and simulations in previously inaccessible conditions (Monchaux, Bourgoin & Cartellier Reference Monchaux, Bourgoin and Cartellier2012; Tenneti & Subramaniam Reference Tenneti and Subramaniam2014; Gustavsson & Mehlig Reference Gustavsson and Mehlig2016; Maxey Reference Maxey2017; Mathai, Lohse & Sun Reference Mathai, Lohse and Sun2020; Brandt & Coletti Reference Brandt and Coletti2022).
One of the most elusive aspects of the problem is the back-reaction of inertial particles on the fluid, referred to as two-way coupling. Turbulence modification by particles is believed to be significant already at modest loadings, i.e. above volume fractions 
$\varPhi _{v}={O}(10^{-6})$ for gas–solid systems (Elghobashi Reference Elghobashi1994). However, the evidence on whether the turbulence is augmented or attenuated, and to which extent, is controversial. A dated but still utilized criterion was proposed by Gore & Crowe (Reference Gore and Crowe1991), who compiled previous data and concluded that turbulence was augmented/attenuated by particles larger/smaller than one-tenth of the integral scale. This implies that sub-Kolmogorov particles (whose diameter 
$d_{p}$ is smaller than the Kolmogorov scale 
$\eta$) will always attenuate turbulence even when it exhibits marginal scale separation, at odds with several experimental observations (e.g. Yang & Shy Reference Yang and Shy2005). Some early experiments focused on wall-bounded flows, some reporting attenuation (e.g. Kulick, Fessler & Eaton Reference Kulick, Fessler and Eaton1994; Paris Reference Paris2001) and others augmentation (e.g. Sato & Hishida Reference Sato and Hishida1996). Tanaka & Eaton (Reference Tanaka and Eaton2008) compiled experimental results for turbulence modification in internal flows. They proposed a criterion based on a non-dimensional parameter combining the particle response time, density and size relative to the flow scales, as well as the flow Reynolds number, signalling the complexity of the coupling. Numerical simulations have struggled to reproduce the measurements, even when laboratory conditions were matched and the resolution was appropriate for direct numerical simulation (DNS) of the single-phase turbulence; see e.g. Vreman (Reference Vreman2015) and Wang et al. (Reference Wang, Fong, Coletti, Capecelatro and Richter2019). This is partly attributable to the limitations of the classic point-particle method, in which particles are treated as material points applying pointwise forcing on the fluid computational grid; see Eaton (Reference Eaton2009), Maxey (Reference Maxey2017) and Brandt & Coletti (Reference Brandt and Coletti2022).
Studying turbulence modification in wall-bounded flows, while highly relevant to several applications, requires disentangling the simultaneous effects on the hierarchy of scales in the boundary layer (Wang & Richter Reference Wang and Richter2019). Moreover, the particle–fluid slip velocity, which plays a crucial role in the momentum coupling, varies significantly with wall distance (Kiger & Pan Reference Kiger and Pan2002; Fong, Amili & Coletti Reference Fong, Amili and Coletti2019; Berk & Coletti Reference Berk and Coletti2020; Baker & Coletti Reference Baker and Coletti2021), complicating the generalization of the results. For example, Righetti & Romano (Reference Righetti and Romano2004) found that turbulence fluctuations were enhanced in the viscous sub-layer but damped further away from the wall. Additionally, at the centre-plane of a channel flow, the turbulence is anisotropic at moderate Reynolds numbers (Andersson, Zhao & Variano Reference Andersson, Zhao and Variano2015). Thus homogeneous isotropic turbulence may be a more suitable terrain to attack the fundamental aspects of the inter-phase coupling. This is the focus of the present investigation.
Creating homogeneous isotropic turbulence in the laboratory is notoriously challenging (Bellani & Variano Reference Bellani and Variano2014). The classic approach is grid-generated turbulence, which was used by Schreck & Kleis (Reference Schreck and Kleis1993), Geiss et al. (Reference Geiss, Dreizler, Stojanovic, Chrigui, Sadiki and Janicka2004) and Poelma, Westerweel & Ooms (Reference Poelma, Westerweel and Ooms2007) to investigate particle–turbulence interaction. The dispersed phase was found to induce large-scale anisotropy, with larger axial velocity fluctuations. Highly inertial particles caused mild turbulence attenuation over the entire spectrum, while particles with a Stokes number 
$St_{\eta } = \tau _{p}/\tau _{\eta } < 1$ energized the small-scale fluctuations and modulated the large-scale ones (where 
$\tau _{p}$ is the particle response time, and 
$\tau _{\eta }$ is the Kolmogorov time scale). The experiments of Poelma et al. (Reference Poelma, Westerweel and Ooms2007) were the first to describe clearly such ‘pivoting’ of the energy spectrum, which had been reported in several previous simulations (as summarized by Poelma & Ooms Reference Poelma and Ooms2006). A similar behaviour had been observed previously in bubble-laden turbulence (Mazzitelli, Lohse & Toschi Reference Mazzitelli, Lohse and Toschi2003a,Reference Mazzitelli, Lohse and Toschib; Rensen, Luther & Lohse Reference Rensen, Luther and Lohse2005).
Grid turbulence, however, presents shortcomings associated with its spatial decay, particularly the weak level of turbulent agitation; in Poelma et al. (Reference Poelma, Westerweel and Ooms2007), the Reynolds numbers based on the Taylor microscale was 
$Re_{\lambda }<30$. Zero-mean flow turbulence chambers (usually stirred by jets) overcome this limitation and can approximate homogeneous isotropic turbulence at Reynolds numbers with substantial scale separation (Bellani & Variano Reference Bellani and Variano2014; Carter et al. Reference Carter, Petersen, Amili and Coletti2016). They were used by Hwang & Eaton (Reference Hwang and Eaton2006a,Reference Hwang and Eatonb) and Tanaka & Eaton (Reference Tanaka and Eaton2010) to show that inertial particles with 
$St_{\eta } = {O}(10^{2})$ attenuate turbulence. On the other hand, Yang & Shy (Reference Yang and Shy2005) reported turbulence augmentation for 
$St_{\eta } = {O}(1)$. In this regime, inertial particles are known to cluster over multi-scale sets (Squires & Eaton Reference Squires and Eaton1991; Monchaux et al. Reference Monchaux, Bourgoin and Cartellier2012; Gustavsson & Mehlig Reference Gustavsson and Mehlig2016; Baker et al. Reference Baker, Frankel, Mani and Coletti2017), which may affect the two-way coupling. A zero-mean flow turbulence chamber was also used by Bellani et al. (Reference Bellani, Byron, Collignon, Meyer and Variano2012) to investigate large, nearly neutrally buoyant particles. They observed moderate turbulence attenuation and pivoting of the energy spectrum. Those experimental findings were enabled by advances in particle image velocimetry (PIV) and particle tracking velocimetry (PTV), in particular their simultaneous application to capture the fluid flow and particle motion, respectively (Kiger & Pan Reference Kiger and Pan2000; Khalitov & Longmire Reference Khalitov and Longmire2002; Poelma, Westerweel & Ooms Reference Poelma, Westerweel and Ooms2006).
Under the action of gravity, particles of different sizes and densities will fall at different speeds. This has far-reaching consequences for the fluid dynamics, as soon as the still-fluid terminal velocity 
$V_{t}=\tau _{p}g$ is comparable to some velocity scale of the turbulence (i.e. when the settling velocity parameter 
$Sv_{\eta } = V_{t}/u_{\eta }$ is not vanishingly small, where 
$u_{\eta }$ is the Kolmogorov velocity scale); see Good et al. (Reference Good, Ireland, Bewley, Bodenschatz, Collins and Warhaft2014) and Petersen, Baker & Coletti (Reference Petersen, Baker and Coletti2019). Particle inertia (i.e. the finite response time to fluid fluctuations, quantified by 
$St_{\eta }$) and gravity (i.e. the drift through the flow with crossing of fluid trajectories, quantified by 
$Sv_{\eta }$) are known to have competing effects on dispersion (Wang & Stock Reference Wang and Stock1993; Berk & Coletti Reference Berk and Coletti2021); but the way their balance affects the turbulence is poorly understood. In a unique microgravity experiment, Hwang & Eaton (Reference Hwang and Eaton2006b) found that turbulence attenuation was stronger than in terrestrial gravity, pointing to the role of the particle potential energy in the fluid energy balance (Hwang & Eaton Reference Hwang and Eaton2006a).
Compared to the few experimental works, many more numerical studies have considered DNS of homogeneous isotropic turbulence two-way coupled with inertial particles (Squires & Eaton Reference Squires and Eaton1990; Elghobashi & Truesdell Reference Elghobashi and Truesdell1993; Boivin, Simonin & Squires Reference Boivin, Simonin and Squires1998; Sundaram & Collins Reference Sundaram and Collins1999; Ferrante & Elghobashi Reference Ferrante and Elghobashi2003; Frankel et al. Reference Frankel, Pouransari, Coletti and Mani2016; Saito, Watanabe & Gotoh Reference Saito, Watanabe and Gotoh2019; among others), often reporting the pivoting effect mentioned above. The cross-over scale at which fluid energy was enhanced/damped by the particles often appeared to depend on 
$St_{\eta }$, and so did the answer to the augmentation/attenuation question Poelma & Ooms (Reference Poelma and Ooms2006). Still, several authors indicated multiple reasons why these results should be interpreted with caution. First, in many cases, steady-state turbulence was forced. Lucci, Ferrante & Elghobashi (Reference Lucci, Ferrante and Elghobashi2010) argued that this is not a correct approach to study two-way coupling, independently on the type of forcing: spectral-space forcing at low wavenumbers causes wide amplitude oscillations of the turbulent kinetic energy, which cannot be distinguished from the particle back-reaction, while physical-space forcing continually modifies the spectrum to maintain a constant turbulent kinetic energy, opposing the pivoting action of the particles. Simulations of decaying turbulence are free from this problem, but suffer limitations similar to grid turbulence experiments, and can be compared only qualitatively against steady-state measurements. Second, the common point-particle forcing on the fluid computational grid presents well-known technical challenges (Eaton Reference Eaton2009; Brandt & Coletti Reference Brandt and Coletti2022). Recent efforts to address the issue show promising results (Capecelatro & Desjardins Reference Capecelatro and Desjardins2013; Subramaniam et al. Reference Subramaniam, Mehrabadi, Horwitz and Mani2014; Gualtieri et al. Reference Gualtieri, Picano, Sardina and Casciola2015; Horwitz & Mani Reference Horwitz and Mani2016, Reference Horwitz and Mani2018; Ireland & Desjardins Reference Ireland and Desjardins2017; Balachandar, Liu & Lakhote Reference Balachandar, Liu and Lakhote2019), but one-to-one comparisons with experiments are lacking. Third, most of those studies neglected the effect of gravity. The few DNS studies considering settling particles have reported contrasting outcomes. For the same range of parameters (density ratios of 
${O}(10^{3})$) and volume fractions of 
${O}(10^{-5})$, Bosse, Kleiser & Meiburg (Reference Bosse, Kleiser and Meiburg2006) found the turbulent kinetic energy to be attenuated; Frankel et al. (Reference Frankel, Pouransari, Coletti and Mani2016) and Rosa, Pozorski & Wang (Reference Rosa, Pozorski and Wang2020) found that it was augmented; and Monchaux & Dejoan (Reference Monchaux and Dejoan2017) found that it was almost unaffected.
Arguably, particle-resolved simulations (Tenneti & Subramaniam Reference Tenneti and Subramaniam2014) are needed to fully capture the particle-turbulence interaction. Driven by progress in numerical strategies and ever-increasing computational power, applications of this approach flourished in recent years and have allowed new insight in the inter-phase dynamics at the particle scale (e.g. Burton & Eaton Reference Burton and Eaton2005; Lucci et al. Reference Lucci, Ferrante and Elghobashi2010; Naso & Prosperetti Reference Naso and Prosperetti2010; Cisse, Homann & Bec Reference Cisse, Homann and Bec2013; Uhlmann & Doychev Reference Uhlmann and Doychev2014; Fornari, Picano & Brandt Reference Fornari, Picano and Brandt2016; Schneiders, Meinke & Schröder Reference Schneiders, Meinke and Schröder2017; Wang, Abbas & Climent Reference Wang, Abbas and Climent2017; Mehrabadi et al. Reference Mehrabadi, Horwitz, Subramaniam and Mani2018). The level of detail gained by such computations can hardly be achieved by experiments, especially because the simultaneous measurements of particle and fluid motion in three dimensions is generally beyond reach (Guala et al. (Reference Guala, Liberzon, Hoyer, Tsinober and Kinzelbach2008) and Ni et al. (Reference Ni, Kramel, Ouellette and Voth2015) being notable exceptions). Still, the range of accessible parameters for fully resolved simulations has been relatively narrow, with particles typically much larger than the Kolmogorov scales, and/or immersed in weak/decaying turbulence.
Considering the above, while the progress in the last decade is undeniable, the conclusion of the review by Balachandar & Eaton (Reference Balachandar and Eaton2010) remains topical: the ‘mechanisms of turbulence modulation and their parametric dependence are poorly understood and are wide open for fundamental investigation’. In a previous review, Poelma & Ooms (Reference Poelma and Ooms2006) recognized the multiplicity of important factors and called for systematic experimental studies in which one parameter at a time is varied; but no such effort has been carried out yet. We recently presented extensive measurements of particle–turbulence interaction in a homogeneous turbulence chamber, spanning a wide range of 
$St_{\eta }$, 
$Sv_{\eta }$, 
$Re_{\lambda }$ and 
$\varPhi _{v}$ (Petersen et al. Reference Petersen, Baker and Coletti2019; Berk & Coletti Reference Berk and Coletti2021). Due to the interdependence of the parameters and the focus on the dispersed phase, the turbulence modification was barely addressed.
Here, we present novel measurements in the same laboratory facility, where we systematically increase the loading of heavy sub-Kolmogorov particles, keeping the other input parameters (forcing of the turbulence and particles’ properties) constant. We repeat that for two values of 
$St_{\eta }$, smaller and larger than unity, keeping the same settling velocity parameter 
$Sv_{\eta }$. We leverage multi-scale, time-resolved imaging of both phases by PIV/PTV to simultaneously resolve all relevant scales of motion, and establish how the increasingly concentrated particles transform the fluid turbulence in which they settle. We report several novel findings. First, in the present range of physical parameters, the turbulent kinetic energy is increased substantially by the presence of falling particles, roughly in proportion to the energy input rate associated with their potential energy being released into the flow. Moreover, the energy-containing eddies become elongated in the direction of gravity, although the horizontal velocity fluctuations are more strongly excited than the vertical ones. Finally, the pivoting of the energy spectrum is due to the particles inducing an enhancement of the turbulence cascade. The paper is organized as follows. Section 2 describes the experimental methodology and qualifies the flow conditions. Section 3 presents the results, detailing how the turbulence is impacted at the different scales. Section 4 summarizes the findings and draws conclusions.
2. Methods
2.1. Experimental apparatus and parameters
The turbulence chamber was introduced in Carter et al. (Reference Carter, Petersen, Amili and Coletti2016) and qualified further in Carter & Coletti (Reference Carter and Coletti2017, Reference Carter and Coletti2018). Here, we only summarize the main features. It consists of a 5 m
$^{3}$ enclosure where two facing walls accommodate jet arrays activated in randomized sequence, following the original concept of Variano and co-workers (Variano, Bodenschatz & Cowen Reference Variano, Bodenschatz and Cowen2004; Variano & Cowen Reference Variano and Cowen2008; Bellani & Variano Reference Bellani and Variano2014). While several other random-jet-array facilities have since been build (Mydlarski Reference Mydlarski2017), the present chamber is the only one using air as working fluid and was specifically designed to study the interaction of inertial particles with homogeneous turbulence. It is especially suitable for this goal because in the central portion of the chamber: (i) the mean flow is much smaller than the root-mean-square (r.m.s.) velocity fluctuations; (ii) the turbulence is statistically homogeneous over a region much larger than the energetic scales; (iii) the spatial gradients of the mean velocity are negligible; and (iv) the attainable Reynolds numbers are sufficient to reach power-law (Kolmogorov Reference Kolmogorov1941) scaling of the velocity structure functions. The intensity and isotropy of the turbulence can be adjusted by varying the average jet firing time and the distance between the jet arrays, and by adding grids in front of the jets (Carter et al. Reference Carter, Petersen, Amili and Coletti2016).
Petersen et al. (Reference Petersen, Baker and Coletti2019) and Berk & Coletti (Reference Berk and Coletti2021) used this facility to investigate clustering and settling of inertial particles, and we adopt a similar set-up. The inertial particles (size-selected spherical glass micro-beads of density 
$\rho _{p} = 2500$ kg m
$^{-3}$) are dropped through a 3 m vertical chute and enter the chamber through a 152 mm circular opening in the ceiling. With respect to those previous studies, we release the particles with an hourglass system instead of a screw-feeder, ensuring continuous feeding and precise control of the mass flow rate. Therefore, we are able to adjust and systematically vary the mass loading. Visual inspection and PIV/PTV imaging confirm that the air turbulence spreads the falling particles over a large fraction of the chamber volume, with no significant variation of spatial concentration over the homogeneous turbulence region or in the field of view (FOV) of the imaging (Petersen et al. Reference Petersen, Baker and Coletti2019). The mean vertical velocity of the air remains much smaller 
$({\lessapprox }10\,\%)$ than its r.m.s. fluctuations for all considered particle volume fractions.
We consider two combinations of particle and turbulence properties (in single-phase realization), summarized in table 1. In both cases, the particles have sub-Kolmogorov size and lie in the 
$St_{\eta }={O}(1)$ range for which clustering is expected. The response time 
$\tau _{p}$ is calculated with the Schiller–Naumann correction (Clift, Grace & Weber Reference Clift, Grace and Weber2005). The Reynolds number based on the still-air terminal velocity 
$Re_{p,t} = d_{p}V_{t}/\nu$ (where 
$\nu$ is the air kinematic viscosity) is smaller than unity, thus negligible particle wakes are expected. The ratio 
$u_{1}'/u_{3}'$ indicates approximate large-scale isotropy. The integral scale and the turbulence dissipation are obtained from two-point velocity correlations and second-order velocity structure functions, as discussed in § 3. To isolate the effect of particle inertia and gravity, the settings are chosen to yield two different 
$St_{\eta }$ and similar 
$Sv_{\eta }$. Under fixed gravity, this necessarily requires varying 
$Re_{\lambda }$ between the cases. While we expect the results to be more sensitive to 
$St_{\eta }$ than 
$Re_{\lambda }$ over the considered range, we will refrain from inferring quantitative Stokes number trends from our data. In fact, the numerical simulations by Tom & Bragg (Reference Tom and Bragg2019) indicate that for the present range of 
$St_{\eta }$, the settling dynamics can be very sensitive to 
$Re_{\lambda }$. Overall, both cases show similar behaviours, and comparisons will rather serve as an indication that the conclusions hold over the present range of parameters. For simplicity, throughout the paper we will refer to both configurations as the 
$St_{\eta }=0.3$ case and 
$St_{\eta }=2.6$ case, respectively. The volume fraction 
$\varPhi _{v}$ and mass fraction 
$\varPhi _{m}=\varPhi _{v} \rho _{p}/\rho _{f}$ (where 
$\rho _{f}$ is the air density) span ranges in which the back-reaction on the turbulence is expected to vary from marginal to substantial (Elghobashi Reference Elghobashi1994; Balachandar & Eaton Reference Balachandar and Eaton2010). As the particles alter the flow, several observables will change, including the fluid velocity fluctuations and the fall speed of the particles themselves. In order to characterize the system, the parameters reported in table 1 are defined for the baseline configuration unaffected by two-way coupling.
Table 1. Properties of the particles and the turbulence (for the unladen flow case) for both investigated configurations. Here, 
$d_{p}$ and 
$\tau _{p}$ are the particle diameter and response time; 
$Re_{p,t}$ is the particle Reynolds number based on the terminal velocity; 
$St_{\eta }$ and 
$Sv_{\eta }$ are the Stokes number and settling velocity parameter based on Kolmogorov scales; 
$u_{1}'$ and 
$u_{3}'$ are the r.m.s. fluid velocity fluctuations in the horizontal and vertical directions, respectively; 
$L_{1,1}$ is the integral scale of the turbulence in the horizontal direction; 
$\epsilon$ is the turbulent dissipation calculated from the second-order velocity structure function; 
$\eta$ and 
$\tau _{\eta }$ are the Kolmogorov length and time scales; 
$Re_{\lambda }$ is the Taylor-microscale Reynolds number of the turbulence; and 
$\varPhi _{v}$ and 
$\varPhi _{m}$ are the volume and mass fraction of the particles, respectively.
2.2. Measurement techniques
In the following, 
$x_{1}$ indicates the horizontal direction parallel to the jet axes, 
$x_{2}$ is horizontal and perpendicular to 
$x_{1}$, and 
$x_{3}$ is vertical; the respective fluid velocity components are 
$U_{1}$, 
$U_{2}$ and 
$U_{3}$. Fluctuating fluid velocities are denoted by 
$u_{i} = U_{i} - \overline {U_{i}}$, where the subscript 
$i$ indicates the 
$i$th component, and the overbar indicates temporal averaging. Angle brackets indicate space–time averaging, and a prime indicates r.m.s fluctuations. The subscript 
$p$ denotes quantities related to the particles.
We perform simultaneous PIV on the fluid phase and PTV on the particle phase. The method is similar to that described in Petersen et al. (Reference Petersen, Baker and Coletti2019) and Berk & Coletti (Reference Berk and Coletti2021); only salient differences are addressed here. Time-resolved imaging is carried out using an Nd:YLF single-pulse laser (Photonics, 30 mJ pulse
$^{-1}$) operated at 4 kHz to illuminate a vertical 
$(x_{1}, x_{3})$-plane at the centre of the chamber. Two synchronized CMOS cameras (Phantom VEO, active sensor of 1280 by 1280 pixels) image the in-plane air motion through DEHS tracers seeding the chamber, and track the inertial particles about 90 cm below the ceiling opening. The images are separated into ‘tracers only’ and ‘particles only’, with a routine that we described in detail in Petersen et al. (Reference Petersen, Baker and Coletti2019) and applied recently to various particle-laden turbulent flows (Fong et al. Reference Fong, Amili and Coletti2019; Berk & Coletti Reference Berk and Coletti2020; Baker & Coletti Reference Baker and Coletti2021).
Figure 1 illustrates a sample instantaneous realization of the two-phase flow. The fluid velocity is measured by PIV via iterative cross-correlation of successive image pairs, with a final interrogation window size of 
$24\times 24$ pixels and 50 % overlap. The imaging and PIV specifications are reported in table 2. To qualify the spatio-temporal scales resolved by the measurements, we refer to the single-phase flow properties in table 1. Both cameras are operated simultaneously and capture two nested fields of view, enhancing the dynamic range of the measurements. One camera mounts a 105 mm Nikon lens (f# 5.6), yielding a FOV large enough to capture the integral length scales. The other mounts a 200 mm Nikon lens (ff# 4) and approximately resolves the Kolmogorov length scales: the vector spacing is about 
$\eta$ and 
$4\eta$ at the lower and higher 
$Re_{\lambda }$, respectively. The acquisition frequency resolves 19 to 96 times the Kolmogorov time scales, again depending on 
$Re_{\lambda }$. Besides considerations about the scales of the flow, the camera resolutions and inter-frame separations are optimized so that the tracer displacements are suitable for velocimetry in both cameras: these are approximately 4 and 9 pixels for the large and small FOV, respectively. For each volume fraction, 10 runs of 43 000 successive images are obtained, for a total measurement time of 100 s or approximately 500 integral time scales (defined below).
Figure 1. (a) Example of nested raw images of the particle-laden flow, captured in the large (left) and small (right) field of view (FOV) for 
$St_{\eta }=2.6$. (b) The corresponding in-plane fluid velocity in m s
$^{-1}$ (colour contours) and particle positions (white dots, not in scale size).
Table 2. Imaging parameters for the small and large field of view (FOV).
The detected particle count in the illuminated volume is used to calculate the volume fraction 
$\varPhi _{v}$; see Petersen et al. (Reference Petersen, Baker and Coletti2019) and Fong et al. (Reference Fong, Amili and Coletti2019), where the approach could be validated in a vertical channel flow with known particle mass loading. Over the present range of 
$\varPhi _{v}$, the fluid turbulence can be characterized with sufficient image quality for successful processing. This is based on the extensive analysis of the phase separation by Petersen et al. (Reference Petersen, Baker and Coletti2019), where similar regimes and volume fractions were reached. In particular, the systematic variation of 
$\varPhi _{v}$ in the present study allows us to verify that the occurrence of non-valid PIV vectors is not correlated with the presence of inertial particles in a given interrogation window. The PIV random error in the small FOV (where the spatial resolution is sufficient to capture the Kolmogorov scales of the single-phase turbulence) is estimated by extrapolating the two-point correlation to vanishing separation (Adrian & Westerweel Reference Adrian and Westerweel2011) and is found to be less than 1 % of the velocity variance. The dominant source of uncertainty on the flow statistics is thus the finite sample size. Convergence tests show that each run is well converged for all considered observables, while run-to-run variation is larger. In the following, where appropriate, the standard deviation of the various runs will be used to indicate error bars.
A specific feature of the present zero-mean flow facility is the spatial homogeneity over scales much larger than the integral scale, for both fluid and particle statistics (Carter et al. Reference Carter, Petersen, Amili and Coletti2016; Carter & Coletti Reference Carter and Coletti2017, Reference Carter and Coletti2018; Petersen et al. Reference Petersen, Baker and Coletti2019). This allows drawing more general conclusions compared to systems in which wall proximity and spatial gradients cause dependence on the boundary conditions. Moreover, the considered particles have reached terminal velocity well before entering the region of interest, as indicated by their negligible mean vertical acceleration (see Berk & Coletti Reference Berk and Coletti2021). Finally, the lack of significant mean flow (especially small in the vertical direction) further limits the possible influence from portions of the volume outside the region of interest. In other words, in the considered central region of the chamber, the integral scale of the turbulence sets the scale over which spatial gradients of the flow properties have an influence. Therefore, the flow condition outside the FOV are not expected to alter the statistics measured in the imaging window, which is of the order of the integral scale and about one metre removed from the chamber walls.
Figure 2 displays the time record of instantaneous 
$\varPhi _{v}$, 
$u_{1}'$ and 
$u_{1,p}'$ during a sample run for the case 
$St_{\eta } = 2.6$ at volume fraction 
$4\times 10^{-5}$. Beyond the large-scale excursions, the plots confirm the steady-state behaviour of both phases. We will use the local particle concentration, velocity and fluid velocity at the particle location (interpolated as in Petersen et al. Reference Petersen, Baker and Coletti2019; Berk & Coletti Reference Berk and Coletti2021), extracted from the small-FOV measurements, to estimate the extra dissipation caused by the particles; see § 3.2. The local concentration is obtained by Voronoi tessellation of the particle field (Monchaux, Bourgoin & Cartellier Reference Monchaux, Bourgoin and Cartellier2010).
Figure 2. Temporal evolution of spatially averaged quantities during one experimental run for the case 
$St_{\eta } = 2.6$ with mean volume fraction 
$4\times 10^{-5}$. The instantaneous volume fraction 
$\varPhi _{v}$, the r.m.s. of the horizontal flow velocity fluctuations, 
$u_{1}'$, and the r.m.s. of the horizontal particle velocity fluctuations, 
$u_{1p}'$, are displayed. The data are from the large FOV.
3. Results
3.1. Large-scale turbulence properties
The results in this section are obtained from the large-FOV measurements. We begin by considering the effect of the particles on the r.m.s. velocity fluctuations of the turbulence. Figure 3(a) shows the laden-to-unladen ratio for the horizontal component, 
$u_{1}'/u_{1,\varPhi =0}'$, for the different volume fractions. A monotonically increasing trend is apparent for both Stokes numbers, with the larger 
$St_{\eta }$ case displaying a stronger increase. The vertical component also increases with volume fraction, but less than the horizontal one. Indeed, the anisotropy ratio 
$u_{1}'/u_{3}'$ grows significantly over the considered range of volume fractions (figure 3b), the growth being again stronger for 
$St_{\eta } = 2.6$ than for 
$St_{\eta } = 0.3$.
Figure 3. (a) The r.m.s. of the horizontal velocity fluctuation as a function of the particle volume fraction, normalized by the one measured in the unladen flow. (b) Ratio between the horizontal and the vertical r.m.s. fluctuations, as a function of the particle volume fraction. Squares indicate the 
$St_{\eta } = 2.6$ case, and circles indicate the 
$St_{\eta } = 0.3$ case. The data are from the large FOV.
The enhancement of the r.m.s. fluid velocity fluctuations in both directions indicates that the turbulent kinetic energy (TKE) is augmented significantly by the presence of the particles, as we will discuss later. Moreover, the fact that the particle-laden turbulence fluctuations are more intense in the horizontal than in the vertical direction is at odds with previous numerical studies: simulations of decaying (Ferrante & Elghobashi Reference Ferrante and Elghobashi2003) and forced (Bosse et al. Reference Bosse, Kleiser and Meiburg2006) homogeneous turbulence found that the addition of settling particles resulted in anisotropy ratios 
$u_{1}'/u_{3}'<1$. Also, the grid turbulence measurements of Geiss et al. (Reference Geiss, Dreizler, Stojanovic, Chrigui, Sadiki and Janicka2004) and Poelma et al. (Reference Poelma, Westerweel and Ooms2007) showed stronger fluid fluctuations in the direction of gravity than in the transverse direction. However, direct comparison with these experiments is hampered by the differences in parameters compared to the present case: the particles in Geiss et al. (Reference Geiss, Dreizler, Stojanovic, Chrigui, Sadiki and Janicka2004) were quasi-ballistic, while those in Poelma et al. (Reference Poelma, Westerweel and Ooms2007) were not small compared to the Kolmogorov scales, and produced wakes with large momentum deficit, which heavily impacted their weak decaying turbulence.
The influence of the settling particles on the large-scale turbulence properties is further illustrated by the integral scales, obtained from the spatial autocorrelation of the velocity fluctuations:
\begin{equation} \rho_{ii}(\boldsymbol{r}) = \frac{\langle u_{i}(\boldsymbol{r})\,u_{i} (\boldsymbol{x}+\boldsymbol{r})\rangle}{\langle u_{i}(\boldsymbol{x})^{2}\rangle}, \end{equation}
where 
$x$ and 
$r$ are the position and separation vectors, respectively, and no index summation is implied. In general, the integral length scales associated with fluctuations 
$u_{i}$ over separations 
$r_{j}$ are defined as
\begin{equation} L_{i,j} = \int_{0}^{r_{0}} \rho_{ii}(r_{j})\,{\rm d} r, \end{equation}
where 
$r_{0}$ is the first zero crossing of 
$\rho _{ii}$. In practice, the extent of the integration is limited by the size of the imaging window, and therefore we take as conventional estimate the separation at which the correlation function drops below 0.5. Other common methods such as extrapolating an exponential fit yield similar results. As in Carter & Coletti (Reference Carter and Coletti2017), we define the integral and transverse scales as 
$L_{L,1}=L_{1,1}$ and 
$L_{T,1}=L_{1,2}$, respectively. The subscript 1 indicates that we use the horizontal component of the velocity, the vertical one returning similar trends.
In figure 4(a), we show the autocorrelation of the horizontal velocity fluctuations along the horizontal separation, 
$\rho _{11}$, for the 
$St_{\eta } = 2.6$ case. The correlation decays faster with increasing volume fraction, which results in a drop of the longitudinal integral scale 
$L_{L,1}$ as shown in figure 4(b). As we already remarked, 
$u_{1}'$ increases with volume fraction, so the shrinking of the integral scale suggests that the falling particles subtract energy from the large scales and inject it at the small scales, with a positive global balance (i.e. a net increase of TKE). This view will be confirmed in the following sections.
Figure 4. (a) Autocorrelation coefficient of the horizontal velocity fluctuations at the different volume fractions for 
$St_{\eta } =2.6$. (b) Corresponding longitudinal integral scale, normalized by the unladen case, as a function of the particle volume fraction. The data are from the large FOV.
To better appreciate the particle influence on the large-scale spatial structure of the turbulence, we plot in figure 5(a) the 0.5 contours of 
$\rho _{11}$ in the 
$(r_{1},r_{2})$-plane. With increasing volume fraction, the horizontal extent of the correlation area shrinks, decreasing 
$L_{L,1}$ as mentioned; instead, its vertical extent grows, increasing the transverse length scale 
$L_{T,1}$. This is quantified in figure 5(b), showing the 
$L_{L,1}/L_{T,1}$ ratio as a function of volume fraction. For the single-phase case, the ratio departs from the canonical value 2 due to incomplete homogeneity and isotropy across the FOV. The ratio shows a clear decreasing trend with increasing particle loading. The 0.5 contours of the vertical velocity autocorrelation 
$\rho _{22}$ (not shown for brevity) display a similar vertical stretch, with an increase of the longitudinal length scale (
$L_{L,2} = L_{2,2}$) and decrease of the transverse length scale (
$L_{T,2}=L_{2,1}$). Taken together, these results indicate that the energy-containing eddies become elongated vertically in the presence of falling particles at increasing concentrations. This contrasts with the above observation that the horizontal fluid fluctuations are excited more than the vertical ones: in single-phase turbulence, the integral scales are stretched in the direction along which fluctuations are more intense (Carter & Coletti Reference Carter and Coletti2017). The presence of the particles clearly alters the way energy is redistributed between different directions and across scales.
Figure 5. (a) Contours of the autocorrelation coefficient 
$\rho _{11}$ in the two-dimensional scale space 
$(r_{1},r_{2})$ for the case 
$St_{\eta } = 2.6$ at different volume fractions. (b) Corresponding ratio between the longitudinal and transverse integral scales, as a function of the particle volume fraction. The data are from the large FOV.
3.2. Turbulent kinetic energy and dissipation rate
As anticipated in the previous subsection, the presence of particles in the present regime augments the turbulent kinetic energy of the carrier fluid. This is evident in figure 6, where we plot TKE calculated as
\begin{equation} {\rm TKE} =\tfrac{1}{2}(2u_{1}^{'2} + u_{3}^{'2}) \end{equation}
from the large- and small-FOV data at different volume fractions. In (3.1), we assume that the out-of-plane velocity variance is equal to the measured horizontal component. This is expected when the baseline single-phase turbulence is approximately isotropic, which is justified by the value of 
$u_{1}'/u_{3}'$ in absence of particles (table 1). The change in TKE compared to the single-phase flow is marginal for 
$\varPhi _{v}\approx 10^{-6}$, in agreement with the order-of-magnitude estimate of Elghobashi (Reference Elghobashi1994). The turbulence is augmented monotonically with increasing volume faction, leading to TKE increases of 70 % at 
$\varPhi _{v} = 4\times 10^{-5}$ for 
$St_{\eta } = 0.3$, and more than 150 % at 
$\varPhi _{v}= 5\times 10^{-5}$ for 
$St_{\eta } = 2.6$. Such a dramatic augmentation of turbulence is in stark contrast with the criterion of Gore & Crowe (Reference Gore and Crowe1991) and with several previous experimental studies: Paris (Reference Paris2001), Hwang & Eaton (Reference Hwang and Eaton2006a,Reference Hwang and Eatonb) and Tanaka & Eaton (Reference Tanaka and Eaton2008) found that solid particles attenuated air turbulence. Those authors, however, considered significantly more inertial particles (
$St_{\eta } \approx 50\unicode{x2013}100$), and in the following we reason that this is the likely cause of the different behaviour.
Figure 6. Turbulent kinetic energy normalized by the unladen case, as a function of the particle volume fraction. Squares indicate the 
$St_{\eta }= 2.6$ case, and circles indicate the 
$St_{\eta } = 0.3$ case. The data are from the small FOV.
The increase in TKE is related to the energy balance in the particle-laden turbulence, where gravitational settling plays a major role: as the particles fall, they transfer potential energy into the flow. Following Hwang & Eaton (Reference Hwang and Eaton2006a), we consider the steady-state energy budget
\begin{equation} E_{j} + E_{g} = \epsilon + \epsilon_{p}, \end{equation}
where the left-hand and right-hand sides represent the energy per unit time, and mass injected and dissipated in the fluid control volume, respectively. Here, 
$E_{j}$ is the forcing from the jets, which is independent of the mass loading and equals the energy dissipation rate in the single-phase cases. The energy per unit time released by the falling particles is 
$E_{g} = \varPhi _{m} \langle U_{3,p} \rangle g$, analogous to estimates for rising bubbles (Riboux, Risso & Legendre Reference Riboux, Risso and Legendre2010; Risso Reference Risso2018). The viscous dissipation of turbulence energy due to the cascade from larger to smaller eddies, 
$\epsilon$, is estimated by evaluating directly the in-plane velocity gradients resolved by PIV from the small FOV (see De Jong et al. Reference De Jong, Cao, Woodward, Salazar, Collins and Meng2009; Carter et al. Reference Carter, Petersen, Amili and Coletti2016). For single-phase turbulence in this facility imaged with similar parameters, such a direct estimate was shown to agree closely with the one based on velocity structure functions (Carter et al. Reference Carter, Petersen, Amili and Coletti2016; Carter & Coletti Reference Carter and Coletti2017, Reference Carter and Coletti2018). The presence of the particles leads to an extra dissipation 
$\epsilon _{p}$ associated with their interaction with the fluid, in particular the boundary layer around the particles, which is far below the PIV resolution limits. This can be estimated with the classic expression proposed by Elghobashi & Abou-Arab (Reference Elghobashi and Abou-Arab1983) and used in several later studies (Rogers & Eaton Reference Rogers and Eaton1991; Kulick et al. Reference Kulick, Fessler and Eaton1994; Hwang & Eaton Reference Hwang and Eaton2006a; Sahu, Hardalupas & Taylor Reference Sahu, Hardalupas and Taylor2016; among others):
\begin{align} \epsilon_{p} &= \frac{\langle C \rangle}{\rho_{f}\tau_{p}}\,(\langle u_{i,fp}u_{i,fp} \rangle - \langle u_{i,fp}u_{i,p} \rangle) + \frac{1}{\rho_{f}\tau_{p}}\,(\langle cu_{i,fp}u_{i,fp} \rangle - \langle cu_{i,fp}u_{i,p} \rangle) \nonumber\\ &\quad +\frac{1}{\rho_{f}\tau_{p}}\,(\langle U_{i,fp} \rangle - \langle U_{i,p} \rangle) \langle cu_{i,fp} \rangle. \end{align}
Here, 
$C$ and 
$c$ correspond to the local particle concentration and the local particle concentration fluctuation, and 
$u_{fp}$ is the fluid velocity fluctuation at the particle location.
Figure 7(a) shows 
$\epsilon _{p}$ as a function of volume fraction (normalized by the unladen turbulence dissipation, 
$\epsilon _{\varPhi =0}$) evaluated in two ways: indirectly from (3.4) as 
$E_{j} + E_{g} - \epsilon$, and directly from (3.5). The results are displayed for 
$St_{\eta } = 2.6$, the outcome for 
$St_{\eta } = 0.3$ being analogous. Both methods show a similar trend, although with discrepancies at higher loadings. We remark that an accurate evaluation of 
$\epsilon _{p}$, following either method, is highly challenging for several reasons. First, the evaluation of 
$\epsilon$ in presence of suspended particles is even more difficult than in unladen turbulence (where it is already notoriously difficult). Second, (3.4) is based on the assumption that the forcing from the jets is independent of the presence of the particles, while the latter could indirectly influence the amount of kinetic energy injected at large scales in the investigated region. Third, and perhaps most important, (3.5) is based on the same assumptions as two-way-coupled point-particle models (Hwang & Eaton Reference Hwang and Eaton2006a), which have well-known shortcomings resulting in epistemic uncertainties that are hard to quantify. Indeed, the error bars representing statistical variability are not displayed for 
$\epsilon _{p}$ based on (3.5), as this is not the main source of error. Despite those limitations, both estimates of 
$\epsilon _{p}$ in figure 7(a) lead to a clear observation: at the larger considered volume fractions, the extra dissipation due to the particles becomes of the order of the baseline dissipation. Likewise, the energy input 
$E_{g}$ associated with gravitational settling, shown in figure 7(b) for the same case, also increases with 
$\varPhi _{v}$ and becomes comparable to 
$\epsilon _{\varPhi =0}$. The increase is due not only to the mass loading itself, but also to the mean settling velocity, which grows by 
$110\,\%$ over the considered range of volume fractions. This may be due to collective drag effects (as indicated by Bosse et al. Reference Bosse, Kleiser and Meiburg2006) or preferential sweeping (which was shown to remain important in recent two-way-coupled simulations by Tom, Carbone & Bragg Reference Tom, Carbone and Bragg2022). The analysis of the settling enhancement is beyond the scope of the present work and will be addressed in detail in a separate study. The important observation is that gravitational settling contributes majorly to the amount of energy injected and dissipated in the fluid per unit time. Whether TKE also increases is related to the time scale over which the dissipation is expressed. Tanaka & Eaton (Reference Tanaka and Eaton2010) pointed out how both the turbulence integral time scale and the particle response time may be appropriate candidates. The latter does not vary with the particle loading, while the former (taken as 
$L_{L}/u'$, averaging between horizontal and vertical components) increases by at most 
${\approx }30\,\%$ for the higher 
$\varPhi _{v}$. Thus, as the time scale is not drastically changed by the presence of the particles, at steady state we can expect the turbulent kinetic energy of the system to rise with particle loading, as observed.
Figure 7. (a) Particle-induced dissipation 
$\epsilon _{p}(a)$, and (b) gravitational energy input rate 
$E_{g}(b)$, both normalized by the turbulent dissipation of the unladen flow, for 
$St_{\eta }= 2.6$. In (a), upward-pointing triangles and downward-pointing triangles represent estimates of 
$\epsilon _{p}$ from (3.5) and (3.4), respectively. The data point from (3.5) at the highest concentration was discarded, due to the uncertainty in estimating the high-order terms. The data are from the small FOV.
This type of reasoning was used by Hwang & Eaton (Reference Hwang and Eaton2006a) to argue that settling contributed to augmenting turbulence also in their system; but that was not sufficient to offset the ability of their quasi-ballistic particles to act as sinks of fluid momentum. Their microgravity experiments (Hwang & Eaton Reference Hwang and Eaton2006b) confirmed such an argument, showing larger turbulence attenuation compared to an equivalent particle-laden turbulence regime under terrestrial gravity. Recently, Saito et al. (Reference Saito, Watanabe and Gotoh2019) used two-way coupled simulations to show that massively inertial point-particles (akin to fixed obstacles in the flow) severely attenuated turbulence, while the effect was much weaker for less inertial particles that more closely follow the flow. According to such a view, the difference between the present observations and those from Eaton and co-workers (Paris Reference Paris2001; Hwang & Eaton Reference Hwang and Eaton2006a,Reference Hwang and Eatonb; where 
$St_{\eta } \approx 50\unicode{x2013}100$) can be attributed to the smaller Stokes number in our system. Indeed, the one previous experimental study with solid particles in air at 
$St_{\eta }={O}(1)$ also found turbulence augmentation (Yang & Shy Reference Yang and Shy2005).
The natural question is why the settling of 
$St_{\eta }={O}(1)$ particles excites turbulence more than the settling of quasi-ballistic ones. A possible explanation is suggested by a comparison with bubble-laden turbulence. Mazzitelli et al. (Reference Mazzitelli, Lohse and Toschi2003a,Reference Mazzitelli, Lohse and Toschib) found numerically that small rising bubbles accumulate in the downward side of turbulent eddies due to the lift force (as recently confirmed experimentally by Salibindla et al. Reference Salibindla, Masuk, Tan and Ni2020), transferring momentum upwards, which enhances their ability to modify the turbulence. Likewise, small heavy particles with 
$St_{\eta }={O}(1)$ form clusters that oversample downward regions of the flow (Wang & Stock Reference Wang and Stock1993; Baker et al. Reference Baker, Frankel, Mani and Coletti2017; Petersen et al. Reference Petersen, Baker and Coletti2019), and therefore they may locally enhance turbulent fluctuations more effectively than ballistic particle that do not cluster. This interpretation is consistent with the larger TKE increase for particles with 
$St_{\eta }>1$, which experience stronger clustering than those with 
$St_{\eta }<1$ under the action of gravity (Matsuda, Onishi & Takahashi Reference Matsuda, Onishi and Takahashi2017; Petersen et al. Reference Petersen, Baker and Coletti2019).
As mentioned in the Introduction, although numerous authors investigated two-way coupling of point-particles, only a few of them included gravity in forced homogeneous turbulence. Bosse et al. (Reference Bosse, Kleiser and Meiburg2006) found turbulence attenuation by particles with 
$St_{\eta } \approx Sv_{\eta }\approx 1$, while Frankel et al. (Reference Frankel, Pouransari, Coletti and Mani2016) considered faster falling particles (
$St_{\eta }=1.6$ and 
$Sv_{\eta }=6.8$) and reported an increase of TKE consistent with our observations. Besides the limitations of the point-particle approach, these comparisons with simulations also support the view that gravitational settling enhances the turbulence activity in the carrier flow.
3.3. Turbulence modification across scales
The change in TKE is the result of turbulence modification over the entire spectrum. Figure 8 shows the longitudinal second-order velocity structure function, which, compared to the Fourier spectrum, is less sensitive to finite-sample biases (Rensen et al. Reference Rensen, Luther and Lohse2005; De Jong et al. Reference De Jong, Cao, Woodward, Salazar, Collins and Meng2009) and to missing data that may result from the presence of particles (Poelma & Ooms Reference Poelma and Ooms2006). We consider the second-order structure function
\begin{equation} D_{ii}(\boldsymbol{r}) = \langle[u_{i}(\boldsymbol{x}+ \boldsymbol{r}) - u_{i}(\boldsymbol{x})]^{2}\rangle, \end{equation}
where 
$\boldsymbol {x}$ is the position vector, and 
$\boldsymbol {r}$ is the separation vector. We focus on the longitudinal structure functions, where the velocity component is parallel to the separation vector, particularly on 
$D_{11}(r_{1})$, 
$D_{22}(r_{2})$ and the transverse structure functions showing the same behaviour. According to Kolmogorov theory (Kolmogorov Reference Kolmogorov1941), (3.6) scales as 
$(\epsilon r)^{{2}/{3}}$ for separations in the inertial range. Although the considered range of 
$Re_{\lambda }$ possesses limited scale separation, Carter et al. (Reference Carter, Petersen, Amili and Coletti2016) and Carter & Coletti (Reference Carter and Coletti2017) showed, and the present data confirm, that the single-phase turbulence in this facility approximately follows Kolmogorov scaling in the range 
$r \approx 20\eta \unicode{x2013} 50\eta$. The addition of particles alters this picture profoundly. In keeping with the pivoting effect, the energy becomes more evenly distributed across the scales as the particle loading increases: the dispersed phase intercepts part of the fluid momentum at the large scales, and injects it back into the fluid at small scales. The effect is already visible for 
$St_{\eta } = 0.3$ and is apparent for 
$St_{\eta } = 2.6$.
Figure 8. Longitudinal second-order structure functions of the 
$u_{1}$ component for the different volume fractions for (a) 
$St_{\eta } = 0.3$ and (b) 
$St_{\eta } = 2.6$. The dashed line indicates the 
$r^{2/3}$ scaling predicted by Kolmogorov (Reference Kolmogorov1941) in the inertial range. The data are from the small FOV.
The cross-over separation 
$r_{c}$, i.e. the separation corresponding to the pivoting point, does not vary significantly with 
$\varPhi _{v}$, and lies in the ranges 
$30\eta$–
$40\eta$ and 
$40\eta$–
$50\eta$ for the lower and higher 
$St_{\eta }$, respectively. This difference could be a 
$St_{\eta }$ effect or a 
$Re_{\lambda }$ effect. The latter would be consistent with the proposal of Sundaram & Collins (Reference Sundaram and Collins1999) that 
$r_{c}$ scales with 
$\lambda$, while the former is at odds with the notion that 
$r_{c}$ decreases with 
$St_{\eta }$ (Poelma & Ooms Reference Poelma and Ooms2006). This conclusion, however, was based on simulations neglecting gravity, which here is crucial to the two-way coupling. Determining a specific trend for 
$r_{c}$ would require further measurements with different 
$Sv_{\eta }$ and a broader range of 
$St_{\eta }$. We do note that the range of scales where the fluctuating fluid energy is increased roughly corresponds to the size of the particle clusters in this regime (Petersen et al. Reference Petersen, Baker and Coletti2019).
Finally, the redistribution of the turbulent energy with increasing particle loading can be appreciated by considering the nonlinear inter-scale energy transfer rate, 
$\varPi$. Starting from the Kármán–Howarth–Monin equation (Monin, Yaglom & Lumley Reference Monin, Yaglom and Lumley1975), we recently performed such an analysis for the single-phase turbulence in this same facility (Carter & Coletti Reference Carter and Coletti2018; see also Lamriben, Cortet & Moisy Reference Lamriben, Cortet and Moisy2011; Gomes-Fernandes, Ganapathisubramani & Vassilicos Reference Gomes-Fernandes, Ganapathisubramani and Vassilicos2015; Portela, Papadakis & Vassilicos Reference Portela, Papadakis and Vassilicos2017; among others). Briefly, given two points 
$\mathbf {x}$ and 
$\mathbf {x}^{\prime \prime }$ associated with velocity fluctuations 
$u_{i}$ and 
$u_{i}^{\prime \prime }$ and separated by 
$\mathbf {r}$, one considers the difference 
$\delta u_{i} = u_{i} - u_{i}^{\prime \prime }$ and the associated energy 
$\delta q_{i} =\delta u_{i}\delta u_{i}$. (Assuming small-scale axisymmetry, as in Carter & Coletti (Reference Carter and Coletti2018), or isotropy leads to similar values and analogous conclusions.) The radial component of the nonlinear inter-scale energy transfer rate, 
$\varPi _{r} = \frac {1}{4}[({\partial }/{\partial r})(\langle \delta u_{i}\,\delta q_{i}^{2}\rangle r_{i}/r) + 1/r^{2} \langle \delta u_{i}\,\delta q_{i}^{2} \rangle r_{i}]$, is the only one making a net contribution in spherical coordinates. Thus, over the inertial range in single-phase turbulence, the average of 
$\varPi _{r}$ across all orientations in scale space is expected to balance the two-point kinetic energy dissipation 
$\epsilon _{r}$ (which in the present homogeneous flow coincides with the standard dissipation 
$\epsilon$). This is confirmed in figure 9, where the ratio 
$\varPi _{r}/\epsilon _{\varPhi =0}$ is plotted as a function of 
$r$ for 
$St_{\eta } = 2.6$. The negative sign of 
$\varPi _{r}$ indicates the expected direction of the energy cascade from larger to smaller scales, and the values of the ratio close to 
$-1$ for 
$r/\eta \approx 20\unicode{x2013}100$ reflect the 
$\varPi \unicode{x2013}\epsilon$ equilibrium in the inertial range of the unladen turbulence. With increasing particle loading, the magnitude of 
$\varPi _{r}$ in the same range grows significantly. This indicates that with increasing loading, the energy injected at the large scales is being transferred more effectively down the spectral pipeline towards the smaller scales. On the one hand, this is consistent with the observation that particles cause the energy to pivot towards a more even distribution across all scales, as illustrated in figure 8. On the other hand, it points to a specific mechanism by which such redistribution is realized, i.e. the enhancement of the energy cascade. The precise dynamics by which this is realized (e.g. by clusters accelerating the breakdown of large eddies, rather than high-concentration sheets causing instabilities; Kasbaoui, Koch & Desjardins Reference Kasbaoui, Koch and Desjardins2019) remains to be clarified, but it is likely specific to the present regime. In fact, particle-resolved DNS (PR-DNS) from Schneiders et al. (Reference Schneiders, Meinke and Schröder2017) found that Kolmogorov-scale particles dampened the inter-scale energy flux in the absence of gravity.
Figure 9. Nonlinear inter-scale energy transfer rate, normalized by the turbulent dissipation and plotted as a function of the scale separation 
$r$, at the different particle volume fractions and for the 
$St_{\eta } = 2.6$ case. The dissipation from the unladen case is used to normalize all line plots; using 
$\epsilon$ from the individual particle-laden cases corroborates the observed trend of inter-scale transfer increasing with volume fraction. The data are from the small FOV.
In this multi-scale and two-way coupled system, the question of the inter-scale energy transfer is complex, and other interpretations of the results are possible. For example, the enhancement of 
$\varPi _{r}$ might be seen as a mere consequence of the energy added to the fluid by the falling particles, rather than by an intensification of the cascade. In such a scenario, however, one would expect a scale-independent shift of the energy spectrum, as opposed to the observed pivoting. Likewise, the increase of fluctuating energy at small scales could be due solely to the local action of the particles. If this were the case, however, such energy augmentation would be limited to scales much smaller than the range that appears affected. In conclusion, the present interpretation appears consistent with the ensemble of the experimental observations, but further data and analysis are required to make it more conclusive.
4. Conclusions
We have performed a systematic experimental study of turbulence modification by sub-Kolmogorov heavy particles. These are in the range 
$St_{\eta }={O}(1)$, 
$\rho _{p}/\rho _{f} = {O}(10^{3})$ that has attracted attention for decades, due to its relevance to countless applications and the rich interaction with the flow. With 
$Re_{\lambda } \approx 150\unicode{x2013}300$, the scale separation of the turbulence is enough to display the hallmarks of Kolmogorov's phenomenology, also possessing a high degree of homogeneity over scales much larger than the energetic eddies. The settling velocity parameter 
$Sv_{\eta }={O}(1)$ is the natural consequence of terrestrial gravity, and separates the present regime from the zero-gravity cases that have formed the bulk of the numerical studies in the literature. By varying incrementally the particle volume fraction, 
$\varPhi _{v} \approx 10^{-6}\unicode{x2013}5\times 10^{-5}$, while keeping all other parameters constant, we have isolated the profound flow modifications caused by the two-way coupling between both phases. We considered both cases 
$St_{\eta } < 1$ and 
$St_{\eta } > 1$, which behave similarly, the latter producing more marked changes in the carrier phase.
We find that the presence of settling inertial particles leads to the horizontal contraction and vertical elongation of the correlation length of the turbulent fluctuations, in agreement with the observation of Ferrante & Elghobashi (Reference Ferrante and Elghobashi2003) that falling particles stretch vortical structures. The r.m.s. velocity fluctuations, however, are more strongly excited in the horizontal than in the vertical direction. Both components are significantly augmented, and overall the TKE is dramatically increased (roughly doubled) at the higher loadings. This is attributed to the role of gravity: falling particles release their potential energy to the fluid, proportionally to their mass fraction and settling velocity. At the largest volume fractions considered, such energy input rate is comparable to the hydrodynamic forcing of the turbulence (in the present case, from actuated jets).
The velocity structure functions confirm that the presence of particles leads to a redistribution of energy across the spectrum: the turbulence activity is intensified at the small scales and dampened at the large scales, with a cross-over scale around 
$40\eta$. The two-point analysis shows further that the inter-scale energy flux is enhanced with increasing particle loading. Thus the pivoting of the energy spectrum is rooted in an intensification of the direct energy cascade. The identification and quantification of this process will be important to devise successful (subgrid-scale and/or point-particle) models able to capture the turbulence modification.
In most previous experimental studies, which focused on massively inertial particles (
$St_{\eta }={O}(10^{2})$), turbulence was found to attenuate TKE. Unlike those cases, the present class of particles exhibits distinct phenomena such as clustering and enhanced settling (as measured in the present facility by Petersen et al. (Reference Petersen, Baker and Coletti2019) and Berk & Coletti (Reference Berk and Coletti2021)). The difference in the net variation of TKE depending on the particle inertia may then be due to the time scale associated to the energy input, which in turn is related to the particle dynamics.
While we mainly attribute to 
$St_{\eta }$ the different TKE modification with respect to previous studies, the particle size may also play a role: in most past experiments that found a reduction of TKE, 
$dp/\eta \approx 1$ or larger (Hwang & Eaton Reference Hwang and Eaton2006a; Tanaka & Eaton Reference Tanaka and Eaton2008; Bellani et al. Reference Bellani, Byron, Collignon, Meyer and Variano2012), while here 
$dp/\eta = {O}(10^{-1})$. To discern and quantify geometric effects on the flow, particle-resolving simulations are needed, although the present regime is especially challenging to PR-DNS.
It is important to point out that the present findings are applicable to the considered type of particles, i.e. sub-Kolmogorov spheres much denser than the carrier fluid. While this class has been the centre of attention in early particle-laden turbulence studies, more recently numerical and experimental works have often considered much larger and less dense particles (Qureshi et al. Reference Qureshi, Bourgoin, Baudet, Cartellier and Gagne2007; Bellani et al. Reference Bellani, Byron, Collignon, Meyer and Variano2012; Uhlmann & Doychev Reference Uhlmann and Doychev2014; Fornari et al. Reference Fornari, Picano and Brandt2016; Baker & Coletti Reference Baker and Coletti2019). Bridging the gap between the two extrema of the parameter space will be important to reach a comprehensive view of the two-way coupling problem (Brandt & Coletti Reference Brandt and Coletti2022). Experimentally, an obvious obstacle on such a path will be the difficult optical access in highly concentrated particle suspensions, especially at density ratios for which refractive index matching is not an option. Recently, Fong & Coletti (Reference Fong and Coletti2022) demonstrated that back-lighting can provide quantitative insight in cluster-induced turbulence at 
$\varPhi _{v}$ approaching 
$10^{-2}$, but did so in a relatively confined configuration and could not access the local fluid velocity. Depending on the media, non-optical techniques such as ultrasound imaging velocimetry (Poelma Reference Poelma2017) and X-ray (Aliseda & Heindel Reference Aliseda and Heindel2021) may provide attractive alternatives.
Funding
This work was funded by the US Department of Defense through the Office of Naval Research and the Army Research Office.
Declaration of interests
The authors report no conflict of interest.
Open access
