## 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 **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 **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 **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

We will first start by presenting the equations of motion under the assumption of inextensibility (i.e. ** ω**=

*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,

**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 (*κ*_{1},*κ*_{2}) and twist (*κ*_{3}) strains in the material frame (table 1).

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 ** κ**.

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 **t**=**d**_{3} [5]. This implies that *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**,

### 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 **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

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

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 **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 *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 *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 **u**
*vice versa*

Conversely, given a rotation matrix **R**, the corresponding rotation vector can be directly computed via the matrix logarithm operator

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

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}|, **t**_{i} is the discrete tangent vector,

Whereas in the continuum setting (§2) all quantities are defined pointwise, in a discrete setting some quantities, and in particular *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

Then, we may discretize the governing equations (2.12)–(2.15) so that they read
^{h} and *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

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 *φ* 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 *ϵ*=∥*φ*−*φ*^{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 *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

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 *ξ*≃*π*^{2}/4 for

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 *b*,*c*,*e*,*f* shows how the proposed numerical method converges to the analytical oscillation period *T** and normalized longitudinal displacement *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 *F* at the free end *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 *x* along the direction **k** can be approximated by the arc-length *a* and the appendix for further details and derivation), hence the use of

If the shear modulus *G* approaches infinity or if the ratio

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*^{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*^{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 (**B** has the finite value

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 *internal* physical and biophysical effects are captured by adding their contributions directly to the internal force

### 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 **f**_{v} and torques *γ*_{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 *N*_{m} control points

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

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 *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 *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}_{⊥}, **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 **v**_{∥}) and rotational (*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 **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 *μ*^{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 (*b*). When the load *F* is increased so that the elongation of the rod approaches *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 *μ*^{f}_{k} is set so that the ratio between inertial and friction forces captured by the Froude number is

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 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 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

**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 **Q** and noting that

## 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 *α* is the stress, *σ*=*e*−1), *C*_{10} is a material property to be determined experimentally and that in the limit of small strains (*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* 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 ^{h} and *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 *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 *ϕ*(0)=0 and *C*_{v} twisting torque and

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*_{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 *a* and angular *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*, *α*_{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 *a* and angular *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*, *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}, *μ*^{b}_{s}, **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.

## Footnotes

Electronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.4111496

- Received October 13, 2017.
- Accepted April 30, 2018.

- © 2018 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.