2.1. Small-oscillation-amplitude limit

The nearly spherical particle executes small-amplitude motion such that the scaled Lagrangian displacements of all material points on its surface are of $O(\delta )$. Consequently, the scaled angular displacements of these material points are also of $O(\delta )$.

Expanding (2.1) to linear order in small $\delta$ gives the (linear) unsteady Stokes equations. This has the consequence that all effects due to nonlinear interactions of either fluid or particle motion (including Lagrangian boundary corrections) are removed, i.e. there are no time-averaged effects or mixing at different frequencies. Due to this linearity, all dependent variables (denoted by $X$) exhibit the explicit time dependence $X(\boldsymbol{r}, t) = {\rm Re}[\tilde {X}(\boldsymbol{r}) \,{\rm e}^{-\mathrm{i} t}]$, where $\boldsymbol{r}$ is the position vector from the particle's geometric centre, ‘$\mathrm{i}$’ is the imaginary unit and $t$ is time. The superfluous $\sim$ notation is omitted for convenience and henceforth all dependent variables refer to their quantities in the frequency domain; note that as per (2.4), $\boldsymbol{U}$ and $\boldsymbol{\varOmega }$ are both quantities in the frequency domain. The unsteady Stokes equations are therefore

(2.3)\begin{equation} \lambda^2 \boldsymbol{u} ={-}\boldsymbol{\nabla} p + \nabla^2 \boldsymbol{u}, \quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \end{equation}

where $\lambda = (1-\mathrm{i})\sqrt {\beta /2}$. In this linear limit, the flow boundary conditions are the velocity of the particle mapped onto the stationary particle surface, i.e.

(2.4)\begin{equation} \boldsymbol{u} = \boldsymbol{U} + \boldsymbol{\varOmega}\times\boldsymbol{r} \quad \text{on}\ \boldsymbol{r}\in S, \end{equation}

where $S$ is the (stationary) surface of the nearly spherical particle, together with decay of the velocity field to zero far from the particle, i.e. (2.2b) holds. Figure 2 shows a cross-section of the particles defined in figure 1, illustrating the difference between the nearly spherical surface, $S$, and the surface of the unit sphere, $\hat {S}$, which is equivalent to the surface of the equivalent-volume sphere after non-dimensionalising by $R_{{eq}}$. Due to linearity of the governing equations/boundary conditions, all results in this study can be trivially extended to arbitrary time-dependent motions of the particle through the use of Fourier/Laplace transforms.

Figure 2.

Cross-sections in the $y_1$–$z_1$, $y_2$–$z_2$ and $y_3$–$z_3$ Cartesian planes of the three particles in figure 1. Solid lines correspond to $S$, the surface of the nearly spherical particle (non-dimensionalised by the radius of an equivalent-volume sphere), and dashed lines correspond to $\hat {S}$, the surface of the unit sphere. See (4.3), (4.11b) and (4.18), with $i=2$, for the shape perturbation functions, $f$ and $g$, that generate $S$ in (a)–(c), respectively. The (infinitesimal) degree of non-sphericity is exaggerated for clarity.

2.2. Nearly spherical particle

The particle is ‘nearly spherical’ and of arbitrary shape with a surface that is specified by the dimensionless radial coordinate:

(2.5)\begin{align} R = 1 + \epsilon f(\theta, \phi) + \epsilon^2 g(\theta, \phi) + O(\epsilon^3). \end{align}

We remind the reader that $R$ and all spatial variables are scaled by the radius of an equivalent-volume sphere, $R_{eq}$. Equation (2.5) defines the shape of the particle and is independent of its motion. The $O(1)$ dimensionless functions $f(\theta, \phi )$ and $g(\theta, \phi )$ define the radial shape perturbation correct to $O(\epsilon ^2)$. The variables $\theta$ and $\phi$ are the usual polar and azimuthal angles, respectively, of the particle's fixed spherical coordinate system (see figure 1) whose origin is at the particle's geometric centre. The parameter $\epsilon$ is the asymptotically small ‘non-spherical parameter’ which is identical to that used by Zhang & Stone (Reference Zhang and Stone1998). The relationship between the two radial shape perturbation functions, $f$ and $g$, is discussed below. Note that for a particle with an inhomogeneous mass density distribution, the particle's geometric centre and centre of mass may not coincide.

The equivalent-volume sphere constraint is

(2.6)\begin{equation} \int_{\hat{S}}\int_0^R \mathrm{d}V = \frac{4{\rm \pi} }{3} , \end{equation}

where $R$ is defined in (2.5). When (2.6) is expanded to first and second order in $\epsilon$, we obtain the respective results

(2.7a)$$\begin{gather} \int_{\hat{S}} f \, \mathrm{d}S = 0, \end{gather}$$

(2.7b)$$\begin{gather}\int_{\hat{S}} g \, \mathrm{d}S ={-}\int_{\hat{S}} f^2 \,\mathrm{d}S , \end{gather}$$

where $\hat {S}$ is the surface of the unit sphere; see § 2.1. Requiring that the origin of the coordinate system is at the particle's geometric centre provides the additional constraint

(2.8)\begin{equation} \int_{\hat{S}}\int_0^R \boldsymbol{r}\, \mathrm{d}V = \boldsymbol{0}, \end{equation}

which when evaluated to first and second order in $\epsilon$ respectively gives

(2.9a)$$\begin{gather} \int_{\hat{S}} f \boldsymbol{n} \, \mathrm{d}S = \boldsymbol{0}, \end{gather}$$

(2.9b)$$\begin{gather}\int_{\hat{S}} g \boldsymbol{n} \, \mathrm{d}S ={-}\frac{3}{2}\int_{\hat{S}} f^2 \boldsymbol{n} \, \mathrm{d}S, \end{gather}$$

where $\boldsymbol{n}$ is the unit normal to $\hat {S}$, directed into the fluid domain. Directives for specifying the functions $f$ and $g$ to ensure satisfaction of (2.7) and (2.9) are given in Appendix B.

We note that the surface perturbation in (2.5) could alternatively be described using $R(\theta,\phi ) = 1+f(\theta,\phi )$. This would generate a different but equivalent asymptotic theory. However, the benefit of using (2.5) is that it naturally enables a spherical harmonic expansion of the surface at each order of $\epsilon$. This is specifically advantageous for particle shapes described using a finite number of terms in their spherical harmonic expansion.

Regime of validity: Importantly, the small-oscillation-amplitude asymptotic expansion to $O(\delta )$ is taken before the near-spherical-particle asymptotic expansion to $O(\epsilon ^2)$. It follows that the results reported in this study apply rigorously when the oscillation amplitude, $\delta$, is much smaller than the height of any surface perturbation on the particle, i.e. $\delta \ll \epsilon$. We note that the opposite limit ($\delta \gg \epsilon$) requires inclusion of nonlinear convective inertial effects in the flow – it will generate nonlinear flow phenomena such as acoustic streaming – which is not the focus of this study; a brief illustration of this effect is reported in § 5.

2.3. Flow generated by a nearly spherical particle

We expand the flow variables in the (infinitesimal) non-spherical parameter $\epsilon$:

(2.10)\begin{equation} (\boldsymbol{u}, p) = (\boldsymbol{u}^{(0)}, p^{(0)}) + \epsilon (\boldsymbol{u}^{(1)}, p^{(1)}) + \epsilon^2 (\boldsymbol{u}^{(2)}, p^{(2)}) + O(\epsilon^3). \end{equation}

Substituting (2.10) into (2.4), and performing a Taylor expansion about $r = 1$, then gives

(2.11)\begin{align} &\boldsymbol{u}^{(0)}|_{r=1} +\epsilon(\boldsymbol{u}^{(1)}|_{r=1} + f \partial_r \boldsymbol{u}^{(0)}|_{r=1} ) \nonumber\\ &\qquad + \epsilon^2\left(\boldsymbol{u}^{(2)}|_{r=1} + g \partial_r \boldsymbol{u}^{(0)}|_{r=1} + f \partial_r \boldsymbol{u}^{(1)}|_{r=1} + \frac{f^2}{2} \partial^2_r \boldsymbol{u}^{(0)}|_{r=1} \right) + O(\epsilon^3) \nonumber\\ &\quad =\boldsymbol{U} + \boldsymbol{\varOmega}\times\boldsymbol{n} + \epsilon f \boldsymbol{\varOmega}\times\boldsymbol{n} + \epsilon^2 g \boldsymbol{\varOmega}\times\boldsymbol{n} + O(\epsilon^3), \end{align}

where $\partial _r$ and $\partial ^2_r$ are the first and second partial derivatives with respect to the radial coordinate, $r$, respectively. Note that the position vector, $\boldsymbol{r}$, is identical to the unit normal vector, $\boldsymbol{n}$, when evaluation occurs on the surface of the unit sphere, $\hat {S}$. Equation (2.11) then specifies the required flow boundary conditions on the nearly spherical particle surface (evaluated at the sphere surface, $r=1$), at each order in $\epsilon$:

(2.12a)$$\begin{gather} \boldsymbol{u}^{(0)}|_{r=1} = \boldsymbol{U} + \boldsymbol{\varOmega}\times\boldsymbol{n}, \end{gather}$$

(2.12b)$$\begin{gather}\boldsymbol{u}^{(1)}|_{r=1} = f \boldsymbol{\varOmega}\times\boldsymbol{n} -f \partial_r \boldsymbol{u}^{(0)}|_{r=1}, \end{gather}$$

(2.12c)$$\begin{gather}\boldsymbol{u}^{(2)}|_{r=1} = g \boldsymbol{\varOmega}\times\boldsymbol{n} -g \partial_r \boldsymbol{u}^{(0)}|_{r=1} - f \partial_r \boldsymbol{u}^{(1)}|_{r=1} - \frac{f^2}{2} \partial^2_r \boldsymbol{u}^{(0)}|_{r=1}, \end{gather}$$

with the individual flows vanishing far from the particle, as per (2.2b).

2.4. Force and torque acting on a nearly spherical particle

For a steady Stokes flow, the resulting force and torque acting on a body can be evaluated by integrating the stress tensor, $\boldsymbol{\sigma }$, on any surface that encloses the body (including any singular points of the flow), because the stress tensor is divergence free (Brenner Reference Brenner1964a). The situation is more involved for unsteady Stokes flows where the stress tensor is not divergence free; the stress tensor must be integrated on the body's surface. Extending the approach of Zhang & Stone (Reference Zhang and Stone1998) to $O(\epsilon ^2)$ gives the following respective expressions for the force and torque exerted by the fluid on the nearly spherical particle:

(2.13a)$$\begin{gather} \boldsymbol{F} = \boldsymbol{F}^{(0)} + \epsilon \boldsymbol{F}^{(1)} + \epsilon^2 \boldsymbol{F}^{(2)} + O(\epsilon^3), \end{gather}$$

(2.13b)$$\begin{gather}\boldsymbol{T} = \boldsymbol{T}^{(0)} + \epsilon \boldsymbol{T}^{(1)} + \epsilon^2 \boldsymbol{T}^{(2)} + O(\epsilon^3), \end{gather}$$

with

(2.14a)$$\begin{gather} \boldsymbol{F}^{(0)} = \int_{\hat{S}} \boldsymbol{n} \boldsymbol{\cdot}\boldsymbol{\sigma}^{(0)}\, \mathrm{d}S, \end{gather}$$

(2.14b)$$\begin{gather}\boldsymbol{F}^{(1)} = \int_{\hat{S}} \boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\sigma}^{(1)} \,\mathrm{d}S, \end{gather}$$

(2.14c)$$\begin{gather}\boldsymbol{F}^{(2)} = \int_{\hat{S}} \boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\sigma}^{(2)} + \frac{\lambda^2 f^2}{2}(\boldsymbol{\varOmega}\times \boldsymbol{n} - \partial_r \boldsymbol{u}^{(0)})\,\mathrm{d}S, \end{gather}$$

and

(2.14d)$$\begin{gather} \boldsymbol{T}^{(0)} = \int_{\hat{S}} \boldsymbol{n}\times (\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\sigma}^{(0)}) \,\mathrm{d}S, \end{gather}$$

(2.14e)$$\begin{gather}\boldsymbol{T}^{(1)} = \int_{\hat{S}} \boldsymbol{n}\times (\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\sigma}^{(1)} + \lambda^2 f\boldsymbol{\varOmega}\times\boldsymbol{n} ) \,\mathrm{d}S, \end{gather}$$

(2.14f)$$\begin{gather}\boldsymbol{T}^{(2)} = \int_{\hat{S}} \boldsymbol{n}\times\left(\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\sigma}^{(2)} + \lambda^2\left[\left(\frac{5}{2}f^2 + g\right) \boldsymbol{\varOmega}\times\boldsymbol{n} - \frac{f^2}{2}\partial_r \boldsymbol{u}^{(0)}\right]\right) \mathrm{d}S, \end{gather}$$

where the individual stress tensors, $\boldsymbol{\sigma }^{(0)}$, $\boldsymbol{\sigma }^{(1)}$, $\boldsymbol{\sigma }^{(2)}$, are defined by the flows generated by the boundary conditions in (2.12a), (2.12b), (2.12c), respectively. As required, all the integrals in (2.14) are specified over $\hat {S};$ see Appendix C for the derivation of the above results.

In principle, (2.14) requires knowledge of the above-mentioned individual stress tensors. Use of the Lorentz reciprocal theorem bypasses this requirement by enabling determination of the integrals in (2.14) using only their respective velocity fields. For an arbitrary unsteady Stokes flow, $(\boldsymbol{u}', \boldsymbol{\sigma }')$, this gives

(2.15a)$$\begin{gather} \int_{\hat{S}} \boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\sigma}' \,\mathrm{d}S ={-}\frac{1}{2}\int_{\hat{S}} \boldsymbol{u}'\boldsymbol{\cdot}[3(1+\lambda)(\boldsymbol{I} - \boldsymbol{n}\boldsymbol{n}) + (3+3\lambda+\lambda^2)\boldsymbol{n}\boldsymbol{n} ] \,\mathrm{d}S, \end{gather}$$

(2.15b)$$\begin{gather}\int_{\hat{S}} \boldsymbol{n} \times (\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\sigma}') \,\mathrm{d}S ={-}\frac{3+3\lambda+\lambda^2}{1+\lambda}\int_{\hat{S}} \boldsymbol{n}\times\boldsymbol{u}' \,\mathrm{d}S , \end{gather}$$

which are derived by substituting this arbitrary flow, and flows generated by the translation and rotation (of unitary velocity magnitude) of a unit sphere, into the Lorentz reciprocal theorem.

3.1. Zeroth-order solution

The zeroth-order solution in $\epsilon$, i.e. $O(1)$, corresponds to the force and torque experienced by a perfect sphere (Stokes Reference Stokes1851). The fluid's velocity field and its radial gradients at the sphere's surface, to be used in the subsequent first- and second-order solutions, are

(3.1a)$$\begin{gather} \boldsymbol{u}^{(0)}|_{r=1} = \boldsymbol{U} + \boldsymbol{\varOmega}\times\boldsymbol{n}, \end{gather}$$

(3.1b)$$\begin{gather}\partial_r \boldsymbol{u}^{(0)}|_{r=1} ={-}\frac{3}{2}(1+\lambda)(\boldsymbol{I} - \boldsymbol{n}\boldsymbol{n}) \boldsymbol{\cdot}\boldsymbol{U} - \frac{2+2\lambda+\lambda^2}{1+\lambda}\boldsymbol{\varOmega}\times\boldsymbol{n}, \end{gather}$$

(3.1c)\begin{gather} \partial^2_r \boldsymbol{u}^{(0)}|_{r=1} = \left(\frac{3}{2}(3+3\lambda+\lambda^2) (\boldsymbol{I}-\boldsymbol{n}\boldsymbol{n}) - 3(1+\lambda)\boldsymbol{n}\boldsymbol{n}\right) \boldsymbol{\cdot}\boldsymbol{U} \nonumber\\ +\, \frac{6+6\lambda + 3\lambda^2 + \lambda^3}{1+\lambda}\boldsymbol{\varOmega}\times\boldsymbol{n}. \end{gather}

The force and torque are evaluated directly using (2.15) and yield the well known results

(3.2a)$$\begin{gather} \boldsymbol{F}^{(0)} ={-}6{\rm \pi} \left(1+\lambda+\frac{\lambda^2}{9}\right)\boldsymbol{U}, \end{gather}$$

(3.2b)$$\begin{gather}\boldsymbol{T}^{(0)} ={-}8{\rm \pi} \left(\frac{1+\lambda+\dfrac{\lambda^2}{3}}{1+\lambda}\right)\boldsymbol{\varOmega}. \end{gather}$$

3.2. First-order solution

The first-order velocity boundary condition in $\epsilon$ is obtained by substituting (3.1b) into (2.12b):

(3.3)\begin{equation} \boldsymbol{u}^{(1)}|_{r=1} = \left(\frac{3}{2}(1+\lambda)(\boldsymbol{I} - \boldsymbol{n}\boldsymbol{n})\boldsymbol{\cdot}\boldsymbol{U} + \frac{3+3\lambda + \lambda^2}{1+\lambda}\boldsymbol{\varOmega}\times\boldsymbol{n}\right) f. \end{equation}

Equations (2.14b), (2.15a) and (3.3) then give

(3.4)\begin{equation} \boldsymbol{F}^{(1)} = \frac{9}{4}(1+\lambda)^2 \int_{\hat{S}} f \boldsymbol{n}\boldsymbol{n} \, \mathrm{d}S \boldsymbol{\cdot} \boldsymbol{U}, \end{equation}

while (2.14e), (2.15b), (3.1a) and (3.3) yield

(3.5)\begin{equation} \boldsymbol{T}^{(1)} = \frac{(3+2\lambda)(3+4\lambda+2\lambda^2)}{(1+\lambda)^2}\int_{\hat{S}} f \boldsymbol{nn} \, \mathrm{d}S\boldsymbol{\cdot} \boldsymbol{\varOmega}. \end{equation}

The formulae in (3.4) and (3.5) coincide with those of Zhang & Stone (Reference Zhang and Stone1998). Expressing $f (\theta, \phi )$ generally as a series involving spherical harmonics (see Appendix B) establishes that only the spherical harmonics of degree $l = 2$ provide a non-zero contribution to both (3.4) and (3.5). That is, the force and torque to $O(\epsilon )$ experienced by a nearly spherical particle are dictated by those of an equivalent ellipsoid (Zhang & Stone Reference Zhang and Stone1998).

Equations (3.4) and (3.5) show that the force and torque (about the particle's geometric centre) are decoupled. The particle experiences no drag force due to its rotational motion, and no torque due to its linear translation. This does not agree with general physical reality and, as such, coupling between translational and rotational motion/forces must occur at higher order in $\epsilon$, as discussed in § 1. Additionally, the $\lambda ^2$ coefficient of $\boldsymbol{T}^{(1)}$ is zero as $|\lambda |\rightarrow \infty$, again in disagreement with reality, i.e. the high-frequency limit mentioned in § 1.

3.3. Second-order solution

We now evaluate the second-order solution in $\epsilon$ for the force and torque, which is the principal aim of this study. Details of this derivation are relegated to Appendix D. These solutions are decomposed into two components: (1) a component involving only $\boldsymbol{u}^{(0)}$ and (2) a component that depends on $\partial _r \boldsymbol{u}^{(1)}$ which is expressed as an infinite series involving the shape perturbation function, $f$. These two components are denoted by subscripts 1 and 2, respectively:

(3.6a,b)\begin{align} \boldsymbol{F}^{(2)} = \boldsymbol{F}^{(2)}_1 + \boldsymbol{F}^{(2)}_2, \quad \boldsymbol{T}^{(2)} = \boldsymbol{T}^{(2)}_1 + \boldsymbol{T}^{(2)}_2, \end{align}

where expressions for each component of $\boldsymbol{F}^{(2)}$ and $\boldsymbol{T}^{(2)}$ follow from (2.12c), (2.14c) and (2.14f), and are detailed in (D4a), (D4b), (D5a) and (D5b). The first components, that involve only $\boldsymbol{u}^{(0)}$, are easily evaluated and give

(3.7)\begin{align} \boldsymbol{F}^{(2)}_1 &= \left[ \frac{3}{8}(1+\lambda)(9 + 9\lambda + 5\lambda^2) \int_{\hat{S}} f^2 (\boldsymbol{I} - \boldsymbol{n}\boldsymbol{n}) \, \mathrm{d}S \right. \nonumber\\ & \quad \left. -\,\frac{3}{4}(1+\lambda)(3 + 3\lambda + \lambda^2) \int_{\hat{S}} f^2 \boldsymbol{n}\boldsymbol{n} \, \mathrm{d}S - \frac{9}{4}(1+\lambda)^2 \int_{\hat{S}} g (\boldsymbol{I} - \boldsymbol{n}\boldsymbol{n}) \, \mathrm{d}S \right] \boldsymbol{\cdot}\boldsymbol{U} \nonumber\\ & \quad - \frac{45 + 90\lambda + 69\lambda^2 + 27\lambda^3 + 5\lambda^4}{4(1+\lambda)}\int_{\hat{S}} f^2 \boldsymbol{n} \, \mathrm{d}S\times\boldsymbol{\varOmega} \end{align}

and

(3.8)\begin{align} \boldsymbol{T}^{(2)}_1 &= \frac{3(3+2\lambda)(6+ 8\lambda + 4\lambda^2 + \lambda^3)}{4(1+\lambda)}\int_{\hat{S}} f^2 \boldsymbol{n} \, \mathrm{d}S \times \boldsymbol{U} \nonumber\\ &\quad + \left[\frac{9 + 18\lambda + 20\lambda^2 + 16\lambda^3 + 7\lambda^4 + \lambda^5}{(1+\lambda)^2}\int_{\hat{S}} f^2 (\boldsymbol{I} - \boldsymbol{n}\boldsymbol{n}) \, \mathrm{d}S \right.\nonumber\\ &\quad \left. -\,\frac{(3+2\lambda)(3+4\lambda+2\lambda^2)}{(1+\lambda)^2}\int_{\hat{S}} g (\boldsymbol{I} - \boldsymbol{n}\boldsymbol{n}) \, \mathrm{d}S \right]\boldsymbol{\cdot}\boldsymbol{\varOmega}. \end{align}

Derivations of the corresponding second components, that depend on $\partial _r \boldsymbol{u}^{(1)}$, are more involved and yield

(3.9)\begin{align} \boldsymbol{F}^{(2)}_2 &= \sum_{l=1}^{\infty}\sum_{m={-}l}^{l} \left[ \kappa_1 \kappa_2 \gamma_l^{(1)} \int_{\hat{S}} f\boldsymbol{R}_{lm}\, \mathrm{d}S \int_{\hat{S}} f \boldsymbol{\boldsymbol{\varTheta}}_{lm}^{*} \, \mathrm{d}S\right. \nonumber\\ & \quad \left. +\, \kappa_2^2 \gamma_l^{(2)} \int_{\hat{S}} f \boldsymbol{\boldsymbol{\varTheta}}_{lm} \, \mathrm{d}S \int_{\hat{S}} f \boldsymbol{\boldsymbol{\varTheta}}_{lm}^{*} \, \mathrm{d}S + \kappa_2^2 \gamma_l^{(3)} \int_{\hat{S}} f \boldsymbol{\boldsymbol{\varPhi}}_{lm} \, \mathrm{d}S \int_{\hat{S}} f \boldsymbol{\boldsymbol{\varPhi}}_{lm}^{*} \, \mathrm{d}S \right] \boldsymbol{\cdot} \boldsymbol{U} \nonumber\\ & \quad + \left[ \kappa_1 \kappa_3 \gamma_l^{(1)} \int_{\hat{S}} f \boldsymbol{R}_{lm} \, \mathrm{d}S \int_{\hat{S}} f \boldsymbol{\boldsymbol{\varPhi}}_{lm}^{*} \, \mathrm{d}S + \kappa_2 \kappa_3 \gamma_l^{(2)} \int_{\hat{S}} f \boldsymbol{\boldsymbol{\varTheta}}_{lm} \, \mathrm{d}S \int_{\hat{S}} f \boldsymbol{\boldsymbol{\varPhi}}_{lm}^{*} \, \mathrm{d}S \right. \nonumber\\ & \quad \left. -\, \kappa_2 \kappa_3 \gamma_l^{(3)} \int_{\hat{S}} f \boldsymbol{\boldsymbol{\varPhi}}_{lm} \, \mathrm{d}S \int_{\hat{S}} f \boldsymbol{\boldsymbol{\varTheta}}_{lm}^{*}\, \mathrm{d}S\right] \boldsymbol{\cdot} \boldsymbol{\varOmega} \end{align}

and

(3.10)\begin{align} \boldsymbol{T}^{(2)}_2 &= \sum_{l=1}^{\infty} \sum_{m={-}l}^{l} \left[ - \kappa_2 \kappa_3 \gamma_l^{(3)} \int_{\hat{S}} f \boldsymbol{\boldsymbol{\varTheta}}_{lm} \, \mathrm{d}S \int_{\hat{S}} f \boldsymbol{\boldsymbol{\varPhi}}_{lm}^{*} \, \mathrm{d}S + \kappa_2 \kappa_3 \gamma_l^{(2)} \int_{\hat{S}} f \boldsymbol{\boldsymbol{\varPhi}}_{lm} \, \mathrm{d}S \int_{\hat{S}} f \boldsymbol{\boldsymbol{\varTheta}}_{lm}^{*} \, \mathrm{d}S \right]\boldsymbol{\cdot}\boldsymbol{U}\nonumber\\ & \quad + \left[ \kappa_3^2 \gamma_l^{(3)} \int_{\hat{S}} f \boldsymbol{\boldsymbol{\varTheta}}_{lm} \, \mathrm{d}S \int_{\hat{S}} f \boldsymbol{\boldsymbol{\varTheta}}_{lm}^{*} \, \mathrm{d}S + \kappa_3^2 \gamma_l^{(2)} \int_{\hat{S}} f \boldsymbol{\boldsymbol{\varPhi}}_{lm} \, \mathrm{d}S \int_{\hat{S}} f \boldsymbol{\boldsymbol{\varPhi}}^{*}_{lm} \, \mathrm{d}S \right]\boldsymbol{\cdot} \boldsymbol{\varOmega}, \end{align}

where $\boldsymbol{R}_{lm}$, $\boldsymbol{\boldsymbol{\varTheta }}_{lm}$, $\boldsymbol{\boldsymbol{\varPhi }}_{lm}$ are vector spherical harmonics defined in Appendix D; the outer product is implied between the pairs of integrals in (3.9) and (3.10); the coefficients are

(3.11) \begin{equation} \left.\begin{array}{c@{}} \displaystyle\gamma_l^{(1)} = \sqrt{l(l+1)},\\ \displaystyle \gamma_l^{(2)} = \dfrac{\lambda^2}{2l+1}\left[1 + \dfrac{\mathrm{h}^{(1)}_{l+1}(\mathrm{i}\lambda)}{\mathrm{h}^{(1)}_{l-1}(\mathrm{i}\lambda)}\right]- 1, \\ \displaystyle \gamma_l^{(3)} = l - \mathrm{i}\lambda\dfrac{\mathrm{h}^{(1)}_{l+1}(\mathrm{i}\lambda)}{\mathrm{h}^{(1)}_{l}(\mathrm{i}\lambda)}, \end{array}\right\} \end{equation}

where $\mathrm{h}^{(1)}_l$ is the spherical Hankel function of the first kind of degree $l$; and

(3.12a––c)\begin{align} \kappa_1 = \frac{1}{2}(3 + 3\lambda + \lambda^2), \quad \kappa_2 = \frac{3}{2}(1+\lambda), \quad \kappa_3 = \frac{3 + 3\lambda + \lambda^2}{1+\lambda}. \end{align}

The (first) infinite sums in (3.9) and (3.10) begin with $l = 1$ because the $l = 0$ terms are identically zero, due to the respective equivalent-volume and geometric-centre requirements in (2.7a) and (2.9a). Equations (3.6)–(3.12) specify the complete second-order solution and require only the shape perturbation functions, $f (\theta, \phi )$ and $g(\theta, \phi )$. In contrast to the first-order solution in the previous section, where only the $l=2$ spherical harmonics of $f(\theta, \phi )$ contribute, all spherical harmonics with $l \ge 2$ can affect the second-order solution; see § B.1.

3.3.1. Infinite series in (3.9) and (3.10)

The number of non-zero terms in the (first) infinite series, over index $l$, in (3.9) and (3.10), is dictated by the nature of the shape perturbation function, $f$. Suppose, $f$ can be expressed as a finite linear combination of $k$ spherical harmonics:

(3.13)\begin{equation} f(\theta, \phi) = \sum_{i = 1}^k f_{l_i m_i} Y_{l_i m_i}(\theta,\phi), \end{equation}

where $l_i$ and $m_i$ are specified integers, i.e. the number of non-zero terms in (B1a) is finite. For each $i$ in (3.13), the only terms in (3.9) and (3.10) that may be non-zero are

(3.14)\begin{equation} l_i - 1 \le l \le l_i + 1, \end{equation}

with

(3.15a,b)\begin{align} {-}|m_i| - 1 \le m \le -|m_i| + 1 \quad \text{or} \quad |m_i| - 1 \le m \le |m_i| + 1. \end{align}

Moreover, let $l = l_{max}$ be the largest degree, $l$, in (3.13). It follows from (3.14) that all terms in (3.9) and (3.10) for $l > l_{max}+1$ are identically zero. The infinite series in (3.9) and (3.10) then become finite series from $l = 1$ to $l_{maxsum}\, (\equiv l_{max} +1)$. This property is used in the Mathematica notebook provided in the supplementary material, with the notebook variable $\textrm{MAXNUM} \equiv l_{maxsum}$; see Appendix A.

3.4. Orientation-averaged solution

It is of practical relevance to study the orientation-averaged force and torque experienced by nearly spherical particles, i.e. the expected value of the force (or torque) exerted by the fluid on a particle that is oscillating along (or around) a random axis. This can arise in experiments when the particle orientation is not controlled, e.g. measurements of the particle's Brownian motion, where all particle orientations are sampled equally. The orientation-averaged force, $\bar {\boldsymbol{F}}$, and torque, $\bar {\boldsymbol{T}}$, are determined by taking a uniformly distributed ensemble of all possible orientations, and calculating the average of the resulting forces and torques, i.e. for given $\boldsymbol{U}$ and $\boldsymbol{\varOmega }$,

(3.16a)$$\begin{gather} \bar{\boldsymbol{F}} = \frac{1}{8{\rm \pi} ^2}\int_{0}^{2{\rm \pi} }\int_{0}^{{\rm \pi} } \int_{0}^{2{\rm \pi} } \boldsymbol{F}|_{(\psi',\theta',\phi')} \sin\theta' \,{\rm d}\psi'\,{\rm d}\theta'\,{\rm d}\phi', \end{gather}$$

(3.16b)$$\begin{gather}\bar{\boldsymbol{T}} = \frac{1}{8{\rm \pi} ^2}\int_{0}^{2{\rm \pi} }\int_{0}^{{\rm \pi} } \int_{0}^{2{\rm \pi} } \boldsymbol{T}|_{(\psi',\theta',\phi')} \sin\theta' \,{\rm d}\psi'\,{\rm d}\theta'\,{\rm d}\phi', \end{gather}$$

where $(\psi ',\theta ',\phi ')$ are a set of Euler angles with respect to the laboratory frame. It can be shown that this procedure gives

(3.17a)$$\begin{gather} \bar{\boldsymbol{F}} = \left[{-}6{\rm \pi} \left(1 + \lambda + \frac{\lambda^2}{9}\right) + \frac{\epsilon^2}{3} \, \mathrm{tr}\, \boldsymbol{A}^{(2)} \right]\boldsymbol{U} + O(\epsilon^3), \end{gather}$$

(3.17b)$$\begin{gather}\bar{\boldsymbol{T}} = \left[{-}8{\rm \pi} \left(\frac{1+\lambda+\dfrac{\lambda^2}{3}}{1+\lambda}\right) + \frac{\epsilon^2}{3}\, \mathrm{tr} \, \boldsymbol{C}^{(2)}\right]\boldsymbol{\varOmega} + O(\epsilon^3), \end{gather}$$

where $\mathrm{tr}\,\boldsymbol{A}^{(2)}$ and $\mathrm{tr}\,\boldsymbol{C}^{(2)}$ are the trace of the tensorial coefficients of $\boldsymbol{U}$ in (3.6a) and $\boldsymbol{\varOmega }$ in (3.6b), respectively. Explicit forms for $\mathrm{tr}\,\boldsymbol{A}^{(2)}$ and $\mathrm{tr}\,\boldsymbol{C}^{(2)}$ are given by (D16a) and (D16b), respectively; see § D.3. As discussed in § 1, particle non-sphericity enters into the orientation averages at $O(\epsilon ^2)$.

4.1. Spheroidal particle

First, we consider a prolate spheroid (illustrated in figures 1a and 2a) with its axis of rotational symmetry aligned along the Cartesian $z$ axis, with basis vector $\hat {\boldsymbol{z}}$. The dimensionless length of the spheroid's major axis is chosen to be $1+\epsilon$. Requiring that its volume is independent of $\epsilon$ and equal to that of an equivalent-volume sphere, i.e. $4{\rm \pi} /3$ as per (2.5) and (2.7), gives a minor-axis length of $1/\sqrt {1+\epsilon }$. The corresponding radial coordinate, for arbitrary $\epsilon$, of the spheroid's surface is then

(4.1)\begin{equation} R = \frac{1+\epsilon}{\sqrt{\cos^2\theta + (1+\epsilon)^3\sin^2\theta}} , \end{equation}

which when expanded in small $\epsilon$ gives

(4.2)\begin{equation} R = 1 + \epsilon f(\theta, \phi) + \epsilon^2 g(\theta, \phi) + O(\epsilon^3), \end{equation}

with

(4.3a)$$\begin{gather} f(\theta, \phi) = 2\sqrt{\frac{{\rm \pi} }{5}}{Y}_{2,0}(\theta, \phi), \end{gather}$$

(4.3b)$$\begin{gather}g(\theta, \phi) ={-}\frac{2\sqrt{{\rm \pi} }}{35}(7 {Y}_{0,0}(\theta, \phi) -4\sqrt{5} {Y}_{2,0}(\theta, \phi) + 9 {Y}_{4,0}(\theta, \phi)), \end{gather}$$

where $Y_{l,m} (\theta, \phi )$ is the scalar spherical harmonic defined in Appendix D. Note that $g(\theta,\phi )$ contains an $(l,m) = (4,0)$ spherical harmonic to exactly represent the spheroid; this does not affect the force and torque, as discussed in § B.1.

4.1.1. First-order solution

The leading-order effects of non-sphericity are given by the first-order solutions in $\epsilon$ for the force and torque, which follow directly from (3.4) and (3.5):

(4.4a)$$\begin{gather} \boldsymbol{F}^{(1)} ={-}\frac{3{\rm \pi} }{5}(1 + \lambda)^2(\boldsymbol{I} - 3\hat{\boldsymbol{z}}\hat{\boldsymbol{z}})\boldsymbol{\cdot} \boldsymbol{U}, \end{gather}$$

(4.4b)$$\begin{gather}\boldsymbol{T}^{(1)} ={-}\frac{4{\rm \pi} (3+2\lambda)(3+4\lambda+2\lambda^2)}{15(1+\lambda)^2}(\boldsymbol{I} - 3\hat{\boldsymbol{z}}\hat{\boldsymbol{z}})\boldsymbol{\cdot} \boldsymbol{\varOmega}. \end{gather}$$

It is evident that the corresponding orientation-averaged first-order solutions are zero for both the force and the torque, as expected.

Force

The first component of the second-order force follows immediately from (3.7):

(4.5)\begin{equation} \boldsymbol{F}^{(2)}_1 = \frac{{\rm \pi} }{35}(1+\lambda)[(111 + 111\lambda + 35\lambda^2)(\boldsymbol{I} - \hat{\boldsymbol{z}}\hat{\boldsymbol{z}} ) + (30+30\lambda+14\lambda^2)\hat{\boldsymbol{z}}\hat{\boldsymbol{z}}]\boldsymbol{\cdot} \boldsymbol{U}. \end{equation}

The radial shape perturbation function, $f$, in (4.3) involves spherical harmonics of degree $l \le l_{max} = 2$. This facilitates exact evaluation (with $l_{maxsum} \equiv l_{max}+1=3$; see § 3.3.1) of the second component of the force in (3.9), giving

(4.6)\begin{align} \boldsymbol{F}^{(2)}_2 &={-}\frac{{\rm \pi} (1+\lambda)}{350(3+3\lambda + \lambda^2)}[(4383 + 8829\lambda + 6039\lambda^2 + 2073\lambda^3 + 350\lambda^4)(\boldsymbol{I} - \hat{\boldsymbol{z}}\hat{\boldsymbol{z}}) \nonumber\\ &\quad + (2952 + 6156\lambda + 4176\lambda^2 + 1272\lambda^3 + 140\lambda^4)\hat{\boldsymbol{z}}\hat{\boldsymbol{z}}]\boldsymbol{\cdot} \boldsymbol{U}. \end{align}

The second-order solution for the force is then obtained by substituting (4.5) and (4.6) into (3.6), yielding the required result:

(4.7)\begin{align} \boldsymbol{F}^{(2)} &={-}\frac{3{\rm \pi} (1+\lambda)}{350(3+3\lambda + \lambda^2)}[(351 + 723\lambda + 183\lambda^2 - 29\lambda^3)(\boldsymbol{I} - \hat{\boldsymbol{z}}\hat{\boldsymbol{z}} )\nonumber\\ &\quad + (684 + 1452\lambda + 852\lambda^2 + 184\lambda^3)\hat{\boldsymbol{z}}\hat{\boldsymbol{z}}]\boldsymbol{\cdot} \boldsymbol{U}. \end{align}

Equation (4.7) is to be compared with known limiting cases from the literature. Namely, Chwang & Wu (Reference Chwang and Wu1975) calculated the second-order force and torque experienced by a prolate spheroid undergoing steady, but arbitrary, translation and rotation. In the steady limit, i.e. $|\lambda | \rightarrow 0$, (4.7) coincides with the exact solution of Chwang & Wu (Reference Chwang and Wu1975). Moreover, for translational motion of the spheroid along its major axis – which extracts the $\hat {\boldsymbol{z}}\hat {\boldsymbol{z}}$ component of (4.7) – the solution here agrees with the corresponding result reported by Lawrence & Weinbaum (Reference Lawrence and Weinbaum1986).

Torque

Repeating this calculation for the torque in (3.6) gives

(4.8)\begin{align} \boldsymbol{T}^{(2)} &={-}\frac{4{\rm \pi} }{105(1+\lambda)^3(15+15\lambda+6\lambda^2+\lambda^3)} [(3861 + 15444\lambda + 27996\lambda^2 + 29391\lambda^3 \nonumber\\ &\quad + 19311\lambda^4 + 8107\lambda^5 + 2130\lambda^6 + 321\lambda^7 + 21\lambda^8)(\boldsymbol{I} - \hat{\boldsymbol{z}}\hat{\boldsymbol{z}}) + 2(1242 + 4968\lambda \nonumber\\ &\quad + 8733\lambda^2 + 8661\lambda^3 + 5196\lambda^4 + 1889\lambda^5 + 390\lambda^6 + 36\lambda^7)\hat{\boldsymbol{z}}\hat{\boldsymbol{z}} ]\boldsymbol{\cdot}\boldsymbol{\varOmega}. \end{align}

We again find that the exact solution of Chwang & Wu (Reference Chwang and Wu1975) is recovered from (4.8) in the steady limit, $|\lambda | \rightarrow 0$, as required.

It is notable that the particle's translation–rotation and its resulting torque–force do not couple, respectively, at first or second order in $\epsilon$. That is, the force depends only on the translational velocity of the particle whereas the torque (about the particle's geometric centre) relies only on its angular velocity. Indeed, this special result is true for all orders in $\epsilon$, and is due to the particle's fore–aft geometric symmetry.

4.1.3. Orientation-averaged force and torque

The orientation-averaged force and torque (defined in § 3.4) are calculated using (3.17), (4.7) and (4.8), giving

(4.9a)\begin{align} &\boldsymbol{\bar{F}} ={-}6{\rm \pi} \boldsymbol{U} \left[1+\lambda+\frac{\lambda^2}{9} + \epsilon^2 \frac{3(1+\lambda)(33+69\lambda+29\lambda^2+3\lambda^3)}{50(3+3\lambda+\lambda^2)}\right] + O(\epsilon^3), \end{align}

(4.9b)\begin{align} & \boldsymbol{\bar{T}} ={-}8{\rm \pi} \boldsymbol{\varOmega} \left[\frac{1+\lambda+\dfrac{\lambda^2}{3}}{1+\lambda}\right. \nonumber\\ &\qquad \left. +\, \epsilon^2 \frac{243 + 972\lambda + 1749\lambda^2 + 1812\lambda^3 + 1167\lambda^4 + 476\lambda^5 + 120\lambda^6 + 17\lambda^7 + \lambda^8}{15(1+\lambda)^3(15+15\lambda+6\lambda^2+\lambda^3)} \vphantom{\boldsymbol{\bar{T}} ={-}8{\rm \pi} \boldsymbol{\varOmega} \left[\frac{1+\lambda+\dfrac{\lambda^2}{3}}{1+\lambda}\right.}\right] + O(\epsilon^3). \end{align}

Interestingly, we observe that particle non-sphericity increases the magnitude of the orientation-averaged force and torque.

4.2. ‘Pear-shaped’ particle

Next, we examine a particle that does not possess fore–aft symmetry which, in principle, enables coupling of the translational force that the particle experiences with its rotational motion. The shape perturbation function, $f$, is specified using a linear combination of $l = 2$ and 3 spherical harmonics, with $m = 0$ for each harmonic; this produces an axisymmetric particle. The result is a particle that resembles a pear; see figures 1(b) and 2(b).

The coefficient of the $l = 2$ spherical harmonic of $f$ is chosen to be identical to that of the spheroid studied in § 4.1. The first-order forces and torques here are identical between the two particles and given by (4.4a) and (4.4b), respectively. The corresponding coefficient of the $l = 3$ harmonic is then chosen to be equal to that of the $l = 2$ harmonic. The second-order shape perturbation function, $g$, is specified using the result in Appendix B with coefficients of the $l=2$ components set to zero. This gives the radial coordinate of the particle's surface:

(4.10)\begin{equation} {R} = 1 + \epsilon f(\theta, \phi) + \epsilon^2 g(\theta, \phi) + O(\epsilon^3), \end{equation}

with

(4.11a)$$\begin{gather} f(\theta, \phi) = 2\sqrt{\frac{{\rm \pi} }{5}}({Y}_{2,0}(\theta, \phi) + {Y}_{3,0}(\theta, \phi)), \end{gather}$$

(4.11b)$$\begin{gather}g(\theta, \phi) ={-}\frac{2\sqrt{{\rm \pi} }}{5}\left(2{Y}_{0,0}(\theta, \phi) + 9\sqrt{\frac{3}{35}}{Y}_{1,0}(\theta, \phi)\right) , \end{gather}$$

where the spherical harmonics are defined in Appendix D.

4.2.1. Second-order solution

The second-order contributions to the force and torque are evaluated from (3.6)–(3.12), again using the Mathematica notebook (supplementary material and Appendix A), giving

(4.12a)$$\begin{gather} \boldsymbol{F}^{(2)} ={-}{\rm \pi} \left[\frac{a_1(\lambda)}{b_1(\lambda)} (\boldsymbol{I} - \hat{\boldsymbol{z}}\hat{\boldsymbol{z}}) + \frac{a_2(\lambda)}{b_2(\lambda)}\hat{\boldsymbol{z}}\hat{\boldsymbol{z}} \right]\boldsymbol{\cdot}\boldsymbol{U} + {\rm \pi}\frac{a_3(\lambda)}{b_3(\lambda)}(\kern 1.5pt\hat{\boldsymbol{y}}\hat{\boldsymbol{x}} - \hat{\boldsymbol{x}}\hat{\boldsymbol{y}})\boldsymbol{\cdot}\boldsymbol{\varOmega}, \end{gather}$$

(4.12b)$$\begin{gather}\boldsymbol{T}^{(2)} = {\rm \pi}\frac{a_3(\lambda)}{b_3(\lambda)}(\hat{\boldsymbol{x}}\hat{\boldsymbol{y}} - \hat{\boldsymbol{y}}\hat{\boldsymbol{x}})\boldsymbol{\cdot}\boldsymbol{U} - {\rm \pi}\left[\frac{a_4(\lambda)}{b_4(\lambda)}(\boldsymbol{I} - \hat{\boldsymbol{z}}\hat{\boldsymbol{z}}) + \frac{a_5(\lambda)}{b_5(\lambda)}\hat{\boldsymbol{z}}\hat{\boldsymbol{z}}\right]\boldsymbol{\cdot} \boldsymbol{\varOmega}, \end{gather}$$

where $a_1 (\lambda ),\ldots, a_5(\lambda ), b_1(\lambda ),\ldots, b_5(\lambda )$ are polynomials specified in Appendix E. Here, the maximum degree of spherical harmonics defining the shape perturbation function, $f$, in (4.11b) is $l_{max}=3$. This gives $l_{maxsum} \equiv l_{max}+1=4$, as per § 3.3.1.

Coupling between translational and rotational motion is evident in (4.12) and appears at $O(\epsilon ^2)$, as anticipated in § 1. Moreover, the derived solution in (4.12) possesses the well-known, and required, symmetry relation in the coupling of force/rotation and torque/translation (Brenner Reference Brenner1964b). Namely, the tensorial coefficient of $\boldsymbol{\varOmega }$ in (4.12a) is the transpose of the coefficient of $\boldsymbol{U}$ in (4.12b).

4.2.2. Orientation-averaged force and torque

The orientation-averaged force and torque (§ 3.4) are derived from (3.17), (4.12a) and (4.12b):

(4.13a)$$\begin{gather} \bar{\boldsymbol{F}} ={-}6{\rm \pi} \boldsymbol{U}\left[1+\lambda+\frac{\lambda^2}{9} + \epsilon^2 \frac{a_6(\lambda)}{b_6(\lambda)}\right] + O(\epsilon^3), \end{gather}$$

(4.13b)$$\begin{gather}\bar{\boldsymbol{T}} ={-}8{\rm \pi} \boldsymbol{\varOmega} \left[\frac{1+\lambda+\dfrac{\lambda^2}{3}}{1+\lambda} + \epsilon^2 \frac{a_7(\lambda)}{b_7(\lambda)}\right] + O(\epsilon^3), \end{gather}$$

where the polynomials $a_6(\lambda )$, $b_6(\lambda )$, $a_7(\lambda )$ and $b_7(\lambda )$ are given in Appendix E. As for the prolate spheroid studied in § 4.1, non-sphericity enhances the magnitude of the orientation-averaged force and torque experienced by the particle.

4.2.3. Steady and high-frequency limits

Due to the complexity of (4.12), we provide its results in the steady and (high-frequency) inviscid limits, i.e. $|\lambda | \rightarrow 0$ and $|\lambda | \rightarrow \infty$, respectively. These may be useful in practice, e.g. in applications involving Brownian motion. The force and torque in the steady limit are

(4.14a)\begin{align} \boldsymbol{F}_{|\lambda|\rightarrow 0} &={-}6{\rm \pi} \left[1 + \epsilon\frac{1}{10}\left(\boldsymbol{I} - 3\hat{\boldsymbol{z}}\hat{\boldsymbol{z}}\right) + \epsilon^2 \left( \frac{501}{700} \left(\boldsymbol{I} - \hat{\boldsymbol{z}}\hat{\boldsymbol{z}}\right) + \frac{24}{35}\hat{\boldsymbol{z}}\hat{\boldsymbol{z}} \right)\right]\boldsymbol{\cdot} \boldsymbol{U} \nonumber\\ & \quad +\epsilon^2\frac{162{\rm \pi} }{5\sqrt{35}}\left(\hat{\boldsymbol{y}}\hat{\boldsymbol{x}} - \hat{\boldsymbol{x}}\hat{\boldsymbol{y}}\right)\boldsymbol{\cdot} \boldsymbol{\varOmega} + O(\epsilon^3), \end{align}

(4.14b)\begin{align} \boldsymbol{T}_{|\lambda|\rightarrow 0} &={-} 8{\rm \pi} \left[1 + \epsilon \frac{3}{10}\left(\boldsymbol{I} - 3\hat{\boldsymbol{z}}\hat{\boldsymbol{z}}\right) + \epsilon^2 \left(\frac{541}{140}(\boldsymbol{I} - \hat{\boldsymbol{z}}\hat{\boldsymbol{z}}) + \frac{227}{175}\hat{\boldsymbol{z}}\hat{\boldsymbol{z}}\right) \right]\boldsymbol{\cdot}\boldsymbol{\varOmega} \nonumber\\ &\quad + \epsilon^2\frac{162{\rm \pi} }{5\sqrt{35}}(\hat{\boldsymbol{x}}\hat{\boldsymbol{y}} - \hat{\boldsymbol{y}}\hat{\boldsymbol{x}})\boldsymbol{\cdot} \boldsymbol{U} + O(\epsilon^3), \end{align}

while the corresponding inviscid results, i.e. $|\lambda | \rightarrow \infty$, are

(4.15a)\begin{align} \boldsymbol{F}_{|\lambda|\rightarrow \infty} &={-}\frac{2{\rm \pi} }{3}\lambda^2\left[1 + \epsilon\frac{9}{10}(\boldsymbol{I} - 3\hat{\boldsymbol{z}}\hat{\boldsymbol{z}}) + \epsilon^2 \left( \frac{111}{140} (\boldsymbol{I} - \hat{\boldsymbol{z}}\hat{\boldsymbol{z}}) + \frac{681}{175}\hat{\boldsymbol{z}}\hat{\boldsymbol{z}} \right)\right]\boldsymbol{\cdot} \boldsymbol{U} \nonumber\\ & \quad +\epsilon^2\lambda^2\frac{6{\rm \pi} }{\sqrt{35}}(\kern 1.5pt\hat{\boldsymbol{y}}\hat{\boldsymbol{x}} - \hat{\boldsymbol{x}}\hat{\boldsymbol{y}})\boldsymbol{\cdot} \boldsymbol{\varOmega} + O(\epsilon^3), \end{align}

(4.15b)\begin{align} \boldsymbol{T}_{|\lambda|\rightarrow \infty} & ={-} 2{\rm \pi} \epsilon^2 \lambda^2 (\boldsymbol{I} - \hat{\boldsymbol{z}}\hat{\boldsymbol{z}})\boldsymbol{\cdot}\boldsymbol{\varOmega} + \epsilon^2\lambda^2\frac{6{\rm \pi} }{\sqrt{35}}(\hat{\boldsymbol{x}}\hat{\boldsymbol{y}} - \hat{\boldsymbol{y}}\hat{\boldsymbol{x}})\boldsymbol{\cdot} \boldsymbol{U} + O(\epsilon^3). \end{align}

The resultant orientation-averaged forces and torques in the steady and inviscid limits, derived from (4.13), are

(4.16a)$$\begin{gather} \bar{\boldsymbol{F}}_{|\lambda|\rightarrow 0} ={-}6{\rm \pi} \boldsymbol{U} \left[1 + \frac{247}{350} \epsilon^2 + O(\epsilon^3) \right], \end{gather}$$

(4.16b)$$\begin{gather}\bar{\boldsymbol{T}}_{|\lambda|\rightarrow 0} ={-}8{\rm \pi} \boldsymbol{\varOmega} \left[1 + \frac{1053}{350} \epsilon^2 + O(\epsilon^3) \right], \end{gather}$$

(4.16c)$$\begin{gather}\bar{\boldsymbol{F}}_{|\lambda|\rightarrow\infty} ={-}\frac{2{\rm \pi} }{3}\lambda^2 \boldsymbol{U} \left[1 + \frac{639}{350} \epsilon^2 + O(\epsilon^3) \right], \end{gather}$$

(4.16d)$$\begin{gather}\bar{\boldsymbol{T}}_{|\lambda|\rightarrow\infty} ={-}\frac{4{\rm \pi} }{3}\epsilon^2\lambda^2 \boldsymbol{\varOmega} + O(\epsilon^3). \end{gather}$$

We observe that the $O(\epsilon ^2)$ term provides the leading-order contribution to the torque in the inviscid limit, as eluded to in § 1. Again, particle non-sphericity enhances the magnitudes of the force and torque.

It remains to be established if the above-mentioned enhancements due to particle non-sphericity, in the magnitude of the orientation-averaged force and torque experienced by the spheroid and ‘pear-shaped’ particles, hold true for all particle shapes. Note that the surface area of a nearly spherical particle always exceeds that of its equivalent-volume sphere, which is consistent with this enhancement.

4.3. ‘Harmonic virion’

Virions are often studied via acoustic techniques (see e.g. Robach et al. Reference Robach, Michels, Cerf, Braunwald and Tripier-Darcy1983; Babincová, Sourivong & Babinec Reference Babincová, Sourivong and Babinec2000; Yang et al. Reference Yang, Lin, Liu, Lu, Hung, Huang, Tsai, Kao, Chen and Sun2015), corresponding to the small-amplitude, unsteady Stokes limit in the present study. The following model provides insight into the hydrodynamics of virions when a spherical geometry cannot be assumed, such as SARS-CoV-2 virions, which exhibit ‘spiky’ protrusions; see figure 3(a). This model may also prove useful in studying the hydrodynamics of SARS-CoV-2 virions in aerosols, which is one of their methods of transmission (Van Doremalen et al. Reference Van Doremalen2020). Rather than consider the true geometry of the virion, its shape is modelled using a single spherical harmonic to facilitate analysis and interpretation. This model is also applicable to other virions and is termed a ‘harmonic virion’; see figure 3(b–d). In this section, we study how the size and number of ‘spiky’ protrusions of the particles affect their hydrodynamics.

Figure 3.

(a) Illustration of a SARS-CoV-2 virion (Eckert & Higgins Reference Eckert and Higgins2020). (b–d) ‘Harmonic virions’ which provide a simple hydrodynamic model for a SARS-CoV-2 virion. The nearly spherical particles are constructed with linear combinations of (b) $(l,m)=(10,\pm 5)$, (c) $(l,m)=(20,\pm 10)$ and (d) $(l,m)=(40,\pm 20)$ for the first-order shape perturbation function $f$ in (4.18a).

Three different ‘harmonic virions’ are studied, representing a progressively rougher and ‘spikier’ virion. The shape perturbation functions for these particles are

(4.17)\begin{equation} {R_n(\theta,\phi)} = 1 + \epsilon f_n(\theta,\phi) + \epsilon^2 g_n(\theta,\phi), \end{equation}

with $n \in \{5,\, 10,\, 20\}$, and

(4.18a)$$\begin{gather} f_n(\theta,\phi) = {Y}_{2n, -n} (\theta,\phi) + ({-}1)^{n} {Y}_{2n, n}(\theta,\phi), \end{gather}$$

(4.18b)$$\begin{gather}g_n(\theta,\phi) ={-}\frac{1}{\sqrt{{\rm \pi} }}{Y}_{0,0}(\theta,\phi). \end{gather}$$

We calculate the orientation-averaged forces and torques experienced by these particles, which may be practically relevant due to the randomising effects of Brownian motion. The resulting analytical expressions are complicated and in this section we focus on the (numericised) asymptotic limits of small and large $|\lambda |$, reported to three significant figures; exact formulae for all $\lambda$ are contained in the Mathematica notebook (supplementary material). The resulting orientation-averaged force is

(4.19)\begin{equation} \bar{\boldsymbol{F}}_n = \bar{\boldsymbol{F}}^{(0)} + \epsilon^2 \bar{\boldsymbol{F}}_n^{(2)} + O(\epsilon^3), \end{equation}

with the subscript $n$ specified in (4.18), and

(4.20a)$$\begin{gather} \bar{\boldsymbol{F}}^{(0)} ={-} \boldsymbol{U}\left\{\begin{array}{@{}ll} 18.8 + 18.8\lambda, & |\lambda| \ll 1 \\ 2.09\lambda^2 + 18.8\lambda, & |\lambda| \gg 1 , \end{array}\right. \end{gather}$$

(4.20b)$$\begin{gather}\bar{\boldsymbol{F}}_{5}^{(2)} ={-} \boldsymbol{U}\left\{\begin{array}{@{}ll} 37.9 + 75.9\lambda, & |\lambda| \ll 1 \\ 12.2\lambda^2 + 174\lambda, & |\lambda| \gg 1 , \end{array}\right. \end{gather}$$

(4.20c)$$\begin{gather}\bar{\boldsymbol{F}}_{10}^{(2)} ={-} \boldsymbol{U}\left\{\begin{array}{@{}ll} 82.7 + 165\lambda, & |\lambda| \ll 1 \\ 27.1\lambda^2 + 654\lambda, & |\lambda| \gg 1 , \end{array}\right. \end{gather}$$

(4.20d)$$\begin{gather}\bar{\boldsymbol{F}}_{20}^{(2)} ={-}\boldsymbol{U}\left\{\begin{array}{@{}ll} 173 + 345\lambda, & |\lambda| \ll 1 \\ 57.1\lambda^2 + 2510\lambda, & |\lambda| \gg 1 . \end{array}\right. \end{gather}$$

The corresponding formula for the orientation-averaged torque is

(4.21)\begin{equation} \bar{\boldsymbol{T}}_n = \bar{\boldsymbol{T}}^{(0)} + \epsilon^2 \bar{\boldsymbol{T}}_n^{(2)} + O(\epsilon^3), \end{equation}

with

(4.22a)$$\begin{gather} \bar{\boldsymbol{T}}^{(0)} ={-} \boldsymbol{\varOmega} \left\{\begin{array}{@{}ll} 25.1, & |\lambda| \ll 1 \\ 8.38\lambda, & |\lambda| \gg 1 , \end{array}\right. \end{gather}$$

(4.22b)$$\begin{gather}\quad \bar{\boldsymbol{T}}_{5}^{(2)} ={-} \boldsymbol{\varOmega} \left\{\begin{array}{@{}ll} 162, & |\lambda| \ll 1 \\ 6.67\lambda^2 + 96.0\lambda, & |\lambda| \gg 1 , \end{array}\right. \end{gather}$$

(4.22c)$$\begin{gather}\bar{\boldsymbol{T}}_{10}^{(2)} ={-} \boldsymbol{\varOmega} \left\{\begin{array}{@{}ll} 342, & |\lambda| \ll 1 \\ 13.3\lambda^2 + 323\lambda, & |\lambda| \gg 1 , \end{array}\right. \end{gather}$$

(4.22d)$$\begin{gather}\bar{\boldsymbol{T}}_{20}^{(2)} ={-} \boldsymbol{\varOmega} \left\{\begin{array}{@{}ll} 702, & |\lambda| \ll 1 \\ 26.7\lambda^2 + 1180\lambda, & |\lambda| \gg 1 . \end{array}\right. \end{gather}$$

Figure 4 gives results for the orientation-averaged forces and torques, as a function of dimensionless frequency, $\lambda$; exact formulae from the Mathematica notebook (supplementary material) are used. The most pronounced difference between these results and that for a perfect sphere is observed for the orientation-averaged torque when $|\lambda | \gg 1$. The torque experienced by a perfect sphere scales as $\lambda$ for $|\lambda | \gg 1$ due to the vanishingly small influence of the no-slip boundary condition; see (4.22a). In contrast, the harmonic virion's orientation-averaged torque varies as $\lambda ^2$ due to its protrusions producing non-zero added mass in this inviscid limit; see (4.22b)–(4.22d). It is also evident that the orientation-averaged force and torque magnitudes increase as the degree and order of the spherical harmonics increase, i.e. the particles become spikier. This is not unexpected due to the increase in particle surface area.

Figure 4.

Orientation-averaged force and torque for the three ‘harmonic virions’ in figure 3. Magnitudes of the second-order orientation-averaged (a) force and (c) torque monotonically increase with $|\lambda |$. Also shown is the ratio of the magnitude of the second-order orientation-averaged (b) force and (d) torque to their leading-order counterparts (i.e. results for a perfect sphere).

The ratio of the second-order to leading-order solutions in figure 4(b,d) reveals intriguing behaviour. This ratio is not constant, as may be naively anticipated, establishing that no perfect sphere can mimic the frequency-dependent – or equivalently, time-dependent – orientation-averaged force and torque experienced by a harmonic virion. Interestingly, this ratio for the force is maximised for intermediate $|\lambda |$, the value of which increases with increasing $(l, m)$. This maximum occurs when the viscous penetration depth, over which vorticity diffuses, matches the spacing of the protrusions. In this case, the force becomes sensitively dependent on the precise shape and spacing of the protrusions. No such maxima occur for the torque due to the different dependencies on $\lambda$ in the high-frequency limit discussed above.

4.4. Comparison with the Brownian motion of anisotropic particles studied by Kraft et al. (Reference Kraft, Wittkowski, Ten Hagen, Edmond, Pine and Löwen2013)

Kraft et al. (Reference Kraft, Wittkowski, Ten Hagen, Edmond, Pine and Löwen2013) studied the motion of anisotropic particles constructed by fusing together spheres of different radii. They measured the friction tensors of the particles by monitoring their Brownian trajectories. For a particle consisting of three fused spheres of different radii (see figure 5), they observed coupling between the translational and rotational motion of the Brownian particles in both experiments and numerical simulations performed in the steady Stokes limit.

Figure 5.

(a) Geometry and numerical results for the resistance of an anisotropic Brownian particle studied by Kraft et al. (Reference Kraft, Wittkowski, Ten Hagen, Edmond, Pine and Löwen2013). The Brownian particle is constructed by fusing together three spheres of radii 2.1, 1.7 and 1.3 $\mathrm{\mu }$m, respectively. The resistance matrix, $\boldsymbol{M}$, relates the force and torque via $(\boldsymbol{F}, \boldsymbol{T}) = \boldsymbol{M}* (\boldsymbol{U}, \boldsymbol{\varOmega })$. Note that the resistance matrix is given in dimensional form where the units of the first, second, third and fourth 3$\times$3 submatrices are $\mathrm{\mu }$m, $\mathrm{\mu }$m$^2$, $\mathrm{\mu }$m$^2$ and $\mathrm{\mu }$m$^3$, respectively. (b) Present theory. Nearly spherical reconstruction of the anisotropic Brownian particle using spherical harmonics with $l \le 5$. The radius of the equivalent-volume sphere is calculated from the combined volumes of the three fused spheres and results in $R_{eq} = 2.5$ $\mathrm{\mu }$m. The resistance matrix is constructed from the theory in § 3 in the limit $\lambda \rightarrow 0$ and then cast into dimensional form using $R_{eq}$ and the dimensions above. This matrix is to be compared with the literature data of Kraft et al. (Reference Kraft, Wittkowski, Ten Hagen, Edmond, Pine and Löwen2013) reported in (a).

We reconstruct this particle using spherical harmonics up to order $l \le 5$, resulting in the shape given in figure 5(b). Using the force and torque formulae in § 3, and taking the limit as $\lambda \rightarrow 0$, we calculate the force and torque reported in figure 5(b). Critically, the only information used from Kraft et al. (Reference Kraft, Wittkowski, Ten Hagen, Edmond, Pine and Löwen2013) is the radii of the three spheres and the origin of the coordinate system as approximated from figure 5(a); no other fitting parameters are used. We observe excellent agreement between our nearly spherical theory and the numerical results of Kraft et al. (Reference Kraft, Wittkowski, Ten Hagen, Edmond, Pine and Löwen2013). Agreement between the components of the coupling tensor are particularly relevant because they arise exclusively from the $O(\epsilon ^2)$ terms reported in this study.