Planar model of friction modulation

We adopt the approach of Hu et al.^13,14 wherein the centerline of a snake of length L is modeled as an inextensible planar curve s^∈[0,L]. The center of mass position and average orientation of the snake are denoted by x¯^(t) and α¯(t), respectively, and the local position and orientation of each point along the snake’s centerline is computed via x^(s^,t^)=x¯^(t^)+I[t(s^,t^)] and α(s^,t^)=α¯(t^)+I[κ^(s^,t^)], respectively, where t(s^,t^)=(cosα(s^,t^),sinα(s^,t^)) is the local tangent vector, κ^(s^,t^) is the local curvature and I[f(s^,t^)]=∫0s^f(s^′,t^)ds^′−1L∫0L∫0s^f(s^′,t^)ds^′ds^ is a mean-zero integration function, which expresses the mathematical machinery that allows us to reconstruct snake’s local positions/orientations from the center of mass, global orientation and curvature information (Fig. 1c). The dimensionless form of κ^(s^,t^) is defined in the main text with all simulations employing ϵ = 7 and k = 1, as in¹³. Differentiating x^(s^,t^) twice with respect to time yields

where n=(−sinα,cosα) is the local normal vector. Writing the snake’s dynamics as a force balance of internal f^ and external F^ forces per unit length yields

where ρ is the line density of the snake. We then scale Eqs. (¹) and (²) by s=s^/L and t=t^/τ to non-dimensionalize the system.

External forces stem entirely from frictional effects captured through the Coulomb friction model, with anisotropy characterized by coefficients in the forward (μ_f), backward (μ_b), and transverse (μ_t) directions. Scaling friction forces such that F=F^/ρgμf allows us to write the friction force as F(s, t) = −N(s, t)μ(s, t) with 𝛍(s,t)=μtμf(u⋅n)n+[H(u⋅t)+μbμf1−H(u⋅t)](u⋅t)t, where u(s)=xt(s)/xt(s) is the unit vector associated with the snake’s local velocity direction, and H = 1/2(1 + sgn(x)) is the Heaviside step function used here to distinguish between forward and backward friction components. Here, μ_b/μ_f = 1.5, in keeping with experimental observations¹³. Previous work¹⁴ and our preliminary investigations found that static friction effects do not appreciably influence the snake’s steady-state behavior, thus we did not consider them here. Moreover, we write the friction modulation wave as N(s,t)=ηN^(s,t), where η=1/∫01N^(s,t)ds is the normalization constant to conserve the overall weight of the snake. Finally, we set Fr = 0.1 throughout, consistent with snakes’ typically low values and without lack of generality.

Assuming the total non-dimensionalized internal forces and torques to be zero (∫01fds=0 and ∫01(x−x¯)×fds=0) yields the snake’s equation of motion

where J=∫01(x−x¯)2ds is the moment of inertia. These equations can then be solved for a prescribed non-dimensional curvature κ(s, t) and friction scaling term N(s, t).

In all cases considered here, Eqs. (³) and (⁴) are numerically solved over 10 undulation periods to allow transient effects from startup to dissipate and the snake to reach steady-state behavior. The snake’s locomotion behavior is then analyzed in terms of the pose angle γ, steering rate θ˙, and effective speed ∣v_eff∣ which are illustrated in Fig. 2a. At steady state, the first trajectory metric that can be computed is the pose angle, which is the angle between the snake’s average orientation t¯=(cosα¯,sinα¯) and its center of mass velocity direction u¯. The average pose angle over one undulation period is defined as γ=∫t0t1arctan2((t¯×u¯)⋅ez,t¯⋅u¯)dt/T where T=∫t0t1dt and e_z is the unit vector of out-of-plane axis. Use of the arctan2 function is required to ensure γ ∈ (−π, π]. Note that T is for a nondimensionalized time period, so over one undulation, T=1. Additional trajectory metrics can be computed by considering the snake’s center of mass as a particle undergoing planar motion in polar coordinates, x¯(t)=(rcosθ,rsinθ), allowing the snake’s trajectory to be quantified in terms of its effective velocity ∣veff∣=∫t0t1rdθdtuθdt/T and steering rate θ˙=∫t0t1dθdtdt/T (see SI Note ¹ for relevant derivations).

For phase space simulations, a simulation grid was defined with 501 equidistant points in both A ∈ [−2, 2] and Φ ∈ [0, 1], leading to 251k simulations for each of the friction ratios considered. Simulations were performed on the Bridges supercomputing cluster at the Pittsburgh Supercomputing Center.