## Abstract

Over the past few decades, oscillating flexible foils have been used to study the physics of organismal propulsion in different fluid environments. Here, we extend this work to a study of flexible foils in a frictional environment. When the foil is oscillated by heaving at one end but is not free to locomote, the dynamics change from periodic to non-periodic and chaotic as the heaving amplitude increases or the bending rigidity decreases. For friction coefficients lying in a certain range, the transition passes through a sequence of *N*-periodic and asymmetric states before reaching chaotic dynamics. Resonant peaks are damped and shifted by friction and large heaving amplitudes, leading to bistable states. When the foil is free to locomote, the horizontal motion smoothes the resonant behaviours. For moderate frictional coefficients, steady but slow locomotion is obtained. For large transverse friction and small tangential friction corresponding to wheeled snake robots, faster locomotion is obtained. Travelling wave motions arise spontaneously, and move with horizontal speeds that scale as transverse friction coefficient to the power 1/4 and input power that scales as the transverse friction coefficient to the power 5/12. These scalings are consistent with a boundary layer form of the solutions near the foil’s leading edge.

## 1. Introduction

Snake locomotion has long been of interest to biologists, engineers and applied mathematicians [1–6], as the lack of limbs distinguishes snake kinematics from other common modes of locomotion including flying, swimming and walking [7–9], and results in distinctive dynamical behaviours [10]. Snakes gain thrust from the surrounding environment with a variety of gaits, including slithering, sidewinding, concertina motion and rectilinear progression [4]. Among these gaits, the slithering or undulatory motion is one of the most common and is used by many different animals, in both fluids and frictional environments. Some examples include swimming at high [8,11] and low Reynold numbers [12–15] and locomotion in granular media [16–19].

The dynamical features and locomotor behaviours of different undulatory organisms depend on how they interact with the environment, and in particular, depend on the type of thrust gained from the environment. In fluids, propulsive forces are obtained as a balance of fluid forces, bending rigidity and possibly body inertia when the body of a flexible swimmer is actuated [20,21]. For snakes, previous work showed how self-propulsion arises through effects such as Coulomb frictional forces and internal viscoelasticity [4,6]. Although there has been extensive work on locomotion in viscous and inviscid fluids, there are very few works discussing the locomotion of a passive flexible body in a dry frictional medium. One of the theoretical challenges compared to propulsion in a low-Reynolds-number fluid is that the frictional force is a nonlinear function of the body motion. Although the frictional force at a point depends only on the direction of the body velocity at that point, here the body moves according to geometrically nonlinear elasticity, so the dynamics have a nonlinear spatial coupling. Owing to Coulomb friction and elasticity, the problem is difficult to analyse theoretically except at small amplitudes. In general, the shape of the body and its velocity need to be solved computationally through a system of nonlinear equations.

Previous work approached snake locomotion by prescribing the motion of the snake in certain functional forms [4,6,22–24]. Guo & Mahadevan [6] studied the effects of internal elasticity, muscular activity and other physical parameters on locomotion for prescribed sinusoidal motions and sinusoidal and square-wave internal bending moments. Hu & Shelley [4] assumed a sinusoidal travelling-wave body curvature, and computed the snake speed and locomotor efficiency with different travelling wave amplitudes and wavelengths. Their results showed good agreement with biological snakes. Jing & Alben [22] found time-periodic kinematics of 2- and 3-link bodies that are optimal for efficiency, including travelling wave kinematics in the 3-link case. Alben [23] found the kinematics of general smooth bodies that are optimal for efficiency for various friction coefficients. Theoretical analysis showed that with large transverse friction, the optimal motion is a travelling wave (i.e. transverse undulation) with an amplitude that scales as the transverse friction coefficient to the −1/4 power. Wang *et al.* [24] studied transverse undulation on inclined surfaces with prescribed triangular and sinusoidal deflection waves. They found numerically and theoretically how the optimal wave amplitude depends on the frictional coefficients and incline angles in the large-transverse-friction-coefficient regime.

In this work, we propose a fully coupled model to solve for the dynamics and locomotion of an elastic foil as a generalization of previous work [23,24]. The elastic foil is a one-dimensional body whose curvature is a function of arclength and time. Body inertia, internal bending rigidity and frictional forces are all included in the model. A sinusoidal heaving motion is prescribed at the leading edge to move the foil, similarly to other recent experiments and numerical studies of swimming at low and high Reynolds numbers [13,25–27]. Passive flexible foils are a useful model system that has been used to study locomotion in fluids and granular media for a few reasons: they require a very simple control, such as harmonic heaving or pitching at the leading edge, they allow for generic body–environment interactions to arise spontaneously, and they mimic the behaviour of passive flexible bodies and tails commonly used for propulsion in swimming and crawling organisms. Alben *et al.* [25] studied the scalings of the swimming speeds of self-propelled flexible foils in inviscid flow and their dependence on fluid–structure resonances. The importance of resonances on propulsion was also highlighted by Quinn *et al.* [28] in their study of heaving rectangular panels in a water channel. Similar systems have been used to study effects of foil shape and stiffness, centre-of-mass oscillations, wall effects, and vortex dynamics for oscillating foils in high-Reynolds-number flows [29–42]. Flexible foils have also been used to investigate propulsion in low-Reynolds-number flows [20,21]. In this work, we will use passive flexible foils to study propulsion in a frictional environment.

Our foil moves according to nonlinear force balance equations which will be solved numerically. To simplify the discussion, we first consider the case where the foil is actuated at the leading edge but is not free to translate (fixed base), and study the dynamics of the foil at different parameter values. Then we allow the foil to move freely in the horizontal direction and consider the kinematics of the locomotion. This approach has been used to study an actuated elastica in a viscous fluid [21], where the fixed-base case was used to derive a scaling law for the propulsive forces and the free translation case was used to calculate the swimming speed. In the fixed-base case, we will study some key dynamical phenomena including resonant vibrations at certain bending rigidities, effects of nonlinearity due to large heaving amplitudes and geometrical nonlinearities, and transitions from periodic to non-periodic states. Similar phenomena have been discussed previously for mechanical vibrations [43] and for swimming in a fluid medium [44–46], but not in a frictional environment. In the freely locomoting case, we will focus on the speed and input power for locomotion, and discuss how they scale with physical parameters in a frictional medium, as a complement to recent studies of similar quantities in low- and high-Reynolds-number swimming systems [21,25,28,47].

We also note that the study of flexible foils in a frictional medium can be applied to various situations where Coulomb friction applies. One example is a granular medium where the resistive force can be modelled using Coulomb friction in the regime of slow movement [19,48]. Peng *et al.* [17] applied this model to investigate the locomotion of a slender swimmer in a granular medium by prescribing a travelling wave body shape, and found the optimal swimming speed and efficiency versus wave number. In a more recent work, Peng *et al.* also considered propulsion in a granular medium under the effects of both elasticity and frictional force. They proposed a model where a rigid rod is connected to a torsional spring under a displacement actuation [18], studied the effects of actuation amplitude and spring stiffness on propulsive dynamics, and found the maximum thrust that could be obtained.

This paper is organized as follows: The foil and friction models and the numerical methods are described in §2. The fixed-base cases are discussed in §3 followed by the free locomotion cases in §4. Conclusions and a discussion are given in §5.

## 2. Modelling

### (a) Foil and friction models

We consider here the motion of a flexible foil in a frictional environment with a prescribed heaving motion at the leading edge. The foil has chord length *L*, mass per unit length *ρ* and bending rigidity *B*. The foil thickness is assumed to be much smaller than its length and width, and therefore we model it as a one-dimensional inextensible elastic sheet.

The instantaneous position of the foil is described as *ζ*(*s*,*t*)=*x*(*s*,*t*)+*iy*(*s*,*t*), where *s* is arclength. Assuming an Euler–Bernoulli model for the foil, the governing equation for *ζ* is
*T*(*s*,*t*) is a tension force accounting for inextensibility, and *f*(*s*,*t*) we denote the frictional force per unit length, and from previous work [4,23,24], it can be described as
*μ*_{f} and *μ*_{t} for motions in the tangential

At the leading edge, the transverse position of the foil is prescribed as a sinusoidal function of time with frequency *ω* and amplitude *A*, and the tangent angle *θ* is set to zero:
*X*_{0}(*t*)≡0. For a locomoting body, *X*_{0}(*t*) is computed by assuming no horizontal force is applied at the leading edge: *T*(0,*t*)≡0. Similar clamped boundary conditions have been used in previous experiments and models to study swimming by flexible foils [25,27,44,45]. We note that other choices of boundary conditions have also been applied. For example, a pitching motion, where the tangent angle is a sinusoidal function of time while the transverse displacement is fixed to be zero, is another popular choice [13,20,21]. A torsional flexibility model, in which a torsional spring is connected to a rigid plate at the leading edge, has also been applied in some works [18,49].

At the trailing edge, the foil satisfies free-end conditions, which state that the tension force, shearing force and bending moment are all zero:

### (b) Non-dimensionalization

We non-dimensionalize the governing equations and boundary conditions (2.1)–(2.4) by the density *ρ*, the chord length *L* and the period of the heaving motion, *τ*=2*πω*. We obtain the following dimensionless parameters:

We note that now the heaving motion always has a period of 1. In the following sections, we will discuss how solutions depend on the four parameters: the bending rigidity *B*, heaving amplitude *A*, transverse frictional coefficient *μ*_{t} and tangential frictional coefficient *μ*_{f}.

### (c) Numerical methods

We couple equations and boundary conditions (2.5)–(2.7) together to solve for the positions of the body *ζ*(*s*,*t*) at each time step. This requires solving a nonlinear system **F**(**x**)=0 at each time step, and Broyden’s method [50] is used to do so, which essentially requires the evaluation of the function **F**(**x**) for a given **x**.

At each time step *t*_{n}, given *κ*_{n}, we first obtain *ζ*_{n}:
*s*=1, the tension force can be computed as
*s*=1, we obtain the curvature:

We use the notation *a*,*b*], if we approximate *u*_{s} and *u*_{n} by linear approximations *As*+*B* and *Cs*+*D*, then the integrals are of the form

## 3. Fixed base

Because little is known about the dynamical behaviour of an elastic body in a frictional medium, we begin with the simplified case in which the base is fixed (*X*_{0}(*t*)≡0), and then study the case of a freely locomoting foil in the following section.

### (a) Zero friction

When the heaving amplitude *A* is small, the nonlinear equation (2.5) can be linearized by assuming *s*≈*x* and *ζ*(*x*,*t*)≈*x*+*iy*(*x*,*t*). Therefore, equation (2.5) becomes
*ϕ*_{i}(*x*) is an eigenfunction and *A*_{i}(*t*) is a time-dependent coefficient (see details in appendix B in the electronic supplementary material). For arbitrary initial conditions, the shape of the foil is non-periodic in time, with a superposition of the heaving frequency 1 and natural frequencies *A* is small, and we validate our numerical scheme by comparing with the linearized model in appendix B.

### (b) Non-zero friction

The frictional force adds damping to the system, and therefore damps out the initial transient, leaving a periodic solution (at small enough *A*) with energy input by heaving and removed by friction. A larger heaving amplitude *A* enhances other vibration modes which introduces more frequencies into the system and eventually results in a non-periodic motion. A more flexible foil (i.e. a smaller *B*) also leads to a non-periodic solution. Therefore, the parameters *A*, *B* and the two friction coefficients will compete to determine whether the vibration will be periodic or not. For simplicity, we only consider homogeneous friction coefficients *μ*=*μ*_{t}=*μ*_{f} in this section. For snakes and snake-like robots, these parameters are generally not the same [4].

We first consider the effect of varying *B* and *μ* with fixed heaving amplitude *A*. In figure 2, we plot a diagram of different dynamical states after 200 periods with *A*=0.3. When *μ* is relatively small (less than 0.7 in this case), the vibration transitions directly from a periodic state (with period 1) to a non-periodic chaotic state as the foil becomes more flexible, as shown in figure 2*a*. The transition value of *B* is almost invariant for a large range of *μ*, becoming slightly smaller as *μ* increases. However, when *μ* is moderate (0.7<*μ*<7 for this case), a transition region is observed between the non-periodic and period-1 states as shown in figure 2*a*,*b*. The transition region has complex dynamical behaviours with a mixture of different periodic and non-periodic states.

We analyse the transition region for *A*=0.3 and *μ*=1 as an example. The complex behaviour of the system is shown in figure 3. In figure 3*a*, we plot both the positive and negative vibration amplitude *A*_{v}, i.e. the positive and negative local maxima of the vertical displacement of the free end within one period, for over 50 periods. The power spectrum density versus *B* for the free end displacement is shown in figure 3*b*, where the power spectrum for each case is normalized by the corresponding maximum *A*_{v}. In figure 4, we choose an example from each periodic and non-periodic state, and plot the phase plot of the free end velocity in the vertical direction ∂_{t}*y* against the vertical displacement *y* and the corresponding snapshots of the foil in 50 periods.

As shown in figures 3 and 4, when the foil is rigid (large *B*), the vibration is periodic and symmetric (region (1) in figure 3) with only a change in the amplitude. As *B* decreases, the symmetry of the vibration about the *x*-axis is broken. Depending on the initial condition, the free end can move more to negative or positive *y*. The vibration goes through several different asymmetric periodic states as *B* continues decreasing, with period 1 (region (2)), period 13 (region (3)), period 10 (region (4)) and period 5 (region (5)), and becomes symmetric again with period 1 in region (6). These six different periodic states can be viewed as a perturbation of the symmetric period 1 vibration, as shown by the shapes of the six foils in figure 4, panels (1)–(6). When *B* becomes even smaller, the foil is flexible enough to allow more complex deformations and can flip over and vibrate in the left half-plane (*x*<0). The vibration goes through several periodic states (region (7) with period 3 and region (9) with period 9) and non-periodic states (region (8)) until it becomes non-periodic and chaotic in region (10) and beyond. The dynamics of the foils become qualitatively different from those in regions (1)–(6).

When *μ* becomes even larger (greater than 7 for *A*=0.3), the transition region disappears again. The vibration changes from period 1 for a more rigid foil to non-periodic for a flexible foil directly, as shown in figure 2*b*. However, at large *μ* and large *B*, the vibration stays asymmetric, which is different from the symmetric motions at moderate *μ* and large *B*.

Next, we consider the effect of *A* on the periodicity of the vibration. We fix the value of *B* at 0.5 and plot the diagram of different states of vibration after 200 periods with various *A* and *μ* in figure 5. When *μ* is small, moderate, or large, the system exhibits different dynamical behaviours as *A* increases. When *μ* is small, the vibration transitions from periodic to non-periodic directly as the heaving amplitude increases. The critical *A* is almost invariant for different *μ*<1. For moderate *μ*, a transition region with different periods and symmetry is observed as *A* increases. When *μ* is large enough, only period-1 vibration is observed in the periodic regime. However, the foil becomes asymmetric first before it becomes non-periodic.

We note that similarly complicated dynamical behaviour was also observed in flapping foil vibrations in fluids. Chen *et al.* [51] identified five different periodic states and three chaotic states when a flag transitions from periodic flapping to a non-periodic state under flow-induced vibrations. Spagnolie *et al.* [46] studied flapping locomotion with passive pitching in a viscous fluid, and observed a bistable regime where the wing can move either forwards or backwards depending on its history. An asymmetric pitching motion was also observed in their work. In different systems, the dynamical behaviour depends on different physical parameters such as the driving frequency, amplitude and Reynolds number. In our work, reduced heaving amplitude, frictional coefficients and bending rigidity are the most important parameters.

### (c) Resonance

When the heaving frequency matches one of the natural frequencies of the vibration, resonance occurs and the amplitude of the oscillation grows linearly with *t*. In the non-dimensionalized model, the heaving period is always 1. In the linearized model, the natural frequency depends on the rigidity of the beam *B*. Therefore, a resonance occurs when
*ω*_{i} are constants given in appendix B in the electronic supplementary material. With the nonlinearities introduced by non-zero frictional forces and larger heaving amplitude *A*, the resonant *B* values will vary.

We first consider the effect of the frictional force on the resonance. In figure 6*a*, we plot the free end amplitude *A*_{v} versus the foil rigidity *B* for a fixed heaving amplitude *A*=0.05 and various *μ*. We plot the free end amplitude in the linearized model with no friction with a dashed line. As *B* decreases, multiple resonances are observed with *μ*=0 according to equation (3.5). We only plot the first three of the infinite sequence of resonances in the panel. As *μ* increases, the amplitude *A*_{v} decreases correspondingly. The resonant *B* values shift to the right as *μ* increases except for the first resonance. As *B* decreases, the shape of the foil changes from the first mode to higher bending modes. In figure 6*b*–*e*, we plot snapshots of the foil in one period for *A*=0.05, *μ*=0.1 and *B*=6, 1, 0.09 and 0.02, respectively, showing the different modes.

Next, we consider the effect of *A* on the resonance. In figure 7*a*, we plot the free end amplitude *A*_{v} versus the bending rigidity *B* with fixed *μ*=0.1 and various *A*. We only consider symmetric vibrations with period 1 in this figure; the curves in figure 7*a* end when the vibration becomes either non-periodic or asymmetric for the particular parameter values.

When *A* is small (*A*<0.1 in this case), the free end amplitude *A*_{v} increases with increasing *A*. Moreover, the resonant *B* value shifts to the right as *A* becomes larger at the second (*B*≈0.1) and third peaks (*B*≈0.01). When *A* becomes larger, we observe bistability near the first resonant value (*B*≈3). In figure 7*b*, we enlarge the region near the first resonant value to show the bistability.

Typical motions in the bistable regime are shown in figure 7*c*,*d*, which show the snapshots of the foil in one period for *μ*=0.1, *A*=0.3 and *B*=3.05. The snapshots corresponding to the lower branch are shown in figure 7*c* and those corresponding to the upper branch are shown in figure 7*d*. For the lower branch, as the leading edge of the foil moves upwards by the heaving motion, the trailing edge moves downwards. For the upper branch, the trailing edge moves in the same direction as the leading edge. This phenomenon is observed for both *A*=0.3 and *A*=0.5, and as *A* increases, the bistability region becomes larger.

The bistability is observed for different values of *μ* as long as *A* is large enough. It is a result of the nonlinearity of the system. Such ‘double-fold’ bifurcations near a resonant peak have been observed in other nonlinear oscillators such as the Duffing oscillator [52].

In general, frictional force increases the damping of the system, increasing the rate of convergence to periodic motions (and their extent in parameter space) and smoothing the resonant peaks. However, decreases in the foil rigidity and increases in the heaving amplitude both increase nonlinear effects and therefore lead to more chaotic behaviour and bistability.

## 4. Free locomotion

We now consider free locomotion of the foil. The foil is still actuated by a periodic vertical heaving motion *x*-direction with position *X*_{0}(*t*). Self-propulsion occurs naturally as a balance between inertia, the internal bending rigidity, tangential friction and transverse friction.

The foil shape *ζ*(*s*,*t*) and the horizontal position *X*_{0}(*t*) will be solved simultaneously by modifying the numerical method from the fixed-base case. To solve for the extra unknown *X*_{0}(*t*), we add one more equation to the system by assuming no tangential force (tension) is applied at the leading edge:

As the foil bends, a horizontal force is obtained from transverse friction, and we expect the foil to move horizontally in general. We define the space and time-averaged horizontal velocity as _{t}*y* is the vertical velocity component and *B*∂_{s}*κ* is the shearing force in the vertical direction at the leading edge.

### (a) Effects of transverse friction

We first consider the effect of the transverse frictional coefficient on *B* for fixed *μ*_{f}=0.01, *A*=0.1, and various *μ*_{t} such that *μ*_{f}≪*μ*_{t}. When *μ*_{t} is small (less than 5), a resonant peak is obtained in *B*≈2 and corresponds to a decrease in the velocity. At larger *μ*_{f} (not shown), the foil has a smaller horizontal speed (unsurprisingly), and stronger resonant-like behaviours. Both features are present in the lower curves in figure 8*a*, and these features are strengthened as *μ*_{f} increases. For a real snake, *μ*_{f}<*μ*_{t} but both have been measured in the range 1–2 (when the snake is anaesthetized) [4,53,54]. In this regime, our passive elastic foil translates slowly (≈0.1 body lengths per period), and for certain values of *B* moves rightwards (towards the free end) at large amplitudes (*A*>0.1). Owing to the slow speed of locomotion, the foil behaviour has strong similarities to the fixed-base case.

The upper curves in figure 8*a* (*μ*_{t}≥5) tend towards the case of wheeled robots with large transverse friction and small tangential friction [5], where the foil has a higher speed and efficiency. We show how the foil motion changes from small to large *μ*_{t} in figure 9*a*–*c*. We plot snapshots of the foil as well as the trajectory of the leading edge in one period with fixed *A*=0.1 and *μ*_{f}=0.01, and various *μ*_{t}=1,10 and 100 and *B*=0.5,2.5 and 10. At the right of each panel, we relocate the foil so that the snapshots at different time instants share the same leading *x* location, to better illustrate the mode shapes. In figure 9*a*, *μ*_{t}=1, and the foil shows two different mode shapes for *B*=0.5 and 10. Near the resonant peak at *B*≈2.5, the foil vibrates in a large-amplitude motion mainly along the transverse direction. There, the *x* velocity decreases significantly, while the input power increases greatly. As *μ*_{t} increases to 10, the foil transitions to a different dynamical regime, where the velocity and the input power vary more smoothly and the resonant peaks are reduced as shown in figure 8. In figure 9*b*, the differences in the mode shapes are reduced compared to figure 9*a*. When *μ*_{t}=100, the foil deflection is much smaller as shown in figure 9*c*. When *μ*_{t} is large enough, transverse friction dominates the system, while the foil is close to a flat plate as *a*, we see that *B* becomes larger, and a maximum speed is obtained at a moderate *B* value. In figure 8*b*, we observe that the input power *B* increases for a fixed value of *μ*_{t}, because the foil converges to a purely vertical oscillation and the work done against transverse friction is independent of *B* for a rigid plate.

Figure 9*c* indicates that, in the large *μ*_{t} regime, the motion of the foil can be approximated by a travelling wave solution, as the snapshots of the foil seem to follow a certain wave track. In one of our previous works [24], we prescribed the shape of the snake as a travelling wave and found the optimal wave amplitude for efficiency. In [23], we showed asymptotically that a travelling wave motion minimizes the cost of locomotion at large *μ*_{t} when neglecting inertia. In the current model, we do not prescribe the shape of the foil as a travelling wave. Instead, the foil shape and the travelling speed arise naturally from the dynamical force balance equation. It is interesting that the flexible foil spontaneously adopts a travelling wave motion at large *μ*_{t}. To clearly illustrate the travelling wave motion, we show a contour plot of *y*(*x*,*t*) versus *x* and *t* within one period in figure 9*d*, for *A*=0.1, *B*=10, *μ*_{f}=0.01 and a large *μ*_{t}=100. For *x* away from the leading edge, we find that the contour curves are close to straight lines, which indicates that *y* is of travelling wave form *g*(*x*−*U*_{w}*t*) where *U*_{w} is a wave speed. Deviations are observed for *t* near 0.25 and 0.75, when *y* reaches its extrema. This is reasonable because the velocity of *y* changes sign at its extrema and the deflection of the foil cannot be characterized as a travelling wave there. We also note that, at the leading edge, the foil is held flat (∂_{x}*y*=0) at all times, while the travelling wave has non-zero slope (∂_{x}*y*≠0). This can be seen in how the contours in figure 9*d* change from zero slope at *x*=0 to non-zero slope (=*U*_{w}) for

Along with a boundary layer structure, the approximate travelling wave solutions at large *μ*_{t} also obey certain scaling laws. Two of the most important quantities are the horizontal speed *μ*_{t}. Here, the solution at large *μ*_{t} spontaneously adopts a travelling wave form, so it is reasonable to expect for the current problem that *μ*_{t}. In figure 10, we give numerical evidence for the power-law scalings with *a* and *b*. By assuming small slopes (|∂_{x}*y*|≪1) and approximate travelling wave solutions outside of boundary layers at the leading edge, we now explain these scaling laws.

When the amplitude *A* is small, we can simplify the foil model in the large limit of *μ*_{t} by approximating the arclength *s* by *x*, and *ζ*(*s*,*t*)≈*x*+*iy*(*x*,*t*). At leading order, the tangent and normal vectors are
_{t}*x* has small variations over arclength and time, and therefore _{t}*ζ*(*x*,*t*)≈*U*+*i*∂_{t}*y*. In the limit of large *μ*_{t}, *U*≫∂_{t}*y*, as we expect the vertical velocity ∂_{t}*y* is proportional to the heaving amplitude *A* while *U* grows with *μ*_{t}. Therefore, the normalized velocity at leading order is *U*<0). The normalized tangential and normal velocity components are, to leading order,
*y*-component of equation (2.5), neglect higher-order terms and obtain a force balance for the vertical motion:

As mentioned already, the leading edge clamped boundary condition is not compatible with a travelling wave solution, so we look for a boundary layer near the leading edge. In figure 11*a*, we plot ∂_{x}*y* versus *x* for *B*=10, *A*=0.1 and *μ*_{f}=0.01 and various *μ*_{t}. We find that ∂_{x}*y* scales as *μ*^{−1/4}_{t} from the numerical simulations. As ∂_{x}*y*=0 when *x*=0, ∂_{x}*y* will increase from 0 to _{t}*y* and acceleration ∂_{tt}*y* are *O*(1) (∼*A*) near the leading edge. Thus, according to equation (4.4), we have the following scalings within the boundary layer:
*α*. Using ∂_{x}*y*∼*μ*^{−1/4}_{t} and *α*=−1/3, and the width of the boundary layer scales as *b*,*c*, we plot ∂_{x}*y* scaled by *μ*^{1/4}_{t}, and *μ*^{−3/4}_{t} versus *xμ*^{1/3}_{t} and find a good collapse for both quantities within the boundary layer, particularly at larger *μ*_{t}.

The time-averaged input power _{x}*y*∼*μ*^{−1/4}_{t} twice with respect to *x*. Each differentiation divides by a factor of *μ*^{−1/3}_{t}, the boundary layer width. Consequently, _{t}*y*∼*A*∼1, *μ*^{5/12}_{t}. This scaling is confirmed by the numerical results in figure 10*b*.

### (b) Effects of heaving amplitude

We now briefly mention the dependence of locomotion on the heaving amplitude *A*. For *a*, we plot the horizontal speed for *μ*_{t}=1. Near the first resonance (*B*≈3), the speed has a local minimum, and as we increase *A* from 0.03 to 0.1, the minimum does not move much, but the trough around the minimum broadens. Increasing *A* further to 0.3 and 0.5, strong nonlinear effects come into play. As the foil mainly moves in the vertical direction and has a small horizontal speed, the motion of the foil is similar to that of the fixed-base case. We find a bistable region for *A*=0.3 near the resonant peak, consistent with the fixed-base case (figure 7*a*,*b*). Insets are shown in figure 12*a*,*b* to enlarge the bistability region (1.5≤*B*≤1.8). In the lower bifurcation branch, the foil reverses its horizontal direction and moves rightwards (*a*,*b*, at *B*=1.6, which shows a reverse (rightward) motion and a regular (leftward) motion. The two flapping modes are clearly quite different. For the reverse motion, the leading and trailing edges mainly move in the same direction vertically, while for the regular motion, they move in the opposite direction. This is consistent with the results we obtained in the fixed-base case (figure 7*c*,*d*). In figure 12*b*, we plot the input power for *μ*_{t}=1 near the first resonance. As *A* increases, the resonant peak broadens and becomes less symmetrical, similarly to the fixed-base case. For *A*=0.3, we also observe a bistability region, which we enlarge in a subplot. For *A*=0.5, a similar bifurcation is expected, but we did not search for it at the same resolution as for the *A*=0.3 case. For foils with large flexibility (*B*<1.8), the locomotion eventually becomes non-periodic, so the time-averaged velocity and power are not defined. Therefore, we neglect these data points here.

We note that a similar bistable behaviour was observed by Spagnolie *et al.* [46] when studying flapping locomotion with passive pitching in a fluid. When the pitching frequency exceeds a critical value, the flapping body can move either forwards or backwards.

For comparison, we plot the horizontal speed and input power with *μ*_{t} increased to 100 in figures 12*c*,*d*. The curves are much smoother as noted previously, and resonance and bistability are no longer observed. As *A* increases from 0.03 to 0.1, the speed increases almost in proportion to *A*. This is consistent with the travelling wave solutions we have found, because the travelling wave amplitude and wavelength are both proportional to *A*. As the foil speed is approximately the wavelength divided by the period (fixed to 1), it too is proportional to *A*. As *A* (0.3 and 0.5), the stronger nonlinearities in the equations lead to nonlinear changes in *c*) and *d*).

## 5. Conclusion and discussion

In this work, we have studied the dynamics and locomotion of an elastic foil in a frictional environment, heaved sinusoidally in time at the leading edge. To understand the basic physics, we began with the case where the base does not locomote horizontally. The foil dynamics depend on three key parameters: the heaving amplitude *A*, the bending rigidity *B* and the frictional coefficients (set to the same constant *μ* for simplicity). With *μ*=0 and small *A*, we obtain the well-known case of an elastic beam in a vacuum, whose motion is a transient term, a superposition of eigenmodes at the natural frequencies set by the initial condition, together with a term set by the applied heaving motion. With non-zero *μ* and small *A*, the transient is damped, leaving a periodic solution. At larger *A*, there is a rich set of nonlinear behaviours that can be summarized succinctly in phase space. As *B* is decreased from large to small, the foil transitions from a motion periodic with the driving period to a non-periodic, chaotic motion. For a band of *μ* values near unity, the transition passes through a set of states that are periodic at various multiples of the heaving period, and with or without bilateral symmetry. As *A* increases from small to large, the same types of transitions occur, again with *N*-periodic states appearing at a finite band of *μ* values. At zero *μ* and small *A*, resonances occur at particular *B* values. The resonant peaks are damped nonlinearly with increasing *μ* and *A*, and bistability is observed at large *A*.

Allowing the base to translate freely in *x*, and allowing different tangential and transverse friction coefficients (*μ*_{f}≠*μ*_{t}), we find many of the same dynamics (e.g. chaotic motions, resonances) together with small horizontal velocities (*μ*_{f} and *μ*_{t} do not differ much in magnitude. In the regime *μ*_{t}>1≫*μ*_{f} corresponding to wheeled snake robots, the foil spontaneously adopts a travelling wave motion with high speed ∼*μ*^{1/4}_{t}, though the input power grows faster, ∼*μ*^{5/12}_{t}. We find that the motion has a boundary layer form near the leading edge in powers of *μ*_{t}, consistent with the speed and input power scalings. The input power scaling in particular depends on the clamped boundary condition, which sets the slope to zero at the leading edge. We hypothesize that other leading edge conditions, such as a pinned leading edge (zero curvature), may lead to different scalings and perhaps higher locomotor efficiency. We leave this for future study, as well as comparisons with experimental work, and inclusion of proprioceptive feedback from the environment into the driving motion, as employed in other locomotion studies [55–57].

We also point out that the foil motions found here are remarkably similar to those which were found to be optimal for locomotor efficiency when the foil motion was fully prescribed at large transverse friction [23,24]. Both motions are approximate travelling waves, and in both cases the foil slope ∂_{x}*y*∼*μ*^{−1/4}_{t}. In this work, the input power grows rapidly with *μ*_{t}, ∼*μ*^{5/12}_{t}, due mainly to the necessary deviation from a travelling wave at the leading edge. For the optimal motion, the leading edge amplitude decays as *μ*^{−1/4}_{t} (versus *O*(1) here), and the input power has the same *μ*_{t}-scaling as the forward speed (both are *O*(1)). For more advantageous choices of leading edge heaving and pitching, the passive elastic foil could approach the optimal foil’s performance with respect to *μ*_{t}.

The oscillating foil with friction is a new example of spontaneously self-organized locomotion of a flexible body. Unlike locomotion in a low-Reynolds-number fluid, chaotic dynamics are generic here. Compared to inviscid models of high-Reynolds-number fluids, the frictional model remains valid for larger-amplitude body motions (up to the point of self-intersection), allowing us to describe the transition to chaotic dynamics at large amplitude and small bending rigidity, across a wide range of frictional parameters. Near resonances we find that bistability can arise, and for the freely locomoting case, speed decreases and power increases. When the foil can locomote freely and transverse friction is large, high speeds arise naturally and resonances are suppressed.

## Ethics

This work does not involve any active collection of snake locomotion data, but only computer simulation of snake behaviour.

## Data accessibility

This work does not have any experimental data. The codes for running the simulation are attached in the electronic supplementary material.

## Authors' contributions

X.W. and S.A. formulated the mathematical models and developed the numerical methods. X.W. wrote and debugged the simulation codes, performed the simulations, and generated the results and figures. X.W. and S.A. contributed to the writing of the manuscript and to the analysis of the results. Both the authors gave their final approval for publication.

## Competing interests

We have no competing interests.

## Funding

This work was supported by NSF-DMS Mathematical Biology grant no. 1329726, a Sloan Research Fellowship to S.A., and the NSF REU program at the University of Michigan Department of Mathematics.

## Acknowledgements

We thank Angelia Wang for performing simulations in a preliminary study of the fixed-body dynamics.

## Footnotes

Electronic supplementary material is available online at http://dx.doi.org/10.6084/m9.figshare.c.3951913.

- Received July 27, 2017.
- Accepted November 29, 2017.

- © 2018 The Author(s)

Published by the Royal Society. All rights reserved.