A Terradynamics of Legged Locomotion on Granular Media

See allHide authors and affiliations

Science  22 Mar 2013:
Vol. 339, Issue 6126, pp. 1408-1412
DOI: 10.1126/science.1229163

Walking on Sand

Studies of objects moving through air or water have provided detailed models for designing objects with better flow dynamics. Examples include aircraft fins and wings, robots used as underwater probes, and even swimsuits to enhance swimmers' competitiveness. Much less is known about the mechanics of moving objects on or in materials that themselves have non-uniform dynamics. For example, when walking on a granular medium like sand, the moving leg and foot may penetrate to varying depths with small changes in material properties. Li et al. (p. 1408; see the Perspective by Hunt) study this system and develop a model for predicting the motion of legged bodies on granular media for a range of leg shapes and motion speeds. Factors that complicate the motion include leg shape and size and movement direction, as well as the size shape, and uniformity of the granular material.


The theories of aero- and hydrodynamics predict animal movement and device design in air and water through the computation of lift, drag, and thrust forces. Although models of terrestrial legged locomotion have focused on interactions with solid ground, many animals move on substrates that flow in response to intrusion. However, locomotor-ground interaction models on such flowable ground are often unavailable. We developed a force model for arbitrarily-shaped legs and bodies moving freely in granular media, and used this “terradynamics” to predict a small legged robot’s locomotion on granular media using various leg shapes and stride frequencies. Our study reveals a complex but generic dependence of stresses in granular media on intruder depth, orientation, and movement direction and gives insight into the effects of leg morphology and kinematics on movement.

The locomotion of animals (1) and devices (24) emerges from the effective interaction of bodies and/or appendages with an environment. For flying in air and swimming in water, there is a history of theoretical predictive models (35) to describe the complex interactions between the locomotor and the surrounding fluids, based on the fundamental equations for fluid flow, the Navier-Stokes equations. These models have not only allowed understanding of the movement of a variety of aerial and aquatic organisms (5) [such as bacteria and spermatazoa (6), insects (7), birds (8), and fish and whales (9)] and their functional morphology, evolution, and ecology (9, 10), but also advanced the engineering design of aircraft (3), marine vehicles (4), and flying (11) and swimming (12) robots. For running and walking on ground, studies using solid ground such as running tracks and treadmills have inspired general models (13, 14); building on these models, researchers have begun to apply dynamical systems theory (15). In these studies, the leg-ground interaction was often approximated as a point contact on a rigid, flat, and nonslip ground (1315).

Many small legged animals (1619) [and increasingly robots (2023)] face the challenges of moving on natural substrates such as sand (16, 17, 21), gravel (16, 20), rubble (20), soil (20, 22), mud (17, 20), snow (18, 20), grass (20, 22), and leaf litter (19, 20, 22), which, unlike solid ground, can flow during movement when a yield stress is exceeded. The complexity of the interactions with such “flowable ground” may rival or even exceed that during movement in fluids. For example, recent studies of legged animals (16) and robots (21) moving on granular media [collections of particles (24)] such as sand and gravel (Fig. 1, A and B) have demonstrated that at an instant of time during a step, each element of a leg moves through the substrate at a specific depth, orientation, and movement direction, all of which can change over time (16, 21). Furthermore, the leg interacts with a material that can display both solid-like and fluid-like features (24) (Fig. 1, C and D). Compared to the theories of aero- and hydrodynamics, predictive models are less well developed for calculating forces and predicting legged locomotion on such flowable ground (16, 21).

Fig. 1

Examples of legged locomotion on flowable ground. (A) A zebra-tailed lizard running on sand (16). (B) A biologically inspired RHex robot (22) walking on dirt [Photo credit: Galen Clark Haynes, Aaron M. Johnson, and Daniel E. Koditschek, University of Pennsylvania]. Dashed boxes in (A) and (B) indicate the regions of leg-ground interaction shown in (C) and (D). Schematic of leg-ground interaction for (C) a hind foot of a zebra-tailed lizard (16) and (D) a c-leg of a RHex robot (21) during a step on granular media. Dashed, solid, and dotted tracings are leg positions at early, mid-, and late stance. Bars and arrows indicate local orientations and movement directions of leg elements. The gray area is the granular substrate.

Research in the field of terramechanics (2) has advanced the mobility of off-road vehicles on flowable ground such as sand and soil. These models were developed for large wheeled and tracked vehicles, which sink only slightly into the substrate (2, 25). Thus, in terramechanical models, interaction with the ground is approximated as the indentation of a horizontal, flat, rectangular plate (2, 25). It was a breakdown of this flat-plate approximation, however, that led to overpredicted speeds for small vehicles such as the Mars rovers, whose small wheels have substantially curved ground contact interfaces (25). Because leg-ground interaction on flowable substrates is a more diverse, complex, and dynamic process (16, 21) than the flat-plate indentation, terramechanics is not expected to apply to legged locomotors on flowable ground (2).

Granular materials such as sand and gravel (24), home to a variety of small desert animals (16, 17, 26), have proved to be a promising model substrate for studying legged locomotion on flowable ground (16, 21). Despite their diversity in particle size, shape, density, friction, polydispersity, and compaction (24), dry granular media are relatively simple as compared to media such as soil and mud, because granular particles interact purely through dissipative, repulsive contact forces and have no cohesion (24). In addition, the penetration resistance of granular media can be repeatably controlled using laboratory apparatus such as an air-fluidized bed (16, 21, 26).

To begin to create a “terradynamics” that allows the prediction of legged locomotion on a flowable ground, we hypothesized that the net forces on a leg (or a body) moving in a granular medium in the vertical plane could be approximated by the linear superposition of resistive forces on infinitesimal leg (or body) elements. Our hypothesis was inspired by our recent success in applying the methods of resistive force theory (6) to predict the forces and movement during the limbless locomotion of a lizard swimming in sand (26) and in describing the drag (26, 27) and lift (27) on simple objects moving in granular media at fixed depths. In these studies, the linear superposition was valid for intruders moving in granular media in the horizontal plane at low enough speeds [for example, ≲0.5 m/s for 0.3-mm glass particles (26)], where intrusion forces are dominated by particle friction (insensitive to speed) and non-inertial (26). However, it was unclear whether linear superposition could apply to legs (or bodies) of complex morphology and kinematics moving in the vertical plane.

To measure resistive forces for leg elements, we moved a thin rigid plate (of area A) in granular media in the vertical plane at 1 cm/s and measured lift fz and drag fx on the plate (in the continuously yielding regime). We determined vertical and horizontal stresses σz,x = fz,x/A as a function of the plate’s depth |z| below the surface, angle of attack β, and angle of intrusion γ (Fig. 2A and movie S1) (28). To test the generality of our resistive force model, we used three dry granular media of various particle size, shape, density, and friction, prepared into flat, naturally occurring, loosely and closely packed states (16, 21, 26) (supplementary text section 2, fig. S3, and table S1). Slightly polydispersed near-spherical glass particles 0.3 and 3 mm in diameter [covering the particle size range of natural dry sand (~0.1 to ~1 mm) (29)] and rounded, slightly kidney-shaped poppy seeds (0.7 mm in diameter) allowed us to probe general principles for naturally occurring granular media of high sphericity and roundness [such as Ottawa sand (30)]. We discuss at the end of the paper possible effects of particle nonsphericity and angularity also found in many natural sands (30).

Fig. 2

Measurement of resistive forces in granular media in the vertical plane using a plate element (28). (A) Lift fz (blue arrow) and drag fx (green arrow) on a thin rigid plate of area A moving in granular media (gray area) in the vertical plane at speed v = 1 cm/s were measured as a function of the plate’s depth |z| below the surface, angle of attack β, and angle of intrusion γ. The granular media were fluidized (and then compacted when a closely packed state was prepared) using an air-fluidized bed (26) before each intrusion (γ 0) and extraction (γ 0). g is gravitational acceleration. (B) Vertical (blue curve) and horizontal (green curve) stresses σz,x = fz,x/A versus |z| for representative intrusion and extraction using (β, γ) = ±(π/6, π/4) (movie S1). Blue and green dashed lines are linear fits with zero intercept over intermediate depths at which the plate was fully submerged and far from the bottom of the container. Horizontal arrows on top indicate the range of leg depths in Figs. 3 and 4. (C) Vertical and (D) horizontal stresses per unit depth αz,x [slopes of dashed fit lines in (B)] versus β and γ. Plate schematics with arrows denote representative orientations and movement directions. Black curves indicate where αz,x = 0. The shaded areas indicate where αz (or αx) is not opposing the plate’s vertical (or horizontal) velocity. Circles indicate αz,x from data shown in (B). Arrows above and below the color bars indicate directions of αz,x for positive and negative values.

In all media tested, we observed that for all attack angles β and intrusion angles γ, σz,x were nearly proportional to depth |z| when the plate was fully submerged and far from the bottom of the container (Fig. 2B and movie S1). This is because friction-dominated forces are proportional to the hydrostatic-like pressure in granular media. Therefore, we modeled the hydrostatic-like stresses asσz,x(|z|,β,γ)={αz,x(β,γ)|z|ifz<00ifz>0 (1)where αz,x are vertical and horizontal stresses per unit depth (slopes of dashed fit lines in Fig. 2B). We found that in all media tested, αz,x depended sensitively on both attack angle β and intrusion angle γ (Fig. 2, C and D, fig. S4, and additional data table S5). αz (or αx) was opposing the plate’s vertical (or horizontal) velocity for most but, counterintuitively, not all β and γ (exceptions are indicated by the shaded regions). For almost all attack angles β, αz,x had larger magnitudes (|αz,x|) for intrusion angle γ ≥ 0 than for γ ≤ 0; i.e., it was harder to push the plate into granular media than to extract it. For all intrusion angles γ except γ = ±π/2, |αz,x| were asymmetric to attack angles β = 0 and β = ±π/2; i.e., only when the plate moved vertically were stress magnitudes the same for vertically or horizontally mirrored orientations (e.g., β = ±π/6). These asymmetries are a result of gravity breaking symmetry in the vertical plane and differ from the case in the horizontal plane (26). Our resistive force measurements are an advance from previous force models based on the flat-plate approximation used in many terramechanical models (2, 25), which capture only the dependence of stresses on intruder depth, but not on its orientation or movement direction (supplementary text section 1 and fig. S1). Despite their different magnitudes and subtle differences in shape, the overall profiles of stresses (per unit depth) αz,x(β, γ) were similar for all media tested. Furthermore, these stress profiles could be approximated (to the first order) by a simple scaling of generic stress profiles (supplementary text section 3, figs. S6 and S7, and tables S2 and S3).

We next tested our hypothesis that forces on a complex intruder moving in granular media in the vertical plane could be approximated by the linear superposition of forces on all intruder elements. We measured the net lift Fz and thrust Fx on thin rigid model legs rotating about a fixed axle [simulating a tethered body (7, 11)] through granular media in the vertical plane at ~1 cm/s (Fig. 3 and movie S2) (28). We then compared them to predictions from the resistive force model by the integration of stresses over the legs (movie S3)Fz,x=Sσz,x(|z|s,βs,γs)dAs=Sαz,x(βs,γs)|z|sdAs(2)where S is the leading surface of the leg; dAs, |z|s, βs, and γs are the area, depth below the surface, angle of attack, and angle of intrusion of infinitesimal leg elements; and αz,xs, γs) are element stresses per unit depth (interpolated from data in Fig. 2, C and D). To test the robustness of our force model, we used three model legs of different geometries [with the same maximal leg length 2R (28)]: a RHex robot’s c-leg (21, 22), a flat leg, and a reversed c-leg (Fig. 3, A to C). In model calculations, each leg was divided into 30 elements.

Fig. 3

The resistive force model predicts forces on intruders of complex morphology and kinematics moving in granular media (28). Three thin rigid model legs of different geometries, (A) a c-leg, (B) a flat leg, and (C) a reversed c-leg, were rotated about a fixed axle at a hip height h through granular media (gray area) in the vertical plane at an angular velocity ω, generating leg speeds of v ~ 1 cm/s, and net lift Fz (blue) and thrust Fx (green) were measured as a function of leg angle θ (movie S2). All three legs had identical maximal length 2R from the axle. g is gravitational acceleration. (D to F) Fz,x versus θ on the three legs measured in the experiment (solid curves) and predicted by the resistive force model (dashed curves) using Eq. 2 (movie S3). Insets: Fz,x versus θ from experiment (solid curves) versus predicted (dotted curves) using Eq. 2 and previous force models in which stresses depended only on depth (supplementary text section 1 and fig. S2).

In all media tested (Fig. 3, D to F, and fig. S5), we observed that for all three legs, the measured net lift and thrust Fz,x as a function of leg angle θ (solid curves) were asymmetric to the vertical downward direction (θ = 0), and were larger during intrusion (θ ≤ 0) than during extraction (θ ≥ 0). Peak Fz,x were largest on the c-leg and smallest on the reversed c-leg. The reversed c-leg experienced significant negative lift (suction force, Fz < 0) during extraction. For all media tested, our resistive force model predicted Fz,x versus θ for all three legs (dashed curves), capturing both the magnitudes and asymmetric profiles. The relative errors of peak forces between data and model predictions were within 10% for the c-leg in four media tested, and within 33% for all three legs in all media tested. The accuracy of our resistive force model was significantly better than that of previous force models in which stresses depended only on depth (Fig. 3, D to F, insets; supplementary text section 1 and fig. S2). Furthermore, our resistive force model revealed that the c-leg generated the largest forces, because its morphology allowed leg elements to not only reach deeper depths but also access larger stress regions in Fig. 2, C and D (particularly for elements at large depths).

Our discovery of the insensitivity of the stress profiles to particle properties (fig. S4) has practical benefits: For granular media of near-monodispersed, near-spherical, rounded particles, as an alternative to measuring αz,x for all attack angles β and intrusion angles γ in the laboratory, one can simply perform a single measurement [of αz(0, π/2), using a horizontal plate penetrating vertically downward] to infer all αz,x(β, γ) by a scaling routine (supplementary text section 3 and fig. S9) and predict forces (with a small loss in accuracy for the c-leg and the flat leg, but a larger loss in accuracy for the reversed c-leg, fig. S8).

We tested the ability of our resistive force model to predict legged locomotion. We chose to study the locomotor performance (speed) of a small RHex-like robot (22) (Fig. 4A, top, and movie S4) moving on granular media (28). The robot’s six legs rotated nearly entirely in the vertical plane during locomotion, and its small size ensured that leg intrusion speeds were low enough for particle inertia to be negligible. We chose poppy seeds as the test granular medium, because the grains were both small enough be prepared in our fluidized bed track (21) and large enough to not jam the robot’s motor and gear trains. The robot’s legs had a similar friction coefficient with poppy seeds to that of the model legs and were sufficiently rigid so that they experienced negligible bending during movement. (28).

Fig. 4

A multibody dynamic simulation using the resistive force model predicts legged locomotion on granular media (28). (A) Side views of a small RHex-like robot (movie S4) at mid-stance during locomotion on granular media, using c-legs (left) and reversed c-legs (right) in the experiment (top, movie S5) and simulation (bottom, movie S6). Arrows in the simulation indicate element forces on one tripod of legs. g is gravitational acceleration. (B) Forward speed vx versus time t from two representative runs using c-legs (red, stride frequency f = 2.0 Hz , curvature 1/r = 1/R ′) and reversed c-legs (blue, f = 2.2 Hz , 1/r = −1/R ′). (C) Average forward speed Embedded Image versus f using legs of seven curvatures 1/r transitioning from reversed c-legs to c-legs (inset), where r is the radius of curvature, 2R′ is the maximal length of the robot legs, and the minus sign denotes reversed legs. (D) Embedded Image versus 1/r at f = 1, 2, 3, and 4 Hz. In (B) to (D), solid and dashed curves indicate the experiment and simulation, respectively. Error bars in (C) denote ±1 SD (<<1 cm/s in the simulation). Experimental data in (D) are interpolated from those in (C) (hence no error bars). Red and blue arrows in (C) and (D) indicate averages from data shown in (B).

Unlike the sand-swimming lizard, which moves within granular media quasistatically (thrust and drag are always roughly balanced) (26), legged locomotion on the surface of granular media is dynamic (forces are not always instantaneously balanced). As a result, the resistive force theory (which solves for speed by balancing forces) (6, 26) cannot be directly applied. Thus, to use our resistive force model to calculate robot speed, we developed a three-dimensional multibody dynamic simulation of the robot (Fig. 4A, bottom) (28). The simulated robot had the same body and leg morphology and used the same alternating tripod gait as the actual robot and had its motion constrained in the vertical plane. We divided each body plate and leg into 30 elements. The velocity v and angular velocity ω of the simulated robot’s body were calculated by{v(t+dt)=v(t)+Fmdtω(t+dt)=ω(t)+NIdt(3)where F and N are the sum of net forces and torques on all the six legs and the body exerted by the granular medium, calculated from our resistive force model by the integration of stresses over each leg and the body using Eq. 2; m and I are the robot’s mass and moment of inertia; and t and dt are time and time step. To test the robustness of our resistive force model and simulation, we used legs of seven geometries with different curvatures 1/r (given maximal leg length 2R′) (Fig. 4C, inset) and varied stride frequency f to up to 5 Hz (28).

We observed similar robot kinematics (Fig. 4A) and forward speed vx versus time t (Fig. 4B) in both the experiment (movie S5) and simulation (movie S6). The robot moved faster and penetrated its legs less deeply during stance using c-legs (Fig. 4A, left; Fig. 4B, red) than using reversed c-legs (Fig. 4A, right; Fig. 4B, blue). Average forward speed v¯x increased with stride frequency f for legs of all curvatures 1/r and was lower at any f using legs of negative curvatures than using legs of positive curvatures (Fig. 4, C and D). The agreement between experiment and simulation in v¯x(f, 1/r) was remarkable: Errors were within 20% for 90% of the f and 1/r tested and within 35% for all the f and 1/r tested. Simulation using our scaling routine also achieved reasonable accuracy (supplementary text section 5 and fig. S13). This was an improvement over simulation using previous force models in which stresses depended only on depth (fig. S13). Our resistive force model and simulation revealed that the robot moved faster using c-legs than using reversed c-legs, because whereas the c-legs penetrated less deeply, their elements accessed larger stress regions in Fig. 2, C and D, resulting in larger leg lift (fig. S14) and smaller body drag. Our model and simulation also allowed the prediction of ground reaction forces on granular media (Fig. 4A, red arrows, and fig. S13), which would be difficult to measure otherwise. Furthermore, our model and simulation predicted that using arc-like legs (given maximal leg length 2R′) of an optimal curvature of 1/r = 0.86/R′, the robot would achieve maximal speed of v¯x = 72 cm/s (≈ 5 body length/s) at 5 Hz. Our approach affords significant reduction in the computational time needed to model movement on granular media. For example, relative to our multiparticle discrete element method (DEM) simulation of movement on granular media (23, 27), our simulation using the resistive force model can achieve a factor of 106 in speed-up (e.g., 10 s versus 30 days using DEM to simulate 1 s of locomotion on a granular bed of 5 × 106 poppy seeds).

We close with a brief discussion of the limitations of our model. We tested the predictive power of our scaling routine (supplementary text section 3 and fig. S9) for two highly polydispersed, nonspherical, highly angular natural sands (supplementary text section 4, figs. S10 to S12, and table S4). We found that the model accuracy for natural sands was slightly worse than found in the glass spheres and poppy seeds (for example, 35% versus 20% error in peak Fz for a rotating c-leg). As was the case for the near-spherical granular media tested, the functional forms of forces on the c-leg and the flat leg were still well captured by our scaling routine. Furthermore, the overestimation was not affected by reducing the polydispersity of the natural sand. This suggests that the nonsphericity and angularity of natural sand particles (30) may be the cause of this overestimation, which may require additional model fitting parameters and scaling factors. Our model is intended for dry sand [~0.1 to ~1 mm in particle diameter (29)] and may not work for dry cohesive powder (≤ ~0.01 mm) (31). We do not expect our model to work if the particle size approaches a characteristic length of the locomotor (for example, ~1-cm particles for our robot of ~1-cm foot size), so that the continuum assumption breaks down and particles effectively become “boulders.” We also do not expect our model to capture wet, cohesive flowable media such as soil and mud.

We have developed a new approach to predicting legged locomotion on granular media. This terradynamics relies on new resistive force measurements and linear superposition (6, 26, 27). The general profiles of these resistive force measurements are insensitive (other than magnitudes) to a variety of granular media composed of slightly polydispersed, approximately spherical, rounded particles. Our terradynamics may not be limited to legged locomotion, because the integration of stresses should in principle work for devices of other morphology and kinematics, such as wheels, tracks, and earth movers moving on granular media (2, 25). For the particle types tested here, an important addition to our model would be to capture three-dimensional effects (16) and spatial and temporal variation in compaction (21) and slopes (17), and test its validity in the high-speed “inertial fluid” regime (when leg intrusion speed is ≳1 m/s, at which particle inertia dominates forces) (23). Our resistive force model also provides opportunities to test and develop new physics theories of dense granular flow (32). Finally, we envision that, in concert with aero- and hydrodynamics (312), a general terradynamics of complex ground will not only advance understanding of how animals move (1) at present (510, 1319, 26) and in the past (17, 33), but also facilitate the development of robots with locomotor capabilities approaching those of organisms (11, 12, 2023).

Supplementary Materials

Materials and Methods

Supplementary Text

Figs. S1 to S14

Tables S1 to S4

References (3438)

Movies S1 to S6

Additional Data Table S5

References and Notes

  1. Materials and methods are available as supplementary materials on Science Online.
  2. Acknowledgments: We thank Y. Ding, P. Umbanhowar, N. Gravish, G. Meirion-Griffith, S. Sharpe, H. Komsuoglu, D. Koditschek, and R. Full for discussions; J. Shen for assistance with robot modification; P. Masarati for multibody dynamic simulator support; S. Sharpe for measuring the angle of repose of 3-mm glass spheres and assistance with photography; P. Umbanhowar and H. Marvi for natural sand collection; and all the members of the Complex Rheology And Biomechanics Lab at Georgia Tech for general assistance. This work was supported by the Burroughs Wellcome Fund, the Army Research Laboratory Micro Autonomous Systems and Technology Collaborative Technology Alliance, the Army Research Office, and the NSF Physics of Living Systems program. C.L. was partially supported by a Miller Research Fellowship from the Miller Institute for Basic Research in Science of the University of California, Berkeley. The authors declare that they have no competing interests. C.L. designed the study, performed resistive force measurements, and performed robot experiments; C.L. and T.Z. performed model calculations; T.Z. performed robot simulation; D.I.G. oversaw the study; and C.L. and D.I.G. wrote the paper.

Stay Connected to Science

Navigate This Article