## Abstract

Soft slender structures are ubiquitous in natural and artificial systems, in active and passive settings and across scales, from polymers and flagella, to snakes and space tethers. In this paper, we demonstrate the use of a simple and practical numerical implementation based on the Cosserat rod model to simulate the dynamics of filaments that can bend, twist, stretch and shear while interacting with complex environments via muscular activity, surface contact, friction and hydrodynamics. We validate our simulations by solving a number of forward problems involving the mechanics of passive filaments and comparing them with known analytical results, and extend them to study instabilities in stretched and twisted filaments that form solenoidal and plectonemic structures. We then study active filaments such as snakes and other slender organisms by solving inverse problems to identify optimal gaits for limbless locomotion on solid surfaces and in bulk liquids.

### 1. Introduction

Quasi one-dimensional objects are characterized by having one dimension, the length *L*, much larger than the others, say the radius *r*, so that *L*/*r*≫1. Relative to three-dimensional objects, this measure of slenderness allows for significant mathematical simplification in accurately describing the physical dynamics of strings, filaments and rods. It is thus perhaps not surprising that the physics of strings has been the subject of intense study for centuries [1–10], and indeed their investigation substantially predates the birth of three-dimensional elasticity.

Following the pioneering work of Galileo on the bending of cantilevers, one-dimensional analytical models of beams date back to 1761 when Jakob Bernoulli first introduced the use of differential equations to capture the relationship between geometry and bending resistance in a *planar elastica*, that is, an elastic curve deforming in a two-dimensional space. This attempt was then progressively refined by Huygens *et al.* [11], until Euler presented a full solution of the planar elastica, obtained by minimizing the strain energy and by recognizing the relationships between flexural stiffness and material and geometric properties. Euler also showed the existence of bifurcating solutions in a rod subject to compression, identifying its first buckling mode, while Lagrange formalized the corresponding multi-modal solution [5]. Non-planar deformations of the elastica were first tackled by Kirchhoff [1,6] and Clebsch [2] who envisaged a rod as an assembly of short undeformable straight segments with dynamics determined by contact forces and moments, leading to three-dimensional configurations. Later, Love [3] approached the problem from a Lagrangian perspective characterizing a filament by contiguous cross sections that can rotate relative to each other, but remain undeformed and perpendicular to the centre line of the rod at all times; in modern parlance this assumption is associated with dynamics on the rotation group SO(3) at every cross section. The corresponding equations of motion capture the ability of the filament to bend and twist, but not shear or stretch. Eventually, the Cosserat brothers [4] relaxed the assumption of inextensibility and cross-section orthogonality to the centre line, deriving a general mathematical framework that accommodates all six possible degrees of freedom associated with bending, twisting, stretching and shearing, effectively formulating dynamics on the full Euclidean group SE(3).

These mathematical foundations [5] prompted a number of discrete computational models [12–16] that allow for the exploration of a range of physical phenomena. These include, for example, the study of polymers and DNA [12,17], elastic ribbons and filaments [14,18,19], botanical applications [20,21], woven cloth [22], and tangled hair and fibres [19,23–25]. Because the scaled ratio of the stretching and shearing stiffness to the bending stiffness for slender filaments is *L*^{2}/*r*^{2}≫1, the assumption of inextensibility and/or unshearability is usually appropriate, justifying the widespread use of the Kirchhoff model in the aforementioned applications.

Fewer studies have considered, in different flavours, the Cosserat model, mostly to take advantage of relaxed extension and shearing constraints to simplify numerical implementations and to facilitate the integration of additional physical effects. For example, specialized models including extension and constitutive laws based on strain rates have been developed for the investigation of viscous threads [15,16,26–30]. Lim and Peskin allow for small numerical shear and axial strain and couple their model with an accurate viscous flow solver to investigate fluid–mechanic interactions of ribbons and flagella [31–36]. The graphics software Corde [13] implements a Cosserat-based fast solver for the rendering of looping systems, accounts for (self) contact, operates in the small extension regime and maintains the unshearability constraint, showcasing a number of visually impressive and physically plausible scenarios. Durville [25] introduced a fibre model specialized to efficiently resolve contact-friction effects in entangled materials. Linn [37] explored an elegant theoretical connection between Cosserat rods and the differential geometry of framed curves, and numerically tested it on the classic Euler’s Elastica and Kirchhoff’s helix problems [38]. Finally, Sonneville [39] presents a geometrically exact finite element formulation on the Euclidean group SE(3) and verifies it on test cases that do not involve stretching or environmental effects.

More recently, there has been a need to generalize the model to explain new experimental phenomena such as the plectoneme–solenoid transition [40], that has been used as the basis for artificial muscles [41–43]. Additional new technologies such as soft robotics [44,45] are also generating applications for highly stretchable and shearable elastomeric structures raising the need for methods able to realistically handle these large strains together with a variety of interface and environmental effects. Moreover, the capability to computationally solve both forward [46–49] and inverse problems is emerging as a crucial paradigm to aid the design of novel, more capable soft-robotic prototypes [50,51], as well as to characterize from an optimality standpoint biophysical phenomena [52–57]. Motivated by these advancements and challenges, we use a versatile implementation of the Cosserat model that we validate in a set of physical simulations, and then deploy in the context of inverse design problems to broaden the spectrum of its potential engineering and biophysical investigations. By taking advantage of the Cosserat formalism, and consistently with the full Euclidean group SE(3), we allow for bend, twist and significant shear and stretch [4], and demonstrate the importance of the latter through an example inspired by artificial muscle actuation [41,42] in which the transition between plectonemes and solenoids is crucially enabled by axial extension. Then, moving beyond the passive mechanics of individual filaments, we account for the interaction between filaments and complex environments with a number of additional biological and physical features, including muscular activity, self-contact and contact with solid boundaries, isotropic and anisotropic surface friction and viscous interaction with a fluid. Finally, we demonstrate the capabilities and the robustness of our solver by embedding it in an inverse design cycle for the identification of optimal terrestrial and aquatic limbless locomotion strategies.

The paper is structured as follows. In §2, we review and introduce the mathematical foundations of the model. In §3, we present the corresponding discrete scheme and validate our implementation against a battery of benchmark problems. In §4, we detail the physical and biological enhancements to the original model, and finally in §5, we showcase the potential of our solver via the study of solenoids and plectonemes as well as limbless biolocomotion. Mathematical derivations and additional validation test cases are presented in the appendix.

### 2. Governing equations

We consider filaments as slender cylindrical structures deforming in three dimensions with a characteristic length *L* which is assumed to be much larger than the radius (*L*≫*r*) at any cross section. Then the filament can be geometrically reduced to a one-dimensional representation, and its dynamical behaviour may be approximated by averaging all balance laws at every cross section [5]. We start with a description of the commonly used Kirchhoff–Love theory that accounts for bend and twist at every cross section but ignores stretch and shear, before moving on to the Cosserat theory that also accounts for these additional degrees of freedom.

#### 2.1. Kirchhoff–Love theory for inextensible, unshearable rods

As illustrated in figure 1, a filament in the Cosserat rod theory can be described by a centre line $\mathbf{\text{r}}:(s\in [0,L]\in \mathbb{R},t\in {\mathbb{R}}^{+})\to {\mathbb{R}}^{3}$ and an oriented frame of reference $\mathbf{\text{Q}}:(s\in [0,L]\in \mathbb{R},t\in {\mathbb{R}}^{+})\to \text{SO}(3)$ equivalent to the orthonormal triad of unit vectors **Q**={**d**_{1},**d**_{2},**d**_{3}}. Here, *s* is the centre line arc-length coordinate in its current configuration and *t* is the time.

Denoting by **x** any generic vector represented in the Eulerian frame and ${\mathbf{\text{x}}}_{\mathcal{L}}$ as the body-convected (Lagrangian) frame of reference allows us to write

**x**in the laboratory canonical basis {

**i**,

**j**,

**k**}, while equation (2.2) expresses the same vector in the body-convected director basis {

**d**

_{1},

**d**

_{2},

**d**

_{3}}. Then, the matrix

**Q**transforms any vector

**x**from the laboratory to the body-convected representation via ${\mathbf{\text{x}}}_{\mathcal{L}}=\mathbf{\text{Q}}\mathbf{\text{x}}$ and, conversely, $\mathbf{\text{x}}={\mathbf{\text{Q}}}^{-1}{\mathbf{\text{x}}}_{\mathcal{L}}={\mathbf{\text{Q}}}^{\mathrm{T}}{\mathbf{\text{x}}}_{\mathcal{L}}$, because

**Q**

^{T}

**Q**=

**Q**

**Q**

^{T}=1. In general, we need to distinguish between the arc-length coordinate

*s*that corresponds to the current filament configuration and the arc-length coordinate $\hat{s}$ associated with the reference configuration of the filament due to stretching (figure 1—throughout we will use hatted quantities to denote the reference configuration).

We will first start by presenting the equations of motion under the assumption of inextensibility (i.e. $s=\hat{s}$, Kirchhoff model), before generalizing them to the stretchable case in the subsequent sections. Denoting the rod angular velocity as ** ω**=

*vec*[(∂

**Q**/∂

*t*)

^{T}

**Q**] and the generalized curvature as

**=**

*κ**vec*[(∂

**Q**/∂

*s*)

^{T}

**Q**], where

*vec*[

**A**] denotes the 3-vector associated with the skew-symmetric matrix

**A**, the following transport identities hold

*ρ*is the constant material density,

*A*is the cross-sectional area,

**v**is the velocity, ${\mathbf{\text{n}}}_{\mathcal{L}}$ and ${\mathit{\tau}}_{\mathcal{L}}$ are, respectively, the internal force and couple resultants,

**f**and

**c**are external body force and torque line densities, and the tensor

**I**is the second area moment of inertia (throughout this study we assume circular cross sections; see the appendix).

To close the above system of equations (2.4)–(2.7) and determine the dynamics of the rod, it is necessary to specify the form of the internal forces and torques generated in response to bend and twist, corresponding to the 3 d.f. at every cross section. The strains are defined as the relative local deformations of the rod with respect to its natural strain-free reference configuration. Bending and twisting strains are associated with the spatial derivatives of the material frame director field {**d**_{1},**d**_{2},**d**_{3}} and are characterized by the generalized curvature. Specifically, the components of the curvature projected along the directors (${\mathit{\kappa}}_{\mathcal{L}}={\kappa}_{1}{\mathbf{\text{d}}}_{1}+{\kappa}_{2}{\mathbf{\text{d}}}_{2}+{\kappa}_{3}{\mathbf{\text{d}}}_{3}$) coincide with bending (*κ*_{1},*κ*_{2}) and twist (*κ*_{3}) strains in the material frame (table 1).

deformation modes | strains | rigidities | loads |
---|---|---|---|

bending about d_{1} |
κ_{1} |
B_{1}=EI_{1} |
${\tau}_{1}={B}_{1}({\kappa}_{1}-{\kappa}_{1}^{o})$ |

bending about d_{2} |
κ_{2} |
B_{2}=EI_{2} |
${\tau}_{2}={B}_{2}({\kappa}_{2}-{\kappa}_{2}^{o})$ |

twist about d_{3} |
κ_{3} |
B_{3}=GI_{3} |
${\tau}_{3}={B}_{3}({\kappa}_{3}-{\kappa}_{3}^{o})$ |

shear along d_{1} |
σ_{1} |
S_{1}=α_{c}GA |
${n}_{1}={S}_{1}({\sigma}_{1}-{\sigma}_{1}^{o})$ |

shear along d_{2} |
σ_{2} |
S_{2}=α_{c}GA |
${n}_{2}={S}_{2}({\sigma}_{2}-{\sigma}_{2}^{o})$ |

stretch along d_{3} |
σ_{3} |
S_{3}=EA |
${n}_{3}={S}_{3}({\sigma}_{3}-{\sigma}_{3}^{o})$ |

We assume a perfectly elastic material so that the *stress–strain* relations are linear. Integration of the torque densities over the cross-sectional area *A* yields the bending and twist rigidities (table 1), so that the resultant *torque–curvature* relations can be generically expressed in vectorial notation as

*B*

_{1}the flexural rigidity about

**d**

_{1},

*B*

_{2}the flexural rigidity about

**d**

_{2}and

*B*

_{3}the twist rigidity about

**d**

_{3}. Here, the vector ${\mathit{\kappa}}_{\mathcal{L}}^{o}$ characterizes the intrinsic curvatures of a filament that in its stress-free state is not straight. We wish to emphasize here that the constitutive laws are most simply expressed in a local Lagrangian form; hence the use of ${\mathit{\kappa}}_{\mathcal{L}}$ and not

**.**

*κ*The Kirchhoff rod is defined by the additional assumption that there is no axial extension or compression or shear strain. Then the arc-length *s* coincides with $\hat{s}$ at all times, and the tangent to the centre line is also normal to the cross section, so that **t**=**d**_{3} [5]. This implies that ${\mathbf{\text{n}}}_{\mathcal{L}}$ serves as a Lagrange multiplier, and that the *torque-curvature* relations of equation (2.8) are linear.

This completes the formulation of the equations of motion for the Kirchhoff rod, and when combined with boundary conditions suffices to have a well-posed initial boundary value problem. For the general stretchable and shearable case, all geometric quantities (A, **I**, ${\mathit{\kappa}}_{\mathcal{L}}$, etc.) must be rescaled appropriately, as addressed in the following sections.

#### 2.2. Cosserat theory of stretchable and shearable filaments

In the general case of soft filaments, at every cross section we also wish to capture transverse shear and axial strains in addition to bending and twisting. As we wish to account for all six deformation modes associated with the 6 d.f. at each cross section along the rod, we must augment the Kirchhoff description of the previous section and add three more constitutive laws to define the local stress resultants ${\mathbf{\text{n}}}_{\mathcal{L}}$. In fact, they are no longer defined as Lagrange multipliers that enforce the condition of inextensibility and unshearability, i.e. we now must allow **t**≠**d**_{3}. We note that this not only enables the model of equations (2.4)–(2.7) to capture a richer set of physical phenomena, but also significantly simplifies its numerical treatment (Lagrange multipliers no longer needed), rendering this model both flexible and relatively straightforward to implement.

The shear and axial strains are associated with the deviations between the unit vector perpendicular to the cross section and the tangent to the centre line, and thus may be expressed in terms of the derivatives of the centre line coordinate **r**. In the material frame of reference, we characterize these strains by the vector ** σ** (figure 1) which then takes the form

**t**is the unit tangent vector.

Whenever the filament undergoes axial stretching or compression the corresponding infinitesimal elements deform and all related geometric quantities are affected. By assuming that the material is incompressible and that the cross sections retain their circular shapes at all times, we can conveniently express the governing equations with respect to the rest reference configuration of the filament (denoted by a hat) in terms of the local dilatation $e(\hat{s},t)$. Then, the following relations hold:

As with the Kirchhoff rod, assuming a linear material constitutive law implies linear *stress–strain* relations. Integration of the stress and couple densities over the cross-sectional area *A* yields both the rigidities associated with axial extension and shear (table 1), so that the resultant *load–strain* relations can be generically expressed in vectorial notation as

*S*

_{1},

*S*

_{2},

*S*

_{3}) is the shear/stretch stiffness matrix, with

*S*

_{1}the shearing rigidity along

**d**

_{1},

*S*

_{2}the shearing rigidity along

**d**

_{2}and

*S*

_{3}the axial rigidity along

**d**

_{3}. Here, as with the Kirchhoff rod, the vector ${\mathit{\sigma}}_{\mathcal{L}}^{o}$ corresponds to the intrinsic shear and stretch, and must be accounted for in the case of stress-free, non-trivial shapes. Although the intrinsic strains ${\mathit{\kappa}}_{\mathcal{L}}^{o}$, ${\mathit{\sigma}}_{\mathcal{L}}^{o}$ are implemented in our solver to account for pre-strained configurations, to simplify the notation in the remaining text, we will assume that the filament is intrinsically straight in a stress-free state, so that ${\mathit{\sigma}}_{\mathcal{L}}^{o}={\mathit{\kappa}}_{\mathcal{L}}^{o}=\mathbf{\text{0}}$.

The rigidities associated with bending, twisting, stretching and shearing are specified in table 1, and can be expressed as the product of a material component, represented by the Young’s (*E*) and shear (*G*) moduli, and a geometric component represented by *A*, **I** and the constant ${\alpha}_{c}=\frac{4}{3}$ for circular cross sections [58]. We note that the rigidity matrices **B** and **S** are assumed to be diagonal throughout this study, although off-diagonal entries can be easily accommodated to model anisotropic materials such as composite elements. In general, this mathematical formulation can be extended to tackle a richer set of physical problems including viscous threads [15,16], magnetic filaments [59], etc., by simply modifying the entries of **B** and **S** and introducing time-dependent constitutive laws wherein ${\mathit{\tau}}_{\mathcal{L}}({\mathit{\kappa}}_{\mathcal{L}},{\mathrm{\partial}}_{t}{\mathit{\kappa}}_{\mathcal{L}})$ and ${\mathbf{\text{n}}}_{\mathcal{L}}({\mathit{\sigma}}_{\mathcal{L}},{\mathrm{\partial}}_{t}{\mathit{\sigma}}_{\mathcal{L}})$, as for example in [15]. We also emphasize here that in the case of stretchable rods, *A* and **I** are no longer constant, rendering the *load–strain* relations nonlinear (equations (2.8) and (2.11)), even though the *stress–strain* relations remain linear. This is consistent with the modelling of hyperelastic materials such as rubbers, silicones and biological tissues and therefore in line with targeted soft robotic applications. Indeed, the combination of linear stress–strain constitutive models with the geometrical rescaling by *e* leads to a reasonable approximation of Neo-Hookean [60] and Mooney–Rivlin models [61,62] (especially tailored to hyperelastic solids) over a compression/extension range up to 30%. See the appendix for a quantitative comparison in the axial stretch case.

Having generalized the constitutive relations to account for filament stretchiness and shearability, we now generalize the equations of motion for this case. Multiplying both sides of equations (2.6) and (2.7) by *ds* and substituting the above identities together with the constitutive laws of equation (2.11) into equations (2.4)–(2.7) yields the final system

**f**and ${\mathbf{\text{c}}}_{\mathcal{L}}$ the force and torque line densities, respectively) so as to account for the dependence on

*e*.

Combined with some initial and boundary conditions, equations (2.12)–(2.15) express the dynamics and kinematics of the Cosserat rod with respect to its initial rest configuration, in a form suitable to be discretized as described in §3.

### 3. Numerical method

To derive the numerical method for the time evolution of a filament in analogy with the continuum model of §2, we first recall a few useful definitions for effectively implementing rotations. We then present the spatially discretized model of the rod, and the time discretization approach employed to evolve the governing equations.

#### 3.1. Rotations

Bending and twisting deformations of a filament involve rotations of its material frame **Q** in space and time. To numerically simulate the rod, it is critical to represent and efficiently compute these nonlinear geometric transformations quickly and accurately. A convenient way to express rotations in space or time is the matrix exponential [63–65]. Assuming that the matrix **R** denotes the rotation by the angle *θ* about the unit vector axis **u**, this rotation can be expressed through the exponential matrix $\mathbf{\text{R}}={\mathrm{e}}^{\theta \mathbf{\text{u}}}:{\mathbb{R}}^{3}\to {\mathbb{R}}^{3\times 3}$, and efficiently computed via the Rodrigues formula [66]

**u**

*vice versa*${[\cdot ]}_{\times}^{-1}:{\mathbb{R}}^{3\times 3}\to {\mathbb{R}}^{3}$.

Conversely, given a rotation matrix **R**, the corresponding rotation vector can be directly computed via the matrix logarithm operator $\mathrm{log}(\cdot ):{\mathbb{R}}^{3\times 3}\to {\mathbb{R}}^{3}$

It is important to note that the rotation axis **u** is expressed in the material frame of reference associated with the matrix **R** (or **Q**). With these tools in hand, we now proceed to outline our numerical scheme.

#### 3.2. Spatial discretization

Drawing from previous studies of unshearable and inextensible rods [13,14,67], we capture the deformation of a filament in three-dimensional space via the time evolution of a discrete set of vertices ${\mathbf{\text{r}}}_{i}(t)\in {\mathbb{R}}^{3},i\in [1,n+1]$ and a discrete set of material frames ${\mathbf{\text{Q}}}_{i}(t)\in {\mathbb{R}}^{3\times 3},i\in [1,n]$, as illustrated in figure 2.

Each vertex is associated with the following discrete quantities:

**v**

_{i}is the velocity ,

*m*

_{i}is a pointwise concentrated mass and

**F**

_{i}is the external force given in equation (2.14).

Each material frame is associated with an edge **ℓ**_{i} connecting two consecutive vertices, and with the related discrete quantities

_{i}=|

**ℓ**

_{i}|, ${\hat{\ell}}_{i}=|{\hat{\mathit{\ell}}}_{i}|$, ${e}_{i}={\ell}_{i}/{\hat{\ell}}_{i}$ are the edge current length, reference length and dilatation factor,

**t**

_{i}is the discrete tangent vector, ${\mathit{\sigma}}_{\mathcal{L}}^{i}$ is the discrete shear/axial strain vector, ${\mathit{\omega}}_{\mathcal{L}}^{i}$ is the discrete angular velocity, ${\hat{A}}_{i}$, ${\hat{\mathbf{\text{J}}}}_{i}$, ${\hat{\mathbf{\text{B}}}}_{i}$, ${\hat{\mathbf{\text{S}}}}_{i}$ are the edge reference cross-section area, mass second moment of inertia, bend/twist matrix and shear/stretch matrix, respectively, and finally ${\mathbf{\text{C}}}_{\mathcal{L}}^{i}$ is the external couple given in equation (2.15).

Whereas in the continuum setting (§2) all quantities are defined pointwise, in a discrete setting some quantities, and in particular ${\mathit{\kappa}}_{\mathcal{L}}$, are naturally expressed in an integrated form over the domain $\mathcal{D}$ along the filament [14,68]. Any integrated quantity divided by the corresponding integration domain length $\mathcal{D}=|\mathcal{D}|$ is equivalent to its pointwise average. Therefore, following the approach of Bergou & Grinspun [14,68], the domain $\mathcal{D}$ becomes the Voronoi region ${\mathcal{D}}_{i}$ of length

*interior*vertices

**r**

^{(int)}

_{i=1,…,n−1}. Each interior vertex is then also associated with the following discrete quantities:

**Q**

_{i}to its neighbour

**Q**

_{i+1}over the segment size ${\hat{\mathcal{D}}}_{i}$ along the rod. Therefore, the relation ${e}^{{\hat{\mathcal{D}}}_{i}{\hat{\mathit{\kappa}}}_{\mathcal{L}}^{i}}{\mathbf{\text{Q}}}_{i}={\mathbf{\text{Q}}}_{i+1}$ holds, so that ${\hat{\mathit{\kappa}}}_{\mathcal{L}}^{i}=\mathrm{log}({\mathbf{\text{Q}}}_{i+1}{\mathbf{\text{Q}}}_{i}^{\mathrm{T}})/{\hat{\mathcal{D}}}_{i}$. Finally, we introduce the bend/twist stiffness matrix ${\hat{\mathbf{\text{B}}}}_{i}$ consistent with the Voronoi representation.

Then, we may discretize the governing equations (2.12)–(2.15) so that they read

^{h}and ${\mathcal{A}}^{h}$ operate on a set of

*N*vectors and returns

*N*+1 vectors, consistent with equations (3.6)–(3.9) (see the appendix for further details).

#### 3.3. Time discretization

The derivation above leads to a system that in general is not Hamiltonian, as this depends on the nature of the external loads acting on the filament as well as on the choice of constitutive laws. Nonetheless, without external loads and under the assumptions of linear stress–strain relations, this derivation amounts to a geometric rescaling through *e* of the classic Cosserat rod model (here directly discretized via standard finite difference and trapezoidal quadrature rule, as outlined in [14,18,69,70]), which is a Hamiltonian system [18,71] with quadratic energy functionals [70]: translational ${E}_{T}=\frac{1}{2}{\int}_{0}^{L}\rho A{\mathbf{\text{v}}}^{\mathrm{T}}\cdot \mathbf{\text{v}}\hspace{0.17em}\mathrm{d}s$; rotational ${E}_{R}=\frac{1}{2}{\int}_{0}^{L}\rho {\mathit{\omega}}_{\mathcal{L}}^{\mathrm{T}}\mathbf{\text{I}}{\mathit{\omega}}_{\mathcal{L}}\hspace{0.17em}\mathrm{d}s$; bending/twist ${E}_{B}=\frac{1}{2}{\int}_{0}^{L}{\mathit{\kappa}}_{\mathcal{L}}^{\mathrm{T}}\mathbf{\text{B}}\mathit{\kappa}{;}_{\mathcal{L}}\hspace{0.17em}\mathrm{d}s$; shear/stretch ${E}_{S}=\frac{1}{2}{\int}_{0}^{L}{\mathit{\sigma}}_{\mathcal{L}}^{\mathrm{T}}\mathbf{\text{S}}{\mathit{\sigma}}_{\mathcal{L}}\hspace{0.17em}\mathrm{d}s$. Therefore, by construction, in the limit of $e\to 1$ (no axial elongation) the equivalence with the Cosserat model holds, and for consistency we opt for an energy-preserving time integrator, and in particular a symplectic, second-order Verlet scheme. We note that, despite the failure of Verlet schemes to integrate rotational equations of motion when represented by quaternions, in our case their use is justified as rotations are represented instead by Euler angles [72].

The second-order position Verlet time integrator is structured in three blocks: a first half-step updates the linear and angular positions, followed by the evaluation of local linear and angular accelerations, and finally a second half-step updates the linear and angular positions again. Therefore, it entails only one right-hand side evaluation of equations (3.8) and (3.9), the most computationally expensive operation (see the appendix for details).

This algorithm strikes a balance between computing costs, numerical accuracy and implementation modularity: by foregoing an implicit integration scheme we can incorporate a number of additional physical effects and soft constraints, even though this may come at the expense of computational efficiency. Indeed, for large Young’s or shear moduli or for very thin rods, the system of equations (2.4)–(2.7) might become stiff, so that small timesteps must be employed to ensure stability. Although we have not derived a rigorous CFL-like condition, throughout this work we employed the empirical relation d*t*=*a* dℓ, with *a*∼10^{−2} s/m, and found it reliable in preventing numerical instabilities. Despite the potential stiffness of this model, as the computational cost per timestep scales linearly with the number of discretization elements *n* (the model is one-dimensional), we could carry out all our computation on a laptop (details and representative timings are summarized in the appendix, figure 20). It is worth noting that because of the condition d*t*=*a* dℓ, a smaller number of elements implies larger timesteps (dℓ=*L*/*n*). Therefore, the time-to-solution may scale between linearly and quadratically, depending on how the timestep is set.

#### 3.4. Validation

We first validate our proposed methodology against a number of benchmark problems with analytic solutions and examine the convergence properties of our approach. Three case studies serve to characterize the competition between bending and twisting effects in the context of helical buckling, dynamic stretching of a loaded rod under gravity, and the competition between shearing and bending in the context of a Timoshenko beam. Further validations reported in the appendix include Euler and Mitchell buckling due to compression or twist, and stretching and twisting vibrations.

##### 3.4.1. Helical buckling instability

We validate our discrete derivative operators beyond the onset of instability (see Euler and Mitchell buckling tests in the appendix) for a long straight, isotropic, inextensible, and unshearable rod undergoing bending and twisting. The filament is characterized by the length *L* and by the bending and twist stiffnesses *α* and *β*, respectively. The clamped ends of the rod are pulled together in the axial direction **k** with a slack *D*/2 and simultaneously twisted by the angle *Φ*/2, as illustrated in figure 3*a*. Under these conditions, the filament buckles into a localized helical shape (figure 3*e*).

The nonlinear equilibrium configuration **r**_{eq} of the rod can be analytically determined [8,73–75] in terms of the total applied slack *D* and twist *Φ*. We denote the magnitude of the twisting torque and tension acting on both ends and projected on **k** by *M*_{h} and *T*_{h}, respectively. Their normalized counterparts *m*_{h}=*M*_{h}*L*/(2*πα*) and *t*_{h}=*T*_{h}*L*^{2}/(4*π*^{2}*α*) can be computed via the ‘semi-finite’ correction approach [74] by solving the system

**r**

_{eq}can be expressed [75] as

**k**is necessary. Following Bergou

*et al.*[14], we rely on the definition of the envelope

*φ*

*θ*is the angular deviation of the tangent

**t**from the axial direction

**k**, and ${\theta}_{max}$ is the corresponding maximum value along the filament. The envelope

*φ*relative to the analytical solution of equation (3.10), and

*φ*

^{n}relative to a numerical model of

*n*discretization elements can be estimated via finite differences. This allows us to determine the convergence order of the solver by means of the norms

*L*

^{1}(

*ϵ*),

*L*

^{2}(

*ϵ*) and ${L}^{\mathrm{\infty}}(\u03f5)$ of the error

*ϵ*=∥

*φ*−

*φ*

^{n}∥.

We simulate the problem illustrated in figure 3 at different space–time resolutions. The straight rod originally at rest is twisted and compressed at a constant rate during the period *T*_{twist}. Subsequently, the ends of the rod are held in their final configurations for the period *T*_{relax} to allow the internal energy to dissipate (according to the model of §(a)) until the steady state is reached. Simulations are carried out progressively refining the spatial discretization *δl*=*L*/*n* by varying *n*=100–3200 and the time discretization *δt* is kept proportional to *δl*, as reported in figure 3.

As can be seen in figure 3*b*,*c*,*e* the numerical solutions converge to the analytical one with second order in time and space, consistent with our spatial and temporal discretization schemes. Moreover, to validate the energy-conserving properties of the solver in the limit of $e\to 1$, we turn off the internal dissipation (figure 3*d*) and observe that the total energy of the filament *E*_{F} is constant after *T*_{twist} and matches its theoretical value *E*_{F}=(*M*_{h}*Φ*+*T*_{h}*D*)/2.

##### 3.4.2. Vertical oscillations under gravity

We consider a system in which a rod hanging from one end and subject to gravity *g* oscillates due to a mass *m*_{p} suspended at the other end, and due to its own mass *m*_{r}, as depicted in figure 4*a*,*d*. This system is analogous to a mass-spring oscillator. The static solution is then obtained by integrating the infinitesimal elongations along the spring due to the local load [76], yielding the total equilibrium extension

*k*is the spring constant,

*ξ*=2 is a constant factor and

*m*

_{eq}=

*m*

_{p}+

*m*

_{r}/

*ξ*is the equivalent mass. Thus, the final equilibrium length of the rod reads $L=\hat{L}+\mathrm{\Delta}{L}^{\ast}$, with $\hat{L}$ being the rest unstretched length.

The dynamic solution is instead characterized by oscillations of period *T** and by the time-varying length *L*(*t*) of the spring

*ξ*depends on the ratio

*m*

_{r}/

*m*

_{p}. In fact it can be shown [76] that

*ξ*≃3 for ${m}_{r}/{m}_{p}\to 0$, and

*ξ*≃

*π*

^{2}/4 for ${m}_{r}/{m}_{p}\to \mathrm{\infty}$.

The analytical results rely on the assumption of *k* being constant in space and time, given a fixed ratio *m*_{r}/*m*_{p}. However, this condition is not met here because *k*(*s*,*t*)=*EA*(*s*,*t*) is a function of space and time, due to dilatation and mass conservation. Nevertheless, as the Young’s modulus $E\to \mathrm{\infty}$, that is as a soft filament becomes stiff, the constant $k\to E\hat{A}$ and our rod model must recover the behaviour of the mass-spring oscillator. Indeed, figure 4*b*,*c*,*e*,*f* shows how the proposed numerical method converges to the analytical oscillation period *T** and normalized longitudinal displacement $(L-\hat{L})/\mathrm{\Delta}{L}^{\ast}$ as *E* increases.

##### 3.4.3. Cantilever beam

We now consider the effect of bend and shear simultaneously by validating our numerical methods against the Timoshenko cantilever of figure 5*a*. Timoshenko’s model accounts for bending elasticity, rotary inertia and shear deformations, building on classical beam theories by Rayleigh (bending elasticity and rotary inertia) and Euler–Bernoulli (bending elasticity only). The model captures the behaviour of short or composite beams in which shear deformations effectively lower the stiffness of the rod [58,77].

We consider a beam clamped at one end $\hat{s}=0$ and subject to the downward force *F* at the free end $\hat{s}=\hat{L}$, as illustrated in figure 5*a*. The static solution for the displacement *y* along the vertical direction **i** of the rod can then be analytically expressed as

**j**=

**k**×

**i**,

*E*and

*G*are, respectively, the Young’s and shear moduli, and ${\alpha}_{c}=\frac{4}{3}$ is the Timoshenko shear factor for circular sections and accounts for the fact that the shear stress varies over the section [58]. Furthermore, the Timoshenko (as well as Rayleigh and Euler–Bernoulli) theory relies on the assumption of small deflections, so that the horizontal coordinate

*x*along the direction

**k**can be approximated by the arc-length $\hat{s}$ (see figure 5

*a*and the appendix for further details and derivation), hence the use of $\hat{s}$ in the above equation.

If the shear modulus *G* approaches infinity or if the ratio $E{\hat{I}}_{1}/({\alpha}_{c}{\hat{L}}^{2}\hat{A}G)\gg 1$, then the Timoshenko model in the static case reduces to the Euler–Bernoulli approximation, yielding

We compare our numerical model with these results by carrying out simulations of the cantilever beam of figure 5*a* in the time–space limit of refinement. As can be noticed in figure 5*b* the discrete solution recovers the Timoshenko one. Therefore, the solver correctly captures the role of shear that reduces the effective stiffness relative to the Euler–Bernoulli solution. Moreover, our approach is shown to converge to the analytical solution in all the norms ${L}^{\mathrm{\infty}}(\u03f5)$, *L*^{1}(*ϵ*), *L*^{2}(*ϵ*) of the error *ϵ*=∥*y*−*y*^{n}∥, where *y*^{n} is the vertical displacement numerically obtained in the refinement limit.

We note that the norms ${L}^{\mathrm{\infty}}(\u03f5)$ and *L*^{1}(*ϵ*) exhibit first-order convergence, while *L*^{2}(*ϵ*) decays with a slope between first and second order. We attribute these results to the fact that while the Timoshenko solution does not consider axial extension or tension, it does rely on the assumption of small deflections ($\hat{s}=x$), therefore effectively producing a dilatation of the rod. On the contrary, our solver does not assume small deflections and does not neglect axial extension, because the third entry of the matrix **B** has the finite value $E\hat{A}$ (see figure 5 for details). This discrepancy is here empirically observed to decrease the convergence order.

These studies, together with the ones reported in the appendix, complete the validation of our numerical implementation and demonstrate the accuracy of our solver in simulating soft filaments in simple settings.

### 4. Including interactions and activity: solid and liquid friction, contact and muscular effects

Motivated by advancements in the field of soft robotics [41,44,45], we wish to develop a robust and accurate framework for the characterization and computational design of soft slender structures interacting with complex environments. To this end, we expand the range of applications of our formalism by including additional physical effects, from viscous hydrodynamic forces in the slender-body limit and surface solid friction to self-contact and active muscular activity. As a general strategy, all new *external* physical interactions are accounted for by lumping their contributions into the external forces and couples **F** and ${\mathbf{\text{C}}}_{\mathcal{L}}$ on the right-hand side of the linear and angular momentum balance equations (2.14) and (2.15). On the other hand, all new *internal* physical and biophysical effects are captured by adding their contributions directly to the internal force ${\mathbf{\text{n}}}_{\mathcal{L}}$ and torque ${\mathit{\tau}}_{\mathcal{L}}$ resultants before integrating equations (2.14) and (2.15).

#### 4.1. Dissipation

Real materials are subject to internal friction and viscoelastic losses, which can be modelled by modifying the constitutive relations so that the internal torques ${\mathit{\tau}}_{\mathcal{L}}({\mathit{\kappa}}_{\mathcal{L}})$ and forces ${\mathbf{\text{n}}}_{\mathcal{L}}({\mathit{\sigma}}_{\mathcal{L}})$ of equation (2.11) become functions of both strain and rate of strain, i.e. ${\mathit{\tau}}_{\mathcal{L}}({\mathit{\kappa}}_{\mathcal{L}},{\mathrm{\partial}}_{t}{\mathit{\kappa}}_{\mathcal{L}})$ and ${\mathbf{\text{n}}}_{\mathcal{L}}({\mathit{\sigma}}_{\mathcal{L}},{\mathrm{\partial}}_{t}{\mathit{\sigma}}_{\mathcal{L}})$. Keeping track of the strain rates increases computational costs and the memory footprint of the solver. However, for the purpose of purely dissipating energy, a simple alternative option is to employ Rayleigh potentials [16,78]. In this case, viscous forces **f**_{v} and torques ${\mathbf{\text{c}}}_{\mathcal{L}}^{v}$ per unit length are directly computed as linear functions of linear and angular velocities through the constant translational *γ*_{t} and rotational *γ*_{r} internal friction coefficients, so that

*γ*, therefore assuming

*γ*

_{t}=

*γ*

_{r}.

##### 4.1.1. Muscular activity

To study limbless biolocomotion on solid substrates and in fluids, we allow our soft filaments to be active, by generating internal forces and torques corresponding to coordinated muscular activity driven, for example, by a central pattern generator [79,80].

Following the approach detailed in [53,56], we express the muscular activity magnitude *A*_{m} as a travelling wave propagating head to tail along the filament

*ϕ*

_{m}is the phase,

*t*is time and

*T*

_{m}and λ

_{m}are, respectively, the activation period and wavelength. The amplitude of the travelling wave is represented by the cubic B-spline $\beta (\hat{s})$ characterized by

*N*

_{m}control points $({\hat{S}}_{i},{\beta}_{i})$ with $i=0,\dots ,{N}_{m}-1$, as illustrated in figure 6. The first and last control points are fixed so that $({\hat{S}}_{0},{\beta}_{0})=(0,0)$ and $({\hat{S}}_{{N}_{m}-1},{\beta}_{{N}_{m}-1})=(\hat{L},0)$, therefore assuming the ends of the deforming body to be free. One of the main advantages of the proposed parametrization is that it encompasses a large variety of patterns with a relatively small number of parameters [56].

A given activation mode can be achieved by multiplying the scalar amplitude *A*_{m} with the appropriate director. For example, if we wish to study earthworm-like locomotion, we may employ a wave of longitudinal dilatation and compression forces, so that

**d**

_{2}and

**d**

_{3}to be coplanar to the ground. These two contributions are directly added to the internal force ${\mathbf{\text{n}}}_{\mathcal{L}}$ and torque ${\mathit{\tau}}_{\mathcal{L}}$ resultants.

In the most general case, all deformation modes can be excited by enabling force and torque muscular activity along all directors **d**_{1}, **d**_{2} and **d**_{3}, providing great flexibility in terms of possible gaits.

##### 4.1.2. Self-contact

To prevent the filament from passing through itself, we need to account for self-contact. As a general strategy, we avoid enforcing the presence of boundaries via Lagrangian constraints as their formulation may be cumbersome [81], impairing the modularity of the numerical solver. We instead resort to calculating forces and torques directly and replacing hard constraints with ‘soft’ displacement–force relations.

Our self-contact model introduces additional forces **F**_{sc} acting between the discrete elements in contact. To determine whether any two cylindrical elements are in contact, we calculate the minimum distance ${d}_{min}^{ij}$ between edges *i*,*j* by parametrizing their centre lines *c*_{i}(*h*)=*s*_{i}+*h*(*s*_{i+1}−*s*_{i}) so that

*r*

_{i}and

*r*

_{j}are the radii of edges

*i*and

*j*, respectively. If

*ϵ*

_{ij}is smaller than zero, then the two edges are not in contact and no penalty is applied. Denoting as ${\mathbf{\text{d}}}_{min}^{ij}$ the unit vector pointing from the closest point on edge

*i*to the closest point on edge

*j*, the self-contact repulsion force is given by

*H*(

*ϵ*

_{ij}) denotes the Heaviside function and ensures that a repulsion force is produced only in case of contact (

*ϵ*

_{ij}≥0). The first term within the square brackets expresses the linear response to the interpenetration distance as modulated by the stiffness

*k*

_{sc}, while the second damping term models contact dissipation and is proportional to the coefficient

*γ*

_{sc}and the interpenetration velocity

**v**

_{i}−

**v**

_{j}.

##### 4.1.3. Contact with solid boundaries

To investigate scenarios in which filaments interact with the surrounding environment, we must also account for solid boundaries. By implementing the same approach outlined in the previous section, obstacles and surfaces are modelled as soft boundaries allowing for interpenetration with the elements of the rod (figure 7). The wall response **F**^{w}_{⊥} balances the sum of all forces **F**_{⊥} that push the rod against the wall, and is complemented by the other two components which help prevent possible interpenetration due to numerics. The interpenetration distance *ϵ* triggers a normal elastic response proportional to the stiffness of the wall, while a dissipative term related to the normal velocity component of the filament with respect to the substrate accounts for a damping force, so that the overall wall response reads

*H*(

*ϵ*) denotes the Heaviside function and ensures that a wall force is produced only in case of contact (

*ϵ*≥0). Here,

**u**

^{w}

_{⊥}is the boundary outward normal (evaluated at the contact point, that is the contact location for which the normal passes through the centre of mass of the element), and

*k*

_{w}and

*γ*

_{w}are, respectively, the wall stiffness and dissipation coefficients.

##### 4.1.4. Isotropic and anisotropic surface friction

Solid boundaries also affect the dynamics of the filament through surface friction, a complex physical phenomenon in which a range of factors are involved, from roughness and plasticity of the surfaces in contact to the kinematic initial conditions and geometric set-up. Here, we adopt the Amonton–Coulomb model, the simplest of friction models.

This model relates the normal force pushing a body onto a substrate to the friction force through the kinetic *μ*_{k} and static *μ*_{s} friction coefficients, depending on whether the contact surfaces are in relative motion or not.

Despite the simplicity of the model, its formulation and implementation may not necessarily be straightforward, especially in the case of rolling motions. Given the cylindrical geometry of our filaments, the effect of surface friction can be decomposed into a longitudinal component associated with purely translational displacements, and a lateral component associated with both translational and rotational motions (figure 8). We use the notation **x**_{⊥}, **x**_{∥}, **x**_{×} to denote the projection of the vector **x** in the directions **u**^{w}_{⊥}, ${\mathbf{\text{u}}}_{\parallel}^{\mathrm{w}}$, **u**^{w}_{×}, as illustrated in figure 8.

The longitudinal friction force **F**_{long} is opposite to either the resultant of all forces **F**_{×} acting on an element (static case) or to the translational velocity **v**_{×} (kinetic case) along the direction **u**^{w}_{×} (figure 8). The Amonton–Coulomb model then reads

**v**

_{×}|≤

*v*

_{ϵ}) and kinetic (|

**v**

_{×}|>

*v*

_{ϵ}) cases. We define

*v*

_{ϵ}in a limit form to accommodate the fact that inequalities are numerically evaluated up to a small threshold value. The static friction force is always equal and opposite to

**F**

_{×}up to a maximum value proportional to the normal force |

**F**

_{⊥}| through the coefficient

*μ*

_{s}. The kinetic friction force is instead opposite to the translational velocity

**v**

_{×}, but does not depend on its actual magnitude and is proportional to |

**F**

_{⊥}| via

*μ*

_{k}. In general

*μ*

_{s}>

*μ*

_{k}, so that it is harder to set a body into motion from rest than to drag it.

The lateral displacement of a filament in the direction ${\mathbf{\text{u}}}_{\parallel}^{\mathrm{w}}={\mathbf{\text{u}}}_{\times}^{\mathrm{w}}\times {\mathbf{\text{u}}}_{\perp}^{\mathrm{w}}$ is associated with both translational (**v**_{∥}) and rotational (${\mathit{\omega}}_{\times}={\omega}_{\times}{\mathbf{\text{u}}}_{\times}^{\mathrm{w}}$) motions, as illustrated in figure 8*b*,*c*. In this case, the distinction between static and kinetic friction does not depend on **v**_{∥}, but on the relative velocity (also referred to as slip velocity) between the rod and the substrate

**v**

_{cont}is the local velocity of the filament at the contact point with the substrate, due to the axial component of the angular velocity

*ω*_{×}.

In the static or no-slip scenario (**v**_{slip}=**0**), the linear momentum balance in the direction **u**^{w}_{∥}, and the angular momentum balance about the axis **u**^{w}_{×} express a kinematic constraint between the linear acceleration *a***u**^{w}_{∥} and angular acceleration ${\mathit{\omega}}_{\times}=({\mathbf{\text{u}}}_{\perp}^{\mathrm{w}}\times a{\mathbf{\text{u}}}_{\parallel}^{\mathrm{w}})/r$, so that

**u**

^{w}

_{×}is

*J*=

*r*

^{2}d

*m*/2, the above system can be solved for the unknown

*a*and

*F*

_{roll}, yielding

Therefore, the lateral friction force **F**_{lat} and the associated torque ${\mathbf{\text{C}}}_{\mathcal{L}}^{\mathrm{lat}}$ can be finally expressed as

*μ*

^{r}

_{s}and

*μ*

^{r}

_{k}are, respectively, the rolling static and kinetic friction coefficients.

So far we have considered isotropic friction by assuming that the coefficients *μ*_{s} and *μ*_{k} are constant and independent of the direction of the total acting forces (static case) or relative velocities (kinetic case). Nevertheless, frictional forces may be highly anisotropic. For example, the anisotropy caused by the presence of scales on the body of a snake crucially affects gaits and performance [82,83].

The Amonton–Coulomb model can be readily extended to account for anisotropic effects by simply assuming the friction coefficients *μ*_{s} and *μ*_{k} to be functions of a given reference direction. The nature of these functions depends on the specific physical problems under investigation. An example of this approach is illustrated in §(b) in the context of limbless locomotion. Isotropic and anisotropic friction validation benchmarks are presented in the appendix.

##### 4.1.5. Hydrodynamics

We also extend our computational framework to address flow–structure interaction problems. In particular, we consider the case in which viscous forces dominate over inertial effects, i.e. we consider systems in which the Reynolds number *Re*=*ρ*_{f}*UL*/*μ*≪1 where *ρ*_{f} and *μ* are the density and dynamic viscosity of the fluid, respectively, and *U* is the characteristic velocity of the rod. Under these conditions, the drag forces exerted by the fluid on our filaments can be determined analytically within the context of slender-body theory [84,85] (for more advanced and accurate viscous flow models we refer the reader to [31–36]). At leading order resistive force line densities scale linearly with the local rod velocities **v** according to

### 5. Applications

We now proceed to illustrate the capabilities of our solver with three different applications. We consider first a static problem in which self-contact, bending and twist give rise to the classic out-of-plane configurations denoted as plectonemes [40], while the addition of stretching and shearing produces a different type of experimentally observed solutions, known as solenoids [40]. Then we turn our attention to two dynamic biophysical problems in which an active filament interacts with a solid and a liquid environment, exhibiting qualitatively different optimal biolocomotion strategies. We emphasize here that one of the main points of these biolocomotion studies, besides verifying the versatility and prediction capabilities of our solver in biophysical settings, is to demonstrate its practical use in an inverse (optimization) design process, which typically represents a demanding stress test for any simulation algorithm, as unusual parameter sets or conditions are explored.

#### 5.1. Plectonemes and solenoids

When an inextensible rod is clamped at one end and twisted a sufficiently large number of times at the other end, it becomes unstable, coils up and generates a characteristic structure known as a plectoneme [87]. While this behaviour has been well characterized both theoretically and experimentally [87], its analogue for highly extensible filaments has been largely ignored [40]. In particular, for large extensional and twisting strains qualitatively different solutions arise, such as those corresponding to tightly packed solenoidal structures [40] whose properties are as yet poorly understood, and whose importance has only recently been recognized in the context of soft robotics and artificial muscles [41–44].

Given the broad scope of our computational framework, we can now study the formation of both solenoids and plectonemes. As illustrated in figure 9*a*, a soft rod of Young’s modulus *E*=10^{6} Pa is clamped at one end, and subject to an axial load *F*, while also being twisted *R* times at the other end. As experimentally and theoretically observed for *F*=0, i.e. in the absence of stretching ($L/\hat{L}\approx 1$), plectonemes are generated (figure 9*b*). When the load *F* is increased so that the elongation of the rod approaches $L/\hat{L}\approx 1.15$, solenoids arise as predicted in [40] and illustrated in figure 9*c*. This test case, therefore, shows the ability of our solver to qualitatively capture different instability mechanisms, driven by the competition between the different modes of deformation of the rod. We leave the details of the explanation of the phase diagram for the formation of plectonemes, solenoids and intermediate structures [40] for a later study.

#### 5.2. Slithering

The mechanics of slithering locomotion typical of snakes has been extensively investigated experimentally [83,88,89], theoretically [82,90,91] and computationally [92,93]. While biological experiments have provided quantitative insights, theoretical and computational models have been instrumental to characterize qualitatively the working principles underlying snake locomotion. Although these models implement different levels of realism, they generally rely on a number of key simplifications. Typically, theoretical models assume planar deformations [82] and/or disregard mechanics by prescribing body kinematics [91]. Computational models offer a more realistic representation, but they have mostly been developed for and tailored to robotic applications [92,93]. For example, snakes are often modelled as a relatively small set of hinges and/or springs representing pointwise localized actuators that connect contiguous rigid segments. Therefore, they do not account for the continuum nature of elastic body mechanics and biological muscular activity. Moreover, in robot replicas the critical feature of friction anisotropy is commonly achieved through the use of wheels [89]. As a consequence computational models often assume only two sources of anisotropy, in the tangential and lateral direction with respect to the body. This is in contrast with biological experiments [83] that highlight the importance of all three sources of anisotropy, namely forward, backward and lateral.

Our approach complements those previous attempts by accounting for physical and biological effects within a continuum framework (equations (2.12)–(2.15)). In this section, we demonstrate the qualitative and quantitative capabilities of our solver by reverse engineering optimal slithering gaits that maximize forward speed.

We consider a soft filament of unit length actuated via a planar travelling torque wave of muscular activity in the direction perpendicular to the ground. The interaction with the substrate is characterized by the ratios ${\mu}_{\mathrm{k}}^{f}:{\mu}_{\mathrm{k}}^{b}:{\mu}_{\mathrm{k}}^{r}=1:1.5:2$ and ${\mu}_{\mathrm{s}}^{f}:{\mu}_{\mathrm{s}}^{b}:{\mu}_{\mathrm{s}}^{r}=1:1.5:2$ with ${\mu}_{\mathrm{s}}^{f}=2{\mu}_{\mathrm{k}}^{f}$, as experimentally observed for juvenile Pueblan milk snakes on a moderately rough surface [83]. The value of the friction coefficient *μ*^{f}_{k} is set so that the ratio between inertial and friction forces captured by the Froude number is $Fr=(L/{T}_{m}^{2})/({\mu}_{\mathrm{k}}^{f}g)=0.1$, as measured for these snakes [83].

In the spirit of [53,56,57], we wish to identify the fastest gaits by optimizing the filament muscular activity. The torque wave generated by the snake is parametrized according to §i and is characterized by *N*_{m}=6 control points and a unit oscillation period *T*_{m}, so that overall we optimize for five parameters, four of which are responsible for the torque profile along the rod (*β*_{1}, *β*_{2}, *β*_{3}, *β*_{4}), while the last one represents the wavenumber 2*π*/λ_{m} (see §i).

These parameters are left free to evolve from an initial zero value, guided by an automated optimization procedure that identifies the optimal values that maximize the snake’s forward average speed *v*^{fwd}_{max} over one activation cycle *T*_{m}. The algorithm of choice is the Covariance Matrix Adaptation-Evolution Strategy [94,95] (CMA-ES) which has proved effective in a range of biophysical and engineering problems, from the optimization of swimming gaits [53], morphologies [56,57] and collective dynamics [55] to the identification of aircraft alleviation schemes [96] or virus traffic mechanisms [52]. The CMA-ES is a stochastic optimization algorithm that samples generations of *p* parameter vectors from a multivariate Gaussian distribution $\mathcal{N}$. Here each parameter vector represents a muscular activation instance, and every generation entails the evaluation of *p*=60 different gaits. The covariance matrix of the distribution $\mathcal{N}$ is then adapted based on successful past gaits, chosen according to their corresponding cost function value *f*=*v*^{fwd}_{max}, until convergence to the optimum.

The course of the optimization is reported in figure 10 together with the kinematic details of the identified fastest gait. As can be noticed in figure 10*e*,*f* the forward scaled average speed approaches *v*^{fwd}_{max}≃0.6, consistent with experimental evidence [97]. Moreover, CMA-ES finds that the optimal wavelength is λ_{m}≃*L* (figure 10*g*), again consistent with biological observations [83,98]. Thus, this value of wavelength strikes a balance between thrust production and drag minimization within the mechanical constraints of the system.

We note that a rigorous characterization of slithering locomotion would require the knowledge of a number of biologically relevant parameters (Young’s and shear moduli of muscular tissue, maximum attainable torques, etc.) and environmental conditions (terrain asperities, presence of pegs, etc.) and goes beyond the scope of the present work. Nevertheless, this study illustrates the robustness, quantitative accuracy and suitability of our solver for the characterization of biolocomotion phenomena.

##### 5.2.1. Swimming

We finally turn to apply the inverse design approach outlined in the previous section to the problem of swimming at low Reynolds numbers where viscous forces dominate inertial effects. We maintain the exact same set-up as in the slithering case, while we change the environment from a solid substrate to a viscous fluid. The flow-filament interaction is then modelled via slender-body theory, as illustrated in §v.

Once again we inverse design planar optimal gaits for forward average speed *f*=*v*^{fwd}_{max} within one activation cycle *T*_{m}, by employing the same muscular activity parametrization as for slithering. To verify *a posteriori* the biological relevance of the identified optimal solution, we consider the case of the sea urchin spermatozoon *Echinus esculentus* [99], which swims by means of helical or planar waves travelling along its flagella of length *L*_{s}≃40 μm. The gait corresponding to planar swimming is characterized by kinematic undulations of wavelength λ_{s}<*L*_{s} and frequency *f*_{s}≃2.8. At *Re*≃10^{−4} the spermatozoon attains the scaled velocity *v*_{s}=*U*_{s}/(*f*_{s}*L*_{s})≃0.08±0.03, where *U*_{s} is the dimensional cruise speed [99]. Although this gait may not be the absolute optimal planar locomotion pattern, the fact that it is replicated in a large number of organisms [85] suggests that it captures some effective features that we expect to qualitatively recover via our numerical optimization.

The course of the optimization is reported in figure 11 together with the kinematic details of the identified fastest gait. As can be noted in figure 11*e*,*f*, the forward average scaled speed and wavelength approach *v*^{fwd}_{max}≃0.055 and λ_{m}≃0.38*L* are qualitatively and quantitatively consistent with experimental evidence [99].

As in the previous section, we note that a rigorous characterization of optimal swimming at low Reynolds numbers would require the knowledge of a number of biologically relevant parameters and environmental conditions, a careful comparison with a rich body of literature [33,35,36,85], and goes beyond the scope of the present work. Nevertheless, this and the previous study illustrate how the interplay between filament mechanics and the surrounding environment crucially affects propulsive gaits, as is biologically evident and automatically recovered via our numerical inverse design approach.

### 6. Conclusions

We have presented a versatile implementation of the Cosserat rod model for the simulation of soft filaments deforming in three-dimensional space. These filaments at any given cross section can experience all deformation modes, namely normal and orthonormal bending, twisting, stretching and shearing. Particular emphasis is placed in realistically accounting for substantial axial extensions and shear strains, targeting emerging soft robotic applications. The solver also handles self-contact, muscular activity, solid boundaries, isotropic and anisotropic friction as well as hydrodynamics. The outcome is a simple algorithm suitable to tackle a plethora of engineering and biophysical phenomena.

We verified this against a battery of benchmark problems entailing different physical aspects and boundary conditions, and used it to solve a range of forward and inverse problems spanning mechanics and biophysics. In particular, we studied the formation of solenoids and plectonemes and the evolutionary optimization of terrestrial limbless locomotion and swimming.

Our results demonstrate the utility of a modelling approach based on Cosserat rods in a wide range of settings that involve active and passive soft filaments, broadening the scope of previous studies. We particularly emphasize that our robust computational approach allows us to approach the solution of an inverse problem by coupling our method with evolutionary strategies successfully, without the problems associated with the (often unforeseen) variety of candidate solutions produced throughout the process. Ongoing work involves its coupling to realistic high Reynolds number flow solvers [100], its integration with sensory feedback models for the characterization of locomotory neural circuitry [80], and its use for the rational design of biohybrid robots [45,101].

### Data accessibility

Our software implementing the study cases presented in the manuscript and in the appendix is publicly available at https://github.com/mattialab/elastica.git.

### Authors' contributions

M.G. and L.M. designed the study. M.G. and A.G.M. implemented the software. M.G. and L.H.D. carried out the analysis of the results. M.G., L.H.D. and L.M. wrote the manuscript.

### Competing interests

The authors declare no competing interests.

### Funding

We thank the Swiss National Science Foundation and the Blue Waters project (OCI-0725070, ACI-1238993), a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications, for partial support (M.G.), and the Harvard MRSEC grant NSF DMR 14-20570, Harvard Biomatter grant NSF DMR 33985 and ARO W911NF-15-1-016 for additional support (L.M.).

## Acknowledgements

We thank Nicholas Charles for a careful reading of and comments on the paper.

## Appendix A. Lagrangian governing equations

The following shows the conversion of governing equations for *unstretchable* filaments from the laboratory to material frame of reference. Note that in the Governing equations and Numerical method sections of the main text, all terms are scaled appropriately by the local dilatation to account for stretching.

The time and space derivatives of the centreline **r**(*s*,*t*) are associated with the velocity **v** and tangent field **t**

**t**|=1, since

*s*is the current arc-length.

Similarly, time and space derivatives of the material frame **Q**, due to the orthonormality of the directors, are associated by definition with the angular velocity ** ω** and generalized curvature

**vectors, so that**

*κ*_{t}

**Q**

^{T}⋅

**Q**=

**×(⋅) and ∂**

*ω*_{s}

**Q**

^{T}⋅

**Q**=

**×(⋅) hold. These kinematic equations combined with the linear and angular momentum balance laws at a cross section [102] yield the governing equations for the Cosserat rod**

*κ**ρ*is the constant material density,

*A*is the cross-sectional area in its current state (so that

*ρA*

**v**is the linear momentum line density),

**h**(

*s*,

*t*) is the angular momentum line density,

**n**(

*s*,

*t*) and

**(**

*τ**s*,

*t*) are the internal force and torque resultants, and

**f**and

**c**are external body force and torque line densities.

The internal force **n** and torque ** τ** resultants depend on the geometric and material properties of the filament. This dependence is embedded via the material or constitutive laws, which must be frame invariant, rendering their definition in a Lagrangian setting most natural. Moreover, we note that the angular momentum line density can be readily expressed in the material frame as $\mathbf{\text{Q}}\mathbf{\text{h}}={\mathbf{\text{h}}}_{\mathcal{L}}=\rho \mathbf{\text{I}}{\mathit{\omega}}_{\mathcal{L}}$, where the tensor

**I**is the second area moment of inertia which, assuming circular cross sections, takes the form

To derive the Lagrangian form of the angular momentum equation, we begin with the definition of the material frame

We can use these definitions to establish relationships between time and space derivatives of the laboratory and material frames.

**I**in the laboratory frame is a function of space and time rendering its use cumbersome, while ${\mathbf{\text{I}}}_{\mathcal{L}}$ is a constant and, in our isotropic case, diagonal matrix in the material frame). We can use equations (A 16) and (A 17) to convert all of the terms in equation (A 7) to the material frame.

**Q**and noting that $\mathbf{\text{I}}\mathit{\omega}={\mathbf{\text{Q}}}^{\mathrm{T}}{\mathbf{\text{I}}}_{\mathcal{L}}{\mathit{\omega}}_{\mathcal{L}}$ yields our final governing equation for the angular momentum, expressed in the material frame

## Appendix B. Axial stretch comparison with Neo-Hookean and Mooney–Rivlin models

We consider the simple example of a beam subject to axial stretching and compare the load response predicted by our approach (that combines linear stress–strain constitutive laws and cross-sectional area rescaling) with Neo-Hookean [60] and Mooney–Rivlin models [61,62] for hyperelastic materials. We assume an initial length $\hat{L}=100\hspace{0.17em}\mathrm{mm}$ and cross section $\hat{A}=10\hspace{0.17em}{\mathrm{mm}}^{2}$. As the beam is stretched the cross-sectional area decreases. In one-dimensional, the Neo-Hookean model reads:

*α*is the stress, $e=L/\hat{L}$ is the elongation factor (related to the strain

*σ*=

*e*−1),

*C*

_{10}is a material property to be determined experimentally and that in the limit of small strains ($e\to 1$) is related to the Young’s modulus

*E*, as above indicated. In one-dimensional, the Mooney–Rivlin model reads:

*C*

_{10}and

*C*

_{01}are material properties to be determined experimentally and that in the limit of small strains ($e\to 1$) are related to the Young’s modulus

*E*as indicated above. In our approach the load reads

*E*=10

^{7}Pa, and that in the Mooney–Rivlin model

*C*

_{10}=

*C*

_{01}, we obtain the load responses illustrated in figure 12. As can be noted our approach reasonably approximates more advanced hyperelastic models within a 30% deformation range.

## Appendix C. Time discretization

To advance the discrete system of equations (3.6)–(3.9), we choose the energy conserving, second-order position Verlet time integration scheme. This allows us to write a full iteration from time *t* to *t*+*δt* as

## Appendix D. Discrete operators

The discrete operators ${\mathrm{\Delta}}^{h}:{\{{\mathbb{R}}^{3}\}}_{N}\to {\{{\mathbb{R}}^{3}\}}_{N+1}$ and ${\mathcal{A}}^{h}:{\{{\mathbb{R}}^{3}\}}_{N}\to {\{{\mathbb{R}}^{3}\}}_{N+1}$ take the form

^{h}and ${\mathcal{A}}^{h}$ operate on a set of

*N*vectors and return

*N*+1 vectors of differences or convert integrated quantities to point-wise, respectively, consistent with equations (3.8) and (3.9).

## Appendix E. Derivation of the vertical displacement of a Timoshenko cantilever beam

We briefly derive here the analytical expression of the vertical displacement *y* of equation (3.15) for the cantilever beam problem of figure 5*a*. To do so we make use of the free body diagram of figure 13, and of the constitutive relations of table 1. Moreover, we disregard axial extension and assume small deformations so that the coordinate *x* (along the direction **k**) is approximated by the arc-length *s*.

Recalling that the bending strain is the space derivative of the bending angle *ψ* (figure 13), we may write

*M*is the bending moment,

*E*the Young’s modulus, and

*I*

_{1}is the second area moment of inertia about the axis

**j**=

**k**×

**i**(figure 5

*a*). The shear angle

*θ*, as illustrated in figure 13, is the difference between the bending angle and the slope of the centreline, so that

*V*is the shear force,

*A*is the cross-sectional area,

*G*is the shear modulus and

*α*

_{c}is the Timoshenko factor [58]. If a point load

*F*is applied downward at

*s*=

*L*, where

*L*is the length of the rod, a free body diagram of the beam yields

*M*=−

*F*(

*L*−

*s*) and

*V*=

*F*, so that

By integrating equation (E 3) with boundary conditions *ψ*=0 at *s*=0, injecting the solution *ψ* into equation (E 4), and integrating again with boundary conditions *y*=0 at *s*=0, we obtain

## Appendix F. Further validations

**F.1. Euler buckling instability**

Euler buckling involves a single straight isotropic, inextensible and unshearable rod subject to an axial load *F*, as depicted in figure 14*a*. The critical axial load *F*_{c} that the rod can withstand before bending can be expressed analytically [3] as

*E*is the modulus of elasticity of the rod,

*I*=

*I*

_{1}=

*I*

_{2}is the area moment of inertia,

*L*is the length and

*b*is a constant which depends on the boundary conditions. If both ends are fixed in space but free to rotate then

*b*=1, thus

*α*=

*EI*is the bending stiffness.

We test our solver against the above analytical solution by simulating the time evolution of an inextensible and unshearable rod. In figure 14*a*, we show the results of our computations in the limit when both ends are free to rotate and all their spatial degrees of freedom are fixed except for one which allows the top end to slide vertically. The inextensible and unshearable conditions are enforced numerically by setting the entries of the matrix **S** to be much larger than those of **B** (details in figure 14).

We explore the phase spaces *F*-*α* and *F*-*L* and determine the probability of a randomly perturbed rod to undergo a buckling event, characterized by the bending energy (defined in the main text as a quadratic functional) exceeding the small threshold value *E*_{B}>*E*_{th}. As can be seen in figure 14*b*,*c*, the obtained probability landscapes exhibit a sharp transition in agreement with the exact solution of equation (F 2).

**F.2. Mitchell buckling instability**

The Euler buckling benchmark allows us to test the capability of our solver to capture the onset of an instability relative to the bending modes. Next, we consider the Mitchell buckling instability that simultaneously accounts for both bend and twist. When the ends of an isotropic, inextensible and unshearable filament are joined together, the resulting configuration is a planar circular shape. If the two ends are twisted by a given angle and glued together, a circular shape with uniform twist density is obtained. When the total twist *Φ*, i.e. the integrated twist line density along the filament, exceeds a critical value *Φ*_{c}, the rod buckles out of plane. This critical total twist can be expressed analytically [103] as

*α*and

*β*are, respectively, the bending and twist rigidities.

We explore the phase space *F*-(*β*/*α*) and determine the probability of a randomly perturbed rod to undergo a buckling event, characterized by the translational energy (defined in the main text as a quadratic functional) exceeding a small threshold value *E*_{T}>*E*_{th}. As can be seen in figure 15, our results show a sharp transition in agreement with the exact solution of equation (F 3).

**F.3. Twist forced vibrations in a stretched rod**

Next, we examine the problem of twisting waves generated in a rod axially stretched, as illustrated in figure 16*a*. The coupling between stretching and the other dynamic modes tests the rescaling in terms of the dilatation factor *e* of equations (2.14) and (2.15).

The rod is clamped at one end, stretched by the quantity $\mathrm{\Delta}L=(e-1)\hat{L}$, and forced to angularly vibrate about the axial direction by applying at the free end the couple ${A}_{v}\mathrm{sin}(2\pi {f}_{v}t)$, where *A*_{v} and *f*_{v} are the corresponding forcing amplitude and frequency. In the case of no internal dissipation and constant circular sections, the induced standing angular wave *θ*(*s*,*t*) can be determined analytically. The dynamics of twisting yields the wave equation for the angular displacement [104]

*G*is the shear modulus and

*ρ*is the density. By assuming a solution of the form $\theta (s,t)=\varphi (s)\mathrm{sin}(2\pi {f}_{v}t)$, and by applying the boundary conditions

*ϕ*(0)=0 and ${d}_{s}\varphi (0)={C}_{v}/({\hat{I}}_{3}G)$ with

*C*

_{v}twisting torque and ${\hat{I}}_{3}$ second area moment of inertia about the axial direction, we can solve equation (F 4) obtaining

As can be seen in figure 16*b*, our numerical method recovers the derived analytical solution for the twist angular displacement along the stretched rod. Moreover, the solver consistently exhibits a second-order convergence in time and space given the error metric *ϵ*=∥*θ*−*θ*^{n}∥, where *θ*^{n} represents our numerical solutions in the limit of refinement (figure 16*c*).

## Appendix G. Further friction validations

**G.1. Validation of rolling friction model**

Here, we validate the friction model introduced in §v on three test cases that can be analytically characterized. In all benchmarks, we consider a rigid, unshearable, inextensible straight rod of mass *m*, length *L*, radius *r*, and axial mass second moment of inertia *J*. The rod interacts with a surface characterized by static and kinetic friction coefficients *μ*_{s} and *μ*_{k}, thus assuming isotropic friction.

In the first test, the rod is initially at rest on a plane inclined of the angle *α*_{s}, as depicted in figure 17*a*. Owing to the gravitational acceleration *g* the rod starts rolling or slipping, depending on *α*_{s}, down the plane. The linear *v* and angular *ω* velocities of the filament, and therefore the corresponding translational *E*_{T}=*mv*^{2}/2 and rotational *E*_{R}=*Jω*^{2}/2 energies, can be analytically determined.

By recalling equation (4.12), the force necessary to ensure rolling without slip takes the form ${F}_{\text{noslip}}=-{F}_{\parallel}/3=-mg\mathrm{sin}({\alpha}_{\mathrm{s}})/3$. Given the maximum static friction force ${F}_{\mathrm{s}}={\mu}_{\mathrm{s}}{F}_{\perp}={\mu}_{\mathrm{s}}mg\mathrm{cos}({\alpha}_{\mathrm{s}})$ in the case |*F*_{noslip}|≤|*F*_{s}|, we have *a*=(*F*_{∥}+*F*_{noslip})/*m*. Therefore, at the time *T* after releasing the rod, linear and angular velocities read, respectively, *v*=*aT* and $\omega =\dot{\omega}T=(a/r)T$, expressing the no slip kinematic constraint between linear *a* and angular $\dot{\omega}$ accelerations. On the contrary, if |*F*_{noslip}|>|*F*_{s}| the rod starts slipping and linear and angular accelerations are no longer coupled, so that after time *T* we have *a*=(*F*_{∥}−*μ*_{k}*F*_{⊥})/*m*, *v*=*aT*, $\omega =\dot{\omega}T=({\mu}_{k}{F}_{\perp}r/J)T$. Therefore, translational and rotational energies as a function of *α*_{s} finally read

*b*, our numerical approach faithfully reproduces the derived analytical solution, accurately capturing the discontinuity at the transition between pure rolling and slipping.

The second test case of figure 17*c*,*d* entails a rod set in motion by the external couple *C*_{s} on a horizontal plane. Depending on the magnitude of the load the filament exhibits pure rolling or slipping motion. By recalling equation (4.12), the force necessary to ensure rolling without slip takes the form *F*_{noslip}=2*C*_{s}/(3*r*). Given the maximum static friction force *F*_{s}=*μ*_{s}*F*_{⊥}=*μ*_{s}*mg*, in the case of |*F*_{noslip}|≤|*F*_{s}| we have *a*=*F*_{noslip}/*m*. Therefore, at the time *T* after releasing the rod, linear and angular velocities read, respectively, *v*=*aT* and $\omega =\dot{\omega}T=(a/r)T$, expressing the no slip kinematic constraint between linear *a* and angular $\dot{\omega}$ accelerations. On the contrary, if |*F*_{noslip}|>|*F*_{s}| the rod starts slipping and linear and angular accelerations are no longer coupled, so that after time *T* we have *a*=*μ*_{k}*F*_{⊥}/*m*, *v*=*aT*, $\omega =\dot{\omega}T={J}^{-1}({C}_{\mathrm{s}}-{\mu}_{\mathrm{k}}{F}_{\perp}r)T$. Therefore, translational and rotational energies as a function of *C*_{s} finally read

*d*, again our numerical approach faithfully reproduces the derived analytical solution, accurately capturing the discontinuity at the transition between pure rolling and slipping.

In the last test case, we consider a rod of initial velocity *V* _{s} in the direction parallel to a horizontal plane and perpendicular to the filament axis (figure 17*e*), and we vary the relative mass moment of inertia ratio λ=2*J*/(*mr*^{2}). Initially, the rod slips on the surface and gradually slows down, due to the effect of the kinetic friction, to a point for which linear *a* and angular velocity *ω* meet the kinematic constraint *v*_{eq}=*rω*_{eq} characteristic of pure rolling motion. The frictional force *F* induces the torque *M*=*rF*, so that over a period of time *T* we have *v*=*V* _{s}−*FT*/*m* and *ω*=*rFT*/*J*. By enforcing the no slip kinematic constraint *V* _{s}−*FT*/*m*=*r*^{2}*FT*/*J*, we have that *FT*=*V* _{s}/(*r*^{2}/*J*+1/*m*). This allows us to directly compute *v*_{eq} and *ω*_{eq} and the corresponding translational and rotational energies, yielding

*f*, our numerical approach accurately reproduces the predicted energies as a function of λ.

**G.2. Validation of anisotropic longitudinal friction model**

After validating our rolling friction model, we turn to its longitudinal counterpart. Here, we consider a straight, rigid, inextensible and unshearable rod which is axially pulled or pushed with force **F** for a fixed period of time *T* (figure 18). Anisotropy is captured through the forward *μ*^{f}_{s}, ${\mu}_{\mathrm{k}}^{\mathrm{f}}$ and backward *μ*^{b}_{s}, ${\mu}_{\mathrm{k}}^{\mathrm{b}}$ static and kinetic coefficients. For **F**⋅**t**≥0, the rod translational energy takes the form

**F**⋅

**t**<0 we have

**G.3. Isotropic versus anisotropic friction driven locomotion**

In this section, we illustrate the effect of symmetry breaking via anisotropic friction in a system constituted by a soft filament interacting via surface friction with a solid substrate. If we consider isotropic friction and a specular muscular activation pattern by setting the control values *β*_{1}=*β*_{4} and *β*_{2}=*β*_{3} and the wave number 2*π*/λ_{m}=0 (see §ii), then the system is symmetric (up to inertial effects that can be neglected for small Froude numbers). Therefore, over any activation cycle the snake centre of mass does not move. On the contrary, if we introduce anisotropy the snake will be able to slowly move (the capability to effectively move is impaired by the absence of the travelling gait component since 2*π*/λ_{m}=0).

As can be observed in figure 19, this prediction is captured by our numerical scheme, which accurately resolves the physical mechanisms at play. Moreover, this test shows once again how our methodology is robust in terms of numerical noise as no spurious displacements or rotations are generated in the symmetric case (figure 19*a*,*b*).

## Appendix H. Representative timings

The proposed algorithm strikes a balance between computing costs, numerical accuracy and implementation modularity: by foregoing an implicit integration scheme we can incorporate a number of additional physical effects and soft constraints, even though this may come at the expense of computational efficiency. Indeed, for large Young’s or shear moduli or for very thin rods, the system of equations (2.4)–(2.7) in the main text might become stiff, so that small timesteps must be employed to ensure stability. Although we have not derived any rigorous CFL condition, throughout this work we employed the empirical relation d*t*=*a* dℓ, with *a*∼10^{−2} s m^{−1}, and found it reliable in preventing numerical instabilities (figure 20).

Despite the potential stiffness of this model, since the computational cost per timestep scales linearly with the number of discretization elements *n* (the model is one-dimensional), we could carry out all our computation on a single MacBook Pro equipped with a 2.3 GHz dual-core Intel Core i5, Turbo Boost up to 3.6 GHz, with 64 MB of eDRAM. Because of the condition d*t*=*a* dℓ, a smaller number of elements implies larger timesteps (dℓ=*L*/*n*). Therefore, the time-to-solution may scale between linearly and quadratically, depending on how the timestep is set. This can be noted by inspecting Euler and helical buckling studies reported in figure 20. In the Euler buckling simulations, the timestep is set to be the same for both case 1 and case 2, and indeed by halving *n*, the time-to-solution and the time per timestep reduce by a factor of 2. In the helical buckling study, the timestep was adaptively selected according to d*t*=*a* dℓ, so that, when *n* is halved, the time-to-solution reduces by a factor of 4, while the time per timestep reduces by a factor of 2.

### References

- 1
Kirchhoff G . 1859 Ueber das gleichgewicht und die bewegung eines unendlich dünnen elastischen stabes.**J. Reine Angew. Math.**, 285–313. (doi:10.1515/crll.1859.56.285) Crossref, Google Scholar**56** - 2
- 3
Love AEH . 1906**A treatise on the mathematical theory of elasticity**. Cambridge, UK: Cambridge University Press. Google Scholar - 4
Cosserat E, Cosserat F . 1909**Théorie des corps déformables**. Ithaca, NY: Cornell University Library. Google Scholar - 5
- 6
Harold DE . 1992 Kirchhoff’s theory of rods.**Arch. Hist. Exact. Sci.**, 1–23. (doi:10.1007/BF00379680) Crossref, ISI, Google Scholar**44** - 7
Langer J, Singer DA . 1996 Lagrangian aspects of the kirchhoff elastic rod.**SIAM Rev.**, 605–618. (doi:10.1137/S0036144593253290) Crossref, ISI, Google Scholar**38** - 8
van der Heijden GHM, Thompson JMT . 2000 Helical and localised buckling in twisted rods: a unified analysis of the symmetric case.**Nonlinear. Dyn.**, 71–99. (doi:10.1023/A:1008310425967) Crossref, ISI, Google Scholar**21** - 9
Goss VGA, Van der Heijden GHM, Thompson JMT, Neukirch S . 2005 Experiments on snap buckling, hysteresis and loop formation in twisted rods.**Exp. Mech.**, 101–111. (doi:10.1007/BF02428182) Crossref, ISI, Google Scholar**45** - 10
Maugin GA, Metrikine AV . 2010**Mechanics of generalized continua**.Advances in mechanics and mathematics, vol. 21 . Berlin, Germany: Springer. Google Scholar - 11
Goss VGA . 2009 The history of the planar elastica: insights into mechanics and scientific method.**Sci. Educ.**, 1057–1082. (doi:10.1007/s11191-008-9166-2) Crossref, ISI, Google Scholar**18** - 12
Yang Y, Tobias I, Olson WK . 1993 Finite element analysis of DNA supercoiling.**J. Chem. Phys.**, 1673–1686. (doi:10.1063/1.464283) Crossref, ISI, Google Scholar**98** - 13
Spillmann J, Teschner M . 2007 Corde: Cosserat rod elements for the dynamic simulation of one-dimensional elastic objects. In*Eurographics/SIGGRAPH Symposium on Computer Animation, San Diego, California*, pp. 63–72. Geneva, Switzerland: Eurographics Association. (doi:10.2312/SCA/SCA07/063-072) Google Scholar - 14
Bergou M, Wardetzky M, Robinson S, Audoly B, Grinspun E . 2008 Discrete elastic rods.**ACM Trans. Graph. (TOG)**, 1–12. Crossref, ISI, Google Scholar**27** - 15
Arne W, Marheineke N, Meister A, Wegener R . 2010 Numerical analysis of cosserat rod and string models for viscous jets in rotational spinning processes.**Math. Models Methods Appl. Sci.**, 1941–1965. (doi:10.1142/S0218202510004738) Crossref, ISI, Google Scholar**20** - 16
Audoly B, Clauvelin N, Brun PT, Bergou M, Grinspun E, Wardetzky M . 2013 A discrete geometric approach for simulating the dynamics of thin viscous threads.**J. Comput. Phys.**, 18–49. (doi:10.1016/j.jcp.2013.06.034) Crossref, ISI, Google Scholar**253** - 17
Wolgemuth CW, Powers TR, Goldstein RE . 2000 Twirling and whirling: viscous dynamics of rotating elastic filaments.**Phys. Rev. Lett.**, 1623–1626. (doi:10.1103/PhysRevLett.84.1623) Crossref, PubMed, ISI, Google Scholar**84** - 18
Fütterer T, Klar A, Wegener R . 2009 An energy conserving numerical scheme for the dynamics of hyperelastic rods.**Int. J. Differ. Equ.**, 718308. (doi:10.1155/2012/718308) Google Scholar**2012** - 19
Bertails F, Audoly B, Cani MP, Querleux B, Leroy F, Leveque JL . 2006 Super-helices for predicting the dynamics of natural hair.**ACM. Trans. Graph.**, 1180–1187. (doi:10.1145/1141911.1142012) Crossref, ISI, Google Scholar**25** - 20
Goriely A, Tabor M . 1998 Spontaneous helix hand reversal and tendril perversion in climbing plants.**Phys. Rev. Lett.**, 1564–1567. (doi:10.1103/PhysRevLett.80.1564) Crossref, ISI, Google Scholar**80** - 21
Gerbode SJ, Puzey JR, McCormick AG, Mahadevan L . 2012 How the cucumber tendril coils and overwinds.**Science**, 1087–1091. (doi:10.1126/science.1223304) Crossref, PubMed, ISI, Google Scholar**337** - 22
Kaldor JM, James DL, Marschner S . 2008 Simulating knitted cloth at the yarn level.**ACM Trans. Graph. (TOG)**, 65. (doi:10.1145/1360612.1360664) Crossref, ISI, Google Scholar**27** - 23
Ward K, Bertails F, Kim TY, Marschner SR, Cani MP, Lin MC . 2007 A survey on hair modeling: styling, simulation, and rendering.**IEEE Trans. Vis. Comput. Graph.**, 213–234. (doi:10.1109/TVCG.2007.30) Crossref, PubMed, ISI, Google Scholar**13** - 24
Daviet G, Bertails-Descoubes F, Boissieux L . 2011 A hybrid iterative solver for robustly capturing coulomb friction in hair dynamics.**ACM. Trans. Graph.**, 139. (doi:10.1145/2070781.2024173) Crossref, ISI, Google Scholar**30** - 25
Durville D . 2005 Numerical simulation of entangled materials mechanical properties.**J. Mater. Sci.**, 5941–5948. (doi:10.1007/s10853-005-5061-2) Crossref, ISI, Google Scholar**40** - 26
Bergou M, Audoly B, Vouga E, Wardetzky M, Grinspun E . 2010 Discrete viscous threads.**ACM Trans. Graph. (TOG)**, 116. (doi:10.1145/1778765.1778853) Crossref, ISI, Google Scholar**29** - 27
Brun PT, Audoly B, Ribe NM, Eaves TS, Lister JR . 2015 Liquid ropes: a geometrical model for thin viscous jet instabilities.**Phys. Rev. Lett.**, 174501. (doi:10.1103/PhysRevLett.114.174501) Crossref, PubMed, ISI, Google Scholar**114** - 28
Brun PT, Ribe NM, Audoly B . 2012 A numerical investigation of the fluid mechanical sewing machine.**Phys. Fluids**, 043102. (doi:10.1063/1.3703316) Crossref, ISI, Google Scholar**24** - 29
Ribe NM . 2003 Periodic folding of viscous sheets.**Phys. Rev. E**, 036305. (doi:10.1103/PhysRevE.68.036305) Crossref, ISI, Google Scholar**68** - 30
Ribe NM, Habibi M, Bonn D . 2012 Liquid rope coiling.**Annu. Rev. Fluid. Mech.**, 249–266. (doi:10.1146/annurev-fluid-120710-101244) Crossref, ISI, Google Scholar**44** - 31
Ko W, Lim S, Lee W, Kim Y, Berg HC, Peskin CS . 2017 Modeling polymorphic transformation of rotating bacterial flagella in a viscous fluid.**Phys. Rev. E**, 063106. (doi:10.1103/PhysRevE.95.063106) Crossref, PubMed, ISI, Google Scholar**95** - 32
Lee W, Kim Y, Olson SD, Lim S . 2014 Nonlinear dynamics of a rotating elastic rod in a viscous fluid.**Phys. Rev. E**, 033012. (doi:10.1103/PhysRevE.90.033012) Crossref, ISI, Google Scholar**90** - 33
Olson SD, Lim S, Cortez R . 2013 Modeling the dynamics of an elastic rod with intrinsic curvature and twist using a regularized stokes formulation.**J. Comput. Phys.**, 169–187. (doi:10.1016/j.jcp.2012.12.026) Crossref, ISI, Google Scholar**238** - 34
Lim S, Peskin CS . 2012 Fluid-mechanical interaction of flexible bacterial flagella by the immersed boundary method.**Phys. Rev. E**, 036307. (doi:10.1103/PhysRevE.85.036307) Crossref, ISI, Google Scholar**85** - 35
Lim S . 2010 Dynamics of an open elastic rod with intrinsic curvature and twist in a viscous fluid.**Phys. Fluids**, 024104. (doi:10.1063/1.3326075) Crossref, ISI, Google Scholar**22** - 36
Lim S, Ferent A, Wang XS, Peskin CS . 2008 Dynamics of a closed rod with twist and bend in fluid.**SIAM. J. Sci. Comput.**, 273–302. (doi:10.1137/070699780) Crossref, ISI, Google Scholar**31** - 37
Linn J . 2016 Discrete kinematics of cosserat rods based on the difference geometry of framed curves. In*The 4th Joint Int. Conf. on Multibody System Dynamics, Montreal, Canada*. Berlin, Germany: Springer. Google Scholar - 38
Linn J, Hermansson T, Andersson F, Schneider F . 2017 Kinetic aspects of discrete cosserat rods based on the difference geometry of framed curves. In*ECCOMAS Thematic Conference on Multibody Dynamics, Prague, Czech Republic*. Berlin, Germany: Springer. Google Scholar - 39
Sonneville V, Cardona A, Bruls O . 2014 Geometrically exact beam finite element formulated on the special Euclidean group SE (3).**Comput. Methods. Appl. Mech. Eng.**, 451–474. (doi:10.1016/j.cma.2013.10.008) Crossref, ISI, Google Scholar**268** - 40
Ghatak A, Mahadevan L . 2005 Solenoids and plectonemes in stretched and twisted elastomeric filaments.**Phys. Rev. Lett.**, 057801. (doi:10.1103/PhysRevLett.95.057801) Crossref, PubMed, ISI, Google Scholar**95** - 41
Haines CS *et al.*2014 Artificial muscles from fishing line and sewing thread.**Science**, 868–872. (doi:10.1126/science.1246906) Crossref, PubMed, ISI, Google Scholar**343** - 42
Haines CS, Li N, Spinks GM, Aliev AE, Di J, Baughman RH . 2016 New twist on artificial muscles.**Proc. Natl Acad. Sci. USA**, 11 709–11 716. (doi:10.1073/pnas.1605273113) Crossref, ISI, Google Scholar**113** - 43
Kim SH *et al.*2017 Harvesting electrical energy from carbon nanotube yarn twist.**Science**, 773–778. (doi:10.1126/science.aam8771) Crossref, PubMed, ISI, Google Scholar**357** - 44
Kim S, Laschi C, Trimmer B . 2013 Soft robotics: a bioinspired evolution in robotics.**Trends. Biotechnol.**, 287–294. (doi:10.1016/j.tibtech.2013.03.002) Crossref, PubMed, ISI, Google Scholar**31** - 45
Park SJ *et al.*2016 Phototactic guidance of a tissue-engineered soft-robotic ray.**Science**, 158–162. (doi:10.1126/science.aaf4292) Crossref, PubMed, ISI, Google Scholar**353** - 46
Sadati SMH, Naghibi SE, Shiva A, Noh Y, Gupta A, Walker ID, Althoefer K, Nanayakkara T . 2017 A geometry deformation model for braided continuum manipulators.**Front. Robot. AI**, 22. (doi:10.3389/frobt.2017.00022) Crossref, Google Scholar**4** - 47
Sadati SMH, Shiva A, Ataka A, Naghibi SE, Walker ID, Althoefer K, Nanayakkara T . 2016 A geometry deformation model for compound continuum manipulators with external loading.*2016 IEEE International Conference on Robotics and Automation (ICRA), Stockholm, Sweden*, pp. 4957–4962. Piscataway, New Jersey: Institute of Electrical and Electronics Engineers. Google Scholar - 48
Trivedi D, Rahn CD, Kier WM, Walker ID . 2008 Soft robotics: biological inspiration, state of the art, and future research.**Appl. Bionics. Biomech.**, 99–117. (doi:10.1080/11762320802557865) Crossref, Google Scholar**5** - 49
Armanini C, DalCorso F, Misseroni D, Bigoni D . 2017 From the elastica compass to the elastica catapult: an essay on the mechanics of soft robot arm.**Proc. R. Soc. A**, 20160870. (doi:10.1098/rspa.2016.0870) Link, Google Scholar**473** - 50
Sadati SMH, Naghibi SE, Shiva A, Walker ID, Althoefer K, Nanayakkara T . 2017 Mechanics of continuum manipulators, a comparative study of five methods with experiments. In*Conference Towards Autonomous Robotic Systems, Guildford, UK*, pp. 686–702. Berlin, Germany: Springer. Google Scholar - 51
Connolly F, Walsh CJ, Bertoldi K . 2017 Automatic design of fiber-reinforced soft actuators for trajectory matching.**Proc. Natl Acad. Sci. USA**, 51–56. (doi:10.1073/pnas.1615140114) Crossref, PubMed, ISI, Google Scholar**114** - 52
Gazzola M, Burckhardt CJ, Bayati B, Engelke M, Greber UF, Koumoutsakos P . 2009 A stochastic model for microtubule motors describes the*in vivo*cytoplasmic transport of human adenovirus.**PLoS. Comput. Biol.**, e1000623. (doi:/10.1371/journal.pcbi.1000623) Crossref, PubMed, ISI, Google Scholar**5** - 53
Gazzola M, van Rees WM, Koumoutsakos P . 2012 C-start: optimal start of larval fish.**J. Fluid. Mech.**, 5–18. (doi:10.1017/jfm.2011.558) Crossref, ISI, Google Scholar**698** - 54
Gazzola M, Hejazialhosseini B, Koumoutsakos P . 2014 Reinforcement learning and wavelet adapted vortex methods for simulations of self-propelled swimmers.**SIAM. J. Sci. Comput.**, B622–B639. (doi:10.1137/130943078) Crossref, ISI, Google Scholar**36** - 55
Gazzola M, Tchieu AA, Alexeev D, de Brauer A, Koumoutsakos P . 2015 Learning to school in the presence of hydrodynamic interactions.**J. Fluid. Mech.**, 726–749. (doi:10.1017/jfm.2015.686) Crossref, ISI, Google Scholar**789** - 56
van Rees WM, Gazzola M, Koumoutsakos P . 2013 Optimal shapes for anguilliform swimmers at intermediate reynolds numbers.**J. Fluid. Mech.**, R3. (doi:10.1017/jfm.2013.157) Crossref, ISI, Google Scholar**722** - 57
van Rees WM, Gazzola M, Koumoutsakos P . 2015 Optimal morphokinematics for undulatory swimmers at intermediate reynolds numbers.**J. Fluid. Mech.**, 178–188. (doi:10.1017/jfm.2015.283) Crossref, ISI, Google Scholar**775** - 58
- 59
Landau LD, Bell JS, Kearsley MJ, Pitaevskii LP, Lifshitz EM, Sykes JB . 1984**Electrodynamics of continuous media**,**vol. 8**. Amsterdam: The Netherlands: Elsevier. Google Scholar - 60
- 61
Mooney M . 1940 A theory of large elastic deformation.**J. Appl. Phys.**, 582–592. (doi:10.1063/1.1712836) Crossref, Google Scholar**11** - 62
Rivlin RS . 1948 Large elastic deformations of isotropic materials iv. further developments of the general theory.**Phil. Trans. R. Soc. Lond. A**, 379–397. (doi:10.1098/rsta.1948.0024) Link, ISI, Google Scholar**241** - 63
Baillieul J, Levi M . 1987 Rotational elastic dynamics.**Phys. D: Nonlinear Phenom.**, 43–62. (doi:10.1016/0167-2789(87)90004-2) Crossref, ISI, Google Scholar**27** - 64
Levi M . 1996 Composition of rotations and parallel transport.**Nonlinearity**, 413–419. (doi:10.1088/0951-7715/9/2/007) Crossref, ISI, Google Scholar**9** - 65
Goriely A, Tabor M . 2000 The nonlinear dynamics of filaments.**Nonlinear. Dyn.**, 101–133. (doi:10.1023/A:1008366526875) Crossref, ISI, Google Scholar**21** - 66
Rodrigues O . 1840 Des lois géométriques qui régissent les déplacements d’un système solide dans l’espace: et de la variation des cordonnées provenant de ces déplacements considérés indépendamment des causes qui peuvent les produire.**Journal de Mathématique Pures et Appliquées**, 380–440. Google Scholar**5** - 67
Sobottka G, Lay T, Weber A . 2008 Stable integration of the dynamic cosserat equations with application to hair modeling.**Journal of WSCG**, 73–80. Google Scholar**16** - 68
Grinspun E, Desbrun M, Polthier K, Schroder P, Stern A . 2006 Discrete differential geometry: an applied introduction.**SIGGRAPH Course Notes**, 1–80. Google Scholar**7** - 69
Lang H, Linn J, Arnold M . 2011 Multi-body dynamics simulation of geometrically exact cosserat rods.**Multibody. Syst. Dyn.**, 285–312. (doi:10.1007/s11044-010-9223-x) Crossref, ISI, Google Scholar**25** - 70
McCormick A . 2013 Discrete differential geometry and physics of elastic curves. PhD thesis, Harvard University. Google Scholar - 71
Ge Z, Kruse HP, Marsden JE . 1996 The limits of hamiltonian structures in three-dimensional elasticity, shells, and rods.**J. Nonlinear Sci.**, 19–57. (doi:10.1007/BF02433809) Crossref, ISI, Google Scholar**6** - 72
Kol A, Laird BB, Leimkuhler BJ . 1997 A symplectic method for rigid-body molecular simulation.**J. Chem. Phys.**, 2580–2588. (doi:10.1063/1.474596) Crossref, ISI, Google Scholar**107** - 73
Coyne J . 1990 Analysis of the formation and elimination of loops in twisted cable.**IEEE J. Ocean. Eng.**, 72–83. (doi:10.1109/48.50692) Crossref, ISI, Google Scholar**15** - 74
Neukirch S, van der Heijden GHM, Thompson JMT . 2002 Writhing instabilities of twisted rods: from infinite to finite length.**J. Mech. Phys. Solids**, 1175–1191. (doi:10.1016/S0022-5096(01)00130-2) Crossref, ISI, Google Scholar**50** - 75
van der Heijden GHM, Neukirch S, Goss VGA, Thompson JMT . 2003 Instability and self-contact phenomena in the writhing of clamped rods.**Int. J. Mech. Sci.**, 161–196. (doi:10.1016/S0020-7403(02)00183-2) Crossref, ISI, Google Scholar**45** - 76
Galloni EE, Kohen M . 1979 Influence of the mass of the spring on its static and dynamic effects.**Am. J. Phys.**, 1076–1078. (doi:10.1119/1.11978) Crossref, ISI, Google Scholar**47** - 77
- 78
Torby BJ . 1984**Advanced dynamics for engineers**. New York, NY: Holt Rinehart and Winston. Google Scholar - 79
Altringham JD, John D, Johnston IA . 1990 Modelling muscle power output in a swimming fish.**J. Exp. Biol.**, 395–402. Crossref, ISI, Google Scholar**148** - 80
Gazzola M, Argentina M, Mahadevan L . 2015 Gait and speed selection in slender inertial swimmers.**Proc. Natl Acad. Sci. USA**, 3874–3879. (doi:10.1073/pnas.1419335112) Crossref, PubMed, ISI, Google Scholar**112** - 81
Riewe F . 1996 Nonconservative lagrangian and hamiltonian mechanics.**Phys. Rev. E**, 1890–1899. (doi:10.1103/PhysRevE.53.1890) Crossref, ISI, Google Scholar**53** - 82
Guo ZV, Mahadevan L . 2008 Limbless undulatory propulsion on land.**Proc. Natl Acad. Sci. USA**, 3179–3184. (doi:10.1073/pnas.0705442105) Crossref, PubMed, ISI, Google Scholar**105** - 83
Hu DL, Nirody J, Scott T, Shelley MJ . 2009 The mechanics of slithering locomotion.**Proc. Natl Acad. Sci. USA**, 10 081–10 085. (doi:10.1073/pnas.0812533106) Crossref, ISI, Google Scholar**106** - 84
Cox RG . 1970 The motion of long slender bodies in a viscous fluid part 1. General theory.**J. Fluid. Mech.**, 791–810. (doi:10.1017/S002211207000215X) Crossref, Google Scholar**44** - 85
Lauga E, Powers TR . 2009 The hydrodynamics of swimming microorganisms.**Rep. Prog. Phys.**, 096601. (doi:10.1088/0034-4885/72/9/096601) Crossref, ISI, Google Scholar**72** - 86
Williams BJ, Anand SV, Rajagopalan J, Saif MA . 2014 A self-propelled biohybrid swimmer at low reynolds number.**Nat. Commun.**, 3081. (doi:10.1038/ncomms4081) Crossref, PubMed, ISI, Google Scholar**5** - 87
Thompson JMT, Champneys AR . 1996 From helix to localized writhing in the torsional post-buckling of elastic rods.**Proc. R. Soc. Lond. A**, 117–138. (doi:10.1098/rspa.1996.0007) Link, ISI, Google Scholar**452** - 88
Gray J, Lissmann HW . 1950 The kinetics of locomotion of the grass-snake.**J. Exp. Biol.**, 354–367. Crossref, ISI, Google Scholar**26** - 89
Transeth AA, Pettersen KY, Liljebäck P . 2009 A survey on snake robot modeling and locomotion.**Robotica**, 999–1015. (doi:10.1017/S0263574709005414) Crossref, ISI, Google Scholar**27** - 90
Gray J . 1946 The mechanism of locomotion in snakes.**J. Exp. Biol.**, 101–120. Crossref, PubMed, ISI, Google Scholar**23** - 91
Alben S . 2013 Optimizing snake locomotion in the plane.**Proc. R. Soc. A**, 20130236. (doi:10.1098/rspa.2013.0236) Link, Google Scholar**469** - 92
Erkmen I, Erkmen AM, Matsuno F, Chatterjee R, Kamegawa T . 2002 Snake robots to the rescue!**IEEE Robot. Autom. Mag.**, 17–25. (doi:10.1109/MRA.2002.1035210) Crossref, ISI, Google Scholar**9** - 93
Tanev I, Ray T, Buller A . 2005 Automated evolutionary design, robustness, and adaptation of sidewinding locomotion of a simulated snake-like robot.**IEEE. Trans. Robot.**, 632–645. (doi:10.1109/TRO.2005.851028) Crossref, ISI, Google Scholar**21** - 94
Hansen N, Ostermeier A . 2001 Completely derandomized self-adaptation in evolution strategies.**Evol. Comput.**, 159–195. (doi:10.1162/106365601750190398) Crossref, PubMed, ISI, Google Scholar**9** - 95
Hansen N, Muller SD, Koumoutsakos P . 2003 Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES).**Evol. Comput.**, 1–18. (doi:10.1162/106365603321828970) Crossref, PubMed, ISI, Google Scholar**11** - 96
Chatelain P, Gazzola M, Kern S, Koumoutsakos P . 2011**Optimization of aircraft wake alleviation schemes through an evolution strategy**.Lecture Notes in Computer Science, vol. 6649 , pp. 210–221. Google Scholar - 97
Shine R, Cogger HG, Reed RR, Shetty S, Bonnet X . 2003 Aquatic and terrestrial locomotor speeds of amphibious sea-snakes (serpentes, laticaudidae).**J. Zool.**, 261–268. (doi:10.1017/S0952836902003242) Crossref, ISI, Google Scholar**259** - 98
- 99
Woolley DM, Vernon GG . 2001 A study of helical and planar waves on sea urchin sperm flagella, with a theory of how they are generated.**J. Exp. Biol.**, 1333–1345. Crossref, PubMed, ISI, Google Scholar**204** - 100
Gazzola M, Chatelain P, van Rees WM, Koumoutsakos P . 2011 Simulations of single and multiple swimmers with non-divergence free deforming geometries.**J. Comput. Phys.**, 7093–7114. (doi:10.1016/j.jcp.2011.04.025) Crossref, ISI, Google Scholar**230** - 101
Pagan-Diaz GJ *et al.*2018 Simulation and fabrication of stronger, larger and faster walking biohybrid machines.**Advanced Functional Materials (advance online publication)**1801145. (doi:10.1002/adfm.201801145) Crossref, ISI, Google Scholar - 102
- 103
Goriely A . 2006 Twisted elastic rings and the rediscoveries of michell’s instability.**J. Elast.**, 281–299. (doi:10.1007/s10659-006-9055-3) Crossref, ISI, Google Scholar**84** - 104
Selvadurai APS . 2000**Partial differential equations in mechanics**. Berlin, Germany: Springer Science and Business Media. Google Scholar