2.1. System parameters

Our model system consists of two spheres submerged in a fluid of density $\rho$ and viscosity $\mu$. We denote the fluid domain $\varOmega$ and the domain of sphere $k$ as $\varOmega _k$ for $k = 1$, 2. We let the radius and density of sphere $k$ be $a_k$ and $\rho _k$, respectively. As shown in figure 1, each of the spheres is oscillating along the line connecting their centres as if connected by a massless rod of variable length applying equal and opposite forces. We refer to the internal force, denoted with superscript $(i)$, applied to sphere 2 as $F_2^{(i)}$ and to sphere 1 as $F_1^{(i)} = -F_2^{(i)}$. While the flow induced by this motion appears time reversible at leading order, over many cycles of oscillation small but finite inertial effects give rise to steady time-averaged drift of the two-sphere system at velocity $U$. We seek to describe the flow velocity and pressure $\boldsymbol {u} = \boldsymbol {u}(\boldsymbol {x},t)$ and $p=p(\boldsymbol {x},t)$ as functions of position $\boldsymbol {x}$ and time $t$ to precisely describe the hydrodynamic mechanisms giving rise to this steady motion and determine the form of $U$ in terms of the system parameters. Some of the symbols used in this work are compiled and presented in table 1 for the reader's reference.

Table 1.

List of symbols used in this paper.

We align the $z$ axis so that it passes through both sphere centres, letting $z_k = z_k(t)$ be the instantaneous position of sphere $k$, and specifying $z_2 > z_1$ and $a_1 < a_2$ so the positive $z$ direction points toward the larger sphere. We adopt the reference frame of the swimmer, so that the time-averaged sphere positions $\smash {{\bar {z}}_k}$ are constants and the time-dependent velocity of each sphere $U_k = \dot {z_k}$ is periodic and zero averaged. In particular, we let

(2.1)\begin{equation} U_k(t) = {\rm i}\omega\hat{z}_k \,{\rm e}^{{\rm i} \omega t}, \end{equation}

where $\omega$ is the frequency of oscillation and $\hat {z}_k$ is the complex oscillation amplitude of sphere $k$. Here and throughout, we implicitly take complex expressions to be equal to their real parts. Generally, for an arbitrary time-dependent field $\psi (t)$, we will denote the zeroth and first Fourier coefficients with a bar $(\bar {\psi })$ and hat $(\hat {\psi })$, respectively.

The positive separation distance $\smash {\zeta = z_2 - z_1}$ is thus an oscillatory function

(2.2a–c)\begin{equation} \zeta(t) = d + \varDelta \,{\rm e}^{{\rm i}\omega t},\quad d = \bar{z}_2-\bar{z}_1,\quad \varDelta = \hat{z}_2-\hat{z}_1, \end{equation}

where we have introduced the time-averaged separation distance $d = \bar {\zeta }$ and oscillation amplitude $\varDelta = \hat {\zeta }$ for ease of notation. Symmetry arguments require that any steady flow must depend only on the magnitude of oscillation $\smash {|\varDelta |}$, not its phase.

Without loss of generality, we assume that the separation distance rate of change $\smash {{{\rm d}\varDelta }/{{\rm d} t}}$ has a Fourier amplitude $\smash {{\rm i}\omega \varDelta \propto 1}$, so that the complex arguments of the velocity Fourier amplitudes (2.1) correspond to phase lags with respect to the sphere separation rate of change.

These amplitudes are given by

(2.3a,b)\begin{equation} \hat{U}_1 = {\rm i} \omega \hat{z}_1 = |\omega \hat{z}_1| \exp({{\rm i}(\omega t + {\rm \pi}+ \phi_1)}),\quad \hat{U}_2 = {\rm i}\omega \hat{z}_2 = |\omega \hat{z}_2| \exp({{\rm i}(\omega t - \phi_2)}), \end{equation}

and we define the phase difference

(2.4)\begin{equation} \phi = \phi_1 + \phi_2 \end{equation}

between the spheres’ extremal positions. The sign of $\phi$ identifies which sphere leads the other: when $\phi > 0$, sphere 1 reaches its maximum velocity before sphere 2; when $\phi = 0$, the spheres are exactly 180$^\circ$ out of phase; and when $\phi < 0$, sphere 2 leads sphere 1. This relationship is illustrated in figure 1(c). Shifting the argument of $\varDelta$ corresponds to shifting both trajectories in figure 1(a) or rotating all of the vectors in figure 1(c) through the same angle, neither of which can have any effect on the time averages obtained by integrating over a period of oscillation. Thus, we expect at leading order $U \sim \smash {\varDelta \varDelta ^*}$, since this is the simplest function of $\varDelta$ that can be constructed that is independent of its phase. We also expect $U$ may depend on $\phi$, $|\hat {U}_1|$ and $|\hat {U}_2|$, as these values are also preserved by the aforementioned shifts.

2.2. Dimensionless system

Introducing characteristic length, velocity, density and stress scales $\smash {L= a_2}$, $\smash {V = L\omega }$, $\rho$ and $\smash {\varSigma = \mu V/L}$, we non-dimensionalise the model system. In the following, variables and parameters should be assumed dimensionless unless otherwise specified. Our choice of the length scale is based on the primary contribution to the Stokes drag of the dimer. We will explicitly represent factors of the dimensionless radius $a_2 = 1$.

The fluid system can be reduced to two non-dimensional parameters. Firstly, we introduce $\smash {M^2 = \rho \omega a_2^2 /\mu }$ relating the time partial derivative in the Navier–Stokes equations to viscous stresses $\varSigma$. This is the same dimensionless number and notation used to describe a single oscillating sphere by Riley (Reference Riley1966), although we let $M \in \mathbb {R}$ while Riley defines it as a complex number. The quantity also appears in Stokes’ second problem as the decay rate of a boundary layer with increasing height above an oscillating surface (Acheson Reference Acheson1990), and it takes the same form as the Womersley number Wo describing pulsing flow through pipes (Nguyen et al. Reference Nguyen, Oza, Abouezzi, Sun, Childress, Frederick and Ristroph2021). In settings similar to the current work, it has also been referred to as a scale parameter $s^2$ (Felderhof Reference Felderhof2016) and oscillatory Reynolds number ${\textit {Re}}_\omega$ (Lauga Reference Lauga2007).

Secondly, we introduce the ratio of oscillation amplitude to system size $\smash {\varepsilon = |\varDelta | / L}$, which is the inverse of the Strouhal number. The scaling of the sphere velocities (2.1) shows the dimensionless sphere velocities $U_k \sim \varepsilon$, so the Reynolds number describing flow about the oscillating spheres is $\smash {\textit {Re} = \varepsilon M^2}$. In the following we will assume that $\smash {|\varDelta | \ll a_1,a_2 \ll d}$, so that the amplitude of oscillation is small and the separation distance is large compared with the sphere sizes. In this regime, the flow is naturally described by the parameters $M^2$ and $\varepsilon$, but at higher oscillation amplitudes, as in the work of Dombrowski et al. (Reference Dombrowski, Jones, Katsikis, Bhalla, Griffith and Klotsa2019) and Dombrowski & Klotsa (Reference Dombrowski and Klotsa2020), ${\textit {Re}}$ and $\varepsilon$ are a more convenient set of independent parameters.

Finally, in consideration of the Stokes case, we define the parameter $S^2 = 2\rho _2 \omega a_2^2/3\mu$, representing the relative magnitude of the spheres’ inertia and viscous stresses in the fluid. The factor of 2/3 is included for comparison between the solid and fluid inertial cases. Using this convention, in the fluid inertial case with $\rho _1=\rho _2=1$, the effective mass (physical mass plus added mass effect) of sphere $k$ is $2 {\rm \pi}M^2 a_k^3$. In the solid inertial case $M^2 = 0$, there is no added mass effect and the mass of sphere $k$ is $2 {\rm \pi}S^2 a_k^3$. Thus, using the respective parameters for these two cases will yield a direct comparison between the effective masses of the spheres.

The flow obeys the Navier–Stokes equations

(2.5a,b)\begin{equation} M^2 \left(\frac{{\partial} \boldsymbol{u}}{{\partial} t}+ \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} \right) ={-}\boldsymbol{\nabla} p + \nabla^2 \boldsymbol{u},\quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{equation}

subject to the boundary conditions

(2.6a,b)\begin{equation} \boldsymbol{u}\left(\boldsymbol{x}_k \right) = U_k \boldsymbol{e}_{z}, \quad \lim_{|\boldsymbol{x}| \to \infty} \boldsymbol{u}(\boldsymbol{x}) ={-}U \boldsymbol{e}_{z}, \end{equation}

where $\boldsymbol {x}_k \in \partial \varOmega _k$ denotes an arbitrary position on the boundary of sphere $k$. The right-hand side of (2.5a,b) can be written as the divergence of the Cauchy stress tensor ${\boldsymbol{\mathsf{T}}}$,

(2.7a,b)\begin{equation} {\boldsymbol{\mathsf{T}}} ={-}p {\boldsymbol{\mathsf{I}}} + 2 \mu {\boldsymbol{\mathsf{E}}},\quad {\boldsymbol{\mathsf{E}}} = \tfrac12\left[\boldsymbol{\nabla} \boldsymbol{u} + \left(\boldsymbol{\nabla} \boldsymbol{u}\right)^\intercal\right], \end{equation}

where ${\boldsymbol{\mathsf{I}}}$ is the identity tensor and ${\boldsymbol{\mathsf{E}}}$ is the rate-of-strain tensor.

As shown in figure 1, the spheres are subjected to a vertical flow-mediated force $F_k^{(f)}$ – including contributions from the Stokes drag, Basset force and added mass effect – in addition to the force pair $F_k^{(i)}$ applied by the central connector, given by

(2.8)\begin{equation} F_k^{(f)} = \int_{\partial\varOmega_k} \boldsymbol{e}_{z} \boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}} \boldsymbol{\cdot} \boldsymbol{e}_{n} \,{\rm d} S, \end{equation}

where $\boldsymbol {e}_{n}$ is the normal vector pointing into the fluid. The net vertical force on sphere $k$ can be calculated from the acceleration as

(2.9)\begin{equation} F_k^{(n)} = {\rm i}M^2 \frac{\rho_k}{\rho} \left(\frac{4}{3} {\rm \pi}a_k^3\right) U_k = {\rm i} 2{\rm \pi} S^2 a_k^3 \frac{\rho_k}{\rho_2} U_k, \end{equation}

and at all times each sphere satisfies

(2.10)\begin{equation} F_k^{(n)} = F_k^{(f)} + F_k^{(i)}. \end{equation}

Since $F_1^{(i)} =-F_2^{(i)}$, the pair of spheres satisfies

(2.11)\begin{equation} F_1^{(n)} - F_1^{(f)} = F_2^{(f)} - F_2^{(n)}. \end{equation}

2.3. Series representation

As mentioned above, the sphere velocities (2.1) scale as $\smash {U_k \sim \varepsilon \ll 1}$, suggesting the advective term in the Navier–Stokes equations $\smash {\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla }\boldsymbol {u} \sim \varepsilon ^2}$ will be small and may, as a product of oscillatory functions, have a steady component. As such, the flow can be described as a double sum over powers of $\varepsilon$ and Fourier modes. See Appendix A for a detailed derivation. The approach follows that of Felderhof & Jones (Reference Felderhof and Jones1994), who applied it to a single-sphere geometry. Ignoring all terms with no effect on the $\textit {O}\left (\varepsilon ^2\right )$ swim speed, we adopt the ansätz

(2.12a,b)\begin{equation} \boldsymbol{u} = \varepsilon \hat{\boldsymbol{u}} \,{\rm e}^{{\rm i}t} + \varepsilon^2 \bar{\boldsymbol{u}},\quad p = \varepsilon \hat{p} \,{\rm e}^{{\rm i}t} + \varepsilon^2 \bar{p}, \end{equation}

so that $(\hat {\boldsymbol {u}}, \hat {p})$ are complex fields describing the Fourier amplitudes of a leading-order oscillatory flow and $(\bar {\boldsymbol {u}},\bar {p})$ the steady flow field generated by inertial effects as described above. A ${\rm \pi}$-periodic flow is also induced at the same order as the steady flow, which it does not influence. For this reason, we omit its description as in the work of Felderhof & Jones (Reference Felderhof and Jones1994). Since we have chosen an inertial reference frame, $U = \varepsilon ^2 \bar {U}$. Similarly, the periodic sphere velocities are decomposed as $U_k(t) = \varepsilon \hat {U}_k \,{\rm e}^{{\rm i}t}$, where $\hat {U}_k$ are the non-dimensionalisation of the complex Fourier amplitudes (2.3a,b), implying the net force is $F_k^{(n)} = \varepsilon \hat {F}_k^{(n)} \,{\rm e}^{{\rm i} t}$. The convention of using hats and bars to denote first and zeroth Fourier amplitudes should now be understood to include this normalisation by $\varepsilon$ and $\varepsilon ^2$, respectively, so that the amplitudes are $\textit {O}\left (1\right )$.

The form of (2.12a,b) shows a potential inconsistency: if $\bar {U}$ grows large at high $M^2$, there may be steady flow at order 1 or even order 0, violating the assumptions under which the ansätz was introduced. However, it is shown in § 4 that the swim speed $\bar {U}$ approaches a finite value $\bar {U}_\infty$ as $M \to \infty$ rather than growing unboundedly. Other similar analyses in the one-sphere geometry (Felderhof & Jones Reference Felderhof and Jones1994, Reference Felderhof and Jones2017) have also shown that the steady flow remains bounded at high inertia. We note that the ansätz analysis does break down at large enough $M^2$ or $\varepsilon$. However, this is due to turbulent effects that are not represented in this laminar description as opposed to an inconsistent set of assumptions.

Substituting the expansion (2.12a,b) into the Navier–Stokes equations (2.5a,b) shows the Fourier amplitudes are described by

(2.13a,b)\begin{equation} \left(\nabla^2 - {\rm i}M^2\right)\hat{\boldsymbol{u}} = \boldsymbol{\nabla} \hat p,\quad \boldsymbol{\nabla} \boldsymbol{\cdot} \hat{\boldsymbol{u}} = 0. \end{equation}

These are a complex version of the Brinkman equations, a combination of the Stokes equations and Darcy flow where both viscous stresses and frictional drag force are non-negligible (Durlofsky & Brady Reference Durlofsky and Brady1987). In this case, the ‘drag term’ $-{\rm i}M^2 \boldsymbol {u}$ is imaginary, arising not from the effects of some porous structure but from the acceleration of the unsteady fluid 90$^\circ$ out of phase. The oscillation gives rise to a boundary layer of width $1/M$ over which the resulting vorticity is diffused (Riley Reference Riley1966).

The steady component of the corresponding advective term can be interpreted as a Reynolds stress $\bar {{\boldsymbol{\mathsf{R}}}}$ driving the second-order steady flow according to the Stokes equations

(2.14a,b)\begin{equation} \nabla^2 \bar{\boldsymbol{u}} = \boldsymbol{\nabla}\bar{p} - \boldsymbol{\nabla} \boldsymbol{\cdot} \bar{{\boldsymbol{\mathsf{R}}}},\quad \bar{{\boldsymbol{\mathsf{R}}}} ={-}\frac{M^2}{2} \hat{\boldsymbol{u}} \otimes \hat{\boldsymbol{u}}^*, \end{equation}

where the asterisk denotes a complex conjugate. The factor of 1/2 in the Reynolds stress arises from considering the real part of a product of complex exponentials, as described in Appendix A.

The sphere boundaries $ {\partial }\varOmega _k$ are moving. To obtain a time-independent system, we derive flow constraints on the time-averaged boundaries $\overline { {\partial }\varOmega }_k$ by Taylor expanding the flow fields about points on this surface and matching terms at each order in the ansätz. This process, described in Appendix A, yields boundary conditions on a static geometry for the Brinkman amplitudes

(2.15a,b)\begin{equation} \hat{\boldsymbol{u}}(\boldsymbol{x})|_{\boldsymbol{x} \in \overline{{\partial}\varOmega}_k} = \hat{U}_k \boldsymbol{e}_{z},\quad \lim_{|\boldsymbol{x}|\to \infty}\hat{\boldsymbol{u}}(\boldsymbol{x}) = \boldsymbol{0} \end{equation}

and steady flow

(2.16a,b)\begin{equation} \bar{\boldsymbol{u}}(\boldsymbol{x}) |_{\boldsymbol{x} \in \overline{{\partial}\varOmega}_k} = \bar{u}_s \boldsymbol{e}_{\theta} := \frac{{\rm i}\hat{U}_k}{2}\frac{{\partial}\hat{\boldsymbol{u}}^*}{{\partial} z},\quad \lim_{|\boldsymbol{x}| \to \infty} \bar{\boldsymbol{u}}(\boldsymbol{x}) ={-}\bar{U} \boldsymbol{e}_{z}. \end{equation}

We have introduced the steady tangential slip velocity $\bar {u}_s = \bar {u}_s(\theta )$ defined on $\overline { {\partial }\varOmega }_{k}$, where $\theta$ is the polar angle from the $z$ axis in spherical coordinates originating at the centre of sphere $k$. We emphasise that this slip velocity is unrelated to the Brinkman boundary layer described above, instead arising as in the ‘swimming sheet’ of Taylor (Reference Taylor1951) from periodic motion of the boundary.

Writing $\hat {\boldsymbol {u}} = \hat {u}_r \boldsymbol {e}_{r} + \hat{u}_\theta \boldsymbol {e}_{\theta }$ and $\bar {\boldsymbol {u}} = \bar {u}_r \boldsymbol {e}_{r} + \bar {u}_\theta \boldsymbol {e}_{\theta }$ (letting $r,\theta$ and their associated unit vectors refer to the spherical coordinates at sphere $k$), we briefly show that $ {\partial }_z \hat {\boldsymbol {u}} \propto \boldsymbol {e}_{\theta }$ on sphere $k$ as claimed. Since $ {\partial }_\theta \hat {\boldsymbol {u}} = 0$ on the sphere surface, $ {\partial }_z\hat {\boldsymbol {u}} = \cos \theta {\partial }_r \hat {\boldsymbol {u}}$, implying $\bar {\boldsymbol {u}}_r \propto {\partial }_r \hat {u}_r$. The divergence-free condition requires $ {\partial }_r \hat {u}_r = -(2 \hat{u}_r + \cot \theta \hat{u}_\theta + {\partial }_\theta \hat{u}_\theta )/r$. Substituting in $\hat{u}_r = \hat {U}_k \cos \theta$ and $\hat{u}_\theta = -\hat {U}_k \sin \theta$ shows this quantity vanishes and $\bar {u}_r = 0$ on the surface. Thus, the order-$\varepsilon ^2$ steady flow field $\bar {\boldsymbol {u}}$ obeys the physical requirement of no flux through sphere surfaces.

The swim speed $\bar {U}$ is unknown, as are the individual sphere oscillation amplitudes $\hat {U}_1$ and $\hat {U}_2$. We introduce

(2.17a,b)\begin{equation} \hat{F}_k^{(f)} := \int_{\overline{{\partial}\varOmega}_k} \boldsymbol{e}_{z} \boldsymbol{\cdot} \hat{{\boldsymbol{\mathsf{T}}}} \boldsymbol{\cdot} \boldsymbol{e}_{n} \,{\rm d} S,\quad \bar{F}_k^{(f)} := \int_{\overline{{\partial}\varOmega}_k} \boldsymbol{e}_{z} \boldsymbol{\cdot} \bar{{\boldsymbol{\mathsf{T}}}} \boldsymbol{\cdot} \boldsymbol{e}_{n} \,{\rm d} S, \end{equation}

where $\hat {{\boldsymbol{\mathsf{T}}}}$ and $\bar {{\boldsymbol{\mathsf{T}}}}$ are the Cauchy stress tensors corresponding to the oscillatory and steady flows $(\hat {\boldsymbol {u}},\hat p)$ and $(\bar {\boldsymbol {u}},\bar {p})$, respectively. The three unknowns are thus fixed by the constraints

(2.18a–c)\begin{equation} \hat{U}_2 - \hat{U}_1 = 1,\quad \hat{F}_{i,1} + \hat{F}_{i,2} = 0,\quad \bar{F}_{f,1} + \bar{F}_{f,2} = 0. \end{equation}

See Appendix A for a complete derivation of the boundary conditions and force conditions.

Due to linearity, the steady solution $\bar {\boldsymbol {u}} = \bar {\boldsymbol {u}}_b + \bar {\boldsymbol {u}}_r$ can be decomposed into flows driven exclusively by the boundary condition ($\bar {\boldsymbol {u}}_b,p_b$) and Reynolds stress ($\bar {\boldsymbol {u}}_r,p_r$), respectively, such that

(2.19a–c)\begin{gather} \nabla^2 \bar{\boldsymbol{u}}_b = \boldsymbol{\nabla} \bar{p}_b,\quad \bar{\boldsymbol{u}}_b(\boldsymbol{x})|_{\boldsymbol{x} \in \overline{{\partial}\varOmega}_k} = u_s \boldsymbol{e}_{\theta},\quad \lim_{|\boldsymbol{x}| \to \infty} \bar{\boldsymbol{u}}_b ={-}\bar{U}_b \boldsymbol{e}_{z}, \end{gather}

(2.20a–c)\begin{gather}\nabla^2 \bar{\boldsymbol{u}}_r = \boldsymbol{\nabla} \bar{p}_r - \boldsymbol{\nabla} \boldsymbol{\cdot} \bar{{\boldsymbol{\mathsf{R}}}},\quad \bar{\boldsymbol{u}}_r(\boldsymbol{x}) |_{\boldsymbol{x} \in \overline{{\partial}\varOmega_k}} = \boldsymbol{0},\quad \lim_{|\boldsymbol{x}| \to \infty} \bar{\boldsymbol{u}}_r ={-}\bar{U}_r \boldsymbol{e}_{z}, \end{gather}

where the mechanism-specific swim speeds $\bar {U}_b$ and $\bar {U}_r$ satisfy $\bar {U} = \bar {U}_b + \bar {U}_r$. In particular, by applying the reciprocal theorem as in Appendix B, the two speeds can be written as

(2.21a,b)\begin{equation} \bar{U}_b = \sum_{k=1}^2 \frac{{\rm i} \hat{U}_k}{2F'} \int_{\overline{{\partial}\varOmega}_k} \boldsymbol{e}_{n} \boldsymbol{\cdot} {\boldsymbol{\mathsf{T}}}^\prime \boldsymbol{\cdot} \frac{{\partial} \hat{\boldsymbol{u}}^*}{{\partial} z} {\rm d} S, \quad \bar{U}_r = \frac{M^2}{2 F'} \int_\varOmega \hat{\boldsymbol{u}} \boldsymbol{\cdot} {\boldsymbol{\mathsf{E}}}^\prime \boldsymbol{\cdot} \hat{\boldsymbol{u}}^* \,{\rm d} V, \end{equation}

where the primed variables correspond to the Stokes flow resulting from towing the dimer at a velocity $U^\prime$ with a force $F^\prime$. The full set of model equations (2.13a,b)–(2.18a–c) may be solved numerically using the finite element method and analytically with the method of reflections. In the next section we use the finite element approach to examine the swim speed as a function of $M^2$ or, when considering a dense swimmer in Stokes flow as in (Hubert et al. Reference Hubert, Trosman, Collard, Sukhov, Harting, Vandewalle and Smith2021), a function of $S^2$. In § 4 we interpret these results analytically.

3.1. Finite element method

The dimensionless Brinkman equations, given a parameter $\alpha$ representing an inverse screening length, are written as

(3.1a,b)\begin{equation} \left(\nabla^2 - \alpha^2\right) \boldsymbol{u} = \boldsymbol{\nabla} p - \boldsymbol{f},\quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0. \end{equation}

They describe flow subject to a body force $\boldsymbol {f} \in L^2(\varOmega )$ in a porous region where a drag force $-\alpha ^2 \boldsymbol {u}$ is of the same magnitude as viscous stresses. If $\alpha = 0$, they are the Stokes equations. Therefore, we develop our numerical solution procedure for arbitrary complex $\alpha$.

If $\boldsymbol {u}$ has Dirichlet boundary conditions, it is well known (Iliev, Lazarov & Willems Reference Iliev, Lazarov and Willems2011) that (3.1a,b) has a unique weak solution $(\boldsymbol {u}, p) \in H^1(\varOmega )^2 \times L_0^2(\varOmega )$ such that, for all test functions $(\boldsymbol {v},q) \in H^1(\varOmega )^2 \times L_0^2(\varOmega )$,

(3.2a,b)\begin{align} \int_\varOmega \left[ \left(\boldsymbol{\nabla} \boldsymbol{v} \boldsymbol{:} \boldsymbol{\nabla} \boldsymbol{u}\right) + \alpha^2 \left(\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{u}\right) - (\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v}) p\right] {\rm d} V = \int_{\varOmega} \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{f} \,{\rm d} V,\quad \int_{\varOmega} q(\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}) \,{\rm d} V = 0. \end{align}

At this point we discretise our axisymmetric domain into two $n_e \times n_e$ grids of curvilinear, quadrilateral elements using the bispherical coordinates $(\xi,\eta )$. The surface of sphere $k$ is given by $\xi = \xi _k$, with $\xi _1 < 0 < \xi _2$, and by symmetry the system is agnostic to the substitution $\eta \to 2{\rm \pi} - \eta$. The grids are periodic over $\eta \in [0,2{\rm \pi} )$, spanning $\xi \in [\xi _1,0]$ and $\xi \in [0,\xi _2]$, respectively. We denote the collection of elements $\mathcal {T}_{n_e}$. On each element, we consider functions defined within the space $Q_m = \text {span}\{\xi ^j \eta ^k : 0 \le j,k \le m\}$, and globally we consider piecewise combinations of these. In particular, we let $\mathcal {C}^1_m(\mathcal {T}_{n_e})$ be the space of such functions that are globally continuous. We seek approximate solutions $\smash {(\boldsymbol {u}^h,p^h) \in C^1_{n_p}(\mathcal {T}_{n_e})^2 \times C^1_{n_p-1}(\mathcal {T}_{n_e})}$, i.e. defined on the so-called generalized Taylor–Hood $\smash {(Q_{n_p}-Q_{n_p-1})}$ element. In terms of basis functions for vector and scalar fields, $\smash {\boldsymbol {\Psi }_k(\boldsymbol {x}) \in C^1_{n_p}(\mathcal {T}_{n_e})^2}$ and $\smash {\psi _k(\boldsymbol {x}) \in C^1_{n_p-1}(\mathcal {T}_{n_e})}$, we write

(3.3a,b)\begin{equation} \boldsymbol{u}^h = \sum_j^{N_u} u_j^h \boldsymbol{\Psi}_j,\quad p^h = \sum_j^{N_p} p_j^h \psi_j, \end{equation}

where $N_u$ and $N_p$ are the number of velocity and pressure degrees of freedom, respectively. Considering test functions $(\boldsymbol {v}^h,q^h) \in C^1_{n_p}(\mathcal {T}_{n_e})^2\times C^1_{n_p-1}(\mathcal {T}_{n_e})$ and substituting into (3.2a,b), the coefficients $u_j^h$ and $p_j^h$ are known to satisfy the linear system

(3.4)\begin{equation} \left(\begin{array}{cc} \boldsymbol{\mathsf{A}}(\alpha^2) & -\boldsymbol{\mathsf{B}}^{\boldsymbol{\mathsf{T}}} \\ \boldsymbol{\mathsf{B}} & 0 \end{array}\right) \left(\begin{array}{c} u^h \\ p^h \end{array}\right) = \left(\begin{array}{c} F \\ 0 \end{array}\right), \end{equation}

where the block matrix elements are given by

(3.5a,b)\begin{equation} {\mathsf{A}}_{ij}(\alpha^2) = \int_\varOmega \left(\boldsymbol{\nabla} \boldsymbol{\Psi}_i \boldsymbol{:} \boldsymbol{\nabla} \boldsymbol{\Psi}_j + \alpha^2 \boldsymbol{\Psi}_i \boldsymbol{\cdot} \boldsymbol{\Psi}_j\right) {\rm d} V,\quad {\mathsf{B}}_{ij} = \int_\varOmega \psi_i \left(\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\Psi}_j\right) {\rm d} V, \end{equation}

and the source term on the right-hand side is

(3.6)\begin{equation} F_i = \int_\varOmega \boldsymbol{\Psi}_i \boldsymbol{\cdot} \boldsymbol{f} \,{\rm d} V. \end{equation}

In bispherical coordinates the velocity boundary conditions can be written as

(3.7a,b)\begin{equation} \boldsymbol{u}(\xi_k,\eta) = V_k \boldsymbol{e}_{z} + v_s(\eta)\boldsymbol{e}_{\theta},\quad \boldsymbol{u}(0,0) ={-}V, \end{equation}

for a set of sphere velocities $V_k$, a boundary slip velocity $v_s$ and the velocity $V$ of the swimmer frame relative to the lab frame. We use $V$ and $v$ to distinguish numerical parameters from system parameters. In all cases below, we directly solve (3.4) using the distributed-memory version of the SuperLU library (Li Reference Li2005) and PETSc library for scientific computation (Balay et al. Reference Balay, Gropp, McInnes and Smith1997, Reference Balay2021a,Reference Balayb), working in C++. The code used to generate the data presented is publicly available at https://github.com/nderr/reciprocal_swimming. The solution procedure is as follows.

3.2. Solution procedures

3.2.1. Brinkman equations

The Brinkman equations (2.13a,b) correspond to the linear system (3.4) with $\smash {\alpha ^2 = iM^2}$, $\smash {\boldsymbol {f} = \boldsymbol {0}}$, $\smash {V_k = \hat {U}_k}$, $\smash {v_s = 0}$ and $V=0$. The unknown numerical solution $(\boldsymbol {u}^h,p^h)$ are complex fields, stemming from the complex symmetric block matrix $A$. Recall the sphere amplitudes $\hat {U}_k$ are unknown and fixed by the constraints (2.18a–c). To find them, we set $\smash {V_k = \tilde {U}_k}$, trial values of the sphere velocity amplitudes, which we consider as a vector in $\mathbb {C}^2$. These yield a corresponding set of force amplitudes $\smash {\tilde F_k^{(i)}\in \mathbb {C}^2}$ upon solution of (3.4).

Since the equations are linear, we have $\smash {\tilde {F}_{k}^{(i)} = {\mathsf{C}}_{kj} \tilde {U}_j}$ for some constant complex matrix $\smash {\boldsymbol{\mathsf{C}} \in \mathbb {C}^{2 \times 2}}$, with $\smash {{\mathsf{C}}_{kj} = {\rm d}\tilde {F}_{k}^{(i)}/{\rm d}\hat {U}_j}$. In our approach, we calculate $\smash {\tilde {F}_{k}^{(i)}(\tilde {U}_j)}$ at $\smash {(\tilde {U}_1,\tilde {U}_2) = (1,0),\ (0,1),}$ and $\smash {(0,0)}$, and we use the results to compute $\smash {{\rm d}\tilde {F}_{k}^{(i)}/{\rm d}\tilde {U}_j}$ at $\smash {(0,0)}$ via finite difference. The constraints (2.18a–c) require $\smash {\hat {F}_{1}^{(i)} = -\hat {F}_{2}^{(i)}}$ and $\smash {\hat {U}_2 - \hat {U}_1 = 1}$. From the problem definition, $\smash {\tilde {F}_{k}^{(i)} = 0}$ at $\smash {(\tilde {U}_1,\tilde {U}_2)=0}$, so we set

(3.8a–c)\begin{equation} \hat{U}_k = \frac{\hat{U}_{k,*}}{\hat{U}_{2,*}-\hat{U}_{1,*}},\quad \hat{U}_{k,*} = C^{{-}1}_{kj} \tilde{F}_{j,*}^{(i)},\quad \tilde{F}_{2,*}^{(i)} ={-}\tilde{F}_{1,*}^{(i)} = 1, \end{equation}

for an intermediate set of force and velocity amplitudes $\smash {\tilde {F}_{k,*}^{(i)}}$ and $\smash {\hat {U}_{k,*}}$. We solve (3.4) once more with this set of velocity amplitudes to obtain a solution that satisfies (2.18a–c).

4.1. Describing flow about two spheres

The motion of two spheres in a Brinkman theorem, especially in an axisymmetric configuration, has been the subject of much study. Kim & Russel (Reference Kim and Russel1985) developed a set of Faxén laws relating the force and moment on a sphere in a Brinkman medium to the background flow at its location. This was accomplished using the well-known method of reflections from low-Reynolds number flow theory, applicable in Brinkman media because of the linearity of the Brinkman equations. An initial approximation to the solution about two spheres is established by considering the Brinkman solution to flow about each sphere in unbounded fluid. The presence of each sphere's flow does not vanish on the other sphere's surface. Consequently, correction flows must be added to correct this violation of the boundary conditions. In the limit of infinite reflections, the series of flows converges to the solution. In practice, the sum can be truncated at a desired precision.

Kim & Russel (Reference Kim and Russel1985) define these correction flows in integral form, a convenient choice for calculating the force and torque on the spheres. They note that the process can be completed using an explicit velocity representation using addition theorems to describe flows about one sphere around the other using spherical harmonics. Recently, Liu & Bhattacharya (Reference Liu and Bhattacharya2020) have employed this approach to devise a solution for Brinkman flow in the presence of two spheres undergoing general relative motion. Here, we have developed a more limited procedure restricted to axisymmetric motion in the absence of torsion, significantly reducing the space of basis functions.

4.2. General solution about one sphere

The general axisymmetric solution to the Brinkman equations can be written in terms of two scalar functions $\chi =\chi (\boldsymbol {x})$ and $\phi =\phi (\boldsymbol {x})$, where the flow fields

(4.1a,b)\begin{equation} \boldsymbol{u} = \boldsymbol{\nabla} \phi + \boldsymbol{\nabla}\times\boldsymbol{\nabla}\times\left(\boldsymbol{x} \chi\right),\quad p ={-}\alpha^2 \mu \phi \end{equation}

are a solution to the Brinkman equations if $\phi$ is harmonic and $\chi$ satisfies the Helmholtz equation,

(4.2)\begin{equation} \left(\nabla^2 - \alpha^2\right)\chi=0. \end{equation}

The presence of the Helmholtz equation significantly alters the structure of the flow and the method of solution as compared with the Stokes equations. First, $\chi = \chi (\alpha r,\theta )$ based on dimensional considerations, suggesting the aforementioned boundary layer effects manifest purely through the function $\chi$. Second, the Helmholtz equation is not separable in bispherical coordinates, so exact solutions in this coordinate system are not possible, unlike in the case of the Stokes equations.

General solutions $f=f(\boldsymbol {x})$ to each of the scalar functions can be written in terms of the Laplace spherical harmonics $Y_l(\theta )$ defined in Appendix C such that

(4.3)\begin{equation} f(\boldsymbol{x}) = \sum_{l=0}^\infty f_l(r) Y_l(\theta), \end{equation}

where the $f_l$ are a series of radially varying functions. Solutions to the Laplace equation are

(4.4a,b)\begin{equation} \phi_l^-= \frac{1}{r^{l+1}},\quad \phi_l^+= r^l, \end{equation}

giving rise to the so-called solid harmonics. The Helmholtz solutions to (4.2) are

(4.5a,b)\begin{equation} \chi^{-}_l = k_l(\alpha r),\quad \chi^+_l = g_l(\alpha r), \end{equation}

where $g_l$ and $k_l$ are the modified spherical Bessel functions of the first and second kind, respectively. The superscripts $+$ and $-$ denote solutions that are regular and singular as $r\to 0$, respectively. These four functions provide four basis functions for Brinkman flows at the $l$th mode of spherical harmonics,

(4.6a–d)\begin{equation} \left.\begin{gathered} \boldsymbol{u}_{l0} = \boldsymbol{\nabla} (\phi_l^- Y_l),\quad \boldsymbol{u}_{l1} = \boldsymbol{\nabla} \times \boldsymbol{\nabla} \times (\boldsymbol{x} \chi_l^- Y_l),\\ \boldsymbol{u}_{l2} = \boldsymbol{\nabla} (\phi_l^+ Y_l),\quad \boldsymbol{u}_{l3} = \boldsymbol{\nabla} \times \boldsymbol{\nabla} \times (\boldsymbol{x} \chi_l^+ Y_l). \end{gathered}\right\} \end{equation}

In this way the scalar fields may be discretised and represented by a series of coefficients. The solution can be obtained to desired accuracy by truncating to $n_l$ terms in the spherical harmonic expansion, representing the scalar field via a vector in $\mathbb {R}^{4n_l}$.

We will denote a set of coefficients defined about sphere $k$ with a superscript label, $c^{(k)}$, describing a flow $\boldsymbol {u}^{(k)} = \boldsymbol {u}^{(-k)} + \boldsymbol {u}^{(+k)}$ such that

(4.7a,b)\begin{equation} \boldsymbol{u}^{({-}k)} = \sum_{l=1}^{n_l} \left(c_{l0}^{(k)} \boldsymbol{u}_{l0} + c_{l1}^{(k)} \boldsymbol{u}_{l1}\right),\quad \boldsymbol{u}^{({+}k)} = \sum_{l=1}^{n_l} \left(c_{l2}^{(k)} \boldsymbol{u}_{l2} + c_{l3}^{(k)} \boldsymbol{u}_{l3}\right). \end{equation}

We similarly denote $c^{(+k)}_{lj}$ and $c^{(-k)}_{lj}$ as the regular and singular coefficients of $c^{(k)}_{lj}$, such that, for the former, $j$ takes the values 0 and 1 and, for the latter, 2 and 3. Note the sum starts from $l=1$ because the $l=0$ flow modes violate the incompressibility condition. For flows about an isolated sphere in unbounded flow, $\boldsymbol {u}^{(+k)}= \boldsymbol {0}$, so the flow vanishes as $r \to \infty$. Note as well that $\boldsymbol {u}_{l1}$ only varies within the boundary layer, since it depends only on $\chi _l(\alpha r)$. Now, we construct the vector of coefficients $c^{(k)}_{lj}$ about sphere $k$ in the presence of another sphere by using the method of reflections.

4.3. Solution procedure about two spheres

As shown by figure 5, we begin by considering a superposition of flows describing motion of each sphere in unbounded flow without regard to the presence of the other. We will iteratively introduce additional flows (‘reflections’) to correct the violation of boundary conditions stemming from existing terms. Let the coefficients $c^{(-k,j)}$ describe the singular coefficients of the flow about sphere $k$ at the $j$th step in the reflection. We let $c^{(k,0)}$ be the ‘zeroth’ reflection, i.e. the flow about an isolated sphere, such that $c^{(+k,0)}_{lm} = 0$. We let $c^{(+k,1)}$ denote the regular modes describing the flow incident at sphere $k$ originating at sphere $m \ne k$. By linearity, there exists an expression

(4.8)\begin{equation} c^{({+}k,j)}_{ln} = {\mathsf{A}}^{mk}_{lnop} c^{({-}m,j-1)}_{op}, \end{equation}

relating the two sets of coefficients via some interaction tensor $\boldsymbol{\mathsf{A}}$. Similarly, each singular reflection flow must be linearly related to the incident flow it corrects, i.e.

(4.9)\begin{equation} c^{({-}k,j)}_{ln} = {\mathsf{B}}^{k}_{lnop} c^{({+}k,j)}_{op} \end{equation}

for some induction tensor $\boldsymbol{\mathsf{B}}$. The functional forms of $\boldsymbol{\mathsf{A}}$ and $\boldsymbol{\mathsf{B}}$ are reported in Appendix C, and together they comprise a total reflection tensor $\boldsymbol{\mathsf{R}}$,

(4.10a,b)\begin{equation} c_{ln}^{({-}k,j)} = {\mathsf{R}}^{mk}_{lnop} c^{({-}m,j-1)}_{op},\quad {\mathsf{R}}^{mk}_{lnop} = {\mathsf{B}}^k_{lnqs} {\mathsf{A}}^{mk}_{pqop}. \end{equation}

Note the coefficients $c^{(k,0)}$ must be linearly related to the velocity $\hat {U}_k$, but that $c^{(k,1)}$ will be proportional to $\hat {U}_m$ as demonstrated in figure 5. In general, $c^{(k,j)}$ contains a factor $\hat {U}_k$ if $j$ is even and $\hat {U}_m$ if $j$ is odd. Introducing the notation

(4.11)\begin{equation} \hat{U}^{(k,j)} = \begin{cases} \hat{U}_k & \text{if } j\text{ mod }2 = 0, \\ \hat{U}_m & \text{otherwise}, \end{cases} \end{equation}

we can rewrite the coefficients above as $c^{(k,j)} \to \tilde {c}^{(k,j)}$ and let $\tilde {c}^{(k,j)} = \hat {U}^{(j,k)} c^{(j,k)}$ for a set of normalised coefficients $c^{(j,k)}$ that are agnostic to the sphere velocity magnitudes. The unknown amplitudes $\hat {U}_1$ and $\hat {U}_2$ can be determined by noting, as before, that $\hat {F}_j^{(i)} = {\mathsf{C}}_{jk} \hat {U}_k$ for a matrix $\boldsymbol{\mathsf{C}}$. Given $n_r$ reflections, the components of $\boldsymbol{\mathsf{C}}$ may be calculated according to

(4.12)\begin{equation} {\mathsf{C}}_{jk} = \begin{cases} {\rm i}M^2 \left(\dfrac{4{\rm \pi} a_k^3}{3} \right) - \displaystyle\sum_{m=0}^{\lfloor n_r/2 \rfloor} c_{1n}^{(j,2m)} F_{nj}, & j=k, \\ \displaystyle-\sum_{m=0}^{\lceil n_r/2 \rceil - 1} c_{1n}^{(j,2m+1)} F_{nj}, & j \ne k, \end{cases} \end{equation}

where $F_{nk}$ is the force applied to sphere $k$ by the $n$th $l=1$ basis flow harmonic $\boldsymbol {u}_{1n}$,

(4.13)\begin{equation} F_{nj} = \int_{\overline{{\partial}\varOmega}_j} \boldsymbol{t}_{1n} \boldsymbol{\cdot} \boldsymbol{e}_{z} \,{\rm d} S, \end{equation}

and $\boldsymbol {t}_{ln}$ is the traction associated with $\boldsymbol {u}_{ln}$. Then $\hat {U}_k$ can be calculated just as in (3.8a–c).

Figure 5.

Schematic of the method of reflections for calculating the Brinkman flow field around two spheres undergoing axisymmetric relative motion. A set of coefficients $c^{(-k,0)}$ represents the exact solution of the flow about sphere $k$ in unbounded fluid, but the presence of this flow violates the boundary conditions at the other sphere $m$. This incident flow can be described in terms of harmonics about sphere $m$ by a set of coefficients $c^{(+m,1)}$, and the spheres’ resistance to deformation creates a reflection flow $c^{(-m,1)}$ that exactly cancels out the incident flow at sphere $m$. However, this induces a new incident flow at sphere $k$ described by $c^{(+k,2)}$. This process can be repeated and truncated after a specified number of iterations. Each set of coefficients are linearly related to those they are induced by. Here, all the red coefficients are proportional to $\hat {U}_2$, and all the blue to $\hat {U}_1$.

Note the swim speeds (2.21a,b) as determined by the reciprocal theorem also depend on the Stokes flow $\boldsymbol {u}'$ corresponding to a static dimer being towed through the medium. This flow can also be represented in terms of a harmonic expansion though the method of reflections, and the process of doing so is similar to the above. Details can be found in, for example, Happel & Brenner (Reference Happel and Brenner2012).

4.4. Asymptotic contributions to swim speed

Using this approach, we construct asymptotic approximations to the swim speed $\bar {U}$ and the contributions $\bar {U}_b + \bar {U}_r = \bar {U}$ related to the steady boundary velocity and Reynolds stress. Because we are interested in the leading-order description of the behaviour, we consider the first two modes’ $(n_l = 2)$ interactions through a single reflection $(n_r = 1)$. In the description of the swim speeds via the reciprocal theorem (2.21a,b), each contribution depends quadratically on the Brinkman amplitude fields $\hat {\boldsymbol {u}}$ and $\hat {\boldsymbol {u}}^*$. Since each term in the harmonic expansion of $\hat {\boldsymbol {u}}$ is proportional to $\hat {U}_1$ or $\hat {U}_2$, each term in the expansions of $\bar {U}_b$ and $\bar {U}_r$ will be proportional to $\hat {U}_1 \hat {U}_2^*$, $|\hat {U}_1|^2$ or $|\hat {U}_2|^2$. Note that since $\hat {U}_1 = |\hat {U}_1|\exp ({\textrm {i}({\rm \pi} +\phi _1)})$ and $\hat {U}_2 =|\hat {U}_2| \exp ({-\textrm {i} \phi _2})$ as in (2.4) and figure 1(c), then

(4.14a,b)\begin{equation} \hat{U}_1 \hat{U}_2^* ={-}|\hat{U}_1| |\hat{U}_2| \,{\rm e}^{{\rm i}\phi},\quad \phi = \phi_1 + \phi_2, \end{equation}

where $\phi$ is the phase difference between the spheres’ extrema. This provides a method of categorising the physical origin of each contribution to the swim speed in terms of flow- or geometry-driven interactions, illustrated in figure 6.

Figure 6.

Schematic of the types of sphere interactions and their relation to products of the complex Fourier amplitudes $\smash {\hat {U}_1}$ and $\smash {\hat {U}_2}$. Sphere interactions can be decomposed into geometric and flow contributions. The presence of the other sphere creates a time-invariant geometric asymmetry akin to an oscillating sphere near a wall (a). Each sphere is also moving within an oscillating background flow generated by the motion of the other (b). By linearity, we can associate products of the velocity amplitudes $\smash {\hat {U}_j \hat {U}_k^*}$ with the different interaction types. The geometric interactions correspond to products of each sphere's amplitude with itself, and the flow interactions to the product of both spheres’ amplitudes.

Contributions to the speeds proportional to $\hat {U}_1 \hat {U}_2^*$ correspond to interactions between one sphere's motion and the oscillating background flow induced by the other. In the presence of a non-zero phase difference $\phi \ne 0$, this interaction is not time reversible, and, thus, can give rise to motion even in the absence of fluid inertia ($M^2 = 0$), since the scallop theorem does not apply. We refer to these terms in the swim speeds as ‘flow interactions’. On the other hand, contributions proportional to $|\hat {U}_1|^2$ and $|\hat {U}_2|^2$ arise from the effects on one sphere's motion caused by the other's presence, which breaks geometric symmetry. We refer to these as ‘geometric interactions’.

4.4.1. Velocity magnitudes and phase difference

Recall from (4.12) that we can relate the fluid force amplitude $\smash {\hat {F}_k^{(f)}}$ to the two-sphere amplitudes $\smash {\hat {U}_k}$ via a mobility matrix defined by the harmonic expansion coefficients. This is sufficient to find the asymptotic dependence of the velocity amplitudes on the fluid inertia $M^2$ (reported in table 2, columns two and three) or, for the Stokes case, the solid inertia $S^2$ (reported in table 2, columns four and five.) In the Stokes limit, we assume the dimensionless ratio $\rho _2/\rho$ becomes large such that $S^2 = 2M^2\rho _2/3\rho$ takes a finite value in the limit $M^2 \to 0$. We introduce a rescaling of the actual mass of sphere $k$,

(4.15)\begin{equation} \tilde{m}_k = a_k^3 \frac{\rho_k}{\rho_2} \end{equation}

such that $2 M^2 m_k = 3 S^2 \tilde{m}_k$. The Stokes scalings in table 2 are reported in terms of this value.

Table 2.

Leading-order scalings of the velocity amplitude magnitudes $|\hat {U}_k|$ and trigonometric functions of the sphere phase difference $\phi$. The second and third columns correspond to the small- and large-inertia limits in the presence of non-zero fluid inertia $(M^2 \ne 0)$, and the fourth and fifth columns are the same limits for Stokes flow $(M^2 = 0)$. For the $M^2 \ne 0$ case, values are given in terms of the effective mass $m_k$ (3.9). For Stokes flow, rescalings of the actual sphere masses $\tilde {m}_k$ (4.15) are used.

In both cases, the velocity magnitudes undergo a transition from depending on the sphere sizes at low-inertia levels to sphere masses at high-inertia levels, and the phase difference $\phi$ vanishes in both limits as the sphere–sphere interactions become dominated by Stokes drag or added mass effects, respectively. However, the rate at which the phase difference changes is not the same, as observed in the functional form of $\sin \phi$. In the fluid inertial case, we observe that at low inertia, the smaller sphere leads the larger sphere since $\phi > 0$ for $a_2 > a_1$ and $\sin \phi \propto (d M)[(a_2-a_1)M]$. This factor of $d$ in the proportionality relation stems from the reduction of the Brinkman boundary layer width as $M$ increases. At high inertia, there is a more complicated comparison for determining the leading sphere which involves both size and mass asymmetries, as $\phi > 0$ if $a_1^2 m_2 > a_2^2 m_1$. Due to the presence of the Basset force that is proportional to $M$, the phase difference falls off as $1/M$. In the Stokes case, there is no boundary layer evolution and no Basset force. The relevant asymmetry is the same in both low- and high-inertial limits, as $\phi > 0$ if $a_1 \tilde {m}_2 > a_2 \tilde {m}_1$, and the dependence on $S$ in the proportionalities is always quadratic since the added mass effect is the only inertial coupling between the spheres.

Below, we report the leading-order contributions for each combination of inertial limit ($M^2 \ll 1$ vs $M^2 \gg 1$), mechanism (boundary velocity versus Reynolds stress), and interaction type (flow versus geometric). See Appendix E for the complete derivation. For ease of notation, we introduce the quantities

(4.16)\begin{equation} \left.\begin{gathered} a_{\smash{{+}}} = a_1 + a_2,\quad a_{\smash{{\hspace{1pt}-}}} = \frac{1}{a_1^{{-}1} + a_2^{{-}1}},\quad \left[a\right] = \frac{a_2 - a_1}{a_1 + a_2}, \\ m_{\smash{{+}}} = m_1 + m_2,\quad m_{\smash{{\hspace{1pt}-}}} = \frac{1}{m_1^{{-}1} + m_2^{{-}1}},\quad \left[m\right] = \frac{m_2 - m_1}{m_1 + m_2}. \end{gathered}\right\} \end{equation}

The quantities $a_{\smash {{+}}}$ and $a_{\smash {{\hspace {1pt}-}}}$ approach the large and small sphere sizes as the difference between them becomes large, and $\left [a\right ]$ is a dimensionless measure of the size asymmetry. The $m$ symbols are the same, with respect to the effective sphere masses $m_k$.

4.4.2. Stokes flow

In the Stokes limit ($M^2 = 0$), the leading-order boundary velocity contribution is

(4.17)\begin{equation} \bar{U}_b ={-}\frac{3a_{\smash{{\hspace{1pt}-}}}}{2d^2} |\hat{U}_1| |\hat{U}_2| \sin \phi, \end{equation}

consistent with previous investigations of inertial dimers in Stokes flow (Hubert et al. Reference Hubert, Trosman, Collard, Sukhov, Harting, Vandewalle and Smith2021). Note that $\phi \ne 0$ is required for motion according to the scallop theorem, that this symmetry is broken by the presence of solid inertia ($S^2 > 0$) and that the dimer swims towards the leading sphere (i.e. towards sphere 1 when $\phi > 0$.) As mentioned, the sign of $\sin \phi$ and, thus, $\bar {U}_b$ is the same in both limits, since it depends in both cases on the quantity $a_1 \tilde {m}_2 - a_2 \tilde {m}_1.$ By definition (2.21a,b) there is no Reynolds stress contribution in the Stokes case, so $\bar {U} = \bar {U}_b$. Substituting in the values for $|\hat {U}_k|$ and $\phi$ in table 2 yields

(4.18)\begin{equation} \bar{U} \approx \begin{cases} \dfrac{a_{\smash{{\hspace{1pt}-}}}}{3 a_{\smash{{+}}}^2 d^2}\left(a_2 \tilde{m}_1 - a_1 \tilde{m}_2\right) S^2, & 0 < S^2 \ll 1 \\ \dfrac{27 a_{\smash{{\hspace{1pt}-}}}}{4\tilde{m}_{\smash{{+}}}^2 d^2}\left(a_2 \tilde{m}_1 - a_1 \tilde{m}_2\right) \dfrac{1}{S^2}, & S^2 \gg 1. \end{cases} \end{equation}

In the case of $\rho _1 = \rho _2$, this implies motion towards the small sphere as observed in our numerical results and the work of Felderhof (Reference Felderhof2016) and Hubert et al. (Reference Hubert, Trosman, Collard, Sukhov, Harting, Vandewalle and Smith2021).

4.4.3. Mechanisms at low inertia

For non-zero fluid inertia $0 < M^2 \ll 1$, the leading-order boundary velocity contribution is

(4.19)\begin{align} \bar{U}_b &= \left[-\frac{3a_{\smash{{\hspace{1pt}-}}}}{2 d^2} \sin\phi + \frac{3 a_{\smash{{+}}} a_{\smash{{\hspace{1pt}-}}}}{4\sqrt{2} d^2} \left[a\right] M\cos \phi\right]|\hat{U}_1||\hat{U}_2| \nonumber\\ &\quad + \frac{3 a_{\smash{{\hspace{1pt}-}}}}{4\sqrt{2}} M\left(a_2 |\hat{U}_2|^2 - a_1 |\hat{U}_1|^2\right). \end{align}

It may appear that the Stokesian term will be small compared with the others as $M$ grows. However, substituting the expressions in table 2 for the velocity magnitudes shows the $\cos \phi$ part of the flow interaction cancels the geometric interaction at leading order,

(4.20)\begin{equation} \bar{U}_b ={-}\frac{3a_{\smash{{\hspace{1pt}-}}}^2}{2 a_{\smash{{+}}} d^2} \sin \phi - \frac{3a_{\smash{{\hspace{1pt}-}}}^2}{2\sqrt{2} d^2} \left[a\right] M \sin^2\left(\frac{\phi}{2}\right), \end{equation}

so that, as before, the dimer swims towards the leading sphere. Upon substitution for $\sin \phi$, the first term dominates and

(4.21a,b)\begin{equation} \bar{U}_b = C^{(b)}_{\smash{{\hspace{1pt}-}}} M^2,\quad C^{(b)}_{\smash{{\hspace{1pt}-}}} ={-}\frac{9a_{\smash{{\hspace{1pt}-}}}^2}{16 d} \left[a\right]. \end{equation}

The boundary velocity contribution always promotes motion towards the smaller sphere, proportional to the dimensionless size asymmetry $\left [a\right ]$. Boundary layer effects are manifested in the scaling for $\sin \phi \propto d M$, so that $C^{(b)}_{\smash {{\hspace {1pt}-}}} \propto 1/d$ rather than $1/d^2$ as in the other coefficients arising from the method of reflections in (4.19). Numerical results for $a_1 \in [0.1,0.9]$ and $d \in \{3,5,10\}$ are shown in the bottom row of figure 7, normalised by the proportionality in (4.21a,b). The agreement is good for $Md < 1$, when each sphere lies within the Brinkman boundary layer surrounding the other.

Figure 7.

Normalised swim speeds at low fluid inertia for an equal-density swimmer with $\rho _1 = \rho _2 = 1$. The swim speed contributions from the boundary velocity $\smash {\bar {U}_b}$ (d–f) and Reynolds stress $\smash {\bar {U}_r}$ (a–c) are plotted with respect to $M^2 \le 1$ for a variety of small sphere radii $a_1$ (denoted by line colour) and separation distances $d$ (separated by column.) In each plot, a black vertical line marks where $M d = 1$. For points to the left of this line, the spheres lie within the others’ Brinkman boundary layer. To the right of the line, the boundary layers are small and do not extend to the other sphere. The speeds are scaled by $\smash {C^{(b/r)}_{-}}$, given by (4.21a,b) and (4.23a,b), and lie parallel to the dashed line of slope 1, showing that in this region $\smash {\bar {U}_{b/r} \approx C^{(b)}_- M^2}$ within the boundary layers.

The leading-order Reynolds stress contribution takes the form

(4.22)\begin{equation} \bar{U}_r = \frac{15 a_{\smash{{+}}}^2 a_{\smash{{\hspace{1pt}-}}}}{32 d^2} \left[a\right] M^2 |\hat{U}_1| |\hat{U}_2| \cos \phi. \end{equation}

There is no contribution from geometric interactions here because it can be shown to vanish at order $M^2/d^2$ (see Appendix E.) Substitution of the values in table 2 yields

(4.23a,b)\begin{equation} \bar{U}_r = C^{(r)}_{\smash{{\hspace{1pt}-}}} M^2,\quad C^{(r)}_{\smash{{\hspace{1pt}-}}} = \frac{15 a_{\smash{{\hspace{1pt}-}}}^2 a_{\smash{{+}}}}{32 d^2} \left[a\right] M^2, \end{equation}

so the Reynolds stress contribution always promotes motion towards the larger sphere. Numerical comparisons to the scaling are shown in the top row of figure 7. As in the case of the boundary velocity contribution, there is good agreement with the computed scaling, although there is less variability with respect to both the small sphere radius and separation distance $d$. As in that case, however, the agreement is better for $Md < 1$.

4.4.4. Mechanisms at high inertias

In the large-inertia limit $M^2 \gg 1$, the boundary velocity takes the form

(4.24)\begin{align} \bar{U}_b &= \left[\frac{1}{2\sqrt{2}} \frac{a_{\smash{{+}}}^2 a_{\smash{{\hspace{1pt}-}}}^2 M}{d^4} \left(\left[a\right] \cos \phi + \sin \phi\right) \right]|\hat{U}_1||\hat{U}_2| \nonumber\\ &\quad + \frac{3}{4\sqrt{2}} \frac{a_{\smash{{\hspace{1pt}-}}} M}{d^2} \left(a_2 |\hat{U}_2|^2 - a_1 |\hat{U}_1|^2\right). \end{align}

Since $Md > 1$ in this limit, the spheres do not exist within each other's boundary layers and, as a result, the flow interactions are weaker $(\propto 1/d^4)$ than in the small-inertial limit $(\propto 1/d^2)$. The geometric interactions remain the same strength $(\propto 1/d^2)$ and, thus, dominate at order $M$. Substituting in the values in table 2 gives

(4.25a,b)\begin{equation} \bar{U}_b \approx C_+ M,\quad C_+= \frac{3}{4\sqrt{2}} \frac{a_{\smash{{\hspace{1pt}-}}}}{d^2 m_{\smash{{+}}}^2} \left(a_2 m_1^2 - a_1 m_2^2\right). \end{equation}

Performing the same analysis for the Reynolds stress shows that at order $M$, the flow and geometric effects exactly cancel. However, there exists an additional geometric contribution at order 1 that is not present in the boundary velocity contribution,

(4.26)\begin{align} \bar{U}_r &={-}\left[\frac{1}{2\sqrt{2}} \frac{a_{\smash{{+}}}^2 a_{\smash{{\hspace{1pt}-}}}^2 M}{ d^4} \left(\left[a\right] \cos \phi + \sin \phi\right) \right]|\hat{U}_1||\hat{U}_2| \nonumber\\ &\quad - \frac{3}{4\sqrt{2}} \frac{a_{\smash{{\hspace{1pt}-}}} M}{d^2} \left(a_2 |\hat{U}_2|^2 - a_1 |\hat{U}_1|^2\right) + \frac{9a_{\smash{{\hspace{1pt}-}}}}{8d^2}\left(|\hat{U}_1|^2 - |\hat{U}_2|^2\right). \end{align}

At leading order, then, the isolated Reynolds stress contribution is

(4.27)\begin{equation} \bar{U}_r \approx{-}C_+ M. \end{equation}

Numerical results comparing to the scalings (4.25a,b) and (4.27) are shown in figure 8. In contrast to the low-inertia case, there is little variation with respect to $a_1$ or $d$. The agreement with the scalings is good, especially at higher $M$ where the effects of the constant term in (4.27) become negligible.

Figure 8.

Normalised swim speeds at high fluid inertia for an equal-density swimmer with $\rho _1 = \rho _2 = 1$. The swim speed contributions from the boundary velocity $\smash {\bar {U}_b}$ (d–f) and Reynolds stress $\smash {\bar {U}_r}$ (a–c) are plotted with respect to $M^2 \ge 1$ for a variety of small sphere radii $a_1$ (denoted by line colour) and separation distances $d$ (separated by column.) The speeds are scaled by the large-$M^2$ coefficient $\smash {C_+}$ from (4.25a,b) and (4.27). The lines’ approaching the 1/2 slope of the dashed line indicate that at high levels of inertia $\bar {U}_{b/r} \approx C_+M$. There is little variation with respect to $d$, as the spheres exist outside of each others’ boundary layers for all of the provided inertial range.

4.4.5. Total swim speed

Using the scalings for each mechanism listed in the previous subsection, we can determine the asymptotic form of the total swim speed in the low- and high-inertial limits. At $0 < M^2 \ll 1$, the boundary velocity scaling (4.21a,b), which scales as $1/d$, dominates that of the Reynolds stress (4.23a,b), at $1/d^2$. As such,

(4.28)\begin{equation} \bar{U} \approx{-}\frac{9 a_{\smash{{\hspace{1pt}-}}}^2}{16 d} \left[a\right] M^2, \end{equation}

so that across the entire parameter space we expect motion towards the small sphere at small-but-finite levels of fluid inertia.

At $M^2 \gg 1$, the total swim speed $\bar {U} = \bar {U}_b + \bar {U}_r$ will have no order-$M$ contribution, since the $C_+ M$ terms in (4.25a,b) and (4.27) will cancel. We are left with the order-1 term in (4.26), which gives a theoretical limiting speed as $M^2 \to \infty$,

(4.29)\begin{equation} \bar{U} \approx \frac{9 a_{\smash{{\hspace{1pt}-}}}}{8 d^2} \left[m\right] =: \bar{U}_\infty, \end{equation}

showing that at high inertia motion is always directed towards the more massive sphere. This confirms the impression formed by figure 3(a,b) that at $M^2 \ll 1$ motion is directed towards the smaller sphere and at $M^2 \gg 1$ towards the heavier sphere.

Numerically calculated swim speeds for several dimers with $a_1 = 1/2$ and $d=3$ (the same as in figures 2 and 3) are plotted in figure 9. Each line corresponds to a different effective mass ratio $m_1/m_2$ and, as before, the densities are chosen such that the dimer is neutrally buoyant. In (a) the raw speeds are plotted as a function of inertia, yielding a range of velocities. In (b) the speeds are normalised by the limiting factor $\bar {U}_\infty$, collapsing for $M^2 \gg 1$ onto a single trajectory.

Figure 9.

Plots of the total swim speed $\bar {U}$ as a function of fluid inertia $M^2$. (a) The swim speed is plotted for a system with $a_1 = 1/2$ and $d = 3$ for a variety of effective mass ratios $m_1/m_2$. Positive values correspond to motion towards sphere 2. At high $M^2$, motion is towards the more massive sphere. (b) The normalised swim speed $\smash {\bar {U} / \bar {U}_\infty }$ is plotted over the same range of inertia. Positive values here correspond to motion in the same direction as the theoretical limiting value $\bar {U}_\infty$ (4.29) as $M^2 \to \infty$. The collapse indicates consistent scaling across mass ratios.

For the single-sphere case, Felderhof & Jones (Reference Felderhof and Jones2017) derived analogous asymptotic limits in terms of the multipole coefficients of the swim stroke. There are a few differences worth noting that stem from the dimer geometry. First, for a single sphere, the boundary is always subject to the Brinkman boundary layer. Thus, the introduction of the separation distance $d$ represents a unique aspect of the parameter system. In particular, it yields two separate transitions. The low-$M^2$ scalings in figure 7 are valid only until $M d > 1$ and the spheres exit each others’ boundary layers. However, the high-$M^2$ scalings in figure 8 do not show good agreement until $M > 1$, when the added mass effect dominates the velocity amplitude ratios in table 2. These multiple transitions give rise to the double direction switch observed in figure 3(b).

Secondly, while multipole coefficients are widely useful in a wide variety of contexts including squirmer analysis (Pedley Reference Pedley2016), the dimer geometry provides a clear and intuitive understanding of the results in terms of different types of asymmetries in the problem. We clearly see the transition from a dependence on size asymmetry (due to the domination of Stokes drag) to mass asymmetry (due to the domination of the added mass effect.) One difference from the one-sphere analysis is that Felderhof & Jones (Reference Felderhof and Jones2021) found swimming could proceed in either direction as a function of the surface deformation stroke. However, in this case the gait (corresponding to the choice of $\hat {U}_1$ and $\hat {U}_2$) is initially unknown and must be solved for as a function of the level of fluid and solid inertia.