Forward and inverse problems in the mechanics of soft filaments

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.

Keywords: Cosserat theory; computational mechanics; soft filaments.

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 {d₁,d₂,d₃}. The corresponding orthogonal rotation matrix Q with row entries d₁, d₂, d₃ transforms a vector x from the laboratory canonical basis {i,j,k} to the material frame of reference {d₁,d₂,d₃} so that ${\mathbf{\text{x}}}_{\mathcal{L}}=\mathbf{\text{Q}}\mathbf{\text{x}}$ and vice versa $\mathbf{\text{x}}={\mathbf{\text{Q}}}^{\mathrm{T}}{\mathbf{\text{x}}}_{\mathcal{L}}$. If extension or compression is allowed, the current filament configuration arc-length s may no longer coincide with the rest reference arc-length $\hat{s}$. This is captured via the scalar dilatation field $e=\mathrm{d}s/\mathrm{d}\hat{s}$. Moreover, to account for shear we allow the triad {d₁,d₂,d₃} to detach from the unit tangent vector t so that d₃≠t (we recall that the condition d₃=t and e=1 correspond to the Kirchhoff constraint for unshearable and inextensible rods, and implies that σ=et−d₃=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.

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 ${\mathbf{\text{Q}}}_i(t)={\{{\mathbf{\text{d}}}_1^i,{\mathbf{\text{d}}}_2^i,{\mathbf{\text{d}}}_3^i,\}}_{i=1,\dots ,N}$. Two consecutive vertices define an edge of length ℓ_i along the tangent unit vector t_i. The dilatation is defined as $e_i={\ell }_i/{\hat{\ell }}_i$, where ${\hat{\ell }}_i$ is the edge rest length. The vector ${\mathit{\sigma }}_i=e_i{\mathbf{\text{t}}}_i-{\mathbf{\text{d}}}_3^i$ represents the discrete shear and axial strains. The mass m_r of the filament is discretized in pointwise concentrated masses $m_{i=1,\dots ,n+1}$ at the locations r_i for the purpose of advecting the vertices in time. For the evolution of Q_i in time, we consider instead the mass second moment of inertia ${\widehat{\mathbf{\text{J}}}}_{i=1,\dots ,n}$ associated with the cylindrical elements depicted in blue.

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 φⁿ(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⁻³δl s. (c) Norms $L^{\mathrm{\infty }}(ϵ)$ (black), L¹(ϵ) (blue) and L²(ϵ) (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 E_F=(M_hΦ+ T_hD)/2 (black dashed line). (e) Equilibrium rod configuration ${\mathbf{\text{r}}}_{\mathrm{eq}}^n$ 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⁻¹, bending stiffness α=1.345 Nm², twisting stiffness β=0.789 Nm², shear/stretch matrix S=10⁵⋅1 N, bend/twist matrix B=diag(α,α,β) Nm², dissipation constant γ= 10⁻² kg (ms)⁻¹, radius r=0.35 m, twisting time T_twist=500 s, relaxation time T_relax=10⁴ s.

Figure 4.

Vertical oscillation under gravity. (a,d) We consider a vertical rod of mass m_r clamped at the top and with a mass m_p attached to the free end. Assuming that the rod is stiff enough (i.e. $k\simeq \hat{A}E=\text{const}$), it oscillates due to gravity around the equilibrium position $\hat{L}+\mathrm{\Delta }{L}^{\ast }$, where ΔL*=g(m_p+m_r/2)/k with a period $T^{\ast }=2\pi \sqrt{(m_p+m_r/\xi )/k}$ with ξ≃3 for m_p≫m_r, and ξ≃π²/4 for m_p≪m_r. Therefore, the rod oscillates according to $L(t)=\hat{L}+[1+\sin (2\pi t/T^{\ast }-\pi /2)]\mathrm{\Delta }{L}^{\ast }$. (a–b) Case m_p≫m_r with m_p=100 kg and m_r=1 kg. (b) By increasing the stiffness E=10⁷, 2×10⁷, 3×10⁷, 5×10⁷, 10⁸, 10¹⁰ Pa, the simulated oscillations (red lines) approach the analytical solution (dashed black line). (c) Convergence to the analytical solution in the norms $L^{\mathrm{\infty }}(ϵ)$ (black), L¹(ϵ) (blue) and L²(ϵ) (red) with ϵ=∥L(t)−L^E(t)∥, where L^E is the length numerically obtained as a function of E. (c–d) Case m_p≪m_r with m_p=0 kg and m_r=1 kg. (e) By increasing the stiffness E=10⁴, 2×10⁴, 3×10⁴, 5×10⁴, 10⁵, 2×10⁵, 10⁹ Pa, the simulated oscillations approach the analytical solution. (f) Convergence to the analytical solution in the norms $L^{\mathrm{\infty }}(ϵ)$, L¹(ϵ) and L²(ϵ) as a function of E. For all studies, we used the following settings: gravity g=9.81 m s⁻², rod density ρ=10³ kg m⁻³, shear modulus G=2E/3 Pa, shear/stretch matrix $\widehat{\mathbf{\text{S}}}=\text{diag}(4G\hat{A}/3,4G\hat{A}/3,E\hat{A})\hspace{0.17em}\mathrm{N}$, bend/twist matrix $\widehat{\mathbf{\text{B}}}=\text{diag}(E{\hat{I}}_1,E{\hat{I}}_2,G{\hat{I}}_3)\hspace{0.17em}{\mathrm{Nm}}^2$, rest length $\hat{L}=1\hspace{0.17em}\mathrm{m}$, rest cross-sectional area $\hat{A}=m_r/(\hat{L}\rho )\hspace{0.17em}{\mathrm{m}}^2$, number of discretization elements n=100, timestep δt=T*/10⁶, dissipation constant γ=0.

Figure 5.

Time–space convergence study for a cantilever beam. (a) We consider the static solution of a beam clamped at one end $\hat{s}=0$ and subject to the downward force F at the free end $\hat{s}=\hat{L}$. (b) Comparison between the Timoshenko analytical y (black lines) and numerical yⁿ (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^{\mathrm{\infty }}(ϵ)$ (black), L¹(ϵ) (blue) and L²(ϵ) (red) of the error ϵ=∥y−yⁿ∥ 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⁻²δl seconds. For all studies, we used the following settings: rod density ρ=5000 kg m⁻³, Young’s modulus E=10⁶ Pa, shear modulus G=10⁴ Pa, shear/stretch matrix $\widehat{\mathbf{\text{S}}}=\text{diag}(4G\hat{A}/3,4G\hat{A}/3,E\hat{A})\hspace{0.17em}\mathrm{N}$, bend/twist matrix $\widehat{\mathbf{\text{B}}}=\text{diag}(E{\hat{I}}_1,E{\hat{I}}_2,G{\hat{I}}_3)\hspace{0.17em}{\mathrm{Nm}}^2$, downward force F=15 N, rest length $\hat{L}=3\hspace{0.17em}\mathrm{m}$, rest radius $\hat{r}=0.25\hspace{0.17em}\mathrm{m}$, dissipation constant γ=10⁻¹ kg (ms)⁻¹, simulation time T_sim=5000 s.

Figure 6.

Muscular activity. Example of muscular activity amplitude profile (solid black line) described by cubic B-spline through N_m=8 control points $({\hat{S}}_i,{\beta }_i)$ with i=0,…,N_m−1. The control points are located along the filament at the positions ${\hat{S}}_i$, and are associated with the amplitude values β_i. The first and last control points are fixed so that $({\hat{S}}_0,{\beta }_0)=(0,0)$ and $({\hat{S}}_{N_m-1},{\beta }_{N_m-1})=(\hat{L},0)$, therefore assuming the ends of the deforming body to be free.

Figure 7.

Contact model with solid boundaries. Obstacles and surfaces (grey) are modelled as soft boundaries allowing for the interpenetration ϵ=r−d with the elements of the filament (blue) characterized by radius r and distance d from the substrate. The surface normal u^w_⊥ determines the direction of the wall’s response F^w_⊥ to contact. We note that F^w_⊥ balances the sum of all forces F_⊥ that push the rod against the wall, and is complemented by the other two components, which allow it to amend to possible interpenetration due to numerics. These components are an elastic one (k_wϵ) and a dissipative one (γ_wv⋅u^w_⊥), where k_w and γ_w are, respectively, the wall stiffness and dissipation coefficients.

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 ${\mathbf{\text{u}}}_{\parallel }^{\mathrm{w}}=\mathbf{\text{t}}\times {\mathbf{\text{u}}}_{\perp }^{\mathrm{w}}$ and a longitudinal one in the direction ${\mathbf{\text{u}}}_{\times }^{\mathrm{w}}={\mathbf{\text{u}}}_{\perp }^{\mathrm{w}}\times {\mathbf{\text{u}}}_{\parallel }^{\mathrm{w}}$. We note that in general u^w_×≠t. The notation x_⊥, x_∥, x_× denotes the projection of the vector x in the directions u^w_⊥, u^w_∥, u^w_×. (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.

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/\hat{L}\approx 1$) and R=4. (c) Formation of a solenoid for F=300 N (leading to the total elongation $L/\hat{L}\approx 1.15$) and R=13. Settings: length L=1 m, radius r=0.025 m, mass m=1 kg, Young’s modulus E= 10⁶ Pa, shear modulus G=2E/3 Pa, shear/stretch matrix $\widehat{\mathbf{\text{S}}}=\text{diag}(4G\hat{A}/3,4G\hat{A}/3,E\hat{A})\hspace{0.17em}\mathrm{N}$, bend/twist matrix $\widehat{\mathbf{\text{B}}}=\text{diag}(E{I}_1,E{I}_2,G{I}_3)\hspace{0.17em}{\mathrm{Nm}}^2$, dissipation constant γ=2 kg (ms)⁻¹, k_sc= 10⁴ kg s⁻², γ_sc=10 kg s⁻¹, discretization elements n=100, timestep δt=0.01δl s, T_twist=75 s, T_relax= 50 s.

Figure 10.

Optimal lateral undulation gait. (a–d) Instances at different times of a snake characterized by the identified optimal gait. (e) Evolution of the fitness function f=v^fwd_max 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 T_m. 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 ρ=10³ kg m⁻³, T_m=1 s, Young’s modulus E=10⁷ Pa, shear modulus G=2E/3 Pa, shear/stretch matrix $\widehat{\mathbf{\text{S}}}=\text{diag}(4G\hat{A}/3,4G\hat{A}/3,E\hat{A})$ N, bend/twist matrix $\widehat{\mathbf{\text{B}}}=\text{diag}(E{I}_1,E{I}_2,G{I}_3)\hspace{0.17em}{\mathrm{Nm}}^2$, dissipation constant γ= 5 kg (ms)⁻¹, gravity g=9.81 m s⁻², friction coefficient ratios ${\mu }_{\mathrm{k}}^f:{\mu }_{\mathrm{k}}^b:{\mu }_{\mathrm{k}}^r=1:1.5:2$ and ${\mu }_{\mathrm{s}}^f:{\mu }_{\mathrm{s}}^b:{\mu }_{\mathrm{s}}^r=1:1.5:2$ with ${\mu }_{\mathrm{s}}^f=2{\mu }_{\mathrm{k}}^f$, friction threshold velocity $v_{ϵ}={10}^{-8}\hspace{0.17em}\mathrm{m}\hspace{0.17em}{\mathrm{s}}^{-1}$, ground stiffness and viscous dissipation k_w=1 kg s⁻² and γ_w=10⁻⁶ kg s⁻¹, discretization elements n=50, timestep δt=2.5⋅10⁻⁵ T_m, 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 |β|^max_i=0,…,5=50 Nm.

Figure 11.

Optimal planar swimming gait at low Reynolds number. (a–d) Instances at different times of a filament swimming according to the identified optimal gait. (e) Evolution of the fitness function f=v^fwd_max 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 T_m. Red lines represent head and tail displacement in time. Settings: Reynolds number Re=10⁻⁴, length L=1 m, radius r=0.025 m, filament density ρ=10³ kg m⁻³, T_m=1 s, Young’s modulus E=10⁷ Pa, shear modulus G=2E/3 Pa, shear/stretch matrix $\widehat{\mathbf{\text{S}}}=\text{diag}(4G\hat{A}/3,4G\hat{A}/3,E\hat{A})$ N, bend/twist matrix $\widehat{\mathbf{\text{B}}}=\text{diag}(E{I}_1,E{I}_2,G{I}_3)\hspace{0.17em}{\mathrm{Nm}}^2$, dissipation constant γ=5 kg (ms)⁻¹, discretization elements n=50, timestep δt=2.5×10⁻⁵ T_m, 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 |β|^max_i=0,…,5=50 Nm.

Figure 12.

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

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 x≃s. 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 d_xy are related to each other so that ψ=θ+d_xy.

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 Nm². 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 $\parallel F_{\mathrm{R}}\parallel \sim \mathcal{U}(0,{10}^{-2})\hspace{0.17em}\mathrm{N}$. The occurrence of an instability is detected whenever the rod bending energy E_B>E_th with E_th=10⁻³ J. Settings: rod’s mass m_r=1 kg, twist stiffness β=2α/3 Nm², shear/stretch matrix S=10⁵×1 N, bend/twist matrix B=diag(α,α,β) Nm², radius r=0.025L m, discretization elements n=100 m⁻¹, timestep δt= 10⁻⁵ s, simulation time T_sim=10 s, dissipation constant γ=0.

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 $\parallel F_{\mathrm{R}}\parallel \sim \mathcal{U}(0,{10}^{-3})\hspace{0.17em}\mathrm{N}$. The occurrence of an instability is detected whenever the rod translational energy E_T>E_th with E_th=10⁻⁴ J. (b–c) Visualization of a Mitchell bucking instability event for Φ=10 and β/α=2. Settings: rod’s mass m_r=1 kg, length L= 1 m, bending stiffness α=1 Nm², shear/stretch matrix S= 10⁵×1 N, bend/twist matrix B=diag(α,α,β) Nm², radius r= 0.01L m, discretization elements n=50 m⁻¹, timestep δt= 10⁻⁵ s, simulation time T_sim=2 s, dissipation constant γ=0.

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 $A_v\sin (2\pi f_{v}t)$ to the free end, and characterized by rest length $\hat{L}$ which is extended to a final length $L=e\hat{L}$. (b) Comparison between analytical θ (black lines) and numerical θⁿ (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=T_v/10 within one loading period $T_v=f_v^{-1}$. (c) Norms $L^{\mathrm{\infty }}(ϵ)$ (black), L¹(ϵ) (blue) and L²(ϵ) (red) of the error ϵ=∥θ−θⁿ∥ 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⁻⁴δl seconds. For all studies, we used the following settings: rod’s density ρ=10 kg m⁻³, Young’s modulus E=10⁶ Pa, shear modulus G=2E/3 Pa, shear/stretch matrix $\widehat{\mathbf{\text{S}}}=\text{diag}(4G\hat{A}/3,4G\hat{A}/3,E\hat{A})\hspace{0.17em}\mathrm{N}$, bend/twist matrix $\widehat{\mathbf{\text{B}}}=\text{diag}(E{\hat{I}}_1,E{\hat{I}}_2,G{\hat{I}}_3)\hspace{0.17em}{\mathrm{Nm}}^2$, forcing amplitude A_v=10³ Nm, forcing frequency f_v=1 s⁻¹, dilatation factor e=1.05, rest length $\hat{L}=\sqrt{E/\rho }/(e{f}_v)\hspace{0.17em}\mathrm{m}$, rest radius $\hat{r}=0.5\hspace{0.17em}\mathrm{m}$, simulation time T_sim=2000 s. We enabled dissipation in the early stages of the simulations, letting γ decay exponentially in time to a zero value.

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 E_T=mv²/2 (analytical solution: black line, numerical solution: blue line) and rotational E_R=Jω²/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=10⁹ Pa, shear modulus G=2E/3 Pa, shear/stretch matrix S=10⁴×1 N, bend/twist matrix B=diag(EI₁,EI₂,GI₃) Nm², dissipation constant γ=10⁻⁶ kg (ms)⁻¹, gravity g=9.81 m s⁻², static and kinetic friction coefficients μ_s=0.4 and μ_k=0.2, friction threshold velocity $v_{ϵ}={10}^{-4}\hspace{0.17em}\mathrm{m}\hspace{0.17em}{\mathrm{s}}^{-1}$, ground stiffness and viscous dissipation k_w=10 kg s⁻² and γ_w=10⁻⁴ kg s⁻¹, discretization elements n=50, timestep δt= 10⁻⁶ s, simulation time T=0.5 s. (c) Rod set in motion by the external couple C_s on a horizontal plane. (d) Translational E_T and rotational E_R energies are plotted against the couple C_s. 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 E_T and rotational E_R 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⁻⁶ m s⁻¹.

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 F^f_long are larger than in the backward direction F^b_long. The plot illustrates the behaviour of the total rod’s translational energy E_T as a function of the force F applied. Blue and solid black lines correspond to, respectively, simulated and analytical E_T for a rod pulled forward. Orange and dashed black lines correspond to, respectively, simulated and analytical E_T for a rod pulled backward. Settings: length L=1 m, radius r= 0.025 m, mass m=1 kg, Young’s modulus E=10⁵ Pa, shear modulus G=2E/3 Pa, shear/stretch matrix S=10⁴×1 N, bend/twist matrix B=diag(EI₁,EI₂,GI₃) Nm², dissipation constant γ=10⁻⁶ kg (ms)⁻¹, gravity g=9.81 m s⁻², static and kinetic forward friction coefficients μ^f_s=0.8 and μ^f_k=0.4, static and kinetic backward friction coefficients μ^b_s=0.4 and μ^b_k=0.2, friction threshold velocity $v_{ϵ}={10}^{-4}\hspace{0.17em}\mathrm{m}\hspace{0.17em}{\mathrm{s}}^{-1}$, ground stiffness and viscous dissipation k_w=10 kg s⁻² and γ_w=10⁻⁴ kg s⁻¹, discretization elements n=50, timestep δt=10⁻⁵ s, simulation time T=0.25 s.

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 T_m 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 T_m 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=10⁷ Pa, shear modulus G=2E/3 Pa, shear/stretch matrix S=10⁵×1 N, bend/twist matrix B=diag(EI₁,EI₂,GI₃) Nm², dissipation constant γ= 10⁻¹ kg (ms)⁻¹, gravity g=9.81 m s⁻², static μ^f_s=0.2, ${\mu }_{\mathrm{s}}^{\mathrm{r}}={\mu }_{\mathrm{s}}^{\mathrm{f}}$, ${\mu }_{\mathrm{s}}^{\mathrm{b}}={\mu }_{\mathrm{s}}^{\mathrm{f}}$ and kinetic μ^f_k=0.1, ${\mu }_{\mathrm{k}}^{\mathrm{r}}={\mu }_{\mathrm{k}}^{\mathrm{f}}$, ${\mu }_{\mathrm{k}}^{\mathrm{b}}={\mu }_{\mathrm{k}}^{\mathrm{f}}$ friction coefficients in the isotropic case, static μ^f_s=0.2, ${\mu }_{\mathrm{s}}^{\mathrm{r}}=2{\mu }_{\mathrm{s}}^{\mathrm{f}}$, ${\mu }_{\mathrm{s}}^{\mathrm{b}}=20{\mu }_{\mathrm{s}}^{\mathrm{f}}$ and kinetic μ^f_k= 0.1, ${\mu }_{\mathrm{k}}^{\mathrm{r}}=2{\mu }_{\mathrm{k}}^{\mathrm{f}}$, ${\mu }_{\mathrm{k}}^{\mathrm{b}}=20{\mu }_{\mathrm{k}}^{\mathrm{f}}$ friction coefficients in the anisotropic case, friction threshold velocity $v_{ϵ}={10}^{-8}\hspace{0.17em}\mathrm{m}\hspace{0.17em}{\mathrm{s}}^{-1}$, ground stiffness and viscous dissipation k_w=1 kg s⁻² and γ_w=10⁻⁶ kg s⁻¹, discretization elements n=100, timestep δt= 10⁻⁵ s, muscular activation period T_m=1 s, wavelength ${\mathrm{\lambda }}_m=\mathrm{\infty }$, phase shift ϕ_m=0, torque B-spline coefficients β_i=0,…,5={0,10,15,15,10,0} Nm.

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.