## Abstract

We integrate biological experiment, empirical theory, numerical simulation and a physical model to reveal principles of undulatory locomotion in granular media. High-speed X-ray imaging of the sandfish lizard, *Scincus scincus*, in 3 mm glass particles shows that it swims within the medium without using its limbs by propagating a single-period travelling sinusoidal wave down its body, resulting in a wave efficiency, *η*, the ratio of its average forward speed to the wave speed, of approximately 0.5. A resistive force theory (RFT) that balances granular thrust and drag forces along the body predicts *η* close to the observed value. We test this prediction against two other more detailed modelling approaches: a numerical model of the sandfish coupled to a discrete particle simulation of the granular medium, and an undulatory robot that swims within granular media. Using these models and analytical solutions of the RFT, we vary the ratio of undulation amplitude to wavelength (*A*/*λ*) and demonstrate an optimal condition for sand-swimming, which for a given *A* results from the competition between *η* and *λ*. The RFT, in agreement with the simulated and physical models, predicts that for a single-period sinusoidal wave, maximal speed occurs for *A*/*λ* ≈ 0.2, the same kinematics used by the sandfish.

## 1. Introduction

Discovering principles of locomotion for organisms in natural environments requires integrating experimental visualization [1–4], theoretical [5–7], computational [8–10] and robotic [11–14] approaches to generate and test hypotheses of interaction between a locomotor and its surroundings [9,15–17]. These approaches have produced an understanding of locomotion biomechanics in a range of terrestrial, aquatic and aerial environments [18–21]. Comparable analysis of locomotion within yielding substrates like sand, soil and debris that display both solid and fluid-like behaviour is less developed [9,22–25]. This is, in part, because visualizing motion in opaque materials is challenging, and validated force models on and within these substrates at the level of those that exist for fluids (see [19,26,27]) do not exist yet.

Attempts to understand burrowing and movement within granular media like desert sand directly confront these issues. Subsurface locomotor behaviours are relevant to a diversity of organisms, such as scorpions, snakes and lizards, that move within the sand to escape heat and predators, and to hunt [28–30]. To answer questions about ecological adaptations [28,31–33] and evolution [34–36] in flowing substrates requires an understanding of how propulsion is generated and how both media interactions and animal kinematics enhance or limit locomotor performance.

Recently, we used high-speed X-ray imaging to visualize the locomotion of a small desert lizard, the sandfish, *Scincus scincus*, within approximately spherical 0.3 mm diameter glass particles, a granular medium representative of desert sands [37]. We found that once it buries rapidly (within one-half of a second) into the sand, the animal swims subsurface at up to 2 bl s^{−1} (bl, body-length is defined as the distance from the snout tip to the base of the tail) by propagating a sinusoidal wave down its body (anterior to posterior) without using its limbs [22]. The animal monotonically increases speed by increasing its oscillation frequency while using the same undulatory kinematics (a single-period sinusoidal travelling wave with an amplitude to wavelength ratio of approximately 0.2) even within media prepared to different initial yield strengths (controlled via volume fraction or ‘media compaction’).

Developing models to describe subsurface locomotion requires understanding resistive forces on intruders and how intruders change the local properties of the media [38]. As a broadly applicable theory of intrusion into granular media does not yet exist in the regime relevant to sand-swimming, in the study of Maladen *et al.* [22] an empirical resistive force theory (RFT) was developed to predict the wave efficiency *η*—the ratio of the forward swimming speed to the wave speed. The granular RFT was inspired by the theory used to predict swimming speeds of micro-organisms [39] in fluids at low Reynolds numbers (*Re*), where inertial effects are small, and balances empirically measured estimates of drag and thrust forces on elements of the animal. The RFT demonstrated that body undulations are sufficient to propel the sandfish forward, and its prediction of *η* was close to the observed *η*.

The RFT predictions are based on a simplified model of the animal kinematics and an empirical model of the granular resistive forces on the body. In this paper, we describe the development of two, more detailed and flexible models of undulatory sand-swimming that reproduce the observed biological sand-swimming and exhibit performance in accord with predictions of the RFT. In the first, a computer simulation approach, we couple a numerical model of the sandfish with a multi-particle discrete element method (DEM) [40] simulation of the granular medium. In the second approach, we measure the locomotor performance of a physical model (a robot). We use all three models to study how the wave kinematics of the rapidly escaping sandfish maximize forward swimming speed within the granular media. More generally, these complementary modelling approaches will allow us to explore and answer questions about undulatory locomotion in complex, yielding media.

## 2 Biological experiments

We use an air-fluidized bed to control the initial preparation of the material. The packing state of the substrate is parametrized by the volume fraction *ϕ*, the ratio of solid to occupied volume [41] of the granular medium. The fluidized bed is a box of area 21.5 cm × 18 cm with a porous base and is filled with granular material 10 cm deep. A continuous flow of air through the base expands the granular medium in the bed to a low *ϕ* state; subsequent pulses of air increase compaction allowing precise control of the initial volume fraction of the bed. All air flow through the bed is turned off before the animal is released onto the medium. We use high-speed X-ray imaging to record the movement of the sandfish, *Scincus scincus*, as it swims subsurface (for further experimental details, see [22]).

In the work of Maladen *et al.* [22], we studied sandfish locomotion within 0.27 ± 0.04 mm diameter approximately spherical glass particles, close in size to the sand found in the animal's natural environment [42]. Here, we study swimming in larger spherical glass particles (diameter: 3.2 ± 0.2 mm, referred to here as ‘3 mm particles’) to reduce the time required to simulate the granular medium (figure 1*a*,*b*). High-speed X-ray imaging experiments reveal that the sandfish (five animals; mass = 16.9 ± 4.0 g; body-length (bl) *L* = 8.8 ± 3.3 cm) swims subsurface in the 3 mm particles with performance comparable to that in the 0.3 mm particles. The sandfish are tested in loosely (*ϕ* = 0.60) and closely (*ϕ* = 0.62) packed media preparations and take 1.2 ± 0.6 s to complete the burial process (time from entry until tail is not visible) independent of media preparation. Interestingly, although swimming performance is comparable in the two bead sizes, the burial time in the 3 mm particles is significantly longer than in the 0.3 mm particles. Once subsurface, the animal places its limbs along its sides and advances by executing an undulatory motion (see electronic supplementary material, movie S1). We quantify the subsurface motion by tracking the mid-line of the animal.

As in the study of Maladen *et al.* [22], all phases of the undulatory motion of the animal midline are well-fit by a posteriorly travelling single-period sinusoidal wave (*R*^{2} > 0.9)
2.1with *y* the displacement from the mid-line of a straight animal, *A* the amplitude, *λ* the wavelength, *f* the wave frequency, *v*_{w} = *f* *λ* the wave speed, *t* the time and *x* the distance along a line joining the endpoints of the animal and parallel to the direction of motion. The *x* and *y* coordinates are set in the laboratory frame. The number of periods along the body of the swimming animal, *ξ*, is approximately 1.0 in experiment. The spatial characteristics of the wave, *A* and *λ*, do not vary significantly within a run or between runs and their ratio is conserved at approximately *A*/*λ* = 0.25 ± 0.07. *A* (average 0.18 ± 0.03 bl, *p* = 0.26) and *λ* (average 0.7 ± 0.1 bl, *p* = 0.13) are not significantly different for loosely and closely packed media preparations. A *p*-value greater than 0.05 indicates that no statistically different relation exists between the measured values for the different conditions tested [43]. As within 0.3 mm particles [22], the animal increases its forward velocity by increasing its undulation frequency (figure 1*d*).

A measure of performance of undulatory locomotion in deformable media is the wave efficiency (also referred to as the ‘amount of slip’; [44]), *η*, the ratio between the average forward velocity of the animal *v*_{x} and the velocity of the wave travelling down its body *v*_{w}. Typical values of *η* for small organisms (such as nematodes) for whom inertial effects are negligible are 0.22 ± 0.02, while swimming in fluids at low *Re* [39] and 0.8–0.9 while creeping along an air–solid interface [44–46]. Because *A*/*λ* for sandfish is independent of *f*, *η* is the slope of the *v*_{x}/*λ* versus *f* curve (figure 1*d*). For 3 mm particles, *η* = 0.54 ± 0.13 (figure 1*d*,*e*), similar to the value of 0.52 ± 0.02 for locomotion in 0.3 mm particles [22]. For both particle sizes, *η* is independent of media preparation.

## 3. Resistive force theory for granular media

We use a previously developed RFT [22] for granular media to predict *η*. The RFT assumes that the animal swims in a horizontal plane, although it actually moves into the medium at an angle relative to the horizontal of approximately 22° in 0.3 mm particles [22] and approximately 29° in the 3 mm particles. We discuss the validity of this assumption in §7. In RFT models, the body of the locomotor is partitioned into infinitesimal elements along its length. When moving relative to the medium, each element is acted on by resistive forces (figure 2*a*) that can be decomposed into thrust and drag components. Resolving these forces into perpendicular (*F*_{⊥}) and parallel (*F*_{‖}) components, the net forward force on an element is
3.1where *θ* is the angle between the direction of the average velocity of the organism and the instantaneous orientation of the infinitesimal element. *θ* increases with the oscillation amplitude and can be obtained by differentiating equation (2.1):
3.2

For granular media, the force on an element is well characterized as a function of only the direction of the velocity relative to its orientation (*ψ*) since in the speed regime relevant to the animal previous studies, Maladen *et al.* [22], Albert *et al.* [47] and Wieghardt [48], showed that granular force is independent of speed. Determining the *y*-component of the velocity from equation (2.1)
3.3allows us to write the angle *ψ* as
3.4The net forward force on the entire body is calculated by integrating the force on each element over the length of the body (snout to tail tip). Since the body propagates a single period wave (*ξ* = 1) and the force on a segment is assumed to be localized, propagation of the wave in essence moves a segment from one end of the body to the other end with the orientation, velocity and force on that segment unchanged. Therefore, the net force, the sum of the forces on all segments, is time-invariant. The drag on the forward surface of the head (element end cap), , is averaged over a cycle. To estimate the average forward speed, the forces on the body and head are summed and set to zero. Assuming that the forces on an infinitesimal body element are proportional to the area of the sagittal cross section (*b**δ**s*, where *b* is the height of the element, and is the arc length) and are functions of *ψ*, the total time-averaged force on the sandfish can be expressed as:
3.5where *P*_{⊥} and *P*_{‖} are stresses perpendicular and parallel to the axis of each element. By assuming a constant average velocity (net forward force along the body equal to zero), the speed of the organism and thus *η* can be determined. Motion and force only in the forward direction are considered.

Unlike in fluids, in an arbitrary granular medium, no validated theory exists to calculate the resistive forces on an element as a function of its angle relative to its velocity direction (*ψ*). Previously, Maladen *et al.* [22] obtained these force laws empirically at different *ψ* in 0.3 mm glass particles—the same granular medium in which the animal was tested in [22]. To estimate thrust and drag on the body and head of the organism (figure 2*a*), empirical fitting functions were used to separate the forces on the rod into terms describing the surface along the length (side) and end-cap of the rod.

Here, to avoid the approximations introduced by using empirical laws to resolve side and end-cap forces, we use a validated simulation of the granular medium to directly measure the surface forces acting on the dragged rod. The granular medium is simulated using a three-dimensional soft sphere discrete element methods (DEM) code [40]. To compute particle–particle and particle–intruder interaction forces, we calculate the normal force [49], *F*_{n}, and the tangential Coulomb friction force, *F*_{s} (see figure 3*c*) at each contact using
3.6where *δ* is the virtual overlap between contacting particles, (see figure 3*c*) is the normal component of relative velocity with *n̂* the unit vector along the line connecting the particles centers, and *k* and *G*_{n} represent the hardness and viscoelastic constants, respectively. *μ* refers to the particle–particle (*μ*_{pp}) or body–particle (*μ*_{bp}) friction coefficients depending on which objects are in contact. For specific values, see electronic supplementary material, table S1. The simulated medium (a 50 : 50 bi-disperse mixture of 3.4 and 3.0 mm glass particles, with density 2.47 g cm^{−3}, which we will refer to as 3 mm particles) is validated by comparing the forces on a cylindrical stainless steel rod (diameter = 1.6 cm, length = 4 cm) dragged horizontally through it with those from experiment (see electronic supplementary material, figure S1). Particle–particle restitution and friction coefficients are experimentally measured, while the hardness is selected such that particle overlap (*δ*) is less than 0.4 per cent of particle radius. The container holds a granular volume 25 × 24 × 17 cm^{3} in extent. We obtain a match between the force measured as the rod is dragged through the granular media in experiment and simulation as a function of time (see electronic supplementary material, figure S1*a*). The mean deviation between experimentally measured average drag force and simulated average drag force as a function of angle is less than 5 per cent (see electronic supplementary material, figure S1*b*).

The simulation allows direct measurement of the forces on the side and end caps of the rod separately from which we obtain *F*_{⊥} and *F*_{‖} (figure 2*b*,*c*). To estimate the effect of head drag, only forces on the leading end cap are considered. We note that *F*_{⊥} and *F*_{‖} are measured on a square cross-section rod (width, height *b* = 1.6 cm, length *d* = 4 cm), as we approximate the cross section of the sandfish as a square. The fitting functions used to obtain continuous functions from the discrete measurements of *F*_{⊥} and *F*_{‖} are reported in the electronic supplementary material. *P*_{⊥} and *P*_{‖} are obtained using the relation *P*_{⊥} = *F*_{⊥}/*bd* and *P*_{‖} = *F*_{‖}/*bd*. We measure the forces within both loosely and closely packed media (see electronic supplementary material, figure S3).

To determine *η* using the RFT, we assume that the animal swims at constant speed, which implies that the net force integrated over the body is zero (i.e. in equation (3.5)). Substituting *F*_{⊥} and *F*_{‖} from the simulation and using the spatial characteristics of the animal, *A*/*λ* = 0.22 (close to the average *A*/*λ* measured in 0.3 [22] and 3 mm particles) and *ξ* = 1, we numerically solve equation (3.5) using an unconstrained nonlinear optimization, which finds the minimum of a scalar function of several variables starting at an initial estimate. With a flat head (maximum drag) this procedure yields *η* = 0.66 in loosely packed and *η* = 0.59 in closely packed media; both values are larger than that measured for the animal (*η* = 0.54 ± 0.13; figure 1*d*,*e*). We discuss this difference in §7.

## 4. Numerical and robotic modelling approaches

To better model the animal and test the predictions of the RFT, we develop two complementary models of the sandfish in which parameters can be readily varied: a numerical simulation that combines a model of the animal with the validated granular medium model, and a sandfish inspired physical robot that is tested in a granular medium. These tools also allow us to investigate the effects of varying parameters like *A*, *λ* and *ξ*, variables which we have no control over in the biological experiment (§2). First we describe the development of the modelling approaches and compare their predictions of *η* to the RFT and the animal experiment. Then using each approach, we vary kinematic parameters to systematically test the RFT predictions as well as show an optimality condition for undulatory swimming in granular media.

### 4.1. Numerical sandfish simulation

The numerical simulation of the sandfish uses the commercial software package Working Model 2D (Design Simulation Technologies). The model is divided into 50 motor-actuated segments along its length (figure 3*a*). The angle of each motor is specified as a function of time such that an approximate sinusoidal wave travels posteriorly from head to tail (figure 3*b*). The angle (*β*_{i}) between segments *i* and *i* + 1 varies as
4.1where *x*_{i} is the arc length measured from the tail tip to the *i*^{th} segment and *l* is the body length. To compare performance of the numerical simulation to the RFT, the morphology of the simulated animal is approximated as a square cross-section tube (figure 3*a*, inset). The model does not incorporate limbs as the sandfish places its limbs along its sides during subsurface swimming (figure 1*b*). No tapering in height of the body along its length is considered, as doing so results in the model rising as it moves forward owing to drag-induced granular lift. This lift results from the vertical component of the normal force on an inclined surface dragged through a granular medium and is discussed in Ding *et al.* [50] and Maladen *et al.* [51].

We perform simulations for a square tube body shape, and to more closely match the animal we perform simulations using a model tapered in the coronal plane (figure 3*a*) with width varying linearly from the snout tip to 1/6 bl, and from 3/5 bl to the tail tip. Model properties such as animal dimensions and density are taken from biological measurements (refer to electronic supplementary material, table S1, for values). The segment height and maximum segment width are both denoted by *b* for the model. The body-particle friction (*μ*_{bp} = 0.27) was determined by measuring the angle at which an anaesthetized sandfish slid snout-first down a monolayer of 3 mm glass particles glued to a flat plate [16]. For simplicity, we use the same normal force parameters for both the particle–particle and body–particle interactions.

The multi-segment model of the animal is combined with the previously described experimentally validated DEM simulation of the granular medium (in a container with comparable dimensions to those in experiment, see electronic supplemental material, table S1) such that Working Model integrates the equations of motion of the coupled links that represent the animal and the DEM calculates the resultant force from both the particle–particle and body–particle interactions. At each time-step, the net force from particles contacting each segment is passed to Working Model, and velocity and position information for each segment is transferred back to the DEM simulation. The animal model kinematics are constrained to prevent rotation of the axis of the travelling wave out of the horizontal plane of motion (roll).

With a flat head and square tubular body and using *A*/*λ* = 0.22 and *ξ* = 1 (from animal experiments), the simulated sandfish moves in loosely packed 3 mm particles with *η* = 0.45, while for the same shape, the RFT predicts *η* = 0.66. With the tapered sandfish body in simulation, the forward velocity of the simulation increases linearly with oscillation frequency, as in the animal experiments and the flat head and square tubular body simulation, and the slope of the relationship, *η* = 0.57 ± 0.01 (figure 1*d*,*e*) is close to *η* = 0.54 + 0.13 observed for the animal. Accounting for tapering along the sandfish body in the RFT is complicated as the taper at the head and tail can become leading and trailing surfaces at different instances. This violates an RFT assumption that the net force can be decomposed into time-invariant body and head terms. We estimated the amount of head drag reduction owing to tapering as approximately 30 per cent (in accord with drag-reduction measurements [52]) by comparing the integrated drag force on the flat head to the tapered head in simulation. The RFT for the reduced head drag over-estimates the wave efficiency, predicting *η* = 0.75. We discuss potential reasons for the differences in *η* between simulation and RFT in §7.

### 4.2. Sandfish-inspired robot

To test the RFT and numerical simulation predictions in a real world environment, we built a physical model, a robot. Robots complement the theoretical and numerical modelling approaches because physical laws that govern the interaction with the granular medium need not be simulated.

The basic mechanical design of our robot is adapted from previously developed snake robots [53], which consist of modules (motors) with a single joint that permits angular excursions in the body plane and that are connected via identical links to form the body. Our design employs six standard size (4 × 3 × 3.7 cm^{3}) servomotors (Hitec, HSR 5980SG) and a passive segment (the head) with the same weight and form factor as a motor with its attached connectors for a total of seven segments. Each servomotor is powered in parallel from a 7.4 V, 30 A power supply. The pulse width-based control signal for each motor is generated using LabVIEW as a single analogue waveform initially determined from equation (4.2), output from a PCI-card (NI-6230) and transmitted via a multiplexer to each motor as control input.

Preliminary tests of the robot revealed that the kinematics prescribed by the sinusoidal travelling wave (equation (4.1)) are accurate only for low amplitudes because of the finite segment length associated with a finite number of segments (fixed total length). To obtain a larger range of amplitudes, we implemented an open loop controller that modulates the angle between adjacent segments using
4.2where *β*_{i}(*t*) is the angle between the *i* and *i* + 1 motor at time *t*, *β*_{0} is the angular amplitude, *ξ* is the number of wavelengths along the body (period) and *N* is the number of motors. We experimentally verified that the motors are able to achieve the prescribed joint excursions by observing the kinematics of the robot subsurface using X-ray imaging.

To reduce the motor torque requirements, we use 5.87 ± 0.06 mm plastic particles with particle density = 1.03 ± 0.04 g cm^{−3} as our granular medium (see electronic supplementary material, table S1, for details); we refer to these particles as ‘6 mm particles’. The granular bed is 110 cm × 40 cm and filled 14 cm deep. To prevent particles from intruding between motor segments and jamming the device, we encase the robot in a double-layer skin; a lycra spandex sleeve with a single seam (located at the top of the device) is slid over a snugly fitting elastic latex sleeve that prevents particles from getting between the motors during phases of the motion when the lycra sleeve is loose. The outer layer reduces the wear on the thin inner layer, allowing greater than 100 subsurface trials without replacement of the inner layer. The control and power wire bundle is run along the top of the motors (under the inner latex sleeve) to a mast on the last segment and is tethered above the container.

Overhead video (100 fps) is collected for each experiment. To track the subsurface position of the robot, a lightweight mast is attached to the first and last motor segment with a marker (visible above ground) (figure 4*a*,*e*). The kinematics of the subsurface motion of the robot are also verified using X-ray video (figure 4*b*) for a representative condition. For each test, the robot is submerged 4 cm (from its top) into the medium and the surface is levelled and smoothed. At least 1.5 cycles of motion are collected for each trial (see electronic supplementary material, movie S3).

Like the sandfish lizard and the RFT and numerical simulation models, the forward velocity of the robot monotonically increases with increasing *f* (figure 5*a*) for *A*/*λ* = 0.2 and a single period wave. However, for the robot, *η* = 0.34 ± 0.03, significantly below the values measured for the animal in experiment, and predicted by the RFT and animal numerical simulation.

We speculate that the *η* of the robot is less than that of the sandfish in biological experiment, RFT and numerical models because the small number of segments do not allow the robot to achieve a smooth sinusoidal profile. To investigate segment number effects, rather than increasing the number of motors in experiment, we adapt the numerical sandfish simulation to match the robot and its granular medium in physical dimension, mass and density (figure 4*c*). We do this in simulation as increasing the number of motors experimentally is challenging: doing so would require an increased number of smaller motors with high torque requirements to maintain the same total length as the seven-segment robot.

The simulated robot, like the actual device, consists of seven segments and the angle between adjacent segments is modulated using equation (4.2). The robot model is tested in a validated granular medium consisting of a 50 : 50 bi-disperse mixture of 5.81 and 5.93 diameter particles; we will refer to this mixture as ‘6 mm plastic particles’. To validate the simulated medium and obtain the values of *μ*_{pp}, *k* and *G*_{n} given above, we dropped an aluminium ball (diameter 6.35 cm and mass 385 g) into the plastic particles with varying impact velocity (0.5–3 m s^{−1}) in both experiment and simulation and set grain interaction parameters to best match the measured and simulated penetration force during the impact collision as a function of time (see electronic supplementary material, figure S2 and table S1). With parameters determined from impact at *v* = 1.4 m s^{−1}, the force profile fits well at other impact velocities. In additional experiments, we directly measured the friction (particle–particle and particle–lycra skin) and restitution (particle–particle) coefficients for the plastic particles and found them to be within 5 and 10 per cent of the fitted values, respectively. The robot simulation results quantitatively match the robot experimental results (figure 5*a*) with velocity linear in *f* and *η* = 0.36 ± 0.02 (see electronic supplementary material, movie S4).

Increasing *N* (for a fixed length device) in the robot simulation causes the robot to advance more rapidly and with greater *η* until *N* ≈ 15, above which *η* remains constant at 0.54 (figure 5*b*), resulting in a performance comparable with that of the animal. *η* obtained from the RFT for a robot with a continuous profile (i.e. *N* = ∞) swimming in 6 mm plastic particles is nearly the same (black horizontal line, figure 5*b*). For the empirical force laws used in the RFT to predict *η* refer to electronic supplementary material, figure S4. For the smooth robot, the animal in 0.3 mm (loosely and closely packed media) [22] and 3 mm particles (§2), and in model predictions (§§3 and 4) in simulated 3 mm glass particles, the values of *η* are all close to 0.5. This suggests that subsurface locomotion in granular media is largely independent of media properties like particle size and density [51].

We postulate that *η* increases with increasing *N* because a smoother body profile facilitates media flow and leads to decreased drag and/or a decreased variation in the spatial form (*A*/*λ*) from the prescribed sinusoidal target (as much as 30% for *A*/*λ* = 0.2 with *N* = 7). Our finding that for *N* > 15, *η* saturates at a value close to that observed in the animal suggests that near-maximal performance occurs when the spatial form of the robot matches well a travelling sinusoidal wave.

## 5. Variation of the sinusoidal kinematics of sand-swimming

To this point, we have limited our investigation of sand-swimming to parameters used by the sandfish, i.e. a sinusoidal travelling wave with fixed kinematics (*A*/*λ* ≈ 0.2 and *ξ* = 1). We now use each of our three models to investigate how swimming depends on variation of the spatial form of the sinusoidal wave. This allows us to systematically test our models and to advance an argument for why the rapidly escaping animal uses only a limited range of *A*, *λ* and *ξ*.

We first discuss the RFT prediction of *η* for a fixed length sandfish model at fixed period (*ξ* = 1) with varying wave amplitude such that 0 < *A*/*λ* < 0.6 (figure 6*a*). Solving equation (3.5) numerically reveals that *η* increases monotonically with increasing *A*/*λ*: slowly for small *A*/*λ*, rapidly for intermediate *A*/*λ* and slowly towards a plateau for large *A*/*λ*.

We compare the functional form of *η* versus *A*/*λ* from the RFT to the results from the numerical model with both constant and tapered cross section (figure 6*a*) and find that while the shape of the curves are similar, the RFT predictions (for the square cross section and for the reduced head drag) systematically overestimate *η*. However, scaling the net force (equation (3.1)) on each element (excluding the head) by a factor of 0.5 in the RFT prediction results in a curve that closely matches the non-tapered simulation (figure 6*a*). We return to this point in §7. For all *A*/*λ*, the simulated sandfish with tapered body progresses with a higher *η* than for the non-tapered sandfish mainly owing to the reduction in head drag. The dependence of *η* on *A*/*λ* from both robot experiment and simulation (figure 6*b*) qualitatively matches that predicted by the RFT and numerical sandfish simulation. As explained in §4.2, the reduced performance of the robot compared with the RFT and numerical sandfish simulation is due to the smoother body profile of the animal model (because of the larger number of segments) when compared with the robot.

### 5.1. Analytical approximate solution for resistive force theory

Since the shapes of the *η* versus *A*/*λ* curves are similar over a wide range of experimental and model conditions, and since performance limits on sand-swimming are determined by the dependence of *η* on *A*/*λ*, it is important to understand how parameters in the model affect this functional relationship. Therefore, we gain insight into this function by developing analytical solutions of the RFT in three regions: low *A*/*λ*, where *η* increases with positive curvature, intermediate *A*/*λ*, where *η* increases significantly and high *A*/*λ*, where *η* is largely independent of *A*/*λ*.

As equation (3.1) shows, thrust is determined by the orientation *θ* of the element, the maximum value of which increases with *A*/*λ*, and the magnitude of the normal force. The increase in the projected area of the element perpendicular to the forward motion of the animal (which increases thrust) explains the increase in *η* as a function of *A*/*λ*. Understanding the form of *η* versus *A*/*λ* in the three regions, however, requires examination of the measured force laws (figure 2*b*,*c*). Because of the non-trivial dependence of the force on *ψ*, we divide the *η* versus *A*/*λ* relationship into three adjacent regions (figure 6*a*) in each of which we approximate the force as either constant or linearly dependent on *ψ*. Below, we give a summary of the calculations; see electronic supplementary material for full details.

The order of the regions in figures 2 and 6 is reversed because *A*/*λ* and *ψ* are inversely related (inset, figure 6*a*). As solving the RFT numerically without head drag does not qualitatively affect the dependence of *η* on *A*/*λ*, we neglect head drag to simplify the analysis.

Region 1 (small *A*/*λ*). For small and increasing oscillation amplitude (and similarly for *A*/*λ* and *θ*), *η* increases from zero with an increasing rate. In this region as *ψ* decreases, the displacement of an element of the body of the animal remains nearly perpendicular to its forward velocity *ψ* ≈ *π*/2 (figure 2*b*,*c*) for a majority of the cycle resulting in *F*_{⊥} remaining nearly constant (at its maximum) while *F*_{‖} increases. Consequently, the stresses perpendicular and parallel to an element can be approximated as:
5.1where *b* and *d* correspond to the height and length of the element dragged through the medium. *C*_{S⊥}^{1} and *C*_{‖}^{1} are fitting parameters for the forces in Region 1 (values of the coefficients are given in electronic supplementary material, table S2). To derive an analytical solution in this region, we examine the numerical solution in figure 6 and find that in this region *η* ≈ 0, *θ* ≈ 0 and *ψ* ≈ *π* /2. Using these approximations and equation (5.1), we solve for *η* in equation (3.5):
5.2

In this region of small *A*/*λ* (or small *η*), as *A*/*λ* increases, *F*_{⊥}/*F*_{‖} decreases rapidly, which competes with increasing sin*θ* and results in the quadratic dependence of *η* on *A*/*λ*.

Region 2 (intermediate *A*/*λ*). For 0.05 < *A*/*λ* < 0.15, *η* increases substantially and rapidly. In this region as *ψ* decreases *F*_{⊥} increases slowly and *F*_{‖} is nearly constant:
5.3where *C*_{⊥}^{2}, *C*_{S⊥}^{2} and *C*_{S‖}^{2} are fitting parameters for the forces in Region 2. Solving for *η* gives
5.4

Physically, since *F*_{⊥}/*F*_{‖} decreases slowly, *η* in this region is more sensitive to variations in sin *θ*, and therefore for increasing amplitude *η* increases rapidly and approximately linearly.

Region 3 (large *A*/*λ*). For increasing *A*/*λ* > 0.15, *η* increases slowly. Here, as *ψ* decreases, *F*_{⊥} decreases rapidly while *F*_{‖} saturates at its maximum as the sides of the dragged rod are nearly parallel to its velocity for a majority of each cycle. This gives
5.5where *C*_{⊥}^{3}, *C*_{S⊥}^{3} and *C*_{S‖}^{3} are fitting parameters for the forces in Region 3. To derive an analytical solution, we assume *η* ≈ 1, *θ* ≈ *π*/2 and *ψ* ≈ 0 which gives,
5.6

Physically, the decrease in *F*_{⊥}/*F*_{‖} competes with the increase of sin*θ* which plateaus for large amplitude (as *θ* approaches *π*/2), which causes *η* to increase slowly.

Overall, the analytical solution of the RFT correctly predicts that *η* increases quadratically for small *A*/*λ*, increases rapidly for intermediate *A*/*λ* and is nearly constant for large *A*/*λ*. For each region, the analytical form of *η* qualitatively agrees with the calculated shape of the *η* versus *A*/*λ* relationship (figure 6*a*). The analytic solutions demonstrate that the difference between swimming in sand and swimming in a low-*Re* Newtonian fluid is predominantly because of the difference in the functional forms of *F*_{⊥} and *F*_{‖} in these two media (in addition to the difference that in granular media, forces are speed independent for low speeds, <40 cm s^{−1} in our study). In low-*Re* fluids [39], *F*_{⊥} and *F*_{‖} are proportional to their respective projected velocities with coefficients in the ratio of 2 : 1. In contrast, in granular media, the functional forms of *F*_{⊥} and *F*_{‖} are more complicated and thus the effective coefficient ratio depends on *A*/*λ*. The larger slope of *F*_{⊥} at smaller *ψ* (small *A*/*λ*) is largely responsible for the increased magnitude of *η* in granular media relative to that in low-*Re* swimmers like nematodes.

## 6. Optimal sand-swimming

All our models predict that the animal increases *η* by increasing *A*/*λ* (figure 6*a*,*b*), but we find that the animal does not operate at high *A*/*λ*. Maladen *et al.* [22] showed using the RFT model that locomotion at large *A*/*λ* (figure 6*c*) comes with a cost: since the animal is of finite length, the distance it travels per cycle decreases with increasing *A*/*λ*. We can see this by expressing body lengths travelled per cycle as *v*_{x}/*fL* = *η**λ*/*L*, as *v*_{x} = *η**f**λ*. While *η* increases with increasing *A*/*λ*, *λ* decreases and thus a maximum in forward progress per cycle can be expected. Figure 7*a*,*b* shows that as predicted by the RFT, both the numerical sandfish model and the physical and simulated robot models display a maximum forward progress per cycle at *A*/*λ* ≈ 0.2. As expected, the simulated sandfish with a tapered body shape moves faster than the square-tube shape but the values of *A*/*λ* that maximize their forward speed are close (figure 7*a*). Similar to the prediction of *η* versus *A*/*λ*, the RFT systematically over-predicts the speed per cycle compared with the numerical simulation, but scaling the net force on each segment (effectively the thrust) by 0.5 in the RFT results in an excellent match between the RFT and numerical simulation of the square tube body. The biological data reside close to the peak of the curve, indicating that the animal is maximizing its sand-swimming speed. This result agrees with the hypothesis that sandfish burial is an escape response [28].

To test the effect of the number of wave periods on performance, we fixed *A*/*λ* = 0.2 and varied 0.3 < *ξ* < 1.6. Testing was possible over only a limited range for the physical robot because for *ξ* < 0.6, the side walls of the test container interfere with the motion owing to large yawing of the robot and for *ξ* > 1.3, the maximum angle required to maintain *A*/*λ* = 0.2 exceeds the range of motion of the servos. The models progress forward fastest for approximately a single period along the body (figure 7*c*,*d*). The animal, as mentioned, operates close to the maximum speed predicted by the models.

We can explain the dependence of speed on *ξ*: for non-integer *ξ* the sandfish model experiences an unbalanced torque that results in a periodic yawing motion (with an amplitude of more than 40° for *ξ* = 0.5), which causes *A* to become effectively smaller and thus results in lower *η*. For fixed *A*/*λ*, as *ξ* increases *λ* decreases owing to the fixed body length. For *ξ* < 1, as *ξ* increases both yaw and *λ* decrease which has competing opposite effects on forward velocity and results in a maximum *v*_{x} at *ξ* < 1. For *ξ* > 1, as *ξ* increases, decreased *v*_{x} mainly results from the decrease in *λ*. We use the RFT to predict *η* only for a single period (*ξ* = 1) because determining *η* for non-integer *ξ* is challenging owing to unbalanced torques in the plane of motion.

## 7. Discussion and conclusion

Motivated by biological experiments of sandfish locomotion, we developed numerical modelling and robotic approaches to study undulatory sand-swimming. We used the models to test the predictions of our previously developed empirical resistive force model and demonstrated that limbless undulatory locomotion with *η* comparable with that of the animal can be achieved within granular media by models using the wave kinematics of the animal. We determined how *η* depends on amplitude and developed approximate analytical solutions that give insight into how the interdependent effects of thrust and drag forces, orientation of body elements and spatial characteristics of the animal influence *η*. With each model, we showed that the competing effects of wave efficiency and wavelength determine forward speed, and that *A*/*λ* ≈ 0.2 and *ξ* ≈ 1 maximize forward swimming speed. The sandfish uses these optimal kinematics to swim rapidly within granular media. The performance measured in experiment and predicted by each model obtained by constraining only kinematics (and not forces) was the same as that of the animal across a range of media properties like particle size and density.

The RFT over-predicts *η* for varying *A*/*λ* but the agreement of the functional form of the relationship compared with the non-tapered simulation is good. The cause of these performance differences may lie in the assumptions of the RFT. First, at larger *A*/*λ*, the assumption that the forces associated with each element are localized may be inaccurate. As the amplitude increases, the region of particles influenced by the adjacent segments may overlap, which could cause a reduction in thrust and thus performance.

Second, the force on the body of the model sandfish is different for the RFT and the numerical simulation. The net force on the body is nearly 50 per cent lower (for *A*/*λ* = 0.22) in simulation than estimated by the RFT because a segment decelerates as it reaches its maximum lateral excursion (*A*) and accelerates after the velocity direction reverses. This results in a varying velocity and a possible hysteresis effect as the lateral displacement reverses, which are not considered in the RFT. Hence, the deviation from the constant steady state velocity assumption of the RFT may be responsible for the over-prediction in *η* when compared with the numerical simulation for the *η* versus *A*/*λ* relation (figure 6*a*).

Third, an effect of possible importance at higher *A*/*λ* is that while the centre of mass (CoM) in both animal and numerical simulation oscillate laterally during locomotion, the RFT assumes that the CoM does not oscillate. Lateral CoM oscillation reduces the effective amplitude and changes the effective orientation (*ψ*) of each element which, according to our force measurements, should decrease *η*. Also yaw motion is not considered in the RFT.

Future comparison of the models will elucidate contributions from these three mechanisms and thus allow for better understanding of the physics of forces in granular media. With advances in our understanding of force generation in granular media, the RFT will be improved, potentially allowing it to accurately predict *η* for more complicated geometries such as the smoothly tapered body of the sandfish or the more discretized shape of the robot.

We now discuss the assumption of planar motion used throughout our modelling. Although the sandfish swims into the medium at angles between 20° and 30° relative to the horizontal, the RFT assumes that locomotion occurs in a horizontal plane; comparison with the numerical simulation and physical robot is made for the same horizontal orientation. Previous calculations showed that since resistive forces increase linearly with depth [47,48], the net thrust and drag generated by the body of the animal are unaffected by the entry angle. This suggests that the dominant effect of an angled swim is an increase in head drag. Based on measurements from side view X-ray data of the animal [22], we estimate that the force per area on the head is 1.5 times larger than on the centre of the body. With increased head drag, the RFT predicts that *η* at 22° should be approximately 15 per cent smaller than that for horizontal swimming. Simulation of a sandfish model constrained in a plane at 22° relative to the horizontal plane shows that swimming speed is approximately 16 per cent slower than in the horizontal case. This finding supports our RFT prediction that the dominant effect of swimming at an angle is an increase in head drag. The fact that the *η* predicted by the numerical model for the tapered body is higher (swimming at a fixed depth) (figure 1*e*) than that measured for the animal (swimming at 22°) is consistent with the increase in head drag for angled dives.

Our complementary modelling techniques have individual strengths and weaknesses. For example, while RFT can be rapidly solved to predict organism performance, predictions in different media and with different preparation require re-determination of the empirical force laws, which is inconvenient. In contrast, although the numerical simulation is slower (several days on a desktop PC), the experimentally validated simulation can give an understanding of these force interactions from the particle perspective and can be used to generate empirical drag laws for input into RFT. The simulation is also a flexible design tool with which parameters like particle–particle friction and animal kinematics can be easily varied to accurately predict performance. The robot is a key complement to the numerical and theoretical approaches as it validates predictions in a real world environment. However, unlike simulations, the robot is not amenable to extensive adjustments and cannot be used conveniently to measure particle forces and flowfields.

The framework of experimental visualization techniques, and theoretical, computational and physical modelling approaches that we have described here will allow us to develop and test mechanical hypothesis of interaction between the locomotor and its surroundings [54]. All models allow variation of parameters that are difficult to vary and control in organisms and can thus address questions concerning mechanical cost of transport, effect of body shape, skin friction, burial/unburial mechanics, manoeuvring, control methodologies [55], sensing and sand-swimming performance. These tools will also allow investigation of the detailed physics that governs the form of the force laws within granular media. The models can be used as design tools to develop the next generation of biophysically inspired sand-swimming robots that will explore complex flowing environments.

## Acknowledgements

We thank Mateo Garcia for characterizing the granular media, and Sarah Steinmetz for help with the animal experiments. We thank Andrew Slatton and Daniel Cohen for their assistance with the numerical sandfish simulation, and we thank Kellar Autumn for helpful discussion. This work was supported by The Burroughs Wellcome Fund Career Award at the Scientific Interface, NSF Physics of Living Systems grant PHY-0749991, and the Army Research Laboratory (ARL) Micro Autonomous Systems and Technology (MAST) Collaborative Technology Alliance (CTA) under cooperative agreement number W911NF-08-2-0004.

- Received November 30, 2010.
- Accepted January 28, 2011.

- This Journal is © 2011 The Royal Society