Unsteady motion of nearly spherical particles in viscous fluids: a second-order asymptotic theory

Abstract

The motion of small non-spherical particles is often studied using the unsteady Stokes equations. Zhang & Stone (J. Fluid Mech., vol. 367, 1998, pp. 329–358) reported an asymptotic treatment for nearly spherical particles, to first order in particle non-sphericity, i.e. $O(\epsilon )$, where $\epsilon$ quantifies the shape deviation from a sphere. Importantly, key physical phenomena are absent at $O(\epsilon )$, including (1) coupling between the torque experienced by the particle and its linear translation, (2) coupling between the force the particle experiences and its rotation and (3) the effect of non-sphericity on the orientation averages of these forces and torques. We present an explicit asymptotic theory to second order in particle non-sphericity, i.e. $O(\epsilon ^2)$, for the force and torque acting on a particle in a general unsteady Stokes flow. The derived analytical formulae apply to particles of arbitrary shape, providing the leading-order asymptotic theory for the three above-mentioned phenomena. The theory is demonstrated for several example nearly spherical particles including a spheroid, a ‘pear-shaped’ particle and a simple model for a SARS-CoV-2 virion. This includes formulae for force and torque as a function of particle orientation and their corresponding orientation averages. Our study reveals that the orientation-averaged forces and torques experienced by a nearly spherical particle cannot be generally represented by a perfect sphere. The reported formulae are validated using finite-amplitude three-dimensional direct numerical simulations of the Navier–Stokes equations. A Mathematica notebook is also provided, facilitating implementation of the theory for particle shapes of the user's choosing.

1. Introduction

The time-dependent motion of small particles immersed in viscous fluids has been studied extensively, providing the foundations for a multitude of applications and technologies. One prominent application is in microrheology where the characterisation of particle dynamics forms the basis for measurement (Squires & Mason 2010). Examples include the optical trapping of small spherical particles (Atakhorrami et al. 2005, 2008) through to monitoring of the rotational diffusion of highly non-spherical microdisks (Cheng & Mason 2003). Particle motion in such microrheological measurements can be induced either passively (e.g. via Brownian fluctuations) or through active forcing (e.g. using an electric/magnetic field).

Another example where unsteady particle motion can be used to advantage is in the autonomous propulsion of non-spherical particles in acoustic fields. It is known that when small particles, with either shape or density asymmetries, are trapped at the pressure node of an acoustic standing wave, they can undergo autonomous (self-induced) propulsion (e.g. see Wang et al. 2012; Rao et al. 2015; Ahmed et al. 2016). The primary mechanism proposed for this experimentally observed propulsion is (nonlinear) acoustic streaming, driven by (linear) small-amplitude oscillatory motion of the particle itself (Nadal & Lauga 2014; Collis, Chakraborty & Sader 2017; Nadal & Michelin 2020). In this application, coupling of translational and rotational motion of the particle drives directed streaming, producing a jet which propels the particle (Collis et al. 2017; Nadal & Michelin 2020; Collis, Chakraborty & Sader 2022). Interestingly, Lippera et al. (2019) proved that no propulsion occurs to first order in the radial shape perturbation of nearly spherical particles. This established that higher-order effects drive propulsion of these particles.

Such applications of unsteady particle motion are not limited to rigid particles. For example, Pelton et al. (2013) studied the small-amplitude (oscillatory) elastic vibration of bipyramidal gold nanoparticles that are performing ultra-high-frequency (gigahertz) elastic vibrations in fluid. The short time scale of the generated unsteady flow (in the picosecond regime) naturally interrogates molecular relaxation processes in simple liquids (Slie, Donfor & Litovitz 1966; O'Sullivan et al. 2019). This allows for direct interrogation of the viscoelastic response of these liquids (Chakraborty et al. 2018), which constitutes an application of nanorheology (Canale et al. 2019). An understanding of the hydrodynamic interaction of small particles is essential for interpreting measurements in all of the applications mentioned above.

Perfect solid spheres are a rare occurrence in nature, yet are used ubiquitously to model real-world phenomena including some of the applications listed above. A canonical example in fluid mechanics is the drag experienced by a perfect solid sphere moving steadily in a viscous fluid at low Reynolds number, i.e. Stokes’ law (Stokes 1851). The resulting formula, and its variants as a function of geometry, have found broad-ranging applications that include the sedimentation of dilute particulate suspensions (Richardson & Zaki 1954) through to the characterisation of optical laser traps (Perkins 2009).

Viscous flows generated by the motion of small particles occur precisely in this low-Reynolds-number regime. Consequently, they are often treated under the framework of the unsteady Stokes equations (that naturally encompass the steady case), which neglect the fluid's convective inertia while retaining its local inertia. This simplified approach is justified in a range of applications, including the small-amplitude oscillatory motion of micro- and nanoparticles. Neglect of the fluid's convective inertia greatly simplifies calculation of the flow field, leading to results that are linear in the particle velocity while generally retaining the dominant physical processes (Batchelor 2000). A well-known exception is ‘Stokes’ paradox’ for the steady, two-dimensional flow generated by a cylinder in an unbounded fluid, where convective inertial effects balance viscous effects in the particle's far field. Importantly, such exceptions do not occur in unsteady flows, flows generated by three-dimensional bodies or in bounded fluid domains.

Because particles are often modelled as perfect spheres, knowledge of the effects of (any degree of) non-sphericity is important for practical implementation. A key contribution to this understanding was reported by Zhang & Stone (1998) who developed an analytical theory that considers the radial degree of non-sphericity to be asymptotically small, quantified by the small non-dimensional parameter $\epsilon$. That previous work examined the effect of non-sphericity on the force and torque experienced by a nearly spherical particle undergoing arbitrary translational and rotational motion.

The theory of Zhang & Stone (1998) is derived to first order in the non-spherical shape parameter $\epsilon$, i.e. it is correct to $O(\epsilon )$. While accounting for much of the fluid physics of nearly spherical particles, this first-order approach does not capture several phenomena exhibited by non-spherical particles. For example, non-spherical particles undergoing pure steady translation experience finite torque. The first-order theory of Zhang & Stone (1998) predicts that this torque (about the particle's geometric centre) is identically zero to $O(\epsilon )$, regardless of the particle's shape and orientation. A prediction of zero torque does not coincide with observations of non-spherical particles. For illustrative purposes, consider the extreme case of a particle consisting of a slender rod attached to a sphere at one of its ends, moving at finite angle of attack to its principal geometric axis. The particle naturally experiences a non-zero torque about its geometric centre. This torque–translation coupling phenomenon is not restricted to highly non-spherical particles (as per this example). We therefore conclude that any such coupling between translational motion and (rotation-inducing) torque for nearly spherical particles must occur at higher order in $\epsilon$, i.e. $O(\epsilon ^2)$. Another illustration of this phenomenon is exhibited by the flagellum of bacteria. Many bacteria swim by rotating a helical flagellum (Berg & Anderson 1973). This conversion of rotation into linear translation is only possible due to rotation–force coupling. A nearly spherical analogue of such a particle will exhibit force–rotation and torque–translation coupling at $O(\epsilon ^2)$.

Another practical situation arises when examining the orientation-averaged force experienced by non-spherical particles in flow. For example, Cheng & Mason (2003) studied the rotational diffusion of highly non-spherical solid microdisks to determine the shear moduli of a surrounding viscoelastic fluid. This required the orientation average of the hydrodynamic force/torque experienced by the disks (Hubbard & Douglas 1993). The above-mentioned $O(\epsilon )$ theory predicts that the orientation average of the force experienced by a nearly spherical particle in its direction of motion is identical to that of a perfect sphere. This demonstrates that the existing first-order theory for nearly spherical particles does not capture the effect of particle non-sphericity in such measurements. A third phenomena is the non-zero torque experienced by a non-spherical particle undergoing oscillatory rotation in the (high-frequency) inviscid flow limit, for which the first-order theory predicts zero torque.

A higher-order theory is required to capture these key phenomena. This forms the aim of the present work, which develops theory for nearly spherical particles to second order in the non-spherical shape parameter $\epsilon$, i.e. to $O(\epsilon ^2)$. The theory is derived for nearly spherical particles of arbitrary shape, by performing a regular perturbation expansion in small $\epsilon$. The first-order flow field is required for the second-order solution, which is obtained using a general solution to the unsteady Stokes equations in vector spherical harmonics. An explicit analytical form for the second-order solution is derived, requiring the shape of the particle surface only. The utility of this general theory is demonstrated through its application to three example nearly spherical particles: (1) a prolate spheroid, (2) a ‘pear-shaped’ particle that exhibits force–rotation and torque–translation coupling and (3) a ‘harmonic virion’. The latter particle is a simple model for a SARS-CoV-2 virion, which is used to study how the number of ‘spiky’ surface protrusions affect the virion's hydrodynamics. In all cases, formulae are presented for generally unsteady, non-axisymmetric particle motion. A Mathematica notebook is provided in the supplementary material available at https://doi.org/10.1017/jfm.2024.1075 (with user notes in Appendix A), facilitating implementation of the derived theory for particle shapes of the user's choosing. An additional comparison is made with an anisotropic Brownian particle studied by Kraft et al. (2013), which demonstrates good agreement of the developed theory with existing results in the case of translation–rotation coupling.

We begin in § 2 by deriving the general theoretical framework for an arbitrary nearly spherical particle immersed in an unbounded and quiescent viscous fluid, to second order in $\epsilon$. This gives the force and torque experienced by the particle in terms of its motion and resulting fluid stress tensors; the latter are specified by the fluid velocity at the particle surface, via the Lorentz reciprocal theorem. Results of the implementation of this framework are given in § 3 where general and explicit analytical expressions for the force and torque, solely in terms of the particle motion, are reported. The above-mentioned example particles are then studied in § 4 using the derived general expressions from § 3. This includes analytical formulae for the force and torque experienced by these example particles, their orientation-averaged expressions and their limits of zero- and high-frequency motion. To demonstrate the validity of the derived formulae, their predictions are compared with three-dimensional finite-amplitude direct numerical simulations (DNS) of the (nonlinear) Navier–Stokes equations in § 5. Details of the theoretical and numerical calculations are relegated to appendices, which also contain information on the Mathematica notebook (supplementary material; information in Appendix A).

2. Theoretical framework

We consider a nearly spherical particle that executes small-amplitude, oscillatory rigid-body motion at a single angular frequency, $\omega$, and amplitude, $a$, in a quiescent and unbounded viscous fluid. The following non-dimensionalisation is performed: all distance variables are scaled by the radius of an equivalent-volume sphere, $R_{eq}$, time by $\omega ^{-1}$, velocity by $\omega a$, angular velocity by $\omega a/R_{eq}$ and pressure by $\omega \mu a / R_{eq}$, where $\mu$ is the fluid's shear viscosity; this specifies the force and torque scales of $\omega \mu a R_{eq}$ and $\omega \mu a R_{eq}^2$, respectively. Henceforth, all variables refer to their non-dimensional counterparts. The non-dimensional oscillation amplitude is $\delta \equiv a/R_{eq}$, and the governing Navier–Stokes equations is

(2.1)\begin{equation} \beta \left( \frac{\partial \boldsymbol{u}}{\partial t} + \delta\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u} \right) ={-}\boldsymbol{\nabla} p + \nabla^2\boldsymbol{u}, \quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \end{equation}

where $\beta = \rho R^{2}_{eq}\omega /\mu$ is the oscillatory Reynolds number of the flow and $\rho$ is the fluid density. The fluid boundary conditions are

(2.2a)$$\begin{gather} \boldsymbol{u} = {\rm Re}[\boldsymbol{U}\,\mathrm{e}^{-\mathrm{i}t}] + {\rm Re}[\boldsymbol{\varOmega}\,\mathrm{e}^{-\mathrm{i}t}]\times\boldsymbol{r} \quad \text{on}\ \boldsymbol{r}\in S_p, \end{gather}$$

(2.2b)$$\begin{gather}\boldsymbol{u} \rightarrow \boldsymbol{0} \quad \text{as}\ |\boldsymbol{r}|\rightarrow\infty, \end{gather}$$

where $S_p$ is the (time-varying) location of the particle surface, the constant vectors $\boldsymbol{U}$ and $\boldsymbol{\varOmega }$ specify the rigid-body (dimensionless) translational and angular velocities of the particle, respectively, and ${\rm Re}$ gives the real part of the expression. Figure 1 gives graphical illustrations of three nearly spherical particles that are studied in § 4 (with exaggerated shape perturbation, for clarity).

Figure 1. Examples of three nearly spherical particles generated using the radial coordinate specified by (2.5). The degree of (infinitesimal) non-sphericity is exaggerated for clarity. (a) Prolate spheroid with a major-axis length of $1+\epsilon$, with radial shape perturbation functions $f$ and $g$ specified in (4.3). (b) ‘Pear-shaped’ particle, with $f$ and $g$ specified in (4.11b). (c) ‘Harmonic virion’ where $f$ is specified by a linear combination of $(l,m) = (20,\pm 10)$ spherical harmonics, see (4.18a) with $n = 10$, and $g$ is specified in (4.18b).

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 (1998). 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 1964a). 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 (1998) 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. Force and torque in terms of the particle motion

We use the general theoretical framework in § 2 to derive analytical expressions for the force and torque on an arbitrary nearly spherical particle. The torque is evaluated with respect to the particle's geometric centre (about which the angular velocity is also specified).

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 1851). 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 (1998). 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 1998).

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. Application to example nearly spherical particles

In this section, we illustrate the utility of the force and torque formulae in § 3 for three particle geometries: (a) a prolate spheroid, (b) a ‘pear-shaped’ particle and (c) a ‘harmonic virion’; see figure 1. These examples serve to demonstrate the emergence of non-zero orientation-averaged forces and torques, and the coupling between translational/rotational particle motion and torque/force, the expressions for which may find use in practice. Analytical formulae for a pear-shaped particle in § 4.2 are compared with independent DNS of the Navier–Stokes equations in § 5. Additionally, a comparison is made with a Brownian motion study by Kraft et al. (2013) to examine the validity of the presented solution for the case of translation–rotation coupling.

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.

4.1.2. Second-order solution

The two components of the second-order solution in $\epsilon$ (defined in § 3.3), denoted by subscripts 1 and 2, are now determined. These are evaluated using the Mathematica notebook (supplementary material and Appendix A), which can be used for particle geometries of the user's choosing. Consequently, the results in this section serve as benchmarks for this future usage.

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 (1975) 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 (1975). 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 (1986).

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 (1975) 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 1964b). 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. 1983; Babincová, Sourivong & Babinec 2000; Yang et al. 2015), 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. 2020). 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 2020). (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. (2013)

Kraft et al. (2013) 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. (2013). 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. (2013) 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. (2013) 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. (2013). 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.

5. Direct numerical simulations at finite amplitude

In this section, we validate the (infinitesimally) small-amplitude theory reported in § 3 by comparison with finite-amplitude, time-dependent DNS of the (nonlinear) incompressible Navier–Stokes equations, (2.1) and (2.2). The DNS are performed using COMSOL Multiphysics. Details are reported in Appendix F.

The ‘pear-shaped’ particle described in § 4.2 is used for this validation, together with a normalised frequency of $\beta = 1$; other frequencies lead to similar conclusions (data not shown). The linear and angular oscillatory velocities are set to $\boldsymbol{U} = \delta \hat {\boldsymbol{z}}$ and $\boldsymbol{\varOmega } = \delta \hat {\boldsymbol{y}}$, respectively. The values of $\delta$ and $\epsilon$ are varied to examine the combined effects of non-sphericity and oscillation amplitude, respectively. For each set of these values, time series of the force and torque experienced by the particle are computed from which their complex Fourier series are determined. This uses a fundamental angular frequency of $\omega = 1$, as dictated by the above non-dimensionalisation. All DNS are converged to at least 98 % with respect to the magnitude of the largest non-zero $\textrm{e}^{-\textrm{i}t}$ Fourier coefficient of the force/torque; this coincides with the particle's oscillatory motion. Results for this fundamental Fourier coefficient of the DNS provide the finite-amplitude counterpart to the small-amplitude theory which is a solution of the (linear) unsteady Stokes equation. This enables direct comparison between the DNS and the small-amplitude theory. Symmetries in the ‘pear-shaped’ particle and oscillatory velocities produce non-zero results for (i) the force in the $x$ direction, denoted $F_x$, (ii) force in the $z$ direction, denoted $F_z$, and (iii) torque about the $y$ axis, denoted $T_y$. Results for these components from the DNS and small-amplitude theory are reported in figure 6. Details of the DNS and their convergence, including validation on a perfect sphere, are relegated to Appendix F.

Figure 6. Force and torque components of the ‘pear-shaped’ particle, with $\boldsymbol{U} = \delta \hat {\boldsymbol{z}}$ and $\boldsymbol{\varOmega } = \delta \hat {\boldsymbol{y}}$. Nearly spherical theory (solid lines); DNS (circles). Results for magnitudes of $F_x$, $F_z$ and $T_y$, for particle motion at $\beta = 1$. (a) Force/torque components at $\beta =1$ (fundamental frequency), for oscillation amplitude, $\delta = 0.1$, as non-spherical parameter $\epsilon$ is varied. (b) Force/torque components at $\beta =1$ (fundamental frequency), for non-spherical parameter, $\epsilon = 0.1$, as oscillation amplitude $\delta$ is varied. (c) Force/torque components at $\beta =0$ (steady streaming), for non-spherical parameter, $\epsilon = 0.1$, as oscillation amplitude $\delta$ is varied.

We first assess the validity of the second-order solution in the developed nearly spherical theory. For this purpose, the dimensionless oscillation amplitude is set to $\delta = 0.1$ and the non-spherical parameter $\epsilon$ is varied; note that $\epsilon =0$ corresponds to a perfect sphere. It is evident from figure 6(a) that the DNS and small-amplitude theory agree for smaller values of $\epsilon$, i.e. $\epsilon \lesssim 0.1$. This provides independent validation of the nearly spherical theory. Specifically, the DNS force in the $x$ direction, $F_x$, coincides precisely with the nearly spherical theory for smaller values of $\epsilon$. A quadratic variation in $\epsilon$ is predicted by the nearly spherical theory and observed in the DNS. A deviation emerges at larger values of $\epsilon$, which is not unexpected given the asymptotic nature of the nearly spherical theory. Similar agreement is observed for $F_z$ and $T_y$ with the second-order solution correctly capturing the DNS trends as the non-spherical parameter $\epsilon$ is increased. Results for the phase of these force and torque components exhibit similar agreement (data not shown). Critically, the DNS shows that the force in the $x$ direction, $F_x$, is coupled to the particle's angular velocity around the $y$ axis. This behaviour is captured correctly by the second-order solution; the first-order theory of Zhang & Stone (1998) predicts $F_x=0 + O(\epsilon ^2)$.

Next, we investigate the validity of the (leading-order) small-$\delta$ expansion employed in the nearly spherical theory, i.e. its foundational use of the unsteady Stokes equations. As noted at the end of § 2.2, this theory is expected to be valid provided the dimensionless oscillation amplitude is much smaller than the non-spherical parameter, i.e. $\delta \ll \epsilon$. For this purpose, we now fix the non-spherical parameter to be $\epsilon = 0.1$ and vary the oscillation amplitude, $\delta$. Results for the magnitude of the force and torque are reported in figure 6(b); similar results are found for the phase (data not shown). Linear dependence on the oscillation amplitude, $\delta$, for $\delta < \epsilon$ is evident in all results, as required. A deviation from linearity emerges for $\delta > \epsilon$ which is again not unexpected. This is due to the significant effect of nonlinear convective inertia in this regime (explored below). Note that the force component in the $x$ direction, $F_x$, is zero for a perfect sphere and depends quadratically on the non-spherical parameter $\epsilon$. The predictions of the nearly spherical theory precisely coincide with the DNS for this purely non-spherical force component. This comparison validates use of the unsteady Stokes equations as the foundation of the nearly spherical theory, and demonstrates its predicted regime of validity, $\delta \ll \epsilon$, given in § 2.2.

Increasing the oscillation amplitude, $\delta$, is expected to produce non-zero steady components of the force and torque. This is due to the effect of nonlinear convective inertia, which is absent from the nearly spherical theory. The DNS data for these steady components are reported in figure 6(c) and exhibit the expected nonlinear variation in oscillation amplitude, $\delta$. These steady components are generally much smaller than the fundamental frequency components reported in figure 6(b). However, we observe that when the steady component of $F_x$ in figure 6(b) is similar in magnitude to the fundamental frequency component in figure 6(c), i.e. $\delta \gtrsim \epsilon = 0.1$, the latter component no longer tracks linearly with the oscillation amplitude, $\delta$. This confirms that the nearly spherical theory – based on the (linear) unsteady Stokes equations – is valid provided the effect of nonlinear convective inertia is negligible. This occurs when $\delta \ll \epsilon$, as predicted in § 2.2. Analytical calculation of this streaming flow is beyond the scope of this study.

The present theory requires the oscillation amplitude, $\delta$, to be smaller than all other length scales, including the non-spherical shape parameter, $\epsilon$. This means that (strictly) the theory is not valid when $\epsilon \ll \delta$. However, the left-hand panel of figure 6(a) appears to show agreement in this regime. This is because the force in the $x$ direction is constrained to (i) zero as $\epsilon \rightarrow 0$ (the result for a sphere) and (ii) the derived non-spherical $O(\epsilon ^2)$ result for $\epsilon \gg \delta$. The intermediate regime of $0< \delta < \epsilon$ smoothly connects these two results, with any deviation from the theory not resolvable with the accuracy of our numerical simulations.

6. Conclusions

We have examined the force and torque experienced by a nearly spherical particle of arbitrary shape moving in a quiescent and unbounded viscous fluid. This builds on the previous work of Zhang & Stone (1998), by providing an asymptotic solution in the small non-sphericity parameter correct to $O(\epsilon ^2)$. In so doing, the derived general formulae enable analytical calculations of important fluid physical phenomena of arbitrary nearly spherical particles. These include hydrodynamic coupling between the particle's translational and rotational motion, and the effect of particle non-sphericity on the orientation-averaged force and torque that it experiences. General unsteady linear motion was considered according to the unsteady Stokes equations. Strictly, the derived general formulae apply when the oscillation amplitude is far smaller than the height of the particle's non-sphericity, i.e. $\delta \ll \epsilon$.

To illustrate the utility of these general formulae, we used them to explore the hydrodynamic forces experienced by three particles: (1) a spheroidal particle that exhibits fore–aft symmetry, (2) a ‘pear-shaped’ particle that does not possess this high degree of symmetry and (3) a ‘harmonic virion’ which models the SARS-CoV-2 virion. Closed-form expressions for the force and torque experienced by these particles were derived correct to $O(\epsilon ^2)$, together with their orientation-averaged formulae. This showed that a perfect sphere cannot be generally used to model the time-dependent, orientation-averaged forces experienced by a nearly spherical particle. Additionally, it was observed that non-sphericity enhances the magnitude of the orientation-averaged force and torque experienced by these three particles, a result whose applicability to general particles is yet to be established. It is envisaged that the general formulae in § 3 will find utility in the analytical solution of flow problems that have previously been inaccessible. This may be useful, for example, in improving the efficiency of computational simulations of particle suspensions that may rely on single particle results.

It is anticipated that the present theory will find use in studying the Brownian motion of nearly spherical particles. Because non-sphericity enters into the orientation-averaged forces at $O(\epsilon ^2)$, non-sphericity is also expected to affect orientationally averaged diffusivities at this order. Unsteady effects may also couple with the effects of non-sphericity to modify the velocity autocorrelation functions relative to those of a perfect sphere.

The derived nearly spherical theory was validated against independent finite-amplitude time-domain simulations of the (nonlinear) Navier–Stokes equations. Agreement with the DNS was observed when $\delta \lesssim \epsilon$, which aligns with the regime of validity of the derived theory; see above.

A Mathematica notebook is included in the supplementary material, which applies the general formulae in § 3 to nearly spherical particles with an arbitrary shape perturbation function. Usage details are discussed in Appendix A.