2.1. Asumptions and simplifications

The starting point for our simplified theory, presented in § 2.2, is a Boltzmann-type equation for particle transport in a dilute polydisperse suspension of non-Brownian, neutrally buoyant particles undergoing unidirectional flow. Only cross-flow variations of particle concentrations are considered. Wall-induced particle migration is neglected under the assumption that the scale of the flow field is large compared with the particle size. Restricting consideration to stationary particle distributions, the collision fluxes of particle species are set to zero.

The collision flux is expanded for smoothly varying particle distributions in § 2.3 to obtain formulas for transport coefficients describing diffusive and drift fluxes of particle species driven, respectively, by gradients of particle concentrations and shear rate. The formulas are appropriately modified to handle regions with vanishing shear rate, as occurs in pressure-driven flows. The results of this section are combined in § 3 to yield an exact formulation for stationary particle distributions in dilute suspensions.

2.2. Boltzmann equation

We consider particle transport in two-dimensional unidirectional flows,

(2.1)\begin{equation} \boldsymbol{v}={v}(X_2)\boldsymbol{e_1} , \end{equation}

where $(X_1,X_2,X_3)$ describes a Cartesian coordinate system and $\boldsymbol {e_k} (k=1,2,3)$ are the corresponding unit vectors. Velocity field (2.1) includes simple shear and planar Poiseuille flow. The particle distribution evolves in the plane perpendicular to the fluid velocity according to the Boltzmann-type equation

(2.2)\begin{equation} \frac{\partial n_i}{\partial t}={-}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{F}_i, \end{equation}

where $n_i(X_2,X_3)$ is the number density of type-$i$ particles ($i=1,2,\ldots, m$). (Explicit time dependence of number density is suppressed here because we investigate stationary particle distributions.) The quantity $F_{ik}(X_2,X_3) (k=2,3)$ is the flux of type-$i$ particles in the $k$-direction resulting from pairwise ‘collisions’ with other particles,

(2.3)\begin{equation} F_{ik}(X_2,X_3)=\sum_{j=1}^m F_{ijk}. \end{equation}

Here, $F_{ijk}$ is the contribution to the flux $F_{ik}$ from collisions with type-$j$ particles ($j=1,2,\ldots, m$) given by the Boltzmann collision integral,

(2.4)\begin{equation} F_{ijk}=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\, {{\rm d}\kern0.7pt x}^{-\infty}_{2}\, {{\rm d}\kern0.7pt x}^{-\infty}_{3} \int_{-{\rm \Delta} X^{ij}_{k}}^0 n_{i}(X^{i,-\infty}_{k})n_{j}(X^{j,-\infty}_{k}) \vert {v}^{ij}\vert \, {\rm d} X^{i,-\infty}_{k}, \end{equation}

where $X^{i,-\infty }_{k} (k=2,3)$ is the initial distance of particle $i$ from the plane $X_{k}$ where the flux is evaluated, ${\rm \Delta} X^{ij}_{k}$ is the displacement of particle $i$ in the $k$-direction due to its binary encounter with particle $j$, and

(2.5)\begin{equation} x^{-\infty}_{k}=X^{j,-\infty}_{k}-X^{i,-\infty}_{k}, \end{equation}

is the relative initial position of the particles in the cross-flow plane, i.e. the trajectory offset, as shown in figure 1. Here, $\vert {v}^{ij} \vert$ is the magnitude of the relative velocity between widely separated particles,

(2.6)\begin{equation} \vert {v}^{ij} \vert = \vert {v}(X^{j,-\infty}_2) -{v}(X^{i,-\infty}_2)\vert, \end{equation}

and ${v}(X_2)$ is the prescribed velocity (2.1). Prior to a pair interaction, particles are widely separated in the $X_1$-direction (i.e. flow direction) with uncorrelated initial positions $(X^{i,-\infty }_{2}, X^{i,-\infty }_{3})$ and $(X^{j,-\infty }_{2}, X^{j,-\infty }_{3})$ in the cross-flow plane.

Figure 1. Schematic showing trajectories of particles undergoing pair interaction projected on the 1–$k$ plane ($k=2,3$); cross-flow displacements ${\rm \Delta} X^{ij}_k$ and ${\rm \Delta} X^{ji}_k$ of particles $i$ and $j$, respectively; $X^{i,-\infty }_k$ and $X^{j,-\infty }_k$ are coordinates of the particles in cross-flow plane ($X_2,X_3$) prior to interaction (dashed lines) and $x^{-\infty }_{k}=X^{j,-\infty }_{k}-X^{i,-\infty }_{k}$ is the initial trajectory offset; particle flux is evaluated at plane of constant $X_k$ (solid line).

Formula (2.4) is obtained using the odd symmetry of particle displacements with respect to trajectory offset,

(2.7)\begin{equation} {\rm \Delta} X^{ij}_{k}(- x^{-\infty}_{k})={-}{\rm \Delta} X^{ij}_{k} (x^{-\infty}_{k}), \end{equation}

for $(k=2,3)$. Particle displacements are, moreover, symmetric with respect to complementary coordinates, i.e.

(2.8a,b)\begin{equation} {\rm \Delta} X^{ij}_2({-}x^{-\infty}_3)={\rm \Delta} X^{ij}_2(x^{-\infty}_3), \quad {\rm \Delta} X^{ij}_3({-}x^{-\infty}_2)={\rm \Delta} X^{ij}_3(x^{-\infty}_2). \end{equation}

Cross-flow convection due to particle migration away from solid boundaries is omitted from (2.2) under the assumption that flow occurs in a channel that is wide compared with particle size. In the absence of wall interactions that can produce cross-swapping trajectories (Zurita-Gotor et al. Reference Zurita-Gotor, Bławzdziewicz and Wajnryb2012; Reddig & Stark Reference Reddig and Stark2013), particle displacements obey the relation

(2.9)\begin{equation} x^{-\infty}_{k} {\rm \Delta} X^{ij}_{k} \leq 0. \end{equation}

The equality holds only for $x^{-\infty }_{k}=0$ and for widely separated particles.

2.2.1. Boundary conditions

Stationary particle distributions are governed by

(2.10)\begin{equation} \boldsymbol{F}_i=0, \end{equation}

under the assumption that the channel walls are impermeable to the particles. In channel flows, a bulk number density, $n_{i\infty }$, can be imposed for each particle type $(i=1,2,\ldots,m)$,

(2.11)\begin{equation} q n_{i\infty}=\int_{S_\perp} {v}(X_2) n_i(X_2, X_3) \, {\rm d} X_2 \, {\rm d} X_3, \end{equation}

where $q$ is the volume flow through a channel with cross-section $S_\perp$.

Diluteness requires low volume fractions,

(2.12)\begin{equation} \phi_{i\infty}\ll 1, \end{equation}

where number densities are related to volume fractions, $\phi _i=n_i {v}_i$, and ${v}_i$ is particle volume. Stationary particle distributions depend nonlinearly on the bulk composition $n_{i\infty }/n_\infty$ given that the particle flux (2.4) is quadratic in number density. However, the total bulk density $n_\infty =\sum _{i=1}^m n_{i\infty }$ can be scaled out of (2.10) and (2.11). Stationary particle density distributions are proportional to the total bulk number density.

2.3. Particle transport

2.3.1. Local transport coefficients

Here, we analyse particle transport in suspensions undergoing flows with velocity (2.1) and derive diffusive and drift fluxes, i.e. fluxes proportional to the local concentration gradients of particle species and a flux associated with a drift velocity of particles from high to low shear rates. Local transport coefficients for these fluxes are obtained. The diffusive fluxes have been derived previously (Zarraga & Leighton Reference Zarraga and Leighton2001) and a drift velocity has been extracted from a Boltzmann collision integral (Rivera et al. Reference Rivera, Zhang and Graham2016). A brief derivation is provided below for completeness and uniformity.

The diffusive and drift fluxes are obtained by evaluating the collision flux (2.4) for perturbative variations in number densities and relative velocities. The number densities and shear rates are expanded up to linear variations

(2.13)\begin{gather} n_i=n_{i 0}+\sum_{k=2}^{3} \frac{\partial n_i}{\partial X_k}X_k+O(X^2_k), \end{gather}

(2.14)\begin{gather}\dot{\gamma}=\dot{\gamma}_0 + \frac{{\rm d}\dot{\gamma}}{{\rm d} X_2}X_2+O(X_2^2), \end{gather}

where $\dot {\gamma }_0$ and $n_{i 0}$, respectively, are the local shear-rate magnitude and number densities of the particles at $X_k=0$ ($k=2,3$); the quantities $d\dot {\gamma }/d X_2$ and $\partial n_i/\partial X_k$ are the corresponding local values of the derivatives. Brownian motion is neglected under the assumption of large Péclet number

(2.15)\begin{equation} \frac{\dot{\gamma}_0 a^2}{\mathcal{D}}\gg 1, \end{equation}

where $\mathcal {D}$ is the molecular diffusivity. In general, the relative velocity (2.6) becomes

(2.16)\begin{equation} \lvert {v}^{ij}\rvert = \lvert x^{-\infty}_2\rvert \left\vert \dot{\gamma}_0+ \frac{{\rm d}\dot{\gamma}}{{\rm d} X_2}\left[X_2^{i, -\infty}+ \frac{1}{2}x^{-\infty}_2\right]\right\vert. \end{equation}

The development presented in this section assumes a non-vanishing shear rate

(2.17)\begin{equation} \dot{\gamma}_0 > \left\vert \frac{{\rm d}\dot{\gamma}}{{\rm d} X_2}\right\vert \left\vert X_2^{i, -\infty}+\frac{1}{2}x^{-\infty}_2\right\vert , \end{equation}

so that the relative velocity (2.16) reduces to

(2.18)\begin{equation} \lvert {v}^{ij}\rvert = \lvert x^{-\infty}_2\rvert \left( \dot{\gamma}_0+ \frac{{\rm d}\dot{\gamma}}{{\rm d} X_2}\left[X_2^{i, -\infty}+ \frac{1}{2}x^{-\infty}_2\right]\right). \end{equation}

The case of vanishing shear rates is analysed in the next section.

Inserting (2.13) and (2.14) into the flux (2.4) and integrating in $X_{k}^{i,-\infty }$ yields the net flux of type-$i$ particles,

(2.19)\begin{equation} F_{ik}=F^{(c)}_{ik}+\delta_{k2} F^{(\dot{\gamma})}_{i2}, \end{equation}

where $F^{(c)}_{ik}$ are the diffusive fluxes

(2.20)\begin{equation} F^{(c)}_{ik}={-}D_{ik}^{s}\frac{{\rm d} n_{i}}{{\rm d} X_{k}}- \sum_{j=1}^m\left( D_{ijk}\frac{{\rm d} n_{j}}{{\rm d} X_{k}}\right), \end{equation}

and $F^{(\dot {\gamma })}_{i2}$ is the drift flux

(2.21)\begin{equation} F^{(\dot{\gamma})}_{i2}={-}n_{i0} V_{i} \frac{{\rm d} \dot{\gamma}}{{\rm d} X_2}. \end{equation}

The Kronecker delta $\delta _{pq}$ in (2.19) indicates that there is a drift flux in the $X_2-$direction only.

The above diffusivities and the drift-velocity coefficient are defined by

(2.22)\begin{gather} D_{ik}^{s}=\dot{\gamma}_0\sum_{j=1}^m n_{j0} I^{A}_{ijk}, \end{gather}

(2.23)\begin{gather}D_{ijk}=\dot{\gamma}_0 n_{i0} ( I^{A}_{ijk}+I^{B}_{ijk}), \end{gather}

(2.24)\begin{gather}V_{i}=\sum_{j=1}^m n_{j0} \left(I^{A}_{ij2}+ \frac{1}{2}I^{B}_{ij2}\right), \end{gather}

where $I^A_{ijk}$, $I^B_{ijk}$ are integrals over the relative cross-flow-plane coordinates (2.5)

(2.25)\begin{gather} I^{A}_{ijk}=\frac{1}{2}\int^{\infty}_{-\infty} \int^{\infty}_{-\infty}\vert x^{-\infty}_{2}\vert ({\rm \Delta} X^{ij}_{k})^2 \, {{\rm d}\kern0.7pt x}^{-\infty}_{2}\, {{\rm d}\kern0.7pt x}^{-\infty}_{3}, \end{gather}

(2.26)\begin{gather}I^{B}_{ijk}={-}\int^{\infty}_{-\infty} \int^{\infty}_{-\infty}|x^{-\infty}_{2}|x^{-\infty}_{k} {\rm \Delta} X^{ij}_{k} \, {{\rm d}\kern0.7pt x}^{-\infty}_{2}\, {{\rm d}\kern0.7pt x}^{-\infty}_{3}. \end{gather}

Both integrals are intrinsically positive given the symmetry (2.8a,b) and sign of particle displacements (2.9).

The self-diffusivity of type-$i$ particles, $D_{ik}^{s}$, defined by (2.22), can be directly obtained as a sum of the rate of mean squared displacements from random encounters with all particle species. Diffusive fluxes occur in both the velocity gradient and vorticity directions and have contributions from concentration gradients of all species; a non-zero diffusive flux of a species with uniform concentration can be generated by a gradient of another species. This formulation of the diffusive fluxes concurs with that presented by Zarraga & Leighton (Reference Zarraga and Leighton2001).

The drift velocities describe the migration of particles from regions of high shear rates. Gradients of the shear-rate magnitude $\dot {\gamma }$ generate an oppositely directed flux. By symmetry, gradients of the shear-rate magnitude do not contribute to the diffusive flux.

The $O(r^{-3})$ far-field velocity gradients of force-free particles perturb the trajectories of deformable particles, where $r$ is distance between the particles normalized by particle size. Integrated along an open trajectory, this perturbation produces pair displacements ${\rm \Delta} X_k^{ij}=O(r^{-\infty })^{-2}$ for widely offset trajectories, $r^{-\infty } \gg 1$, where $r^{-\infty }$ is the radial relative trajectory coordinate, defined in figure 2. Accordingly, integral (2.26) is divergent. The self-diffusivity of deformable particles can be computed from pair interactions because integral (2.25) converges (Loewenberg & Hinch Reference Loewenberg and Hinch1997), but the evaluation of the diffusive and drift fluxes (2.20) and (2.21) requires a numerical cutoff, $r_c^{ij}=r_c^{ji}$, a trajectory offset beyond which interactions between type-$i$ and type-$j$ particles cause cross-flow displacements, i.e.

(2.27a,b)\begin{equation} {\rm \Delta} X^{ij}_k=0, \quad r^{-\infty} > r_c^{ij} . \end{equation}

Figure 2. Cylindrical coordinate system $(r^{-\infty }, \theta )$ for the cross-flow plane $(x_3^{-\infty }, x_2^{-\infty })$, where $x_3^{-\infty }=r^{-\infty }\cos \theta$ and $x_2^{-\infty }=r^{-\infty }\sin \theta$.

With this cutoff, integrals (2.25) and (2.26) become

(2.28)\begin{gather} I^{A}_{ijk}=\frac{1}{2}\int_0^{2{\rm \pi}} \int^{r_c^{ij}}_{0}\vert r^{-\infty}\sin\theta\vert ({\rm \Delta} X^{ij}_{k})^2 r^{-\infty}\, {\rm d} r^{-\infty} \, {\rm d} \theta , \end{gather}

(2.29)\begin{gather}I^{B}_{ijk}={-}\int_0^{2{\rm \pi}} \int^{r_c^{ij}}_{0}|r^{-\infty}\sin\theta |x^{-\infty}_{k} {\rm \Delta} X^{ij}_{k} \, r^{-\infty}\, {\rm d} r^{-\infty} \, {\rm d} \theta, \end{gather}

where $(r^{-\infty }, \theta )$ is the cylindrical coordinate system shown in figure 2 and $x^{-\infty }_{k}$ is defined in the caption.

In general, truncation of the integration domain is an ad hoc procedure, and results depend on $r_c^{ij}$ (Narsimhan et al. Reference Narsimhan, Zhao and Shaqfeh2013). However, boundary condition (2.27a,b) applies rigorously for spherical particles that undergo contact interactions, as discussed in § 4. Such particles have circular upstream collision cross-sections $r_c^{ij}$ that scale with particle size (Zinchenko Reference Zinchenko1984; da Cunha & Hinch Reference da Cunha and Hinch1996).

2.3.2. Particle transport at points of vanishing shear rate

Here, the analysis of particle transport developed in the previous section is extended to regions where the shear rate vanishes linearly. Accordingly, we consider regions where the velocity is locally quadratic,

(2.30)\begin{equation} {v}={v}_0+\frac{\gamma'_2}{2} X_2^2, \end{equation}

where $\gamma '_2$ is the magnitude of the shear-rate gradient and ${v}_0$ is an arbitrary local velocity that can be ignored since only velocity differences are significant. The corresponding magnitude of the shear rate and its gradient are given by

(2.31a,b)\begin{equation} \dot{\gamma}_0 = \gamma'_2\, \lvert X_2\rvert, \quad \frac{{\rm d} \dot{\gamma}}{{\rm d} X_2}=\mbox{sign}(X_2)\gamma'_2. \end{equation}

An example is the velocity field at the centre of planar Poiseuille flow in a channel of half-width $H$, where $\boldsymbol {v}$ is given by (2.1) with

(2.32)\begin{equation} {v}={v}_0\left[1-\left(\frac{X_{2}}{H}\right)^2\right]. \end{equation}

Here, ${v}_0$ is the velocity at centreline, $X_2=0$, and the magnitude of the shear rate is given by (2.31a,b) with $\gamma '_2=2{v}_0 H^{-2}$. Neglecting Brownian motion in regions where the shear rate vanishes requires a more stringent condition; in this case, requirement (2.15) is replaced by

(2.33)\begin{equation} \frac{\gamma'_2 a^3}{\mathcal{D}} \gg 1. \end{equation}

Inserting (2.31a,b) into (2.16) yields the relative velocity magnitude,

(2.34)\begin{equation} \vert {v}^{ij}\vert=\dot{\gamma}'_{2} \vert x^{-\infty}_{2}\vert \vert X_{2}+X^{i,-\infty}_{2}+\tfrac{1}{2}x^{-\infty}_{2}\vert. \end{equation}

For $\lvert X_2\rvert > X_c^{ij}$, the result reduces to the form of (2.18),

(2.35)\begin{equation} \vert {v}^{ij}\vert=\dot{\gamma}'_{2} \vert x^{-\infty}_{2}\vert ( \lvert X_{2}\rvert +{\rm sign}(X_2)[X^{i,-\infty}_{2}+\tfrac{1}{2}x^{-\infty}_{2}]), \quad\lvert X_2\rvert > X_c^{ij}, \end{equation}

where $X_c^{ij}$ is an upper bound of the excluded region $\vert X^{i,-\infty }_{2}+\frac {1}{2}x^{-\infty }_{2}\vert$. Here, $X_c^{ij}$ is given by

(2.36)\begin{equation} X_c^{ij}= {\rm \Delta} X^{ij}_{2, {max}} +\tfrac{1}{2}r_c^{ij}, \end{equation}

where ${\rm \Delta} X^{ij}_{2, {max}}$ is the maximum displacement magnitude of a type-$i$ particle by a pair interaction with a type-$j$ particle and is thus the upper bound for $\vert X^{i,-\infty }_{2}\vert$; the radius of the collision cross-section, $r_c^{ij}$, defined by (2.27a,b), is the upper bound for $\lvert x_2^{-\infty }\rvert$. In general, ${\rm \Delta} X^{ij}_{2, {max}}$, and thus $X_c^{ij}$, are $O(r_c^{ij})$. For spherical particles that undergo contact interactions, $X_c^{ij} \leq r_c^{ij}$, as discussed in § 4.

For $\lvert X_2\rvert > X_c^{ij}$, the local analysis presented in § 2.3.1 applies with shear-rate magnitude and gradient given by (2.31a,b). Accordingly, a linearly varying diffusive flux and constant drift velocity are obtained according to (2.20)–(2.24). However, for $\lvert X_2\rvert < X_c^{ij}$, the relative velocity, ${v}_{ij}$, changes sign within the maximum range of particle displacements, ${\rm \Delta} X^{ij}_{2, {max}}$, that contribute to the particle flux. Within this region, (2.34) must be used to describe the magnitude of the relative velocity. The use of (2.35) is inconsistent and leads to incorrect results. For example, inserting the linearly varying diffusive flux balanced by the constant drift velocity, obtained by the local analysis (2.22)–(2.24) into (2.19) and then (2.10) yields,

(2.37)\begin{equation} X_2 n \frac{{\rm d} n}{{\rm d} X_2}={-}M_1 n^2, \end{equation}

where $M_1$ is a positive constant. A singular distribution, $n\approx \lvert X_2\vert ^{-M_1}$ is thus obtained, as pointed out by Leighton & Acrivos (Reference Leighton and Acrivos1987b).

The region $\lvert X_2\rvert < X_c^{ij}$ defines a boundary layer within which the transport coefficients exhibit a more complex dependence on position $X_2$. By retaining the relative velocity magnitude (2.34) in this region, we obtain the essential dependence of the transport coefficients required to avoid the spurious singularity at $X_2=0$ without the need for any of the ad hoc manoeuvres used in prior studies, as discussed in the introduction.

Particle transport in the vorticity direction ($k=3$) does not have a singular behaviour because there is no drift velocity, so the local analysis presented in § 2.3.1 is uniformly valid.

According to the derivation provided in Appendix A, the particle flux in the velocity-gradient direction is given by

(2.38)\begin{equation} F_{i2}(X_2)=F^{(c)}_{i2}(X_2)+ F^{(\dot{\gamma})}_{i2}(X_2), \end{equation}

where the diffusive and drift fluxes are

(2.39)\begin{gather} F^{(c)}_{i2}(X_2)={-}D_{i2}^{s}(X_2)\frac{{\rm d} n_{i}}{{\rm d} X_{2}}-\sum_{j=1}^m\left( D_{ij2}(X_2) \frac{{\rm d} n_{j}}{{\rm d} X_{2}}\right), \end{gather}

(2.40)\begin{gather}F^{(\dot{\gamma})}_{i2}(X_2)={-}\gamma'_2 n_{i0} V_{i} (X_2), \end{gather}

the diffusivities and drift-velocity coefficient are

(2.41)\begin{gather} D^{s}_{i2}(X_2)=\gamma'_2\lvert X_2\rvert \sum_{j=1}^{m} n_{j0} I_{ij}^{(1)}(X_2) , \end{gather}

(2.42)\begin{gather}D_{ij2}(X_2)=\gamma'_2\lvert X_2\rvert \, n_{i0} I_{ij}^{(2)}(X_2), \end{gather}

(2.43)\begin{gather}V_{i}(X_2)=\sum_{j=1}^m n_{j0} I^{(3)}_{ij}(X_2). \end{gather}

Following the analysis presented in Appendix A, integrals $I_{ij}^{(1)}, I_{ij}^{(2)}, I_{ij}^{(3)}$ are given by

(2.44)\begin{align} I_{ij}^{(1)}(X_2)&=\frac{1}{\vert X_2\vert} \int_{0}^{\rm \pi} \int^{r_c^{ij}}_{0} r^{-\infty}\sin \theta \left[ \frac{1}{2} \vert X_2\vert ({\rm \Delta} X^{ij}_{2})^2 B^2(X_2') +\frac{1}{3} (-{\rm \Delta} X^{ij}_{2})^3(1-B^3(X'_{2})) \right. \nonumber\\ & \quad \left. + \frac{1}{4} r^{-\infty} \sin \theta ({\rm \Delta} X^{ij}_{2})^2 (1-B^2(X'_{2})) \right] r^{-\infty} \, {\rm d} r^{-\infty} \, {\rm d} \theta ,\end{align}

(2.45)\begin{align} I_{ij}^{(2)}(X_2)&=\frac{1}{\vert X_2\vert} \int_{0}^{\rm \pi} \int^{r_{c}^{ij}}_{0} r^{-\infty}\sin \theta \left[ r^{-\infty} \sin \theta \vert X_2\vert (-{\rm \Delta} X^{ij}_{2}) B(X'_{2})+\frac{1}{2}\vert X_2\vert ({\rm \Delta} X^{ij}_{2})^2 B^2(X'_{2}) \right. \nonumber\\ & \quad \left. +\frac{1}{3}(-{\rm \Delta} X^{ij}_{2})^3(1-B^3(X'_{2})) + \frac{3}{4} r^{-\infty} \sin \theta ({\rm \Delta} X^{ij}_{2})^2 (1-B^2(X'_{2})) \right. \nonumber\\ & \quad \left. +\frac{1}{2}(r^{-\infty} \sin \theta)^2 (-{\rm \Delta} X^{ij}_{2})(1-B(X'_{2})) \right] r^{-\infty} \, {\rm d} r^{-\infty} \, {\rm d} \theta , \end{align}

and

(2.46)\begin{align} I_{ij}^{(3)}(X_2)&= \mbox{sign}(X_2)\, \int_{0}^{\rm \pi} \int^{r_{c}^{ij}}_{0} r^{-\infty}\sin \theta \left[ \frac{1}{2}({\rm \Delta} X^{ij}_{2})^2 B^2(X'_{2})+ \frac{1}{2} r^{-\infty}\sin \theta (-{\rm \Delta} X^{ij}_{2}) B(X'_{2}) \right. \nonumber\\ &\quad \left. \vphantom{\frac{1}{2}}+\vert X_2\vert(-{\rm \Delta} X^{ij}_{2}) (1-B(X'_{2})) \right] r^{-\infty} \, {\rm d} r^{-\infty} \, {\rm d} \theta , \end{align}

where $(r^{-\infty },\theta )$ are the cylindrical coordinates defined in figure 2, and

(2.47)\begin{equation} X_2'= \frac{\vert X_2\vert-\frac{1}{2}r^{-\infty}\sin\theta}{-{\rm \Delta} X^{ij}_{2}}. \end{equation}

Here, $B(t)$ is the ramp function sketched in figure 3 and defined by

(2.48)\begin{equation} B(t)=t \varTheta(t)-(t-1)\varTheta(t-1) , \end{equation}

where $\varTheta (t)$ is the Heaviside function. Note that

(2.49)\begin{equation} -{\rm \Delta} X^{ij}_2 > 0\, \quad\mbox{for } 0 < r^{-\infty}< r_{c}^{ij}, \quad 0 < \theta < {\rm \pi}, \end{equation}

by (2.9) and the remark below it. It follows that integrals (2.44)–(2.46), and thus the transport coefficients (2.41)–(2.43), are positive provided $X_2\neq 0$.

Figure 3. Ramp function defined by (2.48).

2.3.3. Discussion

Here, we point out the salient features of the diffusive and drift fluxes in the boundary-layer region.

For $\lvert X_2\rvert \geq X_c^{ij}$, we have $X_2'\geq 1$, according to (2.36) and (2.47), and thus

(2.50)\begin{equation} B(X_2')=1. \end{equation}

Inserting this result into integrals (2.44)–(2.46) reduces them to their local forms

(2.51a–c)\begin{equation} I^{(1)}_{ij}\to I^{A}_{ij2}, \quad I^{(2)}_{ij}\to I^{A}_{ij2}+ I^{B}_{ij2}, \quad I^{(3)}_{ij}\to \mbox{sign}(X_2) (I^{A}_{ij2}+ \tfrac{1}{2}I^{B}_{ij2}), \quad \lvert X_2\rvert \geq X^{ij}_c. \end{equation}

Accordingly, the transport coefficients (2.41)–(2.43) reduce to their corresponding local forms, (2.22)–(2.24), with the shear-rate magnitude and its gradient given by (2.31a,b).

First-derivative discontinuities introduced by the ramp function (2.48) are manifest by second-derivative discontinuities of integrals $I^{(k)}_{ij}(X_2)$ ($k=1,2,3$) and the transport coefficients (2.41)–(2.43) at $\lvert X_2\rvert =X_c^{ij}$. The source of these discontinuities is the discrete interaction range (2.27a,b) for contact interactions between spherical particles.

According to (2.44)–(2.46), integrals $I^{(1)}(X_2)$ and $I^{(2)}(X_2)$ are even functions of $X_2$ and $I^{(3)}(X_2)$ is an odd function. (This is plainly seen in the antecedent (A2)–(A4).) Thus, diffusivities (2.41)–(2.42) are even functions and drift velocities (2.43) are odd functions of $X_2$. Near the centreline, integrals (2.44)–(2.46) reduce to

(2.52)\begin{align} \lim_{X_2\to 0}\lvert X_2\rvert I_{ij}^{(1)}(X_2)&=\int_{0}^{\rm \pi} \int^{r_c^{ij}}_{0} r^{-\infty}\sin \theta \left[ \frac{1}{3} (-{\rm \Delta} X^{ij}_{2})^3 \right. \nonumber\\ & \quad \left. + \frac{1}{4} r^{-\infty} \sin \theta\, ({\rm \Delta} X^{ij}_{2})^2\right] r^{-\infty} \, {\rm d} r^{-\infty} \, {\rm d} \theta , \end{align}

(2.53)\begin{align} \lim_{X_2\to 0}\lvert X_2\rvert I_{ij}^{(2)}(X_2)&=\int_{0}^{\rm \pi} \int^{r_{c}^{ij}}_{0} r^{-\infty}\sin \theta \left[ \frac{1}{3}(-{\rm \Delta} X^{ij}_{2})^3 \right. \nonumber\\ & \quad \left. + \frac{3}{4} r^{-\infty} \sin \theta ({\rm \Delta} X^{ij}_{2})^2 +\frac{1}{2}(r^{-\infty} \sin \theta)^2 (-{\rm \Delta} X^{ij}_{2}) \right] r^{-\infty} \, {\rm d} r^{-\infty} \, {\rm d} \theta , \end{align}

and

(2.54)\begin{equation} \lim_{X_2\to 0} I_{ij}^{(3)}(X_2)= X_2\, \int_{0}^{\rm \pi} \int^{r_{c}^{ij}}_{0} r^{-\infty}\sin \theta(-{\rm \Delta} X^{ij}_{2}) r^{-\infty} \, {\rm d} r^{-\infty} \, {\rm d} \theta . \end{equation}

Inserting these results with (2.49) into (2.41)–(2.43) shows that diffusivities are finite and positive and drift velocities are linear with positive slope for $X_2\to 0$. Inserting these results into (2.38) and (2.10) indicates that particle transport near the centreline is governed by

(2.55)\begin{equation} n \frac{{\rm d} n}\, {{\rm d} y}={-}M_2 y n^2, \quad \lvert y\rvert \to 0, \end{equation}

where $M_2$ is a positive constant, and $y=X_2/X_c^{ij}$ is a dimensionless length variable. By contrast to (2.37), this equation predicts the smooth, non-singular behaviour, $n\approx n_0(1- \frac {1}{2} M_2\, y^2)$, where $n_0$ is the number density at $y=0$. Hence, the source of the classical singular particle distribution, described above, is eliminated. Note that $X_c^{ij}$ introduces a particle-related length scale on the distribution near the centreline, whereas (2.37) describes a distribution without a length scale.

Rivera et al. (Reference Rivera, Zhang and Graham2016) accommodated the non-vanishing diffusivities at the centreline but ignored variations in the drift velocity. This led to the prediction of a linear cusp at the centreline, as seen by inserting (41) from their paper into (2.10) to obtain

(2.56)\begin{equation} n \frac{{\rm d} n}{{\rm d} \tilde y}={-}M_3\, \mbox{sign}(\tilde y) n^2, \quad \lvert y\rvert \to 0, \end{equation}

which yields $n\approx n_0(1-M_3 \lvert \tilde y\rvert )$, consistent with the result shown in figure 10 of their paper. Here, $M_3$ is a positive constant, and $\tilde y=X_2/H$.

The diffusive flux model (Phillips et al. Reference Phillips, Armstrong, Brown, Graham and Abbott1992) also predicts a linear cusp at the centreline, except with $n_0=n_{m}$, corresponding the maximum packing fraction, as shown in Appendix B. Although the diffusive flux model uses local coefficients for diffusivity and drift velocity in the particle transport equation (B1), the singular behaviour (2.37) is avoided through the use of a suspension viscosity model that diverges at a prescribed maximum packing fraction.

3.1. Non-vanishing shear rates

Here, we present an exact stationary solution for the particle distribution in a polydisperse suspension in a flow with a power-law shear-rate magnitude,

(3.1)\begin{equation} \dot{\gamma} = C_1 X_2^\beta, \end{equation}

where $\beta$ and $C_1$ are arbitrary non-zero constants and $X_2>0$ is assumed. Accordingly, the shear rate is non-vanishing (and non-singular).

Inserting (2.19)–(2.24) into (2.10) yields

(3.2)\begin{equation} \sum_{j=1}^m (\dot{\gamma}[I_{ij2}^A\, n_{j} n_{i}' + (I_{ij2}^A+I_{ij2}^B)n_{i} n_{j}' ]+\dot{\gamma}'[n_{i}n_{j}(I_{ij2}^A+\tfrac{1}{2}I_{ij2}^B)])=0, \end{equation}

for $i=1,2,\ldots, m$, where primes are used to denote $X_2$-derivatives.

A power-law particle distribution

(3.3)\begin{equation} n_{i}(X_2)=c_i X_2^{-\beta/2} \end{equation}

is seen to satisfy (3.2) with arbitrary coefficients $c_i$. This is a general result that holds independent of the details of pairwise particle interactions in a given system. There are two features that should be noted here: (i) the effect of particle interactions exactly cancel, i.e. the spatial distribution of each particle species is unaffected by the presence of the others, and (ii) particle size does not affect the particle distribution. These features break down in regions where the shear rate vanishes, as seen below.

3.2. Planar Poiseuille flow

Here, we consider the steady-state particle distribution in quadratic flows (2.30), including regions $X_2\to 0$, where the shear rate vanishes.

3.2.1. Polydisperse suspension

Inserting (2.38)–(2.43) into (2.10) yields the equation governing the stationary particle distribution

(3.4)\begin{equation} \sum_{j=1}^m (\lvert X_2\rvert [I_{ij}^{(1)}(X_2) n_{j} n_{i}' + I_{ij}^{(2)}(X_2) n_{i} n_{j}' ]+I_{ij}^{(3)}(X_2) n_{i}n_{j})=0, \quad\vert X_2\vert < X_c, \end{equation}

for $i=1,2,\ldots, m$. In a polydisperse system, each binary interaction has a distinct boundary-layer half-thickness, $X_c^{ij}$, determined by (2.36), and $X_c$ is the maximum of these, i.e.

(3.5)\begin{equation} X_c=\sup\{X^{ij}_c\}. \end{equation}

For $\lvert X_2\rvert > X_c$, (3.4) reduces to  (3.2), according to (2.51a–c), thus particle distributions are decoupled and obey distribution (3.3) with $\beta =1$, i.e.

(3.6)\begin{equation} n_i(X_2) =n_{ic}X_c^{1/2}\lvert X_2\rvert ^{{-}1/2}, \quad \lvert X_2\rvert > X_c, \end{equation}

where $n_{i c}$ is the number density of species $i$ at $X_2=X_c$.

The spatial distributions of particle species are coupled for $\lvert X_2\rvert < X_c$. By contrast to the scale-free power-law distribution (3.3) that governs particle distributions in regions of non-vanishing shear rate, the coupling that occurs in regions of vanishing shear rate imposes a particle-related length scale $X_c$. A dimensionless coordinate is thus introduced using the length scale $X_c$

(3.7)\begin{equation} y=X_2/X_c. \end{equation}

It will also be useful to define the dimensionless number densities

(3.8a,b)\begin{equation} \bar n_i=n_i/n_{c}, \quad \bar N_i=n_i/n_{ic}, \end{equation}

where $n_{i c}$ is the number density of particle species $i$ at $X_2=X_c$ and $n_c=\sum _{i=1}^m n_{i c}$.

In terms of these variables, (3.4) becomes

(3.9)\begin{equation} \sum_{j=1}^m ( {D}_{ij}^{(1)}(y) \bar n_{j} \bar n_{i}' + {D}_{ij}^{(2)}(y) \bar n_{i} \bar n_{j}' +{V}_{ij}(y)\, \bar n_{i}\bar n_{j})=0, \quad \vert y\vert < 1, \end{equation}

where primes denote $y$-derivatives and dimensionless transport coefficients are defined

(3.10a,b)\begin{equation} {D}_{ij}^{(k)}(y) = \lvert y\rvert I^{(k)}_{ij}(y), \quad k=1,2 ;\quad {V}_{ij}(y) = I^{(3)}_{ij}(y). \end{equation}

Boundary conditions for (3.9) are given by

(3.11)\begin{equation} \bar n_i(1)=x_{i c}, \end{equation}

where $x_{i c}=n_{i c}/n_c$ is the number density fraction at $y=1$. Outside of the coupled boundary-layer region, $\lvert y\rvert > 1$, the dimensionless particle densities are given by

(3.12)\begin{equation} \bar N_i(y) =\lvert y\rvert ^{{-}1/2}, \end{equation}

according to (3.6).

The number densities $n_{ic}$ at $y=1$ are related to the prescribed bulk number densities $n_{i \infty }$, through the conservation constraint (2.11), which in this case becomes

(3.13)\begin{equation} q n_{i \infty}=\int_0^{H} {v}(x)n_i(x) \, {{\rm d}\kern0.7pt x}, \end{equation}

for $i=1,\ldots, m$, where ${v}(x)$ is the velocity (2.32), $q=\frac {2}{3} {v}_0 H$ is the corresponding volume flux per unit depth and $H$ is the half-width of the channel. Using dimensionless variables and (3.12) yields

(3.14)\begin{equation} \bar n_{i \infty}=\tfrac{12}{5}\epsilon^{1/2}x_{ic} (1-\tfrac{5}{8}\epsilon^{1/2}{\rm \Delta} \bar N_i )+O(\epsilon^{5/2}), \end{equation}

where $\bar n_{i \infty }=n_{i \infty }/n_{c}$ and

(3.15)\begin{equation} \epsilon=X_c/H, \end{equation}

is the ratio of the two length scales that affect the particle distributions. Under the assumption that the channel width is large compared with the particle size, $\epsilon \ll 1$. The quantity ${\rm \Delta} \bar N_i$ is the (average) deficit of type-$i$ particle density in the centre region compared with the singular distribution (3.12)

(3.16)\begin{equation} {\rm \Delta} \bar N_i =\int_0^{1} [y^{{-}1/2}- \bar N_i(y)] \, {{\rm d} y}, \end{equation}

where $\bar N_i$ is defined above. This particle deficit results from the modified transport coefficients in the region $\vert y\vert < 1$. The $O(\epsilon ^{5/2})$ error in (3.14) results from neglecting particle–wall interactions. Dividing (3.14) by its sum over all species provides a relation between the composition at $X_2=X_c$ and the bulk composition

(3.17)\begin{equation} x_{i \infty}=x_{i c}\left( \frac{ 1-\dfrac{5}{8}\epsilon^{1/2} {\rm \Delta} \bar N_i}{\displaystyle 1-\frac{5}{8}\epsilon^{1/2} \sum_{j=1}^m x_{j c}{\rm \Delta} \bar N_j}\right) +O(\epsilon^{2}), \end{equation}

for $i=1,\ldots, m-1$, where $x_{i \infty }=n_{i\infty }/n_\infty$ is the prescribed bulk number fraction of species $i$. This result in combination with (3.11) provides the boundary conditions for (3.4). In wide channels, $\epsilon \ll 1$, the composition in the centre region is insensitive to channel width; at leading order, (3.17) yields $x_{i c}\approx x_{i \infty }$ and boundary conditions for (3.4) simplify to

(3.18)\begin{equation} \bar n_i(1)=x_{i \infty}. \end{equation}

Dividing $\bar n_i(y)$ by (3.14) yields the particle distributions normalized by their bulk number densities,

(3.19)\begin{equation} \frac{n_i(y)}{n_{i \infty}}= f_i(\epsilon) \bar N_i(y), \end{equation}

where

(3.20)\begin{equation} f_i(\epsilon)=\frac{\frac{5}{12} \epsilon^{{-}1/2}}{1-\frac{5}{8} \epsilon^{1/2} {\rm \Delta} \bar N_i}. \end{equation}

Outside the centre region, this distribution reduces to

(3.21)\begin{equation} \frac{n_i(\tilde y)}{n_{i \infty}}= \frac{5}{12}\tilde y^{{-}1/2}\left(1-\frac{5}{8}\epsilon^{1/2}{\rm \Delta} \bar N_i \right)^{{-}1}+O(\epsilon^{2}), \quad \lvert y\rvert > 1, \end{equation}

according to (3.12), where $\tilde y=X_2/H$. The result indicates that the channel width sets the length scale of the distribution for $\vert y\vert > 1$.

3.2.2. Monodisperse particle distribution

For a single particle species, (3.9)–(3.11) reduce to the linear boundary value problem

(3.22)\begin{equation} {D}(y)\bar N'+ {V}(y)\bar N =0, \quad \bar N(1)=1, \end{equation}

where $\bar N$ is defined by (3.8b) and the index $i$ is dropped to distinguish the monodisperse case. The transport coefficients are given by

(3.23a,b)\begin{equation} {D}(y)={D}^{(1)}_{11}(y)+{D}^{(2)}_{11}(y), \quad {V}(y)=V_{11}(y). \end{equation}

The solution of boundary value problem (3.22) is

(3.24)\begin{equation} \bar N(y)=\exp\left( \int_{\vert y\vert }^{1} \frac{{V}(t)}{{D}(t)} \, {\rm d} t\right). \end{equation}

For $\lvert y\rvert > 1$, the transport coefficients (3.23a,b) reduce to their local forms

(3.25a,b)\begin{equation} {D}(y)=\lvert y\rvert (2 I^{A}_{112}+I^{B}_{112}), \quad {V}(y)=\mbox{sign}(y) (I^{A}_{112}+\tfrac{1}{2}I^{B}_{112}), \end{equation}

and thus ${V}(y)/{D}(y)=1/2y$, according to (2.51a–c), so that the outer solution (3.12) is recovered. Given that ${D}(y)$ is non-vanishing and ${V}(y)$ vanishes linearly for $\lvert y\rvert \to 0$, as discussed in § 2.3.3, it follows that $\bar N$ has a maximum with ${\rm d}\bar N/{{\rm d} y}=0$ at $y=0$, consistent with the scaling predicted by (2.55). From the continuity of ${D}(y)$ and ${V}(y)$ up to their second derivatives at $\lvert y\rvert =1$, it follows that $\bar N(y)$ is continuous up to its third derivative.

From (3.19), we have

(3.26)\begin{equation} \frac{n(y)}{n_{\infty}}=f(\epsilon) \bar N(y), \end{equation}

where $\bar N(y)$ is given by (3.24), ${\rm \Delta} \bar N$ given by

(3.27)\begin{equation} {\rm \Delta} \bar N=\int_0^1 \left[ y^{{-}1/2}-\exp\left( \int_{\vert y\vert }^{1} \frac{{V}(t)}{{D}(t)} \, {\rm d} t\right) \right]\, {{\rm d} y}, \end{equation}

and

(3.28)\begin{equation} f(\epsilon)=\frac{\frac{5}{12} \epsilon^{{-}1/2}}{1-\frac{5}{8} \epsilon^{1/2} {\rm \Delta}\bar N}. \end{equation}

3.2.3. Superposition approximation

The particle distribution in a polydisperse mixture with weak interactions between particles of different sizes can be approximated by a superposition of monodisperse distributions for each particle size (or size class).

Displacements resulting from collisions between particles of different sizes are usually smaller than displacements resulting from collisions between equal-size particles. This is true for particles that undergo contact interactions, as discussed in § 4.2, because collision rates diminish rapidly with size ratio (Adler Reference Adler1981; Wang, Zinchenko & Davis Reference Wang, Zinchenko and Davis1994; Reboucas & Loewenberg Reference Reboucas and Loewenberg2021b). Thus, the superposition approximation may be expected to hold for particles with vastly different sizes. It may be also expected to hold for similar-size particles because the effect is similar to increasing the total number density and the latter scales out of the equations, as explained in § 2.2.1. The superposition approximation is further supported by the fact that particle distributions are coupled only in the boundary layer; outside the boundary-layer particle distributions are decoupled, as shown in § 3.1.

According to the superposition approximation, distribution (3.19) reduces to

(3.29)\begin{equation} \frac{n_i(y_i)}{n_{i\infty}} \simeq f(\epsilon_i) \bar N(y_i), \end{equation}

for $i=1,\ldots, m$. Here,

(3.30a,b)\begin{equation} y_i=X_2/X_c^{ii}\quad \mbox{and}\quad \epsilon_i=X_c^{ii}/H, \end{equation}

where $\bar N(y_i)$, ${\rm \Delta} \bar N$ and $f(x)$ are given by (3.24) and (3.27) and (3.28); $X_c^{ii}$ is the boundary-layer thickness for a suspension of type-$i$ particles. The result reduces to the monodisperse distribution (3.26) by dropping the index $i$.

Comparing (3.19) and (3.29) indicates that the superposition approximation is equivalent to $\epsilon ^{-1/2}\bar N_i(y_i)\simeq \epsilon _i^{-1/2}\bar N(y_i)$ and $\epsilon ^{1/2}{\rm \Delta} \bar N_i\simeq \epsilon _i^{1/2}{\rm \Delta} \bar N$ and, thus,

(3.31a,b)\begin{equation} \bar N_i(y_i)\simeq \chi_i^{{-}1/2}\bar N(y_i), \quad {\rm \Delta}\bar N_i\simeq \chi_i^{1/2}{\rm \Delta} \bar N, \end{equation}

where $\chi _i$ is the ratio

(3.32)\begin{equation} \chi_i=X_c^{ii}/X_c, \end{equation}

and $\chi _i \leq 1$, according to (3.5).

Under the assumption that $X_c^{ii}$ scales with particle size, the above results predict that particle segregation occurs in the boundary layer region with relative enrichment of smaller particles. The dependence of the centreline particle density on particle size is described by the function $f(\epsilon _i)$. According to (3.28), particle densities at the centreline scale approximately with the inverse square root of particle size, assuming that particles are small compared with channel width, $\epsilon _i\ll 1$. Distributions of smaller particles adhere more closely to the singular outer distribution (3.12) as reflected by their smaller average density deficit (3.16), according to (3.31b).

The error of the superposition approximation, resulting from interactions between particles of different sizes in the boundary layer, is characterized by the polydisperse enrichment

(3.33)\begin{equation} {\rm \Delta} \bar N_{ij}= \int_0^1[\bar N_i(y)-\chi_i^{{-}1/2}\bar N(y/\chi_i )] \, {{\rm d} y}, \end{equation}

where $\bar N_i(y)$ is the exact distribution for type-$i$ particles in a polydisperse suspension and $\chi _i^{-1/2}\bar N(y_i)$ is superposition approximation (3.31a). According to (3.16) and (3.31b), we have the relation ${\rm \Delta} \bar N_i=\chi ^{1/2}{\rm \Delta} \bar N-{\rm \Delta} \bar N_{ij}$.

4.1. Contact interactions

Spherical particles that undergo short-range symmetry-breaking ‘contact’ interactions in shear flow are considered here, specifically, particles with surface roughness, permeable particles and emulsion drops. Pair trajectories of such particles are analytically integrated to yield formulas for binary cross-stream particle displacements ${\rm \Delta} X^{12}_k$, ${\rm \Delta} X^{21}_k$ ($k=2,3$) involving integrals of the standard pair mobility functions for spherical particles. Note that superscripts 1 and 2 are used in this section to refer to particle labels, not particle species, by contrast to the superscripts $i$ and $j$ introduced earlier. Accordingly, ${\rm \Delta} X^{12}_k$ refers to the net displacement of particle-$1$ resulting from its collision with particle-$2$ and vice versa for ${\rm \Delta} X^{21}_k$, whereas ${\rm \Delta} X^{ij}_k$ is the net displacement of any type-$i$ particle resulting from its collision with any type-$j$ particle. Similarly, we will use $r_c$ to denote the collision cross-section for particles 1 and 2, by contrast with the quantity $r_c^{ij}$, defined by (2.27a,b), that refers to the collision cross-section for any type-$i$ and $j$ particles. The particles have radii $a_1$ and $a_2$, and size ratio $\kappa =a_2/a_1$. We assign label 1 to the larger particle thus, only $\kappa \leq 1$ needs consideration.

Recall that the particles are assumed to be non-Brownian and neutrally buoyant. Here, creeping flow conditions are assumed, and particle inertia is neglected

4.1.2. Permeable particles

Particle permeability is another short-range, symmetry-breaking mechanism. The dimensionless permeability is defined $\bar K= k/\bar a^2$, where $k$ is the permeability and Darcy's law is used to describe the intraparticle flow. Weak permeability alleviates the lubrication pressure between particles, allowing particle contact, $h_0=0$, but otherwise has little effect on the pair interaction. Under the assumption that the particles are rigid and not cohesive, the effect of weak permeability closely resembles small-amplitude particle roughness. The following hydrodynamic equivalence between particle roughness and particle permeability was proposed based on the contact time between particles under the action of a constant force (Reboucas & Loewenberg Reference Reboucas and Loewenberg2021a):

(4.1)\begin{equation} \bar\delta \longleftrightarrow 0.7224 \nu^{1/5} \bar K^{2/5}, \end{equation}

where $\bar \delta$ and $\bar K$ are the dimensionless particle roughness and permeability, and $\nu =2\kappa / (1+\kappa )^2$. This relation is based on the assumption that no-slip boundary conditions apply on the particle surfaces. This correlation has been shown to hold accurately for colliding pair trajectories in several types of flow (including shear) for a wide range of size ratios and permeabilities (Reboucas & Loewenberg Reference Reboucas and Loewenberg2021b).

4.2. Trajectories of particles with contact interactions

Relative particle trajectories emanate from $x_1\to -\infty$ for $x_2^{-\infty }>0$ and from $x_1\to +\infty$ for $x_2^{-\infty }<0$. Apart from the contact interactions, trajectories are accurately described using standard pair mobility functions for spherical particles (or drops) in shear flow (Batchelor & Green Reference Batchelor and Green1972a; Zinchenko Reference Zinchenko1978, Reference Zinchenko1980, Reference Zinchenko1983). Accordingly, particles with short-range binary contact interactions have circular upstream collision cross-sections, defined by (2.27a,b), where $r_c$ depends on size ratio, and on the roughness, permeability, or drop- to continuous-phase viscosity ratio, respectively, for rough or permeable particles or drops. Collision cross-sections, or equivalently collision efficiencies $E_{12}=(r_c/(2\bar a))^3$, are available in the literature for rough and permeable particles and drops (Wang et al. Reference Wang, Zinchenko and Davis1994; Reboucas & Loewenberg Reference Reboucas and Loewenberg2021b). Trajectories with offsets $r^{-\infty } > r_c$ are reversible, i.e. ${\rm \Delta} X_k^{12}={\rm \Delta} X_k^{21}=0$.

Pair trajectories with upstream trajectory offsets, $r^{-\infty }< r_c$, reach the contact surface, defined by $s=s^*$, where $s=r/\bar a$ is the centre-to-centre separation, $r$, normalized by the average particle radius. For permeable particles and drops, $s^*=2$; for particles with surface roughness, $s=2+\bar \delta$. On the contact surface, the particles undergo relative rotation through the compressional quadrant of the flow and separate at the equator ($x_1=0$), under the assumption that cohesive forces are absent. The motion on the contact surface is described by a subset of the trajectory equations corresponding to zero relative radial velocity. Examples of contacting trajectories are shown in figures 4 and 5.

Figure 4. Relative (a) and pair (b) particle trajectories in velocity-gradient plane, size ratio $\kappa =a_2/a_1=1/2$, roughness $\bar \delta =d/\bar a=10^{-3}$; initial positions (i), contact surface (ii)–(iii) (dotted lines), final positions (iv). Relative $\boldsymbol {x}=\boldsymbol {X}^{(2)}-\boldsymbol {X}^{(1)}$ and pair $\bar {\boldsymbol {x}} =\boldsymbol {X}^{(1)}+\boldsymbol {X}^{(2)}$ coordinates of particles non-dimensionalized by the average radius $\bar a=\frac {1}{2}(a_1+a_2)$; initial conditions $\boldsymbol {x}=(-\infty,0.1,0)$ and $\bar {\boldsymbol {x}}=(-\infty,0.1,0)$.

Figure 5. Offset relative trajectory shown in cross-flow plane (a) and enlargement (b), size ratio $\kappa =a_2/a_1=1/2$; particles with roughness $\bar \delta = d/\bar a= 5 \times 10^{-3}$, frictionless contact (solid lines), infinite contact friction coefficient (dash-dotted lines); permeable particles $\bar K=6\times 10^{-6}$ from (4.1) (dotted lines); drops with viscosity ratio $\lambda =1$ (dashed lines); initial offset (i), contact surface (ii)–(iii), final offset (iv); collision surface $s^*$ (large circle), collision cross-sections $r_c$ particles (small circle), drops (dashed circle). Relative $\boldsymbol {x}=\boldsymbol {X}^{(2)}-\boldsymbol {X}^{(1)}$ coordinates of particles non-dimensionalized by the average radius $\bar a=\frac {1}{2}(a_1+a_2)$; initial conditions $\boldsymbol {x}=(-\infty,0.25,0.25)$.

Below a critical roughness, the collision cross-section for rough particles vanishes; similarly, there exists a critical permeability below which $r_c=0$. The values of these critical parameters increase with diminishing size ratio (Arp & Mason Reference Arp and Mason1977; Reboucas & Loewenberg Reference Reboucas and Loewenberg2021b). Likewise, drops have a critical viscosity ratio beyond which $r_c=0$, and the critical viscosity ratio decreases with diminishing size ratio (Wang et al. Reference Wang, Zinchenko and Davis1994). Equivalently, there exists a finite critical size ratio ratio, $\kappa ^*$, below which $r_c=0$ for fixed values of particle roughness or permeability, or drop viscosity ratio.

4.2.1. Maximum particle displacements

Maximum particle displacements are important because they determine the thickness of the boundary layer that forms in regions of vanishing shear rate. Maximum displacement magnitudes in the $X_2$-direction are achieved for $r^{-\infty }\to 0$ with $\theta =\pm {\rm \pi}/2$ (i.e. $x_3^{-\infty }=0$). For contact interactions between pairs of inertialess particles in creeping flows, the maximum displacement magnitudes satisfy

(4.2)\begin{equation} {\rm \Delta} X^{12}_{2, {max}}+ {\rm \Delta} X^{21}_{2, {max}} = r_c, \end{equation}

and for equal-size particles,

(4.3a,b)\begin{equation} {\rm \Delta} X^{12}_{2, {max}}={\rm \Delta} X^{21}_{2, {max}}=\tfrac{1}{2}r_c, \quad \kappa=1, \end{equation}

by symmetry.

Bounds for the magnitudes of individual displacements of unequal size particles are given by

(4.4a–c)\begin{equation} 0 < {\rm \Delta} X^{12}_{2, {max}} \leq \tfrac{1}{2}r_c, \quad \tfrac{1}{2}r_c \leq {\rm \Delta} X^{21}_{2, {max}} < r_c , \quad 0 < \kappa \leq 1, \end{equation}

and for extreme size ratios,

(4.5a–c)\begin{equation} {\rm \Delta} X^{12}_{2, {max}}\to 0, \quad {\rm \Delta} X^{21}_{2, {max}}\to r_c, \quad \kappa\to 0. \end{equation}

4.4. Boundary-layer thickness

Given that maximum displacements (4.2) occur for $r^{-\infty }\to 0$ and given ${\rm \Delta} X^{12}_{2}={\rm \Delta} X^{21}_{2}=0$ for $r^{-\infty }=r_c$, suggests that (2.36) can be replaced by the tighter bound

(4.6)\begin{equation} X_c^{ij} = max ({\rm \Delta} X^{ij}_{2, {max}},\tfrac{1}{2}r_c^{ij}), \end{equation}

for particles of types $i$ and $j$ that undergo contact interactions. Our calculations support this relation but our results do not rely on it. The relation is used here to discuss the width of the boundary layer, $X_c$, that forms where the shear rate vanishes.

In a monodisperse suspension, maximum particle displacements are set by the collision cross-section, according to (4.3a,b), and $r_c=r_c^{ii}$, thus the boundary-layer thickness (4.6) reduces to

(4.7)\begin{equation} X_c^{ii}=\tfrac{1}{2} r_c^{ii}. \end{equation}

The formula for the collision cross-section is provided in Appendix C. Below the critical values of particle roughness or permeability or above a critical value of the drop viscosity ratio, $r_c=0$, as seen in figures 8 and 9; particle structuring is not predicted in this case.

In a bidisperse suspension we have,

(4.8)\begin{equation} X_c=max ({\rm \Delta}_{2, {max}}^{21},\tfrac{1}{2} r_{c,{max}}), \end{equation}

given that ${\rm \Delta} X_2^{21}\geq {\rm \Delta} X_2^{12}$, according to (4.4a–c). Here, $r_{c,{max}}$ is the maximum collision cross-section among the permutations of pairs of type-$i$ and $j$ particles, i.e. $r_{c,{max}}=max (r_c^{ii}, r_c^{jj}, r_c^{ij})$. In binary mixtures of particles that differ only in size, the collision cross-section between the larger particles usually controls $r_{c,{max}}$. Consider the case where larger particles are type-$i$ and smaller, type $j$, particles differs only in size and $\kappa < 1$. In this case, we have $r_c^{ii}> r_c^{ij}$ because collision cross-sections increase monotonically with size ratio, and usually $r_c^{ii} > r_c^{jj}$, because collision cross-sections scale with particle size. Accordingly, we have $r_{c,{max}}=r_c^{ii}$. For drops, $r_c^{ii} =\kappa ^{-1} r_c^{jj}$ but this relation is inexact for rough and permeable particles because collision cross-sections increase monotonically with the size-dependent parameters, $d/a$ and $k/a^2$. Away from critical values of these parameters, however, the dependence is sub-linear, as seen in figure 8, thus $r_c^{ii} > r_c^{jj}$ holds and the relation $r_{c,{max}}=r_c^{ii}$ is retained.

The results shown in figures 10 and 11 illustrate relation (4.8), indicating that the boundary-layer thickness is controlled by heterogeneous pair interactions (i.e. ${\rm \Delta} _{2, {max}}^{21} > r_{c,{max}}$) for moderate size ratios, and by the collision cross-section of the larger particles closer to the critical size ratio, $\kappa ^*$. For $\kappa <\kappa ^*$, displacements due to interactions between particles of different sizes vanish identically; in this case, the superposition approximation, described in § 3.2.3, holds exactly.

For particles with contact interactions, pair displacements and the resulting boundary-layer thickness scale with the size of the suspended particles. The particle distributions described in § 3.2 thus depend on particle size.

5.1. Particle distribution in monodisperse suspensions

Results for monodisperse suspensions are presented here, including particle distributions, the average deficit of particle density in the boundary layer compared with the singular outer distribution and the underlying diffusive and drift velocity transport coefficients. Predicted particle distributions are compared with predictions of the diffusive flux model (Phillips et al. Reference Phillips, Armstrong, Brown, Graham and Abbott1992) and the experiments of Koh et al. (Reference Koh, Hookham and Leal1994).

The results in figure 12 show that the drift velocity ${V}$ but not the diffusion coefficient ${D}$ vanishes at the centre of the channel, where the shear rate vanishes, and illustrate the respective odd and even symmetries of these quantities and their behaviour for $X_2\to 0$, as discussed in § 2.3.3. These are the essential features that avoid the classical singularity at points of vanishing shear rate predicted by (2.37) and instead yield the smooth, non-singular behaviour predicted by (2.55). The dashed lines in figure 12 depict the local form of the transport coefficients that apply for $X_2>X_c$. The anticipated second-derivative discontinuity in the transport coefficients at $X_2=X_c$, discussed in § 2.3.3, is evident.

Figure 12. Dimensionless drift velocity (a) and diffusive (b) transport coefficients $\bar {{D}}=D X_c^{-6}$ and $\bar {{V}}=V X_c^{-5}$ (3.23a,b) for monodisperse suspensions of particles with roughness $\delta =d/a$ and drops with viscosity ratio $\lambda$, as indicated (solid lines); outer forms of transport coefficients (3.25a,b) (dotted lines). Boundary-layer thicknesses $X_c/a=0.349$ for rough particles, and $X_c/a=0.778$ for drops.

Examples of particle distributions for rough particles and emulsion drops are shown in figure 13; the two distributions are barely distinguishable. Similar results were obtained for other values of roughness and viscosity ratios. The average deficit of particle density in the boundary layer is insensitive to the particle parameters (roughness and viscosity ratio) and has an almost constant value, $0.711 \leq {\rm \Delta} \bar N\leq 0.714$, over a wide range of parameters, as seen in figure 14. The boundary-layer thickness (4.7) depends on details of the contact interactions between particles, as shown in figures 8 and 9; however, the particle distributions are almost independent of these details when displayed using the rescaled coordinate, $X_2/X_c$.

Figure 13. Distribution (3.24) for monodisperse suspension of particles with roughness $\delta =d/a$ and drops with viscosity ratio, $\lambda$, as indicated (solid lines); fit using (5.5) with ${\rm \Delta} \bar N=0.71$ (dashed line); outer solution (3.12) (dotted lines).

Figure 14. Average deficit of particle density in the boundary layer (3.27) for monodisperse suspensions of particles with roughness $\delta =d/a$ (solid line) and drops with viscosity ratio $\lambda$ (dashed line).

The above observation motivates a polynomial fit of the exact calculations. Enforcing continuity up to second derivatives at $y=\pm 1$, consistent with the expected continuity, discussed below (3.25a,b), and a prescribed value of the particle density deficit ${\rm \Delta} \bar N$ yields

(5.5)\begin{equation} \bar N\doteq \left\{\begin{array}{@{}ll} 1+t+\dfrac{5}{2}t^2+35(3-4{\rm \Delta}\bar N) t^3 & \lvert \tilde y\rvert < 1\\ \lvert \tilde y\rvert^{{-}1/2} & \lvert \tilde y\rvert > 1 \end{array}\right. , \end{equation}

where $t=\frac {1}{4}(1-\tilde y^2)$ and $\tilde y=X_2/H$. Inserting the value of ${\rm \Delta} \bar N$ from figure 14 yields an approximation for $\bar N$ with $0.2\,\%$ maximum pointwise relative error; using the nominal value, ${\rm \Delta} \bar N\doteq 0.71$, yields an approximation error of less than 1 % for most values of roughness or viscosity ratio. Given its accuracy, the polynomial fit is appropriate for use with the (less accurate) superposition approximation.

An outline of the diffusive flux model is provided in Appendix B. For the purpose of comparing the dilute theory with the diffusive flux model (Phillips et al. Reference Phillips, Armstrong, Brown, Graham and Abbott1992), we employ its frequently used approximate solution (Koh et al. Reference Koh, Hookham and Leal1994)

(5.6)\begin{equation} \phi=\frac{\phi_{m}}{1+\alpha (X_2/H)}, \end{equation}

where $\phi _{m}$ is the maximum packing fraction at the centreline of the channel, herein taken as $\phi _{m}=0.68$ following Koh et al. (Reference Koh, Hookham and Leal1994), and $\alpha$ is a numerical coefficient adjusted to enforce the bulk particle concentration (3.13). The prediction $\phi =\phi _{m}$ at the centreline is a consequence of using a suspension viscosity model as explained in Appendix B but is inappropriate for dilute suspensions. Moreover, the cusp-shape profile is unphysical and the corresponding absence of a particle scale is at odds with experiments (Koh et al. Reference Koh, Hookham and Leal1994).

The results in figures 15–17 show that the dilute theory is in approximate agreement with experiments for dilute suspensions and avoids the centreline cusp characteristic of the diffusive flux model. Empirical support for the predicted volume fraction and particle size dependence of the dilute theory is discussed in § 5.3.

Figure 15. Particle distribution in monodisperse suspension with roughness $\delta =d/a=10^{-2}$, bulk volume fraction $\phi _{\infty }=10\,\%$ and channel width $H/a=8.8$. Dilute theory (3.26) (solid line); classical diffusive flux model (5.6) (dashed line); data from Koh et al. (Reference Koh, Hookham and Leal1994), $L_{exp}/L_{ss}=5.0$ ($\bigcirc$).

Figure 16. Same as figure 15 except $H/a=15.6$; data (Koh et al. Reference Koh, Hookham and Leal1994) $L_{exp}/L_{ss}=1.6$ ($\bigcirc$).

Figure 17. Same as figure 15 except $\phi _{\infty }=20\,\%$; data (Koh et al. Reference Koh, Hookham and Leal1994) $L_{exp}/L_{ss}=10.$ ($\bigcirc$).

5.2. Particle distribution in bidisperse suspensions

Results for bidisperse suspensions are presented here, including particle distributions, obtained by exact calculation and by superposition, and a parametric exploration of polydisperse enrichment due to the interactions between particles of different sizes. Predicted particle distributions are compared with the experiments of Lyon & Leal (Reference Lyon and Leal1998b) in bidisperse suspensions.

The bidisperse particle distributions shown in figure 18 show that small particles are enriched relative to large particles at the centreline, consistent with the approximately inverse square-root size dependence of the centreline density, discussed below (3.32). When rescaled using superposition variables (3.30a) and (3.31a), the distributions for both particles approximately collapse onto a universal distribution given by (5.5), as seen in figure 19.

Figure 18. Particle distributions $\bar N_i$ ($i=1,2$) in bidisperse suspensions with bulk volume fractions, $\phi _{1\infty }=\phi _{2\infty }$; rough particles $\delta _1=d/a_1=10^{-3}$, size ratio $\kappa =a_2/a_1=0.6$ (a), $\kappa =0.8$ (b); drops with viscosity ratio $\lambda =1$, $\kappa =1/2$ (c); numerical solution of (3.9) and (3.18) for large (thick solid line) and small (thin solid line) particles; superposition approximation (3.31a) (dashed lines), outer solution (3.12) (dotted lines).

Figure 19. Same as figure 18 except rescaled using (3.30a) and (3.31a). Formula (3.24) (dashed lines).

Polydisperse enrichment (3.33) is a measure of the effect of interactions between particles of different sizes, an effect ignored by the superposition approximation. According to the discussion in § 3.2.3 and the remark at the end of § 4.4, the superposition approximation is expected to apply in bidisperse suspensions for $\kappa \lesssim \kappa ^*$ and for $\kappa \approx 1$, as confirmed by the results presented in figures 20 and 21 which show ${\rm \Delta} \bar N_{ij}\to 0$ in these limits. Particle roughness $\delta _1$ and viscosity ratio $\lambda$ determine the collision cross-section and thus the critical size ratio $\kappa ^*$ below which contact interactions are not predicted.

Figure 20. Polydisperse enrichment (3.33) for large (thick lines) and small (thin lines) particles, $\phi _{2\infty }/\phi _{1\infty }$ as indicated; (a) particles with roughness $\delta _1=d/a_1=10^{-3}$, (b) drops with viscosity ratio $\lambda =1$.

Figure 21. Polydisperse enrichment (3.33) for large (thick lines) and small (thin lines) particles, $\phi _{2\infty }/\phi _{1\infty }=1$; (a) particles, roughness $\delta _1=d/a_1$ as indicated; (b) drops, $\lambda$ as indicated.

The results in figures 20 and 21 demonstrate that the superposition approximation provides an estimate of the particle distribution for all size ratios (examples for intermediate size ratio are shown in figure 18). The weak coupling for all size ratios is likely due to the absence of coupling outside the boundary layer, as discussed in § 3.2.3. For rough particles, the superposition approximation is accurate to within a few per cent for a broad range of parameters; larger errors occur for emulsion drops because the displacements resulting from collisions between unequal drops exceed those for unequal rough particles, as seen by comparing figures 10(a) and 11(a).

Polydisperse enrichment in binary suspensions is most sensitive to size ratio. Polydisperse enrichment of larger particles is predicted for rough particles at all size ratios, as seen in panel (a) of figures 20 and 21, whereas smaller particles are enriched by the presence of much larger particles (small size ratios) and depleted by the presence of moderately larger particles. The particle distributions in panels (a) and (b) of figures 18 and 19 illustrate both regimes. By contrast, a polydisperse depletion of larger drops and enhancement of smaller drops is predicted for emulsions over most of the parameter range, as seen in panel (b) of figures 20 and 21. Panel (c) of figures 18 and 19 illustrates a typical drop distribution and shows the bigger effect of polydispersity.

The bulk composition, $\phi _{2\infty }/\phi _{1\infty }$, modulates the magnitude of polydisperse enrichment and depletion in suspensions of particles and drops but has less effect on whether it occurs, according to the results shown in figure 20. The viscosity ratio $\lambda$ also modulates the magnitude of polydisperse enrichment and depletion as seen in figure 21(b). Particle roughness has a weaker effect on the magnitude of polydisperse enrichment and depletion; figure 21(a) shows that the polydisperse enrichment of large particles is insensitive to roughness.

Figures 22 and 23 show particle distributions in bidisperse suspensions with 30 % volume fraction. The dilute theory quantitatively predicts the distribution of large particles in the experiments of Lyon & Leal (Reference Lyon and Leal1998b), except near the boundary, $\lvert X_2\rvert /H\geq 0.8$, where the data are reported to be unreliable (Lyon & Leal Reference Lyon and Leal1998b). Here, the superposition approximation of the dilute theory is used. Similarly, assuming independent transport of large particles, the authors fit their results for large particles using the suspension balance model for a monodisperse suspension (Morris & Brady Reference Morris and Brady1998).

Figure 22. Particle distribution in bidisperse suspension with size ratio, $\kappa =a_2/a_1=0.3$, roughness $\delta _1=d/a_1=10^{-2}$, bulk volume fractions $\phi _1=22.5\,\%$ and $\phi _2=7.5\,\%$, channel width $H/a_1=11$; dilute theory with superposition approximation (3.29) for large (thick line) and small (thin line) particles; data from Lyon & Leal (Reference Lyon and Leal1998b) large ($\triangle$) and small ($\bigcirc$) particles, $L_{exp}/L_{1_{ss}}=5.0$, $L_{exp}/L_{2_{ss}}=0.14$.

Figure 23. Same as figure 22, except for $\phi _1=\phi _2=15\,\%$; data (Lyon & Leal Reference Lyon and Leal1998b) large ($\triangle$) and small ($\bigcirc$) particles, $L_{exp}/L_{1_{ss}}=3.3$, $L_{exp}/L_{2_{ss}}=0.28$.

By contrast, the experiments depicted in figures 22 and 23 show essentially no flow-induced structuring of small particles; however, this is likely due to insufficient entry length to achieve the stationary distribution of small particles, as suggested by the authors (Lyon & Leal Reference Lyon and Leal1998b). As indicated in the figure captions, the entry length used in the experiments, $L_{exp}$, is at least 3.5 times shorter than the estimated length, $L_{2_{ss}}$, required to establish a stationary distribution of the smaller particles. The slow evolution of particle microstructure is well known (Nott & Brady Reference Nott and Brady1994; Phan-Thien & Fang Reference Phan-Thien and Fang1996; Lyon & Leal Reference Lyon and Leal1998a); however, the hypothesis that particles of different sizes have distinct entry lengths has received less consideration (Semwogerere & Weeks Reference Semwogerere and Weeks2008). The experiments reported by Lyon & Leal (Reference Lyon and Leal1998b), including data not shown here, are consistent with distinct entry lengths predicted by (5.3).

There have been few studies that measure polydisperse enrichment, i.e. deviations from superposition resulting from hydrodynamic interactions between particles of different size. Semwogerere & Weeks (Reference Semwogerere and Weeks2008) observed polydisperse enrichment of large or small particles at the centreline and suggested that particles with a shorter entry length become more enriched in the centre region. Their results are less relevant here because colloidal (Brownian) particles were used in their study.

The experiments shown in figures 22 and 23 have been compared with a bidisperse diffusive flux model (Shauly et al. Reference Shauly, Wachs and Nir1998; Chun et al. Reference Chun, Park, Jung and Won2019). In planar Poiseuille flow, the bidisperse diffusive flux model predicts cusp-shaped distributions of both particle species: an upward cusp for the larger particles with a prescribed maximum volume fraction at the centreline, and a downward cusp for smaller particles with zero concentration at the centreline. By adjusting the physical parameters (size ratio, volume fractions and channel width), it was shown that the distribution of large particles could be made to fit the data but the predictions for the small particles were found to be qualitatively at odds with observations (Lyon & Leal Reference Lyon and Leal1998b; Chun et al. Reference Chun, Park, Jung and Won2019).

5.3. Limitation of dilute theory

Here, we propose a volume fraction criterion for application of the dilute theory. Data aggregated from several experiments that satisfy this criterion and the required entry length (5.4) are used to further validate the dilute theory.

At larger volume fractions, the dilute theory fails by predicting a centreline volume fraction above the maximum packing limit. This observation provides a rough practical upper bound of the bulk particle volume fraction $\phi ^*_\infty$ in a monodisperse suspension and an upper bound $\phi ^*_{i\infty }$ for each particle species in a polydisperse suspension. Using the superposition approximation (3.29) with polynomial fit (5.5), and setting the centreline particle volume fraction equal to the maximum packing limit, $\phi (0)=\phi _{m}$, we obtain

(5.7)\begin{equation} \phi_{i\infty}^*/\phi_{m} \approx 1.61 \epsilon^{1/2}_i, \end{equation}

where $\epsilon _i$ is defined by (3.30b), and dropping the index $i$ for monodisperse suspensions. For particles with $1\,\%$ roughness, $X_c^{ii}\doteq 0.52 a_i$, yielding the upper bound $\phi _{i \infty }^*/\phi _{m} \approx 1.16\, (a_i/H)^{1/2}$. Taking $\phi _{m}=0.68$, we find $\phi _\infty ^* \approx 27\,\%$ for the conditions described in figures 15 and 17, and $\phi _\infty ^* \approx$ 20 % for those in figure 16. For the conditions in figures 22 and 23, we have $\phi _{1 \infty }^* \approx 24\,\%$ and $\phi _{2 \infty }^* \approx 13\,\%$ but a smaller value may be warranted to accommodate a smaller value for $\phi _{m}$.

All experiments from Koh et al. (Reference Koh, Hookham and Leal1994) and Lyon & Leal (Reference Lyon and Leal1998a) that satisfy constraint (5.7) and have the required entry length (5.4) are aggregated in figure 24, omitting data close to the wall, $\lvert X_2\rvert /H\geq 0.8$, as suggested by Lyon & Leal (Reference Lyon and Leal1998b). According to (3.29), these data should collapse onto a single curve when re-plotted using

(5.8)\begin{equation} \bar N_{exp}= \frac{\phi_i/\phi_{i \infty}}{f(\epsilon_i)}, \end{equation}

where $\epsilon _i$ is given by (3.30b) and $f(x)$ is defined by (3.28) with ${\rm \Delta} \bar N\doteq 0.71$. The index $i$ is dropped for monodisperse suspensions. As shown in figure 25, the data approximately collapse according to the dilute theory (and superposition approximation for bidisperse systems) when re-plotted this way.

Figure 24. Data for monodisperse suspensions from Koh et al. (Reference Koh, Hookham and Leal1994) $\phi =10\,\%,\ H/a=8.8, L_{exp}/L_{ss}=5.0$ ($\bigcirc$), $H/a=15.6, L_{exp}/L_{ss}=1.6$ ($\square$); $\phi =20\,\% , H/a=8.8, L_{exp}/L_{ss}=10.$ ($\lozenge$), $H/a=15.6, L_{exp}/L_{ss}=3.2$ ($\triangledown$). Data for large particles in bidisperse suspensions from Lyon & Leal (Reference Lyon and Leal1998b) $H/a_1=11, \phi _1=7.5\,\% , L_{exp}/L_{1_{ss}}=1.7$ ($\blacktriangle$), $\phi _1=10\,\% , L_{exp}/L_{1_{ss}}=2.2$ ($\blacksquare$), $\phi _1=15\,\% , L_{exp}/L_{1_{ss}}=3.3$ ($\blacktriangledown$), $\phi _1=20\,\% , L_{exp}/L_{1_{ss}}=4.4$ ($\mathbin {\blacklozenge }$), $\phi _1=22.5\,\% , L_{exp}/L_{1_{ss}}=5$ ($\bullet$).

Figure 25. Data from figure 24 rescaled using (5.8); theoretical curve $\bar N$ given by (3.24) (solid line).

The collapse of data for large particles in bidisperse suspensions (filled symbols) demonstrates the proportionality to bulk volume fraction because these data differ only in volume fraction. This is also seen in the collapse of data for monodisperse suspensions with 10 % and 20 % volume fractions with the same size particles (compare data points $\bigcirc$ with $\lozenge$ and $\square$ with $\triangledown$). Validation for the size dependence predicted by the dilute theory is seen by comparing data for suspensions with different sized particles at fixed volume fractions (compare $\bigcirc$ with $\square$ and $\lozenge$ with $\triangledown$).

Although the data are seen to collapse when re-plotted according to (5.8) in figure 25, the shape of the universal curve (solid line) is not accurately obtained. We note that the data for large particles from Lyon & Leal (Reference Lyon and Leal1998b) (filled symbols) appear to attain the shape of the universal curve more closely than the data from Koh et al. (Reference Koh, Hookham and Leal1994) (open symbols), consistent with the closer, quantitative agreement seen in figures 22 and 23 compared with that seen in figures 15–17. Systematic errors in the earlier study, discussed by Lyon & Leal (Reference Lyon and Leal1998a), might account for this difference.