Continuum dynamics of suspensions at low Reynolds number

Abstract

The dynamics of suspensions of particles has been an active area of research since Einstein first calculated the leading-order correction to the viscosity of a suspension of spherical particles (Einstein, Proc. R. Soc., vol. A102, 1906, pp. 161–179). Since then, researchers have strived to develop an accurate description of the behaviours of suspensions that goes beyond just leading order in the particle volume fraction. Here, we consider the low-Reynolds-number behaviour of a suspension of spherical particles. Working from the Green's functions for the flow due to a single particle, we derive a continuum-level description of the dynamics of suspensions. Our analysis corrects an error in the derivation of these equations in the work of Jackson (Chem. Engng Sci., vol. 52, 1997, pp. 2457–2469) and leads to stable equations of motion for the particles and fluid. In addition, our resulting equations naturally give the sedimentation speed for suspended particles and correct a separate error in the calculation by Batchelor (J. Fluid Mech., vol. 52, 1972, pp. 245–268). Using the pair-correlation function for hard spheres, we are able to compute the sedimentation speed out to seventh order in the volume fraction, which agrees with experimental data up to 30 %–35 %, and also get higher-order corrections to the suspension viscosity, which agree with experiments up to $\sim$15 %. Then, using the pair distribution for spheres in shear flow, we find alterations to both the first and second normal stresses.

1. Introduction

Mixing small solid particles or phase-separated liquid or gas droplets into a bulk liquid creates a fluid-like suspension that, while behaving somewhat like a simple fluid, also has its own unique properties. Think, for instance, of the motion of mud or paint relative to that of water or honey. Or the more dramatic motion of bubbles in a pint of Guinness, where the small and lighter-than-water nitrogen bubbles flow downward along the edge of the pint glass due to the interplay between the pint geometry and the Boycott effect (Benilov, Cummins & Lee 2011). Suspensions are ubiquitous, from naturally occurring ones, such as sand in water, to a range of various culinary sauces and pastes (dressings, pestos, soups, etc.) to the industrial realm of printer inks and oil slurries. Therefore, understanding their physical properties and having an accurate physical theory of their dynamics are essential. It is not surprising that a vast extent of research over more than a hundred years has sought to make progress in this area.

The first calculation to explore the physical consequences of suspended particles in a fluid was Einstein's classic analysis of the change in the effective viscosity of a fluid due to the presence of a low volume fraction of spheres (Einstein 1906). For a volume fraction $\phi$ of spheres of radius $a$ suspended in a fluid of viscosity $\eta _0$, the leading-order correction to the viscosity was found to be $\eta \ \eta _0 = 1 + 5\phi \ 2$. Since Einstein, other authors have found the same result using differing approaches (Jeffrey 1922; Burgers 1938; Brenner 1958). These calculations all ignore the interactions between the spheres in the suspension. To get to higher-order contributions to the viscosity, these interactions need to be taken into account. Various methods have been used to accomplish this, such as using the method of reflections (Guth & Simha 1936), computation of the added stress from the particles (Batchelor 1977; Brady & Morris 1997) and cage models (Happel & Brenner 1983). All of these methods compute the effective viscosity up to second order in the volume or beyond, but all differ in what value they find for the second-order coefficient. Other groups have used fits to empirical data to determine approximate ad hoc functional relationships for the dependence of suspension viscosity on particle volume fraction (see e.g. Finke et al. 2021; Konijn, Sanderink & Kruyt 2014). Therefore, even in simple suspensions of spherical particles, the proper relationship between the volume fraction and the effective viscosity remains unclear.

Other calculations for suspension flow are equally ambiguous. For example, at low Reynolds number an isolated sphere of radius $a$ and density $\rho _s$ falls through a fluid of viscosity $\eta _0$ and density $\rho _f$ at a speed $U_0 = 2a^2g(\rho _s-\rho _f)\ 9 \eta _0$, where $g$ is the acceleration due to gravity. A suspension of similar spheres, however, falls through the same fluid at a speed $U_S$, which is slower than that for the single sphere, and this hindrance increases nonlinearly with increasing volume fraction. To date, there is no consensus on what the leading-order correction to the sedimentation speed as a function of volume fraction is. The original analysis of this problem was also initiated over a hundred years ago in the work of Smoluchowski (1912). As detailed in Batchelor (1972), the theoretical work on this problem has fallen into three main categories. First are methods using regular geometrical arrangements of particles, such as the work by Hasimoto (1959). The second category involves using ‘cell’ models, where the motion of a single particle within an isolated region is computed. The influence of the other particles is then incorporated through the choice of boundary conditions at the edge of the isolation region. The book by Happel and Brenner describes these models in detail (Happel & Brenner 1983). Both of these types of calculations lead to a first-order effect that scales like $\phi ^{1\ 3}$, since each sphere is confined within a region of volume $a^3\ \phi$. Experimental data, though, do not support that the leading-order hindrance to sedimentation scales so steeply with volume fraction. The third type of calculations that have been performed consider a random distribution of spheres. These calculations generally lead to an order $\phi$ correction to the sedimentation speed. The most notable of these calculations is that by Batchelor (1972). In that paper, Batchelor considered the average effect on the sedimentation speed due to a uniform suspension of spheres under the stipulation that the average total flow (from both the fluid and the particles) is zero. Taking only these effects into account, he computed the relative sedimentation speed to be $U_S\ U_0 = 1 - 5\phi$. Batchelor then considered removing a test sphere embedded within the flow. However, because the test sphere also creates its own disturbance flow, Batchelor reasoned that he needed to correct for the alterations to the boundary conditions on the other spheres, which he did using a two-particle reflection method. With these interactions taken into account, the leading-order correction became $U_S\ U_0 \approx 1 - 6.55\phi$ for monodisperse suspensions and $U_S\ U_0 = 1 - 5.6\phi$ for polydisperse suspensions (Batchelor & Wen 1982). This result, though, does not agree well with data, which has led other groups to posit that hydrodynamic effects may alter the pair-distribution function, which could account for the experimentally observed deviation from Batchelor's result (Cichocki & Sadlej 2005). Many other groups have offered empirically determined functional forms to try to describe the dependence of sedimentation speed on volume fraction (see e.g. Amin, Girolami & Risso 2021).

The aforementioned theoretical calculations on viscosity and sedimentation speed each addressed specific aspects of suspension dynamics. A more general question is whether it is possible to write out a system of equations that accurately describes the physics and dynamics of suspensions. Suspensions fall into the category of multiphase fluids. Two separate methodologies have been employed to study these systems. The first category involves using direct simulations to compute the motions of discrete particles moving within a continuum fluid, such as described in Yin et al. (2019). The second method considers volume elements that are large compared with the size of any of the components within the multiphase fluid and averages over the dynamics of each phase. In this way, each phase becomes an interpenetrating continuum defined by its averaged velocity and volume fraction. This methodology was pioneered by Drew & Segel (1971), and has been an active area of interdisciplinary research since (see e.g. Cogan & Guy 2010; Kim & Mudawar 2015). For a suspension composed of solid particles at volume fraction $\phi$ in a fluid of viscosity $\eta _0$ at low Reynolds number, a force balance equation is needed to describe the net force per volume on the spheres and a separate force balance equation is needed for the force per volume on the fluid. In general, each phase can feel intra-phase forces, inter-phase forces and body forces. Deriving the intra- and inter-phase forces from first principles can be difficult or impossible. Often, some form of approximate closure relationship is proposed. In the context of a dilute suspension of spherical particles at low Reynolds number, Jackson (1997) derived an explicit closure relationship by averaging the Stokes equations and the forces exerted on the particles from the fluid. While these equations were derived directly from well-defined physics, Jackson (1997) found the equations to be ill-posed and to lead to unstable solutions, even for simple one-dimensional fluidization/sedimentation problems.

Here we aim to develop a physically correct continuum theory for two-phase suspensions of spherical particles at low Reynolds number. We begin by considering the Green's function for a sphere in an arbitrary flow field, which can be written as a sum over singular solutions. Each sphere in the suspension can then be described by a similar Green's function. By averaging the Stokes equations over an intermediate length scale and using these Green's functions, we show that if the singularity strengths are known, then the effective stress tensor for the fluid can be defined. In order to compute the singularity strengths, we use the Taylor series expansion of the fluid velocity about the centre of a given sphere to set the singularity strengths on an order-by-order basis. This leads to a set of Faxen law-type equations that define the singularity strengths in terms of the local background flow of the fluid and the disturbance velocities due to all the other spheres. To close the system of equations, we relate the average local fluid velocity about a sphere to the fluid-phase averaged velocity. The resulting solution provides the effective stress tensor for the averaged fluid velocity, correcting the original calculation in Jackson (1997) and leading to stable equations of motion for the suspension.

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 2005)

(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 (1970) and later by Pedley & Kessler (1992) in the context of microorganism swimming. This same term also accounts for the primary difference between our results and the results from Jackson (1997), 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 1985; Chen & Ye 2000), 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 1985; Chen & Ye 2000). 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 (1997),

(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 1935; Guazzelli & Morris 2011), 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 1997), 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 1997). 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 (1972).

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 (1997) 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 (1941); table 3 of Ward & Whitmore (1950); Maron & Fok (1955); and Maron & Levy-Pascal (1955). In this figure, we also compare our equation with the second-order expression calculated by Batchelor & Green (1972) who predicted $\eta _{{eff}}\ \eta _0 = 1 + 2.5\phi + 5.2\phi ^2$ and with the empirically determined equation by Mooney (1951) 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 1990). 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 1977; Brady & Morris 1997), 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 (1941) (black squares), from table 3 of Ward & Whitmore (1950) (other open symbols), from Maron & Fok (1955) (blue filled circles) and from Maron & Levy-Pascal (1955) (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 (1972) 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 1954; Richardson & Zaki 1954; Silva et al. 2015). This exponent does not exactly agree with the aforementioned calculation in Batchelor (1972) 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. 2015).

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 1972; Cichocki & Sadlej 2005), 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 1935; Nijboer & Van Hove 1952) and computed numerically up to $n = 6$ (Ree, Keeler & McCarthy 1966; Kolafa & Labik 2006). We use the analytic expression for $g_1$ (Kirkwood 1935) 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 (1952) 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. (1966) and for $g_4$, $g_5$ and $g_6$ from Kolafa & Labik (2006). 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 (2018) that compiled sedimentation results from 14 different publications and also included three new measurements (figure 2). In that work, Brzinski & Durian (2018) 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 (2002, 2004) 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 (2018).

Figure 2. Comparison of our theory with experimental data on sedimentation speed from the meta-analysis in Brzinski & Durian (2018), 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. (2002, 2004). 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 (1972). 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 1980; Parsi & Gadala-Maria 1987; Brady & Morris 1997). 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 1997):

(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 1987; Brady & Morris 1997).

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.

8. Conclusion

Here we have addressed the problem of determining an accurate physical description of the phase-averaged dynamics in a suspension composed of spherical particles immersed in a viscous fluid. We used a sequence of Green's functions centred about each sphere to account for how the force distribution on the surface of the sphere affects the fluid. This sequence of Green's functions is constructed such that the strength of each term in the sequence for a given sphere corresponds to a specific term in the Taylor series expansion of the background flow about that sphere. Since the Green's function strengths define the net disturbance flow from all of the spheres, the strengths of the singular solutions are coupled. By averaging the flow in a given region, we derived a relationship between the strengths of the Green's functions and the fluid-averaged flow that allows us to derive an averaged equation of motion for the fluid that depends only on the overall average velocity ${\left \langle \boldsymbol {v}\right \rangle }$, its derivative and the body forces that act on the fluid and particles. Truncating the Taylor series expansion of the background flow at second order gives a modified Stokes equation that includes up to fourth derivatives of the velocity and predicts an effective viscosity of $1 + 5\phi \ 2(1-\phi )$ for a uniform suspension described by an isotropic pair-distribution function. The term in the denominator represents a correction to the standard Einstein effective viscosity. The Einstein correction is only considered to be valid up to a few per cent in volume fraction, whereas our expression seems to be valid up to around 15 % volume fractions. We believe that the deviation between our result and experimental data is driven by two factors. First, our expression is partially the result of an isotropic pair-distribution function, which is not expected to be correct in the case of shear flow. Second, the spread in the experimental data at higher volume fractions suggests that experimental conditions (e.g. shear rate, particle size, etc.) may also play significant roles at higher volume fraction.

We used the resulting theory to examine a few representative problems. First, we considered the uniform sedimentation of spheres. We found that the dominant coupling between the spheres comes from an integral over the zeroth-order Green's function and the pair-distribution function for the spheres. We made the simplest assumption that the hard-sphere pair distribution is representative of the actual distribution of spheres in the suspension. Using the analytic and computed density expansion for $g$, we obtained an expression relating the sedimentation speed to a power-series expansion in the volume fraction up to seventh order that agrees surprisingly well with previous experiments and also with the empirically determined Richardson–Zaki equation (6.1) up to 30 %. In addition, we propose a simple approximate expression for the hindered settling function (6.5) that agrees well with experimental observations and our theoretical analysis and also vanishes at a physically meaningful volume fraction. Finally, we examined the case of a non-isotropic pair-distribution function in the context of dilute shear flow. We find that the shear pair-distribution function leads to negligible changes to the effective viscosity, but does produce non-zero first and second normal stress differences.

Theoretical determination of the sedimentation speed and effective viscosity of suspensions has been a long-sought grail in fluid mechanics. In most cases, researchers have sought to address either one of these questions or the other. However, both are the consequence of the same physical system, which suggests that they should be explainable by the same theory. Our approach unifies the computation of these behaviours within the context of the averaged two-phase dynamics of a suspension.

Our averaged equations of motion are fourth order in the overall average velocity, thereby requiring two boundary conditions per dimension at each boundary. This result is convenient, as these equations determine the dynamics of both the fluid and particle velocities; hence, one boundary condition for each component of the velocity. Unfortunately, it remains unclear how to set appropriate boundary conditions within the context of the averaging approach that we are using. While it is possible that the fluid phase maintains a no-slip condition at the boundary, careful consideration of how the averaging should be done near the boundary and the consequences of this for the boundary conditions remains an open question. Even so, the framework that is developed here provides a foundation for exploring a large range of problems in the dynamics of suspensions.