Application and Analysis of an Ionic Liquid Gel in a Soft Robot

Abstract

Due to their light weight, flexibility, and low energy consumption, ionic electroactive polymers have become a hotspot for bionic soft robotics and are ideal materials for the preparation of soft actuators. Because the traditional ionic electroactive polymers, such as ionic polymer-metal composites (IPMCs), contain water ions, a soft actuator does not work properly upon the evaporation of water ions. An ionic liquid polymer gel is a new type of ionic electroactive polymer that does not contain water ions, and ionic liquids are more thermally and electrochemically stable than water. These liquids, with a low melting point and a high ionic conductivity, can be used in ionic electroactive polymer soft actuators. An ionic liquid gel (ILG), a new type of soft actuator material, was obtained by mixing 1-butyl-3-methylimidazolium tetrafluoroborate (BMIMBF₄), hydroxyethyl methacrylate (HEMA), diethoxyacetophenone (DEAP) and ZrO₂ and then polymerizing this mixture into a gel state under ultraviolet (UV) light irradiation. An ILG soft actuator was designed, the material preparation principle was expounded, and the design method of the soft robot mechanism was discussed. Based on nonlinear finite element theory, the deformation mechanism of the ILG actuator was deeply analyzed and the deformation of the soft robot when grabbing an object was also analyzed. A soft robot was designed with the soft actuator as the basic module. The experimental results show that the ILG soft robot has good driving performance, and the soft robot can grab a 105 mg object at an input voltage of 3.5 V.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Traditional robots are typically constructed of rigid motion joints based on hard materials and can perform tasks quickly and accurately. However, such robots have limited movement flexibility and low environmental adaptability and cannot work under space-constrained conditions. These factors limit the application of rigid robots in unstructured and complex environments [1–3].

In nature, mollusks are widely distributed in oceans, rivers, and lakes and on land. After hundreds of millions of years of evolution, this animal has gradually developed the characteristics of a large deformation ability, a light weight, and a high power density ratio and can achieve efficient movement under complex natural environment conditions by changing its body shape. In recent years, researchers have attempted to apply the biological principles of mollusks to the research and design of robotics. A soft robot is composed of a soft material that can withstand large strains, has the capabilities for continuous deformation and drive-structure integration and can arbitrarily change its shape and size over a wide range. Such a robot has strong adaptability to unstructured environments and broad application prospects in military reconnaissance, rescue, and medical surgery [4–6].

A flexible actuator based on an ionic electroactive polymer (EAP) has the advantages of a low driving voltage, a large displacement ability, a light weight, etc. Such actuators have become a research hotspot in the field of bionic robots. In the past few decades, electrochemical actuators, which are substitutes for air- and fluid-derived devices, have been further developed due to their ideal mechanical properties for use in intelligent robots. However, traditional ionic EAP actuators, such as those based on IPMCs, have the characteristics of a short working time in nonwater media, a complex manufacturing process, and a high cost. An ILG is a new type of ionic EAP that can work in air for a long time. Because an ILG offers chemical stability, thermal stability, and simpler ion transport, it is more suitable for the production of soft robot actuators than other EAPs.

In this paper, a soft robot actuator based on an ILG material is proposed and its preparation, driving mechanism, and design method are deeply analyzed. Finally, the validity and rationality of the soft robot are verified by driving performance and grabbing experiments.

2. Design of the ILG Soft Robot

The ILG soft actuator consists of a 5-layer structure, as shown in Figure 1. The middle layer is an electroactive layer composed of the ILG material, which functions to store ionic liquids. The outer two layers that wrap the middle layer are electrode layers, which consist of activated carbon. Activated carbon has a high specific surface area, a high electrical conductivity, a high density, a strong adsorption force, etc., making it very suitable for use as an electrode material [7]. In addition, gold foil is selected as a current collector to cover the surface of the activated carbon layers. When the actuator is working, one end is attached to an external metal electrode, and a wire is led from the external electrode to connect to a power supply.

[figure omitted; refer to PDF]

The soft robot consists of three ILG actuators, a common copper electrode, and three independent copper electrodes, as shown in Figure 2. The input unit consists of a common electrode and an independent electrode. The common electrode is a copper column, with the upper end fixed with a fixture, and the copper electrode is connected to a power source. The soft robot consists of three ILG actuators evenly distributed around the common electrode at 120°. Two copper electrodes sandwich one ILG actuator.

[figure omitted; refer to PDF]

Figure 3 shows the motion cycle of the soft robot to grab an object, which can be described as follows. First, the ILG actuator is naturally stretched, and by applying an input voltage of 3.5 V to the electrodes, the actuator quickly bends outward. Second, when the actuator is close to the target object, the input voltage direction is changed, and the actuator bends inward to clamp the object. Finally, the input signal is changed, and the actuator bends outward to release the object.

[figure omitted; refer to PDF]

2.1. Preparation of the Soft Actuator

The preparation process of an ILG is a process of ionic liquid loading. Ionic liquid loading refers to the immobilization of an ionic liquid to form a solid carrier by physical or chemical means [8]. BMIMBF₄ was selected as the ionic liquid. The carriers of ionic liquids need to have interconnected porous structures, large specific surface area and porosity, and good mechanical strength and electrochemical stability. Therefore, HEMA was selected to prepare the carrier. This material is a colorless and transparent liquid that is soluble in water and has low toxicity. HEMA was polymerized under ultraviolet light to form polyhydroxyethyl methacrylate (PHEMA) with a porous structure.

The mixed solution consists of BMIMBF₄, HEMA, and ZrO₂. Under UV light, the ILG is generated by the polymerization reaction induced by DEAP. The polymer matrix is cross-linked into a porous network structure. ZrO₂ enhances the mechanical strength of the gel, and as the ZrO₂ content increases, the tensile strength of the ILG increases [9, 10]. Figure 4 shows the soft actuator fabrication process. Figure 4(a) shows the ILG sample, which is covered with a layer of BMIMBF₄ due to osmosis. Because the ionic liquid has a certain viscosity, the activated carbon electrode can be directly affixed to both sides of the gel, as shown in Figure 4(b). Finally, the gold foil is evenly affixed to the surface of the activated carbon layers to obtain an ILG soft actuator, as shown in Figure 4(c).

[figures omitted; refer to PDF]

2.2. Soft Actuator Deformation Mechanism

For the ILG actuator, the application of a voltage on the electrode can cause bending deformation, which is the result of the coupling of the electric field, chemical field, flow field, and field force. The average pore diameter of the structure is greater than 10 microns, which is much larger than the diameter of free ions in the ionic liquids (nanoscale). Therefore, BMIM⁺ and BF₄⁻ can move freely inside the carrier. Under the action of the electric field, positive and negative ions accumulate in the electrode layers on both sides of the ILG. Because the cation volume is much larger than the anion volume, the volume of the ILG cathode side expands and the volume of the ILG anode side shrinks [11, 12]. The entire actuator therefore bends toward the anode side, as shown in Figure 5.

[figure omitted; refer to PDF]

When an electric field is applied, the movement of the internal ions of the ILG actuator in the porous medium can be described by the Nernst–Planck equation [13]:$$\begin{array}{}\text{(1)}& j_i=-d_i\left(\nabla {\varpi }_i+\frac{F}{RT}{\varpi }_i{z}_i\nabla \Theta \right),\end{array}$$where i is the ion type (i = 1 represents the cation, and i = 2 represents the anion), j_i is the flow of ions in the PHEMA pores, d_i is the diffusion coefficient of the ions, ${\varpi }_i$ is the ion concentration, $\nabla {\varpi }_i$ is the ion concentration gradient, $\Theta$ is the potential, $\nabla \Theta$ is the potential gradient, F is the Faraday constant, R is the gas constant, T is the absolute temperature, and z_i is the number of ionic charges. The first term on the right-hand side of the formula indicates the contribution of the cation or anion concentration diffusion; the second term indicates the contribution of electromigration.

The continuous equation of concentration changing with time is$$\begin{array}{}\text{(2)}& \frac{\partial {\varpi }_i}{\partial t}=-\nabla j_i,\end{array}$$where $\nabla j_i$ is the flow change of the $i‐\text{th}$ ion in the PHEMA pores and t is the time.

The internal electric field balance equation of the ILG actuator can be expressed as follows:$$\begin{array}{c}\text{(3)}& E=\frac{D}{k_{\text{e}}}=-\nabla \Theta ,\nabla D=\rho =F\left({\varpi }_1-{\varpi }_2\right),\end{array}$$where D is the electric displacement, E is the electric field, ρ is the net charge density, and k_e is the effective dielectric constant of the ILG.

The effect of the electrode layer on the bending deformation of the actuator along the width and the thickness is ignored. Considering only the relationship between the axial deformation of the electrode layer and the net charge density, the axial induced strain can be expressed as follows:$$\begin{array}{}\text{(4)}& {\epsilon }_{\text{c}}\left(s\right)=\rho \left(s\right)\alpha \left(s\right),\end{array}$$where ${\epsilon }_{\text{c}}$ is the induced strain, $\alpha \left(s\right)$ is the induced strain coupling coefficient, and s is the Laplace domain.

Therefore, the axial induced stress is$$\begin{array}{c}\text{(5)}& {\sigma }_{\text{c}}\left(s\right)={\epsilon }_{\text{c}}{E}_{\text{t}}=\rho \left(s\right)d\left(s\right),d\left(s\right)=\alpha \left(s\right)E_{\text{t}},\end{array}$$where E_t is the elastic modulus of the electrode layer and $d\left(s\right)$ is the axial-induced stress coupling coefficient.

The bending moment of the actuator is$$\begin{array}{}\text{(6)}& M\left(s\right)=2{\int }_{b_1}^b{\sigma }_{\text{c}}\left(s\right)xwdx=2{\int }_{b_1}^b\rho \left(s\right)d\left(s\right)xwdx,\end{array}$$where b and b₁ are half the thickness of the actuator and ILG layer, respectively, and $w$ is the width of the actuator.

The ILG actuator is the basic unit of a soft robot and is the source of its deformation. For flexible mechanisms, the degrees of freedom and constraints are determined by the flexibility matrix. The cantilever beam flexibility matrix is as follows [14]:$$\begin{array}{c}\text{(7)}& R=\text{diag}\hspace{0.17em}\hspace{0.17em}\left(R_{\theta x},R_{\theta y},R_{\theta z},R_x,R_y,R_z\right)=\text{diag}\hspace{0.17em}\hspace{0.17em}\left(\frac{l}{E_{\text{a}}{I}_x},\frac{l}{E_{\text{a}}{I}_y},\frac{l}{G_{\text{a}}J},\frac{l^3}{12E_{\text{a}}{I}_y},\frac{l^3}{12E_{\text{a}}{I}_x},\frac{l}{E_{\text{a}}A}\right),\end{array}$$where R_θx, R_θy, and R_θz represent rotational flexibility around the x, y, and z axes, respectively, and R_x, R_y, and R_z represent translational flexibility along the x, y,and z axes, respectively. If the order of magnitude of an element in one direction is significantly greater than that of the element in the other directions, then the mechanism has a degree of freedom in this direction. Using this method, the equivalent constraints or degrees of freedom of the soft robot can be established.

3. Hyperelastic Arruda–Boyce Model

Material nonlinearity is caused by the nonlinearity of the constitutive relationship of the material. The constitutive relationship of a material includes the stress, strain, strain rate, load duration, temperature, etc. After a load is unloaded, residual strain will appear in a nonlinear material, as shown in Figure 6.

[figure omitted; refer to PDF]

The hyperelastic Arruda–Boyce model, typically used for relatively largely deformed materials, is used to simulate ILGs. Hyperelastic materials can be described by the strain energy function W, which is used to define the strain energy density in the material. The Arruda–Boyce strain energy function W is as follows [15–17]:$$\begin{array}{}\text{(8)}& W=G\sum _{k=1}^5\frac{c_k}{{\lambda }_{\text{m}}^{2\left(k-1\right)}}\left(I_1^k-3^k\right)+\frac{1}{D}\left(\frac{J_{\text{el}}^2-1}{2}-\ln \hspace{0.17em}\hspace{0.17em}{J}_{\text{el}}\right),\end{array}$$with $c_1=\left(1/2\right)$, $c_2=\left(1/20\right)$, $c_3=\left(11/1050\right)$, $c_4=\left(19/7000\right)$, and $c_5=\left(519/673750\right)$, where W is the strain energy density, G is the initial shear modulus, λ_m is the deformation rate, J_el is the elastic volume ratio, and D is the temperature-dependent material parameter associated with the bulk modulus. A material is completely incompressible when J_el = 1.

I ₁ is the first strain invariant of the deviatoric strain, and the relationship between I₁ and the main tensile strain rates λ₁, λ₂, and λ₃ is as follows:$$\begin{array}{c}\text{(9)}& I_1={\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2,I_2=\frac{1}{{\lambda }_1^2}+\frac{1}{{\lambda }_2^2}+\frac{1}{{\lambda }_3^2},\\ I_3={\lambda }_1^2\times {\lambda }_2^2\times {\lambda }_3^2,\end{array}$$where λ₁, λ₂, and λ₃ are the main (extension) deformation rates along the x, y, and z directions, respectively. I₁, I₂, and I₃ are the relative changes of the length, the surface area, and the volume of the elastomer, respectively.

The material stress is obtained by differentiating the stretching:$$\begin{array}{}\text{(10)}& {\sigma }_{\text{n}}=\frac{\partial W}{\partial {\lambda }_{\text{n}}},\end{array}$$where ${\sigma }_{\text{n}}$ is the normal stress and λ_n is the stretching in the loading direction.

The isotropic and incompressible deformation process of an ILG is given by$$\begin{array}{}\text{(11)}& \sqrt{I_3}={\lambda }_1{\lambda }_2{\lambda }_3=1.\end{array}$$

When the ILG is uniformly stretched, ${\lambda }_2={\lambda }_3$; equation (11) can be used to calculate$$\begin{array}{}\text{(12)}& {\lambda }_2={\lambda }_3=\frac{1}{\sqrt{{\lambda }_1}}.\end{array}$$

The true stress is obtained by$$\begin{array}{}\text{(13)}& {\sigma }_0={\sigma }_{\text{n}}{\lambda }_{\text{n}}.\end{array}$$

Differentiating equation (8), we obtain$$\begin{array}{}\text{(14)}& {\sigma }_0=G\left[\left({\lambda }_{\text{n}}^2-{\lambda }_{\text{n}}^{-1}\right)+\frac{1}{5{\lambda }_{\text{m}}^2}\left({\lambda }_{\text{n}}^4-{\lambda }_{\text{n}}^{-2}\right)+\frac{11}{175{\lambda }_{\text{m}}^4}\left({\lambda }_{\text{n}}^6-{\lambda }_{\text{n}}^{-3}\right)+\frac{19}{875{\lambda }_{\text{m}}^6}\left({\lambda }_{\text{n}}^8-{\lambda }_{\text{n}}^{-4}\right)+\frac{519}{67375{\lambda }_{\text{m}}^8}\left({\lambda }_{\text{n}}^{10}-{\lambda }_{\text{n}}^{-5}\right)\right].\end{array}$$

4. Model for the Soft Robot

To establish a gripping model scheme as shown in Figure 7, assume that the soft robot is an inextensible, fully flexible elastic rod. Assume that l is the length of the evenly distributed rod, EI is the flexural rigidity, ρ is the mass per unit length, and k_g is the curvature profile when the object is gripped [18, 19].

[figure omitted; refer to PDF]

4.1. Deformation Analysis for Soft Actuators

At $s=s^{\ast }$, the position of any material point can be represented by the following formula [20], as shown in Figure 8.$$\begin{array}{}\text{(15)}& \stackrel{\rightharpoonup }{r}\left(s=s^{\ast }\right)=X\left(s=s^{\ast }\right){\stackrel{\rightharpoonup }{E}}_1+Y\left(s=s^{\ast }\right){\stackrel{\rightharpoonup }{E}}_2,\end{array}$$where, in Cartesian coordinates, X and Y are represented as follows:$$\begin{array}{}\text{(16)}& \left\{\begin{array}{c}X\left(s=s^{\ast }\right)=X\left(s=0\right)+{\int }_0^{s^{\ast }}\cos \left(\theta \left(\eta \right)\right)d\eta ,\hfill \\ Y\left(s=s^{\ast }\right)=Y\left(s=0\right)+{\int }_0^{s^{\ast }}\sin \left(\theta \left(\eta \right)\right)d\eta ,\hfill \end{array}\right.\end{array}$$where η is a dummy variable and θ is the angle between the unit tangent vector and the horizontal direction.$$\begin{array}{}\text{(17)}& \stackrel{\rightharpoonup }{r^{\prime }}=\cos \left(\theta \left(s\right)\right){\stackrel{\rightharpoonup }{E}}_1+\sin \left(\theta \left(s\right)\right){\stackrel{\rightharpoonup }{E}}_2,\end{array}$$where the prime (′) represents the partial derivative with respect to s. Assuming that $\stackrel{\rightharpoonup }{r}$ is continuous, θ and $\stackrel{\rightharpoonup }{r^{\prime }}$ are also continuous.

[figure omitted; refer to PDF]

At point $s=\varsigma$, the jump in any function $\chi =\chi \left(s,\theta \left(s\right),{\theta }^{\prime }\left(s\right)\right)$ is expressed by a compact notation:$$\begin{array}{}\text{(18)}& {⟦\chi ⟧}_{\varsigma }=\chi \left({\varsigma }^{+}\right)-\chi \left({\varsigma }^{-}\right),\end{array}$$where$$\begin{array}{}\text{(19)}& \left\{\begin{array}{c}\chi \left({\varsigma }^{-}\right)=\underset{s\nearrow \varsigma }{\lim }\chi \left(s,\theta \left(s\right),{\theta }^{\prime }\left(s\right)\right),\hfill \\ \chi \left({\varsigma }^{+}\right)=\underset{s\searrow \varsigma }{\lim }\chi \left(s,\theta \left(s\right),{\theta }^{\prime }\left(s\right)\right).\hfill \end{array}\right.\end{array}$$

The jump will be associated with the force applied at the discrete point.

The bending moment M is described by a constitutive equation:$$\begin{array}{}\text{(20)}& \stackrel{\rightharpoonup }{M}=EI\left({\theta }^{\prime }-k_0\right){\stackrel{\rightharpoonup }{E}}_3,\end{array}$$where k₀ represents the curvature. In addition to the gravity per unit length and end load of the rod, the singular force $\stackrel{\rightharpoonup }{F_{\xi }}$ at $s=\xi$ must be considered. Figure 7 shows the force that simulates the contact of an object with the soft robot. The governing equation can be derived from the equilibrium of the linear and angular momentum [21–23]:$$\begin{array}{}\text{(21)}& \left\{\begin{array}{c}\stackrel{\rightharpoonup }{m\prime }-\rho g\stackrel{\rightharpoonup }{E_2}=0,\hfill \\ \frac{\partial }{\partial s}\left(EI\left({\theta }^{\prime }-k_0\right)+m_2\cos \left(\theta \right)-m_1\sin \left(\theta \right)\right)=0,\hfill \\ {⟦\stackrel{\rightharpoonup }{m}⟧}_{\xi }+\stackrel{\rightharpoonup }{F_{\xi }}=0,\hfill \end{array}\right.\end{array}$$where $\stackrel{\rightharpoonup }{m}=m_1{\stackrel{\rightharpoonup }{E}}_1+m_2{\stackrel{\rightharpoonup }{E}}_2$.

4.2. Deformation Analysis Governing Formula for the Soft Robot

The corner of the object is assumed to be in contact with the soft robot at s = ξ. Two new unit vectors are defined as follows, as shown in Figure 9:$$\begin{array}{}\text{(22)}& \left\{\begin{array}{c}{\stackrel{\rightharpoonup }{H}}_1=\cos \left(\theta \left(\xi \right)\right){\stackrel{\rightharpoonup }{E}}_1+\sin \left(\theta \left(\xi \right)\right){\stackrel{\rightharpoonup }{E}}_2,\\ {\stackrel{\rightharpoonup }{H}}_1=-\sin \left(\theta \left(\xi \right)\right){\stackrel{\rightharpoonup }{E}}_1+\cos \left(\theta \left(\xi \right)\right){\stackrel{\rightharpoonup }{E}}_2,\end{array}\right.\end{array}$$where ${\stackrel{\rightharpoonup }{H}}_1$ is the tangent to the centerline of the soft robot and ${\stackrel{\rightharpoonup }{H}}_2$ is perpendicular to the centerline of the soft robot.$$\begin{array}{}\text{(23)}& {\stackrel{\rightharpoonup }{F}}_{\xi }=f_1{\stackrel{\rightharpoonup }{H}}_1+f_2{\stackrel{\rightharpoonup }{H}}_2=N_1{\stackrel{\rightharpoonup }{E}}_1+N_2{\stackrel{\rightharpoonup }{E}}_2,\end{array}$$where f₁ is the normal force and f₂ is the corresponding friction force.

[figure omitted; refer to PDF]

Equation (21) is applied to segment $s\in \left(\xi ,l\right)$, and we obtain the following equation:$$\begin{array}{}\text{(24)}& \left\{\begin{array}{c}\stackrel{\rightharpoonup }{m}\left({\xi }^{+}\right)=-\rho g\left(l-\xi \right){\stackrel{\rightharpoonup }{E}}_2,\\ \stackrel{\rightharpoonup }{m}\left({\xi }^{-}\right)=N_1{\stackrel{\rightharpoonup }{E}}_1+N_2{\stackrel{\rightharpoonup }{E}}_2-\rho g\left(l-\xi \right){\stackrel{\rightharpoonup }{E}}_2.\end{array}\right.\end{array}$$

The expression of the potential energy W of the soft robot is established:$$\begin{array}{}\text{(25)}& W={\int }_0^{\xi }\left\{\frac{EI}{2}{\left({\theta }^{\prime }-k_0\right)}^2+\rho gY\left(s\right)-\stackrel{\rightharpoonup }{m}\left({\xi }^{-}\right)\cdot \stackrel{\rightharpoonup }{r^{\prime }}\right\}ds+{\int }_{\xi }^l\left\{\frac{EI}{2}{\left({\theta }^{\prime }-k_0\right)}^2+\rho gY\left(s\right)\right\}ds,\end{array}$$where$$\begin{array}{}\text{(26)}& Y\left(s\right)={\int }_0^s\sin \left(\theta \left(\xi \right)\right)d\xi .\end{array}$$

Equation (21) can be used to establish boundary value problems.$$\begin{array}{}\text{(27)}& \left\{\begin{array}{cc}EI\left({\theta }^{″}-k_0^{\prime }\right)-\rho g\left(l-s\right)\cos \left(\theta \right)=0,& s\in \left(0,\xi \right),\\ EI\left({\theta }^{″}-k_0^{\prime }\right)-\rho g\left(l-s\right)\cos \left(\theta \right)+N_2\cos \left(\theta \right)-N_1\sin \left(\theta \right)=0,& s\in \left(\xi ,l\right),\end{array}\right.\end{array}$$where the solution $\theta ={\theta }^{\ast }$ satisfies the conditions$$\begin{array}{c}\text{(28)}& \theta \left(0\right)=0,{\theta }^{\prime }\left(l\right)=k_0\left(l\right),\\ {⟦\theta ⟧}_{\xi }=0,\\ {⟦{\theta }^{\prime }⟧}_{\xi }=0,\\ {\int }_0^{\xi }\cos \left(\theta \left(s\right)\right)ds=\frac{c}{2}.\end{array}$$

The soft robot is divided into K segments. The total potential energy expression is approximated as follows:$$\begin{array}{}\text{(29)}& W_{\text{d}}=\frac{EI}{2}\sum _{i=1}^K{\left(\frac{{\theta }_i-{\theta }_{i-1}}{ds}-k_0\right)}^{2}ds+\rho g\left(l-s_0\right)\sin \left({\theta }_0\right)\frac{ds}{2}+\sum _{i=1}^{K-1}\rho g\left(l-s_i\right)\sin \left({\theta }_i\right)ds+\rho g\left(l-s_K\right)\sin \left({\theta }_n\right)\frac{ds}{2}-N_1\cos \left({\theta }_0\right)\frac{ds}{2}-N_1\sum _{i=1}^{I-1}\cos \left({\theta }_i\right)ds-N_1\cos \left({\theta }_I\right)\frac{ds}{2}-N_2\sin \left({\theta }_0\right)\frac{ds}{2}-N_2\sum _{i=1}^{I-1}\sin \left({\theta }_i\right)ds-N_2\sin \left({\theta }_I\right)\frac{ds}{2},\end{array}$$where the constraint function W_c is$$\begin{array}{}\text{(30)}& W_{\text{c}}=\frac{ds}{2}+\sum _{i=1}^{I-1}\cos \left({\theta }_i\right)ds+\cos \left({\theta }_I\right)\frac{ds}{2}-\frac{c}{2}.\end{array}$$

5. Nonlinear Finite Element Method

Geometric nonlinearity arises from the nonlinear relationship between the strain and displacement. Currently, research on geometric nonlinearity mainly focuses on three types of problems: (1) a large displacement with a small strain, (2) a small displacement with a large strain, and (3) a large displacement with a large strain. A geometric nonlinear problem has two main characteristics. First, due to the large deformation of the structure, the strain and displacement of the structure are nonlinear. Second, a balance equation is established at the position after deformation. In the analysis of large deformation, the displacement of the structure changes continuously and appropriate strain, stress, and constitutive relationships should be adopted.

To capture the nonlinear behavior of the structure, the full nonlinear finite element formulas of the truss elements and beam elements are studied.

5.1. Nonlinear FEM for Truss Elements

In the Cartesian coordinate system, the object is displaced to a certain position under the action of external forces, as shown in Figure 10. (x₁, y₁, z₁) and (x₂, y₂, z₂) are the position coordinates of points P₁ and P₂ before deformation, respectively. The object is deformed to a new position under the action of external forces. ($u_1$, $v_1$, $w_1$) and ($u_2$, $v_2$, $w_2$) are the deformation coordinates of points P₁ and P₂ after deformation, respectively.

[figure omitted; refer to PDF]

The formula for the undeformed length of a truss element is as follows:$$\begin{array}{}\text{(31)}& l=\sqrt{{\left(x_2-x_1\right)}^2+{\left(y_2-y_1\right)}^2+{\left(z_2-z_1\right)}^2}.\end{array}$$

The expression of the total Lagrangian strain ε along the deformation axis is$$\begin{array}{}\text{(32)}& \epsilon =\frac{\sqrt{\left(u_2+x_2-u_1-x_1\right)+\left(v_2+y_2-v_1-y_1\right)+\left(w_2+z_2-w_1-z_1\right)}-l}{l}.\end{array}$$

According to the change of the elastic energy, the product of the stiffness matrix and the displacement vector is expressed as follows [2, 24]:$$\begin{array}{}\text{(33)}& \left\{\begin{array}{c}\left[k\right]\left\{u\right\}=AlE\epsilon \left\{\phi \right\},\\ \left\{\phi \right\}={\left\{\frac{\partial \epsilon }{\partial u_1},\frac{\partial \epsilon }{\partial v_1},\frac{\partial \epsilon }{\partial w_1},\frac{\partial \epsilon }{\partial u_2},\frac{\partial \epsilon }{\partial v_2},\frac{\partial \epsilon }{\partial w_2}\right\}}^T,\end{array}\right.\end{array}$$where E is Young’s modulus, A is the cross sectional area, and l is the length of the truss element.

The tangent stiffness matrix is obtained by differentiating equation (33) with respect to the displacement vector:$$\begin{array}{c}\text{(34)}& \left[\stackrel{⌣}{k}\right]=\frac{\partial \left(\left[k\right]\left\{u\right\}\right)}{\partial \left\{u\right\}}=EAl\left\{\phi \right\}\frac{\partial \epsilon }{\partial \left\{u\right\}}+EAl\epsilon \frac{\partial \left\{\phi \right\}}{\partial \left\{u\right\}}=EAl\left\{\phi \right\}{\left\{\phi \right\}}^T+EAl\epsilon \frac{{\partial }^2\left\{\phi \right\}}{\partial {\left\{u\right\}}^2}.\end{array}$$

The mass matrix is the same as that for a linear truss element [25, 26]:$$\begin{array}{}\text{(35)}& \left[m\right]=\frac{\rho Al}{6}\left[\begin{array}{cccccc}2& 0& 0& 1& 0& 0\\ 0& 2& 0& 0& 1& 0\\ 0& 0& 2& 0& 0& 1\\ 1& 0& 0& 2& 0& 0\\ 0& 1& 0& 0& 2& 0\\ 0& 0& 1& 0& 0& 2\end{array}\right].\end{array}$$

5.2. Nonlinear FEM for Beam Elements

The two famous beam theories are the Euler–Bernoulli and Timoshenko beam theories. The Euler–Bernoulli beam theory assumes that the cross section remains planar and normal to the reference line after bending and its stiffness is higher than the actual stiffness. Timoshenko’s beam theory overcomes this problem by introducing shear deformation into the model, which obtains accurate results for thick beam calculations. Because the soft actuators are slender, the shear deformation is negligible. Therefore, nonlinear Euler–Bernoulli beams with von Karman nonlinearity can be used for modeling and analysis of the soft actuators.

In addition to the bending effect, the finite element formula should also include the torsion and tensile effects to reflect the large deformation effect [27, 28]:$$\begin{array}{}\text{(36)}& \left\{\begin{array}{c}{\epsilon }_{11}=e+z{\rho }_2-y{\rho }_3,\\ {\epsilon }_{12}=-z{\rho }_1,\\ {\epsilon }_{13}=y{\rho }_1,\\ e=u^{\prime }+\frac{{\left(v^{\prime }\right)}^2}{2}+\frac{{\left(w^{\prime }\right)}^2}{2},\\ {\rho }_1={\phi }^{\prime },\\ {\rho }_2=-w^{″}\left(1-{\left(w^{\prime }\right)}^2\right),\\ {\rho }_3=-v\left(1-{\left(v\right)}^2\right),\end{array}\right.\end{array}$$where ε_ij is the engineering strain tensor, e is the axial strain, y and z are the coordinates on the cross section, and ρ_i is the deformation curvature. $u$, $v$, and $w$ are the displacements on the cross section.

When a change of energy is applied, the product of the stiffness matrix and the displacement vector is expressed as follows:$$\begin{array}{}\text{(37)}& \left[k\right]\left\{h\right\}={\int }_l{\left[B\right]}^T{\left[\Psi \right]}^T\left[\Phi \right]\left\{\varphi \right\}dx,\end{array}$$where$$\begin{array}{c}\text{(38)}& \left[B\right]=\frac{\partial \left[N\right]}{\partial s},\left[\Psi \right]=\left[\begin{array}{cccccc}1& v^{\prime }& w^{\prime }& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 2w^{\prime }{w}^{″}& 0& 0& -\left(1-{\left(w^{\prime }\right)}^2\right)\\ 0& -2v^{\prime }{v}^{″}& 0& 0& 1-{\left(v^{\prime }\right)}^2& 0\end{array}\right],\\ \left[\Phi \right]=\left[\begin{array}{cccc}EA& 0& 0& 0\\ 0& G{I}_{11}& 0& 0\\ 0& 0& E{I}_{22}& 0\\ 0& 0& 0& E{I}_{33}\end{array}\right],\\ \left\{\varphi \right\}=\left\{\begin{array}{c}e\\ {\rho }_1\\ {\rho }_2\\ {\rho }_3\end{array}\right\}.\end{array}$$

Assume that$$\begin{array}{}\text{(39)}& \left\{h^{\left(i\right)}\right\}=\left\{\overline{h}\right\}+\left\{\Delta h^{\left(i\right)}\right\}.\end{array}$$

By applying Taylor expansion without higher order terms, the tangential stiffness matrix is expressed as follows:$$\begin{array}{}\text{(40)}& \left[\stackrel{⌣}{k}\right]={\int }_l{\left[B\right]}^T\left({\left[\overline{\Psi }\right]}^T\left[\Phi \right]\left\{\Psi \right\}+\left[\Gamma \right]\right)\left[B\right]dx,\end{array}$$where$$\begin{array}{}\text{(41)}& \left[\Gamma \right]=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& EAe-2E{I}_{33}{v}^{″}{\rho }_3& 0& 0& -2E{I}_{33}{v}^{\prime }{\rho }_3& 0\\ 0& 0& EAe+2E{I}_{22}{w}^{″}{\rho }_2& 0& 0& 2E{I}_{33}{w}^{\prime }{\rho }_2\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right].\end{array}$$

The mass matrix is the same as that of a linear beam element.

5.3. Nonlinear Convergence Algorithms

Although the Newton–Raphson algorithm is used for nonlinear static analysis, the tangent stiffness matrix becomes singular at some points, which makes a global equilibrium solution impossible. Riks proposed a process to track the intersection of the normal to the tangent line with the equilibrium path to solve this problem [29], as shown in Figure 11. The Crisfield method uses arcs instead of vertical lines to search for solutions. The increment of the load factor becomes an unknown problem to be solved in the iterative process.

[figure omitted; refer to PDF]

Assume that$$\begin{array}{}\text{(42)}& \left\{F\right\}=\lambda \left\{\stackrel{\rightharpoonup }{F}\right\},\end{array}$$where λ is the load parameter and $\stackrel{\rightharpoonup }{F}$ is the preselected load vector.

The residual force vector of the $i‐\text{th}$ step is expressed as follows:$$\begin{array}{}\text{(43)}& R\left({\left\{\stackrel{\rightharpoonup }{U}\right\}}^i,{\lambda }^i\right)={\left.R\left({\left\{\stackrel{\rightharpoonup }{U}\right\}}^{i-1},{\lambda }^{i-1}\right)+\frac{\partial R}{\partial \left\{\stackrel{\rightharpoonup }{U}\right\}}\right|}_{{\left\{\stackrel{\rightharpoonup }{U}\right\}}^{i-1},{\lambda }^{i-1}}\delta {\left\{\stackrel{\rightharpoonup }{U}\right\}}^i+{\left.\frac{\partial R}{\partial \left\{\lambda \right\}}\right|}_{{\left\{\stackrel{\rightharpoonup }{U}\right\}}^{i-1},{\lambda }^{i-1}}\delta {\lambda }_1^i,\end{array}$$where $\stackrel{\rightharpoonup }{U}$ is the displacement vector.$$\begin{array}{}\text{(44)}& \left[k\right]{\left.\left\{\stackrel{\rightharpoonup }{U}\right\}\right|}_{{\left\{\stackrel{\rightharpoonup }{U}\right\}}^{i-1},{\lambda }^{i-1}}-{\lambda }^{i-1}\left\{\stackrel{\rightharpoonup }{F}\right\}+{\left.\left[k_T\right]\right|}_{{\left\{\stackrel{\rightharpoonup }{U}\right\}}^{i-1},{\lambda }^{i-1}}\delta {\left\{\stackrel{\rightharpoonup }{U}\right\}}_1^i+\left\{\stackrel{\rightharpoonup }{F}\right\}\delta {\lambda }_1^i=\left\{0\right\},\end{array}$$where $\delta {\left\{\stackrel{\rightharpoonup }{U}\right\}}_1^i$ is the first increment of the load parameter in the $i‐\text{th}$ step and $\delta {\lambda }_1^i$ is the first increment of the displacement vector in the $i‐\text{th}$ step.$$\begin{array}{}\text{(45)}& {\left\{\stackrel{\rightharpoonup }{U}\right\}}_1^i={\left\{\stackrel{\rightharpoonup }{U}\right\}}^{i-1}+\delta {\left\{\stackrel{\rightharpoonup }{U}\right\}}_1^i,\text{(46)}& {\lambda }_1^i={\lambda }^{i-1}+\delta {\lambda }_1^i.\end{array}$$

Substituting (45) and (46) into equation (44), we obtain$$\begin{array}{}\text{(47)}& \left[k\right]{\left.\left\{\stackrel{\rightharpoonup }{U}\right\}\right|}_{{\left\{\stackrel{\rightharpoonup }{U}\right\}}_1^i,{\lambda }_1^i}-{\lambda }_1^i\left\{\stackrel{\rightharpoonup }{F}\right\}+{\left.\left[k_T\right]\right|}_{{\left\{\stackrel{\rightharpoonup }{U}\right\}}_1^i,{\lambda }_1^i}\delta {\left\{\stackrel{\rightharpoonup }{U}\right\}}_2^i+\left\{\stackrel{\rightharpoonup }{F}\right\}\delta {\lambda }_2^i=\left\{0\right\}.\end{array}$$

When the second step incremental search path is perpendicular to the normal of the first step incremental path, we obtain$$\begin{array}{}\text{(48)}& \delta {\left\{\stackrel{\rightharpoonup }{U}\right\}}_1^i\cdot \delta {\left\{\stackrel{\rightharpoonup }{U}\right\}}_2^i+\left\{\stackrel{\rightharpoonup }{F}\right\}\delta {\lambda }_1^i\cdot \left\{\stackrel{\rightharpoonup }{F}\right\}\delta {\lambda }_2^i=\left\{0\right\}.\end{array}$$

When solving nonlinear finite element equations by an iterative method, choosing the appropriate convergence criterion is necessary to ensure that the iteration can be terminated. The convergence criterion will directly affect the speed and accuracy of the solution. If the convergence criterion is not appropriately chosen, the calculation will fail.

From equations (47) and (48), we obtain$$\begin{array}{c}\text{(49)}& \delta {\left\{\stackrel{\rightharpoonup }{U}\right\}}_2^i=-{\left({\left.\left[k_T\right]\right|}_{{\left\{\stackrel{\rightharpoonup }{U}\right\}}_1^i,{\lambda }_1^i}\right)}^{-1}\left(\left[k\right]{\left.\left\{\stackrel{\rightharpoonup }{U}\right\}\right|}_{{\left\{\stackrel{\rightharpoonup }{U}\right\}}_1^i,{\lambda }_1^i}-{\lambda }_1^i\left\{\stackrel{\rightharpoonup }{F}\right\}+\left\{\stackrel{\rightharpoonup }{F}\right\}\delta {\lambda }_2^i\right),\delta {\lambda }_2^i=\frac{\delta {\left\{\stackrel{\rightharpoonup }{U}\right\}}_1^i{\left({\left.\left[k_T\right]\right|}_{{\left\{\stackrel{\rightharpoonup }{U}\right\}}_1^i,{\lambda }_1^i}\right)}^{-1}\left(\left[k\right]{\left.\left\{\stackrel{\rightharpoonup }{U}\right\}\right|}_{{\left\{\stackrel{\rightharpoonup }{U}\right\}}_1^i,{\lambda }_1^i}-{\lambda }_1^i\left\{\stackrel{\rightharpoonup }{F}\right\}\right)}{\left(\left\{\stackrel{\rightharpoonup }{F}\right\}\delta {\lambda }_1^i\left\{\stackrel{\rightharpoonup }{F}\right\}\right)-\delta {\stackrel{\rightharpoonup }{U}}_1^i{\left({\left.\left[k_T\right]\right|}_{{\left\{\stackrel{\rightharpoonup }{U}\right\}}_1^i,{\lambda }_1^i}\right)}^{-1}\left\{\stackrel{\rightharpoonup }{F}\right\}}.\end{array}$$

The convergent load parameters and displacement vectors are as follows:$$\begin{array}{c}\text{(50)}& {\lambda }^i={\lambda }^{i-1}+\delta {\lambda }_1^i+\delta {\lambda }_2^i+\delta {\lambda }_3^i+\cdots ,{\left\{\stackrel{\rightharpoonup }{U}\right\}}^i={\left\{\stackrel{\rightharpoonup }{U}\right\}}^{i-1}+\delta {\left\{\stackrel{\rightharpoonup }{U}\right\}}_1^i+\delta {\left\{\stackrel{\rightharpoonup }{U}\right\}}_2^i+\delta {\left\{\stackrel{\rightharpoonup }{U}\right\}}_3^i+\cdots .\end{array}$$

6. Experimental Analysis

6.1. Stress Calculation of the ILG

To verify the accuracy and practicability of the induced stress expression of the ILG, the deformation rate of each point in the uniaxial tensile test is substituted into the stress expression to calculate the stress.

For the Arruda–Boyce model, the parameters G and λ_m are as follows:$$\begin{array}{c}\text{(51)}& G=4.10736,{\lambda }_{\text{m}}=\mathrm{4.62529.}\end{array}$$

A comparison of the uniaxial tensile stress calculation data with the experimental data is given in Table 1.

Table 1

Comparison of the tensile stress calculation data with experimental data obtained in the uniaxial tensile test.

Measurement data		Experimental results		Arruda–Boyce model
Force (mN)	Deformation (mm)	Strain (%)	Stress (kPa)	Calculated stress (σE) (kPa)	Error (%)
4.97	2.0	5.7	0.50	0.54	7.42
8.21	3.5	10	0.86	0.92	6.64
14.88	7.0	20	1.70	1.78	4.70
20.52	10.5	30	2.54	2.64	4.08
25.35	14.0	40	3.38	3.53	4.38
29.47	17.5	50	4.21	4.43	5.21
33.06	21.0	60	5.04	5.34	6.02
36.22	24.5	70	5.86	6.26	6.69
38.94	28.0	80	6.68	7.16	7.31
41.48	31.5	90	7.50	8.06	7.40
43.58	35.0	100	8.30	8.94	7.70
45.50	38.5	110	9.10	9.80	7.65
47.13	42.0	120	9.88	10.63	7.64
48.67	45.5	130	10.65	11.44	7.26
49.98	49.0	140	11.41	12.21	6.93
51.03	52.5	150	12.15	12.96	6.71
51.95	56.0	160	12.86	13.69	6.39
52.75	59.5	170	13.56	14.38	6.02
53.40	63.0	180	14.24	15.05	5.66
53.90	66.5	190	14.89	15.69	5.36
54.32	70.0	200	15.52	16.30	5.01
54.65	73.5	210	16.12	16.88	4.65
54.78	77.0	220	16.70	17.44	4.50
54.86	80.5	230	17.25	17.98	4.31
54.94	84.0	240	17.78	18.50	3.97
54.84	87.5	250	18.28	18.99	3.88
54.72	91.0	260	18.74	19.46	3.72
54.47	94.5	270	19.18	19.91	3.73
54.01	98.0	280	19.57	20.34	4.06
53.61	101.5	290	19.93	20.75	4.22
53.23	105.0	300	20.24	21.15	4.29
52.48	108.5	310	20.50	21.53	5.05
51.83	112.0	320	20.73	21.89	5.59
51.02	115.5	330	20.91	22.23	6.42
50.29	119.0	340	21.04	22.57	7.09
49.21	122.5	350	21.12	22.88	8.51

Table 1, which presents comparisons of the tensile stress calculation data with experimental data, shows that the minimum relative error of the Arruda–Boyce model is 3.72 and the maximum relative error is 8.51. This result shows that the calculation accuracy of the Arruda–Boyce model is high and satisfies the requirement of material performance analysis in soft robot design.

Figure 12 shows that the first half of the stress-strain curve calculated by the Arruda–Boyce model almost coincides with the experimental curve, while for the second half, the relative error between the calculated and experimental stress-strain curves becomes increasingly larger but remains small. The above analyses indicate that the calculation formula of the ILG stress is feasible and reliable.

[figure omitted; refer to PDF]

6.2. Experiment on the Soft Actuator

The square wave of the input signal has a magnitude of 4 V, a period of 30 s, and a duty cycle of 50%. The actuator size is 30 mm × 5 mm × 0.5 mm, the free segment length is 25 mm, and the average thickness of the ILG layer is approximately 0.4 mm.

The interval of the displacement reading is 3 s. Figure 13 shows that the maximum amplitude of the positive axis of the actuator is 4.8 mm, and the maximum amplitude of the negative axis is 5 mm. The positive and negative displacements of the actuator are basically symmetrical and stable, and the attenuation of the long-term working drive performance is very small.

[figure omitted; refer to PDF]

In Figure 14, the simulation results are analyzed: when the input voltage is 3.0 V, the maximum displacement of the soft actuator is 3.7 mm; when the input voltage is 3.5 V, the maximum displacement is 4.3 mm; and when the input voltage is 4.0 V, the maximum displacement is 5.0 mm. The displacement calculation results for the different input voltages are shown in Table 2.

[figures omitted; refer to PDF]

Table 2

Displacement calculation results for the different input voltages.

Force (V)	3.0	3.5	4.0
Deformation (mm)	3.7	4.3	5.0

According to Figure 15, the simulation value curve is almost consistent with the experimental value curve, and the change trends of the displacement curves from the simulation and experiment are almost identical. The relative error between the calculated and experimental displacement curves is small. Therefore, the calculation formula of the soft actuator is feasible and credible.

[figure omitted; refer to PDF]

6.3. Working Performance Test of the Soft Actuator

The long-term stable working performance is an important factor in ionic EAP soft actuators. Since the traditional IPMC actuator contains a water-based electrolyte, maintaining a stable driving performance after operation in air for dozens of seconds is difficult due to the evaporation of water and electrolysis [30]. The ILG soft robot does not contain a water-based electrolyte, which avoids the problem of water molecule loss, and can work stably in air for a long time.

The current response stability of the soft actuator is tested by applying a square wave signal with a magnitude of 4 V, a period of 30 s, and a duty cycle of 50%, and the current response curve is shown in Figure 16. The experimental results show that the current peaks remain unchanged over 30 minutes, and the current response of the ILG actuator is stable. Therefore, the internal electrolyte is not substantially decomposed.

[figure omitted; refer to PDF]

From the above experiments, the ILG actuator can maintain stable current response and bending displacement in air. These results show that the ILG actuator can effectively avoid the defects of traditional EAP actuators and has excellent driving performance.

6.4. Experiment of the Soft Robot Grabbing an Object

First, an input voltage of 3.5 V is applied across the electrodes, and the soft robot gradually opens outward. Next, the operating platform is moved close to the target object. Finally, by changing the direction of the supply voltage, the soft robot reversely bends to clamp the object.

The parameters of the soft robot are shown in Table 3. The values 1∼3 in the first column are the ILG actuator number, corresponding the three sets of execution components on the manipulator, with a total mass of 211 mg.

Table 3

Parameters of the soft manipulator.

Number	Total length (mm)	Mass (g)	Free segment length (mm)	Width (mm)	Total thickness (mm)	ILG layer thickness (mm)
1	35	0.068	30	5	0.5	0.4
2	35	0.070	30	5	0.5	0.4
3	35	0.073	30	5	0.5	0.4

A set of screenshots of the soft robot manipulator’s experimental process is shown in Figure 17. The process of the soft robot grabbing an object is shown in Figure 3. From Figures 17(a) and 17(b), when no object is grabbed, the deformation of the soft actuator is larger, and the deformation is similar to that in Figure 13. In this case, the mass of the object clamped by the soft robot at an input voltage of 3.5 V is approximately 105 mg. The experiments verify the effectiveness of ILG soft robots.

[figures omitted; refer to PDF]

Compared with the soft robot designed by Saito et al., which has a load capacity of approximately 3 mg, the ILG soft robot designed in this paper has a greatly increased load capacity. Since the soft robot consists of three actuators, this ILG soft robot has better adaptability for grabbing complex shape objects than Saito’s soft robot, which consists of two actuators [31].

In addition, the service life is an important performance indicator for EAP robots. The ILG material avoids the defects of IPMCs, and the soft robot proposed in this paper can theoretically work in air for a long time.

7. Conclusions

A soft robot is a new flexible structure, and the materials and mechanism design methods of a soft robot are quite different from those of traditional robots. To avoid the shortcomings of traditional rigid manipulators, such as large energy consumption and complex mechanisms, a soft robot based on a new ILG material is proposed in this paper.

A high-performance ILG actuator based on the principle of light curing was developed, and a modular design method based on motion and constraints was discussed. An ion transport model of porous media based on the Nernst–Planck equation was established, and the deformation mechanism of the ILG actuator was deeply analyzed in combination with the theory of cantilever beam deformation. Based on nonlinear finite element theory, the deformation of the soft robot when grabbing an object was analyzed. A three-finger type soft robot was designed, and the maximum displacement was 5 mm at a 4 V square wave voltage. The soft robot can grab a 105 mg object under a 3.5 V voltage. The results prove that the ILG bionic soft robot has good development prospects.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant nos. 51605334, 51705368, 61825303, and U1713215), Natural Science Foundation of Shanghai (Grant no. 17ZR1441800), Shanghai Sailing Program (Grant no. 17YF1420200), and International Exchange Program for Graduate Students, Tongji University (No. 2018020018).

References

[1] M. A. McEvoy, N. Correll, "Materials that couple sensing, actuation, computation, and communication," Science, vol. 347 no. 6228,DOI: 10.1126/science.1261689, 2015.

[2] B. He, C. Zhang, Y. Zhou, Z. Wang, "A computing method to determine the performance of an ionic liquid gel soft actuator," Applied Bionics and Biomechanics, vol. 2018,DOI: 10.1155/2018/8327867, 2018.

[3] B. He, Z. Wang, M. Li, K. Wang, R. Shen, S. Hu, "Wet adhesion inspired bionic climbing robot," IEEE/ASME Transactions on Mechatronics, vol. 19 no. 1, pp. 312-320, DOI: 10.1109/tmech.2012.2234473, 2014.

[4] D. Trivedi, C. D. Rahn, W. M. Kier, I. D. Walker, "Soft robotics: biological inspiration, state of the art, and future research," Applied Bionics and Biomechanics, vol. 5 no. 3, pp. 99-117, DOI: 10.1155/2008/520417, 2008.

[5] D. Rus, M. T. Tolley, "Design, fabrication and control of soft robots," Nature, vol. 521 no. 7553, pp. 467-475, DOI: 10.1038/nature14543, 2015.

[6] W. Zhao, A. Ming, M. Shimojo, "Development of high-performance soft robotic fish by numerical coupling analysis," Applied Bionics and Biomechanics, vol. 2018,DOI: 10.1155/2018/5697408, 2018.

[7] X. H. Liu, B. He, Z. P. Wang, "Tough nanocomposite ionogel-based actuator exhibits robust performance," Scientific Reports, vol. 3, 2014.

[8] S. J. Zhang, C. M. Xu, X. M. Lü, Ionic Liquid and Green Chemistry, 2009. in Chinese

[9] B. He, C.-H. Zhang, A. Ding, "Finite element analysis of ionic liquid gel soft actuator," Chinese Physics B, vol. 26 no. 12,DOI: 10.1088/1674-1056/26/12/126102, 2017.

[10] S. J. Chang, J. Chen, "Design and fabrication of the large thrust force piezoelectric actuator," Advances in Materials Science and Engineering, vol. 2013,DOI: 10.1155/2013/912587, 2013.

[11] Z. P. Wang, B. He, X. H. Liu, Q. G. Wang, "Development of soft manipulator based on ionogel," Chinese Science Bulletin, vol. 61 no. 23, pp. 2637-2647, DOI: 10.1360/n972015-01459, 2016.

[12] C. Zhang, B. He, A. Ding, S. Xu, Z. Wang, Y. Zhou, "Motion simulation of ionic liquid gel soft actuators based on CPG control," Computational Intelligence and Neuroscience, vol. 2019,DOI: 10.1155/2019/8256723, 2019.

[13] J. Newman, K. E. Thomas-Alyea, Electrochemical Systems, 2012.

[14] X. L. Ding, J. M. Selig, "Screw theoretic view on dynamics of spatially compliant beam (in Chinese)," Applied Mathematics and Mechanics, vol. 31 no. 9, pp. 1118-1132, DOI: 10.1007/s10483-010-1351-9, 2010.

[15] Y. C. Fung, Biomechanics: Mechanical Properties of Living Tissues, 1993.

[16] A. Avanzini, D. Battini, "Integrated experimental and numerical comparison of different approaches for planar biaxial testing of a hyperelastic material," Advances in Materials Science and Engineering, vol. 2016,DOI: 10.1155/2016/6014129, 2016.

[17] J. Bergstrom, M. C. Boyce, "Constitutive modeling of the large strain time-dependent behavior of elastomers," Journal of the Mechanics and Physics of Solids, vol. 46 no. 5, pp. 931-954, DOI: 10.1016/s0022-5096(97)00075-6, 1998.

[18] C. Majidi, R. F. Shepherd, R. K. Kramer, G. M. Whitesides, R. J. Wood, "Influence of surface traction on soft robot undulation," The International Journal of Robotics Research, vol. 32 no. 13, pp. 1577-1584, DOI: 10.1177/0278364913498432, 2013.

[19] C. Majidi, O. M. O’Reilly, J. A. Williams, "On the stability of a rod adhering to a rigid surface: shear-induced stable adhesion and the instability of peeling," Journal of the Mechanics and Physics of Solids, vol. 60 no. 5, pp. 827-843, DOI: 10.1016/j.jmps.2012.01.015, 2012.

[20] X. Zhou, C. Majidi, O. M. O’Reilly, "Soft hands: an analysis of some gripping mechanisms in soft robot design," International Journal of Solids and Structures, vol. 64-65, pp. 155-165, DOI: 10.1016/j.ijsolstr.2015.03.021, 2015.

[21] K. Suzumori, S. Endo, T. Kanda, N. Kato, "A bending pneumatic rubber actuator realizing soft-bodied manta swimming robot," Proceedings 2007 IEEE International Conference on Robotics and Automation, pp. 4975-4980, .

[22] N. Faruksenan, O. Oreilly, T. Tresierras, "Modeling the growth and branching of plants: a simple rod-based model," Journal of the Mechanics and Physics of Solids, vol. 56 no. 10, pp. 3021-3036, DOI: 10.1016/j.jmps.2008.06.005, 2008.

[23] X. Zhou, C. Majidi, O. M. O׳Reilly, "Flexing into motion: a locomotion mechanism for soft robots," International Journal of Non-Linear Mechanics, vol. 74,DOI: 10.1016/j.ijnonlinmec.2015.03.001, 2015.

[24] P. F. Pai, Highly Flexible Structures: Modeling, Computation, and Experimentation, 2007.

[25] K. J. Bathe, Finite Element Procedures, 1996.

[26] H. Shi, H. Salim, Y. Shi, F. Wei, "Geometric and material nonlinear static and dynamic analysis of space truss structures," Mechanics Based Design of Structures and Machines, vol. 43 no. 1, pp. 38-56, DOI: 10.1080/15397734.2014.925808, 2014.

[27] H. Shi, H. Salim, "Geometric nonlinear static and dynamic analysis of guyed towers using fully nonlinear element formulations," Engineering Structures, vol. 99, pp. 492-501, DOI: 10.1016/j.engstruct.2015.05.023, 2015.

[28] X. Guo, D. Liu, W. Zhang, L. Sun, S. Chen, "Nonlinear dynamic analysis of macrofiber composites laminated shells," Advances in Materials Science and Engineering, vol. 2017,DOI: 10.1155/2017/4073591, 2017.

[29] J. N. Reddy, An Introduction to Nonlinear Finite Element Analysis, 2004.

[30] T. Fukushima, K. Asaka, A. Kosaka, T. Aida, "Fully plastic actuator through layer-by-layer casting with ionic-liquid-based bucky gel," Angewandte Chemie International Edition, vol. 44 no. 16, pp. 2410-2413, DOI: 10.1002/anie.200462318, 2005.

[31] S. Saito, Y. Katoh, H. Kokubo, M. Watanabe, S. Maruo, "Development of a soft actuator using a photocurable ionic gel," Journal of Micromechanics and Microengineering, vol. 19 no. 3,DOI: 10.1088/0960-1317/19/3/035005, 2009.