Optimizing snake locomotion on an inclined plane

Optimizing snake locomotion on an inclined plane


We develop a model to study the locomotion of snakes on an inclined plane. We determine numerically which snake motions are optimal for two retrograde traveling-wave body shapes—triangular and sinusoidal waves—across a wide range of frictional parameters and incline angles. In the regime of large transverse friction coefficient, we find power-law scalings for the optimal wave amplitudes and corresponding costs of locomotion. We give an asymptotic analysis to show that the optimal snake motions are traveling-wave motions with amplitudes given by the same scaling laws found in the numerics.

87.19.ru, 87.19.rs, 87.10.Ca, 45.40.Ln

I Introduction

Snake locomotion has recently drawn interest from biologists, engineers and applied mathematicians Astley and Jayne (2007); Hu et al. (2009); Maneewarn and Maneechai (2009). A lack of limbs distinguishes snake kinematics from other common forms of animal locomotions such as swimming, flying, and walking. Snakes propel themselves by a variety of gaits including slithering, sidewinding, concertina motion, and rectilinear progression Hu et al. (2009). They can move in terrestrial Gray and Lissmann (1950); Jayne (1986), aquatic Shine et al. (2003), and aerial Socha (2002) environments. Snake-like robots have shown the potential to realize similar locomotor ability, and have potential applications in confined environments like narrow crevices Erkmen et al. (2002), as well as rough terrain Yim et al. (2003). In such environments the ability to ascend an incline is fundamental, and various studies of snakes and snake-like robots have been carried out on this subject. Maneewarn and Maneechai Maneewarn and Maneechai (2009) examined the pattern of crawling gaits in narrow inclined pipes with jointed modular snake robots and found that high speeds were obtained with short-wavelength motions. Hatton and Choset Hatton and Choset (2010) focused on sidewinding gaits on inclines and presented stability conditions for snakes by comparing sidewinding to a rolling elliptical trajectory and determining the minimum aspect ratio of the sidewinding pattern to maintain stability. Marvi and Hu Marvi and Hu (2012) studied the concertina locomotion of snakes climbing on steep slopes and in vertical crevices. They found that snakes can actively orient their scales and lift portions of their bodies to vary their frictional interactions with the surroundings. Transeth et al. Transeth et al. (2008) considered the obstacle-aided locomotion of snake robots on inclined and vertical planes. Their model included both frictional forces and forces from rigid-body contacts with the obstacles, and their numerical results were consistent with their robotic studies. In this work, we focus on the slithering motion of snakes on an inclined plane by extending a recently-proposed model for motions on a horizontal plane Hu et al. (2009); Hu and Shelley (2012); Jing and Alben (2013); Alben (2013) to motions on an incline. Here the snake is a slender body whose curvature is prescribed as a function of arc length and time. For simplicity, we do not consider elasticity or viscoelasticity in the snake body. The external forces acting on the snake body are Coulomb friction with the ground, and gravity. The model has shown good agreement with biological snakes on a horizontal plane Hu et al. (2009); Hu and Shelley (2012), and previous studies have used the model to find optimally efficient snake motions. Hu and Shelley Hu and Shelley (2012) prescribed a sinusoidal traveling-wave body motion and found the optimally efficient amplitude and wavelength of the traveling wave. Jing and Alben Jing and Alben (2013) used the same model to consider the locomotion of two- and three-link snakes. They found optimal motions analytically and numerically in terms of the temporal function for the internal angles between the links. Alben Alben (2013) considered more general snake shapes and motions by optimizing the curvature as a function of arc length and time with 45 (and 180) parameters, across the space of frictional coefficients. He found that the optimal motions are two types of traveling-wave motions, retrograde and direct waves, for large and small transverse friction coefficients, respectively. In the large transverse friction coefficient regime, he showed analytically that the optimal motion is a traveling wave, and found the scaling laws for the wave amplitude and cost of locomotion with respect to the friction coefficients, both numerically and analytically.

In this paper, we confine our discussion to the regime where the transverse friction coefficient is larger than the tangential friction coefficients. This is the typical regime for biological and robotic snakes Hu and Shelley (2012); Hirose (1993). We prescribe the snake’s motion as a retrograde traveling wave with two shape profiles—triangular and sinusoidal—but with undetermined amplitudes. The triangular wave motion has analytical solutions and embodies many aspects of general traveling-wave motions. We obtain the optimal motions in terms of the amplitudes of the triangular wave in the three-parameter space of transverse frictional coefficient, tangential (forward) frictional coefficient, and the incline angle. We discuss the relative effects of these three parameters and find the upper bound on the incline angle for the snake to maintain an upward motion. We also find the power law scalings for the optimal amplitudes and corresponding costs of locomotion with respect to the three parameters. For the sinusoidal wave motion, we use a numerical method to solve for the position of the snake body from its prescribed curvature. We then obtain the optimal body shape numerically and show that it follows the same scaling laws as the triangular wave motion, providing confirmation of those results. In the last part of the paper, we analytically determine the asymptotic optima for more general snake motions in the regime of large transverse friction coefficient. We obtain the same scalings analytically as were found for the triangular and sinusoidal waves, which confirms and generalizes those results. Our study of snake locomotion can also be generalized to other locomotor systems as long as the same frictional law applies. One example is the undulatory swimming of sandfish lizards in sand Hatton et al. (2013); Ding et al. (2012).

The paper is organized as follows: Section II describes the mathematical model for snake locomotion on an inclined plane. Section III discusses the optima of the triangular wave motion, while Section IV studies the sinusoidal wave motion. The analytical asymptotic discussion is presented in Section V, and the conclusions follow in Section VI.

Ii Model

We use the same frictional snake model as Hu et al. (2009); Hu and Shelley (2012); Jing and Alben (2013); Alben (2013) to describe snake locomotion on an incline. The position of the snake is given by , a function of arc length and time . The unit vectors tangent and normal to the snake body are and respectively. The snake is placed on a plane inclined at angle with respect to the horizontal plane. The - axes are oriented with the axis extending from origin directly up the incline and the axis rotated 90 degrees from it and directed across the incline. Height is constant along the axis. A schematic diagram is shown in figure 1.

Figure 1: Schematic of the snake position on a plane inclined at angle . The arc length is nondimensionalized by the snake length. The tangent and normal vectors are labeled at a point. Forward, backward and transverse velocities are shown with corresponding friction coefficients , and .

The tangent angle and the curvature are denoted and . Given the curvature, the tangent angle and the position of the snake body can be obtained by integrating:


The trailing-edge position and the tangent angle are determined by the force and torque balances for the snake:


Here is the mass per unit length and is the length of the snake. We assume the snake is locally inextensible, and and are constant in time. is the external force per unit length acting on the snake. It includes two parts: the force due to Coulomb friction with the ground Hu et al. (2009), and gravity:


is the Heaviside function, and represents the component of gravity in the downhill () direction. The hats denote normalized vectors and we define to be when the snake velocity is . The friction coefficients are and for motions in the forward , backward and transverse directions. Without loss of generality we take , so the forward direction has the smaller friction if it is unequal in the forward and backward directions.

We prescribe the curvature as a time-periodic function with period and nondimensionalize equations (4)-(6) by length , time and mass . We then obtain:


We neglect the snake’s inertia for simplicity, as for typical steady snake locomotion Hu et al. (2009) and set the left hand sides of (8)-(10) to zero. As discussed by Alben Alben (2013), the model still maintains a good representation of real snake motions. We then obtain the following dimensionless force and torque equations:


and the dimensionless force becomes:


The frictional force tends to 0 as approaches . On a strictly vertically plane, frictional force is unable to balance gravity, so planar locomotion is not obtained by our model in this case (however, snakes can ascend vertical crevices in a non-planar concertina motion Marvi and Hu (2012)).

Given the curvature , we solve the three nonlinear equations (11)-(13) at each time for the three unknowns and . Then we obtain the snake’s position as a function of time by equations (1)-(3). We define the cost of locomotion as


where is the distance traveled by the snake’s center of mass over one period


and is the work done by the snake against frictional forces and gravity over one period


Our objective is to find the curvature that minimizes . We choose the initial orientation of the snake so that its center of mass travels only in the direction (up the incline) over one period of motion.

We briefly mention the case in which the snake moves down the incline, which is equivalent to setting . In this case the snake can slide down the incline with no change of shape, and the work done by gravity and friction are equal in magnitude and opposite in sign. Thus the cost of locomotion is 0 regardless of the frictional parameters. A straight body with experiences purely tangential friction, and achieves the fastest speed among possible body shapes.

By contrast, when the snake ascends the incline, i.e. , we will consider traveling wave motions in which net tangential friction and gravity both act in the direction and transverse friction is necessary to balance the -force equation. In the following discussions we only consider , and look for nontrivial to minimize the cost of locomotion.

Iii Triangular Wave Body Shape

iii.1 Range of

We start our discussion with a triangular wave body motion. It was studied on a level plane in Alben (2013), and the snake’s position, angle, velocity and work can all be obtained analytically. The motion is useful to consider because it embodies many aspects of more general traveling wave motions, while the shape dynamics are easy to understand.

The triangular wave has zero mean deflection:


with the unit tangent and normal vectors:


The force and torque balance equations are satisfied when the snake moves forward with a constant speed . Then the horizontal and vertical speeds are


The net -force and torque for such a motion are identically zero. We determine the horizontal speed by the -force balance equation:


Since , and the tangential velocity is uniformly forward or backward over the whole snake body for the triangular wave, the most efficient motion is obtained when the snake moves forward. Thus does not appear in the frictional force in (21). Notice that the tangential frictional force and gravity both have a component in the direction, and the transverse frictional force is the only source for a balancing force in the direction, up the incline. Solving (21) for we obtain:


The speed of the snake is a function of , , and . We require the speed to be real and nonnegative, and therefore the incline angle must satisfy the following inequality:


Here we use the fact that the amplitude in the triangular wave model and . In figure 2, we plot the upper bound of according to inequality (23).

Figure 2: Lines showing the boundaries of the region on nonnegative forward velocity in the space of and , for different and . On the lines, = 0, and below the lines, . The black dashed line shows , giving a vertical incline.

For a given value of , nonnegative speed is obtained for in the region bounded by a curved line (labeled by ) and the horizontal () axis. As increases, larger transverse friction can be generated for the same , and therefore the range of increases accordingly. As becomes larger, the tangential motion produces a stronger downhill drag which inhibits upward motion, so the corresponding range decreases. When the amplitude varies from 0 to 1, both the transverse and tangential frictional forces vary and their -components have opposite sign. As a result, the range of is non-monotonic with respect to . The largest upper bound is obtained at .

iii.2 Optima and other results

In the triangular wave motion, the velocity and power are both constant over time, so we can simplify the cost of locomotion as


and obtain


We plug the value of (22) into (25) and minimize with respect to . Then in the large- limit we obtain the optimal cost of locomotion and corresponding amplitude as:


We plot the optimal motion with , and over one period in figure 3. We manually offset the body by a constant increment in the direction with every snapshot to clearly show the individual bodies, but we note that for the triangular wave motion, the snake’s center of mass moves purely along the -axis. The peak of the snake shifts to the left in the figure, which indicates the snake moves slower than the traveling triangular wave. The snake slips transversely in the direction to obtain a thrust force in the direction that balances the drag forces due to tangential friction and gravity.

Figure 3: The optimal snake trajectory of the triangular wave body shape over one period obtained with , and .

In figure 4 we plot the optimal and with respect to , , and . Our results extend to below the large- limit, but in this limit the results agree with (26) and (27). For each parameter set, we minimize over using equations (22) and (25).

Figure 4: Scaling laws of the triangle wave optima. (a) vs. with various and fixed . The solid line indicates the scaling . (b) vs. with various and fixed . The solid line shows the scaling . (c) vs. with various and fixed . The solid line is the asymptotic solution , obtained with . (d) vs. with various and fixed . The solid line shows the asymptotic optima with .

We plot versus in figure 4a with fixed and vary the parameter . The asymptotic scaling of is shown with the solid line. The corresponding and the scaling (solid line) are plotted in figure 4b. The optimal magnitude and cost of locomotion both decrease with larger . We vary and in figure 4c and d with fixed , and plot the optimal versus and versus respectively. The analytical solutions of (26) and (27) for are shown with solid lines in both figures, and they agree well with the numerical results at the largest (downward-pointing triangles). The optimal amplitude achieves its minimum at , while the cost of locomotion monotonically increases with as . When goes up, the optimal increases accordingly and its minimum over shifts to larger (figure 4c). But the cost of locomotion varies non-monotonically with . We can rewrite (27) as


The least efficient optimum is obtained when


We call the critical incline angle. The optimal snake moves more efficiently when the incline is either shallower or steeper than the incline at the critical angle. The critical incline angle only depends on the tangential friction coefficient and becomes smaller as increases. On a steeper slope, more work is done against gravity and less work against forward friction, for a given distance travelled. Thus when increases, efficiency can be improved by making the slope steeper (and adjusting the amplitude to achieve the optimum at the new set of parameters).

To better understand the effects of the parameters , and , we plot the costs of locomotion due to transverse friction alone and tangential friction alone versus in figure 5. We decompose the cost of locomotion (25) into three parts:




are the costs due to transverse friction, forward tangential friction, and gravity, respectively.

In figure 5, we vary one of the parameters , and in turn and keep the other two fixed. We use solid lines for transverse friction and dashed lines for tangential friction in all panels. In general, as the amplitude becomes larger, the cost due to transverse friction decreases while the cost due to tangential friction increases. The optimal is obtained when the slopes of the two costs are equal in magnitude and opposite in sign. is independent of so it does not play a role in determining the optimal amplitude with respect to .

Figure 5: Cost of locomotion versus with various , and . (a) vs. with various and fixed and . (b) vs. with various and fixed and . (c) vs. with various and fixed and . The solid lines show the cost due to transverse friction and the dashed lines denote the cost due to tangential friction.

In figure 5a, the sums of the costs of the frictional forces become smaller as increases. Thus the optimal decreases as well as shown in figure 4b. For a given motion (a given A), the slope of the tangential cost is almost unchanged as goes up, while the magnitude of the slope of the transverse cost quickly decays. The point where the two slopes balance shifts to the left at larger . The optimal motion is thus obtained at a smaller amplitude as increases. We show the results only for and in the figure panel, but the same phenomenon holds for all and . Physically, as the transverse coefficient increases, the snake can obtain the same amount of forward force from transverse friction with less deflection of the body and less slipping in the transverse direction, and the cost of the tangential friction is reduced as well due to a shorter path travelled. Thus, the total cost of locomotion decreases with .

We show in figure 5b that the costs due to transverse friction and tangential friction both increase as increases. When the friction coefficient is larger, the snake of the same deflection experiences a stronger downward drag caused by tangential friction, increasing the slipping and consequently the work done against transverse friction as well. The optimal amplitude varies non-monotonically with as shown in figure 4c. The slopes of both costs increase with for given . When , the point where the two slopes are equal in magnitude shifts to smaller as increases. When , the balanced point shifts to larger .

In figure 5c, we fix and , and vary the incline angle . The cost due to gravity is for the triangular wave motion and thus always increases with . Meanwhile, the tangential cost decreases with while the transverse cost increases. The competition of these three costs makes non-monotonic with as shown in figure 4d. For a given motion, the slope of the tangential cost with respect to is nearly unchanged as the incline becomes steeper. However, the magnitude of the slope of the transverse cost increases with , so a point with a given slope of the transverse cost shifts to larger as increases. Therefore, the point where two slopes are equal and opposite shifts to the right as goes up, and the optimal in figure 4a and c grows with .

Iv Sinusoidal Body Shape

iv.1 Numerical Method

We now consider an alternative, sinusoidal snake motion, to check the dependence of our results on the snake shape. We again determine the snake shape which minimizes for a given parameter set . We consider the sinusoidal body shape


where the curvature is prescribed as a sinusoidal function with -period . We fix the wave number in this work and look for the optimal constant to minimize for a given . Later we show that although the wave number affects the optimal value of , it does not change the dependence of the optimal on the other parameters. We fix in this section, since for a horizontal plane, the lowest cost of locomotion is obtained in the limit of large wave numbers according to Alben (2013), and approximates this limit while only requiring a moderate number of grid points in arc length along the snake to discretize the equations accurately. We find the optimal on a sequence of one-dimensional meshes of values with decreasing spacing. Each mesh in the sequence is centered near the minimizer from the previous coarser mesh. The fourth mesh used has a mesh size of . We then use a quadratic curve to fit the data around the minimum point on the fourth mesh and obtain the optimal based on a final, fifth mesh, with mesh width . In Alben (2013) Alben used a BFGS algorithm to minimize the cost of locomotion when the curvature is described by a double-series expansion with 45 parameters. In this work we optimize over only a single parameter (the amplitude), for a given shape. We find that a direct search over the parameter space is typically fast enough since varies smoothly with in the regime of physically admissable, non-overlapping shapes.

The algorithm requires fast routines to evaluate . Here we describe our numerical scheme that solves for the work, distance and cost of locomotion. Given the curvature , we solve the three nonlinear equations (11)-(13) at each time step, over a period, for the three unknowns and . Then we use equations (1)-(3) to compute and obtain and over one period.

We discretize the period interval uniformly with time points: . At each time, we rewrite equations (11)-(13) as equations in unknowns by taking time derivatives on both sides. We solve for over one period and then integrate to obtain . The advantage of replacing and with their time derivatives as variables is that it can reduce the numerical error in computing the discrete time derivatives and decrease the computational complexity by decoupling the large system of equations in unknowns into decoupled systems each containing only equations in unknowns Alben (2013).

We design a time-marching scheme which is second-order in both time and space to solve for at each time level. At time , we set to and solve for using Newton’s method with a finite-difference Jacobian matrix as described in Alben (2013). At each time level , we need to carry out the computation. We obtain an initial guess for the current position and angle by a forward Euler scheme using the previous step solutions and . We then solve for at time . Then we correct by integrating from to . Our method is similar to the prediction-correction algorithm for solving a system of ODEs. We can iterate the same procedure until a certain accuracy is obtained. We find that second-order temporal accuracy is achieved by only performing the above procedure once.

iv.2 Optima and other results

We first consider the parameter regime where . In this region, the tangential motion is purely in the forward direction for the sinusoidal wave motion. Therefore, as for the triangular wave, the term drops out of the force law, and the parameter space is reduced to . We plot the optimal snake trajectory with parameters , and over one period in figure 6. The snake moves from left to right and its center of mass moves mainly along the direction.

Figure 6: The optimal snake trajectory of the sinusoidal wave body shape over one period obtained with , and .
Figure 7: Scaling laws for the sinusoidal wave optima. (a) vs. with various and fixed . The solid line shows the scaling law . (b) vs. with various and fixed . The solid line denotes the scaling . (c) vs. with various and fixed . The solid line is the analytical optimum obtained with (d) vs. with various and fixed . The solid line is the graph of the optimal with

In figure 7, we vary , and , and plot the optimal and cost of locomotion versus these parameters respectively. Some data points are ignored because there is no solution with non-negative -velocity for the corresponding parameter values. We fix and plot the optimal and versus with various in figure 7a and b. In figure 7c and d, we fix and vary from to with different . We find that the optima for the sinusoidal wave motion satisfy essentially the same scaling laws as the triangular wave motion in (26) and (27). The cost of locomotion is the same as (26) and the amplitude is scaled by an extra factor yielding a deflection amplitude that agrees with (27).

Figure 8: Optimal snake shapes at the instant for (a) various from to and fixed and ; (b) various from to and fixed and ; and (c) various from to and fixed and . (d) The deflection of the body rescaled by corresponds to (a). (e) The deflection rescaled by for (b). (f) The deflection rescaled by with the same parameters as (c).

In figure 8, we show the optimal snake body shapes corresponding to different , and in panels a, b and c respectively. We displace the different bodies vertically so that they are easier to distinguish. In figure 8d, e, and f, we rescale the deflections of the optimal shapes with , , and according to the scaling law for the optimal amplitude (27). The centers of mass for all bodies are located at the origin. Here we zoom in on the portion of the body nearest to the center. We find a good collapse of the bodies after rescaling.

In the regime where and are comparable, the sinusoidal body shape model may not yield forward motion. For snake locomotion on a level plane, Alben Alben (2013) found that when , the optimal snake shape is no longer a retrograde traveling wave. He identified a set of locally optimal motions and classified some of them as racheting motions. When the snake climbs uphill and is small, the traveling wave may not be able to provide enough uphill thrust to balance gravity and the drag due to forward motion. The snake may therefore use locomotion modes other than slithering to maintain its position on the incline. For example, concertina locomotion is often observed for snakes moving inside an inclined tunnel. Marvi and Hu Marvi and Hu (2012) found that some snakes can resist sliding by lifting part of the body and reducing the contact with the ground to several localized regions. This shows that a planar model is not sufficient to describe snake locomotion on inclines when the ratio of the transverse-to-tangential friction coefficient is small. Our future work may consider three-dimensional motions to better understand this interesting parameter regime.

V Asymptotic Analysis

v.1 Optimal Shape Dynamics

We now analytically determine how the optimal snake motion depends on the parameters , thus providing theoretical confirmation of our previous results and extending them to general shapes, in the large- regime. We assume that the mean direction of the motion is aligned with the -axis and the deflections along the direction are small, i.e. and all of its temporal and spatial derivatives are for some negative . We also assume that the tangential motion is only in the forward direction, to simplify the derivation. We first expand each term in the force and torque balance equations in powers of and retain only the terms which are dominant at large . A detailed discussion of the expansion of each term can be found in Alben (2013). If we only keep the lowest powers in from each expression, the -force balance equation becomes:


where is the -averaged horizontal velocity. The three terms in the integral represent the drag due to forward friction, the thrust due to transverse friction, and gravity, respectively. We note that both tangential friction and gravity forces act in the direction, and transverse friction is essential to maintain the -force balance. In Alben (2013) it is shown that for large , a minimizer of the cost of locomotion should approximate a traveling wave motion. Therefore we pose the shape dynamics as


which is a traveling wave with a prescribed wave speed . is different from in general, for otherwise the snake moves purely tangentially with no transverse motion. Here is a periodic function with period . We obtain an equation for in terms of and by substituting (36) into (35):


where . As , and no frictional force occurs. Therefore, in the limit , we require such that , i.e, the speed at which tends to depends on the speed at which tends to infinity. This requirement is similar to the upper bound of in the triangular wave motion to obtain forward motion. As , equation (37) holds with . The traveling wave moves backward along the snake at speed while the snake moves forward at a speed slightly less than . Therefore, the snake slips transversely to itself, which provides an uphill thrust to balance gravity and the drag due to tangential friction.

In Alben (2013) it is shown that one must expand the terms in the force balance equation in higher powers of to obtain an optimal motion. We do so, again assuming . Then the force balance equation becomes:


and the cost of locomotion is:


Equation (39) is shown in a simplified form, using the result that as . We substitute equation (38) into (39) and obtain:


If we approximate as constant in time, we obtain


which is minimized for


and the corresponding optimal cost of locomotion is :


In the triangular wave motion, we have


By using (42), we obtain that


which is consistent with the analytical result we obtained in (27).

v.2 Optimal Curvature Analysis

For a more general body shape, to satisfy the -force balance and torque balance, we instead prescribe the curvature and obtain , , and from equations (1)-(3). We now correct the above analysis to satisfy all three equations. We again assume that the deflection from -axis is small and we decompose as:


Prescribing the curvature is the equivalent to prescribing . We set


so is similar to the form given before, with two additional terms: vertical translation and rotation . We determine and by expanding the -force (12) and torque balance (13) equations to leading order in and obtain:


We solve (52) and (53) for and :


We then express the -force balance equation in terms of and and obtain:


Substituting (54) and (55) into (58) we solve for in terms of :


The cost of locomotion then becomes


Following Alben (2013), we expand in a basis of Legendre polynomials for any fixed time . The Legendre polynomials are orthonormal functions with unit weight on . They satisfy the relations


We write as: