A Derivation of the equations of motion and small-amplitude analysis

Optimization of two- and three-link snake-like locomotion


We analyze two- and three-link planar snake-like locomotion and optimize the motion for efficiency. The locomoting system consists of two or three identical inextensible links connected via hinge joints, and the angles between the links are actuated as prescribed periodic functions of time. An essential feature of snake locomotion is frictional anisotropy: the forward, backward and transverse coefficients of friction are different. The dynamics are studied analytically and numerically for small and large amplitudes of the internal angles. Efficiency is defined as the ratio between distance traveled and the energy expended within one period, i.e. the inverse of the cost of locomotion. The optimal set of coefficients of friction to maximize efficiency consists of a large backward coefficient of friction and a small transverse coefficient of friction, compared to the forward coefficient of friction. For the two-link case with a symmetrical motion, efficiency is maximized when the internal angle amplitude is approximately , for transverse coefficient sufficiently large. For the three-link case, the efficiency-maximizing paths are triangles in the parameter space of internal angles.

1 Introduction

Snake locomotion has long been a topic that fascinates researchers, and has recently received a renewed wave of interest in the fields of robotics and control, as well as in organismal biology. Snakes are familiar organisms, but as limbless animals, their locomotion has special features [7]. Terrestrial snakes move using friction between the ground and their belly scales, which have anisotropic frictional properties [8]. It has been proposed that the cost of transport (energetic efficiency) for snake slithering is no greater than that of limbed animals [23, 27].

Some works on modeling snake locomotion are oriented towards wheeled snake robots [3, 4, 6, 9, 11, 19, 21]. These models are typically concerned with motion planning, and assume that the transverse coefficient of friction is high enough to prevent transverse motion, while the forward and backward coefficients of friction are the same. These models work well for wheeled robots and provide valuable insights into the locomotion of biological snakes. In experiments, Hu et al. [12] measured the frictional anisotropy of juvenile milk snakes, and found that the forward, backward and transverse coefficients are different but similar in magnitude. They also studied the effects of the snake’s active modulation of its weight distribution on the ground. Hu & Shelley [13] analyzed the so-called “lateral undulation” motion modeled as a family of sinusoidal traveling-wave shapes and calculated the dependence of speed and efficiency on amplitude and wave length of the kinematics as well as coefficients of friction.

In this work, we adopt the model of [12], but for bodies composed of two and three rigid links. Linked bodies are fundamental for robotic sliding systems. By specializing to bodies with two and three links, we consider the simplest such systems, which nonetheless have nontrivial behaviors. Two- and three-link locomoting bodies have been considered previously as swimmers at zero Reynolds number (Stokes flow). Purcell [22] described the physics of swimming in Stokes flow, and stated the Scallop Theorem: in Stokes flow, net locomotion is not possible if a swimmer deforms in a way that is invariant under the reversal of time [5]. Such is the case for periodic motions of a two-link body (“scallop”). He then proposed a three-link swimmer that moves only one link at a time, in a non-reciprocal motion that results in net locomotion. Subsequent studies have calculated efficiency-optimizing motions for Purcell’s swimmer and similar systems. Becker et al. [2] calculated efficiency-optimizing stroke amplitudes for Purcell’s swimmer, and considered different length ratios of the three links. Tam & Hosoi [25] extended the optimization to arbitrary kinematics (allowing both internal angles to change simultaneously) and arbitrary slenderness ratios. They found an optimal path in the parameter space of internal angles using a Fourier series representation, and showed that the high frequency modes are subdominant to the low frequency modes. Avron & Raz [1] developed a qualitative geometric approach by focusing on the curvature of the local connection matrix to study the same system. Hatton & Choset [10] further developed this technique, and suggested a systematic way of choosing the best body-fixed frame to approximate the inertial frame displacement while accounting for the overall rotation. They calculated optimal motions for other systems such as a three-link fish swimming in infinite Reynolds number (potential flow), which admits a similar formulation. Kanso et al. [15] and Melli et al. [18] gave a geometric formulation for swimming in a potential flow and calculated optimal strokes. Jing & Kanso [14] used this formulation to study the effects of elasticity and body configuration on the stability of passive locomotion.

Here we study the locomotion of two- and three-link bodies not in a fluid but instead on a planar surface with sliding friction. Unlike the fluid studies, the dynamical equations are nonlinear with respect to body velocities, so the Scallop theorem and many other results no longer apply. Like the aforementioned studies, we focus on finding motions which optimize the efficiency of locomotion. The structure of the paper is as follows. In Section 2 we formulate the model for a two-link “snake,” nondimensionalize the system and derive the equations of motion. In Section 3 we use analysis and computation to optimize the two-link model with respect to kinematic parameters and coefficients of friction, for both small- and large-amplitude actuations. In Section 4 we analyze the three-link snake model and compute optimal kinematics of the relative angles at one realistic set of coefficients of friction using a Fourier series representation. In Section 5, we summarize our work and suggest directions for future study.

2 Problem setup for two-link model

Figure 1: Two-link snake model; see text for description.

The snake is modeled in 2D as two identical inextensible line segments (links) connected via a hinge joint as depicted in Figure 1. The total length of the snake is , and for each link is . The snake shape is parametrized by the angle between the tail link (left) and the head link (right), with the positive direction of rotation being counterclockwise. Denoting the snake’s mass per unit length as , the total mass is . We denote as the arc length between any point of the snake and the tip of its tail, so . The tail tip, hinge joint and head of the snake correspond to and respectively.

Motion of the two-link snake can be observed both in an inertial frame with origin at a fixed point , or in a body-fixed frame with origin at the center of mass . The unit vector is parallel to the line connecting the links’ centers, and is rotated by 90 degrees. The position of in the inertial frame is , and the orientation is the angle from to . In the inertial frame, the position of an arbitrary material point on the snake is denoted as , and is the angle between the tangent to the snake at a given point and . The positive direction of rotation is counterclockwise. In the body-fixed frame, the position of the same material point is denoted , and the tangent angle is . For a given material point, we define the configuration variable in the inertial frame , and in the body-fixed frame . Specifically, for , we have and . The relation between the configuration variable in both frames is


where is the transformation matrix. For the tail link, i.e. , the configuration in is given by


where the subscript indicates tail. For the head link, , the configuration is


where represents head.

Figure 2: Schematic of the body-frame velocities for the two-link model. Center-of-mass velocity is , the velocity due to rotation about is , where is the position of a material point in body frame, and the velocity due to shape change consists of a horizontal motion of the link centers and rotations about the link centers with angular velocities for the two links (in this example, and ).

In the inertial frame, the linear and angular velocities of any point on the snake are given by the time derivative of its configuration, that is . In particular, the velocity of in the inertial frame is given by . The velocity of with respect to the inertial frame can also be expressed in the body frame as , where . Similarly, the velocity of any material point with respect to the inertial frame can be expressed in the body frame as


where is due to shape changes and is due to rotation of the snake about (see Figure 2). The relation between the velocity in both frames is given by


Perpendicular to the plane of motion, the forces on the snake are gravity and the supporting force from the ground, which balance each other. Within the plane of motion, the forces on the snake are external friction from the ground and internal forces. Since the coefficients of friction are anisotropic, the frictional force is decomposed into components in different directions. In the body frame, we denote the linear velocity of an arbitrary material point as , the unit vector tangent to the snake as , and the unit normal vector as . Our model for the snake mechanics is essentially the same as that in [12]. The Coulomb frictional force density at a given point on the snake is


where and are the forward, backward and transverse coefficients of friction respectively, , and is the Heaviside function, used to distinguish forward and backward friction [12]. Note that the frictional force density depends on the direction of velocity but not the magnitude. In addition to the external force, there are also forces internal to the snake, and the internal force density is denoted as . The torque density with respect to due to friction is given by , while that due to internal force is . The internal force and torque densities are due to a system of equal and opposite tension and shearing forces acting on adjacent sections of the snake across their interfaces [24], which makes the snake inextensible and enforces its shape. The integrals of internal forces and torques are zero


Now we nondimensionalize the variables. We consider actuations (and resulting snake motions) that are periodic in time with period . Variables are nondimensionalized by scaling by the total length , period and mass . Three important dimensionless numbers for the snake dynamics are


Here is the Froude number, which can be written as a ratio of snake inertia to the forward frictional force acting on it. The other two parameters are friction coefficient ratios. We assume the coefficients of friction are uniform along the snake, and define the forward direction as that with the smaller of the tangential friction coefficients (if ), so by definition. For real snakes, , which means that the snake’s inertia is negligible compared to frictional forces [12, 13].

For simplicity, we now drop the tildes with the understanding that all variables are dimensionless in the remainder of this work. The dimensionless point-wise frictional force density is given in the body frame by


It can be transformed into the inertial frame via the transformation matrix . The governing equations are given by the linear and angular balance laws,


Since the inertia term is negligible, the Froude number is assumed to be zero: . Therefore, substituting (7) into (10), the governing equations are reduced to the integrals of frictional force and torque densities equaling zero,


The above equations give three independent scalar equations for the three unknowns , and the solution depends on the parameters and , as well as the prescribed relative angle . Notice that (11) are intrinsically nonlinear, and the nonlinearity primarily lies in the form of friction given in (9). An explicit derivation of the equations of motion is given in Appendix A.

Since only depends on the direction of the velocity, and no inertia term is present in the equations, the only time scale in the problem is . For a given , if is doubled, the speed of the motion is reduced by half, but the snake will trace exactly the same trajectory, as does Purcell’s three-link swimmer in Stokes flow. This is referred to as the rate independence or time invariance of inertialess systems: if a body undergoes a deformation, the trajectory traveled by the body between two different shape configurations does not depend on the rate of deformation, but only on the sequence of deformation [2, 5, 22]. On the other hand, the Scallop theorem indicates that inverting the shape change sequence corresponds to inverting time for the original shape change sequence (note this is not equivalent to time invariance). A corollary [1, 10] is: in the body-fixed frame, the velocity of the center of mass is proportional to the velocity of the shape change. This is known as the kinematic reconstruction equation in the geometric mechanics literature


where is a shape change vector for systems with multiple degrees of freedom, is the local connection matrix, which is an matrix that relates an shape change velocity to an velocity of in the body frame . For the two-link model, is a scalar and is a vector, hence is a vector. Note that does not depend on .

However, equation (12) does not apply to our snake model. This is due to the anisotropy of coefficients of friction, specifically, the Heaviside function in the force equation (9), that causes the irreversibility of shape change. Instead, for our system, the local connection matrix also depends on the direction (or sign) of , denoted by , but not the magnitude . If there were no Heaviside function in (9), and force were only decomposed into tangential and transverse directions, then the system would become very similar to the multi-link fish problem in a potential flow, for which it is known that with the added inertia decomposed into tangential and transverse directions (12) holds [15]. Consequently, the techniques of analyzing the local connection matrix developed in [1, 10] cannot be directly implemented here since it is based on (12). In general, the modified kinematic reconstruction equation for the snake model can be written as


For the two-link model where is a scalar, the above equation is reduced to . Note when relative angle is prescribed, (13) is a nonlinear algebraic equation rather than a differential equation, and the nonlinearity lies in the form of . The solution of (11) or equivalently (13) is the velocity of expressed in the body frame, . Without loss of generality, assume the snake starts from the origin with zero initial orientation angle, i.e. . The configuration in the inertial frame is then


which is an iterative integral in time since is only known from the integral . The distance traveled by during one period observed in the inertial frame is given by


The work generated by the snake during one period is equal to the energy dissipation due to friction since the system is inertialess, i.e.


here it is more convenient to express in the body frame. The efficiency of locomotion is defined as the ratio between distance and work,


This efficiency , after nondimensionalization, is equivalent to the inverse of the cost of locomotion commonly seen in the animal locomotion literature [5]. Intuitively speaking, is analogous to the concept of “miles-per-gallon”.

(a) vs
(b) vs
Figure 3: (a) Relative angle prescribed as a triangular wave with period and amplitude . (b) Angular velocity for relative angle .

3 Two-link model analysis

We first look at an example of two-link snake locomotion. We prescribe the relative angle as a triangular wave with period and amplitude , and with , i.e.


Figure 3 depicts and within one period. For concreteness, let and , which are taken from the experimental measurement in [13]. The equations of motion (11) are solved numerically using the subroutine fsolve in MATLAB, which implements a Trust-Region-Dogleg algorithm. When , the trajectory of the center of mass is shown in Figure 4(a) with five snapshots of the snake at and 1 overlaid (the head of the snake is represented by ), and the orientation as a function of time is depicted in Figure 4(b). One can see from the trajectory plot that the displacement in the direction is larger than that in the direction. And even for the large amplitude of actuation , the rotation is very small: degrees for all time. The distance traveled during the period is , the work is , and the efficiency is . Figure 5 shows the velocity of the center of mass in the inertial frame. The horizontal velocity is always non-negative. The vertical velocity is non-negative during the first half-period, and non-positive during the second half. The two half periods nearly cancel out and result in nearly zero net displacement in the direction as shown in Figure 4(a). The angular velocity alternates between positive and negative during the period, and is symmetric about . All components of velocity become 0 when and , which corresponds to . From Figures 4 and 5, one can see the resulting velocities depend nonlinearly on . They can only be solved numerically. The nonlinearity arises from both the form of the friction in (9), and the large amplitude of actuation, in this case . In order to focus on understanding the former nonlinearity, we now analyze an actuation with small amplitude.

(a) trajectory of and snapshots
(b) orientation vs.
Figure 4: Motion of two-link snake when is prescribed as a triangular wave given in (18) and : (a) Trajectory of center of mass in the inertial frame with snapshots of the snake at and 1; (b) Orientation as a function of time. The coefficients of friction are and .
(a) vs.
(b) vs.
(c) vs.
Figure 5: Inertial frame velocities of the center of mass for the two-link model with given by (18) and . The coefficients of friction are and .

Small amplitude analysis

For a general kinematics , assume , and consequently . One can show that as follows. Specifically, when the motion of the snake is as depicted in Figure 2, that is and , we show in Appendix A that the integral form of the equations of motion (11) results in


The solutions for equations (19) are given by


where is the solution of the following transcendental equation


Notice depends on and but not on . For , can be approximated by


which is between 0 and -1 for . For details of the derivation of (19) and (20), as well as some numerical results regarding (22), see Appendix A. From (20), one has and since , which is consistent with the observation mentioned before: . Since , one has for a period of time . Hence, the transformation matrix given in (1) is approximately an identity matrix, which means the inertial frame velocity can be approximated by the body frame velocity, i.e. . One can easily obtain the velocities for other combinations of the signs of and by symmetry.

(a) vs.
(b) vs.
(c) vs.
Figure 6: Inertial frame velocities of the center of mass for the two-link model with given by (18) and small amplitude of actuation . The coefficients of friction are and .

As an example, when is given by (18) and , the velocity of is numerically solved from the full nonlinear equations (11), and the result is plotted in Figure 6. The large-amplitude solutions plots (Figure 5) approximately contain those in Figure 6. Specifically, for nearly coincides with the two subintervals where is between for the case with , when time is dilated by a constant factor. Hence, the results in Figure 6 also nearly coincide with the portions around zero in Figure 5. One can observe that the order of magnitude of is much smaller than . For the first quarter period, , and , which is consistent with the analytical solutions in (20). When is close to and , the velocities can no longer be approximated by linear or quadratic functions. This is because can no longer be approximated by a constant since is discontinuous, but the velocities are still bounded. Due to the symmetry in in the four quarter periods, the distance traveled by and the total work during one period are given by


since . Therefore, the efficiency for small-amplitude actuation is approximately


This shows that the efficiency is maximized when and (i.e. . Note that the linearization is only valid when . When is comparable to , which means the body is almost frictionless in the transverse direction, the terms in which are of higher order in cannot be ignored, and these higher-order terms depend on as well. This means that transverse friction provides the leading-order contribution to the work, and tangential (forward and backward) friction is comparable at higher order (see Appendix A for details). However, (24) is valid for with fixed , which indicates , or small amplitude actuation is energetically inefficient, and the efficiency increases linearly with . Thus, maximum efficiency occurs at large amplitude, where geometric nonlinearities play a role.

Large amplitude optimization

We now investigate the case when the amplitude of actuation is not small in general. Consider a periodic actuation of the relative angle that varies in the interval during . Without loss of generality, assume , and . To simplify the analysis, we consider the case that and are the only extrema of during the period. In other words, varies monotonically between the maximum and minimum, i.e. when , and when . In general, if more local extrema exist, one can always divide the period into multiple parts at the extrema, and the following analysis can be modified accordingly. Recall that our system has the modified kinematic reconstruction equation (13). That is, for the two-link problem,


where the exact forms of the components and can be derived from the equations of motion (11). We will show mathematically that the trajectory of only depends on the path of but not on the speed along the path. For a prescribed , one has


where and are integration variables. Once is obtained, the position of can be given by  (14),


where the integration can be evaluated similarly to (26). Therefore, only depends on the path of via the integration limits, and the speed does not explicitly appear in the expression (although its sign does appear). For distance , work and consequently efficiency , the integration is over the whole period from to 1, and hence they only depend on the extrema of during the period: and . This is because the parameter space of shapes is only one dimensional, so the path of is defined by the endpoints and .

Due to the nonlinearity of the system, the calculation of via (26) and (27) is not trivial. However, from the numerical results of and , one can observe that the orientation of the snake is usually very small during the locomotion. This is due to the symmetry in shape about the axis, and the fact that the asymmetry has little effect on rotation. Indeed, even for a large amplitude , one still has degrees for all time. Hence, one can closely approximate the problem by assuming during the period. Therefore, from (4), the velocities of in both frames are approximately the same, i.e. , in which the non-zero components are the linear velocity,


Since the Heaviside function appears in the tangential but not the transverse direction, comparing the force equation in (11) for the same and opposite gives


This is also evident from the results shown in Figure 5(a) and (b): for two instants symmetric about 0.5, which have equal but opposite , the corresponding are nearly equal but are nearly opposite. To simplify the notation, denote , and . They are nonlinear functions of only. From (27) and (28), integrating the velocity in the inertial frame for the whole period yields


where depends on the form of and the integration limit . The distance is given by


since . The power, or rate of work, done by the snake at a given time can be written


Power is identical for the instants with the same value of and same magnitude but opposite sign of . That is, . Therefore, one has . Similarly to our definition of , we define . The total work is given by integrating power over the period,


where . The integrals and only explicitly depend on their upper integration limits (and implicitly on the parameters and ). By definition, they are odd functions, for instance, . For symmetric cases, (e.g. the one shown in Figure 3), and the efficiency is simplified to

(a) vs.
(b) vs.
(c) vs.
Figure 7: Two-link model (a) , (b) and (c) as functions of maximum amplitude of actuation for and .

Figure 7 shows and as functions of the maximum amplitude for parameters and . One can see that is maximized around , but since increases faster at larger , the efficiency is maximized around . We emphasize that, since the parameter space of shapes is only one dimensional, as long as is the same, the efficiency and trajectory traveled by are the same. For example, any function that has the same maximum and minimum (and varies monotonically in between) as the one shown in Figure 3, for example a cosine function with period 1 and amplitude , will also result in the trajectory shown in Figure 5(a) and have the same .

(a) for and various
(b) for and various
Figure 8: Two-link efficiency as a function of actuation amplitude for various and : (a) fixed and various and (b) fixed and various .

We repeat the process above to calculate as a function of for various and , and the results are shown in Figure 8 for (a) and various and (b) and various . In (a), the efficiency-maximizing for , and increases as drops below 1, in which case there is a smaller increase in work for the additional transverse motion associated with a larger amplitude. In (b), one can see that for , the efficiency-maximizing regardless of the value of . Note that when is small, varies almost linearly with , which agrees with the small-amplitude result in (24).

(a) contour
(b) contour
(c) contour
Figure 9: Contour maps of (a) distance , (b) work and (c) efficiency as functions of and (on a log-log scale), for .

We now determine the values of and that maximize . Figure 9(c) shows a contour plot of as a function of and when . Here is evaluated at points on a logarithmic grid with nodes spaced by factors of 1.5, and ranges from to and from to . The highest efficiency occurs at the largest and the smallest in this range. In Figure 9(a), the distance increases when and increase. By contrast, work is mostly independent of and increases when increases, as shown in (b). Note this is consistent with the small amplitude result in (23). The trends shown in Figure 9 hold for a wider range of parameters, which suggests that the efficiency-maximizing parameters are large and small.

4 Three-link model

Figure 10: Three-link snake model; see text for description.

We now consider a snake model with three links (each with length ) connected by hinge joints as depicted in Figure 10. The shape is now given by two relative angles, and . Since the links cannot penetrate each other, the angles are constrained to lie in the set

Figure 11: Shape change parameter plane . The shaded areas are infeasible due to the mutual avoidance of the links. A periodic kinematics is a directional closed path , with the initial state marked by .

In Figure 11 the infeasible regions are shaded at the upper right and lower left corners in the plane. One can describe the motion in the inertial frame or in the body-fixed frame , with the angle between and now given by


Here , and are the orientations of the tail, middle and head links in the inertial frame, respectively. In the body frame, the configuration variable on each of the three links is