Dynamics and optimal actuation of a three-sphere low-Reynolds-number swimmer with muscle-like arms

TL;DR

This paper models a three-sphere microswimmer with muscle-like arms, analyzing its dynamics, control strategies, and optimizing actuation for improved swimming performance at low Reynolds numbers.

Contribution

It introduces a novel model of a microswimmer with muscle-inspired arms, deriving its equations, analyzing solutions, and exploring optimal actuation strategies.

Findings

Analytical and numerical solutions of the swimmer's dynamics.

Optimal actuation strategies for performance enhancement.

Control-theoretic analysis of the swimmer's governing equations.

Abstract

The three-sphere swimmer by Najafi and Golestanian is composed of three spheres connected by two arms. The case in which the swimmer can control the lengths of the two arms has been studied in detail. Here we study a variation of the model in which the swimmer's arms are constructed according to Hill's model of muscular contraction. The swimmer is able to control the tension developed in the active components of the arms. The two shape parameters and the tensions acting on the two arms are then obtained by solving a system of ordinary differential equations. We study the qualitative properties of the solutions, compute analytically their leading order approximation and compare them with numerical simulations. We also formulate and solve some optimisation problems, aimed at finding the actuation strategies maximising performance, for various performance measures. Finally, we discuss the structure of the governing equations of our microswimmers from the point of view of control theory. We show that our systems are control affine systems with drift.

Tables

Table 1. Table 1: Values of the parameters used in the numerical simulations.

$J_s$	$J_p$	$K$	$\alpha$	${\tilde{a}}_1$	${\tilde{b}}_1$	${\tilde{a}}_2$	${\tilde{b}}_2$	$ϵ$
4	3	2	0.1	1	0	0.25	0	0.7

Table 2. Table 2: Values of the parameters used in the numerical simulations.

$J_s$	$J_p$	$K$	$H$	$\alpha$	$\tilde{a}$	$\tilde{b}$	$ϵ$
4	3	2	5	0.1	1	0	0.7

Equations

$$v_1$$

$$v_2$$

$$v_3$$

$$f_1+f_2+f_3=0.$$

$${\dot{L}}_1=v_2-v_1,$$

$${\dot{L}}_2=v_3-v_2.$$

$$x:=\frac{1}{3}(x_1+x_2+x_3).$$

$$\dot{x}=\frac{1}{3}(v_1+v_2+v_3).$$

$$\dot{x}=\Big{(}\frac{1}{L_{1}+L_{2}}-\frac{1}{L_{2}}\Big{)}\frac{f_{1}}{12\pi\mu}+\Big{(}\frac{1}{L_{1}+L_{2}}-\frac{1}{L_{1}}\Big{)}\frac{f_{3}}{12\pi\mu}\,.$$

$$\dot{x}=\frac{a}{6}\Big{[}\Big{(}\frac{\dot{L}_{2}-\dot{L}_{1}}{L_{2}+L_{1}}\Big{)}+2\Big{(}\frac{\dot{L}_{1}}{L_{2}}-\frac{\dot{L}_{2}}{L_{1}}\Big{)}+\frac{\dot{L}_{2}}{L_{2}}-\frac{\dot{L}_{1}}{L_{1}}\Big{]}\,.$$

$$\dot{x}=\frac{a}{6}\Big{[}\Big{(}\frac{\dot{L}_{2}-\dot{L}_{1}}{L_{2}+L_{1}}\Big{)}+2\Big{(}\frac{\dot{L}_{1}}{L_{2}}-\frac{\dot{L}_{2}}{L_{1}}\Big{)}\Big{]}\,,$$

$$\{\begin{array}{ll}& L_1=l_1+U_1\\ & L_2=l_2+U_2\\ & U_i/l_j<<1.\end{array}$$

$$\dot{\xi}=\frac{a}{6}\Big{[}\frac{\dot{U_{2}}-\dot{U_{1}}}{l_{1}+l_{2}}-\frac{(\dot{U_{2}}-\dot{U_{1}})(U_{1}+U_{2})}{(l_{1}+l_{2})^{2}}+\frac{2\dot{U_{1}}}{l_{2}}-\frac{2\dot{U_{1}}U_{2}}{l_{2}^{2}}-\frac{2\dot{U_{2}}}{l_{1}}+\frac{2\dot{U_{2}}U_{1}}{l_{1}^{2}}\Big{]}\,.$$

$$l\Delta x=\frac{a}{6}\Big{[}\frac{1}{l_{1}^{2}}+\frac{1}{l_{2}^{2}}-\frac{1}{(l_{1}+l_{2})^{2}}\Big{]}\int_{0}^{T}(U_{1}\dot{U_{2}}-\dot{U_{1}}U_{2})dt\,.$$

$$-f_1+T_1=0$$

$$-f_2-T_1+T_2=0$$

$$-f_3-T_2=0.$$

$$v_1$$

$$v_2$$

$$v_3$$

$${\dot{L}}_1$$

$${\dot{L}}_2$$

$$T_1$$

$$T_2$$

$${\dot{T}}_1$$

$${\dot{T}}_2$$

$${\Lambda }_1^{\prime }$$

$${\Lambda }_2^{\prime }$$

$${\xi }^{\prime }$$

$${\tau }_1^{\prime }$$

$${\tau }_2^{\prime }$$

$$J_{s,i}:=\frac{k_{s,i}}{\omega B_i},J_{p,i}:=\frac{k_{p,i}}{\omega B_i},K_i:=\frac{k_{s,i}}{\mu \omega l},$$

$${\tau }_1=H({\Lambda }_1-{\lambda }_1),$$

$${\lambda }_1={\lambda }_2=1/2$$

$$K_1=K_2=:K$$

$$J_{s,1}=J_{s,2}=:J_s$$

$$J_{p,1}=J_{p,2}=:J_p.$$

$${\sigma }_i=ϵ{\tilde{\sigma }}_i,$$

$$Y^{\prime }=f(Y)+ϵg(t^{\ast }),$$

Full text

Introduction

The three-sphere swimmer by Najafi and Golestanian [1, 2] is a cornerstone in the literature on low Reynolds number swimming. It is composed of three spheres connected by two arms. The presence of two shape parameters (the lengths of the two arms) allows to perform periodic shape changes which are not invariant under time reversal. This is the key to beat Purcell’s famous scallop theorem [3], which says that low Reynolds number swimmers cannot achieve locomotion through a reciprocal shape change. Another famous example of low Reynolds number swimmer with two shape parameters is the three-link swimmer by Purcell [3]. These two model swimmers have been studied extensively (see for example [4, 5, 6, 7, 8, 9, 10, 11, 12, 13]).

In the works [1, 2] the authors study a three-sphere swimmer which is able to perform periodic shape changes by controlling the lengths of the two arms. Real swimmers however are not capable of directly controlling their shape. Their shape change is generated by a complex interplay between different elements: forces generated by the swimmer through internal mechanisms, elastic properties of the swimmer’s body, and forces due to interactions with the surrounding viscous fluid. The swimmer is able to control only the first of these components. For this reason it interesting to study minimal swimmers which are not capable of directly controlling their shape. An interesting work going in this direction is [14], in which the authors consider a three-sphere swimmer whose arms can generate forces thanks to a system of molecular motors and elastic elements. Another possibility is to assume that one of the two shape parameters can be controlled while the other one is driven by a passive elastic spring. This was done in [15] for the three-link swimmer and in [16] for the three-sphere swimmer.

In the present work we study a three-sphere swimmer whose arms are viscoelastic structures constructed according to Hill’s active state muscle model. This model, which is very useful for computing the mechanical behaviour of muscles, includes passive elastic elements, a mechanism of viscous damping, and an active component. The swimmer is able to control the force generated by the active component. The geometry of the shape change is not known a priori and will have to be computed by solving a system of ODEs. In the first section we introduce our model and derive the ODEs governing the system. In the second section we study the qualitative properties of solutions, compute analytically the leading order term of their asymptotic expansion, present the results of numerical simulations, and study some optimmization problems. In the third section we study a variation of our model, in which one of the two muscle-like arms is replaced by a passive elastic spring.

1 Problem formulation

The aim of this section is to introduce the thee-sphere swimmer with muscle-like arms and to obtain a system of ODEs governing its dynamics. In the first subsection we fix the notations and recall some basic facts about the three-sphere swimmer. In the second subsection we assume that the forces acting through the two arms are known functions of time and shape and obtain the equations governing the evolution of the two shape parameters. In the third subsection we model the arms according to Hill’s muscle model and obtain a system of ODEs governing the forces acting across each arm, the shape of the swimmer and its displacement. In the fourth subsection we write the problem in non-dimensional form. In the fifth and last subsection we consider a variation of the problem in which one of the two muscle-like arms is replaced by a passive elastic spring.

1.1 Three-sphere swimmer

Let us consider a three-sphere swimmer as in Figure 1. Let $a$ be the radius of the spheres, and $L_{1},L_{2}$ the lengths of the two arms. We indicate with $\mu$ the dynamic viscosity of the fluid and with $v_{i}$ and $f_{i}$ respectively the velocity of sphere $i$ and the force that this sphere exerts on the fluid, projected on the unit vector $\mathbf{e}$. By using the Oseen tensor and the approximation $\frac{a}{L_{i}}<<1$ we obtain the following linear relations between forces and velocities

[TABLE]

Due to Newton’s third law of motion, the force exerted by the fluid on the $i$-th sphere is $-f_{i}$. Therefore the force balance equation for the swimmer is

[TABLE]

The geometry of the system implies the following kinematic relations

[TABLE]

Let us indicate with $x_{i}$ the position of the $i$-th sphere on axis corresponding to the unit vector $\mathbf{e}$. We indicate with $x$ the mean point of the three spheres, namely,

[TABLE]

The translational velocity of the swimmer is the velocity of the point $x$. Obviously

[TABLE]

Using equations (1)-(4) we can show that

[TABLE]

Now we can use equations (1)-(6) to express $f_{1}$ and $f_{3}$ in terms of $L_{1},L_{2},\dot{L}_{1},\dot{L}_{2}$. If we plug the resulting expressions into equation (9) and keep only the leading order terms in $a/L_{j}$ we obtain

[TABLE]

Suppose that $L_{1}$ and $L_{2}$ are periodic functions. In this case the terms $\dot{L}_{j}/L_{j}$ average to zero in a full swimming cycle, because they are derivatives of $\log(L_{j})$. If we neglect these terms we obtain the following formula

[TABLE]

which can be used instead of (10) to compute the net displacement in one period when the swimmer performs a periodic shape change.

Now we consider the case of small deformations

[TABLE]

We assume that the deformations are periodic with period $T$. We would like to compute the net displacement in one period to leading order in the amplitude of deformations. Let us set $l:=l_{1}+l_{2}$. $l$ is the body length of the swimmer. If we expand equation (11) and retain only the leading order terms in $U_{i}/l_{j}$ we get

[TABLE]

The terms $\dot{U_{j}},\dot{U_{j}}U_{j},\dot{U_{1}}U_{2}+U_{1}\dot{U_{2}}$ give zero when integrated between [math] and $T$. So the net displacement in one period $l\Delta x$ is given, to leading order in $U_{i}/l_{j}$, by the following formula

[TABLE]

1.2 TSS with assigned tensions

Let us indicate with $T_{1}$ and $T_{2}$ the tension on the tail rod and on the front rod respectively. We assume that $T_{1}$ and $T_{2}$ are known functions of time and shape, namely, $T_{i}=T_{i}(t,L_{1},L_{2})$, $i=1,2$. The equations of force balance for the three spheres are

[TABLE]

Notice that equations (15)-(17) imply that the condition of global force balance (4) is satisfied. Our aim now is to obtain two ODEs for $L_{1}$ and $L_{2}$. From equations (15)-(17) we can easily express the forces as functions of $T_{1}$ and $T_{2}$. If we plug the resulting expressions into equations (1)-(3) we obtain

[TABLE]

These equations, combined with the kinematic relations (5) and (6), yield the following system of ODEs governing the evolution of $L_{1}$ and $L_{2}$

[TABLE]

1.3 TSS with muscle-like arms

Equations (21) and (22) drive the evolution of $L_{1}$ and $L_{2}$ once $T_{1}$ and $T_{2}$ are known functions of shape and time. In this section we introduce a further element: we model the two arms of the swimmer according to Hill’s active state muscle model. A schematicrepresentation of this model is shown in Figure 2. The contractile element is composed of an active component capable of generating a tension $S(t)$ and a linear dashpot with characterisic constant $B$. In addition there are two linear springs, one in parallel and one in series with the contractile element, with elastic constants $k_{p}$ and $k_{s}$ respectively.

We assume that the active components of the two arms can generate tensions $S_{1}(t)$ and $S_{2}(t)$. Let us indicate with $l_{i}$ the sum of the rest lengths of the parallel and series springs in arm $i$ ($i=1,2$). Tensions and lengths are related by the following equations

[TABLE]

By rearranging the terms we find

[TABLE]

By using equations (21) and (22) we can express $\dot{L}_{1}$ and $\dot{L}_{2}$ as functions of $T_{1}$, $T_{2}$, $L_{1}$, and $L_{2}$. Then we can plug the resulting expressions into (10), (25), and (26). As a consequence, we can rewrite (21)-(26) and (10) as a system of five ODEs governing the evolution of $x,L_{1},L_{2},T_{1}$, and $T_{2}$.

1.4 Non-dimensionalization

The aim of this subsection is to write the problem in non-dimensional form. We will discover that the physical parameters influence the behaviour of the dynamics through $6$ non-dimensional numbers. This approach allows to obtain more general and useful results. Let $l=l_{1}+l_{2}$ be the characteristic length of the swimmer and $1/\omega$ a characteristic time. Let us set $t^{*}=\omega t$, $\Lambda_{i}=L_{i}/l$, $\lambda_{i}=l_{i}/l$, $\xi=x/l$, $\alpha=a/l$, $\tau_{i}=T_{i}/(\mu\omega l^{2})$, $\sigma_{i}=S_{i}/(\mu\omega l^{2})$. We will indicate with the prime symbol the derivative with respect to the rescaled time $t^{*}$. Using the new variables our system can be rewritten as

[TABLE]

where the dimensionless parameters are

[TABLE]

for $i=1,2$.

1.5 TSS with one muscle-like arm and one passive elastic arm

Now we replace the tail arm with a passive elastic spring. Let us indicate with $l_{1}$ the rest length and with $h$ the elastic constant of the spring. Equations (27)-(29) and (31) remain valid. Equation (30) is replaced by

[TABLE]

where $\lambda_{1}=l_{1}/l$ and $H=h/(\omega\mu l)$.

2 TSS with two muscle-like arms

In this section we present a detailed study of the TSS with two muscle-like arms. For the sake of simplicity we assume that the two arms are identical. Therefore

[TABLE]

We consider the case in which the tensions developed by the active components are small, namely,

[TABLE]

with $\epsilon<<1$ and $|\tilde{\sigma}_{i}|\leq 1$. We assume that $\sigma_{1}$ and $\sigma_{2}$ are periodic. In the first subsection we study the qualitative properties of the solutions. In the second subsection we compute the solutions analytically to leading order in $\epsilon$. In the third subsection we present the results of numerical simulations and in the fourth one we study some optimization problems.

2.1 Qualitative properties of the solutions

In this subsection we prove that, for small enough values of $\epsilon$, there exists one and only one periodic orbit and that this orbit is asymptotically stable. Our analysis is based on some classical results on periodically perturbed systems which can be found in [17].

We observe that in the ODEs (27)-(31), $\xi$ does not appear in the velocity field. This is due to the fact that the response of the system is invariant under translations. As a consequence, the differential problem for $\xi$ is decoupled from the rest. So we can restrict our attention to the four-dimensional system for $\Lambda_{1},\Lambda_{2},\tau_{1},\tau_{2}$. Instead of $\Lambda_{1}$ and $\Lambda_{2}$ we use the variables $u_{1}=\Lambda_{1}-\lambda_{1}$ and $u_{2}=\Lambda_{2}-\lambda_{2}$. Let us set $Y:=(u_{1},u_{2},\tau_{1},\tau_{2})^{tr}$ and $\tilde{J}:=J_{s}+J_{p}$. We can write our system in the form

[TABLE]

where $f(Y)$ is

[TABLE]

and

[TABLE]

The unperturbed system $Y^{\prime}=f(Y)$ has the equilibrium point $Y=0$. The variational system with respect to this equilibrium point (see [17], page 303) is

[TABLE]

where

[TABLE]

and

[TABLE]

Now we would like to prove that all the eigenvalues of $A$ have a negative real part. If this is true then all the characteristic multipliers of system (42) are in modulus strictly less than one. As a consequence we can apply theorems 6.1.1 and 6.1.3 of [17] and conclude that, for small enough values of $\epsilon$, there exists one and only one periodic orbit. Moreover, this orbit is asymptotically stable. The characteristic polynomial of $A$ is

[TABLE]

Let us indicate with $a_{j}$, $j=0,1,2,3$ the coefficients, so that

[TABLE]

First of all we observe that $a_{j}>0$ for every $j$. This follows from the fact that, if $\alpha$ is small enough, then $Q<0$, $R>0$, $|Q|<R$. Having all the coefficients with the same sign is a necessary condition for $p_{A}(x)$ to have only roots with negative real part. The assumption that $\alpha$ is small is a natural one, because all of our analysis is based on the hypothesis that the radius of the spheres is small compared to the lengths of the arms. The table associated to the polynomial by means of the Routh method is

[TABLE]

where

[TABLE]

Let us study the signs of $b_{3}$ and $c_{2}$. Notice that

[TABLE]

Now we notice that, since $|Q|<R$, we have

[TABLE]

As a consequence,

[TABLE]

Notice that $\text{sign}(c_{2})=\text{sign}(a_{1}b_{3}-a_{3}b_{2})$. By using (56) we see that

[TABLE]

So $b_{3}$ and $c_{2}$ are both positive. Therefore we see from table (48) that there are three roots with negative real part. Since the determinant of $A$ is positive, also the fourth root must have negative real part.

2.2 Asymptotic expansions

In this section we compute analytic expressions for the solutions, to leading order in $\epsilon$. Let us assume that

[TABLE]

with $|\tilde{a}_{j}|+|\tilde{b}_{j}|\leq 1$, $j=1,2$. We express the solutions in power series as follows

[TABLE]

Let us set $X=(u_{1}^{(1)},u_{2}^{(1)},\tau_{1}^{(1)},\tau_{2}^{(1)})^{tr}$. The problem satisfied by $X$ is

[TABLE]

The solution of this problem can be computed with the help of Duhamel’s formula

[TABLE]

However, since the exponential of $A$ is difficult to compute, we choose a different strategy. We know from the previous subsection that the solution converges to a periodic orbit and we would like to study the stationary regime of the system. So we are not really interested in the solution of (64). Rather, we would like to compute the periodic orbit. This problem can be reduced to the problem of solving a linear system. First of all we notice that $g$ can be rewritten as

[TABLE]

where $\hat{g}_{s}=(0,0,J_{s}\tilde{a}_{1},J_{s}\tilde{a}_{2})^{tr}$ and $\hat{g}_{c}=(0,0,J_{s}\tilde{b}_{1},J_{s}\tilde{b}_{2})^{tr}$. We look for a solution in the following form

[TABLE]

with $\hat{X}_{s},\hat{X}_{c}\in\mathbb{R}^{4}$. Notice that $X^{\prime}=-\hat{X}_{c}\sin(t^{*})+\hat{X}_{s}\cos(t^{*})$. So our differential equation becomes

[TABLE]

Since the equality must hold for every time $t^{*}$, we obtain the following linear system

[TABLE]

We can rewrite the problem in compact form

[TABLE]

Now suppose that

[TABLE]

namely, $\hat{X}_{s}=(a_{1},a_{2},c_{1},c_{2})^{tr}$ and $\hat{X}_{c}=(b_{1},b_{2},d_{1},d_{2})^{tr}$. The linear problem (70) corresponds to the following set of equations

[TABLE]

The first four equations allow us to express $a_{1}$, $a_{2}$, $b_{1}$, and $b_{2}$ as linear combinations of $c_{1}$, $c_{2}$, $d_{1}$, and $d_{2}$. If we plug the resulting expressions into the last four equations we obtain

[TABLE]

where

[TABLE]

So now our problem is to check if $M$ is invertible and to compute $M^{-1}$. Let us set

[TABLE]

Straightforward computations show that

[TABLE]

We observe that when $\alpha$ is sufficiently small $|Q|>R$ and so $M_{3}>M_{4}$. It follows that $(M_{3}^{2}-M_{4}^{2})^{2}>0$, thus $\text{det}(M)$ is strictly positive. So we can compute the inverse of $M$, and the result is

[TABLE]

where

[TABLE]

Now we can compute the coefficients $a_{i}$, $b_{i}$, $c_{i}$, $d_{i}$, $i=1,2$. First we observe that

[TABLE]

The remaining coefficients can be computed using (75)-(78). At steady state the shape of the system evolves along a closed loop in the configuration space. Figure 3 shows a comparison between the loops obtained through numerical simulations and the ones corresponding to the leading order approximation.

2.3 Numerical simulations

In this section we present the results of some numerical simulations. We used MATLAB ode45 procedure, which consists of a Runge-Kutta integration scheme with adaptive step size. The values of the dimensionless parameters used in the simulations are shown in Table 1.

In section 2.1 we proved that the system converges to a periodic orbit. This is confirmed by Figure 4 in which we see that, after relaxation, $\Lambda_{1}$ and $\Lambda_{2}$ evolve along a closed loop. Figure 5 shows the evolution of $\xi$. In this simulation $(\tilde{a}_{1},\tilde{b}_{1})$ is proportional to $(\tilde{a}_{2},\tilde{b}_{2})$ (see Table 1). This means that the tensions developed in the active components of the two arms are in phase with each other. If this sinchronization of the active components caused a sinchronization in the length change of the two arms, we would obtain a reciprocal shape change, thus no net motion in view of the scallop theorem. As we can see from figures 4 and 5, this is not the case. This property is very interesting because it implies that the swimmer can move even if the two active components are stimulated at the same frequency. We will go back to this in the following section, where this property of the system will emerge from the formula for the leading order term of the net displacement in one period. In Figure 6 we see the evolution of $\tau_{1}$ and $\tau_{2}$: also in this case the behaviour at steady state is periodic, in agreement with the results of subsection 2.1.

2.4 Optimization

In this subsection we study some optimization problems. We consider three performance measures: net displacement in one period, work per travelled distance, and Lighthill’s efficiency. The value of the performance measures is determined by the dimensionless parameters that appear in our system of ODEs. These parameters depend on the swimmer’s body length $l$, on the actuation frequency $\omega$, on the fluid viscosity $\mu$, and on the physical constants $k_{s},k_{p},h,B$ characterizing the viscoelastic arms. We consider the case in which the fluid and the swimmer are given, so that the only parameter which can be varied is the actuation frequency $\omega$. We would like to study optimization of the different performance measures with respect to $\omega$. Our strategy is to compute the leading order approximation of the performance measures and use the resulting expressions to study our optimization problem. Then we will compare the optimality results obtained through the leading order approximation with the outcome of numerical simulations.

The first performance measure we consider is $\Delta x$, the net displacement per period in units of body length. This quantity can be computed using formula (14). Through a simple change of integration variable we see that

[TABLE]

By using the leading order expressions computed above we find the following leading order approximation for $\Delta x$

[TABLE]

By using equations (75)-(78) we obtain

[TABLE]

Now by using the expressions for $c_{1},c_{2},d_{1}$ and $d_{2}$ computed in subsection 2.2 we obtain

[TABLE]

There are some interesting facts to point out about this formula. First of all we notice that if we change $\sigma_{1}$ with $\sigma_{2}$ and viceversa the displacement changes sign, as expected for symmetry reasons. Secondly, we observe that even in the case in which only one of the two arms is activated we obtain a non-zero displacement. Finally, let us comment on the case in which the two active components are sinchronized. Suppose there exists a constant $\eta>0$ such that

[TABLE]

In this case one might expect a sinchronization of the length change of the two arms, leading to a reciprocal shape change which produces no net motion. However, our asymptotic analysis shows that is not the case provided that $\eta\neq 1$: from formula (99) we obtain that the leading order term of $\Delta x$ is

[TABLE]

Now we compute the mechanical work $\mu\omega l^{3}W$ done by the active components in one period. The power expenditure $\mu\omega^{2}l^{3}\mathcal{P}$ is

[TABLE]

It follows that

[TABLE]

Using the leading order approximations computed in subsection 2.2 we find that

[TABLE]

The other two performance measures we consider are are the mechanical work per travelled distance

[TABLE]

and Lighthill’s efficiency

[TABLE]

Since we have already computed leading order expressions for $\Delta x$ and $W$, we can easily compute leading order expression for $\zeta$ and $\eta$.

Now we would like to optimize the leading order terms of the different performance measures with respect to the actuation frequency $\omega$. First of all we will show that this problem is well-posed, namely, that each performance measure admits an optimal value of $\omega$. Secondly, we will show that the optimality results obtained by studying the leading order approximations are in good agreement with the results of numerical simulations.

The coefficients $a_{i},b_{i},c_{i},$ and $d_{i}$, $i=1,2$, have been computed in subsection 2.2. It is not difficult to check that these coefficients are $O(1/\omega)$ when $\omega\rightarrow+\infty$ and $O(\omega)$ when $\omega\rightarrow 0$. As a consequence, it is easy to verify that

[TABLE]

So the leading order expression for $\Delta x$ vanishes when $\omega$ goes to zero and to $+\infty$. This implies that there exists $\omega_{\Delta x}\in(0,+\infty)$ which optimizes this performance measure. The situation is similar for the other performance measures:

[TABLE]

From these asymptotic regimes we deduce the existence of two optimal frequencies $\omega_{\zeta},\omega_{\eta}\in(0,+\infty)$. Figures 7, 8, and 9 show plots of the different performance measures as functions of $\omega$. Notice the existence of the optimal frequencies and the very good agreement between leading order approximations and numerical simulations.

3 TSS with one muscle-like arm and one passive elastic arm

In this section we study the TSS with one muscle-like arm and one passive elastic arm. We assume that the active component of the muscle-like arm generates a tension

[TABLE]

with $\epsilon<<1$ and $|\tilde{\sigma}|\leq 1$. We assume that $\sigma$ is a periodic functionof time. In the first subsection we study the qualitative properties of the solutions. In the second subsection we compute the the leading order approximation of the solutions in the asymptotic regime $\epsilon<<1$. In the third subsection we present the results of numerical simulations. In the fourth subsection we study some optimization problems.

3.1 Qualitative properties of the solutions

For the three-sphere swimmer with one muscle-like arm and one passive elastic arm $\tau_{1}$ is given by equation (33). So we can restrict our attention to the three-dimensional system of ODEs governing the evolution of the two shape variables and $\tau_{2}$. The state of the system is $Y=(u_{1},u_{2},\tau_{2})^{tr}$ and the problem has the form

[TABLE]

where $f(Y)$ is given by

[TABLE]

and

[TABLE]

The unperturbed system $Y^{\prime}=f(Y)$ has the equilibrium point $Y=0$. The variational system with respect to this equilibrium point (see [17], page 303) is

[TABLE]

with

[TABLE]

where

[TABLE]

The characteristic polynomial of $B$ is

[TABLE]

If we prove that all the eigenvalues of $B$ have negative real part then we can argue as in section 2.1 and conclude that, for $\epsilon$ small enough, there exists one and only one periodic orbit which is asymptotically stable. Let us call $a_{j}$, $j=0,1,2$ the coefficients of the characteristic polynomial, so that

[TABLE]

If $\alpha$ is small enough then $P<0$, $Q<0$, $R>0$, $|P|<R$, and $|Q|<R$. So under this hypothesis all the coefficients of the polynomial are positive. The table associated to the polynomial by means of the Routh method is

[TABLE]

with

[TABLE]

Notice that

[TABLE]

It follows that $b_{2}>0$. So by applying Routh’s criterion we see from table (118) that there are two roots with negative real part. Since the determinant of $B$ is negative, also the third root must have negative real part.

3.2 Asymptotic expansions

Now we would like to compute analytically the leading order term of the solutions. Let us express the solutions in power series as follows

[TABLE]

Let us set $X=(u_{1}^{(1)},u_{2}^{(1)},\tau_{2}^{(1)})^{tr}$. $X(t^{*})$ is the solution of

[TABLE]

The solution of this problem can be computed with the help of Duhamel’s formula. Since we are interested in the periodic behaviour of the solution at steady state, we would like to compute the periodic orbit of the system. We assume that

[TABLE]

We can write $g$ as follows

[TABLE]

where $\hat{g}_{s}=(0,0,J_{s}\tilde{a})^{tr}$ and $\hat{g}_{c}=(0,0,J_{s}\tilde{b})^{tr}$. We look for the periodic solution in the following form

[TABLE]

with $\hat{X}_{s},\hat{X}_{c}\in\mathbb{R}^{3}$. So our differential equation becomes

[TABLE]

Since the equality must hold for every time $t^{*}$, we obtain the following linear system

[TABLE]

We can rewrite the system in compact form

[TABLE]

Now suppose that

[TABLE]

namely, $\hat{X}_{s}=(a_{1},a_{2},c)^{tr}$ and $\hat{X}_{c}=(b_{1},b_{2},d)^{tr}$. Let us rewrite (130) explicitly

[TABLE]

From the first four equations we obtain

[TABLE]

If we plug these expressions into the last two equations we obtain the following problem for $c$ and $d$

[TABLE]

where

[TABLE]

and

[TABLE]

The linear problem (130) has a solution if and only if (185) has a solution. So we need to check if $N$ is invertible. We notice that

[TABLE]

hence $M$ is invertible. The inverse of $M$ is given by

[TABLE]

Now $c$ and $d$ can be computed as follows

[TABLE]

Once we know $c$ and $d$ we can compute $a_{1},a_{2},b_{1},b_{2}$ by using equations (140)-(143). At steady state the two shape variables evolve along a closed loop in the configuration space. Figure 10 shows a comparison between the closed curves obtained through numerical simulations and the ones corresponding to the leading order approximation.

3.3 Numerical simulations

In this section we present the results of some numerical simulations, obtained by choosing the dimensionless parameters as in Table 2.

Figure 11 shows the evolution of the system in the shape space. We observe that, after relaxation, $\Lambda_{1}$ and $\Lambda_{2}$ evolve along a closed loop. In Figure 12 we see the evolution of $\xi$. We notice that the net displacement in one period has a negative sign, which means that the object swims with the passive arm ahead. This is a general property of the TSS with passive elastic tail. Given any periodic deformation of the active arm, it can be shown that the corresponding net displacement in one period has a negative sign [16]. Figure 13 shows a plot of $\tau_{1}$ and $\tau_{2}$ as functions of time. The behaviour at steady state is periodic, as expected from the results of subsection 3.1.

3.4 Optimization

In this subsection we study some optimization problems for the three-sphere swimmer with one muscle-like arm and one passive elastic arm. The optimality measures we consider are the net displacement in one period $\Delta x$, the mechanical work per travelled distance $\zeta$ and Lighthill’s efficiency $\eta$. Like in section 2.4, we assume that the fluid and swimmer are given, so that the only parameter which can be varied is the actuation frequency $\omega$. First we study our optimization problem using the leading order approximation of the different performance measures and then we compare the results with the outcome of numerical simulations.

We start by computing the leading order expressions for the different performance measures. Formula (97) for the leading order approximation of $\Delta x$ remains valid. Let us consider the mechanical work $\mu\omega l^{3}W_{2}$ done by the active component of the muscle-like arm in one period. The power expenditure is

[TABLE]

It follows that

[TABLE]

So by using the results of section 3.2 we obtain

[TABLE]

The work per travelled distance $\eta$ and Lighthill’s efficiency $\eta$ are defined as in section 2.4, and the corresponding leading order approximations can be easily obtained once we know the approximation of $\Delta x$ and $W$.

Using the formulas computed in section 3.2 it is not difficult to see that the coefficients $c,d,a_{i},b_{i}$, $i=1,2$ are $O(1/\omega)$ when $\omega\rightarrow+\infty$ and $O(\omega)$ when $\omega\rightarrow 0$. Therefore, the asymptotic behaviour of the leading order terms of the different performance measures is exactly the same as the one observed in section 2.4. This implies existence of three optimal frequencies $\omega_{\Delta x}$, $\omega_{\zeta}$, and $\omega_{\eta}$. These theoretical predictions are confirmed by Figures 14, 15, and 16. In these figures we observe the existence of the three optimal frequencies $\omega_{\Delta x}$, $\omega_{\zeta}$, $\omega_{\eta}$ and the high accuracy of the leading order approximations.

Conclusions

We studied the dynamics of the three-sphere swimmer with muscle-like arms. We assumed that the forces generated by the swimmer in the active components of the arms have intensity $\epsilon$ and vary periodically with frequency $\omega$. We showed that the two shape parameters and the forces acting across the arms evolve according to a system of ODEs. We proved that the solutions converge to a periodic orbit. Under the assumption that $\epsilon<<1$, we computed the leading order approximation of the solutions at steady state. Then we studied some optimization problems. We considered three different performance measures: net displacement in one period, work per travelled distance, and Lighthill’s efficiency. We studied optimization of these quantities with respect to the actuation frequency $\omega$. We computed leading order approximations of the different performance measures. Using these approximations we showed that each performance measure admits an optimal frequency. In addition we showed that the optimality results obtained through leading order approximations are in very good agreement with the outcome of numerical simulations. Then we introduced the three-sphere swimmer with one muscle-like arm and one passive elastic arm. We studied this model through the same type of analysis done for the three-sphere swimmer with two muscle-like arms.

Appendix A Asymptotic expansions when the forces in the active components are generic periodic functions

A.1 TSS with two muscle-like arms

In section 2.2 we computed the leading order approximation of the solution in the case in which $\tilde{\sigma}_{1}$ and $\tilde{\sigma}_{2}$ are given by equations (58) and (59). Now we would like to consider the case in which $\tilde{\sigma}_{1}$ and $\tilde{\sigma}_{2}$ are two generic period functions with period $2\pi$. We consider the expansion in Fourier series

[TABLE]

Le $X=(u_{1}^{(1)},u_{2}^{(1)},\tau_{1}^{(1)},\tau_{2}^{(1)})^{tr}$ be the leading order term of the solution. $X$ satisfies the differential equation

[TABLE]

where $A$ and $g$ are defined as in (43) and (41) respectively. Notice that $g$ can be written as

[TABLE]

where $\hat{g}_{s,j}=(0,0,J_{s}\tilde{a}_{1,j},J_{s}\tilde{a}_{2,j})^{tr}$ and $\hat{g}_{c,j}=(0,0,J_{s}\tilde{b}_{1,j},J_{s}\tilde{b}_{2,j})^{tr}$. We are interested in computing the periodic orbit of the system, so we look for a solution in the following form

[TABLE]

with $\hat{X}_{s,j},\hat{X}_{c,j}\in\mathbb{R}^{4}$ for every $j$. It is not difficult to see that the ODE (156) is equivalent to the set of linear systems

[TABLE]

for $j\in\mathbb{N}$. Now suppose that

[TABLE]

and

[TABLE]

We may rewrite (159) as follows

[TABLE]

Equations (162)-(165) allow us to express $a_{1,j}$, $a_{2,j}$, $b_{1,j}$, and $b_{2,j}$ as linear combinations of $c_{1,j}$, $c_{2}{j}$, $d_{1,j}$, and $d_{2,j}$. If we plug the resulting expressions into the last four equations we obtain the linear problems

[TABLE]

where

[TABLE]

It is not difficult to check that $M(j)$ is invertible for every $j\in\mathbb{N}$: the computations are similar to the ones done in section 2.2 for the matrix $M$. So we can conclude our computation with the help of the matrices $M(j)^{-1}$, $j\in\mathbb{N}$.

A.2 TSS with one muscle-like arm and one passive elastic arm

In this subsection we would like to generalize the results of subsection 3.2 by considering a more general periodic function $\tilde{\sigma}$. We assume that $\tilde{\sigma}$ is a periodic function with period $2\pi$. We expand this function in Fourier series

[TABLE]

Let $X=(u_{1}^{(1)},u_{2}^{(2)},\tau_{2}^{(2)})^{tr}$ be the leading order term of the solution. $X$ satisfies the ODE

[TABLE]

where $B$ and $g$ are defined by equations (112) and (110) respectively. Notice that $g$ can be written as

[TABLE]

where $\hat{g}_{s,j}=(0,0,J_{s}\tilde{a}_{j})^{tr}$ and $\hat{g}_{c,j}=(0,0,J_{s}\tilde{b}_{j})^{tr}$. We would like to compute the periodic orbit of the system, so we look for a solution in the form

[TABLE]

with $\hat{X}_{s,j},\hat{X}_{c,j}\in\mathbb{R}^{3}$ for every $j$. The ODE (173) is equivalent to the set of linear systems

[TABLE]

for $j\in\mathbb{N}$. Now suppose that

[TABLE]

and

[TABLE]

Let us rewrite (176) explicitly

[TABLE]

Using the first four equations we can express $a_{1,j}$, $a_{2,j}$, $b_{1,j}$, and $b_{2,j}$ as linear combinations of $c_{j}$ and $d_{j}$. Then we can plug the resulting expressions into the last to equations. The result is the following linear problem

[TABLE]

where

[TABLE]

and

[TABLE]

A simple computation shows that $\det(N(j))>0$ for every $j$. So with the help of the matrices $N(j)^{-1}$ we can compute the coefficients of our leading order approximation.

Acknowledgments

Support by the European Research Council through the ERC Advanced Grant 340685-MicroMotility is gratefully acknowledged.