# Mass transfer to freely suspended particles at high Péclet number

John M. Lawson1<EMAIL_ADDRESS>1Aerodynamics and Flight Mechanics
Research Group, University of Southampton, SO17 1BJ, UK

###### Abstract

In a theoretical analysis, we generalise well known asymptotic results to
obtain expressions for the rate of transfer of material from the surface of an
arbitrary, rigid particle suspended in an open pathline flow at large Péclet
number, Pe. The flow may be steady or periodic in time. We apply this result
to numerically evaluate expressions for the surface flux to a freely
suspended, axisymmetric ellipsoid (spheroid) in Stokes flow driven by a steady
linear shear. We complement these analytical predictions with numerical
simulations conducted over a range of $\textrm{Pe}=10^{1}-10^{4}$ and confirm
good agreement at large Péclet number. Our results allow us to examine the
influence of particle shape upon the surface flux for a broad class of flows.
When the background flow is irrotational, the surface flux is steady and is
prescribed by three parameters only: the Péclet number, the particle aspect
ratio and the strain topology. We observe that slender prolate spheroids tend
to experience a higher surface flux compared to oblate spheroids with
equivalent surface area. For rotational flows, particles may begin to spin or
tumble, which may suppress or augment the convective transfer due to a
realignment of the particle with respect to the strain field.

###### keywords:

##  1 Introduction

When small, rigid particles are immersed in a fluid, material (e.g. a solute)
may be transferred away from the surface by convection and diffusion. We shall
refer to this process as mass transfer, although an analogous problem exists
for heat transfer. The engineered and natural world are replete with examples:
planktonic osmotrophs absorb dissolved nutrients (Karp-Boss et al., 1996;
Guasto et al., 2012; Pahlow et al., 1997), bacterial hosts encounter viruses
(Guasto et al., 2012), crystals are grown in agitated suspension to produce
pharmaceuticals and industrial products (Myerson, 2002) and microplastics
leach and sorb pollutants in ocean waters (Law, 2017; Suhrhoff & Scholz-
Böttcher, 2016; Seidensticker et al., 2017). The examples cited above have in
common that the solid particle exchanging material is rarely ever spherical,
is usually small compared to the flow features in which it is embedded, and
the material diffuses slowly such that convection is the dominant mechanism of
mass transfer. The question then naturally follows: what is the rate of mass
transfer from the surface?

This question belongs to a general set of classical problems which have been
well studied (see e.g. works summarised by Clift et al. (1978); Leal (2012);
Michaelides (2003)). The answer depends upon the competing effects of
convection and diffusion, as parametrised by the Péclet number, Pe. When
$\textrm{Pe}\ll 1$, conduction dominates inside an inner region near the
particle and the effects of shape are readily accounted for (Batchelor, 1979;
Acrivos, 1980; Feng & Michaelides, 1997; Pozrikidis, 1997). When
$\textrm{Pe}\gg 1$, convection dominates the mass transfer process and the
surface flux of solute is determined by the particle geometry and the nature
of the relative flow field, parametrised by the Reynolds number Re. In closed
streamline flows, the solute simply recirculates around the particle and the
surface flux approaches a constant as $\textrm{Pe}\rightarrow\infty$ (Frankel
& Acrivos, 1968; Poe & Acrivos, 1976). When the particle is surrounded by a
region of open streamlines (or pathlines), the solution to the convection-
diffusion problem takes the form of a thin concentration boundary layer around
the particle. For this case, Acrivos (1960) and later Goddard & Acrivos (1965)
introduced asymptotic methods to compute the mass transfer rate, which scales
as $\textrm{Pe}^{1/3}$. The flow and the geometry of the problem then
determine the scaling coefficient, which prescribes the surface flux.

For sufficiently small particles, $\textrm{Re}\ll 1$ and the surrounding
relative flow field may be well approximated by a Stokes flow consisting of a
background uniform flow or linear shear, plus a perturbation owing to the
presence of the particle. The available analytical results in the high Péclet
number, low Reynolds number regime are generally limited to simplified flows
or geometries. A general solution to this problem would have great utility,
because solid particles are often neither spherical nor subject to motions as
simple as uniform or axisymmetric flows (Leal, 2012). Specific analytical
results are available for spheres in uniform flow (Acrivos, 1960; Acrivos &
Taylor, 1962) or arbitrary linear shear (Gupalo & Riazantsev, 1972; Poe &
Acrivos, 1976; Batchelor, 1979) and axisymmetric bodies in uniform flow
(Sehlin, 1969; Gupalo et al., 1976; Leal, 2012; Dehdashti & Masoud, 2020). To
our knowledge, equivalent results are not available for arbitrary bodies with
high Péclet numbers in an arbitrary linear shear. Experimental and numerical
studies have focused on the uniform flow case typically with
$\textrm{Pe}=O(\textrm{Re})$ (Masliyah & Epstein, 1972; Clift et al., 1978;
Sparrow et al., 2004; Ke et al., 2018; Kishore & Gu, 2011; Ma & Zhao, 2020);
relatively few studies have examined the high Péclet number, low Reynolds
number regime numerically (Karp-Boss et al., 1996; Pahlow et al., 1997) and
the results available for linear shear flows are limited to a handful of
cases. Therefore, additional empirical data are needed for the general case of
arbitrary linear shear to support a generalisation of asymptotic results.

For freely suspended particles, an additional complication arises which
affects the mass transfer rate. In the absence of body forces, such as the
case of neutrally buoyant particles, the drag and consequently the slip
velocity vanish, such that convection is provided by the linear shear alone.
However, the fluid may exert a couple upon the particle (Jeffery, 1922), which
in the absence of a restoring torque will result in the steady precession or
rotation of the particle about its axis. As a result, the flow field around
the particle and resultant surface flux is in general unsteady (Pahlow et al.,
1997). However, for spherical particles, Batchelor (1979) and Batchelor (1980)
argued that the unsteady contribution to the scalar flux may be neglected at
large Péclet number, so that the average flow field perceived by the body
determines the average mass transfer rate.

In this article, we extend the work of Goddard & Acrivos (1965) and Batchelor
(1979) to consider the mass transfer rate due to convection and diffusion from
an arbitrary body in an unsteady, open pathline flow at high Péclet number. We
then apply this analysis to determine the mass transfer rate from a neutrally
buoyant spheroid freely suspended in an arbitrary linear shear. We obtain
asymptotic expressions for the transfer rate in the resultant average flow
field and compare these to the results of numerical simulations, which provide
a quantitative test of the accuracy of the asymptotic approximation. Our
results allow us to examine the role of particle shape in the mass transfer
process for a very broad class of flows and open up new possibilities in the
numerical simulation of mass transfer in particle suspensions.

The paper is structured as follows. In §2, we review the theoretical
background of the problem and derive a general form for the mass transfer
coefficient for an arbitrary particle in an unsteady, time periodic flow. In
§3, we apply this general expression to the case of a freely suspended
spheroid in Stokes flow. In §4, we introduce a numerical test of the results
given in §3 and discuss the influence of particle shape upon the surface flux.
We present a summary of our results and future perspectives in §5.

##  2 The steady flux at large Péclet number

In this section, we shall extend the analyses of Goddard & Acrivos (1965) and
Batchelor (1979) to derive a general expression for the average solute flux
from the surface of a particle of arbitrary shape in a steady, open-streamline
flow. We begin by introducing the governing equations in dimensionless form,
then examine the case of steady flow. Finally, we generalise the result to
unsteady flow.

###  2.1 Governing equations

The mass transport from the surface of the particle is governed by the
convection-diffusion equation. This may be written in dimensionless form as

$\frac{\partial\theta}{\partial
t}+\bm{u}\cdot\nabla\theta=\frac{1}{\textrm{Pe}}\nabla^{2}\theta$ (1)

where $\bm{u}(\bm{x},t)$ is the velocity field and $\theta(\bm{x},t)$ is the
concentration of the solute (Leal, 2012). The characteristic lengthscale is
$r$, the linear dimension of the particle, so that the spatial coordinate is
non-dimensionalised as $\bm{x}=\bm{x}^{*}/r$. The characteristic timescale is
prescribed by the shear rate $E^{*}$, so that time is non-dimensionalised as
$t=t^{*}E^{*}$. The Péclet number is defined $\textrm{Pe}\equiv
r^{2}E^{*}/\kappa$, where $\kappa$ is the diffusion coefficient of the solute.
Our convention is to write dimensional quantities with a superscript ∗ unless
otherwise stated.

We shall adopt a frame of reference moving with the particle, such that the
boundary of the particle $S_{p}$ is stationary. We impose the boundary
condition

$\begin{matrix}\theta=1\text{ and }\bm{u}=0&\text{for }&\bm{x}\in S_{p}\\\
\theta\rightarrow 0&\text{as }&|\bm{x}|\rightarrow\infty\end{matrix}$ (2)

so that the concentration of solute at the particle surface is uniform. This
boundary condition corresponds to non-dimensionalising the concentration field
$C(\bm{x}^{*},t^{*})$ as $\theta=(C-C_{0})/(C_{1}-C_{0})$, where $C_{1}$ and
$C_{0}$ are the concentration of solute at the surface and infinity in
physical units.

The non-dimensional measure of the surface flux $Q$ is the Sherwood number

$\textrm{Sh}=\frac{Q}{4\pi
r\kappa(C_{1}-C_{0})}=-\frac{1}{4\pi}\iint_{S_{p}}\nabla\theta\cdot\mathrm{d}\mathbf{S}.$
(3)

The denominator in (3) corresponds to the steady state flux from a sphere of
radius $r$ in the absence of flow. A convenient choice for the lengthscale $r$
is $\sqrt{A/4\pi}$, where $A$ is the surface area of the particle. The
Sherwood number is therefore the factor by which convection enhances the mass
transfer relative to the diffusive flux from a spherical particle with the
same surface area.

###  2.2 General solution of the surface flux in steady flow

In steady flow, the scalar transport equation reduces to

$\bm{u}\cdot\nabla\theta=\frac{1}{\textrm{Pe}}\nabla^{2}\theta$ (4)

We seek an asymptotic solution for (4) subject to (2) at large Péclet number.

At large Péclet number, the concentration of solute is small almost everywhere
far from the particle, apart from a thin wake, along which material from the
surface is swept. Near the particle, there is a concentration boundary layer
whose thickness is $O(\delta)$, across which there is a sharp jump in
concentration. For the boundary layer approximation to hold, the pathlines
near the surface of the particle must originate and terminate at infinity,
such that the surface of the particle is exposed to a constant stream of fresh
fluid. Recirculating regions (closed pathlines), which can occur in some
geometries, are not permitted.

We shall construct a general curvilinear coordinate system $(\xi,\eta,\zeta)$
to describe the flow near the surface in terms of the distance from the
surface and the direction of the flow near the surface. This coordinate system
is not necessarily orthogonal and we will find it useful to describe it in
terms of the covariant coordinate vectors

$\bm{h}_{\xi}=\frac{\partial\bm{x}}{\partial\xi}~{},~{}\bm{h}_{\eta}=\frac{\partial\bm{x}}{\partial\eta}~{},~{}\bm{h}_{\zeta}=\frac{\partial\bm{x}}{\partial\zeta}$
(5)

and their pathlines, which we have illustrated in figure 1. We note that the
adoption of a general curvilinear coordinate system, rather than an orthogonal
one, is crucial in the generalisation to an arbitrary flow or body.

Figure 1: Illustration of the curvilinear coordinate system $(\xi,\eta,\zeta)$
defined on the surface ($\xi=0$) of a spheroid in an arbitrary linear shear.
Thick black lines show $\eta$-coordinate pathlines
($\xi=0,\zeta=\textrm{const.}$), which are tangent to the local viscous shear
stress on the surface. Thin grey lines are $\zeta$-coordinate pathlines
($\xi=0,\eta=\textrm{const.}$). Three fixed points of the surface streamlines
are shown: red (source), green (saddle) and blue (sink).

The coordinate $\xi$ is defined as the normal distance from the particle
surface, so that at the particle surface, $\bm{h}_{\xi}$ is a unit vector
normal to the surface. The $\eta$ and $\zeta$ coordinates lie tangent to the
surface. We shall define the $\eta$ coordinate so that at a small distance
$\xi$ above the surface, the component of fluid velocity which is tangent to
the surface is parallel to $\bm{h}_{\eta}$. The curves along the $\eta$
direction are therefore “surface streamlines” which are tangent to the
direction of the surface velocity gradient
$\bm{w}=\partial\bm{u}/\partial\xi|_{\xi=0}$.

We can therefore express the velocity near the surface as

$\displaystyle\bm{u}$
$\displaystyle=u^{\xi}\bm{h}_{\xi}+u^{\eta}\bm{h}_{\eta}+u^{\zeta}\bm{h}_{\zeta}$
(6) $\displaystyle=\bm{w}\xi+O(\xi^{2}).$

so that near the particle surface, the velocity components are of the form

$u^{\xi}=G(\eta,\zeta)\xi^{2}+O(\xi^{3}),u^{\eta}=\xi
F(\eta,\zeta)+O(\xi^{2}),u_{\zeta}=O(\xi^{2})$ (7)

By definition, the velocity at the surface is zero; the tangential component
$u^{\eta}$ is linear in $\xi$ to leading order, whereas the surface-normal
component is quadratic (Leal, 2012). The function $F(\eta,\zeta)$ is therefore
obtained in terms of the velocity gradient at the surface of the particle as

$F=\bm{w}\cdot\nabla\eta$ (8)

and is related to $G(\eta,\zeta)$ through the continuity equation.

To obtain the surface normal velocity component $u^{\xi}$, we express the
continuity equation as (Grinfeld, 2013)

$\rho\nabla_{i}u^{i}=\frac{\partial(\rho
u^{\xi})}{\partial\xi}+\frac{\partial(\rho u^{\eta})}{\partial\eta}=0$ (9)

where $\rho=\rho(\eta,\zeta)\equiv\sqrt{\det{g}}$ and $\det{g}$ is the
determinant of the metric tensor $g_{ij}\equiv\bm{h}_{i}\cdot\bm{h}_{j}$.
Since $g_{\xi\xi}=1$, $g_{\xi\eta}=g_{\xi\zeta}=0$ by definition, the physical
interpretation of $\rho$ is a local area density of the surface streamlines,
such that $\delta A=\rho\delta\eta\delta\zeta$ is the area covered by a small
region on the surface measuring $\delta\eta$ along a surface streamline and
$\delta\zeta$ across adjacent streamlines. We can then define a streamfunction
$\psi$ as

$u^{\xi}=-\frac{1}{\rho}\frac{\partial\psi}{\partial\eta}~{},~{}u^{\eta}=\frac{1}{\rho}\frac{\partial\psi}{\partial\xi}$
(10)

so that

$\psi=\frac{1}{2}\xi^{2}\rho F$ (11)

is a solution which satisfies (7) and the continuity equation (9) to leading
order in $\xi$.

Writing (4) in these curvilinear coordinates, our solution now proceeds
analogously to the works of Goddard & Acrivos (1965) and Batchelor (1979). We
shall formalise Batchelor’s analysis here in our new curvilinear coordinate
system for the purpose of exposition. To leading order in $\xi$, the
convection term is

$u^{i}\nabla_{i}\theta=-\frac{1}{\rho}\frac{\partial\psi}{\partial\eta}\frac{\partial\theta}{\partial\xi}+\frac{1}{\rho}\frac{\partial\psi}{\partial\xi}\frac{\partial\theta}{\partial\eta}+O(\xi^{2})$
(12)

whereas the diffusion term may be written

$\nabla^{2}\theta=\frac{1}{\rho}\left[\frac{\partial}{\partial\xi^{i}}\left(\rho
g^{ij}\frac{\partial\theta}{\partial\xi^{j}}\right)\right]=\frac{1}{\rho}\frac{\partial}{\partial\xi}\left(\rho\frac{\partial\theta}{\partial\xi}\right)+\ldots$
(13)

where $g^{ij}$ is the inverse metric tensor and $\xi^{i}$ indexes the
coordinates $(\xi,\eta,\zeta)$. The terms omitted in the expansion of (13) do
not involve any gradients in $\xi$. By neglecting them, we assume that the
diffusive flux of material across the surface is small in comparison to that
normal to the surface. This essentially requires that the particle surface be
smooth and contain no regions of extreme curvature where this assumption may
break down.

We now rescale our coordinate system and streamfunction so that the surface-
normal concentration gradient is $O(1)$. Rescaling $\xi$ by the boundary layer
thickness $\delta$, we write

$\Xi=\frac{\xi r}{\delta}=\xi\textrm{Pe}^{m}~{},~{}\Psi=\frac{1}{2}\Xi^{2}\rho
F=\textrm{Pe}^{2m}\psi$ (14)

where we suppose $\delta/r$ scales as $\textrm{Pe}^{-m}$. We can rewrite (4)
using (12) and (13)

$-\frac{\partial\Psi}{\partial\eta}\frac{\partial\theta}{\partial\Xi}+\frac{\partial\Psi}{\partial\Xi}\frac{\partial\theta}{\partial\eta}+O(\textrm{Pe}^{-m})=\textrm{Pe}^{3m-1}\frac{\partial}{\partial\Xi}\left(\rho\frac{\partial\theta}{\partial\Xi}\right)$
(15)

The first two terms in (15) are $O(1)$ and thus $m=1/3$, as expected.

We can now solve (15) with the well-known von-Mises transformation. Adopting
the change of variables $(\Xi,\eta,\zeta)\rightarrow(\Psi,\eta,\zeta)$, we
have

$\frac{\partial\theta}{\partial\eta}=\frac{\partial}{\partial\Psi}\left(\rho\frac{\partial\theta}{\partial\Psi}\frac{\partial\Psi}{\partial\Xi}\right)=\frac{\partial}{\partial\Psi}\left(\frac{\partial\theta}{\partial\Psi}(2\rho^{3}F\Psi)^{1/2}\right).$
(16)

Equation (16) now admits a self-similar solution of the form
$\theta=\theta(\chi),\chi=\Psi^{1/2}/\tau^{1/3}$, where the functions
$\theta(\chi)$ and $\tau(\eta)$ must satisfy

$\frac{\mathrm{d}^{2}\theta}{\mathrm{d}\chi^{2}}+\frac{4}{3}\chi^{2}\frac{\mathrm{d}\theta}{\mathrm{d}\chi}=0$
(17)

for

$\frac{\mathrm{d}\tau}{\mathrm{d}\eta}=(2\rho^{3}F)^{1/2}.$ (18)

The solution of (17) satisfying the imposed boundary conditions (2) of
$\theta(0)=1$ and $\theta=0$ as $\chi\rightarrow\infty$ is

$\theta(\chi)=\frac{\Gamma(\frac{1}{3},\frac{4}{9}\chi^{3})}{\Gamma(\frac{1}{3},0)}$
(19)

where $\Gamma$ is the incomplete gamma function. We can obtain $\tau$ by
integrating (18) along the surface streamlines, which must begin and terminate
at critical points $F=0$ on the surface, since the surface shear stress is
continuous. We will therefore choose the constant of integration so that

$\tau(\eta)=\int_{\eta_{0}}^{\eta}(2\rho^{3}F)^{1/2}\mathrm{d}\eta$ (20)

has $\tau(\eta_{0})=0$ at the beginning of the surface streamline
$\eta=\eta_{0}$ and $\tau(\eta_{1})=\tau_{1}$ at the end $\eta=\eta_{1}$.

At the surface, the local flux of material per unit area is

$\displaystyle-\frac{\partial\theta}{\partial\xi}\Big{|}_{\xi=0}$
$\displaystyle=-\textrm{Pe}^{1/3}\frac{\partial\theta}{\partial\Xi}\Big{|}_{\Xi=0}=-\textrm{Pe}^{1/3}\frac{\partial\theta}{\partial\chi}\Big{|}_{\chi=0}\cdot\frac{\partial\chi}{\partial\Xi}\Big{|}_{\Xi=0}$
(21)
$\displaystyle=\textrm{Pe}^{1/3}\frac{12^{1/3}}{\Gamma(\frac{1}{3})}\cdot\frac{\left(\frac{1}{2}\rho
F\right)^{1/2}}{\tau^{1/3}}$

The Sherwood number (3) is therefore given by

Sh
$\displaystyle=-\frac{1}{4\pi}\iint_{S_{p}}\frac{\partial\theta}{\partial\xi}\Big{|}_{\xi=0}\rho\mathrm{d}\eta\mathrm{d}\zeta$
(22)
$\displaystyle=\frac{12^{1/3}\textrm{Pe}^{1/3}}{8\pi\Gamma(1/3)}\iint_{S_{p}}\frac{\left(2\rho^{3}F\right)^{1/2}}{\tau^{1/3}}\mathrm{d}\eta\mathrm{d}\zeta$

which we can integrate by parts by recognising that

$\int_{\eta_{0}}^{\eta_{1}}\frac{\left(2\rho^{3}F\right)^{1/2}}{\tau^{1/3}}\mathrm{d}\tau=\int_{\eta_{0}}^{\eta_{1}}\frac{\mathrm{d}\tau}{\mathrm{d}\eta}\tau^{-1/3}\mathrm{d}\eta=\frac{3}{2}\tau_{1}^{2/3}$
(23)

and thus

$\textrm{Sh}=\frac{0.808\textrm{Pe}^{1/3}}{4\pi}\int\left(\int_{\eta_{0}}^{\eta_{1}}\rho^{3/2}F^{1/2}\mathrm{d}\eta\right)^{2/3}\mathrm{d}\zeta+O(1)$
(24)

As expected, the limiting dimensionless flux scales as
$\alpha\textrm{Pe}^{1/3}$, with a prefactor $\alpha$ which can now be
explicitly computed in terms of the shape of the body and the surface velocity
gradient.

The main point of (24) is that we have generalised the results of Goddard &
Acrivos (1965) and (Batchelor, 1979) to any distribution of surface
streamlines on an arbitrary body, not just those which result in an orthogonal
coordinate system. The caveat remains that the pathlines around the body
should be open for the boundary layer approximation to hold and the boundary
should be smooth and contain no regions of extreme curvature. Furthermore, we
have a natural way of numerically constructing the local area density $\rho$,
as shown graphically in figure 1. By integrating (5) with respect to $\eta$ to
generate a set of surface streamlines over the body, we obtain lines of
constant $\zeta$ at intervals of varying $\eta$. With a sufficiently fine
discretisation, we can numerically approximate the local area density $\rho$
and evaluate the coefficient in (24) numerically. We shall demonstrate this
procedure later in §3.

###  2.3 Unsteady solution at large Péclet number

We now seek to generalise our steady solution to an unsteady, periodic flow.
Our argument is analogous to that proposed by Batchelor (1980), who examined
the scalar flux to spherical particles in turbulent flow.

We shall consider the case where the motion is periodic with period
$T=E^{*}T^{*}$. We seek the time average Sherwood number
$\overline{\vphantom{\big{(}}\textrm{Sh}}$ at the surface, which depends upon
the time average concentration field $\overline{\vphantom{\big{(}}\theta}$.
The time average over a period $T$ is defined simply as

$\overline{\vphantom{\big{(}}\theta}(\bm{x})=\frac{1}{T}\int_{0}^{T}\theta(\bm{x},t)\mathrm{d}t.$
(25)

Applying this averaging procedure to (1), we obtain

$\overline{\vphantom{\big{(}}\bm{u}}\cdot\nabla\overline{\vphantom{\big{(}}\theta}+\overline{\vphantom{\big{(}}\bm{u}^{\prime}\cdot\nabla\theta^{\prime}}=\frac{1}{\textrm{Pe}}\nabla^{2}\overline{\vphantom{\big{(}}\theta}$
(26)

where we have decomposed the velocity field
$\bm{u}=\overline{\vphantom{\big{(}}\bm{u}}+\bm{u}^{\prime}$ and concentration
field $\theta=\overline{\vphantom{\big{(}}\theta}+\theta^{\prime}$ into their
time respective average and fluctuating components. Likewise, the transport
equation for the concentration fluctuation $\theta^{\prime}(\bm{x},t)$ is

$\frac{\partial\theta^{\prime}}{\partial
t}+\bm{u}\cdot\nabla\theta^{\prime}-\overline{\vphantom{\big{(}}\bm{u}^{\prime}\cdot\nabla\theta^{\prime}}-\frac{1}{\textrm{Pe}}\nabla^{2}\theta^{\prime}=\bm{u}^{\prime}\cdot\nabla\overline{\vphantom{\big{(}}\theta}$
(27)

We shall argue that the solution of (27) satisfying the boundary conditions
(2) and

$\theta^{\prime}=0\text{ at }\xi=0\text{ and }\theta^{\prime}=0\text{ as
}\xi\rightarrow\infty$ (28)

is such that $\theta^{\prime}\sim\textrm{Pe}^{-1/3}$ as
$\textrm{Pe}\rightarrow\infty$. As before, we reason that the solution
consists of a slender concentration boundary layer of thickness
$\delta/r\sim\textrm{Pe}^{-1/3}$. Outside the boundary layer $\xi\gg\delta$,
$\overline{\vphantom{\big{(}}\theta}\rightarrow 0$ and
$\theta^{\prime}\rightarrow 0$. Inside the boundary layer, the amplitude of
the concentration fluctuations are
$\theta^{\prime}=O(\textrm{Pe}^{m^{\prime}})$ and the jump in the mean
concentration across the boundary layer is $O(1)$. Clearly, $m^{\prime}\leq 0$
since there is no mechanism in (1) which can amplify scalar fluctuations. If
$m^{\prime}<0$, then scalar fluctuations should decay in amplitude at large
Pe, whereas if $m^{\prime}=0$ they remain comparable in magnitude to the mean
flow.

Let us examine the order of magnitude of the terms in (27). Since scalar
fluctuations are $O(\textrm{Pe}^{m^{\prime}})$ and occur on a timescale of
$O(T)$, the time derivative term scales as

$\frac{\partial\theta}{\partial t}\sim\frac{\textrm{Pe}^{m^{\prime}}}{T}.$
(29)

The convective terms on the left hand side scale as

$\bm{u}\cdot\nabla\theta^{\prime}=u_{\xi}\frac{\partial\theta^{\prime}}{\partial\xi}+\ldots\sim\left(\frac{\delta^{2}}{r^{2}}\right)\left(\frac{r}{\delta}\textrm{Pe}^{m^{\prime}}\right)\sim\textrm{Pe}^{m^{\prime}-1/3}$
(30)

whereas convective term on the right hand side scales as

$\bm{u}^{\prime}\cdot\nabla\overline{\vphantom{\big{(}}\theta}=u_{\xi}^{\prime}\frac{\partial\overline{\vphantom{\big{(}}\theta}}{\partial\xi}+\ldots\sim\left(\frac{\delta^{2}}{r^{2}}\right)\left(\frac{r}{\delta}\right)\sim\textrm{Pe}^{-1/3}.$
(31)

Finally, the diffusion term scales as

$\frac{1}{\textrm{Pe}}\nabla^{2}\theta^{\prime}\approx\frac{1}{\textrm{Pe}}\frac{\partial^{2}\theta^{\prime}}{\partial\xi^{2}}\sim\frac{1}{\textrm{Pe}}\frac{r^{2}}{\delta^{2}}\textrm{Pe}^{m^{\prime}}\sim\textrm{Pe}^{m^{\prime}-1/3}.$
(32)

Annotating (27) with the order of magnitude of each term, we have

$\underbrace{\vphantom{\Big{[}}\frac{\partial\theta^{\prime}}{\partial
t}}_{O(\textrm{Pe}^{m^{\prime}}/T)}\underbrace{\vphantom{\Big{[}}+\bm{u}\cdot\nabla\theta^{\prime}-\overline{\vphantom{\big{(}}\bm{u}^{\prime}\cdot\nabla\theta^{\prime}}-\frac{1}{\textrm{Pe}}\nabla^{2}\theta^{\prime}}_{O(\textrm{Pe}^{m^{\prime}-1/3})}=\underbrace{\vphantom{\Big{[}}\bm{u}^{\prime}\cdot\nabla\overline{\vphantom{\big{(}}\theta}}_{O(\textrm{Pe}^{-1/3})}$
(33)

We see that the solution requires $m^{\prime}=-1/3$, such that the first term
on the left hand side dominates and balances the unsteady convection term on
the right hand side, provided the dimensionless period
$T\ll\textrm{Pe}^{1/3}$. Thus, at large Péclet, we have

$\frac{\partial\theta^{\prime}}{\partial
t}\approx\bm{u}^{\prime}\cdot\nabla\overline{\vphantom{\big{(}}\theta}\sim\textrm{Pe}^{-1/3}$
(34)

which is equivalent to equation (3.9) of Batchelor (1980). However, we have
derived the result on more general grounds. Physically, this result means that
when the timescale of diffusion within the boundary layer becomes large in
comparison to the timescale of the unsteadiness, there is insufficient time
for diffusion to redistribute scalar fluctuations. Instead, scalar
fluctuations inside the boundary layer occur due to the convection of the mean
concentration field by unsteady velocity fluctuations.

Thus, to leading order, the concentration and its time mean are equal and the
mean scalar field is well approximated by

$\overline{\vphantom{\big{(}}\bm{u}}\cdot\nabla\overline{\vphantom{\big{(}}\theta}=\frac{1}{\textrm{Pe}}\nabla^{2}\overline{\vphantom{\big{(}}\theta}+O(\textrm{Pe}^{-2/3})$
(35)

This allows us to apply our general result (24) to unsteady, time periodic
flows where the period of motion is sufficiently small. As before, for the
boundary layer approximation to hold, it is required that the surface of the
particle be exposed to a continuous stream of fresh fluid. Therefore, it is
required that the pathlines of fluid parcels adjacent to the surface be open,
such that the fluid does not simply recirculate around a closed path.

##  3 The steady flux for a spheroid

We shall now apply the results of §2.2 to obtain the scaling coefficient
$\alpha$ for a spheroid in an arbitrary linear shear. Although the motion is
in general unsteady, we have argued in §2.2 that the average mass transfer
rate can be computed for an equivalent spheroid subject to the same average
relative flow field. This flow field can be described by three parameters, but
for a broad class of cases, the flow is reduced to an axisymmetric
configuration described by a single parameter. To proceed, we shall analyse
the relative motion of the particle in §3.1 and obtain expressions for the
average flow field perceived by the particle. We introduce an expression for
the velocity field near the particle in §3.2, then derive the mass transfer
coefficient in the axisymmetric and general cases in §3.3 and §3.4
respectively. Based on these results, we consider the special case of rotation
dominated flow in §3.5 and discuss potential extensions in §3.6.

###  3.1 Motion of a freely suspended spheroid

Let us consider the motion of a spheroidal particle in an arbitrary, linear
shear flow. The background linear shear is $\bm{v}=\bm{Gy}$, where
$\bm{y}=\bm{y}^{*}/r$ is the coordinate system of the fixed laboratory frame
and $\bm{G}=\bm{E}+\bm{W}$ is an arbitrary, steady velocity gradient in this
reference frame. The velocity gradient is composed of a symmetric strain
tensor $\bm{E}=\bm{E}^{*}/E^{*}$ and antisymmetric rotation tensor
$\bm{W}=\bm{W}^{*}/E^{*}$. We shall choose the characteristic shear rate
$E^{*}=(E_{ij}^{*}E_{ij}^{*})^{1/2}$ so that $E_{ij}E_{ij}=1$.

The spheroid is centred at the origin and its semiaxes $a_{i}=(a,c,c)$ are
spanned by an orthogonal set of unit vectors $\bm{p},\bm{q},\bm{r}$. Thus, the
coordinate $\bm{x}$ in the body frame maps to the laboratory frame as
$\bm{y}=\bm{Rx}$, where $\bm{R}=[\bm{p},\bm{q},\bm{r}]$. The unit vector
$\bm{p}$ points along the symmetry axis towards the “pole” of the spheroid,
whilst $\bm{q}$ and $\bm{r}$ point radially outward around the “equator”. The
body rotates with angular velocity $\bm{\Omega}=\bm{\Omega}(t)$, such that
$\dot{\bm{R}}=[\bm{\Omega}]_{\times}\bm{R}$, and the velocity field in the
body frame is

$\bm{u}(\bm{x},t)=\bm{R}^{\textrm{T}}\bm{v}(\bm{R}\bm{x},t)-\bm{\Omega}\times(\bm{Rx}).$
(36)

Thus, the background velocity gradient $\bm{u}=\bm{Ax}$ appears in the body
frame as

$\bm{A}=\bm{R}^{\textrm{T}}(\bm{G}-[\bm{\Omega}]_{\times})\bm{R}$ (37)

and, in general, varies in time as the particle rotates.

Jeffery (1922) showed that, when the couple upon a spheroid is zero, the solid
body rotation rate $\bm{\Omega}$ of the body is given by

$\bm{\Omega}=\frac{1}{2}\bm{\omega}+\frac{\lambda^{2}-1}{\lambda^{2}+1}\bm{p}\times\bm{E}\bm{p},$
(38)

where $\omega_{i}=-\epsilon_{ijk}W_{jk}$ is the background vorticity in the
laboratory frame. Thus, the orientation of the particle $\bm{p}$ then evolves
according to Jeffery’s equation

$\displaystyle\dot{\bm{p}}$
$\displaystyle=\bm{W}\bm{p}+\frac{\lambda^{2}-1}{\lambda^{2}+1}\left(\bm{Ep}-\bm{p}^{\textrm{T}}\bm{Ep}\right)$
(39)

where $W_{ij}=\frac{1}{2}\epsilon_{ijk}\omega_{k}$ is the rate of rotation
tensor.

When the background shear is constant in time, the solution to (39) can be
written in terms of a matrix exponential as

$\bm{p}(t)=\frac{\breve{\bm{p}}}{|\breve{\bm{p}}|}~{},~{}\breve{\bm{p}}(t)=\exp\left(\bm{K}t\right)\bm{p}_{0}~{},~{}\bm{K}=\bm{W}+\frac{\lambda^{2}-1}{\lambda^{2}+1}\bm{E}$
(40)

subject to the initial condition $\bm{p}=\bm{p}_{0}$ at $t=0$ (Szeri, 1993).

We shall now examine the fixed points and stable attractors of (39). This
analysis elaborates upon the previous work by Bretherton (1962). Decomposing
$\bm{K}$ in terms of its eigenvectors $\bm{e}_{i}$ and eigenvalues
$\lambda_{i}$, we can write (40) as

$\breve{\bm{p}}=\sum_{i}(\bm{p}_{0}\cdot\bm{e}_{i})\exp(\lambda_{i}t)\bm{e}_{i}.$
(41)

From (41) we see that the fixed points of (39) coincide with the eigenvectors
$\bm{e}_{i}$. The stability of the fixed points depends upon the eigenvalues
$\lambda_{i}$, and because $\sum_{i}\lambda_{i}=0$, there is always only one
(neutrally) stable fixed point or limit cycle (excepting $\lambda_{i}=0$).
Thus after a finite time, the motion of the particle will always approach a
stable attractor whose nature depends upon the largest eigenvalue(s) of
$\bm{K}$.

We can categorise the stable attractors into five different cases: 1a, 1b, 2a,
2b and 3. These five cases map to four different motions: resting, spinning,
2D (two-dimensional) tumbling and 3D (three-dimensional) tumbling. These are
summarised in table 1. We shall now examine each case.

Case | Eigenvalues | Condition | Motion | $\overline{\vphantom{\big{(}}\bm{A}}$
---|---|---|---|---
1a | $\lambda_{i}\in\mathbb{R},\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}$ | $\Omega_{3}\neq 0$ | Spinning | (44)
1b | ” | $\Omega_{3}=0$ | Resting | (46)
2a | $\lambda_{1}=\lambda_{2}^{\dagger}=\sigma+\mathrm{i}\omega,\lambda_{3}=-2\sigma$ | $\sigma>0$ | Spinning | (44)
2b | ” | $\sigma<0$ | 2D tumbling |
3 | ” | $\sigma=0$ | 3D tumbling |
Table 1: Classification of the motion and time average perceived velocity
gradient experienced by a spheroid.

In cases 1a and 1b, all the eigenvalues
$\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}$ of $\bm{K}$ are real and (40) has
three fixed points $\bm{p}=\bm{e}_{i}$, where $\bm{e}_{i}$ are the
eigenvectors of $\bm{K}$. From (41) we see that only the fixed point
corresponding to the largest eigenvalue $\lambda_{3}$ is stable, so for nearly
any initial condition, the particle will reorient itself to be parallel to
$\bm{e}_{3}$. In this configuration, the particle will in general rotate about
its axis with angular velocity $\bm{\Omega}=\Omega_{3}\bm{e}_{3}$, where
$\Omega_{3}=\bm{\omega}\cdot\bm{e}_{3}/2$ (case 1a). We call this motion
spinning. However, for particular configurations, $\Omega_{3}=0$ and the
particle does not rotate (case 1b). We call this motion resting.

In cases 2a, 2b and 3, $\bm{K}$ has a pair of complex eigenvalues and the
system has one fixed point and a periodic point. In these cases, we can write
the eigenvalues as
$\lambda_{1}=\lambda_{2}^{\dagger}=\sigma+\textrm{i}\omega,\lambda_{3}=-2\sigma$,
where $\dagger$ denotes the complex conjugate. When $\sigma<0$, the fixed
point is stable and the periodic point is unstable (case 2a, spinning). Thus
the particle aligns with $\bm{e}_{3}$ as in cases 1a and 1b, but unlike case
1b, the particle must rotate about this axis because the angular velocity
$\bm{\Omega}=\Omega_{3}\bm{e}_{3}$ is bounded
$\Omega_{3}^{2}\geq\omega^{2}>0$. When $\sigma>0$, the fixed point is unstable
and the periodic point is stable (case 2b). Thus the particle will evolve
towards a planar limit cycle oscillation described by

$\breve{\bm{p}}=\bm{a}\cos(\omega t)+\bm{b}\sin(\omega t)$ (42)

where $\bm{e}_{1}=\bm{e}_{2}^{\dagger}=\bm{a}+\textrm{i}\bm{b}$. We call this
motion 2D tumbling. In case 3, $\sigma=0$ and the motion follows three-
dimensional Jeffery orbits (Jeffery, 1922), where $\bm{p}$ precesses around
the axis $\bm{e}_{3}$ with period $2\pi/\omega$. We call this motion 3D
tumbling. 222Jeffery orbits are in general 3D, but for some initial
conditions, the motion can be reduced to spinning or 2D tumbling.

The solution of (26) depends upon the average flow field in the particle
frame. The instantaneous flow field around the particle consists of a
superposition of the background gradient $\bm{Ax}$ and a perturbation owing to
the presence of the particle. Since this perturbation is linear in $\bm{A}$,
the time average flow field in the particle frame is determined by the average
velocity gradient perceived by the particle. We shall now calculate the
average velocity gradient perceived by the particle in the general limiting
cases identified above.

(a) Spinning

(b) 2D tumbling

(c) 3D tumbling

Figure 2: Motion of a freely suspended spheroid: (2(a)) spinning, (2(b)) 2D
tumbling and (2(c)) 3D tumbling. For 2D tumbling and spinning, the motion is
periodic, whereas in for 3D tumbling only the motion of the symmetry axis is
necessarily periodic.

In case 1a and 2a (spinning), the particle rotates around a fixed axis
$\bm{p}=\bm{e}_{3}$, which is illustrated in figure 2(a). For this spinning
motion, the unit vectors of the body frame rotate as

$\bm{p}=\bm{e}_{3}~{},~{}\bm{q}=\cos(\Omega_{3}t)\bm{q}_{0}+\sin(\Omega_{3}t)\bm{r}_{0}~{},~{}\bm{r}=-\sin(\Omega_{3}t)\bm{q}_{0}+\cos(\Omega_{3}t)\bm{r}_{0}$
(43)

To obtain the time average velocity gradient, we can substitute (43) into (37)
and take the time average over one revolution, whose period is
$T=2\pi/\Omega_{3}$. After a little algebra, we obtain

$\overline{\vphantom{\big{(}}\bm{A}}=\frac{1}{T}\int_{0}^{T}\bm{R}^{\textrm{T}}(\bm{G}-[\bm{\Omega}_{\infty}]_{\times})\bm{R}\mathrm{d}t=E_{3}\begin{bmatrix}1&0&0\\\
0&-\frac{1}{2}&0\\\ 0&0&-\frac{1}{2}\end{bmatrix}$ (44)

where

$E_{3}=\bm{e}_{3}^{\textrm{T}}\bm{E}\bm{e}_{3}=\frac{\lambda^{2}+1}{\lambda^{2}-1}\sigma$
(45)

is the rate of strain perceived by the particle along its fixed axis of
rotation. Thus, the average velocity field experienced by a particle steadily
rotating about its axis is an axisymmetric strain, whose magnitude is given by
the component of strain along the rotation axis.

In case 1b (resting), the particle axis aligns with $\bm{e}_{3}$, but the
rotation rate about this axis vanishes. The average velocity field in the body
frame is then simply

$\overline{\vphantom{\big{(}}\bm{A}}=\bm{RGR}^{T}.$ (46)

The average velocity gradient $\overline{\vphantom{\big{(}}\bm{A}}$ can be
specified with three parameters, up to an arbitrary scaling and rotation about
the symmetry axis of the particle. To see this, we note that $\bm{G}$ has nine
components with eight degrees of freedom, since continuity requires
$\textrm{tr}\left(\bm{G}\right)=0$. The requirement that $\bm{\Omega}=0$
imposes the constraint

$\bm{\omega}=-2\frac{\lambda^{2}-1}{\lambda^{2}+1}\bm{p}\times\bm{E}\bm{p}$
(47)

reducing the number of unknowns by three. The choice of scale introduces
another redundant parameter; our non-dimensionalisation requires
$E_{ij}E_{ij}=1$. Finally, we note that rotations about the particle’s axis
are trivial, since the particle is axisymmetric.

In case 2b (2D tumbling), the motion of the particle is more complex, but
remains periodic, as illustrated in figure 2(b). From (42) it follows that the
symmetry axis $\bm{p}$ rotates in a plane spanned by vectors $\bm{a},\bm{b}$
with period $T=2\pi/\omega$. Thus the rotation
$\bm{Q}(t)=\bm{Q}_{2}\bm{Q}_{1}$ of the body frame
$\bm{R}(t)=\bm{Q}\bm{R}_{0}$ can be composed of as a rotation $\bm{Q}_{1}(t)$
about the axis $\bm{a}\times\bm{b}$ (mapping $\bm{p}_{0}$ onto $\bm{p}$),
followed by a rotation $\bm{Q}_{2}(t)$ about the axis
$\bm{p}=\bm{Q}_{1}\bm{p}_{0}$ by an angle

$\theta_{p}(t)=\int_{0}^{t}\frac{1}{2}\bm{\omega}\cdot\bm{p}(t^{\prime})\mathrm{d}t^{\prime}$
(48)

relative to the plane normal. By inspection of (48) and (42), we see that
$\theta_{p}(T)=0$, since $\bm{p}(t)=-\bm{p}(t+T/2)$. Thus, $\bm{R}(t)$ must be
periodic. The average perceived velocity gradient can then be evaluated
analogously to (44), but the resultant expression is much more complicated. It
suffices to remark that the average strain $\overline{\vphantom{\big{(}}E}$
and vorticity $\overline{\vphantom{\big{(}}\omega}$ components of the average
_perceived_ velocity gradient $\overline{\vphantom{\big{(}}\bm{A}}$ must also
satisfy (47), since in the body frame, the apparent rotation rate of the body
must also be zero. It follows that in case 2b, as in case 1b, the average
perceived velocity gradient may also be described by only three parameters, up
to a trivial rotation and choice of scale.

In case 3 (3D tumbling), the particle motion becomes fully three dimensional,
as illustrated in figure 2(c). Only the motion of its symmetry axis is
necessarily periodic; the equatorial axes $\bm{q}$ and $\bm{r}$ may precess
around and do not necessarily trace a path with the same period. Thus, we
cannot expect the flow in the particle frame to be time periodic, since
$\bm{q}$ and $\bm{r}$ point in a new direction at the start of each cycle, and
we cannot expect to apply (24) in this special case.

We recall that the derivation of (24) required pathlines of the flow to be
open. Far from the particle, pathlines follow $\dot{\bm{y}}=\bm{Gy}$, with
solution $\bm{y}=\exp(\bm{G}t)\bm{y}_{0}$. By a similar argument to that
above, we see that when $\bm{G}$ has eigenvalues $0,\pm\mathrm{i}\omega$, the
pathlines are closed. This provides a necessary condition to apply (24). The
3D tumbling orbits identified by Jeffery (1922) happen to correspond to this
case, which further precludes application of (24) to the case of Jeffery
orbits.

To summarise: by an analysis of the stable attractors of (39) we can construct
five cases of the particle motion categorised by the eigenvalues of $\bm{K}$
(40) as shown in table 1. In cases 1a, 1b and 2a, the particle orients itself
along a fixed axis corresponding to the eigenvector $\bm{K}$ with largest real
eigenvalue. In cases 2b and 3, the particle undergoes a limit cycle
oscillation. In cases 1a, 1b, 2a and 2b, the relative flow field is steady or
periodic. In cases 1a, 2a (spinning) the average perceived velocity gradient
over one period is an axisymmetric strain, because the particle rotates
steadily about its axis. In cases 1b and 2b (resting or 2D tumbling), the
average perceived velocity gradient satisfies (47), which can be specified in
terms of three parameters up to a trivial rotation and choice of scale.

###  3.2 The fluid motion near the particle

In the body frame, the surface at $\xi=0$ is stationary, so to leading order
in $\xi$ the velocity is

$\bm{u}=\xi\frac{\partial\bm{u}}{\partial\xi}\Big{|}_{\xi=0}+O(\xi^{2})$ (49)

After differentiation of Jeffery’s solution for the relative velocity field
for an arbitrary ellipsoid (Jeffery, 1922), rewritten in the more compact
matrix notation used by Kim & Karrila (1991), we find that the velocity
gradient normal to the surface is given by

$\frac{\partial\bm{u}}{\partial\xi}\Big{|}_{\xi=0}=\bm{\Phi}\bm{h}_{\xi}-\bm{h}_{\xi}\left(\bm{h}_{\xi}\cdot(\bm{\Phi}\bm{h}_{\xi}\right))=\bm{w}$
(50)

where $\bm{h}_{\xi}=\bm{Dx}/|\bm{Dx}|$ is the unit vector surface normal,
$\bm{D}$ is a diagonal matrix whose entries are $a_{i}^{-2}$ and

$\bm{\Phi}=\frac{3}{4\pi a_{1}a_{2}a_{3}}\bm{S}.$ (51)

The expression for the stresslet $\bm{S}$ is more involved. We remark here
that it is a symmetric matrix with $\textrm{tr}\left(\bm{S}\right)=0$ whose
elements are a linear combination of the components of the velocity gradient
$\bm{A}$ with coefficients determined by geometry-dependent “resistance
functions”. Complete expressions for $\bm{S}$ can be found in Kim & Karrila
(1991). Equation (50) is valid for the general case of torque-free, tri-axial
ellipsoids. The action of a net torque upon the body adds an additional skew-
symmetric component to (51) which is not treated here.

Equation (50) now defines the surface streamlines, which are everywhere
tangent to the surface velocity gradient $\bm{w}$. The surface streamlines and
their critical points are illustrated for an arbitrary linear shear in figure
1. From (50) we see that critical points occur where the surface normal
$\bm{h}_{\xi}$ coincides with an eigenvector $\pm\bm{q}_{i}$ of the surface
gradient tensor $\bm{\Phi}$. Since $\bm{\Phi}$ is a symmetric matrix, its
eigenvalues are all real and its eigenvectors orthogonal. In general,
$\bm{\Phi}$ has three distinct eigenvalues ordered
$\phi_{1}<\phi_{2}<\phi_{3}$ and there are six critical points. In special
cases, $\bm{\Phi}$ has a repeated eigenvalue and there is a locus of critical
points.

The precise nature of the critical points can be seen by introducing the
surface potential

$\phi(\bm{x})=\bm{h}_{\xi}\cdot(\bm{\Phi}\bm{h}_{\xi})=\frac{\bm{x}^{\textrm{T}}\bm{D\Phi
Dx}}{\bm{x}^{\textrm{T}}\bm{D}^{2}\bm{x}}$ (52)

which has the property that

$\nabla\phi=\frac{2\bm{Dw}}{\bm{x}^{\textrm{T}}\bm{D}^{2}\bm{x}}.$ (53)

The gradient of $\phi$ measured along surface streamlines is
$\bm{w}\cdot\nabla\phi>0$ everywhere on the surface except at the critical
points $\bm{w}=0$ where the surface streamlines terminate. In other words, the
potential $\phi$ always increases as one moves along a surface streamline. The
potential takes a minimum value of $\phi=\phi_{1}$ when
$\bm{h}_{\xi}=\pm\bm{q}_{1}$ (red marker in figure 1), i.e. surface
streamlines originate from a “source” $\phi=\phi_{1}$. It takes a maximum
value of $\phi=\phi_{3}$ when $\bm{h}_{\xi}=\pm\bm{q}_{3}$ (blue marker), i.e.
surface streamlines terminate at a “sink” $\phi=\phi_{3}$. Lastly, when the
eigenvalues are distinct, all streamlines must pass through the contour
$\phi=\phi_{2}$, except at points $\bm{h}_{\xi}=\pm\bm{q}_{2}$, which are
saddle points (green marker). When there is a repeated eigenvalue, either
$\phi_{2}=\phi_{1}$ or $\phi_{2}=\phi_{3}$, so there exists a locus of points
where streamlines originate (or terminate). It should be noted that the
definition of $\phi$ (52) and its properties extend also to torque-free, tri-
axial ellipsoids.

###  3.3 Surface flux under spinning motion

We now proceed to evaluate the surface streamlines about a spheroid aligned
with the axis of an axisymmetric straining flow. This configuration
corresponds to the mean flow field about a spinning spheroid. In the body
frame, the velocity gradient tensor is a diagonal matrix whose elements are
$A_{11}/2=-A_{22}=-A_{33}=E_{3}$. Likewise, the surface velocity gradient
tensor $\bm{\Phi}$ is also diagonal matrix with
$\Phi_{11}/2=-\Phi_{22}=-\Phi_{33}=\beta E_{3}/2$, where the geometry-
dependent prefactor $\beta$ is given by

$\frac{1}{\beta}=\frac{3}{4}\int_{0}^{\infty}\frac{ac^{2}t}{(a^{2}+t)^{3/2}(c^{2}+t)^{2}}\mathrm{d}t.$
(54)

For this configuration, the surface streamlines are axisymmetric, originating
at the “pole” and terminating at the “equator”. Therefore, we can use an
ellipsoidal polar coordinate system to describe the surface, identifying
$\eta$ with the polar angle between the symmetry axis and a point on the
surface and $\zeta$ with the azimuthal angle. We parametrise a point near the
surface as

$\bm{x}=\begin{bmatrix}a\cos\eta\\\ c\sin\eta\cos\zeta\\\
c\sin\eta\sin\zeta\end{bmatrix}+\xi\bm{h}_{\xi}.$ (55)

so that the covariant coordinate vectors at the surface are

$\bm{h}_{\xi}=\frac{1}{h_{\xi}}\begin{bmatrix}c\cos\eta\\\
a\sin\eta\cos\zeta\\\ a\sin\eta\sin\zeta\\\ \end{bmatrix}\\\
,~{}~{}\bm{h}_{\eta}=\begin{bmatrix}-a\sin\eta\\\ c\cos\eta\cos\zeta\\\
c\cos\eta\sin\zeta\\\ \end{bmatrix},~{}~{}\bm{h}_{\zeta}=\begin{bmatrix}0\\\
-c\sin\eta\sin\zeta\\\ +c\sin\eta\cos\zeta\\\ \end{bmatrix}$ (56)

with $h_{\xi}$ chosen to ensure that $\bm{h}_{\xi}$ is a unit vector. Then it
follows that

$\displaystyle
u_{\zeta}=0,~{}~{}u_{\eta}=\frac{3}{2}\xi|\Phi_{11}|\frac{ac\cos\eta\sin\eta}{(a^{2}\sin^{2}\eta+c^{2}\cos^{2}\eta)^{3/2}}=F\xi$
(57)

as required and

$\rho=c\sin\eta(a^{2}\sin^{2}\eta+c^{2}\cos^{2}\eta)^{1/2}.$ (58)

Substituting (57) and (58) into (24), we obtain

$\displaystyle\int_{\eta_{0}}^{\eta_{1}}{\rho}^{3/2}{F}^{1/2}\mathrm{d}\eta$
$\displaystyle=\int_{0}^{\pi/2}\left(\frac{3}{2}\beta|{E}_{3}|{a}{c}^{4}\cos\eta\sin^{4}\eta\right)^{1/2}\mathrm{d}\eta$
(59)
$\displaystyle=\left(\frac{\pi{a}{c}^{4}\beta|{E}_{3}|}{6}\right)^{1/2}\frac{\Gamma(\frac{7}{4})}{\Gamma(\frac{9}{4})}$

and thus

Sh $\displaystyle=0.566({a}{c}^{4}\beta)^{1/3}|E_{3}|^{1/3}\textrm{Pe}^{1/3}$
(60) $\displaystyle=\alpha_{||}\left(\frac{|E_{3}|^{*}r}{\kappa}\right)^{1/3}$

where $\alpha_{||}(\lambda)=0.566({a}{c}^{4}\beta)^{1/3}$ represents the
dependence of the surface flux upon geometry and $|E_{3}|^{*}r/\kappa$ is the
Péclet number based on the strain perceived along the axis of rotation.

The geometry dependent prefactor $\alpha_{||}$ has a relatively weak
dependence upon the aspect ratio of the spheroid, varying between $0.762$ to
$1.042$ over the interval $1/20\leq\lambda\leq 20$. In contrast, from the
definition of the characteristic shear rate $E^{*}$, $|E_{3}|$ may vary over
the interval $0<|E_{3}|\leq 2/\sqrt{6}$, depending on the particular
configuration of the background strain and the axis of the particle.
Therefore, the alignment of the particle with the background velocity gradient
plays a significant role in determining the surface flux for spinning
particles. This is the phenomenon of the partial suppression of convection by
rotation first identified by Batchelor (1979).

Some caution must be taken in utilising the result of (60). In our derivation
of (24), we have assumed that the boundary layer thickness remains small in
comparison to radius of curvature of the surface, so that cross-surface
diffusion and higher-order convection terms are negligible. Yet, when the body
is made infinitely slender, this assumption is violated. We expect that higher
order corrections to (24) are required for slender bodies at moderate Péclet
number.

###  3.4 Surface flux under tumbling and resting motion

Under tumbling and resting motion, convenient expressions for the surface
streamlines are no longer easy to find by hand. To proceed in this general
case, we adopt a numerical approach to parametrise the surface in terms of
coordinates $\bm{x}=\bm{x}(0,\eta,\zeta)$. Essentially, the task is to draw a
set of $j=1\ldots n_{\zeta}$ surface streamlines covering the body, label each
streamline with a unique $\zeta=\zeta_{j}$ and evaluate the position at
$i=1\ldots n_{\eta}$ different locations $\eta=\eta_{i}$ along the streamline.
Then we can create a $n_{\eta}\times n_{\zeta}$ mesh of points
$\bm{x}_{i,j}=\bm{x}(0,\eta_{i},\zeta_{j})$ upon the surface to numerically
approximate the metric $\rho$, which can be used to numerically integrate
(24).

We shall now construct a suitable definition of the streamwise coordinate
$\eta$. We require the surface streamlines be tangent to the surface velocity
gradient $\bm{w}$ (50), so $\bm{h}_{\eta}$ must be of the form

$\bm{h}_{\eta}=\frac{\mathrm{d}\bm{x}}{\mathrm{d}\eta}=\frac{\bm{w}}{\bm{w}\cdot\nabla\eta}$
(61)

which has $\nabla\eta\cdot\bm{h}_{\eta}=1$ and
$\bm{h}_{\eta}\cdot\bm{h}_{\xi}=0$, as required. Furthermore, we require that
$\eta$ should increase monotonically along surface streamlines. We have seen
in §3.2 that the surface potential $\phi$ has this property. Therefore, we
identify $\phi$ (52) with the streamwise coordinate $\eta$.

We require a definition of the $\zeta$ coordinate. Since $\zeta$ is constant
along surface streamlines and every surface streamline passes through
$\eta=\phi_{2}$, $\zeta$ can be thought of as a coordinate along the curve
$\bm{x}_{0}(\zeta)=\bm{x}(0,\phi_{2},\zeta)$. Along $\eta=\phi_{2}$, from (52)
we derive that the unit surface normal satisfies

$\bm{h}_{\xi}\cdot\bm{q}_{1}\pm\sqrt{\frac{\phi_{2}-\phi_{3}}{\phi_{2}-\phi_{1}}}(\bm{h}_{\xi}\cdot\bm{q}_{3})=0$
(62)

so the constraint

$\zeta=\bm{h}_{\xi}\cdot\bm{q}_{2}$ (63)

describes a point along $\bm{x}_{0}(\zeta)$. The surface can split into four
quadrants, depending upon the basin of attraction of the surface streamlines.
To wit, the behaviour in the limits $\bm{h}_{\xi}\rightarrow\pm\bm{q}_{1}$ as
$\eta\rightarrow\phi_{1}$ and $\bm{h}_{\xi}\rightarrow\pm\bm{q}_{3}$ as
$\eta\rightarrow\phi_{3}$ forms our classification. Therefore, in each
quadrant, the coordinates $(\eta,\zeta)\in[-\phi_{1},\phi_{3}]\times[-1,1]$
uniquely define a point on the surface. This definition is general for torque-
free, tri-axial ellipsoids.

We now outline the procedure to evaluate (24) numerically. For each surface
quadrant, we choose a set of points
$\bm{x}_{0,j}=\bm{x}(0,\phi_{2},\zeta_{j})$ which lie on the ellipsoid surface
$\xi=0$ and satisfy $\eta(\bm{x})=\phi_{2}$ (52). We numerically integrate
(61) from the initial condition $\bm{x}=\bm{x}_{0,i}$ at $\eta=\phi_{2}$ to
$\eta=\eta_{j}$, which yields a mesh of points
$\bm{x}_{i,j}=\bm{x}(0,\eta_{i},\zeta_{j})$ upon the particle surface.
Measuring the surface area $\delta S_{ij}$ of the region enclosed over
$[\eta_{i},\eta_{i+1}]\times[\zeta_{i},\zeta_{i+1}]$, we approximate the
surface area density for this surface element as

$\rho_{ij}\approx\delta S_{ij}(\eta_{i+1}-\eta_{i})(\zeta_{i+1}-\zeta_{i})$
(64)

which is, roughly speaking, a “cell average” of $\rho$. The velocity component
$F$ (8) is evaluated at $\bm{x}_{i,j}$ from (52) and (50). Then, the integral
(24) is evaluated as a Riemann sum over the surface elements.

A technical point remains that $\eta_{i}$ and $\zeta_{j}$ should be chosen so
that the surface is densely covered in mesh points $\bm{x}_{i,j}$. To achieve
this, we can generate a uniform sampling of seed points on the surface. These
seed points can be integrated (61) numerically towards the curve $\bm{x}_{0}$
to construct $\zeta_{j}$. This guarantees a satisfactory coverage of the
surface by the streamlines. The $\eta_{i}$ can be chosen as

$\eta_{i}=\begin{cases}\phi_{1}\cos^{2}\nu_{i}+\phi_{2}\sin^{2}\nu_{i}&\text{
if }\nu_{i}<0\\\ \phi_{2}\cos^{2}\nu_{i}+\phi_{3}\sin^{2}\nu_{i}&\text{ if
}\nu_{i}>0\end{cases}$ (65)

for $\nu_{i}$ evenly distributed on the interval $[-\pi/2,\pi/2]$.

###  3.5 Surface flux in rotation dominated flow

It is useful to consider the special case $|\bm{\omega}|\rightarrow\infty$
where the vorticity is large and vortex stretching is non-zero
($\bm{\omega}^{\textrm{T}}\bm{E\omega}\neq 0$). In this case, it can be shown
that the eigenvalues of $\bm{K}$ are complex, its eigenvectors lie parallel
and perpendicular to the direction of the vorticity
$\bm{\hat{\omega}}=\bm{\omega}/|\bm{\omega}|$ and the particle rotates with
constant angular velocity $\bm{\Omega}\rightarrow\bm{\omega}/2$. Therefore,
the motion either corresponds to case 2a (spinning) or case 2b (tumbling). In
the spinning case, the particle rotates with its symmetry axis parallel to the
vorticity vector $\bm{p}=\bm{\hat{\omega}}$. In the tumbling case, the
particle rotates with its symmetry axis always orthogonal to the vorticity
vector, e.g. $\bm{r}=\bm{\hat{\omega}}$. The motion of the body frame is
therefore analogous to (43).

The configuration can be inferred from the exact eigenvalue relationship

$\lambda_{1}\lambda_{2}\lambda_{3}=-2\sigma(\sigma^{2}+\omega^{2})=\frac{\gamma^{3}}{3}\textrm{tr}\left(\bm{E}^{3}\right)+\frac{\gamma}{4}\bm{\omega}^{\textrm{T}}\bm{E\omega}$
(66)

where $\gamma=(\lambda^{2}-1)/(\lambda^{2}+1)$ is the shape co-factor. In the
limit $|\bm{\omega}|\rightarrow\infty$, $\sigma\rightarrow-\gamma
E_{\omega}/2$, where
$E_{\omega}=\bm{\hat{\omega}}^{\textrm{T}}\bm{E}\bm{\hat{\omega}}$ is the
strain rate perceived in the direction of vorticity. Therefore, following
§3.1, we identify that the particle is aligned in the parallel configuration
when $\gamma E_{\omega}>0$ and the orthogonal configuration when $\gamma
E_{\omega}<0$. It follows that, as in (44), the average perceived velocity
gradient is

$\overline{\bm{A}}=E_{\omega}\begin{bmatrix}1&0&0\\\ 0&-\frac{1}{2}&0\\\
0&0&-\frac{1}{2}\end{bmatrix}~{}\textrm{when}~{}~{}\gamma
E_{\omega}>0~{}~{}\textrm{or}~{}~{}E_{\omega}\begin{bmatrix}-\frac{1}{2}&0&0\\\
0&-\frac{1}{2}&0\\\ 0&0&1\end{bmatrix}~{}\textrm{when}~{}~{}\gamma
E_{\omega}<0.$ (67)

In both cases, the average perceived velocity gradient is an axisymmetric
strain. However, the alignment of the symmetry axis of the spheroid with this
strain depends upon the sign of the shape co-factor and vortex stretching. As
a result, we evaluate the surface flux as

$\textrm{Sh}=\left\\{\begin{matrix}\alpha_{||}\textrm{Pe}_{\omega}^{1/3}&\gamma
E_{\omega}>0\\\ \alpha_{\perp}\textrm{Pe}_{\omega}^{1/3}&\gamma
E_{\omega}<0\\\ \end{matrix}\right.$ (68)

where $\textrm{Pe}_{\omega}=E_{\omega}^{*}r/\kappa$ is the Péclet number based
on the vortex stretching, $\alpha_{||}(\lambda)$ is given by (60) and
$\alpha_{\perp}(\lambda)$ is obtained using the numerical procedure outlined
in §3.4.

Figure 3: Variation in the mass transfer coefficient of a spheroid of varying
aspect ratio but constant surface area in axisymmetric straining flow (black
lines) and uniform flow (blue line). The black circle shows Batchelor’s result
for a sphere in axisymmetric strain (Batchelor, 1979), whilst the black
diamond shows the equivalent result for a sphere in uniform flow (Acrivos &
Taylor, 1962).

We have plotted the geometry dependent prefactors $\alpha_{||}$ and
$\alpha_{\perp}$ in figure 3. The marker shows Batchelor’s result
$\textrm{Sh}=0.968\textrm{Pe}_{\omega}^{1/3}$ for the case of a sphere in
rotation dominated flow (Batchelor, 1979). We observe that, for the parallel-
aligned configuration $\alpha_{||}$, the variation in the prefactor with
$\lambda$ is relatively modest, corresponding to a variation in the surface
flux of between $-21.4\%$ and $+7.7\%$ relative to that of a sphere with
equivalent surface area. In contrast, there is a stronger variation observed
for the orthogonal configuration $\alpha_{\perp}$, which exhibits an
equivalent variation of $-10.5\%$ to $+53\%$ over the range shown. From (68),
we see that the surface flux is proportional to the lower branches of the two
black curves in vortex stretching $E_{\omega}>0$ and the upper branches in
vortex compression $E_{\omega}<0$. It follows that the _sign_ of the vortex
stretching can influence the surface flux. In particular, the surface flux to
prolate spheroids is significantly increased under vortex compression.

The shape dependence for spheroids in axisymmetric strain may be contrasted to
the shape dependence observed for axisymmetric, uniform flow around a fixed
spheroid shown in figure 3 (Acrivos & Taylor, 1962; Sehlin, 1969; Dehdashti &
Masoud, 2020). Here, the surface flux exhibits the same scaling
$\textrm{Sh}=\alpha_{U}\textrm{Pe}_{U}^{1/3}+O(1)$, where
$\textrm{Pe}_{U}=U_{\infty}^{*}r/\kappa$ and $U_{\infty}^{*}$ is the magnitude
of the relative velocity (slip velocity) between the free-stream and particle.
Whilst it is of limited value to compare the absolute values of $\alpha$, some
insight can be obtained by comparing the relative variation of each function.
We observe that for a spherical particle with fixed translation velocity and
surface area, flattening or elongating the particle results in a reduction in
the surface flux. Likewise, for a spherical particle in rotation-dominated
vortex stretching, flattening or elongating the particle also reduces the
surface flux. However, in rotation-dominated vortex compression, flattening or
elongating increases the surface flux. This observation is relevant to chain
formation by marine diatoms, where each cell in a chain must experience an
increased nutrient flux per unit area to benefit despite increased competition
for nutrients by its neighbours (Pahlow et al., 1997).

One may also observe that rotation dominated straining flow becomes
considerably more effective than uniform flow at transferring solute from the
particle surface at large aspect ratios. Of course, this result must depend
upon the relative orientation of the free-stream and particle axis, which is
not varied here. Nonetheless, similar comparisons have found utility in
determining when sinking or swimming may become a useful strategy for
phytoplankton to increase their nutrient flux in turbulent water (Karp-Boss et
al., 1996).

###  3.6 Further extensions

Two further generalisations of the surface coordinate system are to cases with
non-zero body forces and torques. This would be necessary, for instance, to
capture gyrotactic effects in the nutrient uptake of phytoplankton (Guasto et
al., 2012) or the behaviour of inertial particles. Such an extension can be
accommodated by a suitable modification of the surface coordinate system, e.g.
the inclusion of a body torque adds a skew-symmetric component to $\bm{\Phi}$
in (50). However, in this case, one would also expect the orientation dynamics
(and therefore the mean flow around the particle) to change, which would also
affect the mean surface flux via the convection-suppression effect.

##  4 Comparisons to numerical simulation

As a test of the results in in §2 and §3, we conducted numerical simulations
of the unsteady scalar transport around spheroids freely suspended in a linear
shear. We describe our numerical methods in §4.1 then examine two classes of
arbitrary strain and rotational flows in §4.2.

###  4.1 Methods

The unsteady convection diffusion equation (1) is solved using a second-order
finite volume method, subject to the Dirichlet boundary condition $\theta=1$
on $S_{p}$ and the von-Neumann boundary condition on the outer boundary of the
simulated domain. Equation (1) is discretised on a structured grid in prolate
(69) or oblate (70) spheroidal coordinates $(\mu,\theta,\phi)$, depending upon
the particle aspect ratio. The inner boundary $(\mu_{0},\theta,\phi)$
corresponds to the surface of the spheroid oriented in the Cartesian $x$
direction. The relative velocity field was evaluated in the body frame based
on the coordinate of each cell centre using the expressions given by Kim &
Karrila (1991). The relative velocity field satisfies the impermeable, no slip
boundary condition $\bm{u}=0$ on $S_{p}$.

$\displaystyle\bm{x}_{ijk}$
$\displaystyle=\begin{bmatrix}f\cosh\mu_{i}\cos\theta_{j}\\\
f\sinh\mu_{i}\sin\theta_{j}\cos\phi_{k}\\\
f\sinh\mu_{i}\sin\theta_{j}\sin\phi_{k}\end{bmatrix}~{},~{}\begin{aligned}
\mu_{i}&\in[\mu_{0},\mu_{\infty}]\\\ \theta_{j}&\in[0,\pi]\\\
\phi_{k}&\in[0,2\pi)\\\ \end{aligned}~{},~{}$ $\displaystyle
f=\sqrt{a^{2}-c^{2}}$ (69) $\displaystyle\bm{x}_{ijk}$
$\displaystyle=\begin{bmatrix}f\sinh\mu_{i}\sin\theta_{j}\\\
f\cosh\mu_{i}\cos\theta_{j}\cos\phi_{k}\\\
f\cosh\mu_{i}\cos\theta_{j}\sin\phi_{k}\end{bmatrix}~{},~{}\begin{aligned}
\mu_{i}&\in[\mu_{0},\mu_{\infty}]\\\ \theta_{j}&\in[-\pi/2,\pi/2]\\\
\phi_{k}&\in[0,2\pi)\\\ \end{aligned}~{},~{}$ $\displaystyle
f=\sqrt{c^{2}-a^{2}}$ (70)

The mesh resolution and spacing was chosen following the studies of Pahlow et
al. (1997) and Karp-Boss et al. (1996). The mesh is discretised into
$150\times 64\times 64$ cells in the $(\mu,\theta,\phi)$ directions
respectively, with uniform spacing in the $\theta$ and $\phi$ directions. Due
to the nature of the spheroidal coordinate system chosen, the resolution
varies across the surface of the particle and the outer boundary is very
slightly oblate or prolate. The dimensions of the spheroid $a_{i}=(a,c,c)$ are
chosen such that the surface area is $4\pi$, equivalent in surface area to a
sphere with unit radius. The outer boundary is chosen such that its largest
dimension is $100$ and is very nearly spherical, having an aspect ratio
between $0.999$ and $1.001$. To adequately resolve the thin concentration
boundary layer, which is of thickness $\delta=\textrm{Pe}^{-1/3}$, we employ a
mesh refinement in the $\mu$ direction such that the grid spacing is
$\Delta\mu_{i+1}=R\Delta\mu_{i}$, where $\Delta\mu_{i+1}=\mu_{i+1}-\mu_{i}$ is
the spacing between adjacent cells in the $\mu$ direction. The initial spacing
$\Delta\mu_{1}$ is chosen such that the thickness of the largest cell near the
particle surface is at most $2\times 10^{-4}\max(a,c)$ and the mesh refinement
factor $R$ is chosen accordingly. For the most extreme aspect ratio
$\lambda=16$ ($\lambda=1/16$) at the largest Péclet number tested, this yields
between $27$ (43) and $70$ (87) cells within a distance $\delta$ from the
surface.

The solver is based on the scalarTransportFoam solver of OpenFOAM, which was
modified to solve (1) for a time varying flow field. The convective term in
(1) is discretised using a standard linear upwind Gaussian integration, and
the diffusive term is discretised using a similar linear Gaussian scheme with
an explicit non-orthogonal correction to maintain second order accuracy. Time
stepping is performed using an implicit Euler scheme and a time step of
$\Delta t=0.02$. Simulations were allowed sufficient time for the surface flux
to reach a steady (or periodic) state. Where the particle motion is unsteady,
the cycle average mass flux was evaluated over the interval $9T\leq t<10T$.
Where the particle motion is steady, the steady state mass flux was evaluated
at $t=100$.

###  4.2 Results and discussion

In this section, we shall compare the results of our numerical simulations
against the asymptotic results derived in §2.

####  4.2.1 Pure strain

As our first test, we consider an arbitrary irrotational background velocity
gradient $\bm{G}=\bm{E}$. The particle motion in this case corresponds to case
1b of §3.1 and the particle aligns itself with the most extensional (or
compressive) direction of strain. Thus, the relative velocity gradient field
is of the form

$\bm{A}=\begin{bmatrix}\sigma_{1}&0&0\\\ 0&\sigma_{2}&0\\\
0&0&\sigma_{3}\end{bmatrix}$ (71)

where $\sigma_{i}$ are the eigenvalues of $\bm{E}$. The topology of the
relative flow field can therefore described by a single parameter $-1\leq
s\leq 1$ (Lund & Rogers, 1994)

$s=\frac{-3\sqrt{6}\sigma_{1}\sigma_{2}\sigma_{3}}{(\sigma_{1}^{2}+\sigma_{2}^{2}+\sigma_{3}^{2})^{3/2}}.$
(72)

When $s=1$, the background flow is an axisymmetric extensional flow; $s=0$
corresponds to a two-dimensional strain and $s=-1$ corresponds to an
axisymmetric compression.

(a) Prolate spheroid, $\lambda=4$

(b) Oblate spheroid, $\lambda=1/4$

Figure 4: Mass transfer coefficient for (4(a)) prolate and (4(b)) oblate
spheroids in a pure straining flow, as a function of the strain topology
parameter $s$. Coloured lines with markers show numerical simulations over
$\textrm{Pe}=10^{1}-10^{4}$. Black lines show the expected scaling coefficient
$\alpha(s,\lambda)$ (24). Black circular markers show the axisymmetric result
(60).

We computed the surface flux from prolate and oblate spheroids in an arbitrary
straining flow over a range of Pe. The result is shown in figure 4 for cases
$\lambda=4$ and $1/4$. For both prolate and oblate spheroids, there is a
modest variation ($15-30\%$) in the mass transfer coefficient as the strain
topology is varied, which becomes more pronounced as the Péclet number is
increased. In contrast, the mass transfer coefficient for a spherical body in
a pure straining flow varies by less than $1\%$ over the same range of $s$
(Batchelor, 1979). When the Péclet number is large, the mass transfer rate
approaches the limiting scaling $\textrm{Sh}=\alpha\textrm{Pe}^{1/3}$. The
numerical coefficient $\alpha$ is well predicted by (24) and shows the correct
qualitative dependence upon $s$. At $\textrm{Pe}=10^{4}$, the discrepancy in
the predicted mass transfer rate between theory and numerics is within $2.5\%$
for the prolate case and within $3.1\%$ for the oblate case. This is partly
due to the mesh refinement; additional tests at a higher resolution of
$300\times 128\times 128$ suggest the numerics slightly under-resolve the
surface flux by around 2% for the oblate case, which would improve the
agreement seen here. Nonetheless, the leading order correction term to (24) is
$O(1)$, which corresponds to a discrepancy in $\textrm{Sh}/\textrm{Pe}^{1/3}$
of $\sim 4.6\%$ at $\textrm{Pe}=10^{4}$ and is comparable to the observed
discrepancy.

To examine the role of particle shape, we plot the scaling coefficient
$\alpha(s,\lambda)$ over a range of $s$ and $\lambda$ in figure 5. Markers
show the value of $\textrm{Sh}/\textrm{Pe}^{1/3}$ from our numerical
simulations at $\textrm{Pe}=10^{4}$ whilst the ruled surface shows the result
of (24). We observe that prolate spheroids tend to experience a larger surface
flux than oblate spheroids of equivalent surface area in the same flow.
However, the trend not always clear-cut and is reversed for strain topologies
near $s=1$, where spherical particles experience a larger surface flux.

Figure 5: Mass transfer coefficient for spheroids in a pure straining flow, as
a function of the strain topology parameter $s$ and aspect ratio $\lambda$.
The ruled surface shows the expected scaling coefficient $\alpha(s,\lambda)$
(24). Black markers show the result of numerical simulations conducted at
$\textrm{Pe}=10^{4}$.

In figure 4 and figure 5, we observe that the surface flux is always larger in
axisymmetric compression ($s=-1,\sigma_{1}>0,\sigma_{2}=\sigma_{3}<0$) than
axisymmetric extension ($s=+1,\sigma_{1}<0,\sigma_{2}=\sigma_{3}>0$). This is
remarkable, since both flows correspond to an axisymmetric strain; only the
direction of the flow is reversed. At first glance, this appears to violate
Brenner’s flow reversal theorem (Brenner, 1967; Masoud & Stone, 2019), which
states that for an isothermal body in steady flow, the surface flux is
preserved under flow reversal $\bm{u}\rightarrow-\bm{u}$. However, we recall
that the stable orientation of a free spheroid shifts under flow reversal. In
this example, prolate spheroids align parallel to the Cartesian $1-$direction
in axisymmetric compression ($s=-1$) and orthogonal to it in axisymmetric
extension ($s=+1$). In fact, the flow configuration is identical to the
parallel and orthogonal flow configurations discussed in §3.5 and shown in
figure 3. Thus, the difference in the average transfer rate between $s=\pm 1$
is due to a change in the stable orientation of the spheroid.

####  4.2.2 Spinning and tumbling

As a numerical test of our result for spinning and tumbling spheroids, we
consider a spheroid in a background velocity gradient $\bm{G}=\bm{E}+\bm{W}$
of the form

$\bm{E}=\frac{1}{\sqrt{6}}\begin{bmatrix}2&0&0\\\ 0&-1&0\\\
0&0&-1\end{bmatrix}~{},~{}\bm{W}=\begin{bmatrix}0&-\sin\theta_{\omega}&0\\\
\sin\theta_{\omega}&0&-\cos\theta_{\omega}\\\
0&\cos\theta_{\omega}&0\end{bmatrix}$ (73)

This corresponds to an axisymmetric straining flow with $|\bm{E}|=1$ and
vorticity magnitude $|\omega|=2$. By varying the angle $\theta_{\omega}$
between the background vorticity $\bm{\omega}$ and the $x_{1}$ axis, we survey
different limiting behaviours of particle motion identified in §3.1 and the
relative flow field experienced by the particle depends upon is geometry. For
example, a prolate particle will spin when $\theta_{\omega}=0$ (case 2a)
whereas it will tumble for $\theta_{\omega}=\pi/2$ (case 2b). The situation is
reversed for an oblate particle.

(a) Prolate spheroid, $\lambda=4$

(b) Oblate spheroid, $\lambda=1/4$

Figure 6: Average mass transfer coefficient for (6(a)) prolate and (6(b))
oblate spheroids in an arbitrary shear (73), as an implicit function of the
mean axial strain rate $E_{3}$ (45). Coloured lines with markers show
numerical simulations over $\textrm{Pe}=10^{1}-10^{4}$. Black solid lines show
the expected scaling coefficient for spinning motion (60). Black dashed lines
show the expected scaling coefficient for tumbling motion (60).

Figure 6 shows how the average mass transfer rate varies for prolate
($\lambda=4$) and oblate $(\lambda=1/4)$ spheroids as the parameter
$\theta_{\omega}$ is varied. The dependence is plotted as an implicit function
of the axial strain rate $E_{3}$ (45) for $\textrm{Pe}=10^{1}-10^{4}$. We
first remark that as the Péclet number becomes large, the mass transfer rate
approaches the scaling predicted by (60) and (24). As $\theta_{\omega}$ is
varied, the particle motion switches from spinning to tumbling at $E_{3}=0$
and a pronounced change can be observed in the limiting mass transfer rate.
There is a marked suppression of mass transfer near $E_{3}=0$ as the particle
approaches the transition from spinning to tumbling. This is a demonstration
of the suppression of convection by rotation first identified by Batchelor
(1979). We note, however, that the presence of vorticity does not always
suppress the convective transport. For instance, for the prolate spheroid in
figure 6(a), tumbling induced by the vorticity component can enhance the
convective transfer relative to the equivalent pure strain case (figure 4(a)).

Some remarks are in order. Firstly, the convection suppression/augmentation
effect is essentially a hydrodynamic effect: it occurs as the result of the
motion of a free spheroid and its alignment with the strain field. Secondly,
this effect can be seen even at moderate Péclet number and becomes more
pronounced as Pe increases. This underlines the observation that particle
shape influences two factors in determining the mass transfer rate: it
determines both the boundary condition for the scalar field and the behaviour
of relative flow field.

##  5 Conclusion

In this paper, we have presented a general method to evaluate the average
surface flux of solute from a rigid particle of arbitrary shape immersed in an
arbitrary, open pathline flow. The main restrictions upon the shape of the
particle are that it contains no sharp edges or regions of extreme curvature,
where the thin boundary layer assumption breaks down. The flow may be steady
or time periodic. When the flow is periodic, the average surface flux is
equivalent to that of the same particle embedded in the mean flow field,
provided the dimensionless period of the motion is $T\ll\textrm{Pe}^{1/3}$.
The Sherwood number scales as $\textrm{Pe}^{1/3}$, with a prefactor $\alpha$
which can be readily obtained through numerical integration once the particle
geometry and surface flow field are specified.

We apply this result to compute the surface flux from a small, freely
suspended, spheroid in a steady linear shear. To do so, we compute the
relative flow field experienced by the particle, which may be unsteady due to
the particle motion. Through an analysis of Jeffery’s equation, we identify
four categories of motion: resting, spinning and 2D or 3D tumbling. The
relative flow field is time periodic in the first three cases. In the spinning
case, the average perceived flow field is an axisymmetric strain. In the 2D
tumbling case, the average flow field always corresponds to an equivalent
spheroid in steady flow (resting). We provide a closed form expression (60)
for the surface flux in the spinning case. We outline the numerical procedure
to obtain the surface flux in the resting or 2D tumbling cases. We also
describe a simplification for the case of rotation dominated flow
$|\bm{\omega}|\rightarrow\infty$. In this limit, there is a larger surface
flux under vortex compression when compared to vortex stretching and this
increase becomes significant for prolate bodies. These procedures may serve as
the basis for analysing other, more complex geometries, or more complex rigid
body dynamics including inertia and gravity.

As a test of these analytical results, we have presented numerical simulations
of the scalar transport and surface flux around spheroids in pure straining
and a simplified shear flow. In all cases, we observe good agreement with the
expected scaling law and mass transfer coefficient, up to the accuracy of the
asymptotic approximation. In pure straining flows, the surface flux is steady
and is prescribed by only three parameters: the Péclet number, the particle
aspect ratio and a parameter describing the strain topology. In surveying this
parameter space, we observe that prolate spheroids tend to experience a
greater surface flux than oblate spheroids of equivalent surface area. When
vorticity is present, we observe that the spinning or tumbling of the particle
may suppress or augment the convective transfer, due to a realignment of the
particle with respect to the average strain field. Two additional parameters
are necessary to characterise the surface flux in all rotational flows and a
complete survey of the parameter space is not practical. However, the space is
sufficiently small that it may be readily tabulated for use as a model in a
numerical simulation.

We anticipate that our results may find application in modelling the mass
transfer from spheroidal or arbitrary particles in numerical simulations of
particle laden flows, where models for the interphase mass transfer rate are
required to take particle shape into account. Although the results presented
here pertain to steady and time periodic flows, the results may serve in the
same spirit in which steady flow mass transfer coefficients are employed to
model the transfer rate in unsteady flows (Crowe et al., 2012). Of particular
interest is the mass transfer in turbulent environments, such as turbulent
ocean waters home to planktonic osmotrophs (Karp-Boss et al., 1996) and
microplastics (Law, 2017). Here, the adaptation of shape to maximise surface
flux may help explain the great diversity in the morphology of osmotrophs
(Guasto et al., 2012), or help identify which shapes of microplastics do the
most harm. Another extension would be to include non-zero body torques or
forces. This would allow, for instance, the consideration of gyrotactic
effects in the motion of phytoplankton (Guasto et al., 2012), or inertial
effects in the rigid body dynamics of suspended particles.

##  6 Acknowledgements

This project has received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska-Curie grant
agreement No 846648. The authors acknowledge the use of the IRIDIS High
Performance Computing Facility, and associated support services at the
University of Southampton, in the completion of this work. The authors report
no conflict of interest.

## References

  * Acrivos (1960) Acrivos, Andreas 1960 Solution of the laminar boundary layer energy equation at high Prandtl numbers. Physics of Fluids 3 (4), 657–658.
  * Acrivos (1980) Acrivos, Andreas 1980 A note on the rate of heat or mass transfer from a small particle freely suspended in a linear shear field. Journal of Fluid Mechanics 98 (2), 299–304.
  * Acrivos & Taylor (1962) Acrivos, Andreas & Taylor, Thomas D. 1962 Heat and mass transfer from single spheres in stokes flow. Physics of Fluids 5 (4), 387–394.
  * Batchelor (1979) Batchelor, G.K. 1979 Mass transfer from a particle suspended in fluid with a steady linear ambient velocity distribution. Journal of Fluid Mechanics 95, 369–400.
  * Batchelor (1980) Batchelor, G. K. 1980 Mass transfer from small particles suspended in turbulent fluid. Journal of Fluid Mechanics 98 (3), 609–623.
  * Brenner (1967) Brenner, H. 1967 On the invariance of the heat-transfer coefficient to flow reversal in Stokes and potential streaming flows past particles of arbitrary shape. J. Math. Phys. Sci. 1, 173–179.
  * Bretherton (1962) Bretherton, F. P. 1962 The motion of rigid particles in a shear flow at low Reynolds number. Journal of Fluid Mechanics 14 (2), 284–304.
  * Clift et al. (1978) Clift, M., Grace, J. R. & Weber, M.E. 1978 Bubbles, Drops and Particles. New York: Academic Press, Inc.
  * Crowe et al. (2012) Crowe, C. T., Schwarzkopf, J. D., Sommerfeld, M. & Tsuji, Y. 2012 Multiphase flows with Droplets and Particles. CRC Press.
  * Dehdashti & Masoud (2020) Dehdashti, Esmaeil & Masoud, Hassan 2020 Forced Convection Heat Transfer from a Particle at Small and Large Peclet Numbers. Journal of Heat Transfer 142 (6), 1–9.
  * Feng & Michaelides (1997) Feng, Zhi Gang & Michaelides, Efstathios E. 1997 Unsteady Heat and Mass Transfer from a Spheroid. AIChE Journal 43 (3), 609–614.
  * Frankel & Acrivos (1968) Frankel, Neil A. & Acrivos, Andreas 1968 Heat and mass transfer from small spheres and cylinders freely suspended in shear flow. Physics of Fluids 11 (9), 1913–1918.
  * Goddard & Acrivos (1965) Goddard, J.D. & Acrivos, A. 1965 Asymptotic expansions for laminar forced-convection heat and mass transfer. J. Fluid Mech. 23 (1965), 273–291.
  * Grinfeld (2013) Grinfeld, Pavel 2013 Introduction to tensor analysis and the calculus of moving surfaces. New York: Springer.
  * Guasto et al. (2012) Guasto, Jeffrey S., Rusconi, Roberto & Stocker, Roman 2012 Fluid Mechanics of Planktonic Microorganisms. Annual Review of Fluid Mechanics 44 (1), 373–400.
  * Gupalo et al. (1976) Gupalo, Iu. P., Polianin, A. D. & Riazantsev, Iu. S. 1976 Diffusion to a particle at large Peclet numbers in the case of arbitrary axisymmetric flow over a viscous fluid. Journal of Applied Mathematics and Mechanics 40 (2), 893–898.
  * Gupalo & Riazantsev (1972) Gupalo, Iu. P. & Riazantsev, Iu. S. 1972 Diffusion on a particle in the shear flow of a viscous fluid. Approximation of the diffusion boundary layer. Journal of Applied Mathematics and Mechanics 36 (3), 447–451.
  * Jeffery (1922) Jeffery, G.B. 1922 The Motion o f Ellipsoidal Particles Immersed. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 102 (715), 161 – 179.
  * Karp-Boss et al. (1996) Karp-Boss, L., Boss, E. & Jumars, P. A. 1996 Nutrient fluxes to planktonic osmotrophs in the presence of fluid motion. Oceanography and Marine Biology: an annual review 34, 71–107.
  * Ke et al. (2018) Ke, Chunhai, Shu, Shi, Zhang, Hao, Yuan, Haizhuan & Yang, Dongmin 2018 On the drag coefficient and averaged Nusselt number of an ellipsoidal particle in a fluid. Powder Technology 325, 134–144.
  * Kim & Karrila (1991) Kim, Sangtae & Karrila, Seppo J. 1991 Microhydrodynamics. Elsevier.
  * Kishore & Gu (2011) Kishore, Nanda & Gu, Sai 2011 Momentum and heat transfer phenomena of spheroid particles at moderate Reynolds and Prandtl numbers. International Journal of Heat and Mass Transfer 54 (11-12), 2595–2601.
  * Law (2017) Law, Kara Lavender 2017 Plastics in the Marine Environment. Annual Review of Marine Science 9 (1), 205–229.
  * Leal (2012) Leal, L. Gary 2012 Convection Effects in Low-Reynolds-Number Flows. In Advanced Transport Phenomena, , vol. 135, p. 663. Cambridge University Press.
  * Lund & Rogers (1994) Lund, Thomas S. & Rogers, Michael M. 1994 An improved measure of strain state probability in turbulent flows. Physics of Fluids 6 (5), 1838–1847.
  * Ma & Zhao (2020) Ma, Huaqing & Zhao, Yongzhi 2020 Convective heat transfer coefficient for a rod-like particle in a uniform flow. International Journal of Heat and Mass Transfer 147.
  * Masliyah & Epstein (1972) Masliyah, Jacob H. & Epstein, Norman 1972 Numerical solution of heat and mass transfer from spheroids in steady axisymmetric flow. In Progress in Heat and Mass Transfer, Volume 6. Proceedings of the International Symposium on Two-Phase Systems (ed. G. Hetsroni & S. Sideman), pp. 613–632. Oxford: Pergamon.
  * Masoud & Stone (2019) Masoud, Hassan & Stone, Howard A 2019 The reciprocal theorem in fluid dynamics and transport phenomena .
  * Michaelides (2003) Michaelides, Efstathios E. 2003 Hydrodynamic force and heat/mass transfer from particles, bubbles, and drops - The freeman scholar lecture. Journal of Fluids Engineering, Transactions of the ASME 125 (2), 209–238.
  * Myerson (2002) Myerson, Allan 2002 Handbook of Industrial Crystallization. Boston: Butterworth-Heinemann.
  * Pahlow et al. (1997) Pahlow, Markus, Riebesell, Ulf & Wolf-Gladrow, Dieter A. 1997 Impact of cell shape and chain formation on nutrient acquisition by marine diatoms. Limnology and Oceanography 42 (8), 1660–1672.
  * Poe & Acrivos (1976) Poe, G.G. & Acrivos, A. 1976 Closed Streamline Flows Past Small Rotating Particles: Heat Transfer at High Peclet Numbers. Int. J. Multiphase Flow 2 (4), 365–377.
  * Pozrikidis (1997) Pozrikidis, C. 1997 Unsteady heat or mass transport from a suspended particle at low Péclet numbers. Journal of Fluid Mechanics 334, 111–133.
  * Sehlin (1969) Sehlin, R.C. 1969 Forced-Convection Heat and Mass Transfer at Large Peclet Numbers From an Axisymmetric Body in Laminar Flow: Prolate and Oblate Spheroids. PhD thesis, Carnegie-Mellon University, Pittsburgh, PA.
  * Seidensticker et al. (2017) Seidensticker, Sven, Zarfl, Christiane, Cirpka, Olaf A., Fellenberg, Greta & Grathwohl, Peter 2017 Shift in Mass Transfer of Wastewater Contaminants from Microplastics in the Presence of Dissolved Substances. Environmental Science and Technology 51 (21), 12254–12263.
  * Sparrow et al. (2004) Sparrow, Ephraim M., Abraham, John P. & Tong, Jimmy C.K. 2004 Archival correlations for average heat transfer coefficients for non-circular and circular cylinders and for spheres in cross-flow. International Journal of Heat and Mass Transfer 47 (24), 5285–5296.
  * Suhrhoff & Scholz-Böttcher (2016) Suhrhoff, Tim Jesper & Scholz-Böttcher, Barbara M. 2016 Qualitative impact of salinity, UV radiation and turbulence on leaching of organic plastic additives from four common plastics — A lab experiment. Marine Pollution Bulletin 102 (1), 84–94.
  * Szeri (1993) Szeri, A. J. 1993 Pattern formation in recirculating flows of suspensions of orientable particles. Philosophical Transactions - Royal Society of London, A 345 (1677), 477–506.