Serpentine locomotion through elastic energy release

A model for serpentine locomotion is derived from a novel perspective based on concepts from configurational mechanics. The motion is realized through the release of the elastic energy of a deformable rod, sliding inside a frictionless channel, which represents a snake moving against lateral restraints. A new formulation is presented, correcting previous results and including situations never analysed so far, as in the cases when the serpent's body lies only partially inside the restraining channel or when the body has a muscle relaxation localized in a small zone. Micromechanical considerations show that propulsion is the result of reactions tangential to the frictionless constraint and acting on the snake's body, a counter-intuitive feature in mechanics. It is also experimentally demonstrated that the propulsive force driving serpentine motion can be directly measured on a designed apparatus in which flexible bars sweep a frictionless channel. Experiments fully confirm the theoretical modelling, so that the presented results open the way to exploration of effects, such as variability in the bending stiffness or channel geometry or friction, on the propulsive force of snake models made up of elastic rods.

FDC, 0000-0001-9045-8418; DM, 0000-0002-7375-671X; NMP, 0000-0003-2136-2396; ABM, 0000-0001-8902-9923; DB, 0000-0001-5423-6033 A model for serpentine locomotion is derived from a novel perspective based on concepts from configurational mechanics. The motion is realized through the release of the elastic energy of a deformable rod, sliding inside a frictionless channel, which represents a snake moving against lateral restraints. A new formulation is presented, correcting previous results and including situations never analysed so far, as in the cases when the serpent's body lies only partially inside the restraining channel or when the body has a muscle relaxation localized in a small zone. Micromechanical considerations show that propulsion is the result of reactions tangential to the frictionless constraint and acting on the snake's body, a counter-intuitive feature in mechanics. It is also experimentally demonstrated that the propulsive force driving serpentine motion can be directly measured on a designed apparatus in which flexible bars sweep a frictionless channel. Experiments fully confirm the theoretical modelling, so that the presented results open the way to exploration of effects, such as variability in the bending stiffness or channel geometry or friction, on the propulsive force of snake models made up of elastic rods.

Introduction
Since the beginning of the research on creeping locomotion, four different mechanisms for snake motion have been identified [1][2][3][4]: serpentine, concertina, side-winding and rectilinear. A special feature characterizes the first of the four movements, namely that the snake has to exert normal forces against projections from the ground 1 to induce longitudinal sliding, so that tangential frictional forces operating between the serpent and the terrain have to be minimized, because they simply oppose gliding forward. This model for locomotion, pioneered by Gray [1][2][3], is based on the muscular elasticity of the snake's body and allows us to explain how locomotion can occur in the absence of tangential frictional forces. The model is in agreement with the experimental observation that serpentine motion becomes impossible when a snake is confined to a straight or circular channel where variations in the elastic energy cannot occur. The mechanical model of Gray is made up of a chain of rigid pieces connected to each other with rotational springs (a set-up also employed for undulatory motion in a fluid with obstacles [8]), while a more refined model in which a snake is represented by an elastic rod (Euler's elastica; see, for instance, [9 -11]) has been introduced in an almost unknown article by Kuznetsov et al. [12], based on results presented in another forgotten article by Lavrentiev and Lavrentiev [13]. The model of a deforming elastic rod has been independently considered later in other studies [5,14,15], including sand-swimming [16] and soft robotics [17]. 2 Recently, a way to achieve motion through a release of elastic energy has been related to the generation of configurational or 'Eshelby-like' forces [24][25][26][27] and its understanding has led to the discovery of torsional locomotion [28]. A novel reconsideration of the mechanics of snake motion as modelled by Kuznetsov et al. [12] has evidenced that their solution is inconsistent at some points (to be detailed later) and does not explain the generation of forces along the body of the snake. Furthermore, an experimental direct validation of the rod model has never been attempted. Therefore, a generalization and a rigorous derivation is presented in this article for the mechanics of an elastic rod sliding into a perfectly smooth channel, together with its experimental validation.
The propulsive force is independently derived using two approaches: (i) an energy formulation and (ii) an approach based on the 'micromechanical' derivation of all forces acting on the system, enhanced by the integration of the equations of motion. These two approaches have different merits, so that only the knowledge of both provides a complete picture of serpentine motion through a perfectly smooth channel. In particular, the energy formulation explains how configurational forces provide propulsion from a release of elastic energy. On the other hand, the micromechanical approach shows that the perfectly frictionless channel provides propulsion by means of tangential localized reactions acting against the snake's body. This is a surprising behaviour, counter-intuitive from a mechanical point of view.
Results are also presented for situations in which the snake is only partially inserted in the channel (a case never considered before), while the rest of the body does not have any restraints for propulsion. It is shown that a snake idealized as an elastic rod can always exit, but never enter, a frictionless channel. Finally, circumstances are also considered for the first time in which a discontinuity is present in the channel curvature or in the bending stiffness of the elastic rod simulating the snake. The latter case may be achieved by the snake with a localized muscular relaxation along the body, which is also interesting for snake robots or in the design of a colonoscope or other endoscopes [29]. To provide a direct validation of the theoretical model, a channel has been designed and realized, with the shape of a clothoid spiral, in which friction has been reduced to a negligible amount through the use of roller bearings. A series of experiments carried out in the 'Instability Laboratory' of the University of Trento, Italy, on elastic rods with uniform and variable bending stiffness are shown to fully confirm the theoretical results (see also the video in the electronic supplementary material). Although the presented experiments are limited to a channel in the form of a clothoid spiral and to three different types of elastic rods, the methodology introduced allows exploration of more complex situations, for instance involving friction at the rod/channel contact or sophisticated variability in the snake's bending stiffness or channel geometry.
2. The propulsion of a rod within a channel 2.1. Formulation of the problem An elastic inextensible rod of length l is considered, rectilinear in its undeformed configuration, and with variable bending stiffness B(s), so that s is the curvilinear coordinate along the rod's axis, s [ [0, l ], and represents a non-inertial reference frame. For generality, the bending stiffness will be assumed to be possibly discontinuous at discrete points.
The elastic rod (or a part of it) is moving inside a curved, frictionless (rigid and smooth) channel with curvilinear length L. The position of the rod's left end is singled out by the curvilinear coordinate j(t) [ [ 2 l, L], which is a function of the time t and measures the distance from the left end of the channel (figure 1).
In figure 1, the possible stages of the rod sliding within the channel are sketched at increasing values of the configurational parameter j(t), as defined in table 1: -case (A), where the left part of the rod is unconstrained while the right end moves within the channel; -cases (B1) and (B2), where the rod is entirely within the channel in the former case ( possible only when l , L) and the rod has both ends external to the channel in the latter case ( possible only when l . L); -case (C), which is the opposite configuration of that described in case (A).
Denoting by a superimposed dot the time derivative, _ j(t) and € j(t) represent the velocity and the tangential acceleration of the rod's left end and also of all the points along the rod, because the rod is inextensible. The curvature k of the rod axis depends on the global coordinate s(s, t) ¼ j(t) þ s and corresponds to the curvature x of the channel at each point s of the rod within the channel, while it is assumed null for the part of the elastic rod outside the channel, namely where s 1 and s 2 define, respectively, the initial and final curvilinear coordinate of the part of rod constrained by the channel, so that the four cases shown in figure 1 correspond to the pairs fs 1 (t), s 2 (t)g reported in table 1.
According to the Euler -Bernoulli model, the elastic rod reacts to bending with a moment M linearly proportional to the curvature k, so that M(s, j(t)) ¼ B(s)k(s þ j(t)). Any part of the elastic rod outside the channel is simply unloaded and therefore remains straight.

Energy approach: elastic energy turning into kinetic
At the generic time t of the motion of the inextensible elastic rod inside the channel, the system is characterized by the elastic flexural energy E(t) stored within the rod, and the kinetic energy T (t) of the rod, which (assuming negligible rotational inertia) can be written as where m ¼ Ð l 0 g(s) ds is the total mass of the rod and g(s) is its linear mass density.
Since the rigid channel is frictionless, the total potential energy of the system, 3 which leads to Defining the propulsive force P(t) as Note from equation (2.7) that the propulsive force P is a 'resultant configurational force', which is different from zero only when a change in the elastic strain energy E of the rod occurs through a change in the configurational parameter j, a situation similar to the Eshelby-like forces in elastic systems constrained by a sliding sleeve [26].
The two types of jump discontinuity are now introduced within the part of the rod constrained by the channel. One discontinuity is in the bending stiffness B(s) of the elastic rod at the local coordinate s ¼ l 1 [ (s 1 , s 2 ) and the other discontinuity is in the curvature of the channel at the global coordinate where the dependence on time has been omitted for simplicity (and also in the following). It is noted that the superscripts þ and 2 define, respectively, the one-sided limit of the relevant quantity from the negative or positive direction.
Considering that the elastic energy (2.3) can be rewritten as which is an integral with moving limits s 1 (j) and s 2 (j), and that the derivatives of the curvature field x(j þ s) with respect to j and s are identical, BðsÞxðj þ sÞx 0 ðj þ sÞ ds; ð2: 13Þ  Table 1. Ranges of the configurational parameter j and values of the coordinates s 1 and s 2 , together with their respective derivatives, for the four cases.
or, equivalently, after integration by parts, as Equations (2.13) and (2.14) are fully equivalent and provide the same value of the propulsive force P for the given geometry and stiffness properties. However, the two formulae are expressed in a way that only one type of discontinuity appears. Indeed, the former (the latter) expression highlights the influence of only jump discontinuities in the squared channel curvature ½½x 2 (in the bending stiffness ½½B), while the discontinuity in the stiffness (in the channel curvature) does not explicitly appear. Specifying equation (2.14) for the four cases shown in figure  1, and considering the values of the derivatives @s j (j)/@j ( j ¼ 1, 2) as computed in table 1, the propulsive force P simplifies to BðsÞxðj þ sÞx 0 ðj þ sÞ ds, ð2:15Þ or, equivalently, ð2:16Þ -case B1 BðsÞxðj þ sÞx 0 ðj þ sÞ ds, ð2:17Þ or, equivalently, BðsÞxðj þ sÞx 0 ðj þ sÞ ds, ð2:19Þ or, equivalently, BðsÞxðj þ sÞx 0 ðj þ sÞ ds, ð2:21Þ or, equivalently, ð2:22Þ Equations (2.15)-(2.22) provide the theoretical evaluation of the propulsive force P that may be analytically or numerically computed for given distributions of bending stiffness B(s) along the rod and curvature of the channel x(s). It is worth mentioning that only the case B1 was previously analysed in [12] under several restrictive hypotheses (continuous channel curvature and rod bending stiffness, null bending moment and axial force at the end of the rod) so that equation (2.18) reduces to the result presented in [12] when ½½B(l 1 ) ¼ 0. Note that equations (2.7) and (2.10) of [12] have incorrect signs. In addition, the normal forces internal to the snake have been erroneously disregarded in equation (2.7) of [30].

Propulsion from Newton's second law and micromechanics
Equations (2.13) and (2.14) provide the propulsive force in a global sense. This force is the result of the channel tangential reactions against the rod (which could be thought to be null at a first-mistaken-glance because the channel is frictionless), the knowledge of which is fundamental in the determination of the stress state in the snake's body during motion. These reactions, which have until now only been used to analyse the special case of the concentrated force at the end of the rod within a channel [5], can be derived with consideration of the dynamics of the rod's element and of the micromechanics of the elastic rod at its ends, at the entrance or the exit of the channel and at discontinuity points in the channel curvature or the rod's bending stiffness. The equilibrium equations for an elementary portion of a planar rod in a curved configuration (figure 2), defined by the curvature k(s), can be written in the non-inertial frame of reference attached to the rod (at its left end, the origin of the curvilinear coordinate s) as where p(s), q(s) and m(s) represent the channel reactions applied to the rod, respectively, in the transverse, axial and rotational directions; -for the part of the rod outside the channel, this is assumed to remain straight (flexural vibrations are disregarded), k(s) ¼ 0, so that the distributed forces are given by s [ ½0,s 1 ðtÞ < ½s 2 ðtÞ, l: ð2:25Þ In this case, considering null boundary conditions at the ends, the integration of the equations of motion leads to It is worth remarking that, in the case of quasi-static motion (namely, negligible inertial density forces), the distributed forces f r (s), f u (s), f s (s) would all be either null or coincident with the channel reactions p(s), m(s), q(s), respectively, for the parts of the rod outside and within the channel.

Micromechanics of a rod inside the channel
The tangential reactions acting at the rod's ends, at the channel entrance/exit, at the points of discontinuity of the channel curvature or of the rod's bending stiffness are derived here. These counter-intuitive reactions are strictly related to the deformability of the constrained rod and to the nonrectilinearity of the channel and have been found in other 'non-classical' situations [26,28,31].
Relaxing for a moment the assumption that the rod is perfectly tight to the channel, a clearance is introduced in the channel constraining the elastic rod as shown in figure 3. Owing to the presence of the gap, assumed small, between the elastic rod and the channel, the former always has a zone of detachment (singled out by the angle F and having arc-length a) from the latter. In the micromechanical analysis developed below it is shown that in the limit when the gap vanishes, so that the detachment region vanishes too, a nonnull tangential reaction is developed by the frictionless channel.
In figure 3, three possible cases of interest are sketched, namely: (i) the rod's end is inside the channel; (ii) the channel curvature has a discontinuity; and (iii) the rod's end is outside the channel. In all cases, the reaction force F is orthogonal to the channel, which is frictionless. Below, these three cases will be analysed together with the further case in which (iv) the rod has a discontinuity in the bending stiffness.

Tangential reaction at the rod end inside the channel
The situation when one end of the elastic rod is inside the frictionless channel with radius of curvature 1/x is sketched in figure 3a.
It is evident that the rod end is subject to the channel reaction F only, orthogonal to the upper contour of the channel. This force provides the moment M ¼ FF/x at the point (at the distance a ¼ F/x from the rod's end) where the elastic rod becomes adherent to the profile, so that M ¼ Bx. At that point the internal axial action can be approximated as F sinF % FF, which becomes B x 2 , the resultant axial force developing at the rod's end within a channel.
With reference to the curvilinear coordinate system s, if the left and the right ends of the rod are inside the channel, s 1 ¼ 0 and s 2 ¼ l, the compressive reaction forces are realized by the channel at the edges of the rod cases B1 and C, and N(l À ) ¼ ÀB(l)x 2 (j þ l), cases A and B, ð2:27Þ which have also been obtained using the principle of virtual power in [5] and are similar to the tangential reaction developing in a perfectly frictionless curvilinear profile constraining a clamped end of the rod [31].

Tangential reaction at a discontinuity point in the channel curvature
In the case when there is a jump in the curvature along the channel (figure 3b), the curvature of the rod takes a value different from the channel curvature because a small gap is present. In particular, the rod has a transition length along which its curvature passes from the channel curvature before to that after the jump. With reference to figure 3b, two concentrated contact reactions F 2 and F þ develop orthogonal to the upper and lower contours of the channel at the rod's two detachment points. Owing to a difference in the angle F in the direction of the two reactions, a tangential resultant is expected. At equilibrium, in the limit of vanishing detachment length a, the two reactions orthogonal to the channel tend to the same value, F þ ¼ F 2 ¼ F, and the tangential resultant is given by F sinF % FF. The bending moments at the two detachment points are defined by the relevant curvature imposed by the channel, so that they are given by M 2 ¼ Bx 2 and M þ ¼ Bx þ . With reference to the linear elastic solution of a clamped rod of length a and loaded at its free end by both a force F and a moment M 2 , the rotation angle at the loaded edge is F ¼ Fa 2 /(2B) þ M 2 a/B. Since it follows by equilibrium that Fa ¼ M þ 2 M 2 , the tangential resultant FF due to the jump in curvature ½½x provides a concentrated tangential reaction equal to B½½x 2 =2. Therefore, with reference to a jump in curvature occurring at the point s ¼ L 1 2 j, the presence of the abovementioned tangential reaction provides a jump in the internal axial force corresponding to Tangential reaction at an exit point from the channel The analysis of the previous case is also relevant to the condition when a rod's end is outside the channel (figure 3c). Indeed, in this case x 2 ¼ 0 and, defining x þ ¼ x, the tangential resultant of the contact forces becomes Bx 2 /2. In the special case of partial insertion of the rod into the channel, s 1 = 0 and s 2 = l, the tangential resultant Bx 2 /2 provides the following jumps in the internal axial force: , cases A and B2, and ½½N(s 2 ) ¼ B(s 2 )x 2 (j þ s 2 ) 2 , cases B2 and C: ð2:29Þ Note the following.
-The presence of the factor 1/2 distinguishes equation (2.29) from equation (2.27), which holds for the case of a rod entirely within the channel, s 1 ¼ 0 and s 2 ¼ l. -The tangential reactions (2.29) have an expression similar to, but different in sign from, the Eshelby-like force generated at a sliding (straight) sleeve providing the reaction moment M ¼ Bx [26].

Tangential reaction at a point of discontinuity in the rod stiffness
At the point s ¼ l 1 , where a jump in the bending stiffness ½½B(l 1 ) is present, a tangential concentrated contact reaction ½½B(l 1 ) x 2 is developed in order to satisfy the rotational equilibrium of the infinitesimal part of the rod taking the centre of the osculating circle as the pivot point. This tangential reaction provides the following jump in the internal axial force: Considering now the possibility of discontinuities in the stiffness and in the curvature, respectively, at the coordinates l 1 and L 1 2 j (equations (2.9) and (2.10)), integration of the longitudinal inertia density (equation (2.31)), over the length of the rod leads to the determination of the propulsive force P as P ¼ Nðl À Þ À ½½Nðs 2 ÞHðl À s 2 Þ À ½½Nðl 1 Þ À ½½NðL 1 À jÞ À ½½Nðs 1 ÞHðs 1 Þ À Nð0 þ Þ þ Bðs 2 Þx 2 ðj þ s 2 Þ 2 ½qðsÞ þ mðsÞxðj þ sÞ ds, ð2:32Þ where Although at first sight the propulsive force P given by equation (2.32) may appear different from that obtained through the energy approach (equation (2.14)), the two expressions are identical. Indeed, as shown through micromechanical concepts, (i) the jumps in the internal axial action ½½N(L 1 À j) and ½½N(l 1 ) are given by equations (2.28) and (2.30), respectively, as well as (ii) the following identities (to be compared with equations (2.27) and (2.29)) hold: and N(l À ) À½½N(s 2 )H(l À s 2 ) ¼ À 1 þ 1 2 @s 2 (j) @j B(s 2 )x 2 (j þ s 2 ), ð2:34Þ The cases in which either the channel curvature or the rod stiffness is constant are considered.
In the case that a rod, with non-constant stiffness B(s), is constrained by a channel with constant curvature x(s) ¼ x, the propulsive force (2.14) reduces to while in the case that a rod with constant stiffness B(s) ¼ B is inside a channel, with non-constant curvature x(s), the propulsive force (2.14) reduces to :

ð3:2Þ
With reference to the possible four cases sketched in figure 1, the propulsive force simplifies as reported in table 2.
The evaluation of the propulsive force P in the different cases shows that: -in the case when the rod has only one end outside the channel, an outward propulsive force is always realized in order to eject the rod from the channel on the side where the rod is initially out, namely € j < 0 for case A and € j . 0 for case C; -in the case of a channel with constant curvature, the propulsive force is null, P ¼ 0, when the rod is completely inserted in the channel (case B1) so that the rod does not move; if both ends of the rod are outside the channel (case B2), an outward propulsive force is realized ejecting the rod away from the exit point of the channel where the rod has the highest bending stiffness; -in the case of a rod with constant bending stiffness, the propulsive force is null, P ¼ 0, when the rod has both ends outside the channel (case B2) so that the rod does not move; in the case when the rod is completely inserted in the channel, the propulsive force moves the rod away from the end deformed with the highest (in absolute value) curvature and towards the end deformed with the smallest one; -in the case of a rod with constant bending stiffness constrained inside a channel with constant curvature the propulsive force is null, P ¼ 0, whenever both the ends are constrained inside the channel (case B1) or both are free (case B2).

The two-piecewise constant rod stiffness inside a two-piecewise constant channel curvature
Let us now consider the case of both rod stiffness and channel curvature described by two-piecewise constant functions as follows: where l 1 þ l 2 ¼ l, L 1 þ L 2 ¼ L and for simplicity it is assumed that fL 1 , L 2 g . l.
Restricting attention to the case of a rod completely inserted inside the channel (case B1 in figure 1), namely j [ (0, L 2 l ), four possible subcases depending on the number of configurational parameters j(t) arise. The four subcases are sketched in figure 4 and listed in table 3 together with the elastic energy E of the system and the propulsive force P evaluated in each case.
It can be concluded that -the propulsive force is null, P ¼ 0, in both subcases (I) and (IV) because a small perturbation in the configuration does not provide any change in the elastic energy of the system; this also occurs in subcases (II) and (III) in the special case when the two curvatures assume the same absolute value, jx 1 j ¼ jx 2 j; -a non-null propulsive force is generated in subcases (II) and (III) in the direction of a decreasing amount of rod constrained by the highest curvature (in absolute value).
With reference to the case jx 1 j , jx 2 j and B 1 , B 2 , the propulsive force P and the strain energy E are reported in figure 4 as functions of the configurational parameter j. The value =2 is introduced to define the strain energy at the transition from subcase II to subcase III. Note that subcases II and III provide a negative propulsive force, so that the motion is characterized by _ j(t) < 0, towards minimization of the elastic energy.

Experiments on snake propulsion with an elastic rod
The experimental validation of the theoretical modelling presented in the previous sections was provided by employing Table 2. Expression for the propulsive force P for the four possible cases (A, B1, B2 and C; figure 1) under the assumption of a channel with constant curvature x (first row) or of a rod with constant bending stiffness B (second row).
rsif.royalsocietypublishing.org J. R. Soc. Interface 14: 20170055 the experimental set-up illustrated in figure 5. A frictionless channel, of length L ¼ 750 mm, was designed and manufactured at the Instability Laboratory in the shape of a clothoid spiral. This spiral, used in road design to provide the best smooth transition between different tracks, was selected to be the ideal candidate for experiments because its curvature x changes linearly with the arc length s, so that x(s) ¼ p s/a 2 , which is the reason why the clothoid is considered the basic form assumed by a snake during its motion [30]. This curve is parametrically expressed through Fresnel integrals as [32] x(s) ¼ a where a represents a scaling parameter. In the design of the smooth channel, a scaling parameter a ¼ 510 mm was assumed, so that the curvature varies between 0 and 0.009 mm 21 .
The two sides of the channel between which the elastic rod slides were realized with pairs of superimposed rollers (Misumi Europe; press-fit straight type, 20 mm diameter and 15 mm length), separated by a 3 mm thick Teflon (PTFE) sheet, which was introduced to reduce the discontinuities created by the gaps between the rollers (having rotation axis spaces of 21.62 mm). Two pairs of 66 and 76 rollers were used, respectively, for the short and long side of the channel, as shown in figure 5.
Different elastic rods were shaped by cutting (with an MDX-540 milling machine; Roland) 'plain' carbon fibre strips (E ¼ 80148 MPa, r ¼ 1620 kg m 23 ), 2.05 + 0.05 mm thick, 40 mm wide and 2000 mm long. A tolerance of 0.5 mm was kept between the carbon rod and the rollers along the channel.

Ejection experiments
Qualitative experiments were performed by inserting (by hand) a rod of uniform bending stiffness into the channel and then suddenly releasing it. The release of elastic energy produces a violent upward ejection of the rod from the channel.
For a carbon-fibre rod of dimensions b ¼ 16.2 mm, t ¼ 1.35 mm, l ¼ 660 mm, a mean exit velocity (calculated considering two subsequent images, taken with a high-speed camera (Sony PXW-FS5) just when the rod was entirely outside the channel) of 8.5 + 0.5 m s 21 was estimated. The theoretical value of this upward velocity was calculated in the absence of dissipation to be equal to 12.1 m s 21 . In addition to the vibration trapped inside the elastic rod after ejection (which is observed to have a small amplitude), there are two sources of dissipation explaining the difference between the data. The first, and less important, is due to the friction inside the channel which is not exactly null, but remains very small, because the rollers perform very well and the elastic rod was oiled. 5 The second source is due to the fact that the passage of the rod induces a fast rotation in the rollers, so that in the channel prototype a non-negligible kinetic energy remains, stored as rotational inertia of the rollers.
A sequence of frames (exported from a high-speed movie) is reported in figure 6 during a test (see the movie in the electronic supplementary material) showing how quickly the energy can be released as well as the efficiency of the propulsion.

Quantitative experiments
For a quantitative evaluation, the propulsive force P was measured by extracting the left edge (s ¼ 0) of the rod from the channel at a constant low speed _ j ¼ À2 mm s 21 , using a loading frame (Midi 10; Messphysik Materials Testing) placed horizontally and simultaneously acquiring the load (Mettler MT1041 load cell RC 200N connected to a NI CompactRio acquisition system).
With the purpose of separately verifying the different contributions to the propulsive force, three rods were tested that differ in the bending stiffness distribution B(s) provided by varying the width b(s), namely: a rod with constant stiffness (see figure 5c), a rod with a localized jump in the stiffness (see figure 5d) and a rod with a linearly varying stiffness (see figure 5e). Rod with a localized jump in stiffness. The propulsive force has been measured when a rod with a finite jump in the bending stiffness (varying from B 1 ¼ 1.985 Â 10 6 mm 4 , provided by b 1 ¼ 34.5 mm, to B 2 ¼ 0.460 Â 10 6 mm 4 , provided by b 2 ¼ 8 mm, at the coordinate l 1 ¼ 1.1L) sweeps the channel (case B2 in figure 1), The propulsive force is theoretically predicted by equation (2.20) assuming B 0 (s) ¼ 0 and taking ½½B( Rod with a linearly varying stiffness. The propulsive force generated during the movement of an elastic rod with a linearly varying bending stiffness (from B 1 ¼ 2.134 Â 10 6 mm 4 , provided by b 1 ¼ 34.5 mm, at the coordinate s ¼ l 1 ¼ 1.1L, to B 2 ¼ 0.112 Â 10 6 mm 4 , provided by b 2 ¼ 2 mm, at the coordinate s ¼ l 1 þ d, with d ¼ 750 mm) has been measured (case B2 in figure 1), The propulsive force can be computed from equation (2.20) assuming and taking ½½B(l 1 ) ¼ 0 as Although the channel is idealized as rigid and frictionless in the theoretical approach, a minimal amount of friction has to be expected in any practical realization. This friction has been measured by pulling through the channel a rod of constant bending stiffness, as shown in scheme B2 in figure 1, and finding that a non-null force was present (which was measured to oscillate between 3 and 4 N), instead of being null as theoretically predicted. This determination of the channel friction is important in our understanding of the experiments. The results of the experiments describing the three cases of constant, discontinuous and linearly varying bending stiffness are, respectively, reported in figure 7a, 7b and 7c.
In figure 7, the propulsive force P (made dimensionless through multiplication by L 2 and division by B 1 ) is reported versus the configurational parameter j (made dimensionless through division by the channel length L). Note that the scale of the vertical axis in figure 7c is different from that in figure 7a and figure 7b.
In all the three cases, the graphs have to be read from right to left, because the tests were performed at decreasing j. The experiments have been designed to initiate from a situation of steady state (initiating from j/L ¼ 0 and terminating at j/L ¼ 2 0.1), in which the propulsive force is theoretically null, P ¼ 0.
In the two cases reported in figure 7a,b, the propulsive force jumps to a finite value when the rod's end (figure 7a) or the rod's discontinuity in B (figure 7b) enters the channel, j/L ¼ 2 0.1. Finally, at the exit from the channel, j/L ¼ 2 1.1, the propulsive force smoothly vanishes.
In the case reported in figure 7c, at decreasing j/L, the propulsive force P increases up to a maximum value attained at j/L ¼ 2 1.1 (and then would decay to zero, if the experiment were continued, which it was not). All theoretical Table 3. Elastic energy E and propulsive force P for the four possible subcases (I, II, III and IV) of a two-piecewise constant bending stiffness rod completely inserted within a two-piecewise constant curvature channel ( figure 4).
rsif.royalsocietypublishing.org J. R. Soc. Interface 14: 20170055 predictions are reported with black curves and all experimental values with blue spots ( joined with orange segments).
Although an immediate qualitative and quantitative agreement between theoretical predictions (black curves) and experimental values (marked in red) can be appreciated from figure 7, the discrepancies can be easily classified and explained as follows.
-The propulsive force displays an oscillation which is the consequence of the gaps between the rollers, an effect only partially alleviated by the Teflon layer. Note that the oscillations are more evident in the cases where either the rod's end (figure 7a), or the localized discontinuity in the rod's thickness (figure 7b) sweeps the channel and hits the rollers (see the detail in figure  5a). In the case when the rod has a continuous change in thickness b (figure 7b), the oscillations are strongly reduced (taking into account the different scale in the vertical axes of the graphs in the figure).
-The propulsive force is slightly reduced with respect to the expected theoretical value, even in the steady-state case, where it is negative rather than being null. This is the direct consequence of the small, but not null, friction in the channel that after normalization yields a spurious force PL 2 /B 1 ¼ 1.5 + 0.5, almost exactly corresponding to the observed shift in the data. -The propulsive force in figure 7a,b approaches zero faster than the theoretical prediction near j/L ¼ 2 1. This is simply due to the tolerance of 0.5 mm between the elastic rod and the rollers, so that the elastic rod reaches its undeformed configuration before it complete exits from the channel.
It can be concluded that the model of snake motion is definitively confirmed by experiments and that the propulsive force can be directly measured with a quasi-static electromechanical testing machine. 5. Conclusions, with a discussion on modelling snake motion A model of serpentine motion, in which an elastic rod (the snake) moves inside a frictionless channel (the lateral restraints provided by projections from the ground or by transverse friction) through a release of (muscular) elastic energy, has been reviewed, corrected and extended, to cover situations occurring when the snake's body lies only partially inside the restraining channel and to explain tangential contact forces (so far ignored, but crucial to the propulsion). The presented theoretical results represent an extension and a clarification based on the concept of configurational force of the mathematical models introduced by Gray [1-3] and developed in several later contributions [5,7,13,15,30]. An experimental set-up has been designed, realized and tested, which allows for the first time direct measurement of the propulsive force in different situations, also involving a jump in the rod's bending stiffness (which a snake can obtain through a relaxation of muscles over a small zone of its body). Experimental results have been found to strongly confirm the theory, so that the developed experimental setup is now available to directly measure the propulsive force and (4.6), respectively). The oscillation of the load is due to the gaps between the rollers; the fact that the experimental values are smaller than the predictions is due to the effect of friction, which is small, but not null. Right: sketch of the initial (j/L ¼ 0) and final (j/L ¼ 2 1.1) configurations for the three cases of variation in the bending stiffness considered for the experiments (for simplicity, only a few rollers aligned on a straight line are reported). (Online version in colour.) on 'prototypical snakes' made up of elastic rods of different materials, stiffnesses and geometries. In this way difficulties related partially to mathematical modelling and partially to the use of living serpents, which are difficult to handle and to test, can be bypassed. The experiments can in fact be extended to situations involving friction (which is difficult to treat mathematically, because the transverse reactions in the channel are not known), or lubrication (the experimental equipment can be submerged in water or other liquids), or different stiffnesses (not only in the snake's body, but also in the channel transverse stiffness), and complex geometries of both the channel and the snake. Imagine, for instance, a snake moving against obstacles in an environment including partially submerged areas: using a physical model of such a situation, the developed experimental set-up would allow direct measurement of the propulsive force.
Data accessibility. Calculations were performed analytically. Experiments are reported in the figures and in the electronic supplementary material. Further data will be provided upon request.