Impact Statement
At typical scales, the velocity of a fluid adjacent to an impermeable solid surface is assumed to be equal to the velocity of the surface. However, in micro and nanoscale applications, a finite ‘velocity slip’ exists between the fluid and solid, which can have a significant influence on the overall flow behaviour. This article overviews an existing technique, based on the reciprocal theorem, for predicting the effects of low levels of slip in low Reynolds number conditions. The contribution of this work is (i) to generalise this technique to allow prediction of a broader range of flow properties; (ii) to apply the generalised method to a range of micro and nano-flow applications not considered before; (iii) to perform thorough numerical verification for the first time; and (iv) to derive general results to gain greater insight into low-slip flows at low Reynolds number. As examples of (iv) it is shown, within the limits of the assumptions, that (a) any particle with any distribution of positive slip length will have lower drag than the same particle with no slip, and (b) for a fixed head loss, flow rate through any straight channel of arbitrary cross-section and positive slip length will be greater than the same channel with no slip.
1. Introduction
In microscale and nanoscale fluid mechanics, the relative tangential motion of a fluid and a solid at an interface, known as velocity slip, is a familiar phenomenon. It is prominent, particularly, in microscale gas flows (Reference Arkilic, Schmidt and BreuerArkilic, Schmidt & Breuer 1997; Reference Gad-El-HakGad-El-Hak 1999; Reference Karniadakis and BeşkökKarniadakis & Beşkök 2002) and in liquid flows at the microscale and nanoscale (Reference Choi, Westin and BreuerChoi, Westin & Breuer 2003; Reference Holt, Park, Wang, Stadermann, Artyukhin, Grigoropoulos, Noy and BakajinHolt et al. 2006; Reference Lauga, Brenner and StoneLauga, Brenner & Stone 2007; Reference Falk, Sedlmeier, Joly, Netz and BocquetFalk et al. 2010; Reference Qin, Yuan, Zhao, Xie and LiuQin et al. 2011; Reference Nicholls, Borg, Lockerby and ReeseNicholls et al. 2012). Typically, these flows are at very low Reynolds number, due to their scale, and this is assumed to be the case throughout this article.
The standard approach to slip modelling is with a Maxwell or Navier slip boundary condition, for gases and liquids, respectively (Reference MaxwellMaxwell 1879; Reference Lockerby, Reese, Emerson and BarberLockerby et al. 2004; Reference Lauga, Brenner and StoneLauga et al. 2007). The two conditions are essentially equivalent, with both relating the velocity slip to the shear stress at the interface. In both cases, the degree of slip can be articulated using a slip length (described in detail later).
The focus of this paper is on low-slip flows, i.e. flows for which the slip length is small relative to a characteristic scale of the flow geometry. While in some applications small amounts of slip can simply be neglected, slip's influence on macroscopic properties such as drag, for which small changes can have significant implications, justifies its study at low levels.
In rarefied-gas dynamics and for common surfaces, the low-slip assumption is equivalent to the geometry being in the so-called ‘slip regime’, for which the adoption of a Maxwell slip condition with the conventional continuum equations is the accepted model. For air at standard-atmospheric pressure, this corresponds to devices and particulate on the scale of microns. Liquid slip in micro and nano geometries is less understood and harder to predict, and so the scale that can be classified as low slip is problem dependent.
Super-hydrophobic surfaces generate slip by trapping pockets of gas within micro or nano structures at the liquid–solid interface (Reference RothsteinRothstein 2010). This creates regions of very high slip (at the gas pockets) adjacent to regions of no slip (at the structures). What is often done is to calculate an ‘effective’ slip length for the heterogeneous surface (Reference Lauga and StoneLauga & Stone 2003; Reference Belyaev and VinogradovaBelyaev & Vinogradova 2010), which can be used with a Navier slip condition on a larger scale than the surface structures to predict their macroscopic impact. Cases where the effective slip length is small relatively to the macro geometry also fall into the category of low slip; and when at very low Reynolds number they are within the scope of this work.
An elegant technique, based on the reciprocal theorem, exists for deriving properties of low-slip Stokes flows in a general and convenient fashion (Reference Ramachandran and KhairRamachandran & Khair 2009; Reference StoneStone 2010); for example, for deriving expressions for the response of particle drag to small amounts of slip. The method has not received much dedicated attention, although it has been previously employed to derive expressions for Navier-slip flows around spheres (Reference StoneStone 2010), spheroids (Reference Masoud and StoneMasoud & Stone 2019) and Janus spheres (Reference Ramachandran and KhairRamachandran & Khair 2009); and non-Navier-slip flow through channels (Reference Michelin and LaugaMichelin & Lauga 2015). Importantly, the approach allows simple analytical results to be derived for low-slip flows, from existing no-slip solutions, for cases where full-slip solutions cannot be obtained or are extremely difficult to derive.
The motivation of this article is to: (i) bring this powerful technique into stronger focus, specifically for Maxwell/Navier slip flows in micro and nano-flow applications; (ii) to generalise the approach to allow prediction of a broader range of flow properties; (iii) to perform thorough numerical verification, which has not been done previously; (iv) to apply the generalised method to a range of micro and nano-flow applications not considered before; and (v) to derive general results to gain greater insight into Maxwell/Navier slip flows. On point (v), for example, we will show that, in low-slip Stokes flow, the drag on a single particle, of any shape, is well approximated by
\begin{equation} D \approx D_0-\frac{l}{ W \mu} \int_{S} {\tau}^2_{0}\, {\rm d}S, \end{equation}
where 
$D_0$ is the no-slip drag result, 
$l$ is the slip length, 
$W$ is the speed of the particle, 
$\mu$ is the viscosity and 
$\tau _0$ is the shear-stress magnitude over the particle surface (
$S$) from the no-slip solution to the same problem. We will also show that the pressure difference (
$\Delta p$) required to drive low-slip flow through straight channels, of any cross-sectional shape, can be calculated (note, very similarly) using
\begin{equation} \Delta p \approx \Delta p_0- \frac{ l {\mathscr{L}}}{Q \mu} \int_{{\mathscr{P}}} {\tau}_{0}^2\, {\rm d} {\mathscr{P}}, \end{equation}
where 
$\Delta p_0$ is the no-slip result, 
${\mathscr {L}}$ is the channel length, 
$Q$ is the volumetric flow rate and now the integral of 
$\tau _0^2$ is over the perimeter of the channel cross-section (
${\mathscr {P}}$).
The paper is structured as follows. In § 2 we will overview the basic theory, the spirit of which follows that contained in Appendix C of Reference Ramachandran and KhairRamachandran & Khair (2009) and in § 2.1 of Reference StoneStone (2010), but is developed here in a formulation allowing for arbitrary moments of the traction force to be estimated. In § 3 we present a number of new examples of its use, relevant to a range of microscale and nanoscale flow applications, including: predicting the mobility of particles with arbitrary slip-length distribution (§ 3.2), assessing the efficiency of a micro journal bearing (§ 3.3); predicting the speed of a self-propelled ‘squirmer’ with Navier-type slip (§ 3.4); and evaluating pressure loss in Navier/Maxwell-type slip flow through straight micro/nano channels (§ 3.5).
The approach adopted in this work also provides a means of numerically estimating the impact of slip in Stokes flow, purely from the post-processing of no-slip solutions or numerical calculations. This is discussed in § 4, alongside other general comments.
2. Theory
In this section we provide an overview, and generalisation, to the techniques presented previously (e.g. Reference Ramachandran and KhairRamachandran & Khair 2009; Reference StoneStone 2010) for deriving estimates of the effect of Navier/Maxwell slip on certain moments of the fluid traction force over a surface (e.g. the effect of slip on drag of a translating particle). In §§ 2.2 and 2.3, we generalise this approach to enable estimation of arbitrary moments of the traction force, which we use later in the article for application to problems not considered previously (in §§ 3.3, 3.4 and 3.5).
In this work, we restrict our attention to very low Reynolds number and steady-state flows, for which the governing equations are the steady Stokes equations
\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \quad \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\sigma} + \boldsymbol{f} = 0, \end{equation}
where 
$\boldsymbol {u}$ is the velocity, 
$\boldsymbol {f}$ is an applied body force and 
$\boldsymbol {\sigma }$ is the stress tensor, composed of both pressure and viscous stresses
\begin{equation} \boldsymbol{\sigma} = {-}p\boldsymbol{\mathsf{I}} + \mu ({\boldsymbol{\nabla} \boldsymbol{u}} + \boldsymbol{\nabla} \boldsymbol{u}^T ), \end{equation}
where 
$p$ is the pressure, 
$\mu$ is the dynamic viscosity, 
$\boldsymbol{\mathsf{I}}$ is the identity tensor and 
${T}$ denotes the transpose. Here we are concerned with solutions to the Stokes equations that satisfy a Navier-type slip condition at the boundary of the domain (
$S$)
\begin{equation} \boldsymbol{u}(\boldsymbol{r}) = \boldsymbol{U}(\boldsymbol{r}) + \frac{\ell(\boldsymbol{r})}{\mu} \boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\sigma}\boldsymbol{\cdot}(\boldsymbol{\mathsf{I}}-\boldsymbol{n}\boldsymbol{n})\quad \mathrm{for} \ \boldsymbol{r} \in S, \end{equation}
where 
$\boldsymbol {r}$ is the position vector, 
$\boldsymbol {u}$ is the velocity of the fluid at the fluid–solid interface, 
$\boldsymbol {U}$ is the velocity of the bounding surface itself (sometimes referred to as ‘the wall’), 
$\boldsymbol {n}$ is a surface normal directed into the fluid and 
$\ell (\boldsymbol {r})$ is a spatially varying slip length
\begin{equation} \ell = l \psi(\boldsymbol{r}), \end{equation}
where 
$l$ is the maximum slip length and 
$0\le \psi (\boldsymbol {r})\le 1$ is a non-dimensional function of position on the boundary. In most of the examples considered in this paper, surface properties are considered to be uniform, 
$\psi =1$, so that 
$\ell \equiv l$. Note, (2.3) is general enough to represent velocity conditions at open boundaries (inlet/outlets) and zero-disturbance far-field conditions, by setting 
$\phi =0$, and specifying 
$\boldsymbol {U}$, accordingly.
The vector 
$\boldsymbol {n}\boldsymbol {\cdot }\boldsymbol {\sigma }$ appearing in (2.3) is the surface traction force, and its inner product with the tensor 
$(\boldsymbol {I}-\boldsymbol {n}\boldsymbol {n})$ removes its surface-normal component, yielding a surface shear-stress vector, 
$\boldsymbol {\tau }$ (which the slip velocity is proportional to). Note, when resolved in the surface-normal direction, (2.3) corresponds to the impermeability condition.
The key parameter of the current work is the non-dimensional maximum slip length (in cases where the distinction is unimportant, we will sometimes refer to this as the ‘non-dimensional slip length’, or just ‘slip length’)
\begin{equation} \xi = \frac{l}{L}, \end{equation}
which expresses the degree of slip relative to a characteristic length scale of the flow in question, 
$L$. In this paper, our attention is restricted to situations where the level of slip is low: 
$\xi \ll 1$.
Equation (2.3) is equivalent to Maxwell's slip boundary condition for isothermal rarefied flows (Reference MaxwellMaxwell 1879; Reference Lockerby, Reese, Emerson and BarberLockerby et al. 2004); where 
$\ell =\lambda (2-\alpha )/\alpha$, 
$\lambda$ is the mean free path and 
$\alpha$ is the accommodation coefficient. For most practical surfaces 
$\alpha \approx 1$, which makes the non-dimensional slip length approximately equal to the Knudsen number
\begin{equation} \xi\approx \frac{\lambda}{L} = \mathit{K\hspace{-0.03cm}n}. \end{equation}
Importantly, in adopting Maxwell's slip model it is already implied that 
$\xi \approx \mathit {K\hspace {-0.03cm}n} \ll 1$. In other words, the assumption of 
$\xi \ll 1$ is consistent with the study of low-speed rarefied-gas flows in the slip regime, for which (2.2) and (2.3) are valid. For higher degrees of rarefaction, Knudsen layers and other gas-kinetic phenomena make their application unsuitable (Reference CercignaniCercignani 1969; Reference SoneSone 2002; Reference Lockerby and ReeseLockerby & Reese 2008; Reference TorrilhonTorrilhon 2016).
2.1 Series expansion for Stokes flow with slip
As in Reference Ramachandran and KhairRamachandran & Khair (2009) and Reference StoneStone (2010), we start by expanding the slip Stokes-flow solution in an infinite power series using the non-dimensional slip length as a small parameter, 
$\xi \ll 1$
\begin{gather} \boldsymbol{u} = \boldsymbol{u}_{0} + \xi\boldsymbol{u}_1 + \xi^2\boldsymbol{u}_2 + \cdots = \sum_{k = 0}^{\infty} \boldsymbol{u}_{k} \xi^k, \end{gather}
\begin{gather}\boldsymbol{\sigma} = \boldsymbol{\sigma}_0 + \xi\boldsymbol{\sigma}_1 + \xi^2\boldsymbol{\sigma}_2 + \cdots = \sum_{k = 0}^{\infty} \boldsymbol{\sigma}_{k} \xi^{k}, \end{gather}
where 
$\boldsymbol {u}_k,\boldsymbol {\sigma }_k$ is a Stokes-flow solution associated with the 
$k$th order of the expansion. Substituting (2.7) and (2.8) into the Navier slip condition (2.3), and equating orders of 
$\xi$, yields the boundary conditions for the successive Stokes-flow solutions in the expansion. The velocity field of the first solution (
$\boldsymbol {u}_{0}$) satisfies
\begin{equation}
\boldsymbol{u}_{0}(\boldsymbol{r}) =
\boldsymbol{U}(\boldsymbol{r}) \quad \mathrm{for} \
\boldsymbol{r} \in S ,
\end{equation}and therefore corresponds to the no-slip solution. Subsequent solutions in the series satisfy a velocity slip proportional to the shear-stress vector of the previous solution
\begin{equation} \begin{aligned}
\boldsymbol{u}_{1}(\boldsymbol{r}) & =
\frac{L}{\mu}\psi\,\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\sigma}_{0}\boldsymbol{\cdot}(\boldsymbol{\mathsf{I}}-\boldsymbol{n}\boldsymbol{n}),
\\ \boldsymbol{u}_{2}(\boldsymbol{r}) & =
\frac{L}{\mu}\psi\,\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\sigma}_{1}\boldsymbol{\cdot}(\boldsymbol{\mathsf{I}}-\boldsymbol{n}\boldsymbol{n}),
\\ & \vdots \\ \boldsymbol{u}_{k}(\boldsymbol{r}) & =
\frac{L}{\mu}\psi\,\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\sigma}_{k-1}\boldsymbol{\cdot}(\boldsymbol{\mathsf{I}}-\boldsymbol{n}\boldsymbol{n})
\quad \mathrm{for} \ \boldsymbol{r} \in S. \end{aligned}
\end{equation}2.2 Moments of the traction force
The main focus of the method is on calculating surface moments of the traction force. The novelty of the development here is in maintaining the generality of the surface moment
\begin{equation} M = \int_{S} \boldsymbol{g}\boldsymbol{\cdot} \boldsymbol{\sigma}\boldsymbol{\cdot} \boldsymbol{n} \,{\rm d}S, \end{equation}
where 
$\boldsymbol {g}(\boldsymbol {r})$ defines the moment in question. For example, if 
$S_p$ is the solid boundary of a particle, then
\begin{equation} \boldsymbol{g}(\boldsymbol{r}) = \left\{\begin{array}{ll} \boldsymbol{i}_x & \quad \mathrm{for}\ \boldsymbol{r} \in S_p, \\ 0 & \quad \mathrm{for}\ \boldsymbol{r} \in S / S_p, \end{array} \right. \end{equation}
produces a moment corresponding to the 
$x$-component of the total hydrodynamic force acting on the particle. If, as another example, the function is of the form
\begin{equation} \boldsymbol{g}(\boldsymbol{r}) = \left\{\begin{array}{ll} \boldsymbol{r}\times\boldsymbol{i}_x & \quad \mathrm{for}\ \boldsymbol{r} \in S_p, \\ 0 & \quad \mathrm{for}\ \boldsymbol{r} \in S / S_p, \end{array} \right. \end{equation}
the moment corresponds to the hydrodynamic torque on the particle about the 
$x$ axis.
Substitution of (2.8) into (2.11), yields a series representation of an arbitrary traction-force moment
\begin{equation} M = M_{0} + \xi M_1 + \xi^2 M_2 + \cdots = \sum_{k = 0}^{\infty} M_{k} \xi^k, \end{equation}where
\begin{equation} M_k = \int_{S} \boldsymbol{g} \boldsymbol{\cdot} \boldsymbol{\sigma}_k \boldsymbol{\cdot} \boldsymbol{n} \, {\rm d}S. \end{equation}For low-slip flows, the first-order approximation is a good one, i.e.
\begin{equation} M \approx M_{0} + \xi M_1, \end{equation}
where 
$M_1$ is the first-order slip-correction coefficient. The motivation of the method is to find a general and convenient way of obtaining 
$M_1$.
2.3 Finding the slip-correction coefficient, 
$M_1$
Let 
$\boldsymbol {u}_0'$, 
$\boldsymbol {\sigma }_0'$ be a no-slip Stokes solution satisfying the boundary condition
\begin{equation} \boldsymbol{u}'_{0}(\boldsymbol{r}) = \boldsymbol{g}(\boldsymbol{r}) \quad \mathrm{for} \ \boldsymbol{r} \in S, \end{equation}
with a body force 
$\boldsymbol {f}'=0$. We will refer to this as the conjugate solution. Substituting (2.17) into (2.15) gives
\begin{equation} M_1 = \int_{S} \boldsymbol{u}'_{0} \boldsymbol{\cdot} \boldsymbol{\sigma}_1 \boldsymbol{\cdot} \boldsymbol{n} \, {\rm d}S. \end{equation}From the reciprocal theorem, (2.18) can be written
\begin{equation} M_1 = \int_{S} \boldsymbol{u}_1 \boldsymbol{\cdot}\boldsymbol{\sigma}_0' \boldsymbol{\cdot} \boldsymbol{n} \, {\rm d}S - \int_{V} \boldsymbol{f} \boldsymbol{\cdot} \boldsymbol{u}'_{0} \, {\rm d}V. \end{equation}
Now, upon substituting the boundary conditions for the first-order solution for velocity, 
$\boldsymbol {u}_1$, from (2.10), we obtain
\begin{equation} M_1 = \frac{L}{ \mu} \int_{S} \psi\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\sigma}_{0}\boldsymbol{\cdot}(\boldsymbol{\mathsf{I}}-\boldsymbol{n}\boldsymbol{n}) \boldsymbol{\cdot}\boldsymbol{\sigma}_0' \boldsymbol{\cdot} \boldsymbol{n} \, {\rm d}S - \int_{V} \boldsymbol{f} \boldsymbol{\cdot} \boldsymbol{u}'_{0} \, {\rm d}V. \end{equation}
Given that 
$(\boldsymbol{\mathsf{I}}-\boldsymbol {n}\boldsymbol {n}) \boldsymbol {\cdot }\boldsymbol {\sigma }_0'=(\boldsymbol{\mathsf{I}}-\boldsymbol {n}\boldsymbol {n}) \boldsymbol {\cdot }\boldsymbol {\sigma }_0' \boldsymbol {\cdot }(\boldsymbol{\mathsf{I}}-\boldsymbol {n}\boldsymbol {n})$, this becomes
\begin{equation} M_1 = \frac{L}{\mu} \int_{S} \psi \boldsymbol{\tau}_{0}\boldsymbol{\cdot}\boldsymbol{\tau}'_{0}\, {\rm d}S- \int_{V} \boldsymbol{f} \boldsymbol{\cdot} \boldsymbol{u}'_{0} \, {\rm d}V, \end{equation}
where 
$\boldsymbol {\tau } = \boldsymbol {n}\boldsymbol {\cdot }\boldsymbol {\sigma }\boldsymbol {\cdot }(\boldsymbol{\mathsf{I}}-\boldsymbol {n}\boldsymbol {n})$ is the tangential shear-stress vector. Importantly, the right-hand side of (2.21) is purely in terms of no-slip solutions.
In all of the examples we will consider in this paper, there is no applied body force, and so we will omit the second term in (2.21), and work with the simpler and normalised expression for the first-order slip-correction coefficient
\begin{equation} \hat{M}_1 = \frac{L}{M_0 \mu} \int_{S} \psi \boldsymbol{\tau}_{0}\boldsymbol{\cdot}\boldsymbol{\tau}'_{0}\, {\rm d}S, \end{equation}
where 
$\hat {M}_1=M_1/M_0$. From hereon, hats denote a dimensionless value normalised with its corresponding no-slip quantity.
Equation (2.22) represents a new, succinct and general expression for evaluating the first-order slip correction for arbitrary moments of the surface traction force.
2.4 The resistive moment, ${\mathscr {R}}$
We now introduce an important special case, which simplifies (2.22). For many applications, the function 
$\boldsymbol {g}$ that defines the moment of the traction force has the same dependence on position as the velocity of the wall, but the opposite sign
\begin{equation} \boldsymbol{U}(\boldsymbol{r}) = {-}k \boldsymbol{g}(\boldsymbol{r}), \end{equation}
where 
$k$ is a constant of proportionality. For example, we might want to calculate the 
$x$-component of drag on a particle (
$\boldsymbol {g}=\boldsymbol {i}_x$) due to its translation in the opposite direction (
$\boldsymbol {U}=-U \boldsymbol {i}_x$). In this case, 
$k=U$, where 
$U$ is the particle speed. Alternatively, we might want to calculate the torque around the 
$x$-axis of a particle (
$\boldsymbol {g}=\boldsymbol {r}\times \boldsymbol {i}_x$) that rotates about the 
$x$-axis in the opposite sense (
$\boldsymbol {U}=-\omega \boldsymbol {r}\times \boldsymbol {i}_x$). Here, 
$k=\omega$, where 
$\omega$ is the magnitude of the angular velocity.
We refer to this particular moment as the resistive moment, using the symbol 
${\mathscr {R}}$ to distinguish it from the general case. The resistive moment can be expanded, to first order in slip length, as previously
\begin{equation} \frac{{\mathscr{R}}}{{\mathscr{R}}_0} = 1 + \hat{{\mathscr{R}}}_1 \xi + {\mathscr{O}}(\xi^2), \end{equation}
where 
$\hat {{\mathscr {R}}}_1={\mathscr {R}}_1/{\mathscr {R}}_0$, and the corresponding no-slip moment is given by
\begin{equation} {\mathscr{R}}_0 = {-}\frac{1}{k} \int_{S} \boldsymbol{u}_{0} \boldsymbol{\cdot} \boldsymbol{\sigma}_0 \boldsymbol{\cdot} \boldsymbol{n} \, {\rm d}S, \end{equation}
which in the absence of external forces is guaranteed to be positive. To find 
${\mathscr {R}}_1$, we substitute (2.9) and (2.17) into (2.23) to give 
$\boldsymbol {u}_0=-k \boldsymbol {u}'_0$, from which it follows that
\begin{equation} \boldsymbol{\tau}_0 = {-}k \boldsymbol{\tau}'_0. \end{equation}Finally, substitution of (2.26) into the general expression (2.22) gives the first-order slip-correction coefficient for the resistive moment
\begin{equation} \hat{{{\mathscr{R}}}}_{1} = {-}\frac{L}{{\mathscr{R}}_0 k \mu} \int_{S} \psi {\tau}^2_{0}\, {\rm d}S, \end{equation}
where 
$\tau _0$ is the magnitude of the tangential shear stress. Note, for the resistive moment there is no additional conjugate no-slip solution.
It is significant that this is necessarily negative for any distribution of positive slip length. Since 
${\mathscr {R}}_0>0$, this tells us that a small amount of slip on any geometry will reduce the resistive moment of the traction force. In short, and for example, low levels of slip will reduce the drag on any translating geometry in Stokes flow (in the absence of applied body forces). Similarly, any distribution of slip length (provided it is small) will reduce the retarding torque on a rotating particle.
3. Application examples and analytical solutions
The remainder of this paper is dedicated to applying the expressions derived above to different flow problems, in order to derive analytical expressions for first-order slip corrections. Generally, these are cases for which full analytical slip solutions either do not exist or are extremely involved to evaluate.
To begin with, however, for the purposes of exposition, we apply the method to a case that has a well-known slip solution, and for which the general expression (2.22) derived in § 2.3 is not required; namely, slip-flow drag of a translating sphere (§ 3.1). In Appendix A, we also derive first-order slip corrections to the drag on prolate and oblate spheroids, first done by Reference Masoud and StoneMasoud & Stone (2019), but which we present in a different form and numerically verify for the first time.
In §§ 3.2–3.5, we consider Navier-slip applications that have not been considered previously and that, for all but one, require the general expression for the traction moment derived in § 2.3: a sphere with arbitrary slip-length distribution (§ 3.2); a micro journal bearing (§ 3.3); a self-propelled particle (§ 3.4); and flow through long, straight micro/nano channels (§ 3.5).
3.1 A sphere
As in Reference StoneStone (2010), we start by verifying the expressions derived in § 2 for a problem that has a simple and well-established analytical treatment for slip flow; namely, translational and rotational motion of a single sphere in free space; due to Reference BassetBasset (1888). Basset's solutions for drag and retarding torque, derived for a uniform slip length (
$\psi =1$), are
\begin{equation} D = D_0\left( \frac{1 + 2 \xi}{1 + 3\xi}\right) \quad \mathrm{and} \quad T = T_0 \left(\frac{1}{1 + 3\xi}\right), \end{equation}
where 
$\xi =l/R$ is the non-dimensional slip length, 
$R$ is the sphere radius, 
$D_0=6\pi \mu W R$ is the no-slip result for drag on a translating sphere with velocity 
$W$ and 
$T_0=8\pi \mu \omega R^3$ is the no-slip result for retarding torque on a rotating sphere with angular velocity 
$\omega$; see figure 1. Expanding Basset's expressions in a Taylor series gives the drag and torque to first order in slip length
\begin{equation} \hat{D} = 1-\xi + {\mathscr{O}}(\xi^2) \quad \mathrm{and} \quad \hat{T} = 1-3 \xi + {\mathscr{O}}(\xi^2), \end{equation}
where 
$\hat {D}=D/D_0$ and 
$\hat {T}=T/T_0$. The first-order slip-correction coefficients are, therefore, 
$\hat {D_1}=-1$ and 
$\hat {T_1}=-3$.
A translating (
$W$) and rotating (
$\omega$) sphere, of radius 
$R$, in spherical coordinates.
A side note: the first-order slip-correction coefficient for drag (
$\hat {D_1}=-1$) was first obtained from kinetic theory for dilute gas flows by Reference EpsteinEpstein (1924), who also demonstrated that Basset's full-slip solution (3.1) was only valid to this order; see Reference Happel and BrennerHappel & Brenner (1983) for more discussion.
We now demonstrate how to obtain the first-order coefficients in (3.3) directly from the corresponding no-slip solutions. For the case of a translating sphere, and in spherical polar coordinates, the distribution of wall shear-stress magnitude in the no-slip solution is
\begin{equation} \tau_0 = \frac{3 \mu W \sin \theta }{2 R}, \end{equation}
where 
$\theta$ is the polar angle (see figure 1). In this case, drag force is the resistive moment, and so the first-order slip-correction coefficient can be obtained directly from (2.27) (with 
${\mathscr {R}}_1= D_1$, 
$L=R$, 
$\psi =1$ and 
$k = W$)
\begin{equation} \hat{D}_1 = {-}\frac{R}{D_0 W \mu} \int_{S} {\tau}^2_{0}\, {\rm d}S = {-} \frac{9\pi \mu W R}{2 D_0} \int^\pi_0 \sin^3\theta\, {\rm d}\theta = {-}1, \end{equation}which agrees with Basset and Epstein's solutions.
For the rotating sphere, the surface shear-stress magnitude with no slip is
\begin{equation} \tau_0 = 3 \mu \omega \sin \theta. \end{equation}
For this case, retarding torque is the resistive moment, and so, again, the first-order slip-correction coefficient is obtained directly from (2.27) (but with 
${\mathscr {R}}_1= T_1$, 
$L=R$, 
$\psi =1$ and 
$k=\omega$)
\begin{equation} \hat{T}_1 = {-}\frac{R}{T_0 \omega \mu} \int_{S} {\tau}^2_{0}\, {\rm d}S = {-} \frac{18\pi \mu \omega R^3}{ T_0} \int^\pi_0 \sin^3\theta\, {\rm d}\theta = {-}3, \end{equation}which, again, agrees with the result of Basset.
3.2 A sphere with varying slip length
We now choose an example for which full-slip solutions, like those due to Basset, do not exist. Consider a rigid sphere, as in figure 1 and § 3.1, but with a non-constant slip length
\begin{equation} \ell = l \psi(\theta,\phi), \end{equation}
where 
$\psi$ is some arbitrary surface function and 
$l$ is the maximum slip length over the surface. The first-order slip-correction coefficient for drag force due to translation comes directly from (2.27) for the resistive moment (with 
${\mathscr {R}}_1= D_1$, 
$L=R$ and 
$k = W$)
\begin{equation} \hat{D}_1 = {-} \frac{R}{D_0 W \mu} \int_{S} \psi {\tau}^2_{0} \, {\rm d}S, \end{equation}
where, as a reminder, 
$\tau _0$ is the surface shear-stress magnitude of a sphere in translation with no slip (3.3). We can express 
$\tau _0^2$ as follows:
\begin{equation} \tau_0^2 = \frac{3\sqrt{\pi} \mu^2 W^2 }{R^2}\left( Y^0_{0}-\frac{1}{\sqrt{5}}Y^0_{2} \right), \end{equation}
where 
$Y^0_0={1}/{2\sqrt {\pi }}$ and 
$Y^0_2= \tfrac {1}{4} \sqrt {{5}/{\pi }} (3 \cos ^2(\theta )-1)$ are spherical harmonics (essentially Legendre polynomials). On substitution into (3.8), and then (2.24), we obtain an expression for the drag on a variable-slip-length sphere
\begin{equation} \hat{D} = 1 -\bar{\psi}\xi + \frac{2\sqrt{5 \pi} }{ 5} \psi^0_{2} \xi + {\mathscr{O}}(\xi^2), \end{equation}
where 
$\bar {\psi }\xi$ is the average slip length, 
$\bar {\psi }= \int _S \psi \, {\rm d}S/(4\pi R^2)$ and 
$\psi ^0_2 =\int _S \psi Y_2^0 \, {\rm d}S/(4\pi R^2)$.
This is a surprising result: for the same average slip length, only variations in slip length that are of the form of the second zonal spherical harmonic (
$Y_2^0$) influence the drag in low-slip flow; we can say this because of the orthogonality of spherical harmonics.
Let us take, for example, a type of ‘Janus particle’: a sphere having a constant slip length on one hemisphere and no slip on the other. Irrespective of the sphere's orientation, the first-order slip-correction coefficient is half that of a sphere with a uniform slip length, since in all orientations 
$\bar {\psi }=\frac {1}{2}$ and 
$\psi ^0_2=0$. In short, in low-slip flow, this type of particle has equal drag in all directions. Note, this case, as well as the case of asymmetric regions of uniform slip length, were considered previously by Reference Ramachandran and KhairRamachandran & Khair (2009) using the approach.
A new example is shown in figure 2, where the slip (or no-slip) region is a central band, again, covering half the surface area. For the case of a central band of constant slip length, figure 2(a), 
$\bar {\psi }=\frac {1}{2}$ and 
$\psi ^0_2=-3 \sqrt {({5}/{\pi }) }/32$. For the case of polar regions of constant slip length, figure 2(b), 
$\psi ^0_0=\frac {1}{2}$ and 
$\psi ^0_2=3 \sqrt {({5}/{\pi }) }/32$. For translation along the polar axis, the drag to first order in slip length is therefore given by
\begin{equation} \hat{F} = 1 -\frac{1}{2}\xi + \frac{3\upsilon}{16}\xi + {\mathscr{O}}(\xi^2), \end{equation}
where 
$\upsilon =-1$ and 
$+1$ for the central-band and polar-regions case, respectively. In both cases, translation in a perpendicular direction to that shown in figure 2 results in a drag force given by (3.11) with 
$\upsilon =0$ (
$\bar {\psi }=\frac {1}{2}$ and 
$\psi ^0_2=0$). As such, for the central-band case, the translation direction with minimum drag is in the direction shown in figure 2(a) (i.e. perpendicular to the band). However, for the sphere with polar regions of slip, the translation direction with lowest drag is perpendicular to that indicated in figure 2(b).
A sphere with: (a) a constant slip length on a central band (dark grey) and no slip on the polar caps (light grey); and (b) vice versa.
An identical analysis can be repeated for the retarding torque due to rotation about the polar axis (see figure 2). In the general case
\begin{equation} \hat{T} = 1 -3 \bar{\psi}\xi + \frac{6\sqrt{5 \pi} }{ 5} \psi^0_{2} \xi + {\mathscr{O}}(\xi^2). \end{equation}As with the drag-force case, the retarding torque in low-slip flow is only affected by the average slip length and slip-length variations in the form of the second zonal spherical harmonic. For the examples illustrated in figure 2
\begin{equation} \hat{T} = 1 -\frac{3}{2}\xi + \frac{9 \upsilon}{16}\xi + {\mathscr{O}}(\xi^2). \end{equation}To verify these analytical results, we perform numerical simulations of Stokes flow with a Navier slip condition, using the method of fundamental solutions. Appendix B gives full details of the numerical methodology and the parameters used.
Tables 1 and 2 compare the numerical simulations with the analytical results for the cases illustrated in figure 2. Other than to verify the derivations for low-slip conditions, the purpose of the comparison is to illustrate the extent to which predictions from the low-slip assumption diverge from the full-slip numerical simulation with increasing 
$\xi$. Note, what is tabulated is the calculation of torque/drag change divided by the slip length (e.g. 
$(\hat {F}-1)/\xi$), which is done so that the numerical results converge to a finite value (the first-order slip-correction coefficient) as 
$\xi \to 0$. As such, because of this division, the observed first-order error of the analytical predictions compared with the numerical results is indicative of the second-order error in 
$\hat {F}$ and 
$\hat {T}$; as expected from (3.11) and (3.13).
Comparison of numerical and analytical predictions of the first-order slip-correction coefficient for drag (
$\hat {D}_1$) via calculation of 
$(\hat {F}-1)/\xi$; for the configuration shown in figure 2.
Comparison of numerical and analytical predictions of the first-order slip-correction coefficient for torque (
$\hat {T}_1$) via calculation of 
$(\hat {T}-1)/\xi$; for the configuration shown in figure 2.
3.3 Journal bearing
A plain journal bearing is shown in figure 3, having a rotating inner shaft of radius 
$r$, with angular velocity 
$\omega$, and an axial offset 
$a$ from a stationary containing sleeve of radius 
$R$. Here, we consider the simplest problem of a single fluid phase between the shaft and sleeve.
Schematic of a journal bearing.
3.3.1 No-slip lubrication analysis
In this subsection, we overview the key assumptions and results of the standard no-slip lubrication analysis, closely following the exposition in Reference ShermanSherman (1990).
The radial clearance of the journal bearing is defined as 
$C=R-r$, and it is assumed that 
$C\ll R$ such that the curvature of the streamlines can be neglected and the local clearance is well approximated by
\begin{equation} H(\theta) = C(1- \eta \cos \theta ), \end{equation}
where 
$\eta =a/C$ is the shaft eccentricity (see figure 3). The no-slip solution predicts the fluid shear stress acting on the shaft
\begin{equation} \tau_{0, {shaft}} = \frac{2 \mu R \omega (2 (\eta ^2 + 2) \eta \cos \theta -5 \eta ^2-1)}{C (\eta ^2 + 2) (\eta \cos \theta -1)^2}, \end{equation}and the sleeve
\begin{equation} \tau_{0, {sleeve}} = \frac{2 \mu R \omega ((\eta ^2 + 2) \eta \cos \theta -4 \eta ^2 + 1)}{C (\eta ^2 + 2) (\eta \cos \theta -1)^2}. \end{equation}The resistive torque (per unit length) acting on the shaft is obtained by integration
\begin{equation} T_0 = {-}R^2 \int^{2\pi}_0 \tau_{0,{shaft}} \, {\rm d} \theta = \frac{4 \pi (2 \eta ^2 + 1) \mu R^3 \omega }{C \sqrt{1-\eta ^2} (\eta ^2 + 2)}. \end{equation}
Note, since 
$C\ll R$ in this lubrication analysis, 
$r\approx R$.
Another important moment of the traction force for the journal bearing is the lift per unit length generated on the shaft (
$\varLambda$); i.e. the net force in the direction of 
$\theta =\pi /2$ (see figure 3)
\begin{equation} \varLambda_0 = \frac{12 \pi \eta \mu R^3 \omega }{C^2 \sqrt{1-\eta ^2} (\eta ^2 + 2)}; \end{equation}
(note, there is a factor of 
$R$ missing in Reference ShermanSherman 1990).
A figure of merit for the journal bearing is given by a dimensionless ratio of the lift to the torque (
${\mathscr {M}}=R \varLambda /T$), providing a measure of the cost of producing lift. The figure of merit in no-slip conditions is
\begin{equation} {\mathscr{M}}_0 = \frac{R \varLambda_0}{T_0} = \frac{3 R \eta }{ C(1 + 2\eta ^2)}, \end{equation}
which shows that the greater the eccentricity (
$\eta$) the greater the efficiency of the bearing (by this measure).
3.3.2 First-order slip corrections
Slip flow in journal bearings has been studied in a variety of contexts (Reference Singh, Rao and MajumdarSingh, Rao & Majumdar 1984; Reference Shahdhaar, Yadawad, Khamari and BeheraShahdhaar et al. 2020; Reference Arif, Kango and ShuklaArif, Kango & Shukla 2022), and is normally modelled using a modified Reynolds equation; i.e. using a lubrication analysis similar to the above, but with slip flow (Reference Li, Chu and ChenLi, Chu & Chen 2006; Reference Zhang, Zhou and MengZhang, Zhou & Meng 2011; Reference Shahdhaar, Yadawad, Khamari and BeheraShahdhaar et al. 2020; Reference Arif, Kango and ShuklaArif et al. 2022). To the author's knowledge, no analytical solution to the slip-modified lubrication model has been presented for the journal bearing.
The first-order slip-correction coefficient to the retarding torque (which is the resistive moment) is obtained directly from (2.27) (with 
${\mathscr {R}}= T$, 
$L=C(1-\eta )$, 
$\psi =1$ and 
$k = \omega$)
\begin{equation} \hat{T}_1 = {-}\frac{C(1-\eta)}{T_0 \omega \mu } \int_0^{2\pi} (\tau_{0, {shaft}}^2 + \tau_{0, {sleeve}}^2 )R\, \mathrm{d}\theta = {-}\frac{-8 \eta ^4 + 22 \eta ^2 + 4}{(\eta + 1) (\eta ^2 + 2) (2 \eta ^2 + 1)}, \end{equation}
where the minimum clearance, 
$C(1-\eta )$, is taken as the characteristic scale of the bearing. Note, integration is over both surfaces, not just the shaft. This analytical solution for low-slip flow (which is equivalent to ‘slip-flow’ conditions in gas bearings), tells us that slip will reduce the retarding torque for all values of eccentricity (
$0\le \eta <1$).
The first-order slip correction for the lift force (which is not the resistive moment) requires evaluation of the more general expression for the slip-correction coefficient (2.22), and requires a conjugate solution: a no-slip solution to shaft translation, in a direction that we wish to evaluate the lift (in the direction 
$\theta =\pi /2$). The shear stress on both shaft and sleeve from a unit translational velocity of the shaft in a direction parallel to 
$\theta =\pi /2$ is
\begin{equation} \tau_0' = {-}\frac{6 \mu R ((\eta ^2 + 2) \cos \theta -3 \eta )}{C^2 (\eta ^2 + 2) (\eta \cos \theta -1)^2}. \end{equation}
Substitution into (2.22) (with 
$M= \varLambda$, 
$L=C(1-\eta )$ and 
$\psi =1$) gives
\begin{equation} \hat{\varLambda}_1 = \frac{C(1-\eta)}{\varLambda_0 \mu } \int_0^{2\pi} \tau'_0(\tau_{0, {shaft}} + \tau_{0, {sleeve}} ) R\, \mathrm{d}\theta = \frac{6 (\eta ^2-2)}{(\eta + 1) (\eta ^2 + 2)}. \end{equation}
Quick inspection of (3.22) reveals that, as for retarding torque, the slip-correction coefficient for lift is negative for all values of eccentricity (
$0\le \eta <1$). In other words, small amounts of slip will always reduce the lift generated by the journal bearing.
The figure of merit can be expanded as follows:
\begin{equation} \hat{{\mathscr{M}}} = \frac{\hat{\varLambda}}{\hat{T}} = \frac{1 + \hat{\varLambda_1}\xi + {\mathscr{O}}(\xi^2)}{1 + \hat{T}_1\xi + {\mathscr{O}}(\xi^2)} = 1 + \hat{{\mathscr{M}}}_1 \xi + {\mathscr{O}}(\xi^2), \end{equation}where
\begin{equation} \hat{{\mathscr{M}}}_1 = \hat{L}_1-\hat{T}_1 = \frac{4 (\eta -1)}{2 \eta ^2 + 1}. \end{equation}The immediate observation is that (3.24) is necessarily negative: small amounts of slip will always lower the bearing's figure of merit. In other words, for low-slip flows, slip reduces the lift force proportionally more than it reduces the resistive torque. However, the impact of slip's negative effect on the figure of merit is reduced for greater eccentricities.
3.3.3 Numerical verification
The governing Reynolds equation for the fluid pressure in the bearing, assuming a uniform slip length 
$\ell$ on both shaft and sleeve, is
\begin{equation} \frac{{\rm d}}{{\rm d}\theta} \left( H^2(H + 6\ell) \frac{{\rm d}p}{{\rm d}\theta} \right) = 6 \mu \omega R^2 \frac{{\rm d} H}{{\rm d} \theta}, \end{equation}
where 
$p$ is the pressure, and which upon integration gives
\begin{equation} \frac{{\rm d}p}{{\rm d}\theta} = 6 \mu \omega R^2 \frac{ H}{H^2(H + 6\ell)} + \frac{A}{H^2(H + 6\ell)}. \end{equation}
The constant of integration, 
$A$, can be found by numerically integrating (3.26) with the condition that the pressure be continuous (
$\int _0^{2\pi } ({{\rm d}p}/{{\rm d}\theta }) \,{\rm d}\theta =0$). Subsequent integration of (3.26) provides the pressure distribution in the bearing. Along with the shear-stress distribution on the shaft
\begin{equation} \tau = {-}\frac{H}{2R}\frac{{\rm d}p}{{\rm d}\theta}-\frac{\mu \omega R}{2\ell + H}, \end{equation}the lift and retarding torque on the shaft can be obtained.
For each of the moments, 
$T$ and 
$\varLambda$, and the figure of merit 
${\mathscr {M}}$, the first-order slip-correction coefficients are estimated from numerical calculations by 
$\hat {M}_1\approx {(M-M_0)}/({M_0 \xi })$, where 
$\xi =l/L$ and 
$L=C(1-\eta )$. Figure 4 compares these numerical results with the analytical results derived above; as expected, as the slip length is reduced, they converge. For larger slip lengths (relative to the minimum clearance), the numerical results differ from the analytical results, but the qualitative variation with changing eccentricity remains similar.
Comparison of numerical and analytical predictions for the first-order slip-correction coefficients. Analytical expressions for 
$\hat {T}_1$, 
$\hat {\varLambda }_1$ and 
$\hat {{\mathscr {M}}}_1$ (solid red lines) from (3.20), (3.22) and (3.24), respectively. Numerical results for varying levels of slip: 
$\xi =l/(C(1-\eta ))=0.1,\ 0.01,\ 0.001$.
3.4 Spherical squirmer
The squirmer, first proposed by Reference LighthillLighthill (1952), is a standard model for a self-propelled particle in Stokes flow. The spherical squirmer (of radius 
$R$) creates axisymmetric surface motions, modelled by tangential and radial surface velocities (
$u_\theta$, 
$u_r$), that in turn generate a translational axial velocity (
$W$); see figure 1 for the spherical coordinate system.
3.4.1 The no-slip solution
Lighthill first developed the no-slip analytical solution for the squirmer, but this was later corrected by Reference BlakeBlake (1971), which we reproduce here in a slightly different form. As an interesting aside, the reciprocal theorem can be used to determine translational and rotational swimming speeds without solving the full boundary-value problem that follows; the interested reader is referred to Reference Stone and SamuelStone & Samuel (1996).
In the laboratory frame, the no-slip Stokes solution around the particle can be decomposed into a part due to translational wall motion (
$\boldsymbol {u}_0^{t}$, 
$\boldsymbol {\sigma }_0^{t}$) and a part due to wall motion relative to that translation, i.e. the squirming motion (
$\boldsymbol {u}_0^{s}$, 
$\boldsymbol {\sigma }_0^{s}$),
\begin{equation} \boldsymbol{u}_0 = \boldsymbol{u}_0^{t} + \boldsymbol{u}_0^{s}, \quad \boldsymbol{\sigma}_0 = \boldsymbol{\sigma}_0^{t} + \boldsymbol{\sigma}_0^{s}, \end{equation}
where 
$\boldsymbol {u}_0$ and 
$\boldsymbol {\sigma }_0$ are the total velocity and stress fields of the no-slip solution, respectively. Both translational (superscript 
$t$) and squirming (superscript 
$s$) components of the solution decay to zero in the far field. The respective no-slip boundary conditions at the particle surface (
$S_1$), in spherical polar coordinates, are
\begin{equation} \left. \begin{array}{ll} {u}^{t}_r(r,\theta) = W_0 \cos \theta, & \displaystyle {u}^{s}_r(r,\theta) = \sum_{n = 1}^{\infty} A_n P_n(\cos \theta) \\ {u}^{t}_\theta(r,\theta) = {-}W_0 \sin \theta, & \displaystyle {u}^{s}_\theta(r,\theta) = \sum_{n = 1}^{\infty} B_n V_n(\cos \theta) \end{array} \right\}\quad \mathrm{for} \ r = R, \end{equation}
where 
$r$ is the radial coordinate, 
$\theta$ is the polar angle, 
$W_0$ is the speed of particle translation along the polar axis (
$\theta =0$), 
$A_n$ and 
$B_n$ are coefficients describing the form and strength of the squirming motion, 
$P_n$ are Legendre polynomials and 
$V_n(\cos \theta )=-2/(n(n+1))\, {\rm d}P_n(\cos \theta )/{\rm d}\theta$. Note, here, we have restricted attention to volume-preserving wall motions (i.e. the particle cannot lose or gain mass).
Blake's no-slip solution for the velocity field, decomposed into the two parts, is
\begin{gather} {u^{t}_r} = W_0 \cos (\theta ) \left(\frac{3 R}{2 r}-\frac{R^3}{2 r^3}\right), \end{gather}
\begin{gather}{u^{t}_\theta} = {-}W_0 \sin (\theta ) \left(\frac{R^3}{4 r^3} + \frac{3 R}{4 r}\right), \end{gather}
\begin{gather}{u^{s}_r} = \frac{1}{2}\sum_{n = 1}^{\infty} \left[ (n A_n -2 B_n)\frac{R^{n}}{r^{n}} + (2B_n-A_n(n-2))\frac{R^{n + 2}}{r^{n + 2}}\right] P_n(\cos \theta), \end{gather}
\begin{gather}{u^{s}_\theta} = \frac{1}{4}\sum_{n = 1}^{\infty} \left[\left((4 -2 n) B_n + A_n n(n-2)\right)\frac{R^{n}}{r^{n}} + \left(2 n B_n-A_n n (n-2)\right)\frac{R^{n + 2}}{r^{n + 2}}\right] V_n(\cos \theta), \end{gather}
which generates the following shear-stress (
$\tau _{r\theta }$) components at the boundary:
\begin{gather} \tau^t_0 = \frac{3 \mu W_0 }{2 R} \sin (\theta ), \end{gather}
\begin{gather}\tau^s_0 = {-}\frac{\mu}{R} \sum _{n = 1}^{\infty}\left(\frac{3}{2} n A_n + (2n + 1)B_n\right) V_n(\cos \theta). \end{gather}
The motile force (of the fluid on the particle) generated by the squirming motion is obtained by integrating the induced traction force over the squirmer surface (
$S_1$)
\begin{equation} F_0 = \boldsymbol{i}_z \boldsymbol{\cdot} \int_{S_1} \boldsymbol{\sigma}^s_0 \boldsymbol{\cdot} \boldsymbol{n}\,{\rm d}S = 2\pi \mu R ( 2B_1-A_1), \end{equation}which must be balanced, if the particle is self-propelled and there is no external force, by the drag generated in translation
\begin{equation} D_0 = {-} \boldsymbol{i}_z \boldsymbol{\cdot} \int_{S_1} \boldsymbol{\sigma}^t_0 \boldsymbol{\cdot} \boldsymbol{n}\,{\rm d}S = 6\pi \mu R W_0, \end{equation}
where 
$\boldsymbol {i}_z$ is a unit vector along the polar axis. Equating (3.36) and (3.37) leads to an expression for the translational speed of the particle
\begin{equation} W_0 = \tfrac{1}{3}(2B_1-A_1), \end{equation}as obtained by Lighthill and Blake. The most significant implication of this result is that it is only the first modes of the radial and tangential surface motions that contribute to particle translation.
3.4.2 First-order slip corrections
Here, we consider the impact of a uniform slip length at the interface between the surface of the squirmer and the suspending fluid. (In some articles, the word ‘slip’ is used to refer to the tangential motion of the squirmer's surface itself – this is not what is meant here.)
Similarly to previous sections, the motile force generated by the squirming motion is expanded to first order in slip length
\begin{equation} \frac{F}{F_0} = 1 + \hat{F}_1 \xi + {\mathscr{O}}(\xi^2), \end{equation}
where 
$\xi =l/R$, and 
$\hat {F}_1=F_1/F_0$ is the first-order slip-correction coefficient that we wish to find. To obtain it, we need the general expression for the first-order slip-correction coefficient, (2.22), and the no-slip shear-stress distribution associated with the squirming motion, 
$\boldsymbol {\tau }_0^s$, from (3.35). Additionally, for the conjugate Stokes flow, we need the no-slip solution for unit translation (
$\boldsymbol {\tau }_0'=\boldsymbol {\tau }_0^t/W_0$), so that the moment of the traction force obtained is in the direction of translation. From (2.22), with 
$M_i=F_i$, 
$L=R$ and 
$\psi =1$, we obtain
\begin{equation} \hat{F}_1 = \frac{R}{F_0 \mu W_0} \int_{S_1} \boldsymbol{\tau}^{t}_0\boldsymbol{\cdot}\boldsymbol{\tau}^{s}_0\, {\rm d}S = \frac{2\pi R^3}{ F_0 \mu W_0} \int^{\pi}_0 {\tau}^{t}_0 {\tau}^{s}_0 \sin{\theta} \,{\rm d}\theta = \frac{3(A_1 + 2B_1)}{A_1-2B_1}. \end{equation}From § 3.1, the slip drag on a translating sphere is shown to be
\begin{equation} \frac{D}{D_0} = 1- \xi + {\mathscr{O}}(\xi^2). \end{equation}
Now, combining (3.37), (3.40) and (3.41) with the condition for self-propulsion (
$F=D$), an expression for the translational velocity is found
\begin{equation} \frac{W}{W_0} = 1 + 4\left(\frac{A_1 + B_1}{A_1-2 B_1}\right)\xi + {\mathscr{O}}(\xi)^2. \end{equation}
This tells us some interesting things about the impact of low levels of slip on the squirmer's swimming speed. For purely tangential squirming motion at the surface (
$A_1=0$), slip hinders swimming (
${W/W_0<1}$). This is because the motile force generated by pure tangential motion is reduced by slip at three times the rate of the translational drag. Conversely, for purely radial wall motion (
$B_1=0$), slip promotes swimming speed (
$W/W_0>1$), because, in this case, motile force is actually increased by low levels of slip.
Table 3 provides numerical verification of (3.42) for three different cases; see Appendix B for numerical details.
Comparison of numerical and analytical predictions of the first-order slip-correction coefficient for translational velocity of the squirmer, calculated via 
$(W/W_0-1)/\xi$.
3.5 Poiseuille flow through arbitrary cross-section channels
In this final example, we consider an internal flow containing inflow and outflow boundaries. Figure 5 shows a long straight channel (length 
${\mathscr {L}}$) with an arbitrary, but constant, cross-section of area 
$A$. The pressure gradient is assumed constant throughout the channel, 
$\boldsymbol {\nabla } p=\boldsymbol {i}(p_{out}-p_{in})/{\mathscr {L}}$, where 
$\boldsymbol {i}$ is a unit vector along the channel length, which generates a volumetric flow rate 
$Q$. The boundary of the fluid domain is separated into two parts: the walls of the channel (
$S_L$), at which there is the potential for slip, which we assume to be constant in the streamwise direction, and the inlet and outlet boundaries (
$S_A$) at which 
$\psi =0$.
A channel flow with arbitrary cross-section and an applied pressure drop (
$\Delta p=p_{in}-p_{out}$).
Our aim in this section is to find the first-order impact of the Navier slip boundary condition on the pressure drop (
$\Delta p=p_{in}-p_{out}$) for a given flow rate. This is closely related to the problem solved by Reference Michelin and LaugaMichelin & Lauga (2015) and Reference Masoud and StoneMasoud & Stone (2019) using reciprocal theorem, except that, there, the tangential slip velocity was prescribed at the channel walls (not the Navier slip condition) and the resulting flow rate determined.
A force balance in the direction of the channel gives us the pressure drop in terms of a traction-force moment over the channel walls
\begin{equation} \Delta p = \frac{1}{A}\int_{S_L} \boldsymbol{i}\boldsymbol{\cdot} \boldsymbol{\sigma}\boldsymbol{\cdot} \boldsymbol{n} \, {\rm d}S, \end{equation}which we expand as previously
\begin{equation} {\Delta p} = {\Delta p}_0 + {\Delta p}_1\xi + \cdots, \end{equation}
where 
$\xi =l/L$ and the characteristic length scale 
$L$ is chosen based on the cross-section.
To evaluate the first-order slip-correction coefficient for this moment we need the general expression from § 2.3, repeated here for convenience
\begin{align} M_1& = \int_{S} \boldsymbol{u}'_0\boldsymbol{\cdot} \boldsymbol{\sigma_1}\boldsymbol{\cdot} \boldsymbol{n} \, {\rm d}S, \end{align}
\begin{align} & = \frac{L}{ \mu} \int_{S} \psi \boldsymbol{\tau}_{0}\boldsymbol{\cdot}\boldsymbol{\tau}'_{0}\, {\rm d}S. \end{align}The conjugate solution to obtain the desired moment is a simple transformation of the original no-slip solution
\begin{equation} \boldsymbol{u}_0' = {-}\frac{{\boldsymbol{u}_0}}{Q} + \frac{\boldsymbol{i}}{A}\, \quad \mathrm{and} \quad \boldsymbol{\sigma}'_0 = {-}\frac{\boldsymbol{\sigma}_0}{Q}. \end{equation}
To demonstrate this choice gives the correct moment, we substitute the conjugate velocity field (3.46a) into (3.45a) (noting that 
$\int _{S_A} \boldsymbol {u}_0' \, {\rm d}A =\boldsymbol {0}$ and 
$\boldsymbol {u}_0=\boldsymbol {0}$ at 
$S_L$)
\begin{equation} M_1 = \frac{1}{A}\int_{S_L} \boldsymbol{i}\boldsymbol{\cdot} \boldsymbol{\sigma_1}\boldsymbol{\cdot} \boldsymbol{n} \, {\rm d}S -\Delta p_1 \boldsymbol{i}\boldsymbol{\cdot} \int_{S_A} \boldsymbol{u}'_0 \, {\rm d}A = \Delta p_1. \end{equation} Finally, substituting the conjugate stress field (3.46b) into (3.45b), and recalling 
$\psi =0$ at 
$S_A$, gives
\begin{equation} {\Delta p}_1 = {-} \frac{L {\mathscr{L}}}{Q \mu} \int_{{\mathscr{P}}} \psi {\tau}_{0}^2\, {\rm d}{\mathscr{P}}. \end{equation} As is intuitive, perhaps, the first-order slip correction is negative for any distribution of positive slip length around the perimeter, 
${\mathscr {P}}$, of any cross-sectional shape: low levels of slip will always reduce the pressure loss in a Poiseuille flow for a given flow rate.
4. Summary and discussion
We have adopted, and generalised (to arbitrary traction moments), a convenient method for deriving analytical solutions to Stokes flows with low levels of slip, and applied this to a range of micro and nano-flow applications. In general, these first-order approximations to slip Stokes flows are both much simpler to derive and much simpler to evaluate than the full-slip solution. Of course, in many situations, the slip length will not be small, and the methods presented here can only be considered approximate. For example, for the drag on a sphere (of radius 
$R$) with slip length (
$l$), the percentage error of the first-order approximation is given by
\begin{equation} \mathrm{\%\,error} = 100\times \frac{3 \xi^2}{2\xi + 1}, \end{equation}
where 
$\xi =l/R$. The error increases from 2.5 % at 
$\xi =0.1$ to nearly 40 % at 
$\xi =0.5$, to greater than 100 % for 
$\xi >1$. Note, however, for rarefied-gas ‘slip flows’, the first-order approximation is, in fact, the only valid one.
The first-order slip-correction coefficient characterises the behaviour of a given property in low-slip conditions (e.g. drag in the slip-flow regime in rarefied-gas dynamics) as a function of slip length. It is therefore more informative than a single prediction at a specific slip length. However, numerically calculating first-order slip-correction coefficients from the gradient of the property in question can be extremely computationally demanding. This is because its evaluation requires calculating the difference between a very low-slip solution and the no-slip solution – a difference that is very small, and thus hard to calculate accurately. The numerical techniques used in this paper to verify the derived analytical expressions (the method of fundamental solutions (MFS), see Appendix B) are very high accuracy for certain classes of geometry. However, in general, the gradient approach to calculating the slip-correction coefficients is not straightforward.
The expressions derived in § 2 offer a more convenient numerical method of obtaining the first-order slip-correction coefficients: one involving numerical integration of the square of the shear-stress magnitude from an analytical/numerical no-slip solution to the flow problem(s). What is required is an accurate means of integrating properties over the bounding surfaces – in the case of an MFS framework a convenient solution exists for closed surfaces (Reference LockerbyLockerby 2022), but a number of techniques can be employed.
An alternative is to calculate 
$\boldsymbol {u}_1, \boldsymbol {\sigma }_1$ directly, by solving the Stokes equations with the boundary conditions taken from the no-slip stress, 
$\boldsymbol {\sigma }_0$, as per (2.10). This allows the construction of the whole slip solution, to first order, using the original expansion: (2.7) and (2.8). In a similar way, higher-order solutions, and their associated traction-force moments, could also be obtained, allowing better predictions at higher 
$\xi$. Of course, at some value of 
$\xi$, the series will diverge, and so additional terms in the expansion will yield diminishing returns.
In this article, general results pertaining to Stokes flow in low-slip conditions have been derived. For example, consider the rotation of a particle of arbitrary geometry driven by an external torque. The addition of any (small) slip length, however it is distributed across the particle surface, will reduce the torque required to maintain the particle's rotational speed. This is because, in general, the correction coefficient for the resistive moment, (2.27), is always negative (for positive slip lengths).
Acknowledgments
The author would like to thank J. Sprittles and H. Stone for very helpful comments on a draft of the manuscript. For the purpose of open access, the author has applied a CC BY public copyright licence to any Author Accepted Manuscript version arising from this submission. This work was supported by the Engineering and Physical Sciences Research Council (EPSRC), UK [grant number EP/V01207X/1].
Declaration of interests
The author reports no conflict of interest.
Data availability statement
The data that support the findings of this study are tabulated within the article.
Appendix A. Prolate and oblate spheroids
Stokes flow along the axis of revolution of a no-slip spheroid has been solved, analytically, with a number of approaches (Reference OberbeckOberbeck 1876; Reference Payne and PellPayne & Pell 1960; Reference Happel and BrennerHappel & Brenner 1983). The problem of slip flow is substantially more complex. Reference Keh and ChangKeh & Chang (2008) dedicated a full article to the derivation, involving an infinite-series form of semi-separation of variables; truncating the series after two terms still requires a page of algebra to define the analytical coefficients.
Here, by contrast, we aim for a short and simple closed-form expression for the drag on a spheroid in low-slip conditions. Note, these first-order results were first derived by Reference Masoud and StoneMasoud & Stone (2019), using the reciprocal-theorem method, but they were presented in a different form and not numerically verified in that work. We repeat them here in the appendix for illustration and completeness.
The implicit equation for the surface of a spheroid in cylindrical polar coordinates is
\begin{equation} \frac{r^2}{b^2} + \frac{z^2}{a^2} = 1, \end{equation}
where 
$b$ is the spheroid's equatorial radius and 
$a$ is the distance from centre to either pole; see figure 6.
A spheroid in axial translation: (a) prolate and (b) oblate.
A.1 Prolate spheroids, $a>b$
The no-slip drag on a prolate spheroid in axial translation is given by (Reference Payne and PellPayne & Pell 1960; Reference ShermanSherman 1990)
\begin{equation} D_0 = \frac{16 \pi \mu {\mathscr{L}} W}{(s^2 + 1) \log \left(\dfrac{s + 1}{s-1}\right)-2 s}, \end{equation}
where 
$W$ is the translational velocity of the spheroid in the direction of 
$z$, 
$s={E}/{\sqrt {E^2-1}}$, 
$E=a/b$ and 
${\mathscr {L}}=\sqrt {a^2-b^2}$ is the focal length.
The shear-stress magnitude from the no-slip solution is given by
\begin{equation} \tau_0 = {-}\frac{4 E \mu W \sin (\eta )}{{\mathscr{L}} (\cos (2 \eta )-2 s^2 + 1) ((s^2 + 1) \log (\coth (\tfrac{1}{2} \cosh^{{-}1}(s)))-s)}, \end{equation}
where 
$\eta$ is a coordinate on the spheroid surface related to cylindrical polar coordinates through 
$z={\mathscr {L}}s \cos (\eta )$ and 
$r={\mathscr {L}} s \sin (\eta )/E$.
The slip-flow drag on the prolate spheroid is approximated to first order by
\begin{equation} D/D_0 = 1 + \hat{D_1} \xi + {\mathscr{O}}(\xi^2), \end{equation}
where 
$\xi =l/b$. For translation, drag is the resistive moment, and so the first-order slip-correction coefficient is obtained by evaluating equation (2.27) (with 
${\mathscr {R}}= D$, 
$L=b$, 
$\psi =1$ and 
$k = W$)
\begin{equation} \hat{D}_1 = {-}\frac{b}{D_0 U \mu} \int_{S} {\tau}^2_{0} \, {\rm d}S = {-}\frac{2(1- s E \tan ^{{-}1}(E/s))}{ E (1-(s + 1/s) \coth ^{{-}1}(s))}. \end{equation}Given the complexity of the full-slip derivation and solution due to Reference Keh and ChangKeh & Chang (2008), (A5) is remarkably simple.
In table 4, numerical calculations using the MFS (see Appendix B) are compared with (A5) for a range of spheroid aspect ratios (
$E$); as expected, the full-slip numerical results converge to the analytical solution as the slip length is reduced. For 
$\xi =10^{-2}$ the analytical result is within 
$3\,\%$ of the numerical solutions; for 
$\xi =10^{-5}$ the analytical result is within 0.02 % of the numerical solutions.
Results for the prolate spheroid. Comparison of the analytical prediction of the first-order slip-correction coefficient for drag (
$\hat {D}_1$) and numerical calculation of 
$(D/D_0-1)/\xi$.
A.2 Oblate spheroids, $a< b$
The no-slip solution for drag on the oblate spheroid is (Reference Payne and PellPayne & Pell 1960; Reference ShermanSherman 1990)
\begin{equation} D_0 = \frac{8 \pi \mu {\mathscr{L}} W}{t-(t^2-1)\cot ^{{-}1}(t)}, \end{equation}
where 
$t=\sinh (\cosh ^{-1}(s))$ and the corresponding shear-stress magnitude distribution is
\begin{equation} \tau_0 = \frac{4 \mu W \sin (\eta )}{{\mathscr{L}} E (t-(t^2-1) \cot ^{{-}1}(t)) (\cos (2 \eta ) + \cosh (2 \sinh ^{{-}1}(t)))}, \end{equation}
where, now, 
$E=b/a$, 
${\mathscr {L}}=\sqrt {b^2-a^2}$ and 
$\eta$ is a coordinate on the surface of the spheroid related to cylindrical polar coordinates through 
$z={\mathscr {L}} t \cos (\eta )$ and 
$r={\mathscr {L}} s \sin (\eta )$.
The drag on the oblate spheroid in slip flow can be expanded, as in (A4), but with 
$\xi =l/a$, and the first-order slip-correction coefficient obtained from (2.27) (with 
${\mathscr {R}}_i= D_i$, 
$L=a$, 
$\psi =1$ and 
$k = W$)
\begin{equation} \hat{D}_1 = {-}\frac{a}{D_0 U \mu} \int_{S} {\tau}^2_{0} \, {\rm d}S = {-}\frac{2( E - t \coth ^{{-}1}(s))}{ E (1-(t-1/t)\cot ^{{-}1}(t))}. \end{equation}As for the prolate case, the simplicity of the derivation and the final result is noteworthy. In table 5, numerical calculations (see Appendix B) verify (A8) and also provide an indication of the loss in accuracy of the low-slip assumption as the slip length increases.
Results for the oblate spheroid. Comparison of the analytical prediction of the first-order slip-correction coefficient for drag (
$\hat {D}_1$) and numerical calculation of 
$(D/D_0-1)/\xi$.
Appendix B. The method of fundamental solutions
For the numerical calculation of the external Stokes flows considered in §§ 3.1, 3.2, A and 3.4, we employ the MFS (Reference Lockerby and CollyerLockerby & Collyer 2016; Reference Cheng and HongCheng & Hong 2020). The MFS uses a superposition of fundamental solutions to the Stokes equations (popularly known as Stokeslets) to construct an analytical solution that approximately satisfies the boundary conditions at the particle surface (a zero-disturbance far-field condition is automatically satisfied by the Stokeslets). The numerical procedure has much in common with the boundary element method, requiring the surface of the particle to be discretised as opposed to the volume of fluid it occupies; this reduces the dimensionality of the numerical calculation and removes the requirement for a finite fluid domain. The primary numerical parameter is the number of ‘boundary nodes’ (sometimes referred to as collocation points) that discretise the particle boundary, 
$N$; as 
$N$ increases the superposed solution becomes a more accurate representation of the true one.
The Stokeslet can be viewed as the Stokes-flow response to a steady-state point forcing in three-dimensional space. The MFS approximates a given flow using a superposition of these Stokeslets, with the location of the point forces (referred to as singularity sites) distributed outside of the fluid domain (e.g. set within the volume of a solid particle, 
$V_p$, see figure 7).
Illustration of site and node arrangement in the MFS applied to external flows around particles.
The governing equations that are solved (exactly) by the MFS are
\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = \boldsymbol{0}, \quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\sigma} = {-}\sum_{s = 1}^{M} \boldsymbol{f}_s \delta(\boldsymbol{r}-\boldsymbol{r}_s), \end{equation}
where 
$\boldsymbol {f}_s$ is the force (located at 
$\boldsymbol {r}_s$) associated with the 
$s$th Stokeslet, and 
$\delta$ is the Dirac delta function. Note, these revert to the Stokes equations, (2.1), within the fluid volume (assuming no external forces). The corresponding analytical solution is
\begin{gather} \boldsymbol{u} = \frac{1}{8 \pi \mu} \sum_{s = 1}^{M}\boldsymbol{f}_s\boldsymbol{\cdot} \left( \frac{\boldsymbol{\mathsf{I}}}{\|\boldsymbol{r}-\boldsymbol{r}_s\|} + \frac{(\boldsymbol{r}-\boldsymbol{r}_s)(\boldsymbol{r}-\boldsymbol{r}_s)}{\|\boldsymbol{r}-\boldsymbol{r}_s\|^3} \right), \end{gather}
\begin{gather}\boldsymbol{\sigma} = {-}\frac{3 }{4 \pi } \sum_{s = 1}^{M}\boldsymbol{f}_s \boldsymbol{\cdot}\left(\frac{(\boldsymbol{r}-\boldsymbol{r}_s)(\boldsymbol{r}-\boldsymbol{r}_s)(\boldsymbol{r}-\boldsymbol{r}_s)}{\|\boldsymbol{r}-\boldsymbol{r}_s\|^5} \right) , \end{gather}
where 
$\boldsymbol{\mathsf{I}}$ is the identity tensor, and 
$\|\cdot \|$ denotes the Euclidean (L2) norm. The primary aim of the MFS is to determine the Stokeslet forces, 
$\boldsymbol {f}_s$, from the boundary conditions of the problem.
Evaluating the velocity and stress field at the boundary nodes gives
\begin{equation} \boldsymbol{u}_{b} = \frac{1}{\mu}\sum_{s = 1}^{M} \boldsymbol{f}_s \boldsymbol{\cdot} \boldsymbol{\mathsf{J}}_{\,bs} \quad \mathrm{and} \quad \boldsymbol{\sigma}_{b} = {-}\sum_{s = 1}^{M} \boldsymbol{f}_s \boldsymbol{\cdot} \boldsymbol{\mathsf{K}}_{bs}, \end{equation}where
\begin{equation} \boldsymbol{\mathsf{J}}_{bs} = \frac{1}{8 \pi}\left( \frac{\boldsymbol{\mathsf{I}}}{\|\boldsymbol{r}_{bs}\|} + \frac{\boldsymbol{r}_{bs} \boldsymbol{r}_{bs}}{\|{\boldsymbol{r}}_{bs}\|^3} \right) \quad \mathrm{and} \quad \boldsymbol{\mathsf{K}}_{bs} = \frac{3}{4 \pi}\left( \frac{ {\boldsymbol{r}}_{bs} {\boldsymbol{r}}_{bs}{\boldsymbol{r}}_{bs} }{\|{\boldsymbol{r}}_{bs}\|^5} \right), \end{equation}
and where 
$\boldsymbol {r}_{bs}=\boldsymbol {r}_b-\boldsymbol {r}_s$, and 
$\boldsymbol {r}_b$ is the position of the 
$b$th boundary node; see figure 7. Substituting (B4a,b) into (2.3) allows the Navier slip boundary condition at node 
$b$ to be written
\begin{equation} \boldsymbol{U}_b = \frac{1}{\mu}\sum_{s = 1}^{M}\boldsymbol{f}_s\boldsymbol{\cdot} \boldsymbol{\mathsf{A}}_{bs}, \end{equation}where
\begin{equation} \boldsymbol{\mathsf{A}}_{bs} = \boldsymbol{\mathsf{J}}_{bs} + \ell_b (\boldsymbol{n}_b \boldsymbol{\cdot} \boldsymbol{\mathsf{K}}_{bs})\boldsymbol{\cdot} (\boldsymbol{\mathsf{I}} -\boldsymbol{n}_b \boldsymbol{n}_b ). \end{equation}
The rank-4 tensor, 
$\boldsymbol{\mathsf{A}}_{bs}=A_{bsij}$, and rank-2 tensors, 
$\boldsymbol {f}_s=f_{sj}$ and 
$\boldsymbol {U}_b=U_{bi}$, can be reshaped with the bijection
\begin{equation} p = 3(b-1) + i \quad \mathrm{and} \quad q = 3(s-1) + j, \end{equation}to obtain
\begin{equation} \tilde{A}_{pq} = A_{bsij}, \quad \tilde{f}_{q} = f_{sj}, \quad \mathrm{and} \quad \tilde{U}_{p} = U_{bi}, \end{equation}
where 
$\tilde {\boldsymbol{\mathsf{A}}}$ has size (
$3N \times 3M$), 
$\tilde {\boldsymbol {f}}$ has size (
$3M\times 1$), 
$\tilde {\boldsymbol {U}}$ has size (
$3N \times 1$), such that the evaluation of the Navier slip condition at all nodes is represented by the matrix equation
\begin{equation} \tilde{\boldsymbol{\mathsf{A}}} \boldsymbol{\cdot} \tilde{\boldsymbol{f}} = \tilde{\boldsymbol{U}}, \end{equation}
which can be solved for the vector of force components (
$\tilde {\boldsymbol {f}}$) by any standard linear-equation solver. In practice, it benefits the numerics to have fewer Stokeslets than boundary nodes, which creates an overdetermined system that can be solved using a linear least-squares method. In this work, 
$M\approx 0.9 N$.
For the simulations presented in this article, the boundary nodes on the surface of each object (spheres or spheroids) are distributed evenly, and found using the Matlab code (DistMesh) written by Reference Persson and StrangPersson & Strang (2004). In the MFS literature the singularity sites are often referred to as ‘source nodes’ and considerable work has been done on deciding how they should be optimally located (Reference KarageorghisKarageorghis 2009; Reference Chen, Karageorghis and LiChen, Karageorghis & Li 2016). However, for the problems considered here, a simple surface-normal projection into the particle works well, that is: 
$\boldsymbol {r}_s= \boldsymbol {r}_b-\alpha _b \boldsymbol {n}_b$, where 
$\alpha _b$ is chosen such that each site is 5 % closer to its respective node than any other. After locating the sites in this way, approximately 10 % of them are deleted, so that 
$M\approx 0.9N$.
Once the forces are known, it is simple to calculate moments of the traction force on a particle via (B.1) and the divergence theorem. For example, the net traction force 
$\boldsymbol {F}$ on the particle surface is simply a summation of the Stokeslet forces
\begin{equation} \boldsymbol{F} = \int_{S} \boldsymbol{\sigma}\boldsymbol{\cdot} \boldsymbol{n} \, {\rm d}S = \int_{V_p} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\sigma} \, {\rm d}V_p = {-} \int_{V_p} \sum_{s = 1}^{M} \boldsymbol{f}_s \delta(\boldsymbol{r}-\boldsymbol{r}_s)\, {\rm d}V_p = {-} \sum_{s = 1}^{M} \boldsymbol{f}_s. \end{equation}B.1 Numerical convergence
To verify the implementation of the MFS for Stokes slip flows, we compare numerical predictions of drag around a sphere with the classic analytic result due to Reference BassetBasset (1888), (3.1). Table 6 shows convergence of the MFS result to the exact drag solution with increasing number of nodes (
$N$) for a range of non-dimensional slip lengths.
Results for drag on a translating sphere in Stokes flow with slip (normalised with no-slip drag). Comparison of an analytical solution (Reference BassetBasset 1888) with the MFS.
Obtaining slip-correction coefficients from numerical simulations can be more demanding. For example, in table 4, the quantity of interest is 
$(\hat {D}-1)/\xi$. To obtain this quantity accurate to 4 significant figures, in the case when 
$\xi =10^{-5}$, requires the calculation of 
$\hat {D}$ accurate to 9 significant figures. The accuracy of the MFS is therefore essential for the purposes of verification.
Table 7 summarises the number of boundary nodes used in the various verification cases of the main article. Halving the number of nodes used in each case (or quartering, in the case of the oblate spheroid) results in a small change in the presented results; see the penultimate column of table 7.
Number of boundary nodes used in MFS verification simulations.
The MFS performs far worse when particles have sharp edges or high aspect ratio. In the example of the variable-slip-length sphere, § 3.2, the MFS has to resolve discontinuities in slip length across the sphere's surface (see figure 2). As such, this represents the most challenging case of the article, and requires a large number of boundary nodes to get accurate results.
Open access
