# Hindered settling of log-normally distributed particulate suspensions:
theoretical models vs. Stokesian simulations

Heng Li1  Lorenzo Botto1<EMAIL_ADDRESS>1Process and Energy Department,
ME Faculty of Mechanical Engineering, TU Delft, 2628 CD Delft, The Netherlands

###### Abstract

Settling velocity statistics for dilute, non-Brownian suspensions of
polydisperse spheres having a log-normal size distribution are analysed by
Stokesian Dynamics, as a function of the total volume fraction and width of
the size distribution. Several hundred instantaneous configurations are
averaged to obtain reliable statistics. Average velocities for each particle
class are compared to the models proposed by Batchelor, Richardson $\&$ Zaki,
Davis $\&$ Gecol, and Masliyah-Lockett-Bassoon (MLB). Batchelor’s model is
shown to give reasonably accurate predictions when the volume fraction is
within 5%. Because of its complexity, this model is however hardly used in
practice, so lower-order models are needed. We found that while the other
hindered settling models can give reasonably accurate predictions of the
velocity of the largest particles, all of them overestimate - in certain cases
by a large margin - the velocity of the smaller particles. By computing the
fluid-particle velocity slip for each particle class and using Batchelor’s
model, we explain why predicting the lower tail of the particle size
distribution is challenging, and propose possible avenues for model
improvement. The analysis of velocity fluctuations suggest quantitative
similarities between velocity fluctuations in monodisperse and polydisperse
suspensions.

###### keywords:

Authors should not enter keywords on the manuscript, as these must be chosen
by the author during the online submission process and will then be added
during the typesetting process (see Keyword PDF for the full list). Other
classifications will be added at the same time.

##  1 Introduction

The prediction of the settling velocity of polydisperse suspensions is
crucially important in applications, such as wastewater treatment (He et al.,
2021), food processing (Wang et al., 2020), nanoparticle sorting (Bonaccorso
et al., 2013), particle size characterisation (Papuga et al., 2021), materials
recycling (Wolf, 2021), and sediment transport modelling (Dorrell & Hogg,
2010). Despite decades of research on polydispersed suspension this field
still offers interesting scientific problems. A central challenge is the
prediction of the settling velocity of each particle class in a polydisperse
system. This information is crucially important, not only because the velocity
of each class dictates the particle concentration profile, but also because
only by knowing the velocity of each class it is possible to separate
particles by size. The accurate prediction of the class-averaged particle
velocity has recently become important because of the need for accurate size
fractionation of micro and nanoparticles (Bonaccorso et al., 2013; Backes,
2016). Furthermore, from the knowledge of the class-averaged particle velocity
one can measure the particle size distribution from the time evolution of the
concentration profile (Papuga et al., 2021), which is the principle underlying
the functioning of an analytical centrifuge (Chaturvedi et al., 2018). The
current work aims to analyse the validity of several models used for the
prediction of the class-averaged velocity, comparing against simulation
results.

For a Stokesian suspension of polydisperse spheres grouped into $m$ distinct
particle classes, the average settling velocity of the $i$-th class can be
written as $\langle u_{i}\rangle=u_{St,i}h_{i}(\boldsymbol{\phi})$, where
$u_{St,i}=\frac{2}{9\mu}a_{i}^{2}(\rho_{p}-\rho_{f})g$ is the single-particle
Stokes velocity of the $i$-th class, $h_{i}(\boldsymbol{\phi})$ is the
hindered settling function of that class, and
$\boldsymbol{\phi}=(\phi_{1},\phi_{2},...,\phi_{m})$ is the vector of volume
fractions (Davis & Acrivos, 1985); $a_{i}$ is the particle radius, $\mu$ is
the fluid viscosity, and $\rho_{p}-\rho_{f}$ is the density difference between
the particles and the fluid. The literature reports several models for
$h_{i}(\boldsymbol{\phi})$, as reviewed by Berres et al. (2005).

Batchelor (1982) showed that in the dilute limit the hindered settling
function can be written as

$h_{i}(\boldsymbol{\phi})=1+\sum_{j=1}^{m}S_{ij}\phi_{j},$ (1)

and using the pair-wise interaction approximation developed analytical
solution for the sedimentation coefficients $S_{ij}$ for different size and
density ratios (Batchelor & Wen, 1982). In equation (1) $S_{ii}=-6.55$, so the
hindered settling function for monodispersed suspensions is recovered for
$m=1$ (Batchelor, 1972).

Davis & Gecol (1994) proposed the following semi-empirical extension of
Batchelor’s formula:

$h_{i}(\boldsymbol{\phi})=(1-\phi)^{-S_{ii}}\left(1+\sum_{j=1}^{m}(S_{ij}-S_{ii})\phi_{j}\right),$
(2)

where $\phi=\sum\phi_{j}$ is the total volume fraction, and the coefficients
$S_{ij}$ are defined as in equation (1). for For $\phi\rightarrow 0$, equation
(2) recovers (1). For $m=1$, equation (2) reduces to a power-law form, with
exponent $-S_{ii}$, similar to the expression of Richardson & Zaki (1954).

The models of Batchelor and Davis $\&$ Gecol have been tested in experiments
and simulations of bidisperse suspensions (Davis & Birdsell, 1988; Al-Naafa &
Selim, 1992; Abbas et al., 2006; Wang & Brady, 2015; Chen et al., 2023).
However, these models are rarely used in practice because they contain a large
number of coefficients. Batchelor’s model furthermore does not smoothly
converge to a form that can handle large solid concentrations, and therefore
cannot be used in one-dimensional simulation where the concentration increases
from dilute in the supernatant to concentrated in the region immediately above
the sediment-supernatant interface. Simpler expressions have therefore been
developed for practical predictions of the hindered settling of polydisperse
suspensions.

Some authors (Davis & Hassen, 1988; Abeynaike et al., 2012; Vowinckel et al.,
2019; Chen et al., 2023) have adapted the model of Richardson & Zaki (1954) to
the velocity of each class by writing

$h_{i}(\boldsymbol{\phi})=(1-\phi)^{n}.$ (3)

This model predicts different velocities for different particle radii $a_{i}$,
because $h_{i}$ contains the single-particle Stokes formula at denominator. In
the current paper, this hindered settling formula will be referred to as
Richardson-Zaki’s model for polydispersed suspensions. The exponent $n$ is
obtained from an empirical fit. For monodispersed suspensions, Richardson &
Zaki (1954) originally proposed $n\approx 5$, but the exact value is still a
subject of debate. For example, Brzinski III & Durian (2018) demonstrated that
the value of $n$ depends on the particle Peclet number, and there are two
branches in the sedimentation curve which are best fitted by $n\approx 5.6$
for monodisperse Brownian spheres, and $n\approx 4.48$ for monodisperse non-
Brownian spheres.

Masliyah (1979) and Lockett & Bassoon (1979) proposed the following hindered
settling formula:

$h_{i}(\boldsymbol{\phi})=(1-\phi)^{n-1}\left(1-\sum_{j=1}^{m}\left(\frac{a_{j}}{a_{i}}\right)^{2}\phi_{j}\right),$
(4)

where $a_{j}/a_{i}$ is the ratio between the radii of the $j$-th and the
$i$-th species. The function
$\left(1-\sum_{j=1}^{m}\left(\frac{a_{j}}{a_{i}}\right)^{2}\phi_{j}\right)$ is
the hindered settling function obtained by including the effect of volume
fraction on the fluid-solid slip velocity, and neglecting the effect of
hydrodynamic interactions on the drag force experienced by each particle. The
prefactor $(1-\phi)^{n-1}$ estimates the effect of hydrodynamic interactions
on the drag force (see Appendix A, where equation (4) is re-derived).

The Masliyah-Lockett-Bassoon (MLB) model is favoured in engineering
applications (Xue & Sun, 2003; Berres et al., 2005; Dorrell & Hogg, 2010). It
is easy to tune, since it has only one fitting parameter. Using the MLB model
for the stability analysis of settling size-polydispersed and equal-density
spheres gives no unphysical lateral segregation and streamers, which are
instead obtained using Davis & Gecol’s model (Bürger et al., 2002; Tory et
al., 2003; Berres et al., 2005). The MLB model has been validated by
comparison of predicted and measured concentration profiles (Xue & Sun, 2003;
Berres et al., 2005). However, this validation is not complete, because the
particle concentration is a convolution of the velocities of the different
size classes. Therefore a reasonably accurate prediction of concentration does
not necessarily imply that the velocity of each class has been predicted
correctly. Direct validation of the settling rates predicted by the MLB model
for all size classes has not been published.

Despite large recent interest in the modelling of suspensions of wide and
continuous size distributions (Pednekar et al., 2018; Gonzalez et al., 2021;
Howard et al., 2022; Malbranche et al., 2023; Lavrenteva et al., 2024), there
is very limited data on the settling of polydisperse suspensions with many
size classes. Most physical experiments have been carried out for bidisperse
or tridisperse suspensions (Lockett & Bassoon, 1979; Davis & Birdsell, 1988;
Al-Naafa & Selim, 1992; Davis & Gecol, 1994; Chen et al., 2023). In these
experiments, the largest size ratio between the two species was around 4, and
the velocity of the largest size class was only measured in the homogeneous
region where a mixture of all size classes was present. Numerical simulations
have been carried out for bidisperse suspensions using Stokesian dynamics
(Revay & Higdon, 1992; Cunha et al., 2002; Wang & Brady, 2015) and the force
coupling method (Abbas et al., 2006), with size ratio up to 4. Simulations of
sedimentation of suspensions with a log-normal distribution have been carried
out by Vowinckel et al. (2019) in a domain bounded by top and bottom walls.
Their objective was to study the effect of cohesive force on the transient
settling process. The velocity of each size class was not quantified.

In this paper, we analyse settling velocity statistics for polydisperse
suspensions of non-Brownian spheres using Stokesian dynamics simulations. The
size distribution is log-normal. Of all the particle size distributions, log-
normal distributions are the most interesting because of their ubiquitous
presence in applications (Vowinckel et al., 2019; Di Vaira et al., 2022;
Rettinger et al., 2022). We vary the volume fraction and ratio between the
standard deviation and mean value of the particle size distribution. All the
particles have the same mass density. The volume fraction ranges from 0.01 to
0.1. The largest size ratio between two classes is 5 and up to 9 classes are
considered. To reduce the large statistical error, following other authors
(Revay & Higdon, 1992; Cunha et al., 2002; Abbas et al., 2006; Wang & Brady,
2015), we produce converged particle velocity statistics by generating a large
number of random fixed particle configurations inside a triply periodic box
and ensemble-averaging over all the configurations.

##  2 Numerical approach and validation

Figure 1: A configuration for volume fraction $\phi=0.05$ and polydispersity
parameter $\alpha=0.4$. The spheres are coloured according to their radii.

Consider a polydisperse suspension of $N$ spheres having the same density but
different radii. The $N$ spheres are divided into $m$ size classes. The radius
of size class $i$ is $a_{i}$. Each sphere in class $i$ is subjected to a force
$\boldsymbol{F}_{i}=\frac{4}{3}\pi
a_{i}^{3}(\rho_{p}-\rho_{f})\boldsymbol{g}$, which includes the particle
weight and buoyancy; $\rho_{p}$ and $\rho_{f}$ are the densities of the
spheres and the fluid, respectively, and $\boldsymbol{g}$ is the gravitational
acceleration. The single-particle Stokes velocity corresponding to each class
is
$\boldsymbol{u}_{St,i}=\frac{2}{9\mu}a_{i}^{2}(\rho_{p}-\rho_{f})\boldsymbol{g}$,
where $\mu$ is the dynamic viscosity of the fluid. In the current work, the
particle velocity statistics are calculated from each configuration of the
spheres by first averaging over the particles in the computational domain and
then ensemble-averaging over statistically identical configurations. Each
configuration is generated by randomly placing the spheres one by one inside
the computational domain, ensuring that each placement gives no overlap
between the spheres (Revay & Higdon, 1992; Wang & Brady, 2015; Cheng & Wachs,
2023). One such configuration is shown in figure 1. In our coordinate system,
gravity is aligned in the $z$ direction, also referred to as vertical
direction in the following. The horizontal direction corresponds to the $x$
and $y$ coordinates.

To calculate the velocities of individual particles, a basic version of the
Stokesian Dynamics method is adopted (Brady et al., 1988; Brady & Bossis,
1988). The velocities of the spheres are calculated by solving the mobility
problem

$\boldsymbol{U}-\langle\boldsymbol{u}\rangle=\mathsfbi{M}\boldsymbol{F},$ (5)

where $\boldsymbol{U}$ is the $3N$ vector containing the velocities of the
spheres, $\boldsymbol{F}$ is the $3N$ vector containing the gravitational
forces acting on the spheres (these forces include the particle weight and the
buoyancy force), and $\mathsfbi{M}$ is the $3N\times 3N$ mobility matrix (Kim
& Karrila, 2013). In equation (5), $\langle\boldsymbol{u}\rangle$ is the
average translational velocity of the suspension. In our simulations
$\langle\boldsymbol{u}\rangle=\boldsymbol{0}$ because of the zero volume flux
condition of batch sedimentation (Berres et al., 2005). Note that in the
current work only velocity-force coupling is considered, i.e. the stresslet
and other force moments are not considered. Brady & Durlofsky (1988) showed
that in a sedimenting suspension the inclusion of the stresslet changes the
settling rate negligibly. Because we work in the relatively dilute limit,
short-range lubrication are also neglected.

The mobility matrix $\mathsfbi{M}$ depends on the positions and radii of the
spheres. We used the Rotne-Prager approximation for this term (Rotne & Prager,
1969; Zuk et al., 2014). This approximation has been shown to give accurate
predictions of the sedimentation velocities of suspensions from dilute to
relatively dense (Brady & Durlofsky, 1988). Triply periodic boundary
conditions are applied to the simulation box. The mobility matrix is
constructed using the Ewald summation technique by splitting the mobility
matrix into a real-space part and a wave-space part (Beenakker, 1986).
Explicit formulae for the mobility matrix for a polydisperse suspension can be
found in Beenakker (1986) and Hase & Powell (2001). As characteristic length
and velocity scales, we choose the mean particle radius $\langle a\rangle$ and
the single particle Stokes velocity corresponding to $\langle a\rangle$. To
make forces non-dimensional we use the effective weight of the mean particle,
$\frac{4}{3}\pi\langle a\rangle^{3}(\rho_{p}-\rho_{f})g$.

Figure 2: Normalized settling velocity vs. volume fraction for a simple cubic
array of monodisperse spheres. The line is the point-force solution of
Hasimoto (1959). Upward triangles are the numerical results of Brady et al.
(1988).

In figure 2 numerical predictions for a single sphere in a triply periodic
cubic box are plotted against Hasimoto’s analytical solution (Hasimoto, 1959)
and the simulation results of Brady et al. (1988). The volume fraction of the
simple cubic array is varied by varying the size of the box. Based on the
point-force assumption, Hasimoto (1959) derived $u/u_{St}=1-1.7601\phi^{1/3}$
for $\phi\ll 1$, where $u_{St}$ is the Stokes velocity of the sphere. Brady et
al. (1988) used Stokesian Dynamics with different approximations for the
mobility matrix. The results of Brady et al. (1988) shown in figure 2
correspond to simulations in the Rotne-Prager approximation. As seen from
figure 2, our results match exactly those of Brady et al. (1988) and converge
to Hasimoto’s solution for $\phi\rightarrow 0$. This test validates our
implementation of the Ewald summation for the periodic boundary conditions.

Figure 3: Normalized relative settling velocity for a pair of spheres as a
function of the centre-to-centre distance for (a) size ratio 2
and (b) size ratio 5. Results of current simulations are shown as symbols, and
analytical results of Wacholder & Sather (1974) are shown as lines.

In figure 3, the normalized relative settling velocity is shown as a function
of the normalized centre-to-centre distance between two unequal spheres with
size ratio 2 and 5, respectively. In our simulations, the radius of the large
sphere is fixed to $a_{l}=2$. The radius of the small sphere is $a_{s}=1$ and
0.4 for size ratio 2 and 5, respectively (these values are chosen because the
largest radius is 2 and the largest size ratio is 5 in the polydisperse
simulations analysed in this paper). The relative settling velocity between
the two spheres is normalized by the Stokes velocity of the large sphere. The
centre-to-centre distance is normalized by the radius of the large sphere. In
figure 3, symbols are results from our simulations, and lines correpond to the
asymptotic solution of Wacholder & Sather (1974), in which only far-field
hydrodynamic interactions were considered. It can be seen that our results
match the analytical solution for both vertically and horizontally aligned
pairs.

The current paper discusses results for bidisperse suspensions and
polydisperse suspensions with more than two classes, also comparing with the
monodisperse case. For the monodisperse case, the radius of the spheres is
$a=1$. For the bidisperse case, two size ratios are considered:
$a_{2}/a_{1}=2$ and 5. The radii of the small size classes are $a_{1}=0.8$ and
0.4 for these two size ratios, respectively. The volume fraction of the small
size class is $\phi_{1}=\frac{3}{11}\phi$ for size ratio 2, and
$\phi_{1}=\frac{1}{76}\phi$ for size ratio 5. These volume fraction ratios are
chosen so that the average radius of the spheres is equal to 1.0 for each
system.

Figure 4: Discrete frequency distributions of particle size for different
values of the polydispersity parameter (symbols). Lines indicate the
continuous log-normal distributions that fit the discrete frequency
histograms.

For the simulations with several size classes, the particle size distribution
follows $p(a)=\frac{1}{a\sigma\sqrt{2\pi}}e^{-(\ln{a}-\mu)^{2}/2\sigma^{2}}$,
where the mean value of the size distribution is $\langle
a\rangle=e^{\mu+\sigma^{2}/2}$ and the standard deviation is $\Delta
a=\sqrt{\left(e^{\sigma^{2}}-1\right)e^{2\mu+\sigma^{2}}}$. We define the
polydispersity parameter as $\alpha=\Delta a/\langle a\rangle$. Four size
distributions are considered, with $\langle a\rangle=1$ and
$\alpha=$0.1,0.2,0.3 or 0.4. Each distribution is cut at the two ends,
resulting in a range $a_{min}\leq a\leq a_{max}$, where $a_{min}$ and
$a_{max}$ are chosen so that at least 95% of the original distribution falls
within this range. The largest size ratio between two spheres is 5. Each
radius range is discretized into between 4 and 9 size classes, with a
difference of 0.2 between the radii of two adjacent size classes.

The discrete number frequency distributions are overlaid on the corresponding
continuous distributions in figure 4 . The frequency of size class $i$ is
calculated as $\frac{p(a_{i})}{\sum_{j=1}^{m}p(a_{j})}$. For each value of
$\alpha$, the volume fraction $\phi$ ranges from 0.01 to 0.10. For each
simulated case, corresponding to a combination of $\alpha$ and $\phi$, a fixed
box size $L=80$ is used and 500 random particle configurations are generated.
The total number of spheres in each case varies from 925 to 12223.

The average velocity of class $i$ is calculated by ensemble-averaging over $M$
configurations as

$\langle\overline{{u}}_{\xi,i}\rangle=\frac{\sum_{k=1}^{M}\overline{u}_{\xi,i}^{k}}{M},$
(6)

where $\xi=1,2,3$ correspond to the three Cartesian coordinates,
$\overline{u}_{\xi,i}$ is the intrinsic volume average of the velocity
component ${u}_{\xi,i}$ within one configuration, and $\left<\cdot\right>$ is
the ensemble-averaging operator. The intrinsic volume average within class $i$
over configuration $k$ is
$\overline{u}_{\xi,i}^{k}=\frac{\sum_{l=1}^{N_{i}}u_{\xi,i,l}^{k}}{N_{i}}$,
where $N_{i}$ is the number of particles in class $i$. The standard deviation
of a certain velocity component within one realisation is calculated as
$\left(u_{\xi,i}^{\prime}\right)^{k}=\sqrt{\frac{\sum_{l=1}^{N_{i}}\left(u_{\xi,i,l}^{k}-\overline{u}_{\xi,i}^{k}\right)^{2}}{N_{i}-1}}$.
Averaging over many realisations gives an improve estimate of the class-
averaged standard deviation. In the bulk of the paper we indicate averages by
the bracket symbol, distinguishing between volume and ensemble average when
necessary.

###  2.1 Relation between the mobility formulation and Batchelor’s formula

In this section we show the connection between the mobility formulation,
equation (5), and Batchelor’s formula, equation (1). For simplicity of
notation, let us consider a specific size class. Without loss of generality we
consider class 1. According to (5) the velocity of the $\alpha$-th sphere in
the 1-st size class is

$\boldsymbol{u}_{\alpha,1}=\sum_{i=1}^{m}\sum_{\beta=1}^{N_{i}}\mathsfbi{M}_{\alpha
1,\beta i}\boldsymbol{F}_{i},$ (7)

where $N_{i}$ is the number of spheres in the $i$-th class, and
$\mathsfbi{M}_{\alpha 1,\beta i}$ is the $3\times 3$ mobility matrix
representing the hydrodynamic interaction between the $\alpha$-th sphere in
the 1-st class and the $\beta$-th sphere in the $i$-th class (Kim & Karrila,
2013). Because $\mathsfbi{M}_{\alpha 1,\alpha 1}=(6\pi\mu a_{1})^{-1}$, (7)
can be written as

$\boldsymbol{u}_{\alpha,1}=\boldsymbol{u}_{St,1}+\sum_{\beta\neq\alpha}\mathsfbi{M}_{\alpha
1,\beta 1}\boldsymbol{F}_{1}+\sum_{i\neq
1}\sum_{\beta=1}^{N_{i}}\mathsfbi{M}_{\alpha 1,\beta i}\boldsymbol{F}_{i}.$
(8)

The average velocity of the 1-st class in this configuration is

$\overline{\boldsymbol{u}}_{1}=\boldsymbol{u}_{St,1}+\frac{1}{N_{1}}\left(\sum_{\alpha}\sum_{\beta\neq\alpha}\mathsfbi{M}_{\alpha
1,\beta 1}\boldsymbol{F}_{1}+\sum_{\alpha}\sum_{i\neq
1}\sum_{\beta=1}^{N_{i}}\mathsfbi{M}_{\alpha 1,\beta
i}\boldsymbol{F}_{i}\right),$ (9)

but because $\boldsymbol{F}_{i}$ is constant within the same size class we can
also write

$\overline{\boldsymbol{u}}_{1}=\boldsymbol{u}_{St,1}+\mathsfbi{s}_{11}\boldsymbol{F}_{1}+\sum_{i\neq
1}\mathsfbi{s}_{1i}\boldsymbol{F}_{i},$ (10)

where $\mathsfbi{s}_{11}$ and $\mathsfbi{s}_{1i}$ describe the intra-class
hydrodynamic interactions (within the 1-st class) and the inter-class
hydrodynamic interactions (between the 1-st and the $i$-th classes),
respectively. These two matrices can be written as
$\mathsfbi{s}_{11}=(N_{1}-1)\overline{\mathsfbi{M}}_{11}$ and
$\mathsfbi{s}_{1i}=N_{i}\overline{\mathsfbi{M}}_{1i}$, where
$\overline{\mathsfbi{M}}_{11}$ and $\overline{\mathsfbi{M}}_{1i}$ are the
average two-sphere mobility matrices. Upon ensemble-averaging, the average
velocity of the 1-st size class can be written as

$\langle\boldsymbol{\overline{u}}_{1}\rangle=\boldsymbol{u}_{St,1}+\langle\mathsfbi{s}_{11}\rangle\boldsymbol{F}_{1}+\sum_{i\neq
1}\langle\mathsfbi{s}_{1i}\rangle\boldsymbol{F}_{i}.$ (11)

The average velocity component in the gravity direction can be written as

$\frac{\langle u_{1}\rangle}{u_{St,1}}=1+\frac{9\mu\langle
s_{11}\rangle}{2a_{1}^{2}n_{1}}\phi_{1}+\sum_{i\neq 1}\frac{9\mu\langle
s_{1i}\rangle}{2a_{1}^{2}n_{i}}\phi_{i},$ (12)

where the formula for the single-particle Stokes velocity is used and $n_{i}$
is the number density of class $i$. The scalar $\langle s_{1i}\rangle$ is the
component of $\langle\mathsfbi{s}_{1i}\rangle$ for the velocity-force coupling
in the gravity direction.

Extending equation (12) to a generic class $i$ yields

$\frac{\langle
u_{i}\rangle}{u_{St,i}}=1+B_{ii}(\boldsymbol{\phi})\phi_{i}+\sum_{j\neq
i}B_{ij}\left(\boldsymbol{\phi},\frac{a_{j}}{a_{i}}\right)\phi_{j}.$ (13)

The dependence of $B_{ii}$ and $B_{ij}$ on the volume fraction vector
$\boldsymbol{\phi}$ comes from the fact that $\langle s_{ij}\rangle$ depends
on the pair distribution functions, and the pair distribution functions in
turn depend on the volume fraction of each class. The dependence of $B_{ij}$
on $a_{j}/a_{i}$ comes from the dependence of the two-sphere mobility matrix
on the size ratio. For $\phi\rightarrow 0$, $B_{ii}$ is a constant and
$B_{ij}=S_{ij}$ is only a function of $a_{j}/a_{i}$. In this limit, equation
(13) recovers Batchelor’s expression (1).

##  3 Hindered settling of monodisperse and bidisperse suspensions

Figure 5: Monodispese case: average settling velocity, normalized by the
single-particle Stokes velocity, versus volume fraction

The normalised average settling velocity $\langle u_{z}\rangle/u_{St}$ for the
monodisperse suspension is plotted in figure 5 as a function of $\phi$. We
include in the plot the Richardson-Zaki correlation $(1-\phi)^{n}$ (Richardson
& Zaki, 1954) for $n=5$, the Batchelor model $1+S\phi$ (Batchelor, 1972)
assuming $S=-6.55$ and the Hayakawa-Ichiki model
$\frac{(1-\phi)^{3}}{1+2\phi+1.429\phi(1-\phi)^{3}}$ (Hayakawa & Ichiki,
1995). The values chosen for the exponent $n$ and the coefficient $S$ here are
typically for non-Brownian particles interacting only hydrodynamically
(Padding & Louis, 2004; Moncho-Jordá et al., 2010).

Our simulation results agree with Batchelor’s model for $\phi$ approximately
smaller than 0.03. For larger volume fractions, the simulation gives larger
values than Batchelor’s model. A similar range of validity for Batchelor’s
model was also found by Abbas et al. (2006) using a force-coupling method. Our
results also agree well with the Hayakawa-Ichiki model for $\phi\leq 0.05$ and
they lie between the predictions of Richardson-Zaki’s correlation and
Hayakawa-Ichiki’s model for $\phi\geq 0.06$. The simulation data for
$\phi=0.01$ is smaller than the values predicted by the three models. This is
expected because of the use of triply-periodic boundary conditions in a domain
of finite size (Phillips et al., 1988).

Figure 6: Bidisperse case: average settling velocity normalized by the single-
particle Stokes velocities for the small (left panels) and large (right
panels) particles. Panels (a) and (b) are for size ratio 2. Panels (c) and (d)
are for size ratio 5.

Normalised average settling velocities for the small and the large particles
in the bidisperse case are plotted as symbols in figure 6 for two size ratios.
The predictions of Batchelor’s model (equation (1)), Davis & Gecol’s model
(equation (2)) and the MLB model (equation (4)) are indicated by lines. It is
seen from figure 6 (a) and (c) that our results for the small particles agree
with the predictions of Batchelor’s model for $\phi\leq 0.05$, and lie between
the predictions from Batchelor’s model and Davis & Gecol’s models for
$\phi\geq 0.06$. For the large particles, our results agree with predictions
from Batchelor’s model for $\phi\leq 0.03$ and lie between the predictions
from Davis & Gecol’s and MLB models for $\phi\geq 0.04$. Stokesian dynamics
calculations by Wang & Brady (2015) that include stresslet and lubrication
contributions also predicted for $\phi$ larger than around 0.05 hindered
settling velocities smaller and larger than those of Davis & Gecol’s for the
small and the large particles, respectively.

In summary our simulation results for the mono- and bi-disperse cases are
comparable to those in the literature.

##  4 Polydisperse suspensions

Figure 7: Configuration for different polydispersity parameters and
$\phi$=0.05, with the spheres coloured according to their settling velocity;
(a) $\alpha$=0, (b) $\alpha$=0.2, and (c) $\alpha$=0.4.

To illustrate the spatial distribution of particle velocities in the
polydisperse particle simulations, in figure 7 we show snapshots of the
simulations with each sphere coloured according to its settling velocity.
Spheres coloured in red have settling velocities in the direction of gravity
whereas spheres coloured in blue have settling velocities opposite to gravity.
Figure 7 (c) shows that the smaller particles can move against gravity, and
have negative velocities comparable in magnitude to the positive settling
velocity of the largest particles.

Figure 8: Probability distribution functions of (a) horizontal
and (b) vertical velocities for different polydispersity parameters and
$\phi$=0.05.

Figure 9: Probability distribution functions of (a) horizontal
and (b) for $\alpha$=0.4 at $\phi$=0.05. In contrast to Fig. 8, here the PDFs
are calculated based on the distribution of velocities within each size class.

Probability distribution functions (PDFs) of horizontal and vertical
velocities, shown in figure 8 for different values of $\alpha$, are
approximately Gaussian, with a variance that increases as $\alpha$ increases.
These PDFs are constructed by considering all the particles in the simulation
domain. However, spheres belonging to the same size class also have a
distribution of settling velocities. Therefore, in figure 9, we show the PDFs
of the horizontal and the vertical velocities of spheres in each size class
for $\alpha$=0.4. For comparison, the PDFs of the velocities of all the
spheres are included in this plot as grey lines. Again, the PDFs are
approximately Gaussian (simulations by Cheng & Wachs (2023) of uniform flow
past fixed polydisperse random arrays indicate also a Gaussian distribution
for the hydrodynamic forces of a given size class). Surprisingly, the PDFs of
the horizontal velocities for different size classes collapse onto a single
curve (figure 9 (a)). From the PDFs of the vertical velocities in 9 (b), it is
seen that the mean velocity increases as the size of the particle increases.
And, different size classes have comparable variances.

Figure 10: Normalized average settling velocity of each size class for
different polydispersity parameters and $\phi$=0.05. The inset is a zoom in
the range $0.8\leq a_{i}\leq 2$. The lines are guides for the eyes.

The average vertical settling velocity of each size class normalized by the
corresponding single-particle Stokes velocity is shown in figure 10 for
different values of $\alpha$. The inset shows a zoom in the range $0.8\leq
a_{i}\leq 2$. Because now the settling velocity is normalised by the single-
particle settling velocity, the information in this plot complements the data
of figure 9 (b). We see that for fixed $\alpha$ the normalized average
settling velocity increases as the particle size increases. This means that
the velocities of small particles are more hindered than the velocities of
large particles. For a given size class, the normalized average settling
velocity decreases as $\alpha$ increases, and decreases faster for small
particles than for large particles.

Figure 11: Normalized average settling velocities of each size class, for
$\alpha$=0.4 and different volume fractions. The inset shows a zoom in the
range $1\leq a_{i}\leq 2$.

In the previous figures, the total volume fraction was fixed, and $\alpha$ was
changed. In figure 11, we instead change $\phi$ for fixed $\alpha$=0.4. This
plot confirms the trend seen in figure 10: for a given volume fraction, the
normalized average settling velocity decreases as the particle size decreases.
The normalized average settling velocity decreases faster with increasing
$\phi$ for small size particles.

To summarise, the smaller particles are more hindered and more affected by
polydispersity than the large ones.

###  4.1 Comparison with hindered settling models

In this subsection, current simulations are compared with predictions from
Batchelor’s (see equation (1)), Davis & Gecol’s (see equation (2)) and MLB
(see equation (4)) models. The accuracy of Richardson-Zaki correlation (see
equation (3)) for polydisperse suspensions is also checked. The values of the
coefficients $S_{ij}$ in Batchelor’s and Davis & Gecol’s models are calculated
from $S_{ij}=-3.50-1.10\lambda-1.02\lambda^{2}-0.002\lambda^{3}$ where
$\lambda=a_{j}/a_{i}$ (Davis & Gecol, 1994). The value of the exponent $n$ in
the MLB model and the Richardson-Zaki correlation is 5.

Figure 12: Comparison between the current simulation results and the
predictions of hindered settling function models for the average velocity of
different size classes, for $\phi$=0.05 and different polydispersity
parameters: (a) $\alpha$=0.1, (b) $\alpha$=0.2, (c) $\alpha$=0.3 and (d)
$\alpha$=0.4.

Hindered settling functions corresponding to different size classes for fixed
$\phi$=0.05 and different $\alpha$ are compared against different theoretical
models in figure 12. The Richardson-Zaki correlation largely overestimates the
hindered settling functions of smaller particles, whereas it gives reasonable
values for larger particles. The discrepancy between the Richardson-Zaki
correlation and the computed hindered settling functions of smaller particles
increases as $\alpha$ increases. For each $\alpha$, the predictions from the
other three models show a consistent trend for each size class. The MLB model
gives the largest settling velocities, Batchelor’s model gives the smallest
settling velocities, and Davis & Gecol’s model gives intermediate values. The
differences between the predictions from these three models get smaller as the
particle size increases.

Figure 13: Relative differences between the average settling velocities from
different models and from the current simulations for fixed $\phi=0.05$: (a)
Batchelor’s model; (b) Davis & Gecol model; (c) MLB model.

Figure 13 shows the normalized relative differences between the computed and
predicted average settling velocities. The Batchelor model and the Davis &
Gecol model predict the average settling velocity of each size class quite
well for all $\alpha$ considered here, with relative errors smaller than 10%,
except for the smallest size class $a_{i}$=0.4 for $\alpha$=0.4 for which the
simulation gives a very small settling velocity. From figure 13 (c), it is
seen that the relative difference between the predictions from the MLB model
and the current simulations decreases as $a_{i}$ increases, or $\alpha$
decreases. For $a_{i}\geq 1$, the MLB model predicts the average settling
velocities quite well, with the relative difference within 10%. For $a_{i}\leq
0.8$, the MLB model starts failing.

Figure 14: Comparison between the current numerical results and the
predictions of different models. The comparison is here evaluated as a
function of $\phi$ for fixed $\alpha$=0.4 . (a) to (f) correspond to size
classes $a_{1}$,$a_{2}$,$a_{4}$,$a_{6}$,$a_{7}$ and $a_{9}$, respectively.

Hindered settling function data for fixed $\alpha$=0.4 and different $\phi$
are compared with the predictions from different models in figure 14. Also,
Richardson-Zaki’s correlation predicts the hindered settling functions of
smaller particles poorly. For the other three models, a similar trend in the
predicted values is observed as the one in the case of varying $\alpha$. The
predictions from the MLB model and the Davis & Gecol model are quite close to
each other for all size classes at each volume fraction, and they are also
quite close to the values from current simulations for larger particles.
However, the Davis & Gecol model slightly underestimates the hindered settling
functions of the larger particles when $\phi>0.06$. For smaller particles, the
predictions from the MLB model and the Davis & Gecol model are larger than the
values from the simulations, and the discrepancies between the predictions
from these two models and the values from current simulations get larger as
$\phi$ increases. The predictions of Batchelor’s model are close to the
simulated values for $\phi$ approximately less than 0.05. As $\phi$ increases,
Batchelor’s model underestimates the hindered settling functions of all size
classes systematically compared to the results of current simulations.

Figure 15: Same as in Fig. 13, but for different volume fractions and fixed
$\alpha=0.4.$

For fixed $\alpha$=0.4, the relative differences between the average settling
velocities predicted by different models and calculated by current simulations
of each size class for different volume fractions are shown in figure 15. For
each size class, the relative difference between the prediction from the
Batchelor model and the current simulations increases as the volume fraction
increases, and it is within 10% when $\phi\leq 0.05$, except for the smallest
size class $a_{i}=0.4$. From figure 15 (b) and (c), it is seen that the
relative differences are quite close for the Davis & Gecol and the MLB models,
with those of the MLB model slightly larger. For larger size classes
($a_{i}\geq 1$), both these two models give quite accurate predictions for all
volume fractions considered, with the relative differences within 10% compared
to the results of current simulations. For smaller size classes ($a_{i}\leq
0.8$), both these two models give predictions with large relative differences
compared to the results of current simulations, and in general the relative
difference gets larger as volume fraction increases or as size of the class
decreases.

We saw that the MLB model, despite its simplicity, gives relatively good
agreement for the large particles. However, it fails for the small particles.
The MLB model is based on a closure relation for the particle-fluid velocity
difference (see Appendix A). Therefore, to understand the limits of validity
of the model, we compute the slip velocity from the simulation data and
compare against the MLB model prediction.

Figure 16: Normalized slip velocity for each size class for $\phi$=0.05 and
different $\alpha$. The line shows the prediction of the MLB model.

Slip velocities for each size class normalized by the corresponding Stokes
velocities are plotted in figure 16 for $\phi$=0.05 and different $\alpha$.
The slip velocity for size class $i$ is defined as the difference between the
average settling velocity of class $i$ and the average velocity of the fluid
phase,

$u_{slip,i}=\langle u_{z,i}\rangle-\langle u_{f}\rangle.$ (14)

The value of $\langle u_{f}\rangle$ is obtained from the the zero volume-flux
condition $\sum\phi_{j}\langle u_{z,j}\rangle+(1-\phi)\langle u_{f}\rangle=0$.
The slip velocity predicted by the MLB model is calculated from (see Appendix
A)

$u_{slip,i}=u_{St,i}(1-\phi)^{n-1}.$ (15)

It is seen that the MLB model does not predict accurately the slip velocities
of the smaller particles. The discrepancy between the MLB model prediction and
the simulation data increases as $\alpha$ increases. The slip velocities of
relatively large particles are reasonably well captured. As the particle size
increases, the simulation data tends to converge to the MLB model prediction.

Figure 17: Normalized slip velocity for each size class for $\alpha$=0.4 and
different values of $\phi$. The dashed line is the prediction of the MLB
model.

The normalized slip velocities for each size class for fixed $\alpha$=0.4 and
varying $\phi$ are plotted in figure 17. It is seen that the prediction of the
MLB model gets increasingly worse as the volume fraction increases for the
small size classes. Predictions for the largest particles are instead
acceptable regardless of the volume fraction.

##  5 Particle velocity fluctuations

Figure 18: Normalized (a) horizontal
and (b) vertical velocity fluctuations for different $\alpha$ and $\phi$=0.05.

Figure 19: Normalized (a) horizontal
and (b) vertical velocity fluctuations for different $\phi$ and $\alpha$=0.4.

Figure 20: Normalized velocity fluctuations of each size class at $\phi=0.05$
for (a) $\alpha=0.2$
and (b) $\alpha=0.4$.

Statistical deviations with respect to the mean particle velocity, as measured
by the root-mean square of the velocity fluctuation, increase as $\alpha$ or
$\phi$ increase (see figures 18 and 19). The normalized horizontal and
vertical velocity fluctuations of each size class are shown for $\alpha=0.2$
and 0.4 with fixed $\phi=0.05$ in figure 20. For a fixed $\alpha$, the
normalized velocity fluctuations decrease as the particle size increases.
Figure 9 seems to suggest that the velocity fluctuations are approximately
independent of the particle radius $a_{i}$. Because the Stokes velocity scales
as $a_{i}^{2}$, it is expected that the velocity fluctuations normalized by
the Stokes velocity scale as $\sim a_{i}^{-2}$. Our data confirm this scaling
(see lines in figure 20): $u_{i}^{\prime}/u_{St,i}=ca_{i}^{-2}$ fits the data
for all the values of $\alpha$ and $\phi$ simulated, as shown in tables 1 and
2 (this scaling is also observed in our simulations of bidisperse suspensions,
as velocity fluctuations for the two classes are similar in magnitude).

$\alpha$ | horizontal direction | vertical direction
---|---|---
0.1 | 0.32 | 1.13
0.2 | 0.38 | 1.33
0.3 | 0.46 | 1.56
0.4 | 0.53 | 1.85
Table 1: Approximate values of the prefactor $c$ in the scalings of the horizontal and the vertical normalized velocity fluctuations for each $\alpha$ at $\phi=0.05$. $\phi$ | horizontal direction | vertical direction
---|---|---
0.01 | 0.28 | 0.95
0.03 | 0.45 | 1.54
0.05 | 0.53 | 1.85
0.08 | 0.61 | 2.08
0.10 | 0.62 | 2.15
Table 2: Same as table 1 but for different $\phi$ at $\alpha=0.4$.

Figure 21: Ratio between vertical and horizontal velocity fluctuation
magnitudes for $a_{i}=1$ and (a) different $\alpha$ at $\phi=0.05$,
or (b) different $\phi$ at $\alpha=0.4$.

The anisotropy ratio between the vertical and the horizontal velocity
fluctuations, plotted in figure 21, is around 3.5 regardless of the values of
$\alpha$ or $\phi$. This value was also observed in the monodisperse and
bidisperse simulations.

Peysson & Guazzelli (1999) measured experimentally the velocity fluctuations
of small and large particles in a dilute bidisperse suspensions with size
ratio 2. They found that the ratio of velocity fluctuations between the small
and the large size classes were around 0.85 and 0.75 in the vertical and
horizontal directions, respectively, close to $(a_{1}/a_{2})^{1/2}\approx
0.71$. To explain their findings, they extended Hinch’s scaling for the
velocity fluctuations in monodisperse suspensions (Hinch, 1988), using two
different correlation lengths for the two size classes. In our simulations,
the correlation length should be the size of the computational box, as found
in other numerical simulations using periodic boxes (Ladd, 1996; Nguyen &
Ladd, 2005) and in other experiments during the initial settling stage of
well-mixed suspensions (Segre et al., 1997; Tee et al., 2002). Using the same
correlation length for the two size classes would predict comparable velocity
fluctuation magnitudes.

##  6 Discussion

Figure 22: Inter-class interaction term appearing on the right hand side of
equation (16) for $a_{i}=0.4$ (“small” particles) and $a_{i}=2$ (“large”
particles). The particle size distribution corresponds to $\alpha=0.4$ and
$\phi=0.05$. Figure 23: (a) Interaction coefficients $S_{ij}-S_{ii}$
and (b) volume fraction distribution corresponding to the interaction term of
figure 22.

The good comparison between Batchelor’s model and the simulation data for
sufficiently small $\phi$ enables us to use this analytical model to
illustrate why the prediction of the velocity of the small particles is highly
dependent on the full particle size distribution, while that of the large
particles is not. Equation (1) can be rewritten as

$h_{i}=1+S_{ii}\phi+\sum_{j=1}^{m}(S_{ij}-S_{ii})\phi_{j},$ (16)

where $S_{ii}=-6.55$ and $\phi$ is the total volume fraction. The influence of
particle class $j$ on particle class $i$ is negligible if
$|(S_{ij}-S_{ii})\phi_{j}|\ll|S_{ii}\phi|$. For $\phi=0.05$, the magnitude of
the intra-class interaction term is $|S_{ii}\phi|=0.33$. Let us compare this
value to the inter-class interaction term for $\phi=0.05$. In figure 22
$(S_{ij}-S_{ii})\phi_{j}$ is shown for $a_{i}=0.4$ (small particles) and
$a_{i}=2$ (large particles), for $\alpha=0.4$. The maximum absolute value of
the inter-class interaction term for the small particle is $0.13$, not
negligible in comparison to $0.33$. The maximum value of the inter-class
interaction term for the largest particles is instead 18 times smaller than
the intra-class interaction term. The question is: why is the inter-class
interaction term small for the large particles? Is this because the
interaction coefficients are small in magnitude? Or because of the
distribution of volume fractions?

To reply to these questions, in figures 23 (a) and (b) we show
$(S_{ij}-S_{ii})$ and $\phi$ separately. It is seen that in our log-normal
distribution the volume fraction corresponding to the small particles is small
in comparison to that of the large particles, and tends to zero as the lower
tail of the particle size distribution is approached. The quantity
$(S_{ij}-S_{ii})$ is on the other hand not diverging for $a_{j}/a_{i}\ll 1$,
and is $O(1)$ in this limit. Therefore, the reason why the lower tail of the
distribution has a small influence on the upper tail is that the volume
fraction corresponding to the lower tail is comparatively small and is
weighted by an interaction term that is not large. This is an insight that
could be applied to other particle size distributions. For example, if the
particle size distribution was such that $\phi$ was comparatively large in the
small particle range, we would expect the settling velocity of the largest
particles to be more affected by the smallest particles than seen in our
simulations.

The analysis above also gives an insights into the condition for which models
parameterised on the total volume fraction can be used as a first, practical
approximation for the prediction of the settling of a dilute polydisperse
suspension. This approximation is reasonable when the inter-class interaction
term is comparatively small. This term is small when either $\phi_{j}\ll 1$
for finite $S_{ij}-S_{ii}$, the case discussed above. Or when the particle
size distribution is narrow so that $|S_{ij}-S_{ii}|\rightarrow 0$, the case
discussed by Davis & Hassen (1988) (see the value of $S_{ij}-S_{ii}$ for
$a_{j}/a_{i}$ approaching 1 in figure 22 (a)). If deviations of $S_{ij}$ from
$S_{ii}$ are small, then it can be seen from Eq. (13) that the hindered
settling function for $\phi\ll 1$ depends only on the total volume fraction.
For a dilute suspension with a narrow size distribution, the use of
Richardson-Zaki’s correlation is for example justified, and is not by chance
that the exponent of the Richardson-Zaki correlation ($n\simeq 5$) is
numerically close to $|S_{ii}|=6.5$. This “lucky coincidence” was also noted
by Davis & Hassen (1988).

The comparison with the simulation data shows that the MLB model tends to
overestimate the hindered settling functions of particles with sizes smaller
than the mean. The MLB model gives more accurate predictions than Richardson-
Zaki’s correlation for comparatively large particles. MLB is also based on
Richardson-Zaki’s formula, but in the MLB model the formula is used to
estimate the particle-fluid slip velocity, not the absolute settling velocity.
For applications where the focus is predicting the sedimentation of the larger
particles (e.g., separation of large particles from a polydisperse mixture),
using the MLB model could be sufficient. For applications where the
stratification in different layers needs to be predicted (e.g. in
sedimentology (Dorrell & Hogg, 2010)), using the MLB model will overestimate
the fraction of the smaller particles in the sediment region. In particle size
fractionation by centrifugation or sedimentation (Bonaccorso et al., 2013;
Backes, 2016)), using the MLB model could give a wrong prediction of the
region where most of the small particles are located, jeopardising the size
fractionation protocol.

Looking at the main assumptions of the MLB model, rederived in the Appendix A,
we can see that the model rests on two key assumptions. The first one is that
the buoyancy force on each sphere in a polydisperse suspension depends on the
suspension density. The conceptual difficulty of modelling the buoyancy force
on a small sphere in a suspension of spheres of a different size has been
addressed in several papers (Gibb, 1991; Di Felice, 1995; Ruzicka, 2006;
Piazza et al., 2013). Experiments show that using the density of the
suspension to evaluate the buoyancy force gives accurate predictions of the
settling velocity of the test sphere when the test sphere is comparable in
size or larger than the other particles, but can lead to inaccuracies when the
test particle is smaller than the other particles (Poletto & Joseph, 1995;
Rotondi et al., 2015). A second assumption in the MLB model is that the Stokes
drag correction for each size class only depends on the total volume fraction,
and we have seen that this cannot be in general a good approximation.

We have not tried to improve the MLB model. However, we can give indications
on possible approaches for improvement. When the MLB model was derived, no
accurate drag correlations for polydisperse sphere arrays were available.
Recently, numerical simulation of flow past fixed sphere arrays (van der Hoef
et al., 2005; Yin & Sundaresan, 2009; Sarkar et al., 2009; Cheng & Wachs,
2023) have shown that data for the drag on polydispersed sphere arrays can be
fitted with low-order polynomials of the the first two moments of the size
distribution. Improved models for poly-disperse suspensions could therefore be
built starting from the MLB model but using the drag force correlation for
polydispersed fixed arrays to model the slip velocity. Such models could also
account for recent studies on the shear viscosity of polydisperse suspensions
(Wagner & Woutersen, 1994; Dörr et al., 2013; Mwasame et al., 2016; Pednekar
et al., 2018), because in certain limits there is a relation between the drag
on a sphere translating through a suspension and the shear viscosity of the
suspension (Squires & Mason, 2010).

For a continuum size distribution, equation (13) which is also valid at non-
negligible volume fractions, reads

$h(a)=1+B_{self}\phi+\phi\int(B_{inter}-B_{self})v\left(a^{\prime}\right)da^{\prime},$
(17)

where $v(a)=\frac{a^{3}p(a)}{\int a^{3}p(a)da}$ is the normalised volume
fraction distribution corresponding to $p(a)$, and the coefficients $B_{self}$
and $B_{inter}$ are averages of the mobility matrices. For a two-parameter
distribution such as the log-normal, $v(a^{\prime})$ is a function of only two
parameters, for example the mean $\langle a\rangle$ and the standard deviation
$\Delta a$. One could therefore postulate a hindered settling function where
the inter-class interaction term and the intra-class interaction terms are
functions of $\phi,\langle a\rangle,\Delta a/\langle a\rangle$. The issue is
that the same functional form should “best fit” a wide range of particle
sizes. This is a search which could benefit from the artificial intelligence
methods (Zhang & Ma, 2020; El Hasadi & Padding, 2023; Wu & Zhang, 2023;
Siddani et al., 2024). The analysis we have provided indicate some constraints
on this search, for example the functional forms expected for $\phi\rightarrow
0$ or $\alpha\rightarrow 0$.

##  7 Conclusions

We quantified numerically the hindered settling function of non-Brownian,
dilute suspensions of polydisperse spheres with a log-normal size
distribution, considering the effects of the polydispersity parameter $\alpha$
and the volume fraction $\phi$.

The average settling velocity $\langle u_{z,i}\rangle$ of each particle size
class is found to decrease as either $\phi$ or $\alpha$ increases. The class-
averaged velocity $\langle u_{z,i}\rangle$ decays with $\phi$ or $\alpha$
faster for the smaller particles than the largest particles. For given
$\alpha$ and $\phi$, the hindered settling function $\langle
u_{z,i}\rangle/u_{St,i}$ decreases as the particle size decreases, indicating
that the settling of the smaller size classes is more hindered compared to
that of the larger particles.

The probability distribution functions of the horizontal and vertical
velocities of each size class tend to follow a Gaussian distribution, with the
probability distributions of horizontal velocities of different size classes
collapsing onto each other for a given $\alpha$ and $\phi$. The magnitude of
the horizontal and vertical velocity fluctuations for each size class
increases as $\phi$ or $\alpha$ increases, and appear to follow the
approximate scaling $u_{i}^{\prime}\sim u_{St,i}(a_{i}/\langle
a\rangle)^{-2}$. Our simulations for the log-normally dispersed system suggest
a value of about 3.5 for the anisotropy ratio between the vertical and the
horizontal velocity fluctuations. This value is comparable to the one observed
in our simulations for monodisperse or bidisperse suspensions.

To verify the accuracy of available polydisperse hindered settling function
models, we compare the predictions of these models with the simulation data.
We found that the Richardson-Zaki correlation, which is often used to model
polydisperse suspensions, largely overestimates the hindered settling
functions of the smaller particles. For $\alpha$=0.4 and $\phi$=0.05, the
value predicted by Richardson-Zaki’s formula for polydisperse suspensions can
be up to seven times larger than the simulated value. Batchelor’s model
(equation (1)) gives quite accurate predictions for all size classes when
$\phi\leq 0.05$, yielding discrepancies of the settling velocities that are
within 10% of the numerical results. The Davis & Gecol model (equation 2) and
the MLB model (equation (4)) give comparable predictions. Both these models
tend to overestimate the hindered settling function of the smaller particles.
The discrepancy between the models and the simulation data increases as
$\alpha$ or $\phi$ increases. In Sec. 6, we use Batchelor’s model to analyse
the conditions under which models parameterised on the total volume fraction
can produce predictions of acceptable accuracy.

Our simulations demonstrate that the modelling of polydispersed suspension is
still challenging even in the dilute limit. Practically usable models such as
Richardson-Zaki or MLB work reasonably well for large particles, but give
significant errors for small particles. This is a major obstacle in the size
fractionation by centrifugation, for example, where one is interested in the
precise estimation of the velocity of each particle class. In these
applications, the small particle fraction is often the most valuable (Backes,
2016).

Our results hold in the Stokes regime. For future work, particle-resolved
simulations based on the solution of the Navier-Stokes equation around each
particle (Uhlmann & Doychev, 2014; Fornari et al., 2016; Willen & Prosperetti,
2019; Yao et al., 2021; Shajahan & Breugem, 2023) could be used to evaluate
the first effect of fluid inertia on our low-Reynolds number observations. In
the presence of fluid or particle inertia, averaging over instantaneous random
configurations cannot be applied, because in the inertial case the particle
velocities depend on the history of the hydrodynamic forces. Thanks to
advances in computational power, it is however now possible to simulate tens
of thousands particles at finite Reynolds numbers (Breugem, 2012; Schwarzmeier
et al., 2023), and with appropriate time averaging it should therefore be
possible to obtain smooth statistics for at least some of the parameters
combinations we explore. The most interesting seem to be the extreme cases,
namely small deviations from uniformity and log-normals with a large variance
(such as our $\alpha=0.4$ case). In the presence of fluid inertia, strong
trapping of small particles in the wake of large particles is expected.

Experimental techniques such as X-ray radiography (Dulanjalee et al., 2020),
magnetic resonance imaging (Boyce et al., 2016), or optical experiments with
fluorescent particles (Snabre et al., 2009) could be used to measure the
velocity of the small particle fraction in polydisperse suspensions. Machine
learning techniques such as symbolic regression (Zhang & Ma, 2020; El Hasadi &
Padding, 2023; Wu & Zhang, 2023) could be used in combination with particle-
resolved simulations (Yao et al., 2022) to extend Batchelor’s model to higher
volume fractions, or to incorporate into MLB’s model information about the
moments of the particle size distribution.

[Acknowledgements]We thank Wim-Paul Breugem and Johan Padding for valuable
discussions. The computations are carried out on the Reynolds cluster in the
Process & Energy department in TU Delft.

[Funding]This project has received funding from the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation program
(Grant Agreement No. 715475, project FLEXNANOFLOW).

[Declaration of interests]The authors report no conflict of interest.

##  Appendix A Derivation of the slip velocity closure in the MLB model

A derivation of the MLB model is provided here to highlight the key
assumptions of the model, which was too concisely described in the original
papers (Masliyah, 1979; Lockett & Bassoon, 1979). Consider a homogeneous
polydisperse suspension with $m$ particulate classes. The radius and density
of the $j$-th class are $a_{j}$ and $\rho_{j}$, respectively, with
$j=1,2,…,m$. The density and dynamic viscosity of the fluid are $\rho_{f}$ and
$\mu$, respectively. Gravity is in the negative $z$ direction. Due to the
differences between the particle and the fluid densities, a macroscopic
pressure gradient $dp/dz$ along the height of the mixture is needed to balance
the excess weight of the particles. This pressure gradient drives the back
flow of the fluid during settling of the particles. Corresponding to this
pressure gradient, each particle experiences a buoyancy force $F_{\nabla
p}=(-dp/dz)V_{p}$, where $V_{p}$ is the volume of that particle. The total
force exerted on each particle by the fluid is given by $F_{\nabla p}$, by the
buoyancy force due to the undisturbed hydrostatic pressure gradient and by the
drag force due to the relative fluid-particle velocity difference.

The steady-state momentum equation for the fluid phase is

$\left(-\frac{dp}{dz}\right)(1-\phi)-\sum_{j=1}^{m}f_{d,j}=0,$ (18)

where $\phi$ is the total volume fraction, and $f_{d,j}$ is the volumetric
drag force density (drag per unit volume) exerted by the $j$-th particle
class. The steady-state particle momentum equation for the $j$-th particle
class is

$\left(-\frac{dp}{dz}\right)\phi_{j}+f_{d,j}-(\rho_{j}-\rho_{f})\phi_{j}g=0,$
(19)

where $\phi_{j}$ is the volume fraction of the $j$-th class. Using equations
(18) and (19) gives

$\frac{dp}{dz}=-\sum_{j=1}^{m}(\rho_{j}-\rho_{f})\phi_{j}g,$ (20)

and

$f_{d,j}=(\rho_{j}-\rho_{susp})\phi_{j}g,$ (21)

where $\rho_{susp}=(1-\phi)\rho_{f}+\sum_{i=1}^{m}\rho_{i}\phi_{i}$ is the
density of the suspension (see e.g. Xia et al. (2022) for the case $m=1$). The
predictive accuracy of equation (21) for small particles immersed in a
suspension of larger particles has been put into question (Poletto & Joseph,
1995; Rotondi et al., 2015).

To calculate the particle velocity, a constitutive equation relating relative
velocity to force must be postulated. The MLB model uses a linear law between
the drag force and the slip velocity $u_{slip,j}$ between the $j$-th particle
class and the average fluid velocity:

$f_{d,j}=-\beta_{j}u_{slip,j}.$ (22)

where $u_{slip,j}$ is defined as in equation (14). The friction coefficient
was calculated as $\beta_{j}=\frac{9\mu\phi_{j}C(\phi)}{2a_{j}^{2}}$. The case
$C=1$ corresponds to no influence of neighbouring particles on the drag force
exerted on a test particle (the factor $\phi_{j}$ is due to the fact that
$f_{d,j}$ is a force per unit volume). To model hydrodynamic interactions on
the drag force, the MLB model assumes $C(\phi)=(1-\phi)^{2-n}$, as for a
monodisperse case at the same total volume fraction (from Richardson-Zaki’s
correlation, the slip velocity in the monodisperse case is $u_{slip}=\langle
u_{p}\rangle-\langle u_{f}\rangle=\frac{\langle
u_{p}\rangle}{1-\phi}=u_{St}(1-\phi)^{n-1}$; equating (21) and (22) using this
slip velocity gives $C(\phi)=(1-\phi)^{2-n}$).

From (21) and (22), the slip velocity for the polydispersed case is

$u_{slip,j}=\frac{2a_{j}^{2}}{9\mu}(1-\phi)^{n-2}(\rho_{j}-\rho_{susp}).$ (23)

If all the particles have the same density,
$\rho_{susp}=(1-\phi)\rho_{f}+\phi\rho_{p}$ and
$\rho_{j}-\rho_{susp}=(1-\phi)(\rho_{p}-\rho_{f})$. In this case the slip
velocity simplifies to

$u_{slip,j}=u_{St,j}(1-\phi)^{n-1},$ (24)

where $u_{St,j}$ is the Stokes velocity of the $j$-th species. Using the
definition of the slip velocity and using mass continuity
$\sum_{j}\phi_{j}\langle u_{j}\rangle+(1-\phi)\langle u_{f}\rangle=0$ yields
equation (4).

It can be seen from the derivation that the main assumptions in MLB’s model
are embedded in equations (21) and (22).

## References

  * Abbas et al. (2006) Abbas, M., Climent, E., Simonin, O. & Maxey, M. R. 2006 Dynamics of bidisperse suspensions under Stokes flows: Linear shear flow and sedimentation. Phys. Fluids 18 (12).
  * Abeynaike et al. (2012) Abeynaike, A., Sederman, A.J., Khan, Y., Johns, M.L., Davidson, J.F. & Mackley, M.R. 2012 The experimental measurement and modelling of sedimentation and creaming for glycerol/biodiesel droplet dispersions. Chem. Eng. Sci. 79, 125–137.
  * Al-Naafa & Selim (1992) Al-Naafa, M.A. & Selim, M. S. 1992 Sedimentation of monodisperse and bidisperse hard-sphere colloidal suspensions. AIChE J. 38 (10), 1618–1630.
  * Backes (2016) Backes, C. et al. 2016 Production of highly monolayer enriched dispersions of liquid-exfoliated nanosheets by liquid cascade centrifugation. ACS Nano 10 (1), 1589–1601.
  * Batchelor (1972) Batchelor, G.K. 1972 Sedimentation in a dilute dispersion of spheres. J. Fluid Mech. 52 (2), 245–268.
  * Batchelor (1982) Batchelor, G.K. 1982 Sedimentation in a dilute polydisperse system of interacting spheres. part 1. general theory. J. Fluid Mech. 119, 379–408.
  * Batchelor & Wen (1982) Batchelor, G.K. & Wen, C.S. 1982 Sedimentation in a dilute polydisperse system of interacting spheres. part 2. numerical results. J. Fluid Mech. 124, 495–528.
  * Beenakker (1986) Beenakker, C.W.J. 1986 Ewald sum of the rotne–prager tensor. J. Chem. Phys. 85 (3), 1581–1582.
  * Berres et al. (2005) Berres, S., Bürger, R. & Tory, E. M. 2005 Applications of polydisperse sedimentation models. Chem. Eng. J. 111 (2-3), 105–117.
  * Bonaccorso et al. (2013) Bonaccorso, F., Zerbetto, M., Ferrari, A. C. & Amendola, V. 2013 Sorting nanoparticles by centrifugal fields in clean media. J. Phys. Chem. C 117 (25), 13217–13229.
  * Boyce et al. (2016) Boyce, C. M., Rice, N.P., Ozel, A., Davidson, J.F., Sederman, A.J., Gladden, L.F., Sundaresan, S., Dennis, J.S. & Holland, D.J. 2016 Magnetic resonance characterization of coupled gas and particle dynamics in a bubbling fluidized bed. Phys. Rev. Fluids 1 (7), 074201.
  * Brady & Bossis (1988) Brady, J. F. & Bossis, G. 1988 Stokesian dynamics. Annu. Rev. Fluid Mech. 20 (1), 111–157.
  * Brady & Durlofsky (1988) Brady, J. F. & Durlofsky, L. J. 1988 The sedimentation rate of disordered suspensions. Phys. Fluids 31 (4), 717–727.
  * Brady et al. (1988) Brady, J. F., Phillips, R. J., Lester, J. C. & Bossis, G. 1988 Dynamic simulation of hydrodynamically interacting suspensions. J. Fluid Mech. 195, 257–280.
  * Breugem (2012) Breugem, W. 2012 A second-order accurate immersed boundary method for fully resolved simulations of particle-laden flows. J. Comput. Phys. 231 (13), 4469–4498.
  * Brzinski III & Durian (2018) Brzinski III, T.A. & Durian, D.J. 2018 Observation of two branches in the hindered settling function at low Reynolds number. Phys. Rev. Fluids 3 (12), 124303.
  * Bürger et al. (2002) Bürger, R., Karlsen, K. H., Tory, E.M. & Wendland, W.L. 2002 Model equations and instability regions for the sedimentation of polydisperse suspensions of spheres. Z Angew Math Mech. 82 (10), 699–722.
  * Chaturvedi et al. (2018) Chaturvedi, S. K., Ma, J., Brown, P. H., Zhao, H. & Schuck, P. 2018 Measuring macromolecular size distributions and interactions at high concentrations by sedimentation velocity. Nat. Commun. 9 (1), 4415.
  * Chen et al. (2023) Chen, H., Jia, X., Fairweather, M. & Hunter, T. N. 2023 Characterising the sedimentation of bidisperse colloidal silica using analytical centrifugation. Adv. Powder Technol. 34 (2), 103950.
  * Cheng & Wachs (2023) Cheng, Z. & Wachs, A. 2023 Hydrodynamic force and torque fluctuations in a random array of polydisperse stationary spheres. Int. J. Multiph. Flow 167, 104524.
  * Cunha et al. (2002) Cunha, F.R., Abade, G.C., Sousa, A.J. & Hinch, E.J. 2002 Modeling and direct simulation of velocity fluctuations and particle-velocity correlations in sedimentation. J. Fluids Eng. 124 (4), 957–968.
  * Davis & Birdsell (1988) Davis, R.H. & Birdsell, K.H. 1988 Hindered settling of semidilute monodisperse and polydisperse suspensions. AIChE J. 34 (1), 123–129.
  * Davis & Acrivos (1985) Davis, R. H. & Acrivos, A. 1985 Sedimentation of noncolloidal particles at low Reynolds numbers. Annu. Rev. Fluid Mech. 17 (1), 91–118.
  * Davis & Gecol (1994) Davis, R. H. & Gecol, H. 1994 Hindered settling function with no empirical parameters for polydisperse suspensions. AIChE J. 40 (3), 570–575.
  * Davis & Hassen (1988) Davis, R. H. & Hassen, M. A. 1988 Spreading of the interface at the top of a slightly polydisperse sedimenting suspension. J. Fluid Mech. 196, 107–134.
  * Di Felice (1995) Di Felice, R. 1995 Hydrodynamics of liquid fluidisation. Chem. Eng. Sci. 50 (8), 1213–1245.
  * Di Vaira et al. (2022) Di Vaira, N. J., Łaniewski-Wołłk, Ł., Johnson, R. L., Aminossadati, S. M. & Leonardi, C. R. 2022 Influence of particle polydispersity on bulk migration and size segregation in channel flows. J. Fluid Mech. 939, A30.
  * Dörr et al. (2013) Dörr, A., Sadiki, A. & Mehdizadeh, A. 2013 A discrete model for the apparent viscosity of polydisperse suspensions including maximum packing fraction. J. Rheol. 57 (3), 743–765.
  * Dorrell & Hogg (2010) Dorrell, R. & Hogg, A. J. 2010 Sedimentation of bidisperse suspensions. Int. J. Multiph. Flow 36 (6), 481–490.
  * Dulanjalee et al. (2020) Dulanjalee, E., Guillard, F., Baker, J., Einav, I. & Marks, B. 2020 Measuring grain size fractions of bidisperse granular materials using x-ray radiography. Opt. Express 28 (20), 29202–29211.
  * El Hasadi & Padding (2023) El Hasadi, Y. M. F & Padding, J. T. 2023 Do logarithmic terms exist in the drag coefficient of a single sphere at high reynolds numbers? Chem. Eng. Sci. 265, 118195.
  * Fornari et al. (2016) Fornari, W., Picano, F. & Brandt, L. 2016 Sedimentation of finite-size spheres in quiescent and turbulent environments. J. Fluid Mech. 788, 640–669.
  * Gibb (1991) Gibb, J. 1991 Pressure and viscous forces in an equilibrium fluidized suspension. Chem. Eng. Sci. 46 (1), 377–379.
  * Gonzalez et al. (2021) Gonzalez, E., Aponte-Rivera, C. & Zia, R. N. 2021 Impact of polydispersity and confinement on diffusion in hydrodynamically interacting colloidal suspensions. J. Fluid Mech. 925, A35.
  * Hase & Powell (2001) Hase, K. R. & Powell, R. L. 2001 Calculation of the Ewald summed far-field mobility functions for arbitrarily sized spherical particles in stokes flow. Phys. Fluids 13 (1), 32–44.
  * Hasimoto (1959) Hasimoto, H. 1959 On the periodic fundamental solutions of the Stokes equations and their application to viscous flow past a cubic array of spheres. J. Fluid Mech. 5 (2), 317–328.
  * Hayakawa & Ichiki (1995) Hayakawa, H. & Ichiki, K. 1995 Statistical theory of sedimentation of disordered suspensions. Phys. Rev. E 51 (5), R3815.
  * He et al. (2021) He, W., Wang, Q., Zhu, Y., Wang, K., Mao, J., Xue, X. & Shi, Y. 2021 Innovative technology of municipal wastewater treatment for rapid sludge sedimentation and enhancing pollutants removal with nano-material. Bioresour. Technol. 324, 124675.
  * Hinch (1988) Hinch, E.J. 1988 Sedimentation of small particles. In Disorder and Mixing: Convection, Diffusion and Reaction in Random Materials and Processes, pp. 153–162. Springer.
  * van der Hoef et al. (2005) van der Hoef, M. A., Beetstra, R. & Kuipers, J.A.M. 2005 Lattice-Boltzmann simulations of low-Reynolds-number flow past mono-and bidisperse arrays of spheres: results for the permeability and drag force. J. Fluid Mech. 528, 233–254.
  * Howard et al. (2022) Howard, A. A., Maxey, M. R. & Gallier, S. 2022 Bidisperse suspension balance model. Phys. Rev. Fluids 7 (12), 124301.
  * Kim & Karrila (2013) Kim, S. & Karrila, S. J. 2013 Microhydrodynamics: principles and selected applications. Butterworth-Heinemann.
  * Ladd (1996) Ladd, A. J. C. 1996 Hydrodynamic screening in sedimenting suspensions of non-Brownian spheres. Phys. Rev. Lett. 76 (8), 1392.
  * Lavrenteva et al. (2024) Lavrenteva, O.M., Smagin, I. & Nir, A. 2024 Shear-induced particle migration in viscous suspensions with continuous size distribution. Phys. Rev. Fluids 9 (2), 024305.
  * Lockett & Bassoon (1979) Lockett, M. J. & Bassoon, K.S. 1979 Sedimentation of binary particle mixtures. Powder Technol. 24 (1), 1–7.
  * Malbranche et al. (2023) Malbranche, N., Chakraborty, B. & Morris, J. F. 2023 Shear thickening in dense bidisperse suspensions. J. Rheol. 67 (1), 91–104.
  * Masliyah (1979) Masliyah, J. H. 1979 Hindered settling in a multi-species particle system. Chem. Eng. Sci. 34 (9), 1166–1168.
  * Moncho-Jordá et al. (2010) Moncho-Jordá, A., Louis, A.A. & Padding, J.T. 2010 Effects of interparticle attractions on colloidal sedimentation. Phys. Rev. Lett. 104 (6), 068301.
  * Mwasame et al. (2016) Mwasame, P. M., Wagner, N. J. & Beris, A. N. 2016 Modeling the effects of polydispersity on the viscosity of noncolloidal hard sphere suspensions. J. Rheol. 60 (2), 225–240.
  * Nguyen & Ladd (2005) Nguyen, N. & Ladd, A. J. C. 2005 Sedimentation of hard-sphere suspensions at low Reynolds number. J. Fluid Mech. 525, 73–104.
  * Padding & Louis (2004) Padding, J. T. & Louis, A. A. 2004 Hydrodynamic and Brownian fluctuations in sedimenting suspensions. Phys. Rev. Lett. 93 (22), 220601.
  * Papuga et al. (2021) Papuga, K., Kaszubkiewicz, J. & Kawałko, D. 2021 Do we have to use suspensions with low concentrations in determination of particle size distribution by sedimentation methods? Powder Technol. 389, 507–521.
  * Pednekar et al. (2018) Pednekar, S., Chun, J. & Morris, J. F. 2018 Bidisperse and polydisperse suspension rheology at large solid fraction. J. Rheol. 62 (2), 513–526.
  * Peysson & Guazzelli (1999) Peysson, Y. & Guazzelli, E. 1999 Velocity fluctuations in a bidisperse sedimenting suspension. Phys. Fluids 11 (7), 1953–1955.
  * Phillips et al. (1988) Phillips, R. J., Brady, J. F. & Bossis, G. 1988 Hydrodynamic transport properties of hard-sphere dispersions. i. suspensions of freely mobile particles. Phys. Fluids 31 (12), 3462–3472.
  * Piazza et al. (2013) Piazza, R., Buzzaccaro, S., Secchi, E. & Parola, A. 2013 On the general concept of buoyancy in sedimentation and ultracentrifugation. Phys. Biol. 10 (4), 045005.
  * Poletto & Joseph (1995) Poletto, M. & Joseph, D. D. 1995 Effective density and viscosity of a suspension. J. Rheol. 39 (2), 323–343.
  * Rettinger et al. (2022) Rettinger, C., Eibl, S., Rüde, U. & Vowinckel, B. 2022 Rheology of mobile sediment beds in laminar shear flow: effects of creep and polydispersity. J. Fluid Mech. 932, A1.
  * Revay & Higdon (1992) Revay, J. M. & Higdon, J. J. L. 1992 Numerical simulation of polydisperse sedimentation: equal-sized spheres. J. Fluid Mech. 243, 15–32.
  * Richardson & Zaki (1954) Richardson, J.F. T. & Zaki, W.N. 1954 Sedimentation and fluidization: part 1. Trans. Inst. Chem. Eng. 32, 35–53.
  * Rotne & Prager (1969) Rotne, J. & Prager, S. 1969 Variational treatment of hydrodynamic interaction in polymers. J. Chem. Phys. 50 (11), 4831–4837.
  * Rotondi et al. (2015) Rotondi, M., Di Felice, R. & Pagliai, P. 2015 Validation of fluid–particle interaction force relationships in binary-solid suspensions. Particuology 23, 40–48.
  * Ruzicka (2006) Ruzicka, M. C. 2006 On buoyancy in dispersion. Chem. Eng. Sci. 61 (8), 2437–2446.
  * Sarkar et al. (2009) Sarkar, S., van der Hoef, M. A. & Kuipers, J. A. M. 2009 Fluid–particle interaction from lattice Boltzmann simulations for flow through polydisperse random arrays of spheres. Chem. Eng. Sci. 64 (11), 2683–2691.
  * Schwarzmeier et al. (2023) Schwarzmeier, C., Rettinger, C., Kemmler, S., Plewinski, J., Núñez-González, F., Köstler, H., Rüde, U. & Vowinckel, B. 2023 Particle-resolved simulation of antidunes in free-surface flows. J. Fluid Mech. 961, R1.
  * Segre et al. (1997) Segre, P. N., Herbolzheimer, E. & Chaikin, P. M. 1997 Long-range correlations in sedimentation. Phys. Rev. Lett. 79 (13), 2574.
  * Shajahan & Breugem (2023) Shajahan, T. & Breugem, W. 2023 Inertial effects in sedimenting suspensions of solid spheres in a liquid. Int. J. Multiph. Flow 166, 104498.
  * Siddani et al. (2024) Siddani, B., Balachandar, S., Zhou, J. & Subramaniam, S. 2024 Investigating the influence of particle distribution on force and torque statistics using hierarchical machine learning. AIChE J. p. e18339.
  * Snabre et al. (2009) Snabre, P., Pouligny, B., Metayer, C. & Nadal, F. 2009 Size segregation and particle velocity fluctuations in settling concentrated suspensions. Rheol. Acta 48, 855–870.
  * Squires & Mason (2010) Squires, T. M. & Mason, T. G. 2010 Fluid mechanics of microrheology. Annu. Rev. Fluid Mech. 42, 413–438.
  * Tee et al. (2002) Tee, S., Mucha, P. J., Cipelletti, L., Manley, S., Brenner, M. P., Segre, P. N. & Weitz, D. A. 2002 Nonuniversal velocity fluctuations of sedimenting particles. Phys. Rev. Lett. 89 (5), 054501.
  * Tory et al. (2003) Tory, E. M., Karlsen, K. H., Bürger, R. & Berres, S. 2003 Strongly degenerate parabolic-hyperbolic systems modeling polydisperse sedimentation with compression. SIAM J. Appl. Math. 64 (1), 41–80.
  * Uhlmann & Doychev (2014) Uhlmann, M. & Doychev, T. 2014 Sedimentation of a dilute suspension of rigid spheres at intermediate Galileo numbers: the effect of clustering upon the particle motion. J. Fluid Mech. 752, 310–348.
  * Vowinckel et al. (2019) Vowinckel, B., Withers, J., Luzzatto-Fegiz, P. & Meiburg, E. 2019 Settling of cohesive sediment: particle-resolved simulations. J. Fluid Mech. 858, 5–44.
  * Wacholder & Sather (1974) Wacholder, E. & Sather, N. F. 1974 The hydrodynamic interaction of two unequal spheres moving under gravity through quiescent viscous fluid. J. Fluid Mech. 65 (3), 417–437.
  * Wagner & Woutersen (1994) Wagner, N. J. & Woutersen, A. T. J. M. 1994 The viscosity of bimodal and polydisperse suspensions of hard spheres in the dilute limit. J. Fluid Mech. 278, 267–287.
  * Wang & Brady (2015) Wang, M. & Brady, J. F. 2015 Short-time transport properties of bidisperse suspensions and porous media: A Stokesian dynamics study. J. Chem. Phys. 142 (9).
  * Wang et al. (2020) Wang, X., Wang, J., Liu, H., Zhao, L., Wang, Y., Wu, X. & Liao, X. 2020 Improving the production efficiency of sweet potato starch using a newly designed sedimentation tank during starch sedimentation process. J. Food Process. Preserv. 44 (10), e14811.
  * Willen & Prosperetti (2019) Willen, D. P. & Prosperetti, A. 2019 Resolved simulations of sedimenting suspensions of spheres. Phys. Rev. Fluids 4 (1), 014304.
  * Wolf (2021) Wolf, A. et al. 2021 Centrifugation based separation of lithium iron phosphate (LFP) and carbon black for lithium-ion battery recycling. Chem. Eng. Process. 160, 108310.
  * Wu & Zhang (2023) Wu, C. & Zhang, Y. 2023 Enhancing the shear-stress-transport turbulence model with symbolic regression: A generalizable and interpretable data-driven approach. Phys. Rev. Fluids 8 (8), 084604.
  * Xia et al. (2022) Xia, Y., Yu, Z., Pan, D., Lin, Z. & Guo, Y. 2022 Drag model from interface-resolved simulations of particle sedimentation in a periodic domain and vertical turbulent channel flows. J. Fluid Mech. 944, A25.
  * Xue & Sun (2003) Xue, B. & Sun, Y. 2003 Modeling of sedimentation of polydisperse spherical beads with a broad size distribution. Chem. Eng. Sci. 58 (8), 1531–1543.
  * Yao et al. (2022) Yao, Y., Biegert, E., Vowinckel, B., Köllner, T., Meiburg, E., Balachandar, S., Criddle, C. S. & Fringer, O. B 2022 Particle-resolved simulations of four-way coupled, polydispersed, particle-laden flows. Int. J. Numer. Methods Fluids 94 (11), 1810–1840.
  * Yao et al. (2021) Yao, Y., Criddle, C. S. & Fringer, O. B. 2021 The effects of particle clustering on hindered settling in high-concentration particle suspensions. J. Fluid Mech. 920, A40.
  * Yin & Sundaresan (2009) Yin, X. & Sundaresan, S. 2009 Fluid-particle drag in low-Reynolds-number polydisperse gas–solid suspensions. AIChE J. 55 (6), 1352–1368.
  * Zhang & Ma (2020) Zhang, J. & Ma, W. 2020 Data-driven discovery of governing equations for fluid dynamics based on molecular simulation. J. Fluid Mech. 892, A5.
  * Zuk et al. (2014) Zuk, P. J., Wajnryb, E., Mizerski, K. A. & Szymczak, P. 2014 Rotne–Prager–Yamakawa approximation for different-sized particles in application to macromolecular bead models. J. Fluid Mech. 741, R5.