A novel slithering locomotion mechanism for a snake-like soft robot

J. Mech. Phys. Solids 99 (2017) 304–320

Contents lists available at ScienceDirect

Journal of the Mechanics and Physics of Solids journal homepage: www.elsevier.com/locate/jmps

A novel slithering locomotion mechanism for a snake-like soft robot ⁎

Yunteng Caoa, Yilun Liub, , Youlong Chena, Liangliang Zhua, Yuan Yana, Xi Chenc,

MARK ⁎

a

International Center for Applied Mechanics, State Key Laboratory for Strength and Vibration of Mechanical Structures, School of Aerospace Engineering, Xi’an Jiaotong University, Xi’an 710049, China State Key Laboratory for Strength and Vibration of Mechanical Structures, School of Aerospace Engineering, Xi’an Jiaotong University, Xi’an 710049, China c Columbia Nanomechanics Research Center, Department of Earth and Environmental Engineering, Columbia University, New York, NY 10027, USA b

A R T I C L E I N F O

ABSTRACT

Keywords: Slithering locomotion Rotational velocity field Traveling wave Poisson's ratio Snake-like soft robot

A novel mechanism for slithering locomotion of a snake-like soft robot is presented. A rectangular beam with an isotropic coefficient of friction of its contact surface with the flat ground can move forward or backward when actuated by a periodic traveling sinusoidal wave. The Poisson's ratio of the beam plays an important role in the slithering locomotion speed and direction, particularly when it is negative. A theoretical model is proposed to elucidate the slithering locomotion mechanism, which is analogous to the rolling of a wheel on ground. There are two key factors of slithering locomotion: a rotational velocity field and a corresponding local contact region between the beam and ground. During wriggling motion of the rectangular beam, a rotational velocity field is observed near the maximum curvature point of the beam. If the beam has a negative Poisson's ratio, the axial tension will cause a lateral expansion so that the contact region between the beam and ground is located at the outer edge of the maximum curvature (the largest lateral expansion point). The direction of the beam's velocity at this outer edge is usually opposite to the traveling wave direction, so the friction force propels the beam in the direction of the traveling wave. A similar scenario is found for the relatively large amplitude of wriggling motion when the beam's Poisson's ratio is positive. Finite element method (FEM) simulation was conducted to verify the slithering locomotion mechanism, and good agreement was found between the FEM simulation results and theoretical predictions. The insights obtained here present a simple, novel and straightforward mechanism for slithering locomotion and are helpful for future designs of snake-like soft robots.

1. Introduction Recently, soft robots have attracted significant interest of both scientists and engineers. Due to their softness and body compliance, some of the limitations of traditional hard robots may be overcome; examples of these limitations include instability when moving in complex environments, limited motion constrained by hard actuators and structures (e.g., metal rods, mechanical joints and electric motors) and difficulty in manipulating unknown or soft objects (Rus and Tolley, 2015). Therefore, soft robots hold great potential for applications in search-and-rescue operations (Tolley et al., 2014), environmental monitoring (Muth et al., 2014;

⁎

Corresponding authors. E-mail addresses: [email protected] (Y. Liu), [email protected] (X. Chen).

http://dx.doi.org/10.1016/j.jmps.2016.11.019 Received 22 July 2016; Received in revised form 24 November 2016; Accepted 30 November 2016 Available online 08 December 2016 0022-5096/ © 2016 Elsevier Ltd. All rights reserved.

J. Mech. Phys. Solids 99 (2017) 304–320

Y. Cao et al.

Vogt et al., 2013), soft-robotic manipulators (Mazzolai et al., 2012) and medical and wearable applications (Park et al., 2014a, 2014b; Polygerinos et al., 2015), to name a few. Inspired by biological systems, many designs of soft robots have been proposed, such as snake-inspired locomotion (Onal and Rus, 2013), caterpillar-inspired locomotion (DeSimone et al., 2015; Lin et al., 2011, 2013), a multi-gait quadruped (Shepherd et al., 2011; Tolley et al., 2014), worm-inspired locomotion (Seok et al., 2010) and jellyfishinspired swimming (Nawroth et al., 2012). Many of these previous studies focused on the design and optimization of the gaits of soft robots to make them easy to fabricate, actuate and control and able to adapt to severe and complex environments. Among limbless animals, the slithering locomotion of snakes is fast and can adapt to complex terrestrial and aquatic environments. Additionally, slithering locomotion is easy to activate; i.e., a slender beam can be activated by a simple periodic traveling wave. Therefore, slithering locomotion underpins the construction of snake-like soft robots. The slithering locomotion mechanism of snakes has long been studied in the literature. In general, the slithering locomotion is represented by the motion of the central line of the contact surface between the snake and ground (Alben, 2013; Goldman and Hu, 2010; Hu et al., 2009). The propulsive force is usually derived from lateral pushing against rocks and branches (Gray, 1946; Gray and Lissmann, 1950; Guo and Mahadevan, 2008; Jayne, 1986; Mahadevan et al., 2004; Mosauer, 1932), and the friction of the snake on the ground (Gray and Lissmann, 1950; Hazel et al., 1999; Hirose, 1993; Hu et al., 2009). The coefficient of friction between the belly snakeskin and the flat ground is anisotropic (Gray and Lissmann, 1950; Hu et al., 2009); that is, the coefficient is largest in the transverse directions of the snake and smallest in the forward direction. Therefore, a net propulsive force is generated by the wriggling motion of the snake. Inspired by snake-slithering locomotion, several snake-like robots have been proposed, including both a snake-like hard robot (Hirose, 1993; Liljebäck et al., 2011; Ma, 2001; Ute and Ono, 2002; Wright et al., 2007) and a snakelike soft robot (Onal and Rus, 2013). The snake-like hard robots consist of several rigid blocks connected by joints and springs; the wriggling motion is actuated by the rotation of the joints. For the snake-like soft robot, the whole body is a soft actuator that actuates the wriggling motion. However, in all of these previous snake-like robots, the propulsive force is generated by the friction anisotropy; this requires passive wheels to sustain the difference in the coefficients of friction between the lateral and forward directions. A major disadvantage is that the rigid encompassing of passive wheels defeats the purpose of a true soft robot. From the soft robotic perspective, a simpler material design is desired. In addition, most previous studies of snake-like slithering locomotion treated the motion as simply that of the central line of the contact surface. However, a rectangular beam actuated by a periodic traveling sinusoidal wave that defines the wriggling motion has a complex velocity distribution. A rotational velocity field is observed near the maximum curvature point of the beam, as shown in Fig. 1. A snake-like robot can also be driven by the well-known rolling mechanism of a wheel on ground if the contact region between the robot and ground is always located at the rotational velocity field region (the maximum curvature point). However, during wriggling motion, the position of the maximum curvature point of the beam changes over time; hence, the difficulty is how to make the contact region between the beam and ground change over time to follow the maximum curvature point of the wriggling motion. New theoretical perspectives need to be established to better control the slithering locomotion. In this work, a new slithering locomotion mechanism is proposed utilizing the rotational velocity field of wriggling motion to drive a snake-like soft robot. This mechanism does not require friction anisotropy or lateral pushing against rocks and branches (as for generating propulsive force in conventional models). Additionally, this type of soft robot has a very simple material structure and is easy to actuate with a simple periodic traveling wave. Thus, the new mechanism enables many potential applications in soft robots. The organization of the paper is as follows. In Section 2, a theoretical model is proposed to elucidate the slithering locomotion mechanism. Based on this mechanism, a new type of snake-like soft robot is developed and an analytical solution of the slithering locomotion speed is obtained. In Section 3, finite element method (FEM) simulations are conducted to verify the slithering locomotion mechanism; good agreement is found between the FEM results and theoretical predictions. The effects on the slithering locomotion speed of a snake-like soft robot of its mechanical properties and geometry, of the loading method and of uneven ground are systematically studied. Two strategies for fabrication of a snake-like soft robot are illustrated in Section 4, followed with discussions of the efficiency of the slithering locomotion in Section 5. Conclusions are given in Section 6.

Fig. 1. Schematic diagram of the wriggling motion of a beam. The upper figures are two snapshots of the beam at two instants of time, where the arrows indicate the velocity of the particles on the beam during wriggling motion. The bottom figures give the velocity field of the beam near the maximum curvature points; the rotational velocity field is observed near maximum curvature points.

305

J. Mech. Phys. Solids 99 (2017) 304–320

Y. Cao et al.

Fig. 2. Schematic diagram of the beam position. The tangent angle and unit vectors tangent and normal to the central line are shown.

2. Slithering locomotion mechanism As previously mentioned, there are two key factors in slithering locomotion: a rotational velocity field and a corresponding local contact region between the object and surface. In this section a theoretical model is proposed to elucidate the slithering locomotion mechanism. First, we derive the velocity field of a rectangular beam actuated by a periodic traveling sinusoidal wave (wriggling motion). Then, the Poisson effect of the beam is utilized to modulate the contact region between the beam and ground. Based on this theoretical model, an analytical solution of slithering locomotion speed is obtained that is compared to the FEM simulation results in Section 3. 2.1. Kinematics of the slithering locomotion The slithering locomotion of a beam can be divided into wriggling motion and locomotion. The displacement of the beam u can be represented as

⎡ ux (s, n , t ) ⎤ ⎡ uCx (t ) + u wx (s, n , t )⎤ u=⎢ ⎥, ⎥=⎢ ⎣ u y (s, n , t )⎦ ⎣ uCy (t ) + u wy (s, n , t ) ⎦

(1)

where ux and uy are the displacements of the beam in the x and y directions, uCx and uCy represent the displacements contributed by the motion of the center of mass of the beam, i.e., the locomotion of the beam, uwx and uwy represent the displacements contributed by the wriggling motion of the beam, t is time and (s, n) is the corresponding coordinates in a curvilinear coordinate l (shown in Fig. 2) following the motion of the beam. s denotes the arc length of the central line of the bottom surface from system ^s -n the left end, n is the distance perpendicular to the central line and θ is the tangent angle of central line. Here, out-of-plane displacement (perpendicular to the ground) is ignored, and the displacement components ux, uy, uwx and uwy are assumed uniform in the height (out-of-plane) direction. The focus is therefore on the displacement and velocity distribution, especially the wriggling motion on the bottom surface of the beam. During locomotion, the acceleration of the beam remains much smaller than the acceleration of gravity g at all times. Therefore, inertia is neglected throughout this paper and the motion of the beam is therefore governed by the condition that the resultant forces and moment exerted on the beam (namely, weight as well as friction and reaction forces from the support) are zero

∬S

fdsdn = 0

∬S

r × fdsdn = 0.

(2)

and (3)

Here, the wriggling motion of the beam is actuated by a periodic sinusoidal traveling wave; i.e., the curvature of the central line is s = s / L , T denotes the period, L κ = −κ0 sin 2π (t − s ), where the dimensionless time t and arc length s are defined as t = t / T , the length of the central line, and κ0 the maximum curvature applied on the beam. The velocity field of wriggling motion on the bottom surface of the beam can therefore be represented as (The detailed derivation can be found in Appendix section.)

⎡

2

∂x (s , n, t ) ∂t

=

L⎢ T⎢

( κ4πL )

⎤

⎥ ( κ2πL cos 2πt ) ⎥ κ L ⎢⎣−[nκ0 L sin 2π (t − s ) + 1]cos ( 2π cos 2π (t − s ) )⎥⎦ ⎡κ L ⎤ κ L cos 2πt − sin ( 2π cos 2πt ) ⎥ ∂y (s , n, t ) L ⎢ 2π Vyw = = ⎥. ∂t T⎢ κ L +[ nκ L sin 2 π ( t − s ) + 1]sin cos 2 π ( t − s ) ⎢⎣ 0 ( 2π )⎥⎦

Vxw =

0

cos 4πt + cos

0

0

0

0

0

(4)

Based on Eq. (4), the velocity fields on the bottom surface of the beam at four different instants, t = 0 , 0.25, 0.5, and 0.75, are plotted in Fig. 3(a). For comparison, the velocity fields obtained from FEM simulation at the same times are also shown in Fig. 3(b). Here, to simulate the wriggling motion alone, the coefficient of friction between the beam and ground was set as zero, so that the center of mass velocity of the beam is zero. Additionally, due to the lack of friction force, the angular velocity of the beam is also zero. The theoretical predictions (Eq. (4)) agree well with the FEM simulation results, except at the ends of the beam due to difference of boundary conditions. In the theoretical model periodic boundary conditions are applied to the two ends of the beam during wriggling motion, while the two ends are free in the FEM simulation, which means the configurations at the two ends cannot match each other 306

J. Mech. Phys. Solids 99 (2017) 304–320

Y. Cao et al.

Fig. 3. Velocity field on the bottom surface of the beam during wriggling motion obtained from a theoretical prediction (Eq. (4)) (a) and from FEM simulation (b). The four snapshots in the top to bottom panels correspond to four different instants, t = 0 , 0.25, 0.5 and 0.75. The instantaneous velocity centers are marked in red points in (a). Each arrow's direction and length indicate the direction and amplitude of the velocity, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

completely. The effects of the boundary conditions and loading method on the slithering locomotion of a snake-like soft robot will be discussed in Section 3. In both the theoretical and FEM numerical results, a rotational velocity field is observed at the maximum curvature region of the beam ( sin 2π (t − s ) = 1). Due to the different boundary conditions applied in the theoretical model and FEM simulation, the locations of the rotational velocity field obtained from the theoretical predictions and from FEM simulation results are somewhat different, especially for t = 0.25 and t = 0.75, as shown in Fig. 3. Inspired by the rolling mechanism of a wheel on ground, if the contact region between the beam and ground is located at the rotational velocity field, the beam can also be driven forward. However, as the rotational velocity field relocates at different points (always the maximum curvature point) of the beam during wriggling motion, the difficulty is how to make the contact region between the beam and ground change over time to follow the maximum curvature point. In the following subsection, the Poisson effect of the beam is considered to modulate the contact region between the beam and ground. 2.2. Modulating the contact region between the beam and ground through the Poisson effect Due to the Poisson effect, the extension or contraction of the beam in the axial direction (x direction) will generate deformations in the transverse directions (y and z direction). Top and side views of the beam during wriggling motion are shown in Fig. 4(a) for a

Negative Poisson’s ratio ν= 0.3

(a)

N

Positive Poisson’s ratio ν=0.3

(b)

N

M

M

M

y

N

z x

y

M

z x

N

z

z

y x

y x

(c)

(d)

z

N

M

M

z

N

(e)

N

M

M

N

x y

x y

(f) y

y

x

x

Fig. 4. The contact regions between the beam and ground during wriggling motion. (a) and (b) are the top and side views of the beam during wriggling motion for a negative and positive Poisson's ratio, respectively. Red color represents positive normal strain in the z direction and blue represents negative normal strain. Note that for the side view figures, the z direction deformation is enlarged 10 times for visibility. (c) and (d) schematically illustrate the deformation of the cross section of the beam at the marked lines M-N and M′-N′. (e) and (f) show the distribution of contact pressure. The red circles indicate the contact areas between the beam and ground. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

307

J. Mech. Phys. Solids 99 (2017) 304–320

Y. Cao et al.

negative Poisson's ratio and in Fig. 4(b) for a positive Poisson's ratio, respectively; the red color represents positive normal strain in the z direction and blue color represents negative normal strain in the z direction. The cross section of the beam during wriggling motion becomes trapezoidal, as schematically illustrated in Fig. 4(c) and (d). As the axial deformation has the largest value at the maximum curvature point, for a negative Poisson's ratio, the contact regions between the beam and ground are always located at the outer edge of the maximum curvature point, i.e., points N and N′ in Fig. 4(c), while for a positive Poisson's ratio the contact regions between the beam and ground are located at the inner edge of the maximum curvature point, i.e., points M and M′ in Fig. 4(d). Through FEM simulation, the distribution of contact pressure of the beam on the ground is shown in Figs. 4(e) and (f) for negative and positive Poisson's ratios, respectively. Indeed, the contact region between the beam and ground is not a geometric point but a small area near points N and N′ for negative Poisson's ratio, or M and M′ for positive Poisson's ratio. In the next section we will discuss the influence of the contact area on the slithering locomotion speed of the beam. Due to the Poisson effect, the contact region between the beam and ground is always located in the region of the rotational velocity field, so that a snake-like soft robot can roll on the contact point; this is analogous to the rolling of a wheel on ground. Recalling the velocity field of the beam during wriggling motion (Eq. (4)), and substituting the relations sin 2π (t − s ) = 1 and n = d/2 into Eq. (4), the velocity at the outer edge of the maximum curvature (points N and N′ in Fig. 4) Vxw,o is

Vxw, o =

2 ⎤ ⎞ ⎛κ L L ⎡ ⎛ κ0 L ⎞ d ⎢⎜ ⎟ cos 4πt + cos ⎜ 0 cos 2πt ⎟ − κ0 − 1⎥ , ⎥⎦ ⎢ ⎠ ⎝ ⎠ ⎝ T ⎣ 4π 2π 2

(5)

and letting n =-d/2, the velocity at the inner edge (the point M and M′ in Fig. 4) Vxw,i is

Vxw, i =

2 ⎤ ⎞ ⎛κ L L ⎡ ⎛ κ0 L ⎞ d ⎢⎜ ⎟ cos 4πt + cos ⎜ 0 cos 2πt ⎟ + κ 0 − 1⎥ , ⎥⎦ ⎠ ⎝ 2π T ⎢⎣ ⎝ 4π ⎠ 2

(6)

where d is the width of the beam. The variation of Vxw,o and Vxw, i with time t results from the approximation solutions of θ0(t), x0(t) and y0(t) in Eq. (A9), as only the leading items of the Taylor series expansion in Eq. (A6) are reserved. In practical slithering locomotion the amplitude of the wriggling motion is relatively small, i.e., κ0L < 2π. Then, by Taylor series expansion of κ L cos( 20π cos 2πt ) to the second term, the velocity Vxw,o and Vxw, i can be expressed as 2 L ⎡⎛ κ L ⎞ d⎤ Vxw, o = − ⎢ ⎜ 0 ⎟ + κ0 ⎥ T ⎢⎣ ⎝ 4π ⎠ 2 ⎥⎦

and

Vxw, i =

⎛ κ L ⎞2 ⎤ L⎡ d ⎢κ 0 − ⎜ 0 ⎟ ⎥ , ⎝ 4π ⎠ ⎥⎦ T ⎢⎣ 2

(7)

which are independent of time t . The velocity at the outer edge of the maximum curvature is negative, while at the inner edge it changes from positive to negative as the maximum curvature κ0 increases; this will be further validated through FEM simulation in Section 3. Here, the slithering locomotion is the superposition of the wriggling motion and the translation of the beam, so according to Eq. (1), the slithering locomotion velocities at the outer edge Vx,o and inner edge Vx, i of the maximum curvature are

Vx, o = Vcom + Vxw, o

and

Vx, i = Vcom + Vxw, i,

(8)

where Vcom=∂uCx/∂t is the center of mass velocity of the beam in the x direction. For the negative Poisson's ratio, the contact region is at the outer edge of the maximum curvature. If the velocity at the contact point is negative (Vx,o < 0, i.e., Vcom < -Vxw,o), the friction force is in the positive x direction, which accelerates the center of mass motion of the beam in the positive x direction (the acceleration stage in Fig. 6(b)). For a positive Poisson's ratio, if the velocity at the inner edge of the maximum curvature is positive (Vx, i > 0, i.e., Vcom > -Vxw, i), the friction force is in the negative x direction, which accelerates the center of mass motion in the negative x direction. Therefore, steady locomotion is obtained when Vx,o=0 for a negative Poisson's ratio and Vx, i=0 for a positive Poisson's ratio, so that the net friction force is zero, i.e. the resultant forces and moment equilibrium equations (Eq. (2) and Eq. (3)) are satisfied. Thus, the steady locomotion speed Vcom for a negative Poisson's ratio is

L ⎡ ⎛ κ0 L ⎞ d⎤ ⎢⎜ ⎟ + κ0 ⎥, T ⎢⎣ ⎝ 4π ⎠ 2 ⎥⎦ 2

Vcom = −Vxw, o =

(9)

and for a positive Poisson's ratio it is

⎛ κ L ⎞2 ⎤ L⎡ d Vcom = −Vxw, i = − ⎢κ0 − ⎜ 0 ⎟ ⎥ . ⎝ 4π ⎠ ⎥⎦ T ⎢⎣ 2

(10)

In what follows, the slithering locomotion of the beam will be explored through FEM simulation, and the effects of the mechanical properties and width of the beam, of the loading method and of uneven ground on the slithering locomotion speed will be systematically studied. 3. FEM simulation of slithering locomotion This section reports our FEM simulation study of the slithering locomotion of a rectangular beam on a flat rigid plane (ground). The periodic sinusoidal traveling wave used to drive the beam propagated in the positive x direction and gravity was in the negative z 308

J. Mech. Phys. Solids 99 (2017) 304–320

Y. Cao et al.

Fig. 5. The slithering locomotion of a rectangular beam on a flat rigid plane. The propagation of the traveling wave is in the positive x direction and gravity is in the negative z direction. The lower right inset is the bottom surface of the beam, where the three red points (Tail, Middle, Head) are trace points to track the displacement and velocity of the beam during slithering locomotion. The beam segments represent the 16 discrete temperature fields applied to actuate the wriggling motion of the beam. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

direction, as shown in Fig. 5. The coefficient of friction of the contact surface between the beam and ground was isotropic and uniform everywhere on the contact surface. Although in practice the wriggling motion of snake-like soft robots may be activated by disparate methods, e.g., shape memory polymers (Bhattacharyya et al., 2010; Monkman, 2000; Schmidt, 2006), electroactive polymers (Bar-Cohen and Zhang, 2011), dielectric elastomers (Carpi et al., 2007; Chen et al., 2014; Fox and Goulbourne, 2008; Mirfakhrai et al., 2007; Pelrine et al., 2000b, 2000c; Suo, 2010; Zhao and Suo, 2010) and pneumatic pressurization (Arroyo and DeSimone, 2014; Ciarletta and Ben Amar, 2012; Ilievski et al., 2011; Martinez et al., 2014; Mosadegh et al., 2014), without losing generality, we describe the wriggling actuation only by the bending curvature function, i.e., κ = −κ0 sin 2π (t − s ). In our FEM simulation, the bending deformation of a beam was achieved by dividing the beam into 16 segments and applying 16 discrete temperature fields to lateral surfaces of the beam, where every temperature field varied with time to simulate actuation by a periodic traveling sinusoidal wave, as shown in the inset of Fig. 5. Note that the sinusoidal curvature of the beam was just piecewise linearly approximated by the discrete temperature field; the slithering locomotion speed was somewhat influenced by the number of the discrete temperature fields. Additionally, the free boundary conditions at two ends of the beam also influenced the slithering locomotion speed. In the following section, the effects of the loading method and boundary conditions on the slithering locomotion of the beam are discussed. The FEM simulation was conducted using the commercial software ABAQUS via the explicit procedure, in which a coupled displacement-temperature 8 node solid element (C3D8T) was used to discretize the rectangular beam and a mesh convergence study was conducted to ensure appropriate mesh density. The penalty method was used to describe the frictional contact between the beam and ground. In the following FEM simulation, the size of the beam was L =1 m, height h =0.05 m, width d =0.05 m, wriggling motion period T=4 s, maximum bending curvature κ0 =4 m−1, coefficient of friction between the beam and ground μ=0.3, and the density and elastic modulus of the beam are ρ=1000 kg/m3 and E =100 MPa, respectively, unless explicitly stated elsewhere. The Poisson effect is one of the key factors affecting slithering locomotion; therefore, various Poisson ratios of the beam ranging from ν=−0.4 to ν=0.4 were considered in the FEM simulation. Further, the contact area between the beam and ground is somewhat influenced by the width, density and elastic modulus of the beam, so that the slithering locomotion speed of the beam is also affected by these factors, as discussed in this section. Although the locomotion of the beam was accelerated by the friction force between the beam and ground, the steady locomotion speed is essentially independent of the friction force. Through FEM simulation, we have found the steady locomotion speed was almost constant if the coefficient of friction between the beam and ground was larger than 0.001. 3.1. Slithering locomotion behaviors We compared the slithering locomotion of beams of different Poisson's ratios ν=−0.3 and ν=0.3. The x-direction displacements and velocities at the trace points (Tail, Middle, Head, as marked in Fig. 5) were tracked during the FEM simulation. As the motion of the beam was the superposition of the wriggling motion and translation motion of the center of mass, the locomotion displacement and speed were defined as the displacement and velocity of the center of mass of the beam to represent the translation motion of the beam. Through FEM simulation, it was found that the rectangular beam actuated by the periodic sinusoidal traveling wave can move forward or backward, depending on the Poisson's ratio of the beam material, as shown in Fig. 6. The slithering locomotion displacements and speeds with the Poisson's ratios ν=−0.3 and ν=0.3 are given in Fig. 6(a) and (b). For ν=−0.3, the locomotion direction was the same as the traveling wave propagation direction (positive x direction), while for the positive Poisson's ratio ν=0.3, the locomotion direction was in the opposite (negative x direction), which is consistent with the theoretical predictions in Section 2.2. In addition, the locomotion speed for the negative Poisson's ratio was much larger than that for the positive Poisson's ratio. For both the negative and positive Poisson's ratios, the acceleration of the beam was completed within the first period of the traveling wave. After that, the locomotion displacement of the beam increased almost linearly and the locomotion speed fluctuated around a constant, which was the steady slithering locomotion stage. The displacements and velocities at the three trace points are also given in Fig. 6(c) and (d) for a negative Poisson's ratio and in Fig. 6(e) and (f) for a positive Poisson's ratio. In the case of ν=−0.3, despite the differences in their motions, all three trace points shared an inclination to move forward accompanied by significant fluctuations, the resultant of the translational motion of the center of mass and the wriggling motion. The scenario was similar for ν=0.3 except for the direction of motion and the locomotion speed. These features of slithering locomotion for different Poisson's ratios were qualitatively predicted in Section 2, where the contact point was suggested as a key. The FEM results of the velocity field of the beam during steady slithering locomotion are shown in Fig. 7, and the corresponding slithering locomotion speed 309

J. Mech. Phys. Solids 99 (2017) 304–320

Y. Cao et al.

Fig. 6. The displacement and velocity of the slithering locomotion of a rectangular beam on a flat rigid plane. (a) and (b) are the slithering locomotion displacement and locomotion speed of the beam for the Poisson's ratio ν=−0.3 (black line) and ν=0.3 (red line), respectively. (c) and (d) show the displacement and velocity at the three trace points (marked in Fig. 3) for ν=−0.3, and (e) and (f) for ν=0.3, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

is graphed in Fig. 8. The red circles in Fig. 7 indicate the contact region between the beam and ground, i.e., the outer edge at the maximum curvature for a negative Poisson's ratio and the inner edge at the maximum curvature point for a positive Poisson's ratio, respectively. The instantaneous velocity center was located near the outer edge for ν=−0.3 and the inner edge for ν=0.3, and the velocity at the contact area was almost zero for both a negative Poisson's ratio v =−0.3 and a positive Poisson's ratio ν=0.3, which further verifies the assumption in Section 2.2, i.e., Vx,o=0 for a negative Poisson's ratio and Vx, i=0 for a positive Poisson's ratio. The slithering locomotion speeds of a rectangular beam with different Poisson's ratios ν and maximum curvature κ0 are shown in Fig. 8(a), where the geometrical parameters of the beam are constant (L=1 m, h=d=0.05 m). Here, based on Eqs. (9) and (10), a dimensionless slithering locomotion speed V is defined as V = Vcom T / L . As shown in Fig. 8(a), the trend of the slithering locomotion speed for a specific Poisson's ratio obtained through FEM simulation was similar to the theoretical prediction based on Eqs. (9) and (10); that is, the relationship between the dimensionless slithering locomotion speed V and the maximum curvature κ0 is quadratic. However, the two results are not exactly the same. To reproduce the slithering locomotion speed of the FEM simulation results, the theoretical predictions must be multiplied by a scaling factor λ. By introducing the scaling factor λ, the theoretical prediction for negative Poisson's ratio (Eq. (9)) agrees well with the FEM simulation results, as shown in Fig. 7(a) for different Poisson's ratios,ν=−0.3, ν=−0.2, ν=−0.1 with λ=0.80, λ=0.80, λ=0.69, respectively. However, for positive Poisson's ratios ν=0.3 and ν=0.1, the theoretical prediction (Eq. (10)) deviates from the FEM simulation results for the maximum curvature κ0 > 5 m−1 even if a scaling factor λ is introduced, as shown in Fig. 8(a). The deviation is caused primarily by the residual error of the Taylor series expansion

310

J. Mech. Phys. Solids 99 (2017) 304–320

Y. Cao et al.

Fig. 7. FEM results of the velocity field of the beam during steady slithering locomotion for the Poisson's ratio ν=−0.3 (a), and ν=0.3 (b). The four snapshots correspond to 4 different instants: t = 2 , 2.25, 2.5 and 2.75 from the top to bottom panels. Each arrow indicates the direction of the velocity, and its length indicates the amplitude of the velocity. The red circles indicate the contact region between the beam and ground. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 8. The slithering locomotion speed of the beam. (a) The relation between the slithering locomotion speed and the maximum curvature κ0 for different Poisson's ratios. The symbols are the FEM results and the dashed lines are the theoretical predictions of Eqs. (9) and (10). The dash dotted line is the theoretical prediction by substituting the average value of Eq. (6) into Eq. (10). (b) The relationship between the scaling factor λ and the Poisson's ratio ν. The red dashed line indicates λ=1, which is the exact theoretical prediction. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) κ L

from Eq. (6) to the second equation of Eq. (7). In such cases the higher order items of cos( 20π cos 2πt ) are no longer negligible compared to Vxw, i causing ineffectiveness of Eq. (10). If the average value of Eq. (6) is substituted into Eq. (10), the theoretical model can describe the FEM simulation results well. As the slithering locomotion speed is significant for negative Poisson's ratios, in the following analysis we focus on a beam with a negative Poisson's ratio. The scaling factor λ increases as the Poisson's ratio decreases. For ν < −0.1, the scaling factor gradually approaches a constant value of approximately 0.80, whereas for −0.1 < ν < 0 it quickly decreases to zero, as shown in Fig. 8(b). Although the theoretical model can predict a trend of the slithering locomotion speed similar to the FEM simulation results, a 311

J. Mech. Phys. Solids 99 (2017) 304–320

Y. Cao et al.

scaling factor λ is required to describe the FEM simulation results quantitatively, especially for the negative Poisson's ratio of the beam. The gap between the theoretical prediction and FEM simulation results is caused by three factors: the small contact area between the beam and ground due to the gravity deformation of the beam, the discrete temperature field applied to actuate the wriggling motion of the beam and the free boundary conditions at the two ends. The following subsections discuss the effects of these factors on the slithering locomotion speed. 3.2. Effects of the mechanical properties of the beam on the slithering locomotion speed The mechanical properties of the beam that may contribute to the slithering locomotion speed, i.e., its density ρ, elastic modulus E and Poisson's ratio ν, are not included in the theoretical predictions (Eqs. (9) and (10)). Moreover, although the Poisson's ratio is expected theoretically to affect the slithering locomotion in an abrupt way (that is, the slithering locomotion speed changes abruptly from Eq. (9) to Eq. (10) as the Poisson's ratio changes from negative to positive), the slithering locomotion speed actually varies smoothly with the Poisson's ratio, as shown in Fig. 8(b). Because the mechanism here deals more with kinematics than dynamics, we analyzed the effects of the mechanical properties qualitatively and verified our analysis via FEM simulation. Eqs. (9) and (10) are based on the assumption that the contact region between the beam and ground is a geometrical point contact. Due to the gravitational deformation of the beam, the contact actually is a small area situated around the outer edge of the maximum curvature point for a negative Poisson's ratio and around the inner edge for a positive Poisson's ratio, as shown in Figs. 4(e) and (f). The instantaneous velocity center of the slithering locomotion hence deviates from the edge of the maximum curvature point. Therefore, the small contact area around the outer edge should decrease the slithering locomotion speed and the theoretical prediction (Eq. (9)) overestimates the slithering locomotion speed for a negative Poisson's ratio. Empirically, the contact area between the beam and ground is larger if the beam is heavier or softer. For a negative Poisson's ratio, the trapezoidal deformation of the cross section is more significant if the Poisson's ratio of the beam is smaller, so that the contact area between the beam and ground is smaller. Therefore, the slithering locomotion speed of the beam decreases with increasing beam density, or decreasing elastic modulus, or increasing negative Poisson's ratio. A dimensionless parameter ρgh/ (E|ν|) can be defined that may describe the effect of the density ρ, elastic modulus E and Poisson's ratio ν on the slithering locomotion speed of the beam, as shown in Fig. 9. The slithering locomotion speed decreases with increasing ρgh/(E|ν|), and it is almost constant for ρgh/(E|ν|) < 2.5×10−5; that is, for the given parameters E=10 MPa, ν=−0.3 and h=0.05 m, it requires ρ < 125 kg/m3. However, even if the density increases ten folds (ρ=1250 kg/m3), the slithering locomotion speed is reduced by only half (VcomT/L=0.08). 3.3. Effect of the beam width on the slithering locomotion speed In addition to the mechanical properties of the beam, the contact area between the beam and ground is also influenced by the width of the beam. The relationship between the slithering locomotion speed and the maximum curvature κ0 for different widths of the beam is given in Fig. 10(a), which shows that the scaling factor λ increases with the width of the beam. Due to the larger lateral expansion at the outer edge of the maximum curvature point, the increasing beam width is beneficial to decreasing the contact area between the beam and ground to increase the scaling factor λ. For ρgh/(E|ν|) < 2.5×10−5, the scaling factor λ first increases, then gradually approaches a constant of 0.90; see the open symbols in Fig. 10(b). Note that even after adjusting the mechanical properties and width of the beam, there is still a gap between the theoretical prediction (Eq. (9)) and FEM simulation results, i.e., the scaling factor λ < 1. This results from the discrete temperature field and the free boundary conditions of both beam ends in FEM simulations, which cannot perfectly simulate the wriggling motion actuated by a periodic traveling sinusoidal wave. 3.4. Effects of the loading method and boundary conditions on the slithering locomotion speed In our FEM simulation, the wriggling motion of the beam was actuated by sixteen discrete temperature segments and free

Fig. 9. The relationship of the slithering locomotion speed to the mechanical properties of the beam. For various densities ρ, elastic moduli E and Poisson's ratios ν, their influence on the slithering locomotion speed can be represented by a dimensionless parameter ρgh/(E|ν|), where g is the gravitational acceleration.

312

J. Mech. Phys. Solids 99 (2017) 304–320

Y. Cao et al.

Fig. 10. (a) The relationship between the slithering locomotion speed and the maximum curvature κ0 for different beam widths with the Poisson's ratio ν=−0.3. (b) The relationship between the scaling factor λ and the beam width for the initial (open squares) and revised (solid squares) models, respectively. The red dashed line indicates λ=1. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 11. (a) The revised model used in the FEM simulation. The beam was divided into thirty two segments and thirty two discrete temperature fields were applied to actuate the wriggling motion of the beam. The periodic boundary conditions were applied to the two ends of the beam. (b) The relationship between the slithering locomotion speed and the maximum curvature κ0 for the initial model (open symbols) and revised model (solid symbols). The width of the beam was d=0.05 m.

boundary conditions were applied to the two ends of the beam. Indeed, during slithering locomotion, the curvature of the central line of the beam does not strictly follow the equation κ = −κ0 sin 2π (t − s ). The velocity field predicted by the theoretical model is therefore to some extent different from the FEM simulation results, especially at the two ends of the beam due to the free boundary conditions in the FEM simulation, as shown in Fig. 3(a) and (b). In the revised FEM model, the beam was divided into thirty-two segments and thirty-two time-varying discrete temperature segments were applied to the beam to refine the actuation of the wriggling motion of the beam, as shown in Fig. 11(a). Apparently, the curvature of the beam during slithering locomotion was closer to the formula κ = −κ0 sin 2π (t − s ) for the actuation with more discrete temperature fields. Additionally, the periodic boundary conditions were applied to the two ends of the beam to force the displacements at the two ends to be equal. The slithering locomotion speed obtained from the revised model was larger than that from the initial model. Meanwhile, the scaling factor λ increased from

313

J. Mech. Phys. Solids 99 (2017) 304–320

Y. Cao et al.

Fig. 12. The slithering locomotion of the beam on uneven ground. (a) The sinusoidal wave-like surface in the x direction and (b) the sinusoid wave-like surface in both directions. (c) The ground with a raised embankment and (d) with a groove. The dashed line in each figure indicates the direction of the slithering locomotion of the beam.

0.80 to 0.93 for ν=−0.3 and from 0.69 to 0.87 for ν=−0.1, as shown in Fig. 11(b). The scaling factor λ for the revised model increased further with the width of the beam and gradually approached a constant of 0.98, as shown in Fig. 10(b). Note that the remaining discrepancy between the FEM results of the revised model and theoretical predictions was mainly caused by the residual error of the Taylor series expansion from Eq. (5) to the first equation of Eq. (7) for negative Poisson's ratios. If the average value of Eq. (5) is substituted into Eq. (9), the scaling factor λ is equal to 1 for the revised model with d > 0.07 m. Therefore, the theoretical slithering locomotion speed (Eq. (9)) is the upper limit of the slithering locomotion speed of the beam with small ρgh/(E|ν|), large width d, refined actuation of the wriggling motion, and periodic boundary conditions at the two ends. However, the refinement of the actuation could increase the complexity of the snake-like soft robot, and the periodic boundary conditions are also difficult to apply to the two ends of a soft robot in practice. Additionally, the materials of the soft robot are constrained by the soft actuator (mostly polymer), so the value of the dimensionless parameter ρgh/(E|ν|) is also limited. Therefore, an optimal design of a snake-like soft robot may be achieved by balancing factors such as efficiency of slithering locomotion, easy fabrication, actuation and control, etc. 3.5. Slithering locomotion on uneven ground In the FEM simulation, the slithering locomotion of the beam occurred on flat ground. Uneven ground, however, is more common in practice. To verify the robustness of the slithering locomotion mechanism, we studied the slithering locomotion of a beam on uneven ground, as shown in Fig. 12. Four representative uneven surfaces were considered: surfaces with a sinusoidal wave in one direction and in both directions, surfaces with a raised embankment and with a sunken channel. Generally, the slithering locomotion of the beam was observed on all of these four uneven surfaces. For the surface with a sinusoidal wave in the x direction, the slithering locomotion occurred on the top of the surface wave, which means the driving mechanism was also rolling at the contact point between the beam and ground, as shown in Fig. 12(a). When on the ground with a sinusoidal wave in both directions, the beam traveled along the valley of the ground, as shown in Fig. 12(b). For this type of slithering locomotion, the driving mechanism is the combination of rolling at the contact point and pushing against the lateral barrier. Therefore, the slithering locomotion mechanism proposed in this work is compatible with other driving mechanisms, such as the lateral pushing against barriers and anisotropic friction. On the rugged surface with a raised embankment, the beam could not cross over the embankment and the slithering locomotion was along one side of it, as shown in Fig. 12(c). In contrast, on the uneven surface with a sunken channel, the beam crossed the channel easily, and the slithering locomotion was almost unaffected, as shown in Fig. 12(d). If an adaptive mechanism can be introduced in the slithering locomotion, such as the head of beam bending upward when it encounters a rise, an embankment may be easily overcome. The effects of uneven surfaces on slithering locomotion will be systematically studied in the future, and a snakelike soft robot with self-adaption to a complex environment will be explored. 4. Two strategies for fabricating a snake-like soft robot Based on the slithering locomotion mechanism mentioned above, we propose two strategies for fabricating a snake-like soft robot: a soft robot made of a dielectric elastomer actuator and a pneumatic soft robot. 314

J. Mech. Phys. Solids 99 (2017) 304–320

Y. Cao et al.

(a)

Body

(b) Inextensible layer

DE actuator

(c)

Fiber Compliant electrodes Dielectric elastomer

Fig. 13. Schematic illustration of the snake-like soft robot made with a dielectric elastomer actuator. (a) The structure of the robot. (b) A unit cell of the auxetic cellular structure. (c) Dielectric elastomer actuator.

4.1. Dielectric elastomer snake-like soft robot A snake-like soft robot with a dielectric elastomer actuator sandwich structure is schematically illustrated in Fig. 13(a). The dielectric elastomer actuator is attached to the lateral sides of the body of the soft robot. An inextensible but flexible layer is inserted in the center of the robot to serve as the neutral plane during bending deformation. The body of the robot is comprised of auxetic cellular structures. In addition to having the required property of a negative Poisson's ratio, this type of material can significantly reduce the weight of the robot due to its porous structure. It has been demonstrated that many natural and manmade cellular materials exhibit negative Poisson's ratios, ranging from macroscopic cellular structures (Masters and Evans, 1996) to micro- and nanoscale materials such as silicates (Alderson and Evans, 2002), zeolites (Grima et al., 2000), and carbon nanotube sheets (Hall et al., 2008). For example, simple cubic auxetic metamaterials (Shen et al., 2014) exhibit a negative Poisson's ratio down to −0.4 and are stable within a compressive strain larger than 0.3. A unit cell of the auxetic cellular structure is illustrated in Fig. 13(b). For a given prestrain of 0.2, the illustrated cellular structure with an 8% imperfection has an approximately constant negative Poisson's ratio for nominal strains ranging from −0.1 to 0.1, which may suit a snake-like soft robot. The actuators of soft robots are usually made of soft active materials due to their excellent deformability and flexibility: dielectric elastomers (Carpi et al., 2007; Chen et al., 2014; Fox and Goulbourne, 2008; Mirfakhrai et al., 2007; Pelrine et al., 2000b, 2000c; Suo, 2010; Zhao and Suo, 2010), electrostrictive polymer (Pelrine et al., 2000a), conducting polymer (Hara et al., 2006; Spinks et al., 2006), etc. The properties of such soft active materials are listed in Table 1 (Brochu and Pei, 2010). The dielectric elastomer actuator is discussed here as an example. It is a thin layer of dielectric elastomer sandwiched between two compliant electrodes, as shown in Fig. 13(c). Fibers are embedded in the lateral direction of the dielectric elastomer to enhance the lateral stiffness. When the dielectric elastomer actuator is subjected to a voltage through its thickness, it reduces its thickness and expands its area; i.e., expansion occurs in both the axial and vertical directions. The axial expansion of the elastomer generates the bending deformation of the soft robot because of the inextensible layer in the robot body, and its vertical expansion is consistent with the lateral deformation of the body caused by axial expansion and Poisson's ratio. Numerous authors have studied the dielectric elastomer actuator, and interested readers are directed to the previous references. Time varying voltages can be applied to the dielectric elastomer actuators in the lateral sides of the robot to actuate the required curvature field, i.e., κ = −κ0 sin 2π (t − s ), to achieve the slithering locomotion.

4.2. Pneumatic snake-like soft robot Another type of snake-like soft robot is based on a pneumatic actuator. Inspired by previous work (Mosadegh et al., 2014), the structure of the pneumatic actuator is schematically illustrated in Fig. 14. The pneumatic snake-like soft robot is connected by several pneumatic actuators that can be actuated separately to achieve the curvature field required for slithering locomotion. In general, each pneumatic actuator consists of two extensible layers and a middle inextensible but flexible layer. The extensible layer Table 1 Properties of soft actuator materials (Brochu and Pei, 2010). Type

Maximum strain (%)

Maximum stress (MPa)

Specific density

Dielectric elastomer (acrylic with prestrain) Dielectric elastomer (silicone with prestrain) Electrostrictive polymer [P(VDF–TrFE)] Shape memory polymer (polyurethane) Conducting polymer (PANI)

380 63 4.3 100 10

7.2 3 43 4 450

1 1 1.8 1 1

315

J. Mech. Phys. Solids 99 (2017) 304–320

Y. Cao et al.

Fig. 14. Schematic illustration of the pneumatic snake-like soft robot. (a) The structure of the pneumatic snake-like soft robot. (b) The inside channels of a unit cell of the pneumatic soft robot. (c) The top (cut from 1 to 1 surface in (a)) and front view (cut from 2 to 2 surface in (a)) of the deformed configuration of the pneumatic soft robot.

contains two rows of discrete chambers, upper and bottom row chambers, connected by a single channel, as shown in Fig. 14(b). When the pneumatic actuator is actuated by pressurized air, all of the chambers expand. Because the chambers in the upper and bottom rows are close to each other, the repulsion between the adjacent expanded chambers causes the actuator to expand in the axial (top view in Fig. 14(c)) and lateral directions (front view in Fig. 14(c)). The axial expansion of the actuator generates the bending deformation of the soft robot. Slithering locomotion is achieved by properly adjusting the phase differences between the expansions of adjacent actuators. Further, the close proximity of a chamber in the upper row and its neighbor in the bottom row causes them to push against each other and results in an elongation. The proposed pneumatic actuator also has an auxetic structure, which is beneficial to obtain high slithering locomotion speed. Indeed, the auxetic property is intrinsic to the nature of a pneumatic actuator because the air pressure is isotropic. Therefore, the magnitude of the negative Poisson's ratio for pneumatic actuator can be modulated over a very wide range by properly adjusting the structure of the pneumatic actuator. In our future work, the impact of the structure of a pneumatic actuator, such as the number of chambers, the wall thickness, and the chamber height, on the mechanical behaviors of the pneumatic actuator will be systematically studied to fine tune the locomotion behaviors of the robot. 5. Performance of slithering locomotion In discussing slithering locomotion performance, it is useful to define the dimensionless energy dissipation η as the ratio of the energy dissipated by friction between the slithering locomotion EL and the sliding motion ES, i.e., η = EL /ES . EL can be derived directly from FEM simulation and ES=μMgS0. μ is the coefficient of friction, Mg is the weight of the beam and S0 is the sliding distance of slithering locomotion. The energy dissipation of the slithering locomotion for different Poisson's ratios is shown in Fig. 15. For a Poisson's ratio smaller than −0.2, the energy dissipation of slithering is smaller than that of sliding. It is interesting

Fig. 15. The energy dissipation of slithering locomotion. The red dashed line indicates η = 1. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

316

J. Mech. Phys. Solids 99 (2017) 304–320

Y. Cao et al.

Fig. 16. Comparison of the movement efficiencies and deformation levels among several soft robots and organisms (Hu et al., 2009; Lin et al., 2013; Maladen et al., 2011; Quillin, 1999; Seok et al., 2010; Shepherd et al., 2011; Tolley et al., 2014).

that although the distance moved (including axial and lateral movements) of the slithering locomotion is much larger than that of the sliding motion, its frictional energy dissipation is smaller. Therefore, for energy consumption, the design of a slithering locomotion soft robot should be better in the region of ν < −0.2. We compare the movement efficiency and the deformation level of the slithering locomotion presented in this work to several other soft robots and organisms (e.g., the quadruped soft robot (Shepherd et al., 2011; Tolley et al., 2014), the caterpillar soft robot (Lin et al., 2013), milk snake (Hu et al., 2009), sandfish lizard (Goldman and Hu, 2010), etc.) as shown in Fig. 16. Here, we define movement efficiency as the ratio of the distance moved per loading cycle (the actuation of the soft robots or soft organisms is usually periodic) to the maximum dimension of the soft robot or the soft organism. The deformation level is defined as the ratio between the maximum dimension and the corresponding minimum dimension during the actuation. For example, the deformation level of slithering locomotion is defined as the ratio of the initial length of the undeformed beam and the length during slithering locomotion (the distance between the head and the tail in Fig. 5). Similar definitions are used for the snake, sandfish lizard, quadruped soft robot and caterpillar soft robot. For the earthworm or the peristaltic soft robot, the deformation level is defined as the ratio between the initial radius and the radius after worm deformation. The deformation level, an important parameter in soft robot design, is usually constrained by the deformability of the soft body, the actuators and the environments. A larger movement efficiency and smaller deformation level represent better performance for soft robots. The slithering locomotion presented in this work has a relatively large movement efficiency and small deformation level, and hence is a promising candidate for soft robots. 6. Conclusions Slithering locomotion, in many cases, is very efficient and adaptive to complex surroundings and has inspired the development of various slithering (mostly snake-like) soft robots. In this work, a novel yet simple slithering locomotion mechanism is proposed for a snake-like soft robot. A rectangular beam with an isotropic and uniform coefficient of friction between its contact surface and ground can move forward or backward when it is actuated by a periodic traveling sinusoidal wave. The mechanism of slithering locomotion is completely different from that of previous snake-like robots, which require friction anisotropy. Both FEM simulation and theoretical analysis have been conducted to explore the mechanism of the slithering locomotion. It was found that there is a rotational velocity field near the maximum curvature point of the beam. For a negative Poisson's ratio, the contact point between the beam and ground is located at the outer edge of the maximum curvature, while for a positive Poisson's ratio, the contact point is located at the inner edge of the maximum curvature. The slithering locomotion mechanism is analogous to the rolling of a wheel on ground. Based on this idea, a theoretical prediction of the slithering locomotion speed is proposed that agrees well with the FEM simulation results. The slithering locomotion presented here has the characteristics of simple structure, easy actuation and control, and high locomotion efficiency, and holds great potential for applications in snake-like soft robot. Acknowledgments Y.L. acknowledges support from the National Natural Science Foundation of China (11302163 and 11572239), and X.C. acknowledges support from the National Natural Science Foundation of China (11572238 and 11372241), ARPA-E (DEAR0000396) and AFOSR (FA9550-12-1-0159). We thank the anonymous reviewers and the editor for their comments on the manuscript.

Appendix. Velocity Field on the bottom surface of the beam actuated by wriggling motion

l (shown in Fig. 2), geometric relations among curvature, coordinates and tangent angle In the curvilinear coordinate system ^s -n

317

J. Mech. Phys. Solids 99 (2017) 304–320

Y. Cao et al.

along the central line are (A1)

∂x (s, 0, t )/∂s = cos θ , ∂y (s, 0, t )/∂s = sin θ , and ∂θ /∂s = κ ,

where x(s, 0, t) and y(s, 0, t) represent the coordinates of the central line in the reference coordinate system x-y. By integrating the curvature κ along the central line (the third equation of Eq. (A1)), the tangent angle θ is

∫0

θ (s , t ) = θ 0 (t ) +

s

κ (s1, t ) ds1

(A2)

During slithering locomotion, the beam is considered to be a Bernoulli-Euler beam; i.e., the cross section of the beam is perpendicular to its central line, which is valid for a slender beam with a shear modulus comparable to Young's modulus. By integrating the first and second equations of Eq. (A1) and using the Bernoulli-Euler beam assumption, the coordinates of a point on the bottom surface of the beam during wriggling motion can be expressed as

x ( s , n , t ) = x 0 (t ) +

∫0

y (s, n, t ) = y0 (t ) +

∫0

s

cos θ (s1, t ) ds1 − n sin θ (s, t )

(A3)

and s

sin θ (s1, t ) ds1 + n cos θ (s, t ).

(A4)

Here, θ0 (t ), x 0 (t ) and y0 (t ) are the tangent angle and coordinates at point (s, n) =(0, 0), which are determined by the boundary conditions of the wriggling motion. Due to the decomposition of the locomotion into the motion of the center of mass and the wriggling motion (Eq. (1)), the velocity of the center of mass caused by the wriggling motion is always zero. Besides, as the initial angular velocity of the beam and the resultant moment exerted on the beam are zero (Eq. (3)), the angular velocity of the beam is also always zero. Since wriggling motion of the beam is actuated by a periodic sinusoidal traveling wave; i.e., the curvature of the central line is κ = −κ0 sin 2π (t − s ), the tangent angle θ is

θ (s , t ) = −

κ0 L κ L cos 2π (t − s ) + θ0 (t ) + 0 cos 2πt . 2π 2π

(A5)

Substituting Eq. (A5) into Eqs. (A3) and (A4), the coordinates of a given point on the bottom surface of the beam can be obtained, from which the velocity field of the beam is derived also. However, it is not possible to obtain an explicit expression of x(s, 0, t) and y(s, 0, t) directly from Eqs. (A3) and (A4), so that it is also not possible to derive a closed-form solution of θ0(t), x0(t) and y0(t) from the boundary conditions. To solve this problem, the two terms cosθ(s1, t) and sinθ(s1, t) in Eqs. (A3) and (A4) are expressed by Taylor's series to the order of θ2, i.e., sin θ = θ + O (θ 3), cos θ = 1 − θ 2 /2 + O (θ 4 ). Substituting θ (s , t ) (Eq. (A5)) into the Taylor's series, the sinθ and cosθ can be represented as

sin θ = −

κ0 L 2π

cos 2π (t − s1) + θ0 (t ) +

and cos θ = 1 −

1 2

[

κ L − 20π

κ0 L 2π

cos 2πt + O (θ 3)

cos 2π (t − s1) + θ0 (t ) +

κ0 L 2π

cos 2πt

2

]

+ O (θ 4 ).

(A6)

Then, substituting Eq. (A6) into Eqs. (A3) and (A4), the coordinates x and y can be expressed as

⎫ n ⎧ κ0 L ⎨ [cos 2πt − cos 2π (t − s )] + θ0 (t ) ⎬ + Ls ⎭ 2 ⎩ 2π κ 2L3 ⎡8πs + 4πs cos 4πt + 8 cos 2πt sin 2π (t − s )⎤ − 0 3⎢ ⎥ ⎦ 64π ⎣− sin 4π (t − s ) − 3 sin 4πt

x (s , n , t ) = x 0 (t ) −

+

Lθ0 (t ) {κ0 L [−2πs cos 2πt − sin 2π (t − s ) + sin 2πt ] − 2π 2sθ0 (t )} 4π 2

and

⎡κ L ⎤ y (s , n , t ) = y0 (t ) + L ⎢ 0 2 (πs cos 2πt − cos π (s − 2t )sin πs ) + sθ0 (t ) ⎥ ⎣ 2π ⎦

+

⎤2 n n ⎡κ L − ⎢ 0 (cos 2πt − cos 2π (t − s )) + θ0 (t ) ⎥ . ⎦ 2 4 ⎣ 2π (A7)

Based on the boundary conditions of the wriggling motion, the approximate solutions of θ0(t), x0(t) and y0(t) are given

θ 0 (t ) = − x 0 (t ) =

κ0 L 2π

κ02 L3 (4π )3

cos 2πt sin 4πt , y0 (t ) =

κ 0 L2 4π 2

sin 2πt .

(A8)

Substituting Eq. (A8) into Eq. (A7), the displacement field on the bottom surface of the beam during wriggling motion is obtained, based on which the velocity field on the bottom surface is calculated as the time derivation of the displacement field,

318

J. Mech. Phys. Solids 99 (2017) 304–320

Y. Cao et al.

⎡

Vxw =

=

L⎢ T⎢

( κ4πL ) 0

cos 4πt + cos

⎤

0

0

0

Vyw =

2

⎥ ( κ2πL cos 2πt ) ⎥ κ L ⎢⎣−[nκ0 L sin 2π (t − s ) + 1]cos ( 2π cos 2π (t − s ) )⎥⎦ ⎡κ L ⎤ κ L cos 2πt − sin ( 2π cos 2πt ) ⎥ ∂y (s , n, t ) L ⎢ 2π = ⎥. ∂t T⎢ κ L +[ nκ L sin 2 π ( t − s ) + 1]sin cos 2 π ( t − s ) ⎢⎣ 0 ( 2π )⎥⎦

∂x (s , n, t ) ∂t

0

0

(A9)

References Alben, S., 2013. Optimizing snake locomotion in the plane. Proc. R. Soc. A, 469. Alderson, A., Evans, K.E., 2002. Molecular origin of auxetic behavior in tetrahedral framework silicates. Phys. Rev. Lett. 89, 225503. Arroyo, M., DeSimone, A., 2014. Shape control of active surfaces inspired by the movement of euglenids. J. Mech. Phys. Solids 62, 99–112. Bar-Cohen, Y., Zhang, Q., 2011. Electroactive polymer actuators and sensors. MRS Bull. 33, 173–181. Bhattacharyya, R., Di Leo, C., Floerkemeier, C., Sarma, S., Anand, L., 2010. RFID tag antenna based temperature sensing using shape memory polymer actuation. IEEE Sens., 2363–2368. Brochu, P., Pei, Q., 2010. Advances in dielectric elastomers for actuators and artificial muscles. Macromol. Rapid Commun. 31, 10–36. Carpi, F., Salaris, C., Rossi, D.D., 2007. Folded dielectric elastomer actuators. Smart Mater. Struct. 16, S300–S305. Chen, B., Bai, Y., Xiang, F., Sun, J.-Y., Mei Chen, Y., Wang, H., Zhou, J., Suo, Z., 2014. Stretchable and transparent hydrogels as soft conductors for dielectric elastomer actuators. J. Polym. Sci., Part B: Polym. Phys. 52, 1055–1060. Ciarletta, P., Ben Amar, M., 2012. Pattern formation in fiber-reinforced tubular tissues: folding and segmentation during epithelial growth. J. Mech. Phys. Solids 60, 525–537. DeSimone, A., Gidoni, P., Noselli, G., 2015. Liquid crystal elastomer strips as soft crawlers. J. Mech. Phys. Solids 84, 254–272. Fox, J.W., Goulbourne, N.C., 2008. On the dynamic electromechanical loading of dielectric elastomer membranes. J. Mech. Phys. Solids 56, 2669–2686. Goldman, D.I., Hu, D.L., 2010. Wiggling through the world: the mechanics of slithering locomotion depend on the surroundings. Am. Sci. 98, 314–323. Gray, J., 1946. The mechanism of locomotion in snakes. J. Exp. Biol. 23, 101–120. Gray, J., Lissmann, H.W., 1950. The kinetics of locomotion of the grass-snake. J. Exp. Biol. 26, 354–367. Grima, J.N., Jackson, R., Alderson, A., Evans, K.E., 2000. Do zeolites have negative Poisson's ratios? Adv. Mater. 12, 1912–1918. Guo, Z.V., Mahadevan, L., 2008. Limbless undulatory propulsion on land. Proc. Natl. Acad. Sci. USA 105, 3179–3184. Hall, L.J., Coluci, V.R., Galvão, D.S., Kozlov, M.E., Zhang, M., Dantas, S.O., Baughman, R.H., 2008. Sign change of Poisson's ratio for carbon nanotube sheets. Science 320, 504–507. Hara, S., Zama, T., Takashima, W., Kaneto, K., 2006. Tris(trifluoromethylsulfonyl)methide-doped polypyrrole as a conducting polymer actuator with large electrochemical strain. Synth. Metals 156, 351–355. Hazel, J., Stone, M., Grace, M.S., Tsukruk, V.V., 1999. Nanoscale design of snake skin for reptation locomotions via friction anisotropy. J. Biomech. 32, 477–484. Hirose, S., 1993. Biologically Inspired Robots: Snakelike Locomotors and Manipulators. Oxford University Press, Oxford. Hu, D.L., Nirody, J., Scott, T., Shelley, M.J., 2009. The mechanics of slithering locomotion. Proc. Natl. Acad. Sci. USA 106, 10081–10085. Ilievski, F., Mazzeo, A.D., Shepherd, R.F., Chen, X., Whitesides, G.M., 2011. Soft robotics for chemists. Angew. Chem. Int. Ed. 50, 1890–1895. Jayne, B.C., 1986. Kinematics of terrestrial snake locomition. Copeia, 915–927. Liljebäck, P., Stavdahl, Ø., Pettersen, K.Y., Gravdahl, J.T., 2011. Two new design concepts for snake robot locomotion in unstructured environments. Paladyn 1, 154–159. Lin, H.T., Leisk, G.G., Trimmer, B., 2011. GoQBot: a caterpillar-inspired soft-bodied rolling robot. Bioinspir. Biomim. 6, 026007. Lin, H.T., Leisk, G.G., Trimmer, B.A., 2013. Soft robots in space: a perspective for soft robotics. Acta Futur. 6, 69–79. Ma, S.G., 2001. Analysis of creeping locomotion of a snake-like robot. Adv. Robot. 15, 205–224. Mahadevan, L., Daniel, S., Chaudhury, M.K., 2004. Biomimetic ratcheting motion of a soft, slender, sessile gel. Proc. Natl. Acad. Sci. USA 101, 23–26. Maladen, R.D., Ding, Y., Umbanhowar, P.B., Kamor, A., Goldman, D.I., 2011. Mechanical models of sandfish locomotion reveal principles of high performance subsurface sand-swimming. J. R. Soc. Interface 8, 1332–1345. Martinez, R.V., Glavan, A.C., Keplinger, C., Oyetibo, A.I., Whitesides, G.M., 2014. Soft actuators and robots that are resistant to mechanical damage. Adv. Funct. Mater. 24, 3003–3010. Masters, I.G., Evans, K.E., 1996. Models for the elastic deformation of honeycombs. Compos. Struct. 35, 403–422. Mazzolai, B., Margheri, L., Cianchetti, M., Dario, P., Laschi, C., 2012. Soft-robotic arm inspired by the octopus: II. From artificial requirements to innovative technological solutions. Bioinspir. Biomim. 7, 025005. Mirfakhrai, T., Madden, J.D.W., Baughman, R.H., 2007. Polymer artificial muscles. Mater. Today 10, 30–38. Monkman, G.J., 2000. Advances in shape memory polymer actuation. Mechatronics 10, 489–498. Mosadegh, B., Polygerinos, P., Keplinger, C., Wennstedt, S., Shepherd, R.F., Gupta, U., Shim, J., Bertoldi, K., Walsh, C.J., Whitesides, G.M., 2014. Pneumatic networks for soft robotics that actuate rapidly. Adv. Funct. Mater. 24, 2163–2170. Mosauer, W., 1932. On the locomotion of snakes. Science 76, 583–585. Muth, J.T., Vogt, D.M., Truby, R.L., Menguc, Y., Kolesky, D.B., Wood, R.J., Lewis, J.A., 2014. Embedded 3D printing of strain sensors within highly stretchable elastomers. Adv. Mater. 26, 6307–6312. Nawroth, J.C., Lee, H., Feinberg, A.W., Ripplinger, C.M., McCain, M.L., Grosberg, A., Dabiri, J.O., Parker, K.K., 2012. A tissue-engineered jellyfish with biomimetic propulsion. Nat. Biotechnol. 30, 792–797. Onal, C.D., Rus, D., 2013. Autonomous undulatory serpentine locomotion utilizing body dynamics of a fluidic soft robot. Bioinspir. Biomim. 8, 026003. Park, Y.-L., Santos, J., Galloway, K.G., Goldfield, E.C., Wood, R.J., 2014a. A soft wearable robotic device for active knee motions using flat pneumatic artificial muscles, IEEE International Conference on Robotics and Automation, Hong Kong, 4805–4810. Park, Y.L., Chen, B.R., Perez-Arancibia, N.O., Young, D., Stirling, L., Wood, R.J., Goldfield, E.C., Nagpal, R., 2014b. Design and control of a bio-inspired soft wearable robotic device for ankle-foot rehabilitation. Bioinspir. Biomim. 9, 016007. Pelrine, R., Kornbluh, R., Joseph, J., Heydt, R., Pei, Q., Chiba, S., 2000a. High-field deformation of elastomeric dielectrics for actuators. Mater. Sci. Eng. C. 11, 89–100. Pelrine, R., Kornbluh, R., Kofod, G., 2000b. High-strain actuator materials based on dielectric elastomers. Adv. Mater. 12, 1223–1225. Pelrine, R., Kornbluh, R., Pei, Q.B., Joseph, J., 2000c. High-speed electrically actuated elastomers with strain greater than 100%. Science 287, 836–839. Polygerinos, P., Wang, Z., Galloway, K.C., Wood, R.J., Walsh, C.J., 2015. Soft robotic glove for combined assistance and at-home rehabilitation. Robot. Auton. Syst. 73, 135–143. Quillin, K.J., 1999. Kinematic scaling of locomotion by hydrostatic animals: ontogeny of peristaltic crawling by the earthworm lumbricus terrestris. J. Exp. Biol. 202 (Pt 6), 661–674.

319

J. Mech. Phys. Solids 99 (2017) 304–320

Y. Cao et al.

Rus, D., Tolley, M.T., 2015. Design, fabrication and control of soft robots. Nature 521, 467–475. Schmidt, A.M., 2006. Electromagnetic activation of shape memory polymer networks containing magnetic nanoparticles. Macromol. Rapid Commun. 27, 1168–1172. Seok, S., Onal, C.D., Wood, R., Rus, D., Kim, S., 2010. Peristaltic locomotion with antagonistic actuators in soft robotics, Proceedings 2010 IEEE International Conference on Robotics and Automation. IEEE, Anchorage, AK, 1228–1233. Shen, J., Zhou, S., Huang, X., Xie, Y.M., 2014. Simple cubic three-dimensional auxetic metamaterials. Phys. Status Solidi B 251, 1515–1522. Shepherd, R.F., Ilievski, F., Choi, W., Morin, S.A., Stokes, A.A., Mazzeo, A.D., Chen, X., Wang, M., Whitesides, G.M., 2011. Multigait soft robot. Proc. Natl. Acad. Sci. USA 108, 20400–20403. Spinks, G.M., Mottaghitalab, V., Bahrami-Samani, M., Whitten, P.G., Wallace, G.G., 2006. Carbon-nanotube-reinforced polyaniline fibers for high-strengthhighstrength artificial muscles. Adv. Mater. 18, 637–640. Suo, Z., 2010. Theory of dielectric elastomers. Acta Mech. Solid. Sin. 23, 549–578. Tolley, M.T., Shepherd, R.F., Mosadegh, B., Galloway, K.C., Wehner, M., Karpelson, M., Wood, R.J., Whitesides, G.M., 2014. A resilient, untethered soft robot. Soft Robot. 1, 213–223. Ute, J., Ono, K., 2002. Fast and efficient locomotion of a snake robot based on self-excitation principle, In: Proceedings of the 7th International Workshop on Advanced Motion Control, IEEE, 532-539. Vogt, D.M., Park, Y.-L., Wood, R.J., 2013. Design and characterization of a soft multi-axis force sensor using embedded microfluidic channels. IEEE Sens. J. 13, 4056–4064. Wright, C., Johnson, A., Peck, A., McCord, Z., Naaktgeboren, A., Gianfortoni, P., Gonzalez-Rivero, M., Hatton, R., Choset, H., 2007. Design of a modular snake robot, Proceedings 2007 IEEE/RSJ International Conference on Intelligent Robots and Systems. IEEE, San Diego, 2609–2614. Zhao, X., Suo, Z., 2010. Theory of dielectric elastomers capable of giant deformation of actuation. Phys. Rev. Lett. 104, 178302.

320