2.1. Model

Consider an axisymmetric microswimmer whose boundary $\varGamma$ is obtained by rotating a generating curve $\gamma$ of length $\ell$ about the $\boldsymbol {e}_3$ axis, as shown in figure 1(a). We adopt the classic envelope model (Lighthill Reference Lighthill1952) and assume that the ciliary tips undergo time-periodic tangential movements along the generating curve. Let $s=\alpha (s_0,t)$ be the ciliary tip arc length coordinate on the generating curve $\gamma$ at time $t$ for a cilium rooted at $s_0$. The tangential slip velocity of this material point in its body-frame is thus

(2.1)\begin{equation} u^{S}(s, t)= u^{S}(\alpha(s_0, t), t)= \partial_t \alpha(s_0, t). \end{equation}

Figure 1.

(a) Schematic of the microswimmer geometry. The shape is assumed to be axisymmetric, obtained by rotating the generating curve $\gamma$ about the $\boldsymbol {e}_3$ axis. The tip of the cilium rooted at $s_0$ at time $t$ is given by $s=\alpha (s_0,t)$. (b) Illustration of the algorithm for computing the slip velocity at the quadrature points $u^{S}(s_q,t)$. We first compute the ‘tip’ position and the corresponding tip velocities (open blue circles) of cilia rooted at the $N_q$ quadrature points $s_q$ (closed blue circles). We then obtain the slip velocities at sample points uniformly distributed along the generating curve (open red squares) by a cubic interpolation. The slip velocity at any arc length (black curve) is then obtained by a high-order B-spline interpolation from the sample points. We have reduced the number of quadrature and sample points in this figure (compared to values used in the numerical experiments) to avoid visual clutter.

In addition to the time-periodic condition, the ciliary motion $\alpha$ needs to satisfy two more conditions to avoid singularity (Michelin & Lauga Reference Michelin and Lauga2010). First, the slip velocities should vanish at the poles

(2.2)\begin{equation} \alpha(0,t) = 0 \quad\text{and} \quad \alpha(\ell,t) = \ell, \quad \forall\, t\in\mathbb{R}^+, \end{equation}

and second, $\alpha$ should be a monotonic function, that is,

(2.3)\begin{equation} \partial_{s_0} \alpha(s_0,t) > 0, \quad \forall\, (s_0,t)\in [0,\ell]\times\mathbb{R}^+. \end{equation}

The last condition ensures the slip velocity is unique at any arc length $s$; in other words, crossing of cilia is forbidden. While in reality, cilia do cross, this condition is enforced to ensure validity of the continuum model.

In the viscous-dominated regime, the flow dynamics is described by the incompressible Stokes equations at every instance of time:

(2.4)\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, and $p$ and $\boldsymbol {u}$ are the fluid pressure and velocity fields, respectively. In the absence of external forces and imposed flow field, the far-field boundary condition is simply

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

The free-swimming microswimmer also needs to satisfy the no-net-force and no-net-torque conditions. Owing to the axisymmetric assumption, the no-net-torque condition is satisfied by construction, and the no-net-force condition is reduced to one scalar equation

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

where $x_1$ is the $\boldsymbol {e}_1$ component of $\boldsymbol {x}$, $\boldsymbol {f}$ is the active force density the swimmer applied to the fluid (negative to fluid traction) and $f_3$ is its $\boldsymbol {e}_3$ component.

Given any ciliary motion $\alpha (s_0,t)$ that satisfies (2.2) and (2.3), there is a unique tangential slip velocity ${u}^{{S}}(s,t)$ defined by (2.1). Such a slip velocity propels the microswimmer at a translational velocity $U(t)$ in the $\boldsymbol {e}_3$ direction, determined by (2.6). 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.7)\begin{equation} \boldsymbol{u}(\boldsymbol{x}(s),t) = {u}^{{S}}(s,t)\boldsymbol\tau(s)+{U}(t)\boldsymbol{e}_3, \end{equation}

where $\boldsymbol \tau$ is the unit tangent vector on $\gamma$. Thereby, the instantaneous power loss $P(t)$ can be written as

(2.8)\begin{align} P(t) &= \int_\varGamma \boldsymbol{f}(\boldsymbol{x},t) \boldsymbol{\cdot} \boldsymbol{u}(\boldsymbol{x}, t) \, \mathrm{d}\varGamma \nonumber\\ &= 2{\rm \pi}\left[ \int_\gamma \boldsymbol{f}(s,t) \boldsymbol{\cdot} \boldsymbol{\tau}(s) u^{S}(s,t)\, x_1 \, \mathrm{d}s + U(t)\int_\gamma \boldsymbol{f}(s,t) \boldsymbol{\cdot} \boldsymbol{e}_3\, x_1 \, \mathrm{d}s \right]. \end{align}

The second term on the right-hand side is zero provided that the no-net-force condition (2.6) is satisfied.

Following Lighthill (Reference Lighthill1952), we quantify the performance of the microswimmer by its swimming efficiency $\epsilon$, defined as

(2.9)\begin{equation} \epsilon = \frac{C_D \langle{U}\rangle^2}{\langle{P}\rangle}, \end{equation}

where ${P}=P(t)$ and ${U}=U(t)$ are the instantaneous power loss and swim speed, $\langle \cdot \rangle$ denotes the time average over one period and $C_D$ is the drag coefficient defined as the total drag force of towing a rigid body of the same shape at a unit speed along the $\boldsymbol {e}_3$ direction. The coefficient $C_D$ depends on the given shape $\gamma$ only; for example, $C_D = 6{\rm \pi} \mu a$ in the case of a spherical microswimmer with radius $a$.

In our simulations, we normalise the radius of the microswimmer to unity, and the period of the ciliary motion to $2{\rm \pi}$. It is worth noting that the swimming efficiency (2.9) is size and period independent, thanks to its dimensionless nature. The Reynolds number of a ciliated microswimmer of radius $100\ \mathrm {\mu } \text {m}$ and frequency 30 Hz submerged in water can be estimated as ${Re} \sim 10^{-4}$, confirming the applicability of the Stokes equations.

2.2. Numerical algorithm for solving the forward problem

Before stating the optimisation problem, we summarise our numerical solution procedure for the governing equations (2.4)–(2.7). By the quasi-static nature of the Stokes equation (2.4), the flow field $\boldsymbol {u}(\boldsymbol {x},t)$ can be solved independently at any given time, and the time averages can be found using standard numerical integration techniques (e.g. trapezoidal rule). We use a BIM at every time step. A similar BIM implementation was detailed in our recent work Guo et al. (Reference Guo, Zhu, Liu, Bonnet and Veerapaneni2021) in which we studied the optimisation of time-independent slip profiles. The main procedures are summarised below, with key equations highlighted in boxes.

We use the single-layer potential ansatz, which expresses the velocity as a convolution of an unknown density function $\boldsymbol {\mu }$ with Green's function for the Stokes equations:

(2.10)\begin{equation} \boldsymbol{u}(\boldsymbol{x}) = \frac{1}{8{\rm \pi}}\int_\varGamma \left(\frac{1}{|\boldsymbol{r}|} \boldsymbol{I} + \frac{\boldsymbol{r} \otimes \boldsymbol{r} }{|\boldsymbol{r}|^3}\right) \, \boldsymbol{\mu}(\boldsymbol{y}) \, \mathrm{d}\varGamma(\boldsymbol{y}), \quad\text{where}\ \boldsymbol{r} = \boldsymbol{x} - \boldsymbol{y}.\end{equation}

The force density can then be evaluated as a convolution of $\boldsymbol {\mu }$ with the (negative of) traction kernel:

(2.11)\begin{equation} \boldsymbol{f}(\boldsymbol{x}) =\frac{1}{2}\boldsymbol{\mu}\left(\boldsymbol{x}\right) + \frac{3}{4{\rm \pi}}\int_{\varGamma} \left(\frac{\boldsymbol{r} \otimes \boldsymbol{r} }{|\boldsymbol{r}|^5}\right) (\boldsymbol{r}\boldsymbol{\cdot}\boldsymbol{n}(\boldsymbol{x}))\boldsymbol{\mu}\left(\boldsymbol{y}\right)\mathrm{d}\varGamma\left(\boldsymbol{y}\right). \end{equation}

We convert these weakly singular boundary integrals into convolutions on the generating curve $\gamma$ by performing an analytic integration in the orthoradial direction, and then apply a high-order quadrature rule designed to handle the log singularity of the resulting kernels (Veerapaneni et al. Reference Veerapaneni, Gueyffier, Biros and Zorin2009). The Stokes flow problem defined at any time $t$ by (2.4)–(2.7) is then recast as the BIM system for the unknowns $\boldsymbol {\mu }$ and $U(t)$ obtained by substituting (2.10) into (2.7) and (2.11) into (2.6). The numerical solution method consists in discretising $\gamma$ into $N_p$ non-overlapping panels, each panel supporting the nodes of a 10-point Gaussian quadrature rule. The single-layer operator is approximated in Nyström fashion, by collocation at the $N_q=10 N_p$ quadrature nodes, while the values of $\boldsymbol {\mu }$ are sought at the same quadrature nodes. The resulting BIM system is

(2.12) \begin{equation} \begin{bmatrix} \mathcal{S} & -\mathcal{B} \\ \mathcal{C} & 0 \end{bmatrix} \begin{bmatrix} \boldsymbol{\mu}\\ {U(t)} \end{bmatrix} = \begin{bmatrix} {\boldsymbol{u}^{{S}}}\\ {0} \end{bmatrix}, \end{equation}

where the vectors $\boldsymbol {\mu }=\boldsymbol {\mu }(s_q,t)$ and $\boldsymbol {u}^{{S}}=\boldsymbol {u}^{{S}}(s_q,t)$ are the unknown density and the given slip velocity at all quadrature nodes $s_q$, $\mathcal {S}$ is the axisymmetric single-layer potential operator (which is fixed for a given shape $\gamma$), $\mathcal {B}$ is the column vector reproducing $\boldsymbol {e}_3$ at each quadrature node and $\mathcal {C}$ is the row vector, such that $\mathcal {C}[\boldsymbol {\mu }] = \int _\varGamma \boldsymbol {f}(\boldsymbol {x})\boldsymbol {\cdot }\boldsymbol {e}_3 \mathrm {d}\varGamma$ is the total traction force in the $\boldsymbol {e}_3$ direction.

The algorithm to obtain the slip velocity at the quadrature nodes at a given time $\boldsymbol {u}^{S}(s_q,t)$ is summarised in figure 1(b). Specifically, we start by computing the corresponding ciliary tip position $s=\alpha (s_q,t)$ and the slip velocity $u^{S}(s,t)$ from (2.1). These tip positions $s$ can be highly non-uniform, depending on the form of $\alpha$, which can be difficult for the forward solver. To circumvent this difficulty and to find a smooth representation of the slip velocities on the quadrature points, we first find the slip velocities at $N_s$ sample points uniformly distributed along the generating curve by interpolating $u^{S}(s,t)$ (we use the routine PCHIP in MATLAB); the slip velocities at the quadrature nodes $u^{S}(s_q,t)$ are then in turn interpolated from the $N_s$ sample points using high-order B-spline bases. An alternative approach could be to follow the position and the slip velocity of each material point. In other words, one can use $\boldsymbol {u}^{S}(s,t)$ directly on the right-hand side of (2.12), which will bypass the interpolation steps mentioned above. However, it requires reassembly of the matrix $\mathcal {S}$ at every time step, significantly increasing the computational cost.

2.3. Optimisation problem

The goal of this work is to find the optimal ciliary motion for a given arbitrary axisymmetric shape; that is, the ciliary motion $\alpha ^{\star }(s_0,t)$ that maximises the swimming efficiency $\epsilon$:

(2.13) \begin{equation} \alpha^{{\star}}=\mathop{\textrm{arg max}}\limits_{\alpha\in\mathcal{A}}\epsilon(\alpha), \end{equation}

where $\mathcal {A}$ is the space of all possible time-periodic ciliary motion satisfying (2.2) and (2.3). It is, however, not easy to define and manipulate finite-dimensional parametrisations of $\alpha$ that remain in that space. To circumvent this difficulty, we follow the ideas in Michelin & Lauga (Reference Michelin and Lauga2010) and represent $\alpha$ in terms of a time-periodic function $\psi (x,t)$, such that

(2.14)\begin{equation} \alpha(s_0,\psi) = \frac{\ell\int_0^{s_0}{[\psi(x,t)]^2\,\mathrm{d}\kern0.04em x}}{\int_0^{\ell}{[\psi(x,t)]^2\,\mathrm{d}\kern0.04em x}}, \end{equation}

where $\ell$ is the total length of the generating curve $\gamma$. Note that $\alpha$ is also (implicitly) a function of time $t$, through $\psi = \psi (x,t)$. It is easy to verify that $\alpha$ given by (2.14) satisfies the boundary conditions (2.2) and the monotonicity requirement (2.3) for any choice of $\psi$. Conversely, for any $\alpha$ satisfying (2.2) and (2.3), there is at least one $\psi$ that provides $\alpha$. As a result, the optimisation problem is recast as finding

(2.15) \begin{equation} \psi^{{\star}} = \mathop{\textrm{arg max}}\limits_{\psi} \epsilon(\psi), \end{equation}

where $\psi (\cdot ,t)$ is only required to be square-integrable over $[0,\ell ]$ for any $t$.

We use a quasi-Newton Broyden–Fletcher–Goldfarb–Shanno (BFGS) method (Nocedal & Wright Reference Nocedal and Wright2006) to optimise the ciliary motion via $\psi$, which requires repeated evaluations of efficiency sensitivities with respect to perturbations of $\psi$. The sensitivities of power loss and swim speed are derived using an adjoint-based method, while the efficiency sensitivity is found using the quotient rule thereafter. The adjoint-based method exhibits a great advantage over the traditional finite-difference method when finding the sensitivities, because regardless of the dimension of the parameter space, the objective derivatives with respect to all design parameters can here be evaluated on the basis of one solve of the forward problem for each given ciliary motion $\alpha$. The derivations are detailed below.

2.4. Sensitivity analysis

We start by finding the sensitivities in terms of the slip profile $u^{S}$. The sensitivities in terms of the auxiliary unknown $\psi$ will be found subsequently by a change of variable. As the concept of an adjoint solution, in general, rests on duality considerations, we recast the forward-flow problem in weak form for the purpose of finding the sought sensitivities of power loss and swim speed, even though the numerical forward-solution method used in this work does not directly exploit that weak form. Specifically, we recast the forward problem (2.4)–(2.7) in mixed weak form (see e.g. Brezzi & Fortin Reference Brezzi and Fortin1991, chapter 6). That is, to find $(\boldsymbol {u}, p, \boldsymbol {f}, U) \in \boldsymbol {\mathcal {V}}\times \mathcal {P}\times \boldsymbol {\mathcal {F}}\times \mathbb {R},$ such that

(2.16a)$$\begin{gather} a(\boldsymbol{u}, \boldsymbol{v}) - b(\boldsymbol{v},p) - b(\boldsymbol{u},q) - \langle\, \boldsymbol{f}, \boldsymbol{v}\rangle_{\varGamma} = 0 \quad \forall\, (\boldsymbol{v},q) \in \boldsymbol{\mathcal{V}}\times\mathcal{P}, \end{gather}$$

(2.16b)$$\begin{gather}\langle \boldsymbol{g}, \boldsymbol{e}_3\rangle_\varGamma U + \langle \boldsymbol{g}, u^{S} \boldsymbol{\tau} \rangle_\varGamma - \langle \boldsymbol{g}, \boldsymbol{u}\rangle_\varGamma = 0 \quad \forall \boldsymbol{g} \in \boldsymbol{\mathcal{F}}, \end{gather}$$

(2.16c)$$\begin{gather}\langle\, \boldsymbol{f}, \boldsymbol{e}_3\rangle_\varGamma = 0, \end{gather}$$

where the bilinear forms $a$ and $b$ are defined by

(2.17a,b)\begin{equation} a(\boldsymbol{u}, \boldsymbol{v}) := \int_\varOmega 2\mu \boldsymbol{D}[\boldsymbol{u}]:\boldsymbol{D}[\boldsymbol{v}]\, \mathrm{d} V, \quad b(\boldsymbol{v},q) := \int_\varOmega q\, \text{div}\,\boldsymbol{v}\, \mathrm{d}V, \end{equation}

and $\boldsymbol {D}[\boldsymbol {u}] := (\boldsymbol {\nabla } \boldsymbol {u}+\boldsymbol {\nabla }^{\textrm {T}}\boldsymbol {u})/2$ is the strain-rate tensor; we use $\langle \cdot , \cdot \rangle _\varGamma$ as shorthand for the inner product on $\varGamma$. For example, $\langle \boldsymbol {f}, \boldsymbol {v}\rangle _\varGamma = \int _\varGamma \boldsymbol {f} \cdot \boldsymbol {v}\, \mathrm {d}\varGamma$. Similarly, with a slight abuse of notation, the power loss functional can be written as $P(u^{S}) := \langle \boldsymbol {f}, u^{S}\boldsymbol {\tau } + U\boldsymbol {e}_3 \rangle _\varGamma$, where $U := U(u^{S})$ is the swim speed functional.

The Dirichlet boundary condition (2.7) is (weakly) enforced explicitly through (2.16b), rather than being embedded in the velocity solution space $\boldsymbol {\mathcal {V}}$, as this will facilitate the derivation of slip-derivative identities; this is in fact our motivation for using the mixed weak form (2.16). Condition (2.16c) is the no-net-force condition (2.6).

First-order sensitivities of functionals at $u^{S}$ are defined as directional derivatives, by considering perturbations of $u^{S}$ of the form

(2.18)\begin{equation} u^{S}_{\eta} = u^{S} + \eta \nu \end{equation}

for some $\nu$ in the slip velocity space and $\eta \in \mathbb {R}$. Then, the directional (or Gâteaux) derivative of a functional $J(u^{S})$ in the direction $\nu$, denoted by $J'(u^{S};\nu )$, is defined as

(2.19)\begin{equation} J'(u^{S};\nu) = \lim_{\eta\to0} \frac{1}{\eta}\left( J[u_{\eta}^{S}]-J[u^{S}] \right). \end{equation}

For the power-loss functional, we obtain (because the derivative of ${u}^{S}$ in the above sense is $\nu$)

(2.20a,b)\begin{equation} P'(u^{S};\nu) = \langle\, \boldsymbol{f}',{u}^{S}\boldsymbol{\tau} + U\boldsymbol{e}_3 \rangle_{\varGamma} + \langle\, \boldsymbol{f},\nu\boldsymbol{\tau} \rangle_{\varGamma} + \langle\, \boldsymbol{f},\boldsymbol{e}_3 \rangle_{\varGamma}U', \end{equation}

where $\boldsymbol {f}'$ and $U'$ are the derivatives of the active force $\boldsymbol {f}$ and swim speed $U$ solving problem (2.16), considered as functionals on the slip velocity $u^{S}$:

(2.21)\begin{equation} \boldsymbol{f}' = \lim_{\eta\to0} \frac{1}{\eta}\left( \boldsymbol{f}[u_{\eta}^{S}]-\boldsymbol{f}[u^{S}] \right), \quad U' = \lim_{\eta\to0} \frac{1}{\eta}\left( U[u_{\eta}^{S}]-U[u^{S}] \right). \end{equation}

Differentiating the weak formulation (2.16) of the forward problem with respect to $u^{S}$ leads to the weak formulation of the governing problem for the derivatives $(\boldsymbol {u}', \boldsymbol {f}', p', U')$ of the solution $(\boldsymbol {u}, \boldsymbol {f}, p, U)$:

(2.22a)$$\begin{gather} a(\boldsymbol{u}',\boldsymbol{v}) - b(\boldsymbol{u}',q) - b(\boldsymbol{v},p') - \langle\, \boldsymbol{f}',\boldsymbol{v} \rangle_{\varGamma} =0 \quad\, \forall(\boldsymbol{v},q)\in\boldsymbol{\mathcal{V}}\times\mathcal{P}, \end{gather}$$

(2.22b)$$\begin{gather}\langle \nu\boldsymbol{\tau},\boldsymbol{g}\rangle_{\varGamma} + U'\langle \boldsymbol{e}_3,\boldsymbol{g}\rangle_{\varGamma} - \langle \boldsymbol{u}',\boldsymbol{g} \rangle_{\varGamma} = 0 \quad \forall\, \boldsymbol{g}\in\boldsymbol{\mathcal{F}}, \end{gather}$$

(2.22c)$$\begin{gather}\langle\, \boldsymbol{f}',\boldsymbol{e}_3 \rangle_{\varGamma} = 0 . \end{gather}$$

Here we have assumed without loss of generality that the test functions in (2.16) verify $\boldsymbol {v}' = \boldsymbol {0}$, $\boldsymbol {g}' = \boldsymbol {0}$, and $q' = 0$, which is made possible by the absence of boundary constraints in $\boldsymbol {\mathcal {V}}$.

At first glance, evaluating $P'(u^{S};\nu )$ in a given perturbation $\nu$ appears to rely on solving the derivative problem (2.22). However, a more effective approach allows us to bypass the actual evaluation of ${\boldsymbol {f}}'$. Let the adjoint problem be defined by

(2.23a)$$\begin{gather} a(\hat{\boldsymbol{u}},\boldsymbol{v}) - b(\hat{\boldsymbol{u}},q) - b(\boldsymbol{v},\hat{p}) - \langle\hat{\boldsymbol{f}},\boldsymbol{v}\rangle_{\varGamma} = 0 \quad \forall(\boldsymbol{v},q)\in\boldsymbol{\mathcal{V}}\times\mathcal{P}, \end{gather}$$

(2.23b)$$\begin{gather}\langle \boldsymbol{e}_3,\boldsymbol{g}\rangle_{\varGamma} - \langle \hat{\boldsymbol{u}},\boldsymbol{g}\rangle_{\varGamma} = 0 \quad \forall\boldsymbol{g}\in\boldsymbol{\mathcal{F}}, \end{gather}$$

i.e. $(\hat {\boldsymbol {u}},\hat {p})$ are the flow variables induced by prescribing a unit velocity $\boldsymbol {e}_3$ on $\varGamma$. For later convenience, we let $F_0$ denote the (non-zero) net force exerted on $\varGamma$ by the adjoint flow:

(2.24)\begin{equation} F_0 := \langle\, \hat{\boldsymbol{f}},\boldsymbol{e}_3 \rangle_{\varGamma}. \end{equation}

Problem (2.23) in strong form is defined by (2.4)–(2.7) with $U=1,\,u^{S}=0$. In fact, $F_0$ takes the same value as the drag coefficient $C_D$ in (2.9).

Then, combining the derivative problem (2.22) with the forward problem (2.16) or the adjoint problem (2.23) with appropriate choices of test functions allows us to derive expressions of $P'(u^{S};\nu )$ and $U'(u^{S};\nu )$ which do not involve the forward solution derivatives.

Specifically, we set the test functions to $(\boldsymbol {v},q,\boldsymbol {g})=(\boldsymbol {u}',p',\boldsymbol {f}')$ in ((2.16a) and (2.16b)) of the forward problem and $(\boldsymbol {v},q,\boldsymbol {g})=({\boldsymbol {u}},{p},{\boldsymbol {f}})$ in ((2.22a) and (2.22b)) of the derivative problem. Then, the combination (2.22a) + (2.22b) $-$ (2.16a) $-$ (2.16b) is evaluated, to obtain

(2.25)\begin{equation} \langle\, \boldsymbol{f}',{u}^{S}\boldsymbol{\tau} + U\boldsymbol{e}_3 \rangle_{\varGamma} = \langle\, \boldsymbol{f},\nu\boldsymbol{\tau} \rangle_{\varGamma} + \langle\, \boldsymbol{f},\boldsymbol{e}_3 \rangle_{\varGamma}U'. \end{equation}

Substituting (2.25) into (2.20), and recalling the no-net-force condition (2.6), we have

(2.26)\begin{equation} \boxed{ P'(u^{S};\nu) = 2\langle\, \boldsymbol{f},\nu\boldsymbol{\tau} \rangle_{\varGamma} ={4{\rm \pi}}\int_\gamma (\boldsymbol{f} \boldsymbol{\cdot} \boldsymbol{\tau}) \, \nu x_1 \, \mathrm{d}s.} \end{equation}

Likewise, setting the test functions to $(\boldsymbol {v},q,\boldsymbol {g})=(\boldsymbol {u}',p',\boldsymbol {f}')$ in the adjoint problem (2.23) and $(\boldsymbol {v},q,\boldsymbol {g})=(\hat {\boldsymbol {u}},\hat {p},\hat {\boldsymbol {f}})$ in ((2.22a) and (2.22b)) of the derivative problem (2.22), then evaluating (2.22a) + (2.22b) $-$ (2.23a) $-$ (2.23b), yields

(2.27)\begin{equation} 0 = \langle\, \hat{\boldsymbol{f}},\nu\boldsymbol{\tau} \rangle_{\varGamma} + \langle\, \hat{\boldsymbol{f}},U'\boldsymbol{e}_3 \rangle_{\varGamma} - \langle\, \boldsymbol{f}',\boldsymbol{e}_3 \rangle_{\varGamma} = \langle\, \hat{\boldsymbol{f}},\nu\boldsymbol{\tau} \rangle_{\varGamma} + F_0U'. \end{equation}

Note that $\langle \boldsymbol {f}',\boldsymbol {e}_3 \rangle _{\varGamma } = 0$ according to (2.22c). Rearranging terms in (2.27), we have

(2.28)\begin{equation} \boxed{ U'(u^{S};\nu) ={-}\frac{1}{F_0}\langle\, \hat{\boldsymbol{f}},\nu\boldsymbol{\tau} \rangle_{\varGamma} ={-} \frac{2{\rm \pi}}{F_0} \int_\gamma (\hat{\boldsymbol{f}}\boldsymbol{\cdot} \boldsymbol{\tau}) \, \nu x_1\, \mathrm{d}s.} \end{equation}

The sensitivity formulas (2.26) and (2.28) are not practically applicable in this form to the current optimisation problem, because the constraints (2.2) and (2.3) are not easy to enforce on parametrisations of the unknown slip profiles $u^{S}$. For this reason, we rewrite the quantities of interest as functionals of $\psi$, and the connection between $\psi$ and $\alpha$ is given by (2.14). Specifically, the slip profile is

(2.29)\begin{equation} u^{S}(s,t) = \partial_{t}\alpha(s_0,\psi) = \partial_\psi\alpha(s_0,\psi;\dot{\psi}) = \partial_\psi\alpha\left( \beta(s,\psi),\psi;\dot{\psi} \right) = v^{S}(s,\psi), \end{equation}

where $\dot {\psi } := \partial _t \psi$ and $\beta (s,\psi )$ is the inverse function of $\alpha$, i.e. $s_0 = \beta (s,\psi )$. The average power-loss and swim-speed functionals are written as

(2.30)\begin{equation} \langle\mathbb{P}\rangle(\psi) := \langle P \rangle(u^{S}), \quad \langle \mathbb{U}\rangle(\psi):= \langle U\rangle(u^{S}) \quad \text{with}\ u^{S}(s,t) = v^{S}(s,\psi). \end{equation}

On applying the change of variables $s = \alpha (s_0, \psi )$ in the integrals (2.26) and (2.28) and averaged over one period, we obtain

(2.31)$$\begin{gather} \boxed{\langle\mathbb{P}\rangle'(\psi;\hat{\psi}) = 2 \int_0^{2{\rm \pi}} \int_\gamma \boldsymbol{f}(\alpha)\boldsymbol{\cdot}\boldsymbol{\tau}(\alpha)\,x_1(\alpha)\, v^{S}{}'(s,\psi;\hat{\psi}) \, \partial_s\alpha \, \mathrm{d}s_0 \,\mathrm{d}t,} \end{gather}$$

(2.32)$$\begin{gather}\boxed{ \langle \mathbb{U}\rangle'(\psi;\hat{\psi}) ={-}\frac{1}{F_0} \int_0^{2{\rm \pi}} \int_\gamma \hat{\boldsymbol{f}}(\alpha)\boldsymbol{\cdot}\boldsymbol{\tau}(\alpha)\,x_1(\alpha)\,v^{S}{}'(s,\psi;\hat{\psi}) \, \partial_s\alpha\, \mathrm{d}s_0 \, \mathrm{d}t,} \end{gather}$$

where $v^{S}{}'(s, \psi ; \hat {\psi })$ is the directional derivative of $u^{S}$ with respect to $\psi$ and in the direction $\hat {\psi }$. Specifically, we can show that

(2.33)$$\begin{gather} v^{S}{}'(s,\psi;\hat{\psi}) \, \partial_s\alpha(s_0,\psi) \mathrm{d}s_0 = \left\{ \partial_s\alpha(s_0,\psi)\,\left[ \partial^2_\psi \alpha \left( s_0,\psi;\hat{\psi},\dot{\psi}\right) + \partial_\psi \alpha \left( s_0,\psi;\dot{\hat{\psi}} \right)\,\right] \right. \nonumber\\ \left. -\partial_{\psi s}\alpha\left( s_0,\psi;\dot{\psi} \right)\; \partial_{\psi}\alpha\left( s_0,\psi;\hat{\psi} \right) \right\} \,\mathrm{d}s_0. \end{gather}$$

The derivation and the explicit expression of each term in (2.33) are given in Appendix A. Finally, the efficiency sensitivity in terms of $\psi$ readily follows by the quotient rule

(2.34)\begin{equation} \epsilon'(\psi;\hat\psi) = C_D\frac{2\langle\mathbb{U}\rangle \langle\mathbb{U}\rangle'\langle\mathbb{P}\rangle - \langle\mathbb{U}\rangle^2\langle\mathbb{P}\rangle'}{\langle\mathbb{P}\rangle^2}. \end{equation}

2.5. Constraints on surface displacement

The unconstrained optimisation problem (2.15) introduced previously has the tendency to converge to unphysical/unrealistic strokes, where each cilium effectively ‘covers’ the entire generating curve. For a more realistic model, we should add a constraint on the length of the cilium. To this end, and again following Michelin & Lauga (Reference Michelin and Lauga2010), we replace the initial unconstrained optimisation problem (2.15) with the penalised optimisation problem

(2.35a,b)\begin{equation} \psi^{{\star}} = \textrm{arg max}_{\psi} E(\psi), \quad E(\psi) = \epsilon(\psi) - C(\psi), \end{equation}

where the (non-negative) penalty term $C(\psi )$, defined as

(2.36)\begin{equation} C(\psi) = \int_0^{\ell} H({A}(\psi) - c) \,\mathrm{d}s_0, \end{equation}

serves to incorporate the kinematic constraint $A(\psi )\leq c$ in the optimisation problem. The functional ${A}(\psi )$ in (2.36) is a measure of the amplitude of the displacement of individual material points for the stroke (through $\alpha$) and $c$ is a threshold parameter to bound $A(\psi )$ (a smaller $c$ corresponding to a stricter constraint). We use $H$ as a smooth non-negative penalty function defined by

(2.37)\begin{equation} H(u) = \varLambda_1\left[1+\tanh \left(\varLambda_2{u}\right)\right]u^2, \end{equation}

which for a large enough $\varLambda _2$ approximates $u\mapsto 2\varLambda _1 u^2 Y(u)$ ($Y$ being the Heaviside unit step function). The multiplicative parameter $\varLambda _1$ then serves to tune the severity of the penalty incurred by violations of the constraint $A(\psi )\leq c$. We use $\varLambda _1 = 10^4$ and $\varLambda _2 = 10^{4}$ in our numerical simulations unless otherwise mentioned. The optimisation results are not sensitive to the choice of $\varLambda _1$ and $\varLambda _2$. A small caveat of the penalty function (2.37) is that it has a (small) bump at $\varLambda _2 u \approx -1.109$. This bump can occasionally trap the optimisations into local extrema that have significantly lower efficiencies, depending on the initial guesses. Perturbing $\varLambda _2$ for such cases helps to alleviate the problem.

The most relevant physical definition of ${A}$ would be the actual displacement amplitude of an individual point, i.e. $\Delta s= [\alpha _{max}(s_0) - \alpha _{min}(s_0)]/2$. However, the strong nonlinearity of this measure is not appropriate for the computation of the gradient. Following Michelin & Lauga (Reference Michelin and Lauga2010), we measure the displacement by its variance in time:

(2.38)\begin{equation} {A}(\psi) = \langle (\alpha(s_0,\psi) - \langle\alpha\rangle(s_0))^2\rangle. \end{equation}

The maximum displacement $\Delta s_{max} = \max _{s_0}(\Delta s)$ will be found post-optimisation for the optimal ciliary motion $\alpha ^{\star }$ to better illustrate our results in § 3.

Like the initial problem (2.15), the penalised problem (2.35a,b) is solvable using unconstrained optimisation methods, and we again adopt a quasi-Newton BFGS algorithm to optimise the ciliary motion. Applying the chain rule to the penalty functional $C(\psi )$, we obtain the derivative of the penalty term in the direction of $\hat \psi$ as

(2.39)\begin{equation} C'(\psi; \hat\psi) = \int_0^\ell H'({A}(\psi) - c) {A}'(\psi; \hat\psi) \,\mathrm{d}s_0. \end{equation}

The derivative of the penalised objective functional $E(\psi )$ is therefore

(2.40)\begin{equation} E'(\psi;\hat{\psi}) = \epsilon'(\psi;\hat{\psi}) - C'(\psi;\hat{\psi}), \end{equation}

where $\epsilon '$ and $C'$ are given by (2.34) and (2.39), respectively.