2.2.1 Kinematic modeling: constant curvature assumption

The robot vertebrae are assumed to locally bend as a constant curvature [Reference Gautreau, Sandoval, Arsicault, Bonnet, Zeghloul, Laribi, Müller and Brandstötter19] meanwhile the curvature is varying from a vertebra to another as they are all different in terms of size and stiffness.

The constant curvature assumption was applied to each subsegment (bending between two vertebrae). Arc geometry formulation was used for vertebra kinematic modeling as described by [Reference Rao, Peyron, Lilge and Burgner-Kahrs23]. Position of $\mathbf{d}(\mathbf{s}_{\mathbf{i}})$ at any point in the arc in function of the base in the $i-1$ frame is given by

(1) \begin{align} \mathbf{d}\left(\mathbf{s}_{\mathbf{i}}\right)=\left[\begin{array}{c} x_{i}\\[3pt] y_{i}\\[3pt] z_{i} \end{array}\right]=\left[\begin{array}{c} r_{i}\cdot \cos \left(\varphi _{i}\right)\cdot \left(1-\cos \left(\theta \left(s_{i}\right)\right)\right)\\[3pt] r_{i}\cdot \sin \left(\varphi _{i}\right)\cdot \left(1-\cos \left(\theta \left(s_{i}\right)\right)\right)\\[3pt] r_{i}\cdot \sin \left(\theta \left(s_{i}\right)\right) \end{array}\right] \end{align}

where $r_{i}=\frac{1}{k_{i}}$ is the radius of curvature. The subsegment curvature $k_{i}$ is defined by Eq. (2), where $, \beta _{i}, \Gamma _{i}$ are, respectively, the curvature of the backbone along $\overrightarrow {x_{i-1}}$ and $\overrightarrow {y_{i-1}}$ . $\varphi _{i}$ is defined by Eq. (3).

(2) \begin{align} k_{i}=\sqrt{\beta _{i}^{2}+\Gamma _{i}^{2}} \end{align}

(3) \begin{align} \varphi _{i}=atan2\left(\Gamma _{i},\beta _{i}\right) \end{align}

Given $\mathbf{d}_{i}$ for $s_{i}=l_{i}$ , the homogeneous transformation matrix from $i-1$ to $i$ is described by [Reference Rao, Peyron, Lilge and Burgner-Kahrs23] as follows:

(4) \begin{align} { }^{i-1}{\mathbf{T}}{}_{i}=\left[\begin{array}{cc} \mathbf{R}_{z}\left(\varphi _{i}\right)\mathbf{R}_{y}\left(\theta _{i}\right)\mathbf{R}_{z}\left(-\varphi _{i}\right) & \ \mathbf{d}_{i}\\ 0 & \ 1 \end{array}\right] \end{align}

The direct geometrical model links the actuator space, giving the pulley angles, to the configuration space and to the cartesian space at the robot tip [Reference Gautreau, Sandoval, Arsicault, Bonnet, Zeghloul, Laribi, Müller and Brandstötter19].

2.2.2 Forces acting on the continuum robot

Given the kinematic modeling, the constant curvature static modeling predicts the robot bending and the robot tip position. The constant curvature static modeling is described by [Reference Yuan, Zhou and Xu21–Reference Rao, Peyron, Lilge and Burgner-Kahrs23]. The four tendons, $m=4$ , are distributed around the backbone at a distance $r_{c}$ from the center of each disk. This distance, $r_{c}$ , can vary as the snake shape (and disks diameters) is varying all along the body [Reference Gautreau, Bonnet, Sandoval, Fosseries, Herrel, Arsicault, Zeghloul and Laribi31]. In this study, $r_{c}$ remain constant. The distribution of tendons on a disk $i$ is described by [Reference Rao, Peyron and Burgner-Kahrs22] as follows:

(5) \begin{align} { }^{i}{\mathbf{p}_{\mathbf{t}}}{}_{i,c}=r_{c,i}\cdot \left[\cos \left(\alpha _{m}\right) \sin \left(\alpha _{m}\right)\ 0\right], with\ \alpha\ \epsilon \left[0,\frac{\pi }{2}, \pi, \frac{3\pi }{2}\right] \end{align}

As defined by [Reference Yuan, Zhou and Xu21], the acting forces in the continuum robot are coming from the gravity, friction, and external forces. The gravity vector defined in the global frame for the mass of the universal compliant joint $(\text{flexible element})i$ , $m_{uc{j_{i}}},$ and the mass of the disk $i$ , $m_{{d_{i}}}$ are given by

(6) \begin{align} { }^{0}{\mathbf{F}_{\mathbf{ucj}}}_{i}=\left[0\ 0\,{m_{ucj}}_{i}\cdot g\ 0\right]^{\mathrm{T}} \end{align}

(7) \begin{align} { }^{0}\mathbf{F}_{\mathbf{d}_{i}}=\left[0\ 0\,{m_{d}}_{i}\cdot g\ 0\right]^{\mathrm{T}} \end{align}

The force $\mathbf{F}_{{\mathbf{d}_{i}}}$ and moments $\mathbf{M}_{{\mathbf{d}}_{i}}$ due to the gravity of the disk are given by

(8) \begin{align} { }^{i-1}\mathbf{F}_{\mathbf{d}_{i}}={ }^{i-1}{\mathbf{T}}{}_{0}\cdot { }^{0}{\mathbf{F}}{}_{i} \end{align}

(9) \begin{align} { }^{i-1}\mathbf{M}_{{\mathbf{d}}_{i}}=\overrightarrow {{{ }^{i-1}{O}{}_{i-1}}{ }^{i-1}{O}{}_{i}}\times { }^{i-1}\mathbf{F}_{\mathbf{d}_{i}} \end{align}

Note that ${ }^{i-1}{\mathbf{T}}{}_{0}$ is the homogeneous transformation matrix. The force ${\mathbf{F}_{\mathbf{ucj}}}_{i}$ and moments ${\mathbf{M}_{\mathbf{ucj}}}_{g}$ due to gravity of the universal compliant joint (flexible element) are given by

(10) \begin{align} { }^{i-1}{\mathbf{F}_{\mathbf{ucj}}}_{i}={ }^{i-1}{\mathbf{T}}{}_{0}\cdot { }^{0}{\mathbf{F}_{\mathbf{ucj}}}_{i} \end{align}

(11) \begin{align} { }^{i-1}{\mathbf{M}_{\mathbf{ucj}}}_{i}=\overrightarrow {{ }^{i-1}{\mathbf{O}}{}_{i-1}{ }^{i-1}{\mathbf{O}_{g}}_{i}}\times { }^{i-1}{\mathbf{F}_{\mathbf{ucj}}}_{i} \end{align}

With ${ }^{i-1}{\mathbf{O}_{\mathbf{ucj}}}_{i}$ the gravity center of the compliant universal joint between two disks. In case of external forces applied at the robot tip, the forces $\mathbf{F}_{\mathbf{ext}}$ and the moment $\mathbf{M}_{\mathbf{ext}}$ induced by the force $\mathbf{F}_{\mathbf{ext}}$ are described by

(12) \begin{align} { }^{v-1}{\mathbf{F}_{\mathbf{ext}}}{}={ }^{v}{\mathbf{T}}{}_{0}\cdot { }^{0}{\mathbf{F}}{}_{e} \end{align}

(13) \begin{align} { }^{v-1}{\mathbf{M}_{\mathbf{ext}}}{}=\overrightarrow {{ }^{v-1}{O}{}_{v-1}{ }^{v-1}{O}{}_{v}}\times { }^{v-1}{\mathbf{F}_{\mathbf{ext}}}{} \end{align}

where ${ }^{0}{\mathbf{F}}{}_{e}$ is the contact force applied to the robot body tip in the frame $R_{0}$ . The driving tendon tensions ${ }^{i-1}{\mathbf{F}}{}_{i,c}$ , triggered by the tendon $c$ at disk $i$ , are the sum of the two driving force vectors $F_{i,c}$ and $F_{i+1,c}$ (Fig. 2(a)) lead by their respective unit vector described as follows and given by [Reference Yuan, Zhou and Xu21] and [Reference Rao, Peyron, Lilge and Burgner-Kahrs23]:

Figure 2. $(a)$ Continuum robot pattern. $(b)$ Kinematic diagram of the continuum robot. $(c)$ Subsegment kinematic diagram.

(14) \begin{align} { }^{i-1}{\mathbf{F}_{i,c}}{}=F_{i,c}\cdot \frac{\overrightarrow {{}^{i-1}{p_{i,c}}{}{}^{i-1}{p_{i-1,c}}{}}}{\overrightarrow {\left\| {}^{i-1}{p_{i,c}}{}{}^{i-1}{p_{i-1,c}}{}\right\| }}+F_{i+1,c}\cdot \frac{\overrightarrow {{}^{i-1}{p_{i,c}}{}{}^{i-1}{p_{i+1,c}}{}}}{\overrightarrow {\left\| {}^{i-1}{p_{i,c}}{}{}^{i-1}{p_{i+1,c}}{}\right\| }} \end{align}

(15) \begin{align} { }^{i-1}{\mathbf{M}_{i,c}}{}={ }^{i-1}{\mathbf{p}_{\mathbf{t}_{i,c}}}{}\times { }^{i-1}{\mathbf{F}_{i,c}}{} \end{align}

Note that when $i=v$ , the force reduces to:

(16) \begin{align} { }^{i-1}{\mathbf{F}_{i,c}}{}=F_{{i_{c}}}\cdot \frac{\overrightarrow {{}^{v-1}{p_{v,c}}{}{}^{v-1}{p_{v-1,c}}{}}}{\overrightarrow {\left\| {}^{v-1}{p_{v,c}}{}{}^{v-1}{p_{v-1,c}}{}\right\| }} \end{align}

The equilibrium equations are defined at the $i-1$ local frame for each subsegment. Then, the resultant moments and forces at vertebra $i=v$ and at point $O_{v-1}$ in the frame $v-1$ are described as follow:

(17) \begin{align} { }^{v-1}{\mathbf{F}}{}_{{O_{v}}}=\sum _{c}^{m}{ }^{v-1}{\mathbf{F}_{v,c}}{}+{ }^{v-1}{\mathbf{F}_{\mathbf{cu}}}{}_{{g_{v}}}+{ }^{v-1}{\mathbf{F}_{\mathbf{d}}}{}_{{g_{v}}} \end{align}

(18) \begin{align} { }^{v-1}{\mathbf{M}}{}_{{O_{v}}}=\sum _{c}^{m}{ }^{v-1}{\mathbf{F}_{v,c}}{}+{ }^{v-1}{\mathbf{M}_{\mathbf{cu}}}{}_{{g_{v}}}+{ }^{v-1}{\mathbf{M}_{\mathbf{d}}}{}_{{g_{v}}} \end{align}

If $i=v$ , the external force $\mathbf{F}_{\mathbf{ext}}$ and moment $\mathbf{M}_{\mathbf{ext}}$ should be included to the equilibrium equations.

2.2.3 Tendons frictions

When the continuum robot is achieving a trajectory, the vertebrae bend according to their stiffness. In the case of discrete interactions (Fig. 2(a)), forces from tendons are locally applied to each disk which led to friction between the tendons and the associated disk. This reflects that the friction coefficient will vary over the continuum robot from the motor to the end effector and during the trajectory execution. To extract the relation between friction coefficient and bending angle, a test bench was specially designed and manufactured as shown in Fig. 3(a).

Figure 3. $(a)$ Friction test bench. $(b)$ Friction coefficient $\mu$ according to deflection angle.

A braided steel-coated cable of 0.5 mm diameter is used for the experiment. The cable goes through a hole of 0.75 mm in the disk. In addition, the cable is guided by two ball bearings to limit the friction and to vary the $\sigma$ angle (Fig. 3(a)). A set of 10 fixed mass $M_{1}$ is set for a range of 12 angles varying from 30 degrees to 85 degrees. The variable mass $M_{2}$ is adjusted for cable balance according to the angle. A linear regression is applied for each angles according to the mass. Thus, a friction coefficient is extracted for each angle (Fig. 3(b)). The friction is described by Coulomb law and defined the coefficient $\mu$ following the Coulomb friction formula as follows:

(19) \begin{align} \mu \left(\sigma \right)=\frac{{F_{f}}_{\sigma }}{{F_{N}}_{\sigma }} \end{align}

With $F_{f}$ the frictional force defined by

(20) \begin{align} {F_{f}}_{\sigma }=\left| M_{1}-M_{2}\right| \end{align}

And the normal force $F_{N}$ defined by

(21) \begin{align} {F_{N}}_{\theta }=M_{1}+M_{2} \end{align}

$M_{1}$ and $M_{2}$ are the fixed mass and the adjusted mass (Fig. 3(a)), respectively. Exponential regression Eq. (22) was the most adapted to describe the friction coefficient $\boldsymbol{\mu }_{i,c}$ in regard of the bending angle ${\unicode[Arial]{x03C3}} .$ Applied to the continuum robot, the coefficient is depicted as following:

(22) \begin{align} \boldsymbol{\mu }_{i,c}=0.689e^{-0.027{\boldsymbol{\sigma }_{i,c}}} \end{align}

$\boldsymbol{\sigma }_{i,c}$ is the local bending angle (Fig. 2(a)) of the vertebrae determined using unit vectors ${ }^{i-1}{\mathbf{U}}{}_{i-1,c}$ and ${ }^{i-1}{\mathbf{U}}{}_{i+1,c}$ :

(23) \begin{align} { }^{i-1}{\mathbf{U}}{}_{i-1,c}=\frac{\overrightarrow {{}^{i-1}{p_{i,c}}{}{}^{i-1}{p_{i-1,c}}{}}}{\overrightarrow {\left\| {}^{i-1}{p_{i,c}}{}{}^{i-1}{p_{i-1,c}}{}\right\| }} \end{align}

(24) \begin{align} { }^{i-1}{\mathbf{U}}{}_{i+1,c}=\frac{\overrightarrow {{}^{i-1}{p_{i,c}}{}{}^{i-1}{p_{i+1,c}}{}}}{\overrightarrow {\left\| {}^{i-1}{p_{i,c}}{}{}^{i-1}{p_{i+1,c}}{}\right\| }} \end{align}

(25) \begin{align} \boldsymbol{\sigma }_{i,c}=\cos ^{-1} \left(\frac{{}^{i-1}{\mathbf{U}}{}_{i-1,c}\cdot {}^{i-1}{\mathbf{U}}{}_{i+1,c}}{\left\| {}^{i-1}{\mathbf{U}}{}_{i-1,c}\right\| \left\| {}^{i-1}{\mathbf{U}}{}_{i+1,c}\right\| }\right) \end{align}

The tension magnitude $F_{i+1,c}$ is updated according to the friction between the disks and the tendons. As depicted in ref. [Reference Yuan, Zhou and Xu21], the resulting tension magnitude can be calculated as

(26) \begin{align} F_{i+1,c}=F_{i,c}-\boldsymbol{\mu }_{i,c }\left| { }^{i-1}{\mathbf{N}_{i,c}}{}\right| \end{align}

where ${ }^{i-1}{\mathbf{N}_{i,c}}{}$ is the normal force applied by tendon tension on the considered disk $i$ defined on the plane $(O_{i},\overrightarrow {x_{i}},\overrightarrow {y_{i}})$ and direction is $(\overrightarrow {p_{ic},O_{i}})$ :

(27) \begin{align} { }^{i-1}{\mathbf{N}_{i,c}}{}={ }^{i-1}{\mathbf{F}_{i,c}}{}-\left({ }^{i-1}{\mathbf{F}_{i,c}}{}\cdot { }^{i-1}{\mathbf{n}}_{i}\right) \end{align}

where ${ }^{i-1}{\mathbf{n}}_{i}$ is the unit normal vector of plane $(O_{i},\overrightarrow {x_{i}},\overrightarrow {y_{i}})$ [Reference Yuan, Zhou and Xu21]. Thus, the net force ${ }^{i-1}{\mathbf{F}_{i,c}}{}$ is updated with the computed tension magnitude $F_{i+1,c}$ . Simulations were performed with the friction coefficient determined experimentally. By applying a range of forces similar to experimental tests ( $F_{0}\in [0,0.7]\,(N)$ ), we show that the friction force is very small and can be neglected during the numerical gait pattern generation (Fig. 4). One note that the friction coefficient is defined in the air without any lubrication. Thus, it is assumed that when the robot is immersed, the coefficient decreases even more due to water lubrication.

Figure 4. Numerical simulation with implementation of friction coefficient. $(a)$ Deformed continuum robot. $(b)$ Comparison of bending angles along the 8 vertebrae (7 angles). $(c)$ Friction coefficient distribution.

2.2.4 Deformation of the compliant vertebra: torsion beam modeling

The continuum compliant robot is designed to reproduce the snake swimming undulations by adopting an elastic behavior similar to that of biological snakes. Thus, the torsion beam diameters are optimally sized in orthogonal directions (both on dorsal/ventral side and lateral side) according to stiffness repartition along biological snakes as described by Gautreau et al. [Reference Gautreau, Bonnet, Laribi and Okada32, Reference Gautreau, Bonnet, Sandoval, Fosseries, Herrel, Arsicault, Zeghloul and Laribi31]. Note that $K_{{y_{i}}}$ is associated to dorsal/ventral stiffness while $K_{{x_{i}}}$ is associated to the lateral stiffness. Thus, the robot bending will vary heterogeneously according to the applied force. This behavior requires static modeling to take into account the tendons tensions and external forces.

(28) \begin{align} \mathbf{K}_{\mathbf{x}}=\left[K_{{x_{1}}},\ldots, K_{{x_{v}}}\right]^{T} \end{align}

(29) \begin{align} \mathbf{K}_{\mathbf{y}}=\left[K_{{y_{1}}},\ldots, K_{{y_{v}}}\right]^{T} \end{align}

The design of each vertebra is then similar to a compliant universal joint where a comprehensive description of the system is provided in ref. [Reference Gautreau, Bonnet, Laribi and Okada32]. The beams compliances are provided by the material flexibility (PA12). The beams are deformed in torsion with small deformations so that each vertebra locally bend to make the continuum robot bending under tendon tensions (Fig, 2(c)). The resulting large deformation of the continuum robot (bending) is allowed by small deformations of each beam in torsions. Small deformations enable to adopt a linear elastic behavior for beams torsion as follows:

(30) \begin{align} { }^{i-1}{\mathbf{M}}{}_{\mathbf{ben}{\mathbf{d}_{i}}}=\left[\begin{array}{c} { }^{i-1}{\tau }{}_{{x_{i}}}\\ { }^{i-1}{\tau }{}_{{y_{i}}}\\ { }^{i-1}{\tau }{}_{{z_{i}}} \end{array}\right] \end{align}

where ${ }^{i-1}{\mathbf{M}}{}_{\mathbf{ben}{\mathbf{d}_{i}}}$ is the torsion moment induced by the torsional stiffness of the torsion beam. ${ }^{i-1}{\tau }{}_{{x_{i}}}$ and ${ }^{i-1}{\tau }{}_{{y_{i}}}$ are the torque obtained from the vertebrae length $l_{i}$ and $\Gamma _{i}$ and $\beta _{i}$ , respectively, the curvature around $y$ axis and $x$ axis. Note that the curvatures $\Gamma _{i}$ and $\beta _{i}$ are the optimized parameters:

(31) \begin{align} { }^{i-1}{\tau }{}_{{x_{i}}}=K_{{x_{i}}}\cdot \theta _{{x_{i}}}=K_{{x_{i}}}\cdot l_{i}\cdot \Gamma _{i} \end{align}

(32) \begin{align} { }^{i-1}{\tau }{}_{{y_{i}}}=K_{{y_{i}}}\cdot \theta _{{y_{i}}}=K_{{y_{i}}}\cdot l_{i}\cdot \beta _{i} \end{align}

(33) \begin{align} { }^{i-1}{\tau }{}_{{z_{i}}}={ }^{i-1}{T}{}_{i}\cdot \left[0\ 0\ G_{i}\cdot I_{{0_{i}}}\cdot \frac{\varepsilon _{i}}{l_{i}}\ 0\right]^{T} \end{align}

where $l_{i}$ is defined by $l_{i}=L_{{11_{i}}}+L_{{12_{i}}}$ . $\varepsilon _{i}$ is the twisting angle of vertebra $i$ . ${ }^{i-1}{\tau }{}_{{z_{i}}}$ is the twist moment induced around z axis described by [Reference Yuan, Zhou and Xu21] and [Reference Rao, Peyron, Lilge and Burgner-Kahrs23]. In the case of the present robot, the twisting moment has not been taken into account in the design phase. Thus, ${ }^{i-1}{\tau }{}_{{z_{i}}}$ is neglected to avoid any twisting moment in numerical simulations.

2.2.5 Archimedes force modeling

The robot is intended to be a locomotor swimming on the $(O_{0},\overrightarrow {x_{0}},\overrightarrow {y_{0}})$ plane. Assuming the skeleton is immersed, the Archimedes force is opposed to gravity and acts on the robot bending angle. For a configuration in which the robot is horizontally setup (Fig. 5), the Archimedes thrust (force and moment) are defined according to the volume of the vertebrae $V_{i}$ , as follows:

Figure 5. Experimental setup of the robot subject to gravity and Archimedes force.

(34) \begin{align} { }^{i-1}{\mathbf{f}}{}_{\mathbf{a},0}=\left(\rho _{w}V_{i}\right)\cdot \left[\begin{array}{cccc} 0 & -1 & 0 & 0 \end{array}\right]^{T} \end{align}

(35) \begin{align} { }^{i-1}{\mathbf{F}}{}_{\mathbf{a},i}={ }^{i-1}{\mathbf{T}_{0}}{}\cdot { }^{i-1}{\mathbf{f}}{}_{\mathbf{a},0} \end{align}

(36) \begin{align} { }^{i-1}{\mathbf{M}}{}_{\mathbf{a},i}=\overrightarrow {{ }^{i-1}{O}{}_{i-1}{ }^{i-1}{O_{g}}_{i}}\cdot { }^{i-1}{\mathbf{F}}{}_{\mathbf{a},i} \end{align}

Thus, the net force and the net moment applied at vertebra $i=v$ are set to:

(37) \begin{align} { }^{v-1}{\mathbf{F}}{}_{{O_{v}}}=\sum _{c}^{m}{ }^{v-1}{\mathbf{F}_{v,c}}{}+{ }^{v-1}{\mathbf{F}_{\mathbf{cu}}}{}_{{g_{v}}}+{ }^{v-1}{\mathbf{F}_{\mathbf{d}}}{}_{{g_{v}}}+{ }^{v-1}{\mathbf{F}}{}_{\mathbf{a},v} \end{align}

(38) \begin{align} { }^{v-1}{\mathbf{M}}{}_{{O_{v}}}=\sum _{c}^{m}{ }^{v-1}{\mathbf{F}_{v,c}}{}+{ }^{v-1}{\mathbf{M}_{\mathbf{cu}}}{}_{{g_{v}}}+{ }^{v-1}{\mathbf{M}_{\mathbf{d}}}{}_{{g_{v}}}+{ }^{v-1}{\mathbf{M}}{}_{\mathbf{a},v} \end{align}

Note that in the case of horizontal static modeling, the gravity is along y axis $.$ A comparison of end effector position between numerical optimization and experiment for a given trajectory will be given in section 4.

2.2.6 Equilibrium equations

Given the resultant forces and net moments at vertebra $i=v$ Eqs. (17)–(18), Eqs. (37)–(38), the equilibrium equations at vertebrae $i$ for $i\in \{0,..,v\}$ are described by [Reference Yuan, Zhou and Xu21] as follow:

(39) \begin{align} { }^{i-1}{\mathbf{F}}{}_{{O_{i-1}}}=\sum _{c}^{m}{ }^{i-1}{\mathbf{F}_{i,c}}{}+{ }^{i-1}{\mathbf{F}_{\mathbf{cu}}}{}_{{g_{i}}}+{ }^{i-1}{\mathbf{F}_{\mathbf{d}}}{}_{{g_{i}}}+{ }^{i-1}{\mathbf{F}}{}_{\mathbf{a},i}+{ }^{i-1}{\mathbf{F}}{}_{{O_{i}}} \end{align}

(40) \begin{align} { }^{i-1}{\mathbf{M}}{}_{{O_{i-1}}}=\sum _{c}^{m}{ }^{i-1}{\mathbf{F}_{i,c}}{}+{ }^{i-1}{\mathbf{M}_{\mathbf{cu}}}{}_{{g_{i}}}+{ }^{i-1}{\mathbf{M}_{\mathbf{d}}}{}_{{g_{i}}}+{ }^{v-1}{\mathbf{M}}{}_{\mathbf{a},i}+{ }^{i-1}{\mathbf{M}}{}_{\mathbf{l}{O_{i}}}+{ }^{i-1}{\mathbf{M}}{}_{{O_{i}}} \end{align}

The discretized forces and moments ${ }^{i-1}{\mathbf{F}}{}_{{O_{i}}}$ and ${ }^{i-1}{\mathbf{M}}{}_{{O_{i}}}$ are the lumped forces and moments at point $O_{i}$ in the frame $i-1$ , and defined as follows:

(41) \begin{align} { }^{i-1}{\mathbf{F}}{}_{{O_{i}}}={ }^{i-1}{\mathbf{T}}{}_{i}\cdot { }^{i}{\mathbf{F}}{}_{{O_{i}}} \end{align}

(42) \begin{align} { }^{i-1}{\mathbf{M}}{}_{{O_{i}}}={ }^{i-1}{\mathbf{T}}{}_{i}\cdot { }^{i}{\mathbf{M}}{}_{{O_{i}}} \end{align}

And ${ }^{i-1}{\mathbf{F}}{}_{{O_{i}}}$ induce a moment ${ }^{i-1}{\mathbf{M}}{}_{\mathbf{l}{O_{i}}}$ in the frame $i-1$ according to the vertebra length $i$ :

(43) \begin{align} { }^{i-1}{\mathbf{M}}{}_{\mathbf{l}{O_{i}}}=\overrightarrow {{ }^{i-1}{O}{}_{i-1}{ }^{i-1}{O}{}_{i}}\cdot { }^{i-1}{\mathbf{F}}{}_{{O_{i}}} \end{align}

As described in the previous subsection, the vertebrae have a known compliance allowing the continuous robot to bend in two planes. This results in the following equilibrium equation as depicted by [Reference Yuan, Zhou and Xu21]:

(44) \begin{align} { }^{i-1}{\mathbf{M}}{}_{{O_{i-1}}}={ }^{i-1}{\mathbf{M}}{}_{\mathbf{ben}{\mathbf{d}_{i}}} \end{align}

The static model of the robot was solved using the numeric method trust-region-recursive algorithm (fsove in MATLAB) [Reference Rao, Peyron, Lilge and Burgner-Kahrs23].

4.1.1 Experimental test bench in air

Two experimental test benches were used for the experiments. For a numerical/experimental coupling, the robot was tested vertically outside of water. The actuator (“neck”) (Fig. 8) was fixed on the test bench, and the tip of the robot (“head”) was mobile and equipped with four force sensors. The motion was recorded through motion capture. The head and the neck were equipped with a cluster of three passive reflective markers (Fig. 8), whereas disks of vertebrae were equipped with only one reflective marker. Velocity and orientation are computed from MOCAP using MATLAB. Tip of the robot is equipped with four force sensors for each of the four actuated tendons. Forces were obtained using Phidget^® control panel at a recording frequency of 128 milliseconds.

Figure 8. Experimental test bench and prototype in air with the schematic set-up.

The work of the combined forces is computed at each position of the imposed trajectory. The general definition of the work $W$ of a force in a plane with an angle $\theta _{w}$ at each point $j$ is described as

(55) \begin{align} W=F\cdot \Delta d\cdot \cos \left(\theta _{w}\right) \end{align}

In the case of the continuum robot, the equivalent force is the combination of the 4 forces measured at the tip of the 4 tendons as shown in Fig. 9.

Figure 9. Equivalent force $\boldsymbol{F}_{\boldsymbol{eq}}$ applied by tendons at the tip of the robot.

The equivalent force computed at the tip of the last disk $O_{v}$ due to tendon force is defined as

(56) \begin{align} \mathbf{F}_{\mathbf{e}{\mathbf{q}_{j}}}=\left[\begin{array}{c} {F_{eq}}_{{{x_{0}}_{j}}}\\[4pt] {F_{eq}}_{{{y_{0}}_{j}}}\\[4pt] {F_{eq}}_{{{z_{0}}_{j}}} \end{array}\right]=\left[\begin{array}{c} \sum _{c=1}^{4}{F_{{x_{0}}}}_{c,j}\\[4pt] \sum _{c=1}^{4}{F_{{y_{0}}}}_{c,j}\\[4pt] \sum _{c=1}^{4}{F_{{z_{0}}}}_{c,j} \end{array}\right] \end{align}

With $c\in [1,4]$ tendons and $j\in [1,n]$ positions on the trajectory. The total works along each axis, according to the equivalent force $\mathbf{F}_{\mathbf{e}{\mathbf{q}_{j}}}$ and the displacement $\Delta \mathbf{d}_{(j,j+1)}$ , result in an equivalent work at each point of the trajectory $\mathbf{W}_{\mathbf{e}{\mathbf{q}_{j}}}$ :

(57) \begin{align} \mathbf{W}_{\mathbf{e}{\mathbf{q}_{j}}}=\left[\begin{array}{c} {F_{eq}}_{{{x_{0}}_{j}}}\\[2pt] {F_{eq}}_{{{y_{0}}_{j}}}\\[2pt] {F_{eq}}_{{{z_{0}}_{j}}} \end{array}\right]\cdot \left[\begin{array}{c} \Delta d_{{{x_{0}}_{j}}}\\[2pt] \Delta d_{{{y_{0}}_{j}}}\\[2pt] \Delta {d_{{z_{0}}}}_{j} \end{array}\right]={F_{eq}}_{{{x_{0}}_{j}}}\cdot \Delta d_{{{x_{0}}_{j}}}+{F_{eq}}_{{{y_{0}}_{j}}}\cdot \Delta d_{{{y_{0}}_{j}}}+{F_{eq}}_{{{z_{0}}_{j}}}\cdot \Delta d_{{{z_{0}}_{j}}} \end{align}

The robot is manually set vertically to begin the imposed trajectory. The starting position is identical to numerical simulation (tip of the robot is at coordinate [0,0], respectively, along $x$ axis and $z$ axis). Each trajectory is repeated three times, and the starting position is not modified between trials. The trajectories N°1 and N°2 were performed on the test bench.

4.1.2 Experimental test bench in water

In the second case, the robot is vertically immersed underwater. The motion was recorded through GoPro^® hero 9 video cameras with the linear mode to avoid image distortion and placed in front of the robot outside of water (Fig. 10). Motion of the robot is extracted using the markerless DeepLabCut software. The software was trained using the resnet50 neuronal network with three labeled videos. In order to allow the deep learning software to track the positions of specific points, we labeled several videos from our experiments. We specified which points to follow for each video. Consequently, any video taken under the same experimental conditions as the labeled ones can be analyzed to track these specific points. Each video is labeled on 40 frames. Velocity and orientation are computed using MATLAB. The force sensors were coated using epoxy to protect strain gauges from water. Forces were obtained using Phidget^® control panel at a recording frequency of 128 milliseconds. The actuator remains outside of water as it is not waterproof. The work was computed based on measured force and displacements $\Delta d_{x},\Delta d_{y}$ , and $\Delta d_{z}$ .

Figure 10. Experimental test bench and prototype in water with the schematic set-up.

In this study, the displacement can only be measured on a plane $(O_{0},\overrightarrow {x_{0}},\overrightarrow {z_{0}})$ . The equivalent work is computed using Eq. (54). All the continuum flexible modules were immersed. Similarly, to the previous experimental test bench and the experimentations, the starting position of the robot is manually set vertically, and the robot position is not modified from a played trajectory to another one. The trajectories N°1 and N°2 were carried out on the test bench.

4.1.3 Experimental test bench in water and horizontally

The last experiments focus on the continuum robot horizontally oriented as the snake robot will swim on the $(O_{0},\overrightarrow {x_{0}},\overrightarrow {z_{0}})$ plane. The Archimedes force is considered for the continuum flexible module.

The actuator is partially immersed, and a waterproof case was design to protect the Dynamixel^® 2XC-430W-250T servomotor. Two GoPro Hero 9 cameras were positioned orthogonally to record the motion. The camera 1 is positioned to record lateral undulations, and camera 2 is positioned to record the dorsal/ventral undulations, as shown in Fig. 11. Two virtual markers were set on the labeled videos (Fig. 11) using DeepLabCut software to extract the 3D movement of the tip of the robot based on the same local frame, by combining the positions of each video. The trajectory N°3 was performed for the corresponding setup.

Figure 11. Prototype set-up dedicated for horizontal experiments.

4.2.1 Planar trajectory

The obtained trajectory represents the location of the snake head undulating in function of the neck on a $(O_{0},\overrightarrow {x_{0}},\overrightarrow {y_{0}})$ plane, corresponding to the lateral side oscillations. Experimental positions of the robot end-effector along $\overrightarrow {x}$ , $\overrightarrow {y}$ , and $\overrightarrow {z}$ are similar to numerical ones as depicted in Fig. 12 and Fig. 13 which confirm the truthfulness of the proposed static modeling for the given continuum robot design.

Figure 12. Position comparison between numerical result, experimental in air, and experimental in water results.

Figure 13. (a) Position errors between numerical and experimental results in air. (b) Position errors between numerical and experimental results in water.

Figure 14. $(a)$ Relative errors between experimental forces in air and numeric forces. $(b)$ Relative errors between experimental forces in water and numeric forces.

Figure 15. Planar trajectory: in the air $(a)$ Tendons forces $(b)$ Displacements $(c)$ Work; Underwater $(d)$ Tendons forces $(e)$ Displacements $(f)$ Work.

Figure 16. Vertical 3D trajectory of the robot tip positions along $x$ , $y$ and $z$ axis.

Figure 17. (a) Position errors between numerical and experimental results in air. (b) Position errors between numerical and experimental results in water.

Figure 18. Computed work for vertical 3D trajectory. $(a)$ Work on the air. $(b)$ Work on water.

Figure 19. Horizontal 3D trajectory end effector positions along $x$ , $y$ and $z$ axis.

Figure 20. Position errors between numerical and experimental results in water.

Remark 1: The error positioning can come from several reasons, the open loop control, the manual initial positioning, as well as the MOCAP analysis, whereas still in ratio to the biological study and mimicry. Snake swimming undulations are approximative and do not require high accuracy (positioning error less than 10 mm is acceptable in this study).

Measured experimental forces present a behavior closely similar to those obtained during numeric optimization. Relative errors in the air between experimental forces and numeric forces are smaller than those obtained underwater (Fig. 14). Errors are slightly high but remains acceptable according to experimental conditions.

In fact, the magnitude force order is quite similar between air, water, and numeric environment results. A comparison of forces (air, water, and numeric environment) was performed to assess the force accuracy. This latter shows that the force error is similar between experimentations in the air and the numeric and the experimental underwater vs the numerical one. In both cases, one notes that during the experimentation, force sensors mounted on nonactuated cables detect a small force variation, which can be explained by the internal friction between cables and disks. These forces are not considered in the static model and as well minimized underwater as expected due to the water effect acting as lubricant.

Remark 2: Force sensors are slightly sensitive to temperature. In experimental conditions, temperature can vary between 20°C and 35°C, which alter the force measurement. Tendons were manually pretensioned, which can lead to a variation of measured forces. Note that these results are preliminary and promising to open the way to new developments.

Given the force and the displacement for a gait pattern, the computed work shows the energy consumption in function of position (Fig. 15(c) (f)). It appears that the work is higher underwater than in air. One assumes water resistance plays a significant role, and water force should be considered in the next experimentations. Work is high in transition phases from an extremum to another. The work is null at the vertical position as it is a resting position. Work shows where and when the robot is consuming energy. This could be a major criterion for robot design to adjust the robot stiffness and tendon routine in the objective of energy consumption minimization.

4.2.2 Vertical 3D trajectory

The 3D trajectory represents the head of the snake undulating in function of the neck on a volume. The positions of the robot tip (end-effector) shown in Fig. 16 follow faithfully the optimized trajectory, and the deviations (Fig. 17) are still acceptable for biomimicry motion (Fig. 7).

The obtained results demonstrate the robot capability to reproduce an imposed gait pattern in a volume. Experimental forces were measured in the air and in water to compute equivalent work (Fig. 18). The work order of magnitude appears higher for underwater motions than for air. The water forces are additional forces that should be considered for fluid–structure interaction.

4.2.3 Horizontal 3D trajectory

This section represents the head of the snake undulating in function of the neck on a volume and horizontally. This configuration was close to biology as the prototype was ¾ immersed. The neck was fixed, and the head was moving. The recorded 3D movement shows the robot achieve faithfully the imposed motion as shown in Fig. 19 and Fig. 20.