Prolate spheroids settling in a quiescent fluid: clustering, microstructures and collisions

In past decades, a series of experimental and numerical studies on the settling of an isolated sphere have been carried out. These earlier studies revealed that the dynamic mode and moving speed of a single falling/rising sphere are determined by two dimensionless parameters, the density ratio $\alpha$ and the Galileo number $Ga$ (Jenny et al., 2004; Horowitz & Williamson, 2010; Zhou & Dušek, 2015; Raaghav et al., 2022). The former measures the inertia of the solid particle, and the later quantifies the ratio between the buoyancy and viscous force acting on the sphere. In this study, our focus is on settling particles so the density ratio is greater than unity. One can also describe this problem by defining a dependent parameter, the Reynolds number $Re_{t}$, based on the a-posteriori settling velocity $V_{t}$ and the diameter of the sphere $D$. In the creeping flow regime ($Re_{t}\ll 1$), the sphere settles vertically with a constant settling velocity by balancing the buoyancy force with the Stokes drag. With the increase of the Reynolds number to a finite value, the introduction of fluid inertia breaks the fore-aft symmetry of the fluid flow around the settling sphere so the wake emerges. As $Re_{t}$ increases, the change of the rear-wake morphology can trigger the path instability of the settling sphere. A variety of settling modes, including vertical, oblique, zigzag, helical and chaotic motions, could be observed with varying inertia of the fluid and particle (Jenny et al., 2004; Horowitz & Williamson, 2010; Ern et al., 2012; Zhou & Dušek, 2015; Raaghav et al., 2022).

Concerning the settling motion of a pair of particles, the hydrodynamic interaction between them must be taken into account. A typical phenomenon is the drafting-kissing-tumbling (DKT) process of a pair of initially vertical-aligned settling spheres in the inertial-flow regime (Fortes et al., 1987; Glowinski et al., 2001). Specifically speaking, in the first stage (drafting stage), the trailing particle accelerates its settling motion as it resides in the wake of the leading particle. The two particles exhibit an attractive relative motion during this stage. Subsequently, the two particles touch and form an elongated body aligning along the vertical direction (kissing stage). However, settling with this configuration is unstable, so the particle pair tumbles and separates under the effect of hydrodynamic interaction (tumbling stage). In this stage, the two particles behave as if they repel each other, and the originally trailing particle becomes the leading one. Hence, the DKT process reflects the complicated hydrodynamic interaction between a pair of settling particles.

When considering the sedimentation of a group of particles, the most well-known phenomenon is the hindered settling motion of dispersed particles in suspensions (Richardson & Zaki, 1954), i.e. the reduction of the mean settling velocity of particles with the increase of the volume fraction $\phi$. The physical explanation of this hindrance effect is as follows. To maintain zero net flux of the whole flow system, a mean upward fluid flow is generated to counteract the downward flux of settling particles. As a result, the upward mean flow increases the hydrodynamic drag acting on the particle phase and thus reduces the mean settling velocity (Di Felice, 1999). An empirical formula of the hindered settling velocity of particles in the creeping-flow regime was also proposed by Richardson & Zaki (1954). More recently, the hindered settling velocity of particles with finite fluid inertia was also observed in the numerical studies by means of the two-way coupling point-particle simulation (Climent & Maxey, 2003) or the particle-resolved direct numerical simulation (PR-DNS) (Yin & Koch, 2007; Zaidi et al., 2014). The empirical expression of the hindered settling velocity has also been improved in a series of studies to incorporate the finite-Reynolds-number correction (Garside & Al-Dibouni, 1977; Di Felice, 1999; Yin & Koch, 2007).

However, over the past two decades, a striking enhancement of the mean particle settling velocity has been observed under certain conditions, thanks to the state-of-the-art PR-DNS of the particle sedimentation (Kajishima & Takiguchi, 2002; Kajishima, 2004; Uhlmann & Doychev, 2014; Zaidi et al., 2014). Investigations into the particle spatial distribution have shown that the enhanced settling velocity is always associated with the formation of column-like particle clusters (Doychev, 2014; Uhlmann & Doychev, 2014; Zaidi et al., 2014). Later on, the experimental work conducted by Huisman et al. (2016) also confirmed these numerical observations. Zaidi et al. (2014) and Moriche et al. (2023) attributed the formation of particle clusters to the DKT-like interactions among settling particles, but this phenomenon can only be observed in dilute suspensions when the Reynolds number $Re_{t}$ is sufficiently high. Conversely, at low Reynolds numbers, the weaker wake-induced attraction among particles results in orderly particle arrangements rather than particle clustering (Yin & Koch, 2007; Zaidi et al., 2014, 2015). While, in dense suspensions, the short distances between particles disrupt particle wakes and thus inhibit the formation of particle clusters as well (Zaidi, 2018b). Readers can refer to Chouippe et al. (2023) for the review of previous studies on this problem.

In practice, the shape of dispersed particles is commonly non-spherical. For instance, ice crystals in clouds, plankton in the marine, and dusts in the atmosphere are usually disk-like or rod-like in shape (Shaw, 2003; Mallios et al., 2020; Slomka & Stocker, 2020). For simplicity, non-spherical particles are often modelled by smooth-surface prolate/oblate spheroids or by polyhedrons with edges. Compared to spherical particles, the orientational behavior and rotational motion of non-spherical particles add complexity to their dynamics in fluid flows (Rahmani & Wachs, 2014; Voth & Soldati, 2017).

As for a single spheroid settling in a quiescent fluid, the particle would maintain its initial orientation in the creeping flow owing to the vanishing hydrodynamic torque. Thus, the settling velocity of the spheroid is orientation-dependent in this regime (Happel & Brenner, 1983). However, when the fluid inertia is taken into account, a non-negligible hydrodynamic torque reorients the settling spheroid to a broad-side-on alignment (Khayat & Cox, 1989; Ardekani et al., 2016; Dabade et al., 2016). Additionally, when the fluid inertia is strong enough to trigger wake instability, the settling motion of the spheroid transitions from the steady vertical falling to complicated unsteady modes, similar to the case of a settling sphere. The velocity and mode (including spiral, zigzag/fluttering, tumbling, and chaotic) of the settling motion are jointly determined by the density ratio, Galileo number and the shape of the spheroid (Chrust et al., 2013; Ardekani et al., 2016; Zhou et al., 2017; Moriche et al., 2021), and can also be altered by the presence of walls (Huang et al., 2014; Yang et al., 2015).

Moreover, according to the numerical work by Ardekani et al. (2016), the DKT process between a pair of settling spheroids is quite different from that of spherical particles. As for a pair of oblate particles with an aspect ratio (the ratio between the polar and equator radius) $\lambda=1/3$, the two particles do not undergo the tumbling stage after they approach and touch. Instead, they fall with a steady pilled-up configuration, as if they are stuck together (Ardekani et al., 2016). Regarding a pair of prolate spheroids with $\lambda=3$, the DKT process is more complicated and dependent on their initial relative angle (Ardekani et al., 2016). If the symmetry axes of the two prolate particles are initially parallel, the DKT process is similar to that of spherical particles. While, if the symmetry axes are perpendicular at the beginning, a stable cross-like configuration is formed and the two prolate particles do not separate for a long time after they touch, similar to the case of oblate particles. In principle, the attraction zone (within which the trailing particle can be attracted by the leading one) is larger, and the interaction time (the time duration for the particle pair to keep in touch) is longer for the DKT process of spheroidal particle pairs, compared to that of spherical ones (Ardekani et al., 2016; Moriche et al., 2023). To conclude, the particle shape plays an important role in the hydrodynamic interaction between settling particles.

As regards the settling of a large number of non-spherical particles, an increased settling velocity of elongated fibres was observed in the creeping flow regime due to the formation of particle streamers aligning in the gravitational direction (Kuusela et al., 2003; Saintillan et al., 2006; Shin et al., 2009). In the finite-fluid-inertia regime, Seyed-Ahmadi & Wachs (2021) numerically studied the settling motion of cubic particles. In contrast to the clustered settling spheres at $Ga=160$ and $\phi=1\%$, the spatial distribution of settling cubes is closer to a random distribution in the same parameter setup, which was attributed to the greater rotational rate of settling cubes. Fornari et al. (2018) simulated the sedimentation of oblate spheroids with $\lambda=1/3$ and $Ga=60$ at different particle volume fractions. They reported appreciable particle clustering for settling oblate particles at a relatively low Reynolds number ($Re_{t}=38.7$), and considerable enhancement of the mean settling velocity up to $\langle V_{s}\rangle\approx 1.33V_{t}$ at $\phi=0.5\%$. Similar results were also reported in the recent work by Moriche et al. (2023), who considered the low-aspect-ratio oblate spheroids with $\lambda=2/3$ and higher Galileo number with $Ga=111$ and 152. As for the case of prolate particles, Lu et al. (2023) simulated the settling of prolate spheroids with $\lambda=2$ and $Ga=41.8$ at $\phi=2.2\%,5.5\%$ and $9.9\%$ using a relatively small periodic computational domain. They reported a decreased mean particle settling velocity and a transition from the hydrodynamic-interaction-dominated regime to the particle-collision-dominated regime with the increase of $\phi$.

The collision rate among dispersed particles plays an important role in the particle coagulation in fluid flows, which are relevant to many industrial and natural processes. In the past, a plenty of work has been carried out to study the collision and coagulation of point-like spherical particles in turbulent flows in the framework of one-way coupling approach (Saffman & Turner, 1956; Sundaram & Collins, 1997; Wang et al., 2000; Ayala et al., 2008). Readers can refer to the reviews by Grabowski & Wang (2013) and Pumir & Wilkinson (2016) for more details. Recently, some researchers extended the work to non-spherical particles, and demonstrated that the orienational behavior of elongated or flattened particles enhances their collision rate in turbulence (Jucha et al., 2018; Siewert et al., 2014; Slomka & Stocker, 2020; Arguedas-Leiva et al., 2022; Grujić et al., 2024). However, when further considering the intricate particle-fluid and particle-particle interactions, the understanding of particle collision rates remains limited. In Wang et al. (2005), the hydrodynamic interactions among particles were addressed by adding the particle-induced disturbance into the background turbulence. These disturbance can either augment or attenuate collision rate, depending on whether they act as the far-field or near-field influence. In the framework of PR-DNS, Chen et al. (2020) simulated the transport of spherical particles in the homogeneous isotropic turbulence, and studied the collision rate of bidispersed inertial particles. Furthermore, Fornari et al. (2019) simulated settling spherical particles in both the quiescent fluid and the turbulent environment, and examined the effect of particle spatial distribution and relative motion on the collision rate. However, to the best of our knowledge, the collision rate of settling non-spherical particles with the full consideration of fluid-particle and particle-particle interactions has not been investigated so far.

According to the above literature review, we are still far from achieving a comprehensive understanding of settling non-spherical particles in suspensions. In the present work, we investigate the sedimentation of prolate particles in an initially quiescent fluid by means of PR-DNS. In particular, considering the significant impact of hydrodynamic interactions on the dynamics of pairwise settling prolate particles (Ardekani et al., 2016), we aim to explore the collective behavior of settling prolate particles at varying volume fractions. This work is motivated by two concerns. First, although the enhancement of particle settling velocity and particle clustering down to $\phi\approx 0.02\%$ have been reported for spherical particles (Doychev, 2014; Huisman et al., 2016), to the best of the authors’ knowledge, little is known about the sedimentation of non-spherical particles with $\phi<0.5\%$ in the inertial-flow regime. Hence, we study the settling motion of prolate spheroids within a wide range of the volume fraction from $\phi=0.1\%$ to $10\%$, with the particle aspect ratio and Galileo number fixed at $\lambda=3$ and $Ga=80$, following the study of a single and a pair of settling prolate particles by Ardekani et al. (2016). The particle-fluid density ratio is set as $\alpha=2$, which is a typical value for a solid-liquid system (Seyed-Ahmadi & Wachs, 2021). Interestingly, we observe a non-monotonic variation of the mean particle settling velocity as $\phi$ increases, so we further investigate the influences of particle clustering, hindrance effect and particle orientation on the settling speed of prolate spheroids. Second, the collision rate among settling non-spherical particles with finite sizes is not well understood so far. Therefore, we also investigate on this issue and scrutinize the particle pair statistics that are essential in determining the particle collision rate.

The remainder of this paper is organized as follows. In section 2, we describe the physical problem and the simulation set-ups of this study. Then, in section 3, we analyze the statistics of particle motions and spatial distributions, followed by the examination of the collision rate of dispersed particles at different particle volume fractions. Finally, we summarize the findings and draw conclusions in section 4.

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}\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\leq 2\%$, with a peak value of $\langle V_{s}\rangle\approx 1.25V_{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.

According to the discussion in section 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 (2007), Zaidi et al. (2015) and Fornari et al. (2018). By definition, $P(\boldsymbol{r})>1$ indicates a higher probability of finding a pair of particles with a separation of $\boldsymbol{r}$ compared to 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.:

$$P\left({r,\varphi}\right)=\frac{1}{{2\pi}}\int_{0}^{2\pi}{P\left({\boldsymbol{r}}\right)d\theta}=\frac{1}{{2\pi}}\int_{0}^{2\pi}{P\left({r,\theta,\varphi}\right)d\theta}.$$

(6)

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\left({r,\varphi}\right)$ 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.:

$$W_{r}^{\left\{{i,j}\right\}}=\frac{{\left({{{\boldsymbol{v}}_{i}}-{{\boldsymbol{v}}_{j}}}\right)\cdot\left({{{\boldsymbol{x}}_{i}}-{{\boldsymbol{x}}_{j}}}\right)}}{{\left\|{{{\boldsymbol{x}}_{i}}-{{\boldsymbol{x}}_{j}}}\right\|}}<0,$$

(7)

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, 2007; Zaidi, 2018b; Seyed-Ahmadi & Wachs, 2021). 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., 2016). Regarding spherical particles, the shear-induced lift force dominates in wake regions, pushing trailing particles outside the wake (Yin & Koch, 2007). 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, 2019, 2021). 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\geq 5\%$.

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, 2007):

$$g\left(r\right)=\frac{1}{2}\int_{0}^{\pi}{P\left({r,\varphi}\right)\sin\varphi d\varphi}.$$

(8)

As shown in figure 15 (a), the spatial correlation of particle pairs is considerable increased with $g(r)>1$ for 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\leq 2\%$, is a monotonically decreasing function of the volume fraction. While, in dense suspensions with $\phi\geq 5\%$, the RDF is close to unity for all separation distance, corresponding to the attenuated particle microstructures in these cases.

In previous studies on settling particles, the degree of particle clustering is usually quantified by the maximum value of RDF (Zaidi et al., 2014; Fornari et al., 2018; Zaidi, 2018a). 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 3 (b)). 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 particles 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 4 (g,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, 2007; Fornari et al., 2018):

$$\left\langle{{P_{2}}}\right\rangle(r)=\frac{{\int_{0}^{\pi}P(r,\varphi){P_{2}}(\cos\varphi)\sin\varphi{\rm{d}}\varphi}}{{\int_{0}^{\pi}P(r,\varphi)\sin\varphi{\rm{d}}\varphi}},$$

(9)

in which ${P_{2}}\left({\cos\varphi}\right)=\left({3{{\cos}^{2}}\varphi-1}\right)/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, 2007). 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 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\geq 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 10D_{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:

	$\displaystyle\langle\zeta\rangle_{y}\left({x,z}\right)$	$\displaystyle=\frac{1}{{{L_{y}}}}\int_{0}^{{L_{y}}}{\zeta\left({x,y,z}\right)dy},$		(10)
	$\displaystyle\langle{{u_{y}}}\rangle_{y}\left({x,z}\right)$	$\displaystyle=\frac{1}{{{L_{y}}}}\int_{0}^{{L_{y}}}{\left[{1-\zeta\left({x,y,z}\right)}\right]{u_{y}}\left({x,y,z}\right)dy}.$		(11)

Here, $\zeta(x,y,z)$ is an indicator function which 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 above-mentioned correlation is weakened, owing to the diminishing particle microstructures in dense suspensions.

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 probability density functions 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 18(a)), 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}^{\prime}={U^{L}_{rel}}^{\prime}+{u_{y}^{f@p}}^{\prime}$, 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, 2014). 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 (2014) and can be attributed to the different local fluid velocity experienced by the particles in the clustering and void regions. Regarding the 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.

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}^{\prime}$ and $u_{y}^{\prime}$, 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, 2014) or rising bubbles (Risso, 2018), is related to 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 notice from figure 18(b) that the skewness of the particle-sampled fluid velocity, $S_{u^{f@p}_{y}}$, is negligible compared to 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.

We finally investigate the collision rate of settling particles in the present flow system. Here, we introduce the collision kernel $\Gamma$ to quantify the collision rate. The dynamic collision kernel, denoted by $\Gamma^{D}$, is defined by (Wang et al., 2000, 2005):

$$\Gamma^{D}=\frac{2\dot{N}_{C}}{n^{2}},$$

(12)

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 $\Gamma^{K}$), i.e. the inward flux of particles across the surface of a sphere with a collision radius $R_{12}$, as (Wang et al., 2000):

$${\Gamma^{K}}=2\pi R_{12}^{2}\langle|W_{r}(R_{12})|\rangle g\left({{R_{12}}}\right).$$

(13)

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 center-to-center 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 (13) is valid only under the flux-balance assumption (Wang et al., 2000), 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 radial relative velocity 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:

$$\langle|W_{r}|\rangle=P_{in}\langle W_{r}^{-}\rangle+(1-P_{in})\langle W_{r}^{+}\rangle,$$

(14)

where $P_{in}$ is the probability of observing negative the radial relative velocity (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., 2000):

$${C_{p}}=\frac{{{P_{in}}\langle W_{r}^{-}\rangle}}{{\left({1-{P_{in}}}\right)\langle W_{r}^{+}\rangle}},$$

(15)

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., 2000). 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., 2000).

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., 2019), and may be related to the hydrodynamic interactions between particle pairs (Ardekani et al., 2016). 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 behavior of settling particles. In the limiting case where all particles settle with the major axis perpendicular the 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 10 (a)).

Now we return to the discussion on the particle collision kernel. When the flux-balance assumption is valid, $\Gamma^{K}$ should be strictly equivalent to $\Gamma^{D}$ for spherical particles with the collision radius being the diameter of the sphere (Wang et al., 2000). However, for spheroidal particles, deriving the exact expression of the kinematic collision kernel $\Gamma^{K}$ is theoretically challenging due to the complexity of their geometry. Alternatively, Siewert et al. (2014) proposed treating the expression (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 $\Gamma^{K}$ underestimates the exact dynamic collision kernel $\Gamma^{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., 2024). Since the flux-balance assumption has been validated, the discrepancy between $\Gamma^{K}$ and $\Gamma^{D}$ should be ascribed to the above-mentioned approximation regarding the geometry of the prolate spheroid when defining the approximate kinematic collision kernel. Even though, $\Gamma^{K}$ still provides a reasonable estimation of the collision kernel, as it qualitatively captures the decreasing trend of $\Gamma^{D}$ with the increasing volume fraction $\phi$.

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 $\Gamma^{K}$ via equation (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 $\Gamma^{K}$. This highlights the significance of the wake-induced particle microstructures (as has been discussed in section 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., 2014; Jucha et al., 2018). 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 19 (a)), and is responsible for the nearly constant value of $\langle|W_{r}(D_{eq})|\rangle$ for different volume fractions.