3.1. Sedimentation of two spheres

Figure 2.

(a) Schematic representation of two spheres falling freely in an infinite domain, with the spheres initially separated either horizontally or vertically. The separation distance $r$ refers to the distance between their centres. (b) Comparison of the settling velocities for a pair of equal spheres separated either horizontally or vertically with the results of Maxey & Patel (Reference Maxey and Patel2001). Here, $U_{{sed}}$ and $U_0$ are the settling velocities for a pair of identical spheres and an isolated sphere, respectively.

We begin our investigation by validating the steady motion of a pair of identical and equal sized spheres settling under gravity in an unbounded domain (i.e. $L \rightarrow \infty$ ). The initial configuration of two spheres is shown in figure 2(a). The spheres are separated either horizontally or vertically to the line of the centres. As pointed out by Batchelor (Reference Batchelor1972), the settling velocity depends on the separation between the particles and the orientation of the line of centres to the vertical. Later, Maxey & Patel (Reference Maxey and Patel2001) introduced a finite-valued force multipole expansion to describe the dynamics of spherical particles sedimenting in Stokes flow. They used the force coupling method to evaluate the settling velocity of two equal spheres in the above-mentioned configurations. If the particles are otherwise isolated and not influenced by other particles or boundaries, then they will have equal settling velocities due to the symmetry of the induced Stokes flow. The mutual influence of the disturbance flow created by each particle leads to a settling velocity larger than that for a single isolated particle. Our numerical result in figure 2(b) shows good agreement with the results of Maxey & Patel (Reference Maxey and Patel2001). The settling velocity $U_{{sed}}$ is scaled by the Stokes velocity $U_0$ of a single sphere, and given in terms of the distance between the centres, so $r/R=2$ corresponds to touching spheres. There is no external torque on the particles and they are free to rotate in response to the local fluid velocity. Figure 2(b) suggests that the particles settle with the same velocity as an isolated sphere in an infinite medium, when they are very far apart from each other (i.e. $r/R \to \infty$ ). However, the mutual influence of the flow disturbance created by each particle results in a settling velocity that is larger than that of a single isolated particle. This effect leads to a velocity ratio $U_{{sed}}/U_0 \approx 1.41$ at $r/R = 2.01$ for horizontally separated spheres.

3.2. Effect of domain size

Figure 3.

Average settling velocity of spherical particles versus the number of particles $M$ at two different solid volume fractions: (a) $\phi = 0.025$ and (b) $\phi = 0.36$ . The system size $L$ is changed according to $M$ to keep the volume fraction constant, the settling velocity is normalized by $U_0$ , and error bars represent the standard deviation of the mean settling velocity. The red symbols represent the value $M=125$ , which we use in the rest of our simulations.

The system size determines how well the system approximates an infinite suspension and whether artificial effects distort the settling dynamics (Climent & Maxey Reference Climent and Maxey2003; Lu et al. Reference Lu, Xu, Zhong, Ni and Tryggvason2023). To ensure that the system size is sufficiently large, we study its influence on the average sedimenting velocity of spherical particles in two cases: a dilute suspension ( $\phi = 0.025$ ) and a moderately dense suspension ( $\phi = 0.36$ ). In both cases, we vary the number of particles $M$ and the system size $L$ , while keeping the volume fraction constant, and measure the average sedimenting velocity. The results indicate that in both cases, a small box size significantly affects the mean settling velocity; see figures 3(a,b). However, as the domain size increases, it soon reaches the asymptotic value. In the rest of the work, we use $M=125$ so the finite size effects are small while the computational cost is still moderate. The red symbols in figures 3(a,b) indicate the selected box size used in our simulations for different solid volume fractions.

3.3. Sedimentation velocity of a monodisperse suspension of spheres and cubes

Figure 4.

Dimensionless mean sedimentation velocity of monodisperse spheres (upper row) and cubes (lower row) at various volume fractions. The error bars indicate the standard error of the mean settling velocity for several independent random configurations.

The microstructure of suspension evolves in time during sedimentation from an initial equilibrium distribution to a new steady-state distribution. Therefore the sedimentation velocity

(3.1) \begin{align} U_{{sed}}(t)= \left\langle \frac {1}{N} \sum _{i=1}^{N} U_i(t) \right \rangle, \end{align}

ensemble averaged over the initial configurations, is time-dependent. To obtain the mean sedimentation velocity, an ensemble average (denoted by $\langle \cdot \rangle$ ) of different initial configurations of the particles is considered. These configurations were randomly sampled from an equilibrium hard-particles distribution function by the Monte Carlo method (Metropolis et al. Reference Metropolis, Rosenbluth, Rosenbluth, Teller and Teller1953). The initial equilibrium averages were taken over $15$ independent initial configurations. The time evolution of $U_{{sed}}(t)$ for randomly distributed suspension of spheres and cubes, at various solid volume fractions from equilibrium to steady state, is shown in figure 4. The time-dependent settling velocity is normalized by the settling velocity of an isolated particle $U_0$ (sphere/cube for corresponding cases) falling under gravity, whereas the time is normalized with $t_s=R/U_0$ . The characteristic length scale $R$ is already mentioned in table 1. The settling velocity of a monodisperse suspension of spheres and cubes fluctuates slightly as time progresses. Figure 4 shows that the mean sedimentation velocity of a suspension of spheres is always slightly higher at all volume fractions compared to cubic suspension. As reported in the study of Seyed-Ahmadi & Wachs (Reference Seyed-Ahmadi and Wachs2019, Reference Seyed-Ahmadi and Wachs2021), anisotropic particles (i.e. cubes, tetrahedrons, etc.) always produce higher rotational displacements due to their sharp edges and rapid variation of hydrodynamic loads. Thus cubic suspension causes a strong upward flux to hinder the settling process as compared to spheres. Figure 4 also shows the deviation of the settling velocity of a random initial suspension from its mean as time progresses (indicated by error bars).

During the settling process the particles take some time to rearrange their configuration, which affects the sedimentation velocity. This relaxation time and the dispersive motion of the particles have been studied by several researchers. The study of Kuusela (Reference Kuusela2005) demonstrated that the relaxation time in the vertical direction for a sphere suspension follows approximately a $\phi ^{-1/2}$ relaxation in the full range of volume fractions. A sharp decrease in the relaxation time is also observed in the investigation of Hamid et al. (Reference Hamid, Molina and Yamamoto2013) at $\phi \gt 0.3$ in the direction of gravity, mainly because of the ordering of the particles, whereas a jump in relaxation time at $\phi \approx 0.5$ could occur because of very-long-range ordering of the particles. At moderate range of volume fractions (i.e. $ 0.1 \leqslant \phi \leqslant 0.4$ ), the intensity of hydrodynamic interactions between particles increases, which in turn enhances the velocity fluctuations in the directions parallel and perpendicular to gravity. As a consequence, a significant time is required for the particles to relax in order to reach an equilibrium state. With the additional increment in $\phi$ , the settling process is strongly hindered, and the settling particles produce a weaker upward flux of the surrounding fluid in a tightly packed space. The velocity fluctuations in both directions – and hence the relaxation time – diminish drastically. As a result, both the suspensions quickly reach the steady state at $\phi \simeq 0.5$ . A detailed analysis of the velocity fluctuations for both suspensions will be discussed in the next subsection.

Figure 5.

Average settling velocity $U_{{sed}}$ of the particles normalized by the Stokes velocity $U_0$ as a function of volume fraction ( $\phi$ ) to highlight the hindered settling of the particles. Solid lines denote our simulation result; dashed lines show various theoretical predictions of Batchelor (Reference Batchelor1972), Hayakawa & Ichiki (Reference Hayakawa and Ichiki1995), Richardson & Zaki (Reference Richardson and Zaki1997), Di Felice (Reference Di Felice1999), and the simulation results of Yamamoto et al. (Reference Yamamoto, Molina and Nakayama2021); the symbols represent the experimental data of Nicolai et al. (Reference Nicolai and Guazzelli1995). Here, the error bars are smaller than the symbols’ size.

In figure 5, we present the normalized average settling velocity $U_{{sed}}/U_0$ as a function of the solid volume fraction, where $U_{{sed}}$ and $U_0$ denote the ensemble average of the suspension settling velocity and the terminal settling velocity of an isolated sphere, respectively. We also compare our computational result for spheres with the existing theoretical predictions (Batchelor Reference Batchelor1972; Hayakawa & Ichiki Reference Hayakawa and Ichiki1995; Richardson & Zaki Reference Richardson and Zaki1997; Di Felice Reference Di Felice1999), experimental study (Nicolai et al. Reference Nicolai and Guazzelli1995) and simulation data (Hamid et al. Reference Hamid, Molina and Yamamoto2013; Yamamoto et al. Reference Yamamoto, Molina and Nakayama2021), which involves many-body hydrodynamic interactions. We observe that the sedimentation velocity decreases with increasing volume fraction. Since both the particles and the fluid are incompressible as particles sediment, an upward flow to balance the flow of mass is generated that hinders the sedimentation velocity. This effect of particle volume fraction on the average sedimentation velocity was first introduced by Batchelor (Reference Batchelor1972). The theoretical prediction of Batchelor assumes a uniform distribution in the separation of pairs of spheres, ignoring the fluid back flow effects, given by $U_{{sed}}/U_0=1-6.55 \phi$ , which is valid only for low volume fractions, where long-range hydrodynamic interactions can be neglected. We find good accuracy in the sedimentation velocity for spheres with Batchelor’s correlation in the range $\phi \leqslant 0.04$ . For larger volume fractions, our simulation result overestimates Batchelor’s model. Hayakawa & Ichiki (Reference Hayakawa and Ichiki1995) also proposed the normalized sedimentation rate in the context of a resistance problem, given by $ {U_{{sed}}}/{U_0} = {(1-\phi )^3}/{ 1+2 \phi + 1.492 \phi (1-\phi )^3 }$ . Our current result also agrees well with the Hayakawa–Ichiki model for $\phi \leqslant 0.05$ . The simplest semi-empirical relation for a hindered settling function has been proposed by Richardson & Zaki (Reference Richardson and Zaki1997) as $U_{{sed}}/U_0=(1- \phi )^n$ , which gives good predictions at higher $\phi$ . Several studies have found different exponents ranging from $n=4.7$ to $n=6.55$ for the Stokes flow, whereas in our study, the best suited exponent is $n=5.9$ for spheres, and $n=6.2$ for cubes. Finally, we observe that the cubes always sediment slower than the sphere suspensions for all volume fractions, except possibly at the very dilute regime.

3.4. Velocity fluctuations

Figure 6.

Dimensionless hydrodynamic velocity fluctuations of particles as a function of volume fraction ( $\phi$ ) are shown for (a) the settling direction ( $\sigma _{\parallel }$ ) and (b) the direction perpendicular to settling ( $\sigma _{\perp }$ ). The insets illustrate the dependence of velocity fluctuations on $\phi ^{1/2}$ for $\phi \leqslant 0.1$ . (c) The anisotropy in velocity fluctuations, $\Gamma = \sigma _{\parallel } / \sigma _{\perp }$ , is plotted as a function of volume fraction ( $\phi$ ).

During settling, particle velocities fluctuate both vertically and horizontally due to inter-particle hydrodynamic interactions. Vertical velocity fluctuations reflect temporal deviations from the mean sedimentation velocity, while horizontal fluctuations, occurring perpendicular to gravity, are characterized by the mean velocities in the $XY$ -plane. The time-averaged standard deviations of particle velocities parallel and perpendicular to gravity are quantified as $\Delta U_z = \sqrt { \langle [U_{iz} - U_{{sed}}]^2 \rangle }$ and $\Delta U_{xy} = \sqrt { \langle U_{ix}^2 + U_{iy}^2 \rangle }$ , respectively. Here, $U_{iz}$ represents the temporal velocity of the $i$ th particle along the gravitational direction, while $U_{ix}$ and $U_{iy}$ are the velocities in the horizontal directions. These standard deviation values are normalized by the corresponding settling velocity of each suspension, commonly referred to as ‘relative fluctuation’ (Nicolai et al. Reference Nicolai and Guazzelli1995; Climent & Maxey Reference Climent and Maxey2003). Figures 6(a) and 6(b) display dimensionless velocity fluctuations of particles in the settling direction ( $\sigma _{\parallel } = \Delta U_z / U_{{sed}}$ ) and in the perpendicular direction ( $\sigma _{\perp } = \Delta U_{xy} / U_{{sed}}$ ), respectively. The fluctuations exhibit similar qualitative behaviour in both directions. At low volume fractions ( $\phi \lt 0.04$ ), $\sigma _{\parallel }$ and $\sigma _{\perp }$ display characteristics of the Stokes regime, with a $\phi ^{1/2}$ dependence (see insets of figures 6 a,b). As the solid volume fraction increases, velocity fluctuations rise gradually up to a moderate range ( $0.04 \leqslant \phi \leqslant 0.2$ ). Beyond this range ( $\phi \gt 0.2$ ), a gradual decline in fluctuations is observed for both directions. Additionally, velocity fluctuations are slightly higher for cubes compared to spheres, with the differences between the two shapes becoming more pronounced as $\phi$ increases. As the volume fraction increases within the range ( $0.04 \leqslant \phi \leqslant 0.2$ ), particles are more crowded and collision events, where steric interactions push the particles apart, become more common. Consequently, momentum is significantly transferred from the axial direction to the transverse direction, which amplifies velocity fluctuations in the transverse direction (Yao et al. Reference Yao, Criddle and Fringer2021). Notably, even a slight reorientation of a cube can generate considerable hydrodynamic torques, leading to rotations that further enhance this effect (Seyed-Ahmadi & Wachs Reference Seyed-Ahmadi and Wachs2019). As a consequence, cubes seem to experience higher transverse velocity fluctuations compared to spheres, thus facilitating a more efficient transfer of momentum from the axial to the transverse direction. With further increment in $\phi$ , the average inter-particle distance, which is inversely proportional to $\phi$ , diminishes, causing particles to become more densely packed (Zaidi Reference Zaidi2018). Beyond $\phi \gt 0.2$ , this dense packing restricts particle motion, resulting in a reduction in velocity fluctuations. In this range, the reduced inter-particle spacing disrupts the fluid structure required for particle clustering, while the frequency of hydrodynamic collisions rises, further suppressing velocity fluctuations.

The increase in velocity fluctuations with rising solid volume fraction is closely associated with the degree of anisotropy in these fluctuations, for both cubic and spherical suspensions. This anisotropy, defined as the ratio of vertical to horizontal velocity fluctuations ( $\Gamma = \sigma _{\parallel } / \sigma _{\perp }$ ), is shown in figure 6(c). In highly dilute conditions, anisotropy is notably high for both particle shapes, primarily due to restricted lateral movement in this regime. For dilute suspensions ( $\phi \leqslant 0.05$ ), $\Gamma$ ranges between 1.7 and 2.1 for both cases, aligning with the experimental findings of Nicolai & Guazzelli (Reference Nicolai and Guazzelli1995) on non-Brownian spheres under very dilute conditions. As the solid volume fraction increases to moderate levels ( $\phi \approx 0.2$ ), velocity fluctuations intensify in both the vertical and horizontal directions, leading to a substantial reduction in anisotropy. In this range, the anisotropy decreases compared to the dilute regime, with $\Gamma$ approaching approximately $1.3$ for cubic suspensions. A similar trend was also observed in the studies of Nicolai & Guazzelli (Reference Nicolai and Guazzelli1995) and Ladd (Reference Ladd1996). In the study of Seyed-Ahmadi & Wachs (Reference Seyed-Ahmadi and Wachs2021) at finite Reynolds numbers, it was observed that higher rotational velocities of cubes contribute to the transfer of momentum and energy from the gravity direction to the transverse direction, which led to a lower velocity fluctuations anisotropy in cubic suspensions than in spherical suspensions. Here, we observe the same effect in the zero Reynolds number regime. As $\phi$ increases further to $0.5$ , the anisotropy continues to decline due to more frequent and multidirectional inter-particle collisions.

3.5. Probability density function of particle sedimentation velocity

Figure 7.

The PDF of vertical velocity is presented for suspensions of (a) spherical and (b) cubic particles across different solid volume fractions. Negative values of $U_i/U_{{sed}}$ indicate that particles are moving vertically upwards, against the direction of gravity. (c) The percentage of particles moving opposite to gravity at any given moment, as a function of solid volume fraction.

To better understand the differences in velocity fluctuations between spherical and cubic particles, we investigate the instantaneous particle velocity during sedimentation. Figures 7(a) and 7(b) illustrate the normalized probability density function (PDF) of particle velocity as a function of bulk concentration. Notably, $U_i/U_{{sed}} \gt 1$ indicates that particles are falling faster than the average sedimentation velocity. Positive values of $U_i/U_{{sed}}$ imply downward movement (i.e. in the direction of gravity), while negative values represent upward movement. Note that as the particles settle, the surrounding fluid generates an upward jet, due to the conservation of mass, that can drag some particles with it. Interestingly, the variance of the velocity in the vertical direction is a non-monotonic function of the volume fraction. In a very dilute regime ( $\phi \approx 0.001$ ), both the spherical and cubic suspensions exhibit a sharp peak, indicating that the velocity variance is small. In this range, the larger inter-particle separation distance results in weaker hydrodynamic interactions, thus deviation from the single-particle dynamics is small. Additionally, the percentage of upward-moving particles ( $N_{U_z}$ ) is calculated by integrating the area under the PDF of the settling velocity. Figure 7(c) illustrates that, in the highly dilute regime, the fraction of upward-moving particles approaches nearly zero for both types of particles.

As the volume fraction increases, velocity fluctuations increase, causing particles to disperse more and collide more frequently. At moderate concentrations, hydrodynamic interactions among nearby settling particles and inter-particle collisions create a strong upward jet, which can lift some particles (both spheres and cubes) along with it. The frequent hydrodynamic interactions and particle rebounds after collisions may also lead to individual particles moving upwards. A closer examination of the PDF of vertical particle velocity reveals that the tails become stretched as $\phi$ increases from $0.025$ to $ 0.125$ , with the left tail slightly longer than the right. Therefore, the fraction of upward-moving particles increases with rising bulk concentration in the low to moderate range. Results in figure 7(c) indicate that in cubic suspensions, approximately $5.3 \pm 0.2 \,\%$ of particles move against the settling direction, compared to only $4.6 \pm 0.1 \,\%$ for spherical particles at $\phi \approx 0.1$ . In this range, a closer inspection of figures 7(a,b) reveals that the peak of the PDF is slightly higher for spherical particles, indicating that a greater number of them attain velocities close to the mean sedimentation velocity. In contrast, cubic particles exhibit a slightly stretched left tail in the negative region, showing a broader distribution of lower velocities; see also supplementary figure 1. The sharp edges and rapid changes in hydrodynamic forces due to variations in orientation make cubes more susceptible to strong translational and rotational displacements, as well as velocity fluctuations. Additionally, cubes tend to acquire slightly higher angular velocities than spheres, resulting in significant rotation-induced forces perpendicular to their direction of movement. Consequently, cubic suspensions generate larger rotation-induced forces due to the higher magnitude of angular rotation, resulting in a greater fraction of particles being entrained by the upward jet compared to spherical suspensions. At moderate $\phi$ , this effect becomes more pronounced, with $6.2 \pm 0.2 \,\%$ of cubic particles moving upwards, compared to only $2.86 \pm 0.08 \,\%$ for spherical particles at $\phi \approx 0.24$ (figure 7 c).

However, as $\phi$ increases further, the system becomes more densely packed, suppressing velocity fluctuations in both directions. Consequently, the tails of the PDF begin to contract for both particle types (figure 7 a,b), leading to a gradual decrease in the percentage of upward-moving particles, as shown in figure 7(c). In the very dense regime ( $\phi = 0.5$ ), only a small fraction of spherical particles ( $0.070 \pm 0.003 \,\%$ ) exhibit instantaneous upward movement, whereas a slightly larger fraction of cubic particles ( $1.3 \pm 0.1 \,\%$ ) continue to move upwards.

3.6. Suspension microstructure

Figure 8.

The RDF $g(r)$ is obtained to characterize the microstructure of (a) spherical and (b) cubic suspensions, evaluated as functions of radial distance across various volume fractions.

The bulk volume fraction is a key factor in determining the spatial distribution and orientation of settling particles. In this subsection, we analyse the microstructure formation of these particles after reaching steady state. The microstructure of a particle suspension reflects the local arrangement of particles relative to each other, which is governed primarily by inter-particle hydrodynamic interactions. To characterize the overall structures within the suspension, we use the radial distribution function (RDF) (Sundaram & Collins Reference Sundaram and Collins1997; Wang et al. Reference Wang, Wexler and Zhou2000), calculated as an ensemble average over multiple independent configurations. The RDF represents the probability of finding a particle at distance $r$ from a reference particle, often expressed as the ratio of number density within a shell at that distance to the bulk number density:

(3.2) \begin{equation} g(r)=\frac {n(r)}{4\pi \rho r^2\, \delta r}. \end{equation}

Here, $n(r)$ represents the number of particles within a shell of radius $r$ , $\rho$ is the average particle density in the computational domain, and $\delta r$ is the shell thickness. Figure 8 displays the time-averaged RDF as a function of solid volume fraction for both suspension types. In the RDF expression given in (3.2), $g(r)$ becomes zero for distances smaller than the diameter of a particle due to steric interactions that prevent particles from overlapping, and it approaches 1 at larger distances as the shell density aligns with the bulk density. Figure 8 also shows that local particle ordering diminishes at shorter distances for lower volume fractions, while increasing volume fractions enhance this ordering. Peaks in the RDF indicate a higher likelihood of finding particles at specific distances, indicating the gradual formation of a particle ‘cage’ around each reference particle.

For dilute suspensions ( $\phi \simeq 2.5\,\%$ ), the RDF exhibits a subtle peak for both spherical and cubic particles, though it is slightly lower for cubic suspensions (figure 8). In this regime, weak hydrodynamic interactions allow particles to form transient pairs, maintaining contact for extended periods. Consequently, a slightly higher local density of spherical particles is observed around the reference particle compared to cubic particles in this regime. As $\phi$ increases, hydrodynamic interactions among multiple particles become more pronounced, leading to the formation of aggregates or clusters in close contact, as reflected in the sharper RDF peaks. Unlike spherical particles, which maintain a fixed contact point at $r/R = 2$ (figure 8 a), cubic particles interact over a broader range of centre-to-centre distances (Seyed-Ahmadi & Wachs Reference Seyed-Ahmadi and Wachs2019). As a result, the RDF profile for cubic particles appears to be more oblique (figure 8 b). At a moderate solid volume fraction ( $\phi \simeq 36\,\%$ ), the particle-pair distribution undergoes a gradual transition. Initially, particle pairs accumulate around $r/R = 2.2$ for spheres and $r/R = 2.33$ for cubes. This enrichment later shifts to a depletion at greater distances, approximately $r/R \simeq 3.4$ for spheres and $r/R \simeq 3.5$ for cubes. The reduced inter-particle spacing compared to dilute conditions increases the likelihood of cluster formation at short distances. This effect is more pronounced for spherical particles, as evidenced by their higher RDF peaks. In contrast, cubic particles exhibit a greater tendency to escape clustering due to their higher rotational velocities, a behaviour also reported by Seyed-Ahmadi & Wachs (Reference Seyed-Ahmadi and Wachs2019). With further increment in $\phi$ , intensified hindered settling results from reduced inter-particle distances, leading to a more isotropic particle structure. This is reflected in the persistent oscillations beyond the first RDF peak for both suspensions at $\phi \simeq 0.5$ . At larger distances from the reference particle ( $r/R \geqslant 5.1$ for both cases), the particle distribution becomes random, with $g(r) \simeq 1$ .

Figure 9.

The pair distribution function in cylindrical coordinates is presented for suspensions of (a) spheres and (b) cubes across varying solid particle volume fractions.

Although the RDF describes the probability of finding particles within a shell regardless of orientation, it is crucial to assess whether particle pairs show directional preferences in anisotropic suspensions, such as those settling under gravity. To analyse orientation preferences, we compute the pair distribution function, $g(z,r)$ , in cylindrical coordinates for both spherical and cubic suspensions. Here, $g(z,r)$ represents the average distribution of particles located vertically at height $z$ and horizontally at distance $r$ from a reference particle. Figure 9 shows the pair distribution function for various solid volume fractions. At low volume fractions, pairwise interactions dominate, while many particle interactions become important at higher concentrations. As a result of inter-particle hydrodynamic interactions, settling particles naturally rearrange themselves away from their initial random configuration, leading to either clustering or dispersion. At low volume fractions ( $\phi \approx 0.025$ ), there is a notable depletion of particles in the vertical direction near the reference particle for both spherical and cubic suspensions (figure 9). Additionally, particle pairs are enriched in the direction perpendicular to gravity, with peaks at $r/R=2.3$ for spheres and $r/R=1.5$ for cubes. However, cubes show a weaker concentration peak in the horizontal direction (figure 8). A similar situation is observed at $\phi =0.125$ , but the anisotropy is weaker.

As the volume fraction increases, short-range hydrodynamic interactions between nearby particles become more pronounced, causing more particle pairs to cluster around the reference particle, and enhancing structural ordering. We observe a transition from a preferential horizontal alignment at lower bulk concentrations ( $ 0.025 \leqslant \phi \leqslant 0.125$ ) to a more concentric, ring-like distribution at moderate concentrations ( $ 0.15 \leqslant \phi \leqslant 0.4$ ). A spherical contour of particle pairs forms at distances $ r/R = 2.05$ to $r/R= 2.65$ for spheres and $ r/R = 1.15$ to $r/R= 1.6$ for cubes at $ \phi = 0.36$ (figure 9), revealing an isotropic ordering. At this volume fraction, a second ring due to second neighbours becomes visible.

The most intriguing case occurs at the highest volume fraction studied, $\phi =0.5$ , where the pair distribution functions of cubes and spheres exhibit striking differences. For spheres, the pair distribution function shows a high concentration of particles at $r/R \approx 2$ in all directions. Beyond this region, the pair distribution function shows strong fluctuations. As the system is very crowded, spheres take longer to explore the configurational space, and it is difficult to reduce the statistical errors. In contrast, the pair distribution function for cubes shows a stronger ordering in the horizontal direction at high density. Distinct vertical lines appear at distances around $r/R \approx 2$ , followed by separations at $r/R \approx 2.7$ and $r/R \approx 4.5$ , indicating a preferred arrangement of cubes in the $xy$ -plane in a square lattice, with lattice spacing $d\approx 2 R \approx 1.25 h$ , where $h= ( 4\pi /3 ) ^{1/3}R$ is the cube side length. The continuous lines along the $z$ -axis indicate that cubic particles lack ordering along the gravity direction. This suggests that cubes form columnar structures capable of sliding past one another within the densely packed suspension; see also supplementary movie 1.

3.7. Segregation of spherical/cubic particles in a bidisperse suspension

Gravitational sedimentation is a fundamental technique for solid–liquid separation, widely applied in polydisperse systems where size segregation and particle sorting are essential. Although monodisperse systems, comprised of particles of uniform size, have been explored extensively (Brady & Bossis Reference Brady and Bossis1988; Phillips et al. Reference Phillips, Brady and Bossis1988; Ladd Reference Ladd1993; Padding & Louis Reference Padding and Louis2008; Guazzelli & Hinch Reference Guazzelli and Hinch2011; Hamid et al. Reference Hamid, Molina and Yamamoto2013), the dynamics of polydisperse suspensions (systems with particles of varying sizes, densities, or shapes) remains a topic of significant interest. Polydisperse systems present complex behaviour, as the distribution of particle sizes and densities critically influences the overall flow response. Predicting the settling velocity of each particle class in a polydisperse suspension is a primary challenge; understanding these velocities is essential, as they shape the concentration profile and enable effective particle separation by size. Here, we focus on bidisperse sedimentation, which serves as a simplified yet valuable model for examining inter-particle interactions within polydisperse systems. In bidisperse suspensions with low total volume fractions, the hindrance effect becomes more complex than in monodisperse suspensions due to differential impacts of fluid displacement on the settling velocities of distinct particle classes. Although previous studies have measured settling velocities in bidisperse suspensions, these primarily focus on the settling velocities of interfaces or isolated particles within distinct zones of the suspension. When translational symmetry in the vertical direction is disrupted, for example, by the presence of horizontal boundaries, distinct concentration regions can form, leading to varied settling dynamics. In such cases, measurements are typically concentrated on the settling velocities at the interfaces between these zones (Davies Reference Davies1968; Lockett & Bassoon Reference Lockett and Bassoon1979; Al-Naafa & Selim Reference Al-Naafa and Selim1992; Chen et al. Reference Chen, Jia, Fairweather and Hunter2023). In bidisperse suspensions, where particles differ in radii, two primary sedimentation fronts generally propagate from top to bottom: one at the boundary between clear fluid and small particles, and another between small particles and a homogeneous suspension containing both particle sizes. The ratios of particle densities and sizes significantly affect size-based segregation in these systems. Numerical simulations using Stokesian dynamics have been conducted on bidisperse suspensions, covering size ratios up to 4 (Cunha et al. Reference Cunha, Abade, Sousa and Hinch2002; Wang & Brady Reference Wang and Brady2015). Recently, Li & Botto (Reference Li and Botto2024) examined the settling velocity statistics for dilute, non-Brownian, polydisperse sphere suspensions using Stokesian dynamics, analysing how total volume fraction and particle size distribution width influence sedimentation behaviour. They also reported that in a polydisperse mixtures of small and large spheres, the normalized average settling velocity decreases more rapidly with increasing $\phi$ for smaller particles.

Figure 10.

Settling velocities of different species (mixture of spheres and cubes) in a bidisperse suspension are shown for different volume fractions ( $\phi$ ) as functions of time. The velocity of each species is normalized by the Stokes settling velocity of an isolated sphere ( $U_0$ ) in a triply periodic box. We usually employ Stokes time $t_s=R/U_0$ as a time scale, where $R$ is taken as the radius of sphere. Average sedimentation velocities of the spheres, cubes and the whole suspension are indicated by green, red and blue solid lines, respectively.

In this subsection, we examine the settling velocity of non-Brownian bidisperse suspensions composed of a mixture of cubes and spheres, spanning a range of particle volume fractions from low to high. We assume that all particles have a uniform mass density and identical volumes. To achieve this, we calculate an equivalent side length for the cubes, ensuring that each cube has the same volume as a sphere with unit radius. In the bidisperse case, we consider a total of $150$ particles, evenly divided between the two particle shapes. We define the total solid volume fraction ( $\phi$ ) as the sum of the volume fractions of each particle type: $\phi = \phi _1 + \phi _2$ , where $\phi _1$ and $\phi _2$ represent the volume fractions of spheres and cubes, respectively. We examine the case of equal volume fractions for both particle types, i.e. $\phi _1 = \phi _2$ , and apply equal downward force ( $f=1$ ) on all the particles. As illustrated in figure 10, the sedimentation velocity of the bidisperse suspension is strongly influenced by the volume fraction of each species. To reduce statistical error, the velocity distributions are averaged across several random initial configurations of particles within a triply periodic simulation box. In the most dilute case ( $\phi = 0.03$ ), the suspension requires a longer time to reach a steady state. The average sedimentation velocity of each species plays a key role in the overall settling dynamics of the suspension. As shown in figure 10(a), spheres settle faster than cubes. This disparity is due in part to the higher angular rotation induced by the sharp edges of the cubes (Seyed-Ahmadi & Wachs Reference Seyed-Ahmadi and Wachs2019, Reference Seyed-Ahmadi and Wachs2021), which enhances the hindrance effect relative to that experienced by spheres. As the volume fraction increases, interparticle hydrodynamic interactions become more pronounced, causing a moderate reduction in the settling velocity of each species. At higher volume fractions ( $\phi = 0.36$ ), the suspension quickly reaches a steady state, with the reduced settling velocity attributed to increased hindrance effects and diminished velocity fluctuations in both perpendicular and parallel directions. Our study shows that spheres settle faster than cubes at all volume fractions, with segregation efficiency improving at lower $\phi$ .

Figure 11.

The average settling velocity of the particles $U_{{sed}}$ , normalized by the Stokes velocity $U_0$ , is shown as a function of volume fraction ( $\phi$ ). In this context, $U_0$ represents the settling velocity of an isolated sphere in the same periodic domain. Here, the error bars are smaller than the size of the symbols. The inset shows the difference between mean sedimentation velocity of the bisdisperse suspension and that of a monodisperse suspension of spheres.

Figure 11 presents the normalized average settling velocity for different particle types in a bidisperse mixture across various bulk concentrations. The hindered settling behaviour follows a trend similar to that observed in monodisperse suspensions (figure 5), but the mean settling velocities of individual species differ significantly. While the overall velocity of the mixture is slightly lower than that of a monodisperse suspension of spheres in the dilute limit, the variation in settling velocities between species is notable. Figure 11 shows that spheres settle at a rate $11\,\%$ faster than cubes in this range (i.e. $U^{{sphere}}_{{sed}} / U^{{cube}}_{{sed}}=1.11$ ). In a mixture of regular and irregular particles, spheres exhibit the most hydrodynamic shape, minimizing drag and enabling faster settling. In contrast, irregular particles such as cubes, with sharp edges and flat faces, experience higher drag, leading to slower descent. As the volume fraction ( $\phi$ ) increases, the generated upward flux of the surrounding fluid increases, which slows the settling rate. In this regime, the relative velocity difference between spheres and cubes increases to $34 \,\%$ at $\phi \approx 0.25$ . With further increment in $\phi$ , hindered settling becomes more pronounced, reducing somewhat the segregation efficiency. At higher concentrations ( $\phi \approx 0.36$ ), the relative velocity difference narrows to $24 \,\%$ , and the overall mixture velocity approaches a magnitude comparable to that of a monodisperse suspension of spheres.

Figure 12.

Dimensionless hydrodynamic velocity fluctuations of particles as functions of volume fraction ( $\phi$ ) are shown for (a) the settling direction ( $\sigma _{\parallel }$ ) and (b) the direction perpendicular to it ( $\sigma _{\perp }$ ). (c) Comparison of anisotropy in velocity fluctuations of bidisperse mixture with the monodisperse sphere suspension.

In a polydisperse suspension, the particle shape also affects the velocity fluctuations. Figures 12(a,b) present the velocity fluctuations of individual species and the entire mixture in directions parallel and perpendicular to gravity. For dilute bidisperse suspensions ( $\phi \approx 0.1$ ), figure 12(a) shows a clear difference in the velocity fluctuations in the vertical direction between the two species. In this regime, velocity fluctuations in both directions are more pronounced than in monodisperse suspensions. Specifically, the fluctuations of the cube velocities in the vertical direction are $26\,\%$ higher in the mixture than in the monodisperse suspension, while for spheres, the relative difference is only $4\,\%$ . This results in an overall rise in vertical velocity fluctuations in the bidisperse mixture compared to monodisperse cases.

As the solid volume fraction increases, reduced interparticle distances lead to more frequent collisions and stronger hydrodynamic interactions. This effect causes a rise in vertical velocity fluctuations ( $\sigma _{\parallel }$ ) at $\phi \approx 0.23$ for both particle types. In this range, cubic particles in a bidisperse mixture show a $40 \,\%$ increase in $\sigma _{\parallel }$ compared to monodisperse cubic suspensions, while spherical particles exhibit a $27 \,\%$ increase. Simultaneously, horizontal velocity fluctuations steadily increase for both species. These trends collectively alter the degree of anisotropy in velocity fluctuations, as shown in figure 12(c). The combined contributions of individual species lead to a net rise in both $\sigma _{\parallel }$ and $\sigma _{\perp }$ for the entire mixture, directly influencing the anisotropy of velocity fluctuations. As the volume fraction increases from low to moderate values, the anisotropy parameter $\Gamma$ decreases (figure 12 c). While $\Gamma$ follows the same qualitative trend as in figure 6(c), a noticeable difference emerges between the whole mixture and monodisperse sphere suspensions. However, the absolute difference diminishes as the solid volume fraction increases, due to a steady rise in horizontal fluctuations. At higher $\phi$ , intensified hydrodynamic interactions suppress velocity fluctuations in both directions, stabilizing particle motion and leading to a gradual decay in $\sigma _{\parallel }$ and $\sigma _{\perp }$ for both species. Consequently, the overall anisotropy of the bidisperse mixture decreases, though it remains slightly higher than in the monodisperse case. The pair distribution function shows a weaker ordering in the bidisperse suspension compared with the monodisperse suspensions; see supplementary figure 2. The weaker ordering could explain the higher-velocity fluctuations observed in bidisperse suspensions.

Figure 13.

The PDFs of the vertical velocity for (a) spherical and (b) cubic particles in a bidisperse suspension are shown across different solid volume fractions. Negative values of $ U_i / U_{{sed}}$ indicate particles moving upwards, opposing gravity. (c) The percentage of particles ( $ N_{U_z}$ ) moving against gravity at any given moment as a function of solid volume fraction. The inset compares $ N_{U_z}$ in a bidisperse suspension to that in a monodisperse system.

In addition to examining the overall settling behaviour and velocity fluctuations of the particle mixture, we also analyse the normalized PDFs of vertical particle velocities (figure 13) for all particles within the simulation domain. While the PDF for bidisperse suspensions follows a trend similar to that of monodisperse suspensions (figure 7), a key distinction emerges in the distribution of vertical velocities between the two particle species in a bidisperse mixture. Across all solid volume fractions, the PDF of vertical velocity exhibits a slightly higher peak for spherical particles than for cubic particles, indicating that a larger fraction of spheres settle with velocities close to their mean sedimentation velocity. Additionally, a small fraction of particles moves instantaneously against gravity (see supplementary movie 2), even in the dilute regime ( $\phi \lt 0.1$ ; figures 13 a,b). The longer tail in the PDF for cubic particles suggests that cubes are more likely to exhibit upward motion than spheres, contributing to the overall velocity distribution of the mixture. As shown in figure 13(c), $7.67 \pm 0.27 \,\%$ of cubic particles display upward motion, compared to $4.94 \pm 0.2 \,\%$ for spherical particles. This results in a total of $6.38 \pm 0.23 \,\%$ of particles in the bidisperse suspension moving instantaneously against gravity, higher than the $5.27 \pm 0.17 \,\%$ observed in monodisperse suspensions (inset of figure 13 c).

At moderate volume fractions, the tail of the PDF extends further into the negative velocity region due to increased velocity fluctuations and particle dispersion in the plane perpendicular to gravity. Consequently, more particles exhibit upward motion in bidisperse suspensions than in monodisperse ones. As shown in the inset of figure 13(c), at $\phi = 0.23$ , $7.39 \pm 0.2 \,\%$ of particles in the bidisperse suspension move upwards, compared to $6.14 \pm 0.17 \,\%$ in the monodisperse case (see also supplementary movies 3 and 4). However, as $\phi$ increases further, intensified hydrodynamic interactions and inter-particle collisions suppress velocity fluctuations in both directions. As a result, at $\phi = 0.36$ , only $6 \pm 0.26 \,\%$ of cubic particles and $3 \pm 0.14 \,\%$ of spherical particles exhibit upward motion, highlighting the impact of enhanced particle interactions on the settling dynamics of the mixture.