2. Statement of the problem

Consider a rigid body immersed in a Newtonian fluid with viscosity $\mu$ at vanishing Reynolds number. Let $V_b \subset \mathbb {R}^3$ be the domain representing the space occupied by the body and $V_f$ the space occupied by the fluid domain. The surface bounding the body is $S_b$, whereas the surface bounding the fluid is $S_b \cup S_w \cup S_\infty$, where $S_w$ is the surface bounding externally the fluid and considered in the proximity of the body, and $S_\infty$ the boundary at infinity, i.e. any surface considered infinitely far from the body. See the schematic representation of the system geometry in figure 1.

Figure 1. Schematic representation of geometry of the system.

The body is immersed in an ambient flow $(\boldsymbol{u}(\boldsymbol{x}), \boldsymbol{\rm \pi} (\boldsymbol{x}))$, which is defined as any flow, regular at the surface of the body $S_b$, satisfying the Stokes equations (Kim & Karrila Reference Kim and Karrila2005). Considering no-slip BCs on $S_w$, the Stokes equations for the ambient flow read

(2.1)\begin{equation} \left. \begin{gathered} - \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\rm \pi} (\boldsymbol{x}) = \mu \Delta \boldsymbol{u} (\boldsymbol{x})- \boldsymbol{\nabla} p(\boldsymbol{x}) = 0,\\ \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} (\boldsymbol{x}) = 0, \quad \boldsymbol{x} \in V_f \cup V_b,\\ \boldsymbol{u}(\boldsymbol{x})= 0, \quad \boldsymbol{x} \in S_w, \end{gathered} \right\} \end{equation}

where $\boldsymbol{u}(\boldsymbol{x})$ represents the velocity field, $p(\boldsymbol{x})$ the pressure field and

(2.2)\begin{equation} \boldsymbol{\rm \pi} (\boldsymbol{x})= p(\boldsymbol{x}) I- \mu (\boldsymbol{\nabla} \boldsymbol{u}(\boldsymbol{x})+\boldsymbol{\nabla} \boldsymbol{u}(\boldsymbol{x})^t) \end{equation}

is the stress tensor. In (2.2) $I$ represents the identity matrix and the superscript ‘$t$’ denotes the transposition operation for a matrix.

Assuming linear homogeneous BCs at the surface of the body $S_b$, expressed by a generic linear operator $\mathcal {L}[\, ]$ acting on the velocity at the surface of the body (this implies that $\mathcal {L}[\boldsymbol{v}(\boldsymbol{x}) ]$ at the point ${\boldsymbol x} \in S_b$ may depend not only on the velocity $\boldsymbol{v}(\boldsymbol{x})$ but also on its derivatives at that point), and no-slip BCs at the surface of the confinement $S_w$, the total (or disturbed) flow ($\boldsymbol{v}(\boldsymbol{x})$, $\boldsymbol{\sigma }(\boldsymbol{x})$) is

(2.3)\begin{equation} \left. \begin{gathered} - \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\sigma} (\boldsymbol{x}) = \mu \Delta \boldsymbol{v} (\boldsymbol{x})- \boldsymbol{\nabla} s(\boldsymbol{x}) = 0,\\ \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v} (\boldsymbol{x}) = 0, \quad \boldsymbol{x} \in V_f,\\ \mathcal{L}[ \boldsymbol{v} (\boldsymbol{x})]=0,\quad \boldsymbol{x} \in S_b,\\ \boldsymbol{v}(\boldsymbol{x})=0, \quad \boldsymbol{x} \in S_w, \end{gathered} \right\} \end{equation}

with the obvious meaning for $\boldsymbol{v} (\boldsymbol{x})$, $s(\boldsymbol{x})$, and $\boldsymbol{\sigma } (\boldsymbol{x})$, representing velocity, pressure and stress tensor fields of the total flow, respectively. Henceforth, we require that the BCs expressed by the linear operator $\mathcal {L}[\, ]$, satisfy the following condition at the surface of the body: for any couple of flows $(\boldsymbol{v}'(\boldsymbol{x}), \boldsymbol{\sigma }'(\boldsymbol{x}))$ and $(\boldsymbol{v}''(\boldsymbol{x}), \boldsymbol{\sigma }''(\boldsymbol{x}))$ solution of (2.3), the following identity holds:

(2.4)\begin{equation} \int_{S_b} (\boldsymbol{v}'(\boldsymbol{x}) \boldsymbol{\cdot} \boldsymbol{\sigma}''(\boldsymbol{x}) - \boldsymbol{v}(\boldsymbol{x}) \boldsymbol{\cdot} \boldsymbol{\sigma}'(\boldsymbol{x})) \boldsymbol{\cdot} \boldsymbol{n}(\boldsymbol{x}) \,{\rm d}S = 0, \end{equation}

where $\boldsymbol{n}(\boldsymbol{x})$ is the unitary normal vector outwardly oriented with respect to the body as shown in figure 1. This condition has been introduced and thoroughly discussed in Procopio & Giona (Reference Procopio and Giona2024) in connection to the concept of the Hinch–Kim dualism, and it is referred to as the condition of BC reciprocity. BCs satisfying (2.4), are therefore called reciprocal BCs. As it will become clear in the remainder, the assumption of BC reciprocity, coupled with the linearity of the BCs expressed by the operator $\mathcal {L}[\,]$, are the necessary prerequisites in the development of the present theory. Linear reciprocal BCs are, for example, no-slip $\boldsymbol{v}(\boldsymbol{x})=0$, complete slip $\boldsymbol{\sigma }(\boldsymbol{x})=0$ and Navier-slip BCs

(2.5)\begin{equation} \boldsymbol{v}(\boldsymbol{x})+\dfrac{\lambda}{\mu} \boldsymbol{n}(\boldsymbol{x})\boldsymbol{\cdot} \boldsymbol{\sigma}(\boldsymbol{x})\boldsymbol{\cdot} \boldsymbol{t}(\boldsymbol{x})=0 \end{equation}

with $\lambda$ being the slip length and $\boldsymbol{t}(\boldsymbol{x})$ an unitary base vectors tangents to the surface of the body. On the other hand, BCs representing the interaction of the Stokes fluid with a non-Newtonian fluid, for example, are not reciprocal because the reciprocity theorem does not hold in the body domain.

In the next paragraph, the solution of the problem (2.3) is expressed in terms of the hydrodynamic solutions of two simpler problems related separately to the confinement of the fluid and to the body: (i) the Green function of the Stokes equations in the domain of the confinement $V_f \cup V_b$ and (ii) the geometrical moments of the body in the unbounded fluid. See the schematic representation in figure 2. For this reason, it is useful to define these solutions and discuss briefly their formal properties, introducing and clarifying in this way the notation that we use throughout this article.

Figure 2. Schematic representation of the decomposition of the main problem (a body in a confined fluid with an ambient flow) in problems separately related to the geometry of the confinement and of the particle: (i) multipoles of the confinement (at the top) centred at the position point (e.g. the barycentric coordinate) in the volume the body; (ii) flows around the body in $n$th-order ambient flows centred at the position of the body provided by the Faxén operators (at the bottom). Green arrows in (a) represent concentrated forces whereas blue arrows in (b) represent velocity fields. The black dot is the position point of the body.

2.1. Green function of the confinement

As discussed in Procopio & Giona (Reference Procopio and Giona2023), the Green function in the confined domain $V_f \cup V_b$ is a bitensorial field, hence a field depending on two points (called field and source points) with entries at both points expressed, in principle, in different coordinate systems.

The Green function $G_{a \beta } (\boldsymbol{x},\boldsymbol{\xi })$ of the confined flow is the solution of the equations

(2.6)\begin{equation} \left. \begin{gathered} -\boldsymbol{\nabla}_b \varSigma_{a b \beta}(\boldsymbol{x},\boldsymbol{\xi})= \Delta G_{a \beta} (\boldsymbol{x},\boldsymbol{\xi})-\boldsymbol{\nabla}_a P_{\beta} (\boldsymbol{x},\boldsymbol{\xi}) ={-}8{\rm \pi} \delta_{a \beta}\delta(\boldsymbol{x}-\boldsymbol{\xi}),\\ \boldsymbol{\nabla}_a G_{a \beta}(\boldsymbol{x},\boldsymbol{\xi}) = 0, \quad \boldsymbol{x}, \boldsymbol{\xi} \in V_f \cup V_b,\\ G_{a \beta} (\boldsymbol{x},\boldsymbol{\xi})= 0, \quad \boldsymbol{x} \in S_w \cup S_\infty, \end{gathered} \right\} \end{equation}

where $G_{a \beta } (\boldsymbol{x},\boldsymbol{\xi })$, $P_{\beta } (\boldsymbol{x},\boldsymbol{\xi })$, $\varSigma _{a b \beta }(\boldsymbol{x},\boldsymbol{\xi })$ are the associated velocity, pressure and stress tensor field. In (2.6) and in the remainder, the notation of the bitensor calculus is applied: (i) Latin letters $a,b,\ldots = 1,2,3$ are used for indices referred to the entries of the tensorial entities at the field point $\boldsymbol{x}$, with Greek letters $\alpha, \beta,\ldots = 1,2,3$ for indices referred to the entries of the tensorial entities at the source points (i.e. the poles of the singularity) $\boldsymbol{\xi }$; (ii) Einstein's summation convention is adopted; and (iii) $\boldsymbol {\nabla }_a$, with the Latin index, is the gradient with respect to the field point $\boldsymbol{x}$, while $\boldsymbol {\nabla }_\beta$, with the Greek index, is the gradient with respect to the source point $\boldsymbol{\xi }$.

It is useful to define also the regular part of the Green function $(W_{a \beta }(\boldsymbol{x},\boldsymbol{\xi }), Q_\beta (\boldsymbol{x},\boldsymbol{\xi }))$ as the bitensorial fields solving the problem

(2.7)\begin{equation} \left. \begin{gathered} -\boldsymbol{\nabla}_b T_{a b \beta}(\boldsymbol{x},\boldsymbol{\xi})= \Delta W_{a \beta} (\boldsymbol{x},\boldsymbol{\xi})-\boldsymbol{\nabla}_a Q_{\beta} (\boldsymbol{x},\boldsymbol{\xi}) = 0,\\ \boldsymbol{\nabla}_a W_{a \beta}(\boldsymbol{x},\boldsymbol{\xi}) = 0, \quad \boldsymbol{x}, \boldsymbol{\xi} \in V_f \cup V_b,\\ W_{a \beta} (\boldsymbol{x},\boldsymbol{\xi})={-}S_{a \beta}(\boldsymbol{x}-\boldsymbol{\xi}), \quad \boldsymbol{x} \in S_w \cup S_\infty, \end{gathered} \right\} \end{equation}

where $S_{a \beta }(\boldsymbol{x}-\boldsymbol{\xi })$ is the Stokeslet (Pozrikidis Reference Pozrikidis1992; Kim & Karrila Reference Kim and Karrila2005), i.e. the bitensorial velocity field of the unbounded Green function

(2.8)\begin{equation} S_{a \beta}(\boldsymbol{x}-\boldsymbol{\xi})=\dfrac{\delta_{a \beta}}{|(\boldsymbol{x}-\boldsymbol{\xi})|}+\dfrac{(\boldsymbol{x}-\boldsymbol{\xi})_a(\boldsymbol{x}-\boldsymbol{\xi})_\beta}{|(\boldsymbol{x}-\boldsymbol{\xi})|^3}. \end{equation}

Therefore, the bounded Green function can be written as the superposition of a regular field $W_{a \beta }(\boldsymbol{x},\boldsymbol{\xi })$ and a singular contribution given by the Stokeslet

(2.9)\begin{equation} G_{a \beta}(\boldsymbol{x},\boldsymbol{\xi})=S_{a \beta}(\boldsymbol{x}-\boldsymbol{\xi})+W_{a \beta} (\boldsymbol{x},\boldsymbol{\xi}). \end{equation}

By differentiating (2.9) at the pole, higher-order singularities in the domain $V_b \cup V_f$ are obtained. For example, the $n$th-order multipole, with $n=1,2,\ldots$, is obtained by

(2.10)\begin{equation} \boldsymbol{\nabla}_{\boldsymbol{\beta}_n} G_{a \beta}(\boldsymbol{x},\boldsymbol{\xi})= \boldsymbol{\nabla}_{\boldsymbol{\beta}_n} S_{a \beta}(\boldsymbol{x}-\boldsymbol{\xi})+ \boldsymbol{\nabla}_{\boldsymbol{\beta}_n} W_{a \beta} (\boldsymbol{x},\boldsymbol{\xi}), \end{equation}

where bold index $\boldsymbol{\beta }_n= \beta _1 \ldots \beta _n$ denotes a multi-index and $\boldsymbol {\nabla }_{\boldsymbol{\beta }_n}= \boldsymbol {\nabla }_{\beta _1} \ldots \boldsymbol {\nabla }_{\beta _n}$ is a compact notation for $n$th-order differentiation.

2.2. Moments on the body and Faxén operators

Let us briefly define the $n$th-order moments, the $(m,n)$th-order geometrical moments and the $n$th-order Faxén operators of a body, addressed in more detail in Procopio & Giona (Reference Procopio and Giona2024).

Consider the same body immersed in an ambient flow ($\boldsymbol{u}(\boldsymbol{x}),\boldsymbol{\rm \pi} (\boldsymbol{x})$) in the unbounded domain. The disturbance flow ($\boldsymbol{w}(\boldsymbol{x})$, $\boldsymbol{\tau }(\boldsymbol{x})$) generated by the body immersed in the ambient flow is solution of

(2.11)\begin{equation} \left. \begin{gathered} - \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\tau} (\boldsymbol{x}) = \mu \Delta \boldsymbol{w} (\boldsymbol{x})- \boldsymbol{\nabla} q(\boldsymbol{x}) = 0,\\ \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{w} (\boldsymbol{x}) = 0, \quad \boldsymbol{x} \in V_f,\\ \mathcal{L}[ \boldsymbol{w} (\boldsymbol{x})]={-}\mathcal{L}[ \boldsymbol{u}(\boldsymbol{x}) ], \quad \boldsymbol{x} \in S_b,\\ \boldsymbol{w}(\boldsymbol{x})=0, \quad \boldsymbol{x} \rightarrow \infty, \end{gathered} \right\} \end{equation}

where $\boldsymbol{w}(\boldsymbol{x})$, $q(\boldsymbol{x})$ and $\boldsymbol{\tau }(\boldsymbol{x})$ are the associated disturbance velocity, pressure and stress tensor fields accounting for the hydrodynamics at the surface $S_b$ of the rigid body due to the interaction with the ambient flow $\boldsymbol{u}(\boldsymbol{x})$.

To begin with, consider the entries $\psi _a(\boldsymbol{x})$ of any force field distribution, with compact support belonging to $V_b$, such that the Ladyzhenskaya (Reference Ladyzhenskaya2014) volume potential reads

(2.12)\begin{equation} \dfrac{1}{8 {\rm \pi}\mu} \int_{V_b} \psi_\beta(\boldsymbol{\xi}) S_{a \beta}(\boldsymbol{x}-\boldsymbol{\xi}) \,{\rm d}V(\boldsymbol{\xi})= {w}_a(\boldsymbol{x}) \quad \boldsymbol{x} \in S_b . \end{equation}

Next, the $n$th-order moments on the body immersed in a generic ambient flow $\boldsymbol{u}(\boldsymbol{x})$, generating a disturbance field at the surface of the body $\boldsymbol{w}(\boldsymbol{x})$, are defined as

(2.13)\begin{equation} M_{\beta \boldsymbol{\beta}_n}(\boldsymbol{\xi})= \int _{V_b} (\boldsymbol{x}-\boldsymbol{\xi})_{\boldsymbol{\beta}_n} \psi_\beta(\boldsymbol{x}) \,{\rm d}V(\boldsymbol{x}) \quad \boldsymbol{\xi} \in V_b , \end{equation}

where $(\boldsymbol{x} - \boldsymbol{\xi })_{\boldsymbol{\beta }_n}=(\boldsymbol{x} - \boldsymbol{\xi })_{\beta _1} \ldots (\boldsymbol{x} - \boldsymbol{\xi })_{\beta _n}$ and where $(\boldsymbol{x}-\boldsymbol{\xi })_{\beta }=g_{\beta a}(\boldsymbol{\xi },\boldsymbol{x})(\boldsymbol{x}-\boldsymbol{\xi })_{a}$ and $\psi _\beta (\boldsymbol{x})= g_{\beta a}(\boldsymbol{\xi },\boldsymbol{x}) \psi _a(\boldsymbol{x})$, $g_{\beta a}(\boldsymbol{\xi },\boldsymbol{x})$ being the transformation matrix between the coordinate systems at the pole and field point (or, more generically, the parallel propagator (Poisson et al. Reference Poisson, Pound and Vega2011), i.e. the bitensor transforming the entries of a vector at the point $\boldsymbol{x}$ into the entries at the point $\boldsymbol{\xi }$).

Consider an $n$th-order polynomial ambient flow $\boldsymbol{u}^{(n)}(\boldsymbol{x},\boldsymbol{\xi })$, centred at a point $\boldsymbol{\xi }' \in V_b$ within the domain of the body, with entries

(2.14)\begin{equation} u_a(\boldsymbol{x})=A_{a \boldsymbol{a}_n}(\boldsymbol{x}-\boldsymbol{\xi}')_{\boldsymbol{a}_n}, \quad \boldsymbol{\xi}' \in V_b \end{equation}

($A_{a \boldsymbol{a}_n}$ being an $(n+1)$-dimensional constant tensor) and its associated disturbance field ${w}^{(n)}_a(\boldsymbol{x},\boldsymbol{\xi }')$, obtained from (2.13) by a force field distribution $\boldsymbol{\psi }^{(n)}(\boldsymbol{\xi }, \boldsymbol{\xi }')$. According to the definition (2.13), the $m$th-order moments on the body immersed in the $n$th-order ambient flow can be expressed as

(2.15)\begin{equation} M^{(n)}_{\beta \boldsymbol{\beta}_m}(\boldsymbol{\xi}, \boldsymbol{\xi}')= \int (\boldsymbol{x}-\boldsymbol{\xi})_{\boldsymbol{\beta}_m} \psi^{(n)}_\beta(\boldsymbol{x}, \boldsymbol{\xi}') \,{\rm d}V(\boldsymbol{x}) \quad \boldsymbol{\xi}, \boldsymbol{\xi}' \in V_b. \end{equation}

By the linearity of the Stokes equations with respect to the constant tensor $A_{a \boldsymbol{a}_n}$, we can define the $(m,n)$th-order geometrical moments $m_{\beta \boldsymbol{\beta }_m \gamma ' \boldsymbol{\gamma }_n'}(\boldsymbol{\xi },\boldsymbol{\xi }')$ by the relation

(2.16)\begin{equation} M^{(n)}_{\beta \boldsymbol{\beta}_m}(\boldsymbol{\xi}, \boldsymbol{\xi}')= 8 {\rm \pi}\mu A_{\gamma \boldsymbol{\gamma}_n} m_{\beta \boldsymbol{\beta}_m \gamma' \boldsymbol{\gamma}_n'}(\boldsymbol{\xi},\boldsymbol{\xi}') , \end{equation}

where the multi-index $\gamma ' \boldsymbol{\gamma }_n'$ is referred to the entries of the field at the point $\boldsymbol{\xi }'$.

Based on the hierarchy of the geometrical moments, the operator

(2.17)\begin{equation} \mathcal{F}_{\beta \gamma' \boldsymbol{\gamma}_n'}= \sum_{m=0}^{\infty} \dfrac{ m_{\beta \boldsymbol{\beta}_m \gamma' \boldsymbol{\gamma}_n'} (\boldsymbol{\xi},\boldsymbol{\xi}') \boldsymbol{\nabla}_{\boldsymbol{\beta}_m} }{m!} \end{equation}

can be introduced. As shown in Procopio & Giona (Reference Procopio and Giona2024), if BC reciprocity holds, $\mathcal {F}_{\beta \gamma ' \boldsymbol{\gamma }_n'}$ represents the $n$th-order Faxén operator of the body. By the assumption that the operator $\mathcal {L}[\boldsymbol{v}(\boldsymbol{x})]$ in the total Stokes system (2.3) belongs to the class of linear homogeneous reciprocal BCs satisfying (2.4), the following relations for the body in the unbounded domain hold (Procopio & Giona Reference Procopio and Giona2023, Reference Procopio and Giona2024):

(2.18)\begin{equation} M_{\beta \boldsymbol{\beta}_n}(\boldsymbol{\xi})= 8 {\rm \pi}\mu \mathcal{F}_{\gamma' \beta \boldsymbol{\beta}_n} u_{\gamma'}(\boldsymbol{\xi}') \end{equation}

and

(2.19)\begin{equation} w_a(\boldsymbol{x})= \sum_{n=0}^{\infty} \dfrac{ \boldsymbol{\nabla}_{\boldsymbol{\gamma}_n'} u_{\gamma'}(\boldsymbol{\xi}') }{n!}\mathcal{F}_{\beta \gamma' \boldsymbol{\gamma}_n'}{S}_{a \beta}(\boldsymbol{x}-\boldsymbol{\xi}). \end{equation}

Furthermore, owing to the property that $\mathcal {F}_{\beta \gamma ' \boldsymbol{\gamma }_n'}$ is a Faxén operator, the disturbance field can be expressed by

(2.20)\begin{align} w_a(\boldsymbol{x}) &= \sum_{m=0}^{\infty} \dfrac{ \mathcal{F}_{\gamma' \beta \boldsymbol{\beta}_m} u_{\gamma'}(\boldsymbol{\xi}') }{m!}\boldsymbol{\nabla}_{\boldsymbol{\beta}_m}{S}_{a \beta}(\boldsymbol{x}-\boldsymbol{\xi})\nonumber\\ &= \dfrac{1}{8{\rm \pi}\mu} \sum_{m=0}^{\infty} \dfrac{ M_{\beta \boldsymbol{\beta}_m}(\boldsymbol{\xi}) }{m!}\boldsymbol{\nabla}_{\boldsymbol{\beta}_m}{S}_{a \beta}(\boldsymbol{x}-\boldsymbol{\xi}) . \end{align}

Finally, it is useful in the remainder to remark that the force exerted by the fluid on the body is $F_\beta = -M_\beta (\boldsymbol{\xi })$, thus, by (2.18),

(2.21)\begin{equation} F_{\beta}={-}8 {\rm \pi}\mu \mathcal{F}_{\gamma \beta} u_{\gamma}(\boldsymbol{\xi}) \end{equation}

whereas the torque $T_\beta = \varepsilon _{\beta \gamma \gamma _1} M_{\gamma \gamma _1}(\boldsymbol{\xi })$, is given by

(2.22)\begin{equation} T_{\beta}= 8 {\rm \pi}\mu \mathcal{T}_{\gamma \beta }\, u_{\gamma}(\boldsymbol{\xi}), \end{equation}

where $\mathcal {T}_{\delta \beta }= \varepsilon _{\beta \gamma \gamma _1} \mathcal {F}_{\delta \gamma \gamma _1}$ and $\varepsilon _{\beta \gamma \gamma _1}$ the Ricci–Levi Civita symbol.

3. The flow due to a body in a confined fluid

3.1. The reflection method

Consider the problem defined by (2.3) providing the total flow in the system in the case of no-slip conditions both on the body surface and on the confinement walls, thus considering the identity matrix as operator $\mathcal {L}[\,]=I$. Owing to the linearity of the equations and of the BCs, we can apply the reflection method (see Happel & Brenner Reference Happel and Brenner1983) to express the solution $(v_a(\boldsymbol{x}),\sigma _{a b}(\boldsymbol{x}))$ as the superposition of a countable system of fields $(v_a^{[k]}(\boldsymbol{x}),\sigma _{a b}^{[k]}(\boldsymbol{x}))$, with $k=0,1,2,\ldots$,

(3.1)\begin{equation} \left. \begin{gathered} v_a(\boldsymbol{x}) = v_a^{[0]}(\boldsymbol{x}) + v_a^{[1]}(\boldsymbol{x})+ \cdots+ v_a^{[k]}(\boldsymbol{x}) + \cdots,\\ \sigma_{a b}(\boldsymbol{x})= \sigma_{ab}^{[0]}(\boldsymbol{x}) + \sigma_{ab}^{[1]}(\boldsymbol{x}) + \cdots+ \sigma_{ab}^{[k]}(\boldsymbol{x}) + \cdots, \end{gathered} \right\} \end{equation}

where

(3.2)\begin{equation} \sigma_{ab}^{[k]}(\boldsymbol{x})=s^{[k]}(\boldsymbol{x})\delta_{ab}-\mu (\boldsymbol{\nabla}_a v^{[k]}_b(\boldsymbol{x})+ \boldsymbol{\nabla}_b v^{[k]}_a(\boldsymbol{x})), \end{equation}

with $s^{[k]}(\boldsymbol{x})$ being the associated pressure, each of which is the solution of the Stokes equations equipped with the following system of BCs:

(3.3)\begin{equation} \left. \begin{gathered} v_a^{[2k+1]}(\boldsymbol{x}) ={-}v_a^{[2k]}(\boldsymbol{x}) , \quad \boldsymbol{x} \in S_b, \\ v_a^{[2k+2]}(\boldsymbol{x}) ={-}v_a^{[2k+1]}(\boldsymbol{x}), \quad \boldsymbol{x} \in S_w. \end{gathered} \right\} \end{equation}

For $k=0$,

(3.4)\begin{equation} v_a^{[0]}(\boldsymbol{x}) = u_a(\boldsymbol{x}), \quad \boldsymbol{x} \in V_b \cup V_f. \end{equation}

Hence, as can be observed from (3.3), for odd $k$ the condition involves the boundary of the body, for even $k$ the walls of the confinement.

3.2. The velocity fields ${v}^{[1]}$ and ${v}^{[2]}$

Let us start by expressing the first velocity fields $v_a^{[1]}(\boldsymbol{x})$ and $v_a^{[2]}(\boldsymbol{x})$ in terms of the Green function of the confinement and the Faxén operator of the body that are supposed to be given.

Comparing (3.3) with (2.11) it is easy to recognise that ${\boldsymbol{v}}^{[1]}(\boldsymbol{x})$ is the disturbance field of the ambient field $\boldsymbol{u}(\boldsymbol{x})$. Therefore, by using (2.19), it is possible to explicit the velocity field with $k=1$ as

(3.5)\begin{equation} v_a^{[1]}(\boldsymbol{x})= \sum_{n=0}^{\infty} \dfrac{ \boldsymbol{\nabla}_{\boldsymbol{\gamma}_n} {u}_\gamma(\boldsymbol{\xi}) }{n!}\mathcal{F}_{\beta \gamma \boldsymbol{\gamma}_n}{S}_{a \beta}(\boldsymbol{x}-\boldsymbol{\xi}). \end{equation}

Alternatively, from (2.20), the first velocity field can be expressed as

(3.6)\begin{equation} v_a^{[1]}(\boldsymbol{x}) = \dfrac{1}{8{\rm \pi}\mu} \sum_{m=0}^{\infty} \dfrac{ {M}_{\beta \boldsymbol{\beta}_m}(\boldsymbol{\xi}) }{m!}\boldsymbol{\nabla}_{\boldsymbol{\beta}_m}{S}_{a \beta}(\boldsymbol{x}-\boldsymbol{\xi}). \end{equation}

Since, by linearity, any $\boldsymbol{v}^{[k]}(\boldsymbol{x})$ is solution of the Stokes equations, equipped with the BCs (3.3), the flow with $k=2$ is the solution of the problem

(3.7)\begin{equation} \left. \begin{gathered} \mu \Delta v_a^{[2]}(\boldsymbol{x})-\boldsymbol{\nabla}_a s^{[2]}(\boldsymbol{x})=0,\\ \boldsymbol{\nabla}_a v_a^{[2]}(\boldsymbol{x})=0, \quad \boldsymbol{x} \in V_b \cup V_f ,\\ v_a^{[2]}(\boldsymbol{x})={-}v_a^{[1]}(\boldsymbol{x}), \quad \boldsymbol{x} \in S_w. \end{gathered} \right\} \end{equation}

By applying the operator

(3.8)\begin{equation} \sum_{n=0}^{\infty} \dfrac{ \boldsymbol{\nabla}_{\boldsymbol{\gamma}_n} {u}_\gamma(\boldsymbol{\xi}) }{n!}\mathcal{F}_{\beta \gamma \boldsymbol{\gamma}_n}, \end{equation}

at a source point $\boldsymbol{\xi } \in V_b$ of the regular part of the Green function defined by the (2.7), and comparing the resulting problem with (3.7), it is easy to conclude, by the uniqueness of the solution of Stokes equations, that

(3.9)\begin{equation} v_a^{[2]}(\boldsymbol{x})= \sum_{n=0}^{\infty} \dfrac{ \boldsymbol{\nabla}_{\boldsymbol{\gamma}_n} {u}_\gamma(\boldsymbol{\xi}) }{n!}\mathcal{F}_{\beta \gamma \boldsymbol{\gamma}_n}{W}_{a \beta}(\boldsymbol{x},\boldsymbol{\xi}) \end{equation}

or, alternatively, by applying the operator

(3.10)\begin{equation} \dfrac{1}{8{\rm \pi}\mu} \sum_{m=0}^{\infty} \dfrac{ {M}_{\beta \boldsymbol{\beta}_m}(\boldsymbol{\xi}) }{m!}\boldsymbol{\nabla}_{\boldsymbol{\beta}_m}, \end{equation}

we obtain the representation

(3.11)\begin{equation} v_a^{[2]}(\boldsymbol{x}) = \dfrac{1}{8{\rm \pi}\mu} \sum_{m=0}^{\infty} \dfrac{ {M}_{\beta \boldsymbol{\beta}_m}(\boldsymbol{\xi}) }{m!}\boldsymbol{\nabla}_{\boldsymbol{\beta}_m}{W}_{a \beta}(\boldsymbol{x},\boldsymbol{\xi}) \end{equation}

and, thus,

(3.12)\begin{align} v_a^{[1]}(\boldsymbol{x})+v_a^{[2]}(\boldsymbol{x}) &= \sum_{n=0}^{\infty} \dfrac{ \boldsymbol{\nabla}_{\boldsymbol{\gamma}_n} {u}_\gamma(\boldsymbol{\xi}) }{n!}\mathcal{F}_{\beta \gamma \boldsymbol{\gamma}_n}{G}_{a \beta}(\boldsymbol{x},\boldsymbol{\xi})\nonumber\\ &= \dfrac{1}{8{\rm \pi}\mu} \sum_{m=0}^{\infty} \dfrac{ {M}_{\beta \boldsymbol{\beta}_m}(\boldsymbol{\xi}) }{m!}\boldsymbol{\nabla}_{\boldsymbol{\beta}_m}{G}_{a \beta}(\boldsymbol{x},\boldsymbol{\xi}) . \end{align}

3.3. The velocity fields ${v}^{[3]}$ and ${ {v}}^{[4]}$

From the BCs (3.3), the velocity field for $k=3$ is the disturbance field of $v_a^{[2]}(\boldsymbol{x})$ and, therefore, by (2.19),

(3.13)\begin{equation} v_a^{[3]}(\boldsymbol{x})= \sum_{\ell=0}^{\infty} \dfrac{ \boldsymbol{\nabla}_{\boldsymbol{\delta}_{\ell}'} v^{[2]}_{\delta'}(\boldsymbol{\xi}') }{\ell!}\mathcal{F}_{{ \gamma \delta'} \boldsymbol{\delta}_{\ell}'}{S}_{a {\gamma}}(\boldsymbol{x}-\boldsymbol{\xi}) \end{equation}

and, equivalently to (3.9),

(3.14)\begin{equation} v_a^{[4]}(\boldsymbol{x})= \sum_{\ell=0}^{\infty} \dfrac{ \boldsymbol{\nabla}_{\boldsymbol{\delta}_{\ell}'} v^{[2]}_{\delta'}(\boldsymbol{\xi}') }{\ell!}\mathcal{F}_{{ \gamma \delta'} \boldsymbol{\delta}_{\ell}'}{W}_{a {\gamma}}(\boldsymbol{x}-\boldsymbol{\xi}). \end{equation}

Substituting (3.11) into (3.13) one obtains

(3.15)\begin{equation} v_a^{[3]}(\boldsymbol{x})= \dfrac{1}{8{\rm \pi}\mu} \sum_{m=0}^{\infty} \dfrac{ {M}_{\beta \boldsymbol{\beta}_m}(\boldsymbol{\xi}) }{m!} \sum_{\ell=0}^{\infty} \dfrac{ \boldsymbol{\nabla}_{\boldsymbol{\delta}_{\ell}'} \boldsymbol{\nabla}_{\boldsymbol{\beta}_m} {W}_{\delta' \beta}(\boldsymbol{\xi}',\boldsymbol{\xi}) }{\ell!}\mathcal{F}_{{\gamma} {\delta}' \boldsymbol{\delta}_{\ell}'}{S}_{a {\gamma}}(\boldsymbol{x}-\boldsymbol{\xi}). \end{equation}

By the equivalence between the two expressions (2.19) and (2.20),

(3.16)\begin{equation} \sum_{\ell=0}^{\infty} \dfrac{ \boldsymbol{\nabla}_{\boldsymbol{\delta}_{\ell}'} \boldsymbol{\nabla}_{\boldsymbol{\beta}_m} {W}_{\delta' \beta}(\boldsymbol{\xi}',\boldsymbol{\xi}) }{\ell!}\mathcal{F}_{{\gamma} {\delta}' \boldsymbol{\delta}_{\ell}'}{S}_{a {\gamma}}(\boldsymbol{x}-\boldsymbol{\xi}) = \sum_{n=0}^{\infty} \dfrac{ \mathcal{F}_{{\delta}' {\gamma} \boldsymbol{\gamma}_{n}} \boldsymbol{\nabla}_{\boldsymbol{\beta}_m} {W}_{\delta' \beta}(\boldsymbol{\xi}',\boldsymbol{\xi}) }{n!} \boldsymbol{\nabla}_{\boldsymbol{\gamma}_{n}}{S}_{a {\gamma}}(\boldsymbol{x}-\boldsymbol{\xi}). \end{equation}

It is useful to introduce the tensor $N_{\alpha \boldsymbol{\alpha }_m \beta \boldsymbol{\beta }_{n}}(\boldsymbol{\xi })$ defined as

(3.17)\begin{equation} N_{\beta \boldsymbol{\beta}_m \gamma \boldsymbol{\gamma}_{n}}(\boldsymbol{\xi})= \left.\mathcal{F}_{\delta' \gamma \boldsymbol{\gamma}_{n}}\boldsymbol{\nabla}_{\boldsymbol{\beta}_m} W_{\delta' \beta}(\boldsymbol{\xi}',\boldsymbol{\xi})\right|_{\boldsymbol{\xi}'=\boldsymbol{\xi}}, \end{equation}

which corresponds to the $n$th-order moment on the body immersed in an ambient field consisting in the regular part of the $m$th derivative of the Green function. The tensor defined in (3.17) is fundamental for the further development of this analysis because, as shown in the following, it completely represents the hydrodynamic interaction between the body and the confinement.

Using the identity (3.16) and the definition (3.17), (3.15) can be expressed as

(3.18)\begin{equation} v_a^{[3]}(\boldsymbol{x})= \dfrac{1}{8{\rm \pi}\mu} \sum_{m=0}^{\infty} \dfrac{ { M}_{\beta \boldsymbol{\beta}_m}(\boldsymbol{\xi}) }{m!} \sum_{n=0}^{\infty} \dfrac{ N_{\beta \boldsymbol{\beta}_m \gamma \boldsymbol{\gamma}_{n}}(\boldsymbol{\xi}) }{n!}\boldsymbol{\nabla}_{\boldsymbol{\gamma}_{n}}{S}_{a {\gamma}}(\boldsymbol{x}-\boldsymbol{\xi}) \end{equation}

and, enforcing the same argument applied above to obtain (3.7)–(3.11), we have

(3.19)\begin{equation} v_a^{[4]}(\boldsymbol{x})= \dfrac{1}{8{\rm \pi}\mu} \sum_{m=0}^{\infty} \dfrac{ {M}_{\beta \boldsymbol{\beta}_m}(\boldsymbol{\xi}) }{m!} \sum_{n=0}^{\infty} \dfrac{ N_{\beta \boldsymbol{\beta}_m \gamma \boldsymbol{\gamma}_{n}}(\boldsymbol{\xi}) }{n!}\boldsymbol{\nabla}_{\boldsymbol{\gamma}_{n}}{W}_{a {\gamma}}(\boldsymbol{x},\boldsymbol{\xi}) \end{equation}

so that

(3.20)\begin{equation} v_a^{[3]}(\boldsymbol{x})+v_a^{[4]}(\boldsymbol{x})= \dfrac{1}{8{\rm \pi}\mu} \sum_{m=0}^{\infty} \dfrac{ {M}_{\beta \boldsymbol{\beta}_m}(\boldsymbol{\xi}) }{m!} \sum_{n=0}^{\infty} \dfrac{ N_{\beta \boldsymbol{\beta}_m \gamma \boldsymbol{\gamma}_{n}}(\boldsymbol{\xi}) }{n!}\boldsymbol{\nabla}_{\boldsymbol{\gamma}_{n}}{G}_{a {\gamma}}(\boldsymbol{x},\boldsymbol{\xi}) . \end{equation}

3.4. The velocity fields ${v}^{[5]}$ and ${v}^{[6]}$

The subsequent velocity fields can be determined following the same procedure used for $\boldsymbol{v}^{[3]}(\boldsymbol{x})$ and $\boldsymbol{v}^{[4]}(\boldsymbol{x})$. In fact, $\boldsymbol{v}^{[5]}(\boldsymbol{x})$ can be considered as the disturbance field of $\boldsymbol{v}^{[4]}(\boldsymbol{x})$ and, thus,

(3.21)\begin{equation} v_a^{[5]}(\boldsymbol{x})= \sum_{\ell=0}^{\infty} \dfrac{ \boldsymbol{\nabla}_{\boldsymbol{\gamma}_{\ell}'} v^{[4]}_{\gamma'}(\boldsymbol{\xi}') }{\ell!}\mathcal{F}_{{\delta} {\gamma'} \boldsymbol{\gamma}_{\ell}'}{S}_{a {\delta}}(\boldsymbol{x}-\boldsymbol{\xi}) \end{equation}

and, equivalently to (3.9),

(3.22)\begin{equation} v_a^{[6]}(\boldsymbol{x})= \sum_{\ell=0}^{\infty} \dfrac{ \boldsymbol{\nabla}_{\boldsymbol{\gamma}_{\ell}'} v^{[4]}_{\gamma'}(\boldsymbol{\xi}') }{\ell!}\mathcal{F}_{{\delta} {\gamma'} \boldsymbol{\gamma}_{\ell}'}{W}_{a {\delta}}(\boldsymbol{x}-\boldsymbol{\xi}). \end{equation}

Enforcing the same argument used previously in (3.15)–(3.20) for $v_a^{[3]}(\boldsymbol{x})$ and $v_a^{[4]}(\boldsymbol{x})$, we obtain

(3.23)\begin{equation} v_a^{[5]}(\boldsymbol{x})= \dfrac{1}{8{\rm \pi}\mu} \sum_{m=0}^{\infty} \dfrac{ {M}_{\beta \boldsymbol{\beta}_m}(\boldsymbol{\xi}) }{m!} \sum_{n=0}^{\infty} \dfrac{ N_{\beta \boldsymbol{\beta}_m \gamma \boldsymbol{\gamma}_{n}}(\boldsymbol{\xi}) }{n!} \sum_{\ell=0}^{\infty} \dfrac{ N_{\gamma \boldsymbol{\gamma}_{n} \delta \boldsymbol{\delta}_{\ell}}(\boldsymbol{\xi}) }{\ell!} \boldsymbol{\nabla}_{\boldsymbol{\delta}_{\ell}}{S}_{a \delta}(\boldsymbol{x}-\boldsymbol{\xi}), \end{equation}

(3.24)\begin{equation} v_a^{[6]}(\boldsymbol{x})= \dfrac{1}{8{\rm \pi}\mu} \sum_{m=0}^{\infty} \dfrac{ {M}_{\beta \boldsymbol{\beta}_m}(\boldsymbol{\xi}) }{m!} \sum_{n=0}^{\infty} \dfrac{ N_{\beta \boldsymbol{\beta}_m \gamma \boldsymbol{\gamma}_{n}}(\boldsymbol{\xi}) }{n!} \sum_{\ell=0}^{\infty} \dfrac{ N_{\gamma \boldsymbol{\gamma}_{n} \delta \boldsymbol{\delta}_{\ell}}(\boldsymbol{\xi}) }{\ell!} \boldsymbol{\nabla}_{\boldsymbol{\delta}_{\ell}}{W}_{a \delta}(\boldsymbol{x},\boldsymbol{\xi}), \end{equation}

so that

(3.25)\begin{equation} v_a^{[5]}(\boldsymbol{x})+v_a^{[6]}(\boldsymbol{x})= \dfrac{1}{8{\rm \pi}\mu} \sum_{m=0}^{\infty} \dfrac{ {M}_{\beta \boldsymbol{\beta}_m}(\boldsymbol{\xi}) }{m!} \sum_{n=0}^{\infty} \dfrac{ N_{\beta \boldsymbol{\beta}_m \gamma \boldsymbol{\gamma}_{n}}(\boldsymbol{\xi}) }{n!} \sum_{\ell=0}^{\infty} \dfrac{ N_{\gamma \boldsymbol{\gamma}_{n} \delta \boldsymbol{\delta}_{\ell}}(\boldsymbol{\xi}) }{\ell!} \boldsymbol{\nabla}_{\boldsymbol{\delta}_{\ell}}{G}_{a \delta}(\boldsymbol{x},\boldsymbol{\xi}). \end{equation}

3.5. The total velocity field

Iterating the same procedure for any $k$, it is possible to generalise the above results in the form

(3.26)\begin{equation} v_a^{[2k+1]}(\boldsymbol{x})\! =\! \dfrac{1}{8{\rm \pi}\mu}\! \sum_{m=0}^{\infty} \dfrac{ {M}_{\beta \boldsymbol{\beta}_m}(\boldsymbol{\xi}) }{m!}\! \sum_{m_1=0}^{\infty}\! \dfrac{ N_{\beta \boldsymbol{\beta}_m \gamma \boldsymbol{\gamma}_{m_1}}(\boldsymbol{\xi}) }{m_1!} \ldots\! \sum_{m_k=0}^{\infty}\! \dfrac{ N_{\delta \boldsymbol{\delta}_{m_{k-1}} \eta \boldsymbol{\eta}_{m_k}}(\boldsymbol{\xi}) }{m_k!} \boldsymbol{\nabla}_{\boldsymbol{\eta}_{m_k}}{S}_{a \eta}(\boldsymbol{x}-\boldsymbol{\xi}) \end{equation}

and

(3.27)\begin{equation} v_a^{[2k+2]}(\boldsymbol{x})\! =\! \dfrac{1}{8{\rm \pi}\mu}\! \sum_{m=0}^{\infty} \dfrac{ {M}_{\beta \boldsymbol{\beta}_m}(\boldsymbol{\xi}) }{m!}\! \sum_{m_1=0}^{\infty}\! \dfrac{ N_{\beta \boldsymbol{\beta}_m \gamma \boldsymbol{\gamma}_{m_1}}(\boldsymbol{\xi}) }{m_1!} \ldots\!\sum_{m_k=0}^{\infty}\! \dfrac{ N_{\delta \boldsymbol{\delta}_{m_{k-1}} \eta \boldsymbol{\eta}_{m_k}}(\boldsymbol{\xi}) }{m_k!} \boldsymbol{\nabla}_{\boldsymbol{\eta}_{m_k}}{W}_{a \eta}(\boldsymbol{x},\boldsymbol{\xi}) ,\end{equation}

so that

(3.28)\begin{align} v_a^{[2k+1]}(\boldsymbol{x})+v_a^{[2k+2]}(\boldsymbol{x})&= \dfrac{1}{8{\rm \pi}\mu} \sum_{m=0}^{\infty} \dfrac{ {M}_{\beta \boldsymbol{\beta}_m}(\boldsymbol{\xi}) }{m!} \sum_{m_1=0}^{\infty} \dfrac{ N_{\beta \boldsymbol{\beta}_m \gamma \boldsymbol{\gamma}_{m_1}}(\boldsymbol{\xi}) }{m_1!} \nonumber\\ &\quad \cdots\sum_{m_k=0}^{\infty} \dfrac{ N_{\delta \boldsymbol{\delta}_{m_{k-1}} \eta \boldsymbol{\eta}_{m_k}}(\boldsymbol{\xi}) }{m_k!} \boldsymbol{\nabla}_{\boldsymbol{\eta}_{m_k}}{G}_{a \eta}(\boldsymbol{x},\boldsymbol{\xi}). \end{align}

Summing all the fields according to (3.1), the total velocity field can be expressed as

(3.29)\begin{align} &v_a(\boldsymbol{x})-u_a(\boldsymbol{x})=\sum_{k=0}^\infty v_a^{[2k+1]}(\boldsymbol{x})+v_a^{[2k+2]}(\boldsymbol{x})\nonumber\\ &\quad =\dfrac{1}{8{\rm \pi}\mu} \sum_{m=0}^{\infty} \dfrac{ {M}_{\beta \boldsymbol{\beta}_m}(\boldsymbol{\xi}) }{m!} \sum_{k=0}^\infty \sum_{m_1=0}^{\infty} \dfrac{ N_{\beta \boldsymbol{\beta}_m \gamma \boldsymbol{\gamma}_{m_1}}(\boldsymbol{\xi}) }{m_1!} \ldots\sum_{m_k=0}^{\infty} \dfrac{ N_{\delta \boldsymbol{\delta}_{m_{k-1}} \eta \boldsymbol{\eta}_{m_k}}(\boldsymbol{\xi}) }{m_k!} \boldsymbol{\nabla}_{\boldsymbol{\eta}_{m_k}}{G}_{a \eta}(\boldsymbol{x},\boldsymbol{\xi}). \end{align}

3.6. Extension to Linear and BC-reciprocal BCs on the body

The extension to more general linear homogeneous reciprocal BCs is straightforward. To this purpose, the expansion (3.1) still applies, but the BCs (3.3) are substituted by the conditions

(3.30)\begin{equation} \left. \begin{gathered} \mathcal{L} [v_a^{[2k+1]}(\boldsymbol{x})] ={-}\mathcal{L}[v_a^{[2k]}(\boldsymbol{x})] , \quad \boldsymbol{x} \in S_b,\\ v_a^{[2k+2]}(\boldsymbol{x}) ={-}v_a^{[2k+1]}(\boldsymbol{x}), \quad \boldsymbol{x} \in S_w, \end{gathered} \right\} \end{equation}

for $k=1,2,\ldots$, keeping for $k=0$

(3.31)\begin{equation} v_a^{[0]}(\boldsymbol{x}) = u_a(\boldsymbol{x}), \quad \boldsymbol{x} \in V_b \cup V_f .\end{equation}

In fact, by applying the operator $\mathcal {L}[\, ]$ to the total field in (3.1) at the surface of the body, and using the linearity of $\mathcal {L}[\, ]$, we have

(3.32)\begin{equation} \left. \begin{aligned} & \mathcal{L}[ v_a(\boldsymbol{x})] = \mathcal{L}[ v_a^{[0]}(\boldsymbol{x}) + v_a^{[1]}(\boldsymbol{x}) + v_a^{[2]}(\boldsymbol{x}) + v_a^{[3]}(\boldsymbol{x}) + \cdots]=\\ & \mathcal{L}[v_a^{[0]}(\boldsymbol{x})] +\mathcal{L}[ v_a^{[1]}(\boldsymbol{x})] + \mathcal{L}[v_a^{[2]}(\boldsymbol{x})] +\mathcal{L}[ v_a^{[3]}(\boldsymbol{x})] + \cdots=0 , \quad \boldsymbol{x} \in S_b, \end{aligned} \right\} \end{equation}

where all the terms on the right-hand side cancel each other for (3.30). Therefore, the Stokes flow provided by (3.1), with BCs (3.30)–(3.31), is a solution of (2.3).

Since the procedure developed in the previous paragraph (3.5)–(3.29) is independent of the BCs at the surface of the body, with the only constraint of the BC reciprocity of $\mathcal {L}[\, ]$, we can conclude that (3.29) is still valid considering the Faxén operators associated with the BCs assumed on the body surface.

4. Matrix representation of the velocity field

In this section, a compact and useful matrix representation of the equations obtained in the § 3 is developed. To this aim, collect the entries of the system of moments ${M}_{\beta \boldsymbol{\beta }_m}(\boldsymbol{\xi })$ in an infinite-dimensional vector (Cooke Reference Cooke1950)

(4.1)\begin{equation} [{M}]= \left[ \begin{array}{c} \boldsymbol{M}_{(0)} \\ \boldsymbol{M}_{(1)} \\ \dfrac{\boldsymbol{M}_{(2)}}{2} \\ \vdots \\ \dfrac{\boldsymbol{M}_{(m)}}{m!} \\ \vdots \end{array} \right], \end{equation}

where $\boldsymbol{M}_{(m)}$ are $3^{m+1}$ dimensional vectors obtained by the vectorisation of the $(m+1)$th-order tensors ${M}_{\beta \boldsymbol{\beta }_m}(\boldsymbol{\xi })$ so that any entry $[M]_i$ corresponds to the entry ${M}_{\beta \boldsymbol{\beta }_m}(\boldsymbol{\xi })$ according the conversion $i \leftrightarrow \beta \boldsymbol{\beta }_m$

(4.2)\begin{equation} i=\sum_{h=0}^m \beta_{h} 3^{m-h}, \end{equation}

where $\beta _0 \equiv \beta$.

The conversion $i \leftrightarrow \beta \boldsymbol{\beta }_m$ for $m=0,1,2$, according to (4.2), is presented in table 1.

Table 1. Conversion according to (4.2) between the index $i$ of the entries of the vector $[{M}]_i$ and the multi-index $\beta \, \boldsymbol{\beta }_m$ of the entries of the $(m+1)$th-order tensors ${M}_{\beta \boldsymbol{\beta }_m}(\boldsymbol{\xi })$.

We use the notation $[{M}_{(n:m)}]$ to indicate the part of the array (4.1) collecting the entries of the tensors with orders going from $n$ to $m$ ($m>n$), i.e.

(4.3)\begin{equation} [{M}_{(n:m)}]= \left[ \begin{array}{c} \dfrac{\boldsymbol{M}_{(n)}}{n!} \\ \vdots \\ \dfrac{\boldsymbol{M}_{(m)}}{m!} \\ \end{array} \right]. \end{equation}

In the same way, the entries of $\boldsymbol {\nabla }_{\boldsymbol{\beta }_m} G_{a \beta }(\boldsymbol{x},\boldsymbol{\xi })$ can be collected in the $3^{m+1} \times 3$ matrices $\boldsymbol{G}_{(0)}, \boldsymbol{G}_{(1)},\ldots, \boldsymbol{G}_{(m)},\ldots$ (with column indices corresponding to the Latin field point indices) to build the $\infty \times 3$ matrix $[G]$ defined by

(4.4)\begin{equation} [G]= \left[ \begin{array}{ccccc} \boldsymbol{G}_{(0)} \\ \boldsymbol{G}_{(1)} \\ \vdots \\ \boldsymbol{G}_{(m)} \\ \vdots \end{array} \right]. \end{equation}

Adopting this representation, (3.12) can be compactly expressed as

(4.5)\begin{equation} \boldsymbol{v}^{[1]}(\boldsymbol{x})+\boldsymbol{v}^{[2]}(\boldsymbol{x})= \dfrac{[{M}]^t[G]}{8 {\rm \pi}\mu} ,\end{equation}

with $[{M}]^t$ being the transpose of $[{M}]$.

It is also possible to define the infinite matrix $[N]$ as

(4.6)\begin{equation} [N]= \left[ \begin{array}{cccccc} \boldsymbol{N}_{(0,0)} & \boldsymbol{N}_{(0,1)} & \cdots & \dfrac{\boldsymbol{N}_{(0,n)}}{n!} & \cdots\\ \boldsymbol{N}_{(1,0)} & \ddots & & \vdots\\ \vdots & & \ddots & \vdots\\ {\boldsymbol{N}_{(m,0)}} & \cdots & \cdots & \dfrac{\boldsymbol{N}_{(m,n)}}{n!} & \cdots\\ \vdots & & & \vdots & \ddots \end{array} \right], \end{equation}

where $\boldsymbol{N}_{(m,n)}$ are $3^{m+1} \times 3^{n+1}$ matrices obtained unfolding the $(m+n+2)$th-order tensors $N_{\beta \boldsymbol{\beta }_m \gamma \boldsymbol{\gamma }_n}(\boldsymbol{\xi} )$ so that the entries $[N]_{i,j}$ are obtained by converting both $i \leftrightarrow \beta \, \boldsymbol{\beta }_m$ and $j \leftrightarrow \gamma \, \boldsymbol{\gamma }_m$ according to (4.2).

Using this representation, (3.20) becomes

(4.7)\begin{equation} \boldsymbol{v}^{[3]}(\boldsymbol{x})+\boldsymbol{v}^{[4]}(\boldsymbol{x})= \dfrac{[{M}]^t[N][G]}{8{\rm \pi}\mu} ,\end{equation}

whereas (3.25) takes the form

(4.8)\begin{equation} \boldsymbol{v}^{[5]}(\boldsymbol{x})+\boldsymbol{v}^{[6]}(\boldsymbol{x})= \dfrac{[{M}]^t[N]^2[G]}{8{\rm \pi}\mu}, \end{equation}

where $[N]^2=[N][N]$. Defining the power of $[N]$ by induction as $[N]^3=[N]^2[N]$, $[N]^k=[N]^{k-1}[N]$ and $[N]^0=[I]$, $[I]$ being the infinite identity matrix, the total velocity field expressed by (3.29) can be compactly represented as

(4.9)\begin{equation} \boldsymbol{v}(\boldsymbol{x})-\boldsymbol{u}(\boldsymbol{x})= \dfrac{1}{8{\rm \pi}\mu}\sum_{k=0}^{\infty}[{M}]^t[N]^k[G]. \end{equation}

Let us consider the sum entering (4.9) truncated up to $k=K$ and multiply it by $([I]- [N])$. It is straightforward to show that

(4.10)\begin{equation} ([I] - [N])\sum_{k=0}^{K}[N]^k= [I]-[N]^{K} \end{equation}

as for the truncated geometric series defined over a scalar field. As shown in Appendix A, the series in (4.10) converges for characteristic distances $\ell _d$ of the body from the nearest walls larger enough than the characteristic length $\ell _b$ of the body itself, since there exists a constant $\varGamma = O(1) > 0$, depending on the geometry of the system, such that

(4.11)\begin{equation} \lim_{K\rightarrow \infty}[N]^K=0, \quad \mbox{for} \ \ell_d > \varGamma \ell_b. \end{equation}

As a consequence

(4.12)\begin{equation} ([I] - [N])\sum_{k=0}^{\infty}[N]^k= [I], \quad \ell_d > \varGamma \ell_b \end{equation}

and, thus,

(4.13)\begin{equation} \sum_{k=0}^{\infty}[N]^k= ([I] - [N])^{{-}1}, \quad \ell_d > \varGamma \ell_b ,\end{equation}

with $([I] - [N])^{-1}$ being the inverse matrix of $([I] - [N])$ (Cooke Reference Cooke1950).

Therefore, the velocity field attains the simple expression

(4.14)\begin{equation} \boldsymbol{v}(\boldsymbol{x})-\boldsymbol{u}(\boldsymbol{x})= \dfrac{ [{M}]^t ([I] - [N])^{{-}1} [G]}{8{\rm \pi}\mu} , \quad \ell_d > \varGamma\, \ell_b \end{equation}

or, alternatively,

(4.15)\begin{equation} \boldsymbol{v}(\boldsymbol{x})-\boldsymbol{u}(\boldsymbol{x})= \dfrac{[{M}]^t [X]}{8{\rm \pi}\mu} , \quad \ell_d > \varGamma \ell_b, \end{equation}

where $[X]$ is the solution of the infinite-matrix equation

(4.16)\begin{equation} ([I] - [N])[X]= [G]. \end{equation}

In the remainder, we consider exclusively the situation $\ell _d > \varGamma \ell _b$, for which (4.14) holds. By (4.14) it is possible to conclude that, for $\ell _d > \varGamma \ell _b$, the solution of a the Stokes flow past a body immersed in a confined fluid can be expressed in terms of the vector $[M]$, depending only on the hydrodynamic interaction of the body with unbounded (external) Stokes flows, the matrix $[G]$, depending only on the hydrodynamic interaction of the confinement with unbounded (internal) Stokes flows and the matrix $[N]$ ($[N]$-matrix, for short), representing the hydrodynamic interaction between the body and the confinement.

5. Force and torque on the particle

By linearity, the force and the torque acting on the particle due to the hydrodynamic interactions with the fluid, are given by the summation of all the forces and torques associated with the terms in (3.1), i.e.

(5.1)\begin{equation} \left. \begin{gathered} \boldsymbol{F}= \boldsymbol{F}^{[0]} + \boldsymbol{F}^{[1]} + \boldsymbol{F}^{[2]}+ \cdots,\\ \boldsymbol{T}= \boldsymbol{T}^{[0]} + \boldsymbol{T}^{[1]} + \boldsymbol{T}^{[2]}+ \cdots, \end{gathered} \right\} \end{equation}

where

(5.2)\begin{equation} \left. \begin{gathered} \boldsymbol{F}^{[k]}={-} \int_{S_p} \boldsymbol{\sigma}^{[k]} (\boldsymbol{x}) \boldsymbol{\cdot} \boldsymbol{n} \,{\rm d}S , \quad k=0,1,2,\ldots,\\ \boldsymbol{T}^{[k]}={-} \int_{S_p} (\boldsymbol{x}-\boldsymbol{\xi})\times \boldsymbol{\sigma}^{[k]} (\boldsymbol{x}) \boldsymbol{\cdot} \boldsymbol{n} \,{\rm d}S , \quad k=0,1,2,\ldots. \end{gathered} \right\} \end{equation}

Since $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol{\sigma }^{[k]}(\boldsymbol{x})=0$, and due to the symmetry of the stress tensors $\boldsymbol{\sigma }^{[k]}(\boldsymbol{x})$, the forces and torques associated to even values of $k$ (i.e. the forces due to regular fields on the boundary of the particle) vanish for the Gauss–Green theorem (Brenner Reference Brenner1962). The only terms contributing to the total force and torque are the terms corresponding to odd value of $k=1, 3, 5,\ldots$

(5.3)\begin{equation} \left. \begin{gathered} \boldsymbol{F}= \boldsymbol{F}^{[1]} + \boldsymbol{F}^{[3]}+ \boldsymbol{F}^{[5]} + \cdots,\\ \boldsymbol{T}= \boldsymbol{T}^{[1]} + \boldsymbol{T}^{[3]}+ \boldsymbol{T}^{[5]} + \cdots. \end{gathered} \right\} \end{equation}

The first contribution $\boldsymbol{F}^{[1]}$ in the sum (5.3) is the force experienced by the body immersed in the unbounded ambient flow $\boldsymbol{u}(\boldsymbol{x})$, therefore it can be obtained by applying the zeroth-order Faxén operator according (2.21)

(5.4)\begin{equation} F_\beta^{[1]}=F_\beta^{[\infty]}={-}{M}_\beta(\boldsymbol{\xi})={-}8 {\rm \pi}\mu \mathcal{F}_{\gamma \beta} {u}_\gamma(\boldsymbol{\xi}), \end{equation}

where we used the notation $F_\beta ^{[\infty ]}$ to remark that the force $F_{\beta }^{[1]}$ is exactly that experienced by the body if the fluid were unbounded.

The other contribution $\boldsymbol{F}^{[3]}+ \boldsymbol{F}^{[5]} + \cdots$ in (5.3) is the force experienced by the body immersed in the ambient flow $\boldsymbol{v}^{[2]}(\boldsymbol{x}) + \boldsymbol{v}^{[4]}(\boldsymbol{x}) + \cdots$. Therefore,

(5.5)\begin{equation} F_\beta^{[3]}+F_{\beta}^{[5]}+ \cdots={-}8 {\rm \pi}\mu \mathcal{F}_{\gamma \beta} ( v_\gamma^{[2]}(\boldsymbol{\xi})+v_{\gamma}^{[4]}(\boldsymbol{\xi})+ \cdots). \end{equation}

Indicating with $[S]$ the $\infty \times 3$-dimensional matrix collecting all the derivatives of the Stokeslet $\boldsymbol {\nabla }_{\boldsymbol{\beta }_m} S_{a \beta }(\boldsymbol{x}-\boldsymbol{\xi })$ (analogously to the definition (4.4) for $[G]$) and with $[W]$ the $\infty \times 3$-dimensional matrix collecting all the derivatives of the regular part of the Green function $\boldsymbol {\nabla }_{\boldsymbol{\beta }_m}W_{a \beta }(\boldsymbol{x},\boldsymbol{\xi })$, according to (2.8) and (2.9), the matrix $[G]$ can be decomposed as

(5.6)\begin{equation} [G]=[S]+[W] \end{equation}

and the sum of the fields $v_\gamma ^{[2]}(\boldsymbol{\xi })+v_{\gamma }^{[4]}(\boldsymbol{\xi })+ \cdots$ with even values of $k$, (4.9), takes the form

(5.7)\begin{equation} \sum_{k=0}^{\infty} \boldsymbol{v}^{[2 k+2]}(\boldsymbol{x})= \dfrac{ [{M}]^t ([I] - [N])^{{-}1} [W] }{8{\rm \pi}\mu}, \end{equation}

whereas the sum of all the fields corresponding to odd values of $k$, associated with the disturbance field due to the body, is given by

(5.8)\begin{equation} \sum_{k=0}^{\infty} \boldsymbol{v}^{[2k+1]}(\boldsymbol{x})= \dfrac{ [{M}]^t ([I] - [N])^{{-}1} [S]}{8 {\rm \pi}\mu}. \end{equation}

Substituting (5.4), (5.5) and (5.7) into (5.3), a compact representation of the force is achieved

(5.9)\begin{equation} \boldsymbol{F}= \boldsymbol{F}^{[\infty]}- [{M}]^t ([I] - [N])^{{-}1} [N_{(:,0)}], \end{equation}

where the matrix

(5.10)\begin{equation} [N_{(:,0)}]= \left[ \begin{array}{c} \boldsymbol{N}_{(0,0)}\\ \boldsymbol{N}_{(1,0)}\\ \vdots\\ \boldsymbol{N}_{(m,0)}\\ \vdots \end{array} \right] \end{equation}

collecting the entries $N_{\alpha \boldsymbol{\alpha }_m \beta \boldsymbol{\beta }_n}(\boldsymbol{\xi })$ for $n=0$, is exactly the $\infty \times 3$ matrix corresponding to the first three columns of the matrix $[N]$.

The same procedure can be applied to obtain an analogous relation for the torque acting on the body. By (2.22), the torque $\boldsymbol{T}^{[1]}=\boldsymbol{T}^{[\infty ]}$ is provided by the operator $\mathcal {T}_{\gamma \beta }=\varepsilon _{\beta \delta \delta _1} \mathcal {F}_{\gamma \delta \delta _1}$ applied at the ambient flow $\boldsymbol{u}(\boldsymbol{x})$, i.e.

(5.11)\begin{equation} T_\beta^{[1]}=T_\beta^{[\infty]}= \varepsilon_{\beta \delta \delta_1}{M}_{\delta \delta_1}(\boldsymbol{\xi})= 8 {\rm \pi}\mu \mathcal{T}_{\gamma \beta} {u}_\gamma(\boldsymbol{\xi}) \end{equation}

and the remaining term in (5.3) is equal to

(5.12)\begin{equation} T_\beta^{[3]}+T_{\beta}^{[5]}+ \cdots=8 {\rm \pi}\mu \mathcal{T}_{\gamma \beta} ( v_\gamma^{[2]}(\boldsymbol{\xi})+v_{\gamma}^{[4]}(\boldsymbol{\xi})+ \cdots). \end{equation}

Therefore, the total torque is compactly expressed by the equation by

(5.13)\begin{equation} \boldsymbol{T}= \boldsymbol{T}^{[\infty]}+[{M}]^t ([I] - [N])^{{-}1} [L] \end{equation}

with

(5.14)\begin{equation} [L]= \left[ \begin{array}{c} \boldsymbol{L}_{(0)}\\ \boldsymbol{L}_{(1)}\\ \vdots\\ \boldsymbol{L}_{(m)}\\ \vdots \end{array} \right], \end{equation}

where $\boldsymbol{L}_{(m)}$ are the $3^{(m+1)} \times 3$-dimensional matrices with entries $\varepsilon _{\beta \delta \delta _1} N_{\gamma \boldsymbol{\gamma }_m \delta \delta _1}(\boldsymbol{\xi })$, thus

(5.15)\begin{equation} [L]=[N_{(:,1)}]\,\boldsymbol{\varepsilon}^t, \end{equation}

where

(5.16)\begin{equation} \boldsymbol{\varepsilon} = \left( \begin{array}{ccccccccccccccccccccccccccc} 0 & 0 & 0 & 0 & 0 & 1 & 0 & -1 & 0 \\ 0 & 0 & -1 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & -1 & 0 & 0 & 0 & 0 & 0 \end{array} \right) \end{equation}

corresponds to the vectorisation of the Ricci–Levi Civita tensor $\varepsilon _{\beta \delta \delta _1}$ with respect to indices $\delta \delta _1$.

This result can be generalised to the moments: the $n$th-order moment $\bar{\boldsymbol{M}}_{(n)}(\boldsymbol{\xi })$ on the particle in a confined fluid is given by

(5.17)\begin{equation} \bar{\boldsymbol{M}}^t_{(n)}(\boldsymbol{\xi})= \boldsymbol{M}_{(n)}^t(\boldsymbol{\xi})+ [{M}]^t ([I] - [N])^{{-}1} [N_{(:,n)}], \end{equation}

where

(5.18)\begin{equation} [N_{(:,n)}]= \left[ \begin{array}{c} \boldsymbol{N}_{(0,n)}\\ \boldsymbol{N}_{(1,n)}\\ \vdots\\ \boldsymbol{N}_{(m,n)}\\ \vdots \end{array} \right]. \end{equation}

6. Error estimate in truncation

The exact results obtained for velocity field, force and torque in (4.14), (5.9) and (5.13) are expressed in terms of infinite matrices. In practical application, it is not possible to take into account all the entries of these matrices. In fact, in the overwhelming majority of cases, only lower order analytical expressions for Faxén operators and multipoles singularities are available and, moreover, there are no recursive relations able to predict higher-order terms even for the simplest geometries. Furthermore, whenever complex geometries are considered, for which no analytical solutions are available, numerical approaches represent the only feasible alternative in order to evaluate moments and multipoles, and approximations or series truncations become necessary.

Therefore, a central issue in the practical applications of reflection methods is the determination of the order of magnitude of the error committed in the approximations/truncations as function of the geometric dimensionless ratio $\ell _b/\ell _d$. Specifically, consider the error deriving by considering only the first moments and multipoles up to the $K$th order in the exact expressions (4.14), (5.9) and (5.13), hence by substituting to $[M]$ its truncated counterpart $[M_{(0:K)}]$ (where, according the notation introduced in (4.3), $[M_{(0:K)}]$ is the vector collecting all the unbounded moments from the zeroth to the $K$th order) and similarly for the other infinite matrices, $[N] \rightarrow [N_{(0:K,0:K)}]$, $[G] \rightarrow [G_{(0:K)}]$,

(6.1)\begin{equation} \left| \boldsymbol{v}(\boldsymbol{x})-\boldsymbol{u}(\boldsymbol{x}) - \dfrac{ [{M}_{(0:K)}]^t ([I_{(0:K,0:K)}] - [N_{(0:K,0:K)}])^{{-}1} [G_{(0:K)}]}{8{\rm \pi}\mu} \right|. \end{equation}

To this aim, consider the (4.5) rewritten in the form

(6.2)\begin{equation} v_a^{[1]}(\boldsymbol{x})+v_a^{[2]}(\boldsymbol{x}) = \dfrac{[M_{(0:K)}]^t[G_{(0:K)}]}{8{\rm \pi}\mu}+ \dfrac{1}{8{\rm \pi}\mu} \sum_{m=K+1}^{\infty} \dfrac{ \left(\boldsymbol{M}_{(m)}\right)^t \boldsymbol{G}_{(m)} }{m!}. \end{equation}

Enforcing the dimensional analysis developed in Appendix A, specifically (A5) and (A7) providing

(6.3)\begin{equation} M_{\beta \boldsymbol{\beta}_{K+1}}(\boldsymbol{\xi})=\mu U_c O\left( \ell_b^{K+2} \right); \quad \boldsymbol{\nabla}_{\boldsymbol{\beta}_{K+1}} G_{a \beta}(\boldsymbol{x},\boldsymbol{\xi})= O\left(1/\ell_f^{K+2} \right), \end{equation}

the leading-order term in the series on the right-hand side of (6.2) is

(6.4)\begin{equation} |\boldsymbol{M}_{(K+1)}^t \boldsymbol{G}_{(K+1)}| = \mu U_c O\left( \dfrac{ \ell_b}{\ell_f} \right)^{K+2}, \end{equation}

with $U_c$ being the characteristic magnitude of the ambient velocity field. Therefore, truncating the series up to the $K$th order, the order of magnitude of the remainder follows

(6.5)\begin{equation} v_a^{[1]}(\boldsymbol{x})+v_a^{[2]}(\boldsymbol{x}) = \dfrac{[M_{(0:K)}]^t[G_{(0:K)}]}{8{\rm \pi}\mu}+ U_c O\left( \dfrac{\ell_b}{\ell_f} \right)^{K+2}. \end{equation}

The velocity fields due to the next reflections, hence for $k=1$, can be written as

(6.6)\begin{align} v_a^{[3]}(\boldsymbol{x})+v_a^{[4]}(\boldsymbol{x}) &= \dfrac{[M_{(0:K)}]^t[N_{(0:K,0:K)}][G_{(0:K)}]}{8{\rm \pi}\mu} \nonumber\\ &\quad + \dfrac{1}{8{\rm \pi}\mu} \sum_{m=K+1}^{\infty} \, \sum_{n=K+1}^{\infty} \dfrac{ \left(\boldsymbol{M}_{(m)}\right)^t \boldsymbol{N}_{(m,n)}\, \boldsymbol{G}_{(n)} }{m!}, \end{align}

where the leading-order term in the series, estimated via (A5), (A7) and (A12), is

(6.7)\begin{equation} |\boldsymbol{M}_{(K+1)}^t \boldsymbol{N}_{(K+1,K+1)} \boldsymbol{G}_{(K+1)}| = \mu U_c O \left( \dfrac{ \ell_b}{\ell_f} \dfrac{ \ell_b}{2 \ell_d} \right)^{K+2}. \end{equation}

Since $\ell _b/(2 \ell _d) < 1$, the term in (6.7) is always smaller than the term in (6.4), thus

(6.8)\begin{equation} |\boldsymbol{M}_{(K+1)}^t \boldsymbol{N}_{(K+1,K+1)} \boldsymbol{G}_{(K+1)}| = \mu U_c O \left( \dfrac{ \ell_b}{\ell_f} \right)^{K+2} \end{equation}

and (6.6) can be rewritten as

(6.9)\begin{equation} v_a^{[3]}(\boldsymbol{x})+v_a^{[4]}(\boldsymbol{x}) = \dfrac{[M_{(0:K)}]^t[N_{(0:K,0:K)}][G_{(0:K)}]}{8{\rm \pi}\mu}+ U_c\, O\left( \dfrac{\ell_b}{\ell_f} \right)^{K+2}. \end{equation}

Reiterating the same procedure for the higher-order reflected velocity fields, hence for $k>1$, we obtain

(6.10)\begin{align} v_a^{[2 k + 1]}(\boldsymbol{x})+v_a^{[2 k +2]}(\boldsymbol{x}) & = \dfrac{ [M_{(0:K)}]^t [N_{(0:K,0:K)}]^k [G_{(0:K)}]}{8{\rm \pi}\mu} +U_c \, O \left( \dfrac{\ell_b}{\ell_f} \right)^{K+2},\nonumber\\ k & =1,2,\ldots.\end{align}

Therefore, the truncation error committed considering only up to $K$th-order terms in the infinity matrix expressions entering (4.14) satisfies the scaling relations

(6.11)\begin{equation} \boldsymbol{v}(\boldsymbol{x})-\boldsymbol{u}(\boldsymbol{x})= \dfrac{ [{M}_{(0:K)}]^t ([I_{(0:K,0:K)}] - [N_{(0:K,0:K)}])^{{-}1} [G_{(0:K)}]}{8{\rm \pi}\mu}+U_c O\left( \dfrac{\ell_b}{\ell_f} \right)^{K+2}. \end{equation}

A similar analysis can be extended to forces and torques, obtaining

(6.12)\begin{equation} \boldsymbol{F}- \boldsymbol{F}^{[\infty]}={-} [{M}_{(0:K)} ]^t ([I_{(0:K,0:K)}] - [N_{(0:K,0:K)}] )^{{-}1} [N_{(0:K,0)}] + F_c O\left( \dfrac{\ell_b}{\ell_d} \right)^{K+2} \end{equation}

and

(6.13)\begin{equation} \boldsymbol{T}- \boldsymbol{T}^{[\infty]}= [{M}_{(0:K)} ]^t ([I_{(0:K,0:K)}] - [N_{(0:K,0:K)}] )^{{-}1} [L_{(0:K)}] + T_c O\left( \dfrac{\ell_b}{\ell_d} \right)^{K+2}, \end{equation}

where $F_c = \mu \ell _b U_c$ and $T_c=\mu \ell _b^2 U_c$.

The scaling analysis of the truncation error addressed previously can be applied to the approximations of the hydromechanical properties addressed in the literature. In the following, we analyse and discuss, from the point of view of the present theory, the classical literature results of hydromechanics of bodies in confined Stokes flow, extending such expressions to more general BCs, geometries and motions of the bodies.

6.1. The approximation for $K=0$ and Brenner's formula

An explicit approximation for the force acting on an arbitrary body translating in a confined fluid has been first derived by Brenner (Reference Brenner1964a) in terms of the resistance matrix of the body in the unbounded fluid and the regular part of the Green function at the position of the body. Within the present formalism, the Brenner (Reference Brenner1964a) formula for the force is

(6.14)\begin{equation} \boldsymbol{F}=\boldsymbol{F}^{[\infty]} \left({ I} + \dfrac{ \boldsymbol{W}_{(0)} \, \boldsymbol{R} }{8{\rm \pi} \mu} \right)^{{-}1} + F_c O\left( \dfrac{\ell_b}{\ell_d} \right)^{2} ,\end{equation}

where $R_{\beta \gamma }=- 8 {\rm \pi}\mu m_{\beta \gamma }$ is the resistance matrix of the body in the unbounded fluid.

For $K=0$ in (6.11)–(6.13), substituting $\boldsymbol{M}_{(0)}^t=-\boldsymbol{F}^{[\infty ]}$, we have the zeroth-order approximation of the velocity field

(6.15)\begin{equation} \boldsymbol{v}(\boldsymbol{x})-\boldsymbol{u}(\boldsymbol{x})={-} \dfrac{ \boldsymbol{F}^{[\infty]} ({I} - \boldsymbol{N}_{(0,0)})^{{-}1} \boldsymbol{G}_{(0)}}{8{\rm \pi}\mu} +U_c O\left( \dfrac{\ell_b}{\ell_f} \right)^{2} ,\end{equation}

the force

(6.16)\begin{equation} \boldsymbol{F}- \boldsymbol{F}^{[\infty]}=\boldsymbol{F}^{[\infty]}({I} - \boldsymbol{N}_{(0,0)} )^{{-}1} \boldsymbol{N}_{(0,0)} + F_c O\left( \dfrac{\ell_b}{\ell_d} \right)^{2} \end{equation}

and the torque

(6.17)\begin{equation} \boldsymbol{T}- \boldsymbol{T}^{[\infty]}={-}\boldsymbol{F}^{[\infty]}({ I} -\boldsymbol{N}_{(0,0)} )^{{-}1} \boldsymbol{L}_{(0)} + T_c O\left( \dfrac{\ell_b}{\ell_d} \right)^{2} ,\end{equation}

with ${I}$ being the $3 \times 3$ identity matrix.

By the dimensional analysis developed in Appendix B.1, the velocity field (6.15) becomes

(6.18)\begin{equation} \boldsymbol{v}(\boldsymbol{x})=\boldsymbol{u}(\boldsymbol{x}) -\boldsymbol{F}^{[\infty]} \left({ I} +\dfrac{ \boldsymbol{W}_{(0)} \boldsymbol{R} }{8{\rm \pi} \mu} \right)^{{-}1} \dfrac{ \boldsymbol{G}_{(0)}}{8{\rm \pi}\mu} +U_c O\left( \dfrac{\ell_b}{\ell_f} \right)^{2} ,\end{equation}

the force (6.16) returns, as expected, the Brenner's formula

(6.19)\begin{equation} \boldsymbol{F}=\boldsymbol{F}^{[\infty]} \left({ I} + \dfrac{ \boldsymbol{W}_{(0)} \boldsymbol{R} }{8{\rm \pi} \mu} \right)^{{-}1} + F_c O\left( \dfrac{\ell_b}{\ell_d} \right)^{2} \end{equation}

and the torque (6.17) is

(6.20)\begin{equation} \boldsymbol{T}= \boldsymbol{T}^{[\infty]} - \boldsymbol{F}^{[\infty]} \left({ I} + \dfrac{ \boldsymbol{W}_{(0)} \boldsymbol{R} }{8{\rm \pi} \mu} \right)^{{-}1} \dfrac{\boldsymbol{W}_{(0)} \boldsymbol{C} }{8 {\rm \pi}\mu} + T_c O\left( \dfrac{\ell_b}{\ell_d} \right)^{2}, \end{equation}

where $C_{\beta \gamma }= 8{\rm \pi} \mu \varepsilon _{\beta \delta \delta _1} m_{\gamma \delta \delta _1}(\boldsymbol{\xi,\boldsymbol{\xi }})$ is the coupling matrix between forces and rotations of the body in the unbounded fluid (Happel & Brenner Reference Happel and Brenner1983; Procopio & Giona Reference Procopio and Giona2024).

Equations (6.18), (6.19) and (6.20), requiring solely the Green function of the confinement and the grand-resistance matrix of the body, provide the first-order terms associated with any flow past a body immersed in a confined fluid. In the case the body translates with velocity $\boldsymbol{U}$ in a quiescent fluid, we can consider the disturbance flow associated with an ambient flow past the still body $\boldsymbol{u}(\boldsymbol{x})=-\boldsymbol{U}$. Therefore, from (6.18) and considering $\boldsymbol{F}^{[\infty ]}=-\boldsymbol{U} \boldsymbol{R}$, the velocity field around the translating body is

(6.21)\begin{equation} \boldsymbol{w}(\boldsymbol{x})= \boldsymbol{U} \boldsymbol{R} \left({ I} +\dfrac{ \boldsymbol{W}_{(0)} \boldsymbol{R} }{8{\rm \pi} \mu} \right)^{{-}1} \dfrac{ \boldsymbol{G}_{(0)}}{8{\rm \pi}\mu} +U_c O\left( \dfrac{\ell_b}{\ell_f} \right)^{2}. \end{equation}

By (6.19), the force is

(6.22)\begin{equation} \boldsymbol{F}={-}\boldsymbol{U} \boldsymbol{R} \left({ I} + \dfrac{ \boldsymbol{W}_{(0)} \boldsymbol{R} }{8{\rm \pi} \mu} \right)^{{-}1} + F_c O\left( \dfrac{\ell_b}{\ell_d} \right)^{2} ,\end{equation}

whereas, by considering that $\boldsymbol{T}^{[\infty ]}=-\boldsymbol{U}\boldsymbol{C}$, from (6.20) the torque is

(6.23)\begin{equation} \boldsymbol{T}={-}\boldsymbol{U} \left(\boldsymbol{C}-\boldsymbol{R} \left(I + \dfrac{ \boldsymbol{W}_{(0)} \boldsymbol{R} }{8 {\rm \pi}\mu} \right)^{{-}1} \dfrac{\boldsymbol{W}_{(0)} \boldsymbol{C} }{8 {\rm \pi}\mu} \right) + T_c O\left( \dfrac{\ell_b}{\ell_d} \right)^{2}. \end{equation}

If the particle rotates (without translating) with angular velocity $\boldsymbol{\omega }$, the force on the particle in the unbounded fluid is given by $\boldsymbol{F}^{[\infty ]}=- \boldsymbol{\omega } \boldsymbol{C}^t$ (Happel & Brenner Reference Happel and Brenner1983) and, therefore, following the same procedure as previously, the velocity field is

(6.24)\begin{equation} \boldsymbol{w}(\boldsymbol{x})=\boldsymbol{\omega} \boldsymbol{C}^t \left({ I} +\dfrac{ \boldsymbol{W}_{(0)} \boldsymbol{R} }{8{\rm \pi} \mu} \right)^{{-}1}\dfrac{ \boldsymbol{G}_{(0)}}{8{\rm \pi}\mu} +U_c O\left( \dfrac{\ell_b}{\ell_f} \right)^{2}, \end{equation}

the force on the body takes the expression

(6.25)\begin{equation} \boldsymbol{F}={-}\boldsymbol{\omega}\,\boldsymbol{C}^t \left({ I} + \dfrac{ \boldsymbol{W}_{(0)} \boldsymbol{R} }{8{\rm \pi} \mu} \right)^{{-}1} + F_c O\left( \dfrac{\ell_b}{\ell_d} \right)^{2} \end{equation}

and the torque, considering $\boldsymbol{T}^{[\infty ]}=-\boldsymbol{\omega } \boldsymbol{\varOmega }$, can be written as

(6.26)\begin{equation} \boldsymbol{T}={-}\boldsymbol{\omega} \left[ \boldsymbol{\varOmega} - \boldsymbol{C}^t \left({ I} + \dfrac{ \boldsymbol{W}_{(0)} \boldsymbol{R} }{8{\rm \pi} \mu} \right)^{{-}1} \,\dfrac{ \boldsymbol{W}_{(0)} \boldsymbol{C}}{8 {\rm \pi}\mu} \right] + T_c O\left( \dfrac{\ell_b}{\ell_d} \right)^{2}, \end{equation}

where $\boldsymbol{\varOmega }$, having entries $\varOmega _{\alpha \beta }=- 8 {\rm \pi}\mu \varepsilon _{\alpha \gamma \gamma _1} \varepsilon _{\beta \delta \delta _1} m_{\gamma \gamma _1 \delta \delta _1}(\boldsymbol{\xi },\boldsymbol{\xi })$, is the angular resistance matrix.

6.2. Extended Swan and Brady's approximation for rigid motion

In obtaining the approximate expressions (6.19) and (6.20), valid to the order $(\ell _b/\ell _d)$, we have neglected the higher-order terms in the zeroth-order Faxén operator for the force and in the first-order Faxén operator for the torque (see the derivation in Appendix B.1). Supposing that these Faxén operators are exactly known, it is possible to obtain more accurate expressions for the force and the torque on a rigid moving body. This has been found by Swan & Brady (Reference Swan and Brady2007, Reference Swan and Brady2010) in the case of confined spherical bodies with no-slip BCs. Specifically, expressing their result in the mobility representation, Swan and Brady provided the velocity $\boldsymbol{U}$ and angular velocity $\boldsymbol{\omega }$ of a sphere under the action of a force $\boldsymbol{F}$ and a torque $\boldsymbol{T}$

(6.27)\begin{equation} \boldsymbol{U}={-} \boldsymbol{F} \boldsymbol{R}^{{-}1} (\boldsymbol{I}-\boldsymbol{R}^{{-}1} \boldsymbol{\varPhi})+ \boldsymbol{T} \boldsymbol{\varOmega}^{{-}1} \boldsymbol{\varPsi}^t \boldsymbol{R}^{{-}1} \end{equation}

and

(6.28)\begin{equation} \boldsymbol{\omega}= \boldsymbol{F} \boldsymbol{R}^{{-}1} \boldsymbol{\varPsi} \boldsymbol{\varOmega}^{{-}1}- \boldsymbol{T} \boldsymbol{\varOmega}^{{-}1} (\boldsymbol{I}-\boldsymbol{\varOmega}^{{-}1} \boldsymbol{\varTheta}), \end{equation}

where

(6.29)\begin{align} \varPhi_{\alpha\beta } &=\left.- 8{\rm \pi} \mu \mathcal{F}_{\gamma \beta} \mathcal{F}_{\gamma' \alpha} W_{\gamma \gamma'}(\boldsymbol{\xi},\boldsymbol{\xi}') \right|_{\boldsymbol{\xi}'=\boldsymbol{\xi}}, \end{align}

(6.30)\begin{align} \varPsi_{\alpha \beta} &=\left. 8{\rm \pi} \mu \mathcal{T}_{\gamma \beta} \mathcal{F}_{\gamma' \alpha} W_{\gamma \gamma'}(\boldsymbol{\xi},\boldsymbol{\xi}') \right|_{\boldsymbol{\xi}'=\boldsymbol{\xi}} \end{align}

and

(6.31)\begin{equation} \varTheta_{\beta \alpha} ={-} 8{\rm \pi} \left.\mu \mathcal{T}_{\gamma \beta} \mathcal{T}_{\alpha' \delta} W_{\gamma \alpha'}(\boldsymbol{\xi},\boldsymbol{\xi}') \right|_{\boldsymbol{\xi}'=\boldsymbol{\xi}}, \end{equation}

with the Faxén operators being those of a sphere with no-slip BCs.

By using the theory developed previously, it is possible to extend the Swan and Brady expressions to bodies with generic shape and generic reciprocal BCs. Furthermore, it is possible to evaluate the error committed by adopting the Swan and Brady approach to estimate the hydrodynamic interaction. In Appendix B.2, we obtain that the force and torque acting on a body rigid moving with velocity $\boldsymbol{U}$ and angular velocity $\boldsymbol{\omega }$ can be expressed as

(6.32)\begin{equation} \boldsymbol{F}={-} \boldsymbol{U} \left(\boldsymbol{R}+\left(\boldsymbol{I}-\boldsymbol{\varPhi}\,\boldsymbol{R}^{{-}1}\right)^{{-}1}\boldsymbol{\varPhi}\right) - \boldsymbol{\omega} \left(\boldsymbol{C}^t+\boldsymbol{\varPsi}^t\left(\boldsymbol{I}-\boldsymbol{\varPhi} \boldsymbol{R}^{{-}1}\right)^{{-}1}\right)+ O\left( \dfrac{\ell_b}{\ell_d} \right)^3 \end{equation}

and

(6.33)\begin{equation} \boldsymbol{T}\!=\! {-}\boldsymbol{U} \left(\!\boldsymbol{C}\!+\!\left(\boldsymbol{I}\! -\! \boldsymbol{\varPhi} \boldsymbol{R}^{{-}1}\right)^{{-}1}\boldsymbol{\varPsi}\! \right)\! -\! \boldsymbol{\omega} \left(\boldsymbol{\varOmega}\!+\! \boldsymbol{\varTheta}\! +\! \boldsymbol{\varPsi}^t \boldsymbol{R}^{{-}1} \left(1\! -\! \boldsymbol{\varPhi}\boldsymbol{R^{{-}1}}\right)^{{-}1}\boldsymbol{\varPsi}\!\right)\! +\! O\left( \dfrac{\ell_b}{\ell_d} \right)^3, \end{equation}

where the Faxén operators entering $\boldsymbol{\varPhi }$, $\boldsymbol{\varPsi }$, $\boldsymbol{\varTheta }$, defined in (6.29)–(6.31), are those of a body with generic shape and reciprocal BCs.

In the case of a spherical body, for which $\boldsymbol{C}=\boldsymbol{0}$, (6.32) and (6.33) become

(6.34)\begin{equation} \boldsymbol{F}={-} \boldsymbol{U} \left(\boldsymbol{R}+\left(\boldsymbol{I}-\boldsymbol{\varPhi}\,\boldsymbol{R}^{{-}1}\right)^{{-}1}\boldsymbol{\varPhi}\right) - \boldsymbol{\omega} \boldsymbol{\varPsi}^t\left(\boldsymbol{I}-\boldsymbol{\varPhi} \boldsymbol{R}^{{-}1}\right)^{{-}1}+ O\left( \dfrac{\ell_b}{\ell_d} \right)^4 \end{equation}

and

(6.35)\begin{align} \boldsymbol{T}={-} \boldsymbol{U} \left(\boldsymbol{I}-\boldsymbol{\varPhi} \boldsymbol{R}^{{-}1}\right)^{{-}1}\boldsymbol{\varPsi} - \boldsymbol{\omega} \left(\boldsymbol{\varOmega}+\boldsymbol{\varTheta}+ \boldsymbol{\varPsi}^t \boldsymbol{R}^{{-}1} \left(1-\boldsymbol{\varPhi}\boldsymbol{R^{{-}1}}\right)^{{-}1}\boldsymbol{\varPsi}\right)+ O\left( \dfrac{\ell_b}{\ell_d} \right)^5. \end{align}

From these equations it is possible to obtain the Swan and Brady ‘mobility’ representation (6.27) and (6.28) by inverting the block grand-resistance matrix (Bhatia Reference Bhatia1997) and considering that, for a spherical particle, $\boldsymbol{\varPhi }=O(\ell _b/\ell _d)$, $\boldsymbol{\varPsi }=O(\ell _b/\ell _d)^2$ and $\boldsymbol{\varTheta }=O(\ell _b/\ell _d)^3$. As shown in Brenner (Reference Brenner1964b), there exists a vast class of bodies for which, choosing the centre of hydrodynamic reaction as centre of rotation, the coupling between translation and rotation vanishes, then $\boldsymbol{C}=\boldsymbol{0}$. Providing that the geometrical moments entering the Faxén operators are evaluated with respect to the centre of hydrodynamic rotation, for all these bodies it is possible to employ (6.34) and (6.35) obtaining, hence, an improved error in the truncation.

Equations (6.32)–(6.33) and (6.34)–(6.35) are fundamental scaling relations permitting to obtain good approximations for the resistance on bodies (especially if possess special symmetries) in confined systems by the knowledge of solely the zeroth- and first-order (for torque) Faxén operators. However, although such corrections improve the approximations in evaluating force and torques on body in confined fluids, the scaling analysis indicates that the approximate relation (6.34) (and, hence, the Swan and Brady expressions) cannot provide correctly the term $O(\ell _b/\ell _d)^4$ for the force. In the Faxén's expressions for a sphere translating near a plane wall and for a sphere translating between two parallel plane walls (see Happel & Brenner Reference Happel and Brenner1983, p. 327), a non-vanishing term $O(\ell _b/\ell _d)^4$ has been found, whereas in Swan & Brady (Reference Swan and Brady2007, Reference Swan and Brady2010) the term $O(\ell _b/\ell _d)^4$ is considered vanishing, but the terms $O(\ell _b/\ell _d)^5$ (which fortuitously agrees with that obtained by Faxén performing high-order reflections) has been obtained. In the next section, we apply the exact expression (6.12) for the force on a sphere translating near a plane wall truncated to the order $O(\ell _b/\ell _d)^5$ and we find that, in accordance with the Faxén result, the term $O(\ell _b/\ell _d)^4$ is not vanishing.

7. A sphere near a plane wall

Next, consider the archetypical hydrodynamic problem of a sphere with radius $R_p$ and centre at distance $h$ from a plane, as depicted in figure 3. Furthermore, consider Navier-slip BCs (2.5) at the surface of the sphere and no-slip BC at the surface of the plane wall.

Figure 3. Schematic representation of a sphere near a plane wall.

For this hydrodynamic system, both the Faxén operator of the body and the Green function of the confinement, required to construct analytically the $[N]$-matrix in (4.6), are available in the literature up to the second order. Specifically, the zeroth-, first- and second-order Faxén operators with pole at the centre of the sphere can be found in Procopio & Giona (Reference Procopio and Giona2024), whereas the Green function and its derivatives in the semi-space have been derived by Blake (Reference Blake1971); Blake & Chwang (Reference Blake and Chwang1974); Procopio & Giona (Reference Procopio and Giona2023). Both the Faxén operators and the multipoles evaluated at the position of the sphere, useful for the computation of the $[N]$-matrix, are reported in the supplementary material available at https://doi.org/10.1017/jfm.2024.651 expressed in the Cartesian coordinate system $(\xi _1,\xi _2,\xi _3)$ represented in figure 3.

According to the definition (3.17), the entries of the $[N]$-matrix, representing the hydrodynamic interaction between the sphere and the plane, are obtained by applying the Faxén operators to the regular part of the multipoles. Therefore, by the knowledge of the Faxén operators up to the second order, it is possible to construct the matrix

(7.1)\begin{equation} [N_{(0:2,0:2)}]= \left[ \begin{array}{cccccc} \boldsymbol{N}_{(0,0)} & \boldsymbol{N}_{(0,1)} & \dfrac{\boldsymbol{N}_{(0,2)}}{2!}\\ \boldsymbol{N}_{(1,0)} & \boldsymbol{N}_{(1,1)} & \dfrac{\boldsymbol{N}_{(1,2)}}{2!}\\ {\boldsymbol{N}_{(2,0)}} & {\boldsymbol{N}_{(2,1)}} & \dfrac{\boldsymbol{N}_{(2,2)}}{2!} \end{array} \right]. \end{equation}

The entries of the submatrix

(7.2) \begin{equation} \boldsymbol{N}_{(0,0)}= \left[ \begin{array}{cccccc} N_{1\ 1} & N_{1\ 2} & N_{1\ 3}\\ N_{2\ 1} & N_{2\ 2} & N_{2\ 3}\\ N_{3\ 1} & N_{3\ 2} & N_{3\ 3} \end{array} \right] \end{equation}

expressed in the Cartesian coordinate system $(\xi _1,\xi _2,\xi _3)$, are

(7.3)\begin{equation} \left. \begin{aligned} N_{1\, 1} & =\left. \mathcal{F}_{\gamma' 1} W_{\gamma' 1}(\boldsymbol{\xi}',\boldsymbol{\xi})\right|_{\boldsymbol{\xi}'=\boldsymbol{\xi}=(0,0,h)}= { \frac{9 R_p(1+2 \hat{\lambda}) }{16 h (1+3 \hat{\lambda} )}- \frac{R_p^3}{16 h^3 (1+3 \hat{\lambda})} },\\ N_{2\, 2} & =\left. \mathcal{F}_{\gamma' 2} W_{\gamma' 2}(\boldsymbol{\xi}',\boldsymbol{\xi})\right|_{\boldsymbol{\xi}'=\boldsymbol{\xi}=(0,0,h)} = { \frac{9 R_p(1+2 \hat{\lambda}) }{16 h (1+3 \hat{\lambda} )}- \frac{R_p^3}{16 h^3 (1+3 \hat{\lambda})} },\\ N_{3\, 3} & =\left. \mathcal{F}_{\gamma' 3} W_{\gamma' 3}(\boldsymbol{\xi}',\boldsymbol{\xi})\right|_{\boldsymbol{\xi}'=\boldsymbol{\xi}=(0,0,h)} = { \frac{9 {R_p} (1+2 \hat{\lambda} ) }{8 h (1+3 \hat{\lambda} )}- \frac{R_p^3}{4 h^3 (1+3 \hat{\lambda})} },\\ N_{1\, 2} & = N_{1\, 3} = N_{2\, 1} = N_{3\, 1}= N_{3\, 2}= N_{2\, 3}=0, \end{aligned} \right\} \end{equation}

with $\hat {\lambda }=\lambda /R_p$ being the dimensionless slip length of the fluid–sphere interface. The entries of the other submatrices entering (7.1) can be obtained by the Faxén operators and multipoles reported in the supplementary material.

7.1. Force and torque on a translating sphere near a plane wall

Consider a spherical body translating with velocity $\boldsymbol{U}$, hence an ambient flow $\boldsymbol{u}(\boldsymbol{x})=-\boldsymbol{U}$. The force acting onto the sphere is given by (6.12). Assuming $K=3$, we have

(7.4)\begin{equation} \boldsymbol{F}= \boldsymbol{F}^{[\infty]}- [{M}_{(0:3)} ]^t ([I_{(0:3,0:3)}] - [N_{(0:3,0:3)}] )^{{-}1} [N_{(0:3,0)}] + F_c O\left( \dfrac{\ell_b}{\ell_d} \right)^{5}, \end{equation}

where

(7.5)\begin{equation} \boldsymbol{F}^{[\infty]}={-}\boldsymbol{M}^t_{(0)}={-}6 {\rm \pi}\mu R_p \left( \dfrac{1+2\hat{\lambda}}{1+3\hat{\lambda}} \right) \boldsymbol{U} \end{equation}

is the well-known Hadamard–Rybczynski force, $\boldsymbol{M}_{1}=\boldsymbol{M}_{(3)}=0$ and, as shown in Appendix C.1,

(7.6)\begin{equation} \boldsymbol{M}_{(2)}^t={-}\dfrac{R_p^2}{3(1+2\hat{\lambda})}\boldsymbol{F}^{[\infty]} \boldsymbol{I}_{(0,2)}, \end{equation}

where $\boldsymbol{I}_{(0,2)}$ is a matrix obtained by the vectorisation of the multi-index $({\beta \beta _1 \beta _2})$ of the tensor $\delta _{\gamma \beta }\delta _{\beta _1 \beta _2}$. As shown in Appendix C.1, since $\boldsymbol{M}_{(3)}=0$, the submatrices $\boldsymbol{N}_{(3,q)}$ and $\boldsymbol{N}_{(q,3)}$ (with $q=0,1,2,3$) entering $[N_{(0:3,0:3)}]$ and $[N_{(0:3,0)}]$ in (7.4) are immaterial and can be assumed vanishing (in point of fact, we show in Appendix C.1, that only the submatrix $\boldsymbol{N}_{(0,p)}$, with $p=0,1,2$, contribute to the force obtained by (7.4)). After substituting the entries of the $N$-matrix in (7.4) and after some algebra we obtain

(7.7)\begin{equation} \boldsymbol{F}= \boldsymbol{F^{[\infty]}} \left( \begin{array}{ccc} \dfrac{1}{1-\beta_{{\parallel}}} & 0 & 0 \\ 0 & \dfrac{1}{1-\beta_{{\parallel}}} & 0 \\ 0 & 0 & \dfrac{1}{1-\beta_{{\perp}}} \end{array} \right) +O \left(\dfrac{R_p}{h} \right)^5, \end{equation}

where

(7.8)\begin{equation} \beta_{{\parallel}}= \dfrac{9 R_p}{16 h} \left( \dfrac{1+2\hat{\lambda}}{1+3\hat{\lambda}} \right) -\dfrac{R_p^3}{8 h^3(1+3 \hat{\lambda})} +\dfrac{45 R_p^4}{256 h^4} \left( \dfrac{(1+2\hat{\lambda})^2}{(1+3 \hat{\lambda})(1+5\hat{\lambda})} \right) \end{equation}

and

(7.9)\begin{equation} \beta_{{\perp}}= \dfrac{9 R_p}{8 h} \left( \dfrac{1+2\hat{\lambda}}{1+3\hat{\lambda}} \right) -\dfrac{R_p^3}{2 h^3(1+3 \hat{\lambda})} +\dfrac{135 R_p^4}{256 h^4} \left( \dfrac{(1+2\hat{\lambda})^2}{(1+3 \hat{\lambda})(1+5\hat{\lambda})} \right). \end{equation}

For $\hat {\lambda }=0$ (i.e. in the no-slip case), existing literature results can be recovered. In fact, the force $F_\perp$ on a sphere translating perpendicularly to the plane wall with velocity $U_\perp$ is in perfect agreement with the Taylor series expansion of the exact solution obtained by Brenner (Reference Brenner1961), which reads

(7.10)\begin{equation} \dfrac{F_{{\perp}}}{6{\rm \pi}\mu R_p U_\perp}=1+\frac{9 {R_p}}{8 h}+\frac{81 {R_p}^2}{64 h^2}+\frac{473 {R_p}^3}{512 h^3}+ \frac{4113 {R_p}^4}{4096 h^4}+O\left(\dfrac{R_p}{h}\right)^5 \end{equation}

and $\beta _{\parallel }$ reduces to the same expression obtained by Faxén (see Happel & Brenner Reference Happel and Brenner1983, p. 327).

As expected by the dimensional analysis performed in § 6, by using the approximate Brenner's relation (6.22), only the first terms $O(R_p/h)$ on the right-hand side of (7.8) and (7.9) are correctly obtained, whereas by using the Swan and Brady approximations (6.34) also the second terms $O(R_p/h)^3$ are correctly evaluated, but the third $O(R_p/h)^4$ term is erroneously vanishing. Figures 4 and 5 depict the results obtained in (7.7)–(7.9) and by the approximated relations developed in §§ 6.1 and 6.2 compared with the exact results (such as those deriving by solving the Goren (Reference Goren1979) equations) and to FEM simulations (in the cases the exact results are not available in the literature). From the visual inspection of figures 4 and 5, it readily follows that (7.7) provides an accurate representation of the force on a spherical body valid for gaps $\delta =(h-R_p) \gtrsim R_p$ or even smaller, if one considers all the terms for $\beta _{\parallel }$ and $\beta _{\perp }$ entering (7.8) and (7.9). Further improvements of the expansion, taking into account higher-order Faxén operators, could increase the range of validity of the theoretical expressions to smaller gap values.

Figure 4. Dimensionless force on a spherical body translating perpendicularly to a planar wall for different slip lengths on the surface of the sphere. Solid black lines represent exact results obtained by solving the equations provided by Goren (Reference Goren1979), red dashed-dotted lines depict (7.7) taking into account all the terms, green dashed lines represent (7.7) where the terms $O(R_p/h)^4$ are neglected (corresponding to the approximated relation (6.34)) and blue dotted lines represent (7.7) where the terms $O(R_p/h)^3$ are neglected (corresponding to the approximated relation (6.22)). For complete slip ($\hat {\lambda }=\infty$), blue dotted and green dashed lines are practically coincident.

Figure 5. Dimensionless force on a spherical body translating parallel to a plane wall for different slip lengths on the surface of the sphere. The solid black line represents the exact no-slip result obtained by solving the equations provided by O'Neill (Reference O'Neill1964), black circles are the results of FEM simulations, red dashed-dotted lines depict (7.7) taking into account all the terms, green dashed lines represent (7.7) where the terms $O(R_p//h)^4$ are neglected (corresponding to the approximated relation (6.34)) and blue dotted lines represent (7.7) where the terms $O(R_p/h)^3$ are neglected (corresponding to the approximated relation (6.22)). For complete slip ($\hat {\lambda }=\infty$), blue and green dashed lines are practically coincident.

It can be observed that there is a larger error in the solution obtained by the Swan and Brady approximation (6.34) to the order $O(R_p/h)^3$ compared with the solution obtained by the Brenner approximation (6.22) to the order $O(R_p/h)$, despite the higher-order reflection considered. This additional error arises from the negative terms in the denominator of the expressions (7.8) and (7.9), associated with the dipole potential of the Stokes solution for the sphere translating in the unbounded fluid, which induce a slight ‘boost’ upon reflection by the plane.

According (6.13), the torque on the translating sphere with velocity $\boldsymbol{U}$ truncated to $K=3$ reads

(7.11)\begin{equation} \boldsymbol{T}= \boldsymbol{T}^{[\infty]}+ [{M}_{(0:3)} ]^t ([I_{(0:3,0:3)}] - [N_{(0:3,0:3)}] )^{{-}1} [L_{(0:3)}] + T_c O\left( \dfrac{\ell_b}{\ell_d} \right)^{5} . \end{equation}

Substituting the expressions for the entries of the $N$-matrix and considering the error estimated in Appendix C.2 (or, equivalently, considering the (7.20) and the symmetry of the grand-resistance matrix), after some algebra one obtains

(7.12)\begin{equation} \boldsymbol{T}= R_p \boldsymbol{F^{[\infty]}} \left( \begin{array}{ccc} 0 & -\gamma & 0 \\ \gamma & 0 & 0 \\ 0 & 0 & 0 \\ \end{array} \right) +O \left(\dfrac{R_p}{h} \right)^6 \end{equation}

or, by (7.5),

(7.13)\begin{equation} \boldsymbol{T}={-}6 {\rm \pi}\mu R_p^2 \left( \dfrac{1+2\hat{\lambda}}{1+3\hat{\lambda}} \right)\boldsymbol{U} \left( \begin{array}{ccc} 0 & -\gamma & 0 \\ \gamma & 0 & 0 \\ 0 & 0 & 0 \\ \end{array} \right) +O \left(\dfrac{R_p}{h} \right)^6, \end{equation}

where

(7.14)\begin{equation} \gamma= \dfrac {1}{8} \left(\dfrac{R_p}{h} \right)^4 \left( {\dfrac{1}{(1+2\hat{\lambda})( 1+ 3 \hat{\lambda})}-\frac{3}{8} \left( \frac{R_p}{h} \right) \frac{(1+5\hat{\lambda} +15 \hat{\lambda}^2)}{(1+3\hat{\lambda})^2 (1+5 \hat{\lambda})} }\right). \end{equation}

As expected from the dimensional analysis performed in § 6.1 (6.23), the $K=0$-order approximation cannot predict the leading-order term $O(R_p/h)^4$ as in (7.13)–(7.14). In fact, since the coupling matrix $\boldsymbol{C}$ of a sphere in the unbounded fluid is vanishing, the torque on a translating sphere near a plane would vanishing according to this approximation. Using the approximation (6.35) it is possible to obtain correctly the leading-order term $O(R_p/h)^4$ in (7.14), but not the term $O(R_p/h)^5$.

The results obtained by (7.14) truncated to the orders $O(R_p/h)^4$ and $O(R_p/h)^5$ are compared in figure 6 with the exact solutions provided by O'Neill (Reference O'Neill1964) for no-slip BCs assumed on the surface of the sphere and with FEM simulations performed for $\hat {\lambda }=0.1, 0.5, 1$.

Figure 6. Dimensionless torque on a spherical body translating parallel to a plane wall for different slip lengths on the surface of the sphere. The solid black line represents the exact no-slip result obtained by solving the equations by O'Neill (Reference O'Neill1964), black circles are the results of FEM simulations, red dashed-dotted lines correspond to (7.14) taking into account all the terms and green dashed lines correspond to (7.14) neglecting the term $O(R_p/h)^5$, corresponding to the approximated relation (6.35).

7.2. Torque and force on a rotating sphere near a plane wall

The torque truncated to $K=4$ on a sphere rotating near a plane wall is obtained by considering

(7.15)\begin{equation} \boldsymbol{T}= \boldsymbol{T}^{[\infty]}+ [{M}_{(0:4)} ]^t ([I_{(0:4,0:4)}] - [N_{(0:4,0:4)}] )^{{-}1} [L_{(0:4)}] + T_c O\left( \dfrac{\ell_b}{\ell_d} \right)^{6}, \end{equation}

where $\boldsymbol{M}_{(0)}=\boldsymbol{M}_{(2)}=\boldsymbol{M}_{(4)}=0$ and, as shown in Appendix C.3,

(7.16)\begin{equation} \boldsymbol{M}_{(1)}=\dfrac{\boldsymbol{T}^{[\infty]} \boldsymbol{\varepsilon}}{2}. \end{equation}

In Appendix C.3, we show that $\boldsymbol{M}_{(3)}$, $\boldsymbol{L}_{(3)}$, $\boldsymbol{N}_{(q,3)}$ and $\boldsymbol{N}_{(3,q)}$ (with $q=0,1,2,3$) are immaterial for this approximation order and, hence, they can be assumed to be zero in (7.15). By substituting the other entries of the $N$-matrix, we obtain

(7.17)\begin{equation} \boldsymbol{T}= \boldsymbol{T}^{[\infty]} \left( \begin{array}{ccc} {1+\nu_{{\parallel}}} & 0 & 0 \\ 0 & {1+\nu_{{\parallel}}} & 0 \\ 0 & 0 & {1+\nu_{{\perp}}} \\ \end{array} \right) +O \left(\dfrac{R_p}{h} \right)^6, \end{equation}

where

(7.18)\begin{equation} \nu_{{\parallel}}= \dfrac{15 R_p^3}{16 h^3} \left( \dfrac{1}{1+3\hat{\lambda}} \right) \end{equation}

and

(7.19)\begin{equation} \nu_{{\perp}}= \dfrac{ R_p^3}{8 h^3} \left( \dfrac{1}{1+3\hat{\lambda}} \right). \end{equation}

The result expressed by (7.17) can be obtained by applying the extended Swan and Brady approximation (6.35). The graph of (7.17) is depicted in figures 7 and 8, compared with the exact expressions for a sphere with no-slip BCs provided by Jeffery (Reference Jeffery1915) and Dean & O'Neill (Reference Dean and O'Neill1963) and with FEM simulations for $\hat {\lambda }=0.1,0.5,1$.

Figure 7. Dimensionless torque on a spherical body rotating with angular velocity perpendicular to a plane wall for different slip lengths on the surface of the sphere. The solid black line represents the exact no-slip result obtained by solving the equations provided by Dean & O'Neill (Reference Dean and O'Neill1963), black circles are the results of FEM simulations and red dashed-dotted lines represent (7.17) which, in the present case, is equivalent to (6.35).

Figure 8. Dimensionless torque on a spherical body rotating with angular velocity parallel to a plane wall for different slip lengths on the sphere. The solid black line represents the exact no-slip result obtained by solving the equations provided by Dean & O'Neill (Reference Dean and O'Neill1963), black circles are the results of FEM simulations and red dashed-dotted lines represent (7.17) which, in the present case, is equivalent to (6.35).

The force acting on a rotating sphere near a plane wall can be deduced from

(7.20)\begin{equation} \boldsymbol{F}= [{M}_{(0:4)} ]^t ([I_{(0:4,0:4)}] - [N_{(0:4,0:4)}] )^{{-}1} [N_{(0:4,0)}] + F_c O\left( \dfrac{\ell_b}{\ell_d} \right)^{6}. \end{equation}

The same result admits a more straightforward derivation by considering the linearity between $\boldsymbol{F}$ and $\boldsymbol{\omega }$ expressed by the coupling matrix $\bar{\boldsymbol{C}}$

(7.21)\begin{equation} \boldsymbol{F}={-} \boldsymbol{\omega} \bar{\boldsymbol{C}} +O\left( \dfrac{R_p}{h}\right)^6 . \end{equation}

Enforcing the symmetry of the grand-resistance matrix (see Appendix A in Procopio & Giona (Reference Procopio and Giona2022) or for the case of a body in a confined fluid considering Navier-slip BCs), we have

(7.22)\begin{equation} \bar{\boldsymbol{C}} = \bar{\boldsymbol{D}}^t, \end{equation}

where $\bar{\boldsymbol{D}}$ is the coupling matrix entering the expression $\boldsymbol{T}=- \boldsymbol{U} \bar{\boldsymbol{D}}$ in (7.13). Therefore,

(7.23)\begin{equation} \boldsymbol{F}={-}6 {\rm \pi}\mu R_p^2 \left( \dfrac{1+2\hat{\lambda}}{1+3\hat{\lambda}} \right)\boldsymbol{\omega} \left( \begin{array}{ccc} 0 & \gamma & 0 \\ -\gamma & 0 & 0 \\ 0 & 0 & 0 \\ \end{array} \right) +O \left(\dfrac{R_p}{h} \right)^6, \end{equation}

where $\gamma$ is given in (7.14). Analogously to (7.12)

(7.24) \begin{equation} \boldsymbol{F}= \dfrac{\boldsymbol{T^{[\infty]}}}{R_p} \left( \begin{array}{ccc} 0 & \delta & 0 \\ -\delta & 0 & 0 \\ 0 & 0 & 0 \\ \end{array} \right) +O \left(\dfrac{R_p}{h} \right)^6, \end{equation}

where

(7.25)\begin{equation} \delta=\left( 1+2\hat{\lambda} \right)\dfrac{3\gamma}{4}. \end{equation}

8. Force on a translating prolate spheroid near a plane wall

In this section, we investigate the effect of the shape and of the orientation of a body in a confined fluid on the hydrodynamic interactions between the body and the confinement in light of the theory developed here. To this aim, we consider a prolate spheroid translating in the Stokes fluid with velocity $\boldsymbol{U}$ near a plane wall, at distance $h$ between its centroid and the plane. Unlike the spherical case addressed in the previous section, two additional geometrical parameters should be introduced: the eccentricity $e$ of the spheroid, accounting for the effects of the shape of the body on its hydromechanics, and the angle $\theta$ between the symmetry axis of the prolate spheroid and the plane (see figure 9), accounting for the effects of the orientation of the body. In order to highlight these effects, the translational motion of the spheroid without rotation is considered, although a full hydromechanic analysis of this system can be developed within the present approach. To obtain the hydrodynamic force acting on the spheroid, we employ the approximate expressions obtained in §§ 6.1 and 6.2, and we show that, owing to these approximations, accurate expressions for the hydromechanics of particles in confined fluid can be derived solely from the knowledge of the unbounded resistance matrix or the lower-order Faxén operators.

Figure 9. Schematic representation of a prolate spheroid near a plane wall. The prolate spheroid translates with velocity $\boldsymbol{U}$, forming an angle $\phi$ with the major axis of the spheroid. Here $\boldsymbol{F}_d$ is the hydrodynamic drag force experienced by the body aligned with $\boldsymbol{U}$ and $\boldsymbol{F}_L$ the hydrodynamic lift force orthogonal to $\boldsymbol{U}$.

We consider a prolate spheroid defined by the surface equation

(8.1)\begin{equation} \dfrac{\zeta_1^2}{a^2}+\dfrac{\zeta_2^2+\zeta _3^2}{b^2}=1, \end{equation}

where $(\zeta _1,\zeta _2,\zeta _3)$ is a Cartesian coordinate system with origin at the centroid, $a$ the semi-length of the major axis and $b$ the semi-length of the minor axis (see figure 9).

The resistance matrix of a prolate spheroid with no-slip BCs in the unbounded fluid, expressed in the coordinate system $(\zeta _1,\zeta _2,\zeta _3)$, is $\boldsymbol{R}=16 {\rm \pi}\mu c \boldsymbol{A}(e)$ (Kim & Karrila Reference Kim and Karrila2005) where ${e}$ is the eccentricity of the spheroid, $c=ae$ the focal length and

(8.2)\begin{align} & \boldsymbol{A}(e)=\nonumber \\ & \left( \begin{array}{ccc} \dfrac{e^2}{(1+e^2)\log{\left(\dfrac{1+e}{1-e}\right)}-2e} & 0 & 0\\ 0 & \dfrac{2e^2}{(1-3e^2)\log{\left(\dfrac{1+e}{1-e}\right)}-2e} & 0\\ 0 & 0 & \dfrac{2e^2}{(1-3e^2)\log{\left(\dfrac{1+e}{1-e}\right)}-2e} \end{array} \right). \end{align}

Substituting the resistance matrix into (6.14), the first-order approximation is straightforward

(8.3)\begin{equation} \boldsymbol{F}= \boldsymbol{F}^{[\infty]} \left(I+2\dfrac{c}{h} { \boldsymbol{w}}\left(\theta \right) \boldsymbol{A}(e) \right)^{{-}1} +O\left(\dfrac{a}{h} \right)^{^2}, \end{equation}

where $\boldsymbol{F}^{[\infty ]}=-16 {\rm \pi}\mu c \boldsymbol{U}\boldsymbol {\cdot } \boldsymbol{A}(e)$ and ${ \boldsymbol{w}}(\theta )$, having as its entries the regular part of the Green function $-{W}_{\alpha \beta }(\boldsymbol{\zeta }/ h,\boldsymbol{\zeta }/h)$ at the position of the centroid $\boldsymbol{\zeta }=(0,0,0)$, reads

(8.4)\begin{equation} w_{\alpha \beta} (\theta)={-} \left( \begin{array}{ccc} \dfrac{3}{8} (3-\cos (2 \theta )) & 0 & \dfrac{3}{8} \sin (2 \theta ) \\ 0 & \dfrac{3}{4} & 0 \\ \dfrac{3}{8} \sin (2 \theta ) & 0 & \dfrac{3}{8} (3+\cos (2 \theta )) \\ \end{array} \right). \end{equation}

To obtain the higher-order terms providing the hydrodynamic force, we can employ the zeroth-order Faxén operator for the prolate spheroid available in the literature (Hasimoto Reference Hasimoto1983; Kim Reference Kim1985; Kim & Karrila Reference Kim and Karrila2005) by using the extended Swan and Brady relations obtained in § 6.2. Specifically, the zeroth-order Faxén operator can be expressed as

(8.5)\begin{equation} \mathcal{F}_{\beta \gamma}= \dfrac{R_{\beta \gamma}}{16 {\rm \pi}\mu c}\int_{\boldsymbol{\varGamma}} d\boldsymbol{\varGamma} \left( 1+\dfrac{(1-e^2)(c^2-\zeta_1^2)\Delta_{\zeta}}{4e^2} \right), \end{equation}

where $\Delta _{\zeta }$ is the Laplacian acting on the coordinates $(\zeta _1,\zeta _2,\zeta _3)$, the curve $\boldsymbol{\varGamma }$ is the line segment between the two foci at $\zeta _1=\pm c$ of the spheroid with parametrisation

(8.6)\begin{equation} \boldsymbol{\varGamma}(s)= \left\{ \begin{array}{c} \zeta_1(s)=s,\\ \zeta_2(s)=0,\\ \zeta_3(s)=0, \end{array} \right. \quad s \in [{-}c,c] \end{equation}

and $\textrm {d}\boldsymbol{\varGamma } \equiv \,\textrm {d}s$ is its measure element.

The geometrical moments entering the Faxén operator (8.5) are expressed with respect to the centroid (being the centre of hydrodynamic reaction) as pole. For the symmetry of the spheroid, a rotation around this point is uncoupled to the translations, thus providing $\boldsymbol{C}=\boldsymbol{0}$. Therefore, we can use (6.34) assuming $\boldsymbol{\omega }=\boldsymbol{0}$, which, after some algebra, attains the more compact form

(8.7)\begin{equation} \boldsymbol{F}={-}\boldsymbol{U} \boldsymbol{R}(\boldsymbol{I}-\boldsymbol{R}^{{-}1}\boldsymbol{\varPhi})^{{-}1}+O\left(\dfrac{a}{h}\right)^4. \end{equation}

From (8.7), in order to obtain the hydrodynamic force acting on the spheroid by (8.7), it is necessary to evaluate the matrix $\boldsymbol{\varPhi }$ defined by (6.29), which, considering the Faxén operator (8.5), reads

(8.8)\begin{align} \varPhi_{\beta \gamma}&={-} \dfrac{R_{\delta' \beta} R_{\delta \gamma}}{32 {\rm \pi}\mu c^2} \left[ \int_{{-}c}^c \,{\rm d}s' \int_{{-}c}^c \,{\rm d}s\, W_{\delta' \delta}\left(\boldsymbol{\zeta}'(s'),\boldsymbol{\zeta}(s)\right) \right.\nonumber\\ &\quad+ \left. \dfrac{1-e^2}{4e^2} \int_{{-}c}^c \,{\rm d}s' \int_{{-}c}^c \,{\rm d}s\, (c^2-s^2) \Delta_{\zeta} W_{\delta' \delta}\left(\boldsymbol{\zeta}'(s'),\boldsymbol{\zeta}(s)\right) \right. \nonumber\\ &\quad+ \left. \dfrac{1-e^2}{4e^2} \int_{{-}c}^c \,{\rm d}s' \int_{{-}c}^c \,{\rm d}s\, (c^2-s'^2) \Delta_{\zeta'} W_{\delta' \delta}\left(\boldsymbol{\zeta}'(s'),\boldsymbol{\zeta}(s)\right) \right]+O(a/h)^4, \end{align}

where the primed subscripts ${\delta }'=1,2,3$ refer to the coordinate system at the point $\boldsymbol{\zeta }'$.

Substituting (8.8) into (8.7), we finally obtain

(8.9)\begin{equation} \boldsymbol{F}= \boldsymbol{F}^{[\infty]} \left(I+\dfrac{c}{2 h} { \boldsymbol{P}}\left(\dfrac{c}{h},\theta \right) \, \boldsymbol{A}(e)+\left(\dfrac{c}{h}\right)^3 \left(\dfrac{1-e^2}{8 e^2} \right) { \boldsymbol{Q}}\left(\dfrac{c}{h},\theta \right) \boldsymbol{A}(e) \right)^{{-}1} +O\left(\dfrac{a}{h} \right)^{^4}, \end{equation}

where the entries of $\boldsymbol{P}({c}/{h},\theta )$ and $\boldsymbol{Q}({c}/{h},\theta )$ are, respectively,

(8.10)\begin{align} P_{\delta' \delta}\left(\dfrac{c}{h},\theta \right) & = \int_{{-}1}^1 \,{\rm d}s'\, \int_{{-}1}^1 \,{\rm d}s\, W_{\delta' \delta}\left(\boldsymbol{\zeta}'(s')\dfrac{c}{h},\boldsymbol{\zeta}(s)\dfrac{c}{h}\right) \end{align}

(8.11)\begin{align} Q_{\delta' \delta}\left(\dfrac{c}{h},\theta \right) & = \left[ \int_{{-}1}^1 \,{\rm d}s' \int_{{-}1}^1 \,{\rm d}s (1-s^2) \Delta_{\zeta} W_{\delta' \delta}\left(\boldsymbol{\zeta}'(s') \dfrac{c}{h},\boldsymbol{\zeta}(s)\dfrac{c}{h}\right) \right. \nonumber\\ & \quad + \left. \int_{{-}1}^1 \,{\rm d}s' \int_{{-}1}^1 \,{\rm d}s (1-s'^2) \Delta_{\zeta'} W_{\delta' \delta}\left(\boldsymbol{\zeta}'(s') \dfrac{c}{h},\boldsymbol{\zeta}(s) \dfrac{c}{h}\right) \right]. \end{align}

To perform the derivatives and integrals along the segment between the foci of the spheroid entering (8.8) (or the normalised integrals entering (8.10) and (8.11)) independently of the orientation of the prolate spheroid, it is fundamental to employ the representation of the Green function invariant both at the field and the pole points, corresponding to the ‘bi-invariant’ form of the regular part of the Green function in the semi-space obtained in Procopio & Giona (Reference Procopio and Giona2023)

(8.12) \begin{align} {W}_{b \beta} (\boldsymbol{x},\boldsymbol{\zeta})&= {S}_{b \beta} (\boldsymbol{x}- \boldsymbol{\zeta}^*) -(\boldsymbol{\zeta}^* -\boldsymbol{\zeta})\boldsymbol{\cdot} \boldsymbol{n}\, J_{\beta \gamma*}\nonumber\\ &\times \left[ n_{\delta^*} \boldsymbol{\nabla}_{\gamma*} {S}_{b \delta^*}(\boldsymbol{x}- \boldsymbol{\zeta}^*) -\dfrac{(\boldsymbol{\zeta}^* -\boldsymbol{\zeta}) \boldsymbol{\cdot} {\boldsymbol n}}{2} \Delta_{\zeta^*} {S}_{a \gamma^*}({\boldsymbol x}- \boldsymbol{\zeta}^*) \right], \end{align}

with $\boldsymbol{x} \equiv \boldsymbol{\zeta }'$ and where $\boldsymbol{n}=(\sin {\theta },0,\cos {\theta })$ is the normal to the plane wall inward to the fluid, $\boldsymbol{J}= I - 2 \boldsymbol{n} \otimes \boldsymbol{n}$ the mirror operator (that in Procopio & Giona (Reference Procopio and Giona2023) has been shown to coincide with the parallel propagator between the fluid space and the reflected space), $\boldsymbol{\zeta }^*=\boldsymbol{J}\boldsymbol {\cdot }\boldsymbol{\zeta }-2 h\, \boldsymbol{n}$ is the reflected point by the plane and the starred subscripts, such as ${\beta ^*}$, refers to the coordinate system at the mirror point $(\zeta ^*_1,\zeta ^*_3,\zeta ^*_3)$.

Solving the integrals (8.10) and (8.11) with the regular part of the Green function (8.11) (normalised accordingly) and expanding the results in Taylor series of $(c/h)$, we obtain

(8.13)\begin{equation} \boldsymbol{F}\!=\! 16 {\rm \pi}\mu c \boldsymbol{U}\boldsymbol{A}(e) \left(\! I\! +\! \dfrac{c}{2 h} { \boldsymbol{p}}_{0} \left(\theta \right) \boldsymbol{A}(e)\! +\! \left(\dfrac{c}{2 h}\right)^3 \left(\! { \boldsymbol{p}}_{2} (\theta)\! +\! \dfrac{1-e^2}{4e^2} \boldsymbol{q}_0( \theta )\! \right) \boldsymbol{A}(e) \right)^{{-}1} \! +\! O\left(\dfrac{a}{h} \right)^{^4}, \end{equation}

where $\boldsymbol{p}_0( \theta )=4 \,\boldsymbol{w}( \theta )$,

(8.14)\begin{equation} \boldsymbol{p}_2= \left( \begin{array}{ccc} \dfrac{60 \cos{(2 \theta )}-9 \cos{(4 \theta)}+5}{96} & 0 & \dfrac{3}{4} \sin \theta \cos ^3\theta \\ 0 & \dfrac{ 9 \cos{(2 \theta)}+1}{24} & 0 \\ \dfrac{3}{4} \sin{\theta} \cos^3{\theta} & 0 & \cos{(2\theta)}+\dfrac{9 \cos (4 \theta )+55}{96} \end{array} \right) \end{equation}

and

(8.15)\begin{equation} \boldsymbol{q}_{0}= \left( \begin{array}{ccc} \dfrac{20}{3}-4 \cos (2 \theta ) & 0 & 8 \sin (\theta ) \cos (\theta ) \\ 0 & \dfrac{8}{3} & 0 \\ 8 \sin (\theta ) \cos (\theta ) & 0 & 4 \cos (2 \theta )+\dfrac{20}{3} \end{array} \right). \end{equation}

From (8.13), we can evince that the second-order term $O(a/h)^2$ in the inverse matrix entering the expression of the force is vanishing for any orientation $\theta$ of the prolate spheroid. Therefore, comparing this result with (8.3), we can conclude that the error committed by using the $K=0$ approximation is $O(a/h)^3$, hence smaller than $O(a/h)^2$ expected from general considerations.

When particle dynamics is considered, it might be useful to express the force in the fixed reference frame of the plane. Specifically, considering the coordinate system with origin on the plane (see figure 9)

(8.16)\begin{equation} \begin{cases} \xi_1 = \zeta_1\, \cos{\theta} -\zeta_3\, \sin{\theta},\\ \xi_2 = \zeta_2,\\ \xi_3 = \zeta_1 \sin{\theta} + \zeta_3 \cos{\theta} -h, \end{cases} \end{equation}

the entries of the matrices $\boldsymbol{p}_0 ( \theta )$, $\boldsymbol{p}_2 ( \theta )$, $\boldsymbol{q}_0 ( \theta )$, expressed in this system, read

(8.17)\begin{gather} \boldsymbol{p}_0(\theta)=\left( \begin{array}{ccc} -3 & 0 & 0 \\ 0 & -3 & 0 \\ 0 & 0 & -6 \end{array} \right), \end{gather}

(8.18)\begin{gather} \boldsymbol{p}_2(\theta)= \left( \begin{array}{ccc} \dfrac{ 11 \cos (2 \theta )+3}{24} & 0 & -\dfrac{ \sin (2 \theta )}{6} \\ 0 & \dfrac{9 \cos (2 \theta )+1}{24} & 0 \\ -\dfrac{ \sin (2 \theta )}{6} & 0 & \dfrac{ 7 \cos (2 \theta )+3}{6} \end{array} \right), \end{gather}

(8.19)\begin{gather} \boldsymbol{q}_0(\theta)=\left( \begin{array}{ccc} \dfrac{8}{3} & 0 & 0 \\ 0 & \dfrac{8}{3} & 0 \\ 0 & 0 & \dfrac{32}{3} \end{array} \right). \end{gather}

An expression for the mobility, conceptually analogous to (8.13), has been proposed by Mitchell & Spagnolie (Reference Mitchell and Spagnolie2015) using a Swan and Brady approach, directly extended to the ellipsoids. However, it can be observed that the entries of the mobility matrix reported in Mitchell & Spagnolie (Reference Mitchell and Spagnolie2015) do not match the entries of the mobility of the sphere obtained by Swan & Brady (Reference Swan and Brady2007) (or obtainable by inverting the resistance matrix provided in § 7.1) in the limit of vanishing eccentricity $e\rightarrow 0$. On the other hand, by inverting the resistance matrix entering (8.13), the limit of spherical bodies is correctly matched. For instance, the velocity of the spheroid $U$ along the plane under the effect of an external force $F^{ ext}$ parallel to the plane reads

(8.20) \begin{gather} \dfrac{6{\rm \pi} \mu a U}{F^{ ext}} = \dfrac{3}{8} \left(\frac{\cos ^2(\theta )}{e A_{11}(e)}+\dfrac{\sin ^2(\theta )}{e A_{33}(e)}\right) -\left(\dfrac{a}{h} \right) \frac{9}{16}+ \left(\dfrac{a}{h} \right)^3 \dfrac{ \left(11 e^2 \cos (2 \theta )-13 e^2+16\right)}{128}, \end{gather}

where $A_{11}(e)$ and $A_{33}(e)$ are the entries of the matrix $\boldsymbol{A}(e)$ in (8.2) in the coordinate system $( \zeta _1,\zeta _2,\zeta _3)$. Since $e A_{11}(e) \rightarrow 3/8$ and $e A_{33}(e) \rightarrow 3/8$ in the limit $e\rightarrow 0$, (8.20) yields the mobility

(8.21)\begin{equation} \dfrac{6{\rm \pi} \mu a U}{F^{ ext}} = 1 -\left(\dfrac{a}{h} \right) \frac{9}{16}+ \left(\dfrac{a}{h} \right)^3 \frac{ 1}{8} \end{equation}

of a sphere in this limit. Equation (8.20) differs from that provided by Mitchell & Spagnolie (Reference Mitchell and Spagnolie2015) for terms $O(a/h)^0$ and $O(a/h)^3$. The discrepancy is likely due to the fact that the relations derived by Mitchell & Spagnolie (Reference Mitchell and Spagnolie2015) were obtained by differentiating only at the field point of the Green function to obtain higher-order terms, also in those cases the derivative at the pole point was necessary (for example to obtain terms $\boldsymbol {\nabla }_\beta \boldsymbol {\nabla }_{\beta '}W(\boldsymbol{\xi },\boldsymbol{\xi }')$), since the pole point cannot be differentiated in the Blake's formulation (Blake Reference Blake1971) of the Green function bounded by a no-slip plane (see Liron & Mochon (Reference Liron and Mochon1976) or Spagnolie & Lauga (Reference Spagnolie and Lauga2012) for a discussion on this point). As regards this technical but important issue, bitensorial calculus represents the proper geometrical setting for addressing the differentiation problems, owing to its bi-invariance both with respect to field and pole point coordinates.

The comparison of the theoretical relations (8.3) and (8.13) and the numerical FEM simulations depicted in figures 10 and 11, shows that the theoretical approximated relations provide accurate results up to $h \approx 2 a$, independently of the orientation of the spheroid. More specifically, figure 10 depicts the drag force $F_d$ (hence, the force aligned with the translatory direction, see figure 9) experienced by the spheroid translating with velocity parallel to the plane along the axis $\xi _1$ and intensity $U$. It can be observed that, for small values of the eccentricity $e$, the approximation (8.13), truncated to the order $O(a/h)^4$, provides a worse approximation with respect to (8.3) truncated to a lower order. This phenomenon could indeed be expected by considering the case of a spherical body addressed in § 7.1, where it is shown that the approximation for the drag force is improved once the next correction term $O(R_p/h)^4$ is accounted for. However, it is possible to observe in figure 10 that, for higher values of the eccentricity $e$, the error decreases with the inclusion of the next correction term $O(a/h)^3$. As shown in figure 11 for the lift force (i.e. the force orthogonal to the translation of the body, see figure 9) experienced by prolate spheroid moving parallel to the plane along the axis $\xi _1$, (8.13) provides a better approximation with respect to the lower-order truncation (8.3) for any value of eccentricity $e$.

Figure 10. Dimensionless drag force acting on prolate spheroids with eccentricity $e=0.3,0.6,0.9$ translating parallel to the plane wall along the axis $\xi _1$ with intensity $U$ as a function of the orientation angle $\theta$ (see figure 9) for different distances $h$ between the centroid and the plane wall. Symbols represent FEM simulations, green dashed lines the $O(a/h)^3$ approximation (8.3) and red lines the $O(a/h)^4$ approximation (8.12).

Figure 11. Dimensionless lift force acting on prolate spheroids with eccentricity $e=0.3,0.6,0.9$ translating parallel to the plane wall along the axis $\xi _1$ with intensity $U$ as a function of the orientation angle $\theta$ (see figure 9) for different distances $h$ between the centroid and the plane wall. The arrows indicate increasing distances $h/a=2, 4, 6,\infty$. Symbols represents FEM simulations, black dashed lines the unbounded case, green dashed lines the $O(a/h)^3$ approximation (8.3) and red lines the $O(a/h)^4$ approximation (8.12).

In contrast to the case of a spheroid in a unbounded flow, the plane wall generates a lift force $\boldsymbol{F}_L$ even when the translation occurs along the principal axes of the spheroid, as shown in figure 12(a). This implies that, if we want to move the spheroid along its principal axes, the external force applied to it cannot be aligned with the axes itself. Specifically, the lift force always opposes the approach of the body to the plane, pointing towards the plane when the direction of motion points out of the plane and vice versa. In figure 12(b–d), it is possible to observe that the hydrodynamic effects of the plane wall on the resistance of the spheroid determine a lift force that may have an opposite sign with respect to the case of the unbounded flow. This occurs when the orientation of the spheroid is $\theta ={\rm \pi} /2$ or $\theta ={\rm \pi} /4$, whereas for $\theta =0$ the spheroid experiences a lift force in the same direction of the unbounded case. This implies that if the spheroid is perpendicular to the plane and we want to move it in a direction forming an angle $\phi ={\rm \pi} /4$ with its axis of symmetry, an external force in a direction forming a smaller angle with the axis should be applied. This is opposite to the unbounded case, where the force direction must have a greater angle.

Figure 12. Dimensionless lift force acting on a prolate spheroid with eccentricity $e=0.5$ at distance $h=3 a$. Panel (a) depicts the lift force acting on a spheroid translating along its principal axis as a function of the orientation angle $\theta$. Panels (b–d) depict the lift force as a function of the angle $\phi$ between the velocity of the body and its symmetry axis (see figure 9) for different orientations of the spheroid. Specifically, in panel (b) the spheroid is normal to the plane, in panel (c) forms an angle $\theta ={\rm \pi} /4$ with the plane and in panel (d) is parallel to the plane. The black dashed lines represent the force on prolate spheroids in the unbounded flow, the green dashed lines the $O(a/h)^3$ approximation (8.3) and the red lines the $O(a/h)^4$ approximation (8.12).