Dynamic Modeling of Underwater Snake Robot by Hybrid Rigid-Soft Actuation

Abstract

For decades, underwater vehicles have been performing underwater operations, which are critical to the development and upgrading of underwater robots. With the advancement of technology, various types of robots have been developed. The underwater robotic snake is a bioinspired addition to the family of underwater robotic vehicles. In this paper, we propose an innovative underwater snake robot actuated by rigid propulsions and soft joints, which can improve the swimming efficiency and flexibility of the robot and reduce the probability of collision leading to damage. Existing math models of robotic snakes typically incorporate only planar motion, rarely considering spatial motion. So, we formulate a complete three-dimensional dynamic model for the robotic snake, which is extended by deriving expressions for the geometric Jacobians. This modeling approach is well suited since it provides compact matrix expressions and easy implementation. We use the constant curvature method to describe the configuration of the soft joint, use the Lagrangian method to obtain its dynamic characteristics, and focus on deriving the visco-hyperelastic mechanical energy of the soft material. Next, the local dynamics of soft members are extended as a nonholonomic constraint form for modeling the snake robot. Finally, the multi-modal swimming behavior of the robot has been verified by simulations, including forward and backward rectilinear motion, yaw turning, pitch motion, and spiral rising motion. The overall results demonstrate the effectiveness and the versatility of the developed dynamic model in the prediction of the robot trajectory, position, orientation, and velocity.

1. Introduction

In the past decades, due to technological innovation, robotic underwater vehicles have increased rapidly, enabling these mechanisms to operate in deep water and harsh seabed environments [1,2]. Today, autonomous underwater vehicles (AUVs) and remotely operated vehicles (ROVs) are widely used subsea for different challenging tasks, such as inspection, surveillance, maintenance, repair, and construction, and they are extensively used in the subsea oil and gas industry [3,4,5].

In addition to the common ROVs and AUVs, the underwater robots also include various types, such as robotic fish [6,7], snake robots [8], etc. The conventional work-class ROVs can work in deep waters and are equipped with a powerful manipulator, which can perform heavy intervention tasks [9]. However, their operation is always expensive, because they need the support of surface ships, professionals for deployment, and constant supervision. In contrast, the inspection-class ROVs are much smaller, easier to deploy, and cheaper to operate [10]. However, they are not powerful enough to perform a heavy intervention. Moreover, ROVs are usually restricted by the attached tether. Some limitations related to ROVs can be solved by using AUVs, such as the need to support surface ships, continuous monitoring, and tether control. Most commercial AUVs are single torpedo-shaped vehicles, which are used for reconnaissance and surveillance, as well as long-term search missions [11]. Although some AUVs possess hovering capabilities and can perform stationary tasks, they are not particularly suitable for intervention operations. In recent years, modular underwater snake robots have developed rapidly, such as ACM [12], AmphiBot [13], and HUMRS [14]. Due to their slender and flexible body structure, they are superior in terms of accessibility and passability, so they can provide access to narrow areas that are difficult to access by other types of underwater vehicles. Furthermore, Pettersen et al. propose the concept of an underwater swimming manipulator [15], which is a promoted snake robot equipped with additional effectors, such as a stern propeller and head gripper. Three generations of products have been released up to now, namely Mamba [16], Eelume [17], and Eely500 [18], undoubtedly, this research is cutting-edge. However, there are two basic problems in the current research on snake robots. The first is that the utilization of articulated joints leads to low degrees of freedom, especially for a robot composed of a small number of modules. The second is that the whole robot is made of rigid materials, which is prone to collision and damage in the complex underwater environment.

Underwater soft robot technology has developed rapidly in recent years, such as octopus tentacle soft arm [19,20], elephant trunk arm [21,22], and deep-sea self-powered soft fish [23], which shows that soft robot has broad prospects in a marine application. This paper, inspired by the locomotion gait of biological snakes, proposes an innovative modular underwater snake robot actuated by rigid propulsions and soft joints, which can improve the swimming efficiency and flexibility of the robot and reduce the probability of collision leading to damage. Based on the above presentation of the current state of underwater robots, and the rigid-soft snake robot presented in this paper, we summarize a variety of properties from size and weight, drag force, accessibility, payload capacity, intervention, collision consequence, and maneuverability, and the author compares these performances from low to large based on our own experience as shown in Figure 1. The proposal of the rigid-soft underwater snake robot provides a novel reference and choice for the diverse development of underwater robots.

At present, the vast majority of existing studies of underwater snake robot modeling have focused on planar motions [24,25,26,27,28,29,30]. Notice that land snakes crawl on the flat surface, while water snakes always swim in three-dimensional motion. Hence, the three-dimensional modeling of an underwater robotic snake will facilitate the understanding of motions in swimming and the design of motion controllers. Several groups have preliminarily studied the three-dimensional modeling of the robotic snake in the literature [31,32,33]. However, these models are typically incomplete in predicting multimodal three-dimensional behaviors. In addition, these studies either involve only simulation or contain limited experimental verification.

In fact, underwater snake robots are a special type of underwater vehicle manipulator system (UVMS) because their body structure is integrated. Therefore, the modeling method of the snake robot in this paper is extended by deriving expressions for the geometric Jacobians, using a method inspired by literature [34]. This modeling approach is well suited since it provides compact matrix expressions, easy implementation, and flexibility that accommodates an arbitrary number of modules, with variable module length. The expressions are also valid if the robot includes soft joints. The geometric Jacobians are then utilized to derive the state-dependent rigid propulsion inertia matrix and the state-dependent propellers configuration matrix. For the soft joint part, we use the Constant Curvature method to describe its configuration [35], use the Lagrangian method to obtain its dynamic characteristics, and focus on deriving the visco-hyperelastic mechanical energy of the soft material. Next, the local dynamics of soft members are extended as a nonholonomic constraint form for modeling the multibody rigid-soft snake robot system. Finally, the multi-modal swimming behavior of the robot has been verified by simulation, including forward and backward rectilinear motion, yaw turning, pitch motion, and spiral rising motion. The overall results demonstrate the effectiveness and the versatility of the developed dynamic model in the prediction of the robot trajectory, position, orientation, and velocity.

2. Methodology

2.1. Notation

The proposed underwater snake robot is composed of rigid propulsion modules and soft joints connected in series, in which rigid propulsions are rigid cylinder shapes and can be installed with multiple groups of propellers, which can be regarded as micro AUV. The soft joints mainly consist of three hydraulic soft actuators that are evenly distributed, which can be bent in all directions by controlling properly. In practical applications, the precise manipulation of soft joints and soft actuators will be very complicated and difficult. In the derivation of the mathematical model, we first used the constant curvature model, assuming that the soft joint bends with a constant curvature when bending, has no torsion, and the section is always perpendicular to the centerline. In addition, when inputting pressure to the joint, there will be interference from a variety of uncertainties. Here, we assume that the output of joint torque is in direct proportion to the water pressure. In addition, when the soft joint is disturbed by the water flow, it will deform, so that the hydrodynamic coefficient is not a fixed value, which greatly increases the difficulty of modeling. For details, see the modeling section below. Figure 2 shows a schematic diagram of the multi-modal motion of the proposed snake robot, taking the single-joint structure as an example, it can perform various motions such as rectilinear, turning, pitching, and spiral motion. Obviously, this performance cannot be accomplished with the rigid joint with only one degree of freedom of a traditional snake robot. We will focus on these four modes in subsequent simulation studies.

The snake robot is modular essentially and can expand the multi-module structure, so in the mathematical modeling we consider that the snake robot consists of n + 1 rigid propulsions connected by n soft joints. Assuming that the rigid propulsions are standard cylindrical shape with their center of mass located at the geometric center, they are numbered from the tail forward and indicated by subscripts i ∈ [0 … n’], where i = 0, the tail rigid propulsion is considered as the base of the snake robot. In addition, each rigid propulsion module is equipped with m propellers numbered j ∈ [1, …, m], including embedded longitudinal propellers and external transverse propellers. The propeller configuration can be described with thrust direction and attack point, which is not affected by joint configuration. The propeller provides forward propulsion and hovering retention by eliminating nonzero buoyancy, water currents, and other disturbances. The soft joints have three degrees of mobility, and they are numbered with i ∈ [1 … n], such that the rigid propulsion i and i − 1 are connected by joint i.

The following reference frames will be referred to when considering the motion description of the snake robot. All frames are displayed in Figure 3 and all frames have y-axes pointing into the plane and x-axes pointing along the tangential direction such that the coordinate systems conform to the right-hand rule.

• Inertial frame F_I is fixed with the earth.

• Rigid propulsion body frames are attached to the centerline on the side closest to the tail, which is represented by F_{i ‘}.

• The rigid propulsion module of the tail is the robot base, so the frame F_{0 ‘} is defined as the base frame F_b.

• The body frames of the soft joints are fixed at the end closest to the tail, which is the same as the definition of a rigid propulsion frame, denoted by F_i.

• There is an additional end frame F_e at the robot head.

• The coordinate of the propeller is described by the frames F_t,i,j, where the subscript i and j are the rigid propulsion index and propeller index, respectively. The origin of the frame F_t,i,j is at the attack point of the propeller, and the thrust direction coincides with the x-axis.

2.2. Kinematic Analysis

Kinematics modeling is the description of the robot’s motion, without considering the forces causing the motion. Forward kinematics can be used to describe the end motion with respect to joints and base. Inverse kinematics is reversed to find the base and joint motion corresponding to the required end motion [36]. Therefore, consistent with industrial manipulators, kinematics includes the conversion between configuration space and task space and vice versa. The kinematic matrix is however invertible only in certain cases. Therefore, differential kinematics may be used to calculate the desired base and joint trajectories instead. This section focuses on the derivation of the forward kinematics and differential kinematics of the snake robot.

(1) Forward Kinematics: The snake robot is an articulated structure that consists of multiple rigid and soft bodies connected to each other, that is, a multibody system. A mathematically convenient and geometrically meaningful method is to use a homogeneous transformation matrix to represent the configuration of each part in the 3D Euclidean space and the transformations among these parts. For example, the pose of F_b with respect to F_I is given by the homogeneous transformation matrix:

(1)${\mathbf{g}}_b^I=\left[\begin{array}{cc}{\mathit{R}}_b^I& {\mathit{p}}_b^I\\ 0_{1\times 3}& 1\end{array}\right]\in SE(3)$

where ${\mathit{R}}_b^I$ ∈ SO(3) is rotation matrix describing the orientation of the base frame and ${\mathit{p}}_b^I$∈ R³ is the vector from the origin of F_I to the origin of F_b. The transformation in (1) indicates a translation by ${\mathit{p}}_b^I$ followed by a rotation given by ${\mathit{R}}_b^I$. We use pose as a general term for position and orientation.

Rigid propulsion is a rigid component that keeps straight all times, so the transformation matrix T from the arbitrary rigid propulsion body frame F_i’ to the soft joint body frame F_i+1 (i.e., one end to the other end of rigid propulsion) is

(2)${\mathbf{T}}_{i+1}^{i^{\prime }}=\left[\begin{array}{cc}{0}_{3\times 3}& {\mathit{p}}_{i+1}^{i^{\prime }}\\ 0_{1\times 3}& 1\end{array}\right]$

A general kinematic model for the soft joint with continuum medium is the constant curvature approach [35] that defines the configuration states as the length L, bending angle θ, and the orientation φ of the deformation for each joint. Furthermore, the configuration state of soft joint i is determined by the three soft actuators’ length state L_i1~3 according to the following relationship (see Figure 4):

(3)$\begin{array}{l}{L}_i=\frac{L_{i1}+L_{i2}+L_{i3}}{3}\\ {\theta }_i=\frac{2\sqrt{L_{i1}^2+L_{i2}^2+L_{i3}^2-L_{i1}{L}_{i2}-L_{i1}{L}_{i3}-L_{i3}{L}_{i2}}}{3r_s}\\ {\phi }_i=\arctan 2(\sqrt{3}(L_{i3}-L_{i2}),L_{i2}+L_{i3}-2L_{i1})\end{array}$

where L_i1~3 = L₀ + q_i1~3, L₀ is the initial length,q_i= [q_i1, q_i2, q_i3] is the actuator space, and r_s is the distance from the actuator position to the center of the soft joint. The pose of soft joint body frame F _i with respect to rigid propulsion body frame F _i’ (i.e., from the bottom to the top of the joint) is given by the homogeneous transformation matrix H,

(4)${\mathbf{H}}_{i^{\prime }}^i(\mathit{q})=\{\begin{array}{c}\left[\begin{array}{cc}{\mathit{I}}_3& {\left[\begin{array}{ccc}{L}_i& 0& 0\end{array}\right]}^T\\ 0& 0\end{array}\right]\theta =0\\ \left[\begin{array}{cc}{\mathit{R}}_{i^{\prime }}^i(\mathit{q})& {\mathit{p}}_{i^{\prime }}^i(\mathit{q})\\ 0& 1\end{array}\right]\theta \ne 0\end{array}$

where I₃ is the third-order identity matrix. If all three actuators have the same length such that the soft joint is not bent but straight, the bending angle θ is zero. Only the lengths among them are different to give the soft joint bending and orientation angle. When θ ≠ 0, the specific expression is:

(5)${\mathbf{H}}_{i^{\prime }}^i(\mathit{q})=\left[\begin{array}{cccc}{c}_{{\phi }_i}^2\left(c_{{\theta }_i}-1\right)+1& c_{{\phi }_i}{s}_{\phi }{}_i\left(c_{{\theta }_i}-1\right)& c_{{\phi }_i}{s}_{{\theta }_i}{L}_i& \frac{c_{{\phi }_i}\left(1-c_{{\theta }_i}{L}_i\right)}{\theta }\\ s_{{\phi }_i}{c}_{{\phi }_i}\left(c_{{\theta }_i}{L}_i-1\right)& c_{{\phi }_i}^2\left(1-c_{{\theta }_i}l\right)+c_{{\theta }_i}{L}_i& s_{{\phi }_i}{s}_{{\theta }_i}{L}_i& \frac{s_{\phi }\left(1-c_{{\theta }_i}{L}_i\right)}{\theta }\\ -c_{{\phi }_i}{s}_{{\theta }_i}{L}_i& -s_{{\phi }_i}{s}_{{\theta }_i}{L}_i& c_{{\theta }_i}{L}_i& \frac{s_{{\theta }_i}{L}_i}{{\theta }_i}\\ 0& 0& 0& 1\end{array}\right]$

where $c_{{\phi }_i}=\cos {\phi }_i,s_{{\phi }_i}=\sin {\phi }_i,c_{{\theta }_i}=\cos {\theta }_i,s_{{\theta }_i}=\sin {\theta }_i$. Therefore, according to (2) and (4), it can be further deduced that the homogeneous transformation matrix from the F_I to F_e of the snake robot is

(6)${\mathbf{g}}_e^I={\mathbf{T}}_b^I{\mathbf{T}}_1^b\underset{{\mathbf{g}}_2^1}{\underbrace{{\mathbf{H}}_{1^{\prime }}^1{\mathbf{T}}_2^{1^{\prime }}}}\cdots \underset{{\mathbf{g}}_n^{n-1}}{\underbrace{{\mathbf{H}}_{(n-1{)}^{\prime }}^{(n-1)}{\mathbf{T}}_n^{(n-1{)}^{\prime }}}}\underset{{\mathbf{g}}_e^n}{\underbrace{{\mathbf{H}}_{n^{\prime }}^n{\mathbf{T}}_e^{n^{\prime }}}}$

where T and H represent the transformation matrix of rigid propulsion and soft joint, respectively, and ${\mathbf{g}}_{i+1}^i$ represents the transformation matrix of their combination.

(2) Differential Kinematics: The differential kinematics describes the velocity relationship between different parts of the snake robot. The motion of the robot can be determined by the speed vector ζ:

(7)$\mathsf{\zeta }=\left[\begin{array}{c}{\mathit{V}}_I^b\\ \dot{\mathit{q}}\end{array}\right]\in {\mathit{R}}^{6+3n},{\mathit{V}}_I^b=\left[\begin{array}{c}{\mathit{v}}_1\\ {\mathit{v}}_2\end{array}\right]$

where V ∈ R⁶ represents the body-fixed velocity vector of the robot base relative to the inertial frame F _I, including linear velocity vector v₁ ∈ R³ and angular velocity vector v₂∈R³, can be obtained by the velocity twist and the transformation matrix ${\mathbf{T}}_b^I$,

(8)$\begin{array}{l}{\dot{\mathbf{T}}}_b^I=({\mathit{V}}_I^b{)}^{\wedge }{\mathbf{T}}_b^I\Rightarrow \\ {({\mathit{V}}_I^b)}^{\wedge }={({\mathbf{T}}_b^I)}^{-1}{\dot{\mathbf{T}}}_b^I=\left[\begin{array}{cc}{\mathit{v}}_2^{\wedge }& {\mathit{v}}_1\\ 0& 0\end{array}\right]\in se(3)\end{array}$

where the operator ^ converts the angular velocity vector from R³ to skew-symmetric matrix R^3×3. The velocity derivation method of any part of the snake robot is similar, here taking the robot end as an example, the velocity twist of the end with respect to the inertial frame F _I can then be expressed in end frame F _e by

(9)$\begin{array}{l}{({\mathit{V}}_I^e)}^{\wedge }={({\mathbf{T}}_e^I)}^{-1}{\dot{\mathbf{T}}}_e^I\\ \hspace{1em}\hspace{1em} ={({\mathbf{T}}_b^I{\mathbf{T}}_e^b)}^{-1}({\dot{\mathbf{T}}}_b^I{\mathbf{T}}_e^b+{\mathbf{T}}_b^I{\dot{\mathbf{T}}}_e^b)\\ \hspace{1em}\hspace{1em} ={\mathbf{T}}_e^b{}^{-1}{({\mathit{V}}_I^b)}^{\wedge }{\mathbf{T}}_e^b+{({\mathit{V}}_b^e)}^{\wedge }\end{array}$

Next, we convert the twist expression in (9) into vector expression

(10)${\mathit{V}}_I^e=A{d}_{{\mathbf{g}}_e^b}^{-1}{\mathit{V}}_I^b+{\mathit{V}}_b^e$

where $A{d}_{{\mathbf{g}}_e^b}$ is the adjoint map for the transformation ${\mathbf{g}}_e^b$. The velocity variable ${\mathit{V}}_b^e$ represents the relative velocity between the robot base and end. When acting on a transformation g, the adjoint map and its inverse are given by

(11)$A{d}_{\mathbf{g}}=\left[\begin{array}{cc}\mathit{R}& {\mathit{p}}^{\wedge }\mathit{R}\\ 0& \mathit{R}\end{array}\right], A{d}_{\mathbf{g}}^{-1}=\left[\begin{array}{cc}{\mathit{R}}^T& -{\mathit{R}}^T{\mathit{p}}^{\wedge }\\ 0& {\mathit{R}}^T\end{array}\right]$

In (10), the inverse adjoint maps transform the body-fixed velocity vector ${\mathit{V}}_I^b$ so that it represents velocities in the end frame F _e. We write the relative end velocity as functions of the joint space velocities according to

(12)${\mathit{V}}_b^e=A{d}_{{\mathbf{g}}_e^b}^{-1}{\mathit{J}}_e\dot{\mathit{q}},{\mathit{J}}_e\triangleq \left[\begin{array}{ccc}A{d}_{{\mathbf{g}}_1^b}{\mathit{J}}_{s,1}& \cdots & A{d}_{{\mathbf{g}}_n^b}{\mathit{J}}_{s,n}\end{array}\right]$

where J_e is the Jacobian matrix of the robot end determined by J_s. J_s,i is Jacobian matrix of each soft joint, and the specific derivation refers to [37]. Finally, by inserting (12) into (10), we can obtain

(13)${\mathit{V}}_I^e={\mathit{J}}_{g,e}\mathsf{\zeta }, {\mathit{J}}_{g,e}\triangleq \left[\begin{array}{cc}A{d}_{{\mathbf{g}}_e^b}^{-1}& A{d}_{{\mathbf{g}}_e^b}^{-1}{\mathit{J}}_e\end{array}\right]$

where J_g,e is the geometric Jacobians, which map the base frame and the joint velocities to the linear and angular velocities of the robot end. See [34] for more theoretical details of this geometric method.

2.3. Dynamic Formulation

The snake robot is controlled by hybrid actuation of rigid propulsions and soft joints. An approach for three-dimensional dynamics modeling decouples motion equations of soft members from those of rigid members. Therefore, simpler sets of dynamics equations will be achieved so that obtained model can be used for control conveniently. In this section, this approach is divided into two sections. Firstly, we model the local dynamics of the soft joint, because the soft joint is different from the general rigid joint, but as an independent module component. Next, the local dynamics of soft members are extended as a nonholonomic constraint form for modeling the multibody rigid-soft snake robot system.

(1) Local Dynamic for Soft Joints: The local dynamic model of a soft joint is obtained by employing the Lagrangian analysis [38]. The equation of motion can be obtained by applying the Euler–Lagrange equation formulated in (14).

(14)$\mathit{F}=\left(\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \dot{\mathit{q}}}\right)-\frac{\partial \mathcal{L}}{\partial \mathit{q}}$

where F is the generalized force, including driving pressure force and load force, L is the Lagrangian defined as

(15)$\mathcal{L}=E_k-E_g-E_m$

where E_k is the kinetic energy, E_g is the gravitational potential energy, and E_m is the mechanical energy caused by the elastic and viscous stresses of soft materials.

Figure 4 depicts the configuration of a soft joint. Taking a similar approach as in literature [39], the soft module is assumed to be made up of infinitesimal slices of thickness Ldξ, where ξ ∈ [0,1] is the normalized position parameter along the axis length from the base to the top, where ξ = 0 is the base. The kinetic energy of the soft joint is

(16)$E_k(\mathit{q},\dot{\mathit{q}})=\frac{1}{2}{{\displaystyle {\int }_0^1({\mathit{V}}_{\xi }^b)}}^T\delta m{\mathit{V}}_{\xi }^{b}d\xi =\frac{1}{2}{\dot{\mathit{q}}}^T\mathit{m}(\mathit{q})\dot{\mathit{q}}$

where $\delta m=m_s\cdot diag\left\{\begin{array}{cccccc}1& 1& 1& \frac{1}{4}{R}_s^2& \frac{1}{4}{R}_s^2& \frac{1}{2}{R}_s^2\end{array}\right\}$ is the mass matrix of the slice, m_s and R_s are the mass and radius of the soft joint, respectively. m is the generalized mass matrix and ${\mathit{V}}_{\xi }^b\in R^6$ is the velocity vector, which are formulated as

(17)$\begin{array}{c}\mathit{m}(\mathit{q})={{\displaystyle {\int }_0^1({\mathit{J}}_{\xi }^b)}}^T\delta m{\mathit{J}}_{\xi }^{b}d\xi \\ {\mathit{V}}_{\xi }^b(\mathit{q})=\left[\begin{array}{c}{\mathit{R}}^T\mathit{p}\\ {({\mathit{R}}^T\mathit{R})}^{\vee }\end{array}\right]={\mathit{J}}_{\xi }^b(\mathit{q})\dot{\mathit{q}}\end{array}$

where ${\mathit{J}}_{\xi }^b\in R^{6\times 3}$ is the Jacobian matrix, according to constant curvature approach, relevant to L, θ, φ. If ξ = 1, R and p in (17) correspond to the rotation matrix and position vector in (4), and ${\mathit{J}}_{\xi }^b$ correspond to J_s in (4). The operator ^∨ is the inverse operation of operator ^, that is, converts the skew-symmetric matrix into a vector.

Next, the gravitational potential energy is

(18)$E_g(\mathit{q})=\mathit{m}(\mathit{q})\cdot {\displaystyle {\int }_0^1{\mathit{p}}^T}{\mathit{e}}_{1}d\xi$

where e₁ = [0.0, 9.8]^T is the gravity vector.

The mechanical stress-related energy takes into account the axial strain energy generated by the tension or compression of an actuator,

(19)$E_m=\frac{1}{2}{\displaystyle \sum _{i=1}^3{V}_{si}}{\sigma }_i{\epsilon }_i$

where i = 1, 2, 3 here represents the three actuators of one soft joint, σ_i is the tensile stress, ε_i is the axial strain and V_si is the chamber volume of an actuator.

Due to the special material properties of soft joints, the soft material models presented in this section are incorporated in the mechanical energy terms of (15) describing the module deformation energy and viscous losses. The sectional view of modeled soft joint in this study is presented in Figure 5, which is composed of multiple material parts. More specifically, each hydraulic soft actuator has a soft silicone chamber of internal radius r_a and silicone wall thickness h_c, surrounded by a fabric-reinforced layer of wall thickness h_s. So, a novel non-Newtonian visco-hyperelastic constitutive model is considered for the soft silicone. First, the circular fabric reinforcement prevents radial expansion in each actuator, and only allows it to elongate when actuated under water pressure. Considering axial tensile stress applied to each actuated chamber, the inverse rule of mixtures as in (20) is applied for determining the mechanical properties of the composite layer as the stress is transverse to the fabric [40].

(20)$\frac{1}{E_c}=\frac{V_f}{E_f}+\frac{V_m}{E_m}$

where E_c is the elastic modulus of the actuator compound body. E_f and E_m are the modulus of the fabric layer and silicone matrix, respectively, and V_f and V_m are the corresponding volumes. Furthermore, it reveals that upon water pressure actuation, in fact, the fabric and the silicone matrix volume fraction change. The formulae for the changing matrix volume fraction V_m and that of the silicone part associated with each of the three channels, V_si in (19), are as follows:

(21)$\begin{array}{l}{V}_{si}=\frac{\pi R_s^2-3\pi {(r_a+h_s+h_c)}^2}{3}+\frac{2\pi (r_a+h_c)h_s}{V_{mi}}+2\pi r_a{r}_s\\ V_{mi}=\frac{{\epsilon }_i+V_{m0}}{{\epsilon }_i+1}\end{array}$

where ${\epsilon }_i=\frac{q_i}{L_0}$ is the strain of the actuator i and V_m0 is the initial volume fraction of the silicone matrix.

The visco-hyperelastic constitutive model proposed for the soft actuator is described by a system of a non-Newtonian viscous damper and a hyperelastic spring in parallel, which is equivalent to the simplified Kelvin–Voigt model [41]. Therefore, σ in (19) includes hyperelastic stress σ_h and viscous stress σ_v.

The standard approach to describe the strain energy of the hyperelastic material is by applying a polynomial energy model [42]. For this study, the simplified polynomial Yeoh hyperelastic model as in (22) has been chosen to describe the energy potential.

(22)$U_h(\mathit{q})={\displaystyle \sum _{i=1}^3{C}_{i0}{({\overline{I}}_1(q_i)-3)}^i}$

where C_10~30 are the constant coefficient of soft material based on the Yeoh model, and its value can be measured by the uniaxial tensile test. ${\overline{I}}_1$ is the first partial strain invariant, depending on the type of load applied. For the case of axial strain under the constant volume hypothesis, the first deviatoric strain invariant for actuator i is described as follows:

(23)${\overline{I}}_1(q_i)={\left(1+\frac{q_i}{L_0}\right)}^2+2{\left(1+\frac{q_i}{L_0}\right)}^{-1}$

The viscous part of the visco-hyperelastic material model [43] is obtained from the following equation:

(24)${\sigma }_v=\eta \frac{d\epsilon (t)}{dt}$

where η is the power law viscosity defined as

(25)$\eta (\dot{\epsilon })={\eta }_0{\left|\frac{d\epsilon (t)}{dt}\right|}^{k-1}$

where η₀ is the consistency of the material and k is the power law index for a non-Newtonian fluid. If k < 1, the power law predicts that the effective viscosity would decrease with increasing deformation rate, which is usually the case for polymeric materials. Thus, the viscous energy potential of the visco-hyperelastic constitutive model is described.

(26)$U_v(\mathit{q},\dot{\mathit{q}})=\frac{1}{2}{\eta }_0\epsilon (t)\frac{d\epsilon (t)}{dt}{\left|\frac{d\epsilon (t)}{dt}\right|}^{k-1}$

The stress σ_h and σ_v can be calculated, respectively according to (22) and (26), then the axial strain energy can be calculated by substituting (21) for (19). Finally, substituting the energy Equations (16), (18), and (19) for (15) and (14) to further deduce the local dynamics of the soft joint. Thus, the state variables q of the soft joints can be obtained.

(2) System Dynamic for Snake Robot: The proposed dynamic approach separates dynamics modeling of the soft elements from the unified overall system. Subsequently, the coupling between the soft members and overall system corresponding sets of equations is realized by the nonholonomic kinematic constraints. According to the local dynamics of the soft joint obtained above, on the basis of the floating frame hypothesis [44], we uniformly use the Lagrange method to calculate the dynamic model of the snake robot. Attaching a floating frame to each rigid member, the motion of the latter rigid propulsion module is specified by a set of configurations of the former soft joint. Finally, the total three-dimensional dynamics model of this rigid–soft multibody system is presented.

The snake robot is considered a free-floating serial chain manipulator [24]. The dynamic formulations are therefore structurally similar to those of underwater vehicle manipulator systems presented in [45]. In this paper, the velocity vector ζ, including the base and all joint components, has been defined in (7). The equations of motion can also be derived from the Lagrangian method [45],

(27)$\mathit{M}(\mathit{q})\dot{\mathsf{\zeta }}+\mathit{D}(\mathit{q},\mathsf{\zeta })\mathsf{\zeta }+\mathit{C}(\mathit{q},\mathsf{\zeta })\mathsf{\zeta }+\mathit{G}(\mathit{q})=\mathsf{\tau }(\mathit{q})$

where M(q) is the inertia matrix, D(q, ζ) is the hydrodynamical drag matrix, C(q, ζ) is the Coriolis and centrifugal force matrix, G(q) is hydrostatic force contained weight and buoyancy, and τ(q) is the driving control force including soft joints and rigid propulsions.

First, M(q) is an inertia matrix containing rigid body and hydrodynamical inertia terms, and D(q, ζ) is a hydrodynamical drag-wrench-related item. For the purpose of simulation, it is extremely challenging to determine the required M(q) and D(q, ζ) to achieve (27). In practice, most are forced to adopt additional simplified assumptions at the expense of model accuracy, which makes this task more manageable. Therefore, it can be assumed that the modules are hydrodynamically decoupled [46]. Under this assumption, there will be no hydrodynamic interaction between modules, and the static surface will confine the fluid. Consequently, the matrices M(q) and D(q, ζ) become block-diagonal and configuration-independent [47].

The inertia matrix M(q) is made up of rigid body and hydrodynamic contributions also called added mass. The rigid body is simply the mass of the body itself, while the added mass forces are hydrodynamic pressure-induced forces developing when it accelerates. As the snake robot consists of many modules in a chain, the inertia matrix contains the sum for each module transformed to the base frame through the Jacobian found in (13),

(28)$\begin{array}{l}\mathit{M}(\mathit{q})={\displaystyle \sum _{i=1}^n{\mathit{J}}_{g,i}^T(\mathit{q})}{\mathit{M}}_i(\mathit{q}){\mathit{J}}_{g,i}(\mathit{q})\\ {\mathit{M}}_i(\mathit{q})={\mathit{M}}_{R,i}(\mathit{q})+{\mathit{M}}_{A,i}(\mathit{q})\end{array}$

where M_R,i = ${\mathit{M}}_{R,i}^T$ is the rigid-body inertia matrix and M_A,i = ${\mathit{M}}_{A,i}^T$ the added inertia matrix. A general expression for M_R,i can be given as,

(29)${\mathit{M}}_{R,i}=\left[\begin{array}{cc}{m}_{R,i}{\mathit{I}}_3& m_{R,i}{({\mathsf{\gamma }}_{g,i}^{\wedge })}^T\\ m_{R,i}{\mathsf{\gamma }}_{g,i}^{\wedge }& {\mathit{I}}_{R,i}\end{array}\right], {\mathit{I}}_{R,i}=diag(I_{xx,i},I_{yy,i},I_{zz,i})$

where m_R,i denotes the mass of rigid propulsion module i, γ_g,ithe position of the center of mass of module i expressed in F _i’, and I_R,i the rigid body rotational inertia matrix, where $I_{yy,i}=I_{zz,i}=\raisebox{1ex}{$\pi R_r^2$}\!\left/ \!\raisebox{-1ex}{$4$}\right.,I_{xx,i}=I_{yy,i}+I_{zz,i}$, R_r is the radius of rigid propulsion consistent with the soft joint.

The added inertia matrix depends on the shape of the module. It must be emphasized that M_A,i appears under the assumption of hydrodynamically decoupled links. A more practical approach in the context of commonly encountered cylindrical geometries is to take the hydrodynamical inertia matrix as follows,

(30)${\mathit{M}}_{A,i}=\rho \pi R_r^2{L}_{R,i}{C}_a\left[\begin{array}{cccccc}{\alpha }_i& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& \frac{1}{2}{L}_{R,i}\\ 0& 0& 1& 0& -\frac{1}{2}{L}_{R,i}& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& \frac{1}{2}{L}_{R,i}& 0& -\frac{1}{3}{L}_{R,i}^2& 0\\ 0& \frac{1}{2}{L}_{R,i}& 0& 0& 0& \frac{1}{3}{L}_{R,i}^2\end{array}\right]$

where ρ is the density of water, L_R,i the length of rigid propulsion module i, C_a is the added mass coefficient, and α_i is a parameter that allows added mass in the surge to be added of each module. These parameters must be determined by experiment or high-quality calculation method.

The drag wrenches are selected as a combination of linear and nonlinear components of the form:

(31)$\begin{array}{l}\mathit{D}(\mathit{q},\mathsf{\zeta })={\displaystyle \sum _{i=1}^n{\mathit{J}}_{g,i}^T(\mathit{q})}{\mathit{D}}_i(\mathit{q},\mathsf{\zeta }){\mathit{J}}_{g,i}(\mathit{q})\\ {\mathit{D}}_i(\mathit{q},\mathsf{\zeta })={\mathit{D}}_{NL,i}(\mathit{q},\mathsf{\zeta })+{\mathit{D}}_{L,i}\end{array}$

where D_NL,i(q, ζ) is the nonlinear drag force contribution and D_L,i is the linear. The nonlinear resistance is calculated by strip theory [48], in which the difficult surface integral is replaced by a line integral over the length of the cylindrical object divided into disk units. The contribution of linear drag influence is given by the following equation,

(32)${\mathit{D}}_{L,i}=\rho \pi L_{R,i}{C}_d{\mathit{V}}_I^i\left[\begin{array}{cccccc}{\beta }_i& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& \frac{1}{2}{L}_{R,i}\\ 0& 0& 1& 0& -\frac{1}{2}{L}_{R,i}& 0\\ 0& 0& 0& {\mathsf{\Omega }}_i{R}_S^2& 0& 0\\ 0& 0& -\frac{1}{2}{L}_{R,i}& 0& \frac{1}{3}{L}_{R,i}^2& 0\\ 0& \frac{1}{2}{L}_{R,i}& 0& 0& 0& \frac{1}{3}{L}_{R,i}^2\end{array}\right]$

where ${\mathit{V}}_I^i$ and C_d are the reference velocity and drag coefficient, respectively. β_i and Ω_i are parameters which allow linear drag effects to be added in surge and roll. It should be noted that it is also a challenging task to determine these parameters, because they depend on other parameters that might not be constant during the motion.

The centripetal and Coriolis forces are described by

(33)$\begin{array}{l}\mathit{C}(\mathit{q},\mathsf{\zeta })={\displaystyle \sum _{i=1}^n\left[{\mathit{J}}_{g,i}^T(\mathit{q}){\mathit{M}}_i\frac{d}{dt}({\mathit{J}}_{g,i}^T(\mathit{q}))-{\mathit{J}}_{g,i}^T(\mathit{q})W_i(\mathit{q},\mathsf{\zeta }){\mathit{J}}_{g,i}(\mathit{q})\right]}\\ W_i(\mathit{q},\mathsf{\zeta })=\left[\begin{array}{cc}{0}_{3\times 3}& {\left\{{\mathit{M}}_i{\mathit{V}}_I^i\right\}}_{v_1}^{\wedge }\\ {\left\{{\mathit{M}}_i{\mathit{V}}_I^i\right\}}_{v_1}^{\wedge }& {\left\{{\mathit{M}}_i{\mathit{V}}_I^i\right\}}_{v_2}^{\wedge }\end{array}\right]\end{array}$

where W_i(q, ζ) is an intermediate variable symbol. ${\left\{{\mathit{M}}_i{\mathit{V}}_I^i\right\}}_{v_1}\in R^3$ and ${\left\{{\mathit{M}}_i{\mathit{V}}_I^i\right\}}_{v_2}\in R^3$ are the first three entries and the last three entries of ${\mathit{M}}_i{\mathit{V}}_I^i$, respectively. Although this term is traditionally referred to as centripetal and Coriolis force, in fact, it takes into account forces that are more complex than centripetal force and Coriolis force alone [49].

When the robot is in the water, hydrostatic force is caused by hydrostatic pressure on the robot. Because the horizontal hydrostatic forces acting on the modules cancel each other out, only the vertical hydrostatic forces contribute to this term. Additionally, this term is given as a sum of the hydrostatic force contribution to each module transformed to the base frame,

(34)$\begin{array}{l}\mathit{G}(\mathit{q})={\displaystyle \sum _{i=1}^n{\mathit{J}}_{g,i}^T(\mathit{q})}{\mathit{G}}_i(\mathit{q}){\mathit{J}}_{g,i}(\mathit{q})\\ {\mathit{G}}_i(\mathit{q})=\left[\begin{array}{c}(\rho \pi R_r^2{L}_{R,i}-m_{R,i}){\mathit{I}}_3\\ \rho \pi R_r^2{L}_{R,i}{\mathsf{\gamma }}_{{}_{b,i}}^{\wedge }-m_{R,i}{\mathsf{\gamma }}_{{}_{g,i}}^{\wedge }\end{array}\right]{({\mathit{R}}_i^I)}^T{\mathit{e}}_1\end{array}$

where γ_b,i and γ_g,i are the location of the center of buoyancy and gravity of rigid propulsion i expressed in F_{i ‘}. Furthermore, the influence of steady (time constant) and uniform (space constant) water currents can also be included without much difficulty.

The driving forces and torques τ that are applied to the snake robot are provided by joint actuators and propellers, and τ can be divided into joint torques τ_J and propulsions torques τ_T,

(35)$\mathsf{\tau }(\mathit{q})=\left[\begin{array}{c}{\mathsf{\tau }}_T(\mathit{q})\\ {\mathsf{\tau }}_J(\mathit{q})\end{array}\right]=\left[\begin{array}{cc}\mathit{B}& 0_{(6+3n)\times 3n}\\ 0_{3n\times m}& \mathit{K}\end{array}\right]\left[\begin{array}{c}{\mathit{u}}_T\\ {\mathit{u}}_J\end{array}\right]$

where B ∈ R^{(6 + 3n) × m} is the propellers’ configuration matrix, u_T ∈ R^m and u_J ∈ R³ⁿ are the control inputs of propellers and joints, respectively. According to our previous research on three-chamber actuated soft manipulators [50], u_J and τ_J is the linear relation and K is the linear coefficient matrix. Contrary to the traditional ROVs or AUVs, the position of the propellers relative to the robot base depends on the joints’ configuration. The propellers also affect the position and orientation of the base as well as the joint configurations. The control input u_T is applied to the propellers to achieve the desired force on the base, while the joints’ configuration is controlled by the joint actuators. The propulsion forces τ_T corresponding to a control input u_T depend on the structural variables of the propellers as well as the fluid density. Therefore, the relationship between the propulsion’s forces τ_T and the control input u_T is highly nonlinear [47]. The relationship can however be simplified to a linear relationship. Each propulsion module has a matrix B_i∈R^6×mi associated with it to describe how the propellers exert forces and torques on that module, where mi is the number of propellers of module i. The matrix B_i is found from the thrust direction ${\mathsf{\beta }}_{t,i,j}$ and the point of attack ${\mathsf{\chi }}_{t,i,j}$ on the propeller j of module i expressed in F_i’. B_i can then be written as:

(36)${\mathit{B}}_i=\left[\begin{array}{ccc}{\mathsf{\beta }}_{t,i,j}& \cdots & {\mathsf{\beta }}_{t,i,m}\\ {\mathsf{\chi }}_{t,i,j}\times {\mathsf{\beta }}_{t,i,j}& \cdots & {\mathsf{\chi }}_{t,i,j}\times {\mathsf{\beta }}_{t,i,m}\end{array}\right]$

The forces and torques on the base and joints from all the propellers can therefore be expressed through the Jacobian matrix as:

(37)${\mathsf{\tau }}_T={\displaystyle \sum _{i=1}^n{\mathit{J}}_{g,i}^T{\mathsf{\tau }}_i}=\left[\begin{array}{cccc}{\mathit{J}}_{g,1}^T{\mathit{B}}_1& {\mathit{J}}_{g,2}^T{\mathit{B}}_2& \cdots & {\mathit{J}}_{g,n}^T{\mathit{B}}_n\end{array}\right]{\mathit{u}}_t$

A flow chart of the snake robot system dynamics of (27), including the local dynamics of soft joints, and joint actuators and propellers’ inputs, can be seen in Figure 6. Finally, we can solve kinematics and dynamics parameters such as robot position, pose, force, linear velocity, and angular velocity by numerical Runge–Kutta method in a MATLAB environment. Therefore, multiple spatial motions of the snake robot will be widely generated in simulations, and then estimated through the following experiments.

3. Results

In this section, we first describe the design prototype of the proposed snake robot and then show four simulation cases to reflect the multi-modal motion of the robot. Finally, the practicability of the model is verified through experiments.

3.1. Robot Prototype Overview

The proposed rigid-soft snake robot is a mechanically robust and easily reconfigurable experimental object, which is developed to support the research on snake robot behavior. The main characteristics of the robot are as follows:

• Modular: The modular design of the robot is adopted, and the rigid propulsion and soft joint modules have universal mechanical and electrical interfaces.

• Flexible: Compared to rigid joint, The soft one can be bent in all directions, which greatly improves flexibility, especially in a few modules.

• Waterproof: All modules of the robot are waterproof at water depths down to at least 5 m.

The types of modules developed for the snake robot include head modules, control modules, soft joint modules, rigid propulsion modules, and tail modules. The snake robot verified in this study is shown in Figure 7a, including 3 rigid modules and 2 soft joints. These rigid modules, whose casings have been 3D printed in plastic using selective laser sintering. The propellers are installed externally and internally to fabricate a rigid propulsion module. Since our soft actuator adopts a fabric-reinforced structure, possessing high strength, the output force and stiffness are also considerable, which can ensure the controllability of the posture during the motion of the snake robot. It must be pointed out that, in the manufacturing process of soft actuators, as multiple steps will inevitably produce errors, we have carried out the static pressure test for three actuators of one joint, and the static bending and force pressure test for one joint, as shown in Figure 8. The error below 0.5 MPa is lower, so the subsequent test is within this interval. Furthermore, in the design and manufacturing process, the mass distribution of the overall robot has been carefully considered. Firstly, the centroid of the robot is calculated by setting the density/mass of each module in Solidworks. Secondly, assuming that the robot has a uniform mass distribution, the buoyancy center of the robot can be obtained by Solidworks. Finally, the density design of the robot is close to water. The center of mass of the robot is slightly lower than the center of buoyancy, but on the same vertical axis as the center of buoyancy. These features contribute to the planar and spatial motion of the snake robot. We have verified the usability of this designed structure in the pool. Taking the horizontal plane turning as an example, the movement process of the robot is shown in Figure 9.

The main controller of the robot adopts an STM32 used for locomotion control, attitude calculation, and data acquisition. Moreover, a credit-card-sized Raspberry Pi acting as an auxiliary processor is used for image processing, localization, and motion planning. Figure 7b shows the implementation diagram of the robot control system. Moreover, the robot integrates several types of sensors onboard. An inertial measurement unit (IMU) is fixed parallel to the body’s principal axes to monitor the motion information of the robot. The motion controller, including 6 propellers and 6 soft actuators, serves as the actuation of the robot system to generate multiple motions such as forward and backward motion, turning, ascending, descending, and spiraling motion. Note that the sonar and depth sensors are not used in this study. The snake robot moves underwater, and the joints do not need to be sealed, so the water medium is a natural driving source. The external connection method was adopted to provide water pressure from the land. The robot has two soft joints, that is, six soft actuators. The driving equipment used is six water pressure cylinders, so the robot is dragging water pipes behind it. However, we will study how to remove the water pipe and integrate the micro water pressure pump into the robot body. Preliminary plan has been made, as shown in Appendix A.

3.2. Numerical Simulation

First of all, for the values of the design parameters described previously, see Table 1. Additionally, the snake robot works in water, characterized by environment parameters summarized in Table 2.

In order to verify the effectiveness and generality of the dynamic model studied in this paper, the multi-modal swimming behavior of the robot has been widely simulated, including forward and backward rectilinear motion, yaw turning, pitch motion, and spiral rising motion.

CASE 1: Rectilinear motion

Rectilinear motion is one of the most basic motion modes of snake robot, which can move quickly in spacious waters. The base frame of the robot in the initial state coincides with the inertial frame, and the body is in the positive direction of the X-axis. Corresponding to the propeller and actuator index in Figure 7a, the control input corresponds to (36). We first set the propeller control parameters u_T(1) = u_T(3) = u_T(4) = u_T(6) = 1, u_T(2) = u_T(5) = 0. It can be seen that the forward motion trajectory of the robot is as shown in Figure 9. Next, set u_T(1) = u_T(3) = u_T(4) = u_T(6) = −1, u_T(2) = u_T(5) = 0 to make the robot move in reverse. The above control input is named condition 1. In addition, in linear motion, one or more groups of propellers can be controlled to change the motion speed, such as u_T(1) = u_T(3) =1, u_T(4) = u_T(6) = u_T(2) = u_T(5) = 0. Such control input is named condition 2. Figure 10 shows the position and velocity of the forward and backward motion under two conditions. In this case, only the propeller is manipulated to perform the rectilinear motion of the robot. It is also possible to only manipulate the joints to perform a winding forward gait, but Pettersen et al. has confirmed that the efficiency of only the propeller being manipulated is higher than that of joints [51].

CASE 2: Yaw motion

Yaw motion is represented by turning motion in a large space. For the snake robot, the turning motion can be realized in different ways. For example, the soft joint input is zero to keep it straight, and then input to the propeller on one side without input or negative signal on the other side, which can make the robot do the overall turning motion like an AUV. Or the input signal of the actuator makes the soft joint bend and then turns on the propeller to make the robot move in circles. We did a simulation comparison test, and the results can be seen in Appendix B. In fact, the motion of implementation achieved by thrusters and joints is the unique advantage of the snake robot, especially in the underwater confined space, the conventional AUVs are difficult or unable to pass, but the snake robot has a flexible body and can be bent locally, which will be easy to pass through. In the simulation, the propeller control parameters are set to u_T(1) = u_T(3) = u_T(4) = u_T(6) = 1, u_T(2) = u_T(5) = 0, and the soft joint control parameters are set to u_J(1) = u_J(4) = 0.3, u_J(2) = u_J(3) = u_J(5) = u_J(6) = 0. The above control input is named condition 1, which can make the robot do yaw and turn in the XOY plane. In order to verify the generality and rationality of the model, the control parameters of the soft joint are reset to u_J(1) = u_J(4) = 0.1, and other parameters remain unchanged. Such a control input is named condition 2. The trajectories of the motion under conditions 1 and 2 are shown in Figure 11, the configuration of soft joints, and the spatial position, linear velocity, orientation, and angular velocity are shown in Figure 11. It is clear that the range of motion increases as the bending angle of the joint decreases, i.e., the radius of gyration increases. Interestingly, the forward speed controlled by the propeller does not change and only changes the bending angle of the joint, which makes the positions of the rotation centers of the two conditions do not coincide and not on the same line. This relationship enables additional considerations for motion controller design and path planning for robotic snakes in practical applications requiring various steering conditions.

CASE 3: Pitch motion

Pitch motion can also be realized in different control modes. For example, the longitudinal propeller can be manipulated to make the robot float or sink in a fixed attitude. However, due to the large side area, the water resistance will be large. The advantage of the snake robot is that it can coordinate the manipulation of the propeller and joints and can quickly reach the predetermined altitude. The simulation shows the pitching motion of the robot on the XOZ plane, in which the first sets of parameters are set as u_T(1) = u_T(3) = u_T(4) = u_T(6) = 1, u_T(2) = u_T(5) = 0, u_J(1) = 0.2, u_J(2) = 0.3, u_J(3) = u_J(4) = u_J(5) = u_J(6) = 0. The second sets of different control parameters are set as the soft joint parameters are changed to u_J(1) = 0.2, u_J(3) = 0.1, u_J(2) = 0, and the others remain unchanged. The pitch and heave motion trajectories of the two conditions are shown in Figure 12. The configuration of soft joints is shown in Figure 12, it should be pointed out that the pitch motion is in the XOZ plane, so the orientation angle of the soft joint is near 90°. Additionally, the spatial position, linear velocity, orientation, and angular velocity are shown in Figure 12. It is necessary to explain the reason for the oscillation at the beginning of the motion because the position of the center of gravity and the center of buoyancy changes in the Z-axis direction, which leads to the change of the propeller configuration matrix, and then becomes smooth as the motion becomes stable.

CASE 4: Spiral motion

Spiral motion is one of the three-dimensional motion modes of snake robot. For spiral motion, the soft joint configuration needs to be adjusted so that the bending angle exists, and the orientation angle is between −π/2~π/2 rad. Therefore, the control parameters are set as follows: the parameters of propellers are set as u_T(1) = u_T(3) = u_T(4) = u_T(6) = 1, u_T(2) = u_T(5) = 0, and the control parameters of soft joints are set as u_J(1) = u_J(2) = 0.2, u_J(3) = u_J(4) = 0.3, u_J(5) = u_J(6) = 0. The trajectory of the motion is shown in Figure 13, and the configuration of soft joints is shown in Figure 13. In this case, we only do one set of simulations, but show the state of the base and end of the robot, and the three-dimensional position, orientation, linear and angular velocity are shown in Figure 13. In the simulation results, we found that the states of the base and the end are inconsistent. The main reason is the existence of joints, especially the soft joints that have less damping and are easier to deform, unlike the integral structure of AUV. Although the combination of soft joints can improve the performance of the snake robot, it cannot be ignored that it will pose a difficult problem for control. In fact, it reveals some general challenges that soft robot technology is faced with, that is, high nonlinearity, serious hysteresis, and poor repeatability, which mainly influence the accuracy.

4. Conclusions

In this study, we have developed a complete three-dimensional dynamic model for an underwater robotic snake actuated by rigid propulsions and soft joints. The modeling method of the snake robot is extended by deriving expressions for the geometric Jacobians. This modeling approach is well suited since it provides compact matrix expressions, easy implementation, and flexibility that accommodates an arbitrary number of modules, with variable module length. The geometric Jacobians are then utilized to derive the state-dependent rigid propulsion inertia matrix and the state-dependent propellers configuration matrix. For the soft joint part, we use the constant curvature method to describe its configuration, use the Lagrangian method to obtain its dynamic characteristics, and focus on deriving the visco-hyperelastic mechanical energy of the soft material. Finally, the local dynamics of soft members are extended as a nonholonomic constraint form for modeling the multibody rigid-soft snake robot. Simulations show that the model can predict multimodal swimming behaviors of the robot, such as forward and backward motion, yaw turning, pitch maneuvers, and even complicated spiraling motions. The results of the robot have demonstrated the effectiveness of the three-dimensional dynamic model in describing multiple swimming modes. In summary, the work presented in this paper provides a versatile dynamic model capable of predicting multiple three-dimensional motions for the hybrid rigid-soft actuated robotic snake.

Author Contributions

Conceptualization, J.Z. and Y.C.; methodology, J.Z.; software, J.Z.; validation, J.Z.; formal analysis, J.Z.; investigation, J.Z.; data curation, J.Z.; writing—original draft preparation, J.Z.; writing—review and editing, Y.C.; visualization, J.Z. and Y.L.; supervision, Y.G.; project administration, Y.G.; funding acquisition, Y.C. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Comparison of various characteristics of common underwater vehicles.

Figure 2. Schematic diagram of multi-modal motion of snake robot: (a) linear motion, (b) turning motion, (c) pitching motion, (d) helical motion.

Figure 3. Schematic diagram of underwater rigid-soft snake robot model, where white parts are rigid propulsions and blue parts are soft joints.

Figure 4. Configuration diagram of a soft joint.

Figure 5. Sectional view and geometric parameters of the soft joint.

Figure 6. Flowchart of the solution for the dynamic equations.

Figure 7. Mechanical and hardware configurations of the snake robot: (a) robot prototype; (b) control system.

Figure 8. Static pressure test results for one soft joint.

Figure 9. Test results of horizontal turning motion of the robot.

Figure 10. Simulation tests for forward and backward rectilinear motion, which contain the trajectories, position and velocity under condition 1, and position and velocity under condition 2. Note that position and velocity are in the inertial reference frame.

Figure 11. Simulation tests for yaw turning motion under different conditions, which contain the trajectories under conditions 1 and 2, configurations of soft joints, position and orientation under condition 1, linear and angular velocity under condition 1, position and orientation under condition 2, and linear and angular velocity under condition 2.

Figure 12. Simulation tests for pitch motion under different conditions, which contain the trajectories under conditions 1 and 2, configurations of soft joints, position and orientation under condition 1, linear and angular velocity under condition 1, position and orientation under condition 2, and linear and angular velocity under condition 2.

Figure 13. Simulation tests for spiral ascending motion under different conditions, which contain the trajectories, configurations of soft joints, position and orientation of robot base, linear and angular velocity of robot base, position and orientation of robot end, and linear and angular velocity of robot end.

Table 1

Design parameters of the prototype snake robot.

Symbol	Meaning	Value
L 0	soft joint length	200 mm
Rs	soft joint radius	90 mm
ra	actuator radius	8 mm
rs	actuator to soft joint center	24 mm
hc	silicone wall thickness	6 mm
hs	fabric layer thickness	2 mm
C 10~30	silicone material parameters	12,563 J/m3−67.784 J/m32.7385 J/m3
η 0	consistency of the material	22,263.8 Pa s0.3
k	power law index	0.3
Em	elastic modulus	2.048 × 105 Pa
ms	soft joint mass	158 g
mR	rigid propulsion mass	563 g
χ	point of attack	76 mm
LR	rigid propulsion length	200 mm
Rr	rigid propulsion radius	90 mm

Table 2

Environment parameters.

Symbol	Meaning	Value
Ca	added mass coefficient	1
αi	added mass in the surge	0.2
CD	drag coefficient	1
βi	linear drag effects in surge	0.2
Ωi	linear drag effects in roll	0.5
ρ	Water density	1000 kg/m3
e 1	Gravitational acceleration	[0.0–9.81]

Appendix A

The appendix Figure A1 shows the preliminary plan of how to remove the water pipe and integrate the micro water pressure pump into the robot body.

Figure A1. Physical diagram of micro water pump integrated system for underwater snake robot.

Appendix B

The appendix Figure A2 shows the comparison of the snake robot turning in straight line state and bending state. In fact, the executive motion realized by thrusters and joints is the unique advantage of the snake-like robot. Especially in the underwater confined space, it is difficult or impossible for the traditional AUV to pass through, but the snake-like robot has a flexible body, which can bend locally and pass easily.

Figure A2. Simulation tests of comparison of turning performance, (a) in a straight state and (b) bending state.

References

1. Zhou, Z.; Liu, J.; Yu, J. A Survey of Underwater Multi-Robot Systems. IEEE/CAA J. Autom. Sin.; 2021; 9, pp. 1-18. [DOI: https://dx.doi.org/10.1109/JAS.2021.1004269]

2. Yoerger, D.R.; Govindarajan, A.; Howland, J.C.; Llopiz, J.K.; Wiebe, P.H.; Curran, M.; Fujii, J.; Gomez-Ibanez, D.; Katija, K.; Robison, B.H. et al. A hybrid underwater robot for multidisciplinary investigation of the ocean twilight zone. Sci. Robot.; 2021; 6, eabe1901. [DOI: https://dx.doi.org/10.1126/scirobotics.abe1901] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34135116]

3. Huvenne, V.A.; Robert, K.; Marsh, L.; Iacono, C.L.; le Bas, T.; Wynn, R.B. ROVs and AUVs. Submarine Geomorphology; Springer: Cham, Switzerland, 2018; pp. 93-108.

4. Teague, J.; Allen, M.J.; Scott, T.B. The potential of low-cost ROV for use in deep-sea mineral, ore prospecting and monitoring. Ocean. Eng.; 2018; 147, pp. 333-339. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2017.10.046]

5. Carreras, M.; Hernández, J.D.; Vidal, E.; Palomeras, N.; Ribas, D.; Ridao, P. Sparus II AUV—A hovering vehicle for seabed inspection. IEEE J. Ocean. Eng.; 2018; 43, pp. 344-355. [DOI: https://dx.doi.org/10.1109/JOE.2018.2792278]

6. Katzschmann, R.K.; DelPreto, J.; MacCurdy, R.; Rus, D. Exploration of underwater life with an acoustically controlled soft robotic fish. Sci. Robot.; 2018; 3, eaar3449. [DOI: https://dx.doi.org/10.1126/scirobotics.aar3449] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33141748]

7. Wang, W.; Dai, X.; Li, L.; Gheneti, B.H.; Ding, Y.; Yu, J.; Xie, G. Three-dimensional modeling of a fin-actuated robotic fish with multimodal swimming. IEEE/ASME Trans. Mechatron.; 2018; 23, pp. 1641-1652. [DOI: https://dx.doi.org/10.1109/TMECH.2018.2848220]

8. Kelasidi, E.; Liljeback, P.; Pettersen, K.Y.; Gravdahl, J.T. Innovation in underwater robots: Biologically inspired swimming snake robots. IEEE Robot. Autom. Mag.; 2016; 23, pp. 44-62. [DOI: https://dx.doi.org/10.1109/MRA.2015.2506121]

9. Huo, X.; Wang, X.; Ge, T. Impulse control method for hydraulic propulsion system used in 3500 m work-class ROV. Appl. Ocean. Res.; 2016; 60, pp. 75-83. [DOI: https://dx.doi.org/10.1016/j.apor.2016.08.008]

10. Capocci, R.; Dooly, G.; Omerdić, E.; Coleman, J.; Newe, T.; Toal, D. Inspection-class remotely operated vehicles—A review. J. Mar. Sci. Eng.; 2017; 5, 13. [DOI: https://dx.doi.org/10.3390/jmse5010013]

11. Ridao, P.; Carreras, M.; Ribas, D.; Sanz, P.J.; Oliver, G. Intervention AUVs: The next challenge. Annu. Rev. Control; 2015; 40, pp. 227-241. [DOI: https://dx.doi.org/10.1016/j.arcontrol.2015.09.015]

12. Yamada, H. Development of Amphibious Snake-like Robot ACM-R5. Proceedings of the the 36th International Symposium on Robotics (ISR 2005); Tokyo, Japan, 29 November–1 December 2005.

13. Porez, M.; Boyer, F.; Ijspeert, A.J. Improved Lighthill fish swimming model for bio-inspired robots: Modeling, computational aspects and experimental comparisons. Int. J. Robot. Res.; 2014; 33, pp. 1322-1341. [DOI: https://dx.doi.org/10.1177/0278364914525811]

14. Zhao, X.; Fei, F. The Development Status and Trend of Commercial Bionic Robots. Proceedings of the 2021 5th International Conference on Electronic Information Technology and Computer Engineering; Xiamen, China, 22–24 October 2021; pp. 1377-1381.

15. Sverdrup-Thygeson, J.; Kelasidi, E.; Pettersen, K.Y.; Gravdahl, J.T. The underwater swimming manipulator—A bioinspired solution for subsea operations. IEEE J. Ocean. Eng.; 2017; 43, pp. 402-417. [DOI: https://dx.doi.org/10.1109/JOE.2017.2768108]

16. Liljebäck, P.; Stavdahl, Ø.; Pettersen, K.Y.; Gravdahl, J.T. Mamba-A waterproof snake robot with tactile sensing. Proceedings of the 2014 IEEE/RSJ International Conference on Intelligent Robots and Systems; Chicago, IL, USA, 14–18 September 2014; pp. 294-301.

17. Liljebäck, P.; Mills, R. Eelume: A flexible and subsea resident IMR vehicle. Proceedings of the OCEANS 2017—Aberdeen Conference; Aberdeen, UK, 19–22 June 2017; pp. 1-4.

18. Zwilgmeyer, G.O. Creating a Synthetic Underwater Dataset for Egomotion Estimation and 3D Reconstruction. Master’s Thesis; Norwegian University of Science and Technology: Trondheim, Norway, 2021.

19. Xu, F.; Wang, H.; Au, K.W.S.; Chen, W.; Miao, Y. Underwater dynamic modeling for a cable-driven soft robot arm. IEEE/ASME Trans. Mechatron.; 2018; 23, pp. 2726-2738. [DOI: https://dx.doi.org/10.1109/TMECH.2018.2872972]

20. Xu, F.; Wang, H.; Wang, J.; Au, K.W.S.; Chen, W. Underwater dynamic visual servoing for a soft robot arm with online distortion correction. IEEE/ASME Trans. Mechatron.; 2019; 24, pp. 979-989. [DOI: https://dx.doi.org/10.1109/TMECH.2019.2908242]

21. Gong, Z.; Fang, X.; Chen, X.; Cheng, J.; Xie, Z.; Liu, J.; Chen, B.; Yang, H.; Kong, S.; Hao, Y. et al. A soft manipulator for efficient delicate grasping in shallow water: Modeling, control, and real-world experiments. Int. J. Robot. Res.; 2021; 40, pp. 449-469. [DOI: https://dx.doi.org/10.1177/0278364920917203]

22. Wang, Z.; Cui, W. For safe and compliant interaction: An outlook of soft underwater manipulators. Proc. Inst. Mech. Eng. Part M J. Eng. Marit. Environ.; 2021; 235, pp. 3-14. [DOI: https://dx.doi.org/10.1177/1475090220950911]

23. Li, G.; Chen, X.; Zhou, F.; Liang, Y.; Xiao, Y.; Cao, X.; Zhang, Z.; Zhang, M.; Wu, B.; Yin, S. et al. Self-powered soft robot in the Mariana Trench. Nature; 2021; 591, pp. 66-71. [DOI: https://dx.doi.org/10.1038/s41586-020-03153-z] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33658693]

24. Kelasidi, E.; Pettersen, K.Y.; Gravdahl, J.T.; Liljebäck, P. Modeling of underwater snake robots. Proceedings of the 2014 IEEE International Conference on Robotics and Automation (ICRA); Hong Kong, China, 31 May–7 June 2014; pp. 4540-4547.

25. Tanaka, M.; Sawabe, H.; Nakajima, M.; Ariizumi, R. Redundant Control of a Planar Snake Robot with Prismatic Joints. Int. J. Control Autom. Syst.; 2021; 19, pp. 3475-3486. [DOI: https://dx.doi.org/10.1007/s12555-020-0607-2]

26. Mohammadi, A.; Rezapour, E.; Maggiore, M.; Pettersen, K.Y. Maneuvering control of planar snake robots using virtual holonomic constraints. IEEE Trans. Control Syst. Technol.; 2015; 24, pp. 884-899. [DOI: https://dx.doi.org/10.1109/TCST.2015.2467208]

27. Ariizumi, R.; Matsuno, F. Dynamic analysis of three snake robot gaits. IEEE Trans. Robot.; 2017; 33, pp. 1075-1087. [DOI: https://dx.doi.org/10.1109/TRO.2017.2704581]

28. Kohl, A.M.; Pettersen, K.Y.; Kelasidi, E.; Gravdahl, J.T. Planar path following of underwater snake robots in the presence of ocean currents. IEEE Robot. Autom. Lett.; 2016; 1, pp. 383-390. [DOI: https://dx.doi.org/10.1109/LRA.2016.2517827]

29. Ouyang, W.; Liang, W.; Li, C.; Zheng, H.; Ren, Q.; Li, P. Steering motion control of a snake robot via a biomimetic approach. Front. Inf. Technol. Electron. Eng.; 2019; 20, pp. 32-44. [DOI: https://dx.doi.org/10.1631/FITEE.1800554]

30. Zhang, J.; Chen, Y.; Yang, S.; Li, J.; Song, F.; Gong, Y. Dynamic modeling and analysis of underwater swimming snake robot with soft joint. Proceedings of the 31st International Ocean and Polar Engineering Conference; Rhodes, Greece, 20–25 June 2021.

31. Billah, M.M.; Khan, M.R.; Shafie, A.A. Autonomous flexible snake robot for 3D motion. Proceedings of the 2015 IEEE International Symposium on Robotics and Intelligent Sensors (IRIS); Langkawi, Malaysia, 18–20 October 2015; pp. 105-110.

32. Fu, Q.; Gart, S.W.; Mitchel, T.W.; Kim, J.S.; Chirikjian, G.S.; Li, C. Lateral oscillation and body compliance help snakes and snake robots stably traverse large, smooth obstacles. Integr. Comp. Biol.; 2020; 60, pp. 171-179. [DOI: https://dx.doi.org/10.1093/icb/icaa013] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32215569]

33. Malayjerdi, M.; Akbarzadeh, A. Analytical modeling of a 3-D snake robot based on sidewinding locomotion. Int. J. Dyn. Control; 2019; 7, pp. 83-93. [DOI: https://dx.doi.org/10.1007/s40435-018-0441-z]

34. From, P.J.; Gravdahl, J.T.; Pettersen, K.Y. Vehicle-Manipulator Systems; Springer: Berlin/Heidelberg, Germany, 2016.

35. Rolf, M.; Steil, J.J. Constant curvature continuum kinematics as fast approximate model for the Bionic Handling Assistant. Proceedings of the 2012 IEEE/RSJ International Conference on Intelligent Robots and Systems; Vilamoura-Algarve, Portugal, 7–12 October 2012; pp. 3440-3446.

36. Jazar, R.N. Theory of Applied Robotics: Kinematics, Dynamics, and Control; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2010.

37. Webster, R.H., III; Jones, B.A. Design and kinematic modeling of constant curvature continuum robots: A review. Int. J. Robot. Res.; 2010; 29, pp. 1661-1683. [DOI: https://dx.doi.org/10.1177/0278364910368147]

38. Mustaza, S.M.; Elsayed, Y.; Lekakou, C.; Saaj, C.; Fras, J. Dynamic modeling of fiber-reinforced soft manipulator: A visco-hyperelastic material-based continuum mechanics approach. Soft Robot.; 2019; 6, pp. 305-317. [DOI: https://dx.doi.org/10.1089/soro.2018.0032]

39. Godage, I.S.; Medrano-Cerda, G.A.; Branson, D.T.; Guglielmino, E.; Caldwell, D.G. Dynamics for variable length multisection continuum arms. Int. J. Robot. Res.; 2016; 35, pp. 695-722. [DOI: https://dx.doi.org/10.1177/0278364915596450]

40. Peters, S.T. Handbook of Composites; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2013.

41. Banks, H.T.; Hu, S.; Kenz, Z.R. A brief review of elasticity and viscoelasticity for solids. Adv. Appl. Math. Mech.; 2011; 3, pp. 1-51. [DOI: https://dx.doi.org/10.4208/aamm.10-m1030]

42. Rackl, M. Curve fitting for Ogden, Yeoh and polynomial models. ScilabTEC Conf.; 2015; 2015, pp. 1-11.

43. Winter, H. Viscous dissipation term in energy equations. Calc. Meas. Tech. Momentum Energy Mass Transf.; 1987; 7, pp. 27-34.

44. Witteveen, W.; Pichler, F. On the relevance of inertia related terms in the equations of motion of a flexible body in the floating frame of reference formulation. Multibody Syst. Dyn.; 2019; 46, pp. 77-105. [DOI: https://dx.doi.org/10.1007/s11044-018-09662-0]

45. Antonelli, G. Modelling of underwater robots. Underwater Robots; Springer: Berlin/Heidelberg, Germany, 2018; pp. 33-110.

46. Kanso, E.; Marsden, J.E.; Rowley, C.W.; Melli-Huber, J.B. Locomotion of articulated bodies in a perfect fluid. J. Nonlinear Sci.; 2005; 15, pp. 255-289. [DOI: https://dx.doi.org/10.1007/s00332-004-0650-9]

47. Schmidt-Didlaukies, H.M.; Sørensen, A.J.; Pettersen, K.Y. Modeling of articulated underwater robots for simulation and control. Proceedings of the 2018 IEEE/OES Autonomous Underwater Vehicle Workshop (AUV); Porto, Portugal, 6–9 November 2018; pp. 1-7.

48. Faltinsen, O. Sea Loads on Ships and Offshore Structures; Cambridge University Press: Cambridge, UK, 1993.

49. Spong, M.W.; Hutchinson, S.; Vidyasagar, M. Robot Modeling and Control; John Wiley & Sons: New York, NY, USA, 2006.

50. Chen, Y.-L.; Li, W.-S.; Guo, Q.; Gong, Y.-J. Kinematics Modeling of Soft Manipulator Interacting with Environment Using Segmented Variable Curvature Method. Int. J. Control Autom. Syst.; 2022; 20, pp. 255-267. [DOI: https://dx.doi.org/10.1007/s12555-020-0559-6]

51. Kelasidi, E.; Pettersen, K.Y.; Gravdahl, J.T. Modeling and propulsion methods of underwater snake robots. Proceedings of the 2017 IEEE Conference on Control Technology and Applications (CCTA); Maui, HI, USA, 27–30 August 2017; pp. 819-826.