The locomotion of long slender crawlers arises from the interaction of muscles and ligaments and external forces such as Coulomb friction and viscous drag from fluids. The mechanism of movement for long thin organisms has been the subject of many previous investigations (e.g., Refs. [1][2][3]). Long slender organisms across length scales and various media have adopted a variety of locomotion gaits. Terrestrial snake locomotion, for example, can be classified into one of four main locomotion modes: lateral undulation, concertina, sidewinding, and rectilinear [1]. It is known that snakes switch between these modes subject to their environment as snakes are motile in diverse media such as granular, aquatic, and aerial domains [4][5][6]. From a robotics perspective, the hyper-redundant body morphology of snakes, which provides mastery of such varied terrain and stability, makes them excellent subjects for biologically inspired designs of snake robots, surgical devices, and targeted drug delivery mechanisms.
Many animals use lateral undulation to move in all sorts of environments. Examples include crawling snakes in grass, eels in water, and sandfish in sedimentary environments such as sand [7]. Furthermore, some animals combine lateral undulation with other motility techniques, such as fracture and crack propagation [8], burrowing by sediment removal, and forms of two-anchor crawling utilizing peristaltic waves [9]. In addition to the utility of lateral undulation in various environments, it also occurs across a wide range of length scales, from reticulated pythons all the way down to nematode worms, such as beet eelworms, which propagate waves of curvature from head to tail down their body to move in sand and soft mud [10]. Since lateral undulation relies on a sliding interaction with the ground, friction is both essential and destructive to the crawler. Without friction, there can be no net forces driving locomotion; yet the friction results in energy dissipation at the crawler's skin and affects their cost of transport. Yet, the transport cost of snakelike animals moving by sliding across the ground is comparable to that of birds and limbed creatures. Because of the applicability of lateral undulation in various environments, across length scales and for various body morphologies, it forms the focus of this study [11,12].
Lateral undulation has also been employed in a new class of snake robots that are soft and compliant in nature [13][14][15]. These soft robots hold out the promise of highly customized robots whose body motion and dynamics can be programed into the material of the body. Recent soft robotic snakes can mimic their animal prototypes and use curvature waves for locomotion [16]. These robots can be used for complex search operations, inspection, and maintenance in congested and chaotic environments.
Different techniques have been proposed to model the kinematics and dynamics of crawling robots [17]. These models represent the interface between the crawler and its environment as either no-slip or frictional slip conditions. For the no-slip condition, the ground interaction is modeled as a nonholonomic constraint at every point on the body's surface. This constraint treats a snake belly as a flexible ice skate. In the slipping condition, on the other hand, the environmental forces are represented with specific frictional or resistive models, and the body is allowed to move laterally relative to the environment. Such a model has been frequently used for studying the slithering snake motion, notably in Refs. [18,19].
Slithering snakes distribute their weight away from the points of maximum curvature along their bodies in a locomotion strategy known as sinus lifting. Hu et al. [18] advanced the understanding of slithering locomotion by introducing a model explaining the critical role of weight distribution in snake locomotion. This sinus-lifting action has also been leveraged to improve the functionality of snake robots without wheels [20,21].
Researchers have used frictional models to study the optimal locomotion problem, studying the interaction of a crawler's body with a frictional environment [18,22,23]. Alben studied a planar crawler subject to sliding friction in the forward, transverse, and backward directions and used a quasi-Newton method to find the time and space curvature trajectory that optimizes locomotion efficiency [22]. Among other things, Alben discovered that the optimal direction of locomotion changes as a function of environmental friction. Later, it was shown that the optimal gait is mainly a function of the transverse and forward friction coefficient ratio, and optimal crawling performance is in the direction of the curvature wave propagation when transverse friction is low and opposite that when the transverse friction is large [23]. Alben referred to these as "direct" and "retrograde" curvature waves. Guo and Mahadevan [19] investigated the coupled effect of internal and external parameters, including the viscous and elastic moment of a viscoelastic crawler with the nonholonomic, no-slip velocity constraint. They found that parameters in the external environment mainly determine a crawler's body shape, and the speed of the crawler is primarily determined by internal parameters (exogenous and endogenous parameters, respectively). McMillen et al. found that damping plays a critical role in the swimming of lampreys [24]. McMillen et al. showed that muscle contraction in slender swimmers happens slightly after muscle activation, and this delay increases progressively from head to tail. McMillen et al. related this observation to the viscoelastic property of the tissue, the geometry of the swimmer, and the muscle dynamics, more so than the effect of fluid forces. This work indicates that damping is important in lateral undulation in general and significant for designing embodied intelligence into robotic systems.
The role of body viscoelastic damping is especially important for pneumatically actuated soft robots. Damping is inherent to these robots, which are composed of materials with very low Young's moduli and include fluidic channels for actuation [25,26]. By tuning the damping behavior of soft robot actuators, it is also possible to improve their performances [27]. However, most of the prior research focused on the importance of internal compliance within the body of a crawler [28,29]. It has been demonstrated experimentally that by leveraging the robot's flexibility, it is possible to increase the snake robot's ability to traverse two-and three-dimensional obstacles [30,31]. However, little is known about how internal damping affects the locomotion of long slender organisms and robots relative to their stiffness. The notable exception is the work done by Guo and Mahadevan, which investigated the role of internal viscous damping on a crawler with the nonholonomic constraint [19]. However, this nonholonomic constraint does not apply to the general locomotion case and particularly does not apply to soft bioinspired robots where substantial lateral slipping motion is seen [32]. In this paper, we explore the effects of internal damping and its connection with the lateral friction coefficient on the crawling performance of a modeled slender system that can serve as a guide for making decisions about the effect of damping in soft robot designs and also provide insight into locomotion. Figure 1 illustrates the concept of modeling damping effects internal to long slender organisms and robots with a nonlinear Euler-Bernoulli beam.
We study the locomotion of long slender organisms using a three-dimensional nonlinear Euler-Bernoulli beam formulation originally proposed by Tjavaras [33]. This model has been employed before to study the dynamics of animal locomotion from swimming fish to climbing lampreys [34][35][36]. In the model, an Euler-Lagrangian dual coordinate system is employed to formulate the motion, namely, the global coordinate system (x, y, z) and a local reference system with unit vectors τ, n, and b in the tangential, normal, and binormal directions of the centerline of a long cylindrical crawler, respectively. Figure 2 shows the free-body diagram of a long slender crawler, coordinate frames, external forces, and inter-FIG. 2. Free-body diagram of the local coordinate system, environmental force density, and internal torque density. The body mechanics are parametrized by arc length s. Forces and velocities are written in a local coordinate frame n, τ, and b. Muscle actuation is modeled using a traveling wave of bending moment depicted with red dashed arrows Q(s, t ). Environmental force is modeled as a linear force density F e depicted with black dashed arrows. nal torques. The two coordinate systems are related by
where i, j, and k are unit vectors in the x, y, and z directions, respectively. C is an orthogonal rotation tensor representing the change of orientation from the global to the local reference frames. C is constructed by using the singularity-free method of Euler parameters β 0 , β 1 , β 2 , and β 3 to describe the rotation based on a single angle ϑ and a principle direction as
where
Dynamic equations are derived by imposing the conservation of translational and angular momenta of an infinitesimal segment of the body. We have
and
where m is the mass per unit length. s is the distance from this point to the tail end along the unstretched crawler.
are the translational and angular velocities, respectively.
) is the axial strain, and (s, t ) = τ τ + n n + b b is the Darboux vector measuring the material torsion and curvatures of the body centerline. F e is the frictional force density, and Q is the internal moment density generated by the body, both to be defined later.
One compatibility relation is required to ensure a continuous body configuration in space and time,
The internal forces and moments are related to the strain ε and the Darboux vector through the constitutive relations
By definition, A is the cross-sectional area. EI and GJ are the bending and torsional stiffness, respectively. Aggregate elastic Young's modulus E el and viscoelastic coefficient D = αE el represents the combined contribution from all the passive elements during bending. This, in real animals, includes a viscoelastic spine, skin, white muscle, and the inactive part of red muscles. Mathematically, linear viscoelasticity is introduced by defining E = E el (1 + α ∂/∂t ). The final system of equations is
where
A detailed description of the matrices Hand P is provided by Tjavaras [33]. Overall, 13 boundary conditions are required at the two ends of the body. At the tail end, we apply seven boundary conditions: β 0 = 0; T = 0; and M = 0. Six additional boundary conditions are applied at the head: i β 2 i = 1; T = 0; and M n = 0 and M b = 0. To solve Eq. ( 7), we divide the body into N p -1 segments, each with length s, by uniform distribution of points s k (k = 1, . . . , N p ) along the unstretched length of the body. An implicit box method is then applied to integrate the system of equations from the time step t i-1 to t i = t i-1 + t [33]. The method is second-order accurate in space and first-order accurate in time. A lower-order time integration scheme is used for better numerical stability. Following the convergence study, N p = 50 is found to be sufficient for the current study.
Snakes and other crawling animals possess body features such as scales that give their bodies a higher coefficient of friction when sliding sideways relative to when they are sliding forward or backward [18,22,37]. To replicate this property in the current mathematical model, an anisotropic Coulomb friction model is used. In this model, the friction force is calculated in the body-attached Lagrangian coordinate system as a function of the body velocity V as
where Vτ and Vn are the regularized direction of velocity vectors along the body in the tangential and normal directions, respectively. They are defined as
where is a small regularized parameter chosen to ensure the friction forces change continuously with the body's velocity [23]. This parameter is chosen from sensitivity analysis; so the results do not change more than 1%.
It is known that body actuation for snakes, and limbless lizards that employ lateral undulation for locomotion, is caused by muscle actuation propagating along their bodies from head to tail [3]. This internal body actuation is modeled as a traveling wave of muscle torque given by
where L is the body length, n w is the number of waves along the body, and c is the wave propagation speed. For this study, n w is set to unity based on the previous finding that a curvature wavelength close to unity results in an optimized planar crawling locomotion [22].
In this computational study, the crawler starts from a straight-line initial configuration and from the rest condition. Internal torque produces shape changes, which, in turn, generates environmental friction interactions and net forces that accelerate the body of the crawler until the time-averaged thrust and drag forces are balanced. Once the body reaches this limit cycle gait, the phase-averaged values of different performance metrics are calculated, as will be discussed next.
Equations ( 7) and ( 8) are scaled by the body length L, the body mass per unit length m, and the wave propagation speed of the muscle torque. The following nondimensional parameters quantify the effects of viscoelastic damping, body stiffness, and muscle torque amplitude:
The frictional environmental force is quantified with two nondimensional parameters:
Another important parameter is the body curvature . The nondimensional curvature is defined as ˜ = L. From now on, to simplify the notation, we refer to ˜ as . As nominal values for the study, we select a material stiffness ẼI = 13 and torque amplitude Q = 66. These numbers are chosen to produce locomotion (for the low-damping case) with a traveling wave form with amplitude-to-wavelength ratio of ≈0.2 consistent with organisms that move in granular media as well as efficient organisms at low Reynolds number [7,38]. We note that the resulting curvature is approximately 2, which is lower than what is observed for terrestrial snakes whose same nondimensional body curvature typically reaches approximately 7 [39]. We explore the role of κ 1 while fixing the value of κ 2 = 7.8 consistent with the friction values of snake scales traveling at normal speeds based on measured values from previous studies [18,19].
As noted, various taxa, including terrestrial snakes that utilize undulation for locomotion, propagate waves of curvature from head to tail down their bodies. As in anguilliform locomotion, these body waves often increase in amplitude from the head to the tail. We quantify this increased amplitude with a dimensionless ratio d tail /d head shown in Fig. 3. Here, d head (d tail ) is the difference between the maximum and minimum lateral excursions of the head (tail).
At any given moment in time, a part of the friction density contributes to thrust, and the other part contributes to drag. To FIG. 3. In this paper, the shape ratio is the ratio of the tail-to-head lateral excursion. separate these, we defined two operators, H + and H -, as
and define thrust and drag forces acting on the swimmer along the locomotion direction as
where t v is the unit vector in the mean locomotion direction.
In addition, the muscle power output at any given time can be calculated from
It is assumed that the crawler cannot harvest power from the environment. Hence the contribution to the power expenditure through the length is always non-negative. We refer to the average of the thrust and power outputs as F thrust and P muscle , respectively. The distance traveled in one gait period is named the stride length λ, and the mean velocity is given by |v|.
Previous studies have defined efficiency as the ratio of work done by the force thrust on the environment to either distance traveled [22] or to work done to drag a straight crawler in the frictional medium [18]. Here, since we use the internal muscle torque as the actuation mechanism, it is possible to directly quantify the efficiency of the crawler by comparing the rate of work done by the thrust force to the mean power expenditure by the muscles over one gait cycle and define the efficiency η as
Figure 4 shows computed steady locomotion results for different internal damping α values. Here the environmental friction is transverse dominant with κ 1 = 0. We find that for low internal damping, each cycle produces motion opposite the muscle moment wave propagation direction. At the value of α ≈ 0.5, the crawler undulates in place with no net locomotion. Damping values beyond unity produce first increasing stride length in the direction of torque propagation and then quickly drop off to zero. The crawler achieves the maximum efficiency of ≈0.024 at the critical value of α ≈ 6.
In Fig. 4, the shape ratio, d tail /d head , is around 0.8 for cases with low internal damping. The shape ratio increases FIG. 4. The dependence of the locomotion performance on the internal damping value in a friction environment where the normal friction force dominates the tangential friction (κ 1 = μ τ /μ n = 0). The vertical axis is stride length λ scaled by body length L, efficiency η in percent, and lateral excursion shape ratio d tail /d head , depending on the curve. It is observed that the longest stride displacement is associated with a relatively larger damping value and is accompanied by increased efficiency of η ≈ 2.4%. to its maximum of approximately 4 near the critical damping condition before returning to approximately 1 at larger α. In other words, for sufficiently low or high damping cases, the head deviates approximately as much from the direction of travel as the tail. However, the tail deviates approximately four times as much as the head at the critical damping point. The results are expanded to a larger range of the environmental and damping effects. Figure 5 shows how the stride displacement modifies over a large 40 × 40 parameter space of damping and environmental frictional ratio, μ τ /μ n . We find that the stride displacement crosses zero in the parameter space at μ τ /μ n ≈ 0.5, and for α ≈ 0.5. These zero crossings are denoted with a thick dashed line. Additionally, there is a saddle point where these zero crossings intersect. Moreover, locomotion that proceeds in the direction opposite to the torque wave propagation is referred to as "forward" in the figure and highlighted in light green. Similarly, locomotion displacement in the same direction as the torque wave propagation is labeled as "reverse" and highlighted in darker green. These features divide the parameter space into four quadrants. The global maximum of the stride length is obtained in the bottom right quadrant with higher damping and transverse dominant friction The second optimal stride length is obtained in the bottom left quadrant with low damping and transverse dominant friction condition. The top left quadrant has a much smaller local maximum for stride length that corresponds to high tangential friction and moderate damping. High damping Figure 6 shows how the shape ratio d tail /d head changes with respect to α and μ τ /μ n . Taking Figs. 5 and6 together, we see that a high shape ratio is associated with regimes with small stride lengths and often accompanied by no net locomotion. Conversely, a low shape ratio over the whole parameter space corresponds to longer stride lengths, except for the unique case of high tangent friction and high damping.
Figure 7 depicts the locomotion over one phase at the points of interest [regions (a)-(d)] from Fig. 5. Locomotion is slow in the tangentially dominated regions shown in Fig. 7 for regions (a) and (b); a salient difference is that the crawler's body exhibits more lateral motion for the lower-damping case FIG. 6. Lateral excursion ratio as a function of environmental friction and damping. We quantify that effect with the shape ratio or lateral excursion ratio d tail /d head (as defined in Fig. 3).
of region (a). In the case of region (b), in contrast, the crawler undergoes approximate rigid body pitching rotation.
For regions (c) and (d), Fig. 7 illustrates the gait changes that cause the crawler to change direction as a function of damping. In the top row graph for region (c), the crawler moves from left to right (dark to light), opposite to the moment propagation direction. A different direction of motion is observed in the top row graph for region (d), where it is shown that the crawler is moving from right to left in the same direction as the moment wave. A gait transition happens between left-to-right and right-to-left motion as the coordinate in the friction-damping parameter space crosses over the zero displacement line. We note that in this region, the tail has a significant deviation outside the path of the head. Region (d) in Fig. 7 offers the highest stride length gait. Here, the amplitude of the crawler is significantly attenuated and characterized by a lower curvature amplitude and higher wavelength.
The curvature ( b ) contour plots in Fig. 7 show that in both the highly damped cases, i.e., regions (b) and (d), isoclines are increasing from left to right, showing that curvature waves are propagating from tail to head. However, we see the opposite isocline slope in regions (a) and (c).
The thrust distribution plots, F e (s, t ) • t v , in Fig. 7 show that region (a) is uniquely dysfunctional in that very high thrust is being generated but is being dissipated entirely by corresponding negative forces. Moreover, we can see that even though the stride of region (d) is longer than that of region (c), the locomotion generates lower forces. This is explained by the reduced gait wave-form amplitude, which exposes the crawler to less drag force.
The previous section showed the dependence of a crawler's gait on two parameters identified to be most influential for locomotion, the damping parameter α and the environmental friction ratio μ τ /μ n . In Fig. 8, the left panel provides a summary of the damping-friction observed response regimes, and the right panel shows the dependence of stride length on the remaining dimensionless groups: μ n , Q, and ẼI. This sensitivity analysis reveals that the fastest forward gait shown in region (c) is not sensitive to variation in environmental friction, internal stiffness, or internal actuation. Likewise, the reverse gait shown in region (d) is not sensitive to either internal torque or material stiffness. However, this fastest gait is sensitive to the normal friction coefficient, such that low normal friction can cause the crawler to come to a complete stop in this regime and, similarly, a 50% increase can cause the stride length to double. For the transverse-dominated friction regimes, regions (a) and (c) in Fig. 8, an increase of bending moment or normal friction, or reduction of stiffness causes an increase in the stride length. This is the same trend that agrees with the effect of reducing the μ τ /μ n ratio in Fig. 5.
This work using the viscoelastic beam equations shows that lateral undulation is effective when normal or lateral friction is very high relative to tangential friction. This finding is consistent with the prevalence of wheels on snake robots that move using mainly lateral undulation since wheels provide a low value for μ τ /μ n . This is also consistent with the findings of a numerical study by Hu and Shelley [39] and the observation of Jayne that terrestrial snakes in isotropic friction environments tend to favor other gaits such as concertina over lateral undulation [2].
Considering these results in light of biological organisms, the stride length minima in the parameter space at μ τ /μ n ≈ 0.5 is roughly on the edge of the friction regimes measured for the scales of unconscious snakes (μ τ /μ n ≈ 0.6) [39]. This places snake-type crawlers on the edge of the transition between forward and backward with very low global stride displacement. In this friction-damping regime, the crawling is similar to region (a) in Fig. 7 and transitions to the uniform traveling wave similar to region (c) in Fig. 7 when damping is reduced, and the environmental friction becomes dominant in the transverse direction. As noted before, pure lateral undulation (that is, without weight redistribution) is known to be ineffective for snakes on isotropic surfaces. Moreover, in a comparable friction regime, Alben [22] showed an optimal gait to be very different from a uniform traveling wave (see Fig. 10 of Ref. [22] with μ t = 1.5 and μ t = 2). However, damping can still steer forward and backward motions in this regime, albeit with small stride displacements. More work is needed here since the effect of scale angle activation is thought to be significant for terrestrial snakes.
Comparing the effect of damping in Fig. 7, especially in regions (c) and (d), we see that internal damping stretches out the body of an undulating crawler. From a locomotion standpoint, this has the benefit of increasing stride length by reducing the curvature. To the extent that these crawlers' interactions with the environment are highly anisotropic, and more specifically the friction force in n is dominant, one may think of them as being similar to a skate on ice. In the same way that an ice skater can modulate their skate angle to achieve optimal motion, a crawler with large internal damping and anisotropic frictional forces can benefit from the reduced amplitude of the body wave to expose its body to less drag as defined in Eq. ( 15). This drag reduction explains the increased performance in this region and further suggests that the modulation of a robotic crawler's body curvature amplitude can be leveraged to achieve higher locomotion efficiency, perhaps in a burst-and-coast locomotion pattern similar to what is employed by some fishes [40].
Viscoelasticity in biological organisms varies throughout the gait cycle due to muscle activity and has been quantified by several studies. Measurements of the viscoelasticity of silver lamprey, Ichthyomyzon unicuspis, bodies show that it has a damping value of 0.22 kPa s. The same study suggests that the American eel, Anguilla rostrata, has damping 50 times greater [41]. A study on dog sharks, Squalus acanthias, reported a damping value on the order of 10 -6 N m 2 s [42]. Although these aquatic animals have distinct body structures that may give them internal damping properties different from terrestrial crawlers, we use these measurements to give context about where biological organisms are situated in the damping continuum studied. To do so, we scale these values according to our nondimensionalization, which places these animals in the α range 10 -3 -10 0 , roughly in the middle of our damping axis. Previous numerical studies have used damping (scaled by our nondimensionalization) of approximately 0.6 [43].
The damping also changes the crawler's shape wave form by changing the crawling direction from opposite the torque wave propagation to moving in the same direction as the torque wave. Since the transverse dominant friction regime is similar to wheeled snake robots, this phenomenon could be used to control the direction of a robotic crawler by changing the internal damping during runtime, perhaps through a modulated aperture inside the crawler's body. Many mechanisms and fluidic circuitry can now be built into the material of a soft robot through the use of specialized manufacturing processes [44].
Transverse dominant friction is also comparable to organisms moving in sedimentary media. Some organisms on various length scales, such as razor clams and polychaete worms, can move forward and backward in sedimentary media [45]. Recently, biologically inspired robots have been developed that mimic this locomotion [46]. Using internal damping as a control strategy to change the gait of a crawler in transverse dominant locomotion regimes might be useful in these soft robots to alter their burrowing direction. Considering the nondimensional parameters of Eq. ( 11), damping is a more significant factor at smaller length scales. Thus this locomotion modification strategy is especially relevant for small robots.
The effect of viscoelastic properties on the locomotion properties of a flexible crawler is studied here. A traveling wave is employed as the body actuation, and the swimmer's response at different environmental friction conditions is quantified. The study reveals an optimal damping value for the crawlers with transverse dominant friction forces. This condition can be leveraged to improve the performance of certain bioinspired robots, such as wheeled snake robots, and systems that can adjust the normal friction coefficient through surface interface manipulation. The results reveal that the internal damping reduces the amplitude of a crawler's undulatory wave form and changes the direction of curvature propagation with an accompanying change in locomotion direction. The unique functionality of the internal damping in switching the direction of the locomotion can be used in future robotic design to enable a seamless transition between forward to backward motion by controlling internal damping.