3.3.1. The spherical-shell spectral analysis

Figure 6.

Spherical-shell energy spectra for single-phase flow (black), particle-laden flows at ${\textit{Ga}}=0$ (cyan), ${\textit{Ga}}=40$ (red), ${\textit{Ga}}=80$ (green) and ${\textit{Ga}}=120$ (blue). The energy spectra are normalized by the total kinetic energy $E_k = \sum_{\kappa = 0}^{{\kappa _{\textit{max}}}}E(\kappa )$ . The vertical dashed and dotted lines depict the cutoff wavenumber $\kappa _{\!f}$ and the wavenumber corresponding to the particle scale, $\kappa _d=2\pi /D$ , respectively.

We start with the energy spectrum of the fluid flow. In figure 6 we present the energy spectra of the single-phase and particle-laden turbulent flows as functions of the spherical-shell wavenumber $\kappa = \|\boldsymbol{\kappa }\|$ . As we can see, in comparison with the single-phase flow, the energy spectra rise at small scales whereas the energy contained in the large scales is attenuated in the particle-laden turbulent flows, an observation also made in earlier studies on particle-laden turbulence (Wang et al. Reference Wang, Ayala, Gao, Andersen and Mathews2014; Bhattacharjee et al. Reference Bhattacharjee, Tom, Carbone and Bragg2024). This phenomenon of energy depletion at large scales and reinjection into small scales has been attributed to two main mechanisms (Peng et al. Reference Peng, Sun and Wang2023): (i) the promotion of energy transfer by particles due to the enhanced dissipation rate in the vicinity of the particle surfaces (Tanaka & Eaton Reference Tanaka and Eaton2010; Botto & Prosperetti Reference Botto and Prosperetti2012); (ii) the enhanced vorticity stretching caused by the particles (Peng et al. Reference Peng, Sun and Wang2023). As a consequence of the particle-induced energy transfer, a $\kappa ^{-3}$ scaling of the energy spectra is observed in the particle-laden flows, an observation widely reported in the case of bubble/particle-induced pseudoturbulence (Lance & Bataille Reference Lance and Bataille1991; Martínez Mercado et al. Reference Martínez Mercado, Chehata Gómez, Van Gils, Sun and Lohse2010; Pandey et al. Reference Pandey, Ramadugu and Perlekar2020; Innocenti et al. Reference Innocenti, Jaccod, Popinet and Chibbaro2021; Jiang et al. Reference Jiang, Brandt, Xu and Zhao2025).

Moreover, we also notice from figure 6 that the pivot wavenumber $\kappa ^\bigstar$ , the intersection point of the energy spectra of the single-phase and particle-laden flows, is close to the particle scale ( $\kappa ^\bigstar \approx \kappa _d$ ) for the case without or with weak particle sedimentation velocity (at ${\textit{Ga}}=0$ or 40), a similar observation also reported in previous studies of particle-laden turbulent flows (Wang et al. Reference Wang, Ayala, Gao, Andersen and Mathews2014; Peng et al. Reference Peng, Sun and Wang2023; Chiarini & Rosti Reference Chiarini and Rosti2024). However, it is noteworthy that the pivot wavenumber $\kappa ^\bigstar$ shifts towards lower wavenumbers with the increase of the Galileo number, indicating an enhanced particle-induced turbulence modulation as the particle sedimentation becomes faster. Accordingly, the $\kappa ^{-5/3}$ scaling of the energy spectrum, which characterizes the classic energy cascade in the inertial range of the single-phase turbulent flow, becomes hardly visible in the suspensions with ${\textit{Ga}}=80$ and 120. Note that a Kolmogorov scaling $\kappa ^{-5/3}$ of the energy of bubble-induced turbulence was identified at the sub-bubble scales in the recent experimental work by Ma et al. (Reference Ma, Tan, Ni, Hessenkemper and Bragg2025) (these authors report the equivalent $r^{2/3}$ scaling of the second-order structure function, with $r$ being the separation distance) at the high Galileo number of ${\textit{Ga}} \sim O(1000)$ . We, however, did not observe such a Kolmogorov scaling at small scales in the present study, probably due to the fact that the Galileo number under consideration is not high enough. Whether a similar scaling can be observed in the settling-particle pseudoturbulence is worth future explorations by experiments/simulations of particle-laden flows at high ${\textit{Ga}}$ . We also note that the recent numerical work by Chiarini et al. (Reference Chiarini, Tandurella and Rosti2025) reported a $\kappa ^{-4}$ scaling of the energy spectra for the HIT laden with (non-settling) Kolmogorov-size particles. The turbulence modulation by such small particles with gravity sedimentation also needs further examination.

Figure 7.

Cumulative SBS energy budget for (a) case 1 with ${\textit{Ga}}=0$ and (b) case 4 with ${\textit{Ga}}=120$ . The total viscous dissipation $\mathcal{D}^\varSigma _{\textit{tot}}=-\sum _{\kappa = 0}^{{\kappa _{\textit{max}}}} {\mathcal{D}( \kappa )}$ is used to normalize the data.

We now examine the SBS energy budget of the fluid flow in each suspension. We first compare in figure 7 the cumulative spherical-shell SBS energy budget for the flow with no particle settling and with the highest Galileo number under study. The cumulative contribution to the budget, denoted by $\mathcal{J}^\varSigma (\kappa )$ , is defined as the sum of the spherical-shell spectrum up to the given wavenumber $\kappa$ , i.e.

(3.10) \begin{align} \mathcal{J}^\varSigma (\kappa )=\sum \limits _{\kappa = 0}^{{\kappa _{\textit{max}}}} \mathcal{J}\left ( \kappa \right ). \end{align}

As shown in figure 7(a), the external EP forces inject energy into the flow system at the scales larger than that of the cutoff wavenumber $\kappa _{\!f}$ , effectively serving as an energy source; on the other hand, the energy is dissipated by the viscous stresses. The nonlinear interactions transfer energy from large to small scales with zero net energy production/dissipation in the flow (Pope Reference Pope2000), i.e. ${\varPi }^\varSigma =$ 0 at the highest wavenumber.

It is noteworthy that the contribution from the IB forces plays a similar role of transferring energy from large to small scales (indicated by the negative slope of ${\varPi }^\varSigma (\kappa )$ at low wavenumbers and positive slope at high wavenumbers). However, its contribution to the interscale energy transfer is less important than that of the nonlinear effect at small $\kappa$ , while it directly balances dissipation at large $\kappa$ . Note additionally that the particle–fluid interactions have a net contribution of energy dissipation (indicated by the negative value of $\mathcal{F}_{\textit{ib}}^\varSigma$ at the highest wavenumber $\kappa _{max}$ ) in the case with ${\textit{Ga}}=0$ , a fact which deserves an explanation.

Let us therefore discuss the non-zero energy dissipation due to the particle–fluid interactions. To this end, we resort to the budget of the total fluid kinetic energy $E_k=(u_x^2+u_y^2+u_z^2)/2$ in the physical space (see also Chouippe & Uhlmann Reference Chouippe and Uhlmann2015):

(3.11) \begin{equation} \frac {{d}}{{\mathrm{d}t}}{\left \langle {{E_k}} \right \rangle _\varOmega } = - {\langle \epsilon \rangle _\varOmega } + {\left \langle {{{\boldsymbol{f}}_{\!\textit{ext}}} \boldsymbol{\cdot }{\boldsymbol{u}}} \right \rangle _\varOmega } + {\left \langle {{{\boldsymbol{f}}_{\!\textit{ib}}} \boldsymbol{\cdot }{\boldsymbol{u}}} \right \rangle _\varOmega }. \end{equation}

Here, $\langle \boldsymbol{\cdot }\rangle _\varOmega$ represents the spatial average. The terms on the right-hand side of (3.11) represent the viscous dissipation ( $- {\langle \epsilon \rangle _\varOmega }$ ), the energy input by external forces (denoted by $\varPsi _{\textit{ext}}=\left \langle {{{\boldsymbol{f}}_{\!\textit{ext}}} \boldsymbol{\cdot }{\boldsymbol{u}}} \right \rangle _\varOmega$ ) and the contribution from the particle–fluid interactions (denoted by $\varPsi _{p}=\left \langle {{{\boldsymbol{f}}_{\!\textit{ib}}} \boldsymbol{\cdot }{\boldsymbol{u}}} \right \rangle _\varOmega$ ). These terms are related to the energy budget terms in the spectral space by

(3.12) \begin{align} \langle \epsilon \rangle _\varOmega =-\sum \limits _{\boldsymbol{\kappa }} \mathcal{D}(\boldsymbol{\kappa }),\quad \varPsi _{\textit{ext}}=\left \langle {{\boldsymbol{f}}_{\!\textit{ext}}} \boldsymbol{\cdot }{\boldsymbol{u}} \right \rangle _\varOmega = \sum \limits _{\boldsymbol{\kappa }} \mathcal{F}_{\textit{ext}}(\boldsymbol{\kappa }),\quad \varPsi _{p}=\left \langle {{\boldsymbol{f}}_{\!\textit{ib}}} \boldsymbol{\cdot }{\boldsymbol{u}} \right \rangle _\varOmega =\sum \limits _{\boldsymbol{\kappa }} \mathcal{F}_{\textit{ib}}(\boldsymbol{\kappa }). \end{align}

Considering the particle energy budget, the term $\varPsi _p$ can be related to the particle motions as follows:

(3.13) \begin{equation} \begin{aligned} {\varPsi _p} & \approx \underbrace { - \frac {\phi }{{\rho _{\!p}}{V_{\!p}}} \left ( {\frac {{{\rho _{\!p}}}}{{{\rho _{\!f}}}} - 1} \!\right )\!\left [ {{ \rho _{\!p}V_{\!p} \frac {{d}}{{\textrm {d}t}}{{\left \langle {\frac {1}{2}{\boldsymbol{v}} \boldsymbol{\cdot }{\boldsymbol{v}}} \right \rangle }_p}} + {I_{\!p}}\frac {{d}}{{\mathrm{d}t}}{\left \langle {\frac {1}{2}{{\boldsymbol{\omega }}} \boldsymbol{\cdot }{{\boldsymbol{\omega }}}} \right \rangle _p}} \right ]}_{\varPsi _p^{acc}} + \underbrace {\phi \left ( {\frac {{{\rho _{\!p}}}}{{{\rho _{\!f}}}} - 1} \right ){{\left \langle {\boldsymbol{v}} \right \rangle }_p} \boldsymbol{\cdot }{\boldsymbol{g}}}_{\varPsi _p^{\textit{gra}}} \\ &\quad + \underbrace {\frac {\phi }{{{\rho _{\!f}}{V_{\!p}}}}\left [ {{{\left \langle {{\boldsymbol{v}} \boldsymbol{\cdot }{{\boldsymbol{F}}_{\!C}}} \right \rangle }_p} + {{\left \langle {{\boldsymbol{\omega }} \boldsymbol{\cdot }{{\boldsymbol{T}}_C}} \right \rangle }_p}} \right ]}_{\varPsi _p^{col}} \underbrace { - \frac {\phi }{{{V_{\!p}}}}{{\left \langle {{\boldsymbol{v}} \boldsymbol{\cdot }\int _{{\varOmega _p}} {{{\boldsymbol{f}}_{\!\textit{ext}}}{\rm d}V} + {\boldsymbol{\omega }} \boldsymbol{\cdot }\int _{{\varOmega _p}} {{\boldsymbol{r}} \times {{\boldsymbol{f}}_{\!\textit{ext}}}{\rm d}V} } \right \rangle }_p}}_{\varPsi _p^{\textit{ext}}}. \end{aligned} \end{equation}

The derivation of the above equation, which was originally performed in Chouippe & Uhlmann (Reference Chouippe and Uhlmann2015), is provided in Appendix A for completeness. According to (3.13), the total work related to the particle–fluid interactions $\varPsi _p$ is composed of the contributions from the particle acceleration $\varPsi _p^{acc}$ , the particle sedimentation $\varPsi _p^{\textit{gra}}$ , interparticle collisions $\varPsi _p^{col}$ and the term related to the external forces $\varPsi _p^{\textit{ext}}$ . Note that the term inside the square bracket of $\varPsi _p^{acc}$ can be reformulated as

(3.14) \begin{align} {{ \rho _{\!p}V_{\!p} \frac { {d}}{{\textrm {d}t}}{{\left \langle {\frac {1}{2}{\boldsymbol{v}} \boldsymbol{\cdot }{\boldsymbol{v}}} \right \rangle }_p}} + {I_{\!p}}\frac {{d}}{{\mathrm{d}t}}{\left \langle {\frac {1}{2}{{\boldsymbol{\omega }}_p} \boldsymbol{\cdot }{{\boldsymbol{\omega }}_p}} \right \rangle _p}} = \frac {1}{2}\frac {{d}}{{\textrm {d}t}}{\left \langle {{\rho _{\!p}}{V_{\!p}}{{\boldsymbol{v}}^2} + {I_{\!p}}{{\boldsymbol{\omega }}^2}} \right \rangle _p}, \end{align}

namely the time derivative of the particle-averaged kinetic energy. Hence, $\varPsi _p^{acc}$ can only have a non-zero contribution to the energy budget for a non-equilibrium flow; for example, it works either as a energy source or sink, depending on the particle size, in particle-laden decaying turbulence (Lucci et al. Reference Lucci, Ferrante and Elghobashi2010). In the steady state of sustained turbulence, however, $\varPsi _p^{acc}$ is equal to zero. Additionally, $\varPsi _p^{\textit{gra}}$ is equal to zero when the particles do not settle ( ${\textit{Ga}}=0$ ). As regards the other two terms in (3.13), the contribution from particle–particle collisions was reported to be negligible at the volume fraction $\phi \sim O(0.1)$ (Chiarini & Rosti Reference Chiarini and Rosti2024), which is also confirmed in the present simulations. Hence, from our simulation results, the last term $\varPsi _p^{\textit{ext}}$ is mainly responsible for the negative value of $\varPsi _p$ in the case of ${\textit{Ga}}=0$ , corresponding to a weak but non-negligible net energy dissipation associated with the particle–fluid interactions (see figure 7 a).

As shown in Appendix B, the net dissipation originated from the IB forces at ${\textit{Ga}}=0$ is related to the non-zero EP external forces applied also inside the rigid particles. However, as shown in the same appendix, this contribution becomes almost negligible when the external forcing is excluded from the solid region; then the residual contribution is merely from interparticle collisions ( $\varPsi _p^{col}$ ). Note, also, from the definition of $\varPsi _p^{\textit{ext}}$ that a negative value of this term indicates a positive correlation between the particle velocities and the external forces sustaining the turbulence, which is attributed to the tendency of neutrally buoyant particles to follow local fluid flows produced by the external forces. Although it might have been preferable to force the turbulence only in the fluid phase, we show that this does not change the physics discussed in this work, especially above the wavenumbers of the external forcing, which acts at the large scale of the computational domain.

Now let us move back to the discussion on figure 7(b), which provides the cumulative SBS budget of the fluid kinetic energy in the suspension with ${\textit{Ga}}=120$ . Compared with the non-sedimentation case ( ${\textit{Ga}}=0$ in figure 7 a), it is obvious that the total work done by the IB forces is significantly augmented. The positive value of $\mathcal{F}_{\textit{ib}}^\varSigma (\kappa _{max})$ indicates a net energy input by the particle–fluid interactions at ${\textit{Ga}}=120$ . This is because the contribution from the particle sedimentation, i.e. $\varPsi _p^{\textit{gra}}$ in (3.13), does positive work in the flow system, which overwhelms the energy dissipation caused by collisions and external forces inside the particles (negative values of $\varPsi _p^{col}$ and $\varPsi _p^{\textit{ext}}$ ); the work done by the settling particles, instead of the energy production by the external EP forces, serves as the major source of the energy input. As a consequence, to close the energy budget, the viscous dissipation becomes greater than $\varPsi _{\textit{ext}}$ , different from the physics in the single-phase turbulence.

To have a clearer comparison between different cases, we show the cumulative contributions to the energy budget in figure 8: in particular, we display the nonlinear transfer ${\varPi }^\varSigma$ , viscous dissipation $\mathcal{D}^\varSigma$ and the particle–fluid interactions $\mathcal{F}_{\textit{ib}}^\varSigma$ at different ${\textit{Ga}}$ . These comparisons illustrate the modulation of the energy budget in the particle-laden turbulence with varying particle settling velocity. We first observe from figure 8(a) that the total energy input by the IB forces, i.e. the value of $\mathcal{F}_{\textit{ib}}^\varSigma$ at the highest wavenumber, increases monotonically with the Galileo number. This enables us to provide an estimate of the ratio between the energy input by the gravity sedimentation and the external forces that sustain the single-phase turbulence:

(3.15) \begin{equation} \varLambda = \frac {{{\varPsi _p^{\textit{gra}}}}}{{{\varPsi _{\textit{ext}}}}} \approx \frac {{{\varPsi _p^{\textit{gra}}}}}{{{\epsilon _0}}} = \frac {{\phi \left ( {{\rho _{\!p}}/{\rho _{\!f}} - 1} \right )g\left \langle {{V_s}} \right \rangle }}{{{\epsilon _0}}} = \frac {{\phi G{a^2}{{\textit{Re}}_s}}}{{5832\epsilon _0^*}}. \end{equation}

Here, ${\textit{Re}}_s=\langle V_s \rangle D/\nu$ is the Reynolds number based on the particle mean settling velocity $\langle V_s \rangle$ , and $\epsilon _0^*=\epsilon \tau _p^3/D^2$ is the dimensionless viscous dissipation of the single-phase turbulence provided in table 2. Note that we use the single-phase viscous dissipation $\epsilon _0$ to estimate the energy production by the external forces $\varPsi _{\textit{ext}}$ , which is only strictly valid for the single-phase flow. However, the difference between $\epsilon _0$ and $\varPsi _{\textit{ext}}$ is also found to be small (up to approximately $10\,\%$ ) in the particle-laden cases under consideration. As the dimension of the particle is an order of magnitude smaller than the energy-input scale (i.e. $\kappa _d/\kappa _{\!f}=12.8 \gg 1$ ), the energy input by the large-scale external forces is not significantly affected by the presence of particles.

In table 4 we list the value of $\varLambda$ computed by (3.15) with the mean particle settling velocity obtained from the simulations. Not surprisingly, $\varLambda$ increases with the Galileo number because of the increased particle settling velocity. Also, it is noteworthy that the ratio is as low as $\varLambda =0.07$ for the case with ${\textit{Ga}}=40$ , indicating the negligible energy input by gravity sedimentation compared with the external forces. This is in agreement with the marginal difference regarding the SBS budgets between the cases with ${\textit{Ga}}=0$ and ${\textit{Ga}}=40$ (the red and cyan lines in figure 8). Nevertheless, increasing the Galileo number to ${\textit{Ga}}=80$ , the energy injection due to the particle sedimentation becomes comparable to that of the external forces with $\varLambda =0.94$ ; the former contribution becomes the predominant mechanism of energy production at ${\textit{Ga}}=120$ when the contribution from the particle sedimentation exceeds the EP forces by a factor approximately 4 ( $\varLambda =3.92$ ).

Table 4.

Ratio $\varLambda$ between the energy input by the particle settling and by the external forces maintaining turbulence for the different cases under examination. See text for the definition of $ \varLambda$ in (3.15).

Figure 8.

Comparison of the cumulative contributions to the SBS energy budget for the different Galileo number under consideration: (a) IBM forcing term; (b) nonlinear transfer; (c) viscous dissipation. Here we only present the data for $\kappa \gt \kappa _{\!f}$ .

We also observe from figure 8(a) that the large-scale energy transfer by the particle–fluid interactions, indicated by the negative slope of $\mathcal{F}_{\textit{ib}}^\varSigma (\kappa )$ in the suspensions with ${\textit{Ga}}=0$ and 40, becomes negligible at higher Galileo numbers. As discussed below, the decrease with ${\textit{Ga}}$ of the net energy removal by the settling particles at large scales is associated with a cancellation between the additional injection of energy into vertically elongated scales and the quenching of the horizontal flow correlations, so at the origin of the flow anisotropy. Furthermore, the range of scales on which the energy is injected by the particle–fluid interactions (region with positive slope of $\mathcal{F}_{\textit{ib}}^\varSigma (\kappa )$ ) is enlarged when increasing the Galileo number.

The different behaviour of the particle–fluid interactions also alters the energy transport by the nonlinear terms. As shown in figure 8(b), this energy transfer is weakened as the Galileo number increases. For the greatest particle settling velocity at ${\textit{Ga}}=120$ , the contribution from the nonlinear transfer to the energy budget becomes subdominant compared with those from the particle–fluid interaction and the viscous dissipation (see figure 7 b). In other words, the energy input by the particle sedimentation is almost in direct balance with the viscous dissipation within a wide wavenumber range, i.e. $\mathcal{F}_{\textit{ib}} (\kappa ) \approx -\mathcal{D}_{\textit{ib}} (\kappa )$ , at ${\textit{Ga}}=120$ . This observation supports the local-balance assumption proposed by Lance & Bataille (Reference Lance and Bataille1991) regarding the pseudoturbulence induced by rising bubbles. In our recent study on the particle sedimentation in an initially quiescent fluid (Jiang et al. Reference Jiang, Brandt, Xu and Zhao2025), we also observed this local balance in the particle-induced pseudoturbulence, which leads to the $\kappa ^{-3}$ scaling of the energy spectra. Note, however, that this local balance is only observed at the smallest scales at lower ${\textit{Ga}}$ .

As regards the spherical-shell spectra of the viscous dissipation, we observe from figure 8(c) that $\mathcal{D}^\varSigma (\kappa ) \sim -\log (\kappa )$ on the scales close to the particle size, corresponding to the $E(\kappa ) \sim \mathcal{D}(\kappa )/\kappa ^2 \sim \kappa ^{-3}$ scaling of the energy spectra for the particle-laden flows shown in figure 6. Note, interestingly, that this scaling is not clearly observed in the particle-laden turbulence with ${\textit{Ga}}\leqslant 40$ (in absence of or with weak gravity sedimentation), when the local-balance mechanism only holds at the smallest scales.

To gain more insights on the small-scale dynamics, we further examine the energy budget at scales smaller than the particle size. Specifically, we compute the integrals of the nonlinear transfer term, IB forcing term and the viscous dissipation term at $\kappa \geqslant \kappa _d$ , denoted by

(3.16) \begin{align} {\varPi }^\lt =\sum \limits _{\kappa \geqslant \kappa _d} {\varPi }\left ( \kappa \right ), \quad \mathcal{F}^\lt _{\textit{ib}}=\sum \limits _{\kappa \geqslant \kappa _d} \mathcal{F}_{\textit{ib}}\left ( \kappa \right ), \quad \mathcal{D}^\lt =\sum \limits _{\kappa \geqslant \kappa _d} \mathcal{D}\left ( \kappa \right ). \end{align}

The physical interpretations of these three quantities are: (i) ${\varPi }^\lt$ the total nonlinear energy transfer by inertia from larger scales to subparticle scales; (ii) $\mathcal{F}^\lt _{\textit{ib}}$ the total energy contribution from the particle–fluid interactions at subparticle scales; (iii) $\mathcal{D}^\lt$ the total energy dissipation at subparticle scales. The values of ${\varPi }^\lt _{\textit{ib}}$ , $\mathcal{F}^\lt _{\textit{ib}}$ and $\mathcal{D}^\lt$ are shown in figure 9 for all cases being examined, normalized by the single-phase dissipation rate $\epsilon _0$ . Note first that the balance of the energy budget suggests $\mathcal{D}^\lt = -({\varPi }^\lt + \mathcal{F}^\lt _{\textit{ib}})$ , which is confirmed by the data shown in figure 9. In contrast to the single-phase turbulence with only the nonlinear transfer as energy source, the particle–fluid interaction always provides a net energy input at subparticle scales in the particle-laden cases. In other words, the IB forces do positive work at small scales, regardless of the Galileo number. As the Galileo number increases, the energy input by particle sedimentation is enhanced due to the greater particle settling velocity; the energy brought to small scales by nonlinear effects decreases in particle-laden flows with ${\textit{Ga}}\leqslant 40$ and then increases, remaining, however, lower than the input by the sedimentation at high ${\textit{Ga}}\geqslant 80$ . It is also interesting to note that the absolute value of ${\varPi }^\lt$ is larger than that of the single-phase turbulence at ${\textit{Ga}}=120$ . This indicates that the energy injection by the particle–fluid interaction, which starts at scales larger than the particle size (see the positive slope of $\mathcal{F}^\varSigma _{\textit{ib}}$ at $\kappa \lesssim \kappa _d$ in figure 8 a), would be partially transferred to small scales by inertial effects. The results shown in figures 8 and 9 suggest that more energy is available at $\kappa = \kappa _d$ at high ${\textit{Ga}}$ , which is sufficient to trigger nonlinear velocity fluctuations and to force an extension of the dissipative range towards smaller scales.

Figure 9.

The energy budget at small scales (smaller than the particle size, $\kappa \geqslant \kappa _d$ ) for the different suspension flows under investigation. The numbers on the horizontal axis represent the Galileo number, except for ‘SP’ standing for the single-phase turbulence.

The analysis of the SBS budget indicates how the combined effects of the particle-enhanced energy cascade (evident at ${\textit{Ga}}=0$ or 40) and the energy input by particle sedimentation (more evident at ${\textit{Ga}}=80$ or 120) alter the energy transfer path in the particle-laden turbulent flows. In particular, the energy transfer induced by the particle–fluid interactions quenches the classic inertia-driven energy cascade (characterized by $E(\kappa ) \sim \kappa ^{-5/3}$ in the inertial range), leading to a robust $\kappa ^{-3}$ scaling of the energy spectra. When settling is null or relatively weak, the local balance between fluid–structure interactions and dissipation is seen only at small scales. As soon as the settling speed (or equivalently the Galileo number) increases, the $\kappa ^{-3}$ scaling of the energy spectra appears on wider scale ranges, even at the largest flow scales, completely masking the classic $\kappa ^{-5/3}$ spectra for the highest ${\textit{Ga}}$ under the current investigation. Indeed, we notice the coexistence of the $\kappa ^{-5/3}$ and $\kappa ^{-3}$ scalings of the energy spectra at different scales in the suspensions with ${\textit{Ga}}=0$ and ${\textit{Ga}}=40$ (see the cyan and red lines in figure 6); in these flows, the modulation of the energy transfer by the particle–fluid interactions does not reach the inertial scales and remains limited to structures smaller than the particle size. A similar observation regarding the coexistence of inertia-dominated and fluid–structure-dominated dynamics underlying the energy spectrum was also reported in bubble-laden and particle-laden turbulent flows (Pandey, Mitra & Perlekar Reference Pandey, Mitra and Perlekar2023; Chiarini & Rosti Reference Chiarini and Rosti2024).

3.3.2. The spectral analysis on the two-dimensional plane

Figure 10.

Energy spectral density in the $\kappa _h-\kappa _y$ plane. The notation $\tilde {E}^*$ denotes the normalized energy spectral density, i.e. $\tilde {E}^*=\tilde {E}/(E_k D^3)$ . The dashed grey contour lines depict the result for the particle-laden turbulence without particle sedimentation ( ${\textit{Ga}}=0$ ), while the coloured solid lines illustrate the result for (a,e) the single-phase flow, and the particle-laden flows with (b, f) ${\textit{Ga}}=40$ , (c,g) ${\textit{Ga}}=80$ and (d,h) ${\textit{Ga}}=120$ . The linear scale is used for the axes in (a–d) while the log scale is used in (e–h).

To further understand the directional information of the turbulence modulated by the settling particles, we analyse the energy budget on the two-dimensional $\kappa _h-\kappa _y$ plane, corresponding to horizontal and vertical wavenumbers. First, we present in figure 10 the two-dimensional distribution of the energy spectral density. Here, we display the data with linear scale in figure 10(a–d) to better illustrate the isotropy/anisotropy and with log scale on the bottom row to clarify the large-scale information. In each panel, we display with dashed lines the results pertaining to the suspension with ${\textit{Ga}}=0$ .

Figure 10(a,e) compare the single-phase turbulence with the particle-laden flow of neutrally buoyant particles. As expected, the energy distribution retains isotropy in the absence of particle sedimentation, as shown by the concentric-circle contour lines of $\tilde {E}^*$ with the linear scale of $\kappa _h$ and $\kappa _y$ (figure 10 a), or the diagonal-symmetric contour lines with the log-scale axes (figure 10 e). The presence of particles induces a slower decay of $\tilde {E}^*$ along the radial direction, i.e. an increase of energy at the smallest scales, which is in agreement with the raised spherical-shell energy spectrum below the particle scale shown in figure 6.

For the cases with non-zero Galileo numbers, we observe from figure 10(b–d, f–h) that the particle sedimentation breaks the isotropy of the energy distribution in the spectral space. The contour lines of $\tilde {E}^*$ are elongated along the $\kappa _h$ -direction on the scales close to or moderately smaller than the particle size ( $\kappa _h \gtrsim \kappa _d$ ), an observation becoming more prominent as the Galileo number increases. Specifically, for a given solid contour line at intermediate scales in figures 10(c) or 10(d), which represents flow structures with equal energy, we see that its intercept on the horizontal axis occurs at larger wavenumber than that on the vertical axis. This indicates the larger physical dimension of this flow structure along the vertical direction than the horizontal, and can be related to the local flow structures around the settling particles: thin shear layers and elongated wakes appear in the lateral and vertical directions, respectively (see also Fornari et al. Reference Fornari, Picano and Brandt2016a ). Note, however, that the contour lines of $\tilde {E}^*$ tend to recover symmetry far from the origin of the $\kappa _h-\kappa _y$ plane, indicating the recovery of the flow isotropy at small scales. As concerns the suspension at ${\textit{Ga}}=40$ , however, the anisotropy is hardly visible (figure 10 b, f). This can be ascribed to the marginal effect of particle sedimentation as the particle mean settling velocity is significantly smaller than the turbulent velocity fluctuations ( $\langle V_s \rangle / u'=0.21$ ) in this case.

Figure 11.

Spectra of the IB forcing work in the $\kappa _h-\kappa _y$ plane for the particle-laden flows under investigation. The data is normalized by $\tilde {\mathcal{F}}_{\textit{ib}}^*=\tilde {\mathcal{F}}_{\textit{ib}}/(\mathcal{D}_{\textit{tot}}^\varSigma D^3)$ , and only positive values of $\tilde {\mathcal{F}}_{\textit{ib}}^*$ are illustrated by contour lines (denoted by the subscript ‘+’) . The dashed grey lines correspond to non-settling particles ( ${\textit{Ga}}=0$ ), while the coloured solid lines illustrate the results for (a) ${\textit{Ga}}=40$ , (b) ${\textit{Ga}}=80$ and (c) ${\textit{Ga}}=120$ . Log-scale axes are used. The symbols and indicates the region with negative values of $\tilde {\mathcal{F}}_{\textit{ib}}^*$ and the position of the positive peak of $\tilde {\mathcal{F}}_{\textit{ib}}^*$ .

Figure 12.

Same as figure 11 but for the two-dimensional spectra of the nonlinear term $\tilde {{\varPi }}$ .

Next, we present the two-dimensional contours of the spectral density associated with the IB forcing in figure 11 and with the nonlinear term in figure 12. In each panel, the results pertaining to a non-zero settling velocity are compared with the case of ${\textit{Ga}}=0$ . As shown in figure 11(a) with the results at ${\textit{Ga}}=0$ and 40, the particle–fluid interactions mainly transfer energy from large to small scales along the radial direction in the spectral space so the flow anisotropy is hardly visible; in particular, the directional bias of $\tilde {\mathcal{F}}_{\textit{ib}}$ is null for the case with ${\textit{Ga}}=0$ and very weak for ${\textit{Ga}}=40$ . However, when the effect of particle sedimentation becomes significant with the increase of ${\textit{Ga}}$ to 80 and 120, we observe a clear positive peak of $\tilde {\mathcal{F}}_{\textit{ib}}$ around $\kappa _h \approx 0.3-0.6\kappa _d$ for low values of $\kappa _y$ , indicated by the symbol in figure 11(b,c), a signature of the energy injection by the settling particles. We therefore demonstrate that the energy input by the gravitational sedimentation is anisotropic in the spectral space, unlike the isotropic large-scale energy injection by the external EP forces. This mechanism is responsible for the anisotropic energy distribution shown in figure 10(c,d,g,h). Note also that, the region of the spectral space characterized by negative values of $\tilde {\mathcal{F}}_{\textit{ib}}$ (marked by in figure 11), indicating energy removal by the particle–fluid interactions, loses its isotropy at large scales as ${\textit{Ga}}$ increases. This is associated with the aforementioned positive peak of $\tilde {\mathcal{F}}_{\textit{ib}}$ as well. The increased energy injection by particle sedimentation partially cancels out the particle-induced energy transfer from large to small scales, which leads to a decreased energy removal by the IB forces at low wavenumbers (see figure 8 a). Additionally, the anisotropic energy input by the particle sedimentation becomes more pronounced for higher ${\textit{Ga}}$ , along with the positive value of $\tilde {\mathcal{F}}_{\textit{ib}}$ expanding to lower wavenumbers. To sum up, the settling particles induce flow anisotropy by forcing vertically elongated structures while breaking the flow coherence in the horizontal direction; this only leads to the weakened energy injection/removal at large scales when examining the spherical-shell energy spectra, yet leaves an imprinting on the flow (its anisotropy).

To conclude this section, let us elaborate on the two-dimensional spectrum of the energy transfer by nonlinear terms, with the results shown in figure 12. Similar to figure 10(a,e) and figure 11(a), the two-dimensional plot of $\tilde {{\varPi }}$ retains isotropy on the $\kappa _h-\kappa _y$ plane for the particle-laden flows with ${\textit{Ga}}=0$ and ${\textit{Ga}}=40$ (figure 12 a). Increasing the particle settling speed, $\tilde {{\varPi }}$ exhibits anisotropy in the cases at ${\textit{Ga}}=80$ and 120 (figures 12 b and 12 c), which can be seen as a consequence of the anisotropic energy input by the particle–fluid interactions (the IB forcing term).

It is interesting to note the negative peak of $\tilde {{\varPi }}$ on the horizontal axis at $\kappa _h \approx 0.5\kappa _d$ for the ${\textit{Ga}}=120$ case, highlighted by the symbol in figure 12(c). This indicates the removal of energy by the nonlinear effect at these intermediate scale, as also indicated by the negative slope of ${\varPi }^\varSigma (\kappa )$ for the associated case at $\kappa \approx 0.3\kappa _d-0.5\kappa _d$ in figure 8(b). Notice also that the location of the negative peak of $\tilde {{\varPi }}$ in figure 12(c) coincides with the positive peak of $\tilde {\mathcal{F}}_{\textit{ib}}$ shown in figure 11(c). Therefore, for the flow with a substantial energy supply from the particle sedimentation, part of the energy input by the particle–fluid interaction is transferred outwards at intermediate scales by the nonlinear effects. This is analogy to the energy extraction by inertial effects at large scales for the single-phase turbulence.

3.4.1. Voronoi analysis

In this section, we examine the spatial clustering and collision rate of the dispersed particles in the turbulent flows. To begin with, we employ the Voronoi analysis to investigate the particle spatial distribution (Monchaux, Bourgoin & Cartellier Reference Monchaux, Bourgoin and Cartellier2010). With this analysis tool, the domain is partitioned into Voronoi tessellations of the same number of particles, such that any point in one tessellation is closest to the associated particle. In figure 13(a) we present the probability distribution of the volumes of the Voronoi tessellations in each particle-laden case under study. The symbol $\overline {V}_{Voro}$ denotes the Voronoi volume ${V}_{Voro}$ normalized by the mean value for the same particle concentration, $\langle {V}_{Voro} \rangle = V_{\varOmega }/N_{\!p}$ . Note that we also show the data for the randomly distributed non-overlapping particles at the same volume fraction $\phi =10\,\%$ for the sake of comparison. The statistics for the randomly distributed particles are collected by repeatedly generating non-overlapping particles with a Poisson process (Uhlmann Reference Uhlmann2020). The data in the figure reveal that the PDF of $\overline {V}_{Voro}$ is more intermittent in HIT than for the random spatial distribution, indicating the occurrence of particle clustering (Monchaux et al. Reference Monchaux, Bourgoin and Cartellier2010).

To quantify the degree of particle clustering, we compute the standard deviation of the normalized Voronoi volumes, $\sigma (\overline {V}_{Voro})$ , and define the clustering indicator $C$ as the ratio $\sigma (\overline {V}_{Voro})$ to that of the randomly distributed particles (Tagawa et al. Reference Tagawa, Roghair, Prakash, Van Sint Annaland, Kuipers, Sun and Lohse2013):

(3.17) \begin{equation} C=\sigma (\overline {V}_{Voro})/\sigma (\overline {V}_{Voro}^{rand}). \end{equation}

As shown in figure 13(b), the value of $C$ is greater than unity in each case but does not change much with the Galileo number, which is consistent with the PDFs shown in figure 13(a). In previous studies, heavy, finite-size particles (without sedimentation) were reported to exhibit clustering in HIT (Fiabane et al. Reference Fiabane, Zimmermann, Volk, Pinton and Bourgoin2012; Uhlmann & Chouippe Reference Uhlmann and Chouippe2017; Chiarini & Rosti Reference Chiarini and Rosti2024), unlike neutrally buoyant particles (Fiabane et al. Reference Fiabane, Zimmermann, Volk, Pinton and Bourgoin2012). Nevertheless, particles settling in an otherwise quiescent fluid (without background turbulence) were not observed to cluster in both dilute and non-dilute suspensions at ${\textit{Ga}}=80$ (Jiang et al. Reference Jiang, Brandt, Xu and Zhao2025, Reference Jiang, Zhao and Brandt2026). Therefore, turbulence transport, instead of particle sedimentation, plays an essential role in the clustering of finite-size particles; the so-called ‘sweep-stick’ mechanism (Goto & Vassilicos Reference Goto and Vassilicos2008), i.e. particles preferentially accumulate in weak-fluid-acceleration regions (sticky points), is considered responsible for the particle clustering in most scenarios (Uhlmann & Chouippe Reference Uhlmann and Chouippe2017; Chiarini & Rosti Reference Chiarini and Rosti2024). This may explain the weak variation of the particle clustering intensity with the particle settling velocity, as shown by the data presented in figure 13(b). Moreover, in figure 13(b) we also quantitatively compare the clustering indicator from the present data with earlier studies. The discrepancy suggests the sophisticated dependence of the particle clustering on the particle properties, in agreement with the conclusion in Chiarini & Rosti (Reference Chiarini and Rosti2024). In particular, the effect of particle sedimentation on the clustering of finite-size particles in HIT remains to be further explored. Note, additionally, that the clustering was found to contribute to the enhancement of particle settling (or bubble rising) velocity in dilute scenarios (Jiang et al. Reference Jiang, Xu and Zhao2024b ; Hessenkemper et al. Reference Hessenkemper, Bragg, Lucas and Ma2025). This effect appears to be negligible in the present set-up given the relatively high-volume fraction; the particle mean settling velocity is always lower than that of an isolated particle as discussed in § 3.2.

Figure 13.

(a) Probability density function (PDF) of the normalized Voronoi volumes of the dispersed particles. (b) The clustering indicator $C$ versus the Galileo number ${\textit{Ga}}$ . The present results (indicated by $\bullet$ ) are compared with previous studies, ( $\blacksquare$ ) heavy particles in HIT ( $\rho _{\!p}/\rho _{\!f}=1.5$ , $D/\eta = 11$ , $\phi =0.5\,\%$ ) (Uhlmann & Chouippe Reference Uhlmann and Chouippe2017); ( $\blacktriangle$ ) heavy particles in HIT ( $\rho _{\!p}/\rho _{\!f}=1.29$ , $D/\eta = 16$ , $\phi =7.9\,\%$ ) (Chiarini & Rosti Reference Chiarini and Rosti2024); ( $\blacklozenge$ ) settling particles in the otherwise quiescent fluid ( $\rho _{\!p}/\rho _{\!f}=2$ , $\phi =5\,\%$ ) (Jiang et al. Reference Jiang, Zhao and Brandt2026).

3.4.2. Collision kernel and the particle pair statistics

Another issue we are concerned with is the collision rate of the dispersed particles, which is relevant to a variety of physical processes in environmental and industrial applications (Pumir & Wilkinson Reference Pumir and Wilkinson2016). In a variety of existing studies, the collision kernel $\varGamma$ is employed to quantify the collision rate in a particulate system (Smoluchowski Reference Smoluchowski1917; Saffman & Turner Reference Saffman and Turner1956; Pumir & Wilkinson Reference Pumir and Wilkinson2016; Wang Reference Wang, Subramaniam and Balachandar2023). For the monodisperse system, the collision kernel is defined as the rate of collision events normalized by the square of the particle number density (Wang, Wexler & Zhou Reference Wang, Wexler and Zhou2000),

(3.18) \begin{equation} \varGamma = 2\dot {N}_C/n^2. \end{equation}

Here, $\dot {N}_C$ denotes the number of collision events per unit volume and per unit time, and $n=N_{\!p}/V_{\varOmega }$ is the particle number density. To obtain $\varGamma$ , we count the number of collision events within the statistically steady state during the simulations and then compute the collision kernel by its definition (3.18). As shown in figure 14, the collision kernel $\varGamma$ , normalized by the single-phase r.m.s. velocity $u'_0$ , first decreases slightly with the introduction of particle sedimentation at ${\textit{Ga}}=40$ , and then increases when further increasing ${\textit{Ga}}$ .

To better examine the combined effect of sedimentation and turbulent transport on the collision rate, we compare the present results with the particle collision kernel in the absence of background turbulence. In our previous study on the particle sedimentation in an initially quiescent fluid (Jiang et al. Reference Jiang, Zhao and Brandt2026), the collision kernel, normalized by the single-particle settling velocity, was found to be $\varGamma /(V_{s0}D^2)=0.29$ at ${\textit{Ga}}=80$ – a value obtained at a volume fraction $\phi =5\,\%$ but found to be insensitive to the particle volume fraction. If we similarly normalize the collision kernels by $V_{s0}$ for the present cases, we can obtain the results of $\varGamma /(V_{s0}D^2)$ as listed in table 5. For the same ${\textit{Ga}}=80$ the collision kernel is increased by approximately 4.75 times in the present turbulent flow when compared with that in the quiescent fluid (Jiang et al. Reference Jiang, Zhao and Brandt2026). Additionally, the value of $\varGamma /(V_{s0}D^2)$ decreases considerably as ${\textit{Ga}}$ increases due to the increase of $V_{s0}$ , which is in contrast to the weak variation of $\varGamma /(u'_0 D^2)$ with ${\textit{Ga}}$ . Thus, we can conclude that the turbulent transport is the dominant effect for particle collisions in the present flow configuration.

Table 5.

Collision kernel normalized by the single-particle settling velocity.

Figure 14.

Collision kernel normalized by the r.m.s. velocity fluctuation of the single-phase turbulence as a function of ${\textit{Ga}}$ .

Finally, we attempt to interpret the variations in the collision kernel with the Galileo number shown in figure 14. To this end, we introduce the kinematic description of the collision kernel (Wang et al. Reference Wang, Wexler and Zhou2000; Wang Reference Wang, Subramaniam and Balachandar2023),

(3.19) \begin{equation} \varGamma ^K=2\pi R^2 \langle |W_r(R)| \rangle g (R). \end{equation}

Here, $R=D$ is the collision radius, and $\langle |W_r(R)| \rangle$ and $g(R)$ represent the mean absolute RRV and the RDF for particles at the separation distance $r=R$ . Readers can refer to Wang et al. (Reference Wang, Wexler and Zhou2000) for the definitions of RRV and RDF presented in (3.19). To provide a physical interpretation, the kinematic collision $\varGamma ^K$ is the flux of particles across a collision sphere centred at the reference particle with a collision radius $R=D$ . The RRV provides the information of the relative motion among particles, whereas the RDF reflects the particle spatial distribution by measuring the probability of finding particles pairs with a given separation distance. Theoretically, the kinematic collision kernel $\varGamma ^K$ is equivalent to the dynamic collision kernel $\varGamma$ defined in (3.18) (Wang et al. Reference Wang, Wexler and Zhou2000).

Figure 15.

(a) Radial distribution function (RDF) and (b) relative radial velocity (RRV) of the dispersed particles as function of the separation distance $r$ . The inset panels are the zoom-in view close to $r=D$ . The RDF of randomly distributed non-overlapping particles (the hard-sphere distribution, see Liu (Reference Liu2024)) is displayed in panel (a) with dashed line for comparison.

Figure 16.

(a,b) The particle pair distribution function and (c,d) relative particle velocity at (a,c) ${\textit{Ga}}=0$ and (b,d) ${\textit{Ga}}=120$ as function of the distance $r$ and the angle $\psi$ .

In figure 15, we display the RDF and RRV as functions of the particle separation distance $r$ . As shown in figure 15(a), the RDF exhibits a pronounced peak at $r=D$ ; this peak indicates the enhanced probability of finding touching particle pairs in the turbulent flow compared with randomly distributed particles. As the separation distance $r$ increases, the RDF quickly decreases to $g(r)=1$ , indicating a rapid decorrelation of particle locations. As regards the RRV, $\langle |W_r| \rangle$ is directly related to the particle velocity fluctuations for the largest separation distances (Jiang et al. Reference Jiang, Zhao and Brandt2026), so the value of $\langle |W_r| \rangle$ at $r \gtrsim 3D$ decreases with the Galileo number (shown in figure 15 b) owing to the reduced particle velocity fluctuations (see figure 5 a). At small separations, instead, the particle relative motion is dominated by the interparticle hydrodynamic interactions. The RRV decreases rapidly as the separation distance approaches $r \approx D$ owing to the significant lubrication effects among nearby particles (Jiang et al. Reference Jiang, Xu and Zhao2024b ). In addition, it is noteworthy that the value of $\langle |W_r (r=D)| \rangle$ is highest at ${\textit{Ga}}=120$ among all cases, which explains the highest collision kernel shown in figure 14 for the highest- ${\textit{Ga}}$ case, although the value of $g(r=D)$ is slightly reduced at ${\textit{Ga}}=120$ (see the inset in figure 15 a).

At last, to further examine the effect of gravity on the particle-pair dynamics, we compute the particle-pair correlation function $g(r,\psi )$ and the absolute RRV $|W_r(r,\psi )|$ as two-dimensional functions of the independent variables $r$ (the radial separation) and $\psi$ (the polar angle relative to the vertical direction). This allows us to further understand the directional information about the particle spatial distribution and relative motions. Note that one can obtain the one-dimensional functions $g(r)$ and $\langle |W_r(r)|\rangle$ by integrating $g(r,\psi )$ and $|W_r(r,\psi )|$ over the polar angle $\psi$ as (Yin & Koch Reference Yin and Koch2007)

(3.20) \begin{align} g(r) &= \frac {1}{2}\int _0^\pi {g(r,\psi )\sin \psi {\rm d}\psi }, \end{align}

(3.21) \begin{align} \langle |W_r(r)|\rangle &= \frac {1}{2}\int _0^\pi {|W_r(r,\psi )| g(r,\psi ) \sin \psi {\rm d}\psi }. \end{align}

The functions $g(r,\psi )$ and $|W_r(r,\psi )|$ are depicted in figure 16 for the suspensions at ${\textit{Ga}}=0$ (no gravity) and ${\textit{Ga}}=120$ (strongest particle sedimentation). As shown in figure 16(a,c), $g(r,\psi )$ and $|W_r(r,\psi )|$ only vary with the radial separation $r$ for non-sedimenting particles ( ${\textit{Ga}}=0$ ), and do not exhibit any directional preference. In contrast, at ${\textit{Ga}}=120$ , we can clearly observe a directional bias in the spatial distribution and relative motion of particle pairs due to the sedimentation. Specifically, the data in figure 16(b) show that the probability of finding particle pairs along the horizontal direction is higher than that in the vertical direction, which could be related to the wake-induced hydrodynamic interactions among settling particles. This observation reminds us of the horizontally aligned particles settling in a quiescent fluid at a low Galileo number ${\textit{Ga}}\approx 18$ (Yin & Koch Reference Yin and Koch2007), but is in contrast to the column-like particle structures observed in settling-particle suspensions at high ${\textit{Ga}}$ (Uhlmann & Doychev Reference Uhlmann and Doychev2014). The different behaviours may be ascribed to the complexity of the DKT-like interactions at different particle Reynolds numbers or in the turbulent/non-turbulent environments.

Regarding the particle relative motion (cf. figure 16 d), the value of $|W_r(r,\psi )|$ is found to be greater for the horizontally aligned particle pairs ( $\psi =90^\circ$ ) than for the vertically aligned ones ( $\psi =90^\circ$ ) with a given radial separation distance $r$ . Note that the directional preference of $|W_r(r,\psi )|$ observed in the current study is opposite to that observed for settling particles in a quiescent fluid: $|W_r(r,\psi )|$ was found to be greater along the vertical direction than horizontal direction for settling particles in the absence of background turbulence (see e.g. figure 8e in Jiang et al. Reference Jiang, Zhao and Brandt2026). This difference is reminiscent of the opposite variation of the particle vertical velocity fluctuations as ${\textit{Ga}}$ increases in turbulent and quiescent fluids (see the discussion in § 3.2). The decrease of $\sigma (v_y)$ at ${\textit{Ga}}=120$ (see figure 5 a) is related to the smaller particle relative velocities along the vertical direction shown in figure 16(d), since the particle vertical velocity fluctuations can promote relative motion between the vertical-aligned particles. Moreover, considering the relationship (3.21), the large values of $g(r,\psi )$ and $|W_r(r,\psi )|$ in the region $(r,\psi ) \approx (D,90^{\circ })$ explain the higher value of $\langle |W_r(r=D)|\rangle$ in the suspension with ${\textit{Ga}}=120$ than for the non-sedimenting particles; this also results in an enhanced collision rate in the presence of gravitational sedimentation.