1. Introduction
The sedimentation of bodies in viscous fluids is important in many natural settings and industrial processes, from paper manufacturing (Steenberg & Johansson Reference Steenberg and Johansson1958) to blood circulation (Caro Reference Caro2012) to the settling of contaminant particles through oil in internal combustion engines (Guazzelli Reference Guazzelli2006). The scientific study of these viscous sedimentation processes encompasses both analytical and numerical treatments dating back to the work of Stokes (Reference Stokes1851) on the flow past a sphere in an unbounded fluid. A number of exact solutions have since been derived to describe the dynamics of sedimenting bodies of simple shape and symmetric orientation. A generalization known as the Faxén law gives formulae for the force and torque on a sphere placed in an arbitrary background flow (Faxén Reference Faxén1922, Reference Faxén1924). Stimson & Jeffery (Reference Stimson and Jeffery1926) considered two sedimenting spheres of equal density and radius, with one placed directly above the other, and showed that the settling speed is increased by their interaction through the fluid. Exact series solutions for two arbitrarily oriented identical spheres were then derived (Goldman, Cox & Brenner Reference Goldman, Cox and Brenner1966). Single non-spherical bodies were also treated classically by Oberbeck (Reference Oberbeck1876), Edwardes (Reference Edwardes1892), Jeffery (Reference Jeffery1922) and Lamb (Reference Lamb1932), who found the force and torque on a triaxial ellipsoid in a linear flow field in terms of ellipsoidal harmonics. Later, Chwang & Wu (Reference Chwang and Wu1975, Reference Chwang and Wu1976) gave a simpler solution of the same problem using the singularity method, in which fundamental solutions of the Stokes equations are placed internal to the body surface with coefficients selected so as to satisfy the no-slip boundary condition. For particle–wall interactions, an exact solution for a sphere translating and rotating near a plane wall was obtained by Brenner (Reference Brenner1961), O’Neill (Reference O’Neill1964) and Goldman, Cox & Brenner (Reference Goldman, Cox and Brenner1967a ), and in a shear flow by Goldman, Cox & Brenner (Reference Goldman, Cox and Brenner1967b ). The study of other body types generally requires methods of approximation such as exploiting particle slenderness, as in the various slender-body theories (Batchelor Reference Batchelor1970; Cox Reference Cox1970; Tillett Reference Tillett1970; Keller & Rubinow Reference Keller and Rubinow1976; Lighthill Reference Lighthill1976; Johnson Reference Johnson1980; Blake, Tuck & Wakeley Reference Blake, Tuck and Wakeley2010) or weak particle flexibility (Li et al. Reference Li, Manikantan, Saintillan and Spagnolie2013).
A widely employed strategy for incorporating the hydrodynamic effect of a plane wall is the method of reflections, an iterative solution procedure where boundary conditions are alternately enforced on the particle and the wall, leading to asymptotically valid representations of fluid forces and particle velocities. Problems involving the extreme cases of spheres and rods have been more frequently investigated than spheroids of intermediate eccentricity; a prominent exception is work by Wakiya (Reference Wakiya1959) wherein the mobility of a general ellipsoid near a wall is approximated using reflections of Lamb’s general solution (Lamb Reference Lamb1932). Various investigations using slender bodies include that of Russel et al. (Reference Russel, Hinch, Leal and Tieffenbruck1977), who derived an asymptotic expression for the rotation of a slender cylinder sedimenting near a plane wall, ignoring end effects; they observed two types of trajectories which they termed glancing and reversing. A related result using matched asymptotic expansions is due to Katz, Blake & Paveri-Fontana (Reference Katz, Blake and Paveri-Fontana1975), who gave an analytical solution of the mobility problem for a slender rod near a wall. Later, Barta & Liron (Reference Barta and Liron1988) used resistive force theory to find simplified integral equations for the rigid motion of a (possibly non-straight) slender body near a wall, including the special case of a slender prolate body. Yang & Leal (Reference Yang and Leal1983) studied the more general problem of motion near an interface between fluids of two different viscosities, giving asymptotically valid ordinary differential equations describing the trajectories of slender rods.
Many generalizations of the problem of a sphere moving near a wall in a viscous fluid have been investigated. A series of papers (Bossis, Meunier & Sherwood Reference Bossis, Meunier and Sherwood1991; Cichocki et al. Reference Cichocki, Jones, Kutteh and Wajnryb2000; Swan & Brady Reference Swan and Brady2007) addressed the problem of constructing the (positive-definite) grand mobility tensor for many spherical particles above a plane wall. The extensive literature on the subject includes treatments of the electrophoresis of a single charged non-conducting sphere near a wall (Keh & Anderson Reference Keh and Anderson1985), a deformable droplet moving between two parallel plates (Shapira & Haber Reference Shapira and Haber1988) and a colloidal sphere translating perpendicularly between two parallel plane walls (Keh & Wan Reference Keh and Wan2008). Similar efforts have been extended to the study of the trajectories of swimming micro-organisms near surfaces (see Spagnolie & Lauga Reference Spagnolie and Lauga2012 and references provided therein).
Particle dynamics in the presence of surfaces have also been studied numerically. Hsu & Ganatos (Reference Hsu and Ganatos1989, Reference Hsu and Ganatos1994) computed the solution to the resistance problem for a spheroid of arbitrary aspect ratio near a plane wall using a combined single- and double-layer representation of the flow. In addition to confirming the glancing and reversing scenarios found by Russel et al. (Reference Russel, Hinch, Leal and Tieffenbruck1977), these authors studied the sedimentation problem for inclined walls and observed trajectories in which the particle escapes from the wall as well as a stable steady solution where the particle translates parallel to the wall without rotation. Huang, Hu & Joseph (Reference Huang, Hu and Joseph1998) and Swaminathan, Mukundakrishnan & Hu (Reference Swaminathan, Mukundakrishnan and Hu2006) considered the sedimentation of a prolate body in a circular or rectangular cylinder at finite Reynolds number, Pozrikidis (Reference Pozrikidis2007) considered a sphere near the interface of two fluids of varying viscosity and Kutteh (Reference Kutteh2010) treated flows containing several non-spherical particles by modelling irregular particles as rigid ensembles of spheres. Boundary integral methods are commonly used to solve particle–wall interaction problems numerically (see Pozrikidis Reference Pozrikidis1992), but novel numerical methods have also been developed, including the method of regularized Stokeslets with images by Ainley et al. (Reference Ainley, Durkin, Embid, Boindala and Cortez2008).
In this paper we consider the problem of sedimenting prolate and oblate spheroids of arbitrary aspect ratio in a wall-bounded Stokes flow. We present a unified description of the particle dynamics, combining and generalizing the trajectories observed by Russel et al. (Reference Russel, Hinch, Leal and Tieffenbruck1977) and Hsu & Ganatos (Reference Hsu and Ganatos1994), along with several novel trajectory types. Using the method of images, we derive asymptotically valid ordinary differential equations to describe the body dynamics, accurate up to
$O(h^{-4})$
in the translational velocity and
$O(h^{-5})$
in the rotational velocity, where
$h$
is the distance from the particle centroid to the wall. The resulting system is further reduced to yield analytical solutions for the complete particle trajectory in many cases, and the predictions are found to agree very well with full numerical simulations. These numerical simulations are carried out using a novel double-layer boundary integral formulation, a method of stresslet images. We describe various types of trajectories that can arise during sedimentation near a wall, from glancing and reversing to periodic tumbling orbits, as well as the sliding trajectory that can arise if the wall is tilted relative to gravity. We also generalize previously published work by treating arbitrary particle orientations instead of laterally symmetric configurations. This leads to fully three-dimensional dynamics, which may result in periodic wobbling, three-dimensional glancing and three-dimensional reversing trajectories.
The paper is organized as follows. In § 2 we describe the geometry of the problem and the equations of motion and we present the method of stresslet images. In § 3 we conduct a numerical survey of the particle dynamics and present a qualitative overview of the zoology of trajectory types. In § 4 we apply the method of images to reduce the sedimentation problem for arbitrarily oriented spheroids to a two- or three-dimensional system of ordinary differential equations. In § 5 we analyse this system of equations and provide closed-form results describing particle trajectories in many special cases; for example, a simple inequality indicates whether or not a particle of a given shape and initial data will escape from the wall. We conclude with a discussion of applications and possible directions for future work in § 6.
2. Equations of motion and numerical method
Consider a spheroid of uniform density sedimenting through a viscous fluid near an infinite plane wall located along the
$xy$
-plane. Its surface, denoted by
$S^{\ast }$
, is described by
$$\begin{eqnarray}S^{\ast }=\left\{a(x,y,h)^{\text{T}}+a\unicode[STIX]{x1D64D}_{{\it\phi}}\boldsymbol{\cdot }\unicode[STIX]{x1D64D}_{{\it\theta}}\boldsymbol{\cdot }(X,Y,Z)^{\text{T}}:X^{2}+(a/b)^{2}\,Y^{2}+(a/c)^{2}\,Z^{2}=1\right\},\end{eqnarray}$$
where
$a$
is a length,
$a\boldsymbol{x}_{0}=a(x,y,h)$
is the position of the centroid and
$$\begin{eqnarray}\unicode[STIX]{x1D64D}_{{\it\theta}}=\left(\begin{array}{@{}ccc@{}}\cos {\it\theta} & 0 & -\!\sin {\it\theta}\\ 0 & 1 & 0\\ \sin {\it\theta} & 0 & \cos {\it\theta}\end{array}\right),\quad \unicode[STIX]{x1D64D}_{{\it\phi}}=\left(\begin{array}{@{}ccc@{}}\cos {\it\phi} & -\!\sin {\it\phi} & 0\\ \sin {\it\phi} & \cos {\it\phi} & 0\\ 0 & 0 & 1\end{array}\right),\end{eqnarray}$$
are rotation operators (see figure 1), with
${\it\theta}\in (-{\rm\pi}/2,{\rm\pi}/2]$
and
${\it\phi}=[0,2{\rm\pi})$
. The semi-axis lengths satisfy
$a>b=c$
for prolate spheroids and
$a=b>c$
for oblate spheroids. The body eccentricity is given by
$e=\sqrt{1-c^{2}/a^{2}}\in [0,1]$
in both cases, and the vectors
$\boldsymbol{x}=(x,y,h)^{\text{T}}$
and
$\boldsymbol{X}=(X,Y,Z)^{\text{T}}$
in (2.1) are dimensionless; contact with the wall occurs if
$h=(\sin ^{2}{\it\theta}+(1-e^{2})\cos ^{2}{\it\theta})^{1/2}$
. The body is subject to a gravitational force
$\boldsymbol{F}^{\ast }={\rm\Delta}{\it\rho}gV\hat{\boldsymbol{x}}$
, where
${\rm\Delta}{\it\rho}$
is the density difference between the body and the fluid,
$g$
is gravitational acceleration and
$V$
is the body volume; the study of particle sedimentation near an inclined wall is achieved by considering a gravitational force at an angle
${\it\beta}$
relative to
$\hat{\boldsymbol{x}}$
. In response, the body moves with translational velocity
$\boldsymbol{U}^{\ast }$
and rotational velocity
${\it\bf\Omega}^{\ast }$
which depend on the particle position and orientation. The system is made dimensionless by scaling lengths upon
$a$
and defining the dimensionless translational and rotational velocities
$\boldsymbol{U}=(6{\rm\pi}{\it\mu}a)({\rm\Delta}{\it\rho}gV)^{-1}\boldsymbol{U}^{\ast }$
and
${\it\bf\Omega}=(6{\rm\pi}{\it\mu}a^{2})({\rm\Delta}{\it\rho}gV)^{-1}{\it\bf\Omega}^{\ast }$
. In the case of a spherical particle in an unbounded fluid this results in a dimensionless body of radius and sedimentation speed both equal to unity, the well-known Stokes drag law (Stokes Reference Stokes1851). The surface
$S$
denotes the dimensionless body surface (the surface
$S^{\ast }$
scaled by the length
$a$
).
Figure 1. Schematic of a prolate spheroid near an infinite plane wall located along the
$xy$
-plane. The body, with centroid a distance
$h$
from the wall, is rotated through an angle
${\it\phi}$
about the
$\hat{\boldsymbol{z}}$
axis and is pitched at an angle
${\it\theta}$
about its lateral axis. The dimensionless force due to gravity acts at an angle
${\it\beta}$
relative to the wall,
$\boldsymbol{F}=\cos {\it\beta}\hat{\boldsymbol{x}}-\sin {\it\beta}\hat{\boldsymbol{z}}$
, with
${\it\beta}$
set to zero in the illustration above.
In the theoretical limit of zero-Reynolds-number flow, appropriate for modelling the flow generated by small or slow moving particles, or larger particles in highly viscous fluids, the equations of fluid motion are the Stokes equations,
$$\begin{eqnarray}-\boldsymbol{{\rm\nabla}}p^{\ast }+{\it\mu}{\rm\nabla}^{2}\boldsymbol{u}^{\ast }=\mathbf{0},\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}^{\ast }=0,\end{eqnarray}$$
where
$p^{\ast }$
is the pressure,
$\boldsymbol{u}^{\ast }$
is the fluid velocity and
${\it\mu}$
is the viscosity (see Batchelor Reference Batchelor2000). The dimensionless fluid velocity is also scaled as above,
$\boldsymbol{u}=(6{\rm\pi}{\it\mu}a)({\rm\Delta}{\it\rho}gV)^{-1}\boldsymbol{u}^{\ast }$
, and is assumed to satisfy no-slip boundary conditions on the particle surface,
$\boldsymbol{u}(\boldsymbol{x}\in S)=\boldsymbol{U}+{\it\bf\Omega}\times (\boldsymbol{x}-\boldsymbol{x}_{0})$
, and on the wall,
$\boldsymbol{u}(z=0)=\mathbf{0}$
. The fluid velocity is assumed to decay to zero as the distance from the body tends towards infinity. In the inertialess limit, the integrated fluid stress on the particle surface must balance the external force due to gravity, and there must be zero net fluid torque on the body. These six conditions, along with (2.3), close the system of equations for the fluid velocity and pressure and the particle’s instantaneous translational and rotational velocities,
$\boldsymbol{U}$
and
${\it\bf\Omega}$
. The position and orientation of the body at any time are described by the state variable
${\it\bf\Phi}=(x,y,h,{\it\theta},{\it\phi})$
. Since time does not appear explicitly in the Stokes equations, any solution of the mobility problem for general
${\it\bf\Phi}$
yields an autonomous system of ordinary differential equations,
$\dot{{\it\bf\Phi}}=\boldsymbol{{\mathcal{F}}}(h,{\it\theta},{\it\phi})$
, describing the trajectory of the particle sedimenting under the influence of a constant gravitational force.
2.1. Fundamental singularities and image systems
The linearity of the Stokes equations opens the door to numerous analytical and numerical approaches to solving fluid–body interaction problems that rely on fundamental singularities or Green’s functions. In one particularly useful and clean approach to solving such problems, the singularities are placed internal to an immersed body and their strengths are chosen so as to match the boundary conditions on the surface (Chwang & Wu Reference Chwang and Wu1975). For instance, the Stokes flow around a no-slip spherical boundary in an unbounded flow may be represented as a linear combination of a Stokeslet singularity,
$$\begin{eqnarray}\unicode[STIX]{x1D642}(\boldsymbol{x},\boldsymbol{x}_{0})=\frac{\unicode[STIX]{x1D644}}{|\boldsymbol{x}-\boldsymbol{x}_{0}|}+\frac{(\boldsymbol{x}-\boldsymbol{x}_{0})(\boldsymbol{x}-\boldsymbol{x}_{0})^{\text{T}}}{|\boldsymbol{x}-\boldsymbol{x}_{0}|^{3}},\end{eqnarray}$$
with
$\unicode[STIX]{x1D644}$
the identity operator, and a potential source dipole (see Kim & Karrila Reference Kim and Karrila1991).
The effect of a wall on the trajectory of a moving body can be studied using image singularity systems (Blake Reference Blake1971; Blake & Chwang Reference Blake and Chwang1974). Image systems for Stokeslets of varying orientation relative to a no-slip wall were presented by Blake (Reference Blake1971). As an example, consider an
$x$
-directed Stokeslet
$\boldsymbol{G}_{x}(\boldsymbol{x},\boldsymbol{x}_{0})=\unicode[STIX]{x1D642}(\boldsymbol{x},\boldsymbol{x}_{0})\boldsymbol{\cdot }\hat{\boldsymbol{x}}$
located in the fluid at a point
$\boldsymbol{x}_{0}=(0,0,h)$
. The image system cancels the fluid velocity on the surface
$z=0$
when placed at the image point
$\boldsymbol{x}^{\ast }=(0,0,-h)$
, and is given by
$$\begin{eqnarray}\boldsymbol{G}_{x}^{\ast }(\boldsymbol{x},\boldsymbol{x}^{\ast })=-\boldsymbol{G}_{x}(\boldsymbol{x},\boldsymbol{x}^{\ast })-2h\frac{\partial }{\partial x}\boldsymbol{G}_{z}(\boldsymbol{x},\boldsymbol{x}^{\ast })+2h^{2}\frac{\partial }{\partial x}\boldsymbol{U}_{\!P}(\boldsymbol{x},\boldsymbol{x}^{\ast }),\end{eqnarray}$$
where
$\boldsymbol{G}_{z}=\unicode[STIX]{x1D642}\boldsymbol{\cdot }\hat{\boldsymbol{z}}$
is a
$z$
-directed Stokeslet and
$$\begin{eqnarray}\boldsymbol{U}_{\!P}(\boldsymbol{x},\boldsymbol{x}_{0})=\frac{\boldsymbol{x}-\boldsymbol{x}_{0}}{|\boldsymbol{x}-\boldsymbol{x}_{0}|^{3}}\end{eqnarray}$$
is a potential flow point source. Similarly,
$$\begin{eqnarray}\boldsymbol{G}_{z}^{\ast }(\boldsymbol{x},\boldsymbol{x}^{\ast })=-\boldsymbol{G}_{z}(\boldsymbol{x},\boldsymbol{x}^{\ast })+2h\frac{\partial }{\partial z}\boldsymbol{G}_{z}(\boldsymbol{x},\boldsymbol{x}^{\ast })-2h^{2}\frac{\partial }{\partial z}\boldsymbol{U}_{\!P}(\boldsymbol{x},\boldsymbol{x}^{\ast }).\end{eqnarray}$$
Image systems for derivatives of the Stokeslet may be determined by careful manipulation of the image systems of Blake (Reference Blake1971), though some care must be taken as the image system of the derivative is not in general the derivative of the image system (see Kim & Karrila Reference Kim and Karrila1991 and Spagnolie & Lauga Reference Spagnolie and Lauga2012). Note for instance that the coefficients in (2.5) are
$h$
-dependent. Since
$\boldsymbol{{\rm\nabla}}_{\boldsymbol{x}}\unicode[STIX]{x1D642}(\boldsymbol{x},\boldsymbol{x}_{0})=-\boldsymbol{{\rm\nabla}}_{\boldsymbol{x}_{0}}\unicode[STIX]{x1D642}(\boldsymbol{x},\boldsymbol{x}_{0})$
, we may use Blake’s result to construct the image system for a difference quotient. As an example, the image system of
$\partial ^{2}\,\boldsymbol{G}_{x}/\partial x\partial z$
is given by
$$\begin{eqnarray}\displaystyle \left[\frac{\partial ^{2}}{\partial x\partial z}\boldsymbol{G}_{x}\right]^{\ast }(\boldsymbol{x},\boldsymbol{x}^{\ast }) & = & \displaystyle -\frac{\partial ^{2}}{\partial x\partial z}\boldsymbol{G}_{x}(\boldsymbol{x},\boldsymbol{x}^{\ast })-2h\frac{\partial ^{3}}{\partial x^{2}\partial z}\boldsymbol{G}_{z}(\boldsymbol{x},\boldsymbol{x}^{\ast })+2\frac{\partial ^{2}}{\partial x^{2}}\boldsymbol{G}_{z}(\boldsymbol{x},\boldsymbol{x}^{\ast })\nonumber\\ \displaystyle & & \displaystyle +\,2h^{2}\frac{\partial ^{3}}{\partial x^{2}\partial z}\boldsymbol{U}_{\!P}(\boldsymbol{x},\boldsymbol{x}^{\ast })-4h\frac{\partial ^{2}}{\partial x^{2}}\boldsymbol{U}_{\!P}(\boldsymbol{x},\boldsymbol{x}^{\ast }).\end{eqnarray}$$
2.2. Numerical method: the method of stresslet images
The fundamental singularities of Stokes flow may be used to derive a representation formula for the fluid flow in terms of singular boundary integrals (see Power & Miranda Reference Power and Miranda1987 and Pozrikidis Reference Pozrikidis1992). The presence of a nearby wall has been incorporated into various forms of the boundary integral formulation, for instance using regularized Stokeslets with their images (Ainley et al. Reference Ainley, Durkin, Embid, Boindala and Cortez2008). A well-posed double-layer form of the boundary integral formulation may be adapted for use near an infinite wall using image singularities of the stresslet, as suggested by Spagnolie & Lauga (Reference Spagnolie and Lauga2012). In this double-layer formulation with stresslet images, the fluid velocity is given by (see Pozrikidis Reference Pozrikidis1992)
$$\begin{eqnarray}\boldsymbol{u}(\boldsymbol{x})=-\int _{S}\boldsymbol{q}(\boldsymbol{y})\boldsymbol{\cdot }(\unicode[STIX]{x1D61B}(\boldsymbol{x},\boldsymbol{y})+\unicode[STIX]{x1D61B}^{\ast }(\boldsymbol{x},\boldsymbol{y}^{\ast }))\boldsymbol{\cdot }\hat{\boldsymbol{n}}(\boldsymbol{y})\,\text{d}S+\frac{1}{8{\rm\pi}}\left(\unicode[STIX]{x1D642}(\boldsymbol{x},\boldsymbol{x}_{0})+\unicode[STIX]{x1D642}^{\ast }(\boldsymbol{x},\boldsymbol{x}^{\ast })\right)\boldsymbol{\cdot }\boldsymbol{F},\end{eqnarray}$$
where
$\hat{\boldsymbol{n}}$
is the unit normal vector pointing into the fluid,
$\boldsymbol{y}$
is an integration variable over the body surface,
$\boldsymbol{q}(\boldsymbol{y})$
is an unknown density,
$$\begin{eqnarray}\unicode[STIX]{x1D61B}(\boldsymbol{x},\boldsymbol{y})=-6\frac{(\boldsymbol{x}-\boldsymbol{y})(\boldsymbol{x}-\boldsymbol{y})(\boldsymbol{x}-\boldsymbol{y})}{|\boldsymbol{x}-\boldsymbol{y}|^{5}}\end{eqnarray}$$
is the stresslet singularity, a third-order tensor, and
$\unicode[STIX]{x1D61B}^{\ast }(\boldsymbol{x},\boldsymbol{y}^{\ast })$
is the associated image system, which is singular at the image point
$\boldsymbol{y}^{\ast }$
inside the wall and is given by the formula
$$\begin{eqnarray}\unicode[STIX]{x1D61B}_{ijk}^{\ast }(\boldsymbol{x},\boldsymbol{y}^{\ast })=\frac{6\hat{X}_{i}X_{j}X_{k}}{|\boldsymbol{X}|^{5}}+12x_{3}\frac{{\it\beta}_{ik}y_{3}X_{j}+{\it\beta}_{ij}y_{3}X_{k}-{\it\delta}_{jk}x_{3}{\it\beta}_{i\ell }X_{\ell }}{|\boldsymbol{X}|^{5}}-60x_{3}y_{3}{\it\beta}_{i\ell }\frac{X_{j}X_{k}X_{\ell }}{|\boldsymbol{X}|^{7}},\end{eqnarray}$$
where
${\it\beta}_{ij}={\it\delta}_{ij}-2{\it\delta}_{3i}{\it\delta}_{3j}$
is the reflection operator,
$\boldsymbol{y}^{\ast }={\bf\beta}\boldsymbol{y}$
(and
$\boldsymbol{y}={\bf\beta}\boldsymbol{y}^{\ast }$
),
$\boldsymbol{X}={\bf\beta}(\boldsymbol{x}-\boldsymbol{y}^{\ast })$
and
$\hat{\boldsymbol{X}}=\boldsymbol{x}-\boldsymbol{y}$
. The expression (2.11) is the result of applying the Lorentz reflection (see Kim & Karrila Reference Kim and Karrila1991 or Kuiken Reference Kuiken1996) to the original stresslet, with some manipulation. The dimensionless force due to gravity acts at an angle
${\it\beta}$
relative to the wall,
$\boldsymbol{F}=\cos {\it\beta}\hat{\boldsymbol{x}}-\sin {\it\beta}\hat{\boldsymbol{z}}$
(the wall is parallel to gravity when
${\it\beta}=0$
), and
$\unicode[STIX]{x1D642}^{\ast }(\boldsymbol{x},\boldsymbol{x}^{\ast })\boldsymbol{\cdot }\boldsymbol{F}=\cos {\it\beta}\,\boldsymbol{G}_{x}^{\ast }(\boldsymbol{x},\boldsymbol{x}^{\ast })-\sin {\it\beta}\,\boldsymbol{G}_{z}^{\ast }(\boldsymbol{x},\boldsymbol{x}^{\ast })$
.
In the limit as the point
$\boldsymbol{x}$
tends towards a point on the boundary,
$\boldsymbol{x}\in S$
, the no-slip boundary condition on the body surface provides an integral equation to be solved for
$\boldsymbol{q}$
,
$$\begin{eqnarray}\displaystyle \boldsymbol{U}+{\it\bf\Omega}\times (\boldsymbol{x}-\boldsymbol{x}_{0}) & = & \displaystyle -\int _{S}\left(\boldsymbol{q}(\boldsymbol{y})-\boldsymbol{q}(\boldsymbol{x})\right)\boldsymbol{\cdot }(\unicode[STIX]{x1D61B}(\boldsymbol{x},\boldsymbol{y})+\unicode[STIX]{x1D61B}^{\ast }(\boldsymbol{x},\boldsymbol{y}^{\ast }))\boldsymbol{\cdot }\hat{\boldsymbol{n}}(\boldsymbol{y})\,\text{d}S\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{8{\rm\pi}}\left(\unicode[STIX]{x1D642}(\boldsymbol{x},\boldsymbol{x}_{0})+\unicode[STIX]{x1D642}^{\ast }(\boldsymbol{x},\boldsymbol{x}^{\ast })\right)\boldsymbol{\cdot }\boldsymbol{F}.\end{eqnarray}$$
The integrand is finite with a jump at the singular point,
$\boldsymbol{x}=\boldsymbol{y}$
. Further investigation of the integral operator leads to relations between the velocities and the density
$\boldsymbol{q}$
, closing the system:
$$\begin{eqnarray}\boldsymbol{U}=-\frac{4{\rm\pi}}{S_{A}}\int _{S}\boldsymbol{q}(\boldsymbol{x})\,\text{d}S,\quad {\it\bf\Omega}=-\mathop{\sum }_{m=1}^{3}\frac{4{\rm\pi}}{A_{m}}\boldsymbol{e}_{m}\left(\boldsymbol{e}_{m}\boldsymbol{\cdot }\int _{S}(\boldsymbol{x}-\boldsymbol{x}_{0})\times \boldsymbol{q}(\boldsymbol{x})\,\text{d}S\right),\end{eqnarray}$$
where
$S_{A}$
is the surface area of the particle,
$\boldsymbol{e}_{m}$
is the
$m$
th Cartesian unit vector and
$A_{m}=\int _{S}|(\boldsymbol{e}_{m}\times (\boldsymbol{x}-\boldsymbol{x}_{0}))|^{2}\,\text{d}S$
(see Pozrikidis Reference Pozrikidis1992).
To solve the integral equations above we use a collocation scheme, enforcing the equations at the nodes of the quadrature rule used to approximate the surface integrals. The body surface is parameterized using a spherical coordinate system. Integration is performed with respect to the zenith angle using Gaussian quadrature with
$N_{{\it\phi}}$
points and with respect to the azimuthal direction using the trapezoidal rule with
$N_{{\it\theta}}$
points, with
$N_{{\it\theta}}$
varying with the dimensionless ring circumference
$R\in (0,1)$
so as to achieve a roughly uniform distribution over the surface; namely,
$N_{{\it\theta}}$
is taken to be the greatest integer not exceeding
$N_{{\it\phi}}R$
in the prolate case or
$2.5N_{{\it\phi}}R^{2}$
in the oblate case. The integrand in (2.12) is set to zero at the jump discontinuity, a convenient method pioneered by Power & Miranda (Reference Power and Miranda1987), resulting in a quadrature scheme that is second-order accurate in the grid spacing. Where time stepping is required we employ a second-order Runge–Kutta method. The grid spacing and time-step size are chosen based on the stiffness of the problem under investigation, and changes to the numerical results presented in the paper are negligible when compared with simulations with much finer resolution. A convergence study, comparisons with previously published numerical results, and other methods of validation are included as appendix B.
A major benefit of the double-layer boundary integral method with stresslet images is that the equation for the density
$\boldsymbol{q}$
is a Fredholm integral equation of the second kind and is therefore well-posed, unlike approaches built upon a single-layer or combined formulation (Pozrikidis Reference Pozrikidis1992; Stakgold & Holst Reference Stakgold and Holst2011). Moreover, the flow can be computed using adaptive quadrature for near-wall interactions without suffering from poor conditioning problems, and other subtle issues like the careful selection of a regularization parameter may be avoided. The method of stresslet images used in this paper is more accurate than the popular method of regularized Stokeslets with images by Ainley et al. (Reference Ainley, Durkin, Embid, Boindala and Cortez2008) in the tests performed in appendix B, and does not rely on an extra regularization parameter. A more complete discussion on the numerical method will be presented elsewhere.
3. A zoology of particle trajectories
We begin by conducting a numerical survey of the trajectories exhibited by sedimenting prolate and oblate bodies near a vertical or tilted wall. Figures 2 and 3 show a selection of representative body dynamics, which depend on the body shape, initial data and wall inclination angle. We will discuss each trajectory type in turn.
Figure 2. Trajectories of prolate spheroidal bodies sedimenting near walls, as determined by full numerical simulation: (a) symmetric glancing; (b) symmetric reversing; (c) periodic tumbling (the distance travelled along the wall has been scaled by a factor of 100); (d) stable sliding near an inclined wall. The white markers on the tumbling body illustrate the rotation. Initial data and body shapes in each case are given in appendix C.
Figure 3. Three-dimensional glancing (a,b) and reversing (c,d) of prolate (a,c) and oblate (b,d) bodies near a vertical wall. The black rectangle in the background of each frame represents a strip of the wall,
$\{(x,y,0):-2\leqslant y\leqslant 2\}$
. Gravity is parallel to the wall, i.e. vertical on the page; the horizontal axis is the
$y$
direction. The lateral movements are plotted to scale, while the movements in the
$x$
direction have been greatly reduced for visualization purposes. Animations of these four trajectories are included as supplementary data available at http://dx.doi.org/10.1017/jfm.2015.222. The initial data used to generate these trajectories are given in appendix C. The movies, along with movies of periodic tumbling and wobbling of nearly spherical prolate and oblate bodies, are included as supplementary data.
3.1. Glancing, reversing and tumbling near a vertical wall
Our investigation begins in the simplest setting, where the wall is parallel to gravity,
${\it\beta}=0$
, and the geometry has symmetry through the
$xz$
-plane,
${\it\phi}=0$
. In this case, all trajectories can be described by tracking the distance
$h$
from the particle centre to the wall together with the angle
${\it\theta}$
measuring rotation in the
$xz$
-plane. Figure 2(a) shows the glancing dynamics of a slender prolate spheroid of eccentricity
$e=0.98$
, placed initially at a distance
$h=3$
from the wall and at an orientation angle
${\it\theta}=-20^{\circ }$
. Due to the drag anisotropy of slender bodies in viscous flows, the body initially drifts towards the wall. Hydrodynamic interactions with the surface then cause the particle to rotate until
${\it\theta}=0$
, at which point, in accordance with the time-reversal symmetry of the Stokes equations, the body continues to rotate and migrates away from the surface along a trajectory symmetric with its initial approach. As the particle escapes from the wall, its rotation rate diminishes so that
${\it\theta}$
tends towards a constant value. The same body follows a markedly different trajectory if the initial orientation angle is larger, as shown in figure 2(b), a reversing trajectory. In this example the initial orientation is
${\it\theta}\approx -70^{\circ }$
, the body rotates in the opposite direction, and the leading edge becomes the trailing edge after the closest approach to the wall.
These glancing and reversing dynamics were explored numerically, analytically and experimentally by Russel et al. (Reference Russel, Hinch, Leal and Tieffenbruck1977) for very slender particles. In the trajectories described in that work, a particle released far from the wall with
${\it\theta}\in (-{\rm\pi}/2,0)$
always approaches the wall, rotates, and then escapes from the wall, just as shown in figure 2(a,b). Russel et al. (Reference Russel, Hinch, Leal and Tieffenbruck1977) distinguished glancing from reversing trajectories according to the orientation of the particle at closest approach to the wall; in the former the particle is oriented parallel to the wall at closest approach while in the latter the particle is oriented normal to the wall at closest approach. The distinction between glancing and reversing can also be described in terms of a vector
$\boldsymbol{d}$
aligned with the axis of body symmetry:
$\hat{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{d}$
does not change sign in prolate glancing and oblate reversing orbits, while it does change sign in oblate glancing and prolate reversing orbits.
A third trajectory type prevails for a nearly spherical body released near the wall, as shown in figure 2(c), where we take
$e=0.15$
. On releasing the body at
$h=3$
with
${\it\theta}=0$
, we observe a new type of dynamics, a slow periodic tumbling motion (the distance travelled along the wall is scaled by a factor of 100 in figure 2
c). This behaviour can be understood as a perturbation of the well-known trajectory of a sphere near a vertical wall, i.e. translation in the direction of gravity together with a rolling-type rotation due to the torque induced by the presence of the wall. The slight eccentricity of the body causes a drag anisotropy which in turn leads to a migration velocity of the particle either towards or away from the surface, depending on the orientation angle. This slight migration, in concert with the rolling-type rotation, leads to a periodic tumbling trajectory. The dynamics is similar to the periodic tumbling of two identical non-spherical bodies placed side by side, as studied by Kim (Reference Kim1985, Reference Kim1986) and Jung et al. (Reference Jung, Spagnolie, Parikh, Shelley and Tornberg2006), with an important difference: here the rotation of the body is in the opposite direction,
$\dot{{\it\theta}}<0$
(relating to the opposing orientation of the Stokeslet in the image system in (2.5)).
3.2. Sliding along an inclined wall
Another type of particle trajectory arises when the bounding wall is not parallel to gravity. Figure 2(d) shows the dynamics of a prolate spheroid with
$e=0.98$
near a wall that is tilted at an angle
${\it\beta}=9.17^{\circ }$
(the initial data for the cases shown are included as appendix C). Here, a behaviour appears that does not exist for sedimentation near a vertical wall, which we term sliding. The body settles into a steady motion with a fixed orientation angle and distance from the wall. In this quasi-steady equilibrium the horizontal velocity induced by the particle orientation exactly balances the approach of the wall as the particle falls, and the rotation of the body due to the interaction with the wall has vanished; this type of trajectory was observed numerically by Hsu & Ganatos (Reference Hsu and Ganatos1994). This quasi-steady state owes its existence to the breaking of gravity–wall symmetry in connection with the choice
${\it\beta}>0$
, which weakens the consequences of time reversibility on the dynamics.
Assuming
${\it\phi}=0$
, the dynamics remains two-dimensional and it is natural to ask whether the glancing, reversing and periodic tumbling trajectories found near a vertical wall still may be found near an inclined wall, and if so how they share the phase space with the sliding trajectory. Although not shown in figure 2, the combination of small wall inclination angle
${\it\beta}$
and large eccentricity
$e$
allows both glancing- and reversing-like trajectories to occur in the full numerical simulations. However, these trajectories are less symmetric in that the limiting orientation angle after the wall encounter is no longer the opposite of the value before the wall encounter for the same distance to the wall; instead, the wall interaction tends to focus the orientations of escaping particles into a narrow band of escape angles. As
${\it\beta}$
increases or
$e$
decreases, it becomes increasingly difficult for the particle to escape from the wall, and the concentration of escape angles increases until it yields an attracting fixed point, namely the sliding trajectory discussed above. For still larger
${\it\beta}$
the equilibrium particle–wall gap size becomes extremely small, resulting in excessive computational costs and possible wall impact, and we do not study this regime. On the other hand, a careful tuning of
${\it\beta}$
against particle eccentricity can produce geometries where the fixed point is arbitrarily far from the wall, yet finite. In § 5 we derive analytical results quantifying this phenomenon, illustrated in figure 8.
The periodic tumbling orbits mentioned earlier no longer exist with
${\it\beta}>0$
. Instead, a nearly spherical body is found to rotate in nearly periodic orbits, but with a slow drift towards the wall (for
${\it\beta}>0$
) until eventually the orbit approaches the wall very closely. These initially near-periodic trajectories may be of more mathematical interest than practical application, since the region in parameter space where they arise is so limited.
3.3. Three-dimensional glancing, reversing and wobbling
In general, a particle near a surface will undergo lateral translations and out-of-plane rotations, leading to a fully three-dimensional trajectory. Consider a particle falling near a vertical wall,
${\it\beta}=0$
, but with no lateral symmetry,
${\it\phi}\neq 0$
. Four such trajectories are shown in figure 3, where prolate and oblate spheroids with aspect ratio
$a/c=2$
have been released with non-zero values of both
${\it\theta}$
and
${\it\phi}$
. A movie depicting these trajectories is included as supplementary data. Depending on the initial data, we again observe glancing- and reversing-like trajectories wherein the particle approaches and then escapes from the wall; in three dimensions, glancing and reversing trajectories can be classified just as in the two-dimensional case, in terms of a vector
$\boldsymbol{d}$
aligned with the axis of body symmetry. Once again,
$\hat{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{d}$
does not change sign in prolate glancing and oblate reversing orbits, while it does change sign in oblate glancing and prolate reversing orbits. Unlike in the two-dimensional case, the particle can drift laterally, subject to the constraint of symmetry about the moment of closest approach of the centroid to the wall. Importantly, the lateral behaviour depends on whether the particle is prolate or oblate. A glancing prolate body and a reversing oblate body continue to drift laterally without a change in the direction of drift (figure 3
a,d), while a glancing oblate body and a reversing prolate body return in the direction from which they came at the point of closest approach (a sign change in
${\dot{y}}$
, figure 3
b,c).
The periodic tumbling trajectory also has fully three-dimensional analogues. Included in the supplementary data is a movie showing the nearly spherical prolate tumbling and oblate tumbling trajectories which have a periodic lateral wobble with zero net lateral drift. These tumbling orbits, rotated away from the two-dimensional dynamics previously described, are now found to undergo periodic lateral motions in the
$y$
direction. As in the two-dimensional case, the trajectory can be understood as a combination of spherical rolling and reversing. The difference in prolate and oblate lateral drift in figure 3(c,d) also emerges in the three-dimensional tumbling orbits, so that the body changes lateral direction at the point of closest approach in the prolate case and at the point farthest from the wall in the oblate case.
In the more general setting with
${\it\phi}\neq 0$
and
${\it\beta}>0$
we have observed in numerical simulations that for small
${\it\beta}$
(small wall tilt angle) the wall interactions induce a concentration of the three-dimensional dynamics (escape angles tend towards a narrower band). For larger values of
${\it\beta}$
we see the emergence of an attracting fixed point. The wall inclination damps
${\it\phi}$
towards 0, and the fixed point is the same as in the case of lateral symmetry, as illustrated in figure 2(d). We will return to all of these trajectory types once we have developed analytical expressions with which to study them.
4. The method of images for wall-bounded Stokes flow
The numerical investigations described above are somewhat computationally intensive; for each time step in a trajectory a large linear system representing the discrete version of the surface integral equation must be inverted. At the same time, the dynamics can be fully described by tracking two or three scalar parameters, and the derivation of ordinary differential equations describing their dynamics would be of considerable value. To obtain an explicit system of differential equations which can be rapidly integrated or further studied analytically, we will apply the method of images and the method of reflections to an arbitrarily oriented prolate or oblate spheroid near a vertical or inclined wall.
The method of reflections takes an especially convenient form when the flows are constructed from systems of fundamental singularity solutions of the Stokes equations. The flow due to the motion of a spheroidal body in an infinite fluid may be represented by a collection of singularities placed at points interior to the body surface; image systems are then placed at the reflections of these points inside the wall to enforce the no-slip boundary condition on the wall. A generalization of Faxén’s law then gives the effect of this auxiliary velocity field on the body as a first-order correction of the trajectory due to the wall. The process may be continued to develop higher-order approximations of the effect of the wall on the body trajectory (see Kim & Karrila Reference Kim and Karrila1991). Wakiya (Reference Wakiya1959) carried out a similar procedure using Lamb’s solutions in ellipsoidal coordinates, producing expressions for the force and torque on a body moving with lateral symmetry near a wall.
The flow field associated with a spheroidal body in an unbounded fluid may be represented by an integrated distribution of image singularities on the centreline (prolate case) or a circular disk (oblate case). However, far from the particle,
$r=|\boldsymbol{x}-\boldsymbol{x}_{0}|\gg 1$
, this velocity field may be written as a multipole expansion of singularities placed at the body centroid. As shown in Kim & Karrila (Reference Kim and Karrila1991), the dimensionless fluid velocity far from the body is given by
$$\begin{eqnarray}\boldsymbol{u}^{(0)}(\boldsymbol{x})=\frac{3}{4}\left(\frac{\sinh (D)}{D}\right)\unicode[STIX]{x1D642}(\boldsymbol{x},\boldsymbol{x}_{0})\boldsymbol{\cdot }\boldsymbol{F},\end{eqnarray}$$
where
$\boldsymbol{F}=\cos {\it\beta}\hat{\boldsymbol{x}}-\sin {\it\beta}\hat{\boldsymbol{z}}$
is the external gravitational force on the body,
$\unicode[STIX]{x1D642}$
is the Stokeslet singularity given in (2.4) and
$$\begin{eqnarray}\frac{\sinh (D)}{D}=\mathop{\sum }_{n\geqslant 0}\frac{1}{(2n+1)!}D^{2n}=1+\frac{1}{6}D^{2}+\cdots \,,\end{eqnarray}$$
with
$D^{2}=\partial _{XX}+(b/a)^{2}\partial _{YY}+(c/a)^{2}\partial _{ZZ}$
. Truncation of the series after the two terms shown above results in errors in the flow (from (4.1)) that scale as
$r^{-5}$
as
$r\rightarrow \infty$
. It will prove useful to transform the differential operator
$D^{2}$
, which is diagonalized in the coordinate system of the spheroidal body axes, into the usual coordinate system oriented with the wall at
$\{z=0\}$
. Given the definitions of
${\it\theta}$
and
${\it\phi}$
from (2.1), the second derivatives can be written as linear combinations of derivatives along the standard axes,
$$\begin{eqnarray}\displaystyle \partial _{XX} & = & \displaystyle \cos ^{2}{\it\theta}\cos ^{2}{\it\phi}\,\partial _{xx}+\cos ^{2}{\it\theta}\sin ^{2}{\it\phi}\,\partial _{yy}+\sin ^{2}{\it\theta}\,\partial _{zz}\nonumber\\ \displaystyle & & \displaystyle +\,\sin (2{\it\theta})\left(\cos {\it\phi}\,\partial _{xz}+\sin {\it\phi}\,\partial _{yz}\right)+\cos ^{2}{\it\theta}\sin (2{\it\phi})\,\partial _{xy},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \partial _{YY} & = & \displaystyle \sin ^{2}{\it\phi}\,\partial _{xx}+\cos ^{2}{\it\phi}\,\partial _{yy}-\sin (2{\it\phi})\,\partial _{xy},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \partial _{ZZ} & = & \displaystyle \sin ^{2}{\it\theta}\cos ^{2}{\it\phi}\,\partial _{xx}+\sin ^{2}{\it\theta}\sin ^{2}{\it\phi}\,\partial _{yy}+\cos ^{2}{\it\theta}\,\partial _{zz}\nonumber\\ \displaystyle & & \displaystyle -\,\sin (2{\it\theta})\left(\cos ({\it\phi})\,\partial _{xz}+\sin {\it\phi}\,\partial _{yz}\right)+\sin ^{2}({\it\theta})\sin (2{\it\phi})\,\partial _{xy}.\end{eqnarray}$$
$x$
- and
$z$
-directed Stokeslets and selected second derivatives. To obtain an image system, we employ Blake’s image system for the Stokeslet and expressions such as (2.8) for the required second derivatives. This leads to an analytical expression for the reflection flow,
$\boldsymbol{u}^{(1)}(\boldsymbol{x})$
. Since we have truncated the series in (4.2), neglecting terms of size
$D^{4}(r^{-1})$
, the error in this reflected flow scales as
$(r^{\ast })^{-5}$
for
$r^{\ast }\rightarrow \infty$
, where
$r^{\ast }=|\boldsymbol{x}-\boldsymbol{x}^{\ast }|$
and
$\boldsymbol{x}^{\ast }=(0,0,-h)$
.
Finally, the effect of the reflected flow on the original spheroidal particle is given by the mobility relations between the particle motion
$(\boldsymbol{U},{\it\bf\Omega})$
and the external force and torque
$(\boldsymbol{F},\boldsymbol{T})$
through a generalized Faxén law (see Kim & Karrila Reference Kim and Karrila1991),
$$\begin{eqnarray}\displaystyle -\boldsymbol{F} & = & \displaystyle (X^{A}\boldsymbol{d}\boldsymbol{d}^{\text{T}}+Y^{A}(I-\boldsymbol{d}\boldsymbol{d}^{\text{T}}))\left(\frac{\sinh (D)}{D}\boldsymbol{u}^{(1)}(\boldsymbol{x}_{0})-\boldsymbol{U}\right),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle -\boldsymbol{T} & = & \displaystyle \frac{2}{3}(X^{C}\boldsymbol{d}\boldsymbol{d}^{\text{T}}+Y^{C}(I-\boldsymbol{d}\boldsymbol{d}^{\text{T}}))\left(\left.\frac{3}{D}\frac{\partial }{\partial D}\left(\frac{\sinh (D)}{D}\right)\boldsymbol{{\rm\nabla}}\times \boldsymbol{u}^{(1)}\right|_{\boldsymbol{x}_{0}}-2{\it\bf\Omega}\right)\nonumber\\ \displaystyle & & \displaystyle -\,\frac{4}{3}Y^{H}\left\{\frac{3}{D}\frac{\partial }{\partial D}\left(\frac{\sinh (D)}{D}\right)\boldsymbol{\cdot }\boldsymbol{E}^{(1)}(\boldsymbol{x}_{0})\boldsymbol{\cdot }\boldsymbol{d}\right\}\times \boldsymbol{d},\end{eqnarray}$$
$\boldsymbol{E}^{(1)}=(\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{(1)}+[\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{(1)}]^{\text{T}})/2$
is the symmetric rate-of-strain tensor and
$\boldsymbol{d}$
is a unit vector oriented along the particle’s axis of symmetry, with 
$$\begin{eqnarray}\boldsymbol{d}=\left\{\begin{array}{@{}ll@{}}(\cos {\it\theta}\cos {\it\phi},\cos {\it\theta}\sin {\it\phi},\sin {\it\theta})\quad & \text{for prolate bodies,}\\ (-\!\sin {\it\theta}\cos {\it\phi},-\!\sin {\it\theta}\sin {\it\phi},\cos {\it\theta})\quad & \text{for oblate bodies}.\end{array}\right.\end{eqnarray}$$
The constants
$X^{A},Y^{A},X^{C},Y^{C}$
and
$Y^{H}$
depend only on the eccentricity
$e$
and whether the particle is prolate or oblate, and are included in table 1. The dimensionless torque
$\boldsymbol{T}$
is the result of scaling upon
$a|\boldsymbol{F}^{\ast }|$
. The prolate and oblate problems are solved together in a single calculation through the introduction of a parameter that is equal to 1 for prolate bodies and
$-1$
for oblate bodies. Each component of
$\boldsymbol{d}\boldsymbol{d}^{\text{T}}$
can then be rewritten using appropriate trigonometric identities; for example,
$d_{1}d_{1}$
reduces to
$\cos ^{2}{\it\phi}(1\pm \cos (2{\it\theta}))/2$
.
Table 1. Parameter definitions for the prolate and oblate cases in the Faxén laws (4.6) and (4.7) and the general system (A 1)–(A 10), from Kim & Karrila (Reference Kim and Karrila1991). Of these, only
$X^{A}$
and
$Y^{A}$
appear directly in (A 1)–(A 10) because we have exploited the relation
$Y^{H}/Y^{C}=\pm e^{2}/(2-e^{2})$
. The quantity
$\pm (X^{A}-Y^{A})$
is negative for both body types.
By setting
$\boldsymbol{T}=0$
in (4.7) to solve the sedimentation problem of interest and truncating the differential operators in (4.6) and (4.7) at second order using (4.2) and
$$\begin{eqnarray}\frac{3}{D}\frac{\partial }{\partial D}\left(\frac{\sinh (D)}{D}\right)=1+\frac{D^{2}}{10}+\cdots \!,\end{eqnarray}$$
we obtain linear equations for
$\boldsymbol{U}$
and
${\it\bf\Omega}$
. Solution of those equations and extraction of
$\dot{{\it\theta}}$
and
$\dot{{\it\phi}}$
from
${\it\bf\Omega}$
results in a set of ordinary differential equations governing the dynamics of the body. The full system, for an inclined wall, is given in appendix A together with the formula for
${\it\bf\Omega}$
. When the wall is vertical, the system reduces to an especially tidy form:
$$\begin{eqnarray}\displaystyle {\dot{x}} & = & \displaystyle \frac{(2-\cos ^{2}{\it\phi}(1\pm \cos (2\,{\it\theta})))(X^{A}-Y^{A})+2Y^{A}}{2\,X^{A}Y^{A}}-\frac{9}{16\,h}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{2\,e^{2}(\cos (2\,{\it\theta})\pm 1)\cos ^{2}{\it\phi}+18\,e^{2}\cos ^{2}{\it\theta}-e^{2}(17\pm 7)+16}{128\,h^{3}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\dot{y}} & = & \displaystyle -\frac{\sin (2{\it\phi})(1\pm \cos (2\,{\it\theta}))(X^{A}-Y^{A})}{4X^{A}Y^{A}}+\frac{e^{2}\sin (2{\it\phi})(\cos (2\,{\it\theta})\pm 1)}{128\,h^{3}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\dot{h}} & = & \displaystyle \frac{\pm \cos {\it\phi}\sin (2{\it\theta})(Y^{A}-X^{A})}{2X^{A}Y^{A}}-\frac{e^{2}\cos {\it\phi}\sin (2{\it\theta})}{32h^{3}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \dot{{\it\theta}} & = & \displaystyle \frac{9e^{2}\cos {\it\phi}\cos (2{\it\theta})}{32(2-e^{2})h^{2}}-\frac{3\cos {\it\phi}}{64(2-e^{2})h^{4}}\nonumber\\ \displaystyle & & \displaystyle \times \,[4-10e^{2}+(7\pm 1)e^{4}+e^{2}\cos ^{2}{\it\theta}(9e^{2}\cos ^{2}{\it\theta}-(15\pm 2)e^{2}+12)],\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \dot{{\it\phi}} & = & \displaystyle \left\{\begin{array}{@{}ll@{}}\displaystyle \frac{3\sin {\it\phi}\tan {\it\theta}}{64(2-e^{2})}\left(\frac{-6e^{2}}{h^{2}}+\frac{3e^{4}\cos ^{2}{\it\theta}-8e^{4}+10e^{2}-4}{h^{4}}\right)\quad & (\text{prolate}),\\ \displaystyle \frac{3\sin {\it\phi}\cot {\it\theta}}{64(2-e^{2})}\left(\frac{-6e^{2}}{h^{2}}-\frac{3e^{4}\sin ^{2}{\it\theta}-2e^{4}-2e^{2}-4}{h^{4}}\right)\quad & (\text{oblate}),\end{array}\right.\end{eqnarray}$$
$\pm$
signs should be replaced with
$+$
in the prolate case and
$-$
in the oblate case, and the constants
$X^{A}$
and
$Y^{A}$
have different definitions in the two cases, as indicated in table 1. Importantly, the derivatives of the particle–wall distance
$h$
and of the angles
${\it\theta}$
and
${\it\phi}$
are independent of the positions
$x$
and
$y$
, so that the system is fundamentally three-dimensional; the positions
$\{x(t),y(t)\}$
may be determined directly once
$\{h(t),{\it\theta}(t),{\it\phi}(t)\}$
have been found. The errors in the expressions above and in the general setting (in appendix A) are
$O(h^{-4})$
in the translational velocity and
$O(h^{-5})$
in the rotational velocity for
$h\gg 1$
.
The full expression of the rotational velocity
${\it\bf\Omega}$
could be used to deduce the rotation of the body about its axis of symmetry, a third angle that in addition to
${\it\theta}$
and
${\it\phi}$
prescribes the precise history of the body as it evolves in time. However, the translational and rotational velocities computed at any moment are invariant to rotations about the axis of symmetry, so that this third angle may be determined after solving for
$\{h(t),{\it\theta}(t),{\it\phi}(t)\}$
just as may be done for the drift positions
$x(t)$
and
$y(t)$
.
5. Analysis of particle trajectories
The ordinary differential equations describing the body dynamics can be integrated numerically, and in some cases analytically, to produce approximate trajectories for the sedimentation problem in the general setting. In this section we will derive analytical formulae for the particle trajectory in various special cases, beginning with the assumption of two-dimensional dynamics and then proceeding to the general case.
5.1. Analysis of glancing, reversing and tumbling dynamics
Consider first the case of two-dimensional motion,
${\it\phi}=0$
, near a vertical wall,
${\it\beta}=0$
. The evolution of the particle position and orientation is governed by the reduced system (from (4.12) and (4.13))
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{{\it\theta}}=\frac{\cos (2{\it\theta})}{h^{2}}\left[A-\frac{B}{h^{2}}-C\frac{\cos (2{\it\theta})}{h^{2}}\right]-\frac{D}{h^{4}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\dot{h}}=\sin (2{\it\theta})\left[E-\frac{F}{h^{3}}\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle A=\frac{9e^{2}}{32(2-e^{2})},\quad B=\frac{3e^{2}(6-(3\pm 1)e^{2})}{64(2-e^{2})},\quad C=\frac{27e^{4}}{256(2-e^{2})},\\ \displaystyle D=\frac{48-48e^{2}+21e^{4}}{256(2-e^{2})},\quad E=\pm \frac{Y^{A}-X^{A}}{2X^{A}Y^{A}},\quad F=\frac{e^{2}}{32}.\end{array}\right\}\end{eqnarray}$$
The limiting case of a spherical body is dramatically simpler, with
$A=B=C=E=F=0$
and
$D=3/32$
. For a general particle eccentricity
$e$
the system has a fixed point at
${\it\theta}=0$
and a particle–wall distance that satisfies
$$\begin{eqnarray}h^{2}=\frac{B+C+D}{A}=\frac{4+2e^{2}-(-1\pm 1)e^{4}}{6e^{2}}.\end{eqnarray}$$
This unstable fixed point corresponds to a particle aligned with the wall and falling vertically without rotating, and can occur only at a specific particle–wall distance where the competing dynamics that give rise to glancing and reversing precisely balance each other. The corresponding particle–wall gap size
$h-\sqrt{1-e^{2}}$
from (5.4) is plotted against particle eccentricity in figure 4, along with the values computed using the full numerical simulation, shown as symbols. The equilibrium distance is unbounded as
$e\rightarrow 0$
(a sphere always rotates in the same direction at any finite distance from the wall). As particle eccentricity increases, the numerically determined gap size decreases and then vanishes, and the accuracy of the estimate from (5.4) is poor for large particle eccentricity where the particle equilibrium distance is very close to the wall.
Figure 4. The gap size for which a body with
${\it\theta}=0$
does not rotate as a function of particle eccentricity, which distinguishes the transition from glancing to periodic tumbling. The results of full numerical simulations are shown as symbols and those according to (5.4) as lines. The results for prolate and oblate bodies are remarkably similar, with the accuracy of the analytical estimates for both deteriorating for large particle eccentricity where the equilibrium distance is very close to the wall.
A reduced but analytically tractable approximation of the system above is found in the limit of large particle distances from the wall and upon inspection of the coefficients. It is appealing to keep only the terms of size
$O(h^{-2})$
in (5.1) and of size
$O(1)$
in (5.2), but higher-order terms become dominant when
$\cos (2{\it\theta})=0$
. Moreover, for nearly spherical particles,
$e\approx 0$
, a more appropriate comparison of terms involves the ratio
$e/h^{2}$
; for instance,
$B/h^{4}\sim e^{2}/h^{4}\ll D/h^{4}$
for
$e\ll 1$
and
$h\gg 1$
. Neglecting the terms with coefficients
$B,C$
and
$F$
results in the reduced system
$$\begin{eqnarray}\dot{{\it\theta}}=\frac{A\cos (2{\it\theta})}{h^{2}}-\frac{D}{h^{4}},\quad {\dot{h}}=E\sin (2{\it\theta}).\end{eqnarray}$$
This system is autonomous, and the two derivatives can be divided to obtain a single first-order equation for
$\text{d}{\it\theta}/\text{d}h$
. The transformations
${\it\gamma}=1/h$
and
${\it\eta}=-\!\cos (2{\it\theta})/2$
then yield a linear differential equation,
$$\begin{eqnarray}\frac{\text{d}{\it\eta}}{\text{d}{\it\gamma}}=\frac{2A}{E}{\it\eta}+\frac{D}{E}{\it\gamma}^{2}.\end{eqnarray}$$
Multiplying by
$\exp (-2A{\it\gamma}/E)$
and integrating leads to
$$\begin{eqnarray}{\it\eta}\exp \left(-\frac{2A}{E}{\it\gamma}\right)=\frac{D}{E}\exp \left(-\frac{2A}{E}{\it\gamma}\right)\left[\frac{-E}{2A}{\it\gamma}^{2}-\frac{E^{2}}{2A^{2}}{\it\gamma}-\frac{E^{3}}{4A^{3}}\right]+c_{0},\end{eqnarray}$$
where
$c_{0}$
is a constant of integration. For each trajectory the relation
$$\begin{eqnarray}2c_{0}=\exp \left(-\frac{2A}{E}{\it\gamma}\right)\left(2{\it\eta}+\frac{D}{A}\left({\it\gamma}^{2}+\frac{E}{A}{\it\gamma}+\frac{E^{2}}{2A^{2}}\right)\right)\end{eqnarray}$$
holds, and therefore each trajectory must follow a level set of the function
$$\begin{eqnarray}{\it\Psi}(h,{\it\theta})=\exp \!\left(-\frac{2A}{Eh}\right)\left(-\!\cos (2{\it\theta})+\frac{D}{A}\left(h^{-2}+\frac{E}{Ah}+\frac{E^{2}}{2A^{2}}\right)\right)\end{eqnarray}$$
in the
${\it\theta}h$
-plane.
Figure 5. Two-dimensional trajectories of prolate spheroids sedimenting near a vertical wall are depicted by plotting the particle–wall distance
$h$
against the orientation angle
${\it\theta}$
, for (a)
$e=0.02$
, (b)
$e=0.15$
, (c)
$e=0.30$
and (d)
$e=0.80$
. Unphysical regions with
$h\leqslant \sqrt{\sin ^{2}{\it\theta}+(1-e^{2})\cos ^{2}{\it\theta}}$
, corresponding to body penetration through the wall, are shaded. The results of the full numerical simulations are shown as red symbols. For small
$e$
we observe periodic orbits near the wall (circles). As
$e$
increases, the periodic trajectories are replaced by reversing (squares) and glancing (triangles) trajectories. Arrows indicate the direction of time. The contours of the scalar function
${\it\Psi}$
from (5.9), shown as black lines, give accurate predictions of the full numerical results.
Figure 5 shows the level sets of (5.9) as black lines, together with the results of the full numerical simulations (see § 2.2) as red symbols, for prolate spheroids of four different eccentricities:
$e\in \{0.02,0.15,0.3,0.8\}$
. Periodic tumbling orbits are indicated by circles, reversing trajectories by squares and glancing trajectories by triangles. Arrows indicate the direction of time. For
$e=0.02$
the particle is nearly spherical and we see the periodic orbits described earlier. In these orbits the particle is farthest from the wall when
${\it\theta}=0$
, i.e. when the major axis is parallel to the wall. As noted in § 3, the period of the tumbling orbit is extremely long; using (5.1) for the limiting case of a sphere,
$e\rightarrow 0$
, we find the period of full rotation
$T=64{\rm\pi}h_{0}^{4}/3$
, where
$h_{0}$
is the constant distance from the wall. The sedimentation distance
$X$
travelled during one period of a tumbling orbit, using (4.10), is
$$\begin{eqnarray}X=\frac{64{\rm\pi}h_{0}^{3}}{3}\left(1-\frac{9}{16h_{0}}+O(h_{0}^{-3})\right).\end{eqnarray}$$
For eccentricities
$e=0.15$
and 0.3 in figure 5, the periodic orbits are restricted to a narrower region near the wall, while glancing and reversing trajectories approach ever nearer to the wall before turning back. Finally, with
$e=0.8$
, we see only glancing- and reversing-type trajectories, in concurrence with the slender-body work in Russel et al. (Reference Russel, Hinch, Leal and Tieffenbruck1977); figure 5(a) can be compared directly with figure 4 in that work, bearing in mind that the
${\it\theta}$
defined therein is the same as our
$-{\it\theta}$
.
These contours concur with the numerical survey in classifying the trajectories in this symmetric version of the problem into glancing, reversing and tumbling types. The quantitative agreement with the numerical solutions is generally very good, though imperfect in some cases where the body comes very close to the wall. Our experience suggests that the results of the reduced equations should be used with some caution for
$h<2$
, though the dynamics generally remains qualitatively sound for much smaller values of
$h$
. Replacement of the contours of
${\it\Psi}$
with trajectories determined by numerical integration of the fourth-order system (5.1) and (5.2) results in only minor changes.
From the analytical picture above we are in a position to predict solely from initial conditions whether or not a particle will escape, and if so to predict the final orientation it will assume far from the wall. The initial data
$h=h_{0}$
and
${\it\theta}={\it\theta}_{0}$
determine a level set of
${\it\Psi}$
, and if the particle escapes this contour must have an asymptote with
$h\rightarrow \infty$
. In this limit we obtain from (5.9) the relation
$$\begin{eqnarray}\cos (2{\it\theta})=\frac{DE^{2}}{2A^{3}}-{\it\Psi}(h_{0},{\it\theta}_{0}).\end{eqnarray}$$
If the right-hand side of (5.11) has magnitude greater than one, there is no solution and a periodic orbit is predicted. Otherwise, (5.11) predicts the limiting orientation angle that the particle takes once it has escaped from the wall. Among escaping particles we can distinguish the glancing from the reversing trajectories by examining this asymptotic orientation angle more closely. That is, we consider the problem of determining the angle
${\it\theta}^{\ast }$
that divides glancing from reversing trajectories at a given eccentricity far from the wall. This can be done analytically by solving the equation
${\it\Psi}(\sqrt{D/A},0)=\lim _{h\rightarrow \infty }{\it\Psi}(h,{\it\theta}^{\ast })$
, since
$(h=\sqrt{D/A},{\it\theta}=0)$
is the fixed point in the reduced model used to derive
${\it\Psi}$
and the glancing–reversing separatrix passes through this fixed point. This gives an expression for
${\it\theta}^{\ast }$
in terms of the constants
$A,D,E$
defined in (5.3):
$$\begin{eqnarray}{\it\theta}^{\ast }=\frac{1}{2}\arccos \left(2{\it\kappa}^{-2}\left(1-\frac{{\it\kappa}+1}{\exp ({\it\kappa})}\right)\right),\quad \text{where}~{\it\kappa}=\frac{2A^{3/2}}{E\sqrt{D}}.\end{eqnarray}$$
Making the substitutions in (5.3) and in table 1, we can write
${\it\kappa}$
explicitly in terms of the eccentricity
$e$
:
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\kappa}_{prolate}=\frac{12\sqrt{6}e^{6}}{(2-e^{2})\sqrt{16-16e^{2}+7e^{4}}\left((3-e^{2})\log ((1+e)(1-e)^{-1})-6e\right)}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\kappa}_{oblate}=\frac{6\sqrt{6}e^{6}}{(2-e^{2})\sqrt{16-16e^{2}+7e^{4}}\left((3-2e^{2})\text{arccot}(\sqrt{1-e^{2}}/e)-3e\sqrt{1-e^{2}}\right)}.\qquad & \displaystyle\end{eqnarray}$$
To assess the accuracy of this analytical result we consider prolate spheroids and determine
${\it\theta}^{\ast }$
numerically for several values of
$e$
by computing trajectories that start from a fixed large
$h$
and various initial angles
${\it\theta}\in (-{\rm\pi}/2,0)$
and continue until
${\it\theta}$
reaches either 0 (glancing) or
$-{\rm\pi}/2$
(reversing). These numerical results are not obtained by integrating (4.12) and (4.13) but by solving the full Stokes equations using the numerical method described in § 2.2. The two values of the initial angle where the outcome changes from glancing to reversing determine an interval containing
${\it\theta}^{\ast }$
. These results are shown in figure 6. For
$e\ll 1$
the interval of uncertainty is quite small, and we report
${\it\theta}^{\ast }$
to an accuracy of
$0.1^{\circ }$
. For
$e\approx 1$
the trajectories of interest pass extremely close to the wall, resulting in excessive computational costs, so that the intervals to which we are able to constrain
${\it\theta}^{\ast }$
are large enough to be visible in figure 6. In an experimental setting, trajectories with initial angles within these ranges of uncertainty may result in wall collisions due to imperfections in the particle or wall geometry. The stars in figure 6, corresponding to an eccentricity
$e=0.99986$
, are the values reported by Russel et al. (Reference Russel, Hinch, Leal and Tieffenbruck1977) as the result of numerical work and experiments with aluminium wires of aspect ratio 60. The results compare well; for
$e\in \{0.1,0.3,0.5,0.7,0.9\}$
the formula (5.12) gives a result within one degree of the numerical range, and for
$e=0.999861$
the result is within one degree of the range reported by Russel et al. (Reference Russel, Hinch, Leal and Tieffenbruck1977). The theory predicts for oblate bodies a transition angle
${\it\theta}^{\ast }$
slightly greater than in the prolate case, but the difference is within one degree for
$e\leqslant 0.866$
.
Figure 6. The transition angle between glancing and reversing trajectories for prolate bodies as a function of eccentricity. The solid curve is from the explicit formula (5.12). Numerical results are shown as red circles; at low eccentricity the transition angle is well resolved but at greater eccentricities we report only a range of possible values. The two stars on the right indicate the maximum glancing angle and the minimum reversing angle reported by Russel et al. (Reference Russel, Hinch, Leal and Tieffenbruck1977), who in addition to a numerical study used aluminium wires of aspect ratio
$a/c\approx 60$
or
$e=0.999861$
. Between these angles they reported wall impacts. For oblate bodies, the theory predicts a slightly larger value of
${\it\theta}^{\ast }$
than in the prolate case, but the difference is less than one degree for
$e<0.866$
.
5.2. Analysis of three-dimensional dynamics near a vertical wall
The fully three-dimensional equations, while more complicated, can still be investigated analytically. With
${\it\phi}\neq 0$
(and
${\it\beta}=0$
as before), the discussion in the preceding subsection remains relevant because
$\dot{{\it\theta}}$
and
${\dot{h}}$
depend on
${\it\phi}$
only through the common factor
$\cos {\it\phi}$
in (4.12) and (4.13), and may be divided as before to make
$h$
the independent variable. This division now implicitly assumes that
$\cos {\it\phi}$
does not vanish on any open time interval, so we must note separately the steady solutions at
${\it\theta}=0$
,
${\it\phi}={\rm\pi}/2$
(prolate) and
${\it\theta}={\rm\pi}/2$
,
${\it\phi}={\rm\pi}/2$
(oblate). The argument then proceeds in the same way and results in an incomplete but still valuable description of the three-dimensional orbit: the projection of the trajectory in
$(h,{\it\theta},{\it\phi})$
-space onto the
$h{\it\theta}$
-plane must lie on a single level set of
${\it\Psi}$
in (5.9) as determined by the initial condition.
A difference between two- and three-dimensional trajectories that is visible in the
$h{\it\theta}$
-plane is that the projection of a three-dimensional trajectory may traverse only part of a contour of
${\it\Psi}$
instead of all of it. In particular, an orbit may be periodic even though the contour of
${\it\Psi}$
has an asymptote for some finite value of
${\it\theta}$
. This is possible because the periodic orbit repeatedly traverses a subset of the contour, doubling back on itself at regular intervals. To explain this behaviour, we note that
$\cos {\it\phi}>0$
implies that a particle with
${\it\theta}<0$
is moving towards the wall, whereas for
$\cos {\it\phi}<0$
the opposite holds. The result is that initial conditions with
${\it\phi}=0$
for which (5.11) predicts periodicity also lead to periodic trajectories when the initial
${\it\phi}$
is modified.
For a trajectory where
${\it\phi}\neq 0$
, (4.11) implies that
$\dot{y}\neq 0$
, i.e. the particle will move laterally. In the case of a periodic trajectory, the body wobbles periodically as it falls, drifting laterally back and forth along the wall as described in § 3.3. Plots of
$h(t)$
,
${\it\theta}(t)$
and
${\it\phi}(t)$
for such trajectories are shown in figure 7(a), where we have set
$e=0.05$
(prolate),
${\it\beta}=0^{\circ }$
, and initially
$(h,{\it\theta},{\it\phi})=(5,-50^{\circ },20^{\circ })$
. An animation of such a trajectory is also provided as supplementary material. The dynamics is akin to a three-dimensional reversing trajectory that fails to escape from the wall. When the body is closest to the wall,
${\it\theta}$
is nearly
$\pm {\rm\pi}/2$
, so that a small body rotation leads to a rapid change in
${\it\phi}$
from nearly zero to nearly
${\rm\pi}$
, or vice versa, and the body begins to drift laterally back in the direction from which it came. However, since the body is nearly spherical, the rotation induced by the wall is sufficient to rotate the major axis fast enough to redirect the body towards the wall yet again, and another reversing-type interaction ensues.
Figure 7. Three-dimensional trajectories of sedimenting prolate spheroids, from full numerical simulations (symbols) and integration of the fourth-order accurate ordinary differential equations (ODEs) from § 4 (lines). (a) Three-dimensional periodic tumbling orbit of a nearly spherical particle near a vertical wall:
$e=0.05$
,
${\it\beta}=0^{\circ }$
. The dynamics is akin to a three-dimensional reversing trajectory that fails to escape from the wall. The analytical prediction captures the shapes and amplitudes of the particle–wall distance
$h$
and the orientation angles
${\it\theta}$
and
${\it\phi}$
, but with an error in the frequency of oscillation. (b) A reversing trajectory of a more eccentric particle near a vertical wall:
$e=0.7$
,
${\it\beta}=0^{\circ }$
. The body visits the wall and departs, settling to a constant final orientation in both
${\it\theta}$
and
${\it\phi}$
as
$h\rightarrow \infty$
. (c) The same particle as in (b) near a slightly tilted wall converges to the stable sliding trajectory:
$e=0.7$
,
${\it\beta}=0^{\circ }$
. The particle initially rotates while approaching the wall and then recedes towards a limiting separation distance and orientation, with
${\it\phi}\rightarrow {\rm\pi}$
(the dynamics tends towards the two-dimensional sliding trajectory).
Figure 7(b) shows the trajectory for a more eccentric particle, with
$e=0.7$
, with the same initial condition,
$(h,{\it\theta},{\it\phi})=(5,-50^{\circ },20^{\circ })$
, which results in a complete three-dimensional reversing trajectory (also depicted in figure 3
c). Just as in the previous case, when the body reaches the point nearest to the wall there is a rapid rotation in
${\it\phi}$
, but in this case the body then ceases to rotate and escapes, settling to a constant orientation. This limiting orientation has an interesting relationship to the dynamics near an inclined wall, to which we now turn.
5.3. Analysis of the fully general problem and the sliding trajectory
We now consider the most general version of the problem, with an inclined wall (
${\it\beta}>0$
) and fully three-dimensional sedimentation dynamics (
${\it\phi}\neq 0$
). A non-zero wall inclination angle reduces the symmetry in the problem and weakens the constraints of time reversibility on the dynamics. One consequence is that there are no longer periodic orbits of the form discussed in the previous section. Another is that the three-dimensional dynamics can be driven towards the two-dimensional state, with
${\it\phi}$
reducing in magnitude to zero as
$t\rightarrow \infty$
. While the complete system presented in appendix A resists analytical treatment, the existence and stability of a sliding trajectory may still be investigated as follows.
Neglecting terms of
$O(h^{-3})$
in (A 3), (A 9) and (A 10) in appendix A, a fixed point (the sliding trajectory as depicted in figure 2
d) may be found explicitly. Taking
${\it\phi}=0$
(the two-dimensional laterally symmetric case) gives
$\dot{{\it\phi}}=0$
, and then the reduced expression for
$\dot{{\it\theta}}$
vanishes when
${\it\theta}={\it\theta}_{0}$
, where
$$\begin{eqnarray}{\it\theta}_{0}={\textstyle \frac{1}{2}}\,\tan ^{-1}\left({\textstyle \frac{2}{3}}\,\cot ({\it\beta})\right),\end{eqnarray}$$
an equation previously derived by Hsu & Ganatos (Reference Hsu and Ganatos1994) using expressions from Wakiya (Reference Wakiya1959). Finally, with
${\it\theta}={\it\theta}_{0}$
and
${\it\phi}=0$
the reduced expression for
${\dot{h}}=U_{z}$
vanishes when
$h=h_{0}$
, where
$$\begin{eqnarray}h_{0}=\frac{9\,X^{A}Y^{A}}{8Y^{A}\pm 4(X^{A}-Y^{A})\left((3+2\cot ^{2}{\it\beta})(9+4\cot ^{2}{\it\beta})^{-1/2}\pm 1\right)},\end{eqnarray}$$
where the
$\pm$
signs should be replaced by
$+$
for the prolate case and by
$-$
for the oblate case, and where
$X^{A}$
and
$Y^{A}$
also have different definitions, as indicated in table 1. On linearizing the reduced system about
$(h=h_{0},{\it\theta}={\it\theta}_{0},{\it\phi}=0)$
, this fixed point is found to be stable to arbitrary small perturbations for
${\it\beta}>0$
(the eigenvalues associated with the linear system may be shown to always be negative). As an example of a body that is attracted to this stable sliding trajectory, we consider again a particle of eccentricity
$e=0.7$
, but near a wall tilted at an angle
${\it\beta}=2.5^{\circ }$
. Details for the body trajectory are shown in figure 7(c), using once again the initial condition
$(h,{\it\theta},{\it\phi})=(5,-50^{\circ },20^{\circ })$
. Unlike in the vertical-wall case, the body approaches the wall and rotates in a reversing-type manner, but then settles to a finite wall separation distance as
$t\rightarrow \infty$
. Meanwhile, the rotations in
${\it\theta}$
and
${\it\phi}$
(and the distance
$h$
) are no longer symmetric about the time at which the centroid is closest to the wall. Since
$h$
remains bounded, the interaction with the wall continues to influence the rotational velocity of the body, and the body continues to rotate into the plane of the two-dimensional dynamics (
${\it\phi}\rightarrow {\rm\pi}$
).
More generally, the equilibrium particle–wall separation
$h_{0}$
is a function of the particle eccentricity and the wall inclination angle. The positive contours of
$h_{0}$
for prolate bodies are plotted in the
$e{\it\beta}$
-plane in figure 8. Numerical simulations indicate that the contours with
$h<2$
may overestimate the height of the fixed point, but the higher contours, near the boundary of the escaping trajectories, are reliable.
Figure 8. A prolate body near a sufficiently inclined wall cannot escape; it either approaches the wall so closely that particle or wall imperfections or other physics become important, or it assumes a stable orientation at a constant separation distance. In the latter case, the asymptotic separation
$h_{0}$
is a function of the wall inclination angle and the particle eccentricity, with contours plotted above. Near the boundary of the escaping trajectories in the
$e{\it\beta}$
-plane one can find arbitrarily large values of
$h_{0}$
.
A sliding trajectory exists if the equilibrium particle–wall separation
$h_{0}$
in (5.16) is positive and finite. To ascertain whether an eccentricity
$e$
and wall inclination angle
${\it\beta}$
result in a sliding trajectory, we set the denominator on the right-hand side of (5.16) to zero. Solving for
${\it\beta}$
in terms of
$e$
gives a critical wall inclination
${\it\beta}^{\ast }(e)$
beyond which the sliding trajectory arises. For
$0<{\it\beta}<{\it\beta}^{\ast }(e)$
, the wall is sufficiently vertical that the particle may escape if the initial condition is suitable. In this case
$h_{0}$
is negative and the fixed point does not describe physical behaviour. The glancing and reversing trajectories are similar to their vertical-wall counterparts shown in figure 3, except in that the orientation of the particle after the wall encounter need not be symmetric with its orientation beforehand. The rare cases of tumbling-type particles with positive
${\it\beta}$
exhibit perhaps the richest and most complex dynamics due to the low symmetry constraints. These trajectories are not perfectly periodic; there is a gradual approach to the wall together with increased rotation rate until, possibly after many full revolutions, the particle collides with the wall.
Meanwhile, if
${\it\beta}>{\it\beta}^{\ast }(e)$
, escape is impossible and a particle beginning from any initial condition instead approaches the sliding fixed point whose coordinates are given (according to the
$O(h^{-2})$
theory) above. In a numerical study, Kutteh (Reference Kutteh2010) reported a second critical value of
${\it\beta}$
beyond which the sliding trajectory disappears and ‘the particle monotonically approaches the wall until it makes contact’. The analytical results presented here simply indicate a small particle–wall equilibrium distance (the small contour heights in figure 8), but the underlying assumptions are not suitable for modelling close particle–wall interactions.
The long formula for
${\it\beta}^{\ast }(e)$
(not shown here) reproduces the four values given in table 9 of Hsu & Ganatos (Reference Hsu and Ganatos1994) for prolate and oblate bodies of aspect ratios
$c/a\in \{0.1,0.5\}$
, providing a useful check on the method in a situation where the particle is far from the wall. In the limit of very slender prolate particles, as
$e\rightarrow 1$
, we find
$$\begin{eqnarray}{\it\beta}^{\ast }(e\rightarrow 1)=\cos ^{-1}\sqrt{{\textstyle \frac{3}{11}}(5-\sqrt{3})}\approx 19.25^{\circ }.\end{eqnarray}$$
For walls inclined more steeply than this angle, a prolate body of any eccentricity cannot escape. The drag anisotropy in the limit
$e\rightarrow 1$
is not nearly as significant in the oblate case as it is in the prolate case (see Happel & Brenner Reference Happel and Brenner1983), which implies that escape from the wall is more difficult; in fact a wall inclination greater than 11.48° is sufficient to prevent escape for all oblate bodies.
6. Discussion
In this paper we studied the sedimentation of rigid prolate and oblate spheroids in a highly viscous fluid near a vertical or tilted wall. A system of ordinary differential equations governing the fully three-dimensional trajectories was derived. In numerous special cases, the system of equations yielded approximate analytical results for particle trajectories. The analytical predictions were compared with the results of full numerical simulations of the Stokes equations using a novel double-layer boundary integral scheme, the method of stresslet images. These two approaches were used to investigate a wide array of trajectory types for bodies of arbitrary eccentricity, and near a vertical or inclined wall. When the wall is vertical, a nearly spherical body may undergo a periodic tumbling motion. In three-dimensional versions of the tumbling trajectory, a periodic lateral wobbling arises. For more eccentric particles, three-dimensional glancing and reversing trajectories appear, with the body approaching the wall only once before receding back into the bulk fluid. When the wall is tilted, the symmetry in the system is weakened. As a consequence, new trajectory types appear, while the periodic tumbling orbit vanishes. Glancing- and reversing-type behaviour is still possible, but a sliding trajectory emerges for many combinations of particle eccentricity and wall inclination angle. The sliding trajectory was found to be asymptotically stable to small translational and rotational perturbations in the far-field hydrodynamic theory. Critical wall inclination angles distinguishing sliding from either escaping or colliding with the wall were also presented.
Improvement of the analytical predictions given in this paper might be challenging. For instance, the inclusion of lubrication effects would be beneficial for understanding the near-wall interactions but would require other techniques similar to those employed for sphere–sphere interactions by Durlofsky, Brady & Bossis (Reference Durlofsky, Brady and Bossis1987) (see also Brady & Bossis Reference Brady and Bossis1988). At the same time, the techniques we used here could be extended with no conceptual (but perhaps some algebraic) difficulty to deal with general triaxial ellipsoids. The mobility problem for imposed torques can also be solved in a similar manner, which could be used to obtain the solution of the general resistance problem to the same level of accuracy. Body deformability, multiple-body interactions and the inclusion of a background shear flow may be considered in future work.
Acknowledgements
We are grateful to M. Graham, J. Meiss and J.-L. Thiffeault for helpful comments. We also acknowledge the authors of the freely available 3D and 2D visualization software packages Mayavi and Matplotlib (Hunter Reference Hunter2007; Ramachandran & Varoquaux Reference Ramachandran and Varoquaux2011). This research was supported in part by NSF grant DMS 1147523.
Supplementary data
Supplementary data is available at http://dx.doi.org/10.1017/jfm.2015.222.
Appendix A. Approximate sedimentation dynamics of spheroids in the general setting
The solution of the mobility problem for a sedimenting prolate or oblate body near a wall with inclination angle
${\it\beta}$
is discussed in § 4. The components of the translational velocity,
$\boldsymbol{U}=U_{x}\hat{\boldsymbol{x}}+U_{y}\hat{\boldsymbol{y}}+U_{z}\hat{\boldsymbol{z}}$
, are given by
$$\begin{eqnarray}\displaystyle U_{x} & = & \displaystyle \frac{(2\cos {\it\beta}-(1\pm \cos (2{\it\theta}))\cos {\it\beta}\cos ^{2}{\it\phi}\pm \cos {\it\phi}\sin {\it\beta}\sin (2{\it\theta}))(X^{A}-Y^{A})+2Y^{A}\cos {\it\beta}}{2X^{A}Y^{A}}\nonumber\\ \displaystyle & & \displaystyle -\,\frac{9\cos {\it\beta}}{16h}+\frac{1}{128h^{3}}\bigg[4e^{2}\cos {\it\phi}\sin {\it\beta}\sin (2{\it\theta})\nonumber\\ \displaystyle & & \displaystyle +\,(2e^{2}(\cos (2\,{\it\theta})\pm 1)\cos ^{2}{\it\phi}+18\,e^{2}\cos ^{2}{\it\theta}-(17\pm 7)e^{2}+16)\cos {\it\beta}\bigg],\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle U_{y} & = & \displaystyle \frac{\sin {\it\phi}(\pm \!\sin {\it\beta}\sin (2\,{\it\theta})-(1\pm \cos (2\,{\it\theta}))\cos {\it\beta}\cos {\it\phi})(X^{A}-Y^{A})}{2X^{A}Y^{A}}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{e^{2}\sin {\it\phi}(2\sin {\it\beta}\sin (2\,{\it\theta})+(\cos (2{\it\theta})\pm 1)\cos {\it\beta}\cos {\it\phi})}{64h^{3}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle U_{z} & = & \displaystyle \frac{9\sin {\it\beta}}{8h}-\frac{2Y^{A}\sin {\it\beta}\pm (\cos {\it\beta}\cos {\it\phi}\sin (2{\it\theta})+(\cos (2{\it\theta})\pm 1)\sin {\it\beta})(X^{A}-Y^{A})}{2X^{A}Y^{A}}\nonumber\\ \displaystyle & & \displaystyle -\,\frac{e^{2}\cos {\it\beta}\cos {\it\phi}\sin (2{\it\theta})-(14e^{2}\sin ^{2}{\it\theta}+(1\pm 5)e^{2}-16)\sin {\it\beta}}{32h^{3}}.\end{eqnarray}$$
The three components of the rotational velocity,
${\it\bf\Omega}={\it\Omega}_{x}\hat{\boldsymbol{x}}+{\it\Omega}_{y}\hat{\boldsymbol{y}}+{\it\Omega}_{z}\hat{\boldsymbol{z}}$
, are given by
$$\begin{eqnarray}\displaystyle {\it\Omega}_{x} & = & \displaystyle \frac{9e^{2}\sin {\it\phi}\left((1\pm 1-2\sin ^{2}{\it\theta})\cos {\it\beta}\cos {\it\phi}-3\sin {\it\beta}\sin (2{\it\theta})\right)}{64(2-e^{2})h^{2}}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{6\cos {\it\beta}\cos {\it\phi}\left(6e^{2}\sin ^{4}{\it\theta}-8\sin ^{2}{\it\theta}+(1\pm 1)(4-e^{2}-2e^{2}\sin ^{2}{\it\theta})\right)}{128(e^{2}-2)h^{4}}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{3e^{2}\sin {\it\beta}\sin (2{\it\theta})\sin {\it\phi}\left(12e^{2}\sin ^{2}{\it\theta}+(3\pm 5)e^{2}-18\right)}{128(e^{2}-2)h^{4}},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\it\Omega}_{y} & = & \displaystyle \frac{27e^{2}\sin {\it\beta}\cos {\it\phi}\sin (2{\it\theta})+9e^{2}\cos {\it\beta}\left(2-4\cos ^{2}{\it\theta}+(2\cos ^{2}{\it\theta}\pm 1-1)\sin ^{2}{\it\phi}\right)}{64(2-e^{2})h^{2}}\nonumber\\ \displaystyle & & \displaystyle -\,\frac{3e^{2}\sin {\it\beta}\cos {\it\phi}\sin (2{\it\theta})}{128(2-e^{2})h^{4}}(12e^{2}\cos ^{2}{\it\theta}-5e^{2}(3\pm 1)+18)-\frac{6\cos {\it\beta}}{128(2-e^{2})h^{4}}\nonumber\\ \displaystyle & & \displaystyle \times \,\Big[e^{4}\cos ^{4}{\it\theta}(6\sin ^{2}{\it\phi}-9)+2e^{2}\cos ^{2}{\it\theta}\sin ^{2}{\it\phi}(4\pm e^{2}-5e^{2})+e^{2}\sin ^{2}{\it\phi}(4-3e^{2})\nonumber\\ \displaystyle & & \displaystyle \times \,(\pm 1-1)-e^{2}\cos ^{2}{\it\theta}(12-(15\pm 2)e^{2})-(7\pm 1)e^{4}+10e^{2}-4\Big],\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\it\Omega}_{z} & = & \displaystyle \cos {\it\beta}\sin {\it\phi}\sin (2{\it\theta})\left[\frac{-9e^{2}}{64(2-e^{2})h^{2}}+\frac{3e^{2}(6e^{2}\cos ^{2}{\it\theta}-e^{2}(7\pm 5)+8)}{256(2-e^{2})h^{4}}\right].\end{eqnarray}$$
$\pm$
signs should be replaced with
$+$
in the prolate case and
$-$
in the oblate case. The constants
$X^{A}$
and
$Y^{A}$
also have different definitions in these two cases, as indicated in table 1. The errors in the expressions above are of size
$O(h^{-4})$
in the translational velocity and
$O(h^{-5})$
in the rotational velocity for
$h\gg 1$
(verified by comparison with full numerical simulations in figure 9 for a test problem with
${\it\theta}={\rm\pi}/5$
,
${\it\phi}={\rm\pi}/7$
,
${\it\beta}={\rm\pi}/100$
and
$e=\sqrt{3}/2$
).
Figure 9. The boundary integral method and the differential equations agree with the expected rate of convergence in the distance from the wall, with errors scaling as
$O(h^{-4})$
in the translational velocity and as
$O(h^{-5})$
in the rotational velocity. Using
${\it\theta}={\rm\pi}/5$
,
${\it\phi}={\rm\pi}/7$
,
${\it\beta}={\rm\pi}/100$
and
$e=\sqrt{3}/2$
, the differences between the computed and analytical translational and rotational velocity vectors in the
$\Vert \cdot \Vert _{\infty }$
norm are shown. Both prolate and oblate bodies are considered, using 1972 quadrature nodes for the prolate body and 2281 for the oblate body. This may be viewed either as a validation of the full numerical scheme or as a check of the differential equations.
The time derivatives of
${\it\phi}$
and
${\it\theta}$
can be obtained from
${\it\bf\Omega}$
through the relations
$$\begin{eqnarray}\displaystyle & \dot{{\it\theta}}=-{\it\Omega}_{y}\cos ({\it\phi})+{\it\Omega}_{x}\sin ({\it\phi}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \dot{{\it\phi}}=\left\{\begin{array}{@{}ll@{}}{\it\Omega}_{z}-\tan ({\it\theta})[{\it\Omega}_{x}\cos ({\it\phi})+{\it\Omega}_{y}\sin ({\it\phi})]\quad & (\text{prolate}),\\ {\it\Omega}_{z}+\cot ({\it\theta})[{\it\Omega}_{x}\cos ({\it\phi})+{\it\Omega}_{y}\sin ({\it\phi})]\quad & (\text{oblate}).\end{array}\right. & \displaystyle\end{eqnarray}$$
${\it\phi}$
is indeterminate (
${\it\theta}=0$
for oblate bodies and
${\it\theta}={\rm\pi}/2$
for prolate bodies), these formulae require the choice
${\it\phi}=0$
; this degeneracy requires no extra bookkeeping in our problem since particles must already have
${\it\phi}=0$
when passing through these indeterminate positions. Equations (A 7) and (A 8) can then be simplified to give 
$$\begin{eqnarray}\displaystyle \dot{{\it\theta}} & = & \displaystyle \frac{18e^{2}\cos {\it\beta}\cos {\it\phi}\cos (2{\it\theta})-27e^{2}\sin {\it\beta}\sin (2{\it\theta})}{64(2-e^{2})h^{2}}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{3}{256\,(2-e^{2})h^{4}}\Big[4e^{2}\cos {\it\theta}\sin {\it\beta}\sin {\it\theta}(18-(9\pm 5)e^{2}+6\,e^{2}\cos (2\,{\it\theta}))\nonumber\\ \displaystyle & & \displaystyle -\,\cos {\it\beta}\cos {\it\phi}(16-16e^{2}+7e^{4}+e^{2}\cos (2{\it\theta})(24-(12\pm 4)e^{2}+9\,e^{2}\cos (2\,{\it\theta})))\Big],\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \dot{{\it\phi}} & = & \displaystyle \left\{\begin{array}{@{}ll@{}}\displaystyle \frac{3\cos {\it\beta}\sin {\it\phi}\tan {\it\theta}}{64(2-e^{2})}\left(\frac{-6e^{2}}{h^{2}}+\frac{3e^{4}\cos ^{2}{\it\theta}-8e^{4}+10e^{2}-4}{h^{4}}\right)\quad & (\text{prolate}),\\ \displaystyle \frac{3\cos {\it\beta}\sin {\it\phi}\cot {\it\theta}}{64(2-e^{2})}\left(\frac{-6e^{2}}{h^{2}}-\frac{3e^{4}\sin ^{2}{\it\theta}-2e^{4}-2e^{2}-4}{h^{4}}\right)\quad & (\text{oblate}),\end{array}\right.\end{eqnarray}$$
$(x,y,h,{\it\theta},{\it\phi})$
, as discussed in § 4.
Appendix B. Verification and validation of the numerical method
In this section we give additional details on the accuracy of the method of stresslet images described in § 2.2. This includes a convergence study and direct comparisons with previously published results for an inclined oblate body and exact solutions for the motion of a sphere.
We begin by describing some checks on the time-stepping trajectory computations in this work. For the trajectories with lateral symmetry of figures 2 and 5, the computed out-of-plane motion provides a simple quantitative estimate of the accumulated error; for these trajectories we rejected results for which the ratio of out-of-plane to in-plane translations or rotations exceeded
$10^{-4}$
in any time step. For the three-dimensional trajectories of figures 3 and 7 this simple check is not available, so we reversed the direction of gravity after each computation and evolved back to the original position. In all of the results presented here the initial position was recovered with a relative error of less than
$10^{-3}$
in the
$\Vert \cdot \Vert _{\infty }$
-norm. Typical trajectory runs used
${\sim}500$
nodes and
${\sim}250$
time steps, i.e. required the inversion of 500 dense matrices of size 1500.
B.1. Convergence tests for an oblate body near a wall
We now benchmark our numerical method against the results obtained by Hsu & Ganatos (Reference Hsu and Ganatos1989). For this test we consider an oblate body with unit radius and aspect ratio 2 or 10 (in our notation
$e=\sqrt{3}/2$
or
$e=3\sqrt{11}/10$
), and with position and orientation described by the parameters
$h=1.1$
or
$h=1.5$
,
${\it\theta}=75^{\circ }$
and
${\it\phi}=0$
. We impose a wall-normal velocity
$\boldsymbol{U}=(0,0,1)$
and zero rotation and solve the resistance problem, reporting the
$z$
-component of drag normalized by the value that would prevail in the absence of a wall. These normalizations are known exactly, 15.084358 for
$a/c=2$
and 11.862466 for
$a/c=10$
. The results are indicated for various discretization levels
$N$
in table 2 together with the values computed by Hsu & Ganatos (Reference Hsu and Ganatos1989). The numerical method of that work was optimized to treat the case of an axisymmetric body and the quoted results were accordingly obtained more cheaply, with
$N\leqslant 100$
. In three of the four cases our results agree, but with
$h=1.1$
and
$a/c=2$
there is a small difference which is probably due to insufficient grid resolution in the previously published work.
Table 2. Convergence tests for a resistance problem with an inclined oblate spheroid of aspect ratio 2 or 10 and with centroid at height 1.1 or 1.5, with
$N$
the number of nodes used to discretize the body surface.
Table 3. Drag on a sphere of radius 0.1 translating parallel to a wall at six distances using four discretizations, normalized by the value predicted by Stokes’ law for an unbounded fluid. The shaded values were found using the method of stresslet images, as described in § 2.2. These are accompanied by results obtained via regularized Stokeslets, quoted from table 1 in Ainley et al. (Reference Ainley, Durkin, Embid, Boindala and Cortez2008).
B.2. Comparison with exact solutions and regularized Stokeslets for a sphere near a wall
The resistance problem for a spherical body translating without rotation in a fluid bounded by a plane wall was solved exactly using bispherical coordinates by O’Neill (Reference O’Neill1964) for motion parallel to the wall and by Brenner (Reference Brenner1961) for motion normal to it; see Goldman et al. (Reference Goldman, Cox and Brenner1967a ) for a summary of these results. More recently, Ainley et al. (Reference Ainley, Durkin, Embid, Boindala and Cortez2008) solved this problem numerically using the method of regularized Stokeslets. In this section we tabulate the results of the method of stresslet images against these exact and numerical solutions.
The geometry of the problem is determined by two parameters, the sphere radius
$a$
and the bispherical parameter
${\it\alpha}$
, which satisfies
$\cosh ({\it\alpha})=d/a$
, where
$d$
is the distance from the particle centre to the wall. The gap size
$d$
–
$a$
is the distance from the particle surface to the wall. Following Ainley et al., we consider
$a=0.1$
and
${\it\alpha}=10,3,2,1,0.5,0.3$
. The integrals over the sphere are discretized using
$N$
-point quadrature rules, for
$N=468,812,1486,2718$
. We then calculate the drag
$F$
and torque
$T$
when the sphere translates at speed 1 with no rotation; these values are non-dimensionalized by the drag and torque that would occur in the absence of the wall.
In table 3 we give the component of drag in the same direction as the translation when the particle moves parallel to the wall, normalized by the drag predicted by Stokes’ law for an unbounded fluid. Table 4 gives this drag correction factor when the direction of translation is normal to the wall. In both cases, the method of stresslet images gives more accurate results than the method of regularized Stokeslets at large and moderate gap sizes, and without the need of a regularization parameter; at the smallest gap size the performance of the two methods is similar. The exact solutions were recalculated following equation (2.19) in the work of Brenner (Reference Brenner1961) and equation (26) in the work of O’Neill (Reference O’Neill1964). These values appear in the rightmost columns of tables 3 and 4.
Table 4. Drag on a sphere of radius 0.1 translating perpendicular to a wall at six distances using four discretizations, normalized by the value predicted by Stokes’ law for an unbounded fluid. The shaded values were found using the method of stresslet images, as described in § 2.2. These are accompanied by results obtained via regularized Stokeslets, quoted from table 2 in Ainley et al. (Reference Ainley, Durkin, Embid, Boindala and Cortez2008).
