Royal Society Open Science
Open AccessResearch article

Forward and inverse problems in the mechanics of soft filaments

M. Gazzola

M. Gazzola

Department of Mechanical Science and Engineering, and National Center for Supercomputing Applications, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA

Google Scholar

Find this author on PubMed

,
L. H. Dudte

L. H. Dudte

John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA

Google Scholar

Find this author on PubMed

, and
L. Mahadevan

L. Mahadevan

John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA

Department of Physics, Harvard University, Cambridge, MA 02138, USA

[email protected]

Google Scholar

Find this author on PubMed

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 [110], 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 [1216] 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,2325]. Because the scaled ratio of the stretching and shearing stiffness to the bending stiffness for slender filaments is L2/r2≫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,2630]. 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 [3136]. 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 [4143]. 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 [4649] 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 [5257]. 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 (Lr) 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 r:(s[0,L]R,tR+)R3 and an oriented frame of reference Q:(s[0,L]R,tR+)SO(3) equivalent to the orthonormal triad of unit vectors Q={d1,d2,d3}. Here, s is the centre line arc-length coordinate in its current configuration and t is the time.

Figure 1.

Figure 1. The Cosserat rod model. A filament deforming in the three-dimensional space is represented by a centre-line coordinate r and a material frame characterized by the orthonormal triad {d1,d2,d3}. The corresponding orthogonal rotation matrix Q with row entries d1, d2, d3 transforms a vector x from the laboratory canonical basis {i,j,k} to the material frame of reference {d1,d2,d3} so that xL=Qx and vice versa x=QTxL. If extension or compression is allowed, the current filament configuration arc-length s may no longer coincide with the rest reference arc-length s^. This is captured via the scalar dilatation field e=ds/ds^. Moreover, to account for shear we allow the triad {d1,d2,d3} to detach from the unit tangent vector t so that d3t (we recall that the condition d3=t and e=1 correspond to the Kirchhoff constraint for unshearable and inextensible rods, and implies that σ=etd3=0). The dynamics of the centre line and the material frame is determined at each cross section by the internal force and torque resultants, n and τ. External loads are represented via the force f and couple c line densities.

Denoting by x any generic vector represented in the Eulerian frame and xL as the body-convected (Lagrangian) frame of reference allows us to write

laboratory:x=x¯1i+x¯2j+x¯3k2.1
and
body-convected:xL=x1d1+x2d2+x3d3,2.2
where equation (2.1) expresses x in the laboratory canonical basis {i,j,k}, while equation (2.2) expresses the same vector in the body-convected director basis {d1,d2,d3}. Then, the matrix Q transforms any vector x from the laboratory to the body-convected representation via xL=Qx and, conversely, x=Q1xL=QTxL, because QTQ=QQT=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 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=s^, Kirchhoff model), before generalizing them to the stretchable case in the subsequent sections. Denoting the rod angular velocity as ω=vec[(∂Q/∂t)TQ] and the generalized curvature as κ=vec[(∂Q/∂s)TQ], where vec[A] denotes the 3-vector associated with the skew-symmetric matrix A, the following transport identities hold

Qxt=xLt+ωL×xLandQxs=xLs+κL×xL.2.3
Using the above equations (full derivation in the appendix), we can express the advection of the rod positions and local frames, as well as the linear and angular momentum balance in a convenient Eulerian–Lagrangian form
rt=v, 2.4
djt=(QTωL)×dj,j=1,2,3 2.5
(ρAv)t=(QTnL)s+f 2.6
and(ρIωL)t=τLs+κL×τL+Qrs×nL+(ρIωL)×ωL+cL, 2.7
where ρ is the constant material density, A is the cross-sectional area, v is the velocity, nL and τ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 {d1,d2,d3} and are characterized by the generalized curvature. Specifically, the components of the curvature projected along the directors (κL=κ1d1+κ2d2+κ3d3) coincide with bending (κ1,κ2) and twist (κ3) strains in the material frame (table 1).

Table 1.Constitutive laws. The generalized curvature κL is associated with the bending κ1, κ2 about the principal directions (d1, d2) and the twist κ3 about the longitudinal one (d3), while σL=Q(etd3) is associated with the shears σ1, σ2 along the principal directions (d1, d2) and the axial extensional or compression σ3 along the longitudinal one (d3). The material properties of the rod are captured through the Young’s (E) and shear (G) moduli, while its geometric properties are accounted for via the cross-sectional area A, the second moment of inertia I and the constant αc=43 for circular cross sections [58]. The diagonal entries of the bending/twist BR3×3 and shear/stretch SR3×3 matrices are, respectively, (B1, B2, B3) and (S1, S2, S3). Pre-strains are modelled via the intrinsic curvature/twist κLo and shear/stretch σLo.

deformation modes strains rigidities loads
bending about d1 κ1 B1=EI1 τ1=B1(κ1κ1o)
bending about d2 κ2 B2=EI2 τ2=B2(κ2κ2o)
twist about d3 κ3 B3=GI3 τ3=B3(κ3κ3o)
shear along d1 σ1 S1=αcGA n1=S1(σ1σ1o)
shear along d2 σ2 S2=αcGA n2=S2(σ2σ2o)
stretch along d3 σ3 S3=EA n3=S3(σ3σ3o)

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

τL=B(κLκLo),2.8
where BR3×3=diag(B1,B2,B3) is the bend/twist stiffness matrix, with B1 the flexural rigidity about d1, B2 the flexural rigidity about d2 and B3 the twist rigidity about d3. Here, the vector κLo 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 κ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 s^ at all times, and the tangent to the centre line is also normal to the cross section, so that t=d3 [5]. This implies that nL 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, κ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 nL. In fact, they are no longer defined as Lagrange multipliers that enforce the condition of inextensibility and unshearability, i.e. we now must allow td3. 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

σL=Q(rs^d3)=Q(etd3).2.9
Here, the scalar field e(s^,t)=ds/ds^ expresses the local stretching or compression ratio (figure 1) relative to the rest reference configuration (s^) and 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(s^,t). Then, the following relations hold:

ds=eds^,A=A^e,I=I^e2,B=B^e2,S=S^eandκL=κ^Le.2.10

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

nL=S(σLσLo),2.11
where SR3×3=diag(S1, S2, S3) is the shear/stretch stiffness matrix, with S1 the shearing rigidity along d1, S2 the shearing rigidity along d2 and S3 the axial rigidity along d3. Here, as with the Kirchhoff rod, the vector σLo 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 κLo, σLo 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 σLo=κLo=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 αc=43 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 τL(κL,tκL) and nL(σL,tσ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

rt=v, 2.12
djt=(QTωL)×dj,j=1,2,3 2.13
dm2rt2=s^(QTS^σLe)ds^shear/stretch internal force+Fexternal force 2.14
anddJ^eωLt=s^(B^κ^Le3)ds^+κ^L×B^κ^Le3ds^bend/twist internal couple+(Qt×S^σL)ds^shear/stretch internal couple+(dJ^ωLe)×ωLLagrangian transport+dJ^ωLe2etunsteady dilatation+CLexternal couple, 2.15
where dm=ρA^ds^=ρAds is the infinitesimal mass element, and dJ^=ρI^ds^ is the infinitesimal mass second moment of inertia. We note that the left-hand side and the unsteady dilatation term of equation (2.15) arise from the expansion of the original rescaled angular momentum dJ^t(ωL/e) via the chain rule. We also note that the external force and couple are defined as F=efds^ and CL=ecLds^ (with f and cL 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 [6365]. 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 R=eθu:R3R3×3, and efficiently computed via the Rodrigues formula [66]

eθu=11+sinθU+(1cosθ)UU.3.1
Here, UR3×3 represents the skew-symmetric matrix associated with the unit vector u
U=[u]×=(0u3u2u30u1u2u10),u=[U]×1=(U3,2U3,1U2,1),
where the operator []×:R3R3×3 allows us to transform a vector into the corresponding skew-symmetric matrix, and vice versa []×1:R3×3R3.

Conversely, given a rotation matrix R, the corresponding rotation vector can be directly computed via the matrix logarithm operator log():R3×3R3

θu=log(R)={0ifθ=0θ2sinθ[RRT]×1ifθ0,θ(π,π),θ=arccos(trR12).

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 ri(t)R3,i[1,n+1] and a discrete set of material frames Qi(t)R3×3,i[1,n], as illustrated in figure 2.

Figure 2.

Figure 2. Discretization model. A discrete filament is represented through a set of vertices r(t)i=1,…,N+1 and a set of material frames Qi(t)={d1i,d2i,d3i,}i=1,,N. Two consecutive vertices define an edge of length ℓi along the tangent unit vector ti. The dilatation is defined as ei=i/^i, where ^i is the edge rest length. The vector σi=eitid3i represents the discrete shear and axial strains. The mass mr of the filament is discretized in pointwise concentrated masses mi=1,,n+1 at the locations ri for the purpose of advecting the vertices in time. For the evolution of Qi in time, we consider instead the mass second moment of inertia J^i=1,,n associated with the cylindrical elements depicted in blue.

Each vertex is associated with the following discrete quantities:

ri=1,,n+1vi=rit,mi,Fi,3.2
where vi is the velocity , mi is a pointwise concentrated mass and Fi 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

Qi=1,,ni=ri+1ri,i=|i|,^i=|^i|,ei=i^i,ti=iiσLi=Qi(eitidi3),ωLi,A^i,J^i,B^i,S^i,CLi,}3.3
where ℓi=|i|, ^i=|^i|, ei=i/^i are the edge current length, reference length and dilatation factor, ti is the discrete tangent vector, σLi is the discrete shear/axial strain vector, ωLi is the discrete angular velocity, A^i, J^i, B^i, S^i are the edge reference cross-section area, mass second moment of inertia, bend/twist matrix and shear/stretch matrix, respectively, and finally CLi 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 κL, are naturally expressed in an integrated form over the domain D along the filament [14,68]. Any integrated quantity divided by the corresponding integration domain length D=|D| is equivalent to its pointwise average. Therefore, following the approach of Bergou & Grinspun [14,68], the domain D becomes the Voronoi region Di of length

Di=i+1+i2,3.4
which is defined only for the interior vertices r(int)i=1,…,n−1. Each interior vertex is then also associated with the following discrete quantities:
ri=1,,n1(int)Di,D^i,Ei=DiD^i,κ^Li=log(Qi+1QiT)D^iandB^i=B^i+1^i+1+B^i^i2D^i,3.5
where D^i is the Voronoi domain length at rest and Ei is Voronoi region dilatation factor. Recalling that the generalized curvature expresses a rotation per unit length about its axis, the quantity D^iκ^Li naturally expresses the rotation that transforms a material frame Qi to its neighbour Qi+1 over the segment size D^i along the rod. Therefore, the relation eD^iκ^LiQi=Qi+1 holds, so that κ^Li=log(Qi+1QiT)/D^i. Finally, we introduce the bend/twist stiffness matrix B^i consistent with the Voronoi representation.

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

rit=vi,i=[1,n+1], 3.6
di,jt=(QiTωLi)×di,j,i=[1,n],j=1,2,3, 3.7
mivit=Δh(QiTS^iσLiei)+Fi,i=[1,n+1] 3.8
andJ^ieiωLit=Δh(B^iκ^LiEi3)+Ah(κ^Li×B^iκ^LiEi3D^i)+(Qiti×S^iσLi)^i+(J^iωLiei)×ωLi+J^iωLiei2eit+CLi,i=[1,n], 3.9
where Δh:{R3}N{R3}N+1 is the standard discrete difference operator and Ah:{R3}N{R3}N+1 is the averaging operator (trapezoidal quadrature rule) to transform integrated quantities over the domain D to their point-wise counterparts [14,69]. We note that Δh and Ah 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 ET=120LρAvTvds; rotational ER=120LρωLTIωLds; bending/twist EB=120LκLTBκ;Lds; shear/stretch ES=120LσLTSσLds. Therefore, by construction, in the limit of e1 (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 dt=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 dt=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 3a. Under these conditions, the filament buckles into a localized helical shape (figure 3e).

Figure 3.

Figure 3. Time–space convergence study for localized helical buckling. (a) We consider a rod originally straight whose ends are pulled together in the axial direction k with a slack D/2 and simultaneously twisted by the angle Φ/2. (b) Comparison between the analytical envelope function φ(s) and numerical approximations φn(s) at different levels of time–space resolution. Here, the time discretization δt is slaved by the spatial discretization δl=L/n according to δt=10−3δl s. (c) Norms L(ϵ) (black), L1(ϵ) (blue) and L2(ϵ) (red) are plotted against the number of discretization elements n. (d) Time evolution of the total energy of a rod (n=800, here the energy is computed assuming quadratic functionals, a suitable representation for an inextensible rod) simulated assuming no dissipation γ=0 (red line) versus the theoretical total energy EF=(MhΦ+ ThD)/2 (black dashed line). (e) Equilibrium rod configuration reqn numerically obtained given the discretization n=800, and assuming dissipation. For all studies, unless specified otherwise, we used the following settings: length L= 100 m, twist Φ=27⋅2π, slack D=3 m, linear mass density ρ= 1 kg m−1, bending stiffness α=1.345 Nm2, twisting stiffness β=0.789 Nm2, shear/stretch matrix S=105⋅1 N, bend/twist matrix B=diag(α,α,β) Nm2, dissipation constant γ= 10−2 kg (ms)−1, radius r=0.35 m, twisting time Ttwist=500 s, relaxation time Trelax=104 s.

The nonlinear equilibrium configuration req of the rod can be analytically determined [8,7375] 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 Mh and Th, respectively. Their normalized counterparts mh=MhL/(2πα) and th=ThL2/(4π2α) can be computed via the ‘semi-finite’ correction approach [74] by solving the system

DL=4π2th(1mh24th),Φ=2πmhβ/α+4arccos(mh2th).
Then, the analytical form of req can be expressed [75] as
req=L[12πth4thmh2 sech(πs¯4thmh2)sin(mhπs¯)]iL[12πth4thmh2 sech(πs¯4thmh2)cos(mhπs¯)]j+L[s¯12πth4thmh2 tanh(πs¯4thmh2)]k,3.10
where s¯=s/L0.5 is the normalized arc-length 0.5s¯0.5. Here, we make use of equation (3.10) to investigate the convergence properties of our solver in the limit of refinement. To compare analytical and numerical solutions, a metric invariant to rotations about k is necessary. Following Bergou et al. [14], we rely on the definition of the envelope φ
φ=cosθcosθmax1cosθmax,θ=arccos(tk)3.11
where θ is the angular deviation of the tangent t from the axial direction k, and θ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 L1(ϵ), L2(ϵ) and L(ϵ) 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 Ttwist. Subsequently, the ends of the rod are held in their final configurations for the period Trelax 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 3b,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 e1, we turn off the internal dissipation (figure 3d) and observe that the total energy of the filament EF is constant after Ttwist and matches its theoretical value EF=(MhΦ+ThD)/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 mp suspended at the other end, and due to its own mass mr, as depicted in figure 4a,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

ΔL=gmeqk=g(mp+mr/ξ)k,3.12
where k is the spring constant, ξ=2 is a constant factor and meq=mp+mr/ξ is the equivalent mass. Thus, the final equilibrium length of the rod reads L=L^+ΔL, with L^ being the rest unstretched length.
Figure 4.

Figure 4. Vertical oscillation under gravity. (a,d) We consider a vertical rod of mass mr clamped at the top and with a mass mp attached to the free end. Assuming that the rod is stiff enough (i.e. kA^E=const), it oscillates due to gravity around the equilibrium position L^+ΔL, where ΔL*=g(mp+mr/2)/k with a period T=2π(mp+mr/ξ)/k with ξ≃3 for mpmr, and ξπ2/4 for mpmr. Therefore, the rod oscillates according to L(t)=L^+[1+sin(2πt/Tπ/2)]ΔL. (ab) Case mpmr with mp=100 kg and mr=1 kg. (b) By increasing the stiffness E=107, 2×107, 3×107, 5×107, 108, 1010 Pa, the simulated oscillations (red lines) approach the analytical solution (dashed black line). (c) Convergence to the analytical solution in the norms L(ϵ) (black), L1(ϵ) (blue) and L2(ϵ) (red) with ϵ=∥L(t)−LE(t)∥, where LE is the length numerically obtained as a function of E. (cd) Case mpmr with mp=0 kg and mr=1 kg. (e) By increasing the stiffness E=104, 2×104, 3×104, 5×104, 105, 2×105, 109 Pa, the simulated oscillations approach the analytical solution. (f) Convergence to the analytical solution in the norms L(ϵ), L1(ϵ) and L2(ϵ) as a function of E. For all studies, we used the following settings: gravity g=9.81 m s−2, rod density ρ=103 kg m−3, shear modulus G=2E/3 Pa, shear/stretch matrix S^=diag(4GA^/3,4GA^/3,EA^)N, bend/twist matrix B^=diag(EI^1,EI^2,GI^3)Nm2, rest length L^=1m, rest cross-sectional area A^=mr/(L^ρ)m2, number of discretization elements n=100, timestep δt=T*/106, dissipation constant γ=0.

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

T=2π(mp+mr/ξ)k,L=L^+[1+sin(2πtTπ2)]ΔL.3.13
In this case, unlike the static solution, the factor ξ depends on the ratio mr/mp. In fact it can be shown [76] that ξ≃3 for mr/mp0, and ξπ2/4 for mr/mp.

The analytical results rely on the assumption of k being constant in space and time, given a fixed ratio mr/mp. 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, that is as a soft filament becomes stiff, the constant kEA^ and our rod model must recover the behaviour of the mass-spring oscillator. Indeed, figure 4b,c,e,f shows how the proposed numerical method converges to the analytical oscillation period T* and normalized longitudinal displacement (LL^)/ΔL 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 5a. 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].

Figure 5.

Figure 5. Time–space convergence study for a cantilever beam. (a) We consider the static solution of a beam clamped at one end s^=0 and subject to the downward force F at the free end s^=L^. (b) Comparison between the Timoshenko analytical y (black lines) and numerical yn (with n=400, red dashed lines) vertical displacements with respect to the initial rod configuration. As a reference we report in blue the corresponding Euler–Bernoulli solution. (c) Norms L(ϵ) (black), L1(ϵ) (blue) and L2(ϵ) (red) of the error ϵ=∥yyn∥ at different levels of time–space resolution are plotted against the number of discretization elements n. Here, the time discretization δt is slaved by the spatial discretization n according to δt=10−2δl seconds. For all studies, we used the following settings: rod density ρ=5000 kg m−3, Young’s modulus E=106 Pa, shear modulus G=104 Pa, shear/stretch matrix S^=diag(4GA^/3,4GA^/3,EA^)N, bend/twist matrix B^=diag(EI^1,EI^2,GI^3)Nm2, downward force F=15 N, rest length L^=3m, rest radius r^=0.25m, dissipation constant γ=10−1 kg (ms)−1, simulation time Tsim=5000 s.

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

y=FαcA^Gs^FL^2EI^1s^2+F6EI^1s^3,3.14
where L^ is the length of the rod, A^ is the constant cross-sectional area, I^1 is the area second moment of inertia about the axis j=k×i, E and G are, respectively, the Young’s and shear moduli, and αc=43 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 s^ (see figure 5a and the appendix for further details and derivation), hence the use of s^ in the above equation.

If the shear modulus G approaches infinity or if the ratio EI^1/(αcL^2A^G)1, then the Timoshenko model in the static case reduces to the Euler–Bernoulli approximation, yielding

y=FL^2EI^1s^2+F6EI^1s^3,3.15
as the shear term of equation (3.15) becomes negligible.

We compare our numerical model with these results by carrying out simulations of the cantilever beam of figure 5a in the time–space limit of refinement. As can be noticed in figure 5b 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(ϵ), L1(ϵ), L2(ϵ) of the error ϵ=∥yyn∥, where yn is the vertical displacement numerically obtained in the refinement limit.

We note that the norms L(ϵ) and L1(ϵ) exhibit first-order convergence, while L2(ϵ) 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 (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 EA^ (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 CL 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 nL and torque τ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 τL(κL) and forces nL(σL) of equation (2.11) become functions of both strain and rate of strain, i.e. τL(κL,tκL) and nL(σL,tσ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 fv and torques cLv 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

fv=γtv4.1
and
cLv=γrωL.4.2
This approach does not model the physical nature of viscoelastic phenomena, although it does dissipate energy, effectively mimicking overall material friction effects. In the context of our numerical investigations, we did not observe any appreciable difference between the two outlined methods, so that, for the sake of simplicity and computational efficiency, we opted for the second one. Throughout the remainder of the text we will then employ equations (4.1) and (4.2) with a single dissipation constant γ, 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 Am as a travelling wave propagating head to tail along the filament

Am=βm(s^)sin(2πTmt+2πλms^+ϕm),4.3
where ϕm is the phase, t is time and Tm and λm are, respectively, the activation period and wavelength. The amplitude of the travelling wave is represented by the cubic B-spline β(s^) characterized by Nm control points (S^i,βi) with i=0,,Nm1, as illustrated in figure 6. The first and last control points are fixed so that (S^0,β0)=(0,0) and (S^Nm1,βNm1)=(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].
Figure 6.

Figure 6. Muscular activity. Example of muscular activity amplitude profile (solid black line) described by cubic B-spline through Nm=8 control points (S^i,βi) with i=0,…,Nm−1. The control points are located along the filament at the positions S^i, and are associated with the amplitude values βi. The first and last control points are fixed so that (S^0,β0)=(0,0) and (S^Nm1,βNm1)=(L^,0), therefore assuming the ends of the deforming body to be free.

A given activation mode can be achieved by multiplying the scalar amplitude Am 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

nLm=Q(Amd3).4.4
Similarly, if we wish to investigate a slithering snake characterized by a planar kinematic wave, we may consider a torque activation of the form
τLm=Q(Amd1),4.5
assuming d2 and d3 to be coplanar to the ground. These two contributions are directly added to the internal force nL and torque τL resultants.

In the most general case, all deformation modes can be excited by enabling force and torque muscular activity along all directors d1, d2 and d3, 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 Fsc acting between the discrete elements in contact. To determine whether any two cylindrical elements are in contact, we calculate the minimum distance dminij between edges i,j by parametrizing their centre lines ci(h)=si+h(si+1si) so that

dminij=maxh1,h2[0,1]ci(h1)cj(h2).4.6
If dminij is smaller than the sum of the radii of the two cylinders, then they are considered to be in contact and penalty forces are applied to each element as a function of the scalar overlap ϵij=(ri+rjdminij), where ri and rj 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 dminij 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
Fsc=H(ϵij)[kscϵijγsc(vivj)dminij]dminij,4.7
where 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 ksc, while the second damping term models contact dissipation and is proportional to the coefficient γsc and the interpenetration velocity vivj.

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

Fw=H(ϵ)(F+kwϵγwvuw)uw,4.8
where H(ϵ) denotes the Heaviside function and ensures that a wall force is produced only in case of contact (ϵ≥0). Here, uw 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 kw and γw are, respectively, the wall stiffness and dissipation coefficients.
Figure 7.

Figure 7. Contact model with solid boundaries. Obstacles and surfaces (grey) are modelled as soft boundaries allowing for the interpenetration ϵ=rd with the elements of the filament (blue) characterized by radius r and distance d from the substrate. The surface normal uw determines the direction of the wall’s response Fw to contact. We note that Fw balances the sum of all forces F that push the rod against the wall, and is complemented by the other two components, which allow it to amend to possible interpenetration due to numerics. These components are an elastic one (kwϵ) and a dissipative one (γwvuw), where kw 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 uw, uw, uw×, as illustrated in figure 8.

Figure 8.

Figure 8. Surface friction. (a) The forces produced by friction effects between an element of the rod and the substrate are naturally decomposed into a lateral component in the direction uw=t×uw and a longitudinal one in the direction u×w=uw×uw. We note that in general uw×t. The notation x, x, x× denotes the projection of the vector x in the directions uw, uw, uw×. (b,c) Kinematic and dynamic quantities at play at any cross section in the case of (b) rolling and slipping and (c) pure rolling motion. Red arrows correspond to forces and torques, green arrows correspond to velocities, and black arrows correspond to geometric quantities.

The longitudinal friction force Flong 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 uw× (figure 8). The Amonton–Coulomb model then reads

Flong={max(|F×|,μs|F|)F×|F×|if |v×|v,μk|F|v×|v×|if |v×|>v,
where vϵ0 is the absolute velocity threshold value employed to distinguish between static (|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 uw=u×w×uw is associated with both translational (v) and rotational (ω×=ω×u×w) motions, as illustrated in figure 8b,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

vslip=v+vcont,vcont=ruw×ω×,4.9
where vcont 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 (vslip=0), the linear momentum balance in the direction uw, and the angular momentum balance about the axis uw× express a kinematic constraint between the linear acceleration auw and angular acceleration ω×=(uw×auw)/r, so that

(F+Froll)uw=dmauw4.10
and
T×u×wruw×Frolluw=Juw×auwr,4.11
where F=Fuw and T×=T×u×w are the forces and torques acting on the local element, and Froll=Frolluw is the rolling friction force at the substrate–filament interface necessary to meet the no-slip condition. By recalling that a disk mass second moment of inertia about uw× is J=r2 dm/2, the above system can be solved for the unknown a and Froll, yielding
Froll=rF2T×3ruw.4.12

Therefore, the lateral friction force Flat and the associated torque CLlat can be finally expressed as

Flat={max(|Froll|,μsr|F|)Froll|F roll|if |vslip|vϵμkr|F|v slip|v slip|if |vslip|>vϵ,CLlat=Q(Flat×ruw),
where μrs and μrk 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=ρfUL/μ≪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 [3136]). At leading order resistive force line densities scale linearly with the local rod velocities v according to

fH=4πμln(L/r)(I12tTt)v.4.13
We note that the matrix (I12tTt) introduces an anisotropic effect for which
fH=2πμln(L/r)|(vt)t|,fH=4πμln(L/r)|v(vt)t|,4.14
where fH=(fHt)t and fH=fHfH are, respectively, tangential and orthogonal viscous drag components. The coupling of liquid environment, filament mechanics and muscular activity provides a flexible platform to characterize biological locomotion at the microscopic scale (bacteria, protozoa, algae, etc.) and to design propulsion strategies in the context of artificial micro-swimmers [45,86].

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 [4144].

Given the broad scope of our computational framework, we can now study the formation of both solenoids and plectonemes. As illustrated in figure 9a, a soft rod of Young’s modulus E=106 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/L^1), plectonemes are generated (figure 9b). When the load F is increased so that the elongation of the rod approaches L/L^1.15, solenoids arise as predicted in [40] and illustrated in figure 9c. 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.

Figure 9.

Figure 9. Formation of plectonemes and solenoids. (a) We consider a soft rod clamped at one end, subject to a constant vertical load F and twisted R times at the other end. (b) Formation of a plectoneme for F=0 (leading to the total elongation L/L^1) and R=4. (c) Formation of a solenoid for F=300 N (leading to the total elongation L/L^1.15) and R=13. Settings: length L=1 m, radius r=0.025 m, mass m=1 kg, Young’s modulus E= 106 Pa, shear modulus G=2E/3 Pa, shear/stretch matrix S^=diag(4GA^/3,4GA^/3,EA^)N, bend/twist matrix B^=diag(EI1,EI2,GI3)Nm2, dissipation constant γ=2 kg (ms)−1, ksc= 104 kg s−2, γsc=10 kg s−1, discretization elements n=100, timestep δt=0.01δl s, Ttwist=75 s, Trelax= 50 s.

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 μkf:μkb:μkr=1:1.5:2 and μsf:μsb:μsr=1:1.5:2 with μsf=2μkf, as experimentally observed for juvenile Pueblan milk snakes on a moderately rough surface [83]. The value of the friction coefficient μfk is set so that the ratio between inertial and friction forces captured by the Froude number is Fr=(L/Tm2)/(μkfg)=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 Nm=6 control points and a unit oscillation period Tm, 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 vfwdmax over one activation cycle Tm. 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=vfwdmax, 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 10e,f the forward scaled average speed approaches vfwdmax≃0.6, consistent with experimental evidence [97]. Moreover, CMA-ES finds that the optimal wavelength is λmL (figure 10g), 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.

Figure 10.

Figure 10. Optimal lateral undulation gait. (ad) Instances at different times of a snake characterized by the identified optimal gait. (e) Evolution of the fitness function f=vfwdmax as a function of the number of generations produced by CMA-ES. Solid blue, solid red and dashed black lines represent, respectively, the evolution of f corresponding to the best solution, the best solution within the current generation, and the mean generation value. (f) Scaled forward (red) and lateral (blue) centre of mass velocities versus normalized time. (g) Gait envelope over one oscillation period Tm. Red lines represent head and tail displacement in time. Settings: Froude number Fr=0.1, length L=1 m, radius r= 0.025 m, density ρ=103 kg m−3, Tm=1 s, Young’s modulus E=107 Pa, shear modulus G=2E/3 Pa, shear/stretch matrix S^=diag(4GA^/3,4GA^/3,EA^) N, bend/twist matrix B^=diag(EI1,EI2,GI3)Nm2, dissipation constant γ= 5 kg (ms)−1, gravity g=9.81 m s−2, friction coefficient ratios μkf:μkb:μkr=1:1.5:2 and μsf:μsb:μsr=1:1.5:2 with μsf=2μkf, friction threshold velocity vϵ=108ms1, ground stiffness and viscous dissipation kw=1 kg s−2 and γw=10−6 kg s−1, discretization elements n=50, timestep δt=2.5⋅10−5 Tm, wavelength λm=0.97L, phase shift ϕm=0, torque B-spline coefficients βi=0,…,5={0,17.4,48.5,5.4,14.7,0} Nm, bounds maximum attainable torque |β|maxi=0,…,5=50 Nm.

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=vfwdmax within one activation cycle Tm, 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 Ls≃40 μm. The gait corresponding to planar swimming is characterized by kinematic undulations of wavelength λs<Ls and frequency fs≃2.8. At Re≃10−4 the spermatozoon attains the scaled velocity vs=Us/(fsLs)≃0.08±0.03, where Us 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 11e,f, the forward average scaled speed and wavelength approach vfwdmax≃0.055 and λm≃0.38L are qualitatively and quantitatively consistent with experimental evidence [99].

Figure 11.

Figure 11. Optimal planar swimming gait at low Reynolds number. (ad) Instances at different times of a filament swimming according to the identified optimal gait. (e) Evolution of the fitness function f=vfwdmax as a function of the number of generations produced by CMA-ES. Solid blue, solid red and dashed black lines represent, respectively, the evolution of f corresponding to the best solution, the best solution within the current generation, and the mean generation value. (f) Scaled forward (red) and lateral (blue) centre of mass velocities versus normalized time. (g) Gait envelope over one oscillation period Tm. Red lines represent head and tail displacement in time. Settings: Reynolds number Re=10−4, length L=1 m, radius r=0.025 m, filament density ρ=103 kg m−3, Tm=1 s, Young’s modulus E=107 Pa, shear modulus G=2E/3 Pa, shear/stretch matrix S^=diag(4GA^/3,4GA^/3,EA^) N, bend/twist matrix B^=diag(EI1,EI2,GI3)Nm2, dissipation constant γ=5 kg (ms)−1, discretization elements n=50, timestep δt=2.5×10−5 Tm, wavelength λm=2.6L, phase shift ϕm=0, torque B-spline coefficients βi=0,…,5={0,50,50,50,50,0} Nm, bounds maximum attainable torque |β|maxi=0,…,5=50 Nm.

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

v=rtandt=rs,A 1
with |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

djt=(QTej)t=QTtej=QTtQdj=ω×dj,j=1,2,3A 2
and
djs=(QTej)s=QTtej=QTsQdj=κ×dj,j=1,2,3,A 3
where the equivalences ∂tQTQ=ω×(⋅) and ∂sQTQ=κ×(⋅) 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
rt=v, A 4
djt=ω×dj,j=1,2,3, A 5
(ρAv)t=ns+f A 6
andht=τs+rs×n+c, A 7
where ρ is the constant material density, A is the cross-sectional area in its current state (so that ρAv 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 Qh=hL=ρIωL, where the tensor I is the second area moment of inertia which, assuming circular cross sections, takes the form

I=(I1000I2000I3)=A24π(100010002)=A24πdiag(1,1,2).

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

dj=QTej,j=1,2,3A 8
and definitions of time and space derivatives of its frame directors
djt=ω×dj,j=1,2,3A 9
and
djs=κ×dj,j=1,2,3.A 10

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

xt=(QTxL)t A 11
=(j=13djxLj)t A 12
=j=13t(djxLj) A 13
=j=13(tdj)xLj+j=13dj(txLj) A 14
=j=13(ω×dj)xLj+j=13dj(txLj) A 15
xt=ω×(QTxL)+QTxLt, A 16
andxs=κ×(QTxL)+QTxLs. A 17
Because the second moment of inertia is defined in the material frame, the governing equation for the angular momentum is also most naturally expressed in the material frame (I in the laboratory frame is a function of space and time rendering its use cumbersome, while IL 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.
(ρIω)t=QT(ρIω)Lt+ω×(QT(ρIω)L)(LHS) A 18
τs=QTτLt+ω×(QTτL)(RHS) A 19
rs×n=QT(rs)L×QTnL(RHS) A 20
andc=QTcL(RHS). A 21
Multiplying by Q and noting that Iω=QTILωL yields our final governing equation for the angular momentum, expressed in the material frame
(ρILωL)t=τLs+κL×τL+Qrs×nL+(ρILωL)×ωL+cL.A 22
Aside from the balance of angular momentum, converting equations (A 5) and (A 6) to the material frame simply requires a direct application of the definition of the material frame to the RHS quantities.

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 L^=100mm and cross section A^=10mm2. As the beam is stretched the cross-sectional area decreases. In one-dimensional, the Neo-Hookean model reads:

lime1E=6C10,α=2C10(e1e2)eandF=αA=αA^e=2C10A^(e1e2),
where α is the stress, e=L/L^ is the elongation factor (related to the strain σ=e−1), C10 is a material property to be determined experimentally and that in the limit of small strains (e1) is related to the Young’s modulus E, as above indicated. In one-dimensional, the Mooney–Rivlin model reads:
lime1E=6(C10+C01),α=2C10(e21e)+2C01(e21e)1e,andF=αA=αA^e=A^e[2C10(e21e)+2C01(e21e)1e],
where C10 and C01 are material properties to be determined experimentally and that in the limit of small strains (e1) are related to the Young’s modulus E as indicated above. In our approach the load reads
F=AEσ=AE(e1)=A^eE(e1).
Assuming that the beam is made of rubber with E=107 Pa, and that in the Mooney–Rivlin model C10=C01, 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.
Figure 12.

Figure 12. Comparison with hyperelastic models. Load response in a stretched rubber beam for different models.

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

ri(t+δt2)=ri(t)+δt2vi(t),i=[1,n+1] C 1
Qi(t+δt2)=exp[δt2ωLi(t)]Qi(t),i=[1,n] C 2
vi(t+δt)=vi(t)+δtdvidt(t+δt2),i=[1,n+1] C 3
ωLi(t+δt)=ωLi(t)+δtdωLidt(t+δt2),i=[1,n] C 4
ri(t+δt)=ri(t+δt2)+δt2vi(t+δt2)i=[1,n+1] C 5
andQi(t+δt)=exp[δt2ωLi(t+δt2)]Qi(t+δt2).i=[1,n]. C 6
The time integrator is structured in three blocks: a first half step position update (equations (C 1) and (C 2)), followed by the evaluation of local accelerations (equations (C 4) and (C 5)), and finally, the second half step position update (equations (C 5) and (C 6)). Therefore, the second-order position Verlet scheme entails only one right-hand side evaluation of equations (3.8) and (3.9), the most computationally expensive operation. We also emphasize that for the numerical integration of equations (C 2) and (C 6) the Rodrigues formula is employed, implying the direct use of ωL instead of ω=QTωL. By foregoing an implicit integration scheme we can quickly incorporate a number of additional physical effects and soft constraints. This algorithm strikes a balance between computing costs, numerical accuracy and implementation modularity.

Appendix D. Discrete operators

The discrete operators Δh:{R3}N{R3}N+1 and Ah:{R3}N{R3}N+1 take the form

yj=1,,N+1=Δh(xi=1,,N)={x1if j=1,xjxj1if 1<jNxNif j=N+1,,yj=1,,N+1=Ah(xi=1,,N)={x12if j=1,xj+xj12if 1<jNxN2if j=N+1.,
We note that Δh and Ah 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 5a. 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.

Figure 13.

Figure 13. Free body diagram of an infinitesimal beam element. The sketch represents an infinitesimal element undergoing bending and shear deformations, assuming small deflections so that xs. The bending moment is denoted as M, the shear force as V , the vertical displacement as y. The bending angle ψ, the shear angle θ and the slope of the vertical displacement dxy are related to each other so that ψ=θ+dxy.

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

ψs=MEI1,E 1
where M is the bending moment, E the Young’s modulus, and I1 is the second area moment of inertia about the axis j=k×i (figure 5a). The shear angle θ, as illustrated in figure 13, is the difference between the bending angle and the slope of the centreline, so that
ψys=VαcAG,E 2
where 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(Ls) and V =F, so that
ψs=F(Ls)EI1E 3
and
ψys=FαcAG.E 4

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

y=FαcA^Gs^FL^2EI^1s^2+F6EI^1s^3.E 5

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 14a. The critical axial load Fc that the rod can withstand before bending can be expressed analytically [3] as

Fc=π2EI(bL)2,F 1
where E is the modulus of elasticity of the rod, I=I1=I2 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
Fc=π2αL2,F 2
where α=EI is the bending stiffness.
Figure 14.

Figure 14. Euler buckling instability. (a) Inextensible and unshearable rod vertically initialized and subject to the axial load F. (b) Probability P of observing an instability event as a function of the compression force F and the bending rigidity α for a fixed length L=1 m. The corresponding analytical solution is overlaid as reference (white line). (c) Probability P of observing an instability event as a function of the compression force F and the length L for a fixed bending rigidity α=1 Nm2. The corresponding analytical solution is overlaid as reference (white line). The probability P is determined by performing 10 simulations for each pair of parameters (F,α) or (F,L). Each simulation is initially perturbed by applying to every discretization node a small random force sampled from a uniform distribution, such that FRU(0,102)N. The occurrence of an instability is detected whenever the rod bending energy EB>Eth with Eth=10−3 J. Settings: rod’s mass mr=1 kg, twist stiffness β=2α/3 Nm2, shear/stretch matrix S=105×1 N, bend/twist matrix B=diag(α,α,β) Nm2, radius r=0.025L m, discretization elements n=100 m−1, timestep δt= 10−5 s, simulation time Tsim=10 s, dissipation constant γ=0.

We test our solver against the above analytical solution by simulating the time evolution of an inextensible and unshearable rod. In figure 14a, 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 EB>Eth. As can be seen in figure 14b,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

Φc=23πβ/α,F 3
where α 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 ET>Eth. As can be seen in figure 15, our results show a sharp transition in agreement with the exact solution of equation (F 3).

Figure 15.

Figure 15. Mitchell buckling instability. (a) Probability P of observing an instability event as a function of the total twist Φ and the ratio between bending and twist rigidities β/α. The corresponding analytical solution is overlaid as reference (white line). The probability P is determined by performing 10 simulations for each pair of parameters (Φ,β/α). Each simulation is initially perturbed by applying to every discretization node a small random force sampled from a uniform distribution, such that FRU(0,103)N. The occurrence of an instability is detected whenever the rod translational energy ET>Eth with Eth=10−4 J. (bc) Visualization of a Mitchell bucking instability event for Φ=10 and β/α=2. Settings: rod’s mass mr=1 kg, length L= 1 m, bending stiffness α=1 Nm2, shear/stretch matrix S= 105×1 N, bend/twist matrix B=diag(α,α,β) Nm2, radius r= 0.01L m, discretization elements n=50 m−1, timestep δt= 10−5 s, simulation time Tsim=2 s, dissipation constant γ=0.

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 16a. 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).

Figure 16.

Figure 16. Time-space convergence study for twist forced vibrations in a stretched rod. (a) We consider a rod clamped at one end, forced to vibrate by applying the periodic couple Avsin(2πfvt) to the free end, and characterized by rest length L^ which is extended to a final length L=eL^. (b) Comparison between analytical θ (black lines) and numerical θn (red dashed lines) angular displacement with respect to the reference configuration along a stretched rod. Each red (numerical simulation) and black (analytical solution) line corresponds to the angular displacement along a rod discretized with n=1600 elements, and sampled at regular intervals Δt=Tv/10 within one loading period Tv=fv1. (c) Norms L(ϵ) (black), L1(ϵ) (blue) and L2(ϵ) (red) of the error ϵ=∥θθn∥ at different levels of time-space resolution are plotted against the number of discretization elements n. Here, the time discretization δt is slaved by the spatial discretization n according to δt=10−4δl seconds. For all studies, we used the following settings: rod’s density ρ=10 kg m−3, Young’s modulus E=106 Pa, shear modulus G=2E/3 Pa, shear/stretch matrix S^=diag(4GA^/3,4GA^/3,EA^)N, bend/twist matrix B^=diag(EI^1,EI^2,GI^3)Nm2, forcing amplitude Av=103 Nm, forcing frequency fv=1 s−1, dilatation factor e=1.05, rest length L^=E/ρ/(efv)m, rest radius r^=0.5m, simulation time Tsim=2000 s. We enabled dissipation in the early stages of the simulations, letting γ decay exponentially in time to a zero value.

The rod is clamped at one end, stretched by the quantity ΔL=(e1)L^, and forced to angularly vibrate about the axial direction by applying at the free end the couple Avsin(2πfvt), where Av and fv 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]

2θs2=1cs2θt2,F 4
where cs=G/ρ is the shear wave velocity, G is the shear modulus and ρ is the density. By assuming a solution of the form θ(s,t)=ϕ(s)sin(2πfvt), and by applying the boundary conditions ϕ(0)=0 and dsϕ(0)=Cv/(I^3G) with Cv twisting torque and I^3 second area moment of inertia about the axial direction, we can solve equation (F 4) obtaining
θ(s,t)=Avcs2πfvG(I^3/e2)sin(2πfves^/cs)cos(2πfveL^/cs)sin(2πfvt).F 5

As can be seen in figure 16b, 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 16c).

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 17a. 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 ET=mv2/2 and rotational ER=2/2 energies, can be analytically determined.

Figure 17.

Figure 17. Rolling static and dynamic friction. (a) A rod initially at rest on a plane inclined of the angle αs rolls or slips down with linear and angular velocities v and ω due to its own weight mg. (b) Translational ET=mv2/2 (analytical solution: black line, numerical solution: blue line) and rotational ER=2/2 (analytical solution: dashed black line, numerical solution: orange line) energies are plotted against the angle αs. Settings: length L=1 m, radius r= 0.025 m, mass m=1 kg, Young’s modulus E=109 Pa, shear modulus G=2E/3 Pa, shear/stretch matrix S=104×1 N, bend/twist matrix B=diag(EI1,EI2,GI3) Nm2, dissipation constant γ=10−6 kg (ms)−1, gravity g=9.81 m s−2, static and kinetic friction coefficients μs=0.4 and μk=0.2, friction threshold velocity vϵ=104ms1, ground stiffness and viscous dissipation kw=10 kg s−2 and γw=10−4 kg s−1, discretization elements n=50, timestep δt= 10−6 s, simulation time T=0.5 s. (c) Rod set in motion by the external couple Cs on a horizontal plane. (d) Translational ET and rotational ER energies are plotted against the couple Cs. Colour scheme and settings are identical to those of panel (b) except for the simulation time T=1 s. (e) Rod with initial velocity V s slows down due to kinetic friction until the no slip condition is reached. (f) Translational ET and rotational ER energies are plotted against the relative mass moment of inertia ratio λ/2. Colour scheme and settings are identical to those of (b) except for the simulation time T=2 s, and the friction threshold velocity vϵ=10−6 m s−1.

By recalling equation (4.12), the force necessary to ensure rolling without slip takes the form Fnoslip=F/3=mgsin(αs)/3. Given the maximum static friction force Fs=μsF=μsmgcos(αs) in the case |Fnoslip|≤|Fs|, we have a=(F+Fnoslip)/m. Therefore, at the time T after releasing the rod, linear and angular velocities read, respectively, v=aT and ω=ω˙T=(a/r)T, expressing the no slip kinematic constraint between linear a and angular ω˙ accelerations. On the contrary, if |Fnoslip|>|Fs| the rod starts slipping and linear and angular accelerations are no longer coupled, so that after time T we have a=(FμkF)/m, v=aT, ω=ω˙T=(μkFr/J)T. Therefore, translational and rotational energies as a function of αs finally read

ET={2mg2T2sin2(αs)9if |Fnoslip||Fs|mg2T2[sin(αs)μkcos(αs)]22if |Fnoslip|>|Fs|,
and
ER={2Jg2T2sin2(αs)9r2if |Fnoslip||Fs|μk2m2g2r2T2cos2(αs)2Jif |Fnoslip|>|Fs|.
As can be noted in figure 17b, 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 17c,d entails a rod set in motion by the external couple Cs 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 Fnoslip=2Cs/(3r). Given the maximum static friction force Fs=μsF=μsmg, in the case of |Fnoslip|≤|Fs| we have a=Fnoslip/m. Therefore, at the time T after releasing the rod, linear and angular velocities read, respectively, v=aT and ω=ω˙T=(a/r)T, expressing the no slip kinematic constraint between linear a and angular ω˙ accelerations. On the contrary, if |Fnoslip|>|Fs| the rod starts slipping and linear and angular accelerations are no longer coupled, so that after time T we have a=μkF/m, v=aT, ω=ω˙T=J1(CsμkFr)T. Therefore, translational and rotational energies as a function of Cs finally read

ET={2T2Cs29mr2if |Fnoslip||Fs|,mμk2g2T22if |Fnoslip|>|Fs|
and
ER={2JT2Cs29r4m2if |Fnoslip||Fs|,(Csμkmgr)2T22Jif |Fnoslip|>|Fs|.
As can be noted in figure 17d, 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 17e), and we vary the relative mass moment of inertia ratio λ=2J/(mr2). 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 veq=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 sFT/m and ω=rFT/J. By enforcing the no slip kinematic constraint V sFT/m=r2FT/J, we have that FT=V s/(r2/J+1/m). This allows us to directly compute veq and ωeq and the corresponding translational and rotational energies, yielding

ETeq=mVs221(1+λ/2)2andEReq=mVs22λ/2(1+λ/2)2.G 1
As depicted in figure 17f, 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 μfs, μkf and backward μbs, μkb static and kinetic coefficients. For Ft≥0, the rod translational energy takes the form

ETf={0if |F||μsfmg|T22m(|F|μkfmg)2if |F|>|μsfmg|,
while for Ft<0 we have
ETb={0if |F||μsbmg|T22m(|F|μkbmg)2if |F|>|μsbmg|.
As can be seen in figure 18, our numerical method reproduces the above theoretical predictions.
Figure 18.

Figure 18. Anisotropic static and kinetic longitudinal friction. A rod initially at rest on a horizontal plane is pulled longitudinally with force F. The system is characterized by anisotropic friction so that resistance forces in the forward direction Fflong are larger than in the backward direction Fblong. The plot illustrates the behaviour of the total rod’s translational energy ET as a function of the force F applied. Blue and solid black lines correspond to, respectively, simulated and analytical ET for a rod pulled forward. Orange and dashed black lines correspond to, respectively, simulated and analytical ET for a rod pulled backward. Settings: length L=1 m, radius r= 0.025 m, mass m=1 kg, Young’s modulus E=105 Pa, shear modulus G=2E/3 Pa, shear/stretch matrix S=104×1 N, bend/twist matrix B=diag(EI1,EI2,GI3) Nm2, dissipation constant γ=10−6 kg (ms)−1, gravity g=9.81 m s−2, static and kinetic forward friction coefficients μfs=0.8 and μfk=0.4, static and kinetic backward friction coefficients μbs=0.4 and μbk=0.2, friction threshold velocity vϵ=104ms1, ground stiffness and viscous dissipation kw=10 kg s−2 and γw=10−4 kg s−1, discretization elements n=50, timestep δt=10−5 s, simulation time T=0.25 s.

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 19a,b).

Figure 19.

Figure 19. Isotropic versus anisotropic friction driven locomotion. (a) Gait envelope computed over the 10th muscular activation cycle in the case of isotropic friction. The blue triangle denotes the location of the snake’s centre of mass at time t=0, reported as reference. (b) Lateral (blue) and forward (red) velocities as functions of time normalized by the activation period Tm in the case of isotropic friction. (c) Gait envelope computed over the 10th muscular activation cycle in the case of anisotropic friction. The blue triangle denotes the location of the snake’s centre of mass at time t=0, reported as reference. (d) Lateral (blue) and forward (red) velocities as functions of time normalized by the activation period Tm in the case of anisotropic friction. Settings: length L=0.5 m, radius r= 0.025 m, mass m=1 kg, Young’s modulus E=107 Pa, shear modulus G=2E/3 Pa, shear/stretch matrix S=105×1 N, bend/twist matrix B=diag(EI1,EI2,GI3) Nm2, dissipation constant γ= 10−1 kg (ms)−1, gravity g=9.81 m s−2, static μfs=0.2, μsr=μsf, μsb=μsf and kinetic μfk=0.1, μkr=μkf, μkb=μkf friction coefficients in the isotropic case, static μfs=0.2, μsr=2μsf, μsb=20μsf and kinetic μfk= 0.1, μkr=2μkf, μkb=20μkf friction coefficients in the anisotropic case, friction threshold velocity vϵ=108ms1, ground stiffness and viscous dissipation kw=1 kg s−2 and γw=10−6 kg s−1, discretization elements n=100, timestep δt= 10−5 s, muscular activation period Tm=1 s, wavelength λm=, phase shift ϕm=0, torque B-spline coefficients βi=0,…,5={0,10,15,15,10,0} Nm.

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 dt=a dℓ, with a∼10−2 s m−1, and found it reliable in preventing numerical instabilities (figure 20).

Figure 20.

Figure 20. Timings for various representative simulations. The time-to-solution values reported correspond to an average over 10 measurements for the same study case.

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

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.

References

Comments