Optimal slip velocities of micro-swimmers with arbitrary axisymmetric shapes

Abstract

This article presents a computational approach for determining the optimal slip velocities on any given shape of an axisymmetric micro-swimmer suspended in a viscous fluid. The objective is to minimize the power loss to maintain a target swimming speed, or equivalently to maximize the efficiency of the micro-swimmer. Owing to the linearity of the Stokes equations governing the fluid motion, we show that this PDE-constrained optimization problem reduces to a simpler quadratic optimization problem, whose solution is found using a high-order accurate boundary integral method. We consider various families of shapes parameterized by the reduced volume and compute their swimming efficiency. Among those, prolate spheroids were found to be the most efficient micro-swimmer shapes for a given reduced volume. We propose a simple shape-based scalar metric that can determine whether the optimal slip on a given shape makes it a pusher, a puller or a neutral swimmer.

1. Introduction

The squirmer model (Lighthill 1952; Blake 1971) has been widely adopted by mathematicians and physicists over the past decades to model ciliated micro-swimmers such as Opalina, Volvox and Paramecium (Lauga & Powers 2009). On a high level, this continuum model, sometimes referred to as the envelope model, effectively tracks the motion of the envelope formed by the tips of the densely packed cilia, located on the swimmer's body, while neglecting the motion below the tips. Individual and collective ciliary motions could be mapped to travelling waves of the envelope on the surface. Assuming no radial displacements of the surface and a time-independent tangential velocity led to the simpler steady squirmer model (see Pedley 2016), wherein, a prescribed slip velocity on the boundary propels the squirmer. While the model was originally designed for spherical shapes, it has since been adapted to more general shapes and has recently been shown to capture realistic collective behaviour of suspensions (Kyoya et al. 2015).

Shape is also a key parameter in the design of artificial micro-swimmers for promising applications such as targeted drug delivery. In particular, the squirmer model is often employed to study the propulsion of phoretic particles, which are micro- to nano-metre sized particles that propel themselves by exploiting the asymmetry of chemical reactions on their surfaces (Anderson 1989; Golestanian, Liverpool & Ajdari 2007). A classical example is the Janus sphere (Howse et al. 2007), which consists of inert and catalytic hemispheres. When submerged in a suitable chemical solution, the asymmetry between the chemical reactions on the two hemispheres creates a concentration gradient. The gradient creates an effective steady slip velocity on the surface via osmosis that naturally suits the squirmer model. Besides the classical Janus spheres and bi-metallic nanorods (Paxton et al. 2004), more sophisticated shapes have also been proposed recently, such as two spheres (Valadares et al. 2010; Palacci et al. 2015), spherocylinder (Uspal et al. 2018), matchsticks (Morgan et al. 2014) and microstars (Simmchen et al. 2017). Interestingly, Uspal et al. (2018) showed that special shapes of phoretic particles exhibit novel properties such as ‘edge following’ when placed close to chemically patterned surfaces.

Studying the efficiency of biological micro-swimmers is pivotal to understanding natural systems and designing artificial ones for accomplishing various physical tasks. The mechanical efficiency (Lighthill 1952) of the spherical squirmer can be directly computed as its rate of viscous energy dissipation, or power loss, can be written in terms of the modes of the squirming motion. Michelin & Lauga (2010) found the optimal swimming strokes of unsteady spherical squirmers by employing a pseudo-spectral method for solving the Stokes equations that govern the ambient fluid and a numerical optimization procedure. Their approach, however, does not readily generalize to arbitrary shapes. On the other hand, Leshansky et al. (2007) analytically investigated the efficiency of micro-swimmers with prolate spheroid shapes with a time-independent ‘treadmilling’ slip velocity and found that the optimal efficiency increases unboundedly with the aspect ratio. Vilfan (2012) optimized the steady slip velocity and the shape at the same time, with constraints on its volume and maximum curvature. The work considered power loss inside the squirmer surface, which could be an order of magnitude higher than the outside power loss (Keller & Wu 1977; Ito, Omori & Ishikawa 2019). However, it assumed that the tangential force on the squirmer surface is linear to its local slip velocity, which is not always the case for micro-swimmers.

In this paper, we address the following broader question: Given an axisymmetric shape of a steady squirmer, what is the slip velocity that maximizes its swimming efficiency? The optimization problem, being quadratic, is reduced to a linear system of equations solved by a direct method, while forward exterior flow problems are solved using a boundary integral method. Those combined features produce a simple and efficient solution procedure. We introduce the optimization problem and our numerical solver in § 2, present the optimal solution for various shape families, summarize the correlations between the shapes and the optimal slip velocities and propose a shape-based scalar metric to predict whether the optimized swimmer would be a pusher or a puller in § 3, followed by conclusions and a discussion on future research directions in § 4.

2. Problem formulation and numerical solution

2.1. Model

Consider an axisymmetric micro-swimmer whose boundary $\varGamma$ can be obtained by rotating a curve $\gamma$ about the $\boldsymbol {e}_3$, axis as shown in figure 1(a). Using the arc length $s\in [0,\ell ]$ to parameterize the generating curve, its coordinate functions can be written as $\gamma (s) = (x_1(s), 0, x_3(s))$. Here, we restrict our attention to shapes of spherical topology, therefore, all shapes considered satisfy the conditions $x_1(0) = x_1(\ell ) = 0$ and $x_1(s)>0, \ \forall \, s\in (0,\ell )$. We assume that the micro-swimmer is suspended in an unbounded viscous fluid domain. The governing equations for the ambient fluid in the vanishing Reynolds number limit are given by the Stokes equations

(2.1a,b)\begin{equation} -\mu\nabla^{2}\boldsymbol{u} + \boldsymbol{\nabla} p = \boldsymbol{0},\quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \end{equation}

where $\mu$ is the fluid viscosity, $p$ and $\boldsymbol {u}$ are the pressure and flow field respectively. In the absence of external forces and imposed flow fields, the far-field boundary condition simply is

(2.2)\begin{equation} \lim_{\boldsymbol{x}\rightarrow\infty}\boldsymbol{u}(\boldsymbol{x}) = \boldsymbol{0}. \end{equation}

A tangential slip ${u}^{{S}}$ defined on $\gamma$ propels the micro-swimmer forward with a translational velocity $U$ in the $\boldsymbol {e}_3$ direction. Its angular velocity as well as the translational velocities in the $\boldsymbol {e}_1$ and $\boldsymbol {e}_2$ directions are zero by symmetry. Consequently, the boundary condition on $\gamma$ is given by

(2.3)\begin{equation} \boldsymbol{u} = {u}^{{S}}\boldsymbol\tau+{U}\boldsymbol{e}_3, \end{equation}

where $\boldsymbol {\tau }$ is the unit tangent vector on $\gamma$. Note that, in order to avoid singularities, the slip must vanish at the endpoints

(2.4)\begin{equation} u^{{S}}(0) = u^{{S}}(\ell) = 0. \end{equation}

Due to the axisymmetry of $\varGamma$, the required no-net-torque condition on the freely suspended micro-swimmer is automatically satisfied while the no-net-force condition reduces to one scalar equation

(2.5)\begin{equation} \int_\varGamma \boldsymbol{f}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{e}_3 \, \mathrm{d}S=2{\rm \pi}\int_\gamma {f}_3(\boldsymbol{x})\, x_1\mathrm{d}s= 0, \end{equation}

where $\boldsymbol {f}$ is the active force density on the micro-swimmer surface (negative to fluid traction) and $f_3$ is its $\boldsymbol {e}_3$ component.

Figure 1. $(a)$ Schematic of the micro-swimmer geometry. The shape is assumed to be axisymmetric, obtained by rotating the generating curve $\gamma$ about the $\boldsymbol {e}_3$ axis. $(b)$ Biological swimmers (Lynn 2008, chapter 4, figure 4.6). $(c)$ Scanning electron microscope (SEM) image of a single half-coated Janus particle; inset: the dark blue shows the location of the platinum (Pt) cap (Choudhury et al. 2017). $(d)$ The SEM image of a phototactic swimmer, which consists of a haematite particle extruded from a colloidal bead (Aubret & Palacci 2018).

We quantify the performance of the micro-swimmer with slip velocity $u^{{S}}$ by its power loss while maintaining a target swimming speed $U$. The power loss is defined by

(2.6)\begin{equation} P = \int_\varGamma\boldsymbol{f}\boldsymbol{\cdot}\boldsymbol{u} \,\mathrm{d}S = 2{\rm \pi}\int_\gamma \boldsymbol{f}\boldsymbol{\cdot}({u}^{{S}}\boldsymbol\tau+U\boldsymbol{e}_3) x_1\,\mathrm{d}s. \end{equation}

Note that $P$ can be made arbitrarily small by lowering the swimming speed $U$. It is therefore necessary to compare the power loss of different swimmers that have the same swimming speed $U$. We note that a lower $P$ with a fixed shape and swimming speed $U$ corresponds to a higher efficiency, $\eta = C_DU^{2}/P$, as defined by Lighthill (1952), where $C_D$ is the drag coefficient of the given swimmer.

2.2. Boundary integral method for the forward problem

Before stating the optimization problem, we summarize our numerical solution procedure for (2.1a,b)–(2.3). Again, we fix the swimming speed $U$, referred to from here onwards as the ‘target swimming speed’, and assume that the tangential slip $u^{S}$ is given. In general, an arbitrary pair of ${u}^{{S}}$ and $U$ does not satisfy the no-net-force condition (2.5). This condition will be treated as a constraint in our optimization problem. Therefore, the goal is to find the active force density $\boldsymbol {f}$ given the velocity on the boundary $\gamma$ as in (2.3). We use the single-layer potential ansatz (Pozrikidis 1992), which expresses the velocity as a convolution of an unknown density function $\boldsymbol {\mu }$ with the Green's function for the Stokes equations $G$, from which the force density can be determined by convolution with the traction kernel $T$

(2.7a,b)\begin{align} \boldsymbol{u}(\boldsymbol{x}) = \int_\varGamma G(\boldsymbol{x}\,{-}\,\boldsymbol{y}) \, \boldsymbol{\mu}(\boldsymbol{y}) \, \textrm{d}\varGamma(\boldsymbol{y}), \enspace \boldsymbol{f}(\boldsymbol{x}) ={-}\tfrac{1}{2}\boldsymbol{\mu}\left(\boldsymbol{x}\right)+\boldsymbol{n}\left(\boldsymbol{x}\right)\int_{\varGamma}T\left(\boldsymbol{x}\,{-}\,\boldsymbol{y}\right)\boldsymbol{\mu}\left(\boldsymbol{y}\right)\,\textrm{d}\varGamma\left(\boldsymbol{y}\right), \end{align}

where $\boldsymbol {n}$ is the unit normal vector pointing into the fluid. We can solve for $\boldsymbol {\mu }$ by taking the limit of $\boldsymbol {x} \rightarrow \varGamma$ in the above ansatz and substituting in (2.3). The boundary integrals in (2.7a,b) become weakly singular on $\varGamma$, requiring specialized quadrature rules. Here, we use the approach of Veerapaneni et al. (2009), which performs an analytic integration in the $\theta -$direction reducing the integrals to convolutions on the generating curve and applies a high-order quadrature rule designed to handle the log-singularity of the resulting kernels. More details on the numerical scheme are provided in appendix B.

2.3. Optimization problem and its reformulation

The goal is to find a slip profile $u^{{S}*}(s)$ that minimizes the power loss $P$ while maintaining the target swimming speed $U$ of a given axisymmetric micro-swimmer. Let ${J}$ be the objective function, here equated to $P$ defined in (2.6), and ${F}$ be the net force functional

(2.8a,b)\begin{equation} {J}(u^{{S}}) := 2{\rm \pi}\int_\gamma \boldsymbol{f}(u^{{S}})\boldsymbol{\cdot}(u^{{S}}\boldsymbol\tau + U\boldsymbol{e}_3)\, x_1 \mathrm{d}s, \quad {F}(u^{{S}}) :=2{\rm \pi} \int_\gamma \boldsymbol{f}(u^{{S}}) \boldsymbol{\cdot} \boldsymbol{e}_3\, x_1 \mathrm{d}s. \end{equation}

They are slip velocity functionals as their values are completely determined by $u^{{S}}$. The optimization problem can now be stated as follows:

(2.9)\begin{equation} u^{{S*}}=\mathop{{\textrm{arg\,min}}}\limits_{u^{{S}}\in \mathcal{U}} {J}(u^{{S}})\quad \text{subject to } {F}(u^{{S}}) = 0, \end{equation}

with $\mathcal {U}$ being the space of the all possible slip velocities satisfying (2.4). Notice that the no-net-force condition (2.5) is added as a constraint here.

By (2.3) and linearity of the Stokes equation (2.1a,b), the forward solution $\boldsymbol {u}$ and the net force ${F}$ are affine in $u^{{S}}$ ($\boldsymbol {u}$ is linear in $u^{{S}}$ if $F=0$). Consequently, ${J}(u^{{S}})$ is a quadratic functional and (2.9) is inherently a quadratic optimization problem. To make it more explicit, consider a discretized version of the slip optimization problem where $u^{{S}}$ is sought in the form

(2.10)\begin{equation} u^{{S}}(\boldsymbol{x}) = \sum_{k=1}^{m} U\xi_k\, u^{{S}}_k(s), \end{equation}

for some set of $m$ basis functions $u^{{S}}_k$ satisfying (2.4). We adopt a B-spline formulation for these basis functions (see appendix A for more details). Let $(\boldsymbol {u}_0,p_0,\boldsymbol {f}_0)$ and $(\boldsymbol {u}_k,p_k,\boldsymbol {f}_k)$ (with $1\le k\le m$) denote the solutions of the forward problem (2.1a,b) with $\boldsymbol {u}=\boldsymbol {e}_3$ and $\boldsymbol {u}=u_k^{{S}}\boldsymbol {\tau }$ being their boundary conditions on $\gamma$, respectively.

The net force ${F}(u^{{S}})$ is then given by ${F}(u^{{S}})= 2{\rm \pi} U\mathcal {F}(\boldsymbol \xi )$, where

(2.11)\begin{equation} \mathcal{F}(\boldsymbol\xi):= \int_\gamma \left(\boldsymbol{f}_0 + \sum_{k = 1}^{m} \xi_k\,\boldsymbol{f}_k \right)\boldsymbol{\cdot} \boldsymbol{e}_3\, x_1 \,\mathrm{d}s= F_0 + \boldsymbol{F}^{\mathrm{T}}\boldsymbol{\xi}. \end{equation}

Here, $\boldsymbol {\xi }=(\xi _1,\ldots ,\xi _m)^{\mathrm {T}}$, $\boldsymbol {F} = (F_1,\ldots ,F_m)^{\mathrm {T}}$ and $F_k = \int _\gamma \boldsymbol {f}_k \boldsymbol {\cdot } \boldsymbol {e}_3\, x_1\, \mathrm {d}s$ for $k=0,1,\ldots ,m$.

Similarly, we have ${J}(u^{{S}}) = 2{\rm \pi} U^{2}\mathcal {J}(\boldsymbol \xi )$, where

(2.12)\begin{equation} \mathcal{J}(\boldsymbol\xi):= \int_\gamma \left(\boldsymbol{f}_0 + \sum_{k = 1}^{m} \xi_k\,\boldsymbol{f}_k \right)\boldsymbol{\cdot}\left(\boldsymbol{e}_3 + \sum_{j=1}^{m} \xi_j u^{{S}}_j\boldsymbol\tau\right)\, x_1 \mathrm{d}s = \boldsymbol{\xi}^{\mathrm{T}}\boldsymbol{A}\boldsymbol{\xi} + 2\boldsymbol{\xi}^{\mathrm{T}}\boldsymbol{F} + F_0. \end{equation}

The elements of the $m\times m$ matrix $\boldsymbol {A}$ are given by $A_{kj} = \int _\gamma \boldsymbol {f}_k \boldsymbol {\cdot } u^{{S}}_j \boldsymbol \tau \, x_1 \mathrm {d}s$. We have used the fact that $\int _\gamma \boldsymbol {f}_0 \boldsymbol {\cdot } u^{{S}}_k \boldsymbol \tau \, x_1 \mathrm {d}s = \int _\gamma \boldsymbol {f}_k \boldsymbol {\cdot } \boldsymbol {e}_3\, x_1 \mathrm {d}s$ for the linear term by the reciprocal theorem (Happel & Brenner 1973). We note that $\boldsymbol {A}$ is symmetric, also by the reciprocal theorem. Physically speaking, $\boldsymbol \xi ^{\mathrm {T}}\boldsymbol {A}\boldsymbol \xi$ represents the scaled power loss of the swimmer being held still with its slip velocity parametrized by $\boldsymbol \xi$, implying that $\boldsymbol {A}$ is positive–definite; $\boldsymbol {\xi }^{T} \boldsymbol {F}$ is the scaled power loss of the active force along the swimming direction; $F_0$ is the scaled power loss of towing a rigid body with the same shape as the micro-swimmer at unit speed.

Now, the discretized optimization problem becomes

(2.13)\begin{equation} \min_{\boldsymbol{\xi}\in\mathbb{R}^{m}} \mathcal{J}(\boldsymbol{\xi}) \quad \text{subject to } \mathcal{F}(\boldsymbol{\xi})=0. \end{equation}

Introducing the Lagrangian $L(\boldsymbol {\xi },\lambda ):= \mathcal {J}(\boldsymbol {\xi })-2\lambda \mathcal {F}(\boldsymbol {\xi })$, the slip optimization problem is reduced to solving the first-order stationarity equations for $L$ given by

(2.14)\begin{equation} \begin{bmatrix} \boldsymbol{A} & -\boldsymbol{F} \\ -\boldsymbol{F}^{\mathrm{T}} & 0 \end{bmatrix} \, \begin{bmatrix} \boldsymbol{\xi} \\ \lambda \end{bmatrix} = \begin{bmatrix} -\boldsymbol{F} \\ F_0 \end{bmatrix}. \end{equation}

Note that forming the matrix requires $(m+1)$ solves of the forward problem (2.1a,b) with appropriate boundary conditions. Since the micro-swimmer is assumed to be rigid, the single-layer potential operator as well as the traction operator, required for forming $\boldsymbol {A}$ and $\boldsymbol {F}$, are both fixed for a given shape. Therefore, we only need to form them once.

3. Results

We tested the convergence of our numerical solvers rigorously; the boundary discretization for all the numerical examples presented here is chosen so that at least 6-digit solution accuracy is attained (determined via self-convergence tests). The optimal slip velocity for a particular prolate spheroid tested against the (truncated) analytical solution given by Leshansky et al. (2007) is shown in figure 2. Our numerical solution is indistinguishable against the analytical solution at their finer truncation level $L=10$. Additional validation results can be found in appendix B.

Figure 2. Optimal slip velocity compared to Leshansky et al. (2007, figure 4). The aspect ratio of the prolate spheroid is $(1+2.5^{2})^{1/2}$. Our numerical optimization is depicted in black solid curve, while dash curves represent analytical solutions at different truncation levels $L=4$ (red) and $L=10$ (blue).

Here, we focus on analysis of the optimal solutions for various micro-swimmer shape families. Let $V$ be the volume enclosed by the swimmer. We normalize lengths by the radius of a sphere of equivalent volume i.e. by $R=(3V/4{\rm \pi} )^{1/3}$, and velocities by the swimming speed $U$. A simple calculation shows that, for a micro-swimmer submerged in water of size $R=5\ \mathrm {\mu }\text {m}$ and the speed of one body length per second, the Reynolds number $(Re) \approx 5\times 10^{-5}$; thereby, confirming the validity of the Stokes equation (2.1a,b). We will use the dimensionless reduced volume, defined by $\nu = 6\sqrt {{\rm \pi} }V/A^{3/2}$ where $A$ is the surface area of the given shape, to characterize each shape family. The largest possible value of $\nu$, attained by spheres, is $\nu =1$, while for example $\nu$ decreases monotonically for spheroids as the aspect ratio is increased.

We first consider six different micro-swimmer shapes and plot their optimal slip profiles obtained by solving (2.14) in figure 3. In each case, we also show the flow fields in both the body and laboratory frames. The optimal slip velocities plotted against the arc length, measured from north pole to south pole, are shown in the insets. In the case of a sphere (figure 3a), we recover the standard result that the optimal profile is a sine curve (Michelin & Lauga 2010). The optimal slip velocity of the prolate swimmer, shown in figure 3(b), ‘flattens’ the sine curve in the middle while that of the oblate swimmer, shown in figure 3(c), ‘pinches’ the sine curve. Additionally, the peak value of the optimal slip velocity is low for the prolate swimmer, and high for the oblate swimmer, compared to the spherical swimmer.

Figure 3. Flow fields and the optimal slip velocity for a few swimmers with typical shapes: $(a)$ sphere, $(b)$ prolate spheroid, $(c)$ oblate spheroid, $(d)$ wavy, $(e)$ spherocylinder, $(f)$ stomatocyte. Insets show the optimal slip velocities as functions of arc length along the generating curve. The optimization is performed using 21 control points on the generating curve to represent the slip velocity. The colour map holds for both the slip velocity and the flow fields.

Next, we consider three shapes corresponding to different shape families. In figure 3(d), we consider the ‘wavy’ configuration obtained by adding high-order axisymmetric modes to the spherical shape. The optimal slip velocity follows the general trend for that of figure 3(a), while lower slip velocities are observed at the troughs, qualitatively consistent to those obtained in Vilfan (2012). The spherocylinder (figure 3e) resembles closely the prolate spheroid of figure 3(b) with the same aspect ratio, its optimal slip velocity being nearly the same (albeit with a slightly narrower plateau and higher peak slip velocity). Finally, we investigate the optimal slip velocity of the stomatocyte shape (figure 3f), which is the only non-convex shape among those considered here. Similar to that of the oblate swimmer, the general slip velocity is like a pinched sine wave. However, one distinguishing feature is that slip velocity is nearly zero over part of its surface, namely the cup-like region in its posterior.

The optimal slip velocity strongly depends on the local geometry of the micro-swimmer. Generally speaking, the optimal slip velocity is high if the material point is far away from the axis of symmetry. This could be seen most clearly in the cases of spheroids, see figure 3(a–c). Specifically, the peak value of the optimal slip velocity is the highest for the oblate spheroid and lowest for the prolate spheroid among the three. Intuitively, an object that has a larger radius would endure a higher fluid drag compare to one with a smaller radius when moving in the same speed. Thus extra effort, in the form of a slip velocity, would need to be put in to balance the drag. Additionally, the slip velocity is high when the orientation of the generating curve aligns with the swimming direction (axis of symmetry), and low otherwise. This is understandable as the slip velocity is constructed to be tangential to the generating curve, and a slip velocity perpendicular to the swimming direction generates little swimming velocity at the cost of additional power loss. This could be seen most clearly in the wavy shape of figure 3(d). Specifically, comparing the two points A and B marked in the panel, although point B has a larger radius than point A, the slip velocity of point B is lower because the orientation of the generating curve is almost perpendicular to the swimming direction.

Additionally, we note that the optimal slip velocity is proportional to the target swimming speed $U$ due to linearity of the Stokes equations. As a consequence, while the results only showcase micro-swimmers propelling themselves in the positive $\boldsymbol {e}_3$ direction, the optimal solution $u^{{S}*}$ for swimming in the opposite direction is merely a change of sign.

Micro-swimmers can be loosely classified as pushers that repel fluid from the body along the axis of symmetry, pullers that draw fluid to the body along the axis of symmetry or neutral swimmers that do not repel or draw fluid along the axis of symmetry (Lauga & Powers 2009). At first sight, the flow fields for all optimal swimmers studied here seem to be neutral swimmers. A closer look into the stresslet tensor $\boldsymbol {S}$, however, reveals a more interesting story. For an axisymmetric swimmer whose swimming direction is $\boldsymbol {e}_3$, the stresslet tensor could be simplified to $\boldsymbol {S} = S(\boldsymbol {e}_3\boldsymbol {e}_3 - \boldsymbol {I}/3)$, where $\boldsymbol {I}$ is the identity matrix. The sign of $S$ characterizes whether the swimmer is a pusher ($S<0$) or a puller ($S>0$).

It is easy to prove by contradiction that the optimal ‘front–back symmetric’ swimmers cannot be pushers nor pullers: flipping the direction of the slip velocity would make a pusher into a puller of the same shape with an equal (minimal) power loss, contradicting the unique solution guaranteed by the quadratic nature of the problem. However, the contradiction does not apply for ‘front–back asymmetric’ swimmers as flipping the swimming direction would essentially change the shape of the swimmer. In fact, the optimal ‘front–back asymmetric’ swimmers are not always neutral. For example, the stomatocyte shown in figure 3(f) is a puller where the stagnation point in the laboratory frame's flow field is in front of the micro-swimmer.

Conventionally, pusher and puller particles have been associated with ‘tail-actuated’ swimmers (e.g. spermatozoa) and ‘head-actuated’ swimmers (e.g. Chlamydomonas reinhardtii) respectively (Saintillan & Shelley 2015). It is, however, not immediately clear whether a micro-swimmer should be a pusher (tail actuated) or a puller (head actuated) to optimize its efficiency when given an arbitrary shape. Here, capitalizing on our earlier observation on the dependence of local geometry and optimal slip velocity, we propose a shape-based scalar metric $\mathbb {A}$ that can be used to predict whether the optimal swimmer for a given shape is a pusher or puller without the need of optimization. Simply speaking, $\mathbb {A}$ quantifies the relative ‘nominal actuation’ of the ‘head’ part and the ‘tail’ part of the swimmer based solely on the swimmer shape

(3.1)\begin{equation} \mathbb{A} = \log\left(\frac{\int_{\gamma_h} \boldsymbol{\tau}\boldsymbol{\cdot}\boldsymbol{e}_3\, x_1^{2}\, \mathrm{d}s/\int_{\gamma_h} x_1\, \mathrm{d}s}{\int_{\gamma_t} \boldsymbol{\tau}\boldsymbol{\cdot}\boldsymbol{e}_3\, x_1^{2}\, \mathrm{d}s/\int_{\gamma_t} x_1\, \mathrm{d}s}\right), \end{equation}

where the generating curve $\gamma$ is divided into two curves $\gamma = \gamma _h \cup \gamma _t$; $\gamma _h$ represents the generating curve of the head part and $\gamma _t$ represents the generating curve of the tail part. The numerator and denominator inside the logarithm function are the surface averages of the nominal actuation for the head and tail parts, respectively. The nominal actuation is stronger if the generating curve aligns with the swimming direction better (larger $\boldsymbol \tau \boldsymbol {\cdot }\boldsymbol {e}_3$), or if the material point is farther away from the axis of symmetry (larger $x_1$). For front–back symmetric shapes, we naturally divide $\gamma$ in the middle thus $\mathbb {A}\equiv 0$; for front–back asymmetric shapes, we divide $\gamma$ at the arc length where $x_1$ is the largest along the generating curve $s^{*}={\textrm {arg\,max}}_{s\in \gamma } x_1(s)$, or the average $s^{*}$ if ${\textrm {arg\,max}}$ returns more than one $s^{*}$. Positive $\mathbb {A}$ corresponds to a shape whose head part actuates stronger than its tail part, which indicates that the micro-swimmer is likely to be a puller; similarly negative $\mathbb {A}$ indicates that the micro-swimmer is likely to be a pusher.

The predictions based on $\mathbb {A}$ for various families of asymmetric shapes are shown in figure 4. Specifically, most of the shapes are correctly predicted as they lie in the first and the third quadrants; the ones that are misclassified, on the other hand, have close to zero $\mathbb {A}$ and $S$, which means the head and tail are similarly actuated and the optimal swimmers are close to neutral.

Figure 4. The value of $\mathbb {A}$ provides a simple prediction of the swimmer type. Swimmers with $\mathbb {A}<0$ are predicted to be pushers ($S<0$), and swimmers with $\mathbb {A}>0$ are predicted to be pullers ($S>0$). Swimmers in the first and third quadrants are correctly predicted. Shape families are shown in figure 6 and the generating curves are given in appendix C.

Next, we study the optimal active force density $\boldsymbol {f}$ corresponding to the same shapes. Its normal and tangential components are plotted in figure 5. We note that by the no-net-force condition (2.5), the power loss reduces to $P= 2{\rm \pi} \int _\gamma \boldsymbol {f}\boldsymbol {\cdot }({u}^{{S}}\boldsymbol \tau ) x_1\,\mathrm {d}s,$ implying that only the tangential component contributes to the power loss. The change in tangential forces as a function of arc length loosely resembles that of the optimal slip velocity, mediated by the local curvature of the generating curve. Qualitatively, a low local curvature suppresses the traction relative to the slip velocity, and a high local curvature amplifies it. Slip velocities scaled by their local curvatures are shown in black dotted curves for a reference.

Figure 5. Active force density on the swimmer surface as a function of arc length along the generating curve. Normal and tangential components of the force densities are depicted by blue and orange curves. Scaled optimal slip velocities $2u^{S*}\kappa R/U$ are shown in dotted curves, where $\kappa$ is the local curvature of the generating curve. Insets are the shapes of the corresponding swimmers.

In figure 6, we plot the minimal power loss as a function of the reduced volume for various shape families. The power loss is scaled by the minimal power loss of a spherical swimmer with the same volume $J_o = 12{\rm \pi} \mu RU^{2}$ with $R = (3V/4{\rm \pi} )^{1/3}$. The minimal power loss for prolate spheroids monotonically decreases as the shape gets more slender; in contrast, it is well known that the shape with the minimal fluid drag is one with approximately $2\,{:}\,1$ aspect ratio (Pironneau 1973). By slender body theory, the power loss of a prolate spheroid scales as $\sim \mu \alpha ^{2/3} U^{2}$, where $\alpha$ is the aspect ratio (see Leshansky et al. 2007). On the other hand, the minimal power loss for oblate spheroids grows rapidly as the reduced volume is increased. Shapes of the spherocylinder family behave similarly to the prolate spheroids, and converge to the spherical case when the length of the cylinder reduces to 0, as expected. It is, however, worth pointing out that spherocylinder costs more power loss than prolate spheroids with the same reduced volume; this relates to the fact that the peak slip velocity for spherocylinder is higher than that of the prolate spheroids (figure 3b,e). The stomatocyte family is constructed by ‘pulling’ the rim of the shape, effectively making the shape ‘taller’ and curls deeper and deeper inside. We find that ‘taller’ shapes require lower power loss for this shape family, which is qualitatively consistent with the spheroid family. Finally, we note that the power loss of the snowman family (two spheres attaching with each other) is quite robust to the relative sizes of the two spheres. The power loss is only approximately $25\,\%$ higher than that of a single sphere in the limit case where the two spheres are of the same size.

Figure 6. Scaled minimal power loss of different shape families, plotted against the reduced volume $\nu$. Example shapes are colour coded by the optimal slip velocity. The dotted line shows the approximation of power loss given by the slender body theory $P\sim \mu \alpha ^{2/3}U^{2}$ (Leshansky et al. 2007).

A few other examples that take more generic shapes are also shown in figure 6. The optimal slip velocities are coloured on their surfaces while their power loss is shown in the form of scatter points. The generating curves of these shapes are formed by spherical harmonics. We note that the optimal performance of shapes that appear similar can be very different. For example, the difference in power loss between examples 6 and 8 is approximately $150\,\%$ of the spherical swimmer, or $60\,\%$ of example 6. This result is a strong indicator that the slip velocity of the artificial swimmer, as well as its shape, must be carefully designed to achieve good performance.

We note that the minimal power loss for all the shape families considered here are bounded from below by the curve for prolate spheroids. However, since the current work does not optimize shape, whether the prolate spheroids are universally optimal remains to be tested.

4. Conclusions

In this work, we provided a solution procedure for the PDE-constrained optimization problem of finding the optimal slip profile on an axisymmetric micro-swimmer that minimizes the power loss required to maintain a target swimming speed. While it can be extended to other objective functions, we exploited the quadratic nature of the power loss functional in the control parameters to simplify and streamline the solution procedure. In the general case, an adjoint formulation and iterative optimization algorithms can be employed. Regardless of the formulation, however, the use of the boundary integral method to solve the Stokes equations greatly reduces the computational cost due to dimensionality reduction. Solving any of the examples presented in this work, for example, required only a few seconds on a standard laptop. Extending our procedure to fully three-dimensional (non-axisymmetric) shapes is straightforward; the key technical challenge is incorporating a high-order boundary integral solver, for which open-source codes are now available (e.g. see Gimbutas & Veerapaneni 2013).

Based on our numerical results, we came up with a heuristic metric that can classify the optimal swimming pattern for a given shape. It measures relative actuation of the ‘head’ and the ‘tail’ of the swimmer and predicts whether the optimal swimmer is head actuated (puller) or tail actuated (pusher). This metric could inform the early design of optimal slip for a given shape without the need for carrying out numerical optimization.

The optimization procedure developed in this work can directly be employed in the design pipeline of autophoretic particles. For example, in the case of diffusiophoresis, the computed optimal slip profile for a given shape can be used to formulate the chemical coating pattern of the phoretic particles. We acknowledge that the cost function for such optimization may need to be modified accordingly to reflect the chemical nature of the problem (Sabass & Seifert 2012). Another natural extension of this work is to relax the steady slip assumption and consider time-periodic squirming motion as done in Michelin & Lauga (2010). This would be particularly useful for studying the ciliary locomotion of micro-organisms with arbitrary shapes. Furthermore, building on the recent work of Bonnet, Liu & Veerapaneni (2020), we are developing solvers for the shape optimization problem of finding the most efficient micro-swimmer shapes under specified area, volume or other physical constraints.