3.1. Particle settling velocity

The first observable we discuss is the mean settling velocity of dispersed particles. Here, the settling velocity is defined as the component of particle velocity along the gravitational direction, i.e. $V_s=\boldsymbol {v}\boldsymbol {\cdot } \boldsymbol {e}_g=-v_y$. As shown in figure 2, the mean settling velocity $\langle V_s \rangle$ exhibits a non-monotonic variation with the increase of particle volume fraction from $\phi =0.1\,\%$ to $\phi =10\,\%$. Specifically, the mean settling velocity is greater than the settling velocity of an isolated particle when $\phi \le 2\,\%$, with a peak value of $\langle V_s \rangle \approx 1.25 V_t$ at $\phi =1\,\%$, and decreases to less than $V_t$ when the particle volume fraction exceeds $5\,\%$. In the following, we look into the particle clustering, hindrance effect and the particle orientation to interpret the non-monotonic variation of the particle mean settling velocity with the change of the volume fraction.

Figure 2. Mean settling velocity of dispersed particles at different volume fractions. The empirical correlation of the hindered settling velocity (Richardson & Zaki Reference Richardson and Zaki1954) (depicted by the red dashed line) is included for comparison.

3.1.1. Particle clustering

In previous studies it has been reported that the enhancement of the particle mean settling velocity is highly related to the formation of particle clustering (Kajishima Reference Kajishima2004; Uhlmann & Doychev Reference Uhlmann and Doychev2014; Fornari et al. Reference Fornari, Ardekani and Brandt2018). Thus, we first examine the spatial distribution of dispersed particles in the present simulations using the Voronoi analysis (Monchaux, Bourgoin & Cartellier Reference Monchaux, Bourgoin and Cartellier2010). In this analysis method the entire computational domain is partitioned into $N_p$ cells, i.e. Voronoi tessellations. The partitioning rule ensures that a given spatial point inside the $i$th tessellation is closest to the centroid of the $i$th particle among all particles. Accordingly, the spatial distribution of dispersed particles can be quantified by the statistics of the normalized volume of Voronoi tessellations, $\bar {V}_{Voro}(i)={V}_{Voro}(i) N_p/V_{tot}$, where ${V}_{Voro}(i)$ is the volume of the $i$th Voronoi tessellation and $V_{tot}$ is the volume of the whole domain. If particles are orderly distributed in the space (like molecules in a crystal), the entire domain would be evenly partitioned so that $\bar {V}_{Voro}\equiv 1$ and the standard deviation is $\sigma (\bar {V}_{Voro})=0$. In contrast, in a system where particle clustering arises, the prevalence of particle accumulations (represented by small values of $\bar {V}_{Voro}$) and voids (represented by large values of $\bar {V}_{Voro}$) would increase the intermittency of the probability density function (p.d.f.) of the standard deviation of $\bar {V}_{Voro}$ (Monchaux et al. Reference Monchaux, Bourgoin and Cartellier2010).

To quantify the degree of particle clustering, Tagawa et al. (Reference Tagawa, Roghair, Prakash, Van Sint Annaland, Kuipers, Sun and Lohse2013) defined the clustering indicator $C$ based on the standard deviation of the normalized Voronoi volume of dispersed particles as

(3.1)\begin{equation} C = \sigma \left( {{{\bar{V} }_{Voro}}} \right)/\sigma \left( {{{\bar{V} }_{Voro,rand}}} \right), \end{equation}

where the subscript ‘$rand$’ represents the assembly of particles with a random spatial distribution. According to this definition, particles are considered to form clusters when the clustering indicator exceeds unity, and a higher level of clustering is identified by a greater value of $C$. As for point-like particles, the statistics of $\bar {V}_{Voro,rand}$ conforms to a gamma distribution, yielding a standard deviation of $\sigma (\bar {V}_{Voro,rand})=0.447$ (Ferenc & Néda Reference Ferenc and Néda2007). However, the value of $\sigma (\bar {V}_{Voro,rand})$ is a decreasing function of $\phi$ for non-overlapping finite-size particles (Uhlmann Reference Uhlmann2020). In the present work, we generate randomly distributed prolate particles with random orientations following the method proposed by Anoukou et al. (Reference Anoukou, Brenner, Hong, Pellerin and Danas2018) and compute $\sigma (\bar {V}_{Voro,rand})$ accordingly. To ensure the convergence of $\sigma (\bar {V}_{Voro,rand})$ with sufficient samples, we repeat the generating process 1000 times for $\phi =0.1\,\%$, 200 times for $\phi =0.5\,\%$ and 100 times for other volume fractions.

The statistical distributions of the normalized Voronoi volume in the current work are illustrated in figure 3(a). We can clearly observe a raised tail for the p.d.f. of $\bar {V}_{Voro}$ in the cases with $0.5\,\% \le \phi \le 2\,\%$. In contrast, the distribution of $\bar {V}_{Voro}$ is narrowed in the densest suspension at $\phi =10\,\%$. Furthermore, the variation of the clustering indicator $C$ as the function of the volume fraction is illustrated in figure 3(b). Interestingly, the value of $C$ varies non-monotonically with a peak at $\phi =1\,\%$, which coincides with the highest mean settling velocity of particles as is observed in figure 2. For the cases with a lower or higher volume fraction, the clustering indicator is reduced, although its value remains to be greater than unity. Thus, the particle clustering becomes less pronounced in more dilute or denser suspensions.

Figure 3. Results of the Voronoi analysis at different volume fractions. (a) The p.d.f. of the normalized volume of Voronoi tessellations. (b) Clustering indicator $C$ at different volume fractions $\phi$.

In the previous studies of particle sedimentation the DKT-like interactions among settling particles are regarded as the essential mechanism in the formation of particle clusters (Kajishima Reference Kajishima2004; Zaidi et al. Reference Zaidi, Tsuji and Tanaka2014; Fornari et al. Reference Fornari, Ardekani and Brandt2018; Moriche et al. Reference Moriche, Hettmann, García-Villalba and Uhlmann2023). In particular, during the drafting stage of a DKT event, a pair of particles can attract each other, reducing the distance between them. Furthermore, if these interacting particles attract additional particles before they separate, the number of accumulated particles can increase progressively, which eventually results in particle clustering in the suspension (Moriche et al. Reference Moriche, Hettmann, García-Villalba and Uhlmann2023). Compared with settling spheres, spheroidal particles are more likely to be drawn into the wake of a leading particle and tend to have a longer interaction time in the DKT process (Ardekani et al. Reference Ardekani, Costa, Breugem and Brandt2016). This explains the occurrence of clustered prolate particles in this study, similar to the behaviour of settling oblate particles (Fornari et al. Reference Fornari, Ardekani and Brandt2018; Moriche et al. Reference Moriche, Hettmann, García-Villalba and Uhlmann2023), in contrast to the absence of particle clustering of settling spheres at a comparable Reynolds number (Zaidi et al. Reference Zaidi, Tsuji and Tanaka2014).

In figure 4 we provide the visualization of the flow system at three typical volume fractions $\phi =0.1\,\%$, $1\,\%$ and $10\,\%$. In the most dilute case with $\phi =0.1\,\%$ (figure 4a,d), the particles are sparsely distributed in the space without appreciable particle clustering. However, individual DKT events induced by the wake-related hydrodynamic interactions can still be found (see figure 4g,h). As for the case with $\phi =1\,\%$, while, we can evidently observe locally accumulated particles and some regions devoid of particles (see figure 4e). Meanwhile, large-scale flow structures formed by the interconnected particle wakes are also illustrated in figure 4(b). These structures exhibit the footprint of the column-like particle clusters meandering along the vertical direction (with more details provided in § 3.2). While, in another limit of dense suspension, particles are crowded in the space as shown in figure 4(c, f). The too small distance between neighbouring particles frequently perturbs particle wakes, so as to inhibit the formation of particle clustering (Zaidi et al. Reference Zaidi, Tsuji and Tanaka2014). Incidentally, the slight growth of the clustering indicator $C$ from $\phi =5\,\%$ to $10\,\%$ (see figure 3b) may be related to the more ordered arrangement for the randomly distributed particles, which reduces $\sigma (\bar {V}_{Voro,rand})$ in (3.1).

Figure 4. (a–c) Snapshots of the instantaneous flow field at (a) $\phi =0.1\,\%$, (b) $\phi =1\,\%$ and (c) $\phi =10\,\%$. Dispersed prolate particles are depicted in grey. The background contour represents the vertical fluid velocity $u_y$ (normalized by $V_t$), with isosurfaces of $u_y=-\langle V_s \rangle$ shown in red. (d–f) Vertical sections with a thickness of $D_{eq}$ along the $z$ direction taken from panels (a–c). Dispersed prolate particles are shown in blue. Isosurfaces of the Q-criterion at $Q=0.2V_t^2/D_{eq}^2$ are coloured by the magnitude of vorticity $\|\boldsymbol {\varOmega }\|$ (normalized by $V_t/D_{eq}$). (g,h) Zoom-in views of panel (a) to illustrate the touching particle pairs at $\phi =0.1\,\%$.

Moreover, we also look into the statistics of the nearest neighbour distance (NND) of particles in the present simulations. Here, the NND distance of the $i$th particle, denoted by $d_{NN}(i)$, is defined by (Zaidi et al. Reference Zaidi, Tsuji and Tanaka2014)

(3.2)\begin{equation} d_{NN} (i) = \mathop {\min }_{j = 1,2,\ldots,{N_p},j \ne i} \left\| {{{\boldsymbol{x}}_i} - {{\boldsymbol{x}}_j}} \right\|, \end{equation}

where ${\boldsymbol {x}}_i$ is the centroid position of the $i$th particle. In figure 5 the p.d.f.s of $d_{NN}$ at different volume fractions are provided, with a comparison with that of the randomly distributed particles. It is observed that the statistical distribution of $d_{NN}$ shifts to the side of smaller values in all cases. This observation indicates that the dispersed particles become more locally crowded than the randomly distributed particles, which is consistent with the clustering nature of particles manifested by the Voronoi analysis (see figure 3). In addition, the probability of finding touching particles with a centre-to-centre distance $d_{NN}=2b$ is increased, which can be ascribed to the DKT-like interactions between nearby particles. Especially, there is a noticeable secondary peak of the p.d.f. of $d_{NN}$, evaluated at $d_{NN}=2b$, in the most dilute suspension at $\phi =0.1\,\%$, revealing the prevalence of touching particle pairs undergoing the kissing stage of the DKT process. The touching particle pairs, however, would become less stable with the intensified hydrodynamic disturbances and more frequent inter-particle collisions as $\phi$ increases. As a result, this secondary peak of the p.d.f. of $d_{NN}$ becomes invisible at higher volume fractions. Moreover, we also calculate the ensemble average of $d_{NN}$, denoted by $\langle d_{NN} \rangle$, of each case (not presented here). The results show that in the most dilute case at $\phi =0.1\,\%$, the value of $\langle d_{NN} \rangle$ is $4.0D_{eq}$, considerably larger than that of the case with the strongest particle clustering (i.e. $\langle d_{NN} \rangle =1.9 D_{eq}$ at $\phi =1\,\%$). To make a fair comparison, we normalize $\langle d_{NN} \rangle$ by that of a particle assembly with a random spatial distribution (denoted by $\langle d_{NN} \rangle _{rand}$), and obtain $\langle d_{NN} \rangle /\langle d_{NN} \rangle _{rand}=0.93$ for $\phi =0.1\,\%$ and 0.70 at $\phi =1\,\%$. This indicates that particles are closer to the random distribution at the lowest volume fraction, consistent with the attenuation of the particle clustering as $\phi$ decreases from 1 % to 0.1 % obtained by the Voronoi analysis. However, this observation is in a qualitative disagreement with the intensified clustering trend for settling spherical particles at $Ga=178$ with the volume fraction decreasing from $\phi =0.5\,\%$ to $\phi =0.05\,\%$ (Doychev Reference Doychev2014). We speculate that the disagreement is related to the weaker fluid inertia effect (characterized by the lower value of $Ga$ and also $Re_t$) in the present study. The wake-induced hydrodynamic interactions, which are essential for the attraction among settling particles, may not be strong enough to make sparsely distributed particles form clusters at a very low volume fraction. This argument, however, needs to be further examined by the simulations with the volume fraction lower than 1 %.

Figure 5. The p.d.f. of the NND of particles at volume fractions of (a) $\phi =0.1\,\%, 0.5\,\%, 1\,\%$ and (b) $\phi =2\,\%, 5\,\%, 10\,\%$. The solid lines are the results of the present simulations while the dashed lines represent the results of randomly distributed prolate spheroids with random orientation. Each line starts from $d_{NN}=2b$, which is the smallest centre-to-centre distance between two finite-sized prolate particles.

Then, we further investigate the relationship between the clustering and settling velocity of particles in the present flow system. Moriche et al. (Reference Moriche, Hettmann, García-Villalba and Uhlmann2023) reported the positive correlation between the standard deviation of Voronoi volumes and the mean settling velocity of low-aspect-ratio oblate particles. Here, we compute the averaged settling velocity conditioned on the Voronoi volume, denoted by $\langle V_s \rangle _{V_{Voro}}$, and present the results in figure 6. It is shown that particles with smaller Voronoi tessellations tend to settle faster, irrespective of the volume fraction. This correlation can be explained by the so-called ‘swarm effect’ (Koch & Hill Reference Koch and Hill2001; Wang et al. Reference Wang, Zhang, Li, Xiao and Yang2022), which suggests that a cluster of settling particles experiences lower total drag compared with the same number of individual particles. To gain further understanding, we compute the fluid velocity sampled by the particle, denoted by $\boldsymbol {u}^{f@p}$, by averaging the local fluid velocity on the surface of a sphere centred at the particle centroid with a radius of $1.5 D_{eq}$ (Kidanemariam et al. Reference Kidanemariam, Chan-Braun, Doychev and Uhlmann2013). Figure 7 shows the joint p.d.f. of the Voronoi volume and the sampled vertical fluid velocity of the particle. In this figure we observe a positive correlation between $\bar {V}_{Voro}$ and $u_y^{f@p}$, revealing that the particles in the clustering regions (represented by the small value of $\bar {V}_{Voro}$) are prone to sample downward fluid flows, while in the void zones (represented by the large value of $\bar {V}_{Voro}$) particles tend to experience stronger upward flows. This observation can be attributed to the fact that the clustered particles are more likely to reside in the wake of other particles, where the downward flow is dominated. On the contrary, the fluid moves upwards in the void regions so as to decelerate the particle settling motion. However, this correlation becomes less pronounced at $\phi \geq 5\,\%$, seemingly due to the disruption of particle wakes and the diminished distinction between the ‘wake region’ and ‘void region’ in dense suspensions.

Figure 6. Averaged settling velocity of dispersed particles conditioned on the Voronoi volume at different volume fractions.

Figure 7. Joint p.d.f. (scaled by its maximum value) of the normalized Voronoi volume $\bar {V}_{Voro}$ and the vertical fluid velocity sampled by particles $u_y^{f@p}$ at (a) $\phi =0.1\,\%$, (b) $\phi =0.5\,\%$, (c) $\phi =1\,\%$, (d) $\phi =2\,\%$, (e) $\phi =5\,\%$ and ( f) $\phi =10\,\%$. The horizontal and vertical dashed lines represent the mean value of $u_y^{f@p}$ and $\bar {V}_{Voro}$, respectively.

Furthermore, we also examine the relationship between the translational velocity of the particle and the local fluid velocity seen by the particle. In figure 8 we present the joint p.d.f. of the vertical component of the particle velocity and the particle-sampled fluid velocity. The results demonstrate that the particles tend to settle rapidly when experiencing vertical downward flows (with the negative value of $\langle u_y^{f@p} \rangle$), and vice versa. The correlations presented in figures 7 and 8 altogether can account for the decreased settling velocity of particles with larger Voronoi volumes shown in figure 6. In addition, we also compute the average value of the fluid vertical velocity sampled by particles, denoted by $\langle u_y^{f@p} \rangle$, and compare it with the ensemble averaged fluid vertical velocity, $\langle u_y \rangle _f$, in figure 9(a). Interestingly, $\langle u_y^{f@p} \rangle$ is always smaller than $\langle u_y \rangle _f$, regardless of the volume fraction. This observation reveals the preferential sampling of downward fluid flows by dispersed particles, which was also reported by Uhlmann & Doychev (Reference Uhlmann and Doychev2014) for settling spheres at $Ga=178$. Moreover, the difference between $\langle u_y^{f@p} \rangle$ and $\langle u_y \rangle _f$ is largest at $\phi =1\,\%$, corresponding to the strongest particle clustering with varying volume fraction. Therefore, the preferential sampling of downward flows, which is most significant for the strongest particle clustering, is the underlying mechanism of the aforementioned swarm effect to enhance the particle mean settling velocity.

Figure 8. Joint p.d.f. (scaled by its maximum value) of the vertical fluid velocity sampled by particles, $u_y^{f@p}$, and the particle vertical velocity, $v_y$, at (a) $\phi =0.1\,\%$, (b) $\phi =0.5\,\%$, (c) $\phi =1\,\%$, (d) $\phi =2\,\%$, (e) $\phi =5\,\%$ and ( f) $\phi =10\,\%$. The horizontal and vertical dashed lines represent the mean value of $v_y$ and $u_y^{f@p}$, respectively.

Figure 9. (a) Mean vertical fluid velocity sampled by the particles, $\langle u_y^{f@p} \rangle$, and the ensemble averaged velocity of the fluid flow, $\langle u_y \rangle _f$. (b) Averaged relative velocity between the particle motion and the mean flow, $U_{rel}$, and between the particle motion and the local fluid flow, $U_{rel}^L$.

Additionally, we also look into the relative motion between the particle and fluid phases. Here we define $U_{rel}=v_y-\langle u_y \rangle _f$ as the relative vertical velocity between the particle motion and the mean flow, and $U_{rel}^L=v_y-u_y^{f@p}$ as the local relative velocity. In figure 9(b) we provide the ensemble average of these two relative velocities, denoted by $\langle U_{rel} \rangle$ and $\langle U_{rel}^L \rangle$, at different volume fractions. It is observed that the variation of $\langle U_{rel}^L \rangle$ with increasing $\phi$ is alleviated compared with that of $\langle U_{rel} \rangle$. Especially in dilute cases with $\phi \leq 1\,\%$, the difference of $\langle U_{rel}^L \rangle$ among different cases is less than $4\,\%$, similar to the observation for settling spheres in dilute suspensions (Doychev Reference Doychev2014). Therefore, the variation of the global relative particle–fluid velocity $\langle U_{rel} \rangle$ as $\phi$ changes can be substantially attributed to the different level of particle clustering and the preferential sampling of the fluid velocity. More discussion about the variation of $\langle U_{rel}^L \rangle$ is provided in § 3.1.3.

3.1.2. Hindrance effect

Let us now turn to the reduced particle mean settling velocity (i.e. $\langle V_s \rangle < V_t$) in dense suspensions at $\phi \geq 5\,\%$ (see figure 2). The ensemble averaged fluid velocity, which can be calculated by the flux conservation of the whole system as $\langle u_y \rangle _f=\phi /(1-\phi ) \langle V_s \rangle$ (Yin & Koch Reference Yin and Koch2007), is enhanced as $\phi$ increases (see figure 9a). The enhanced upward fluid flow has an opposite effect upon the sampling of downward flows by particles. In the meantime, as $\phi$ increases, particle clustering is attenuated (with the clustering indicator $C$ decreasing) and the preferential sampling of downward flows becomes less pronounced (with $\langle u_y^{f@p} \rangle$ approaching $\langle u_y \rangle _f$). Consequently, the value of $\langle u_y^{f@p} \rangle$ decreases in magnitude when $\phi >1\,\%$ and even becomes positive at the highest volume fraction $\phi =10\,\%$. As a result, the hindrance effect becomes predominant when the volume fraction exceeds approximately 5 %, leading to the reduction of the mean settling velocity in this regime.

As for the sedimentation of spherical particles, Richardson & Zaki (Reference Richardson and Zaki1954) proposed the well-known empirical formula of the hindered settling velocity as a function of the particle volume fraction, i.e.

(3.3)\begin{equation} \langle V_s \rangle/V_t=(1-\phi)^n. \end{equation}

The exponent $n$ in (3.3) was found to be an decreasing function of the settling Reynolds number $Re_t$ and can be fitted by (Garside & Al-Dibouni Reference Garside and Al-Dibouni1977)

(3.4)\begin{equation} \frac{5.1-n}{n-2.7}=0.1Re_t^{0.9}. \end{equation}

In figure 2 we also depict the empirical hindered settling velocity as a function of $\phi$ given by (3.3), with the exponent $n=3.17$ obtained by substituting $Re_t=61.8$ in (3.4). It is shown that the reduced settling velocity observed in the present simulations at $\phi \geq 5\,\%$ approaches the prediction by the empirical formula (3.3). The remaining discrepancy can be ascribed to the weak effect of clustering and the change of orientation (see the discussion on figure 10 in the following) of settling prolate particles.

Figure 10. The statistics of the orientation of settling prolate particles at different volume fractions. (a) The p.d.f. of the cosine value of the pitch angle $\psi$. (b) Mean value of $|\cos \psi |$ as a function of the volume fraction.

3.1.3. Particle orientation

At last, we would like to study the orientation of settling prolate particles and its influence on the particle settling motion. First, we compute the statistics of particle orientation in the present simulations and display the results in figure 10. It is shown that in the cases with low volume fractions the broad-side-on orientation (corresponding to $|\cos \psi |=0$) of settling prolate spheroids still prevails. This is the stable orientation of an isolated settling prolate spheroid under the effect of the fluid inertia torque (Ardekani et al. Reference Ardekani, Costa, Breugem and Brandt2016; Dabade et al. Reference Dabade, Marath and Subramanian2016). However, with the increasing volume fraction, the orientation of particles progressively shifts towards a random distribution, demonstrated by the flattening of the p.d.f. of $|\cos \psi |$ and the corresponding increase in the average value, $\langle |\cos \psi |\rangle$. This observation manifests the overwhelming effect of particle–particle interactions to perturb the stable orientation of settling prolate spheroids in dense suspensions. Then, we also examine the correlation between the settling velocity and the orientation of particles in the present flow system by computing the joint p.d.f. of these two quantities. As shown in figure 11, prolate particles tend to settle faster as their orientation deviates more from the broad-side-on alignment, irrespective of the volume fraction, just as the case of an isolated settling prolate spheroid.

Figure 11. Joint p.d.f. (scaled by its maximum) of the absolute cosine of the particle pitch angle and the vertical velocity of the particle at (a) $\phi =0.1\,\%$, (b) $\phi =0.5\,\%$, (c) $\phi =1\,\%$, (d) $\phi =2\,\%$, (e) $\phi =5\,\%$ and ( f) $\phi =10\,\%$. The solid line represents the averaged settling velocity conditioned on the pitch angle (the data with ${\rm p.d.f.} (| \cos \psi |) <0.01$ are masked). The dashed line represents the variation of an isolated settling prolate spheroid with an artificially fixed pitch angle.

Furthermore, we examine the influence of particle orientation on the settling motion. As for a single spheroid in a uniform flow, the drag coefficient depends on the attack angle between the symmetry axis and the incoming flow (Zastawny et al. Reference Zastawny, Mallouppas, Zhao and van Wachem2012). To examine the case in the suspension, we compute the mean drag coefficient of the dispersed particles, $\overline {C_d}$, in each case under consideration. The definition of the mean drag coefficient is based on the mean hydrodynamic drag force and the mean local particle–fluid relative velocity, i.e.

(3.5)\begin{equation} \overline{C_d}=\frac{\langle F_y^H\rangle}{\dfrac{1}{2}\rho_f\langle U_{rel}^L\rangle^2A}, \end{equation}

where $F_{y}^H$ is the vertical component of the hydrodynamic force acting on the particle and $A={\rm \pi} D_{eq}^2/4$ the characteristic frontal area of the particle. Note that the time average of $F_{y}^H$ is in balance with the buoyancy of the particle in the statistically steady state, i.e. $\langle F_y^H\rangle ={\rm \pi} (\rho _p-\rho _f)g D_{eq}^3/6$. Thus, according to the definition (3.5), the variation of $\langle U_{rel}^L\rangle$ shown in figure 9(b) is determined by $\overline {C_d}$ for different $\phi$. As shown in figure 12, $\overline {C_d}$ first decreases slightly as $\phi$ increases from $0.1\,\%$ to $0.5\,\%$, which could be attributed to the increased deviation from the broad-side-on orientation of settling prolate particles. However, $\overline {C_d}$ then increases monotonically with the volume fraction when $\phi > 0.5\,\%$. This cannot be explained by the change of particle orientation since the dispersed particles deviate more from the broad-side-on orientation as $\phi$ continues to grow (see figure 10). Therefore, we would like to conclude that the change of particle orientation plays a minor role in the mean settling velocity of dispersed particles, although the settling velocity of individual particles is still orientation dependent. Incidentally, the increasing trend of the mean drag coefficient with the increasing particle volume fraction was also reported in the studies of the flow past a fixed or mobile assembly of spherical particles (Tenneti, Garg & Subramaniam Reference Tenneti, Garg and Subramaniam2011; Tavanashad, Passalacqua & Subramaniam Reference Tavanashad, Passalacqua and Subramaniam2021). This trend can be interpreted by the change of the local flow conditions in the vicinity of the particle. Specifically speaking, with the increase of the volume fraction, the fluctuations of the particle and fluid velocities grow with the intensified particle–particle and particle–fluid interactions (see figure 18a for more details). As a consequence, the dispersed particles exposed to turbulent local flows would experience increased unsteady and nonlinear drags (Fornari, Picano & Brandt Reference Fornari, Picano and Brandt2016), which contributes to the drag enhancement of settling particles.

Figure 12. Mean drag coefficient as a function of the particle volume fraction.

3.2. Particle microstructures and particle–fluid interactions

According to the discussion in § 3.1.1, the spatial distribution of settling particles is non-uniform in the suspensions. Therefore, it is of interest to examine the particle microstructures. First, we calculate the particle pair distribution function $P(\boldsymbol {r})$, which provides the information about the probability of finding another particle relative to a reference particle with a separation vector $\boldsymbol {r}$. The definition of $P(\boldsymbol {r})$ is referred to Yin & Koch (Reference Yin and Koch2007), Zaidi et al. (Reference Zaidi, Tsuji and Tanaka2015) and Fornari et al. (Reference Fornari, Ardekani and Brandt2018). By definition, $P(\boldsymbol {r}) > 1$ indicates a higher probability of finding a pair of particles with a separation of $\boldsymbol {r}$ compared with the uniform distribution of particles. In addition, as $P(\boldsymbol {r})$ is axisymmetric about the direction of gravity, we calculate the average of $P(\boldsymbol {r})$ over the isotropic horizontal plane (i.e. $x-z$ plane), and obtain the pair distribution function as a function of the separation distance $r=\| \boldsymbol {r} \|$ and the polar angle $\varphi$ (the angle between the positive $y$ direction and the vector $\boldsymbol {r}$), i.e.

(3.6)\begin{equation} P\left( {r,\varphi } \right) = \frac{1}{{2{\rm \pi} }}\int_0^{2{\rm \pi} } {P\left( {\boldsymbol{r}} \right)\,{\rm d}\theta } = \frac{1}{{2{\rm \pi} }}\int_0^{2{\rm \pi} } {P\left( {r,\theta ,\varphi } \right)\,{\rm d}\theta }. \end{equation}

Here, the variable of integration $\theta$ is the azimuth angle between the positive $x$ direction and the projection of vector $\boldsymbol {r}$ onto the horizontal plane.

In figure 13 we display the particle pair distribution function at different volume fractions. It is observed that the value of $P( {r,\varphi } )$ is significantly greater than unity along the vertical direction in dilute suspensions. Therefore, we can infer that the dispersed particles prefer to form column-like structures. To gain more information, we also compute the probability of observing attracting particle pairs with the separation vector $\boldsymbol {r}$, denoted by $P_{in}(\boldsymbol {r})$. Here, an attracting particle pair is identified if the relative radial velocity $W_r$ between the two particles is negative, i.e.

(3.7)\begin{equation} W_r^{\left\{ {i,j} \right\}} = \frac{{\left( {{{\boldsymbol{v}}_i} - {{\boldsymbol{v}}_j}} \right) \boldsymbol{\cdot} \left( {{{\boldsymbol{x}}_i} - {{\boldsymbol{x}}_j}} \right)}}{{\left\| {{{\boldsymbol{x}}_i} - {{\boldsymbol{x}}_j}} \right\|}}<0, \end{equation}

where $i$ and $j$ represent the indices of the two particles. Same as the particle pair function, we compute the average of $P_{in}(\boldsymbol {r})$ over the isotropic horizontal directions and obtain the two-dimensional function of $P_{in}(r,\varphi )$, as shown in figure 14. In dilute suspensions, particle pairs exhibit a strong tendency to attract each other along the vertical direction (indicated by $P_{in}>0.5$), but behave in a repulsive manner along the horizontal direction. We attribute these observations to the DKT-like interactions among settling prolate particles. As a result, the attraction and entrapment of particles in the wake of leading ones give rise to the column-like particle microstructures, consistent with the results shown in figure 13. In contrast to our result, it was reported that settling spheres tend to form particle deficits along the vertical direction, while the spatial distribution of cubic particles is more uniform under similar conditions ($Re_t<70$ and $\phi \approx 1\,\%$) (Yin & Koch Reference Yin and Koch2007; Zaidi Reference Zaidi2018b; Seyed-Ahmadi & Wachs Reference Seyed-Ahmadi and Wachs2021). These differences highlight the effect of particle shape on the wake-induced hydrodynamic interactions among settling particles. For prolate spheroids, the attraction between particle pairs in the DKT-like events is strong enough to entrap particles in the wake regions (Ardekani et al. Reference Ardekani, Costa, Breugem and Brandt2016). Regarding spherical particles, the shear-induced lift force dominates in wake regions, pushing trailing particles outside the wake (Yin & Koch Reference Yin and Koch2007). While, settling cubic particles are found to have greater rotational rates, which generate lateral Magnus forces that help them escape from the wake of leading particles. Hence, cubic particles are less likely to form obvious microstructures (Seyed-Ahmadi & Wachs Reference Seyed-Ahmadi and Wachs2019, Reference Seyed-Ahmadi and Wachs2021). Additionally, due to the disruption of particle wakes, the probability of observing attracting/repelling particle pairs is reduced, and the regions where these events dominate diminish with the increasing volume fraction (see figure 14). As a result, the column-like microstructures are gradually attenuated and eventually becomes negligible at $\phi \ge 5\,\%$.

Figure 13. Pair distribution function of particles at (a) $\phi =0.1\,\%$, (b) $\phi =0.5\,\%$, (c) $\phi =1\,\%$, (d) $\phi =2\,\%$, (e) $\phi =5\,\%$ and ( f) $\phi =10\,\%$. Here, the horizontal and vertical directions are denoted by $r_x=r \sin \varphi$ and $r_y=r \cos \varphi$, respectively.

Figure 14. Probability of observing attracting particle pairs as a function of the separation $\boldsymbol {r}$. Results are shown for (a) $\phi =0.1\,\%$, (b) $\phi =0.5\,\%$, (c) $\phi =1\,\%$, (d) $\phi =2\,\%$, (e) $\phi =5\,\%$, ( f) $\phi =10\,\%$.

Furthermore, we quantify the radial characteristics of particle microstructures by computing the radial distribution function (RDF) of the dispersed particles. The RDF, denoted by $g(r)$, is calculated from the two-dimensional pair distribution function $P(r,\varphi )$ by (Yin & Koch Reference Yin and Koch2007)

(3.8)\begin{equation} g\left( r \right) = \frac{1}{2}\int_0^{\rm \pi} {P\left( {r,\varphi } \right)\sin \varphi \,{\rm d}\varphi }. \end{equation}

As shown in figure 15(a), the spatial correlation of particle pairs is considerably increased with $g(r)>1$ for a small separation distance $r$ in dilute suspensions. As the separation distance $r$ grows, the RDF gradually decays to $g(r) \sim 1$, indicating the recovery to the uniform distribution for particle pairs with long-distance separations. The peak value of $g(r)$, which is evaluated at the separation distance $r=2b$ at $\phi \le 2\,\%$, is a monotonically decreasing function of the volume fraction. While, in dense suspensions with $\phi \ge 5\,\%$, the RDF is close to unity for all separation distances, corresponding to the attenuated particle microstructures in these cases.

Figure 15. (a) Radial distribution function and (b) order parameter of particle pairs as functions of the radial separation distance $r$.

In previous studies on settling particles the degree of particle clustering is usually quantified by the maximum value of the RDF (Zaidi et al. Reference Zaidi, Tsuji and Tanaka2014; Fornari et al. Reference Fornari, Ardekani and Brandt2018; Zaidi Reference Zaidi2018a). However, in the present work, we find that the monotonic decrease of the maximum value of $g(r)$ with the increasing volume fraction does not align with the non-monotonic variation of the clustering indicator based on the Voronoi analysis (see figure 3b). To interpret this inconsistency, which primarily occurs for $\phi <1\,\%$, we recall the statistics of the NND provided in figure 5. By comparing these two statistics, the RDF and the NND, we attribute the peak of the RDF at $r=2b$ to the prevalence of touching particle pairs in the suspensions, corresponding to the rise of the p.d.f. of $d_{NN}$ at $d_{NN}=2b$. In other words, the RDF actually reflects the pairwise information of the dispersed particles instead of the particle clustering, which is a concept in a global sense. The subtle difference between these two quantities becomes especially evident in the most dilute case with $\phi =0.1\,\%$, where the presence of touching particles involved in individual DKT events (see figure 4g,h) significantly increases the value of $g(r)$ at $r=2b$ and gives rise to the secondary peak of the p.d.f. of $d_{NN}$ at the same distance. Therefore, we shall argue that the maximum value of the RDF is not a proper criterion to measure particle clustering in the present flow system.

In addition, we also compute the order parameter $\langle P_2 \rangle (r)$ to measure the orientational feature of the particle microstructures. The order parameter is defined as the angular average of the second Legendre polynomial (Yin & Koch Reference Yin and Koch2007; Fornari et al. Reference Fornari, Ardekani and Brandt2018)

(3.9)\begin{equation} \left\langle {{P_2}} \right\rangle (r) = \frac{{\int_0^{\rm \pi} P (r,\varphi ){P_2}(\cos \varphi )\sin \varphi \,{\rm{d}}\varphi }}{{\int_0^{\rm \pi} P (r,\varphi )\sin \varphi \,{\rm{d}}\varphi }}, \end{equation}

in which ${P_2}( {\cos \varphi } ) = ( {3{{\cos }^2}\varphi - 1} )/2$. The value of $\langle P_2 \rangle (r)$ would be equal to 1 if all particle pairs with the separation $r$ are vertically aligned, 0 for the isotropic arrangement and -1/2 if the particle pairs are horizontally aligned (Yin & Koch Reference Yin and Koch2007). In figure 15(b) we depict the order parameter as a function of the radial separation distance $r$ at different volume fractions. In all cases under consideration, the order parameter is greater than zero at $r=2b$, indicating the tendency of nearby particle pairs to align vertically. This phenomenon again manifests the DKT-like hydrodynamic interactions among settling particles in the present system, which also make a difference even at $\phi \geq 5\,\%$ with weak particle clustering. While, the peak value of $\langle P_2 \rangle$ decreases as the volume fraction increases, indicating the weakening influence of particle wakes. Additionally, $\langle P_2 \rangle (r)$ decays with the separation distance $r$ rapidly to $\langle P_2 \rangle \sim 0$, indicating the recovery to an isotropic particle arrangement, in the suspensions with $\phi \ge 2\,\%$. However, for the lower volume fractions, the order parameter does not completely decay to zero with a small residual positive value even for the long-distance particle separation of $r\sim 10 D_{eq}$. This indicates that the wake-induced hydrodynamic interactions have a long working distance for sparsely distributed particles in dilute suspensions.

Moreover, to demonstrate the fluid–particle interactions in the present flow system, we examine the instantaneous vertical average of the particle concentration (denoted by $\langle \zeta \rangle _y$) and fluid vertical velocity (denoted by $\langle {{u_y}} \rangle _y$) at different volume fractions. The definitions of $\langle \zeta \rangle _y$ and $\langle {{u_y}} \rangle _y$ are as follows:

(3.10)$$\begin{gather} \langle \zeta \rangle_y \left( {x,z} \right) = \frac{1}{{{L_y}}}\int_0^{{L_y}} {\zeta \left( {x,y,z} \right)\,{{\rm d} y}}, \end{gather}$$

(3.11)$$\begin{gather}\langle {{u_y}} \rangle_y \left( {x,z} \right) = \frac{1}{{{L_y}}}\int_0^{{L_y}} {\left[ {1 - \zeta \left( {x,y,z} \right)} \right]{u_y}\left( {x,y,z} \right)\,{{\rm d} y}}. \end{gather}$$

Here, $\zeta (x,y,z)$ is an indicator function that is equal to 1 if a spatial point $(x,y,z)$ locates inside a particle, otherwise $\zeta (x,y,z)=0$. As shown in figure 16, the distribution of $\langle \zeta \rangle _y$ on the horizontal plane is far from uniform in dilute suspensions, indicating the formation of column-like particle microstructures. Also, we can evidently observe spatial correlation between the high value of $\langle \zeta \rangle _y$ and the extreme negative value of $\langle {{u_y}} \rangle _y$, and vice versa. In the regions where particles accumulate, the fluid flow moves downwards due to the drag by settling particles, but moves upwards in the regions devoid of particles to keep zero net flux of the whole system. However, as the volume fraction increases, the structures for $\langle \zeta \rangle _y$ and $\langle {{u_y}} \rangle _y$ become fragmented and the abovementioned correlation is weakened, owing to the diminishing particle microstructures in dense suspensions.

Figure 16. Instantaneous particle concentration $\langle \zeta \rangle _y$ (scaled by its maximum and represented by the background coloured contour) and fluid vertical velocity $\langle {{u_y}} \rangle _y$ (represented by the contour lines) averaged over the vertical direction at different volume fractions. Results are shown for (a) $\phi =0.1\,\%$, (b) $\phi =0.5\,\%$, (c) $\phi =1\,\%$, (d) $\phi =2\,\%$, (e) $\phi =5\,\%$, ( f) $\phi =10\,\%$.

At last, we look into the statistics of the velocity fluctuations of the particle and fluid phases under the effect of the complicated particle–fluid interactions. In figure 17 we illustrate the p.d.f.s of the fluctuation velocity of the particle and fluid vertical motions. In general, the magnitudes of the velocity fluctuation of the two phases, which are quantified by the standard deviations $\sigma _{v_y}$ and $\sigma _{u_y}$ (depicted in figure 18a), increase with the volume fraction. Since the velocity fluctuation of the particle motion can be decomposed to the contributions from the particle-sampled fluid flow and the local particle–fluid relative motion, i.e. $v_y'={U^L_{rel}}'+{u_y^{f@p}}'$, we have $\sigma _{v_y}^2 \approx \sigma _{u^L_{rel}}^2+ \sigma _{u_y^{f@p}}^2$ (with a negligible residual due to the crossing term of these two contributions) (Uhlmann & Doychev Reference Uhlmann and Doychev2014). Therefore, to gain further understanding, we also present the information about the second and third moments of $U^L_{rel}$ and $u_y^{f@p}$ in figure 18. As shown in figure 18(a), the contributions from $\sigma _{u_y^{f@p}}$ and $\sigma _{u^L_{rel}}$ to the particle velocity fluctuation are almost equal at $\phi =0.1\,\%$, 0.5 % and 10 %. However, interestingly, the former contribution dominates the latter one in the other three cases, and this tendency is most significant at $\phi =1\,\%$, corresponding to the strongest particle clustering. This observation was also discussed in Uhlmann & Doychev (Reference Uhlmann and Doychev2014) and can be attributed to the different local fluid velocity experienced by the particles in the clustering and void regions. Regarding the local particle–fluid relative motion, we find that $\sigma _{u^L_{rel}}$ increases monotonically as $\phi$ increases. This increase of fluctuation is presumably due to the influence of the greater turbulence intensity (manifested by the increase of $\sigma _{u^y}$) on the particle–fluid relative motion in dense suspensions.

Figure 17. The p.d.f. of the vertical component of (a) the particle velocity fluctuations, $v_y'=v_y-\langle v_y \rangle$, and (b) the fluid velocity fluctuations, $u_y'=u_y-\langle u_y \rangle _f$.

Figure 18. (a) Standard deviation (normalized by $V_t$) and (b) skewness of the vertical component of the particle velocity $v_y$, the fluid velocity $u_y$, the fluid velocity sampled by the particle $u^{f@p}_y$ and the relative velocity between the particle motion and the local fluid flow $U^L_{rel}$.

In addition, the asymmetric distribution of the particle and fluid velocity fluctuations shown in figure 17 is worthy of further discussion. First, the asymmetry of the p.d.f. of both $v_y'$ and $u_y'$, manifested by the higher negative tails, is appreciable in dilute cases. The skewed distribution of the fluid vertical velocity fluctuation, which has been reported in the flow induced by settling particles (Doychev Reference Doychev2014) or rising bubbles (Risso Reference Risso2018), is related to the dominant effect of wakes. With the increase of the volume fraction, the p.d.f.s of velocity fluctuations gradually recover to be symmetric, which is confirmed by the decrease of the skewness in magnitude shown in figure 18(b) and indicates the attenuated effect of the particle wakes. Second, we note from figure 18(b) that the skewness of the particle-sampled fluid velocity, $S_{u^{f@p}_y}$, is negligible compared with that of the local fluid–particle relative velocity, $S_{U^L_{rel}}$. Therefore, the asymmetric distribution of the particle vertical velocity at low volume fractions seems to be mainly related to the local particle–fluid motion. We speculate that the considerable variation of $S_{U^L_{rel}}$ with varying volume fraction may be caused by the change of the flow condition in the vicinity of the dispersed particles, which is worthy of further investigation.

3.3. Particle-particle collisions

We finally investigate the collision rate of settling particles in the present flow system. Here, we introduce the collision kernel $\varGamma$ to quantify the collision rate. The dynamic collision kernel, denoted by $\varGamma ^D$, is defined by (Wang et al. Reference Wang, Wexler and Zhou2000, Reference Wang, Ayala, Kasprzak and Grabowski2005)

(3.12)\begin{equation} \varGamma^D = \frac{2\dot{N}_C}{n^2}, \end{equation}

where $\dot {N}_C$ represents the number of collision events per unit volume per unit time and $n=N_p/V_{tot}$ is the number density of particles in the suspension. In the meantime, the collision kernel can also be described in a kinematic form (namely the kinematic collision kernel $\varGamma ^K$), i.e. the inward flux of particles across the surface of a sphere with a collision radius $R_{12}$, as (Wang et al. Reference Wang, Wexler and Zhou2000)

(3.13)\begin{equation} {\varGamma ^K} = 2{\rm \pi} R_{12}^2 \langle |W_r(R_{12})| \rangle g\left( {{R_{12}}} \right). \end{equation}

We notice that the particle pair statistics are involved in the above definition. Specifically, $\langle |W_r(R_{12})| \rangle$ represents the average absolute radial relative velocity (RRV) of particle pairs with a centre-to-centre distance $R_{12}$, and $g(R_{12})$ is the RDF evaluated at the collision radius $R_{12}$. It is important to note that the definition (3.13) is valid only under the flux-balance assumption (Wang et al. Reference Wang, Wexler and Zhou2000), which requires the equality between the inward and outward fluxes of the particle motion across the surface of the collision sphere. For clarity, we define the averaged magnitude of the negative and positive RRV as the inward and outward velocity $\langle W_r^ -\rangle$ and $\langle W_r^ +\rangle$, respectively. Then, the average absolute RRV can be expressed as

(3.14)\begin{equation} \langle |W_r| \rangle=P_{in}\langle W_r^ -\rangle +(1-P_{in})\langle W_r^ +\rangle, \end{equation}

where $P_{in}$ is the probability of observing negative the RRV (as has been shown in figure 14). Accordingly, the flux-balance assumption is valid when the ratio between the inward and outward flux, $C_p$, defined by (Wang et al. Reference Wang, Wexler and Zhou2000)

(3.15)\begin{equation} {C_p} = \frac{{{P_{in}}\langle W_r^ -\rangle }}{{\left( {1 - {P_{in}}} \right)\langle W_r^ +\rangle }}, \end{equation}

is equal to unity at the separation distance $r=R_{12}$.

In figure 19 we present the particle pair statistics as functions of the radial separation distance $r$. To begin with, figure 19(a) displays the relative radial velocities at different volume fractions, leading to the following key observations. First, for a large separation distance $r$, the relative radial velocities generally increase with the particle volume fraction $\phi$, which is ascribed to the intensified particle–fluid and particle–particle interactions as $\phi$ grows. However, the slight decrease of $\langle |W_r| \rangle$ from $\phi =5\,\%$ to $\phi =10\,\%$ may be related to the weakening effect of particle wakes, reminiscent of the decrease of the particle and fluid velocity fluctuations shown in figure 18(a). Second, a notable peak value of $\langle |W_r| \rangle$ is observed at the separation distance $r\approx 1.5D_{eq}$ in the most dilute case of $\phi =0.1\,\%$. This reflects the influence of wake-induced DKT-like interactions on the particle relative motion. However, this peak value diminishes as $\phi$ increases due to the disruption of particle wakes. Third, as the separation distance approaches $r=2b$, the relative radial velocities of particles decrease rapidly. The deceleration of the approaching/separating motions between nearby particles should be attributed to the lubrication effect. Last, the mean inward velocity $\langle W_r^ -\rangle$ and outward velocity $\langle W_r^ +\rangle$ exhibit slight discrepancies, indicating that the particle velocity field is compressible (Wang et al. Reference Wang, Wexler and Zhou2000). While, as shown in figure 19 (b), the ratio $C_p$ between the inward and outward fluxes is equal to unity (with some fluctuations due to the statistical error), irrespective of the volume fraction. This indicates that the flux-balance assumption remains to be valid in the present four-way coupling particle–fluid system, same as in the one-way coupling simulations of point-particle-laden turbulent flows (Wang et al. Reference Wang, Wexler and Zhou2000).

Figure 19. Particle pair statistics as functions of the separation distance $r$. Each line starts from $r=2b$ because of the geometric constraint for finite-size particles. The vertical black dotted line in each panel represents the separation distance $r=D_{eq}$. (a) Relative radial velocity. Dashed lines: average inward velocity $\langle W_r^- \rangle$; dotted lines: average outward velocity $\langle W_r^- \rangle$; solid lines: average absolute velocity $\langle |W_r| \rangle$. (b) The ratio between the inward and outward flux. (c) Average relative angle between particle pairs. The horizontal dash-dotted lines represent $\langle \theta _{rel} \rangle =45^\circ$ and $\langle \theta _{rel} \rangle =60^\circ$.

Moreover, we also illustrate the average relative angle $\langle \theta _{rel} \rangle$ between the symmetry axes of particle pairs with the varying separation distance in figure 19(c). Interestingly, $\langle \theta _{rel} \rangle$ seems to converge to approximately $\langle \theta _{rel} \rangle \approx 45^\circ - 50^\circ$ with the approaching of the particle–particle distance, irrespective of the volume fraction. This observation is quite different from the alignment of passive directors in the turbulence (Zhao et al. Reference Zhao, Gustavsson, Ni, Kramel, Voth, Andersson and Mehlig2019), and may be related to the hydrodynamic interactions between particle pairs (Ardekani et al. Reference Ardekani, Costa, Breugem and Brandt2016). Additionally, as the separation distance increases, the value of $\langle \theta _{rel} \rangle$ saturates to a certain value dependent on the particle volume fraction: the saturated value of $\langle \theta _{rel} \rangle$ is an increasing function of $\phi$, which can be explained by the orientational behaviour of settling particles. In the limiting case where all particles settle with the major axis perpendicular to gravity, the random orientation of particles on the two-dimensional horizontal plane will lead to $\langle \theta _{rel} \rangle =45^\circ$. The case with $\phi =0.1\,\%$ is close to this limit as the broad-side-on orientation dominates therein. While, in another limit where the particle orientations are totally random in the three-dimensional space, the average relative angle should be $\langle \theta _{rel} \rangle =60 ^\circ$. This interprets the increase of $\langle \theta _{rel} \rangle$ as $\phi$ increases, in consideration of the randomization of the particle orientation in the dense suspension (see figure 10a).

Now we return to the discussion on the particle collision kernel. When the flux-balance assumption is valid, $\varGamma ^K$ should be strictly equivalent to $\varGamma ^D$ for spherical particles with the collision radius being the diameter of the sphere (Wang et al. Reference Wang, Wexler and Zhou2000). However, for spheroidal particles, deriving the exact expression of the kinematic collision kernel $\varGamma ^K$ is theoretically challenging due to the complexity of their geometry. Alternatively, Siewert et al. (Reference Siewert, Kunnen and Schröder2014) proposed treating the expression (3.13) as an approximate kinematic collision kernel for spheroidal particles by using the equivalent diameter of the spheroid as the collision radius (i.e. $R_{12}=D_{eq}$).

In figure 20(a) we illustrate the dynamic and approximate kinematic collision kernel of settling prolate particles at different volume fractions. The results demonstrate that $\varGamma ^K$ underestimates the exact dynamic collision kernel $\varGamma ^D$ in the present flow system, similar to the case of settling spheroids in a quiescent fluid with the negligence of particle–particle interactions (Jiang et al. Reference Jiang, Huang, Xu and Zhao2024). Since the flux-balance assumption has been validated, the discrepancy between $\varGamma ^K$ and $\varGamma ^D$ should be ascribed to the abovementioned approximation regarding the geometry of the prolate spheroid when defining the approximate kinematic collision kernel. Even though, $\varGamma ^K$ still provides a reasonable estimation of the collision kernel, as it qualitatively captures the decreasing trend of $\varGamma ^D$ with the increasing volume fraction $\phi$.

Figure 20. (a) Dynamic and kinematic collision kernels and (b) the RDF and RRV evaluated at $r=D_{eq}$ at different particle volume fractions.

Furthermore, we depict the variation of $g(D_{eq})$ and $\langle | W_r (D_{eq}) | \rangle$, both of which contribute to the kinematic collision kernel $\varGamma ^K$ via (3.13), with the change of the volume fraction in figure 20(b). On the one hand, $g(D_{eq})$ decreases monotonically as $\phi$ increases, which plays an essential role in the reduction of $\varGamma ^K$. This highlights the significance of the wake-induced particle microstructures (as has been discussed in § 3.2) to the collision rate for settling particles. Accordingly, we remark that although the maximum value of RDF is not an appropriate criterion to quantify the particle clustering, the RDF is particularly relevant to the collision efficiency of particles as it directly contributes to the kinematic collision kernel. On the other hand, the average absolute RRV, $\langle | W_r (D_{eq}) | \rangle$, exhibits a minor degree of variation with the change of the volume fraction. Previous studies on settling spheroidal particles in turbulence, which ignored particle–particle interactions, reported the enhancement in the average RRV owing to the dispersion of the settling velocity of randomly oriented spheroidal particles (Siewert et al. Reference Siewert, Kunnen and Schröder2014; Jucha et al. Reference Jucha, Naso, Lévêque and Pumir2018). However, this mechanism does not apply to the present flow system, although particle orientations become more randomized as $\phi$ increases. With the fully resolved particle–particle hydrodynamic interactions herein, we attribute the abovementioned difference to the predominate effect of lubrication on the relative motion among nearby particles. Under the lubrication effect, the relative radial velocity of particle pairs reduces with the separation distance approaching $r=D_{eq}$ (see figure 19a), and is responsible for the nearly constant value of $\langle | W_r (D_{eq}) | \rangle$ for different volume fractions.