2.1. Static theoretical model of WHAM

The force characteristics of the artificial muscles are related to the wall thickness, the geometrical features of end plugs, the friction, and other factors. In the previous modeling approach based on energy conservation, the effects of wall thickness and friction were ignored and the rubber hose was assumed to be an ideal cylinder. Considering the effect of water depth on artificial muscles, the output force of the artificial muscles in the previous theoretical model $N_0$ can be expressed as [Reference Ching-Ping and Hannaford21]

(1) \begin{equation}{N_0} = \frac{{\pi {D^2}(p - {p_e})(3{{\rm{cos}}^2}\alpha - 1)}}{4}\end{equation}

To take the factor of elastic force produced by the rubber hose into account, the following assumptions are made [Reference Hocking and Wereley22,Reference Kothera, Jangid, Sirohi and Wereley23]:

Figure 1.

Force analysis of WHAM considering elasticity.

Based on of the above assumptions, the force analysis of WHAM considering elastic force is shown in Fig. 1, and the force equilibrium equations are deduced as

(2) \begin{equation}\left\{ \begin{array}{l}{F_g} + \pi {R^2}(p - {p_e}) = {\sigma _z}A + T\cos \alpha \\ \\[-7pt] (p - {p_e})lR = {\sigma _z}tl + nT\sin \alpha\end{array} \right.\end{equation}

Considering the impact of friction between the rubber hose and the braided net and between the fibers, the empirical friction coefficient $K_f$ are introduced into the model [Reference Woods, Kothera and Wereley24]. The output force of the improved theoretical model N taking into account elasticity and friction can be expressed as:

(3) \begin{equation}\begin{array}{l}N = (1 - {k_f})(\dfrac{{({p} - {p_e})(3{l^2} - {b^2})}}{{4\pi{n^2}}} + \dfrac{{VE}}{l}(\dfrac{l}{{{l_0}}} - 1)+\dfrac{{{l^2}Et}}{{n\sqrt {{b^2} - {l^2}} }} - \dfrac{{{l^2}Et}}{{2\pi {n^2}{R_0}}})\end{array}\end{equation}

Under the same input pressure, the comparison of the relationship between output force and contraction rate ( $\frac{{l_0}-l}{{{l_0}}}$ ) of WHAM in the previous model and the improved model is shown in Fig. 2. When the contraction rate is 0, the output force is smaller in the improved model. And when the output force is zero, the contraction rate in the previous model is 0.330, and that in the improved model is 0.268.

Figure 2.

Comparison of the relationship between output force and contraction rate of WHAM in the previous model and the improved model under the same input pressure.

2.2. Kinematic model of the continuous manipulator

Figure 3 shows the structural model of the continuous manipulator studied in this paper. Each rigid-flexible segment consists of a bending section and two non-deformable sections, including three vertical WHAMs and three horizontal ones. The paper refer to the two non-deformable sections connecting the bending sections of the upper and the lower two segments as a connecting plate. In each connecting plate, three horizontal WHAMs and 3D-printed components make up a variable-diameter gripper. When the transverse muscles are pressurized, the diameter of the gripper decreases to grasp objects. The manipulator can bend in different positions and attitudes by inputting different combinations of pressure to the vertical muscles. Table I shows the changes in the position and attitude of the manipulator under different input pressures of WHAMs. And by the cooperation of the vertical and horizontal muscles, the manipulator can swallow and disgorge objects in imitation of the action of snakes. The characteristics of the rigid-flexible structural make the manipulator keep flexible and have specific stiffness to improve its load capacity. See more details about the manipulator in this paper [Reference Wang, Xu, Yu, Li, Yang, Chen, Chen, Zhang and Zhou20].

Figure 3.

Structural model of the manipulator.

Table I.

Manipulator input pressure and the corresponding position and attitude.

The PCC model can be applied to manipulators with non-deformable sections between the bending sections [Reference Gong, Cheng, Chen, Sun, Fang, Hu, Xie, Wang and Wen25]. In the PCC model, the bending section composed of WHAMs can be considered as joined segments of constant curvature (CC). The following two assumptions are made to predict the bending deformation and spatial position based on the PCC model [Reference Jones and Walker26]:

For the continuous manipulator in this paper with n segments, the kinematic model is established as shown in Fig. 4(a). According to the PCC model, the paper set up reference coordinate systems in various segments to describe the changes of the position and attitude through the homogeneous transformation matrices. All reference coordinate systems are established by the same standard, their x-axis all passes through a WHAM on the same side. The homogeneous transformation matrices are represented by three configuration space parameters, which are curvature radius $r_i$ , bending angle $\theta_i$ , and the angle between the bending plane and the positive $X_i$ -axis $\phi_i$ (i =1,2…,n). The three WHAMs in each segment are defined as WHAM1, WHAM2, and WHAM3. Then for the WHAMj (j=1,2,3) of Segmenti, its length, input pressure, output force, and constraint reaction force are defined as $l_{ij}$ , $p_{ij}$ , $N_{ij}$ , and $T_{ij}$ respectively. Take all the structure between $O_i$ disk and $O_n$ disk as a whole part, Fig. 4(b) shows its force analysis when underwater working with load force. The whole part is affected by the load force $F_N$ , the gravity of each connecting plate $G_i$ under water (considering buoyancy), the output forces of the WHAMs $N_{ij}$ , and the constraint reaction forces of the WHAMs $T_{ij}$ . The force balance equations can be listed as follows:

(4) \begin{equation}\left\{ \begin{array}{l}({F_N} + \sum\limits_{k = i}^n {{G_k}} )\cos ({\theta _{ii}}) = \sum\limits_{j = 1}^3 {{N_{ij}}} \\({F_N} + \sum\limits_{k = i}^n {{G_k}} )\sin ({\theta _{ii}}) = \sum\limits_{j = 1}^3 {{T_{ij}}} \\\end{array} \right.\end{equation}

The ends of the WHAMs are installed perpendicular to the connecting plates, so the output force $N_{ij}$ acts perpendicular to connecting plate $O_i$ and the constraint reaction force $T_{ij}$ acts parallel to connecting plate $O_i$ . The component forces of the gravity and load force perpendicular to the connecting plates are balanced with the output forces. Since $N_{ij}$ is determined by the input pressure and the length of the WHAM, the first force balance equation determines the length of the WHAMs when the input pressure is given. The component forces of gravity and load force parallel to the connecting plates are balanced with the constrained reaction forces. Because the forces are perpendicular to the ends of the WHAMs, the second force balance equation does not affect the length of the WHAMs within a specific range. The bending angles $\theta_i,\theta_{ii}$ are also the minimum dihedral angles between the connecting plates. The value of the minimum dihedral angles can be expressed by the geometric relationship with the length of the WHAMs. While $\theta_{ii}$ is the minimum dihedral angle between the $O_i$ plate and the $O_0$ plate on the horizontal plane, when $0 \le {\theta _{ii}} \le 0.5\pi$ , the range of $\theta _{ii}$ can be expressed as

(5) \begin{equation}{\theta _{ii}} \ge {\theta _1} + {\theta _2} + ... + {\theta _i}\end{equation}

Substitute Eq. (5) into Eq. (4). When the input pressure is given, the solution interval of the length of the WHAMs can be expressed as

(6) \begin{equation}\left\{ \begin{array}{l}{N_{11}} + {N_{12}} + {N_{13}} = ({F_N} + \sum\limits_{k = 1}^n {{G_{k}}} )\cos ({\theta _1})\\{N_{i1}} + {N_{i2}} + {N_{i3}} \le ({F_N} + \sum\limits_{p = i}^n {{G_{p}}})\cos (\sum\limits_{q = 1}^n {{\theta _q}} )\\\end{array} \right.\\(i = 2,3,...,n)\end{equation}

where ${\theta _i} = f({l_{ij}}){\rm{\;\;\;\;\;\;\;\;\;\;\;\;\;}}$

$\qquad\!\!\!\!{N_{ij}}= g({l_{ij}},{p_{ij}}){\rm{\;\;\;\;\;\;}}(i = 1,2,...,n;\,j = 1,2,3)$

Figure 4.

Kinematic model and force analysis of the underwater manipulator with load.

3.1. Forward kinematic solution

For linear elastic statically indeterminate systems with unknown loading processes, the principle of minimum strain energy is valid as long as the system is in a stable state under the final load. The system has the condition that the strain energy of the whole system is minimized. Therefore, the minimum elastic strain energy of WHAMs in the process of pressurization and loading is set as the optimization function to obtain the solution of the length of the WHAMs. According to the strain energy function of silicon rubber material described by the neo-Hookean model, the strain energy function of WHAMs can be expressed as [Reference Sheikhi, Shojaeifard and Baghani27]

(7) \begin{equation}{U_{ij}} = {C_{10}}(\frac{{{l_{ij}}^2}}{{l_0^2}} + \frac{{{b^2} - {l_{ij}}^2}}{{4{\pi ^2}{n^2}R_0^2}} - 2){\rm{\;\;\;}}(i = 1,2,...,n;j = 1,2,3)\end{equation}

The original length of WHAMs before pressurized is specified as the rated maximum working length. The optimization algorithm can be expressed as

(8) \begin{equation}\begin{array}{l}\textbf{X} = {\left[ {{l_{11},l_{12},l_{13},l_{21},...,l_{ij}}} \right]^T}{\rm{\;\;\;\;}}{\rm{\;\;\;}}(i = 1,2,...,n;j = 1,2,3)\\ \\[-7pt] \min f(\textbf{X}) = \sqrt[2]{{\dfrac{{\sum\limits_{i{\rm{ = }}1}^n {\sum\limits_{j{\rm{ = }}1}^3 {{U_{ij}}} } }}{n}}}\\ \\[-7pt] {\rm{s}}.{\rm{t}}.{\rm{\;\;}}h(\textbf{X}) = ({F_N} + \sum\limits_{k = 1}^n {{G_k}} )\cos ({\theta _1}) - \sum\limits_{j = 1}^3 {{N_{1j}}} {\rm{ = }}0\\ \\[-7pt] {\rm{\;\;\;\;\;\;\;}}g^{1}(\textbf{X}) =\sum\limits_{j = 1}^3 {{N_{ij}}} - ({F_N} + \sum\limits_{p = i}^n {{G_p}} )\cos (\sum\limits_{q = 1}^i {{\theta _q}} ) \le 0\\ \\[-7pt] {\rm{\;\;\;\;\;\;\;}}g^{2}(\textbf{X}) = {N_{ij}}\ge 0{\rm{\;\;\;\;}}\end{array}\end{equation}

To solve the configuration space parameters, the bending section of each segment is divided into m infinitesimals as shown in Fig. 4(c). The length of three WHAMs in infinitesimals is ${l_{ij}^s}=\frac{l_{ij}}{m}$ . The cross-sectional radius of each infinitesimal $d_i ^s$ varies gradually with s and can be expressed by the length of the transverse muscles. According to the geometric relationship, the configuration space parameters in each infinitesimal $r_i ^s,\theta_i ^s,\phi_i ^s$ can be expressed as

(9) \begin{equation}\left\{ \begin{array}{l}r_i^s = \dfrac{{(l_{i1}^s + l_{i2}^s + l_{i3}^s)d_i^s}}{{2\sqrt {\sum\limits_{j = 1}^3 {{{\left( {l_{ij}^s} \right)}^2}} - l_{i1}^sl_{i2}^s - l_{i2}^sl_{i3}^s - l_{i1}^sl_{i3}^s} }}\\[10pt]\theta _i^s = \dfrac{{2\sqrt {\sum\limits_{j = 1}^3 {{{\left( {l_{ij}^s} \right)}^2}} - {{\left( {l_{i1}^s} \right)}^2}{{\left( {l_{i2}^s} \right)}^2} - {{\left( {l_{i2}^s} \right)}^2}{{\left( {l_{i3}^s} \right)}^2} - {{\left( {l_{i1}^s} \right)}^2}{{\left( {l_{i3}^s} \right)}^2}} }}{{3d_i^s}}\\[10pt]\phi _i^s = \left\{ \begin{array}{l}{\arctan}\dfrac{{\sqrt 3 (l_{i2}^s + l_{i3}^s - 2l_{i1}^s)}}{{3(l_{i2}^s - l_{i3}^s)}}(l_{i2}^s \ne l_{i3}^s)\\[10pt]- 0.5\pi(l_{i2}^s = l_{i3}^s < l_{i1}^s)\\[5pt]0.5\pi(l_{i2}^s = l_{i3}^s > l_{i1}^s)\end{array} \right.\\\end{array} \right.(i=1,2,...,n;\ j=1,2,3;\ s=1,2,...,m)\end{equation}

The coordinate transformation matrix of each infinitesimal ${}_i^sT$ can be expressed by the configuration space parameters as

(10) \begin{equation}\begin{array}{l}{}_i^sT = {\rm{Rot}}(z,\phi _i^s){\rm{Rot}}(y,\theta _i^s){\rm{Trans}}(r_i^s(1 - {\rm{c}}\theta _i^s),0,r_i^s{\rm{s}}\theta _i^s){\rm{Rot}}(z, - \phi _i^s)\\ \\[-7pt] {\rm{ = }}\left[ {\begin{array}{*{20}{c}}{{{\rm{cos}}^2}\phi _i^s({\rm{cos}}\theta _i^s - 1) + 1}&{{\rm{sin}}\phi _i^s{\rm{cos}}\phi _i^s({\rm{cos}}\theta _i^s - 1)}&{{\rm{cos}}\phi _i^s{\rm{sin}}\theta _i^s}&{r_i^s{\rm{cos}}\phi _i^s(1 - {\rm{cos}}\theta _i^s)}\\ \\[-7pt] {{\rm{sin}}\phi _i^s{\rm{cos}}\theta _i^s({\rm{cos}}\theta _i^s - 1)}&{{{\rm{cos}}^2}\phi _i^s(1 - {\rm{cos}}\theta _i^s) + {\rm{cos}}\theta _i^s}&{{\rm{sin}}\phi _i^s{\rm{sin}}\theta _i^s}&{r_i^s{\rm{sin}}\phi _i^s(1 - {\rm{cos}}\theta _i^s)}\\ \\[-7pt] { - {\rm{cos}}\phi _i^s{\rm{sin}}\theta _i^s}&{ - {\rm{sin}}\phi _i^s{\rm{sin}}\theta _i^s}&{{\rm{cos}}\theta _i^s}&{r_i^s{\rm{sin}}\theta _i^s}\\ \\[-7pt] 0&0&0&1\end{array}} \right]\\ \\[-7pt] {\rm{\;\;\;}}(i = 1,2,...,n;\ s = 1,2,...,m)\end{array}\end{equation}

The coordinate transformation matrices of the non-deformable section at both ends of each segment ${{}_a^iT,{}_c^iT}$ are represented by transformation matrices. Then the terminal position and attitude matrix $^nT$ can be expressed as

(11) \begin{equation}\begin{array}{l}^nT = \prod\limits_{i = 1}^n {{}_a^iT} {}_b^iT{}_c^iT\\\end{array}\end{equation}

(12) \begin{equation}\left\{ \begin{array}{l}{}_b^iT = \mathop \prod \limits_{s = 1}^m {}_i^sT{\rm{\;\;\;\;\;\;\;}}\\{}_a^iT = {}_c^iT = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}1&0&0&0\\ \\[-7pt]0&1&0&0\\ \\[-7pt]0&0&1&{{l_w}}\\ \\[-7pt]0&0&0&1\end{array}} \right]\end{array} \right.(i = 1,2,...,n)\end{equation}

For the terminal position and attitude matrix

\[\begin{array}{l}^nT = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{{e_1}}&{{e_4}}&{{e_7}}&{{e_{10}}}\\ \\[-7pt]{{e_2}}&{{e_5}}&{{e_8}}&{{e_{11}}}\\ \\[-7pt]{{e_3}}&{{e_6}}&{{e_9}}&{{e_{12}}}\\ \\[-7pt]0&0&0&1\end{array}} \right]\\ \\[-7pt]\end{array}\]

where ${e_k} = g({l_{ij}})\qquad\qquad\!\!\!\! (k = 1,2,...,12)$

${\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;}{l_{ij}} = f({p_{ij}},{N_{ij}})\qquad(i = 1,2,...,n;\ j = 1,2,3)$

Its terminal position and attitude ${}^n\gamma$ can be expressed as

(13) \begin{equation}{}^n\gamma{\rm{ = }}\left[ {\begin{array}{*{20}{c}}x\\ \\[-7pt]y\\ \\[-7pt]z\\ \\[-7pt]\delta \\ \\[-7pt]\beta \\ \\[-7pt]\psi\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}{{e_{10}}}\\ \\[-7pt]{{e_{11}}}\\ \\[-7pt]{{e_{12}}}\\ \\[-7pt]{{\rm{atan2}}({e_6},{e_9})}\\ \\[-7pt]{{\rm{atan2}}( - {e_3},\sqrt {e_1^2 + e_2^2} )}\\ \\[-7pt]{{\rm{atan2}}({e_2},{e_1})}\end{array}} \right]\end{equation}

$\delta,\beta,\psi$ are roll angle, pitch angle, and yaw angle, respectively. For the continuum manipulator with three segments studied in this paper, the solution process is shown in Fig. 5. $p_{h1}, p_{h2}, p_{h3}$ is the input pressures of the horizontal WHAMs in Segment1, Segment2, Segment3.

Theoretically, the larger the value of the infinitesimal number m, the more accurate the model will be, but the computational time will be significantly longer. With the increase of m, the computed result changed significantly and then tended to be stable. When m is greater than 30, the change of computed result becomes very small. Considering both the computational time and the accuracy of the model, m is set as 100.

3.2. Forward kinematics case study

This section presents five groups of calculation results as a case study for comparison, as shown in Table II. The first three groups correspond to the results caused by proportional changes of input pressure under the same load force. The last three groups correspond to the results of different load forces under the same input pressure.

Table II.

Input pressure and corresponding terminal position and attitude.

Figure 5.

Forward kinematic solution process.

3.3. Inverse kinematic solution

The inverse kinematic solution process is to solve input pressures of the manipulator according to the given terminal position and attitude and load force. Also, the input pressure of the variable-diameter grippers should be given as a known condition. The solution process of inverse kinematics can be expressed as

(14) \begin{equation}{}^n\gamma {\rm{ = }}\left[ {\begin{array}{*{20}{c}}x\\ \\[-7pt]y\\ \\[-7pt]z\\ \\[-7pt]\delta \\ \\[-7pt]\beta \\ \\[-7pt]\psi\end{array}} \right];F_N = \left[ {F_N} \right];{}^HP = \left[ {\begin{array}{c@{\quad}c}{{p_{h1}} - {p_e}}\\ \\[-7pt]{{p_{h2}} - {p_e}}\\ \\[-7pt]{{p_{h3}} - {p_e}}\end{array}} \right] \Rightarrow {}^VP{\rm{ = }}\left[ {\begin{array}{c@{\quad}c@{\quad}c}{{p_{11}} - {p_e}}&{{p_{12}} - {p_e}}&{{p_{13}} - {p_e}}\\ \\[-7pt]{{p_{21}} - {p_e}}&{{p_{22}} - {p_e}}&{{p_{23}} - {p_e}}\\ \\[-7pt]{{p_{31}} - {p_e}}&{{p_{32}} - {p_e}}&{{p_{33}} - {p_e}}\end{array}} \right]\end{equation}

Due to the complexity of calculation and many variables and matrix operations, it is very complex and difficult to solve inverse kinematics by inverse derivation. A DNN is designed to solve the inverse kinematics in this paper. In the DNN, the terminal position and attitude $x,y,z,\delta,\beta,\psi$ and the input pressures of the grippers $^HP$ are input, and the input pressures of the vertical muscles $^VP$ in the three segments are output. The TensorFlow library of Python is used to build the DNN. The DNN has 10 input nodes: the terminal position and attitude, load force, and input pressures of the grippers in the three segments $x,y,z,\delta,\beta,\psi, F_N, p_{h1},p_{h2},p_{h3}$ . To unify the dimensions, training data are normalized before input. The output node has nine nodes: the input pressures of the vertical muscles in the three segments $p_{11}, p_{12}, p_{13}, p_{21}, p_{22}, p_{23}, p_{31}, p_{32}, p_{33}$ . The DNN can be expressed as

(15) \begin{equation}^VP{\rm{ = DNN\_Net(}}^3\gamma {\rm{,}}{^HP}{\rm{)}}\end{equation}

The network has 6 hidden layers, each with 200 nodes. ADAM optimization algorithm is selected as the training method. RELU activation function is selected as the activation function of the hidden layers and linear activation function is selected in the output layer. To avoid overfitting and ensure generalization ability, Batch Normalization (BN) layers are constructed and the net is trained in batches. Mean Square Error (MSE) is selected as the loss function. Ninety-thousand forward kinematic calculation results in the working space are randomly chosen as the training set. At training epoch 12,000, the MSE dropped to $7e^{-4}$ . Another 10,000 forward kinematic results are input to the model as the testing set to test the network’s performance, and the correlation coefficient between the output and target samples reaches 0.945.

4.1. Characteristics of working space

According to the calculation process as shown in Fig. 5, the terminal position and attitude of the underwater manipulator with three segments studied in this paper are calculated. By loop assignment of input pressure and load force, the terminal coordinates are plotted, and the working space of the manipulator is figured out. As shown in Fig. 6, the working place of the manipulator is a bowl-shaped entity with a symmetrical center and an undulating edge. Point A, B, and C are convex points obtained with the maximum input pressure on one side and the minimum input pressure on the other two sides. Point D, E, and F are concave points obtained with the maximum input pressure on two sides and the minimum input pressure on the other side. Point G is the lowest point of the depression on the surface and is obtained when all input pressures are at the maximum value.

4.2. Effect of load force on working space

Figure 7 shows the effect of the load force $F_N$ on the working space under the same maximum of input pressure ${p_{\rm{max}}}$ . The load forces are selected according to the actual application of the manipulator. For easy understanding, cross-section figures are added to the right side of the three-dimensional figures. With the increase of the load force, the working space of the manipulator becomes smaller. The paper compute the volume of the workspaces and compared the effects of impact factors in proportion, as shown in Table III.

Table III.

Effect of load force on working space.

Figure 6.

Working space with special points.

Figure 7.

Effect of the load force on working space.

4.3. Effect of maximum input pressure and water depth on working space

For the input pressure matrix of the underwater manipulator

\[^IP=\left[ {\begin{array}{c@{\quad}c@{\quad}c}{{p_{11}} - {p_e}}&{{p_{12}} - {p_e}}&{{p_{13}} - {p_e}}\\{{p_{21}} - {p_e}}&{{p_{22}} - {p_e}}&{{p_{23}} - {p_e}}\\{{p_{31}} - {p_e}}&{{p_{32}} - {p_e}}&{{p_{33}} - {p_e}}\\{{p_{h1}} - {p_e}}&{{p_{h2}} - {p_e}}&{{p_{h3}} - {p_e}}\end{array}} \right]\]

When ${p_{\rm{max}}}$ increases, the working space of the underwater manipulator also increases. At the same time, as the water depth increases, the external pressure ${p_e}$ increases, which is equivalent to reducing the effective input pressure of the WHAMs. Each increase in a water depth of 10 m is approximately equal to a reduction in the maximum input pressure of 0.1 MPa. Figure 8 and Table IV show the effect of ${p_{\rm{max}}}$ on the working space under the same load force.

Table IV.

Effect of maximum input pressure on working space.

Figure 8.

Effect of the maximum input pressure on working space.

4.4. Maximum load force and the rate of no solution points

As shown in Fig. 9(a), with the increase of $p_{\rm{max}}$ , the maximum load force $F_{N{\rm{max}}}$ that the manipulator can bear also increases. Under the maximum input pressure of 1.0 MPa, $F_{N{\rm{max}}}$ reaches 1598 N. With the change of the given load force and the maximum input pressure, the rate of no solution points in all input data sets also changes. As shown in Fig. 9(b), the greater the maximum input pressure, the less the no solution points under the same load force.

Figure 9.

(a) Relationship between maximum load force and input pressure; (b) Effects of load force and input pressure on no solution points.

Figure 10.

(a) Experimental system setups; (b) Software interface; (c) Experimental subsystems.

Figure 11.

Position and attitude of the manipulator under different maximum input pressures.

Figure 12.

Comparison of experimental results and computational results.