2.1. Design of soft actuators

The main body of the modular bionic soft robot in Fig. 1 is made of silicone rubber. The whole robot consists of three driving modules, each of which is divided into four parts: two soft actuators juxtaposed above, a top elastomer layer, a bottom elastomer layer and the strain limiting layer located in the middle [Reference Wu, Jiangbei and Yanqiong30]. One soft actuator has 11 air chambers and each of them is connected by a tunnel located at the bottom of the air chambers (Fig. 1). Each driving module has two air pipes, each of which is inserted into the middle air chamber of the actuator. During inflation, the air passing through the inflatable pipe reaches the middle air chamber of the soft actuator and then reaches each air chamber through the tunnel. With the increase of air chamber pressure, each air chamber expands to drive the driving module to bend. As shown in Fig. 2(b), the partition walls of chambers bulge to contact each other after inflation. The manufacturing parameters of the air chamber are shown in Fig. 2(a) and (c) and their values are shown in Table I. They are determined by model and experience. The establishment of the model is described in Section 2.2.

Figure 1.

The spatial view of the modular bionic soft robot.

Figure 2.

Sectional view of air chambers. (a) Before inflating view and (b) after inflating view. (c) Top sectional view.

Table I.

Manufacturing parameters of soft actuators.

The nonlinear deformation process of each soft actuator can be divided into four stages: (1) driving module 1 inflates. In the initial process, since module 2 and module 3 of the modular bionic soft robot are flat on the ground, the friction between the front and the rear of module 1 and the ground is different, so the rear of the modular bionic soft robot is rolled forward after module 1 is inflated. (2) Driving module 1 deflates while driving module 2 inflates. During the deflation of driving module 1, due to the deformation and bending of the modular bionic soft robot after inflation, the rear of the modular bionic soft robot is in inclined contact with the ground, resulting in large rear end friction, which makes the modular bionic soft robot extend forward during the deflation process. When driving module 2 is inflated at the same time, the rear end of driving module 2 will stretch driving module 1 forward. (3) Driving module 3 inflates while driving module 2 deflates. The process is similar to stage 2. (4) Driving module 3 deflates. Since driving module 1 and driving module 2 of the modular bionic soft robot are flat on the ground, the friction between the front and the rear of driving module 3 and the ground is different, so the front of the modular bionic soft robot extends forward when driving module 3 deflates.

During the locomotion of the modular bionic soft robot, the bottom friction of the modular bionic soft robot is different when different modules inflate and deform. As shown in Fig. 3, the distance of the robot moving forward in a modular bionic soft is S through the cyclic inflating and deflating behavior of module 1, module 2 and module 3.

Figure 3.

Locomotion process of the modular bionic soft robot in one period.

2.2. Analysis of nonlinear deformation characteristics

The nonlinear deformation analysis of the soft actuator is carried out to determine the variation of displacement h of the expansion wall of the air chamber with the inflation pressure P. As shown in Fig. 2(a) and (c), the cavity in the exhausted state is rectangular with a gap k. When inflated, the partition walls of the air chamber expand and contact each other, as shown in Fig. 2(b). We assume that the partition wall expands into a cylindrical surface (i.e., an arc on the sectional view) and contacts each other in a tangent manner. Compared with the contact wall, the deformation of the thicker top and side walls is much smaller, so it is assumed that the chamber height b is constant during expansion. In addition, when the expansion is extruded by an external force, the top and side walls will protrude. Through experiments, we have observed that the convexity of the top and side walls is much smaller than that of the contact walls. In this study, in order to simplify the model, only the deformation of the contact wall is considered, and the deformation is h. Assuming that the silicone rubber material is incompressible (I ₃ = 1) and isotropic, the deformation of silicone rubber can be analyzed by strain energy density function [Reference Jianlong, Guangjuan and Zhengwei31].

(1) \begin{align}W = W\left( {{I_1},{I_2},{I_3}} \right)\end{align}

where W is the strain energy density, and I ₁ , I ₂ and I ₃ are the three invariants of the deformation tensor.

(2) \begin{align}\left\{ \begin{array}{l}{I_1} = \lambda _1^2 + \lambda _2^2 + \lambda _3^2\\[5pt]{I_2} = \lambda _1^{ - 2} + \lambda _2^{ - 2} + \lambda _3^{ - 2}\\[5pt]{I_3} = \lambda _1^2\lambda _2^2\lambda _3^2 = 1\end{array} \right.\end{align}

Among Eq. (2), λ _i is the main elongation ratio in three directions.

According to the relationship between stress and strain of silicone rubber material, the relationship between principal stress and principal elongation ratio of silicone rubber material can be obtained:

(3) \begin{align}{\sigma _1} = \frac{2}{{{\lambda _1}}}\left({\lambda _1}^2 - \frac{1}{{\lambda _1^2\lambda _2^2}}\right)\left(\frac{{\partial W}}{{\partial {I_1}}} + \lambda _2^2\frac{{\partial W}}{{\partial {I_2}}}\right)\end{align}

(4) \begin{align}{\sigma _2} = \frac{2}{{{\lambda _2}}}\left({\lambda _2}^2 - \frac{1}{{\lambda _1^2\lambda _2^2}}\right)\left(\frac{{\partial W}}{{\partial {I_1}}} + \lambda _1^2\frac{{\partial W}}{{\partial {I_2}}}\right)\end{align}

Among Eqs. (3) and (4), ${\sigma _1},{\sigma _2}$ are the principal stresses in both directions of the tensile plane.

With Yeoh model, the binomial parameter form [Reference Yeoh32] of the strain energy function is

(5) \begin{align}W = {C_1}({I_1} - 3) + {C_2}{({I_2} - 3)^2}\end{align}

Among Eq. (6), ${C_1},{C_2}$ are material coefficients, which are determined by tensile tests (BYES, 5000 N, 700 mm, China). According to the standard ASTMD412 and then curve fitting, which results in C ₁ = 0.112 and C ₂ = 0.019. o is located at the center of a single air chamber as the coordinate origin, and Y direction is the expansion direction of the air chamber wall. In the process of inflation, any point on the air chamber wall moves in the Y direction, while the X and Z coordinates remain unchanged (Fig. 2(b)). The displacement of any point on the left air chamber wall is U, as shown in Eq. (8). By integrating the left air chamber wall, the corresponding strain energy function is obtained, as shown in Eqs. (4) and (5).

(6) \begin{align}E = {\int {Wdv = \int_{{{c \over 2}}}^{{{c \over 2}} + d} {\int_{ - {{b \over 2}}}^{{{b \over 2}}} {\int_{ - {{a \over 2}}}^{{{a \over 2}}} {{C_1}({I_1} - 3) + {C_2}({I_2} - 3)} } } } ^2}dxdzdy\end{align}

(7) \begin{align}{I_1} = tr\left[\left(I + \frac{{\partial (U)}}{{\partial (x,y,z)}}\right){\left(I + \frac{{\partial (U)}}{{\partial (x,y,z)}}\right)^T}\right]\end{align}

(8) \begin{align}U = \left({U_x},{U_y},{U_z}\right) = \left(0,hf\left(\frac{{2z}}{a}\right)f\left(\frac{{2x}}{b}\right)\right)\end{align}

where ν is the volume of a single left air chamber wall, I is a unit matrix and f(τ) is called the profile function. a,b,c and d are manufacturing parameters of the air chamber (Fig. 2). According to Eqs. (6)–(8), strain energy can be obtained.

(9) \begin{align}E = 4{C_2}\left[{A_0}{A_1}\left(\frac{{{a^2}}}{{{b^2}}} + \frac{{{b^2}}}{{{a^2}}}\right) + 2{C^2}\right]\frac{{{h^4}d}}{{ab}} + 4{C_1}{B_0}{B_1}\left(\frac{a}{b} + \frac{b}{a}\right){h^2}d\end{align}

Among Eq. (9), ${A_0},{A_1},{B_0},{B_1},C$ are constants, which are determined by function $f(\tau )$ :

(10) \begin{align}\left\{ \begin{array}{l}{A_0} = \int_0^1 {{{[f(\tau )]}^4}d\tau } \\[4pt]{A_1} = \int_0^1 [f'(\tau )]^4 d\tau \\[4pt]{B_0} = \int_0^1 [f(\tau )]^2 d\tau \\[4pt]{B_1} = \int_0^1 [f'(\tau )]^2 d\tau \\[4pt]C = \int_0^1 [f(\tau )]^2 [f^{\prime}(\tau )]^2 d\tau \end{array} \right.\end{align}

Generally, f(τ) should comply with the actual deformation of the left air chamber wall as following conditions, $f(0) = 1\,,\,f( - 1) = f(1) = 0\,,\,f'(0) = 0$ .

To simplify the function, we construct the function $f(\tau ) = 1 - {\tau ^2}$ . Then, the constants are calculated according to f(τ).

By using the principle of virtual work, the expansion pressure P is connected with the maximum displacement h, and Eq. (11) is obtained. The left side is the virtual change of the total strain energy of the left air chamber wall with the expansion height, and the right side is the virtual work of the left air chamber wall under the expansion pressure.

(11) \begin{align}\frac{{dE(h)}}{{dh}}\delta h = P\frac{{dV(h)}}{{dh}}\delta h\end{align}

(12) \begin{align}V(h) = hab\int_0^1 {{{[f(\tau )]}^2}d\tau } \end{align}

Among Eq. (12), V(h) is the increased volume inside the air chamber after inflating. Among Eqs. (9)–(12), the relationship between the pressure P inside the air chamber and the maximum displacement h of the air chamber expansion wall can be obtained.

(13) \begin{align}P = \left[0.2340\left(\frac{{{a^2}}}{{{b^2}}} + \frac{{{b^2}}}{{{a^2}}}\right) + 0.0334\right]\frac{{{h^3}d}}{{{a^2}{b^2}}}\end{align}

According to Eq. (10), the pressure P is proportional to the cubic of the deformation h of the air chamber expansion wall. With the increase of the pressure in the air chamber, the expansion speed of the air chamber slows down. The manufacturing parameters of the air chamber are shown in Table I. Parameters a, b and d are obtained from Eq. (13). Other parameters have little influence on the deformation of the air chamber during inflation. Therefore, other parameters are determined by the manufacturing requirements.

When the pressure P increases, the deformation of the air chamber increases. The bending angle of the driving module will also increase. However, in the inflation experiment, when the pressure P exceeds 90 kPa, the structure of the driving module will be damaged. Considering the characteristics of silicone rubber and the safety of the driving module, the air pressure of the driving module is set at 70 kPa.

3.1. Type A: module 2 is the main motion module

When module 2 is the main motion module, the modular bionic soft robot model is shown in Fig. 4, then the overall motion height of the modular bionic soft robot is $H = {h_0} + {h_1} + {h_2}$ .

Figure 4.

Module 2 is the main motion module, (a) the model of the modular bionic soft robot and (b) the actual modular bionic soft robot.

It is easy to get h ₀ from the manufacturing parameters of a single driving module. The bending angle $\phi$ of module 2 can be obtained in Eq. (14). Equation (14) shows that the bending angle $\phi$ is related to the inflation pressure P of module 2.

(14) \begin{align}\varphi = \dfrac{{2(N - 1)\left\{ \sqrt[3]{{\dfrac{{P{a^2}{b^2}}}{{d\left[0.234\left(\dfrac{{{a^2}}}{{{b^2}}} + \dfrac{{{b^2}}}{{{a^2}}}\right) + 0.0334\right]}}}} - \dfrac{k}{2}\right\} }}{{\dfrac{a}{2} + c}}\end{align}

As shown in Fig. 4(a), L is the original length of the main motion module, R is the radius of curvature when the structure is bent and $\phi$ is the central angle after deformation.

(15) \begin{align}L = (f + 2d)N + (N - 1)k + 2n = R\varphi \end{align}

(16) \begin{align}{h_1} = R - R\cos \frac{\varphi }{2}\end{align}

After obtaining h ₁ , the force analysis of the auxiliary motion module 1 is carried out, as shown in Fig. 4(b). f is the frictional force on the auxiliary motion module 1. F ₀ is the ground support force of the auxiliary motion module 1. F ₂ is the force exerted by the main motion module on the auxiliary motion module 1.

(17) \begin{align}{F_2}\sin \theta + mg = {F_0}\end{align}

(18) \begin{align}{F_2}\cos \theta = \mu {F_0}\end{align}

After the lifting angle θ of auxiliary motion module 1 is obtained, it is necessary to determine the lifting height when auxiliary motion module 1 is not inflated. At the beginning of the movement, the inflation pressure of module 2 increases, the bending angle $\phi$ of module 2 increases with the increase of inflation pressure P, and the lifting angle θ of auxiliary motion module 1 increases from 0°. As shown in Fig. 5(a), the lifting height h of auxiliary motion module 1 is a piecewise function f(θ) which varies with the lifting angle of auxiliary motion module 1, and the expression of f(θ) is as follows:

(19) \begin{align}f{(\theta )_1} = k\sin \frac{\theta }{2} + (2d + f + n)\sin \theta ,\,0 \le \theta \le \pi - 2{\cos ^{ - 1}}\frac{k}{{2a}}\end{align}

(20) \begin{align}f{\left(\theta \right)_2}\, =\, & k\left[\sin \left(\theta - \frac{\pi }{2} + {\cos ^{ - 1}}\frac{k}{{2a}}\right) + \sin \left(\frac{\theta }{2} - \frac{\pi }{2} + {\cos ^{ - 1}}\frac{k}{{2a}}\right)\right] + \left(2d + f + n\right)\sin \theta \nonumber\\[3pt]& + \left(2d + f\right)\sin \left(\theta - \pi + 2{\cos ^{ - 1}}\frac{k}{{2a}}\right),\quad \pi - 2{\cos ^{ - 1}}\frac{k}{{2a}} \le \theta \le 2\left(\pi - 2{\cos ^{ - 1}}\frac{k}{{2a}}\right) \end{align}

(21) \begin{align}f{\left(\theta \right)_3} \,=\, & k\left[\sin \left(\theta - \frac{\pi }{2} + {\cos ^{ - 1}}\frac{k}{{2a}}\right) + \sin \left(\theta - \frac{{3\pi }}{2} + 3{\cos ^{ - 1}}\frac{k}{{2a}}\right) + \sin \left(\frac{\theta }{2} - \pi + 2{\cos ^{ - 1}}\frac{k}{{2a}}\right)\right]\nonumber\\[4pt] & + \left(2d + f + n\right)\sin \theta + \left(2d + f\right)\left[\sin \left(\theta - \pi + 2{\cos ^{ - 1}}\frac{k}{{2a}}\right) \right.\nonumber\\[4pt] & \left. + \sin \left(\theta - 2\pi + 4{\cos ^{ - 1}}\frac{k}{{2a}}\right)\right] ,2\left(\pi - 2{\cos ^{ - 1}}\frac{k}{{2a}}\right) \le \theta \le 3\left(\pi - 2{\cos ^{ - 1}}\frac{k}{{2a}}\right)\end{align}

(22) \begin{align}f{\left(\theta \right)_4} =\, & k\left[\sin \left(\theta - \frac{\pi }{2} + {\cos ^{ - 1}}\frac{k}{{2a}}\right) + \sin \left(\theta - \frac{{3\pi }}{2} + 3{\cos ^{ - 1}}\frac{k}{{2a}}\right) + \sin \left(\theta - \frac{{5\pi }}{2} + 5{\cos ^{ - 1}}\frac{k}{{2a}}\right)\right.\nonumber\\[5pt]& + \left. \sin \left(\frac{\theta }{2} - \frac{{3\pi }}{2} + 3{\cos ^{ - 1}}\frac{k}{{2a}}\right)\right] + \left(2d + f + n\right)\sin \theta + \left(2d + f\right)\left[\sin \left(\theta - \pi + 2{\cos ^{ - 1}}\frac{k}{{2a}}\right)\right.\nonumber\\[5pt]& \left. + \sin \left(\theta - 2\pi + 4{\cos ^{ - 1}}\frac{k}{{2a}}\right) + \sin \left(\theta - 3\pi + 6{\cos ^{ - 1}}\frac{k}{{2a}}\right)\right] ,3\left(\pi - 2{\cos ^{ - 1}}\frac{k}{{2a}}\right) \le \theta\nonumber\\&\quad \le 4\left(\pi - 2{\cos ^{ - 1}}\frac{k}{{2a}}\right)\end{align}

Figure 5.

(a) When the modular bionic soft robot moves, the displacement of the auxiliary motion module 1 changes with time, (b) when module 2 is the main motion module, the force analysis diagram of module 1 during inflation, (c) when module 2 is the main motion module, the force analysis diagram of auxiliary motion module 1.

Through Eqs. (19)–(22), the deformed height of the actuator in other states can be deduced. After getting the elevation height of auxiliary motion module 1 in the initial state, the change of h ₂ is analyzed when its stiffness changes (i.e., driving module 1 inflates). When the expansion wall displacement $h \le k/2$ , driving module 1 bends downward. When the expansion wall displacement $h \gt k/2$ , driving module 1 expands upward.

1. (1) Expansion wall displacement $h \le k/2$ .

   When the inflation pressure of module 2 is 70 kPa, module 1 is lifted by six chambers in the initial state, and the geometric analysis is carried out, as shown in Fig. 5(c).

   (23) \begin{align}\left\{ \begin{array}{l}{L_0} = 7(2d + f) + 12h = (R + m + 12h)\alpha \\[5pt]{L_1} = 7(2d + f) + 6k = (R + c + m + a)\alpha \end{array} \right.\end{align}

   It can be solved:

   (24) \begin{align}\alpha = \frac{{6(k - 2h)}}{{c + a - 12h}}\end{align}

   So we can get the height h2,

   (25) \begin{align}{h_2} = l\sin \frac{\alpha }{2}\end{align}

2. (2) Expansion wall displacement $h \gt k/2$ .

   The inflation pressure of module 2 is 70 kPa. Module 1 begins to expand and bend, and the geometric analysis is carried out, as shown in Fig. 5 (b).

   (26) \begin{align}\left\{ \begin{array}{l}{L_1} = 11(2d + f) + 10k = R\alpha \\[5pt]{L_0} = 11(2d + f) + 20h = (R + c + \dfrac{a}{2})\alpha \end{array} \right.\end{align}

   It can be solved:

   (27) \begin{align}\alpha = \frac{{10(2h - k)}}{{c + \dfrac{a}{2}}}\end{align}

   It can be obtained from geometric relations.

   (28) \begin{align}\gamma = \frac{\alpha }{2} - \beta \end{align}
   (29) \begin{align}{h_2} = l\sin \gamma \end{align}

Therefore, by combining h, the relationship curves between the overall height H of the modular bionic soft robot and the inflation pressure P of the module 1 and module 3 can be obtained.

The curves in Fig. 6 show that when the inflation pressure P of module 1 and module 3 is less than 20.07 kPa, the stiffness of module 1 and module 3 increases, the bending of module 1 and module 3 decreases and the overall height H of module decreases. When the inflation pressure of module 1 and module 3 is greater than 20.07 kPa, the stiffness of module 1 and module 3 continues to increase, resulting in obvious expansion deformation, which makes the overall height H of the module increase rapidly. Therefore, when the inflation pressure of module 1 and module 3 equals 20.07 kPa, the overall height of the modular bionic soft robot is the lowest, and the height of the narrow terrain through which it can pass is the lowest. At this time, the stiffness of module 1 and module 3 is the best.

Figure 6.

Relationship between inflation pressure P of the modular bionic soft robot and the whole height H of the modular bionic soft robot.

By comparing the experimental data, as shown in Fig. 6, the solid line is the theoretical curve and the dotted line is the experimental acquisition curve. By comparing the experimental data with the theoretical data, the error is found, and the error is relatively stable. The main reason is that the thickness of the bottom plate of the modular bionic soft robot is slightly thicker than that of the design, which makes the overall height of the actual module slightly higher than that of the theoretical data.

As shown in Fig. 7, with the inflation pressure of module 2 unchanged, the inflation pressure of module 1 and module 3 increases, which changes the stiffness of the module. The overall height of the modular bionic soft robot first decreases and then increases. When the inflation pressure is 20 kPa, the overall height of the modular bionic soft robot is the lowest, and the height of the narrow terrain through which it could pass is the lowest. This height is 42.8 mm.

Figure 7.

Schematic diagram of height change of the modular bionic soft robot.

3.2. Type B. Module 1 or module 3 is the main motion module.

At this time, the model of the modular bionic soft robot is shown in Fig. 8, and the overall height of the modular bionic soft robot is $H = {h_0} + {h_1}$ .

Figure 8.

When module 1 or module 3 is the main motion module, (a) the model of the modular bionic soft robot and (b) the actual modular bionic soft robot.

Module 1 needs to be analyzed when calculating h ₁ , as shown in Fig. 8(a). Assuming that the center of the circle corresponding to module 1 is O, and the distance between the center of the circle and the highest point of module 1 is R, the vertical distance between the point O and the ground can calculate the overall height of the three modules.

(30) \begin{align}\theta = {\sin ^{ - 1}}\left[\frac{{{h_2}}}{{2R\sin \left(\frac{\varphi }{2}\right)}}\right]\end{align}

(31) \begin{align}\alpha = \frac{\pi }{2} - \frac{\varphi }{2} - \theta \end{align}

(32) \begin{align}{y_0} = R\sin \varphi \end{align}

(33) \begin{align}{h_1} = R - {y_0}\end{align}

Comparing the data of type B with the data of type A, it can be found that under the same conditions, the overall height of the modular bionic soft robot is higher when module 2 is the main motion module, which is the highest height that the modular bionic soft robot can pass through obstacles.

After adjusting the stiffness of different modules, it can be observed in the experiment that the motion height of the modular bionic soft robot decreases, the vibration amplitude decreases and the motion speed increases. The optimum stiffness of the auxiliary motion module will also change under different air pressure of the main motion module. The optimum air pressure of the deflation module can be calculated by Eqs. (25) and (29).

When the manufacturing size of the module remains unchanged, the minimum height H of the modular bionic soft robot varies with the change of the pressure of the main motion module P, as shown in Fig. 9(a). The relationship between module dimensions a, b, d and inflation pressure P of deflation module is shown in Fig. 9(b).

Figure 9.

(a) The relationship between the main motion module P and the minimum height H of the modular bionic soft robot (b) the relationship between the module size a, b, d and the inflation pressure P of the auxiliary motion module.