2.1. Model systems

We model the dynamics of a suspension of hard colloidal spheres of varying sizes $a_i$ suspended in a Newtonian solvent of density $\rho$ and viscosity $\eta$, all enclosed in a hard spherical cavity of size $R$. Here, $i$ refers to particles of a given size. The motion of the particles sets the surrounding fluid into motion, where the small size of the colloids produces a vanishingly small Reynolds number $Re=\rho U a/\eta$, where $U$ is the characteristic particle velocity. As a result, fluid motion is governed by the Stokes equations, and because the particles are freely mobile, their motions are coupled to each other and with the enclosing cavity by many-body reflected hydrodynamic and lubrication interactions. Accounting for Brownian motion, hydrodynamic interactions and other deterministic forces in colloidal suspensions can be done analytically for a pair of particles, but requires computational methods beyond the semi-dilute limit. Properly accounting for Stokesian particle dynamics in such systems requires two challenges to addressed: first, representing near-field hydrodynamic interactions near contact, which become divergent lubrication forces; the second challenge is representing infinitely reflected fluid-mediated interactions between many particles simultaneously, i.e. many-body hydrodynamic interactions. One approach for addressing these challenges is the Stokesian dynamics algorithm (Durlofsky et al. Reference Durlofsky, Brady and Bossis1987; Bossis & Brady Reference Bossis and Brady1989; Sierou & Brady Reference Sierou and Brady2001; Banchio & Brady Reference Banchio and Brady2003), developed for unconfined suspensions, which leverages the linearity of Stokes flow to solve the near-field and many-body problems separately and then superimpose them for a complete solution.

CSD expands that framework by incorporating a spherical enclosure, enabling measurement of its effect on particle microstructure and hydrodynamic interactions for monodisperse (Aponte-Rivera & Zia Reference Aponte-Rivera and Zia2016; Aponte-Rivera et al. Reference Aponte-Rivera, Su and Zia2018) and polydisperse (Gonzalez et al. Reference Gonzalez, Aponte-Rivera and Zia2021) suspensions. To do so, we expanded mobility functions derived from Ladyzhenskaya's integral equation (Ladyzhenskaya Reference Ladyzhenskaya1963) which we combined with Faxén relations to represent the enclosure. These functions were superimposed with previously developed near-field theory (Jeffery Reference Jeffery1915; Stimson & Jeffery Reference Stimson and Jeffery1926; Dean & O'Neill Reference Dean and O'Neill1963; O'Neill Reference O'Neill1964; Majumdar Reference Majumdar1967; Cooley & O'neill Reference Cooley and O'neill1968; O'Neill & Majumdar Reference O'Neill and Majumdar1970a; Jeffrey & Onishi Reference Jeffrey and Onishi1984; Jeffrey Reference Jeffrey1992; Jones Reference Jones2009). We used the resulting CSD model to study equilibrium diffusion in concentrated confined monodisperse and polydisperse suspensions (Aponte-Rivera et al. Reference Aponte-Rivera, Su and Zia2018; Gonzalez et al. Reference Gonzalez, Aponte-Rivera and Zia2021). This model was limited to equilibrium conditions however, because some elements of the hydrodynamic functions were not yet developed: the coupling between particle motion and the hydrodynamic stresslet, the symmetric first moment of the hydrodynamic traction on particle surfaces. We compute these new hydrodynamic functions in the present work. We get started by defining coordinate reference frames that facilitate solving the near-field and many-body problems for confined suspensions. Two coordinate systems are needed and are described here because they will be used throughout subsequent sections to solve the Stokes equations and relate these solutions to the rheology of confined suspensions.

For pair interactions (particle–particle or particle–cavity), the reference frame is identical to an unbound pair (centred on a colloidal particle). Placing the origin at the centre of a particle facilitates transformations between Cartesian, spherical, cylindrical and bispherical coordinates needed to construct and solve the relevant equations of motion and boundary conditions which requires going back and forth between coordinate systems. To model just the pair-level hydrodynamic interactions between a colloid and the cavity, the origin of a Cartesian and curvilinear coordinate systems (cylindrical, spherical and bispherical) is located at the centre of the colloid, as shown in figure 1(a). Here, $\boldsymbol{r}$ denotes the centre-to-centre vector between the cavity and the colloid with the origin at the colloid centre. For pair theory, the geometry of the system is fully specified by the particle position and the particle-to-cavity size ratio $\lambda _c=a/R$ and $\boldsymbol{r}$, equivalently called ‘confinement ratio’ throughout this article.

Figure 1. Conceptual sketch of model systems (a) for a confined colloid in a spherical cavity with origin at the centre of the colloid and (b) for a confined polydisperse suspension in a spherical cavity with origin at the centre of the cavity, where the domain is discretised in concentric spherical bins of constant width.

In contrast, many-body interactions require placement of the origin at the centre of the confining cavity (figure 1b), as though it were itself the reference particle, which facilitates the discretisation of the volume inside the cavity throughout which particles interact. The centre-to-centre vector between the cavity and any colloid is opposite in sign to that used for particle–particle interactions. For monodisperse confined suspensions of $N$ particles, the composition of the system is fully specified by the particle-to-cavity size ratio $\lambda _c=a_1/R$ and the volume fraction $\phi = \lambda _c^3 N$. However, for a polydisperse confined suspension with $M$ different particle sizes $a_1< a_2< \cdots < a_M$, one must also specify $M-1$ particle size ratios $\lambda _{p(i)} = a_i/a_1$ and $M-1$ volume compositions $\phi _i/\phi$, where $\phi = \sum _i^M\phi _i$ and $\phi _i = 4/3{\rm \pi} a_i^3n_i$ and $n_i$ the number density of particles with radius $a_i$. To facilitate the analysis of the suspension dynamics, particles are tracked by the position of their particle centres, $\boldsymbol{r}$, and located in one of 100 concentric bins of constant width $\Delta r/R$. The edge of the outermost bin corresponds to the location of the centre of the smallest particle when it is at contact with the cavity wall and each bin is identified by its radial midpoint $r_m$ (between the innermost and outermost edges), see figure 1(b).

2.2. Mobility and resistance formulations

Viscosity, viscoelastic moduli, osmotic pressure and normal-stress differences are all elements of the stress tensor. Thus, obtaining the stress tensor is one of the first steps in characterising the complete rheology of a material. In a colloidal suspension, the stress differs from that of the suspending fluid and of the material properties of the suspended particles, as demonstrated by bulk rheological experiments (Mewis & Wagner Reference Mewis and Wagner2011). Batchelor (Reference Batchelor1970) showed that the average suspension stress can be obtained as an average over the entire volume, solvent plus particles:

(2.1)\begin{equation} \langle\boldsymbol{\sigma}\rangle =-\langle\, p\rangle_f\boldsymbol{I} + 2\eta\langle \boldsymbol{e} \rangle + \langle \boldsymbol{\varSigma}\rangle_{PP}. \end{equation}

The angle brackets signify an average over the entire suspension $\langle {\cdot }\rangle$, the fluid only $\langle {\cdot }\rangle _f$, or the particle phase only $\langle {\cdot }\rangle _{PP}$. The first term on the right-hand side gives the dynamic pressure $p$ averaged throughout the fluid phase; here, $\boldsymbol {I}$ is the identity tensor. In the second term the strain rate $e$ is averaged over the fluid and particles, for all configurational arrangements. The stress contributed by the particle phase, $\langle \boldsymbol{\varSigma }\rangle _{PP}$, is averaged over all particle configurations. Batchelor (Reference Batchelor1970) showed that the configuration average in homogeneously deformed suspensions is equivalent to a volume average. Related to this result, the particle-phase stress emerges from the interface with the fluid, producing a surface integral that depends on the configuration of the particles:

(2.2)\begin{equation} \langle\boldsymbol{\varSigma}\rangle_{PP} = \frac{1}{V}\sum_{i}^N\int_{S_{P,i}}\left[\boldsymbol{\sigma}\boldsymbol{\cdot}\boldsymbol{n}\,\boldsymbol{r} -\eta(\boldsymbol{vn}+\boldsymbol{nv}) \right] \,\mathrm{d}S. \end{equation}

Here $V$ is the total volume of the suspension, $N$ is the total number of particles, $S_{P,i}$ is the surface of particle $i$ and $\boldsymbol{n}$ is a unit normal vector pointing outward from the particle surface. We evaluate the fluid stress $\boldsymbol{\sigma }$ and velocity $\boldsymbol{v}$ at the particle surface. The first term in (2.2) is the first moment of the surface traction $\boldsymbol{t}=\boldsymbol{\sigma }\boldsymbol{\cdot}\boldsymbol{n}$ exerted by the fluid on a particle, which induces stress. The second term accounts for particle deformability that can reduce stress. The first moment of the traction can be separated into its symmetric and antisymmetric parts,

(2.3)\begin{equation} \int_{S_{P,i}}\boldsymbol{\sigma}\boldsymbol{\cdot}\boldsymbol{n}\,\boldsymbol{r} \, \mathrm{d}S= \frac{1}{2} \int_{S_{part}} \left[ \left( \boldsymbol{\sigma}\boldsymbol{\cdot}\boldsymbol{n}\,\boldsymbol{r} + \boldsymbol{r}\,\boldsymbol{\sigma}\boldsymbol{\cdot}\boldsymbol{n}\right) + \left( \boldsymbol{\sigma}\boldsymbol{\cdot}\boldsymbol{n}\,\boldsymbol{r} - \boldsymbol{r}\,\boldsymbol{\sigma}\boldsymbol{\cdot}\boldsymbol{n}\right)\right] \, \mathrm{d}S. \end{equation}

The antisymmetric part is the hydrodynamic torque tensor and exerts no dynamical influence on the suspension stress for a suspension of hard spheres. The symmetric part defines the hydrodynamic stresslet on particle surfaces:

(2.4)\begin{equation} \frac{1}{V}\sum_{i}^N \int_{S_{P,i}}\left[\tfrac{1}{2} \left( \boldsymbol{\sigma}\boldsymbol{\cdot}\boldsymbol{n}\boldsymbol{r}+\boldsymbol{r}\boldsymbol{\sigma}\boldsymbol{\cdot}\boldsymbol{n} \right) -\eta(\boldsymbol{vn}+\boldsymbol{nv}) \right]\mathrm{d}S \equiv \frac{1}{V}\sum_{i}^N\boldsymbol{S}_i^H \equiv \boldsymbol{\varSigma}^{PP}. \end{equation}

The deformability term is zero for hard particles, but we retain it in the expression for the stresslet because it is useful for imposing an infinitesimal deformation to solve the Stokes equations for a straining flow, an essential step in developing the missing hydrodynamic functions. This mathematical tool of imposing a perturbative flow to obtain hydrodynamic functions has been used previously in similar stresslet calculations for an unbound pair of hard spheres (Jeffrey Reference Jeffrey1992). The foregoing expression shows that calculation of the average suspension stress and suspension rheology thus requires calculation of the hydrodynamic stresslet in its general, configuration-dependent form and must capture the effect of hydrodynamic interactions on the total fluid stress. In Appendix A, we provide a brief recap of the moments of the surface traction that produce hydrodynamic force, torque and stresslet; for an $N$-particle suspension, the method of multipole Taylor expansion of Ladyzhenskaya's integral equation (Ladyzhenskaya Reference Ladyzhenskaya1963) for flow velocity that relates these couplings to fluid flow; and implementation of Faxèn formulae that relate particle motion to these hydrodynamic traction moments through a grand mobility matrix (GMM), $\mathcal {M}$, and, in turn, a grand resistance matrix (GRM), $\mathcal {R}$. These matrices are discussed next.

As outlined in Appendix A, the microhydrodynamics formalism provides a compact expression of the linearity of Stokes flow of a suspension relating the surface traction exerted by the fluid on particle surfaces to particle motion:

(2.5) \begin{equation} \left( \begin{array}{@{}c@{}} {\boldsymbol{F}}^H \\ {\boldsymbol{L}}^H \\ {\boldsymbol{S}}^H \\ \vdots\end{array} \right) =-\left( \begin{array}{@{}cccc@{}} {\boldsymbol{R}}^{FU} & {\boldsymbol{R}}^{F\varOmega} & {\boldsymbol{R}}^{FE} & \cdots \\ {\boldsymbol{R}}^{LU} & {\boldsymbol{R}}^{L\varOmega} & {\boldsymbol{R}}^{LE} & \cdots \\ {\boldsymbol{R}}^{SU} & {\boldsymbol{R}}^{S\varOmega} & {\boldsymbol{R}}^{SE} & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{array}\right) \boldsymbol{\cdot} \left( \begin{array}{@{}c@{}} {\boldsymbol{U}} - {\boldsymbol{U}}^\infty \\ {}\boldsymbol{\varOmega} - \boldsymbol{\omega}^\infty\\ {\boldsymbol{E}}-{\boldsymbol{E}}^\infty \\ \vdots \end{array} \right). \end{equation}

Here, ${\boldsymbol {F}}^H$, ${\boldsymbol {L}}^H$ and ${\boldsymbol {S}}^H$ are vectors and tensors containing the hydrodynamic force, torque and stresslet on the surfaces of each of the $N$ particles in the suspension. Each is linear in the boundary conditions, i.e. particle motion relative to the fluid: ${\boldsymbol {U}} - {\boldsymbol {U}}^\infty$, $\boldsymbol{\varOmega } - \boldsymbol{\varOmega }^\infty$, ${\boldsymbol {E}}-{\boldsymbol {E}}^\infty$, translation, rotation and straining motion, respectively. In the present case, there is no imposed flow, thus ${\boldsymbol {U}}^\infty = \boldsymbol{\varOmega }^\infty = \textbf {0}$ and ${\boldsymbol {E}}^\infty = \textbf {0}$. The particles can also move owing to Brownian motion and interparticle forces. Disregarding the origin of particle motion, the GRM couples the traction moments to motion. More specifically, GRM contains tensors that couple force to translation (${\boldsymbol {R}}^{FU}$), rotation (${\boldsymbol {R}}^{F\varOmega }$), straining motion (${\boldsymbol {R}}^{FE}$) and an infinitude of higher-order couplings as indicated by the ellipses (see Appendix A). Each of these block matrices within the GRM contains the geometrical transformation for all $N$ suspended particles that relate a particle motion to a surface traction moment. A similar linearity relation defines how prescribed motion produces hydrodynamic traction moments:

(2.6) \begin{equation} \left( \begin{array}{@{}c@{}} \boldsymbol{U} \\ \boldsymbol{\varOmega}\\ \boldsymbol{E} \\ \vdots \end{array} \right)=- \left( \begin{array}{@{}cccc@{}} {\boldsymbol{M}}^{UF} & {\boldsymbol{M}}^{U L} & {\boldsymbol{M}}^{US} & \cdots \\ {\boldsymbol{M}}^{\varOmega F} & {\boldsymbol{M}}^{\varOmega L} & {\boldsymbol{M}}^{\varOmega S} & \cdots \\ {\boldsymbol{M}}^{EF} & {\boldsymbol{M}}^{EL} & {\boldsymbol{M}}^{ES} & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{array} \right) \boldsymbol{\cdot} \left( \begin{array}{@{}c@{}} {\boldsymbol{F}}^H \\ {\boldsymbol{L}}^H \\ {\boldsymbol{S}}^H \\ \vdots \end{array} \right) . \end{equation}

The GRM $\mathcal {R}$ and the GMM $\mathcal {M}$ hold a reciprocal relationship

(2.7) \begin{equation} \mathcal{M}= \left( \begin{array}{@{}cccc@{}} \boldsymbol{M}^{UF} & \boldsymbol{M}^{U L} & \boldsymbol{M}^{US} & \cdots \\ \boldsymbol{M}^{\varOmega F} & \boldsymbol{M}^{\varOmega L} & \boldsymbol{M}^{\varOmega S} & \cdots \\ \boldsymbol{M}^{EF} & \boldsymbol{M}^{EL} & \boldsymbol{M}^{ES} & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{array} \right) = \left( \begin{array}{@{}cccc@{}} \boldsymbol{R}^{FU} & \boldsymbol{R}^{F\varOmega} & \boldsymbol{R}^{FE} & \cdots \\ \boldsymbol{R}^{LU} & \boldsymbol{R}^{L\varOmega} & \boldsymbol{R}^{LE} & \cdots \\ \boldsymbol{R}^{SU} & \boldsymbol{R}^{S\varOmega} & \boldsymbol{R}^{SE} & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{array}\right)^{-1} =\mathcal{R}^{-1}, \end{equation}

so generating either provides access to the other.

To extract motion along and transverse to the line of centres $\boldsymbol {\hat {r}}$ connecting interacting particles, each tensor within the GRM can be projected through orthogonal bases, $\boldsymbol {\hat {r}}\boldsymbol {\hat {r}}$ and $\boldsymbol {I} -\boldsymbol {\hat {r}}\boldsymbol {\hat {r}}$ that set the direction of hydrodynamic coupling, each scaled by a configuration-dependent scalar function that sets the decay of the coupling with separation distance. We illustrate this result with the familiar couplings between translation and rotation that induce hydrodynamic translation and rotational drag:

(2.8)$$\begin{gather} R_{ij}^{FU}=6{\rm \pi}\eta a\left[ X^A(\mathcal{C})\hat{r}_i\hat{r}_j+ Y^A(\mathcal{C})\left(\delta_{ij}-\hat{r}_i\hat{r}_j\right) \right], \end{gather}$$

(2.9)$$\begin{gather}R_{ij}^{LU}=\left(R_{ij}^{F\varOmega}\right)^{\rm T}= 6{\rm \pi}\eta a^2\left[ Y^B(\mathcal{C})\epsilon_{ijk}\hat{r}_k\right], \end{gather}$$

(2.10)$$\begin{gather}R_{ij}^{L\varOmega} = 6{\rm \pi}\eta a^3\left[ X^C(\mathcal{C})\hat{r}_i\hat{r}_j+ Y^C(\mathcal{C})\left(\delta_{ij}-\hat{r}_i\hat{r}_j\right)\right]. \end{gather}$$

Here, the coefficients $X(\mathcal {C})$ and $Y(\mathcal {C})$ are scalar resistance functions that encode the configuration-dependent coupling of motion along and transverse to the line of centres $\hat {\boldsymbol{r}}=\boldsymbol{r}/r$ between a particle and the cavity. The superscripts $A$, $B$ and $C$ refer to the force/translation, torque/translation and torque/rotation couplings, respectively. The variable $\mathcal {C}$ defines system configuration, including particle size and particle-to-cavity size ratio. For a single confined particle, $\mathcal {C} = \mathcal {C}(r, \lambda _c)$. For concentrated monodisperse suspensions, $\mathcal {C} = \mathcal {C}(r, \lambda _{c},\phi )$ ($\lambda _p$ does not matter). For concentrated polydisperse suspensions of $K$ different particle sizes, $\mathcal {C} = \mathcal {C}(r, \lambda _{c(1)}, \lambda _{p(2)},\ldots,\lambda _{p(K)},\phi, \phi _i/\phi,\ldots, \phi _K/\phi )$, where $\lambda _{c(1)}\equiv a_1/R$ is the smallest particle-to-cavity ratio, $\lambda _{p(i)}\equiv a_i/a_1>1$ is the particle-to-particle size ratio defined with respect to the smallest particle, and $\phi /\phi _i$ is the ${\textit {i}}{\rm th}$ partial volume fraction. Here, the dimensions of the resistance tensors are captured by the prefactors $6{\rm \pi} \eta a^n$ which also render the hydrodynamic functions dimensionless. Finally, $\delta _{ij}$ is the Kronecker delta, and $\epsilon _{ijk}$ is the Levi-Civita or alternating tensor.

The goal of the present work is to obtain the total suspension stress via the hydrodynamic stresslet; we thus focus specifically on those resistance tensors relating the stresslet to particle motion. The relevant resistance tensors from the GRM are

(2.11)$$\begin{align} R_{ijk}^{SU} &= \left( R_{ijk}^{FE} \right)^{\rm T} =6{\rm \pi}\eta a^2\left[ X^{G}(\mathcal{C})\left( \hat{r}_i\hat{r}_j - \tfrac{1}{3}\delta_{ij} \right)\hat{r}_k+ Y^{G}(\mathcal{C})\left( \hat{r}_i\delta_{jk}+\hat{r}_j\delta{ik}-2\hat{r}_i\hat{r}_j\hat{r}_k \right)\right], \end{align}$$

(2.12)$$\begin{align}R_{ijk}^{S\varOmega} &= \left(R_{ijk}^{LE} \right)^{\rm T} = 6{\rm \pi}\eta a^3\left[ Y^H(\mathcal{C})\left(\hat{r}_i\epsilon_{jkm}\hat{r}_m +\hat{r}_j\epsilon_{ikm}\hat{r}_m\right)\right], \end{align}$$

(2.13) \begin{align} R_{ijkl}^{SE} &= 6{\rm \pi}\eta a^3\left[ X^M(\mathcal{C})\tfrac{3}{2} \left( \hat{r}_k\hat{r}_l - \tfrac{1}{3}\delta_{kl} \right)\right.\nonumber\\ &\quad + \left. Y^M(\mathcal{C})\tfrac{1}{2}\left( \hat{r}_i\delta_{jl}\hat{r}_k+\hat{r}_j\delta_{il}\hat{r}_k+\hat{r}_i\delta_{jk}\hat{r}_l+\hat{r}_j\delta_{ik}\hat{r}_l-4\hat{r}_i\hat{r}_j\hat{r}_k\hat{r}_l \right) \right.\nonumber\\ &\quad + \left. Z^M(\mathcal{C})\tfrac{1}{2}\big(\delta_{ik}\delta_{jl}+\delta_{jk}\delta_{il}-\delta_{ij}\delta_{kl}+\hat{r}_i\hat{r}_j \delta_{kl} +\delta_{ij}\hat{r}_k\hat{r}_l+\hat{r}_i\hat{r}_j\hat{r}_k\hat{r}_l-\hat{r}_i \delta_{jl}\hat{r}_k-\hat{r}_j\delta_{il}\hat{r}_k \right. \nonumber\\ &\quad \vphantom{\left(\hat{r}_i\hat{r}_j - \tfrac{1}{3}\delta_{ij} \right)}- \left. \hat{r}_i\delta_{jk}\hat{r}_l-\hat{r}_j\delta_{ik}\hat{r}_l\big)\vphantom{\Big( \hat{r}_k\hat{r}_l - \tfrac{1}{3}\delta_{kl} \Big)}\right]. \end{align}

The superscripts, $G$, $H$ and $M$ indicate the stresslet/translation, stresslet/rotation and stresslet/straining coupling.

Each of the hydrodynamic functions (2.8)–(2.13) encodes the influence of the specific flow associated with a given surface traction, each being chosen to reveal particle motion and rheology of interest. The function $X^A$ (2.8) corresponds to the hydrodynamic force on the colloid associated with axisymmetric translation along the line of centres. This line of centres is shown for a particle inside a sphere in figure 2(a) as the z-axis. The function $Y^A$ plays the same role for translation transverse to the line of centres, where the axis is chosen here for convenience to be along the x-axis in figure 2(b). We could equivalently choose motion along the y-axis. The stresslet sought here can arise from uniaxial, shear and planar-elongational flows, and the corresponding couplings of superscript $G$, $H$ and $M$ each have longitudinal and transverse components $X^G$, $X^M$, $Y^G$, $Y^H$, $Y^M$ and $Z^M$. Here, the $Z$ coupling corresponds to transverse-planar straining motion, related to uniaxial, shear and planar-elongational flow depicted in figures 2(a)–2(c).

Figure 2. Conceptual sketch of model straining, translational and rotational flows for (a) axisymmetric, (b) transverse and (c) transverse-planar motions of a colloid inside a spherical cavity.

The expressions (2.8)–(2.10) and (2.11)–(2.13) have been previously derived for unbound dilute suspensions (Jeffrey & Onishi Reference Jeffrey and Onishi1984; Jeffrey Reference Jeffrey1992) and that framework subsequently leveraged to derive corresponding hydrodynamic functions for concentrated, unconfined suspensions (Durlofsky et al. Reference Durlofsky, Brady and Bossis1987; Ladd Reference Ladd1990; Sierou & Brady Reference Sierou and Brady2001). In these unconfined conditions, all entries shown in (2.7) for the GMM and GRM in are available. In contrast, the corresponding pair-level functions for spherical confinement, which are necessary for developing the confined, concentrated framework, are incomplete. Specifically, there is no existing theory for hydrodynamic functions that represent the stresslet for any flow, nor functions to represent straining flows, for a particle inside a spherical cavity. These resistance tensors are essential both for the near-field interactions between pairs and for many-body interactions; these are, in turn, required to calculate rheological properties of crowded suspensions by means of the confined Stokesian algorithm. Earlier models of the dynamics of a confined particle truncate the GRM, using only the forces and torques produced from translational and rotational flows (2.8)–(2.10) (Jeffery Reference Jeffery1915; Stimson & Jeffery Reference Stimson and Jeffery1926; Dean & O'Neill Reference Dean and O'Neill1963; O'Neill Reference O'Neill1964; Majumdar Reference Majumdar1967; Cooley & O'neill Reference Cooley and O'neill1968; O'Neill & Majumdar Reference O'Neill and Majumdar1970a; Jones Reference Jones2009), which is sufficient for modelling equilibrium diffusion of a confined suspension (Aponte-Rivera & Zia Reference Aponte-Rivera and Zia2016; Aponte-Rivera et al. Reference Aponte-Rivera, Su and Zia2018). The strategy used in these earlier studies was to solve the Stokes equations for a translating sphere and a rotating sphere, and then calculate the surface traction moments $\boldsymbol{F}^H$ and $\boldsymbol{L}^H$.

For a particle inside a spherical cavity, solution of Stokes’ equations for the pressure and velocity fields follows solutions of Laplace's equation with boundary conditions on both surfaces (Appendix B). Jeffery (Reference Jeffery1912) developed a solution using bispherical coordinates, elegantly describing a pair of spheres either internal or external to each other (Appendix B). Jeffery's approach underlies all subsequent bispherical-coordinates solutions to Stokes equations for a pair of spheres (internal or external to each other) (Jeffery Reference Jeffery1915; Stimson & Jeffery Reference Stimson and Jeffery1926; Dean & O'Neill Reference Dean and O'Neill1963; O'Neill Reference O'Neill1964; Majumdar Reference Majumdar1967; Cooley & O'neill Reference Cooley and O'neill1968; O'Neill & Majumdar Reference O'Neill and Majumdar1970a; Jones Reference Jones2009). A brief overview of the development of this body of solutions is as follows. Jeffery (Reference Jeffery1915) contributed the methodology of solving Stokes’ equations in bispherical coordinates and with that methodology, he contributed the flow and pressure due to axisymmetric rotation of a pair of spheres and, from it, calculated the hydrodynamic torque. Next, Stimson & Jeffery (Reference Stimson and Jeffery1926) solved the problem of axisymmetric translation and the corresponding hydrodynamic force on the sphere. In that work, they simplified the mathematics substantially by representing the Stokes’ stream function in bispherical coordinates.

The problems of a sphere rotating and translating transverse to an enclosing cavity were solved by O'Neill & Majumdar (Reference O'Neill and Majumdar1970a) using a method developed by Dean & O'Neill (Reference Dean and O'Neill1963) for a sphere near a flat wall, which used a separation of variables method. Examining the problem in cylindrical coordinates revealed that the $\theta$ component of the velocity and pressure fields is independent of the two other coordinates ($\rho,z$) and that those two coordinates can be equivalently described in any coordinate system of revolution (e.g. cylindrical, spherical and bispherical coordinates, among others). Here Dean & O'Neill (Reference Dean and O'Neill1963) proposed an ansatz in cylindrical coordinates $\{\rho,z,\theta \}$ where the angular dependence of the pressure and velocity fields is independent of their spatial dependencies in the radial and axial directions. After utilising this ansatz in the equation of motion, the resultant second order differential equations were solved in bispherical coordinates as a special case of the solution to Laplace's equation reported by Jeffery, which is discussed in detail in § 3.3. This separation of variables approach was applied to the problems of a sphere translating and rotating transverse to the enclosing spherical cavity by O'Neill & Majumdar (Reference O'Neill and Majumdar1970a). The authors reported expressions for the force and torque on the colloid for both velocity fields. The methods based on bispherical coordinates also give solutions to Stokes’ equations and the surface traction moments for two external spheres, which preceded and are equivalent to the work based on twin multipole methods (Jeffrey & Onishi Reference Jeffrey and Onishi1984; Kim & Karrila Reference Kim and Karrila1991; Jeffrey Reference Jeffrey1992; Jeffrey, Morris & Brady Reference Jeffrey, Morris and Brady1993). In summary, this work produced the force and torque entries in the GRM for a particle in a spherical cavity. In the present work, we utilise these methods to derive additional entries in the GRM.

In the next section, we present the methods for obtaining solutions to Stokes’ equations in bispherical coordinates for three mutually orthogonal straining flows (shown in figure 2) for a spherical colloid inside a cavity. In addition, we present the new expressions derived here for calculating the stresslet on the colloid, and the corresponding hydrodynamic functions for translational, rotational and straining flows ((2.11), (2.12) and (2.13)).

3.1. Solution to the Stokes equations for the axisymmetric straining flow

Axisymmetric straining flow, also known as uniaxial straining flow, has axial symmetry about the line of centres connecting the particle and the cavity (z-axis in figure 2a), yielding a flow that is independent of the azimuthal angle $\theta$, with strain rate $\boldsymbol{E}$ given by

(3.3)\begin{equation} \boldsymbol{E} =E\left[\begin{array}{@{}ccc@{}} -\tfrac{1}{3} & 0 & 0\\ 0 & -\tfrac{1}{3} & 0\\ 0 & 0 & \tfrac{2}{3} \end{array} \right]_{(x,y,z)}. \end{equation}

This flow produces boundary conditions on the surface of the colloid, $S_a$, and on the cavity, $S_{cav}$ given by

(3.4a)$$\begin{gather} \boldsymbol{v} = \boldsymbol{E}\boldsymbol{\cdot}\boldsymbol{r} = \left(-\tfrac{1}{3}x\boldsymbol{i}_x -\tfrac{1}{3}y\boldsymbol{i}_y + \tfrac{2}{3}z\boldsymbol{i}_z \right)E = \left(- \tfrac{1}{3}\rho\boldsymbol{i}_\rho + \tfrac{2}{3}z\boldsymbol{i}_z\right)E \quad \text{at } S_a, \end{gather}$$

(3.4b)$$\begin{gather}\boldsymbol{v} = \boldsymbol{0} \quad \text{at } S_{cav}. \end{gather}$$

Here, $\boldsymbol{i}$ is a unit normal vector pointing in the direction of the subindex. In (3.4a), the boundary conditions are given in both Cartesian and cylindrical coordinates. The former is useful for algebra in index notation whereas the latter facilitates evaluation of the boundary conditions on coordinate systems of revolution (cylindrical, spherical and bispherical coordinates).

By symmetry, neither $p$ nor $\boldsymbol{v}$ depend on $\theta$; hence, in cylindrical coordinates $\{\rho, z, \theta \}$ (we remark that the traditionally defined cylindrical coordinates $\{ \rho, \theta, z^*\}$ have equivalences with bispherical coordinates $\{\sigma, \xi, \theta \}$ centred between the limiting points $\pm c$, whereas the cylindrical coordinates $\{ \rho, \theta, z\}$ mentioned previously are centred at the particle centre. The axial coordinate transformation is $z = z^*-c\coth (\alpha )$ (see also (3.40)); this relationship makes the differential operators on $z$ and $z^*$ equivalent $\partial / \partial z = \partial / \partial z^*$; (3.5)–(3.9) reflect this connection) the Stokes equations simplify to

(3.5a)$$\begin{gather} \frac{\partial p}{\partial \rho} = \eta\,\biggl( \frac{\partial^2 v_\rho}{\partial \rho^2} + \frac{\partial^2 v_\rho}{\partial z^2} + \frac{1}{\rho}\frac{\partial v_\rho}{\partial \rho} - \frac{v_\rho}{\rho^2} \biggr), \end{gather}$$

(3.5b)$$\begin{gather}\frac{\partial p}{\partial z} = \eta\left( \frac{\partial^2 v_z}{\partial \rho^2} + \frac{\partial ^2 v_z}{\partial z^2} + \frac{1}{\rho}\frac{\partial v_z}{\partial \rho} \right), \end{gather}$$

(3.5c)$$\begin{gather}\frac{1}{\rho}\frac{\partial }{\partial \rho}\left( \rho v_\rho \right) +\frac{\partial v_z}{\partial z} = 0, \end{gather}$$

where $v_\rho$ and $v_z$ are the velocity components in the $\rho$ and $z$ directions. Symmetry in the azimuthal angular direction is also present in the flow resulting from translation of the colloid along the $z$-axis, which was solved by Stimson & Jeffery (Reference Stimson and Jeffery1926) utilising the Stokes flow stream function $\varPsi$. In this work, we adopt a similar approach to solve Stokes’ equations for axisymmetric straining motion.

Following this programme, in terms of $\varPsi$, the velocity components are defined by

(3.6a,b)\begin{equation} v_\rho = \frac{1}{\rho}\frac{\partial \varPsi}{\partial z}, \quad v_z=-\frac{1}{\rho}\frac{\partial \varPsi}{\partial \rho}, \end{equation}

which automatically satisfy continuity. Using $\varPsi$, the equations of motion (3.5a,b) simplify to

(3.7)\begin{equation} E^4[\varPsi] = 0, \end{equation}

where the differential operator $E^4$ is given by

(3.8)\begin{equation} E^4 [{\cdot}]= \left[ E^2\right]^2 [{\cdot}]=\left[ \rho\frac{\partial }{\partial \rho} \left( \frac{1}{\rho}\frac{\partial}{\partial \rho} \right) + \frac{\partial ^2}{\partial z^2} \right]^2[{\cdot}], \end{equation}

and the boundary conditions (3.4) simplify to

(3.9a)$$\begin{gather} \frac{\partial \varPsi}{\partial z}=-\frac{1}{3}\rho^2E \quad \text{and}\quad \frac{\partial \varPsi}{\partial \rho} =-\frac{2}{3}\rho zE \quad \text{at } S_a, \end{gather}$$

(3.9b)$$\begin{gather}\frac{\partial \varPsi}{\partial z}=0\quad \text{and}\quad \frac{\partial \varPsi}{\partial \rho}=0 \quad \text{at } S_{cav}. \end{gather}$$

To find a solution to the operator $E^4$ and evaluate the boundary conditions on the surfaces of the particle and the cavity, we use the equivalence between coordinate systems of revolution (see Appendix C), which allows one to easily transform between cylindrical and bispherical coordinates.

A key advantage of using bispherical coordinates is that if the line of centres is aligned with the $z$-axis of the cylindrical coordinate system (figure 3), one can simultaneously locate two spheres of any relative size and separation (external or internal to each other) with a single set of coordinates. This approach is indispensable for evaluation of boundary conditions and equations of motion for a pair of off-centred spheres. The variable change between cylindrical and bispherical coordinates is

(3.10a,b)\begin{equation} z = \frac{c\sinh(\xi)}{\cosh(\xi)-\sigma}-c\coth(\alpha) , \quad \rho=\frac{c(1-\sigma^2)^{{1}/{2}}}{\cosh(\xi)-\sigma}, \end{equation}

whereas the azimuthal angle $\theta$ is the same in both coordinates and $-\infty <\xi <+\infty$ and $-1<\sigma <+1$, where $\xi$ and $\sigma$ are bispherical coordinates equivalent in spatial representation to the cylindrical coordinates variables $\rho$ and $z$. Here, $c$ is a characteristic scale given by the size and separation of the two spheres (refer to figure 3)

(3.11)\begin{equation} c=\frac{R}{2}\frac{1}{r/R}\left[\left( 1-\lambda _c^2\right)^2 -2\left( 1+\lambda _c^2\right)(r/R)^2+(r/R)^4\right]^{{1}/{2}}. \end{equation}

Rotating curves of either constant $\xi$ or constant $\sigma$ forms bodies of revolution. For example, a sphere is given by surfaces of constant $\xi$ according to

(3.12)\begin{equation} z^2+\rho^2 = c^2\text{csch}^2(\xi). \end{equation}

In this work, $\xi =\alpha$ for the colloid and $\xi =\beta$ for the cavity, such that

(3.13a,b)\begin{equation} a=c\,\text{csch}(\alpha),\quad R=c\,\text{csch}(\beta). \end{equation}

Utilising these definitions, the operator $E^2$ in bispherical coordinates is given by

(3.14)\begin{align} E^2[{\cdot}]&=\frac{\cosh(\xi)-\sigma}{c^2}\left\{ \frac{\partial }{\partial \xi}\left[ \left( \cosh(\xi) -\sigma \right)\frac{\partial}{\partial \xi} \right] \right.\nonumber\\ &\quad + \left. \left(1-\sigma^2 \right)\frac{\partial}{\partial \sigma} \left[ \left( \cosh(\xi) -\sigma \right) \frac{\partial}{\partial \sigma} \right] \right\}[{\cdot}]. \end{align}

Owing to the convenience of working in bispherical coordinates for two-sphere problems, a solution to the operator $E^4$ in bispherical coordinates was reported by Stimson & Jeffery (Reference Stimson and Jeffery1926), for the axisymmetric translation of a sphere along the $z$-axis. This canonical solution is also valid for axisymmetric straining flow, since both flows possess the same azimuthal symmetry. The expression for $\varPsi$ that satisfies (17) is given by

(3.15a)\begin{equation} \varPsi(\sigma,\xi) = \left( \cosh(\xi)-\sigma \right)^{-{3}/{2}}\sum_{n=1}^{\infty}U_n\left( \xi \right)V_n\left(\sigma \right), \end{equation}

where

(3.15b)\begin{align} U_n\left(\xi\right) &=c^3 E \left\{ A_n\cosh\left[(n-\tfrac{1}{2})\xi\right]+B_n\sinh\left[(n-\tfrac{1}{2})\xi\right] +C_n\cosh\left[(n+\tfrac{3}{2})\xi\right]\right\}\nonumber\\ &\quad + \left. D_n\cosh\left[(n+\tfrac{3}{2})\xi\right]\right\} \end{align}

and

(3.15c)\begin{equation} V_n\left( \sigma \right)=P_{n-1}(\sigma)-P_{n+1}(\sigma). \end{equation}

Here, $P_n(\sigma )$ are Legendre polynomials. Determination of the constants $A_n$, $B_n$, $C_n$ and $D_n$ requires evaluation of the four boundary conditions on the particle and the cavity surfaces (3.9). To this end, we use again the equivalence between coordinate systems of revolution (see Appendix C) and change the boundary conditions to bispherical coordinates $\{ \sigma, \xi,\theta \}$ to obtain

(3.16a)\begin{gather} \left.\begin{gathered} \varPsi(\sigma, \xi) =-\frac{a^3}{3}E\left.\left[ \frac{\sinh^2(\xi)}{\cosh(\xi)-\sigma}-\cosh(\xi) \right]\right|_{\xi = \alpha}\\ \text{and}\\ \frac{\partial \varPsi}{\partial \xi}=a^2E\left[ \frac{\sinh^2(\xi)}{\cosh(\xi)-\sigma}-\cosh(\xi) \right] \left.\frac{c\sinh^2(\xi)(1-\sigma^2)}{(\cosh(\xi)-\sigma)^3}\right|_{\xi = \alpha} \quad\text{at } S_a, \end{gathered}\right\} \end{gather}

(3.16b)\begin{gather} \left.\varPsi=0\right|_{\xi = \beta}\quad \text{and}\quad \left.\frac{\partial \varPsi}{\partial \xi }=0 \right|_{\xi = \beta}\quad \text{at } S_{cav}. \end{gather}

Then, utilising Legendre polynomials generating functions such as

(3.17)\begin{equation} \left[\cosh(\xi)-\sigma\right]^{-{1}/{2}}=\sqrt{2}\sum_{n=0}^{\infty}\exp({-(n+1/2)\xi})P_n(\sigma), \end{equation}

and considering the orthogonality of the $V_n(\sigma )$ functions lead to a set of equations that are only a function of $\xi$, evaluated as $\xi =\alpha$ and $\xi = \beta$ at the colloid and the cavity, respectively, which can be manipulated to get closed expressions for the coefficients $A_n$, $B_n$, $C_n$ and $D_n$. Using this approach, we obtained the following explicit expressions:

(3.18a) \begin{align} A_n &= 2 \sqrt{2} n (n+1)\left[(-4 n^3+3 n+1)\exp({2 \alpha (n+3)+2 \beta n})\right. \nonumber\\ &\quad +(4 n^3-7 n+3)\exp({2 (\alpha (n+2)+\beta (n+1))}) \nonumber\\ &\quad +2 (2 n^2+n-1)\exp({3 \alpha + 2 \beta (n+1)})\nonumber\\ & \quad+(-4 n^2+2 n+2)\exp({3 \alpha +2 \beta (n+2)})\nonumber\\ &\quad +(-4 n^2-6 n+4)\exp({\alpha +2 \beta (n+1)}) \nonumber\\ &\quad +(4 n^2-2 n-2)\exp({\beta +2 \alpha (n+3)}) +(4 n^2+6 n-4)\exp({\beta +2 \alpha (n+1)})\nonumber\\ &\quad +(-8 n^2-4 n+2)\exp({\beta +2 \alpha (n+2)}) \nonumber\\ &\quad +(-4 n^3-8 n^2+n+2)\exp({2 (\alpha +\alpha n+\beta n)})\nonumber\\ &\quad +(4 n^3+8 n^2-n-2)\exp({2 (\beta +n (\alpha +\beta ))}) \nonumber\\ &\quad +(-8 n^3-8 n^2+4 n+1)\exp({2 (n+1) (\alpha+\beta )})\nonumber\\ &\quad +(8 n^3+8 n^2-1)\exp({2 \alpha (n+2)+2 \beta n}) \nonumber\\ &\quad +4 (n-1) \exp({3 \alpha +3 \beta +4 \beta n})-4 n \exp({\alpha +3 \beta +4 \beta n}) \nonumber\\ &\quad \left.+\left.\! 2 n (2 n+1) \exp({\alpha +2 \beta (n+2)})\vphantom{(-4 n^3+3 n+1)}\right]\right/ \left[\vphantom{(-8 n^2-8 n+6)}3 (\exp({2 \alpha })-1)(2 n-1) (2 n+1) \right.\nonumber\\ &\quad \times ((-8 n^2-8 n+6)\exp({2 (n+1) (\alpha +\beta )})\nonumber\\ &\quad +(2 n+1)^2 \exp({2 \alpha (n+2)+2 \beta n}) +(2 n+1)^2 \exp({4 \beta +2 n (\alpha +\beta)}) \nonumber\\ &\quad - \left.\vphantom{(-8 n^2-8 n+6)}4 \exp({3 \alpha +\beta +4 \alpha n})-4 \exp({\alpha +3 \beta +4 \beta n}))\right], \end{align}

(3.18b)\begin{align} B_n&=-2 \sqrt{2} n (n+1) \left[(-4 n^3+3 n+1) \exp({2 \alpha (n+3)+2 \beta n})\right.\nonumber\\ &\quad +(4 n^3-7 n+3) \exp({2 (\alpha (n+2)+\beta (n+1))}) \nonumber\\ &\quad -2 (2 n^2+n-1) \exp({3 \alpha +2 \beta (n+1)})+(-4 n^2+2 n+2) \exp({\beta +2 \alpha (n+3)})\nonumber\\ &\quad +(-4 n^2-6 n+4) \exp({\beta +2 \alpha (n+1)}) +(4 n^2-2 n-2) \exp({3 \alpha +2 \beta (n+2)})\nonumber\\ &\quad +(4 n^2+6 n-4) \exp({\alpha +2 \beta (n+1)}) +(8 n^2+4 n-2) \exp({\beta +2 \alpha (n+2)}) \nonumber\\ &\quad +(-4 n^3-8 n^2+n+2) \exp({2 (\alpha +\alpha n+\beta n)})\nonumber\\ &\quad +(4 n^3+8 n^2-n-2) \exp({2 (\beta +n (\alpha +\beta ))}) \nonumber\\ &\quad +(-8 n^3-8 n^2+4 n+1) \exp({2 (n+1) (\alpha +\beta )})\nonumber\\ &\quad +(8 n^3+8 n^2-1) \exp({2 \alpha (n+2)+2 \beta n}) \nonumber\\ &\quad +4 (n-1) \exp({3 \alpha +3 \beta +4 \beta n})-4 n \exp({\alpha +3 \beta +4 \beta n}) \nonumber\\ &\quad -\!\left.\left.\vphantom{(-4 n^3+3 n+1)}2 n (2 n+1) \exp({\alpha +2 \beta (n+2)})\right]\right/ \left[\vphantom{(-4 n^3+3 n+1)}3 (\exp({2 \alpha })-1) (2 n-1) (2 n+1) \right.\nonumber\\ &\quad \times((-8 n^2-8 n+6) \exp({2 (n+1) (\alpha +\beta )})\nonumber\\ &\quad +(2 n+1)^2 \exp({2 \alpha (n+2)+2 \beta n})+(2 n+1)^2 \exp({4 \beta +2 n (\alpha +\beta )}) \nonumber\\ &\quad -\left.\vphantom{(-8 n^2-8 n+6)}4 \exp({3 \alpha +\beta +4 \alpha n})-4 \exp({\alpha +3 \beta +4 \beta n}))\right], \end{align}

(3.19a)\begin{align} C_n&=-2 \sqrt{2} n (n+1) \left[2 (2 n^2+n-3) \exp({\beta +2 \alpha (n+2)})\right.\nonumber\\ &\quad -2 (2 n^2+n-3) \exp({3 \alpha +2 \beta (n+1)})+2 (2 n^2+5 n+2)\exp({\beta +2 \alpha n}) \nonumber\\ &\quad -2 (2 n^2+5 n+2) \exp({\alpha +2\beta n})-2 (4 n^2+6 n+1) \exp({\beta +2 \alpha (n+1)}) \nonumber\\ &\quad +(4 n^2+6 n+2) \exp({3 \alpha +2\beta n})\nonumber\\ &\quad +(-4 n^3-4 n^2+5 n+3) \exp({2 (\alpha (n+3)+\beta (n+1))}) \nonumber\\ &\quad +(4 n^3+4 n^2-5 n-3) \exp({2 (n+2) (\alpha +\beta )})\nonumber\\ &\quad -(4 n^3+12 n^2+5 n-6) \exp({2 (n+1) (\alpha +\beta )}) \nonumber\\ &\quad +(8 n^3+16 n^2+4 n-3) \exp({2 (\alpha (n+2)+\beta (n+1))})\nonumber\\ &\quad -(8 n^3+16 n^2+8 n+1) \exp({2 (\alpha (n+1)+\beta (n+2))}) \nonumber\\ &\quad +(n+2) (2 n+1)^2 \exp({4 \beta +2 n (\alpha +\beta )})\nonumber\\ &\quad +4 (n+1) \exp({3 \alpha +3\beta +4 \beta n})-4 (n+2) \exp({\alpha +3 \beta +4 \beta n})\nonumber\\ &\quad +\!\left.\left.\vphantom{(2 n^2+n-3)}2 n (2 n+3) \exp({\alpha +2 \beta (n+1)})\right]\right/\left[3 (\exp({2 \alpha })-1) (4 n^2+8 n+3)\right.\nonumber\\ &\quad\times ((-8 n^2-8 n+6) \exp({2 (n+1) (\alpha +\beta )})+(2 n+1)^2 \exp({2 \alpha (n+2)+2 \beta n}) \nonumber\\ &\quad +(2 n+1)^2 \exp({4 \beta +2 n (\alpha +\beta )})\nonumber\\ &\quad -\left.\vphantom{(2 n^2+n-3)}4 \exp({3 \alpha +\beta +4 \alpha n}) -4 \exp({\alpha +3 \beta +4 \beta n}))\right], \end{align}

(3.19b) \begin{align} D_n&= 2 \sqrt{2} n (n+1)\left[2 (2 n^2+n-3) \exp({\beta +2 \alpha (n+2)})\right.\nonumber\\ &\quad -2 (2 n^2+n-3) \exp({3 \alpha +2 \beta (n+1)})+2 (2 n^2+5 n+2) \exp({\beta +2 \alpha n}) \nonumber\\ &\quad -2 (2 n^2+5 n+2) \exp({\alpha +2 \beta n})-2 (4 n^2+6 n+1)\exp({\beta +2 \alpha (n+1)}) \nonumber\\ &\quad +(4 n^2+6 n+2) \exp({3 \alpha +2 \beta n})+(-4 n^3-4 n^2+5 n+3) \exp({2 (n+2) (\alpha +\beta )}) \nonumber\\ &\quad +(4 n^3+4 n^2-5 n-3) \exp({2 (\alpha (n+3)+\beta (n+1))})\nonumber\\ &\quad +(4 n^3+12 n^2+5 n-6) \exp({2 (n+1) (\alpha +\beta )}) \nonumber\\ &\quad -(8 n^3+16 n^2+4 n-3) \exp({2 (\alpha (n+2)+\beta (n+1))})\nonumber\\ &\quad +(8 n^3+16 n^2+8 n+1) \exp({2 (\alpha (n+1)+\beta(n+2))}) \nonumber\\ &\quad +(n+2) (2 n+1)^2 (-\exp({4 \beta +2 n (\alpha +\beta )}))-4 (n+1) \exp({3 \alpha +3\beta +4 \beta n})\nonumber\\ &\quad +4 (n+2) \exp({\alpha +3 \beta +4 \beta n})\nonumber\\ &\quad \left.+\left.\! 2 n (2 n+3) \exp({\alpha +2 \beta (n+1)})\vphantom{(2 n^2+n-3)}\right]\right/\left[ \vphantom{(2 n^2+n-3)}3 (\exp({2 \alpha })-1) (2 n+1) (2 n+3)\right.\nonumber\\ &\quad \times ((8 n^2+8 n-6) \exp({2 (n+1) (\alpha +\beta )})+(2 n+1)^2 (-\exp({2 \alpha (n+2)+2 \beta n}))\nonumber\\ &\quad -(2 n+1)^2 \exp({4 \beta +2 n (\alpha +\beta )})+4 \exp({3 \alpha +\beta +4 \alpha n})\nonumber\\ &\quad +\left.\vphantom{(2 n^2+n-3)}4 \exp({\alpha +3 \beta +4 \beta n}))\right]. \end{align}

Insertion of these coefficients into the appropriate expressions above gives the fluid velocity and pressure outside a sphere immersed in an straining flow that is axisymmetric to the line connecting that sphere to the centre of an enclosing sphere (the cavity).

Figure 3. Sketch of the bispherical coordinate system $(\xi,\sigma,\theta )$ about a vertical $z$-axis such that the plane of $z=0$ splits the domain in half. The points $z=\pm c$ are the limiting points of the coordinate system. The coordinate $\theta$ is the same azimuthal angle as in cylindrical or spherical coordinates, $0<\theta <2{\rm \pi}$. The usefulness of the coordinates $(\xi,\sigma )$ stems from the type of bodies of revolution that are formed by rotating about the $z$-axis curves of either constant $\xi$ or $\sigma$. The transformation circle related to $\xi$ can be infinite or point-sized, $-\infty <\xi <+\infty$. Curves of constant $\xi$ are half-circles: a value of $\xi = \pm \infty$ are points of zero radius located at the limiting points $\pm c$ (see curve $\xi _1$), whereas values $\xi < +\infty$ ($\xi > -\infty$) are half-circles centred in the upper (lower) half of the $z$-axis (see curves $\xi _2, \xi _3$ and $-\xi _2$), and values of $\xi \rightarrow 0$ from either $+\infty$ or $-\infty$ limit describe circles of infinite radius whose perimeter lays on the $\pm c$ planes. Thus, rotating curves of constant $\xi$ about the $z$-axis generates spheres. Two spheres resulting from $\xi$ of different signs will be external to each other, i.e. one centred above $+c$ and another centred below $-c$. Relevant for this work is the scenario of two spheres with the same $\xi$ sign, where one sphere is inside the other. The coordinate $\sigma$ ranges between $-1<\sigma <+1$, where $\sigma = 0$ at both limiting points $\pm c$. Curves of constant $\sigma$ are arcs of circles: values $0<\sigma <1$ are greater in length than semicircles (see curve $\sigma _1$), a value $\sigma = 0$ is a circle centred between the limiting points $\pm c$ (see curve $\sigma _2$) and values $-1<\sigma <0$ form are less in length than semicircles (see curve $\sigma _3$). Thus, rotating curves with constant $\sigma$ generate spindle toruses.

3.2. Calculation of the stresslet and hydrodynamic functions for axisymmetric flows

We now calculate the traceless part of the hydrodynamic stresslet on the colloid induced by the axisymmetric straining flow for which we just obtained the velocity and pressure. We start by calculating the integrand of (3.2). First, $\boldsymbol{\sigma }$ and $\boldsymbol{v}$ are written in cylindrical coordinates, $(r, z, \theta )$, with $\theta$ the azimuthal angle. The expression simplifies because $\boldsymbol{\sigma }$ and $\boldsymbol{v}$ are independent of $\theta$. It is further simplified by enforcing the boundary conditions on the integrand (3.4). Next, $\boldsymbol{n}$ and $\boldsymbol{r}$ are written in spherical coordinates $(r, \theta _1, \theta )$, where $\theta _1$ is the polar angle. Finally, all vectorial elements are transformed to Cartesian coordinates to facilitate execution of inner and outer products, resulting in a $3\times 3$ matrix in Cartesian coordinates. The scalar coefficients of each element in the matrix now appear as a combination of cylindrical and spherical coordinates. Finally, the elements and limits of integration are expressed in spherical coordinates. The first integration is taken over the azimuthal angle $\theta \in [0,2{\rm \pi} ]$ to yield

(3.20)\begin{align} \boldsymbol{S}^H&= \left[ \boldsymbol{i}_x\boldsymbol{i}_x+\boldsymbol{i}_y\boldsymbol{i}_y \right] \int_0^{\rm \pi}{\rm \pi} a^3 \left[-p\left( \frac{1}{3} -\cos^2(\theta_1) \right)+\eta\left( 2\frac{\partial v_\rho}{\partial \rho} \sin^2(\theta_1) \right.\right.\nonumber\\ &\quad + \!\left.\left. \left(\frac{\partial v_z}{\partial \rho}+\frac{\partial v_\rho}{\partial z} \right) \cos(\theta_1)\sin(\theta_1) -2\frac{v_\rho}{a}\sin(\theta_1) \right)\right]\sin(\theta_1)\,\text{d}\theta_1 \nonumber\\ &\quad - \boldsymbol{i}_z\boldsymbol{i}_z\int_0^{\rm \pi} 2 {\rm \pi}a^3 \left[-p\left( \frac{1}{3} -\cos^2(\theta_1) \right)-\eta \left(2\frac{\partial v_z }{\partial z }\cos ^2\theta_1 \right.\right.\nonumber\\ &\quad + \!\left.\left. \left( \frac{\partial v_z}{\partial \rho}+ \frac{\partial v_\rho}{\partial z} \right)\cos(\theta_1)\sin(\theta_1)+\frac{v_z}{a}\cos(\theta_1) \right) \right] \sin(\theta_1)\,\text{d}\theta_1 . \end{align}

Owing to the linearity of Stokes’ equations, only diagonal components are present in (3.20), and these obtain the same algebraic relationship as the diagonal elements in the rate of strain imposed, (3.3), i.e. the solution is linear in the boundary data. As expected, by continuity, the trace term in the integrand of (3.20), $\boldsymbol{x}\boldsymbol{\cdot}\boldsymbol{\sigma }\boldsymbol{\cdot}\boldsymbol{n}$, contributes only a pressure component.

To make further progress, the velocity components are rewritten in terms of the stream function (3.6a,b), followed by a three-step simplification procedure: first the partial derivatives are transformed to spherical coordinates; next, we integrate by parts for terms with partial derivatives in $\theta _1$; and third, we enforce the equations of motion (3.5a,b) to eliminate the pressure term. All together leads to a compact integrand that is a function solely of the stream function and its derivatives:

(3.21)\begin{align} \boldsymbol{S}^H &=\int _0^{\rm \pi} {\rm \pi}a^3\left\{ \left[ p\left( \cos^2(\theta_1) -\frac{1}{3}\right)+\eta\frac{\cos(\theta_1)\sin(\theta_1)}{r}E^2\left[ \varPsi \right] \right]\sin(\theta_1)\,\text{d}\theta_1 \right\}\nonumber\\ &\quad\times\left[ \boldsymbol{i}_x\boldsymbol{i}_x+\boldsymbol{i}_y\boldsymbol{i}_y -2\boldsymbol{i}_z\boldsymbol{i}_z\right]. \end{align}

This integral is performed in bispherical coordinates utilising the stream function expression from (3.15a–c), the definition of the operator $E^2$ from (3.14), as well as other coordinate-transformation related expressions in Appendix D, (D1)–(D4). All together these operations yield the final expression for the axisymmetric stresslet,

(3.22)\begin{align} \boldsymbol{S}^H&=-6{\rm \pi}\eta a^2\mathcal{K}_1\mathcal{K}_2\left( \frac{\sqrt{2}}{3}{\rm e}^{-\alpha}\sinh(\alpha) \sum_{n=1}^{\infty}(1+2n)\left\{ \left[n-{\rm e}^{2\alpha}(n-1) \right]\left( A_n+B_n \right) \right.\right.\nonumber\\ &\quad + \!\left.\left. \left[(n+2)- {\rm e}^{2\alpha}(n+1) \right]\left( C_n+D_n \right) \right\} \right)\left(\frac{1}{3}\right)\left[ \boldsymbol{i}_x\boldsymbol{i}_x+\boldsymbol{i}_y\boldsymbol{i}_y-2\boldsymbol{i}_z\boldsymbol{i}_z \right], \end{align}

where the constants $\mathcal {K}_1= a\sinh (\alpha )$ and $\mathcal {K}_2=E$. The stresslet expressed by (3.22) is valid for axisymmetric flows arising not only for straining flows but also for flow arising from particle translation. The difference in the stresslet in each flow emerges from the boundary conditions, which are encoded in their respective expressions of $A_n$, $B_n$, $C_n$ and $D_n$, as well as the values of the constants $\mathcal {K}_1$ and $\mathcal {K}_2$. For axisymmetric translation, $\mathcal {K}_1 = 1$ and $\mathcal {K}_2=U$ and expressions for $(A_n + B_n)$ and $(C_n+D_n)$ are reported in Appendix E.

To retrieve the specific expressions of the coefficients $X^M$ and $X^G$, we must consider the tensor associated with each coefficient, given by (2.11) and (2.13), respectively. We do so, obtaining

(3.23)\begin{align} X^M &= \frac{\sqrt{2}}{3}{\rm e}^{-\alpha}\sinh^2(\alpha)\sum_{n=1}^\infty (1+2n)\left\{ \left[{\rm e}^{2\alpha}(n-1)-n \right]\left( A_n+B_n \right) \right.\nonumber\\ &\quad + \left. \left[ {\rm e}^{2\alpha}(n+1)-(n+2) \right]\left( C_n+D_n \right) \right\}, \end{align}

and

(3.24)\begin{align} X^G &= \frac{\sqrt{2}}{3}{\rm e}^{-\alpha}\sinh(\alpha)\sum_{n=1}^\infty (1+2n)\left\{ \left[{\rm e}^{2\alpha}(n-1)-n \right]\left( A_n+B_n \right) \right.\nonumber\\ &\quad + \left. \left[ {\rm e}^{2\alpha}(n+1)-(n+2) \right]\left( C_n+D_n \right) \right\}. \end{align}

A basic validation for both the solution to Stokes’ equations for axisymmetric straining and the stresslet expression for axisymmetric flows is to calculate the force on the colloid produced by axisymmetric straining where, by symmetry of the GRM (§ 2.2), the hydrodynamic coefficient $X^{\tilde {G}}$ from $\boldsymbol{R}^{FE}$ must be the same as $X^{G}$ from $\boldsymbol{R}^{SU}$ (see (2.11)). The hydrodynamic force on the colloid is given by

(3.25)\begin{equation} \boldsymbol{F}^H=\int_{S_a}\boldsymbol{\sigma}\boldsymbol{\cdot}\boldsymbol{n}\,\text{d}S. \end{equation}

Utilising a similar strategy as the one taken to evaluate the stresslet, we obtain

(3.26)\begin{equation} \boldsymbol{F}^H =-6{\rm \pi}\eta a^2E \left(\frac{\sqrt{2}}{3}\right) \sinh^2(\alpha)\sum_{n=1}^{\infty}\left[ (2n+1)\left(A_n + B_n + C_n+D_n\right) \right]\boldsymbol{i}_z. \end{equation}

As expected, the hydrodynamic force exerted by the fluid on a sphere placed in an axisymmetric straining flow is very similar to the expression reported for the force on a colloid undergoing axisymmetric translation (Stimson & Jeffery Reference Stimson and Jeffery1926; Jones Reference Jones2009). Considering the axisymmetric orthogonal tensor in (2.11) we obtain

(3.27)\begin{equation} X^{\tilde{G}}={\left(\frac{\sqrt{2}}{3}\right) \sinh^2(\alpha)}\sum_{n=1}^{\infty}{\left[ (2n+1)\left(A_n + B_n + C_n+D_n\right) \right]}. \end{equation}

As shown in figure 4, coefficients $X^{\tilde {G}}$ and $X^{G}$ match perfectly over the entire cavity and for all the levels of confinement reported here. This agreement validates our methodology for solving the Stokes equations for axisymmetric straining flow § 3.1, and our methodology to obtain the stresslet expression for axisymmetric flows, i.e. (3.22).

Figure 4. Comparison of hydrodynamic resistance functions for a colloid inside a cavity related to either stresslet (open symbols) or to straining flows (filled symbols) for (a) axisymmetric and (b,c) transverse motions.

For general linear flows, we shall need the complete representation of the stresslet, which requires additional hydrodynamic functions. To this end, next we obtain the functions $Y^G$, $Y^H$, $Y^M$ and $Z^M$ associated with transverse motions to the line of centres.

3.3. General solution method to the Stokes equations in bispherical coordinates and the explicit solution for transverse planar straining flow

At the beginning of this section we put forth a plan to solve the Stokes equations for three flows, to present a new methodology for obtaining such solutions, and to deduce all the hydrodynamic functions associated with the traceless stresslet. Having just obtained the solution for one flow, here we simultaneously present our new generalised method and the solution for the second flow.

We present now our general method for solving Stokes’ equations for a pair of spheres in bispherical coordinates based on a method of separation of variables. The foundations of the method were proposed by Dean & O'Neill (Reference Dean and O'Neill1963) in their solution of Stoke's equations for the problem of transverse rotation of a hard sphere near a flat wall. The same method was later utilised by O'Neill & Majumdar (Reference O'Neill and Majumdar1970a) to obtain the relative translation and rotation transverse to the line of centres between a pair of spheres. To avoid repetition we briefly summarise the technique, emphasising the novel steps in our approach that make it general for both axisymmetric and transverse flows for bodies described by bispherical coordinates.

We start with the ansatz that separates the azimuthal angular dependence ($\theta$) from the spatial dependence ($\rho,z$) on the velocity and pressure fields according to

(3.28a)$$\begin{gather} p(\rho,\theta, z) = \eta \frac{\mathcal{K}_{bc}}{c}P(\rho,z)P_{\theta}(\theta), \end{gather}$$

(3.28b)$$\begin{gather}\boldsymbol{v} =\left[ \begin{array}{@{}c@{}} v_\rho \\ v_\theta \\ v_z \end{array} \right]= \left[ \begin{array}{@{}c@{}} \mathcal{K}_{bc}U(\rho,z)U_\theta(\theta) \\ \mathcal{K}_{bc}V(\rho,z)V_\theta(\theta)\\ \mathcal{K}_{bc}W(\rho,z)W_\theta(\theta) \end{array} \right]_{(\rho,\theta,z)}, \end{gather}$$

where $\mathcal {K}_{bc}$ is a constant with units of velocity (length/time) given by the boundary conditions and the characteristic length $c$ from bispherical coordinates ((3.11) and figure 3). The premise of the ansatz is the amenability of homogeneous and linear partial differential equations (Habberman Reference Habberman2004) to solution by the method of separation of variables. The pressure satisfies Laplace, as does the homogeneous part of the velocity field under the decomposition (Appendix B)

(3.29)\begin{equation} \boldsymbol{v} = \boldsymbol{v}^p+\boldsymbol{v}^h, \end{equation}

where $\boldsymbol{v}^h$ is the homogeneous solution of (3.1a), i.e. it satisfies the Laplacian of a vector field

(3.30)\begin{equation} \nabla^2 \boldsymbol{v}^h =\boldsymbol{0}. \end{equation}

The particular solution $\boldsymbol{v}^p$ is

(3.31)\begin{equation} \boldsymbol{v}^p=\frac{1}{2\eta}\boldsymbol{x}p, \end{equation}

and satisfies

(3.32)\begin{equation} \boldsymbol{\nabla} p = \eta \nabla^2\boldsymbol{v}^p. \end{equation}

The method utilises the particular solution (3.31) and the boundary conditions in the Stokes equations (3.1a) and (3.1b), first, to obtain the functionalities of all the $\theta$-dependent functions $P_\theta, U_\theta, V_\theta,$ and $W_\theta$, which are of the form

(3.33a)$$\begin{gather} P_\theta = U_\theta = W_\theta = \gamma(m\theta), \end{gather}$$

(3.33b)$$\begin{gather}V_\theta = \frac{2\gamma(m\theta)'}{m(m-1)+2}, \end{gather}$$

(3.33c)$$\begin{gather}V_\theta = \frac{\gamma(m\theta)'}{m!}, \end{gather}$$

with $m\in \mathbb {Z}$ and, second, to pose the differential equation for $P(\rho,z)$

(3.34)\begin{equation} L_m^2[P] = 0, \end{equation}

where the operator $L_m^2$ is defined by

(3.35)\begin{equation} L_m^2[{\cdot}] =\left( \frac{\partial^2}{\partial \rho^2}+ \frac {1}{\rho}\frac{\partial}{\partial \rho} -\frac{m^2}{\rho^2}+\frac{\partial^2}{\partial z^2}\right)[{\cdot}]. \end{equation}

By construction, both velocity profiles $\boldsymbol{v}^h$ and $\boldsymbol{v}^p$ hold the same angular dependence; therefore, the Laplacian of $\boldsymbol{v}^h$ (3.30) gives in cylindrical coordinates

(3.36a)$$\begin{gather} r: \quad 0= \mathcal{K}_{bc}U_\theta\left( \frac{1}{\rho}\frac{\partial U^h}{\partial \rho} + \frac{\partial ^2 U^h}{\partial \rho^2} -(m^2 + 1)\frac{U^h}{\rho^2}+ 2m\frac{V^h}{\rho^2} + \frac{\partial ^2 U^h}{\partial z^2}\right), \end{gather}$$

(3.36b)$$\begin{gather}\theta: \quad 0=\mathcal{K}_{bc}V_\theta \left( \frac{1}{\rho}\frac{\partial V^h}{\partial \rho} + \frac{\partial^2 V^h}{\partial \rho^2} - (m^2 + 1)\frac{V^h}{\rho^2} +2m\frac{U^h}{\rho^2} + \frac{\partial^2 V^h}{\partial z^2} \right), \end{gather}$$

(3.36c)$$\begin{gather}z: \quad 0= \mathcal{K}_{bc}W_\theta \left( \frac{1}{\rho}\frac{\partial W^h}{\partial \rho} + \frac{\partial^2 W^h}{\partial \rho^2 } - m^2\frac{W^h}{\rho^2} + \frac{\partial ^2 W^h}{\partial z^2} \right). \end{gather}$$

It is convenient to define the variables

(3.37a–c)\begin{equation} \psi = U^h + V^h, \quad \chi =U^h-V^h \quad \text{and} \quad \varphi = W^h, \end{equation}

where $\psi$ and $\varphi$ naturally arise by adding and subtracting (46a) and (46b). Altogether this yields the set of partial differential equations governing the flow:

(3.38)\begin{equation} L^2_{|m-1|}[\psi] = L^2_{m}[\varphi] = L^2_{m}[P] = L^2_{m+1}[\chi] = 0, \end{equation}

with respect to the same differential operator. The variables in (3.37a–c) define the velocity and pressure profiles according to

(3.39a)$$\begin{gather} p(\rho,z,\theta) = \eta \frac{\mathcal{K}_{bc}}{c}P\gamma(m\theta), \end{gather}$$

(3.39b)$$\begin{gather}v_\rho =\mathcal{K}_{bc} U(\rho,z)U_\theta=\frac{\mathcal{K}_{bc}}{2}\left[\frac{r}{c}P + (\psi+\chi)\right]\gamma(m\theta), \end{gather}$$

(3.39c)$$\begin{gather}v_\theta = \mathcal{K}_{bc}V(\rho,z)V_\theta= \frac{\mathcal{K}_{bc}}{2}\left[ \psi-\chi\right]\left[ \frac{\gamma(m\theta)'}{m!} \right], \end{gather}$$

(3.39d)$$\begin{gather}v_z = \mathcal{K}_{bc}W(\rho,z)W_\theta = \frac{\mathcal{K}_{bc}}{2}\left[\frac{z^*}{c}P+2\varphi \right]\gamma(m\theta). \end{gather}$$

The variable $z^*$ is defined as (the variable $z^*$ corresponds to the cylindrical coordinate system that is centred between the limiting points of bispherical coordinates $\pm c$, see figure 3, whereas the variable $z$ corresponds to the cylindrical coordinate system that is centred at the colloid)

(3.40)\begin{equation} z^* = z+c\coth{\alpha}=\frac{c\sinh(\xi)}{\cosh(\xi)-\sigma}. \end{equation}

There are two clear advantages associated with this technique. The first advantage was developed previously and demonstrated for flows induced by transverse translation and rotational motion (Dean & O'Neill Reference Dean and O'Neill1963; O'Neill Reference O'Neill1964; Majumdar Reference Majumdar1967; Cooley & O'neill Reference Cooley and O'neill1968; O'Neill & Majumdar Reference O'Neill and Majumdar1970a; Jones Reference Jones2009): the problem of solving Stokes’ equations is reduced to finding the solution of only one differential operator $L^2_m$ (3.35).

The second advantage is an outcome of this work: we show here that the technique is general, beyond the two simple flows analysed previously. The method is applicable to all the Stoke's flow problems involving two spheres either external or internal to each other. (In fact, this general solution to the Stokes equations is valid for all geometries that can be described in bispherical coordinates, a particle near a flat wall, and even a single sphere or a spindle.) The generality of the method emerges from its basis in the fundamental properties of the Stokes equations, namely that the pressure is harmonic and, by linearity, the velocity can be split into a particular and a homogeneous solution, where the earlier is related to the pressure and once more the latter is also harmonic. Our recognition of these facts inspired us to pursue a more abstract representation of the solution method that recovers all previous solutions and produces the new solutions we present here. This general solution was actually captured in (3.33)–(3.39), where the differences among flow problems are encapsulated in the value of the constant $m$, the functional form of $\gamma (m\theta )$, and the value of the constant $\mathcal {K}_{bc}$. This is summarised as follows:

(3.41a)$$\begin{gather} \gamma(m\theta) = \sin(m\theta) \quad m=0 \quad \text{for axisymmetric rotation}, \end{gather}$$

(3.41b)$$\begin{gather}\gamma(m\theta) = \cos(m\theta) \quad m=0 \quad \text{for axisymmetric translation and straining}, \end{gather}$$

(3.41c)$$\begin{gather}\gamma(m\theta) = \cos(m\theta) \quad m=1 \quad \text{for transverse rotation, translation and straining}, \end{gather}$$

(3.41d)$$\begin{gather}\gamma(m\theta) = \cos(m\theta) \quad m=2 \quad\text{for transverse planar straining}, \end{gather}$$

and

(3.42a)$$\begin{gather} \mathcal{K}_{bc} = \varOmega c \quad \text{for rotation}, \end{gather}$$

(3.42b)$$\begin{gather}\mathcal{K}_{bc} = \mathcal{U} \quad \text{for translation and} \end{gather}$$

(3.42c)$$\begin{gather}\mathcal{K}_{bc} = E c \quad\text{for straining}. \end{gather}$$

Finally, this method directly produces an explicit expression for the pressure (3.39a). This is particularly advantageous for the solution of axisymmetric translation and straining through this method, which is intractable using the previous method based on the stream function § 3.1. We exploit this advantage to calculate the near-field hydrodynamic functions required for the osmotic pressure, which we report in a separate work.

The general solution to the operator $L^2_m$ can be obtained from the solution to Laplace's equation to a scalar field $f(\xi, \sigma, \theta )$ in bispherical coordinates reported by Jeffery (Reference Jeffery1912). He employed a method of separation of variables similar to that explained above with the azimuthal and spatial dependencies split according to $f(\xi, \sigma, \theta ) = \bar {f}(\xi, \sigma )\cos (m\theta )$. Following this programme, we can rewrite the Laplacian in bispherical coordinates as

(3.43)\begin{equation} \nabla^2 f (\xi, \sigma, \theta) = \cos(k\theta)L^2_m\left[\bar{f}(\xi, \sigma)\right]. \end{equation}

This strategy yields a solution of the form

(3.44)\begin{equation} f(\xi, \sigma, \theta) = \cos(k\theta)\sqrt{\cosh(\xi)-\sigma}\sum_{n=k}^\infty \left\{ g_n^+ \exp({(n+1/2)\xi}) +g_n^-\exp({-(n+1/2)\xi}) \right\}P_n^m(\sigma), \end{equation}

where $g_n^+$ and $g_n^-$ are constants that depend on the boundary conditions, and $P_n^m(\sigma )$ are the associated Legendre polynomials defined by (some authors (Habberman Reference Habberman2004) define the associated Legendre polynomials as the negative of the definition of (3.45), thus care must be taken while utilising recursive relations or definitions related to the $P_n^m(\sigma )$ functions)

(3.45)\begin{equation} P_n^m(\sigma) = (1-\sigma^2)^{m/2}\frac{\mathrm{d}^mP_n(\sigma)}{\mathrm{d}\sigma^m}. \end{equation}

From c the general solution to the operator $L^2_m$ is easily inferred:

(3.46)\begin{align} \bar{f}(\xi, \sigma) = \sqrt{\cosh(\xi)-\sigma}\sum_{n=m}^\infty \left\{ g_n^+\exp({(n+1/2)\xi}) +g_n^-\exp({-(n+1/2)\xi})\right\}P_n^m(\sigma). \end{align}

In the following, we use the generalised separation of variables method to solve the Stokes equations for transverse planar flow, where the boundary conditions are given by

(3.47a)$$\begin{gather} \boldsymbol{v} = \boldsymbol{E}\boldsymbol{\cdot}\boldsymbol{x} = \left(\rho\cos(\theta)\boldsymbol{i}_x -\rho\sin(\theta)\boldsymbol{i}_y \right)E = \left( \rho\cos(2\theta)\boldsymbol{i}_\rho -\rho\sin(2\theta) \boldsymbol{i}_\theta\right)E \quad \text{at }S_a, \end{gather}$$

(3.47b)$$\begin{gather}\boldsymbol{v} = \boldsymbol{0} \quad \text{at } S_{cav}, \end{gather}$$

which reflect the strain rate $\boldsymbol{E}$ characteristic of planar elongation (see figure 2c)

(3.48)\begin{equation} \boldsymbol{E} = E\left[\begin{array}{@{}ccc@{}} 1 & 0 & 0\\ 0 & -1 & 0\\ 0 & 0 & 0 \end{array} \right]_{(x,y,z)}. \end{equation}

Utilising the ansatz of (3.28) alongside the boundary conditions of (3.47) yields a solution to Stokes’ equations of the form

(3.49a)$$\begin{gather} p(\rho,z,\theta) = \eta EP\cos(2\theta), \end{gather}$$

(3.49b)$$\begin{gather}v_\rho= Ec U U_\theta = \frac{Ec}{2}\left[\frac{\rho}{c}P + (\psi+\chi)\right]\cos(2\theta), \end{gather}$$

(3.49c)$$\begin{gather}v_\theta = Ec VV_\theta = \frac{Ec}{2}\left[ (\psi-\chi)\right]\left[ -\sin(2\theta)\right], \end{gather}$$

(3.49d)$$\begin{gather}v_z = Ec W W_\theta = \frac{Ec}{2}\left[\frac{z^*}{c}P+2\varphi \right]\cos(2\theta), \end{gather}$$

where the functions $P$, $\psi$, $\chi$ and $\varphi$ satisfy the differential operator $L_m^2$ according to

(3.50)\begin{equation} L_1^2[\varphi] = L_2^2[\psi]= L_2^2[P]= L_3^2[\chi]=0, \end{equation}

with solutions given by

(3.51a)$$\begin{gather} \psi(\sigma,\xi) = \sqrt{\cosh(\xi)-\sigma}\sum_{n=1}^\infty \left\{ d_n^+\exp({(n+1/2)\xi}) +d_n^-\exp({-(n+1/2)\xi} )\right\}P_n^1(\sigma), \end{gather}$$

(3.51b)$$\begin{gather}P(\sigma,\xi) = \sqrt{\cosh(\xi)-\sigma}\sum_{n=2}^\infty \left\{ b_n^+\exp({(n+1/2)\xi}) +b_n^-\exp({-(n+1/2)\xi}) \right\}P_n^2(\sigma), \end{gather}$$

(3.51c)$$\begin{gather}\varphi(\sigma,\xi) = \sqrt{\cosh(\xi)-\sigma}\sum_{n=2}^\infty \left\{ a_n^+\exp({(n+1/2)\xi}) +a_n^-\exp({-(n+1/2)\xi}) \right\}P_n^2(\sigma), \end{gather}$$

(3.51d)$$\begin{gather}\chi(\sigma,\xi) = \sqrt{\cosh(\xi)-\sigma}\sum_{n=3}^\infty \left\{ f_n^+\exp({(n+1/2)\xi}) +f_n^-\exp({-(n+1/2)\xi})\right\}P_n^3(\sigma). \end{gather}$$

To obtain the values of the $8n$ coefficients $a_n^+$, $a_n^-$, $b_n^+$, $b_n^-$, $d_n^+$, $d_n^-$, $f_n^+$ and $f_n^-$, we utilised the boundary conditions and enforced continuity. Rewriting the boundary conditions (3.47) in terms of the functions $P$, $\psi$, $\chi$ and $\varphi$ gives

(3.52a)$$\begin{gather} P=-\frac{2c}{z^*}\varphi, \quad \chi = \frac{\rho}{z^*}\varphi, \quad \psi = \frac{2\rho}{c} +\frac{\rho}{z^*}\varphi, \quad \text{at } S_a, \end{gather}$$

(3.52b)$$\begin{gather}P=-\frac{2c}{z^*}\varphi, \quad \chi = \frac{\rho}{z^*}\varphi, \quad \psi = \frac{\rho}{z^*} \varphi, \qquad \text{at } S_{cav}, \end{gather}$$

and the same strategy for the equation of continuity (3.1b) yields

(3.53)\begin{equation} \left( 3+ \rho\frac{\partial}{\partial \rho}+z^*\frac{\partial}{\partial z} \right)P + \left( c \frac{\partial}{\partial \rho}-\frac{c}{\rho} \right)\psi +\left( c\frac{\partial}{\partial \rho}+ 3\frac{c}{\rho} \right) \chi + 2c\frac{\partial \varphi}{\partial z} = 0. \end{equation}

Utilising the definitions of $P$, $\psi$, $\chi$ and $\varphi$ (3.51a–d) in the boundary conditions given in 3.52 alongside recurrence formulas of the associated Legendre polynomials gives a set of equations that describes the coefficients $b_n^-$, $d_n^+$, $d_n^-$, $f_n^+$ and $f_n^-$ as a function of instances of the coefficients $a_n^+$ and $a_n^-$:

(3.54) \begin{align} &\sinh(\xi)\left[ b_n^+\exp({(n+1/2)\xi})+b_n^-\exp({-(n+1/2)\xi})\right] \nonumber\\ &\quad =-2\cosh(\xi)\left[ a_n^+\exp({(n+1/2)\xi})+a_n^-\exp({-(n+1/2)\xi})\right] \nonumber\\ & \qquad + 2\left[ a_{n-1}^+\exp({(n-1/2)\xi})+a_{n-1}^-\exp({-(n-1/2)\xi}) \right]\left(\frac{n-2}{2n-1}\right)\nonumber\\ & \qquad + 2\left[ a_{n+1}^+\exp({(n+3/2)\xi})+a_{n+1}^-\exp({-(n+3/2)\xi}) \right]\left(\frac{n+3}{2n+3}\right), \quad {\text{with}}\ \xi \in (\alpha, \beta), \end{align}

(3.55) \begin{align} &\sinh(\xi)\left[ f_n^+\exp({(n+1/2)\xi})+f_n^-\exp({-(n+1/2)\xi}) \right] \nonumber\\ &\quad =+ \left[ a_{n-1}^+\exp({(n-1/2)\xi})+a_{n-1}^-\exp({-(n-1/2)\xi})\right] \left(\frac{1}{2n-1}\right) \nonumber\\ &\qquad -\left[ a_{n+1}^+\exp({(n+3/2)\xi})+a_{n+1}^-\exp({-(n+3/2)\xi})\right] \left(\frac{1}{2n+3}\right)\quad{\text{with }\xi \in (\alpha, \beta)}, \end{align}

(3.56a)\begin{align} &\sinh(\alpha)\left[ d_n^+\exp({(n+1/2)\alpha})+d_n^-\exp({-(n+1/2)\alpha})\right]&\nonumber\\ &\quad = 2\sqrt{2}\sinh(\alpha)\exp({-(n+1/2)\alpha}) \nonumber\\ &\qquad -\left[ a_{n-1}^+\exp({(n-1/2)\alpha})+a_{n-1}^-\exp({-(n-1/2)\alpha})\right] \frac{(n-2)(n-1)}{2n-1} \nonumber\\ &\qquad + \left[ a_{n+1}^+\exp({(n+3/2)\alpha})+a_{n+1}^-\exp({-(n+3/2)\alpha}) \right] \frac{(n+2)(n+3)}{2n+3}, \end{align}

(3.56b)\begin{align} &\sinh(\beta)\left[ d_n^+\exp({(n+1/2)\beta})+d_n^-\exp({-(n+1/2)\beta})\right]\nonumber\\ &\quad =-\left[ a_{n-1}^+\exp({(n-1/2)\beta})+a_{n-1}^-\exp({-(n-1/2)\beta})\right] \frac{(n-2)(n-1)}{2n-1} \nonumber\\ &\qquad + \left[ a_{n+1}^+\exp({(n+3/2)\beta})+a_{n+1}^-\exp({-(n+3/2)\beta})\right] \frac{(n+2)(n+3)}{2n+3}. \end{align}

This strategy also gives two extra equations from the equation of continuity

(3.57) \begin{align} &{\mp}2(n-2)a_{n-1}^{\pm} \pm 2(2n+1)a_{n}^\pm\mp2(n+3)a_{n+1}^\pm- (n-2)b_{n-1}^\pm+5b_n^\pm+ (n+3)b_{n+1}^\pm \nonumber\\ &\quad -d_{n-1}^\pm+2d_n^\pm-d_{n+1}^\pm+ (n-2)(n-3)f_{n-1}^\pm-2(n-2)(n+3)f_n^\pm\nonumber\\ &\quad +(n+3)(n+4)f_{n+1}^\pm=0, \end{align}

where each equation follows either the upper or the lower signs. Together (3.54)–(3.57) form a set of $8n$ equations to get the $8n$ coefficients $a_n^+$, $a_n^-$, $b_n^+$, $b_n^-$, $d_n^+$, $d_n^-$, $f_n^+$ and $f_n^-$, confirming that the problem is fully specified. Utilising the definitions of the coefficients $b_n^-$, $d_n^+$, $d_n^-$, $f_n^+$ and $f_n^-$ in terms of $a_n^+$ and $a_n^-$ in (3.57) generate a set of recursive relations for the coefficients $a_n^+$ and $a_n^-$ given by

(3.58a) $$\begin{gather} M_{nn-1}^{+{+}}a_{n-1}^++ M_{nn}^{+{+}}a_{n}^{+}+M_{nn+1}^{+{+}}a_{n+1}^{+}+ M_{nn-1}^{+{-}}a_{n-1}^{-}+M_{nn}^{+{-}}a_{n}^{-}+ M_{nn+1}^{+{-}}a_{n+1}^-=L_n^+, \end{gather}$$

(3.58b) $$\begin{gather}M_{nn-1}^{-{+}}a_{n-1}^{+}+ M_{nn}^{-{+}}a_{n}^{+}+M_{nn+1}^{-{+}}a_{n+1}^{+}+ M_{nn-1}^{-{-}}a_{n-1}^{-}+M_{nn}^{-{-}}a_{n}^{-}+ M_{nn+1}^{-{-}}a_{n+1}^-=L_n^-, \end{gather}$$

where $n \in \mathbb {Z}^{+}$. The expressions for the $M$-coefficients are reported in Appendix F, whereas the coefficients $L_n^+$ and $L_n^-$ are

(3.59a)\begin{align} L_n^+&=4\sqrt{2}\left[\frac{1}{\exp({(2n-1)\alpha})-\exp({(2n-1)\beta})}+\frac{1}{\exp({(2n+3)\alpha}) -\exp({(2n+3)\beta})}\right.\nonumber\\ &\quad - \left.\frac{2}{\exp({(2n+1)\alpha})-\exp({(2n+1)\beta})}\right], \end{align}

(3.59b)\begin{align} L_n^-&=4\sqrt{2}\left[ \frac{\exp({(2n+3)\beta})}{\exp({(2n+3)\beta})-\exp({(2n+3)\alpha})}+\frac{\exp({\alpha+2n\beta})} {\exp({\alpha+2n\beta})-\exp({\beta+2n\alpha})}\right.\nonumber\\ &\quad + \left. 2\frac{\exp({\beta+2n\beta})}{\exp({\alpha+2n\alpha})-\exp({\beta+2n\beta})} \right], \end{align}

(3.59c)\begin{align} L_n^-&=4\sqrt{2}\left[ \frac{1}{1-\exp({(2n+3)(\alpha-\beta)})}+\frac{1}{1-\exp({(2n-1)(\alpha-\beta)})}\right.\nonumber\\ &\quad - \left. 2\frac{1}{1-\exp({(2n+1)(\alpha-\beta)})}\right], \end{align}

and derive from the first term on the right-hand side of the boundary condition for $\varphi$ in (3.52a). Thus, solving 3.58 gives access to coefficients $a_n^+$ and $a_n^-$ which can then be used in (3.54)–(3.56) to get the coefficients $b_n^+$, $b_n^-$, $d_n^+$, $d_n^-$, $f_n^+$ and $f_n^-$.

In this work, we also obtained the solution to Stokes’ equations for the transverse straining flow associated with a simple shear flow, and the details of the Stokes equations solution following the method of separation of variables is given in Appendix G. We have thus completed solution of the Stokes flow problem for fluid velocity and pressure. Next we apply the results to compute the stresslet from which we will deduce additional hydrodynamic functions.

3.4. Calculation of the traceless stresslet for transverse-planar flows

In this section, we outline the methodology to obtain the hydrodynamic stresslet on the colloid for the transverse-planar straining flow. A similar procedure applies to obtain the stresslet for transverse flows involving translation, rotation or straining motion (Appendix H).

The general strategy to calculate the hydrodynamic stresslet exerted on the colloid by the fluid (3.2) starts by posing the surface of integration, $S_a$, in spherical coordinates centred at the colloid $(r, \theta _1, \theta )$, and then writing the integrand considering the velocity and pressure profiles (3.49) according to the separation of variables ansatz in cylindrical coordinates $(\rho, z, \theta )$, but projected to Cartesian coordinates. Then, the traceless stresslet expression is simplified by applying the boundary conditions at the colloid and doing an integration over the azimuthal angle $\theta \in (0,2{\rm \pi} )$ to yield

(3.60)\begin{align} \boldsymbol{S}^{H,TL} &= \left\{ \frac{\rm \pi}{2}\eta a \left[ \int_0^{\rm \pi} a^2 \sin(\theta_1)\, \mathrm{d}\theta_1 \left\{-P\sin^2(\theta_1)+ 2c\frac{\partial U}{\partial \rho}\sin^2(\theta_1)\right.\right.\right.\nonumber\\ &\quad +c\left(\frac{\partial W}{\partial \rho}+\frac{\partial U}{\partial z}\right)\cos(\theta_1)\sin(\theta_1) + c\left(\frac{\partial V}{\partial \rho} +\frac{2U}{\rho} - \frac{V}{\rho}\right)\sin^2(\theta_1)\nonumber\\ &\quad + \!\left.\left.\left. c\left( \frac{2W}{\rho}+\frac{\partial V}{\partial z}\right)\cos(\theta_1)\sin(\theta_1) \right\} \right]-\frac{8}{3}{\rm \pi}\eta a^3 \right\}E\left[ \boldsymbol{i}_x\boldsymbol{i}_x-\boldsymbol{i}_y\boldsymbol{i}_y \right]. \end{align}

To make progress towards the evaluation of the integral in (3.60), first, redundant terms are eliminated with the use of the equation of continuity written in terms of the variables $U$, $V$ and $W$, according to

(3.61)\begin{equation} {c \left( \frac{U}{\rho} + \frac{\partial U}{\partial \rho}-\frac{2V}{\rho}+\frac{\partial W}{\partial z} \right) = 0}, \end{equation}

and, second, the equivalence between coordinate systems of revolution (see Appendix C) is utilised by executing a variable change from cylindrical coordinates $(\rho, z, \theta )$ to spherical coordinates $(r,\theta _1 , \theta )$, where terms that have partial derivatives with respect to $\theta _1$ are eliminated upon integration. Altogether these steps give the expression for the integral as

(3.62a)\begin{align} & \frac{\rm \pi}{2}\eta a \int_0^{\rm \pi} a^2\sin(\theta_1)\, \mathrm{d}\theta_1\left[-P\sin^2(\theta_1)+\sin(\theta_1)c\left( \frac{\partial U}{\partial r}+\frac{\partial V}{\partial r} +\frac{U}{r}+\frac{V}{r}\right) \right] \end{align}

(3.62b)\begin{align} &\quad = \frac{\rm \pi}{4}\eta a \int_0^{\rm \pi} a^2\sin^2(\theta_1)\, \mathrm{d}\theta_1 \left[ \sin(\theta_1)r\frac{\partial P}{\partial r} + \frac{2c}{r}\frac{\partial \rho\psi}{\partial r} \right] \end{align}

(3.62c)\begin{align} &\quad =- \frac{\rm \pi}{4}\eta a^3 \int_{-1}^{+1} \, \mathrm{d}\sigma \left[ \frac{\sinh^3(\alpha)(1-\sigma^2)}{(\cosh(\alpha)-\sigma)^3}\left(\frac{\partial P}{\partial \xi } \right) \right.\nonumber\\ &\qquad + \left.2 \frac{\sinh^3(\alpha)(1-\sigma^2)^{1/2}}{(\cosh(\alpha)-\sigma)}\frac{\partial}{\partial \xi}\frac{ \psi}{(\cosh(\xi)-\sigma)} \right]. \end{align}

Here, progression from the first to the second expression relies on use of the definitions of $U$ and $V$ (3.49). The entire integral is transformed to bispherical coordinates ($\sigma,\xi,\theta$) to arrive at (3.62c). One could evaluate (3.62c) numerically, but instead one can generate an analytical solution by utilising the definitions of the functions $P$ and $\varPsi$, the recurrence relations of the associated Legendre polynomials, and the orthogonality of $P_n^2(\sigma )$ for the integration in $\sigma$:

(3.63)\begin{equation} \int_{-1}^{+1}P_m^2(\sigma)P_n^2(\sigma)\, \mathrm{d}\sigma = \frac{2(n+2)(n+1)n(n-1)}{2n+1}\delta_{mn}. \end{equation}

The final expression for the traceless stresslet is given by

(3.64) \begin{align} \boldsymbol{S}^{H,TL}&=-6{\rm \pi}\eta a^3 E\biggl[\frac{1}{45} \sqrt2 \sinh(\alpha)\sum_{n=2}^\infty\left\{n(n + 1)\left[\left( 3(n + 2)(n - 1)b_n^++ 5d_n^+\right)\right.\right.\nonumber\\ &\quad - \!\left.\left. 2\exp({-(1 + 2n)\alpha})\left((n + 2)(n - 1)b_n^-+ 5d_n^-\right)\right]\right\}+\frac{4}{9}\biggr]\left[ \boldsymbol{i}_x\boldsymbol{i}_x-\boldsymbol{i}_y\boldsymbol{i}_y \right]. \end{align}

The hydrodynamic coefficient $Z^M$ is inferred from this expression by recalling the tensorial representation given in (2.13), giving

(3.65)\begin{align} Z^M&=\frac{4}{9} + \frac{1}{45} \sqrt2 \sinh(\alpha) \sum_{n=2}^\infty \left\{n(n + 1)\left[\left( 3(n + 2)(n - 1)b_n^++ 5d_n^+\right)\right.\right.\nonumber\\ &\quad - \!\left.\left. 2\exp({-(1 + 2n)\alpha})\left((n + 2)(n - 1)b_n^-+ 5d_n^-\right)\right]\right\}. \end{align}

A similar procedure followed to generate the stresslet from transverse translation, rotation and straining; the associated equations are provided in Appendix H. These transverse flows share the same functional form for the stresslet:

(3.66)\begin{align} \boldsymbol{S}^{H,TL}&=-6{\rm \pi} \eta a^2\left\{ K_0+\frac{1}{45\sqrt{2}}\sinh^2(\alpha)\right.\nonumber\\ &\quad \times\sum_{n=0}^\infty \left[-2\left\{n(n+1)\left[\left(4(2n+1)-7\coth(\alpha)\right)b_n^++10a_n^+\right]\right.\right.\nonumber\\ &\quad+ \left. 5 \left[2n+1-\coth(\alpha)\right]d_n^+\right\} +\left\{n(n+1)\left\{\left[2(2n+1)-\coth(\alpha)\right]b_n^-+10a_n^-\right\}\right.\nonumber\\ &\quad + \!\left.\left.\left. 5\left[2n+1-\coth(\alpha)\right]d_n^-\right\}\exp({-(2n+1)\alpha})\right]K_1 \vphantom{\frac{1}{45\sqrt{2}}}\right\}K_2 \left[ \boldsymbol{i}_x\boldsymbol{i}_z-\boldsymbol{i}_z\boldsymbol{i}_x \right] . \end{align}

Here the value of the constants $K_0$, $K_1$ and $K_2$ encode the flow-appropriate expression and units:

(3.67)\begin{equation} \left.\begin{array}{llll} K_0 = 0 & K_1 = 1 & K_2 = U & \text{for transverse translation,}\\ K_0 = 0 & K_1 = \sinh(\alpha) & K_2 = a\varOmega & \text{for transverse rotation and} \\ K_0 =\dfrac{2}{9} & K_1 =\sinh(\alpha) & K_2 = aE & \text{for transverse straining.} \end{array}\right\} \end{equation}

In (3.66), the constants $a_n^+$, $a_n^-$, $b_n^+$, $b_n^-$, $d_n^+$ and $d_n^-$ are different for each flow profile, with values set by the boundary conditions specific to each problem. The expressions for the coefficients for translation and rotation can be found in the work of O'Neill & Majumdar (Reference O'Neill and Majumdar1970a), with corrections by Jones (Reference Jones2009). We derived in the present work the functional forms for the constants for transverse straining, associated with simple shear motion. We report these results in Appendix G.

The expressions for the hydrodynamic coefficients $Y^G$, $Y^H$ and $Y^M$ related to transverse flows are obtained by mapping the stresslet expression (3.66) in terms of the resistance formulation with the tensorial representation given in (2.11), (2.12) and (2.13) for transverse translation, rotation and straining, respectively. Doing so yields

(3.68)\begin{align} Y^G &= \frac{1}{45\sqrt{2}}\sinh^2(\alpha)\sum_{n=0}^\infty \left[-2\left\{n(n+1)\left[\left(4(2n+1)-7\coth(\alpha)\right)b_n^++10a_n^+\right]\right.\right. \nonumber\\ &\quad +\left. 5\left[2n+1-\coth(\alpha)\right]d_n^+\right\}+ \left\{n(n+1)\left\{\left[2(2n+1)-\coth(\alpha)\right]b_n^-+10a_n^-\right\}\right.\nonumber\\ &\quad + \!\left.\left. 5\left[2n+1-\coth(\alpha)\right]d_n^-\right\}\exp({-(2n+1)\alpha})\right], \end{align}

(3.69)\begin{align} Y^H&= \frac{1}{45\sqrt{2}}\sinh^3(\alpha)\sum_{n=0}^\infty\left[-2\left\{n(n+1)\left[\left(4(2n+1)-7\coth(\alpha) \right)b_n^++10a_n^+\right]\right.\right. \nonumber\\ &\quad + \left. 5\left[2n+1-\coth(\alpha)\right]d_n^+\right\} +\left\{n(n+1)\left\{\left[2(2n+1)-\coth(\alpha)\right]b_n^-+10a_n^-\right\}\right.\nonumber\\ &\quad + \!\left.\left. 5\left[2n+1-\coth(\alpha)\right]d_n^-\right\}\exp({-(2n+1)\alpha})\right], \end{align}

and

(3.70)\begin{align} Y^M&= \frac{2}{9}+\frac{1}{45\sqrt{2}}\sinh^3(\alpha) \sum_{n=0}^\infty\left[-2\left\{n(n+1)\left[\left(4(2n+1)-7\coth(\alpha)\right)b_n^++10a_n^+\right]\right.\right. \nonumber\\ &\quad + \left. 5\left[2n+1-\coth(\alpha)\right]d_n^+\right\} +\left\{n(n+1)\left\{\left[2(2n+1)-\coth(\alpha)\right]b_n^-+10a_n^-\right\}\right.\nonumber\\ &\quad + \!\left.\left. 5\left[2n+1-\coth(\alpha)\right]d_n^-\right\}\exp({-(2n+1)\alpha})\right]. \end{align}

A simple way to verify the solutions to Stokes’ equations for transverse straining (Appendix G) and the stresslet expression for transverse flows (3.66) is to appeal to the symmetry of the resistance matrix (§ 2.2). Using this approach, we leverage the fact that the expressions for the hydrodynamic coefficients $Y^{\tilde {G}}$ and $Y^{\tilde {H}}$ come from transverse straining motion that corresponds to the force and torque on the colloid, respectively. To do so, we first compute the hydrodynamic force and torque resulting from transverse straining motion:

(3.71)\begin{equation} \boldsymbol{F}^H =-6{\rm \pi} \eta a^2 \sinh^2(\alpha)\frac{-1}{3\sqrt{2}}\sum_{n=0}^\infty\left[n(n+1)b_n^++d_n^+ \right]E\,\hat{\boldsymbol{i}}_x, \end{equation}

and

(3.72) \begin{align} \boldsymbol{L}^H &=-6{\rm \pi}\eta a^3\sinh^3(\alpha)\left(\frac{-\sqrt{2}}{9}\right) \sum_{n=0}^\infty\left\{2\left[n(n+1)\left(2a_n^++b_n^+\coth(\alpha)\right)\right.\right.\nonumber\\ &\quad -\!\left.\left.\left(2n+1-\coth(\alpha)\right)d_n^+\right]- \left[n(n+1)\left(2a_n^-+b_n^-\coth(\alpha)\right)\right.\right.\nonumber\\ &\quad -\!\left.\left.\left(2n+1-\coth(\alpha)\right)d_n^-\right] \exp({-(2n+1)\alpha})\right\}E \,\hat{\boldsymbol{i}}_y. \end{align}

These expressions have the same functional form as those reported by O'Neill & Majumdar (Reference O'Neill and Majumdar1970a), but the specific values of the constants $a_n^+$, $a_n^-$, $b_n^+$, $b_n^-$, $d_n^+$ and $d_n^-$ (derived in Appendix H) are set by the boundary conditions of the transverse straining flow. The coefficients $Y^{\tilde {G}}$ and $Y^{\tilde {H}}$, are deduced directly from the expressions for $\boldsymbol{F}^H$ and $\boldsymbol{L}^H$ given the coupling tensors $\boldsymbol{R}^{FE}$ and $\boldsymbol{R}^{LE}$, which are projected transverse to the line of centres between the colloid and the cavity, $\hat {\boldsymbol{r}}$:

(3.73)\begin{equation} Y^{\tilde{G}} = { -\sinh^2(\alpha)\frac{1}{3\sqrt{2}}}\sum_{n=0}^\infty{\left[n(n+1)b_n^++d_n^+ \right]}, \end{equation}

and

(3.74)\begin{align} Y^{\tilde{H}} &=-\sinh^3(\alpha)\frac{\sqrt{2}}{9}\sum_{n=0}^\infty \left\{2\left[n(n+1)\left(2a_n^++b_n^+\coth(\alpha)\right)-\left(2n+1-\coth(\alpha)\right)d_n^+\right]\right. \nonumber\\ &\quad -\left.\left[n(n+1)\left(2a_n^-+b_n^-\coth(\alpha)\right)-\left(2n+1-\coth(\alpha)\right)d_n^-\right] \exp({-(2n+1)\alpha})\right\}. \end{align}

Verification of our generalised approach is shown in figures 4(b) and 4(c) by the equality of the coefficients $Y^G = Y^{\tilde {G}}$ and $Y^H = Y^{\tilde {H}}$ over all positions $r/R$ in the cavity, and for several values of confinement $\lambda _c$.

Finally, the agreement of the hydrodynamic functions complies with the symmetry of the GRM, which is the fundamental validity test for the correctness of the stresslet-related functions reported here, because the symmetry emerges from the minimum dissipation theorem.

Next, we further stress-test the method by recovering previously reported hydrodynamic functions from lubrication theory and far-field expressions for a particle in a spherical cavity, as well as for the problem of a particle near a flat wall.

4.1. Validation via lubrication theory

The lubrication theory for a particle inside a spherical cavity (one sphere internal to another) is equivalent to that of a pair of unbound spheres (two spheres external to each other), where a variable change inverts the curvature of the bigger sphere from convex (external spheres) to concave (an internal sphere). This equivalence was utilised by O'Neill & Majumdar (Reference O'Neill and Majumdar1970b) to obtain the lubrication theory for transverse translation and rotation for a pair of spheres internal and external to each other. In fact, the hydrodynamic coefficients from lubrication theory are only a function of the separation between the surfaces in close approach and the size ratio of the spheres.

We utilise the lubrication equations for the stresslet of a pair of external spheres reported in Jeffrey (Reference Jeffrey1992) to generate those for a particle inside a spherical cavity. To do so, we transform the variables that describe the surface-to-surface distance and the size ratio of the spheres. We start with $\xi$, the surface-to-surface separation between the two external spheres (following notation from Jeffrey Reference Jeffrey1992, pp. 17 and 18),

(4.1)\begin{equation} \xi = {\dfrac{2\lambda_c}{\lambda_c-1}\epsilon }, \end{equation}

where $\epsilon$ is the surface-to-surface separation between a sphere and the cavity wall, normalised by the colloid radius $a$:

(4.2)\begin{equation} \epsilon = \dfrac{R-a-r}{a} = \frac{1-\lambda_{c}-{r}/{R}}{\lambda_{c}}. \end{equation}

The particle size ratio $\lambda$ with respect to the smallest sphere (following notation from Jeffrey Reference Jeffrey1992, p. 17) is transformed as

(4.3)\begin{equation} \lambda = {-\dfrac{1}{\lambda _c}}, \end{equation}

where the negative sign accounts for the fact that when one sphere resides inside the other, they must both lie along either the positive or the negative $z$-axis. In contrast, two external spheres must reside on opposite axes (see figure 3).

Following this programme, the lubrication equations describing the stresslet hydrodynamic functions for a particle inside a cavity are

(4.4)\begin{align} X^G &= \frac{2}{3}\left[ \frac{3}{2 (\lambda_c -1)^2 \epsilon } +\frac{3 \left(\lambda_c ^2-12 \lambda_c -4\right) \log \left[\epsilon\right]}{10 (\lambda_c -1)^3}\right.\nonumber\\ &\quad + \left.\frac{\left(-5 \lambda_c ^4+181 \lambda_c ^3+453 \lambda_c ^2+566 \lambda_c +65\right) \epsilon \log \left[\epsilon\right]}{70 (\lambda_c -1)^4}\right] + O(1), \end{align}

(4.5)\begin{align} X^M &= \frac{10}{9}\left[\frac{3}{5 (\lambda_c -1)^2 \epsilon }+\frac{3 \left(\lambda_c ^2-17 \lambda_c -9\right) \log \left[\epsilon\right]}{25 (\lambda_c -1)^3}\right.\nonumber\\ &\quad + \left.\frac{\left(-5 \lambda_c ^4+272 \lambda_c ^3+831 \lambda_c ^2+1322 \lambda_c +415\right) \epsilon \log \left[\epsilon\right]}{175 (\lambda_c -1)^4}\right] + O(1), \end{align}

(4.6)\begin{align} Y^G &= \frac{2}{3}\left[\frac{\left(4 \lambda_c ^2+\lambda_c +7\right) \log \left[\epsilon\right]}{10 (\lambda_c -1)^3}- \frac{\left(32 \lambda_c ^4+179 \lambda_c ^3+532 \lambda_c ^2+356 \lambda_c +221\right) \epsilon \log \left[\epsilon\right]}{250 (\lambda_c -1)^4} \right]\nonumber\\ &\quad+ O(1), \end{align}

(4.7)\begin{align} Y^H &= \frac{4}{3} \left[\frac{(2 \lambda_c +1) \log \left[\epsilon\right]}{10 (\lambda_c -1)^2} -\frac{\left(16 \lambda_c ^3+61 \lambda_c ^2+180 \lambda_c -2\right) \epsilon \log \left[\epsilon\right]}{250 (\lambda_c -1)^3}\right] + O(1),\end{align}

(4.8)\begin{align} Y^M &= \frac{10}{9}\left[ \frac{6 \left(\lambda_c ^2+\lambda_c +4\right) \log \left[\epsilon\right]}{25 (\lambda_c -1)^3}-\frac{6 \left(8 \lambda_c ^4+67 \lambda_c ^3+294 \lambda_c ^2+394 \lambda_c +197\right) \epsilon \log \left[\epsilon\right]}{625 (\lambda_c -1)^4} \right]\nonumber\\ &\quad + O(1), \end{align}

(4.9)\begin{align} Z^M &= \frac{10}{9}\left[\frac{3 \left(\lambda_c ^2+1\right) \epsilon \log \left[\epsilon\right]}{5 (\lambda_c -1)^4} \right] + O(1). \end{align}

Here, $O(1)$ constants are a function of the confinement ratio $\lambda _c$ and are obtained by matching the first derivative of lubrication equations to that of the complete solution. We compare the lubrication theory to the full solution in figure 5. As expected, the lubrication theory equations (4.4)–(4.9) have the same functional form as the full-solution hydrodynamic functions developed in §§ 3.2 and 3.4 ((3.23), (3.24), (3.65), (3.68), (3.69) and (3.70)) when the particle is close to the wall. The agreement holds for several levels of confinement, over the narrow region where lubrication effects dominate. This test validates our complete solutions over the near-contact part of the domain. As expected, lubrication theory does not predict the hydrodynamic functions away from contact; to validate our solutions in region far from contact, we next compare our work with the far-field approximation.

Figure 5. Comparison of the resistance hydrodynamic functions obtained in this work (solid lines) with lubrication theory (hollow circles) for a colloid inside a cavity in close proximity to the cavity wall. Shaded regions are a guide to the eye to emphasise the zones where both solutions match.

4.2. Validation via far-field approximation

The far-field approximation is a methodology developed in the framework of the Stokesian dynamics algorithm (Durlofsky et al. Reference Durlofsky, Brady and Bossis1987) that allows the calculation of a GRM that captures many-body hydrodynamic interactions between particles and with the wall (Aponte-Rivera et al. Reference Aponte-Rivera, Su and Zia2018). This methodology is accurate for surface-to-surface distances between particles or between particles and the cavity wall large compared with colloid size. The elements of the far-field GMM for a spherically confined particle were reported by Zia and coworkers (Aponte-Rivera & Zia Reference Aponte-Rivera and Zia2016; Gonzalez et al. Reference Gonzalez, Aponte-Rivera and Zia2021) and its inversion yields the far-field approximation to the many-body GRM.

We compare the coefficients for the far-field confined hydrodynamic stresslet via the far-field approximation with our complete many-body solution ((3.23), (3.24), (3.65), (3.68), (3.69) and (3.70)) in figure 6 for several confinement ratios $\lambda _c$. Agreement is excellent for particles located two or more radii from the cavity wall. As expected, the hydrodynamic coefficients from off-diagonal resistance matrices ($X^G$, $Y^G$ and $Y^H$) tend to zero at the centre of the cavity owing to spherical symmetry, whereas the diagonal elements ($X^M$, $Y^M$ and $Z^M$) tend to a single-like particle value (${\sim }O(1)$).

Figure 6. Comparison of the resistance hydrodynamic functions obtained in this work (solid lines) with the far-field resistance matrix (hollow circles) for a colloid inside a cavity.

We conclude that our theoretical model accurately captures the hydrodynamic coefficients for finite particle-to-cavity size ratios. The final limiting regime to test is that of an infinitely large cavity.

4.3. Validation via the solution of a particle near a flat wall

Here we demonstrate that our model correctly predicts the behaviour of a particle near a flat wall. The model system of a particle in a spherical cavity approaches that of a particle near a flat wall when the cavity becomes large, $R/a\to \infty$. This limiting behaviour is easily captured in bispherical coordinates (our coordinate system of choice), where the body of revolution at constant $\xi$ goes from a point to a sphere to a plane as $\xi$ varies from infinity to zero. Therefore, with the theory developed in this work, we can approximate the behaviour of a particle near a flat wall in the limit of very weak confinement, here chosen as $\lambda _c = 10^{-7}$. To facilitate the comparison between our model system and the geometry for a particle near a flat wall, in this section, the position of the confined particle is described with respect to its distance from the wall ($h/a$) rather than the centre-to-centre distance ($x/R$) (see figure 7).

Figure 7. Comparison of the stresslet resistance hydrodynamic functions obtained in this work (solid line) in the limit of a large cavity with solutions for a colloid near a flat wall (circles).

To the best of the authors’ knowledge, there are only a few exact expressions for the stresslet hydrodynamic functions $X^{G,{flat \ wall}}$, $Y^{G,flat \ wall}$, $Y^{H,flat \ wall}$ and $Y^{M,flat \ wall}$ (Cichocki & Jones Reference Cichocki and Jones1998; Swan Reference Swan2010), although Goldman, Cox and Brenner much earlier reported the hydrodynamic drag and torque on a sphere in shear flow near a wall (Goldman, Cox & Brenner Reference Goldman, Cox and Brenner1967). We considered differences in notation and scaling while comparing our theory to these prior results for a particle near a flat wall; the detailed equivalences are reported in Appendix I. Our results are compared with previous approaches in figure 7.

We start the analysis with figure 7(a), the hydrodynamic coefficient coupling axisymmetric translation to the stresslet, $X^G$, for which there are results from the method of Stokes eigenfunction expansions and lubrication theory (Cichocki & Jones Reference Cichocki and Jones1998). We find excellent agreement between our solution (red solid line) and the eigenfunction solution (orange diamonds). Very close to contact, our theory gives excellent prediction of the divergent growth predicted by lubrication theory (blue star). The eigenfunction solution of Cichocki and Jones deviates slightly from ours close to contact because Cichocki and Jones published only the first 32 terms of the infinite sum. In contrast, close to the wall our complete solution matches the lubrication equation indistinguishably.

Next, in figures 7(b) and 7(c), we plot the hydrodynamic coefficients $Y^{G,flat \ wall}$ and $Y^{H,flat \ wall}$ coupling the stresslet to transverse translation and rotation, respectively. We show complete solutions from both methods (the Stokes eigenfunction expansions (Cichocki & Jones Reference Cichocki and Jones1998) and the separation of variables (Swan Reference Swan2010) in black circles), as well as lubrication theory (Cichocki & Jones Reference Cichocki and Jones1998) and far-field approximations (Swan Reference Swan2010) (green diamonds). These complete solutions for a particle near a flat wall all agree with each other and agree with our results: they all exhibit similar divergence near contact, matching lubrication theory. The functions also tend to zero as the distance between the flat wall and the particle increases, in agreement with the far-field approximation.

In figures 7(d)–7(f) we consider couplings for straining flows. The coupling $Y^{M, flat \ wall}$ (related to shear strain rate) is special in that it is actually the only hydrodynamic function for a particle near a flat wall relating the stresslet to a straining flow for which there is a prior complete solution. In figure 7(d), the black filled circles are an incorrect solution obtained by Swan (Reference Swan2010), who actually generated two separate solutions, one for near contact via a complete solution by separation of variables and a far-field solution. Whereas the Swan's complete solution (filled black circles) recovers divergence near the wall, it disagrees with the far-field solution presented in their same work (green diamonds). The error in Swan's theory was in the expression for the stresslet integral, after the integrand is expressed in bispherical coordinates. In this work, we obtained a complete solution by the method outlined in § 3.4 (red curve), for confinement $\lambda _c = 10^{-7}$. Our complete solution shows excellent agreement with Swan's far-field approximation yet also predicts the correct divergence near the wall, providing a complete solution accurate over the entire domain.

In figures 7(e) and 7(f) we plot our complete solutions for the coefficients $X^{M,flat \ wall}$ and $Z^{M,flat \ wall}$. These predict the correct behaviour over the entire domain, recovering the far-field solutions away from the wall divergence near contact.

In this article we have reported the theoretical relations and plots of the flat-wall results predicted by our work. Since our theory recovers and expands prior theory for a particle near a flat wall, it will be useful to make the information available in tabulated form. Thus, our complete solutions of the coefficients $Y^{flat \ wall}$, $X^{M,flat \ wall}$ and $Z^{M,flat \ wall}$ for a particle near a flat wall are reported in figure 7, alongside tabular data for the complete set of stresslet hydrodynamic functions for a flat wall (table 1 in Appendix J).

In a final validation, we compare our theory to the calculation reported by Goldman et al. (Reference Goldman, Cox and Brenner1967) of the force and torque on a colloid located near a flat wall while embedded in a simple shear flow (symbols in figure 8). To recover the force and torque using our approach requires the input of the hydrodynamic functions $Y^{G,flat \ wall}$ and $Y^{H,flat \ wall}$, respectively, as well as other previously reported hydrodynamic functions (see (I3a) and (I3b) in Appendix I). The large-cavity approximation applied to the hydrodynamic functions derived in the present work produce the force and torque on a colloid embedded in a simple shear flow (solid lines) and show excellent agreement with the solution of Goldman, Cox and Brenner.

Figure 8. Comparison of the force and torque on a colloid embedded in a simple shear flow (inset), calculated with theory developed for a particle near a flat wall (Goldman et al. Reference Goldman, Cox and Brenner1967) (filled symbols) and with the theory developed in this work for a particle inside a spherical cavity with a large radius (open symbols and line).

In summary, we validated the stresslet hydrodynamic functions derived in this work by comparison with the three limiting behaviours presented in this section (lubrication theory, far-field approximation and a particle near a flat wall). We also demonstrated recovery of the symmetry relations dictated by the GRM shown in §§ 3.2 and 3.4. In the next section these hydrodynamic functions will be implemented into the CSD algorithm (Aponte-Rivera et al. Reference Aponte-Rivera, Su and Zia2018; Gonzalez et al. Reference Gonzalez, Aponte-Rivera and Zia2021) to generate the many-body contribution to the suspension stress, which gives access to rheological data such as the high-frequency viscosity.

5.1. Many-body stresslet

The CSD algorithm captures the hydrodynamic coupling of particles with each other and with the cavity wall to model the dynamics and rheology of spherically confined concentrated colloidal suspensions (Aponte-Rivera & Zia Reference Aponte-Rivera and Zia2016; Aponte-Rivera et al. Reference Aponte-Rivera, Su and Zia2018; Gonzalez et al. Reference Gonzalez, Aponte-Rivera and Zia2021). These hydrodynamic interactions are represented via the grand resistance tensor $\boldsymbol{\mathcal {R}}$, as discussed in § 2.2, and depend only on system geometry, including particle and confinement shape and size, as well as the relative separations between objects. Analytical expressions can fully capture such interactions for any separation only for very simple systems. In suspensions in which many-body interactions matter and for which computational approaches are required, the challenge of handling both singular near-contact lubrication interactions and far-field many-body couplings requires special methods. The Stokesian dynamics algorithm solves the two problems separately and then combines them into a full solution by superposition (Bossis & Brady Reference Bossis and Brady1984; Durlofsky et al. Reference Durlofsky, Brady and Bossis1987; Aponte-Rivera et al. Reference Aponte-Rivera, Su and Zia2018; Gonzalez et al. Reference Gonzalez, Aponte-Rivera and Zia2021). As a result, the total resistance tensor is built up by superimposing a many-body far-field resistance tensor $\boldsymbol{\mathcal {R}}_{ff}$ and a pair-wise near-field resistance tensor $\boldsymbol{\mathcal {R}}_{nf}$. Detailed explanations of the method for building the grand resistance tensor $\boldsymbol{\mathcal {R}}$ as well as $\boldsymbol{\mathcal {R}}_{ff}$ and $\boldsymbol{\mathcal {R}}_{nf}$ are well-established and can be found elsewhere (Aponte-Rivera & Zia Reference Aponte-Rivera and Zia2016; Aponte-Rivera et al. Reference Aponte-Rivera, Su and Zia2018). The resulting GRM can be used to produce particle motion as well as to relate the motion and surface tractions to rheological quantities (cf. (A9)). Here we focus on obtaining the rheology from the stresslet.

Non-deformable particles embedded in a flow develop hydrodynamic surface tractions as they resist deforming with the fluid streamlines. This stresslet gives rise to suspension stress (Batchelor Reference Batchelor1970). The stresslet developed on the surface of particles that are free of imposed forces and torques but are embedded in a linear flow is

(5.1)\begin{equation} \boldsymbol{S}^H=-\boldsymbol{R}^{SU}\boldsymbol{\cdot} \left( \boldsymbol{U} - \boldsymbol{u}^{\infty} \right) - \boldsymbol{R}^{S\varOmega}\boldsymbol{\cdot} \left( \boldsymbol{\varOmega} - \boldsymbol{\omega}^{\infty} \right) + \boldsymbol{R}^{SE}:\boldsymbol{E}^{\infty} + \boldsymbol{S}^B +\boldsymbol{S}^{P,dis}. \end{equation}

Note that Brownian motion produce disturbance flows and interparticle forces entrain particles to produce dissipative interactions, both inducing hydrodynamic stresslets on particle surfaces, $\boldsymbol{S}^B$ and $\boldsymbol{S}^{P,dis}$ (Chu & Zia Reference Chu and Zia2016, Reference Chu and Zia2017; Huang & Zia Reference Huang and Zia2021). We recall that the hydrodynamic stresslet comprises several parts (Huang & Zia Reference Huang and Zia2021): a component arising directly from externally imposed flow acting on no-slip surfaces, plus a contribution due to Brownian disturbance flows and dissipative particle interactions, $\boldsymbol{S}^{H} \equiv \boldsymbol{S}^{H,ext} + \boldsymbol{S}^B + \boldsymbol{S}^{P,dis}$.

In Stokesian dynamics simulations, the ‘purely hydrodynamic’ part $\boldsymbol{S}^{H,ext}$ is calculated in two parts, by decomposing the GRM into far- and near-field parts to give a near-field stresslet $\boldsymbol{S}_{nf}$ and a many-body far-field contribution $\boldsymbol{S}_{ff}$:

(5.2)\begin{equation} \boldsymbol{S}^{H,ext} =-\boldsymbol{R}^{SU}_{nf}\boldsymbol{\cdot} \left( \boldsymbol{U} - \boldsymbol{u}^{\infty} \right)-\boldsymbol{R}^{S\varOmega}_{nf}\boldsymbol{\cdot} \left( \boldsymbol{\varOmega} - \boldsymbol{\omega}^{\infty} \right) + \boldsymbol{R}^{SE}_{nf}:\boldsymbol{E}^{\infty}+\boldsymbol{S}_{ff}. \end{equation}

The theory developed in the present work (§§ 3.2, 3.4 and 4.1) provides the near-field particle–cavity interactions that were previously missing from the literature. These new functions ((3.23), (3.24), (3.65), (3.68), (3.69) and (3.70)) now complete the near-field resistance tensors in (5.2), allowing us to compute the full hydrodynamic stresslet in a confined suspension. The purely hydrodynamic stresslet $\boldsymbol{S}^{H,ext}$ (5.2) can be used to calculate the high-frequency viscosity, described next.

5.2. High-frequency viscosity for spherically confined suspensions

Consider a quiescent incompressible Newtonian fluid of viscosity $\eta _{1}$ enclosed by a no-slip spherical cavity of volume $V_c$. That fluid can dissipate energy $\mathcal {E}^{(1)}$ viscously under a weak perturbation as given by

(5.3)\begin{equation} \mathcal{E}^{(1)}=\int_{V_c}\boldsymbol{\nabla}\boldsymbol{\cdot}\left( \boldsymbol{\sigma}^{(1)}\boldsymbol{\cdot}\boldsymbol{v}^{(1)} \right) \mathrm{d}V = 2\eta_{1} \boldsymbol{E}^{(1)}:\boldsymbol{E}^{(1)}V_c, \end{equation}

where $\boldsymbol{v}^{(1)}= \boldsymbol{E}^{(1)}\boldsymbol{\cdot}\boldsymbol{r}$ at the enclosure.

Defining the viscosity as the rheological quantity that relates stress $\boldsymbol{\sigma }$ to strain rate $\boldsymbol{E}$, then the constitutive equation of a generalised incompressible Newtonian fluid of viscosity $\eta$ is given by

(5.4)\begin{equation} \boldsymbol{\sigma} =-p\boldsymbol{I} + 2\eta \boldsymbol{E}. \end{equation}

Now, instead suppose that the fluid is a colloidal suspension with bulk viscosity $\eta _2$ (also Newtonian for weak flow), which is enclosed in the same spherical cavity where the same rate of strain $\boldsymbol{E}^{(1)}$ is imposed and maintained at the enclosure. Then it follows that

(5.5)\begin{equation} \frac{\mathcal{E}^{(2)}}{\mathcal{E}^{(1)}} = \frac{\eta_2}{\eta_1}. \end{equation}

Here $\mathcal {E}^{(2)}$ is given by

(5.6)\begin{equation} \mathcal{E}^{(2)} =\int_{V_c}\boldsymbol{\nabla} \boldsymbol{\cdot} \left( \boldsymbol{\sigma}^{(2)}\boldsymbol{\cdot}\boldsymbol{v}^{(2)} \right) \, \mathrm{d}V, \end{equation}

where

(5.7a)\begin{equation} \boldsymbol{v}^{(2)}=\boldsymbol{E}^{(1)}\boldsymbol{\cdot}\boldsymbol{r}, \end{equation}

at the cavity and

(5.7b)\begin{equation} \boldsymbol{v}^{(2)}=\boldsymbol{U}^{(i)}+ \boldsymbol{\varOmega}^{(i)}\times \boldsymbol{r}^{(i)}, \end{equation}

on each sphere centred at $\boldsymbol{r}^{(i)}$. Utilising the divergence theorem, (5.6) can be written as

(5.8)\begin{align} \mathcal{E}^{(2)} &=\int_{S_c}\tilde{\boldsymbol{n}}\boldsymbol{\cdot}\boldsymbol{\sigma}^{(2)}\boldsymbol{\cdot}\boldsymbol{v}^{(1)}\, \mathrm{d}S + \sum_{i}^N\int_{S_i}\tilde{\boldsymbol{n}}\boldsymbol{\cdot}\boldsymbol{\sigma}^{(2)} \boldsymbol{\cdot}\boldsymbol{v}^{(2)}\, \mathrm{d}S \nonumber\\ &=\int_{S_c + \sum S_i}\tilde{\boldsymbol{n}}\boldsymbol{\cdot}\boldsymbol{\sigma}^{(2)}\boldsymbol{\cdot}\boldsymbol{v}^{(1)}\, \mathrm{d}S +\sum_{i}^N\int_{S_i}\tilde{\boldsymbol{n}}\boldsymbol{\cdot}\boldsymbol{\sigma}^{(2)} \boldsymbol{\cdot}\left(\boldsymbol{v}^{(2)}-\boldsymbol{v}^{(1)}\right)\, \mathrm{d}S, \end{align}

where $N$ is the total number of suspended particles and the normal vectors $\tilde {\boldsymbol{n}}$ are directed outwards from the cavity into the fluid, and is antiparallel to the unit-surface normal pointing outwards from a particle. We utilise the Lorentz reciprocal theorem to simplify the first integral and rearrange terms to obtain

(5.9)\begin{align} \mathcal{E}^{part} = \mathcal{E}^{(2)}- \mathcal{E}^{(1)} = \sum_{i}^N\left[-\int_{S_i}\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\sigma}^{(1)}\boldsymbol{\cdot}\boldsymbol{v}^{(2)}\, \mathrm{d}S- \int_{S_i}\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\sigma}^{(2)}\boldsymbol{\cdot} \left(\boldsymbol{v}^{(2)}-\boldsymbol{v}^{(1)}\right)\, \mathrm{d}S\right],\end{align}

where $\mathcal {E}^{part}$ is the extra dissipation arising from the presence of the particles, and $\boldsymbol{n}$ is directed outwards from a particle into the fluid. Utilising the divergence theorem, the first integrand vanishes because there are no body forces imposed on the pure solvent nor is there a rate of strain inside the particles. To evaluate the second set of integrals, we recognise that particles are undergoing solid-body motion (5.7b) and that the velocity profile of the fluid $\boldsymbol{v}^{(1)}$ evaluated at the centre of the particles, $\boldsymbol{r}^{(i)}$, is to leading order a linear flow:

(5.10)\begin{equation} \boldsymbol{v}^{(1)}(\boldsymbol{r}^{(i)}) \approx \boldsymbol{U}^{(1)}(\boldsymbol{r}^{(i)})+\boldsymbol{\varOmega}^{(1)}(\boldsymbol{r}^{(i)})\times\boldsymbol{r}^{(i)}+\boldsymbol{E}^{(1)}\boldsymbol{\cdot}\boldsymbol{r}^{(i)}, \end{equation}

where $\boldsymbol{U}^{1}(\boldsymbol{r}^{(i)}$ and $\boldsymbol{\varOmega }^{1}(\boldsymbol{r}^{(i)}$ are the translational and rotational velocities of particles embedded in a flow of the same strength as ‘case 1’, a fluid-only filled cavity, at the positions where the particles under evaluation are located. This motion is augmented by disturbance flows due to the presence of the particles, the last term in (5.10).

Insertion of the surface velocity into (5.9) yields

(5.11)\begin{equation} \mathcal{E}^{part} = \sum_{i}^N\left[ \left( \boldsymbol{U}^{(1)} - \boldsymbol{U}^{(i)} \right)\boldsymbol{\cdot} \boldsymbol{F}^{H,i} + \left( \boldsymbol{\varOmega}^{(1)} - \boldsymbol{\varOmega}^{(i)} \right)\boldsymbol{\cdot} \boldsymbol{L}^{H,i} + \boldsymbol{E}^{(1)}:\boldsymbol{S}^{H,i} \right], \end{equation}

where $\boldsymbol{F}^{H,i}$, $\boldsymbol{L}^{H,i}$ and $\boldsymbol{S}^{H,i}$ are the hydrodynamic force, torque and stresslet on particle $i$. For a suspension of freely mobile particles, inserting (5.11) into (5.5) gives (eliminating the label ‘1’ on solvent viscosity for simplicity):

(5.12)\begin{equation} \frac{\eta_2}{\eta} = 1 + \frac{{\displaystyle \sum_i^N\boldsymbol{E}:\boldsymbol{S}^{H,i}}}{2\eta V_c\boldsymbol{E}:\boldsymbol{E}}. \end{equation}

This result is general for any given imposed strain rate and the resulting stresslet induced on the confined particles. One need only select the flow conditions to determine the corresponding viscosity.

We now use this result to stress-test the model and also to deduce the contribution to the high-frequency viscosity $\eta '_\infty$ from the confined particles, beyond that due to the no-slip cavity (the extra dissipation due to the presence of the particles) isolating the purely hydrodynamic contribution from that due to Brownian motion. That is, the purely hydrodynamic part of $\eta _2$, denoted $\eta '_\infty$, is obtained from (5.12) because of the hydrodynamic stresslet.

The hydrodynamic stresslet on the surface of particle $i$ in a spherically confined suspension needed in (5.12) is given by

(5.13)\begin{equation} \boldsymbol{S}^{H,i}(\mathcal{C},\boldsymbol{r}) =-6{\rm \pi}\eta a^3\left[ \tfrac{3}{2}X^{M,i}(\mathcal{C})\boldsymbol{d}^{X}(\hat{\boldsymbol{r}})+\tfrac{1}{2}Y^{M,i}(\mathcal{C})\boldsymbol{d}^{Y}(\hat{\boldsymbol{r}}) + \tfrac{1}{2}Z^{M,i}(\mathcal{C})\boldsymbol{d}^{Z}(\hat{\boldsymbol{r}})\right]:\boldsymbol{E}, \end{equation}

where $\mathcal {C} = (\phi,\lambda _c,r_i)$, the fourth-order tensors $\boldsymbol{d}^{X}, \boldsymbol{d}^{Y}$ and $\boldsymbol{d}^{Z}$ are defined as in (2.13), and the hydrodynamic functions $X^{M,i}, Y^{M,i}$ and $Z^{M,i}$ reflect the many-body coupling between particles and the enclosure. Considering $\eta _2 = \eta '_\infty$ and combining (5.12) and (5.13) yields an expression for the high-frequency viscosity of a suspension confined in a spherical cavity

(5.14)\begin{equation} \frac{\eta'_\infty}{\eta} = 1 + \frac{5}{2}\phi\left[ \frac{9}{50} \left( \langle X^M \rangle + 2 \langle Y^M \rangle +2 \langle Z^M\rangle\right)\right], \end{equation}

where

(5.15a)$$\begin{gather} \langle X^M \rangle = \frac{1}{N_{conf}} \frac{1}{N}\sum _j^{N_{conf}} \sum _i^N X^{M,(i,j)}, \end{gather}$$

(5.15b)$$\begin{gather}\langle Y^M \rangle = \frac{1}{N_{conf}}\frac{1}{N} \sum _j^{N_{conf}} \sum _i^N Y^{M,(i,j)}, \end{gather}$$

(5.15c)$$\begin{gather}\langle Z^M \rangle = \frac{1}{N_{conf}}\frac{1}{N} \sum _j^{N_{conf}} \sum _i^N Z^{M,(i,j)}. \end{gather}$$

In these equations, $\langle {\cdot } \rangle$ denotes an ensemble average over $N_{conf}$ particle configurations in the confined domain. The effect of crowding is embedded both the confinement ratio and the volume fraction. The latter, $\phi$, appears explicitly in (5.14) and implicitly in the hydrodynamic functions. The volume-fraction-dependent values of the hydrodynamic functions were obtained via our statistical sampling algorithm, which is explained in Appendix K.

To ascertain the effect of confinement on the high-frequency viscosity (5.14), we compare its value at a range of volume fractions and cavity sizes with that of an unconfined suspension at the same volume fractions in figure 9. Any increase in $\eta '_\infty$ reflects the increase in viscous dissipation generated by having more solid surfaces over which the fluid streamlines must bend. Comparing the red curve for $\lambda _c=0.10$ with the black curve for an unconfined suspension reveals that a confined suspension can be equally as dissipative as an unconfined suspension at a lower volume fraction. A representative particle in the suspension ‘sees’ only other particles in the unbound domain, but when confined, it sees the cavity as well. One might thus expect that adding more and more particles to a confined suspension at a given confinement ratio (increasing the volume fraction), the particle contribution will dominate over that of the cavity, driving the viscosity towards the unbound value at the same volume fraction. However, we find the opposite: as volume fraction increases at fixed confinement ratio (moving along a curve), the viscosity moves farther and farther away from the unconfined value. Tighter confinement increases this effect even further, as shown by the increase in $\eta '_\infty$ as confinement increases at a fixed value of volume fraction. This strongly suggests that direct hydrodynamic interactions with the cavity have a dominant effect on the viscosity.

Figure 9. High-frequency viscosity for confined colloidal suspensions (5.14) as a function of volume fraction compared with the $\eta '_\infty$ values of an unconfined suspension.

To understand the root of the increment in dissipation with confinement in crowded suspensions, we will eliminate from the $\eta '_\infty$ expression the trivial contributions from the pure solvent and the Einstein (single-particle) contribution, and focus only on that due to hydrodynamic interactions given by

(5.16)\begin{equation} \frac{2}{5}\left( \frac{\eta'_\infty}{\eta} - 1 \right)\frac{1}{\phi}= \frac{9}{50} \left( \langle X^M \rangle + 2 \langle Y^M \rangle +2 \langle Z^M\rangle\right). \end{equation}

This quantity approaches unity for small $\phi$ and as $\lambda _c \rightarrow 0$, recovering the Einstein contribution (as shown in Appendix L, figure 15).

As expected, the particle–particle contribution in an unbound suspension increases with volume fraction, and the added particle–cavity contribution drives viscosity up over all volume fractions. We split this total combined contribution into its far-field and a near-field components, as suggested by (5.2). The far-field contribution is initially higher at lower volume fractions, but barely grows with increased crowding. In contrast, the near-field contribution vanishes at low volume fraction but grows markedly with increased volume fraction. This behaviour reveals that near-field interactions drive the pronounced increase in dissipation with increased crowding (and could arise from particle–particle or particle–cavity interactions). Although we expect the same trends in an unconfined suspension, the increment in near-field particle dissipation outweighs the total unconfined value. This behaviour can be understood by recognising that, even at equal volume fraction, the entropic restriction of the cavity demands more configurations with closely spaced particles, which drive stronger hydrodynamic interactions. We next test this idea that confinement-driven changes in particle–particle interactions drive increased dissipation.

We examine the detailed effect of confinement on near-field interactions by splitting it into the particle–particle and particle–cavity contributions. The near-field contribution is dominated by two-body interactions and, for confined suspensions, emerges both from particle–particle and particle–cavity interactions. These two contributions are each plotted in figure 10. The particle–particle interactions contribute the most to the near-field; indeed, nearly $60\,\%$ of the viscosity comes from particles located less than two radii away from the wall, regardless of volume fraction, simply because particles always locate preferentially near the wall to increase configurational entropy in the bulk of the suspension. Nearly $70\,\%$ of particles can be found in that location (see Appendix L, figure 16). In comparison, particle–cavity near-field interactions contribute less at all volume fractions. Notably, the growth of particle–cavity interactions with increased volume fraction is weak, because the annular volume near the wall is finite and fixed. Occupancy of this annulus saturates. One can observe the direct contribution of confinement, the interaction between the cavity and particles, graphically, as the distance between the open crossed orange squares and the open green squares, which has the same value as the open red squares with Xs.

Figure 10. Evaluation of various contributions to $\eta '_\infty$ (5.16), comparing unconfined with confined colloidal suspension at $\lambda _c = 0.2$ as a function of volume fraction. The confined data are split into the far-field and near-field contributions, and the near-field contribution is further split into particle–cavity vs particle–particle interactions.

Individual pairs of particle–cavity interactions are actually two- to four-fold stronger than the interactions between a pair of particles (Appendix L, figure 16), but the latter are far more numerous. As such, strong particle–cavity near-field interactions are stronger but still contribute less than the particle–particle near-field interactions.

Figure 9 shows that tighter confinement increases viscous dissipation. Now, a renormalisation of the volume fraction on the horizontal axis with the random close packing volume fraction (Man et al. Reference Man, Donev, Stillinger, Sullivan, Russel, Heeger, Inati, Torquato and Chaikin2005) collapses all the viscosity data into a single curve (figure 11), a normalisation routine in the evaluation of colloidal suspension rheological data (Solomon & Boger Reference Solomon and Boger1998). Physically, the master curve suggests that the divergence of the viscosity with respect to volume fraction depends on the proximity of the given suspension to $\phi _{rcp}$. Shown here for the first time in spherically confined suspensions, our results agree with past results for unbound suspensions in both experiments (D'Haene Reference D'Haene1992; Man et al. Reference Man, Donev, Stillinger, Sullivan, Russel, Heeger, Inati, Torquato and Chaikin2005; Shewan & Stokes Reference Shewan and Stokes2014) and simulations (Mewis & Wagner Reference Mewis and Wagner2011; Roquier Reference Roquier2016).

Figure 11. High-frequency viscosity for confined colloidal suspensions as a function of volume fraction normalised by the random close packing volume fraction, $\phi _{rcp}$, of the different confinement levels (Man et al. Reference Man, Donev, Stillinger, Sullivan, Russel, Heeger, Inati, Torquato and Chaikin2005). The values of an unconfined suspension are also normalised and showed for reference.

In summary, by deriving and analysing the confined stresslet, we have shown that confinement increases viscous dissipation over all volume fractions, more so as packing increases, by gathering particles near the wall and, in so doing, promoting particle–particle near-field interactions. We find that the vast majority of particles can be found within two radii of the cavity wall and interactions between them (not with the wall) provide nearly $60\,\%$ of the viscous dissipation beyond the solvent and Einstein correction. Near-field interactions between particles dominate this effect, more important than the far field and also exceeding the total unconfined viscosity.

5.3. Effect of variable particle size

We examined the interrelated effects of confinement, packing fraction and size bidispersity, building on our previous work (Gonzalez et al. Reference Gonzalez, Aponte-Rivera and Zia2021). We studied here bidisperse suspensions at moderate volume fractions ($\phi \in [0.35,0.50]$), varying the mixture from the limit of all large particles to the limit of all small particles, as shown in figure 12. In the range $\phi \in [0.35,0.45]$ (figure 12a), the viscosity of bidisperse suspensions decreases smoothly between the two monodisperse limits. But this changes for suspensions close to maximum packing. (The geometrical restriction for the maximum volume fraction that particles can attain inside the cavity is the maximum random close packing, which is a function of the confinement ratio, $\lambda _c$.) For confinement ratio $\lambda _c=0.2$, where maximum packing is $\phi _{max}\sim 0.53$ (Man et al. Reference Man, Donev, Stillinger, Sullivan, Russel, Heeger, Inati, Torquato and Chaikin2005), the viscosity drops markedly when the composition changes from monodisperse large to bidisperse composition. Indeed, even a slight degree of bidispersity fluidises the suspension, as shown in figure 12(b).

Figure 12. High-frequency viscosity of bidisperse suspension with size ratio $\lambda _p = 2$ and confinement of the largest particle of $\lambda _{c_2}=0.2$ as a function of the percentage volume of small particles for several levels of crowding: (a) $\phi =0.35$, 0.40 and 0.45 and (b) $\phi =0.50$, 0.51 and 0.53.

We confirm the validity of the high-frequency viscosity of confined suspensions predicted via stresslet measurements (5.14) by a comparison with the same quantity measured via two-point microrheology (Aponte-Rivera & Zia Reference Aponte-Rivera and Zia2021), an independent alternative measurement of $\eta '_\infty$. We show this comparison in figures 13(a) and 13(b) for two values of confinement. Both data follow the same qualitative trends: $\eta '_\infty$ increases with tighter confinement (figure 13c) and increased packing fraction (figure 13a,b). The quantitative difference between the two measurements of $\eta '_\infty$ stems from the difference in perturbation utilised to access the rheological response: the former is based on changes in the stress field generated by fluid motion (Stresslet) whereas the latter is based changes in particle configuration generated by forces (Stokeslet) (Zia Reference Zia2018). This quantitative difference between the two measurements is also present in the $\eta '_\infty$ data for unconfined suspensions (Appendix L, figure 17). For confined suspensions, the two approaches are complementary: the two-point microrheology approach is more readily accessible through experiments whereas the stresslet approach can be expanded with little effort to predict the viscosity of confined suspensions out of equilibrium.

Figure 13. High-frequency viscosity $\eta '_\infty$ (5.14) for confined suspensions predicted in silico from the suspension stresslet (empty circles) and from two-point microrheology (empty squares) for (a) $\lambda _c=0.10$ and for (b) $\lambda _c=0.15$. The $\eta '_\infty$ from two-point microrheology (c) follows the same qualitative trend as that predicted from the suspension stresslet: $\eta '_\infty$ increases with both confinement and crowding.