A limiting model for a low Reynolds number swimmer with N passive elastic arms

TL;DR

This paper analyzes a low Reynolds number swimmer with an active arm and multiple passive elastic arms, proving the convergence of a discrete model to a continuous limit and exploring its swimming performance through numerical experiments.

Contribution

It introduces a generalized model with N passive springs, proves its convergence to a continuous model as N increases, and evaluates its swimming efficiency numerically.

Findings

Convergence of discrete to continuous model established

Swimmer's displacement depends on oscillation frequency and amplitude

Numerical results demonstrate performance variations with parameters

Abstract

We consider a low Reynolds number artificial swimmer that consists of an active arm followed by $N$ passive springs separated by spheres. This setup generalizes an approach proposed in Montino and DeSimone, Eur. Phys. J. E, vol. 38, 2015. We further study the limit as the number of springs tends to infinity and the parameters are scaled conveniently, and provide a rigorous proof of the convergence of the discrete model to the continuous one. Several numerical experiments show the performances of the displacement in terms of the frequency or the amplitude of the oscillation of the active arm.

Tables

Table 1. Table 1: Values of the parameters used in the numerical simulations, matching those of [ 14 ] . We have taken for μ 𝜇 \mu the dynamic viscosity of water at 25 ∘ ​ C superscript 25 𝐶 25^{\circ}C .

$\tilde{a}$	$1\cdot {10}^{-5}m$
$a_1$	$1\cdot {10}^{-5}m$
$\Lambda$	$4\cdot {10}^{-4}m$
$L$	$3\cdot {10}^{-5}m$
$\tilde{k}$	$1\cdot {10}^{-8}N{m}^{-1}$
$\mu$	$8.9\cdot {10}^{-4}Pas$

Equations

$$\left\{\begin{array}[]{ll}f_{j}^{F}=-6\pi\mu aV_{j}\text{ for }j\geq 3,\vspace{6pt}\\ f_{j}^{F}=-6\pi\mu a_{1}V_{j}\text{ for }j=1,2.\end{array}\right.$$

$$\left\{\begin{array}[]{ll}f_{2}^{R}&\displaystyle=k(L_{1}-h)=k\frac{\ell_{1}}{N}\vspace{6pt}\\ f_{j}^{R}&=k\big{(}(L_{j-1}-h)-(L_{j-2}-h)\big{)},\vspace{6pt}\\ &\displaystyle=k\frac{\ell_{j-1}-\ell_{j-2}}{N}\quad\text{ for }3\leq j\leq N+1\vspace{6pt}\\ f_{N+2}^{R}&\displaystyle=-k(L_{N}-h)=-k\frac{\ell_{N}}{N}.\end{array}\right.$$

$$\left\{\begin{array}[]{ll}f_{1}^{F}+\dots+f_{N+2}^{F}=0,\vspace{6pt}\\ f_{j}^{R}+f_{j}^{F}=0\text{ for }j\geq 3.\end{array}\right.$$

$$\left\{\begin{array}[]{ll}f_{1}^{F}-f_{2}^{F}=6\pi\mu a_{1}(V_{2}-V_{1})=6\pi\mu a_{1}\dot{L}_{0},\vspace{6pt}\\ f_{1}^{F}+f_{2}^{F}=f_{3}^{R}+\dots+f_{N+2}^{R}=-k\ell_{1}/N,\end{array}\right.$$

$$\left\{\begin{array}[]{ll}f_{1}^{F}=\frac{1}{2}(+6\pi\mu a_{1}\dot{L}_{0}-\tilde{k}\ell_{1}),\vspace{6pt}\\ f_{2}^{F}=\frac{1}{2}(-6\pi\mu a_{1}\dot{L}_{0}-\tilde{k}\ell_{1}).\end{array}\right.$$

$${\dot{\ell }}_j={\Lambda }^{2}K\frac{{\ell }_{j-1}-2{\ell }_j+{\ell }_{j+1}}{h^2},2\le j\le N,$$

$${\ell }_{N+1}=0,$$

$$h{\dot{\ell }}_1={\Lambda }^{2}K\frac{{\ell }_2-{\ell }_1}{h}-\frac{\Lambda K\tilde{a}}{2a_1}{\ell }_1-\frac{\Lambda }{2}{\dot{L}}_0.$$

$${\ell }_j(t)=({\alpha }_d{\gamma }_{+}^{j-1}+{\beta }_d{\gamma }_{-}^{j-1})e^{i\omega t},$$

$${\gamma }_{\pm }=\frac{i/(K_{\omega }{N}^2)+2\pm \sqrt{\Delta }}{2}$$

$$\left\{\begin{array}[]{ll}\displaystyle\ell_{1}(t)=b_{d}e^{i\omega t}=e^{i\omega t}(\alpha_{d}+\beta_{d}),\vspace{6pt}\\ \displaystyle\ell_{N+1}(t)=e^{i\omega t}(\alpha_{d}\gamma_{+}^{N}+\beta_{d}\gamma_{-}^{N})=0\,,\end{array}\right.$$

$${\alpha }_d=\frac{-{\gamma }_{-}^N{b}_d}{({\gamma }_{+}^N-{\gamma }_{-}^N)},{\beta }_d=\frac{{\gamma }_{+}^N{b}_d}{({\gamma }_{+}^N-{\gamma }_{-}^N)}.$$

$$b_d=-\frac{\epsilon i/2}{i/N+N{K}_{\omega }(1-z_d)+K_{\omega }\frac{\tilde{a}}{2a_1}},$$

$${\partial }_t\ell (y,t)=K{\Lambda }^2{\partial }_{yy}\ell (y,t),\forall (y,t)\in [0,\Lambda ]\times {\mathbb{R}}_{+}^{\star }.$$

$${\Lambda }^{2}K{\partial }_y\ell (0,t)-\Lambda K\frac{\tilde{a}}{2a_1}\ell (0,t)=\frac{\Lambda }{2}{\dot{L}}_0(t),\forall t>0.$$

$$\left\{\begin{array}[]{l}\displaystyle\partial_{t}\ell(y,t)-\Lambda^{2}K\partial_{yy}\ell(y,t)=0,\vspace{6pt}\\ \displaystyle\Lambda^{2}K\partial_{y}\ell(0,t)-\Lambda K\frac{\tilde{a}}{2a_{1}}\ell(0,t)=\frac{\Lambda}{2}\dot{L}_{0}(t),\vspace{6pt}\\ \displaystyle\ell(\Lambda,t)=0.\end{array}\right.$$

$$\mathcal{V}=\left\{u\in H^1((0,\Lambda ))|u(\Lambda )=0\right\},$$

$$\frac{d}{dt}{\int }_0^{\Lambda }\ell vdy+{\Lambda }^{2}K{\int }_0^{\Lambda }{\partial }_y\ell {\partial }_{y}vdy$$

$$+\frac{\Lambda K\tilde{a}}{2a_1}\ell (0,t)v(0)=-\frac{\Lambda }{2}{\dot{L}}_0(t)v(0)$$

$$\kappa :(u,v)\mapsto {\Lambda }^{2}K{\int }_0^{\Lambda }{\partial }_{y}u(y){\partial }_{y}v(y)dy+\frac{\Lambda K\tilde{a}}{2a_1}u(0)v(0),$$

$$i\underline{\ell }={\Lambda }^2{K}_{\omega }{\partial }_{yy}\underline{\ell }.$$

$$\underline{\ell }(y)=\alpha e^{ry}+\beta e^{-ry},$$

$$\left\{\begin{array}[]{ll}\displaystyle-(\alpha+\beta)\frac{\tilde{a}}{2a_{1}}+\Lambda r(\alpha-\beta)=\frac{i\varepsilon}{2K_{\omega}},\vspace{6pt}\\ \displaystyle\alpha e^{r\Lambda}+\beta e^{-r\Lambda}=0,\end{array}\right.$$

$$\left\{\begin{array}[]{ll}\displaystyle\alpha=\frac{i\varepsilon}{\displaystyle 2K_{\omega}\big{(}\frac{\tilde{a}}{2a_{1}}(e^{2r\Lambda}-1)+\Lambda r(e^{2r\Lambda}+1)\big{)}},\vspace{6pt}\\ \displaystyle\beta=-e^{2r\Lambda}\alpha.\end{array}\right.$$

$$\frac{d(M_h{L}_h)}{dt}+K_h{L}_h=\tilde{f}(t),$$

$$(K_{h})_{ij}=\left\{\begin{array}[]{ll}-\Lambda^{2}K/h&\mbox{for }|i-j|=1,\\ 2\Lambda^{2}K/h&\mbox{for }i=j\geq 2\,,\\ \Lambda^{2}K/h+\Lambda K\tilde{a}/(2a_{1})&\mbox{for }i=j=1\,,\end{array}\right.$$

$$(M_{h})_{ij}=\left\{\begin{array}[]{ll}h/6&\mbox{for }|i-j|=1,\\ 2h/3&\mbox{for }i=j\geq 2\,,\\ h/3&\mbox{for }i=j=1\,.\end{array}\right.$$

$${\tilde{M}}_h=\left(\begin{array}{ccc}h& & 0\\ & \ddots & \\ 0& & h\end{array}\right).$$

$$\frac{d({\tilde{M}}_h{L}_h)}{dt}+K_h{L}_h=\tilde{f}(t),$$

$${\ell }_h(0)={\ell }_{0,h}\in {\mathcal{V}}_h.$$

$${\bar{M}}_h=\left(\begin{array}{cccc}h/2& & & 0\\ & h& & \\ & & \ddots & \\ 0& & & h\end{array}\right)$$

Full text

1 Introduction

As stated by Purcell’s Scallop Theorem [3], reciprocal shape changes in a swimmer never leads to a net displacement of the system in a low Reynolds number setting. Indeed, a microscopic scallop opening and closing its valve would be completely unable to swim, due to negligible inertial forces in this situation [4]. Several simple mechanisms have then been introduced (see e.g. [5]) to overcome this obstruction, most of them using two degrees of freedom in order to create closed curves with nonzero surface in the shape space of the swimmer.

One of the simplest mechanisms introduced in the literature is probably Najafi and Golestanian’s three-sphere swimmer [6], which consists in three spheres linked by two extensible arms of negligible thickness, moving along a single direction. This model is much simpler than Purcell’s original three-link swimmer [3], or Purcell’s rotator [7], as there is no rotational motion involved. This swimmer has two degrees of freedom, activated periodically in time with a phase lag in order to produce the loop. Both Purcell’s and Najafi and Golestanian’s swimmers have been extensively studied in [8, 9, 10, 11, 12, 13].

As an extension of this three-sphere swimmer, Montino and DeSimone then introduced a three-sphere swimmer with a passive elastic arm [14]. This adaptation has only one degree of freedom, which is the length of the non-elastic arm. Thanks to a resonant effect at natural frequency of the system (which depends on the viscosity of the fluid, the masses and the spring constant), an out-of-phase oscillation of the spring is created, which ultimately leads to a net motion of the swimmer. However, at very low or very high frequency, no net motion is possible after a stroke. Having this passive elastic arm also confines net motion to only one direction on the swimming axis, swimming direction is thus limited, and the swimmer can only move with its passive arm ahead. This was also denoted by Passov [15], when looking at Purcell’s three-link swimmer with a passive elastic tail.

In this paper, Montino and DeSimone’s swimmer is extended by adding a large number $N$ of passive elastic arms to their one-dimensional swimmer, thus turning it into an $(N+2)$-sphere swimmer. This simple swimmer then leads to a limit model with an elastic tail resembling a one-dimensional flagella along which compressive waves propagate.

The paper is organized as follows. In sect. 2, we describe the $N$-spring swimmer, and its equations of motion, before looking at the limit model, when the number of springs tends to infinity, in sect. 2.3. We prove the convergence of the discrete model to the continuous one in sect. 3, using the fact that it is found to be a non-conventional mass lumping discretization of the limit model. Sect. 4 introduces two formulas in order to compute the net displacement of both swimmers, discrete and continuous. Finally, in sect. 5 we study numerically the movement and displacement of our swimmer depending on various system parameters, in order to find optimal swimming parameters to obtain the largest net displacement possible.

2 Problem’s formulation and study: $N$-spring discrete model and its continuous limit

The swimmer studied in this paper is an extension of the three-sphere swimmer with a single passive elastic arm [14], to a swimmer with $N+2$ spheres and $N$ passive elastic arms, presented in figure 1. The first arm of this artificial swimmer is a rod of negligible thickness, surrounded by two spheres of radius $a_{1}$. This arm has a prescribed periodic movement around a length at rest $L$, of the form $L_{0}(t)=L(1+\tilde{\varepsilon}\cos(\omega t))$ where $\tilde{\varepsilon}\in[0,1)$ is a non-dimensional parameter. $\tilde{\varepsilon}<1$ so that the active arm always has a positive length. We define $\varepsilon$ as $\varepsilon=L\tilde{\varepsilon}$. The rest of the swimmer has a total length at rest $\Lambda$ that does not depend on $N$. In order to keep a constant length and have an elastic force that does not depend on $N$, all the other spheres have a radius $a=\tilde{a}/N$, the springs have each a rest length $h=\Lambda/N\gg a$, and an elastic constant $k=\tilde{k}N$, with $\tilde{k}$ and $\tilde{a}$ prescribed and independent of $N$.

If the swimmer is able to control the length of the front rod with the prescribed periodic function $L_{0}(t)$, the length of the $N$ remaining springs are governed by the balance of viscous and elastic forces. At any time $t$, the length $L_{j}(t)$ of the $j$-th arm, $j\geq 1$ is written as $L_{j}(t)=\frac{\ell_{j}(t)}{N}+h$. Let us then denote by $\mu$ the fluid viscosity, $f_{j}^{F}$ and $f_{j}^{R}$ the hydrodynamic and elastic forces on the $j$-th sphere. We also call $x_{j}$ the coordinate of its center, so that $V_{j}=\dot{x}_{j}$ is the velocity of the $j$-th sphere. The geometry of the system entails $\dot{L}_{j}=V_{j+2}-V_{j+1}$ for all $j=0,\,\dots,\,N$.

In order to effectively swim, our $N$-spring swimmer undergoes periodic harmonic but non-reversible deformations, just like the original swimmers from Najafi and Golestanian [6], and Montino and DeSimone [14]. However, due to the geometry, we expect a wave to propagate along the tail. This is the behaviour of this wave that we aim at describing in the remainder of the paper.

2.1 First approximations

In a first approximation, we consider the case where the hydrodynamic force on the $j$-th sphere only depends on the speed of that same sphere, and neglect interactions between spheres. This leads to the following set of equations on (fluid) forces and velocities:

[TABLE]

The elastic forces on each sphere can be written as:

[TABLE]

At low Reynolds number, inertial forces are negligible. This, together with the fact that the artificial swimmer is self-propelled, gives:

[TABLE]

Using (1), (2) and (3), we obtain the expression of fluid forces on each sphere with respect to the length of the adjacent arms. In particular, for the first two spheres:

[TABLE]

which finally leads to:

[TABLE]

2.2 Movement of the spheres

In order to write the equations governing the system, we use equations (1-5) to find ODEs on the elongation $l_{j}(t)$ of the $j$-th arm, for $j\geq 1$. We first consider the case $j\geq 2$. Writing $\dot{L}_{j}=V_{j+2}-V_{j+1}=\frac{1}{6\pi\mu a}(f^{R}_{j+2}-f^{R}_{j+1})$, one deduces

[TABLE]

where we have added a fictitious variable

[TABLE]

and with $K=\displaystyle\frac{\tilde{k}}{6\pi\mu\tilde{a}}$ .

To determine the equation for the first elastic arm, we use the fact that $\displaystyle\dot{L}_{1}=V_{3}-V_{2}=-\frac{1}{6\pi\mu a}f^{F}_{3}+\frac{1}{6\pi\mu a_{1}}f^{F}_{2}$ to obtain, using equations (2) and (5):

[TABLE]

We can easily verify that the ODE problem (6,7,8) is well-posed using Cauchy-Lipschitz theorem, and provides a unique solution $(\ell_{j}(t))_{1\leq j\leq N+1}$ for any initial configuration.

Seeking for periodic (complex) solutions to equation (6) leads to

[TABLE]

where $\alpha_{d},\beta_{d}\in\mathbb{C}$ and

[TABLE]

and $\displaystyle\Delta=\frac{-1}{K_{\omega}^{2}N^{4}}+\frac{4i}{K_{\omega}N^{2}}$, where $\displaystyle K_{\omega}=\frac{K}{\omega}=\frac{\tilde{k}}{6\pi\mu\tilde{a}\omega}$ is an adimensional number. Notice that $|\gamma_{+}|>1$ while $|\gamma_{-}|<1$. The constants $\alpha_{d}$ and $\beta_{d}$ may be determined through the boundary conditions. Namely assuming, from the linearity of the problem, $\ell_{1}=b_{d}e^{i\omega t}$, with $b_{d}\in\mathbb{C}$ and recalling $l_{N+1}=0$ enables us to write

[TABLE]

to finally obtain

[TABLE]

Then, we use (8) to determine $b_{d}$:

[TABLE]

where $\displaystyle z_{d}=\frac{\gamma_{+}^{N}\gamma_{-}-\gamma_{-}^{N}\gamma_{+}}{\gamma_{+}^{N}-\gamma_{-}^{N}}$ .

2.3 Limit model with an infinite number of springs

As we increase the number of springs in our swimmer, a limit model arises, with an elastic-like tail, as shown in figure 2. This elastic tail compresses and dilates itself in the same way that the springs do, following the active arm, in order to create a global displacement of our swimmer.

Equations (6) and (8) can be viewed as a finite element discretization of a PDE, which describes the continuous version of our swimmer. Limit expressions for this PDE model are formally derived throughout this section while the convergence of the $N$-spring model to the continuous model will be proven in Sect. 3.

First, as $h\to 0$ ($N\to\infty$), $\displaystyle\frac{\ell_{j-1}-2\ell_{j}+\ell_{j+1}}{h^{2}}$ formally converges to a second order derivative. More precisely, we introduce a new space variable $y_{j}=(j-1)h$ for $1\leq j\leq N+1$. The points $y_{j}$ are equally spaced and thus different from the previous $x_{j}$. Since $y_{1}=0$, the $y$ variable can be seen as a local space coordinate attached to the second sphere, and we assume $\ell(y_{j})=\ell_{j}$ for a smooth enough function $\ell$. Passing to the formal limit in (6) leads to a heat equation:

[TABLE]

Concerning the boundary conditions, we first notice that $\ell_{N+1}=0$ leads to $\ell(\Lambda,t)=0$ for all $t>0$. As $h\to 0$, the equation (8) on $\ell_{1}$ formally becomes a Fourier-type boundary condition:

[TABLE]

Therefore, we finally obtain the following continuous problem:

Find $\ell\in\mathcal{C}^{2}([0,\Lambda]\times\mathbb{R}^{*}_{+})$ such that $\forall(y,t)\in(0,\Lambda)\times\mathbb{R}^{\star}_{+},$

[TABLE]

2.4 Well-posedness of the problem

Equation (15) belongs to the class of problem for which the classical theory of parabolic equations applies. Namely, calling

[TABLE]

which is a Hilbert space with the scalar product $(u,v)_{\mathcal{V}}=\int_{0}^{\Lambda}\partial_{y}u\,\partial_{y}v\,dy$, the variational formulation reads:

Let $T>0$, find $\ell(y,t)\in L_{t}^{\infty}(0,T;L_{y}^{2}((0,\Lambda)))\cap L_{t}^{2}(0,T;\mathcal{V})$ such that for all $t\in(0,T)$ and for all $v\in\mathcal{V}$

[TABLE]

with $\ell(y,0)=\ell_{0}(y)\in L^{2}((0,\Lambda))$ a given initial data.

Defining the bilinear form $\kappa$ in $\mathcal{V}\times\mathcal{V}$ as:

[TABLE]

which is symmetric and coercive on $\mathcal{V}$, well-posedness of the problem (16) follows from standard results on parabolic equations (see e.g. [16]). Moreover, it is well known that the solution $\ell(\cdot,t)$ is of class $\mathcal{C}^{\infty}([0,\Lambda])$ for any time $t>0$.

2.5 Analytical periodic solutions

Let us now solve the system (15) using the following ansatz $\ell(y,t)=\underline{\ell}(y)e^{i\omega t}$. From (14) we deduce the following equation:

[TABLE]

The characteristic polynomial associated to (18) has two roots, $\displaystyle r:=\frac{1+i}{\Lambda\sqrt{2K_{\omega}}}$ and $-r$, which leads to the following solutions:

[TABLE]

with $\alpha,\,\beta\in\mathbb{C}$.

We then determine $\alpha$ and $\beta$ using boundary conditions:

[TABLE]

i.e.,

[TABLE]

We notice that $\displaystyle r\Lambda=\frac{1+i}{\sqrt{2K_{\omega}}}$ only depends on $K_{\omega}$.

3 Convergence of the discrete model towards the continuous one

We first notice that the discrete problem (6) is a kind of non conventional mass-lumped version of a finite element discretization of the continuous one (15). In order to clarify this statement, we introduce the finite element setting. Let $\mathcal{V}_{h}\subset\mathcal{V}$ the space of continuous, piecewise linear functions $g$ on the one-dimensional partition $T_{h}=\{y_{1},\,\cdots,\,y_{N+1}\}$ of $(0,\Lambda)$, and that verify the Dirichlet boundary condition $g(\Lambda)=0$. Let $\{\Phi_{j}\}_{j=1,\,N}$ be the standard basis for $\mathcal{V}_{h}$ consisting of the hat functions defined by $\Phi_{j}(y_{k})=\delta_{j,k}$ for $1\leq j,k\leq N$.

Let $\ell_{h}\in\mathcal{V}_{h}$ be the continuous, piecewise linear function such that for $1\leq j\leq N+1$, $t>0$, $\ell_{h}(y_{j},t)=\ell_{j}(t)$. Using the basic semi-discrete Galerkine method would lead to the discretization of (21) in the matrix form:

[TABLE]

with $L_{h}(t)=(\ell_{1}(t),\,\cdots,\,\ell_{N}(t))^{T}$. Similarly, $\tilde{f}=(-\frac{\Lambda}{2}\dot{L}_{0},\,0,\,\cdots,\,0)$, $(M_{h})_{i,j}=\int\limits_{0}^{\Lambda}\Phi_{i}(y)\Phi_{j}(y)dy$ and $(K_{h})_{i,j}=\kappa(\Phi_{i},\Phi_{j})$, where $\kappa$ is defined in equation (17).

Computing explicitly the coefficients of the matrices $K_{h}$ and $M_{h}$ gives

[TABLE]

and

[TABLE]

The key observation is that Eqs. (6) and (8) are nothing but a mass-lumped discretization of (15) where the mass matrix $M_{h}$ has been replaced by the diagonal version

[TABLE]

Hence, $\ell_{h}$ actually solves

[TABLE]

together with the initial condition

[TABLE]

The classical mass-lumped method, on the other hand, would have consisted in replacing the tridiagonal mass matrix $M_{h}$ by a diagonal matrix $\bar{M}_{h}$ using an integration formula on the vertices of the partition. Namely, using the trapezoidal formula \int\limits_{0}^{\Lambda}g\sim\big{(}\frac{1}{2}g(y_{1})+\sum\limits_{j=2}^{N}g(y_{j})+\frac{1}{2}g(y_{N+1})\big{)}h=\big{(}\frac{1}{2}g(y_{1})+\sum\limits_{j=2}^{N}g(y_{j})\big{)}h, for a function $g$ satisfying $g(\Lambda)=0$ leads to the mass-lumped matrix

[TABLE]

which differs from $\widetilde{M}_{h}$.

We shall then study the ODE (22), (23) using the method presented in [17] which provides us with a convergence result for the mass-lumped method with $\bar{M}_{h}$.

We introduce the two following inner products on $\mathcal{V}_{h}$ associated with $\bar{M}_{h}$ and $\widetilde{M}_{h}$ respectively. Namely, for $(u_{h},v_{h})\in\mathcal{V}_{h}$

[TABLE]

and

[TABLE]

We also call $\|\cdot\|_{h}$ the norm associated to $(\cdot,\cdot)_{h}$, while the $L^{2}$ norm and inner products are denoted by $\|\cdot\|$ and $(\cdot,\cdot)$ respectively. Gerschgörin Theorem applied to $M_{h}$ shows the equivalence of the norms $\|\cdot\|$ and $\|\cdot\|_{h}$ on $\mathcal{V}_{h}$ uniformly in $h$, and, more precisely, we have the estimate, valid for all $v_{h}\in\mathcal{V}_{h}$

[TABLE]

from which we also deduce

[TABLE]

Finally, we introduce, for $u_{h},v_{h}\in\mathcal{V}_{h}$, $\delta_{h}(u_{h},v_{h})=(u_{h},v_{h})_{h}-(u_{h},v_{h})$ the quadrature error.

Lemma 3.1*.*

Let $u_{h},v_{h}\in\mathcal{V}_{h}$. We have, for $h$ sufficiently small:

[TABLE]

for a constant $C$ that does not depend on $u_{h}$, $v_{h}$ or $h$.

Proof.

In all what follows, $C$ denotes a constant that may vary from line to line, being always independent of $h$. Let $u_{h},v_{h}\in\mathcal{V}_{h}$. We write $|\delta_{h}(u_{h},v_{h})|\leq|(u_{h},v_{h})_{h}-\langle u_{h},v_{h}\rangle_{h}|+|\langle u_{h},v_{h}\rangle_{h}-(u_{h},v_{h})|$. Thomée [17] provides us with an estimate of the error between $\langle u_{h},v_{h}\rangle_{h}$ and $(u_{h},v_{h})$, namely,

[TABLE]

and

[TABLE]

for some constant $C>0$ that does not depend on $u_{h}$, $v_{h}$ or $h$. The latter estimate is obtained by an inverse inequality.

It remains to estimate the term $\tilde{\delta}_{h}(u_{h},v_{h})=(u_{h},v_{h})_{h}-\langle u_{h},v_{h}\rangle_{h}$.

We notice that:

[TABLE]

Similarly, (28) together with (25) gives:

[TABLE]

This yields (26) and (27). ∎

Theorem 3.1*.*

If $\ell$ and $\ell_{h}$ are solution to (16) and (22), (23) respectively, and $\ell_{0}\in H^{2}((0,\Lambda))$, we have, for all $t\geq 0$,

[TABLE]

Proof.

Let $R_{h}$ be the Ritz projector from $\mathcal{V}$ on $\mathcal{V}_{h}$, associated with $\kappa(\cdot,\cdot)$. Namely, for $g\in\mathcal{V}$, $R_{h}g$ is defined by

[TABLE]

for all $v_{h}\in\mathcal{V}_{h}$. We write $\ell_{h}-\ell=\left(\ell_{h}-R_{h}\ell\right)+\left(R_{h}\ell-\ell\right)=\theta_{h}+\rho$ (Notice that $\theta_{h}\in\mathcal{V}_{h}$). Standard estimations show that $\rho(t)$ satisfies $\|R_{h}\ell-\ell\|\leq Ch^{2}\|\partial_{yy}\ell\|$. In order to estimate $\theta_{h}$, we write, for all $\chi_{h}\in\mathcal{V}$

[TABLE]

Setting $\chi_{h}=\theta_{h}$, we obtain

[TABLE]

Here, we have at once, using Cauchy-Schwarz and Poincaré inequalities:

[TABLE]

Using the first equation of Lemma 3.1, and the fact that $\|\partial_{y}R_{h}u\|\leq C\|\partial_{y}u\|$, we also obtain

[TABLE]

from which we deduce that

[TABLE]

using the coercivity of $\kappa(\cdot,\cdot)$ on $\mathcal{V}$. We therefore infer

[TABLE]

We now recall that $\|\cdot\|_{h}$ and $\|\cdot\|$ are equivalent norms on $\mathcal{V}_{h}$, uniformly in $h$, and hence

[TABLE]

Here $\|\theta_{h}(0)\|=\left\|\ell_{0,h}-R_{h}\ell_{0}\right\|$ and

[TABLE]

whence $\theta_{h}(t)$ is bounded as desired. ∎

Theorem 3.2*.*

If $\ell$ and $\ell_{h}$ are solution to (16) and (22), (23) respectively we have, for $t\geq 0$,

[TABLE]

Proof.

We now set $\chi_{h}=\partial_{t}\theta_{h}$ in equation (31) for $\theta_{h}$ to obtain:

[TABLE]

Here, as in the proof of Theorem 3.1,

[TABLE]

Further, by the second line of Lemma 3.1,

[TABLE]

Using again the equivalence between the norms $\|\cdot\|_{h}$ and $\|\cdot\|$ on $\mathcal{V}_{h}$, we conclude:

[TABLE]

so that, after integration, and using the coercivity of $\kappa(\cdot,\cdot)$ on $\mathcal{V}$

[TABLE]

This, together with the standard estimate for $\partial_{y}\rho(t)$ completes the proof. ∎

We proved the convergence of the discrete $N$-spring swimmer to the continuous model we formally derived in the previous section. Note that we obtain only a first-order (resp. half order) convergence in $L^{2}$ norm (resp. $H^{1}$ norm) while the standard estimations for the mass-lumping method leads to a second-order (resp. first order) convergence . This is due to the Fourier-type boundary condition at [math] which differs from the classical Dirichlet boundary condition used in [17].

4 Mathematical expression of the displacement

4.1 Net displacement of the $N$-spring swimmer

We seek the swimmer’s displacement by looking at the displacement of the first of the largest sphere, meaning we only compute $V_{1}=\dot{x}_{1}$, and integrate over a period $(0,2\pi/\omega)$.

Taking into account the hydrodynamic interactions due to the $i^{\text{th}}$-sphere with $i\in\{2,\cdots,N+2\}$ on the first sphere, we have

[TABLE]

Using expressions (2) and (5), we obtain

[TABLE]

where we recall that, by convention, $\ell_{N+1}=0$.

Finally, by integrating over one period, and noticing that both $\ell_{2}$ and $\dot{L}_{1}/L_{1}$ have a vanishing time-average, we obtain, for any value of $h=\Lambda/N$, the displacement of the corresponding $N$-spring swimmer:

[TABLE]

4.2 Net displacement of the limit model

We may find an expression for the displacement of the limit model as $h$ tends to 0, by passing to the limit in the preceding expression.

Indeed, for $h$ and $y$ given, we define $j_{h}(y)$ the unique integer such that $j_{h}(y)h\leq y\leq(j_{h}(y)+1)h$. Then, defining $\chi_{h}$ the function

[TABLE]

we may write

[TABLE]

where $\ell_{h}$ is the piecewise linear function defined in the previous section.

Finally, the displacement $\Delta_{h}x_{1}$ of the $N$-spring swimmer during one time period can be rewritten as

[TABLE]

Now, using the fact that $j_{h}(y)h\to y$ when $h\to 0$, together with the $L^{2}$ and $H^{1}$ convergence of $\ell_{h}$ to $\ell$, we obtain that, for any $y$ and $t$,

[TABLE]

Moreover $0\leq\chi_{h}(y,t)\leq\max_{t}\frac{1}{L_{0}(t)}=\frac{1}{L(1-\tilde{\varepsilon})}$, shows that $\chi_{h}$ is uniformly bounded.

Therefore, using dominated convergence theorem, we deduce that $\chi_{h}$ converges to $\chi$ in $L^{2}(0,2\pi/\omega;(0,\Lambda))$ as $h$ tends to 0.

Using the convergence theorems proven in the preceding section, we may pass to the limit $h\rightarrow 0$ in $\Delta_{h}x_{1}$, and obtain the following expression for the displacement during one period for the limit model

[TABLE]

5 Numerical experiments

In this section, we numerically study the discrete model’s convergence towards the continuous one. Then, we investigate the influence of the two parameters $\omega$ and $\tilde{\varepsilon}$ on the system and on its displacement, while the rest of the swimmer is determined by the values in table 1. All simulations are achieved using Matlab. We consider here that the default length $L$ of the active arm is small compared to the rest of the swimmer. The first sphere thus acts like the head of a sperm cell, and the active arm like a link between the head and the flagella, which gives a signal so that the rest of the system oscillates.

5.1 Convergence of the discrete models to the continuous one

We investigate numerically the convergence estimations obtained in section 3. We recall that the continuous solution $\ell$ solves the heat equation PDE with the Fourier-type boundary conditions (15). We consider, in this section, periodic forcing for which explicit solutions are given by (19, 20).

5.1.1 Convergence of the $N$-spring discrete model

We recall that the discrete solution $\ell_{h}$ is the $P^{1}$ discrete function based on the $(\ell_{i})_{i}$ solution to the $N$-spring ODE system (6,7,8). This discrete system corresponds to a semi-discretization in space of the continuous model, based on a non conventional mass-lumping method. The solution $(\ell_{i})_{i}$ of the discrete problem in the periodic setting is given in equations (9,10,12,13).

The space step $h$ (or equivalently the number of springs $N$) being given, the discrete error is defined as the error between $\ell_{h}$ and the $P^{1}$ interpolation of $\ell$. We plot in figure 3, the $L^{2}$ (resp. $H^{1}$) error denoted by $e_{h,L^{2}}$ (resp. $e_{h,H^{1}}$).

We observe that the $L^{2}$ error converges with order one, as expected from theorem 3.1. Concerning the $H^{1}$ error, we observe a superconvergence phenomenon: as the $L^{2}$ error, it converges at order 1, while theorem 3.2 predicts a convergence at order 1/2. This can be explained by the regularity of the considered periodic solution.

5.1.2 Influence of mass-lumping

As mentioned earlier, the $N$-spring model turns out to be a discretization in space of the continuous problem (15), based on an unconventional mass-lumping method. The convergence proof that we proposed in section 3 is based on the results of Thomée[17]. He shows that, for a standard mass-lumping discretization, the usual order of convergence for finite elements is obtained: convergence of order 2 for the $L^{2}$ error and 1 for the $H^{1}$ error.

We investigate here the influence of the space discretization, by comparing the $N$-spring model (22), solved numerically this time, to the classical mass-lumping method (24) and the standard Galerkin finite element method (21). Again we consider the periodic framework for which the exact solution is available. The time discretization of the three ODE systems is achieved using a Crank-Nicolson scheme for which the time step is chosen to be small enough so that the error due to the time discretization is negligible.

The corresponding $L^{2}$ (resp. $H^{1}$) error is given on figure 4 (resp. figure 5). We can see that, as expected, the $L^{2}$ error converges at order 1 for the $N$-spring model, while it converges at order 2 for both the classical mass-lumping method and the standard Galerking discretization. Again, due to the regularity of the solutions, a super-convergence phenomenon of the $H^{1}$ error is observed for all three methods: as the $L^{2}$ error, it converges at order 1 for the $N$-spring model and order 2 for the other two discretizations.

5.2 Swimming strokes

In this section, we investigate the swimming ability of the $N$-spring swimmer. The stroke being periodic, we use the explicit solutions given in section 2.2. The computations are achieved for $N=2\,000$ springs.

5.2.1 Deformation of the swimmer

Figure 6 shows a full stroke of the swimmer, in which we notice that a wave is propagating along its tail. Remember that this wave is a contraction wave along the horizontal tail. This tail appears to be oscillating fairly efficiently for the side close to the head, while the amplitude of the contraction decays considerably on the second half of the tail.

The movement shown corresponds to the stretch of $\ell_{j}$, and not to the actual deformation which would be $\ell_{j}/N$, for all $1\leq j\leq N$. We thus remark that this deformation is relatively small compared to the size of the artificial swimmer, which matches the approximation of small deformations that we made in the first place.

5.2.2 Displacement

In this section, we study the influence of the parameters $\tilde{\varepsilon}$ and $K_{\omega}$ on the swimmer’s displacement  (35), in order to maximize its absolute value.

In figure 7, we plot the displacement of the swimmer as a function of time, for different values of $\tilde{\varepsilon}$. The displacement is computed through numerical integration of equation (35). The graph shows that the swimmer globally swims backwards, and we recognize the back and forth motion which is characteristic of low Reynolds number artificial swimmers. A larger amplitude $\tilde{\varepsilon}$ of the forcing leads to a larger displacement and we observe (see figure 8), that $\Delta x_{1}$ is proportional to $\tilde{\varepsilon}^{2}$, which is what is expected in theory (similar behaviors are observed, e.g., in [2, 7, 1] and explained as the surface of loops in the space of shapes [11]).

As we want to maximize $\Delta x_{1}$ while having $\tilde{\varepsilon}<1$, we choose a fixed value $\tilde{\varepsilon}=0.7$ which, although arbitrary, allows for an easier comparison to Montino and DeSimone’s results [14], as they made a similar parameter choice.

Figure 9 shows $\Delta x_{1}$ depending on $K_{\omega}$, for different values of $\tilde{\varepsilon}$. At any fixed $K_{\omega}$, we observe once again that larger $\tilde{\varepsilon}$ leads to larger $\Delta x_{1}$. We first observe that, if $K_{\omega}\to\infty$, the net displacement of the swimmer vanishes. According to the expression of $K_{\omega}$, this is the case for example when $\omega\to 0$: the oscillation disappears, immobilizing the artificial swimmer. This can also happen when $\tilde{k}\to\infty$: the springs become so rigid that the tail of the swimmer can no longer deform. In that case, the swimmer has only one degree of freedom left to deform and faces Purcell’s Scallop theorem’s obstruction. Similarly, letting $K_{\omega}\to 0$ immobilizes the swimmer. An optimal value $K_{\omega}^{\text{opt}}$ for the non-dimensional parameter is reached between these two limiting cases, in order to maximize the displacement on one time period. According to the figure, $K_{\omega}^{\text{opt}}\simeq 0.3765$. A complete mathematical expression of $K_{\omega}^{\text{opt}}$ does not seem available, due to the largely nonlinear nature of the problem contrarily to the final expression obtained in [14]. A pair of optimal values for $\omega$ and $k$ to obtain this $K_{\omega}^{\text{opt}}$ are $\omega=1\,rads^{-1}$ and $\tilde{k}\simeq 6.207\cdot 10^{-8}Nm^{-1}$. Moreover, the expression of $K_{\omega}$ guarantees that $\omega$ must vary proportionally to $\tilde{k}$ for the pair $(\tilde{k},\omega)$ to remain at the optimum.

Indeed, the softer the spring, the slower the first arm needs to oscillate in order to generate a large movement.

Looking at the other parameters separately, we can also clearly see from equation (35), that the displacement depends linearly on $\tilde{a}$, which is predictable. However, this parameter has a direct consequence on the size of the artificial swimmer and must stay in a reasonable range (in our case no more than $1e-5\,m$) so that the swimmer stays at microscopic scale.

Finally, we notice that the value of $\Lambda$ and the ratio $a_{1}/\tilde{a}$ has little to no influence on the previous analysis. We therefore keep for those parameters values that seem coherent with the scale we are working at, and that match with numerical experiments provided in [14].

6 Conclusion

We analyzed the dynamics of two low-Reynolds-number swimmers. The first one, which is an extension of [14], is made of $N$ passive springs, and the second one is the corresponding limit model with an elastic tail. Both are activated by an active arm that elongates and retracts periodically with amplitude $\varepsilon$ and angular frequency $\omega$.

Noting that the $N$-spring swimmer is a non-conventional mass lumping discretization of the limit model, we proved its convergence, when $N$ tends to infinity, to the continuous model, by extending the results of Thomée [17] to the case of a Fourier-type boundary condition.

For both swimmers, a phase difference between the oscillations of the active arm and the tail is created by the interaction between elastic and hydrodynamic forces. Then, both models undergo non-reciprocal shape changes and thus circumvent Scallop Theorem’s obstruction [3]. Numerical simulations indeed show a wave propagating along the swimmers’ tails. Similarly to what was shown in [14], our models are able to swim but there is no control over the swimming direction.

Then, we focused on computing the net displacement over a time period of the swimmer in both cases, in view of its optimization. We obtain explicit formulae for this displacement as a function of the local elongation during the stroke. We numerically recover the classical back and forth swimming and the second-order scaling of the displacement as a function of the maximum elongation of the forcing active arm. Moreover, we highlight a dimensionless parameter $K_{\omega}$, driving the movement of the swimmer when its geometry ($\Lambda$, $a$, $a_{1}$) is given. Some optimal values for this parameter can be estimated by numerical experiments.

Lastly, we noticed that, similarly to the behavior of Machin’s swimming rod [18], the deformations of both our swimmers is rapidly attenuating along their passive parts, which suggests that some form of activation is needed in order to mimic the type of behavior observed in nature.