Locomotion of Self-Excited Vibrating and Rotating Objects in Granular Environments

: Many reptiles, known as ‘sand swimmers’, adapt to their speciﬁc environments by vibrating or rotating their body. To understand these type of interactions of active objects with granular media, we study a simpliﬁed model of a self-excited spherical object (SO) immersed in the granular bed, using three-dimensional discrete element method (DEM) simulations. Modelling the vibration by an oscillatory motion, we simulate the longitudinal locomotion of the SO in three modes: transverse vibration, rotation around different axes, and a combination of both. We ﬁnd that the mode of oscillation in y direction coupled with rotation around x -axis is optimal in the sense that the SO rises fastest, with periodic oscillations, in the z direction while remaining stable at the initial x position. We analyze the physical mechanisms governing the meandering up or down and show that the large oscillations are caused by an asynchronous changes between the directions of oscillation and rotation. We also observed that the SO’s rising rate is sensitive to three parameters: the oscillation amplitude, the oscillation frequency, f , and the rotation angular velocity, Ω . We report the following results. 1. When the frequencies of the rotation and transverse motion are synchronised, SO rises when Ω < 0 and sinks when Ω > 0; the average rising/sinking rate is proportional to | Ω | . 2. The rising rate increases linearly with the oscillation amplitude. 3. There exists a critical oscillation frequency, above and below which the rising mechanisms are different. Our study reveals the range of parameters that idealized ‘swimmers’ need to use to optimize performance in granular environments.


Introduction
Granular Materials (GMs) are common in nature. A GM is defined as an assembly of many discrete macroscopic solid particles that interact via a hard-core potential and dissipate energy on collisions [1]. GMs can behave like a solid, supporting loads, or flow like a fluid [2]. However, full continuum theories for the solid and fluid [3] behaviors are still to emerge. Therefore, understanding GM in general is still a challenging problem.
Some sand lizards, such as Scincus scincus, are well known to 'swim' through sand beds by effecting a solid-to-fluid and back transition of the media around them [4,5]. Such lizards can both dive into desert sand, to hide from heat and predators, and travel beneath the dry sand surface with considerable speed over reasonable distances [6][7][8]. These capabilities are enabled: by physiological attributes: the spatula-shaped snout, the streamlined body shape, the smooth integument, and fringed digits [4,9,10], which are well adapted to life in the desert, and by the rapid lateral oscillations that loosen the sand surrounding the swimmer's body, leading to vertical burrowing. Another mechanism of moving in GM is by helical rotation, such as employed by Pelargonium and Erodium seeds for digging into cohesionless soils [11,12]. This mechanism is reminiscent of a common propulsion method of many bacterial species, such as E. coli, that propel themselves by rotating bundles of helical filaments [13,14]. Thus, discovering mechanism of locomotion for organisms in GMs is of great engineering significance, not least to the development of walking robots in sandy terrains.
Recently, Maladen et al. [7] have studied the swimming of sandfish Scincus scincus in sand using high-speed X-ray imaging. Takashi et al. [15] have constructed a simplified model of periodic contraction and extension of large intruders in granular beds, using an event-driven simulation. A number of researchers applied the resistive force theory (RFT) [16,17] in fluid mechanics to model the sandfish swimming dynamics [18,19] and concluded that RFT predicts reasonably well the motion on granular surfaces [20]. Legged sand walking robots and sand swimming robots have been developed, based on experiments and simulations [21,22]. However, the aforementioned studies focused mainly on the in-plane motion the walking or swimming, less is known about the mechanics of large-scale vertical motion within the GM.
The main aim of this work is to investigate the vertical behavior of a self-energized SO within a dense granular medium. In this paper, we report three-dimensional DEM simulations of the locomotion of vibrating and/or rotating SOs, studying their trajectories under various conditions. To the best of our knowledge, this is the first work to study such a coupled mode of motion, which could prove very efficient for subsurface locomotion. Specifically, we studied the effect of three parameters: A, Ω, and f , on the rise/sink of SO in this mode of motion. The results of this study have the potential to advance the design of a bio-inspired robot to rapidly rise to the surface or dive into sand to a significant depth. The rest of the paper is organized as follows. In Section 2 we describe the numerical simulation method and the system. In Section 3 we present all the simulation results and includes our analysis and discussion. In Section 4, we present our conclusions and suggest future work.

Simulation Method
Numerical simulation of GMs are used extensively both as a cheap alternative to experiments and because they enable good resolution of grain-scale analysis. Advances in computer processing speed continue making this approach ever more attractive, as well as the inherent difficulties in constructing continuum theoretical models for the dynamics of, and in, non-cohesive GMs [23]. In particular, the DEM [24][25][26] has become a powerful and reliable tool for numerical modelling of granular materials. The rationale underlying DEM simulation of granular dynamics is based on the dynamics of collisions between particles. Specifically, the method requires a time-discretized form of equations of motion for particles displacements and rotations, as well as a force-displacement relation for particle interactions [27].
In the DEM method, the simulation starts from an initial state of a collection of particles within a designated region. The particles have mathematically well-defined shapes, which can be arbitrary, in principle. Particle properties are used to calculate interaction forces between particles on collision [28]. The (intergranular) interactions are often modelled as linear springs, dashpots, and joints [28,29], whose normal and tangential components (with respect to the overlap) are given by, using Hertz contact model, where the δ is the normal overlap, α is the tangential displacement of the colliding particles from the initiation of the contact, µ is the particle-particle friction coefficient, and v n and v t are the normal and tangential component of the relative velocity at the contact points, respectively. The parameters k n , k t , η n , and η t represent, respectively, the normal and tangential hardness and viscoelastic coefficients, which are functions of material properties: The material characteristics in (3), are calculated as follows, where Y is Young's modulus, ν Poisson's ratio, e the restitution coefficient, G the shear modulus, R the particle's radius, and m its mass. Using the forces, determined from Equtions (1) and (2), in Newton's second law, a particle's acceleration, a, is determined and, from it, the particles velocity and displacement increments from time t to time t + ∆t, using the Velocity-Verlet integration method: This is done at each time step for each particle.

Simulation Software
We used the open-source software LIGGGHTS [30,31] to simulate numerically a SO moving within a granular medium. This software is particularly convenient for simulating spherical objects [32]. For simplicity and to avoid using too many parameters, we represent the self-excited object as a sphere, larger than the bed particles. The spherical object (SO) can vibrate, rotate or do both simultaneously. We investigated in the simulations the dynamics of the SO, in particular the effects on its vertical locomotion, as well as of the granular medium surrounding it. To compare with existing experiments, we parameterized the particle material properties to imitate SO dynamics in sand or glass bead beds, for which standard values are used in Molecular Dynamic simulations [33]. The material properties of the SO are adjusted to achieve the largest possible rise and sink rates. Table 1 summarized the values of the material properties of the bed particles and SO used in our simulations.The positions of the particles are updated in discrete time steps [34]. No unit 0.9 (bed particle-bed particle) 0.5 (SO-bed particle) Coefficient of Friction No unit 0.45 (bed particle-bed particle) 0.5 (SO-bed particle) Poisson Ratio No unit 0.25 0.25 Gravity Acceleration m/s 2 9.81 9.81 To model the physical phenomenon, the larger the number of simulated bed particles the better. Yet, this number os constrained by running time and memory resources, partly mitigated by parallel processing on high performance supercomputers. We used the supercomputer, Sunway-TaihuLight, in Wuxi. The service consists of a Sunway26010 massively parallel processor (MPP) distributed memory system. It has 128 main memory per 24 core processor, using Non-Uniform Memory Architecture (NUMA). In our simulation, we specified processors with a factorization P = 6 × 6 × 6, i.e., 6 processors in each of the x, y, and z directions.

The Simulated System
The simulations took place in a container of dimensions 0.5 m × 1.0 m × 1.3 m. We used periodic boundary conditions in the x and y directions, with the dimensions in these direction deemed sufficiently large to avoid boundary effects, as well as interactions of the SO with its periodic image. The bottom of the container was fixed at z = −0.5 m. Earth gravity was in the −z direction. We generated bed particles with diameter d 0 = 0.02 m and density ρ 0 = 2500 kg/m 3 at the top of the container and let them fall under gravity and fill the box up to z = 0.0 m. At this stage, a SO of diameter d i = 0.12 m and density ρ i = 2500 kg/m 3 was placed as close as possible to position (0, 0, 0). Then, more bed grains were poured into the box up to the level of z = 0.5 m. In total, a bed of 8 * 10 5 particles was generated, filling a volume of 0.5 m × 1.0 m × 1.0 m, with the SO buried in its center ( Figure 1). Before starting each formal simulation, we ran an initial equilibration simulation to stabilize the system structure. The system was considered in equilibrium when the kinetic energy dropped to 10 −4 of its initial value. Starting from the equilibrated initial state, we carried out two sets of simulations. In the first, the SO was driven by excitation in a composite mode of oscillation and constant speed velocity rotation. The oscillation is simulated by a periodic oscillation in the y direction, y = Asin(2π f t), with amplitude A = 0.03 m and frequency f = 3 Hz, using the values in [35]. We set the SO rotations around x-axis, y-axis, z-axis, and an xz-axis at an oscillating angular velocity of amplitude ω = 80 rad/s, using the value in [36], which was shown to lead to a significant rise of SOs. In the second set of simulations, the SO was driven by both directional and rotational oscillations, y = A sin(2π f t) and ω = Ω sin(2π f t), respectively. We varied the three parameters, A, Ω, and f in a series of DEM simulations, described in Section 2.1, and analyzed their effect on the motion of the SO. The time step in all the simulations was 10 −5 s. Altogether, we ran 58 simulations, each with different parameters, 56 of which for 4 × 10 5 time steps, corresponding to = 4 s physical time, and 2 simulations for 5 × 10 5 time steps, corresponding to = 5 s physical time. One simulation of 4 s physical time took about 2 h to run in CPU time.

Coupling Oscillatory Translation and Constant Rotation Speed
We set the oscillation amplitude and frequency to A = 0.03 m and f = 3 Hz, and the rotation angular velocity to ω = 80 rad/s, with ω > 0 denoting henceforth a counterclockwise rotation of the SO. We ran five types of simulations: only vibration, vibration and rotation about the x-, y-, z-, or xz-axis. In each run, we monitored the dynamics and tracked the trajectory of the SO. Observing that its movement in the y direction is relatively small, we focus on the displacement in x and z directions.
Our observations, shown in Figure 2, can be summarized as follows.
• The displacement in the x direction oscillates strongly at frequency f whenever ω has a component about the z-axis, while otherwise it is relatively smooth. • The x displacement is vanishingly small except when the rotation is about the y-axis, in which case the SO moves steadily in the −x direction. • The displacement in the z direction oscillates strongly at frequency f when ω has a component about the x-axis, otherwise it is relatively smooth. • The displacement in the z direction is always upward, with rotation enhancing the rising rate. To understand the dynamics under coupled vibration and rotation, consider, for example, rotation about the x-axis. As the SO moves in the +y direction, the pressure from the medium on its forward moving (+y) surface is larger than the pressure on the rear surface. The constant counter-clockwise rotation around the x-axis then induces a friction force in the −z-direction on the forward moving surface, as illustrated in Figure 3, and a lower friction force in the +z direction on the rear surface. Therefore, the combined action of vibration and rotation leads to an overall downward force during this one-quarter of the cycle. On changing direction, the pressure is higher on the surface facing the −y direction and, since the rotation direction is constant, the overall friction force is upwards. Thus, the friction force oscillates between up and down with the translational motion. Nevertheless, the displacement in the +z direction is larger than in the −z direction because the layer supporting the SO is slightly less fluid than the medium above it, being deeper in the medium. This leads to an overall rise with time, consistently with the displacement marked in the red line in Figure 2b. Since the displacements in the x direction are negligible compared to those in the y − z plane, shown in Figure 2, we focus in the following on the rising and sinking in the z direction.

Coupling Oscillatory Translation and Oscillatory Rotation Speed
The above understanding of the oscillations in the z displacement during one cycle suggests changing the rotation direction in phase with the translational one, with Ω the amplitude of the rotational angular velocity and f its oscillating frequency, which was the same as that of the translational vibration in the y-direction. The idea was that the direction of the displacement, would then be either always negative or always positive, resulting in a steady rise or sink during the cycle. The results for the z displacement of the SO for A = 0.03 m, Ω = ±80 rad/s, and f = 3 Hz, shown in Figure 4, confirmed the above expectation-the SO sank steadily with positive Ω and rose steadily with negative Ω.
We tracked the trajectories and show in Figure 5 examples of slices of the system in the y − z plane through the SO at times 0 s, 2.5 s, and 5 s for both the rising and sinking. In Tables 2 and 3 we summarize the observed dynamics during a cycle for Ω = −80 rad/s.  The motion in the z direction, both up and down, is due to the imbalance in the friction forces on the SO induced by the rotation. When Ω = −80 rad/s, the SO is subject to a steady upward friction force, given in Tables 2 and 3, and consequently rises steadily, as shown in Figure 4 (black line). When Ω = 80 rad/s, the friction force on the SO is downward, leading to steady sinking, also shown in Figure 4 (red line). As can be observed in Figure 4, the steady rising and sinking were accompanied by fluctuations that were much smaller than those shown in Figure 2. These fluctuations, which appeared also at frequency f , were caused by a combined effect of particles falling into the cavity, which the moving SO left behind, and of the stress release in the compressed medium ahead of the SO as it changes direction. Thus, the in-phase coupling of the rotation to the vibration eliminates the large oscillations, observed the constant rotation mode, and the steady rise and sink are faster in this mode. It follows that this mode of locomotion is a more efficient technique for desert robots and most likely for burrowing desert animals as well.
To determine the conditions for the fastest locomotion, we turn next to investigate the effects of A, Ω, and f on the rising or sinking rates in this mode. Table 2. The corresponding values of vibration displacement and rotation velocity at one period.

The Effect of Rotation Amplitude, Ω
We fixed the vibration amplitude at A = 0.03 m, the frequency at f = 3 Hz, and varied the angular velocity amplitude Ω = ±10 × 2 n rad/s, with n = 0, 1, 2, ..., 7. Plots of the SO displacement vs. time for negative and positive Ω are shown in Figure 6a,b. The mean rate of displacement, U z = d(∆z)/dt depends sensitively on the rotation angular velocity. Figure 6c,d summarize our following observations. (i) As established in Section 3.2, the SO rose and sank for Ω negative or positive, respectively. (ii) The rising and sinking rates increased with | Ω |. (iii) When | Ω | exceeded 320 rad/s, the SO either broke the top surface or stabilized very close to the bottom of the container. The relationship between the rising and sinking rates and Ω can be fitted well by where the C 1 and C 2 are different constants for the rising and sinking processes. These fits are shown in Figure 6c,d.

The Effect of the Translational Amplitude, A
To test the effect of the translational amplitude on the vertical displacement, we carried out a series of simulations in which we fixed rotation angular velocity at Ω = ±80 rad/s and the frequency at f = 3 HzA, and varied A from 0.01 m to 0.10 m. The direction of ∆z for negative and positive Ω, respectively, remained unaffected, but the rates of rising and sinking did vary, as shown in Figure 7a,b. In Figure 7c,d, we show that the mean rising rate increased linearly with A until the SO gets close to the top surface. In contrast, the sinking rate was linear only down to about A > −0.06 m. The reason for plateauing at this specific depth is probably because of the increased solidity of the supporting layer as a result of both the increased depth and the proximity to the bottom layer. With the exception of this effect, the rising rate and the amplitude appears to be linear:  Figure 8a,b and the rates are shown in Figure 8c,d. In this case, the rates display a richer behavior then before. While the SO sinks faster with frequency, there appears to exist a critical frequency, f c = 5 Hz, up to which the rising rate increases and above which it decreases. We understand this latter behavior as follows. When f < f c , the rising mechanism is as described above, namely, because of upwards force generated by friction. This force persists throughout every stroke. The temporal z displacement is determined only by this force and the amplitude A and is therefore constant at each stroke. As a result, the rising rate is proportional to the the number of strokes per unit time, i.e., frequency.
The same mechanism persists also when f > f c , but in this regime, the high kinetic energy of the SO starts to de-solidofy the layer of bed particles supporting it and it also sinks slightly every stroke. These two mechanisms compete and while the combined effect is still upward, the rate of rising is reduced. The higher the frequency the higher the SO's kinetic energy and the less solid the supporting layer.
The sinking rate with frequency, seen in Figure 8d, is also the result of the same two mechanisms. In this case, the friction force acts downwards, forcing the SO to sink against its supporting layer. At low frequencies, this layer is relatively solid and it becomes increasingly difficult to fluidize it with depth. As the frequency increases there is a combined effect of the force downward and the partial de-solidifcation of the support layer. However, in our system there is a third effect -the approach to the rigid bottom boundary. The latter prevents de-solidification of the support layer, thus slowing the rate of sinking, as observed in Figure 8d.

Conclusions and Future Research
Inspired by the biological experiments of sand lizard diving into desert sand and Escherichia coli propelling themselves forward by rotating helically in a viscous fluid, we used DEM numerical simulations to design a rising-sinking device in granular materials, driven by a combination of vibration and rotation. The simulations allowed us to test the vertical motion of the SO in response to periodic same-frequency rotational and translational transverse oscillations. The transverse oscillation was in the y direction, y = Asin(2π f t), and the rotation was around the x axis, ω = Ωcos(2π f t). We have found that the SO rose to the bed surface when Ω < 0 and sank to approximately Z = −0.3 m when Ω > 0. The saturation in the sinking is likely an effect of the bottom boundary. To understand the effects of the key SO actuation parameters: A, f , and Ω, on its rising and sinking rates, we varied these parameters systematically. Our results show that the rising rate is linearly proportional to the A and Ω and that there is an optimal operation region for f . These findings can inform locomotion biology in understanding how animals appear to move effectively across a diversity of complex substrates. This work also has the potential to facilitate the design, control and development of a new type of bio-physically inspired robots on or within complex seabed or desert terrains.
This research can be extended in several directions, some of which we intend to explore. The range of motion and material parameters can be broadened, including the density, size and friction coefficient, so as to mimic living organisms, which often operate optimally in variable environments [7]. For example, the cross-sectional shape of the sandfish has been shown to aid rapid burial into granular media [37]. This also suggests that it would be interesting to explore the effect of this morphology, combined with head tapering on performance. Other intruder shapes should also be investigated. Another interesting direction would be to simulate non-spherical, e.g., polyhedral, bed particles, to understand the effect on the SO dynamics. We believe that in such environments the SO vertical motion would be slower, but this remains to be investigated. Another important direction is theoretical-a mathematical predictive model for this phenomenon does not yet exist. Such a theoretical model could be then tested quantitatively against our numerical results.