2. Green's functions and the averaged Stokes equations

At low Reynolds number, the fluid component of a suspension obeys the Stokes equations:

(2.1)\begin{equation} \left.\begin{gathered} \eta_0\nabla^2{\boldsymbol{v}} - \boldsymbol{\nabla} P+ {\boldsymbol{b}}_f = 0,\\ \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}= 0, \end{gathered}\right\}\end{equation}

where ${\boldsymbol {v}}$ is the velocity, $P$ is the pressure and ${\boldsymbol {b}}_f$ is the body force on the fluid. In the presence of suspended spheres, these equations are only valid in the fluid region. Since the Stokes equations are linear, the effect due to the presence of the spheres can be accounted for by a Green's function contribution from each sphere. We, therefore, decompose the fluid velocity into a background flow ${\boldsymbol {v}}_0$ and a disturbance flow from the spheres ${\boldsymbol {u}}$, such that the velocity at point ${\boldsymbol {x}}$ is ${\boldsymbol {v}}({\boldsymbol {x}}) = {\boldsymbol {v}}_0({\boldsymbol {x}}) + {\boldsymbol {u}}({\boldsymbol {x}})$, with the total disturbance flow being a sum over the contribution from each individual sphere:

(2.2)\begin{equation} {\boldsymbol{u}}({\boldsymbol{x}}) = \sum_n {\boldsymbol{u}}_n({\boldsymbol{x}}). \end{equation}

Consider the sphere located at $\boldsymbol {x}_m$. This particle experiences an external flow given by the background flow and the disturbance flows from all the other particles, excluding itself, which we will denote as $\boldsymbol {u}_{-m}$. That is, the external flow about this particle is $\boldsymbol {v}_{-m} = \boldsymbol {v}_0 + \boldsymbol {u}_{-m}$. This flow field is not the same as the velocity field that would be present if the sphere at $\boldsymbol {x}_m$ were removed, as is described in more detail in § 4. In general, the disturbance flow due to a spherical object in an arbitrary background flow at low Reynolds number can be written as an expansion in singular solutions (i.e. Green's functions) centred at the sphere centre. For example, if the external flow is constant, the disturbance flow at location $\boldsymbol {x}$ due to the sphere centred at $\boldsymbol {x}_m$ is

(2.3)\begin{equation} \boldsymbol{u}_m(\boldsymbol{x}) = \frac{{\hat{\boldsymbol{F}}^0}}{8{\rm \pi}\eta_0}\boldsymbol{\cdot}\left(1 + \frac{a^2}{6}\nabla^2\right){\hat{\boldsymbol{S}}}\left(\boldsymbol{x},\boldsymbol{x}_m\right), \end{equation}

where ${\hat {\boldsymbol {F}}^0}$ is the net force that the particle exerts on the fluid and ${\hat {\boldsymbol {S}}}$ is the ‘Stokeslet’ or Oseen tensor:

(2.4)\begin{equation} {\hat{S}}_{ij} = \frac{\delta_{ij}}{\left|\boldsymbol{x} - \boldsymbol{x}_m\right|} + \frac{\left(x-x_m\right)_i\left(x - x_m\right)_j}{\left|\boldsymbol{x} - \boldsymbol{x}_m\right|^3}. \end{equation}

Similarly, for a sphere in pure linear flow, the disturbance flow is given by

(2.5)\begin{equation} \boldsymbol{u}_m(\boldsymbol{x}) = \frac{{\hat{\boldsymbol{F}}^1}}{8{\rm \pi}\eta_0}:\left(1 + \frac{a^2}{10}\nabla^2\right)\boldsymbol{\nabla}{\hat{\boldsymbol{S}}}\left(\boldsymbol{x},\boldsymbol{x}_m\right), \end{equation}

where ${\hat {\boldsymbol {F}}^1}$ is a rank-2 tensor with its antisymmetric part giving the net torque that the particle exerts on the fluid and the symmetric part giving the added fluid stress due to the presence of the particle. For a general arbitrary external flow, the disturbance velocity can be written as an expansion of singular solutions as

(2.6)\begin{equation} \boldsymbol{u}_m(\boldsymbol{x}) = \frac{1}{8{\rm \pi}\eta_0}\sum_{k=0}^\infty{\hat{\boldsymbol{F}}^k}\odot{\hat{\boldsymbol{S}}^k}. \end{equation}

Here ${\hat {\boldsymbol {F}}^k}$ are rank $k+1$ tensors that are determined from boundary conditions at the surface of the sphere. Tensors ${\hat {\boldsymbol {S}}^k}$ are rank $k+2$ tensors given by (Kim & Karilla Reference Kim and Karilla2005)

(2.7)\begin{equation} {\hat{\boldsymbol{S}}^k}\left(\boldsymbol{x},\boldsymbol{x}_m\right) = \left(1+\frac{a^2}{2\left(2k+3\right)}{\nabla}^2_{\boldsymbol{x}_m}\right)\boldsymbol{\nabla}^k_{\boldsymbol{x}_m} {\hat{\boldsymbol{S}}}\left(\boldsymbol{x},\boldsymbol{x}_m\right), \end{equation}

and we are using notation that $\odot$ denotes the full contraction of ${\hat {\boldsymbol {F}}^k}$ with the 2 to $k+2$ indices of ${\hat {\boldsymbol {S}}^k}$. For example, for $k=2$,

(2.8)\begin{equation} {\hat{\boldsymbol{F}}^2}\odot{\hat{\boldsymbol{S}}^2}={\hat{F}^2_{jkl}}\hat{S}^2_{ijkl}, \end{equation}

where we are using the Einstein summation convention that repeated indices are summed over. Given the magnitudes of the singularity strengths ${\hat {\boldsymbol {F}}^k}$ for each sphere, the general solution for the fluid velocity is

(2.9)\begin{equation} \boldsymbol{v}\left(\boldsymbol{x}\right) = \boldsymbol{v}_0\left(\boldsymbol{x}\right) + \frac{1}{8{\rm \pi}\eta_0}\sum_n\sum_k{\hat{\boldsymbol{F}}^k_n}\odot{\hat{\boldsymbol{S}}^k}\left(\boldsymbol{x},\boldsymbol{x}_n\right), \end{equation}

where the sum on $n$ is over the particle positions $\boldsymbol {x}_n$. Note that the velocity defined here is only valid in the fluid region. In the region surrounding the particle centred at $\boldsymbol {x}_m$, the velocity is given by $\boldsymbol {v}(\boldsymbol {x}) = \boldsymbol {v}_{p,m} + {\boldsymbol {\omega }}_{p,m}\times (\boldsymbol {x} - \boldsymbol {x}_m)$, where $\boldsymbol {v}_{p,m}$ and ${\boldsymbol {\omega }}_{p,m}$ are the translational and rotational velocities of the particle, respectively.

The corresponding Stokes equations for the fluid velocity are then given by

(2.10)\begin{equation} \eta_0\nabla^2\boldsymbol{v} - \boldsymbol{\nabla} P + \boldsymbol{b}_f + \sum_n\sum_k{\hat{\boldsymbol{F}}^k_n}\odot\boldsymbol{\nabla}_{\boldsymbol{x}_n}^k\left(1+\frac{a^2}{2\left(2k+3\right)}\nabla^2_{\boldsymbol{x}_n}\right) \delta^3\left(\boldsymbol{x}-\boldsymbol{x}_n\right)=0. \end{equation}

Note that the $\odot$ operator now denotes a full contraction of the singularity strength ${\hat {\boldsymbol {F}}^k}$ with the $k$th power of the gradient operator, which acts on the coordinate $\boldsymbol {x}_n$, and $\delta ^3 (\boldsymbol {x})$ is the three-dimensional Dirac delta function. Equation (2.10) is the effective dynamic equation for the fluid valid at all points in space and defines the fluid velocity in all regions outside the particles. To get to a continuum description of the suspension, we average (2.10) over a mesoscopic scale that is large compared with the size of the spheres, but small compared with the system size. We define an averaging kernel $K(|\boldsymbol {x}-\boldsymbol {y}|)$ that satisfies the conditions

(2.11)\begin{equation} \int K\left(\left|\boldsymbol{x}-\boldsymbol{y}\right|\right) {\rm d}^3 y = 1; \quad K\left(\left|\boldsymbol{x}-\boldsymbol{y}\right|\right) = 0 \quad {\rm as} \left|\boldsymbol{x}-\boldsymbol{y}\right|\rightarrow \infty, \end{equation}

and assume that $K(x)$ is a smooth function that is differentiable as often as required by our manipulations. We also define $R$ to be the effective radius of the averaging kernel. When computations do not involve derivatives of the averaging kernel, it is sometimes convenient to choose a step function where $K(x) = 3\ 4{\rm \pi} R^3$ for $x \leq R$ and $K(x)=0$ for $x > R$, which enables direct calculation of integrals over this kernel. Multiplying (2.10) by the averaging kernel and integrating over all space, we find

(2.12)\begin{equation} \eta_0\nabla^2\left\langle\boldsymbol{v}\right\rangle - \boldsymbol{\nabla}\left\langle P\right\rangle + \langle\boldsymbol{b}_f\rangle + \frac{3}{4{\rm \pi} a^3}\sum_k\left({-}1\right)^k\left( 1 + \frac{a^2}{2\left(2k+3\right)}\nabla^2\right)\boldsymbol{\nabla}^k\odot(\phi\langle{\hat{\boldsymbol{F}}^k}\rangle) = 0, \end{equation}

where

(2.13)\begin{equation} \left\langle\boldsymbol{v}\right\rangle\left(\boldsymbol{x}\right) = \int K\left(\left|\boldsymbol{x}-\boldsymbol{y}\right|\right) \boldsymbol{v}\left(\kern0.7pt\boldsymbol{y}\right) {\rm d}^3 y, \end{equation}

with a similar definition for $\left \langle P\right \rangle$, and where

(2.14)\begin{equation} \frac{3\phi}{4{\rm \pi} a^3}\langle {\hat{\boldsymbol{F}}^k}\rangle \left(\boldsymbol{x}\right) = \sum_n K\left(\left|\boldsymbol{x}-\boldsymbol{x}_n\right|\right){\hat{\boldsymbol{F}}^k_n}, \end{equation}

with the volume fraction defined by

(2.15)\begin{equation} \phi\left(\boldsymbol{x}\right) = \sum_n\int_{\left|\boldsymbol{x}\right|'\leq a} K(|\boldsymbol{x}-\boldsymbol{x}_n-\boldsymbol{x}'|){\rm d}^3x'. \end{equation}

Since the averaged velocity in (2.13) is over regions containing particles and fluid,

(2.16)\begin{equation} \left\langle \boldsymbol{v}\right\rangle = \left(1-\phi\right)\boldsymbol{v}_f + \phi\boldsymbol{v}_p, \end{equation}

where $\boldsymbol {v}_f$ is the fluid-phase averaged velocity and $\boldsymbol {v}_p$ is the particle-phase averaged velocity.

In order to define the effective stress tensor for the fluid, we can rewrite (2.12) as

(2.17)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\sigma^{e} + \langle \boldsymbol{b}_f\rangle + \frac{3\phi\langle {\hat{\boldsymbol{F}}^0}\rangle }{4{\rm \pi} a^3} = 0,\end{equation}

with the general form for the averaged two-phase fluid stress tensor $\sigma ^{e}$ given by

(2.18)\begin{align} \sigma^{e}_{ij} &= \eta_0\left(\frac{\partial\left\langle v\right\rangle _i}{\partial x_j} + \frac{\partial\left\langle v\right\rangle _j}{\partial x_i}\right) - P\delta_{ij} + \frac{1}{8{\rm \pi} a}\frac{\partial}{\partial x_j}(\phi\langle {\hat{F}^0}\rangle _i)\nonumber\\ &\quad +\frac{3}{4{\rm \pi} a^3}\sum_{k=1}^\infty\left({-}1\right)^k\left(1 + \frac{a^2}{2\left(2k+3\right)}\nabla^2\right)\boldsymbol{\nabla}^{k-1}\odot(\phi\langle {\hat{\boldsymbol{F}}^k}\rangle). \end{align}

This is the first major result of our paper. Given the singularity strengths for the spheres, this equation completely defines the effective stress tensor for the two-phase fluid. An interesting aspect of this stress tensor is that it is not necessarily symmetric. For example, the last term in the first line includes the $j$th component of the gradient and the $i$th component of the force ${\hat {\boldsymbol {F}}^0}$. In single-phase fluids, the stress tensor should be symmetric to prevent self-torquing. However, consider a suspension of sedimenting particles with a horizontal gradient in the volume fraction. In this case, any given volume element has more spheres on one side than it does on the other. Since the spheres are acted on by gravity, they exert a net force on the fluid, which is larger on one side of the volume element than on the other. Consequently, the spheres produce a net torque on the volume element, which is accounted for in this stress tensor. The possibility of an antisymmetric stress tensor was previously noted by Batchelor (Reference Batchelor1970) and later by Pedley & Kessler (Reference Pedley and Kessler1992) in the context of microorganism swimming. This same term also accounts for the primary difference between our results and the results from Jackson (Reference Jackson1997), which derived a similar term:

(2.19)\begin{equation} \frac{\partial}{\partial x_j}\left( \frac{3\eta_0\phi}{4}(v_{p,i} - v_{f,i})\right). \end{equation}

For non-interacting particles, Faxen's first law gives ${\hat {\boldsymbol {F}}^0} = 6{\rm \pi} \eta _0 a(\boldsymbol {v}_p - (1 + (a^2\ 6)\nabla ^2)\boldsymbol {v}_f)$ (Kim Reference Kim1985; Chen & Ye Reference Chen and Ye2000), in which case the last term in the first line of (2.18) is

(2.20)\begin{equation} \frac{\partial}{\partial x_j}\left( \frac{3\eta_0\phi}{4}\left(v_{p,i} - v_{f,i} - \frac{a^2}{6}\nabla^2 v_{f,i}\right)\right). \end{equation}

In the case of neutrally buoyant particles for which ${\hat {\boldsymbol {F}}^0} = 0$, our term does not contribute to the stress tensor, whereas Jackson's term becomes

(2.21)\begin{equation} \frac{\partial}{\partial x_j}\left( \frac{3a^2\phi}{24}\nabla^2 v_{f,i} \right), \end{equation}

which gives a force per volume that is positively proportional to the bi-Laplacian of the velocity and causes short-wavelength perturbations, on lengths smaller than the size of the individual particles, to be unstable in Jackson's formulation. As we show, our formulation leads to a term in the force per volume that is negatively proportional to the bi-Laplacian of the velocity, thereby providing a stabilizing effect to the dynamics.

In what follows, we consider a truncation of the sum on $k$ such that it only includes terms up to $k=2$. In that case, the stress tensor simplifies to

(2.22) \begin{align} \sigma^{e}_{ij} &= \eta_0\left(\frac{\partial\left\langle v\right\rangle _i}{\partial x_j} + \frac{\partial\left\langle v\right\rangle _j}{\partial x_i}\right) - P\delta_{ij} + \frac{1}{8{\rm \pi} a}\frac{\partial}{\partial x_j}(\phi\langle {\hat{F}^0}\rangle _i)\nonumber\\ &\quad -\frac{3}{4{\rm \pi} a^3}\left(\phi\langle {\hat{F}^1}\rangle _{ij} + \frac{a^2}{10}\nabla^2(\phi\langle {\hat{F}^1}\rangle _{ij}) - \frac{\partial}{\partial x_k}(\phi\langle {\hat{F}^2}\rangle _{ijk})\right). \end{align}

We show that the leading-order contribution to each singularity strength ${\hat {\boldsymbol {F}}^k}$ is proportional to $a^{2k+1}$ times the $k$th derivative of the average velocity field. Since the contribution from each ${\hat {\boldsymbol {F}}^k}$ comes into the stress as $\boldsymbol {\nabla }^{k-1}\odot {\hat {\boldsymbol {F}}^k}$, truncating at $k=2$ omits terms that are proportional to the sixth derivative of the velocity field and are of order $(a/L)^4$ and smaller in the force per volume, where $L$ is a natural length scale for the system. To be consistent, we have, therefore, left off the term in (2.22) that depends on the Laplacian of ${\hat {\boldsymbol {F}}^2}$ as it is of order $(a/L)^4$.

3. Determining the singularity strengths

In order to complete the closure of the stress tensor, we need to relate the singularity strengths (i.e. ${\hat {\boldsymbol {F}}^k}$) to the flow profile. This is handled by enforcing the no-slip boundary condition on the surface of the spheres. We begin by considering the Taylor series expansion of the external flow, $\boldsymbol {v}_{-m}$, around a sphere with centre located at $\boldsymbol {x}_m$ that is translating with velocity $\boldsymbol {v}_p$ and rotating with angular speed ${\boldsymbol {\omega }}_p$:

(3.1)\begin{equation} v_{{-}m,i} = v_{{-}m,i}\left(\boldsymbol{x}_m\right) + \left.{\hat{x}}_j\frac{\partial v_{{-}m,i}}{\partial x_j}\right|_{\boldsymbol{x}_m} + \left.\frac{1}{2}{\hat{x}}_j{\hat{x}}_k\frac{\partial^2 v_{{-}m,i}}{\partial x_j\partial x_k}\right|_{\boldsymbol{x}_m} + \cdots , \end{equation}

where ${\hat {\boldsymbol {x}}}$ is the vector displacement measured from the centre of the sphere. Comparing the velocity defined by the Green's functions evaluated at the surface of the sphere with the terms in the Taylor series expansion, we find the following relationships:

(3.2)$$\begin{gather} {\hat{F}^0_i} - \frac{3}{14 a^2}(5{\hat{F}^2_{ijj}} - 2{\hat{F}^2_{jij}} - 2{\hat{F}^2_{jji}}) = 6{\rm \pi}\eta_0 a(v_{p,i} - v_{{-}m,i}), \end{gather}$$

(3.3)$$\begin{gather}4{\hat{F}^1_{ij}} - {\hat{F}^1_{ji}} ={-}20{\rm \pi}\eta_0 a^3\left(\frac{\partial v_{{-}m,i}}{\partial x_j} - \epsilon_{ijk}\omega_{p,k}\right), \end{gather}$$

(3.4)$$\begin{gather}6{\hat{F}^2_{ijk}} - {\hat{F}^2_{jik}} - {\hat{F}^2_{kji}} - {\hat{F}^2_{llk}}\delta_{ij}-({\hat{F}^2_{ljl}} + {\hat{F}^2_{jll}})\delta_{ik} ={-}\frac{14{\rm \pi}\eta_0 a^5}{3}\frac{\partial^2 v_{{-}m,i}}{\partial x_j\partial x_k}. \end{gather}$$

As mentioned previously, ${\hat {\boldsymbol {F}}^k}$ are proportional to the $k$th derivative of the velocity. They also depend on the ${\hat {\boldsymbol {F}}^{k+2}}$ terms. Indeed, it is a general result that the $k$th term in the Taylor expansion only depends on ${\hat {\boldsymbol {F}}^k}$ and ${\hat {\boldsymbol {F}}^{k+2}}$, since the functional form of the lower-order Green's functions only includes powers of $\boldsymbol {x}^k$ when evaluated at the surface of the sphere, and because the bi-Laplacian of the velocity field is zero, terms of higher order than $k+2$ also do not contribute. To be consistent with our previous truncation in the stress, we also truncate the Taylor series at $k=2$. Because the truncation ignores the higher-order derivatives in the velocity field, we are effectively not considering short-range lubrication forces or contact forces between the spheres. We have also used that $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {v}_{-m}=0$ (which also implies that ${\hat {F}^1_{ii}} = 0$). Contracting (3.4) on each pair of its indices, we find

(3.5a–c)\begin{align} {\hat{F}^2_{ijj}} = \frac{2{\rm \pi}\eta_0 a^5}{3}\nabla^2 v_{{-}m,i}, \quad {\hat{F}^2_{jij}} = \frac{20{\rm \pi}\eta_0 a^5}{9}\nabla^2 v_{{-}m,i}, \quad {\hat{F}^2_{jji}} = \frac{16{\rm \pi}\eta_0 a^5}{9}\nabla^2 v_{{-}m,i}, \end{align}

which we can use to simplify (3.2) as

(3.6)\begin{equation} {\hat{F}^0_i} = 6{\rm \pi}\eta_0 a\left( v_{p,i} - v_{{-}m,i} - \frac{a^2}{6}\nabla^2 v_{{-}m,i} \right).\end{equation}

Equations (3.3) and (3.6) are equivalent to Faxen's laws, with the imposed fluid velocity replaced by $\boldsymbol {v}_{-m}$ (Kim Reference Kim1985; Chen & Ye Reference Chen and Ye2000). Note, however, that (3.3) does not contain the Laplacian term that is often included in Faxen's third law (which sets the symmetric part of ${\hat {\boldsymbol {F}}^1}$). This term comes in at order $k=3$, which would be inconsistent to include here, without adding the $k=3$ term. It is interesting to note that the fluid stress tensor (2.18) does include the Laplacian term acting on ${\hat {\boldsymbol {F}}^1}$. The solution to the third equation, (3.4), is found to be

(3.7)\begin{align} {\hat{F}^2_{ijk}} &={-}\frac{{\rm \pi}\eta_0 a^5}{6}\left( 5\frac{\partial^2 v_{{-}m,i}}{\partial x_j\partial x_k} + \frac{\partial^2 v_{{-}m,j}}{\partial x_i\partial x_k} + \frac{\partial^2 v_{{-}m,k}}{\partial x_i\partial x_j}\right)\nonumber\\ &\quad +\frac{{\rm \pi}\eta_0 a^5}{18}(3\nabla^2 v_{{-}m,i}\delta_{jk} + 11\nabla^2 v_{{-}m,j}\delta_{ik} + 7\nabla^2 v_{{-}m,k}\delta_{ij}). \end{align}

4. Relating the local flow about a sphere to the fluid-phase averaged velocity

Equations (3.2)–(3.4) determine the singularity strengths in terms of the external flow about a sphere, $\boldsymbol {v}_{-m}$. This velocity field represents the local background flow that a sphere located at $\boldsymbol {x}_m$ experiences. That is, it is the total flow minus the disturbance flow from this sphere:

(4.1)\begin{equation} \boldsymbol{v}_{{-}m}\left(\boldsymbol{x}\right) = \boldsymbol{v}\left(\boldsymbol{x}\right) - \frac{1}{8{\rm \pi}\eta_0}\sum_{k=0}^\infty{\hat{\boldsymbol{F}}^k_m}\odot{\hat{\boldsymbol{S}}^k}\left(\boldsymbol{x},\boldsymbol{x}_m\right). \end{equation}

As mentioned previously, this is not the same as the velocity field that would be present if the sphere at $\boldsymbol {x}_m$ were removed, since removing the sphere would alter the singularity strengths on the other spheres. The specific values of the singularity strengths are important, as they are defined in such a way as to provide the correct boundary conditions on the surfaces of the spheres, up to the order of our truncation. As defined, then, $\boldsymbol {v}_{-m}$ is not measurable. To develop a physically meaningful continuum description, we need to relate $\boldsymbol {v}_{-m}$ to an experimentally measurable quantity, such as the fluid-phase averaged velocity, $\boldsymbol {v}_f$. Since $\boldsymbol {v}_f$ is an averaged velocity, we consider first the average of $\boldsymbol {v}_{-m}$ over the same averaging region as is defined by $K(x)$, using a similar point averaging as was done for ${\hat {\boldsymbol {F}}^k}$ in (2.14):

(4.2) \begin{align} \left\langle \boldsymbol{v}_{{-}m}\right\rangle \left(\boldsymbol{x}\right) = \frac{4{\rm \pi} a^3}{3\phi}\sum_m \left(\vphantom{\rule{1pt}{19pt}}\boldsymbol{v}_0\left(\boldsymbol{x}_m\right) + \frac{1}{8{\rm \pi}\eta_0}\smash{\sum_k\sum_{n\ne m}}{\hat{\boldsymbol{F}}^k_n}\odot{\hat{\boldsymbol{S}}^k} \left(\boldsymbol{x}_m,\boldsymbol{x}_n\right)\right)K\left(\left|\boldsymbol{x}-\boldsymbol{x}_m\right|\right), \end{align}

where the sum on $n\ne m$ is equivalent to subtracting the contribution from the sphere at $\boldsymbol {x}_m$.

The fluid-phase averaged velocity is defined by integrating only over the fluid volume $V_f$ contained within the averaging region. Following the definition given in Jackson (Reference Jackson1997),

(4.3)\begin{equation} \boldsymbol{v}_f\left(\boldsymbol{x}\right) = \frac{1}{\left(1-\phi\right)} \int_{V_f}\left( \boldsymbol{v}_0\left(\kern0.7pt\boldsymbol{y}\right) + \frac{1}{8{\rm \pi}\eta_0}\sum_k\sum_n{\hat{\boldsymbol{F}}^k_n}\odot{\hat{\boldsymbol{S}}^k}\left(\kern0.7pt\boldsymbol{y},\boldsymbol{x}_n\right)\right)K\left(\left|\boldsymbol{x}-\boldsymbol{y}\right|\right){\rm d}^3y. \end{equation}

Using these two definitions and taking the continuum limit (as shown in Appendix A), we find

(4.4)\begin{align} \left\langle \boldsymbol{v}_{{-}m}\right\rangle \left(\boldsymbol{x}\right) &= \boldsymbol{v}_f\left(\boldsymbol{x}\right) + \frac{3\phi}{32{\rm \pi}^2\eta_0 a^3\left(1-\phi\right)}\sum_k \left[ {\hat{\boldsymbol{F}}^k}\left(\boldsymbol{x}\right)\odot \int_{r = a}^{r = 2R} \left( g\left(\boldsymbol{r}\right) - 1 \right){\hat{\boldsymbol{S}}^k}\left(\boldsymbol{r}\right){\rm d}^3r\right.\nonumber\\ &\quad + \left. \frac{a^2\phi}{10} {\hat{\boldsymbol{F}}^k}\left(\boldsymbol{x}\right)\odot \int_{r = a}^{r = 2R} g\left(\boldsymbol{r}\right)\nabla^2{\hat{\boldsymbol{S}}^k}\left(\boldsymbol{r}\right){\rm d}^3r\right], \end{align}

where $r = |\boldsymbol {r}|$, $g(\boldsymbol {r})$ is the pair-distribution function, which gives the probability of finding a particle at location $\boldsymbol {r}$ with respect to the centre of another particle (Kirkwood Reference Kirkwood1935; Guazzelli & Morris Reference Guazzelli and Morris2011), and we have dropped the brackets on the averages of ${\hat {\boldsymbol {F}}^k}$ for simplicity.

For isotropic $g=g(r)$, only even values of $k$ give non-zero integrals. In addition,

(4.5)\begin{equation} \int_{r=a}^{r=2R} g\left(r\right)\nabla^2{\hat{\boldsymbol{S}}^k} \,{\rm d}^3r = 0, \end{equation}

since the flow rate through any closed surface enclosing one of the singularities of the Green's functions is zero (Pozrikidis Reference Pozrikidis1997), and it can be shown that only the $k=0$ term contributes to the first integral. Therefore, when the pair-distribution function is isotropic,

(4.6)\begin{equation} \left\langle \boldsymbol{v}_{{-}m}\right\rangle \left(\boldsymbol{x}\right) = \boldsymbol{v}_f\left(\boldsymbol{x}\right) + \frac{3\phi}{32{\rm \pi}^2\eta_0 a^3\left(1-\phi\right)} {\hat{\boldsymbol{F}}^0}\left(\boldsymbol{x}\right)\odot {\hat{\boldsymbol{Q}}^0},\end{equation}

where we have defined

(4.7) \begin{equation} {\hat{\boldsymbol{Q}}^0}= \int_{r = a}^{r = 2R} \left( g\left(\boldsymbol{r}\right) - 1 \right){\hat{\boldsymbol{S}}^0}\left(\boldsymbol{r}\right){\rm d}^3r. \end{equation}

In general, for non-isotropic $g$, the other integrals can contribute to the relationship between $\left \langle \boldsymbol {v}_{-m}\right \rangle$ and $\boldsymbol {v}_f$.

Using similar calculations, it is straightforward to compute $\left \langle \boldsymbol {\nabla }\boldsymbol {v}_{-m}\right \rangle$. We find that

(4.8)\begin{align} \left\langle \boldsymbol{\nabla}\boldsymbol{v}_{{-}m}\right\rangle \left(\boldsymbol{x}\right) = \left\langle \boldsymbol{\nabla}\boldsymbol{v}\right\rangle _f\left(\boldsymbol{x}\right) + \frac{3\phi}{32{\rm \pi}^2\eta_0 a^3\left(1-\phi\right)} {\hat{\boldsymbol{F}}^1}\left(\boldsymbol{x}\right)\odot \int_{r = a}^{r = 2R} \left( g\left(\boldsymbol{r}\right) - 1 \right)\boldsymbol{\nabla}{\hat{\boldsymbol{S}}^1}\left(\boldsymbol{r}\right){\rm d}^3r, \end{align}

where $\left \langle \boldsymbol {\nabla }\boldsymbol {v}\right \rangle _f$ is the fluid-phase average of the gradient of $\boldsymbol {v}$, which is not the same as $\boldsymbol {\nabla }\boldsymbol {v}_f$. Indeed, from the definitions of the averaging procedure,

(4.9)\begin{equation} \left(1-\phi\right)\left\langle \boldsymbol{\nabla}\boldsymbol{v}\right\rangle _f + \phi\left\langle \boldsymbol{\nabla}\boldsymbol{v}\right\rangle _p = \boldsymbol{\nabla} [\left(1-\phi\right)\boldsymbol{v}_f + \phi\boldsymbol{v}_p], \end{equation}

with $\left \langle \right \rangle _p$ denoting the particle-phase average (see Jackson Reference Jackson1997). For isotropic $g(\boldsymbol {r})$, the integral in (4.8) is zero, and we find

(4.10)\begin{equation} \left\langle \boldsymbol{\nabla}\boldsymbol{v}_{{-}m}\right\rangle = \frac{1}{\left(1-\phi\right)}\boldsymbol{\nabla}[\left(1-\phi\right)\boldsymbol{v}_f + \phi\boldsymbol{v}_p] - \frac{\phi}{\left(1-\phi\right)}\left\langle \boldsymbol{\nabla}\boldsymbol{v}\right\rangle _p. \end{equation}

For cases where the pair-distribution function is not isotropic, such as in the presence of shear, the gradient of $\boldsymbol {v}_{-m}$ can also depend on the values of ${\hat {\boldsymbol {F}}}$. We discuss this case later.

Using the same analysis, we also find

(4.11)\begin{equation} \left\langle \boldsymbol{\nabla}\boldsymbol{\nabla}\boldsymbol{v}_{{-}m}\right\rangle = \left\langle \boldsymbol{\nabla}\boldsymbol{\nabla}\boldsymbol{v}\right\rangle _f.\end{equation}

Taking into account the discontinuity in the gradient of $\boldsymbol {v}$ at the surface of the spheres,

(4.12)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\nabla}\left\langle \boldsymbol{v}\right\rangle = \left(1-\phi\right)\left\langle \boldsymbol{\nabla}\boldsymbol{\nabla}\boldsymbol{v}\right\rangle _f + \sum_m\frac{1}{V}\int_{S_p}\left(\boldsymbol{\nabla}\boldsymbol{v}\right)\boldsymbol{n}\,{\rm d}S_p,\end{equation}

where the integral is over the surfaces $S_p$ of the sphere included in the averaging region and $\boldsymbol {n}$ is the outward normal to the surface. Combining the results from (4.11) and (4.12),

(4.13) \begin{equation} \boldsymbol{\nabla}_k\boldsymbol{\nabla}_j\left\langle v_i\right\rangle = \langle \boldsymbol{\nabla}_k\boldsymbol{\nabla}_j v_{{-}m,i}\rangle -\frac{\phi}{20{\rm \pi}\eta_0 a^3}(4{\hat{F}_i^0}\delta_{jk} - {\hat{F}_j^0}\delta_{ik} - {\hat{F}_k^0}\delta_{ij}).\end{equation}

The Laplacian of $\boldsymbol {v}_{-m}$ is then

(4.14)\begin{equation} \langle \nabla^2\boldsymbol{v}_{{-}m}\rangle = \nabla^2[\left(1-\phi\right)\boldsymbol{v}_f + \phi\boldsymbol{v}_p] + \frac{\phi}{2{\rm \pi}\eta_0 a^3}{\hat{\boldsymbol{F}}^0}. \end{equation}

Note that both the gradient and the double gradient of $\boldsymbol {v}_{-m}$ only depend on the full average velocity and therefore

(4.15)\begin{equation} \left.\begin{gathered} \left\langle \boldsymbol{\nabla}\boldsymbol{v}_{{-}m}\right\rangle = \frac{1}{\left(1-\phi\right)}\boldsymbol{\nabla}{\left\langle\boldsymbol{v}\right\rangle} - \frac{\phi}{\left(1-\phi\right)} \left\langle \boldsymbol{\nabla}\boldsymbol{v}\right\rangle _p,\\ \langle \nabla^2\boldsymbol{v}_{{-}m}\rangle = \nabla^2{\left\langle\boldsymbol{v}\right\rangle} + \frac{\phi}{2{\rm \pi}\eta_0 a^3}{\hat{\boldsymbol{F}}^0}. \end{gathered}\right\} \end{equation}

The second term on the right-hand side of the Laplacian equation above is equivalent to the $\phi \ 2$ correction that was found by Batchelor (Reference Batchelor1972).

5. The dynamic equations for low-Reynolds-number suspensions

The results from the previous sections can now be used to write out the full dynamic equations for the two-phase suspension. We begin with the averaged equations of motion for the spheres. These equations are dictated by the condition that the sums of the forces and torques on the spheres are equal to zero. These two conditions are directly related to Faxen's first and second laws, (3.6) and (3.3), and set the values for ${\hat {\boldsymbol {F}}^0}$ and the antisymmetric part of ${\hat {\boldsymbol {F}}^1}$. In the absence of body forces on the spheres, the net force that the sphere exerts on the fluid must be zero, which gives that ${\hat {\boldsymbol {F}}^0} = 0$. When there are body forces on the fluid and the spheres, such as when the spheres are not neutrally buoyant, then ${\hat {\boldsymbol {F}}^0}$ needs to be related to the body forces. We begin by integrating (2.10) over the volume of a single sphere, from which we find

(5.1)\begin{equation} -\int_{S_p} {\boldsymbol{\sigma}}\boldsymbol{\cdot}\boldsymbol{n} \,{\rm d}S_p + \frac{4{\rm \pi} a^3}{3}\boldsymbol{b}_f+ {\hat{\boldsymbol{F}}^0} = 0, \end{equation}

where $\boldsymbol {n}$ is the outward normal to the sphere and ${\boldsymbol {\sigma }}$ is the pure fluid stress tensor (not the effective stress tensor). The integral of the stress tensor dotted with the normal is the force that the fluid exerts on the sphere. In the presence of body forces $\boldsymbol {b}_p$ on the particles, this force must be equal to $-4{\rm \pi} a^3\boldsymbol {b}_p\ 3$, from which we get that

(5.2)\begin{equation} {\hat{\boldsymbol{F}}^0} = \frac{4{\rm \pi} a^3}{3}(\boldsymbol{b}_p - \boldsymbol{b}_f).\end{equation}

When the body forces are from gravity, ${\hat {\boldsymbol {F}}^0}$ is just equal to the difference between the gravitational force on the sphere and the buoyant force from the fluid. If the net torque on the particles is zero, then the antisymmetric piece of ${\hat {\boldsymbol {F}}^1}$ must be zero. From (3.3), (3.6) and the relations between $\boldsymbol {v}_{-m}$ and ${\left \langle \boldsymbol {v}\right \rangle }$, the equations of motion for the spheres are then

(5.3)$$\begin{gather} \boldsymbol{v}_p - {\left\langle\boldsymbol{v}\right\rangle} - \frac{a^2\left(1-\phi\right)}{6}\nabla^2{\left\langle\boldsymbol{v}\right\rangle} = \frac{a^2}{9\eta_0}\left((2-\phi-\phi^2){\boldsymbol{\mathcal{I}}}+ \frac{9\phi}{8{\rm \pi} a^2}{\hat{\boldsymbol{Q}}^0}\right)\boldsymbol{\cdot}(\boldsymbol{b}_p - \boldsymbol{b}_f), \end{gather}$$

(5.4)$$\begin{gather}\left(1 + \phi\right)\epsilon_{ikj}\omega_{p,k} = \frac{\partial\left\langle v\right\rangle _i}{\partial x_j} - \frac{\partial\left\langle v\right\rangle _j}{\partial x_i}, \end{gather}$$

with ${\boldsymbol {\mathcal {I}}}$ the identity tensor. One difference between this result and the results from Jackson (Reference Jackson1997) is the presence of the overall average velocity in the Faxen force in the first equation. In addition, this equation includes the dependence of the particle velocity on the volume fraction of the particles. Indeed, if the full pair-distribution function is known, then (5.3) should give the correct relationship between the particle velocities, the fluid velocity and the body forces for any volume fraction. Another interesting aspect is the presence of the factor of $1+\phi$ between the magnitudes of the angular velocity of the particles and the vorticity of the average fluid velocity. It is important to note that the spheres still rotate with the local vorticity; the factor of $1+\phi$ arises because the average fluid vorticity is not the same as the local fluid vorticity.

The force balance equation for the fluid is given by (2.17) with ${\hat {\boldsymbol {F}}^0}$ as in (5.2) and the effective stress tensor

(5.5) \begin{align} \sigma^{e}_{ij} &= \eta_0\left(\frac{\partial\left\langle v\right\rangle _i}{\partial x_j} + \frac{\partial\left\langle v\right\rangle _j}{\partial x_i}\right) - P\delta_{ij} + \frac{a^2}{6}\frac{\partial}{\partial x_j}(\phi(b_{p,i}-b_{f,i}))\nonumber\\ &\quad +\frac{5\eta_0}{2\left(1-\phi\right)}\left(1+\frac{a^2}{10}\nabla^2\right)\left(\phi\left(\frac{\partial\left\langle v\right\rangle _i}{\partial x_j} + \frac{\partial\left\langle v\right\rangle _j}{\partial x_i}\right) \right) - \frac{\partial}{\partial x_k}\left(\frac{3\phi}{4{\rm \pi} a^3}\langle {\hat{F}^2}\rangle _{ijk}\right), \end{align}

where ${\hat {\boldsymbol {F}}^2}$ is given by (3.7) with the double gradient of $\boldsymbol {v}_{-m}$ given in (4.13). This stress tensor has the nice property that it only depends on the overall average velocity, and does not depend independently on either $\boldsymbol {v}_f$ or $\boldsymbol {v}_p$. The first term in the second line reflects not only the Einstein correction to the viscosity, but also an additional correction that goes as $1\ (1-\phi )$. We also find terms of third order in gradients of the velocity. For neutrally buoyant particles at constant volume fraction, the fluid equation of motion simplifies to

(5.6)\begin{equation} \eta_{eff}\nabla^2{\left\langle\boldsymbol{v}\right\rangle} - \frac{a^2\eta_0\phi}{4}\nabla^4{\left\langle\boldsymbol{v}\right\rangle} - \boldsymbol{\nabla} P = 0, \end{equation}

where the effective viscosity is

(5.7)\begin{equation} \eta_{eff} = \eta_0\left(1 + \frac{5\phi}{2\left(1-\phi\right)}\right).\end{equation}

As mentioned, the bi-Laplacian term in the dynamic equation is negative, which provides stability to the equations.

The force and torque balance on the spheres, (5.3) and (5.4), along with the effective stress tensor (5.5) and the averaged Stokes equation (2.17) are the dynamic equations for the flow of the two-phase suspension. These equations of motion are closed by the conservation of volume condition:

(5.8)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}{\left\langle\boldsymbol{v}\right\rangle} = 0, \end{equation}

which defines the pressure. For cases where the volume fraction is not constant, (5.6) will include terms that involve derivatives of the volume fraction, and the continuity equation

(5.9)\begin{equation} \frac{\partial\phi}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\phi\boldsymbol{v}_p) = 0 \end{equation}

needs to be solved for the dynamics of $\phi$ with $\boldsymbol {v}_p$ given by (5.3).

In order to close this system of equations and solve for the dynamics of a suspension in a specified situation, boundary conditions need to be defined. However, our averaging method defines the average fluid and particle velocities at position $\boldsymbol {x}$ using a spherically symmetric averaging kernel of effective size $R$ centred about this location. Therefore, these averages are only well defined in regions where $\boldsymbol {x}$ is sufficiently far away from any boundaries (at least a distance $R$ or greater). Because of this limitation, we are not currently able to define the proper boundary conditions for the suspension, but we are actively working on the proper way to do this. Even so, it is obvious that our model predicts an additional correction to the viscosity beyond the Einstein result, since our effective viscosity $\eta _{{eff}}$ includes an additional factor of $1\ (1-\phi )$. In figure 1, we plot this function against the following experimental data: Eilers (Reference Eilers1941); table 3 of Ward & Whitmore (Reference Ward and Whitmore1950); Maron & Fok (Reference Maron and Fok1955); and Maron & Levy-Pascal (Reference Maron and Levy-Pascal1955). In this figure, we also compare our equation with the second-order expression calculated by Batchelor & Green (Reference Batchelor and Green1972) who predicted $\eta _{{eff}}\ \eta _0 = 1 + 2.5\phi + 5.2\phi ^2$ and with the empirically determined equation by Mooney (Reference Mooney1951) that estimated $\eta _{{eff}}\ \eta _0 = \exp (K_1\phi \ (1-K_2\phi ))$, where $K_1$ is typically between 2 and 5 and $K_2$ is between 0.6 and 2 (Pal Reference Pal1990). The effective viscosity that we predict works well out to 15 % volume fractions, but then underestimates the effective viscosity. It is likely that the deviation between our result and the experimental data at higher volume fractions is due to our $k=2$ truncation, which ignores lubrication forces that become more important as the suspended particle density increases. It is interesting to note that this is also the point where the experimental data begin to spread out more, suggesting that the experimental conditions may have a strong effect on the measured viscosity. There is also another, potentially more important, reason that our results might underestimate the effective viscosity at higher volume fractions. Our analysis assumes that the pair-distribution function is isotropic, which leads to the integral in (4.8) being zero. Consequently, there are no additional corrections to the viscosity. In the presence of shear, however, it is expected that the pair-distribution function is not isotropic (Batchelor Reference Batchelor1977; Brady & Morris Reference Brady and Morris1997), which could lead to additional corrections to the viscosity. The consequences of a shear-induced, non-isotropic distribution function are described in § 7.

Figure 1. Comparison of experimental measurements of suspension viscosity from Eilers (Reference Eilers1941) (black squares), from table 3 of Ward & Whitmore (Reference Ward and Whitmore1950) (other open symbols), from Maron & Fok (Reference Maron and Fok1955) (blue filled circles) and from Maron & Levy-Pascal (Reference Maron and Levy-Pascal1955) (filled green circles) with the results of our theory for the effective viscosity $\eta _{{eff}}$ (red curve). As a comparison with other theories, we have also included the Batchelor & Green (Reference Batchelor and Green1972) equation (black line) and the Mooney equation with $K_1 = 2.5$ and $K_2 = 1.0$ (blue line).

6. The speed of sedimenting particles

A long-standing problem in suspension dynamics is the question of the correct theoretical description of the speed for sedimenting particles. To date, one of the most widely cited expressions relating the settling speed $U_S$ to the volume fraction is the empirically determined Richardson–Zaki equation:

(6.1)\begin{equation} U_s = \left(1-\phi\right)^n U_0, \end{equation}

where $n$ is an empirically determined parameter that is found to approach 4.65 as the Reynolds number goes to zero (Richardson Reference Richardson1954; Richardson & Zaki Reference Richardson and Zaki1954; Silva et al. Reference Silva, Garcia, Faia and Rasteiro2015). This exponent does not exactly agree with the aforementioned calculation in Batchelor (Reference Batchelor1972) that derived $U_S \approx (1 - 6.55\phi )U_0$. Some authors also use a modified form of the Richardson–Zaki equation, where the volume fraction is scaled by the maximum random packing fraction; however, this introduces an extra parameter that can vary between theories and experiments (Silva et al. Reference Silva, Garcia, Faia and Rasteiro2015).

For the case of sedimentation in a stationary fluid where $\left \langle \boldsymbol {v}\right \rangle = 0$, (5.3) can be used to determine the settling speed predicted by our analysis for a given pair-distribution function. While there has been significant debate in the literature about what the appropriate distribution function should be (see e.g. Batchelor Reference Batchelor1972; Cichocki & Sadlej Reference Cichocki and Sadlej2005), we make the simplest assumption and use the hard-sphere pair-distribution function. This choice is convenient because the pair distribution is isotropic and the functional dependence of the distribution function on the volume fraction has been well studied. For an isotropic $g(r)$, the sedimentation speed is found from (5.3) to be

(6.2)\begin{equation} U_S = \left( 1 - \frac{1}{2}\phi - \frac{1}{2}\phi^2 + \frac{3\phi}{a^2}\int_{a}^{2R} r \left(g\left(r\right) - 1\right) {\rm d}r\right)U_0, \end{equation}

where we have ignored terms of order $a\ R$ and smaller. In general, the pair-distribution function is zero for $r<2a$ and can be written as an expansion in $\phi$:

(6.3)\begin{equation} g\left(r\geqslant2a\right) = 1 + \sum_{n=1}^\infty \phi^n g_n\left(r\right). \end{equation}

The functions $g_n$ in this expansion have been determined analytically for $n = 1,2$ (Kirkwood Reference Kirkwood1935; Nijboer & Van Hove Reference Nijboer and Van Hove1952) and computed numerically up to $n = 6$ (Ree, Keeler & McCarthy Reference Ree, Keeler and McCarthy1966; Kolafa & Labik Reference Kolafa and Labik2006). We use the analytic expression for $g_1$ (Kirkwood Reference Kirkwood1935) to compute the first-order correction to the integral in (6.2). For $g_2$, the analytic expression is sufficiently complex that we just use the tabulated values given in Nijboer & Van Hove (Reference Nijboer and Van Hove1952) and use the midpoint rule to integrate $g_2$ numerically. Similar numerical integrations are done using the tabulated values of $g_3$ from Ree et al. (Reference Ree, Keeler and McCarthy1966) and for $g_4$, $g_5$ and $g_6$ from Kolafa & Labik (Reference Kolafa and Labik2006). Our analysis then predicts that the hindered settling function is

(6.4)\begin{align} \frac{U_S}{U_0}& = 1 - 5\phi + \frac{127}{10}\phi^2 - 25.76\phi^3 + 45.53\phi^4 - 73.06\phi^5 + 108.65\phi^6 - 151.70\phi^7\nonumber\\ &\quad + {{O}}(\phi^8). \end{align}

In order to compare the sedimentation speed predicted from (6.4) with experiments, we use the meta-analysis data in Brzinski & Durian (Reference Brzinski and Durian2018) that compiled sedimentation results from 14 different publications and also included three new measurements (figure 2). In that work, Brzinski & Durian (Reference Brzinski and Durian2018) discovered two separate branches in the hindered settling function, depending on whether the Péclet number was greater or less than ${{O}}(10^8)$.The agreement between our model and the high-Péclet-number data is quite good out to around 30 %, beyond which the truncated expansion in powers of $\phi$ breaks down. This equation also compares extremely well with the empirical Richardson–Zaki equation with $n=4.65$ (figure 2). It is somewhat surprising that this analysis works so well using the hard-sphere pair-distribution function. It may be that due to random movements of the spheres, they behave sufficiently like a gas of particles that this pair distribution is reasonable. In order to explore how well the pair-distribution function does at higher volume fractions, we used the Monte Carlo and molecular dynamics simulations data from Kolafa, Labik & Malijevsky (Reference Kolafa, Labik and Malijevsky2002, Reference Kolafa, Labik and Malijevsky2004) that compute the pair-distribution functions at various volume fractions ranging from 10 % to 53 %. These functions were then integrated numerically to determine $U_S/U_0$ at specific values of $\phi$ (figure 2). From these results, we find that our results with the hard-sphere pair-distribution function agree well with the Richardson–Zaki equation out to around 40 %, at which point our results drop to zero near 50 %, whereas the Richardson–Zaki equation approaches zero as $\phi \rightarrow 1$. Above 40 %, both the Richardson–Zaki and the hard-sphere pair-distribution function results fall below the data from Brzinski & Durian (Reference Brzinski and Durian2018).

Figure 2. Comparison of our theory with experimental data on sedimentation speed from the meta-analysis in Brzinski & Durian (Reference Brzinski and Durian2018), for high Péclet numbers ($Pe > {{O}}(10^8)$) (black triangles) and lower Péclet numbers ($Pe <{{O}}(10^8)$). The red curve shows (6.4) and the yellow squares are computed using our equation with the discrete values of $g(r)$ from the Monte Carlo data in Kolafa et al. (Reference Kolafa, Labik and Malijevsky2002, Reference Kolafa, Labik and Malijevsky2004). We also compare these data with the Batchelor results (black solid and dashed lines) and the Richardson–Zaki equation with $n=4.65$ (green circles). Our theory agrees well with the high-Péclet-number data. The inset shows a log–linear plot of the same comparisons, with the Richardson–Zaki function omitted for clarity.

A weakness of our hindered settling function is that the polynomial expansion is cumbersome and has a limited range of validity. While the Richardson–Zaki equation is often used and agrees well with our analytic expression, one fault with this function is that it does not go to zero at the correct volume fraction. To address these concerns, we found that our analytic expression and the hindered settling function computed with the numerically determined hard-sphere pair distribution are very well fitted by

(6.5)$$\displaystyle{{U_S} \over {U_0}} = \displaystyle{{1-6\phi {\rm /}\pi } \over {{\left( {1 + \phi } \right)}^3}}.$$

The term in the numerator goes to zero at the cubic packing limit, $\phi \approx 0.52$. It should be noted that on a term-by-term basis, this function is not the same as our polynomial expansion; nevertheless, it closely follows (6.4) and provides a simple expression that can be used up to the cubic packing limit.

It is worthwhile to discuss the difference between our results and those from Batchelor (Reference Batchelor1972). Indeed, our leading-order result of $1 - 5\phi$ is identical to the Batchelor result before he corrected for the flows to satisfy boundary conditions and for exactly the same reasons: (i) in the dilute limit, when $g(r) = 1$ for $r<2a$ (which is what Batchelor used), the integral in (6.2) is simply over the single-sphere disturbance flow from $a< r<2a$, which contributes a factor of $-4.5\phi$ to the hindered settling function; (ii) the restriction that the net flow is zero leads to an additional factor of $-\phi$; and (iii) the average Laplacian term from the Faxen force contributes a term $\phi \ 2$. In Batchelor's procedure, though, he considered removing an arbitrarily chosen test sphere from the suspension and replacing it with a uniform velocity pure fluid. He argued that he then needed to correct for the changes to the boundary conditions on the other spheres to account for the missing test sphere. Using the two-particle sedimentation interaction, this procedure led to a flow correction of $-1.55\phi$, and a total linear-order hindered settling function of $1 - 6.55\phi$. However, the test sphere is an idealization for which its removal should leave the average statistical distribution of the suspension unaffected. That is, the remaining spheres should ‘see’ the same distribution of particles with or without the test sphere. Consequently, the removal of the test sphere should not alter the flow experienced by the remaining spheres, and, therefore, no corrections to the sedimentation speed need be considered. By accounting for alterations due to the test sphere, the statistical distributions between the two scenarios are not the same, and the system including the test sphere must have a larger volume fraction. Thus, Batchelor's calculation overestimates the leading order of the hindered settling function.

In our analysis, rather than removing or introducing a test sphere, we considered the background flow, $\boldsymbol {v}_{-m}$, that a given sphere experiences due to the presence of all the other spheres. As mentioned already, this is not the flow that would exist if the sphere were removed, but is rather the real background flow at the current sphere, caused by the disturbance flows from all the other spheres. We define these disturbance flows using Green's functions that automatically satisfy the boundary conditions on the surfaces of the spheres, when the background flow is known in a Taylor-series-type expansion. Hence, we are considering the average real flow about an average sphere in the suspension, with no need to correct the forces on the other spheres to account for the presence of this sphere. Its effect is already included.

7. The effect of shear flows on neutrally buoyant suspensions

Shear flows lead to anisotropy in the pair distribution (Gadala-Maria & Acrivos Reference Gadala-Maria and Acrivos1980; Parsi & Gadala-Maria Reference Parsi and Gadala-Maria1987; Brady & Morris Reference Brady and Morris1997). In the context of the work presented here, this anisotropy does not alter the averaged Stokes equation (2.12) nor does it affect the Faxen laws (3.2)–(3.4); however, it does change the integrals that relate the average local background flow about the sphere $\left \langle \boldsymbol {v}_{-m}\right \rangle$ to the average fluid velocity $\boldsymbol {v}_f$. To examine the consequences in a specific system, we consider a suspension of neutrally buoyant spheres in shear flow directed along the $x$ axis with shear rate ${\dot {\gamma }}$ along the $z$ axis. Our choice of a neutrally buoyant suspension sets ${\hat {\boldsymbol {F}}^0} = 0$. In addition, the symmetry of the pair-distribution function under shear, $g(\boldsymbol {x}) = g(-\boldsymbol {x})$, maintains that odd and even $k$ terms do not affect one another. Therefore, to the order of our truncation, the only modification to our results comes from the integral in (4.8), for which we need $g(\boldsymbol {r})$. Using the diffusion coefficient for the spheres $D = k_BT\ 6{\rm \pi} \eta _0 a$, where $k_BT$ is thermal energy, we define a Péclet number $Pe = {\dot {\gamma }}a^2\ 2D$. In the dilute limit, the pair-distribution function is then defined by the two-particle Smoluchowski equation (Brady & Morris Reference Brady and Morris1997):

(7.1)\begin{equation} \nabla^2 g - Pe z \frac{\partial g}{\partial x} = 0, \end{equation}

with boundary conditions

(7.2)\begin{equation} \left.\begin{gathered} \left.\frac{\partial g}{\partial r}\right|_{r=2} = 2 Pe \sin\theta\cos\theta\cos\phi g(r=2),\\ g(r=\infty) = 1, \end{gathered}\right\}\end{equation}

where $\theta$ is the azimuthal angle from the $z$ axis, $\phi$ is the polar angle and we have non-dimensionalized length by the sphere radius. We solve this equation numerically using a finite volume spatial discretization on a grid in spherical coordinates. Due to the symmetry of the flow, we solve on a hemispherical domain defined by $0\leq \theta \leq {\rm \pi}$ and $-{\rm \pi} /2 \leq \phi \leq {\rm \pi}/2$. Since the distance over which $g$ varies depends on $Pe$, we use a hemispherical radius $2\leq r\leq R_{max}$, where $R_{max} = 16$ for $Pe \leq 1$ and $R_{max} = 2.15 + 30\ (1+Pe)$ otherwise, where this range was empirically determined by requiring that the integral of $g-1$ over the domain had sufficiently converged. We discretized angular coordinates with 45 nodes, the radial coordinate with 200 nodes and then solved (7.1) over a range of Péclet numbers from 0.1 to 1000. Our results for the pair-distribution function show the expected increase in $g$ at the boundary in the compressional quadrants and decrease in the extensional quadrants and also agree with previously reported solutions of (7.1) (figure 3; Koo & Hess Reference Koo and Hess1987; Brady & Morris Reference Brady and Morris1997).

Figure 3. Numerical solution of the pair-distribution function at four values of the Péclet number: $Pe = 0.1$ (a), $1.0$ (b), $10$ (c) and $100$ (d). The cross-section shows the pair distribution in the $x$–$z$ plane. Only half the radial domain is shown in order to highlight the distribution near the boundary.

Using the numerically determined value of $g$, we then computed

(7.3)\begin{equation} {\hat{\boldsymbol\varGamma}} = \int_{r = 2a}^{r = R_{max}} \left( g\left(\boldsymbol{r}\right) - 1 \right)\boldsymbol{\nabla}{\hat{\boldsymbol{S}}^1}\left(\boldsymbol{r}\right){\rm d}^3r, \end{equation}

by integrating over the full spherical domain, where the lower limit on the integral is $r=2a$ because in the range $a\leq r\leq 2a$ where $g(\boldsymbol {r}) = 0$ the integral is zero. Note that ${\boldsymbol \varGamma }$ is a rank-4 tensor. From the symmetry of the velocity gradient, components that are odd in the $y$ direction are identically zero. Using (3.3) and (4.8) leads to the linear system of equations that defines the $k=1$ singularity strengths:

(7.4)\begin{equation} \left( \mathbb{I} + \frac{5\phi}{16{\rm \pi}\left(1-\phi\right)}{\hat{\boldsymbol{\varGamma}}} \right)\odot {\hat{\boldsymbol{F}}^1} ={-}\frac{10{\rm \pi}\eta_0 a^3{\dot{\gamma}}}{3\left(1-\phi\right)}( {\hat{\boldsymbol{x}}}{\hat{\boldsymbol{z}}} + {\hat{\boldsymbol{z}}}{\hat{\boldsymbol{x}}} ), \end{equation}

where $\mathbb {I}$ is the rank-4 identity tensor. The only non-zero components of ${\hat {\boldsymbol {F}}^1}$ are the diagonal, $xz$ and $zx$ components, and the effective stress tensor is

(7.5)\begin{equation} \sigma^{e}_{ij} = \eta_0\left(\frac{\partial\left\langle v\right\rangle _i}{\partial x_j} + \frac{\partial\left\langle v\right\rangle _j}{\partial x_i}\right) - P\delta_{ij} -\frac{3\phi}{4{\rm \pi} a^3}{\hat{F}^1_{ij}}. \end{equation}

The effective viscosity is then given by $\eta = \sigma _{xz}\ \dot {\gamma }$ and the first and second normal stress differences are $\sigma _{xx}-\sigma _{zz}$ and $\sigma _{xx} - \sigma _{yy}$, respectively. For the effective viscosity, we compare the results with our previous results using the isotropic pair distribution, $\eta _{{eff}} = 1+ 5\phi \ 2(1-\phi )$. As shown in figure 4(a), the effect is largest at moderate values of $Pe$ and at larger values of $\phi$, but over the full range of $Pe$ and for $\phi < 0.3$, the viscosity varies by no more than 1 % from the result predicted with the isotropic pair-distribution function. The absence of a significant alteration to the viscosity could be due to the fact that we are calculating the pair distribution in the dilute limit. We do, however, find that the shear pair-distribution function leads to non-negligible first and second normal stress differences (figure 4b,c), with both being negative. Similar to the viscosity, the normal stress differences are largest at moderate values of $Pe$, with a peak near $Pe\approx 8$ and increase with $\phi$. Both also go to zero as $Pe \rightarrow 0$ and asymptote to a non-zero negative value as $Pe \rightarrow \infty$. Fitting the high-$Pe$ results to a quadratic polynomial, we estimate that $\sigma _{xx} - \sigma _{zz} \approx -\eta _0{\dot \gamma }(0.057\phi +0.37\phi ^2)$ and $\sigma _{xx}-\sigma _{yy} = -\eta _0{\dot {\gamma }}(0.32\phi +2.0\phi ^2)$ at large $Pe$.

Figure 4. Fluid stresses as a function of volume fraction and Péclet number. (a) The ratio of predicted viscosity to $\eta _{eff}$. The shear flow pair-distribution function only alters the viscosity by small amounts at moderate values of $Pe$ and larger volume fractions. The non-isotropic, shear flow pair distribution produces both first (b) and second (c) normal stress differences. Both are negative and are largest at moderate values of $Pe$. The colourbars for (b,c) are the same.