The Role of the Double-Layer Potential in Regularised Stokeslet Models of Self-Propulsion

Abstract

The method of regularised stokeslets is widely used to model microscale biological propulsion. The method is usually implemented with only the single-layer potential, the double-layer potential being neglected, despite this formulation often not being justified a priori due to nonrigid surface deformation. We describe a meshless approach enabling the inclusion of the double layer which is applied to several Stokes flow problems in which neglect of the double layer is not strictly valid: the drag on a spherical droplet with partial-slip boundary condition, swimming velocity and rate of working of a force-free spherical squirmer, and trajectory, swimmer-generated flow and rate of working of undulatory swimmers of varying slenderness. The resistance problem is solved accurately with modest discretisation on a notebook computer with the inclusion of the double layer ranging from no-slip to free-slip limits; the neglect of the double-layer potential results in up to 24% error, confirming the importance of the double layer in applications such as nanofluidics, in which partial slip may occur. The squirming swimmer problem is also solved for both velocity and rate of working to within a small percent error when the double-layer potential is included, but the error in the rate of working is above 250% when the double layer is neglected. The undulating swimmer problem by contrast produces a very similar value of the velocity and rate of working for both slender and nonslender swimmers, whether or not the double layer is included, which may be due to the deformation’s ‘locally rigid body’ nature, providing empirical evidence that its neglect may be reasonable in many problems of interest. The inclusion of the double layer enables us to confirm robustly that slenderness provides major advantages in efficient motility despite minimal qualitative changes to the flow field and force distribution.

1. Introduction

In our contribution to this Special Issue, we discuss the double-layer potential in the regularised stokeslet boundary integral equation, focussing on the circumstances in which it can be formally eliminated before describing a method for its ‘meshless’ implementation in the manner of the almost ubiquitous single-layer regularised stokeslet method. This approach is then applied to three representative resistance and swimming problems in microscale biological fluid mechanics in which the double-layer potential may be important. In the Introduction below, we give a brief overview of the context of the work before describing the mathematical background to the method of regularised stokeslets and its computational implementation, also discussing the nearest-neighbor approach that is used in this manuscript. In Materials and Methods and supporting Appendices, we describe our implementation of the double-layer potential with a meshless point cloud and its application to the problems considered herein. In Results, we discuss the convergence and efficiency of the method and demonstrate the calculation of flow fields, swimming velocity/trajectory, and rate of working, focussing on the effect of including or neglecting the double-layer term in situations in which this is not formally justified. Finally in the Discussion, the findings are summarised, and potential future applications and areas of further methodological development are considered.

1.1. Literature Review

Microscale biological flow is typically dominated by viscous forces relative to inertia; in the absence of non-Newtonian effects, the fluid dynamics may be approximated by the Stokes flow equations,

$$-\nabla p+\mu {\nabla }^2\mathit{u}=\mathbf{0},\nabla ·\mathit{u}=0,$$

(1)

where $\mathit{u}=\mathit{u}(\mathit{x},t)$ is fluid velocity, $p=p(\mathit{x},t)$ is pressure, and $\mu$ is dynamic viscosity. In systems associated with microscale propulsion and pumping, e.g., by cilia and flagella, these equations are typically accompanied by the no-slip, no-penetration boundary condition $\mathit{u}(\mathit{X},t)=\dot{\mathit{X}}\left(t\right)$, for points $\mathit{X}\left(t\right)$ on the boundary. Due to the Stokes flow equations being linear and having no explicit time dependence, mathematical complexity in microscale biological flow arises from the shape and movement of the boundaries, for example, the movement of the cell body through the fluid, beating of cilia and flagella, and peristaltic contractions of bounding walls.

Prior to 2001, the majority of modeling studies either used local force-drag approximations [1], slender body theory [2,3], perhaps accompanied by image systems [4], or the standard boundary element method [5,6,7]. The latter approach is distinguished by both accuracy and efficiency; however, the learning curve in their application is steep due to both the need to compute singular integrals and requirement to generate a ‘true surface mesh’, i.e., a set of surface points which are assigned to the elements of a surface partitioning. Methods based on volumetric discretisations (e.g., finite difference/element/volume methods) have typically been less popular in the field due to the relatively high computational cost associated with constructing and moving a body-fitted volumetric mesh. A significant exception is the immersed boundary method [8,9,10], which exploits a regular discretisation and moreover enables the solution of the nonlinear equations arising from finite Reynolds number and viscoelasticity—although much of the application of this method in microscale flow has been for 2D flow problems rather than fully 3D.

Since 2001, the Method of Regularised Stokeslets, proposed by R. Cortez [11] and further developed in three dimensions [12,13], has become a popular approach for these systems due to its ease of implementation, particularly in its original ‘Nyström’ form, which discretises the boundary integral equation using a single quadrature rule, then applies collocation at each quadrature point. The method removes the need for a volumetric mesh as well as the requirements to generate a true surface mesh or evaluate weakly singular integrals, as would be needed by the standard boundary element method. In practical terms, it is necessary only to generate a list of points covering the surfaces in the flow, for example, the bodies and flagella of swimming cells, beating cilia, and epithelial surfaces. The problem of computing drag on a body undergoing prescribed motion (resistance problem) is then reduced to inverting a linear system to find a vector of unknown stokeslet strengths (forces) at each discretisation point. The problem of computing translation and rotation of a body due to movements prescribed only in a frame moving with the body (swimming problem) can be expressed similarly as inverting a linear system for the force augmented with additional unknowns for the velocity and angular velocity and additional constraints on the total force and moment [14].

Recent examples of the success of the method of regularised stokeslets implemented via a Nyström discretisation include: elucidating the role of bacteria flagellar polymorphism in bundling/unbundling and hence run-and-tumble behaviour [15], estimation of hydrodynamic forces triggering protist contraction [16], simulating the dynamics of elongated fibres in shear [17,18], modeling of flows due to sperm-templated microrobots [19], and modelling eukaryotic flagella synchronisation [20]. Other key examples over the last two decades include applications to cilia-driven flows [13], flows due to confinement of flagellated bacteria [21], modeling patchy deposition in filtration of bacteria [22], sperm motility [23,24,25], explaining and quantifying the ‘wiggling’ trajectories of multiflagellated bacteria [26], assessing Stokes flow theory of helical propulsion against experiment [27], and modeling sperm interaction with epithelium [28]. The above gives a necessarily incomplete sample but—along with the other articles in this Special Issue—should indicate the breadth of application and impact of the method. These studies employed a mixture of geometric representations for swimmers and other biological entities, including particles, lines, and most commonly, surface discretisations, the latter emphasising the continuing relevance of the boundary integral equation framework.

The significant advantage of the regularised stokeslet method in implementational ease does, however, come at an additional computational cost. This issue was clear in the initial analysis of the method in three dimensions [12], which derived both the linear dependence of the regularisation error on the regularisation parameter $ϵ$ and inverse dependence of the quadrature error on this parameter (in Reference [12], the quadrature error was given as $O\left({ϵ}^{-3}\right)$, although this can be improved to $O\left({ϵ}^{-1}\right)$ [29]). The reduction of the regularisation error therefore necessitates the use of more quadrature points, and within a Nyström discretisation, more degrees of freedom for the unknown force distribution and hence a larger (dense) linear system. Several strategies have been developed to address this. Smith [30] applied a combination of constant panel boundary elements, and slender body theory–style line integrals to the method of regularised stokeslets, removing the dependence of the force discretisation error and hence linear system size on the regularisation parameter, and moreover enabling adaptive quadrature of the stokeslet integral. Slender body theory and approaches based on regularised stokeslet line integrals have since been developed in much more depth (to give an example, refs. [31,32,33]); the boundary element approach has since found an application to embryonic nodal cilia [34,35], and with higher order force discretisation, C. elegans swimming [36], phoretic propulsion [37,38], and bacteria shape evolution [39].

Nevertheless, and particularly for systems in which nonslender bodies play a significant role, the boundary element discretisation still requires the construction of a true mesh, making this approach less technically straightforward than the original Nyström discretisation. Barrero-Gil [40] described a method based on augmenting the discretisation of the ‘inner’ integral with additional stokeslet points in order to weaken the dependence of the error on the regularisation parameter. Smith [41] and Gallagher [42] formulated a method based on a coarse discretisation for force and finer discretisation for quadrature, with the force being interpolated from the coarse to fine discretisations by nearest-neighbor interpolation, leading to the nearest-neighbor regularised stokeslet method. This approach, analysed in Reference [29], removes the dependence of the linear system size on the regularisation parameter and moreover (for disjoint force and quadrature discretisations) weakens the dependence of the quadrature error on the regularisation parameter. The resulting method which is quite simple to implement and analyse, possesses the accuracy and efficiency necessary to model systems with greater than 100 cilia and surrounding domain [43], and moreover lends itself to implementation through basic linear algebra operations, exploiting underlying hardware and software parallelisation [44].

Other key areas of methodological development for the method of regularised stokeslets include the development of image systems for plane boundaries [13,45,46]; extension to Brinkman/oscillatory Stokes flow [47], triply [48], doubly [49,50] and singly [51] periodic boundary conditions; the use of radial basis functions to represent force distributions [52]; improvement to the near-field regularisation error [53], far-field regularisation error [54]; Richardson extrapolation in the regularisation parameter [55]; and perhaps most powerfully, methods based on kernel-independent fast multipole method [56,57] and treecode [58]. Regularisation at the level of the boundary integral equation itself combined with asymptotic corrections (a related but different perspective from regularisation of the equation from which the fundamental solution is derived) has been shown to provide highly accurate results for problems including the near contact of droplets without the need for specialised quadrature [59,60].

An aspect which has received relatively less attention in the regularised stokeslet literature is the role of the double-layer potential, which strictly is necessary for the boundary integral representation of any problem in which the boundary motion cannot be decomposed into rigid body motions. Elimination of the double layer is possible in certain situations, for example, flows around nonrigid volume-conserving bodies; however, a modification to the single-layer density is required, which complicates the calculation of the total force and moment (and solution of swimming problems which have nonzero total force and moment) and may affect the calculation of the rate of working. Notable exceptions include Spagnolie and Lauga’s application of the method of regularised stokeslets with both single- and double-layer potential to model slip-velocity squirmers interacting with a plane boundary [61] and Montenegro-Johnson, Lauga and colleagues’ work on both phoretic [37] and undulating [36] swimmers.

Given that many problems of interest, particularly those involving bending flagella, undulating cell bodies, effective slip flows produced by cilia [62], and nanofluidics applications [63] cannot strictly be represented by combinations of rigid body motions, we continue this line of research in the present manuscript. We describe how to implement the double-layer potential in the context of the nearest-neighbor regularised stokeslet method, before testing the method for two problems with known analytical solutions: the calculation of the drag on a translating sphere with free-slip or partial-slip boundary condition and the calculation of the translational velocity and rate of working of a slip-velocity spherical squirmer in an infinite fluid. Finally, we apply the method to undulatory slender and nonslender swimmers (to assess whether slenderness affects the surface deformation in a way that influences the importance of the double-layer potential) both in infinite fluid and between plane boundaries and assess the effect of the neglect of the double-layer potential on the calculation of rate of working.

In the remainder of the Introduction, we briefly recapitulate the regularised stokeslet boundary integral equation and the circumstances in which the elimination of the double-layer potential is formally justified, before describing the methodology for meshless discretisation of the double-layer potential and some example applications.

1.2. Mathematical Background of Regularised Stokeslets and the Double-Layer Potential

The linearity of Equation (1) enables boundary conditions to be satisfied by taking discrete sums, line or surface integrals of fundamental solutions. The classical (singular) stokeslet or Oseen tensor is the solution to

$$\begin{array}{cc}\hfill -\nabla p+\mu {\nabla }^2\mathit{u}+\delta (\mathit{x}-\mathit{y}){\widehat{\mathit{e}}}_j& =\mathbf{0},\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \nabla ·\mathit{u}& =0,\hfill \end{array}$$

(3)

where $\delta$ is the three dimensional Dirac delta function and ${\widehat{\mathit{e}}}_k$ is the unit basis vector pointing in the k-direction. Defining,

$$\begin{array}{cc}\hfill P_j(\mathit{x},\mathit{y})& =2\frac{x_j-y_j}{{|\mathit{x}-\mathit{y}|}^3},\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill S_{ij}(\mathit{x},\mathit{y})& =\frac{{\delta }_{ij}}{|\mathit{x}-\mathit{y}|}+\frac{(x_i-y_i)(x_j-y_j)}{{|\mathit{x}-\mathit{y}|}^3},\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill T_{ijk}(\mathit{x},\mathit{y})& =-\frac{6(x_i-y_i)(x_j-y_j)(x_k-y_k)}{{|\mathit{x}-\mathit{y}|}^5}.\hfill \end{array}$$

(6)

the pair of tensors $(S_{ij},P_j)$ then provide the solutions $\mathit{u}={\left(8\pi \mu \right)}^{-1}(S_{1j},S_{2j},S_{3j})$ and $p={\left(8\pi \right)}^{-1}{P}_j$ to Equation (3), with ${\sigma }_{ij}={\left(8\pi \mu \right)}^{-1}{T}_{ijk}$ the corresponding stress.

The regularised stokeslet [11,12] is the exact solution to the Stokes flow equations with spatially smoothed point force,

$$\begin{array}{cc}\hfill -\nabla p+\mu {\nabla }^2\mathit{u}+{\varphi }_{ϵ}(\mathit{x}-\mathit{y}){\widehat{\mathit{e}}}_j& =\mathbf{0},\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \nabla ·\mathit{u}& =0.\hfill \end{array}$$

(8)

where ${\varphi }_{ϵ}\left(\mathit{x}\right)$ is a family of “blob” functions, which approximates $\delta \left(\mathit{x}\right)$ as $ϵ\to 0$ and $k=1,2,3$. With the most frequently used choice of

$${\varphi }_{ϵ}\left(\mathit{x}\right)=\frac{15{ϵ}^4}{{8\pi \left(\right|\mathit{x}|}^2+{ϵ}^2{)}^{7/2}},$$

(9)

the regularised stokeslet pressure, velocity and stress fields are defined by

$$\begin{array}{cc}\hfill P_j^{ϵ}(\mathit{x},\mathit{y})& =\frac{x_j-y_j}{{(|\mathit{x}-\mathit{y}|}^2+{ϵ}^2{)}^{5/2}}{\left(2\right|\mathit{x}-\mathit{y}|}^2+5{ϵ}^2),\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill S_{ij}^{ϵ}(\mathit{x},\mathit{y})& =\frac{{\delta }_{ij}{(|\mathit{x}-\mathit{y}|}^2+2{ϵ}^2)+(x_i-y_i)(x_j-y_j)}{{(|\mathit{x}-\mathit{y}|}^2+{ϵ}^2{)}^{3/2}},\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill T_{ijk}^{ϵ}(\mathit{x},\mathit{y})& =-\frac{6(x_i-y_i)(x_j-y_j)(x_k-y_k)}{{(|\mathit{x}-\mathit{y}|}^2+{ϵ}^2{)}^{5/2}}\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill & -\frac{3{ϵ}^2[(x_i-y_i){\delta }_{jk}+(x_j-y_j){\delta }_{ki}+(x_k-y_k){\delta }_{ij}]}{{(|\mathit{x}-\mathit{y}|}^2+{ϵ}^2{)}^{5/2}}.\hfill \end{array}$$

(13)

Hence,

$$p={\left(8\pi \right)}^{-1}{P}_j^{ϵ},u_i={\left(8\pi \mu \right)}^{-1}{S}_{ij}^{ϵ},\mathrm{and}{\sigma }_{ik}=T_{ijk}^{ϵ}={(8\pi )}^{-1}(-P_j^{ϵ}{\delta }_{ik}+\mu ({\partial }_k{S}_{ij}^{ϵ}+{\partial }_i{S}_{kj}^{ϵ}))$$

define the pressure, velocity and stress due to the spatially smoothed point force.

The method of regularised stokeslets is then based on formulating a boundary integral equation formulated in terms of the near-singular fundamental solutions (10)–(13). Following [12], a ‘regularised stokeslet version’ of the Lorentz reciprocal theorem can be derived for any Stokes flow $(\mathit{u},p)$ and point $\mathit{y}$,

$$\frac{1}{8\pi \mu }\frac{\partial }{\partial x_k}\left(S_{ij}^{ϵ}(\mathit{x},\mathit{y}){\sigma }_{ik}\left(\mathit{x}\right)-\mu u_i\left(\mathit{x}\right)T_{ijk}^{ϵ}(\mathit{x},\mathit{y})\right)=u_j\left(\mathit{x}\right){\varphi }_{ϵ}(\mathit{x}-\mathit{y}).$$

(14)

In Equation (14) and throughout, repeated $i,j,k$ indices imply summation over $\{1,2,3\}$.

Equipped with Equation (14), we now consider Stokes flow in the unbounded three-dimensional space external to a volume D with smooth, orientable surface $\partial D$. This volume could represent, for example, a sedimenting rigid body, a droplet or a swimming cell. We take D to be topologically open, so that its boundary $\partial D\subset {\mathbb{R}}^3\backslash D$. Define $B_R$ to be a ball of radius R sufficiently large to contain D, and consider the region consisting of the boundary $\partial D$ and exterior volume of fluid which is given by ${\mathsf{\Omega }}_R=B_R\backslash D$. Integrating over the volume ${\mathsf{\Omega }}_R$, we have for any point $\mathit{y}$ either exterior to D or on its surface,

$$\frac{1}{8\pi \mu }{\int \int \int }_{{\mathsf{\Omega }}_R}\frac{\partial }{\partial x_k}\left(S_{ij}^{ϵ}(\mathit{x},\mathit{y}){\sigma }_{ik}\left(\mathit{x}\right)-\mu u_i\left(\mathit{x}\right)T_{ijk}^{ϵ}(\mathit{x},\mathit{y})\right)d{V}_{\mathit{x}}={\int \int \int }_{{\mathsf{\Omega }}_R}{u}_j\left(\mathit{x}\right){\varphi }_{ϵ}(\mathit{x}-\mathit{y})d{V}_{\mathit{x}}.$$

(15)

Denoting the outward pointing unit normal to ${\mathsf{\Omega }}_R$ by $\mathit{n}$ (inward pointing to D on $\partial D$), applying the Divergence Theorem, and finally taking the limit as $R\to \infty$ then yields the (exact) boundary integral equation

$$\begin{array}{cc}\hfill & -\frac{1}{8\pi \mu }{\int \int }_{\partial D}{S}_{ij}^{ϵ}(\mathit{x},\mathit{y})f_i\left(\mathit{x}\right)d{S}_{\mathit{x}}-\frac{1}{8\pi }{\int \int }_{\partial D}{u}_i\left(\mathit{x}\right)T_{ijk}^{ϵ}(\mathit{x},\mathit{y})n_k\left(\mathit{x}\right)d{S}_{\mathit{x}}\hfill \\ \hfill & ={\int \int \int }_{{\mathbb{R}}^3\backslash D}{u}_j\left(\mathit{x}\right){\varphi }_{ϵ}(\mathit{x}-\mathit{y})d{V}_{\mathit{x}},\hfill \end{array}$$

(16)

where $f_i:=-{\sigma }_{ik}{n}_k$ is the traction, i.e., the force per unit area exerted on the body by the fluid. For the standard singular boundary integral equation, the right-hand side of Equation (16) would involve a Dirac delta function, thereby evaluting precisely to $c\left(\mathit{y}\right)u_j\left(\mathit{y}\right)$, where $c\left(\mathit{y}\right)=1$ if $\mathit{y}$ is interior to the fluid region ${\mathbb{R}}^3\backslash D$. In the case that $\mathit{y}\in D$, then this coefficient is adjusted based on the solid angle exterior pointing into the fluid at $\mathit{y}$; in the usual case that the surface is smooth, then $c\left(\mathit{y}\right)=1/2$.

In the regularised case, the right-hand side of Equation (16) can only be written approximately in terms of $u_j\left(\mathit{y}\right)$ and an error which is linear in $ϵ$ in the vicinity of $\partial D$ and quadratic otherwise. Given that calculations typically involve boundary collocation, we use the worst case of linear error throughout, yielding the regularised approximate boundary integral equation, valid for $\mathit{y}\in {\mathbb{R}}^3\backslash D$, i.e., on the boundary $\partial D$ or exterior to D,

$$\begin{array}{cc}\hfill -\underset{{\mathrm{SLP}}_j^{ϵ}(\mathit{y};\mathit{f};\partial D)}{\underbrace{\frac{1}{8\pi \mu }{\int \int }_{\partial D}{S}_{ij}^{ϵ}(\mathit{x},\mathit{y})f_i\left(\mathit{x}\right)d{S}_{\mathit{x}}}}-\underset{{\mathrm{DLP}}_j^{ϵ}(\mathit{y};\mathit{u};\partial D)}{\underbrace{\frac{1}{8\pi }{\int \int }_{\partial D}{u}_i\left(\mathit{x}\right)T_{ijk}^{ϵ}(\mathit{x},\mathit{y})n_k\left(\mathit{x}\right)d{S}_{\mathit{x}}}}& =c\left(\mathit{y}\right)u_j\left(\mathit{y}\right)\hfill \\ \hfill & +O\left(\overline{\kappa }\left(\mathit{y}\right)ϵ\right),\hfill \end{array}$$

(17)

where $\overline{\kappa }\left(\mathit{y}\right)$ is the mean curvature of $\partial D$ at $\mathit{y}$. The dependency on curvature is specific to the full equation with double layer not eliminated and associated approximation of ${\int \int \int }_{{\mathbb{R}}^3\backslash D}{u}_j{\varphi }_{ϵ}dV$, as opposed to ${\int \int \int }_{{\mathbb{R}}^3}{u}_j{\varphi }_{ϵ}dV$ (see Reference [38], which extends the analysis of Cortez et al. [12]).

The left-hand sides of Equations (16) and (17) consist of two terms: a surface integral of regularised stokeslets $S_{ij}^{ϵ}$ multiplying the traction, referred to (by analogy with electric potential theory) as the single-layer potential (SLP), and a surface integral of the regularised stokeslet stress $T_{ijk}^{ϵ}{n}_k$ multiplying the surface velocity, referred to as the double-layer potential (DLP). A similar equation can be derived for situations involving completely or partially bounded Stokes flow (Appendix A).

The ‘more nearly singular’ ($O\left({ϵ}^{-2}\right)$ as $\mathit{x}\to \mathit{y}$) behaviour of the stress tensor $T_{ijk}^{ϵ}$ has motivated most studies using the method of regularised stokeslets to eliminate this term, leaving a single-layer equation. There are (at least) two situations in which this elimination can be mathematically justified through identical arguments to singular boundary integral theory—one of which involves an important caveat.

Suppose that D is undergoing rigid body motion, equivalent to the flow throughout the interior and surface $D\cup \partial D$ being given by $u_i\left(\mathit{x}\right)=U_i+{\epsilon }_{ijk}{\mathsf{\Omega }}_j{x}_k$ where $\mathit{U}$ is velocity and $\mathsf{\Omega }$ is angular velocity. Then, the double-layer potential can be simplified as (Appendix B),

$$\begin{array}{cc}\hfill {\mathrm{DLP}}_j^{ϵ}(\mathit{y};\mathit{u};\partial D)& ={\int \int \int }_D{u}_j\left(\mathit{x}\right){\varphi }_{ϵ}(\mathit{x}-\mathit{y})d{V}_{\mathit{x}},\mathrm{for}\mathrm{all}\mathit{y}\in {\mathbb{R}}^3\backslash D.\hfill \end{array}$$

(18)

Substituting into Equation (16), we then have, for $\mathit{y}$ on the surface of, or exterior to the body,

$$\begin{array}{cc}\hfill -{\mathrm{SLP}}_j^{ϵ}(\mathit{y};\mathit{f};\partial D)& ={\int \int \int }_D{u}_j{\varphi }_{ϵ}(\mathit{x}-\mathit{y})d{V}_{\mathit{x}}+{\int \int \int }_{{\mathbb{R}}^3\backslash D}{u}_j\left(\mathit{x}\right){\varphi }_{ϵ}(\mathit{x}-\mathit{y})d{V}_{\mathit{x}},\hfill \\ \hfill & ={\int \int \int }_{{\mathbb{R}}^3}{u}_j\left(\mathit{x}\right){\varphi }_{ϵ}(\mathit{x}-\mathit{y})d{V}_{\mathit{x}}\mathrm{for}\mathrm{all}\mathit{y}\in {\mathbb{R}}^3\backslash D.\hfill \end{array}$$

(19)

The right-hand side can then be approximated [12] by $u_j\left(\mathit{y}\right)+O\left(ϵ\right)$ (note the lack of dependence on $\overline{\kappa }$ in this case).

By similar arguments, Equation (18) can also be approximated, up to the regularisation error,

$${\mathrm{DLP}}_j^{ϵ}(\mathit{y};\mathit{u};\partial D)=c\left(\mathit{y}\right)u_j\left(\mathit{y}\right)+O\left(\overline{\kappa }\left(\mathit{y}\right)ϵ\right),\mathrm{for}\mathit{y}\in \partial D.$$

(20)

This identity, which is an approximate form of a well-established result for the standard boundary integral equation [64], is useful for solving the swimming problem with the double-layer potential.

The second situation in which the double-layer potential may be eliminated refers to bodies which are not necessarily rigid but which do not change volume—for example, flexible swimming cells or certain models of ciliates involving a surface tangential slip velocity, i.e.,

$${\int \int }_{\partial D}{u}_i\left(\mathit{x}\right)n_i\left(\mathit{x}\right)d{S}_{\mathit{x}}=0.$$

(21)

This condition implies the existence of a (unique) solution to the Stokes flow equations in the region interior to D, with velocity $u_i^{*}\left(\mathit{x}\right)$ and stress tensor ${\sigma }_{ij}^{*}\left(\mathit{x}\right)$ satisfying $u_i^{*}\left(\mathit{y}\right)=u_i\left(\mathit{y}\right)$ for all $\mathit{y}\in D$ (a ‘fictitious solution’). Applying Equations (14) and (16) leads to (Appendix C) the single-layer approximate boundary integral equation,

$$-{\mathrm{SLP}}_j^{ϵ}(\mathit{y};\mathit{q};\partial D)=u_j\left(\mathit{y}\right),\mathrm{all}\mathit{y}\in \partial D,$$

(22)

where for brevity here and below the regularisation error associated with the approximation of the right-hand side by $u_j\left(\mathit{y}\right)$ is omitted but implied.

The key issue for practical application is that Equation (22) is not formulated in terms of the physical traction $\mathit{f}$, which may present difficulties if the total force and moment must be determined or are part of the specification of the problem. For example, the resistance problem involves specifying $u_j\left(\mathit{y}\right)$ for $\mathit{y}\in \partial D$ and determining the force and moment,

$$\mathcal{F}:={\int \int }_{\partial D}\mathit{f}\left(\mathit{x}\right)d{S}_{\mathit{x}},\mathcal{M}:={\int \int }_{\partial D}\mathit{x}\wedge \mathit{f}\left(\mathit{x}\right)d{S}_{\mathit{x}}.$$

(23)

If the adjusted density $\mathit{q}$ but not the physical traction $\mathit{f}$ is known, then neither $\mathcal{F}$ nor $\mathcal{M}$ can be determined directly; therefore, the resistance problem cannot in general be solved from Equation (22) accurately if the body is not undergoing a rigid body motion. This issue is exhibited in Section 3.1.

The swimming problem involves specifying the force and moment $\mathcal{F}$ and $\mathcal{M}$, along with the surface velocity relative to a body frame moving with an unknown rigid body motion $\mathit{U}+\mathbf{\Omega }\wedge \mathit{x}$ and then determining the surface traction $\mathit{f}$ along with $\mathit{U}$ and $\mathbf{\Omega }$. Again, the fact that Equation (A6) involves only $\mathit{q}$ is problematic. In the particular case $\mathcal{F}=\mathbf{0}=\mathcal{M}$, the characteristic of inertialess and neutrally buoyant swimming, the fictitious flow stress can be shown to result in zero total force and moment (Appendix D), and hence, the conditions in terms of the adjusted single-layer density,

$${\int \int }_{\partial D}\mathit{q}\left(\mathit{x}\right)d{S}_{\mathit{x}}=\mathbf{0},{\int \int }_{\partial D}\mathit{x}\wedge \mathit{q}\left(\mathit{x}\right)d{S}_{\mathit{x}}=\mathbf{0},$$

(24)

are equivalent to the conditions in terms of the physical density $\mathit{f}$,

$${\int \int }_{\partial D}\mathit{f}\left(\mathit{x}\right)d{S}_{\mathit{x}}=\mathbf{0},{\int \int }_{\partial D}\mathit{x}\wedge \mathit{f}\left(\mathit{x}\right)d{S}_{\mathit{x}}=\mathbf{0},$$

(25)

The zero-force and moment, volume conserving swimming problem can therefore be formulated in terms of the single-layer potential only, giving (up to numerical error) equivalent results for the swimming velocity $\mathit{U}$ and angular velocity $\mathbf{\Omega }$—although quantities derived from $\mathit{f}$ such as the rate of working ${\int \int }_{\partial D}\mathit{u}·\mathit{f}dS$ may not in general be given accurately by the equivalent calculation from $\mathit{q}$. This issue is evident in the model of the squirming swimmer (Results Section 3.2).

1.3. Numerical Discretisation

The final aspect to introduce is the spatial discretisation of the force and quadrature as typically implemented for the Nyström or nearest-neighbor discretisations. First, consider the resistance problem, i.e., the calculation of the total force and moment due to rigid body motion $u_j\left(\mathit{y}\right)=U_j+{ϵ}_{jk\ell }{\mathsf{\Omega }}_k{y}_{\ell }$. From Equation (19), we have (subject to regularisation error),

$$-{\mathrm{SLP}}_j(\mathit{y};\mathit{f};\partial D)=U_j+{\epsilon }_{jk\ell }{\mathsf{\Omega }}_k{y}_{\ell }\mathrm{for}\mathrm{all}\mathit{y}\in \partial D.$$

(26)

The Nyström discretisation involves approximating the integral in the single-layer potential by a quadrature rule with abscissae $\left\{\mathit{x}\right[1],\dots ,\mathit{x}[N\left]\right\}$ and weights $\left\{w\right[1],\dots ,w[N\left]\right\}$. Then, if the surface metric is $dS\left(\mathit{x}\right[n\left]\right)$, the area associated to each node is given by $A\left[n\right]:=w\left[n\right]dS\left(\mathit{x}\right[n\left]\right)$, and hence, Equation (26) has discrete approximation,

$$-\frac{1}{8\pi \mu }\sum _{n=1}^N{S}_{ij}^{ϵ}(\mathit{x}\left[n\right],\mathit{x}\left[m\right])f_i\left(\mathit{x}\left[n\right]\right)A\left[n\right]=U_j+{\epsilon }_{jk\ell }{\mathsf{\Omega }}_k{x}_{\ell }\left[m\right],\mathrm{for}\mathrm{all}m=1,\dots N.$$

(27)

Denoting $F_i\left[n\right]:=-f_k\left(\mathit{x}\left[n\right]\right)A\left[n\right]$ then yields a system of $3N$ equations in the $3N$ unknowns $F_i\left[1\right],\dots ,F_i\left[N\right]$. The simplicity of Equation (27) and absence of the need for a true mesh with local geometry, or indeed singular quadratures, constitute significant practical advantages over the standard boundary integral method.

Nevertheless, this simplicity comes at the cost of the matrix system size depending on the regularisation parameter. Recall that Equation (19) possesses an error which is linear in the regularisation parameter. Moreover, it can be shown that the quadrature error of Equation (27) is $O(h^2/ϵ)$, where h is the quadrature spacing [29]. Therefore, reducing $ϵ$ with the aim of reducing the regularisation error then demands a finer spatial discretisation in order to control the quadrature error. The size of the resulting linear system then grows, rapidly escalating the computational cost. As discussed above, the strategies to address this issue include improving the order of the regularisation error [53,55] and boundary elements/regularised stokeslet segments [30,33].

The approach we focus on is nearest-neighbor interpolation, which preserves the simple meshless nature of the method of regularised stokeslets while significantly improving efficiency and accuracy. This method [41,42] utilises coarse force and fine quadrature discretisations, with nearest-neighbor interpolation being used to project from the fine-to-coarse discretisation. This approach removes the coupling between the system size and regularisation parameter, and for nonoverlapping force/quadrature sets, it also removes the dependence of the quadrature error on the regularisation parameter. Moreover, the method can be implemented in terms of basic linear algebra operations to exploit associated hardware and software optimisations [44].

The basis for the nearest-neighbor method is the use of two point cloud discretisations, one coarse set $\left\{\mathit{x}\right[1],\dots ,\mathit{x}[N\left]\right\}$ which is sufficient to capture the variation of the traction, and which dictates the size of the linear system, and another finer quadrature discretisation set $\left\{\mathit{X}\right[1],\dots ,\mathit{X}[Q\left]\right\}$, which has sufficient resolution to capture the rapidly varying kernel $S_{ij}^{ϵ}(\mathit{x},\mathit{y})$ for $\mathit{y}$ in the vicinity of $\mathit{x}$. The slowly varying force per unit area at a quadrature point $\mathit{f}\left(\mathit{X}\right[q\left]\right)$ is approximated by its value at the nearest point on the coarse force set $\mathit{f}\left(\mathit{x}\left[\widehat{n}\right]\right)$; defining $\nu [q,\widehat{n}]=1$ and $\nu [q,n]=0$ for $n\ne \widehat{n}$, and a similar approximation is made for the surface metric $dS$. Hence, the discrete problem takes the slightly modified form,

$$u_j\left(\mathit{x}\left[m\right]\right)=\frac{1}{8\pi \mu }\sum _{n=1}^N\left(\sum _{q=1}^Q{S}_{ij}^{ϵ}(\mathit{X}\left[q\right],\mathit{x}\left[m\right])\nu [q,n]\right)F_i\left[n\right],$$

(28)

where $F_i\left[n\right]:=-f_i\left(\mathit{x}\left[n\right]\right)dS\left(\mathit{x}\left[n\right]\right)$ can be interpreted as the force on each quadrature point neighbouring $\mathit{x}\left[n\right]$. The error analysis of this method is given in Reference [29] and discussed further in Reference [55]; the most important aspects are that with quadrature spacing $h_q$, the quadrature error is $O(h_q^2/ϵ)$ at collocation points which are contained in the quadrature set, and $O(h_q^3/{\delta }^2)$ at collocation points which have minimum distance $\delta$ from the quadrature set (hence, no asymptotic dependence as $ϵ\to 0$); the force discretisation error and hence linear system size have no intrinsic dependence on $ϵ$. In the following section, we focus on the discretisation of the double-layer potential within this framework.

2. Materials and Methods

The additional complications to including the double-layer potential in the method of regularised stokeslets are:

• Evaluation of the higher order near-singularity ($T_{ijk}^{ϵ}=O(1/{ϵ}^2)$ for $|\mathit{x}-\mathit{y}|\to 0$).
• The need to calculate the surface normal $\mathit{n}\left(\mathit{y}\right)$ without a true mesh/local geometry.
• The fact that the surface metric $dS\left(\mathit{y}\right)$ must be calculated explicitly for the implementation of the double-layer term, as opposed to being absorbed into the force density, as is the case for the single-layer term.

Issue 1. is straightforward to solve, using (the regularised version of) a well-known technique from the boundary element literature. Equation (18) with the specific choice of rigid body motion $U_j={\delta }_{jk}$, ${\mathsf{\Omega }}_j=0$ and $\mathit{y}\in \partial D$ yields (up to regularisation error),

$$\begin{array}{c}\hfill \frac{1}{8\pi }{\int \int }_{\partial D}{T}_{ijk}^{ϵ}(\mathit{x},\mathit{y})n_k\left(\mathit{y}\right)d{S}_{\mathit{x}}={\int \int \int }_D{\delta }_{ij}{\varphi }_{ϵ}(\mathit{x}-\mathit{y})d{V}_{\mathit{x}}={\delta }_{ij}c\left(\mathit{y}\right),\mathrm{for}\mathit{y}\in \partial D.\end{array}$$

(29)

Recall that if $\partial D$ is smooth in the neighborhood of $\mathit{y}$ then $c\left(\mathit{y}\right)=1/2$. This identity is a mathematical fact about the double-layer potential that can be exploited to construct an approximation for the inner quadrature even when the surface motion is not rigid, as follows: splitting the surface integral into an inner component on $\partial D_{\delta }\left(\mathit{y}\right)=\{\mathit{x}\in \partial D:|\mathit{x}-\mathit{y}|<\delta \}$ sufficiently small that the variation in the velocity is negligible, and by an outer component on $\partial D\backslash \partial D_{\delta }\left(\mathit{y}\right)$, we may approximate the double-layer potential by,

$$\begin{array}{cc}\hfill {\mathrm{DLP}}_j^{ϵ}(\mathit{y};\mathit{u};\partial D)& =\frac{1}{8\pi }{\int \int }_{\partial D_{\delta }\left(\mathit{y}\right)}{u}_i\left(\mathit{x}\right)T_{ijk}^{ϵ}(\mathit{x},\mathit{y})n_k\left(\mathit{x}\right)d{S}_{\mathit{x}}\hfill \\ \hfill & +\frac{1}{8\pi }{\int \int }_{\partial D\backslash \partial D_{\delta }\left(\mathit{y}\right)}{u}_i\left(\mathit{x}\right)T_{ijk}^{ϵ}(\mathit{x},\mathit{y})n_k\left(\mathit{x}\right)d{S}_{\mathit{x}},\hfill \\ \hfill & \approx \frac{1}{8\pi }{u}_i\left(\mathit{y}\right){\int \int }_{\partial D_{\delta }\left(\mathit{y}\right)}{T}_{ijk}^{ϵ}(\mathit{x},\mathit{y})n_k\left(\mathit{x}\right)d{S}_{\mathit{x}}\hfill \\ \hfill & +\frac{1}{8\pi }{\int \int }_{\partial D\backslash \partial D_{\delta }\left(\mathit{y}\right)}{u}_i\left(\mathit{x}\right)T_{ijk}^{ϵ}(\mathit{x},\mathit{y})n_k\left(\mathit{x}\right)d{S}_{\mathit{x}}\hfill \\ \hfill & \approx u_i\left(\mathit{y}\right)\left(\frac{{\delta }_{ij}}{2}-\frac{1}{8\pi }{\int \int }_{\partial D\backslash \partial D_{\delta }\left(\mathit{y}\right)}{T}_{ijk}^{ϵ}(\mathit{x},\mathit{y})n_k\left(\mathit{x}\right)d{S}_{\mathit{x}}\right)\hfill \\ \hfill & +\frac{1}{8\pi }{\int \int }_{\partial D\backslash \partial D_{\delta }\left(\mathit{x}\right)}{u}_i\left(\mathit{x}\right)T_{ijk}^{ϵ}(\mathit{x},\mathit{y})n_k\left(\mathit{y}\right)d{S}_{\mathit{x}}.\hfill \end{array}$$

(30)

Consider the application of the Nyström discretisation, with collocation at $\mathit{y}=\mathit{x}\left[m\right]$. Taking $\delta$ as the radius of the part of the surface belonging to $\mathit{x}\left[m\right]$, the double-layer operator evaluated at a collocation point can then be approximated by,

$$\begin{array}{cc}\hfill {\mathrm{DLP}}_j^{ϵ}(\mathit{x}\left[m\right];\mathit{u}[·];\partial D)& =\left(\frac{{\delta }_{ij}}{2}-\frac{1}{8\pi }\sum _{n=1,n\ne m}^N{T}_{ijk}^{ϵ}(\mathit{x}\left[n\right],\mathit{x}\left[m\right])n_k\left(\mathit{x}\left[n\right]\right)A\left[n\right]\right)u_i\left(\mathit{x}\left[m\right]\right)\hfill \\ \hfill & +\frac{1}{8\pi }\sum _{n=1,n\ne m}^N{T}_{ijk}^{ϵ}(\mathit{x}\left[n\right],\mathit{x}\left[m\right])n_k\left(\mathit{x}\left[n\right]\right)A\left[n\right]u_j\left(\mathit{x}\left[n\right]\right),\hfill \end{array}$$

(31)

where $A\left[n\right]$ is the surface metric at $\mathit{x}\left[n\right]$ (this term may also be multiplied by a quadrature weight if a higher-order quadrature rule is used).

Via the nearest-neighbor discretisation, the double-layer operator can be approximated by,

$$\begin{array}{cc}\hfill & {\mathrm{DLP}}_j^{ϵ}(\mathit{x}\left[m\right];\mathit{u}[·];\partial D)\hfill \\ \hfill & =\left(\frac{{\delta }_{ij}}{2}-\frac{1}{8\pi }\sum _{n=1}^N\left(\sum _{q=1,q\ne m}^Q{T}_{ijk}^{ϵ}(\mathit{X}\left[q\right],\mathit{x}\left[m\right])\nu [q,n]\right)n_k\left(\mathit{X}\left[q\right]\right)A\left[q\right]\right)u_j\left(\mathit{x}\left[m\right]\right)\hfill \\ \hfill & +\frac{1}{8\pi }\sum _{n=1}^N\left(\sum _{q=1,q\ne m}^Q{T}_{ijk}^{ϵ}(\mathit{X}\left[q\right],\mathit{x}\left[m\right])\nu [q,n]n_k\left(\mathit{X}\left[q\right]\right)A\left[q\right]\right)u_j\left(\mathit{x}\left[n\right]\right).\hfill \end{array}$$

(32)

If evaluating the double-layer potential at any point $\mathit{y}\notin \partial D$, then the simpler approximation

$${\mathrm{DLP}}_j^{ϵ}(\mathit{x}\left[m\right];\mathit{u}[·];\partial D)=\frac{1}{8\pi \mu }\sum _{n=1}^N\sum _{q=1}^Q{T}_{ijk}^{ϵ}(\mathit{X}\left[q\right],\mathit{x}\left[m\right])\nu [q,n]n_k\left(\mathit{X}\left[q\right]\right)A\left[q\right]u_j\left(x\left[n\right]\right)$$

(33)

may be used. We briefly note that in the MATLAB^® implementation accompanying this paper, the identity $T_{ijk}^{ϵ}(\mathit{x},\mathit{y})=-T_{ijk}^{ϵ}(\mathit{y},\mathit{x})$ is used, such that matrix rows correspond to collocation points and matrix columns to quadrature points.

Issues 2 and 3 require an approximate reconstruction of the local geometry from the quadrature discretisation. Below, we describe a method to accomplish this using principal component analysis, implemented via the accompanying MATLAB^® (Mathworks, Natick, MA, USA) code that can be applied to any point cloud discretisation.

Given a quadrature discretisation $\left\{\mathit{X}\right[1],\dots ,\mathit{X}[Q\left]\right\}$, and any point $\mathit{y}\in B$, the discrete neighborhood of the closest J points to $\mathit{y}$ is denoted ${\mathcal{N}}_J(\mathit{y};B)=\{\mathit{X}\left[q_1\right],\dots ,\mathit{X}\left[q_J\right]\}$. The choice $J=9$ has been found to work in practice and is used for the results in this paper (Figure 1). Performing principal component analysis on the set $\{(\mathit{X}\left[q_1\right]-\mathit{y}),\dots ,(\mathit{X}\left[q_J\right]-\mathit{y})\}$ yields a triplet of eigenvalues ${\lambda }_1>{\lambda }_2>{\lambda }_3$ and corresponding orthogonal directions ${\mathit{v}}_1,{\mathit{v}}_2,{\mathit{v}}_3$, with the first two vectors approximating the tangent space of B at $\mathit{y}$ and the third vector the normal. These vectors are normalised to produce an approximate local basis triple. To ensure that the normal direction is consistently inward pointing, at the discretisation stage, one or more points in the interior of the volume bounded by B are supplied. Utilising the nearest interior point, the sense of ${\mathit{v}}_3$ can be determined and, hence, the inward-pointing normal $\mathit{n}\left(\mathit{y}\right)=\pm {\mathit{v}}_3$.

Finally, it is necessary to approximate the surface metric. For a given point $\mathit{X}\left[q\right]$, in the quadrature discretisation, the nearest (nonself) point $\mathit{X}\left[\overline{q}\right]$ can be found, generating a basis vector for the tangent space ${\mathit{X}}_a\left[q\right]:=\mathit{X}\left[\overline{q}\right]-\mathit{X}\left[q\right]$. Vectors from $\mathit{X}\left[q\right]$ to the other neighboring points are constructed and then projected onto ${\mathit{X}}_a\left[q\right]$. Finding the neighbor $\mathit{X}\left[\check{q}\right]$ with the minimum component along this direction, the vector ${\mathit{X}}_b\left[q\right]:=\mathit{X}\left[\check{q}\right]-\mathit{X}\left[q\right]$ almost certainly contains a significant component orthogonal to ${\mathit{X}}_a\left[q\right]$. The surface metric and normal can then be approximated using a similar formula to the boundary element method [65],

$$\begin{array}{cc}\hfill A\left[q\right]& \approx |{\mathit{X}}_a\left[q\right]\wedge {\mathit{X}}_b\left[q\right]|,\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill \mathit{n}\left[q\right]& \approx \frac{{\mathit{X}}_a\left[q\right]\wedge {\mathit{X}}_b\left[q\right]}{|{\mathit{X}}_a\left[q\right]\wedge {\mathit{X}}_b\left[q\right]|}.\hfill \end{array}$$

(35)

The accuracy of this method in approximating the surface area of a sphere for various values of point spacing h is shown in Figure 1, showing approximately linear convergence for each of $J=4,6,9$.

2.1. Problem 1: Drag on a Translating Sphere with No-Slip, Partial-Slip, or Free-Slip Boundary Conditions

Consider the problem of calculating the drag force $\mathcal{F}={\int \int }_B\mathit{f}\left(\mathit{x}\right)d{S}_{\mathit{x}}$ on a sphere B translating with unit velocity $\mathit{U}={(1,0,0)}^T$, with three choices of boundary condition,

N

no-slip, no-penetration, i.e., $\mathit{u}\left(\mathit{x}\right)=\mathit{U}·\mathit{n}$ for all $\mathit{x}\in B$.

P

partial-slip, no-penetration, i.e., $\mathit{u}\left(\mathit{x}\right)·\mathit{n}\left(\mathit{x}\right)=\mathit{U}·\mathit{n}$ and $\lambda (\mathit{I}-\mathit{n}(\mathit{x}\left)\mathit{n}\right(\mathit{x}\left)\right)·\mathit{f}\left(\mathit{x}\right)=\mu (\mathit{I}-\mathit{n}(\mathit{x}\left)\mathit{n}\right(\mathit{x}\left)\right)·\left(\mathit{u}\right(\mathit{x})-\mathit{U})$ for all $\mathit{x}\in B$.

F

free-slip, no-penetration, i.e., $\mathit{u}\left(\mathit{x}\right)·\mathit{n}\left(\mathit{x}\right)=\mathit{U}·\mathit{n}$ and $(\mathit{I}-\mathit{n}(\mathit{x}\left)\mathit{n}\right(\mathit{x}\left)\right)·\mathit{f}\left(\mathit{x}\right)=\mathbf{0}$ for all $\mathit{x}\in B$.

Problems [N] and [F] can be viewed as special cases of problem [P] in the limits as the slip length $\lambda \to 0$ and $\lambda \to \infty$, respectively. The solution to problem [N] is well known as the Stokes law, $\mathcal{F}=6\pi a\mathit{U}$ where a is the radius of the sphere, and the solution to problem [F] is $\mathcal{F}=4\pi a\mathit{U}$ [66]. Problem [P] has a solution [63],

$$\mathcal{F}=6\pi a\frac{1+2\lambda /a}{1+3\lambda /a},$$

(36)

which approaches the limits for [N] and [F] as $\lambda \to 0$, ∞, respectively.

The resistance problem can be formulated as a second kind Fredholm integral equation by moving the double-layer potential in Equation (17) across to the right-hand side,

$$\begin{array}{c}\hfill -{\mathrm{SLP}}_j^{ϵ}(\mathit{y};\mathit{f};\partial D)={\mathrm{DLP}}_j(\mathit{y};\mathit{u};\partial D)+c\left(\mathit{y}\right)u_j\left(\mathit{y}\right).\end{array}$$

(37)

For problem [N], the velocity $u_i\left(\mathit{y}\right)={\delta }_{j1}$ on $\partial D$, and so, we formulate the following problem for the traction $\mathit{f}$,

$$\begin{array}{cc}\hfill -{\mathrm{SLP}}_j^{ϵ}(\mathit{y};\mathit{f};\partial D)& ={\mathrm{DLP}}_j^{ϵ}(\mathit{y};(1,0,0);\partial D)+\frac{{\delta }_{j1}}{2}\hfill \\ \hfill & ={\delta }_{j1},\hfill \end{array}$$

(38)

for all $\mathit{y}\in \partial D$, the last line following by Equation (20).

For problem [F], only the normal component of surface velocity is known, i.e., $u_i\left(\mathit{x}\right)n_i\left(\mathit{x}\right)={\delta }_{i1}{n}_i\left(\mathit{x}\right)=n_1\left(\mathit{x}\right)$; the system is completed by the requirement that the unknown traction satisfies the free-slip condition $f_i{t}_i^1=0=f_i{t}_i^2$. Hence, we have the augmented system (with all unknowns on the left-hand side),

$$\begin{array}{cc}\hfill -{\mathrm{SLP}}_j^{ϵ}(\mathit{y};\mathit{f};\partial D)-\frac{u_j\left(\mathit{y}\right)}{2}-{\mathrm{DLP}}_j^{ϵ}(\mathit{y};\mathit{u};\partial D)& =0,\hfill \\ \hfill u_j\left(\mathit{y}\right)n_j\left(\mathit{y}\right)& =n_1\left(\mathit{y}\right),\hfill \\ \hfill f_j\left(\mathit{y}\right)t_j^{\alpha }\left(\mathit{y}\right)& =0,\alpha =1,2,\hfill \end{array}$$

(39)

for all $\mathit{y}\in \partial D$. Equation (39) involves solving for both $\mathit{f}$ and $\mathit{u}$ on the surface; the single and double-layer potentials are discretised as in Equation (32).

Equation (39) can be generalised to provide the following formulation of problem P:

$$\begin{array}{cc}\hfill -{\mathrm{SLP}}_j^{ϵ}(\mathit{y};\mathit{f};\partial D)-\frac{u_j\left(\mathit{y}\right)}{2}-{\mathrm{DLP}}_j^{ϵ}(\mathit{y};\mathit{u};\partial D)& =0,\hfill \\ \hfill u_j\left(\mathit{y}\right)n_j\left(\mathit{y}\right)& =n_1\left(\mathit{y}\right),\hfill \\ \hfill \lambda f_j\left(\mathit{y}\right)t_j^{\alpha }\left(\mathit{y}\right)& =u_j\left(\mathit{y}\right)t_j^{\alpha }\left(\mathit{y}\right),\alpha =1,2,\hfill \end{array}$$

(40)

for all $\mathit{y}\in \partial D$. A key step in the numerical implementation is to utilise the metric $dS\left(\mathit{X}\right[q\left]\right)$ and nearest-neighbor matrix $\nu [q,n]$ in order to relate the $F_j\left[n\right]$ (force) unknowns of the discrete system to the $f_j\left(\mathit{y}\right)$ (force per unit area) variable of the continuous problem. For the purposes of comparison, we also evaluate the effect of (without formal justification) eliminating the double-layer potential, replacing the double-layer term of Equation (40) by $u_j\left(\mathit{y}\right)/2$.

Once the discrete surface force solution $\mathit{F}\left[n\right]:=-\mathit{f}\left(\mathit{x}\right[n\left]\right)A\left[n\right]$ and surface velocity ${\mathit{U}}^s\left[n\right]:=\mathit{u}\left(\mathit{x}\left[n\right]\right)$ have been found for $n=1,\dots ,N$ through numerical solution of Equation (39), the total drag on the sphere is calculated numerically from the weighted sum,

$$\mathcal{F}=\sum _{n=1}^N\sum _{q=1}^Q\nu [q,n]\mathit{F}\left[n\right];$$

(41)

the velocity field at any point $\mathit{y}$ in the fluid exterior to D is calculated from the discrete surface force and velocity solutions via

$$\begin{array}{cc}\hfill u_j\left(\mathit{y}\right)& =\frac{1}{8\pi \mu }\sum _{n=1}^N\sum _{q=1}^Q{S}_{ij}^{ϵ}(\mathit{x}\left[n\right],\mathit{y})F_i\left[n\right]\nu [q,n]\hfill \\ \hfill & -\frac{1}{8\pi \mu }\sum _{n=1}^N\sum _{q=1}^Q{U}_i^s\left[q\right]T_{ijk}^{ϵ}(\mathit{X}\left[q\right],\mathit{y})n_k\left(\mathit{X}\left[q\right]\right)dS\left(\mathit{X}\left[q\right]\right).\hfill \end{array}$$

(42)

2.2. Slip-Velocity Squirmer

The slip-velocity spherical squirmer is one of the canonical models of zero Reynolds number propulsion [62,67,68]; this model provides a coarse-grained representation of microorganisms which produce a surface flow due to the action of many beating cilia. By scaling coordinates so that the radius of the squirmer is 1, we denote the angle relative to the north pole of the squirmer as $\mathsf{\Theta }\left(\mathit{y}\right)$ (so that $0⩽\mathsf{\Theta }\left(\mathit{y}\right)⩽\pi$) for $\left|\mathit{y}\right|=1$; we also define ${\mathit{t}}^{\mathsf{\Theta }}$ to be the basis vector in the tangent space pointing toward the north pole. Taking a single squirming mode, the velocity boundary condition in a frame of reference moving with the squirmer is,

$${\left(u_j\left(\mathit{y}\right)\right)}^{\mathrm{bf}}{n}_j\left(\mathit{y}\right)=0,({\delta }_{ij}-n_i\left(\mathit{y}\right)n_j\left(\mathit{y}\right)){\left(u_j\left(\mathit{y}\right)\right)}^{\mathrm{bf}}=sin\mathsf{\Theta }\left(\mathit{y}\right)t_i^{\mathsf{\Theta }}\left(\mathit{y}\right),\mathrm{for}\mathrm{all}\left|\mathit{y}\right|=1,$$

(43)

where the superscript ‘$\mathrm{bf}$’ denotes a frame of reference in which the swimmer is stationary, referred to as the body frame. If the squirmer frame of reference has translational velocity $\mathit{U}$ and angular velocity $\mathbf{\Omega }$ about the origin, then

$$u_j\left(\mathit{y}\right)=U_j+{\epsilon }_{jk\ell }{\mathsf{\Omega }}_k{y}_{\ell }+{\left(u_j\left(\mathit{y}\right)\right)}^{\mathrm{bf}}.$$

(44)

Because ${\int \int }_{\partial D}{u}_j{n}_{j}dS=0$, the double-layer potential can be eliminated in a similar manner to Equation (26) to provide a single-layer boundary integral equation in terms of the adjusted density $\mathit{q}$,

$$-{\mathrm{SLP}}_j^{ϵ}(\mathit{y};\mathit{q};\partial D)-U_j-{\epsilon }_{jk\ell }{\mathsf{\Omega }}_k{y}_{\ell }={\left(u_j\left(\mathit{y}\right)\right)}^{\mathrm{bf}}.$$

(45)

The left-hand side of the equation is a linear operator on the unknowns $\mathit{q}\left(\mathit{x}\right)$, $\mathit{U}$, and $\mathbf{\Omega }$. As discussed in Appendix D, to determine the swimming velocity $\mathit{U}$ and angular velocity about the origin $\mathbf{\Omega }$, the force- and moment-free conditions can be expressed in terms of $\mathit{q}$, as given in Equation (24). Equations (24) and (43)–(45) therefore provide a complete model from which the swimming velocity can be determined, although the physical density $\mathit{f}$ and, hence, derived quantities such as the rate of working are not available.

Alternatively, equipped with the double-layer implementation, we may avoid the elimination of the double layer and work directly with the physical traction $\mathit{f}$. Again by making use of (20), the integral equation is replaced by,

$$-{\mathrm{SLP}}_j^{ϵ}(\mathit{y};\mathit{q};\partial D)-\left(U_j+{ϵ}_{jk\ell }{\mathsf{\Omega }}_k{y}_{\ell }\right)=\frac{{\left(u_j\left(\mathit{y}\right)\right)}^{\mathrm{bf}}}{2}+{\mathrm{DLP}}_j^{ϵ}\left(\mathit{y};{\mathit{u}}^{\mathrm{bf}};\partial D\right)\mathrm{for}\mathrm{all}\mathit{y}\in \partial D.$$

(46)

Equations (25), (43), (44) and (46) form a complete model for the squirming swimmer in terms of the physical traction $\mathit{f}$. Once the traction is computed, the rate of working is given by the integral,

$$W_f={\int \int }_{\partial D}{u}_j\left(\mathit{y}\right)f_j\left(\mathit{y}\right)d{S}_{\mathit{y}}.$$

(47)

An analogous quantity $W_q$ can be computed from the adjusted density for the purpose of comparison.

The numerical model can be compared directly to the analytical results of Blake [62]. For a single squirming mode,

$${\left(u_{\mathsf{\Theta }}\right)}^{\mathrm{bf}}\left(\mathit{y}\right)=\frac{3}{2}sin\mathsf{\Theta },$$

(48)

the exact swimming velocity and rate of working are given by,

$$\mathit{U}=(0,0,1),W_f=12\pi \mu .$$

(49)

2.3. Undulating Swimmer

The undulating swimmer, resembling a worm or headless eukaryotic flagellum, is described by an origin ${\mathit{X}}^0\left(t\right)$ (in this case, the front tip of the swimmer) and body frame $\mathbf{B}\left(t\right)=({\mathit{b}}_1\left(t\right),{\mathit{b}}_2\left(t\right),{\mathit{b}}_3\left(t\right))$, where ${\mathit{b}}_j={(b_{1j},b_{2j},b_{3j})}^T$ are orthogonal unit vectors. Then, a point ${\mathit{y}}^{\mathrm{bf}}$ in body frame coordinates has laboratory frame coordinates,

$$y_i\left(t\right)=X_i^0\left(t\right)+b_{i\ell }{y}_{\ell }^{\mathrm{bf}}.$$

(50)

By denoting the translational velocity $\mathit{U}={\dot{\mathit{X}}}^0$ and angular velocity $\mathbf{\Omega }$ such that ${\dot{\mathit{b}}}^j=\mathbf{\Omega }\wedge {\mathit{b}}^j$, the kinematics of the swimmer can be expressed by,

$$u_i\left(\mathit{y}\right)=U_i+{\epsilon }_{ijk}{\mathsf{\Omega }}_j(y_{\ell }-X_{\ell }^0)+b_{i\ell }{\dot{y}}_i^{\mathrm{bf}}.$$

(51)

The prescribed waveform swimming problem involves specifying a time-varying function $y_i^{\mathrm{bf}}\left(t\right)$ at each material point on the surface of the swimmer, hence prescribing the body frame velocity ${\dot{y}}_i^{\mathrm{bf}}$. At any instant in time, given the current position ${\mathit{X}}^0$ and orientation $\mathbf{B}$, the aim is to solve for $(\mathit{U},\mathbf{\Omega },\mathit{f})$.

To model an undulatory swimmer in an unbounded fluid, the single-layer boundary integral equation in terms of the adjusted density $\mathit{q}$ is

$$-{\mathrm{SLP}}_j^{ϵ}(\mathit{y};\mathit{q};\partial D)-\left(U_j+{\epsilon }_{jk\ell }{\mathsf{\Omega }}_k(y_{\ell }-X_{\ell }^0)\right)=b_{i\ell }{\dot{y}}_i^{\mathrm{bf}}.$$

(52)

Alternatively, working with the full boundary integral model in terms of the physical traction $\mathit{f}$, we have,

$$-{\mathrm{SLP}}_j^{ϵ}(\mathit{y};\mathit{f};\partial D)-\left(U_j+{\epsilon }_{jk\ell }{\mathsf{\Omega }}_k(y_{\ell }-X_{\ell }^0)\right)=\frac{b_{j\ell }{\dot{y}}_{\ell }^{\mathrm{bf}}}{2}+{\mathrm{DLP}}_j^{ϵ}\left(\mathit{y},\mathbf{B}·{\dot{\mathit{x}}}^{\mathrm{bf}};\partial D\right).$$

(53)

The notation ${\dot{\mathit{x}}}^{\mathrm{bf}}$ in the final term is chosen intentionally to denote a variable of integration (as opposed to the free variable ${\dot{\mathit{y}}}^{\mathrm{bf}}$). Equation (53) can be adapted to include the effect of stationary boundaries; details are given in Appendix E.

The waveform of the centreline of the cell is prescribed in the body frame; because the present study is a proof of principle rather than a biologically focussed analysis of any specific species, we make use of the activated beat of a sperm flagellum by Dresdner and Katz [69]

$$\tilde{y}(\tilde{x},t)=(0.1087\tilde{x}+0.0543)sin(2\pi \tilde{x}-t).$$

(54)

To construct the centreline of the waveform in the body frame, the waveform $(\tilde{x},\tilde{y})$ is translated and rotated so that, the first point of the cell is at ${\mathit{y}}^c=\mathbf{0}$ and the flagellum is aligned with the $y_1^{\mathrm{bf}}$-axis.

The centreline waveform is parameterised as ${\mathit{y}}^c(s,t)$ where $s\in [0,1]$ is the arclength; the normal $\widehat{\mathit{n}}(s,t)$ and binormal $\widehat{\mathit{b}}(s,t)$ are then constructed from the centreline tangent $\widehat{\mathit{t}}={\partial }_s{\mathit{y}}^c$. The surface of the swimmer is then defined by

$${\mathit{y}}^{\mathrm{bf}}(s,\theta ,t):={\mathit{y}}^c(s,t)+a(s)(cos(\theta )\widehat{\mathit{n}}(s,t)+sin(\theta )\widehat{\mathit{b}}(s,t)),$$

(55)

where the thickness function is defined to be

$$a\left(s\right):=a_0\left(1-0.95\left[e^{-200s^2}+e^{-200{(1-s)}^2}\right]\right).$$

(56)

Note that due to the use of the arclength parameterisation, the length of the swimmer is ensured to be constant throughout.

The trajectory of the swimmer is formulated as a pair of ordinary differential equations for ${\mathit{X}}^0\left(t\right)$, ${\mathit{b}}^1\left(t\right)$ and ${\mathit{b}}^2\left(t\right)$ (with ${\mathit{b}}^3\left(t\right):={\mathit{b}}^1\left(t\right)\wedge {\mathit{b}}^2\left(t\right)$),

$$\begin{array}{cc}\hfill {\dot{\mathit{X}}}^0\left(t\right)& =\mathit{U}\left(t\right),\hfill \end{array}$$

(57)

$$\begin{array}{cc}\hfill {\dot{\mathit{b}}}^{\alpha }\left(t\right)& =\mathbf{\Omega }\left(t\right)\wedge {\mathit{b}}^{\alpha }\left(t\right),\alpha =1,2.\hfill \end{array}$$

(58)

Equations (57) and (58) combined with the initial position and orientation are solved with the fourth-order Runge–Kutta method with timestep $\Delta t=0.02$, the inner problem of finding the instantaneous rates of change $\mathit{U}$ and $\mathsf{\Omega }$ through the solution of Equations (25) and (53).

All computations were carried out on a midspecification notebook computer with Intel^® Core${}^{\mathrm{TM}}$ i5-1135G7 CPU and 8 GB RAM.

2.4. Rate of Working

The instantaneous rate of working ${\mathit{W}}_{*}\left(t\right)$ is calculated as Equation (47) using the physical traction $*=\mathit{f}$ for the full boundary integral model and the adjusted density $*=\mathit{q}$ for the single-layer model. This quantity is then averaged over one beat cycle to yield, ${\overline{W}}_{*}:={\int }_0^1{W}_{*}\left(t\right)dt.$

3. Results

In this section, we present computational experiments with and without the inclusion of the double-layer potential for: (1) the resistance problems for a sphere with no-slip, free-slip, and partial-slip boundary conditions, comparing results against known analytical results; (2) the swimming problem for a squirming sphere, calculating the flow field and comparing the velocity and rate of working against analytical results; (3) undulatory swimming, examining the effect of slenderness on velocity and rate of working, and examining the effect of nearby plane boundaries on the flow field.

3.1. Resistance Problem with No-Slip, Free-Slip, and Partial-Slip Boundary Conditions

Figure 2 shows the discretisation and computed velocity fields for the resistance problems [N] and [F] around a translating sphere. While the flow fields are superficially similar, it can be seen that in the free-slip case, fluid to the left and right of the sphere is dragged downwards less strongly. Convergence results are given in Appendix F, Figure A1.

Results for the partial-slip resistance problem as a function of slip length $\lambda$ are given in Figure 3, showing very close correspondence to the exact solution over five orders of magnitude (red dots compared with solid blue line). These results also confirm that unjustified neglect of the double-layer potential for the resistance problem for nonrigid boundary conditions results in a significant error, which grows as the slip length increases (black circles compared with solid blue line) to around 24%.

3.2. Slip Velocity Squirmer

The flow field for the slip velocity squirmer is shown in Figure 4a; the source dipole nature of the exact solution is apparent. The walltime required for difference choices of nearest-neighbor discretisation (DoF and $h_q$) is shown in Figure 4b. The convergence of both the swimming velocity $\mathit{U}$ and rate of working $\mathcal{W}$ with quadrature spacing is shown in Appendix F, Figure A2. With 7200 DoF, a relative error in $\mathit{U}$ below 1% is found with all values of quadrature spacing used, whereas the relative error in $\mathcal{W}$ is rather larger, with values of 1–3% being achievable with $h_q⩽0.033$. Taking these results together, the results accurate to within a few percent error for both velocity and rate of working can be computed within around a minute without specialist computing hardware.

With the single-layer-only formulation in terms of the adjusted density $\mathit{q}$, the velocity was accurate to within $1.2$% of the exact value with the most refined discretisation (7200 DoF and $h_q=0.0167$). The calculation of rate of working in terms of $\mathit{q}$ yielded a value of $W_q=95.8$ as compared with the exact value $W=12\pi =37.7$, a 254% relative error.

3.3. Undulatory Swimmer

The results in this section utilise a discretisation with 324 DoF, and quadrature spacing $h_q=0.00245$ for the most slender swimmer and $h_q=0.00516$ for the least slender; computed instantaneous velocity and angular velocity were found to be within 3% of results with the force and quadrature discretisations four times finer (Appendix G,

Table A1. Convergence of the numerical solution of the swimming problem for three discretisation levels, for the most slender swimmer with $a_0=0.02$.

DoF	${\mathit{h}}_{\mathit{q}}$	U	$\mathsf{\Omega }$
324	$0.00938$	($-0.15495$, $0.35808$, $0.00000$)	($0.00106$, $0.00016$, $0.65005$)
1284	$0.00484$	($-0.15731$, $0.35873$, $0.00008$)	($0.00024$, $0.00014$, $0.64233$)
5112	$0.00245$	($-0.15777$, $0.36000$, $0.00007$)	($0.00011$, $0.00027$, $0.63180$)

and

Table A2. Convergence of the numerical solution of the swimming problem for three discretisation levels, for the least slender swimmer with $a_0=0.06$.

DoF	${\mathit{h}}_{\mathit{q}}$	U	$\mathsf{\Omega }$
324	$0.01975$	($-0.10700$, $0.38590$, $0.00023$)	($0.00093$, $0.00080$, $0.62631$)
1284	$0.01018$	($-0.10424$, $0.38616$, $0.00004$)	($0.00004$, $0.00008$, $0.62500$)
5112	$0.00516$	($-0.10402$, $0.38621$, $0.00004$)	($0.00001$, $0.00014$, $0.62060$)

). A discretisation of the swimmer showing the outward pointing normal vector $-\mathit{n}(\mathit{y},t)$ is given in Figure 5.

Figure 6 and Figure 7 show the velocity fields (white arrows) and stokeslet distributions (red arrows) for the slender and nonslender undulating swimmers at five time points during a beat cycle, along with the trajectory of the leading tip of the swimmer shown relative to the initial position for scale. The velocity field and force distributions are qualitatively very similar, despite the significant difference in slenderness, exhibiting the counter-rotating vortices known since the work of Gray [70]. However, the trajectories differ very significantly in the magnitude of progression per beat, a 68% reduction from $0.050$ body lengths per beat to $0.016$ body lengths per beat. Further visualisation of the three dimensional flow field via sections and instantaneous streamlines is given in Figure 8 (slender) and Figure 9 (nonslender). Each panel shows one of four approximately equally spaced intervals during the beat cycle, with instantaneous streamlines shown to highlight the fluid motion. The velocity fields along each plane ($(x_1,x_2)$, $(x_1,x_3)$, and $(x_2,x_3)$) with associated quiver plots showing direction are projected onto the sides and bottom of each panel. There is a strong qualitative similarity between the flow fields generated by the slender and nonslender swimmers at each step in all three plane sections and at each time point. The vertical sections in particular evidence radial streamlines, consistent with a stresslet-type flow field. To demonstrate application of the method to the common situation of inclusion of nearby solid boundaries (as implemented in Appendix E), Figure 10 and Figure 11 show flow fields with stationary solid plane boundaries measuring $2L\times 2L$ equidistant $0.25L$ above and below the swimmer, with in both cases vortical behaviour being evident both above and below the body due to boundary effects.

Finally, we consider the calculation of mean rate of working over a full beat cycle, as shown in Figure 12. From a methodological perspective, comparing the calculation with the physical potential ${\overline{\mathcal{W}}}_{\mathit{f}}$ to the single-layer-only calculation with the adjusted density ${\overline{\mathcal{W}}}_{\mathit{q}}$ yields a remarkably small discrepancy of only $1.6$–$2.0$%. Focussing on physical interpretation, as the slenderness ratio varies from $0.02$ to $0.06$, the mean rate of working varies almost linearly, increasing by 72% from $0.61$ to $1.05$. Combined with the dramatic reduction in swimming velocity, the slender swimmer is much more efficient than the nonslender swimmer.

4. Discussion

This manuscript reviewed the derivation of the method of regularised stokeslets, focussing on the formal justification, or otherwise, of eliminating the double-layer potential, before describing a method based on principal components analysis to enable its implementation within the meshless framework of both the original implementation and the nearest-neighbor method. Two test cases with known analytical solutions were then evaluated: the resistance problem for a translating sphere with no-slip, free-slip, and partial-slip boundary conditions, and the problem of finding the swimming velocity and rate of working for a slip velocity squirmer. By taking the double-layer potential into account, the resistance problem was solved with excellent accuracy and reasonable efficiency across five orders of magnitude in slip length; omitting the double-layer potential results in an error which grows with slip length, asymptoting to 24% in the free-slip limit. The squirming swimmer problem was solved with high accuracy for swimming velocity (better than 1% error within around a minute of walltime) although showed a slightly larger error of 1–3% in the calculation of rate of working. The double-layer term was confirmed to be critical in the calculation of rate of working, with the error being above 250% when the double layer was neglected.

Finally, the method was assessed in the evaluation of the swimming problem for an undulating worm-like organism, focussing on the effect of slenderness on progressive velocity and rate of working. The full boundary integral implementation and single-layer-only implementations differed in predicted values of rate of working by only $1.6$–$2.0$%, a value comparable to the regularisation and discretisation errors. Increasing from a peak radius of $0.02$ body lengths to $0.06$ body lengths, the progressive velocity was reduced by around 68%, and the mean rate of working increased by 72%, confirming the major propulsive advantage of slenderness for undulating swimmers. These results follow a similar trend to the classic boundary element study of Phan-Thien et al. [6], which showed that for rotating helical flagella attached to a spheroidal cell body, efficiency increased as the flagellum is made thinner, relative to the cell body. A similar effect has also been reported in the context of self-motile ribbons, with wider swimmers moving more slowly than narrow [71]. Given that the inclusion of the double layer incurs additional computational cost, determining when it may be eliminated is an important question. In stark contrast to the squirming swimmer case, the neglect of the double-layer potential resulted in a rather small discrepancy in mean rate of working, on the order of the numerical error. The latter finding provides evidence that, while not formally justified, the neglect of the double-layer potential may be acceptable for undulating swimmers such as worms and sperm. We postulate that the character of undulating motion—i.e., the relative rotation of adjacent segments may mean that the swimmer is behaving in a similar way to a collection of rigid bodies—and hence, arguments for the formal neglect of the double-layer potential may transfer approximately; the latter may form a topic for future work. This manuscript also did not assess the role of the regularisation parameter $ϵ$ and regularisation error; the effect of reducing $ϵ$, utilizing higher order blobs [53], or Richardson extrapolation [55] in problems involving the double-layer potential with the aim of reducing the regularisation error may also form a topic for future investigation.

This study focussed on extending the range of application of the regularised stokeslet method, so that its convenient implementation can be applied reliably to problems involving slip boundary conditions and moreover assessing the effect of neglecting the double-layer potential for undulating swimmers. The aim is not to supplant the classical boundary element method in accuracy or efficiency but rather to provide a relatively simple and accurate approach that avoids the need to generate a true surface mesh or evaluate singular quadratures. This ability may be particularly relevant when working with biological data in which geometry can only be extracted with limited precision, but rapid analysis is valuable. In addition to applications in microscale biological flow, the ability to accurately and rapidly model flows and bodies with partial-slip boundary conditions may be valuable in the field of nanofluidics [72]. On the scales of tens of nanometers, object surface roughness properties, material hydrophobicity, and a range of other phenomena entail that the no-slip condition is not appropriate and partial slip provides a more realistic description of the boundary–fluid interface behaviour [73]. Generally, the simulation tools used in the nanofluidics community are relatively costly particle-based methods such as molecular dynamics [74] and dissipative particle dynamics [75]. While inherently deterministic, Stokes flow is generally considered to give a valid mean description of fluid flow down to the 10 nm scale [76] and thereby provides potential advantages in efficiency.

In summary, the method of regularised stokeslets has proved a highly versatile, accessible, and effective method in microscale biological fluid dynamics. The double-layer potential can be included in this method without significant additional implementational complexity, enabling an accurate solution of velocity field, swimming velocity, and mean rate of working for problems involving surface slip; for undulating swimmers, our results provide confidence in the neglect of the double-layer potential.