Open Access Article
This Open Access Article is licensed under a Creative Commons Attribution-Non Commercial 3.0 Unported Licence

Designing optimal elastic filaments for viscous propulsion

Mariia Dvoriashyna *ab and Eric Lauga *a
aDepartment of Applied Mathematics and Theoretical Physics, University of Cambridge, Wilberforce Rd, Cambridge CB3 0WA, UK. E-mail: m.dvoriashyna@ed.ac.uk; e.lauga@damtp.cam.ac.uk
bSchool of Mathematics and Maxwell Institute for Mathematical Sciences, University of Edinburgh, Peter Guthrie Tait Rd, Edinburgh EH9 3FD, UK

Received 22nd November 2024 , Accepted 20th March 2025

First published on 8th April 2025


Abstract

The propulsion of many eukaryotic cells is generated by flagella, flexible slender filaments that are actively oscillating in space and time. The dynamics of these biological appendages have inspired the design of many types of artificial microswimmers. The magnitude of the filament's viscous propulsion depends on the time-varying shape of the filament, and that shape depends in turn on the spatial distribution of the bending rigidity of the filament. In this work, we rigorously determine the relationship between the mechanical (bending) properties of the filament and the viscous thrust it produces using mathematical optimisation. Specifically, by considering a model system (a slender elastic filament with an oscillating slope at its base), we derive the optimal bending rigidity function along the filament that maximises the time-averaged thrust produced by the actuated filament. Instead of prescribing a specific functional form, we use functional optimisation and adjoint-based variational calculus to formally establish the link between the distribution of bending rigidity and propulsion. The optimal rigidities are found to be stiff near the base, and soft near the distal end, with a spatial distribution that depends critically on the constraints used in the optimisation procedure. These findings may guide the optimal design of future artificial swimmers.


1 Introduction

The physics of motile microorganisms, such as bacteria, microalgae and spermatozoa has recently been a subject of active research at the intersection between physics, mathematics and biology.1–5 Studying cell motility is important not only for our general understanding of biological and biophysical phenomena,6,7 but also for the development of biomimetic, synthetic microrobots in a myriad of potential applications, including drug delivery, smart surgery, sensing and detoxification.8–11

Inspired by the actuation and geometry of motile microorganisms, different designs of swimming microrobots have been developed,12,13 coming in a variety of shapes (e.g. rigid helical filaments, straight or helical elastic filaments, collection of beads), actuation mechanisms (e.g. rotation of a helical shape, planar or circular oscillation of an elastic tail, prescribed shape changes) and materials (cobalt composite on helical lipids, nickel or iron-filled carbon nanotubes, magnetic beads linked with DNA strands).10,14–18 For example, recently developed undulatory multilink nano-swimmers are controlled by a planar oscillating magnetic field, which sets elastic tails into oscillatory motion and creates forward propulsion, similar to the swimming of spermatozoa.17 Another example includes artificial bacterial flagella, in which a rigid helical shape attached to a thin soft-magnetic head rotated by an external rotating magnetic field, resulting in a corkscrew-type propulsion similar to that of flagellated bacteria.18

Biological microorganisms swim on small scales, with relevant Reynolds numbers ranging between [scr O, script letter O](10−5) and [scr O, script letter O](10−1).19 In this regime, inertia of both the swimmer and the surrounding fluid are negligible, and swimming strategies have to rely on viscous (Stokesian) forces to generate propulsion. One of the simplest swimming modes capable of creating propulsive forces at zero Reynolds number is the classical ‘flexible oar’, initially proposed by Purcell in his famous lecture on ‘life at low Reynolds numbers’.20 In this setup, a flexible tail is made to passively oscillate in a viscous fluid, generating travelling waves along the filament (whose amplitude tends to decay exponentially if the forcing is localised along the filament) and resulting in propulsion in the direction opposite to that of the wave propagation.21,22

This mode of swimming has been used experimentally in multiple artificial swimmers. For instance, Guo et al.23 developed a (cm)-scale fish-like microbot, with a main body made of wooden and styrol materials placed on a permanent magnet, and the fin made of a polyimide film sheet. Directly inspired by spermatozoa, Williams et al.24 developed a (mm)-scale biohydrid swimmer made of PDMS and actuated by contractile cardiomyocyte cells lining the base of the filament. Magdanz et al.25,26 developed IRONSperm, which comprised of bovine sperm cells (≈60 μm long) coated with a suspension of 100-nm rice grain-shaped maghemite nanoparticles and controlled by a magnetic field. A related device was designed by Celi et al.,27 composed of a superparamagnetic head and a flexible gold/polypyrrole (Au/PPy) flagellum, and actuated by a magnetic field. Mushtaq et al.28 developed a nano-eel, a multifunctional piezoelectric tailed nanorobot, which consisted of a smart flexible tail made of a polyvinylidene fluoride-based copolymer linked to a PPy nanowire head, decorated with nickel (Ni) rings for magnetic actuation. Theoretical modelling can help understanding physics underlying propulsion of these artificial microswimmers and propose ways of enhancing it.

From a theoretical standpoint, and following the pioneering contributions of Machin,21 the dynamics and propulsion of an elastic filament with constant thickness oscillating in a viscous fluid were modelled physically and mathematically in classical studies.22,29 These original models relied on a small-amplitude approximation, which allows linearisation of the equations of motion for the shape of the elastic filament. The propulsion is generated as a result of the quasi-steady balance of viscous and elastic forces. These linear models have been validated experimentally, using a (cm)-scale robot with an elastic tail actuated at its base, immersed in high viscosity silicone oil.30 Beyond the small-amplitude approach, the full elastohydrodynamic problem of an oscillating filament was studied numerically using the lattice-Boltzmann method,31 asymptotic coarse-grained elastohydrodynamics,32 and regularised Stokeslets,33 while further work also addressed the non-linear dynamics of biological filaments resulting from their molecular structure.34,35

With a view on applications in micro-scale robotics, designing a hydrodynamically optimal synthetic swimmer, i.e. one that can self-propel with the highest swimming velocity for a given energy expenditure, is a problem of fundamental interest. Biologically, this optimisation problem has received a lot of attention in the context of motile eukaryotic cells. It was shown that spermatozoa have an optimal flagellum-to-head length ratio, with the most efficient swimming mode involving symmetrical travelling waves in the flagellum.36 The travelling waves were shown to lead to the optimal motion in an active filament,37–39 while the mathematically optimal shape of the wave was shown to be sawtooth-like,40 with shape singularities that may be regularised by additional biophysical constraints.41,42 Recently, by making use of the new BOSO-Micro database,43 which includes a comprehensive collection of experimental data from the literature on the swimming speed and morphological characteristics of unicellular organisms, it was shown that amplitude-to-wavelength aspect ratios of flagellar eukaryotes are very close to the hydrodynamic optimum.44 Additionally, optimal head-to-tail ratios were studied for a model swimmer equipped with elastic filament actuated at the base,45 and the optimal shape of the head was investigated for a synthetic swimmer with magnetic head and an elastic tail.46

Focusing on the specific case of actuated elastic filaments, past work demonstrated that tapering an elastic filament may lead to an increase in viscous propulsion.47 For a filament with a circular cross-section of radius r(s), which depends on the arc-length s along the filament (see e.g.Fig. 1), the bending rigidity of the filament, A, scales as Ar4. Thus, the tapered filaments are more rigid near the actuation point than the free end. Further studies considered specific functional forms for the filament radius r or the bending rigidity A, such as linear, quadratic, and exponential functions of s.47–50 In most cases, the propulsive force was shown to increase if the filament is stiffer at the base and softer at the distal end. A recent experimental investigation suggested that a microrobot equipped with DNA-based filament wrapped in a tapered bundle swims faster than the one with a bundle of constant thickness.51 These studies reveal a nontrivial relationship between the bending rigidity of the filament (in particular, how it cross-section varies) and the magnitude of the propulsion generated.


image file: d4sm01388c-f1.tif
Fig. 1 Sketch of the model problem optimised in this paper: an elastic filament of variable thickness is made to oscillate about its hinged base. At each instant, its shape is obtained as the balance between the viscous forces ([f with combining macron]vis) and the elastic forces ([f with combining macron]elastic). We use variational calculus, along with an adjoint field, to determine the bending rigidity leading to maximal propulsion.

In this work, we propose to rigorously determine the relationship between the mechanical properties of the filament and the produced thrust. By considering a specific model system consisting of a beating elastic filament with an oscillating slope at its base, we determine the optimal bending rigidity function that maximises the time-averaged thrust produced by the actuated filament. Instead of prescribing a specific functional shape, we use functional optimisation and variational calculus to establish formally the link between the bending rigidity function and propulsion for different sets of constraints. We find that the optimal bending rigidity for this model swimmer is invariably stiff at the base of the filament and soft at the distal end.

This paper is structured as follows. We first formulate the model problem of a beating elastic filament with an oscillating slope at its base in Section 2.1, and simplify it in Section 2.2. We then pose in Section 2.3 the optimisation problem consisting in finding the cross-sectional shape that maximises the propulsive force generated by the filament. It is solved using an adjoint approach, which is implemented numerically in a steepest ascent algorithm in Section 2.4. We present the results of the optimisation procedure in Sections 3 and 4 and conclude with a discussion in Section 5.

2 Model, methods and optimisation

2.1 Model system: oscillating slender elastic filament

We study in this paper a well-defined model system, which will allow us to solve the optimisation problem rigorously.
2.1.1 Geometry. We consider the case of an oscillating slender elastic filament whose slope at its base is made to oscillate in time at a prescribed frequency (see sketch in Fig. 1). Using bars to denote dimensional variables, the filament is aligned on average with the [x with combining macron] axis and its base is located at [x with combining macron] = 0. With [s with combining macron] denoting the arclength along the filament, oscillations of the slope of the filament at [x with combining macron] = [s with combining macron] = 0 result in a planar wave propagating between the filament base and its free end located at [s with combining macron] = L. We denote the position of the centreline of the filament at time [t with combining macron] by ȳ([s with combining macron],[t with combining macron]), while the circular cross-section of the filament has radius [r with combining macron]([s with combining macron]). The filament is assumed to be slender, so [r with combining macron]([s with combining macron]) ≪ L everywhere. The goal of the paper is to determine the function [r with combining macron]([s with combining macron]) which maximises the propulsion generated by the filament.
2.1.2 Equations of motion. In order to tackle the problem analytically, we consider a linearised version of the elastohydrodynamic problem, assuming that the displacements of the filament are always small, so that [s with combining macron][x with combining macron]. This approach captures the physics of the problem, and is known to remaining quantitatively accurate even beyond the strict small-amplitude regime of validity.30

In the absence of inertia, the shape of the filament is determined by a quasi-steady balance between hydrodynamic and elastic forces. To compute the viscous forces, we exploit the slenderness of the filament and use resistive force theory (RFT), which relates the hydrodynamic force density at any point along the filament to its local velocity via drag coefficients in the directions orthogonal and parallel to the local tangent to the filament, ζ = 4πμ/(ln(L/[r with combining macron]) + 1/2), ζ = 2πμ/(ln(L/[r with combining macron]) − 1/2),52 where μ denotes the dynamic viscosity of the fluid. Within the assumptions of RFT, the viscous force density per unit length of the filament at each point can be written as

 
[f with combining macron]vis = −ζ[v with combining macron] − (ζζ)t(t·[v with combining macron]),(1)
where t is the local tangent vector, and [v with combining macron] is the local velocity distribution along the swimmer.19 For small-amplitude beating, the viscous force occurs primarily in the direction orthogonal to the flagellum axis (i.e. the ȳ direction), with density
 
[f with combining macron]vis = −ζ[v with combining macron],(2)
where image file: d4sm01388c-t1.tif is the instantaneous velocity in ȳ direction.

To compute the elastic force density, [f with combining macron]elastic, we make use of classical elastic beam theory.22,53 For displacement of the filament written as ȳ([x with combining macron]), the leading-order elastic force in the ȳ direction is given by

 
image file: d4sm01388c-t2.tif(3)
where Ā([x with combining macron]) = ([x with combining macron]) is the bending stiffness, a product of the material's Young's modulus, E, (assumed to be constant) and the second moment of cross-sectional area, Ī([x with combining macron]) = π[r with combining macron]([x with combining macron])4.

Imposing a quasi-steady balance between the viscous and elastic forces, eqn (2) and (3), we obtain the governing (classical) hyper-diffusion equation that describes the position of the centreline of the filament in time, in the linearised limit, as

 
image file: d4sm01388c-t3.tif(4)

Note that the inextensibility of the filament comes in at higher order in the oscillation amplitude and, hence, does not appear in the linearised version above.

2.1.3 Boundary conditions. In the specific model problem considered in this paper, we prescribe that the slope of the filament at its base oscillates with frequency ω and impose that the position of the filament is fixed at the base. At the free end, we adopt force and torque-free conditions. The four boundary conditions (BCs) therefore read
 
ȳ(0,[t with combining macron]) = 0 fixed base,(5a)
 
image file: d4sm01388c-t4.tif(5b)
 
image file: d4sm01388c-t5.tif(5c)
 
image file: d4sm01388c-t6.tif(5d)

Therefore, in the linearised limit ε ≪ 1, the leading-order position of the filament centreline satisfies eqn (4) along with the BCs in eqn (5).

2.1.4 Propulsive force. In order to compute the propulsive force, we need to calculate the component of the viscous force experienced by the filament in the [x with combining macron] direction, i.e. perpendicular to the beating direction. Although the drag coefficients ζ and ζ appearing in eqn (1) depend on the local cross section of the filament, their dependence is logarithmic; we may thus treat them as constant in the slender limit.

Within the small-amplitude approximation we can write t ≈ (1,∂ȳ/∂[x with combining macron]) and [v with combining macron] = (0,∂ȳ/∂[t with combining macron]). The local force in the [x with combining macron] direction is then classically given by image file: d4sm01388c-t7.tif (please see ref. 19 and 54 for detailed derivation). Since waves are expected to propagate in the positive [x with combining macron] direction, and propulsion to occur in the opposite direction, the magnitude of the time-averaged propulsive force experienced by the filament, [F with combining macron], is obtained by taking a double integral in time and space of −[f with combining macron]x, so that

 
image file: d4sm01388c-t8.tif(6)
With this sign convention, the magnitude of the force is positive (and the net force acts in the negative [x with combining macron] direction).

2.2 Model simplification

2.2.1 Non-dimensionalisation. The first step in simplifying the model consists in non-dimensionalising the problem. The relevant time scale is given by the inverse frequency of oscillations, while the relevant length scales are the length of the filament L along the [x with combining macron] direction, and the small-amplitude motion of magnitude εL along ȳ. A relevant scale A0 for the rigidity can be defined using the length of the filament, A0 = L4ωζ, so that Ā = A0A. This choice corresponds to a relevant dimensionless Sperm number, defined as Sp = L(ζω/A0A)1/4, of one when A = 1.19

We thus define the dimensionless variables

 
y = ȳ/εL, t = [t with combining macron]ω, x = [x with combining macron]/L, A = Ā/A0,(7)
and the dimensionless version of (4) is now
 
image file: d4sm01388c-t9.tif(8)
The time averaged propulsive force [F with combining macron] in eqn (6) now scales with relevant magnitude (ζζ)ε2L2ω and the dimensionless force, F, is given by
 
image file: d4sm01388c-t10.tif(9)

2.2.2 Separation of variables. To further simplify the mathematical problem, we exploit the linear dynamics and write the time-varying amplitude of the dimensionless filament as
 
y(x,t) = f(x)cos[thin space (1/6-em)]t + g(x)sin[thin space (1/6-em)]t.(10)
Substituting the expression for y in eqn (10) into the force balance in eqn (8), we obtain two coupled equations for f(x) and g(x) as
 
f + (A(x)g′′)′′ = 0,(11a)
 
g + (A(x)f′′)′′ = 0.(11b)
We will denote with [capital script C]1 and [capital script C]2 the left hand sides of eqn (11a) and (11b), respectively.

Using eqn (10), the propulsive force from eqn (9) can be then simplified as

 
image file: d4sm01388c-t11.tif(12)

Finally, the non-dimensional boundary conditions (5) transform into

 
f(0) = f′(0) = f′′(1) = f′′′(1) = 0,(13a)
 
g′(0) = 1, g(0) = g′′(1) = g′′′(1) = 0.(13b)

2.3 Optimisation of bending rigidity

In this section we outline the mathematical formalism used to determine the bending rigidity of the filament leading to maximal propulsion.
2.3.1 Dimensional optimisation problem. The dimensional optimisation problem consists in maximising the magnitude of the time-averaged propulsive force [F with combining macron] (which is positive) over all possible distributions of bending rigidities, Ā, subject to the equations for the filament dynamics, i.e.
 
image file: d4sm01388c-t12.tif(14)
2.3.2 Dimensionless optimisation problem. In dimensionless variables, the optimisation problem in eqn (14) can be rewritten as
 
image file: d4sm01388c-t13.tif(15)

The functional optimisation problem in eqn (15) in general needs to be solved numerically.

2.3.3 Special case: constant bending rigidity. If we consider the special case of a spatially homogeneous A, the optimisation can be carried out analytically, as shown in previous work.22 In Fig. 2, we plot the dependence of the resulting dimensionless propulsive force on the value of the bending rigidity (A). There is almost no propulsion for very flexible (small A) or very rigid filaments (large A), and an optimal value is A ≈ 0.0794. In the inset we also display the time-varying position of the centreline of the filament, y(x,t), for this optimal choice of A.
image file: d4sm01388c-f2.tif
Fig. 2 Time-averaged (dimensionless) propulsive force in the case of a constant bending rigidity (A). Maximum propulsion for a hinged filament is F ≈ 0.06 at the optimal value A ≈ 0.0794. Inset: Time-varying position of the hinged filament y(x) at different points in time for this optimal value of A: t = 0 is depicted with a thick black line and the following time steps are shown with decreasing greyscale.
2.3.4 Functional derivative and adjoint functions. To solve the functional optimisation problem (15) without introducing any assumptions on the specific shape of A, we will use variational calculus and will compute the functional derivative of dF/dA with the help of two adjoint functions.

To derive the adjoint problem we consider small perturbations, AA + δA, ff + δf and gg + δg, of both [capital script C]1(11a) and [capital script C]2(11b). Keeping only the first-order terms we obtain

 
δf − (δAg′′)′′ − (Aδg′′)′′ = 0,(16a)
 
δg + (δAf′′)′′ + (Aδf′′)′′ = 0.(16b)
We next write eqn (16a) in weak form, multiplying by a test function h and integrating over the domain
 
image file: d4sm01388c-t14.tif(17)
Integrating by parts four times the second term of the left-hand side, and two times the right-hand side, and applying boundary conditions (13) leads to
 
image file: d4sm01388c-t15.tif(18)

We will specify the boundary terms later.

Following a similar calculation, the weak form of eqn (16b) with a test function j can be written as

 
image file: d4sm01388c-t16.tif(19)

Taking the difference between eqn (18) and (19) we obtain

 
image file: d4sm01388c-t17.tif(20)

The variation of the propulsive force FF + δF in eqn (12) at linear order reads

 
image file: d4sm01388c-t18.tif(21)
which after integration by parts becomes
 
image file: d4sm01388c-t19.tif(22)

In order to compute the functional derivative, we need to make an explicit link between δA and δF. To do this, we equate the left-hand side of (20) with the right-hand side of (22), by choosing the adjoint functions h and j to satisfy the field equations

 
h − (Aj′′)′′ = g′,(23a)
 
j + (Ah′′)′′ = f′,(23b)
associated with the boundary conditions
 
h(0) = h′(0) = h′′(1) = 0, (Ah′′)′(1) = f(1)/2,(24a)
 
j(0) = j′(0) = j′′(1) = 0, (Aj′′)′(1) = −g(1)/2.(24b)
With this choice of h and j, we obtain an explicit equation linking δF with δA and eqn (22) becomes
 
image file: d4sm01388c-t20.tif(25)
The functional derivative of the propulsion F in the space of f and g functions is given by
 
dF/dA := f′′j′′ + g′′h′′,(26)
with the four functions satisfying eqn (11) (physical problem) and (23) (adjoint problem) with boundary conditions in eqn (13) and (24).

2.4 Constraints on bending and numerical algorithm

2.4.1 Singular solutions. To compute the solution to the optimisation problem, we use the functional derivative obtained in eqn (26) and the analogue of steepest ascent in functional space to find the optimal solution. In other words, we set up an iteration procedure, and at every iteration step follow the functional derivative to increase the value of F.55 Without additional constraints, however, computations show that the solution always ends up being singular, with the bending rigidity (and hence the thickness of the filament) becoming zero at some point along the filament. Note that for sufficiently soft filaments (high Sp numbers) the assumptions of our linearised model may become invalid. Thus, additional constraints must be considered to regularise the problem and to keep the validity of our approach. Here, we consider two such constraints, motivated by the design and/or manufacturing of the microswimmers, and we use a modified steepest ascent algorithm to account for these constraints.
2.4.2 Constraints C1 and C2. In the first constraint C1, we assume that the function A(x) is bounded below and above by constant fixed values, a and b respectively.

In the second constraint C2, A is bounded below (by a) and the filament has a fixed volume. Recalling that Ā = Eπ[r with combining macron]4, we can rewrite it in dimensionless form as A = Ẽr4 with image file: d4sm01388c-t21.tif, where r0 is a characteristic thickness of the filament. The dimensionless volume of the filament is image file: d4sm01388c-t22.tif. The constraint for fixed volume, V, can thus be written as image file: d4sm01388c-t23.tif.

The two sets of constraints we consider are therefore

• (C1) aA(x) ≤ b; or

• (C2) aA(x) and image file: d4sm01388c-t24.tif.

2.4.3 Numerical solution. To solve the optimization problem with constraints C1 or C2 numerically, we need to modify the steepest ascent algorithm. For C1, a projection function is introduced to ensure that the solution at each iteration remains within the fixed bounds a and b. To address C2, the projection function is adapted to enforce the condition aA(x), while the objective function is adjusted to incorporate the fixed volume constraint. Detailed descriptions of the algorithms used to handle both constraints are provided in Appendix A.

3 Optimal shapes with constraint C1 (bounded bending rigidity)

In this section, we use numerical simulations and theoretical analysis to determine the filament shapes that lead to optimal propulsion. Here we consider the shapes that follow the constraint C1 and thus whose bending rigidity must remain within a prescribed interval.

3.1 Optimal shapes have piece-wise constant rigidity

We first run numerical simulations for several chosen values of [a,b] and for an initial smooth profile satisfying the boundary conditions A1(0) = b, A1(1) = a with image file: d4sm01388c-t25.tif. An example where the shape is initially chosen to have a cubic profile is shown in Fig. 3 as dashed line in the case [a,b] = [0.01,0.1], an interval that contains the optimal constant value from Fig. 2.
image file: d4sm01388c-f3.tif
Fig. 3 Typical numerical solution of the optimisation problem when the dimensionless rigidity is required to remain within the interval A ∈ [0.01,0.1] (bounded constraint C1). Initial condition for the function A(x): dashed line; final optimal function A(x): solid line. Inset: Time-varying position of the centreline of the filament y(x,t) at different times (greyscale) for the optimal A(x).

We applied the algorithm described in Section A.1 setting the tolerances to δ1 = δ2 = 10−5. Independently of the shape taken for the initial condition, we systematically find that the final solution for the optimisation problem is a step-like function, illustrated with solid line in Fig. 3, that takes value A = b at the proximal end, the value of A = a at the distal end and undergoes a sudden jump between a and b at a dimensionless point x0; for example, we have x0 ≈ 0.47 for [a,b] = [0.01,0.1]. When the filament is built of the same material throughout, a step-like function for the bending rigidity corresponds to a piecewise-constant radius for the filament. The corresponding spatial positions of the filament y(x,t) at different time points are shown in the inset of the Fig. 3. The resulting propulsive force is F ≈ 0.116, which is almost twice higher than the optimal force in the case of uniform rigidity (Fig. 2).

To test whether it is the only optimal solution for given [a,b], we run the algorithm for 20 random initial conditions (represented as a truncated Fourier series with random modes) and always obtain the same final step-like solution. For different choices of a and b, the optimal solutions are also invariably step functions, and the values of the jump points x0 depend on the values of both a and b.

3.1.1 Interpretation of the optimal solution. An intuitive mathematical explanation for this optimal solution can be gained by examining the shape of the functional derivative in the case of spatially uniform bending profile, A. When A is constant, the adjoint problem in eqn (23) and (24) can be solved analytically and we can calculate the functional derivative of F with respect to A exactly; we show the results in Fig. 4 for different dimensionless values of A. Clearly, the derivative dF/dA is always positive at x = 0, so in the optimisation steps we add positive values to A around x = 0 until we reach the upper bound b. On the other hand, near x = 1 the functional derivative is always negative, and this drives the repeated reduction of the bending stiffness A near the distal end of the filament with each iteration step until we reach the lower bound a.
image file: d4sm01388c-f4.tif
Fig. 4 Normalised functional derivative dF/dA in the case of a constant dimensionless bending rigidity A (three different values, see inset).

Physically, our solution indicates that the front of the optimal filament tends to be stiff, which creates a finite amplitude for the beat (akin to a filament with clamped end). The distal part, in turn, is much softer, which leads to an asymmetry in the stroke and generates propulsion.

3.2 Analytical solution for piece-wise constant A

The numerical results suggest that the optimal bending profile takes the form of a piece-wise constant value. For this type of profile, we can in fact compute a fully analytical solution for the filament shape y(x,t).

Assuming that A takes form A = b + (ab)H(xx0), where H is a Heaviside function (i.e. H(x) = 0 for x ≤ 0 and 1 for x > 0), we can split the problem in eqn (8) into two separate equations on the left and right of the jump point x0

 
image file: d4sm01388c-t26.tif(27a)
 
image file: d4sm01388c-t27.tif(27b)

The boundary conditions at the point of discontinuity of A are continuity of function y and its derivative, i.e.

 
image file: d4sm01388c-t28.tif(28)
and continuity of the elastic torque and force,
 
image file: d4sm01388c-t29.tif(29)
These are accompanied by the boundary conditions on the right and left boundaries as in eqn (5).

Looking for solutions of the form y = Re(ỹeit), these equations are simplified to

 
iỹ1 = −bỹ(iv)1, xx0,(30a)
 
iỹ2 = −aỹ(iv)2, x > x0.(30b)

The general solution to eqn (30) is

 
1 = C1eαx + C2eαx + C3eiαx + C4eiαx,(31a)
 
2 = D1eβx + D2eβx + D3eiβx + D4eiβx,(31b)
where image file: d4sm01388c-t30.tif and image file: d4sm01388c-t31.tif. Applying the boundary conditions from eqn (28) and (29), we obtain a linear system for the values of Ci and Di, which is easily inverted numerically.

The propulsive force (9) can then be computed analytically using the dynamics in eqn (27), boundary conditions (28) and (29) and integration by parts,

 
image file: d4sm01388c-t32.tif(32)
where we used stars to denote complex conjugate.

3.3 Maximum propulsive force for different bounds [a,b]

The existence of this analytical solution makes it straightforward to mathematically study the impact of varying the values of the bounding constraints a and b on the propulsive force. In Fig. 5a we show the dependence of the maximum propulsive force F on different dimensionless values of a and b, both varying from 10−3 to 10. The corresponding value for the jump point x0 is displayed in Fig. 5b.
image file: d4sm01388c-f5.tif
Fig. 5 Optimisation of a filament with piecewise constant bending rigidity A. (a) Maximum dimensionless F (colours) as a function of dimensionless interval limits a and b. Inset: Maximum value of F for fixed b and variable a. The lines correspond to the lines of the same colour in the main plot: red is for b = 1, cyan for b = 0.1 and magenta for b = 0.01. (b) Corresponding values of the optimal transition point x0. Inset: The position of the centreline of the filament y(x) at different times for a = 10−3 and b = 1 (denoted by a red dot on the plot); grey scale corresponds to different times during the period of oscillation (starting with thicker black line and fading to lighter grey colours).

We first observe from Fig. 5a that F grows systematically when b is increasing and a is decreasing. Increasing the value of b further than 1, however, results in a very slow growth of F. Indeed, for a = 10−3, increasing b from b = 1 to b = 10, results in increasing F by only about 2%. Our results show that F reaches a maximum of ≈0.26, which is over four times higher than the maximum value of 0.06 obtained with a spatially uniform rigidity (see Fig. 2).

3.3.1 Limiting behaviour for large b and small a. In the limit of large b and low a (i.e. the bottom right corners of Fig. 5, and thus close to the unbounded limit), the optimal values for x0 appear to all converge to x0 ≈ 0.67, as shown in Fig. 5b. To understand the physical origin of this limiting value, we may consider separately the contribution of the rigid and soft parts of the filament to the propulsion.

The contribution of the distal end of the filament to the propulsion (x > x0) is captured by the term

 
image file: d4sm01388c-t33.tif(33)
in eqn (32) (where we have used the condition (29) for the second equality), and the contribution from the proximal end is then defined as F1 := FF2. In Fig. 6 we plot the values of F1 and F2 for a = 0.001 and b = 1 (shown as red dots in Fig. 5) as functions of x0. For intermediate values of x0 (away from 0 and 1), the distal part of the filament (x > x0, F2, red line) gives the dominant contribution to the propulsive force, while the almost rigid proximal part (x < x0, F1, blue line) is always much smaller. This suggests that the optimal value of x0 is determined by the dynamics of the soft (distal) part of the filament.


image file: d4sm01388c-f6.tif
Fig. 6 Dimensionless propulsive force as a function of x0 for a = 0.001 and b = 1. F1 and F2 correspond to the contribution of the proximal and the distal ends of the filament, respectively, with F1 := FF2. Fend corresponds to a simplified model, where the proximal end is considered rigid. The vertical line corresponds to the optimal location x0 ≈ 0.67.

To understand physically why the optimal value of x0 is set by the flexible end, let us consider a simplified setup and focus solely on the distal end portion, x ∈ [x0,1], in the limit where it is set into motion by an oscillating proximal end that is straight and rigid; in that case, by reversibility, the proximal end does not contribute to any propulsion. The problem in eqn (27) then reduces to,

 
image file: d4sm01388c-t34.tif(34)
The position and slope of this filament at x = x0 have to match with those of the rigid (proximal) piece of the filament (29). Thus, the position will follow the displacement at point x0, y(x0,t) = x0[thin space (1/6-em)]sin[thin space (1/6-em)]t and the slope will coincide with the actuation at the base of the rigid filament yx(x0,t) = sin[thin space (1/6-em)]t; the actuation of the flexible portion of the filament is therefore a combination of oscillating position and slope.

The solution to eqn (34) can be found analytically, as it is a problem for a filament with constant bending rigidity.22 The resulting propulsive force, termed Fend, is shown in Fig. 6 as a function of x0 for the values a = 0.001 and b = 1. The optimal value for x0 is seen to be the same for this simplified problem as for the full optimisation, demonstrating that the propulsion is indeed governed by the distal portion of the filament. In this simplified problem, the location of x0 determines both the length of the soft distal part of the filament (1 − x0) and the displacement at its base. The intermediate value of x0 is set by a balance between the requirement of having a finite-sized length for the oscillating distal part (so x0 not too close to 1) and the requirement of having a finite oscillating amplitude (so x0 not too close to 0).

4 Optimal shapes with constraint C2 (fixed filament volume)

The previous section demonstrated that, in the case of filaments whose bending rigidities are constrained to remain within a prescribed interval (constraint C1), the optimal filament shape was always piece-wise constant. We now consider the optimisation problem with the fixed-volume constraint (C2).

4.1 Optimal shapes have smoothly-decreasing bending rigidity, followed by a constant value

We use the algorithm in Appendix A.2 to solve the optimisation problem for given values of a and V0 starting with a random initial condition. The resulting optimal bending functions A for different values of V0 and fixed a = 0.01 are shown in Fig. 7a; the inset shows the time-varying shape of the filament for the optimal profile in the specific case where V0 = 0.5 (purple curve in main figure).
image file: d4sm01388c-f7.tif
Fig. 7 Optimisation of a filament with of fixed volume. (a) Numerical solution for optimal A for various values of fixed volume and a = 0.01; the colours correspond to different values of V0. Inset: Position of the filament y(x) at different time points for the optimal A for V0 = 0.5 (purple curve in the main figure). Black thick line corresponds to t = 0 and further times are shown with fading greyscale. (b) The dimensionless propulsive force F corresponding to each optimisation as a function of V0 (x-axis) and with a = 0.01; the colours are the same as in (a).

Similarly in the bounded case, the value of A is seen to be constant and equal to the lower boundary a beyond a certain location x0. However, in contrast with the previous case, closer to the base (i.e. for x < x0) the bending rigidity A is no longer constant but is always a continuously decreasing function of x.

In Fig. 7b we show the resulting propulsive force. We find that larger prescribed volumes V0 lead to systematically higher propulsive force. However, as seen in the figure, the increase slows down; from V0 = 0.4 to V0 = 0.9, the value F increases only by about 4%, while from V0 = 0.7 to V0 = 0.9 only by 0.5%.

In the case where the material has a constant Young's modulus, our result for the optimal bending rigidity means we have a solution for optimal radius of the filament, given by r = (A/)1/4. An example of the optimal r for V = 0.5 and a = 0.01 is shown with solid red line in the inset of Fig. 8.


image file: d4sm01388c-f8.tif
Fig. 8 Maximum propulsive force (colours) for A given by the empirical ansatz in eqn (35) as a function of the parameters V0 and a. Inset: The spatially variable thickness of the filament image file: d4sm01388c-t36.tif for V0 = 0.5 and a = 0.01 (i.e. we assumed = 1). The solid curve is the optimal shape obtained numerically, while the dashed curve is the approximated optimal shape from eqn (35).

4.2 Approximated optimal A

In order to further study the dependency of F on values of a and V0, we introduce the empirical approximation for the optimal A as
 
image file: d4sm01388c-t35.tif(35)
where the constant a0 is chosen to satisfy the constant volume condition and where the transition point x0 is fixed. The radius of the filament, given by r = (A/)1/4 with A from the model in eqn (35) is shown in the inset of Fig. 8 in black dashed line for the choices V0 = 0.5, a = 0.01 and = 1. The relative error for this approximation compared to the optimal case (solid red line), is below 1% in terms of the propulsive force F.

Using this model, we may vary the values of both a and V0 and compute the maximum value for F over possible values of x0, as shown in Fig. 8. The results are very similar to those obtained with the optimisation under the constraint C1. Indeed, the propulsive force is systematically larger if a is small and if V0 is large, as was seen in Fig. 5a, while the values of the optimal transition point x0, after which the function A, is a constant are very similar to Fig. 5b (not shown).

The physical interpretation of the results in this section is very similar to the one we proposed in the previous section and for both constraints A is stiff at the base and soft at the distal end with x0 governing the length of the soft distal end.

5 Discussion

5.1 Summary

Enhancing the propulsion of elastic filaments through optimised design is a significant challenge within the realm of micro-robotics.8 One promising approach involves the optimisation of the bending rigidity along the length of the filament, here denoted by A in its dimensionless version. In this study, we focus on the optimisation of A with the aim of maximising the propulsive force generated by a model system: an elastic slender filament with an oscillating slope at its base. To accomplish this, we employ techniques of functional optimisation, first by computing the functional derivative of A using an adjoint problem, and then by using this derivative within a steepest ascent algorithm operating in the functional space. Our optimisation process considers two sets of additional constraints: ensuring that A remains within defined bounds, or maintaining a fixed volume for the filament.

Our main finding demonstrates that, in all tested cases, the optimal filament must be stiff at the base and soft at the distal end for optimal propulsion, which we were able to interpret physically. In the case where A is constrained to remain bounded between prescribed values a and b, the optimal solution is invariably piecewise constant: A assumes the upper bound near the filament's base and the lower bound at the distal end. When b is substantially larger than a the propulsive force increases by more than a factor of 4 above the optimal value for a spatially uniform A. The transition point, x0, where the rigidity abruptly goes from the upper bound to the lower bound is located about 2/3 along the filament's length, which we demonstrated is governed by the propulsive physics at the distal end. For the fixed volume constraint, the optimal filament continues to have higher stiffness at the base and becomes softer towards its distal end, but in this case the shape decreases smoothly near the base.

5.2 Practical considerations

The two constraints we considered indicate possible practical ways of designing an optimal elastic filament. The bending rigidity of a filament, A, is determined by the product of the Young's modulus of the material and the moment area of inertia, proportional to r4, where r is the filament's radius. Changes in A may thus be obtained either by changing the Young's modulus or by modifying the filament radius. In the case of bounded A, the optimal piecewise-constant rigidity could be achieved by merging two segments with different Young's moduli. For the fixed volume case, on the other hand, the Young's modulus is fixed and spatially variable radius of the filament is given by r = (A/)1/4, shown for example in the inset of Fig. 8.

To further illustrate the advantage of using a filament with varying bending rigidity, consider the (mm)-scale PDMS filament from Williams et al.24 with length L = 1.5 mm, radius r = 5 μm and Young's modulus of 3.86 MPa. Using the viscosity μ = 1.15 mPa s24 and frequency ω/2π = 3.6 Hz, we obtain the estimate A ≈ 0.03. Small oscillations with angle ε = 0.4 rad of this filament yields a propulsive force of [F with combining macron] ≈ 0.3 nN. If instead the distal part of the bending rigidity is reduced by one order of magnitude compared to the proximal part, the force increases by almost one order of magnitude to about [F with combining macron] ≈ 1.7 nN. Similarly, for a PPy filament with L = 9 μm, r = 100 nm,17 a Young's modulus of 100 MPa56 and oscillating with frequency ω/2π = 20 Hz,17 we estimate A ≈ 1 in 65% glycerol solution (μ = 15 mPa s57). Under an oscillation slope of ε = 0.4 rad, we obtained a small propulsive force of [F with combining macron] ≈ 0.22 pN; however, lowering the bending rigidity of the distal part by two orders of magnitude increases the force by a factor of about 18, to [F with combining macron] ≈ 4 pN.

5.3 Comparison with past work

Our results demonstrate rigorously the advantages of having a stiff base and a soft distal end for enhancing propulsion, compared to a uniform filament with optimal constant rigidity, which is consistent with past work.47–50 Singh and Yadava48 reported a 22% increase in propulsion for a linearly tapered elastic filament whose slope at the base was oscillating, while Peng et al.50 found up to 12.5% increase in propulsion using an exponential taper for a clamped oscillating filament. In contrast, the optimal bending rigidity reported in our work (also for imposed oscillating slope at the base) results in a more substantial (four-fold) increase in propulsion.

In Peng et al.50 the authors studied the optimal propulsion of a filament with two segments: one rigid and one soft. They showed that for the case of a clamped oscillating filament, the maximum thrust is achieved when the base is rigid and the distal end is soft, with a transition point of 0.82 of the total length. In contrast, for a hinged filament, the propulsion is maximised if the base of the filament is soft and the distal end is stiff, highlighting the importance of the actuation mechanism in determining the optimal bending rigidity.

5.4 Biological relevance

We can draw parallels between our findings and the structure of a human spermatozoon. The flagellum of a spermatozoon cell is tapered,6 with an estimated proximal-to-distal ratio in bending rigidity of 35,58 and after about 0.65 of the length of the filament, the distal end maintains a constant level of stiffness. Moreover, a sperm flagellum is forced by a continuous internal actuation via molecular motors distributed along the whole flagellum54 except for the passive soft region at the distal end (the so-called ‘end piece’). In the study by Neal et al.58 the authors noted that this passive end piece can result in a remarkable increase in swimming speed, up to 72%, and a substantial boost in hydrodynamic efficiency, up to 438%. This result is consistent with our findings, where the presence of a soft distal end is crucial for enhancing the propulsion of the filament. We note, however, that the comparison with biological cells is limited, owing to the qualitatively different actuation mechanisms.

5.5 Outlook

The methodology used in this study could be extended to tackle optimisation under other actuation mechanisms. This could allow comparison with existing artificial swimmers employing rotating51 or oscillating magnetic heads.17 When applied to the case of a clamped filament with an oscillating base, we found much more complex optimal shapes than for a hinged filament, hinting at a very rich optimal landscape.

Additionally, our study focused on maximising solely the propulsive force, and it will be valuable to explore other objectives such as maximising swimming speed or hydrodynamic efficiency. In our simplified model, maximising swimming speed for a swimmer with a passive filament would lead to the same optimal solution. Indeed, in that case the steady horizontal constant velocity V of the swimmer is obtained by balancing drag and propulsion as image file: d4sm01388c-t37.tif, where Dh is drag coefficient for the head.48,54 Since, at small amplitude, only the propulsive force depends on filament dynamics, the optimisation for V is identical to that of [F with combining macron].

Our study exhibits certain limitations that could be subject to future research. Firstly, our assumption of linearised elastohydrodynamics means that our model is not suitable for describing large displacements of the filament. Secondly, we rely on resistive-force theory, which simplifies our analysis by neglecting hydrodynamic interactions between different parts of the filament. While this is a reasonable approximation for small amplitude deformations, it may not hold true when the filament undergoes significant curvature. Furthermore, the model does not account for the change in drag coefficients with r, which appear in the force balance eqn (4) and in the expression of the propulsive force (6). In the present work, we looked at changes in bending stiffness up to 3–4 orders of magnitude, which would correspond to changes in r by about one order of magnitude. This will result in changing drag coefficients by about a factor of 2, which may alter the precise quantitative nature of the optimal solution.

Data availability

All results and figures in the article were generated using Matlab R2024a. The Matlab codes can be found at the following DOI: 10.5281/zenodo.15111675.

Conflicts of interest

There are no conflicts to declare.

Appendix

A.1 Algorithm for constraint C1 (bounded A)

For a bounded A(x), we have to make sure that at each optimisation step, it remains within the chosen domain [a,b]. To ensure this, we introduce a projection operator P(A) = max(a,min(A,b)) and we adapt the projected gradient method, which consists of the following steps:55

• Choose the initial condition A1, aA1b.

• For k = 1,2,3,…, repeat the following steps until convergence criteria are achieved

(1) Solve for filament dynamics, eqn (11)–(13), and the adjoint problem, eqn (23) and (24), with given Ak. We use Matlab solver bvp4c to solve the ODE system;

(2) Set sk = dF/dA|Ak, where the functional derivative is computed using eqn (26);

(3) Choose step σk by projected Armijo rule (see below) such that F(P(Ak + σksk)) > F(Ak);

(4) Set Ak+1 := P(Ak + σksk).

Defining the errors

 
err1 = |F(Ak) − F(Ak−1)|/F(Ak−1), err2 = max|AkAk−1|,(36)
the convergence criteria are picked to be err1 < δ1 and err2 < δ2, where and δ1 and δ2 are prescribed tolerances.

To choose the optimisation step σk, we use the projected Armijo rule and choose maximum σk ∈ {1,1/2,1/4,…} for which

 
image file: d4sm01388c-t38.tif(37)
where ‖·‖ refers to the L2 norm. As shown in classical work, this algorithm converges.55

A.2 Algorithm for constraint C2 (fixed volume)

In the case of constraint (C2), we use a similar algorithm to the one reported in Section 2.4 but modify the objective function to account for constant volume constraint,
 
image file: d4sm01388c-t39.tif(38)
Here, the second term ρ is a penalty parameter for the volume constraint image file: d4sm01388c-t40.tif. The optimisation algorithm is then similar to the one in Section A.1, replacing F with FV and considering the projection function only for lower boundary, P(A) = max(a,A). Since we do not want the final solution to depend on the value of ρ, we set the penalty parameter ρ to be an increasing sequence ρm = ρ02m, where ρ0 is an initial choice. For each value of ρm we solve the minimisation problem reported in Section A.1 until convergence is reached for given tolerances δm1 and δm2. For convenience, we work in terms of the radius r instead of A, and thus the errors at each step k are err1 = |F(Ak) − F(Ak−1)|/F(Ak−1) and err2 = max|rkrk−1|. We then consider the next value of ρm and follow the procedure again, until the error err3 = |F(Am) − F(Am−1)|/F(Am−1) + max|rmrm−1| gets smaller than a set tolerance δ3 or a maximum number of iterations is reached. In our solution we set min[thin space (1/6-em)]δm1 = 10−4, min[thin space (1/6-em)]δm2 = 10−4 and δ3 = 5 × 10−5.

Acknowledgements

The authors would also like to thank M. Tătulea-Codrean for useful discussions and feedback. This project has received funding from the European Research Council under the European Union's Horizon 2020 Research and Innovation Programme (Grant No. 682754 to E. L.).

References

  1. L. J. Fauci and R. Dillon, Annu. Rev. Fluid Mech., 2006, 38, 371–394 CrossRef.
  2. D. L. Koch and G. Subramanian, Annu. Rev. Fluid Mech., 2011, 43, 637–659 CrossRef.
  3. R. E. Goldstein, Annu. Rev. Fluid Mech., 2015, 47, 343 CrossRef PubMed.
  4. N. Wadhwa and H. C. Berg, Nat. Rev. Microbiol., 2022, 20, 161–173 CrossRef CAS PubMed.
  5. E. Lauga and T. R. Powers, Rep. Prog. Phys., 2009, 72, 096601 CrossRef.
  6. E. A. Gaffney, H. Gadêlha, D. J. Smith, J. R. Blake and J. C. Kirkman-Brown, Annu. Rev. Fluid Mech., 2011, 43, 501–528 CrossRef.
  7. E. Lauga, Annu. Rev. Fluid Mech., 2016, 48, 105–130 CrossRef.
  8. B. J. Nelson, I. K. Kaliakatsos and J. J. Abbott, Annu. Rev. Biomed. Eng., 2010, 12, 55–85 CrossRef CAS PubMed.
  9. K. E. Peyer, A. W. Mahoney, L. Zhang, J. J. Abbott and B. J. Nelson, Microbiorobotics, 2012, 165–199 Search PubMed.
  10. K. E. Peyer, L. Zhang and B. J. Nelson, Nanoscale, 2013, 5, 1259–1272 Search PubMed.
  11. J. Li, B. E.-F. de Ávila, W. Gao, L. Zhang and J. Wang, Sci. Rob., 2017, 2, eaam6431 CrossRef PubMed.
  12. X.-Z. Chen, M. Hoop, F. Mushtaq, E. Siringil, C. Hu, B. J. Nelson and S. Pané, Appl. Mater. Today, 2017, 9, 37–48 CrossRef.
  13. H. Zhou, C. C. Mayorga-Martinez, S. Pané, L. Zhang and M. Pumera, Chem. Rev., 2021, 121, 4999–5041 Search PubMed.
  14. S. Schuerle, S. Pané, E. Pellicer, J. Sort, M. D. Baró and B. J. Nelson, Small, 2012, 8, 1498–1502 Search PubMed.
  15. R. Dreyfus, J. Baudry, M. L. Roper, M. Fermigier, H. A. Stone and J. Bibette, Nature, 2005, 437, 862–865 CrossRef CAS PubMed.
  16. H. C. M. Sun, P. Liao, T. Wei, L. Zhang and D. Sun, Micromachines, 2020, 11, 404 CrossRef PubMed.
  17. B. Jang, E. Gutman, N. Stucki, B. F. Seitz, P. D. Wendel-Garca, T. Newton, J. Pokki, O. Ergeneman, S. Pané and Y. Or, et al. , Nano Lett., 2015, 15, 4829–4833 CrossRef CAS PubMed.
  18. L. Zhang, J. J. Abbott, L. Dong, B. E. Kratochvil, D. Bell and B. J. Nelson, Appl. Phys. Lett., 2009, 94, 064107 CrossRef.
  19. E. Lauga, The fluid dynamics of cell motility, Cambridge University Press, 2020, vol. 62 Search PubMed.
  20. E. M. Purcell, Am. J. Phys., 1977, 45, 3–11 CrossRef.
  21. K. Machin, J. Exp. Biol., 1958, 35, 796–806 CrossRef.
  22. C. H. Wiggins and R. E. Goldstein, Phys. Rev. Lett., 1998, 80, 3879 CrossRef CAS.
  23. S. Guo, Q. Pan and M. B. Khamesee, Microsyst. Technol., 2008, 14, 307–314 CrossRef.
  24. B. J. Williams, S. V. Anand, J. Rajagopalan and M. T. A. Saif, Nat. Commun., 2014, 5, 3081 CrossRef PubMed.
  25. V. Magdanz, I. S. Khalil, J. Simmchen, G. P. Furtado, S. Mohanty, J. Gebauer, H. Xu, A. Klingner, A. Aziz and M. Medina-Sánchez, et al. , Sci. Adv., 2020, 6, eaba5855 CrossRef CAS PubMed.
  26. V. Magdanz, A. Klingner, L. Abelmann and I. S. Khalil, J. Micro Bio Rob., 2022, 18, 49–60 CrossRef.
  27. N. Celi, D. Gong and J. Cai, Sci. Rep., 2021, 11, 21728 CrossRef CAS PubMed.
  28. F. Mushtaq, H. Torlakcik, M. Hoop, B. Jang, F. Carlson, T. Grunow, N. Läubli, A. Ferreira, X.-Z. Chen and B. J. Nelson, et al. , Adv. Funct. Mater., 2019, 29, 1808135 CrossRef.
  29. C. H. Wiggins, D. Riveline, A. Ott and R. E. Goldstein, Biophys. J., 1998, 74, 1043–1060 CrossRef CAS PubMed.
  30. T. S. Yu, E. Lauga and A. Hosoi, Phys. Fluids, 2006, 18, 091701 CrossRef.
  31. T.-H. Wu, R.-S. Guo, G.-W. He, Y.-M. Liu and D. Qi, J. Theor. Biol., 2014, 349, 1–11 CrossRef PubMed.
  32. C. Moreau, L. Giraldi and H. Gadêlha, J. R.Soc., Interface, 2018, 15, 20180235 Search PubMed.
  33. M. Jabbarzadeh and H. C. Fu, J. Comput. Phys., 2020, 418, 109643 CrossRef.
  34. A. Hilfinger, A. K. Chattopadhyay and F. Jülicher, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2009, 79, 051918 CrossRef PubMed.
  35. D. Oriola, H. Gadêlha and J. Casademunt, R. Soc. Open Sci., 2017, 4, 160698 CrossRef PubMed.
  36. D. Tam and A. Hosoi, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2011, 83, 045303 CrossRef PubMed.
  37. O. Pironneau and D. Katz, J. Fluid Mech., 1974, 66, 391–415 Search PubMed.
  38. O. Pironneau, D. Katz, T. Wu, C. Brokaw and C. Brennen, Swimming and Flying in Nature, 1975, vol. 1, pp. 161–172 Search PubMed.
  39. E. Lauga, Phys. Rev. Fluids, 2020, 5, 123101 CrossRef.
  40. S. J. Lighthill, Mathematical biofluiddynamics, SIAM, 1975 Search PubMed.
  41. S. E. Spagnolie and E. Lauga, Phys. Fluids, 2010, 22, 455 Search PubMed.
  42. E. Lauga and C. Eloy, J. Fluid Mech., 2013, 730, R1 CrossRef.
  43. M. F. Velho Rodrigues, M. Lisicki and E. Lauga, PLoS One, 2021, 16, e0252291 CrossRef CAS PubMed.
  44. M. Lisicki, M. F. V. Rodrigues and E. Lauga, J. Fluid Mech., 2024, 978, R1 CrossRef.
  45. E. Lauga, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2007, 75, 041916 CrossRef PubMed.
  46. H. Gadêlha, Regul. Chaotic Dyn., 2013, 18, 75–84 CrossRef.
  47. J. S. Rathore, R. Majumdar and N. N. Sharma, IEEE Trans. Nanotechnol., 2012, 11, 1117–1121 Search PubMed.
  48. T. S. Singh and R. Yadava, Biomed. Phys. Eng. Express, 2017, 4, 015019 CrossRef.
  49. R. Kotesa, J. Rathore and N. Sharma, Bionanoscience, 2013, 3, 343–347 CrossRef.
  50. Z. Peng, G. J. Elfring and O. S. Pak, Soft Matter, 2017, 13, 2339–2347 RSC.
  51. A. M. Maier, C. Weig, P. Oswald, E. Frey, P. Fischer and T. Liedl, Nano Lett., 2016, 16, 906–910 CrossRef CAS PubMed.
  52. R. Cox, J. Fluid Mech., 1970, 44, 791–810 CrossRef.
  53. A. Öchsner, Classical beam theories of structural mechanics, Springer, 2021, vol. 42 Search PubMed.
  54. O. S. Pak and E. Lauga, Fluid–Structure Interactions in Low-Reynolds-Number Flows, The Royal Society of Chemistry, 2015, pp. 100–167 Search PubMed.
  55. M. Hinze, R. Pinnau, M. Ulbrich and S. Ulbrich, Optimization with PDE constraints, Springer Science & Business Media, 2008, vol. 23 Search PubMed.
  56. J. D. Madden, Electroactive Polymers for Robotic Applications: Artificial Muscles and Sensors, Springer, 2007, ch. 5, pp. 121–152 Search PubMed.
  57. J. B. Segur and H. E. Oberstar, Ind. Eng. Chem., 1951, 43, 2117–2120 CrossRef CAS.
  58. C. V. Neal, A. L. Hall-McNair, J. Kirkman-Brown, D. J. Smith and M. T. Gallagher, Phys. Rev. Fluids, 2020, 5, 073101 CrossRef.

This journal is © The Royal Society of Chemistry 2025
Click here to see how this site uses Cookies. View our privacy policy here.