Dynamic Finite Element Modeling and Simulation of Soft Robots

Abstract

Soft robots have become important members of the robot community with many potential applications owing to their unique flexibility and security embedded at the material level. An increasing number of researchers are interested in their designing, manufacturing, modeling, and control. However, the dynamic simulation of soft robots is difficult owing to their infinite degrees of freedom and nonlinear characteristics that are associated with soft materials and flexible geometric structures. In this study, a novel multi-flexible body dynamic modeling and simulation technique is introduced for soft robots. Various actuators for soft robots are modeled in a virtual environment, including soft cable-driven, spring actuation, and pneumatic driving. A pneumatic driving simulation was demonstrated by the bending modules with different materials. A cable-driven soft robot arm prototype and a cylindrical soft module actuated by shape memory alley springs inspired by an octopus were manufactured and used to validate the simulation model, and the experimental results demonstrated adequate accuracy. The proposed technique can be widely applied for the modeling and dynamic simulation of other soft robots, including hybrid actuated robots and rigid-flexible coupling robots. This study also provides a fundamental framework for simulating soft mobile robots and soft manipulators in contact with the environment.

1 Introduction

Soft robots have tremendous potential for application in various fields owing to their safety and flexibility embedded at the material level. The design and fabrication of soft robots represents a considerable challenge, and their effective behavior often arises from complex interactions among the controllers, morphology, and environment. This usually necessitates design, fabrication, and multiple designs, which require considerable time and resources [1]. Additionally, soft sensors and soft skin are other challenges for soft robots in the field of large deformation measurements and environmental perception [2]. To alleviate this problem, it is important to develop and apply effective and high-fidelity physical simulation tools [3]. However, it is arduous to establish a mechanical model for soft robots owing to the hyper-redundant degrees of freedom (DOFs), hyper-elasticity, and nonlinearity of their soft structures [4]. The strong nonlinearity and complex geometries of soft actuators hinder the development of analytical models for describing their motion. The nonlinear effects imply that extensive computational processes must be employed for simulations with high fidelity. Moreover, the dynamics of soft materials are difficult to simulate for the same reason. Generalized soft robots, particularly cable-driven continuum robots [5], have unique advantages for in-situ equipment repair/maintenance [6, 7] and medical applications [8]. However, great challenges remain in the modeling, control, and simulation of continuum robots interacting with environments or humans.

Various efficient simulation software tools are available in the field of traditional rigid robotics. The software currently used include ADAMS [9], Process Simulate [10], V-REP [11], GAZEBO [12], ROBOTRAN, and VEROSIM [13]. Commercial software such as ABAQUS [14, 15], ANSYS [16], COMSOL [17], and MARC [18] can be used for soft robot simulation. These software packages have quasi-static limitations and are mostly used in the simulation of structures or single bending modules. A simulation environment that can be applied for the testing, analysis, and optimization of soft robots still needs to be developed. The known tools and methods of robotic simulation have not yet been able to cope with various aspects in the field of soft robotics. For example, soft robotic simulation requires capabilities ranging from simulation of a single soft body over itself and multibody interaction to advanced robot/environment interactions [19]. Despite many difficulties, a number of researchers have made many contributions to this field, as discussed below.

Finite element (FE) modeling is a powerful numerical method for performing a piecewise approximation continuously with prior knowledge of the material properties, which provides an effective solution for predicting performance and optimizing soft actuator designs [20]. Unlike most analytical models, FE models can be easily adapted to different geometries and the deformations of a component can be readily visualized, leading to a better understanding of the influence of local strain on global actuator performance [15]. Moreover, FE models can provide deeper insights into the internal interactions inside a part, such as the interactions between layers of different materials. This rapid and efficient design framework reduces the cost and development time.

The FE method has been employed for the modeling and real-time control of soft robots using the open-source framework “Simulation Open Framework Architecture” with the “SOFT ROBOTS” plugin [21, 22]. However, it is limited to quasi-static conditions and requires the linearization of structural elasticity. Lipson et al. [1, 23] developed an open-source simulator (VoxCAD) with a GUI to simplify the modeling of the robots. However, it is impossible to approximate some geometrical shapes without an exaggerated number of voxels, which significantly increases the computation time [24]. Recently, Grazioso et al. [25] presented a geometrically exact model for soft continuum robots, and developed a dynamic simulation environment “SimSOFT” based on the FE method. Gazzola [26] developed a software “Elastica” to simulate the dynamics of soft filaments based on the Cosserat rod model. Similarly, Min et al. [27] developed a novel framework “SoftCon” for simulating and controlling soft-bodied animals based on the deep reinforcement learning algorithm. Hu et al. developed a physical simulator “ChainQueen” based on “Taichi” language for potential use in soft robot simulations[28]. Medvet et al. [29] presented a simulation tool for the optimization of 2-D voxel-based soft robots, “2D-VSR-Sim”, which can be used in VSRs by researchers from different disciplines.

Most previous simulation studies concentrated on soft robot fingers, manipulators, or single actuators. The multi-drive dynamic simulation of complex soft robots remains unfulfilled. In this study, the dynamic software RecurDyn was used to evaluate the deformation of soft robots. Various actuation methods have been used to drive soft robots in virtual environments. Superior to traditional FE analysis, the adopted multi-flexible-body dynamics technology can be used not only for system motion (position, velocity, and acceleration) and force analysis, but also for system control/optimization. To the best of the authors’ knowledge, this is the first time that dynamic FE simulations have been used in soft robot motion simulations. Our simulations and analysis can help understand the locomotion mechanics of soft robots and aid in the design, optimization, and damage prediction of complex multi-driven soft robots.

The remainder of this paper is organized as follows. In Section 2, the procedure for dynamic modeling and simulation using different actuation methods is presented in detail. Section 3 presents the hyperelastic constitutive models for soft silicon materials. Section 4 presents an experiment based on a soft-robot arm and spring-driven module study to corroborate the simulation model. Finally, conclusions and future work are presented in Section 5.

2 Modeling and Simulation Method

2.1 Cable-Driven Soft Robot Arm

A cable-driven soft robot arm inspired by an octopus was manufactured using silicone (EcoFlex 00-30) to verify the simulation model. As shown in Figure 1, four inextensible cables were embedded in the body of the silicone conical arm, and there were 13 image marked points fixed by silicone glue on the soft robot arm at intervals of 20 mm. The parameters indicated in Figure 1 are L = 260 mm, R_max = 17 mm, R_min = 8 mm, a = 10 mm, and b = 5 mm. Where L is the total length of the robot arm, and R_max and R_min are the radii of the base and tip anchorage cross sections, respectively. a and b are the distances of the cable from the midline at the base level and at the fixed tip level, respectively.

Figure 1

Structure of the soft robot arm prototype

2.2 Flexible-Body Modeling and Simulation

2.2.1 Modeling and Mesh

First, the users should establish a new model and set the units. One method is to draw the geometry in RecurDyn, the other is to import the CAD geometry from other 3D modeling software, such as PTC Creo, SolidWorks, or CATIA. The flexible bodies are then selected for meshing and the material is assigned to each of these parts. For the hyperelastic material, Solid8 elements (hexahedron elements) were used with auto or advanced mesh method in the Mesher tool to mesh the raw model. Hexahedral elements are more preferable to tetrahedral elements because they can better resemble the hyperelastic stress-strain relationship of a realistic material under large deformations [30]. The minimum average element size was defined by the structure size, which affects the calculation time; therefore, the element size and quantity are very important for the simulation.

2.2.2 Insert Boundary Conditions and Loads

Second, the flexible-surface contact option was used for the interaction area. Here, we take a cable-driven soft robot arm as an example. The cable was free to move into and out of the robot arm when pulled. Contact analysis was performed to predict the distribution of the stresses at the moment when the cable touched the soft material and avoid thrusting the cable into the rest of the robot. The friction between the cable and arm can be set as a constant, corresponding to the actual situation. Additionally, the base cross-section of the arm is fixed with the connection object, such as the prototype fixed on the testbed or robot body. The material properties and actuation methods are discussed separately in the next section.

2.2.3 Dynamic Analysis Settings

Certain parameters must be set before the simulation calculation analysis. One of the most important things is gravity, then set the “End time” and “Step”; other parameters include: “Maximum time step” (default as 0.01), “Initial time step” (default as 1E-6), “Error tolerance” (default as 5E-3), and take the value from 0 to 1 as “Numerical damping”. Because there are many nodes, the total amount of output data would be too large if the data were automatically saved for every node. Instead, the output data are saved only for nodes that are selected by users.

2.2.4 Submit Job and Postprocessing

Before the final calculation, a pre-analysis should be performed to check if the constraints are redundant and to find improper connections and unconstrained components. The simulation job is then submitted to solve for the results. The results are evaluated and processed after the simulation is successful. The final results include displacement, velocity, acceleration of the preset nodes, contact force and torque of the contact surface or joint, strain and stress of the flexible bodies, and the change in the load applied to the actuator. These data can be used to verify the analytical models or further optimize the design of soft robots. The overall simulation procedure is illustrated in Figure 2.

Figure 2

Simulation analysis flow chart in RecurDyn

2.3 Actuation Modeling

Rigid robots are primarily driven by motors or hydraulics. There are various actuators or materials for soft robots [31, 32], such as pneumatic, smart materials (including shape memory materials, hydrogel, liquid crystal elastomer, dielectric elastomer), chemical reaction drives, physical fields (magnetic field or light field), and motor or cables can also be used in soft robots. Here, we demonstrate the modeling methods for the cable drive, spring drive, and pneumatic drive.

2.3.1 Cable-driven Modeling Methods

Cables are a driving method for soft robot actuation and soft manipulation. This is becoming increasingly important in soft robot deformation and locomotion. For example, the cable-driven soft robot arms [33], the biped “Flippy” robot [34], and a soft starfish-like robot with a motor-tendon actuator [35]. Unlike rods, cables do not present shear strains because of their low cross-section thickness, and it is difficult to simulate cables in a virtual environment as they cannot bear shear force and axial pressure, not to mention interacting with soft materials. Here we present three methods for simulating the flexible cables. As shown in Figure 3, in the first method, the cables are modeled by thin quadrangular rods using Solid4 (tetrahedron) elements with an orthotropic elastic material [4]. The Young's modulus in the direction of the length was 1000 MPa, whereas the Young's modulus perpendicular to the length direction was 1 MPa. The second method used a beam element. This element supports isotropic materials only. Young's modulus was set to 1000 MPa. It is worth noting that the solid element with orthotropic material is not supported in the new version of V9R3, and the beam element modeling method in V9R3 is more stable than the old versions. The last approach to modeling flexible cables is using piecewise multiple rigid bodies to subdivide a cable into small pieces. Each rigid body can be approximated by a cylinder or prism, while each segment is connected by spherical joints. It has been used in wire rope modeling, but this method is time consuming in multi-flexible body simulations because a large number of rigid-flexible contacts are required. Consequently, the quadrangular rod using a solid element and the line using beam element are used to simulate flexible cables in this paper.

Figure 3

Modeling methods of flexible cable

2.3.2 Spring Modeling

Shape memory alloys (SMAs) have been widely employed in soft robots owing to their high work density, large recovery stress and strain, silent operation, low stimulating voltage, corrosion resistance, and biocompatibility[36]. They have been used in soft robot arm inspired by octopuses [37], the “GoQBot” [38], and starfish robot with multi-gaits [39]. However, there is little research on SMA spring-drive simulations. Using the spring force function in RecurDyn, we can model SMA springs to actuate the soft robot arm elongation or bending by contracting the springs in a radial direction, similar to octopus’ transverse muscles [37]. First, the base point and action point determine the location of the spring. Then, the spring diameter, coil diameter, and number of coils define the spring geometry characteristics. Finally, the stiffness coefficient, damping coefficient, and free length were used to define the mechanical properties of the spring.

A mechanical model is necessary to predict of the performance of the SMA springs. The linear-elastic model is based on Castigliano’s second theorem and gives the spring rate K (for austenite and martensite) as the output [40, 41]:

$$K = \frac{{G_{a} d^{4} }}{{8nD^{3} }},$$

(1)

where $G_{a} = E/[2(1 + \nu )]$ is the shear modulus, ν is Poisson’s ratio, E is Young’s modulus, d is the wire diameter, D is the average spring diameter, and n is the number of coils.

2.3.3 Pneumatic Driving

Pneumatic actuation is the most common driving method, for pneumatic soft manipulators, soft fingers/hands, and soft grippers [20, 42], as well as for various pneumatic locomotion robots [43, 44]. Here, we used the pressure load function in the flexible toolkit to generate pneumatic actuation. First, we selected surfaces with internal cavities on the meshed model and set them as “Patches”. Second, we defined the pressure load using the pre-input expression. Subsequently, the pressure load on the defined patches was added. It is noteworthy that users should change the “tolerance angle” set (0 – 90) to select the complex inner cavity surface and change the pressure direction up or down according to the real model.

3 Hyperelastic Material Models

The silicon used to make the soft robot is a hyperelastic material, and its stress-strain relationship is nonlinear. The material properties of the soft arm have a significant influence on the simulation results. Four constitutive models of hyperelastic materials are provided in the RecurDyn software, the Arruda-Boyce model, the Neo-Hookean model, the Ogden model [45], and the Mooney-Rivlin model [46]. The Neo-Hookean and Mooney-Rivlin model are based on the linear approximations of the strain invariants from the Ogden model. Although they may be accurate in these low-strain regimes and improve the solving speed, the accuracy of these simplified models is limited at higher strains [47]. The Arruda-Boyce model is based on thermodynamic statistics that require a variety of experimental tests to determine the properties of the material [48]. The constitutive model expression for the Ogden model is as follows:

$$U = \sum\limits_{i = 1}^{N} {\frac{{\mu_{i} }}{{\alpha_{i} }}} (\lambda_{1}^{{\alpha_{i} }} + \lambda_{2}^{{\alpha_{i} }} + \lambda_{3}^{{\alpha_{i} }} - 3),$$

(2)

where μ_i and α_i are the primary material constants, N is the number of terms, λ_i is the principal stretches. The six-parameter model (N = 3) is the most commonly used for large strain problems, that is, at or above 400%.

When N = 1 and α = 2, the Ogden model transforms into the Neo-Hookean model:

$$U = C_{{1{0}}} (I_{1} - 3)= \frac{\mu }{{2}}(I_{1} - 3),$$

(3)

where μ is the shear modulus.

When N = 2, α₁ = 2, and α₂ = −2, the Ogden model transforms into the Mooney-Rivlin model, which is expressed as

$$U = C_{1} (I_{1} - 3) + C_{{2}} (I_{2} - 3),$$

(4)

where C₁ and C₂ are material constants determined by experiments or experience, I₁ and I₂ are the first and second invariants of the stress tensor, $I_{1} = \lambda_{1}^{2} + \lambda_{2}^{2} + \lambda_{3}^{2}$, $I_{2} = \lambda_{1}^{2} \lambda_{2}^{2} + \lambda_{2}^{2} \lambda_{3}^{2} + \lambda_{1}^{2} \lambda_{3}^{2}$. This model is used for moderate deformations, that is, lower than 200%.

The Yeoh model is another frequently used model for large-strain problems [20], which is similar to the Mooney-Rivlin model, but it does not include the second tensor invariant, making it simpler than the Mooney-Rivlin model. The second-order Yeoh model is expressed as

$$U = C_{1} (I_{1} - 3) + C_{2} (I_{{1}} - 3)^{{2}} .$$

(5)

Silicone as EcoFlex 00-30 is an incompressible nonlinear hyperelastic material that is the most commonly used in soft robot structure, with a Poisson ratio of 0.5, and its main elongation satisfies λ₁λ₂λ₃ = 1. The material properties of the three-term Ogden model are as follows: μ₁ = 0.024361, μ₂ = 6.6703 × 10^-5, μ₃ = 4.5381 × 10^-4, α₁ = 1.7138, α₂ = 7.0679, α₃ = −3.3659, all the μ terms have the units of MPa, and all the α terms are dimensionless [47]. In the actual simulation, if μ₃ is set to be negative to make μ_i·α_i greater than 0, the stability of the algorithm can be ensured. According to the relationship between these models, the material properties of the Mooney-Rivlin model are as follows: C₁ =0.01218 MPa, and C₂ = −3.33515 × 10^-5 MPa. The material property of the Neo-Hookean model is C₁₀ = 0.01218 MPa. Additionally, the Ecoflex 00-50 material using the three-term Ogden model parameters are as follows: μ₁ = 0.1079, μ₂ = 2.147 × 10^-5, μ₃ = −0.0871, α₁ = 1.55, α₂ = 7.86, α₃ = −1.91 [14], all the μ terms have the units of MPa, and all the α terms are dimensionless.

In this study, all simulations were performed using an Intel^® Core^TM i7-8700 central processing unit (CPU) at 3.20 GHz with a GTX1060 graphics processing unit (GPU) with 16.0 GB random-access memory.

4 Simulation and Experimental Results

4.1 Simulation of Pneumatic Soft Modules

FE simulations can be used to optimize the structure design of different materials, and dynamic simulations can be used to analyze the deformation and motion performance. As shown in Figure 4(a) and (b), a circular cross-section was modeled with an external diameter of 25 mm and 65 mm in length, 60 mm in chamber length, 3.7 mm in chamber radius, and a distance between the chamber wall and the outer wall of 1.5 mm. Here, the influence of gravity was not considered. Figure 4(c)–(e) shows the simulation results using the Ecoflex 00-30 material, where θ is defined as the bending angle of the module, and Figure 4(f)–(h) shows the simulation results using the Ecoflex 00-50 material. The input pressure was 0.01 × t (MPa) for the 00-30 module and 0.1 × t (MPa) for the 00-50 module, respectively. Figure 4(e) and h shows the semi-sectional view of the chamber, where we can check the internal deformation and predict where damage might occur according to the stress distribution.

Figure 4

3-Dimensional model of circular cross-section and the simulation results.

Figure 5 presents the time change of the bending angle and the end-center point motion velocity, which corresponds to Figure 4. It is clear that the softer material Ecoflex 00-30 has better bending properties and requires less pressure to deform than Ecoflex 00-50 does. The expansion of the air chamber using Ecoflex 00-50 is more distinct. According to previous studies, using fibers or different material layers can further reduce the bulging effect and achieve good bending capability [14, 15].

Figure 5

Time change of the bending angle and the end-center point motion velocity

4.2 Axial Contraction and Elongation

Using the cable-driven method, we can simulate the longitudinal muscles that are arranged along the arm to make the soft arm bend in space or contract in the axial direction. It is difficult to achieve pure longitudinal contraction in experiments owing to manufacturing and actuating errors, and these difficulties can be overcome in a simulation environment. When two opposite cables or four cables were driven with the same tension, the soft arm contracted along the axial direction. Figure 6(a) shows the initial state of the simulated model. The results for T ₁ = T₃ = 2 N, and T₂ = T₄ = 0 N are shown in Figure 6(b) and (c), and the results for T₁ = T₂ = T₃ = T₄ = 2 N are shown in Figure 6(d) and (e), where T_i (i = 1, 2, 3, 4) is the tension force applied to the cables that are distributed clockwise at π/2. Figure 6(b) and d shows the results with gravity, and Figure 6(c) and (e) shows the results without gravity. The pull force of each cable was increased from 0 to 2 N over time (T = 2 × t) when t ≤ 1 s, else T_i = 2 N.

Figure 6

Simulation results of two cables driven and four cables driven constriction simulation

Figure 7 presents the dynamic changing process of the arm contractional length corresponding to the cases shown in Figure 6. It can be observed that the relationship between the contraction length of the arm and tensile force is approximately linear in the two cases. The FE model could quantitatively predict the contraction capabilities of a real prototype. When gravity is added to the simulation, oscillations are present at the beginning, while it still approximates a linear relationship in the main contraction deformation.

Figure 7

Time change of contraction length in simulation

To verify the accuracy of the spring drive method, a cylindrical deformation module was constructed using smooth-on 00-30 silicon. Figure 8(a) and (b) shows the front and top views of the initial state of the deformation module. The main characteristics of the module included an external diameter of 40 mm, an inside diameter of 20 mm, and an altitude of 22 mm. There are two SMA springs in the vertical distribution: one 8 mm below the upper surface and the other 6 mm above the lower surface. For the SMA spring, the wire diameter was 0.5 mm, and the average spring diameter was 2.4 mm. The number of coils was 10, and the maximum effective working length was 8 mm. The shear modulus G_a was approximately 14.1 GPa for the martensitic phase and 24.7 GPa for the austenitic phase [49]. The corresponding stiffness rate calculated by Eq. (1) is 0.8 N/mm for the martensitic phase and 1.4 N/mm for the austenitic phase. The springs were driven by a direct current power supply of 12V/1A. Figure 8(c) and (d) shows the front and top views of the post-deformation state of the deformation module; the height increased to 28 mm, the minimum outside diameter was 26 mm, and the maximum outside diameter was 47 mm.

Figure 8

SMA spring actuation experiment and simulation

Figure 8(e)–(h) represents the FE model and simulation results. The FE model had the same structure as the soft actuator prototype. Because the spring requires a base body and an action body, each spring is divided into two parts and connected by a tiny point of mass. The stiffness is set as a linear function that changes with the spring contraction length. As the speed of the spring contraction in the experiment depends on the current, the influence of spring damping was not studied in the simulation, and the damping coefficient was set to the default value of 1. The simulation result is as follows: the height increased to 28.96 mm, the minimum outside diameter was 25.55 mm, the maximum diameter was 46.92 mm. The relative error in the height change between the simulation and experimental results was 4.4%.

4.3 Experiments of the Soft Robot Arm

4.3.1 Free Bending of the Cantilever Arm

When a soft arm is placed horizontally in air like a cantilever beam, it bends freely under the action of gravity until the initial potential energy is completely consumed. The free bending of the soft robot arm presents a damped oscillation effect and is a simple method to confirm the accuracy of the FE models, particularly, to describe the different responses of various hyperelastic material models. As shown in Figure 9, the simulation assigned different material models, including the Ogden model, the Mooney-Rivlin model, and the Neo-Hookean model, to the cantilever model. The base section is fixed with ground coordinate and only gravity acts on the arm. The simulation duration was set to 30 s, and the time step was 0.01 s.

Figure 9

Simulation results of the cantilever arm free bending

Figure 10(a) shows the experimental results obtained using the soft arm prototype. Figure 10(b) shows a comparison between the FE simulation and experiment. The tip point errors with the experiment and time consumption of different material models are listed in Table 1. It is obvious that the Ogden model is the most accurate and time-consuming among the three models. In this study, the Ogden model was used by default, unless otherwise mentioned. Additionally, the relative error is defined as the ratio of the absolute distance between two endpoints to the arm length.

Figure 10

Experimental results of the cantilever arm and comparison with the simulation

Table 1 Tip point errors with experiment and time-consumption of different material models

4.3.2 Plane Bending of the Soft Arm

When only one cable is driven, the soft arm exhibits planar bending. As shown in Figure 11, the cable was modeled using a beam element, and the contour color on the arm represents the Von Mises stress.

Figure 11

Simulation results of one-cable-driven arm using beam element

Figure 12(a) shows snapshots of the experimental results obtained using the soft arm prototype with different tensile forces in underwater conditions. Figure 12(b) shows the experimental and simulation results using different cable modeling methods in one picture, where the labeled tension represents the experimental value. It was found that the measured tension of the cable in the experiment was significantly greater than the value in the FE simulation because of cable friction, as we did not consider friction in the simulation. Another phenomenon is that it generates obvious contraction in the FE simulation owing to the soft material hyperelasticity. Therefore, the boundary of the workspace was smaller than that in the experimental results as shown in Figure 12(b). Moreover, because part of the tension force was used to make the arm yield contraction at the beginning of the simulation, the tension force did not change continuously in the simulation. The average relative error of the tip point between the solid element cable method and experiment was 2.90%, and that for the beam element method was 5.97%. Here, the tip position errors between the simulation and experimental results are defined as

$$Error = (|p_{s} - p_{e} |/L) \times 100\% ,$$

(6)

Figure 12

Experimental results of one-cable-driven arm and comparison with simulation

where p_s and p_e are the tip point coordinates in the simulation and experiment, respectively. L is the length of the soft arm.

As shown in Figure 13, when two adjacent cables are driven equally, the soft arm also yields planar bending in the 3D space. The input tension is T = 2 × t for each cable. The contour color in Figure 13(a) represents the Von Mises strain and the contour color in Figure 13(b) represents the Von Mises stress.

Figure 13

Simulation results of two-cable-driven soft robot arm

Figure 14(a) shows snapshots of the experimental results obtained using the soft arm prototype with different tensile forces in underwater conditions. We simultaneously obtained pictures from two perpendicular directions, and then calibrated the two groups marked points to obtain the approximate curve of the soft arm middle line. As shown in Figure 14(b), the experimental and simulation results obtained using different cable modeling methods are drawn in one photograph. The simulation coordinate results also used the average values of the two labeled edges. The average relative error of the tip point between the solid element cable method and experiment was 1.67%, and that for the beam element method was 1.63%.

Figure 14

Experimental results of two-cable-driven arm and comparison with simulation

When the soft arm was underwater, we assumed that buoyancy and gravity were equivalent; therefor, gravity was not considered in the plan bending simulation. The computation time of different simulation modeling methods is still worth noting. In the one-cable-driven simulation for 1.5 s, the solid element cable model used 21 h and 1.22 min, whereas the beam element cable model used 41 h and 48.67 min. Additionally, the Mooney-Rivlin model and Neo-Hookean can be applied in the case of small deformations (for example, when tension is smaller than 3–4 N in one-cable-driven experiment) to find a trade-off between the accuracy and time consumption, as it can save more than half of the calculation time under the same conditions.

5 Conclusions and Future Work

In this study, a multi-flexible body dynamic simulation technique was innovatively applied to soft robots. A dynamic FE simulation process for soft robots was proposed systematically, considering a hyperelastic constitutive model for soft materials.

• Three driving modeling methods (cable-driven, spring actuate, and pneumatic) for soft robots are proposed and applied. The pneumatic method was shown using a bending actuator. Springs were used on the cylindrical module. Cables using solid elements and beam elements were used in the cable-driven soft robot arm.

• A cable-driven soft robot arm and cylindrical SMA spring actuated module inspired by an octopus were manufactured. The soft arm exhibited free bending and cable-driven plane bending. The SMA module exhibited elongation ability. Relative FE simulation models of the prototypes were established. The experimental results show that the accuracy of the simulation is acceptable, which demonstrates the potential applications of the FE simulation.

• In the future, we will strive to build a hybrid actuated soft-legged robot and study its mobility and trafficability in a virtual environment.