The Effectiveness of Resistive Force Theory in Granular Locomotion A)

Resistive force theory (RFT) is often used to analyze the movement of microscopic organisms swimming in fluids. In RFT, a body is partitioned into infinitesimal segments , each of which generates thrust and experiences drag. Linear superposition of forces from elements over the body allows prediction of swimming velocities and efficiencies. We show that RFT quantitatively describes the movement of animals and robots that move on and within dry granular media (GM), collections of particles that display solid, fluid, and gas-like features. RFT works well when the GM is slightly polydisperse, and in the " frictional fluid " regime such that frictional forces dominate material inertial forces, and when locomotion can be approximated as confined to a plane. Within a given plane (horizontal or vertical) relationships that govern the force versus orientation of an elemental intruder are functionally independent of the granu-lar medium. We use the RFT to explain features of locomotion on and within granular media including kinematic and muscle activation patterns during sand-swimming by a sandfish lizard and a shovel-nosed snake, optimal movement patterns of a Pur-cell 3-link sand-swimming robot revealed by a geometric mechanics approach, and legged locomotion of small robots on the surface of GM. We close by discussing situations to which granular RFT has not yet been applied (such as inclined granular surfaces), and the advances in the physics of granular media needed to apply RFT in such situations. C 2014 AIP Publishing LLC.


I. INTRODUCTION
][3][4][5][6][7][8][9][10] However, compared to the research efforts to understand swimming 11 (Fig. 1(a)), and flight 12 (Fig. 1(b)) in fluids, movement in and on terrestrial substrates (Figs.1(c)-1(f)) is less understood.In aquatic and aerial locomotion, the lift and drag forces during propulsive movements can, in principle, be calculated by solving the Navier-Stokes equations (often numerically). 12,13 hile traditionally considered rigid, many terrestrial environments are composed of rheologically complex "soft" materials that flow upon interaction (Fig. 2); examples include sand, mud, snow, and leaf debris.Such "flowable" terrestrial materials are not yet modeled at the level of Newtonian fluids.Modeling locomotion in such environments is therefore an open and challenging problem.
While little is known about the mechanics by which organisms move on or within flowable terrestrial substrates, there has been considerable effort applied to modeling the interaction of wheeled/tracked vehicles with terrain.This approach is called terramechanics, 14 and it models soil as (roughly speaking) a nonlinear damped spring. 15At each contact point on the wheel or tread, normal and shear stresses are computed, often using many fitting parameters (depending on the stress model 16,17 ).The stress projected to the vertical direction is then balanced with the load of the system to determine the static and dynamic "sinkage" of the wheel during slip.Terramechanics has advanced the off-road mobility of vehicles in the past 50 years on highly deformable soil.However, since the method was primarily developed for large, heavy vehicles (such as tanks), its applicability towards small, light-weight robots is unknown.Also, the interaction of geometries other than wheels and tracks (e.g., legs) with soil is not well modeled. 18here have been scattered efforts to monitor and model the movement of small organisms in soil, sand, and mud.Charles Darwin's final book focused on the actions of worms on soil. 19In the 1960s using capacitive sensing, Norris and Kavanau studied subsurface locomotion of four species of limbless reptiles (three snakes and one lizard). 20Trueman investigated the burrowing of mole crabs 21 and the razor clam Ensis 22 in soft substrates.Recently, Winter et al. studied the burrowing mechanism of razor clams in mud and found that the animal fluidizes the substrate to reduce drag. 23roups led by Jung 24 and Arratia 25 studied the swimming of nematodes in wet granular media and enhanced mobility was found when compared with fluids.The digging of polychaete under muddy sediments was studied by Che et al., 26 and they explained their results using fracture mechanics.Gidmark et al. 27 used x-ray to study the kinematics of the sandlance fish's burrowing process.
Our group has focused on movement in dry granular media, conducting studies of animal and robot locomotion on and within GM. [28][29][30][31] Such locomotion generally involves a body/limb of complex geometry moving in nontrivial patterns in a granular substrate.For example, S. scincus (Fig. 1(e)) has a wedge-shaped snout and a smooth body which generates a traveling wave pattern during swimming. 28This complex interaction is beyond the scope of traditional terramechanics (which assumes purely plastic interactions).Experiments and simulations have indicated that fluidlike aspects of the GM are important in the regimes relevant to sand locomotion (Fig. 3). 28,32 n such regimes particles are continuously pushed and stirred by the moving body, causing the media to behave like a dense frictional fluid.
To model our locomotion results, we were inspired by the approach developed by pioneers in fluid mechanics who did not solve (and, without computational fluid dynamics, could not) the full Navier-Stokes equations for complex flows in the presence of complex moving boundaries.Instead, simpler approximations were made in viscous fluids.The best known of these, called Resistive Force Theory (RFT), 33 assumes that the deforming body can be partitioned into segments, each experiencing drag, and that the flow/force fields from these segments are hydrodynamically decoupled and do not influence the fields of other segments.Therefore, the normal and tangential forces on a small element (Fig. 4) depend only on the local properties, namely, the length of the element ds, the velocity v and the orientation t (the fluid is homogeneous so position dependence is dropped).The net force on the swimmer is then computed from the integral where the functional forms of f ⊥ and f can in principle be determined from Stokes equations; however, experiments are often needed because analytical solutions to Stokes equations cannot be easily obtained even in simple cases (e.g., a finite length cylinder moving in the axial direction 34 ).In GM where there is no constitutive law, f ⊥ and f can only be determined from experiment or DEM (Discrete Element Method 32 ) simulation (details in Sec.II).
In viscous fluids, the RFT approximation has had some success when compared to experiment [35][36][37][38] .However, in the physiologically relevant regimes for which RFT was originally formulated (the movement of flagella for bacteria, for example), discrepancies between RFT and experiment have recently been observed; 39 to remedy this, more comprehensive theories have been developed, and these account for higher order corrections (see Ref. 39 for discussion).As we will show, RFT with empirically determined f ⊥ and f can accurately model (and predict aspects of) locomotion of animals and physical models of animals (robots) on and within dense GM. 28-31, 40, 41 This has been somewhat surprising, since we initially worried that interaction with a medium that can display features of a solid, 42 a compressible fluid, 43,44 and/or a gas would require models which accounted for the physics of each regime as well as transitions between them. 45

II. GRANULAR RFT, FORCE RELATIONS, AND SOLVING FOR LOCOMOTION
Granular RFT is simple to implement once input forces are specified.To apply RFT, we partition an arbitrary intruder into either cylindrical or plate-like elements; we then drag/penetrate (at constant speed) these elements in granular media, and measure the resistance force to determine the RFT force relations.In experiment, the intrusion speed is low enough (e.g., v 50 cm/s for 0.3 mm glass particles) that the reaction forces from GM are insensitive to speed 46 and increase with depth (described below).The forces increase with increasing prepared volume fraction, which ranges from 0.58 to 0.62 in our laboratory apparatus in 0.3 mm glass particles.
We have empirically determined the force relations for two cases: drag in the horizontal and vertical planes (with respect to gravity) (Figs. 5 and 6).We now describe some details associated with these experiments and the RFT integrals in the two different regimes.v direction of movement FIG. 4. Illustration of the basic idea behind Resistive Force Theory (RFT) for 2D swimming.The body propagates a wave of prescribed shape in the negative x-direction and experiences reaction forces which propel it forward.Each infinitesimal element ds on the swimmer's body is characterized by its tangent direction t (or normal direction n) and its velocity v direction; each element experiences a force dF ⊥, .In true fluids, these forces can be described by Stokes' law, while for granular media, they are measured in experiment.

A. Horizontal plane
We first introduced RFT to model subsurface swimming of the sandfish lizard S. scincus (see Sec. III B). 28 In the spirit of the Gray and Hancock formulation of the RFT in fluids, 33 the body was divided into segments approximated by cylinders.Since no constitutive force law was available for GM, f ⊥ , and f in Eq. ( 1) were determined from drag measurements (Figs.5(a)-5(c)).The RFT force integral for the horizontal plane can then be written with the empirical form: where kρg|z| is the effective hydrostatic-like pressure at the depth |z| (the pressure due to the weight of material with effective density ρ (defined to be the particle density times the volume fraction) above |z|, with k ≈ 2.5 for loosely packed 0.3 mm glass particles), r is the radius of the cylinder.
The anistropic factor of the normal force has the form with C = 1.8 and γ 0 = 13.8 • for the media considered.The locomotion modeling using RFT is discussed in Sec.III B.

B. Vertical plane
Recently we developed RFT for the vertical plane 31,40 .We now describe the penetration experiments by which we obtained the vertical plane force relations 31 .A thin rigid plate of area A was moved in GM in the vertical x − z plane at 1 cm/s.We measured the lift F z and drag F x on the plate in the continuously yielding regime (Fig. 6).By the RFT assumptions, the stresses σ z, x = F z, x /A are functions of the penetration depth |z| below the granular surface, the angle of movement direction γ and the angle of orientation β of the plate (Fig. 6(a)).In all the lab GM substrates (∼1 mm kidney-shaped poppy seeds and slightly poly-dispersed approximately spherical glass particles of 0.3 and 3 mm diameter) tested, we found that σ z, x increased approximately linearly with |z|, for all β and γ in the regime when the plate was fully submerged and far from the container boundaries (Fig. 6(b), shadowed region).We note that our measurements include the special case (β, γ ) = (0, π /2) of a plate pushed directly downwards. 47he stress may be written as Although little is understood about α z, x , the functional forms for all tested GM have similar characteristics (Fig. 7).After applying a discrete Fourier transform to α z,x (in its configuration space (β, γ )), the Fourier components (Fig. 7(b)) for all different substrates could be normalized by a single scaling factor ζ i , which indicated the "strength" (resistance) of the media.With the generic Fourier coefficients M 0 (Fig. 7

7(d))
, which apply to most dry GM (ζ i can be determined from a vertical penetration test such that (β, γ ) = (0, π/2)).The force predictions were compared with experiments and are discussed in Sec.III A.

C. Predicting motion
To model locomotion, we integrated the force relations with various models of animals and robots; in all studies, we prescribed the kinematics of the self-deformation of a body.The simplest models were those that were described analytically and for which we could ignore body inertia.Assuming a balance of forces and torques at all instants, we determined the translation of the center of mass and/or rotation of the body.The force balance assumption is a reasonable approximation (see Ding et al. 29 for discussion of this assumption) for subsurface locomotors.However, it is not a good approximation for the above-ground locomotors, since for the robots we studied, body inertia was non-negligible.Therefore, we also integrated the RFT relations into a multibody simulation environment, MBDyn, which solves equations of motion under constraints using external forces and torques from the RFT relations.The latter approach (called "terradynamics") achieved good accuracy in predicting legged locomotions of lightweight robots moving on homogeneous GM. 31

III. APPLICATIONS OF GRANULAR RFT
Below we list examples where RFT has been used to predict forces and model locomotion in dry, homogeneous (no variation in volume fraction in space), flat, pristine (not previously disturbed) GM.While this may appears restrictive, we will show that RFT quantitatively describes many seemingly complex granular locomotion situations, and optimal solutions obtained in the RFT are in accord with observed (and hypothesized optimal) biological locomotion strategies.Following this section, we give short descriptions of regimes to which RFT could be applied (but has not been thus far).

A. Lift and drag
There have been many studies of drag 46,[48][49][50] and penetration 47,51,52 forces in dry GM and several empirical force models have been proposed. 47,51 owever, such studies involved intruders of simple geometries (e.g., cylinder and sphere) undergoing simple horizontal and vertical movements in GM.Though DEM can predict forces in more complex situations (e.g., translation and/or rotation of complex-shaped legs inside GM during locomotion), 32 it is presently impractical to simulate large-scale problems where there are more than a few million particles.RFT, on the other hand, is a non-PDE based continuous model, and therefore does not suffer from issues associated with such computation.
The validity of the vertical plane RFT was tested by Ding et al. 40 in a study of drag induced lift forces in 3 mm diameter glass particles.An experimentally validated DEM was used to investigate the leading surface stress distributions on the intruders.The normal σ and tangential τ stresses on a cylinder at a particular surface orientation angle β were nearly equal to those on a plate with the same β at the same |z| (Fig. 8(a)).By summing the local stresses, the predicted net lift F z quantitatively matched direct measurements in simulation and experiment for rods of several different cross-sectional shapes. 40enetration forces for complex shapes (e.g., a c-shaped leg) undergoing non-trivial movements in the vertical plane (e.g., rotation in the vertical plane about a fixed axis (Fig. 8(b), left panel) were studied by Li et al. 31 and compared to RFT calculations.At any instant during the rotation, the net forces from vertical RFT are written as where α, β and γ follow the same convention in Section II B, ds is the area of the surface element, and the integral is over the leading surface S of the considered geometry.The RFT-predicted forces F z,x vs the rotation angle θ were in accord with experimental measurements in both magnitudes and asymmetric functional forms (Fig. 8(b), right panel).The relative error of peak forces between the data and model predictions were within 10% for a c-shaped leg in poppy seeds and also compared well in the other substrates tested (Fig. 7(a)).

B. Swimming: Biomechanics
We now explain how we used horizontal plane RFT to study the undulatory swimming mechanism of the sandfish lizard S. scincus (Fig. 1(e)). 28X-ray videos revealed that the body undulation pattern of the animal was well approximated by a planar sinusoidal wave that propagated from the head to the tail, while the center of mass moved forward (Fig. 9(a)) in approximately the horizontal plane.The RFT prediction of body wave curvature to achieve optimal undulation efficiency (the forward progress per cycle normalized by wavelength) (Fig. 9(c), red lines) fell close to the kinematics used by the animals (Fig. 9(c), red cross).In later work, 32  confirming that the animal swam in a frictional fluid.However, relative to DEM, RFT overestimated the performance of the swimmer by up to 30%.When force distributions along the swimmer's body were compared between RFT and DEM, 29 the major discrepancy was revealed to be due to transient phenomena in the initiation of drag in GM (see Sec. IV A).Recently, RFT was applied to model another sand-swimmer, the Mojave shovel-nosed snake C. occipitalis (Fig. 1(f)).X-ray measurements revealed that the animal also used a traveling wave of body curvature to burrow and swim in a granular substrate of 0.3 mm glass spheres.Compared with S. scincus, little slip was observed during its subsurface movement (Fig. 9(a), 9(b)).Here slip can be visualized as deviation from movement within a tube created by the swimmer.However, the animal did not strictly move within a tube (as would be the case for zero slip) (Fig. 9(b)).Since the slip was small but finite, in collaboration with Dr. Stephan Koehler we applied RFT to the snake (accounting for its lower skin friction and its different body undulation number ξ , defined as the average number of waves along the body).Like the sandfish the snake optimized its undulatory pattern to achieve optimal forward progress per cycle (Fig. 9(c), blue cross).Our calculations revealed that both swimmers also minimized their energy use (here the mechanical cost of transport).

C. Swimming: Neuromechanics
4][55] This is referred to as neuromechanical phase lags (NPL).To explain the origin of NPL, at least two approaches have been taken: A "bottom-up" approach that models many biocomponents (such as muscles and tendons of the body), analyzes the complex interactions among them, and integrates these models with external forces from the environment (e.g., fluid partial differential equations).A second, "top-down" approach attempts to guide the study from reduced order ("template" 5 ) models that seek general principles of the system's behavior.
We studied the NPL during kinematic sand swimming using the "top-down" approach. 41The nervous/musculoskeletal system was abstracted by assuming that passive body forces were small and that internal torque was synchronized with neural activation timing.As before (see Sec. III B), the body motion of the animal in the center of mass frame was prescribed as a sinusoidal wave in the horizontal plane.Because accelerations are assumed negligible in the over-damped environment (see Ref. 29 for discussion of this assumption), the center of mass motion ( Ṙ, θ ) can be numerically determined by solving the steady state force-torque equations: which use input forces from RFT (R and θ were dropped due to translational and rotational symmetry, respectively, because of homogeneity of the media in the horizontal plane).
Assuming that the torque on each portion of the body is balanced at any instant of time, i.e., τ muscle + τ ext = 0, the local torque distribution can be solved at any point of interest along the body.Because we assume that the time interval between neural activation and muscle force development is small compared with the undulation period of S. scincus, the non-zero value of τ muscle directly corresponds to neural activity; positive τ ext (or negative τ muscle ) corresponds to the onset of muscle contraction on right/left side of body (Fig. 10(a)).Therefore, only the signs of τ are used to predict muscle activation (Fig. 10(b)).
The muscle activation over the whole body displays a traveling wave pattern (Fig. 10(c), solid green lines).The phase speed (visualized as the inverse of the slope of the onset/offset curve) of the activation is smaller than the curvature wave speed (the inverse of the slope of the shaded region), i.e., phase lag naturally arises from the RFT model.Without corrections or passive body forces, the average phase difference between the beginning and ending of positive τ in the model compared with electromyography onset and offset in experiments is less than 5%, (Fig. 10(c) dashed lines and solid green lines).The agreement between experiment and model implies that NPL is intrinsic to swimming in continuous media because individual segments must act to generate torques that move all segments through the medium.

D. Swimming: Geometric mechanics
Geometric mechanics 56,57 takes advantage of symmetries to reduce the effective dimensionality of the dynamical system and facilitate interpretation of dynamics in terms of geometric concepts, such as areas, lengths, and curvatures.1][62] While the method was developed and tested for movement (including swimming in fluids and the movement of "snakeboards") where analytically describable linear dynamics could be derived, only recently was the concept extended to swimming in GM, with the help of new tools in geometric mechanics 62 and the granular RFT. 30 We demonstrate the key idea of applying geometric theory to locomotion in GM using the so-called "3-link swimmer" (Fig. 11(a)), originally introduced by Purcell to illustrate features of locomotion in low Reynolds number environments. 63The core assumption in the application of geometric mechanics to locomotion is that one can describe the center of mass motion only in the manifold of local variables (shape velocities, etc.).For any given shape with joint angles specified by the vector α = (α 1 , α 2 ) (Fig. 11(b)), the swimmer's body velocity ξ is assumed to be linearly proportional to its shape velocity α, such that the relationship between shape configuration, shape velocity, and body velocity can be expressed as where A(α) is referred to as the local connection (or Jacobian) matrix.
In GM, the local connection matrix was determined numerically using RFT.By integrating the RFT forces along the swimmer's body for combinations (α, α, ξ ) and solving for the implicit force equilibria F(α, α, ξ ) = 0, we found that for a given α, the components of ξ were nearly planar (linear) functions of the shape velocity α.This planarity means that, despite the lack of an analytical linear relationship between shape and body velocities, there exists a linear relationship between body and shape velocity, and thus a local connection can be generated.Using the best fit plane (least square), each row of A(α) was obtained (which can be visualized as a vector field on the shape space (α 1 , α 2 ), Fig. 11(c)).
A swimmer's gait is then visualized as a closed loop S in the shape space (Fig. 11(c), red circle; though in general the loop can be of arbitrary path) and the net displacement of the body can be obtained from integrals in the shape space), which is finally written as a "surface" integral of the curl over the enclosed area S. The derivation here is not strict, because Strokes theorem in principle only applies to systems where the integrations are commutative.However, in the 3-link swimmer rotation and translation do not commute.Thus, the non-Abelian Stokes' theorem was used in calculations.The curl in the non-Abelian case was replaced by the so called Curvature Constraint Functions (CCF, Fig. 11(d)) 64,65 in which an augmenting Liebracket was added to account for the noncommutativity, see Ref. 30 for further details.Agreement between theory, experiment, and simulation is excellent at small amplitudes and good at larger amplitudes (Fig. 11(e)).However, while the CCFs qualitatively predict forward movement at the highest amplitudes, quantitative agreement is lacking.From previous work, 62 we hypothesize that the cause of this error is a kinematic singularity (for certain joint configurations the local connection matrix is singular and the center of mass speed can only be maintained if the joint velocity goes to infinity).

E. Walking on the surface
Inspired by our success in modeling subsurface swimming, we next tested if RFT could model locomotion on the surface of a granular medium where the movement of legs are approximately confined to a vertical plane.We tested the approach on a small hexapedal robot (Xplorer, a ∼150g, ∼13 cm RHex-like robot) on ≈1 mm diameter poppy seeds. 31We observed similar robot kinematics (Fig. 12(a)) and forward speed v x versus time t (Fig. 12(b)) in experiment and simulation.To test the robustness of the theory, we further equipped Xplorer with 3D-printed legs of different curvatures (Fig. 12(c)).For low stride frequencies (f ≤ 5 Hz when the linear speed of the leg during locomotion is small and the substrate behaves like a frictional fluid) vertical-plane RFT accurately captured the performance of Xplorer for all the legs tested (Fig. 12(d)).Recently, we used the vertical-plane RFT (which we now term "terradynamics" to emphasize its utility as a companion to terramechanics for smaller locomotors with intrusions which are more complex than wheels or tracks) to predict the locomotion performance of the larger SandBot (∼2.3kg, ∼30 cm RHex-like robot, Fig. 12(e)) also walking on a bed of "aerated" 66 ≈1 mm diameter poppy seeds.Though air flow potentially adds more complications to the system, the vertical penetration resistance α(0, π/2, v air ) in experiments decreased approximately linearly with increasing air flow speed v air (below the onset of fluidization 67 ).We posited that the shape of the angular dependence map α z,x was preserved under the air flow.The decrease of granular resistance was then modeled by a single scaling factor (like ζ in the generic stress calculation).The average forward speed vx vs. v air in modified terradynamic simulations (Fig. 12(f), solid lines) and experiments (Fig. 12(f), symbols) agreed well for low air speed v air ≈ 0. Though in experiments vx decreased faster with increasing v air than in RFT (especially at the highest frequency ω = 8 rad/s), the model qualitatively captured the experimentally observed failure of the locomotor at large v air .

IV. CONCLUSION AND FUTURE DIRECTIONS
In this review, we have discussed examples from our group of the efficacy of Resistive Force Theory in modeling aspects of locomotion on and within granular media.The studies above have focused on application of RFT to dry, homogeneous, level, pristine GM (slightly polydispersed in particle sizes).We expect that the RFT approach should prove useful in other dry granular settings and below we give a brief summary of a few situations and regimes.[70] This article is copyrighted as indicated in the article.Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions.

A. Transient material response
The current formulation of granular RFT assumes that the forces acting on elements are described by the steady state drag forces.This approximation is valid at low Reynolds numbers, as upon a sudden movement of a body to a constant speed, the accompanying flow fields are developed rapidly (the so-called "instantaneity" assumption 71 ) and are thus time independent for such an interaction.However, in granular media, even for low constant intrusion speed, the material can change state over timescales relevant to the intrusion event.The ∼ 30% over-prediction of sand-swimming speeds in RFT compared to the DEM in Ref. 32 was attributed to such a transient effect.To understand this more clearly, experiment and simulation were designed and conducted to study transient effects of granular response to reciprocal motions. 29Oscillation of a rod in the horizontal plane revealed discrepancies between the simulation and the empirical force relations.In RFT, the oscillatory motion generates a drag force which is described by a square wave.However, in the DEM simulation, from the instant the rod direction is reversed (or started from rest), a significant portion of the cycle (≈1/5 of the amplitude) was required for the force to increase to its steady state value.A similar transient weakening effect was observed and studied in a cyclically sheared granular medium 72 ; the underlying physics of this effect may be related to changes in geometric structure among particles due to preparation methods or previous disturbances. 72,73

B. Surface effects: Disturbed ground and inclines
Several aspects of interaction with surfaces were not incorporated into the granular RFT.For example, unlike fluids, soft materials like sand can display memory effects: for example, GM can form sand craters upon impact 74,75 which will not be removed (as in a fluid) by the spontaneous fluctuations of particles or by further flow (since the material solidifies post-impact).We have discovered that such effects in which surface properties change due to disturbances can play significant roles in locomotor performance.In studies of the locomotion of a flipper-based robot, drag and penetration forces were found to decrease substantially if an intrusion occurred close to previously disturbed ground. 76Such memory effects of "disturbed ground" were posited to be the cause of locomotor failure 76,77 during surface movement.
Movement on sandy inclines presents another interesting challenge for RFT, in particular since (during uphill movement) regions of the locomotor can come into contact with flowing material and/or the region of disturbed material can be large.Recently, to understand the movement of sidewinder snakes on inclines, we performed the first measurements of the drag force on a plate moving downwards on a tilted granular surface. 78The component of gravity along the surface of the medium weakens the drag force along that direction.At a tilting angle of 20 • the drag force was reduced by approximately 50%.This phenomena is quite different from penetration into a flat granular surface at the same angle, where grain flow happens at both the bottom and top of the intruder.Currently less is understood about how the granular penetration force is affected by this additional anisotropy, so it is unclear how to incorporate this effect into an RFT formulation.

C. Inertial effects in the medium
In granular media, RFT has been tested in the friction dominated regime where the penetration force is independent of the speed of the intruder.When the intrusion speed is large, the inertia of grains contribute significantly to the forces, 44,52,79 and we have discovered that lightweight locomotors on GM can utilize such forces 80 to maintain high performance.Since the force fields are not necessarily local, it is unclear if the RFT superposition approximation will be useful.At the minimum, speed and acceleration effects must be incorporated into the force relations.

101308- 4 TFIG. 3 .
FIG.3.Multi-particle Discrete Element Method simulation (DEM) reveals how material flows around a swimming sandfish simulation.(a) A DEM simulation of the sandfish lizard S. scincus in a medium of approximately 3 × 10 5 approximately 3 mm diamater "glass" spheres.The snapshots of the animal's body are rendered in red.(b) In DEM simulation, when particles (of radius r) with velocities v 1 and v 2 collide with each other, the inter-particle contact forces F n,s are computed from a "virtual" overlap δ which represents the elastic deformation at the contact area.(c) Top view of the subsurface swimming in DEM (particles above the plane of interest have been rendered transparent).The movement of the animal's body (indicated by the white region) mobilizes the particles around it.Colored stripes are used for visualization.Colors in the stripe correspond to the initial y positions (decrease from red to blue) of particles before the medium is disturbed.Image created by Yang Ding.

FIG. 5 .FIG. 6 .
FIG. 5. Resistive Force Theory (RFT) relations in the horizontal plane.RFT force relations were determined from rod drag experiments.(a) Experimental setup.The substrate was fluidized (and then compacted) to a controlled volume fraction before each drag experiment was executed; air flow is off during drag.(b) The measured drag forces (in 0.3 mm glass particles) plotted as functions of time for the attack angle ψ = 60 • (defined at the inset of (c)) at v = 10 cm/s.The rod-particle friction μ ≈ 0.15.Dotted lines correspond to F and solid lines to F ⊥ .The average was taken in the shadowed region when the forces had reached asymptotic values and the rod was still far from the boundary.Adapted with permission from R. D. Maladen, Y. Ding, C. Li, and D. I. Goldman, Science 325, 314-8 (2009).Copyright 2009 by the Authors.Reprinted with permission.LP/CP stands for loosely packed/closely packed granular medium (volume fractions of approximately 0.58 and 0.62, respectively.(c) Averaged F ⊥, on the cylinder wall vs. angle of attack ψ.F ⊥ is enhanced in GM compared with a low Re fluid.Adapted with permission from Y. Ding, S. S. Sharpe, A. Masse, and D. I. Goldman, PLoS Comput.Biol.8, e1002810 (2012).Copyright 2012 by the Authors.Reprinted with permission.

FIG. 7 .
FIG. 7. Resistive Force Theory (RFT) relations in vertical plane for different granular materials.(a) GM tested in the study.(Top) regular images.(Bottom) microscope images.The length of each scale bar is 1 mm.(b) Coefficients from discrete Fourier transforms can be scaled and approximated by a generic coefficient curve M 0 (thick black line with stars).Convention of the coefficients is given in the supplementary materials of Ref. 31.LP/CP stands for loosely/closely packed media, respectively.(c) and (d) Generic stress profiles (per unit depth) constructed from M 0 .Color indicates magnitude of stress.Adapted with permission from C. Li, T. Zhang, and D. I. Goldman, Science 339, 1408-12 (2013).Copyright 2013 by the Authors.Reprinted with permission.
FIG. 8. RFT force predictions compared with experiment in the vertical plane during drag (a) and rotational penetration (b).(a) Comparison of DEM calculated normal σ and tangent τ stresses on a circular cylinder (translating left to right) with data from plate drag experiments.Here the orientation angle β is defined such that it ranges from [0, π ].Reprinted with permission from Y. Ding, N. Gravish, and D. I. Goldman, Phys.Rev. Lett.106, 028001 (2011).Copyright 2011 American Physical Society.(b) Rotational penetration forces in experiment and RFT versus the rotation angle θ .The radius of the c-leg R = 3.8 cm.The angular speed ω = 2π rad/s; doubling ω does not affect the profiles of forces.Solid lines denote experimental measurements, dashed lines the RFT calculations.Adapted with permission from C. Li, T. Zhang, and D. I. Goldman, Science 339, 1408-12 (2013).Copyright 2013 by the Authors.Reprinted with permission.

FIG. 9 .
FIG. 9. Application of horizontal plane RFT to biomechanics.(a) Subsurface swimming pattern of S. scincus.Color indicates the "snapshot" time.(b) Movement of C. occipitalis in a substrate of 0.3 mm glass particles.(c) RFT predicated undulation efficiency (the forward progress per cycle normalized by wavelength) versus mean relative body curvature κλ s .Cross marks refer to animal swimming data.LP/CP stands for loosely packed/closely packed granular medium (volume fractions of approximately 0.58 and 0.62, respectively.ξ denotes the undultion number, the average number of waves on the body.

101308- 11 T
FIG. 10.RFT in neuromechanics.(a) RFT model of Neuromechanical Phase Lags (NPL).Magenta arrows represent body velocities and green arrows represent RFT-predicted forces from the media.Negative (−) muscle torque τ muscle indicates muscle contraction on the right side of the body (red thick line) where electromyography signals were measured.A positive rate of change in curvature κ(+) indicates that the body segment is currently bending towards the left.ψ has the same convention as in Fig. 5(c).(b) RFT-predicted temporal profile of κ (blue curve) and torque (green curve) at the normalized horizontal position x = 0.71π (black circular marker in (a)).The time span of muscle activation on the right side of body is indicated by the dashed vertical lines.The shaded region indicates κ(−).(c) Predicted initiation and cessation of muscle activation from RFT (green lines) along the body compared with electromyography measurements (black points with error bars).Dashed lines are a visual guide.Time range 0 − 2π is the normalized period of undulation.Figure modified from Ref. 41.

101308- 12 T
FIG. 11.Application of RFT to geometric mechanics.(a) The three-link swimmer resting on a medium of 6 mm plastic particles.(b) Configuration of the three-link swimmer.The body frame x b − y b is attached to the center of mass of the robot and its orientation θ is the weighted average of the links.(c) A x the x-row of the local connection matrix (see details in the text) and (d) The Curvature Constraint Function (see text for detail) of A x .(e) Body displacement from RFT-enabled CCF model (red curve) compared with experiments (dots with error bars) and DEM simulations (open circles), as functions of the stroke amplitude (the radius of the circle in the shape space).The undulating frequency of the swimmer is set at 0.5 Hz (thus in the non-inertial regime).The dashed black curve in the inset represents a "butterfly" gait.For the butterfly gait, in the main panel, the dashed line corresponds to the CCF prediction and the dotted line to DEM calculation.Adapted with permission from R. Hatton, Y. Ding, H. Choset, and D. I. Goldman, Phys.Rev. Lett.110, 078101 (2013).Copyright 2013 by American Physical Society.

101308- 13 T
FIG.12.RFT model for locomotion on a granular surface.(a) Xplorer on 1 mm diameter poppy seeds.Photo credit: Chen Li.(b) Forward speed v x (t) of Xplorer in RFT model (dashed line) compared with experiment (solid line).31(c) 3D-printed robot legs with different curvatures.All legs have the same projected length 2R.The curvature 1/r is considered "positive" if the leg is bending towards the left; negative towards the right.The rotation direction of the leg is given by ω.(d) The average speed vx of Xplorer as a function of the relative curvatures R/r of the legs.Dashed lines are simulation results from RFT. Markers correspond to data from experiment.Color indicates leg rotation frequency (f = 1, 2, 3, 4 Hz, from blue to red).The diagram of the leg being used is shown (at the position corresponding to the data markers) to facilitate visualization.Adapted with permission from C. Li, T. Zhang, and D. I. Goldman, Science 339, 1408-12 (2013).Copyright 2013 by the Authors.Reprinted with permission.(e) SandBot with cylindrical wheels.Photo credit: Feifei Qian.(f) Average forward speed vx of SandBot on poppy seeds as a function of air flow speed for different leg rotation frequencies (ω = 2, 4, 6, 8 rad/s, from blue to red).Solid curves correspond to RFT model predictions.Markers correspond to experiment results.