A material point-based simulation method for soft robots with free boundary interactions

Abstract

Soft robots show an advantage when conducting tasks in complex environments due to their enormous flexibility and adaptability. However, soft robots suffer interactions and nonlinear deformation when interacting with soft and fluid materials. The reason behind is the free boundary interactions, which refers to undetermined contact between soft materials, specifically containing nonlinear deformation in air and nonlinear interactions in fluid for soft robot simulation. Therefore, we propose a new approach using material point method (MPM), which can solve the free boundary interactions problem, to simulate soft robots under such environments. The proposed approach can autonomously predict the flexible and versatile behaviors of soft robots. Our approach entails incorporating automatic differentiation into the algorithm of MPM to simplify the computation and implement an efficient implicit time integration algorithm. We perform two groups of experiments with an ordinary pneumatic soft finger in different free boundary interactions. The results indicate that it is possible to simulate soft robots with nonlinear interactions and deformation, and such environmental effects on soft robots can be restored.

1. Introduction

Soft robotics has become a popular research field. Compared to traditional robots, soft robots can achieve large deformations while ensuring interaction safety, making them advantageous in areas such as flexible grasping [1, 2], medical rehabilitation [3–5], underactuated finger [6, 7], and operations in unstructured environments [8, 9].

However, since the actuation of soft robots is often intrinsic to their structure, their kinematics result from a coupling of geometric relationships and complex mechanical mechanisms, making it difficult to directly solve for design parameters based on the desired motion. During the design process of soft robots, the researchers need to use simulations to evaluate the performance and behavior of the prototype, so that the soft robot works appropriately in the consistency of expectation. As soft materials such as rubber and hydrogel provide soft robots with enormous flexibility and adaptability, the behaviors of soft robots exhibit nonlinear interactions and deformation in complex environments and are hard to predict. Thus, the physical simulation based on continuum mechanics is adopted. Simulation of soft robots can be regarded as mathematics of physics. It contains continuum mechanical, geometrical, discrete, and surrogate models [10]. Continuum body assumption is made that material is continuously distributed in the occupied space of an object and can be divided into infinitesimal elements of identical properties [11]. Based on continuum body assumption, solving the partial differential equation derived from constitutive equations such as conservation of mass and energy leads to physical simulation based on continuum mechanics. Finite element analysis (FEA) is a widely used Galerkin method for solving PDEs numerically and has been applied to continuum mechanics successfully for a long time. A soft finger for hand rehabilitation was developed, and its deformation was evaluated with FEA [11]. Another work Comber et al. [12] shows that a soft robot can drive the surgical needle, and they inspected the actuator displacement with FEA. An iterative algorithm was developed to optimize the structural design of pneumatic soft robots [13], and FEA was used to quantify the deformation of the intermediate optimized result. These approaches show that the physical simulation based on continuum mechanics supports and promotes the development of soft robots.

However, these approaches do not consider the effect of complex environments. The shortcoming of FEA, which is the incapability to handle free boundary interactions, does not support such consideration. For FEA, the two-step contact formulation procedure [14] and slave–master body technique [15] are used to handle naive fixed boundary coupling by fixing elements position and applying contact force on the contact boundary. A two-step scheme is adopted to solve solid and fluid systems separately and alternately [16] when the contact problem comes to a scenario with both fluid and solid materials. Aforementioned methods fail in such interactions and deformation scenarios because the contact cannot be determinant before simulation. Thus, FEA is incapable of simulating the scenario of soft robots in complex environments, and a new approach is required.

Meshless methods, that is, material point method (MPM) [17] and smoothed particle hydrodynamics (SPH) [18], are significantly different from FEA because they can autonomously solve the coupling of free boundary. The reason is how the shape function is supported. In FEA, a shape function is supported by a set of grouped nodes (which is an element). When boundary contact happens, two elements will go across as they will not contribute to the shape functions of each other. While in the meshless methods, all nodes (in favor of using the terminology of particles) within the support range will contribute to the same support function, which facilitates the information exchange between different phases and objects of materials. Though the shape function of SPH and MPM are supported similarly, SPH requires extra computation to search supporting particles compared to MPM. Research shows that MPM takes advantage of SPH in both computational efficiency and accuracy [19]. For such reasons, we use MPM to develop our approach for the simulation of soft robots. The proposed approach can handle nonlinear interactions and nonlinear deformation of soft robots.

Recently, many methods [20] and frameworks have been proposed to address the interaction problems between free boundaries and the environment. Differentiable simulators have also been integrated with gradient-based optimization algorithms specifically developed for soft robotics applications. DiffPD [21], a differentiable soft-body simulator based on finite element method (FEM), is designed for efficient soft-body learning and control and can be used to simulate an underwater starfish robot [22]. SoftCon [23], another differentiable soft-body simulator based on FEM, serves as a versatile framework for designing and controlling underwater soft-bodied organisms. ChainQueen [24, 25] and its follow-up work DiffTaichi [26] introduce differentiable physics-based soft-body simulators using the MPM. There are also other differentiable frameworks, such as Elastica [27]. In addition, many researchers use fluid–structure interaction (FSI) to simulate the deformation of soft robots in fluid environments. For instance, Xavier et al [28] employ bidirectional FSI to study the deformation of soft robots underwater impact, while Soomro et al [29] use FSI to investigate the interaction between a soft biomimetic frog robot and water.

These methods have made significant contributions to soft robot simulation with free boundary interactions. However, MPM has advantages over other methods in handling free boundary interactions in soft robotics. MPM naturally handles large deformations and interactions, which are common in soft robotics but are often not modeled in mesh-based approaches due to computational expense. As mentioned earlier, by utilizing a mesh-free approach, MPM handles complex boundary problems more effectively compared to mesh-based methods like FEM [21, 23] and Coupled Eulerian–Lagrangian (CEL). Compared to other mesh-free methods like SPH, MPM offers higher accuracy and computational.

Although ChainQueen has already applied MPM to soft robot simulation, our method offers greater stability. The explicit time integration in their method requires small timesteps to maintain numerical stability and leads to significant numerical errors when the timestep is large. In contrast, our MPM framework uses implicit time integration, which ensures numerical stability and guarantees higher numerical accuracy. Additionally, we have introduced automatic differentiation to relatively improve computational efficiency and accuracy.

In this study, the main contribution in this article can be described in two aspects:

• • We address the open problem of simulation of soft finger with free boundary interactions. The proposed approach can be applied automatically without complicated algorithm design, especially in nonlinear deformation and interaction cases.

• • We integrate automatic differentiation and implicit time integration into the modeling framework, which makes it more computationally efficient and stable, respectively.

Our key strengths include the ability to manage large deformation problems and complex boundary conditions with excellent stability. Although this process involves several simplifications when extending to multiphysical phenomena, such as Rayleigh damping and Coulomb friction, it can be theoretically proven that their impact is negligible in most cases. We have also experimentally verified that the simulation results remain highly accurate even after disregarding these factors.

The rest of the paper is organized as follows. In section 2, we will discuss the dynamic model principles for simulation of soft robot. It contains MPM framework and our modification in detail. Then, we illustrate the implemented algorithm in section 3. The simulation and experiment results will be shown in section 4.

Figure 1. An illustration of MPM update framework. From left to right: particle to grid, grid update, grid to particle and reset grid state. [17]

2. Dynamic modeling for simulation of soft robot

2.1. MPM based modeling framework

When conducting simulation via MPM, the space is discretized by background grids, while the object is discretized by particles. The physical information is stored on particles, but grids enable such physical information to be updated. The MPM can perform both dynamic and static simulations, but for the simulation of soft robots, the dynamic simulation is more important since inertia affects the robot’s behavior. Each time step of a dynamic MPM simulation involves three steps [24]. In the first stage, the physical properties of particles are transferred to the adjacent grids (spatial points at which the shape functions are centered). Mass, momentum, and external forces are transferred from particles to grids by $m_{i}=\sum _{p}\alpha _{ip}m_{p}$ , $m_{i}\boldsymbol{v}_{i}^{n}=\sum _{p}\alpha _{ip}m_{p}\boldsymbol{v}_{p}^{n}$ , and $\boldsymbol{f}_{i}^{ext}=\sum _{p}\alpha _{ip}\boldsymbol{f}_{p}^{ext}$ , where variables with subscript $p$ is the physical information on particles and variables with subscript $i$ is the intermediate variables on grids, $\alpha _{ip}$ is the weight between specific particle $p$ and grid $i$ . In the second stage, the internal force $\boldsymbol{f}_{i}^{{int}^{n}}$ is computed and intermediate variables on grids are to be updated using explicit time integration, $\boldsymbol{v}_{i}^{n+1}=\boldsymbol{v}_{i}^{n}+{\Delta} t\boldsymbol{a}_{i}^{n}$ , where $\boldsymbol{a}_{i}^{n}=(\boldsymbol{f}_{i}^{{int}^{n}}+\boldsymbol{f}_{i}^{ext})/m_{i}$ . In the third stage, the intermediate variables on grids are transferred back to update the particles as $\boldsymbol{v}_{p}^{n+1}=\sum _{i}\alpha _{ip}\boldsymbol{v}_{i}^{n+1}$ and $\boldsymbol{x}_{p}^{n+1}=\boldsymbol{x}_{p}^{n}+{\Delta} t\boldsymbol{v}_{p}^{n+1}$ . These stages are combined and referred as update framework of MPM. After these main stages, intermediate variables on grids are to be reset. An illustration is shown in Figure 1.

However, numerous MPM implementations have modified those above equations. Not all MPM frameworks are suitable for soft robot simulation. The most representative frameworks of MPM are fluid implicit particle (FLIP) and particle in cell (PIC). However, PIC and FLIP have flaws that cannot be applied directly in the simulation of soft robots. The PIC technique does not guarantee the conservation of angular momentum because it overlooks the local rotational motion, which severely dissipates energy. As soft robots often suffer from nonlinear deformation of shear and rotation, PIC method cannot handle the simulation correctly. FLIP was proposed to solve such a flaw in the simulation of fluids, but angular momentum is not conserved unless using a full mass matrix [30]. Thus, FLIP is not a proper framework for the simulation of soft robots.

Affine particle in cell (APIC) introduced affine velocity field tensor to store the angular momentum on particles, making APIC free from energy dissipation when simulating deformable objects. As the affine velocity field is added into the framework, it needs to be computed and updated for numerical computation. In the third stage, the updated affine velocity field $\boldsymbol{C}_{p}^{n+1}$ is computed as $\boldsymbol{C}_{p}^{n+1}=4\sum _{i}\boldsymbol{v}_{i}^{n+1}(\boldsymbol{x}_{i}-\boldsymbol{x}_{p})^{T}\alpha _{ip}/h^{2}$ , it is used to update the deformation gradient $\boldsymbol{F}_{p}^{n+1}=(\boldsymbol{I}+{\Delta} t\boldsymbol{C}_{p}^{n+1})\boldsymbol{F}_{p}^{n}$ , where I represents the identity matrix and which will be used to compute internal force. The proof of exact angular momentum conservation has been reported in a previous work [31]. APIC is more stable when compared with PIC and FLIP [32]. Moving least squares material point method (MLS-MPM) [24, 33] simplifies the computation of APIC framework by introducing the moving least squares function into interpolation. Our MPM soft robot simulation approach is based on APIC and MLS-MPM.

2.2. The nonlinear material models with automatic differentiation

Automatic differentiation is a technique to compute the derivatives of functions specified by a program [34]. When a mathematical function is computed, it executes a consecutive series of basic operation and functions. By using chain rule, the overall derivative of the given function is the product of all derivatives during each operation. Moreover, this idea makes sure the result is purely analytical when compared to the finite difference method.

We mainly focus on soft solids and Newtonian fluids as they are ubiquitous in soft robotic systems. Soft robots are often fabricated by materials like hydrogel, silicon-based rubber, and some sorts of polymers. Hyperelastic models can restore the nonlinear stress-strain relationship exhibited by those materials, which fails in linear elastic models. We use Mooney-Rivlin, Ogden, and Yeoh [35] model to exhibit the simulation of soft robot. Their constitutive equations are listed in Table 1 (In the table, d refers to the dimensionality of the simulation). The simulation results based on three material models are shown in Figure 2 (visualized using Matplotlib [36]).

Table 1. Constitutive equations for solids.

Figure 2. The simulation results of a bending soft finger based on three material models: (a) Monney-Rivlin. (b) Ogden. (c) Yeoh.

For Newtonian fluids [37], the authors mentioned that using fixed corotational model (which is a SEDF) is capable to simulate Newtonian fluid as $\Psi =\lambda /2(J-1)^{2}$ , where $\lambda$ is bulk modulus for the fluid and J is the determinant of the deformation gradient. Furthermore, taking the derivative of this function yields a relation identical to the equation of state (EOS) used in SPH simulations [38]. For fluids simulation, the deformation gradient needs to be reset as $\boldsymbol{F}=J^{\frac{1}{d}}\boldsymbol{I}$ .

To incorporate these material models into the MPM framework. The relation exists in MLS-MPM [33] as

(1) \begin{align}\boldsymbol{f}_{i}^{{int}^{n}}=-4/h^{2}\sum _{p}V_{p}^{0}\partial \Psi /\partial \boldsymbol{F}_{p}^{n}\left(\boldsymbol{F}_{p}^{n}\right)^{T}\left(\boldsymbol{x}_{i}-\boldsymbol{x}_{p}\right)\alpha _{ip} \end{align}

Therefore, the derivative $\partial \Psi /\partial \boldsymbol{F}$ is required. A manual calculation is usually performed to obtain this value. However, for some complex SEDF, such a calculation is not feasible. In these cases, the strain energy of the local volume of the deformed soft robot can be computed over the entire volume of the soft robot to obtain the total strain energy.

(2) \begin{align}e=\int _{\Omega _{0}}\Psi (\boldsymbol{F}(\boldsymbol{X}))d\boldsymbol{X}=\sum _{p}V_{p}^{0}\Psi \left(\boldsymbol{F}_{p}\right) \end{align}

Accordingly, the strain-based internal force is the derivative of total strain energy to the location of space grid $\boldsymbol{x}_{i}$ :

(3) \begin{align}\boldsymbol{f}_{i}^{{int}^{n}}\left(\hat{\boldsymbol{x}}_{i}\right)=-\partial e/\partial \boldsymbol{x}_{i} \end{align}

where the location of space grid $\boldsymbol{x}_{i}$ is the function of space grid velocity $\boldsymbol{v}_{i}$ , which is

(4) \begin{align}\boldsymbol{f}_{i}^{{int}^{n}}\left(\hat{\boldsymbol{x}}_{i}\left(\boldsymbol{v}_{i}\right)\right)=-\left(\partial \boldsymbol{v}_{i}/\partial \boldsymbol{x}_{i}\right)\left(\partial e/\partial \boldsymbol{v}_{i}\right)=-\left(1/\Delta t\right)\left(\partial e/\partial \boldsymbol{v}_{i}\right) \end{align}

As we measure the internal force before the grid updating, thus the grid velocity is $\mathbf{0}$ , as it is an intermediate variable. We can get the internal force by computing $e$ and its derivative using automatic differentiation. This modification is used in explicit time integration framework.

2.3. The implicit time integration

As mentioned previously, the computation of MPM framework relies only on the current states of the system (known as the explicit time integration). As a result of explicit time integration, numerical errors accumulate, reducing the simulation’s stability and accuracy in three aspects:

1) the material parameters are significant in numbers (from MPa to GPa for common materials), which increases this error.

2) the interval between time steps needs to be small enough to compensate the accumulation error.

3) there is always a tradeoff between material parameters and accuracy to get a robust simulation.

A solution to these problems is implicit time integration. Unlike the explicit framework, implicit time integration requires solving an equation involving both the current time step of the system and the next time step. Usually, it is solved by an optimization algorithm with finite difference method providing gradient information at the cost of accuracy and running time. We incorporate automatic differentiation to overcome the flaws of finite difference method using implicit time integration. The difference in stability between explicit and implicit time integration is illustrated using a simulation case of a cube freely falling under gravity with a large time step. The results, shown in Figure. 3, demonstrate that MPM with an explicit framework fails when using a large time step, whereas MPM with an implicit framework and automatic differentiation proves to be more robust, improving the stability of the simulation.

Figure 3. The difference of explicit and implicit time frameworks’ simulation results: (a) MPM using explicit framework fails when the simulation with large time step and large material parameters. (b) MPM using implicit framework and automatic differentiation improves the robustness of simulation.

A previous work [39] describes the implicit time integration scheme of APIC. However, our approach is innovative in three aspects. First, we use grid velocity as independent variable rather than intermediate grid position, which is more accurate and precise. Second, our method is based on MLS-MPM. Lastly, the adoption of automatic differentiation makes our approach to compute derivatives much faster and accurate.

The only difference between implicit and explicit time integration is the update of grid velocity. For an implicit time integration, the update of grid intermediate variable $\boldsymbol{v}_{i}$ becomes [40]:

(5) \begin{align}\boldsymbol{v}_{i}^{n+1}=\boldsymbol{v}_{i}^{n}+{\Delta} t\boldsymbol{a}_{i}^{n+1} \end{align}

According to kinematics, $\boldsymbol{a}_{i}^{n+1}=(\boldsymbol{v}_{i}^{n+1}-\boldsymbol{v}_{i}^{n})/{\Delta} t$ , and dynamics, $\boldsymbol{a}_{i}^{n+1}=(\boldsymbol{f}_{i}^{{int}^{n+1}}+\boldsymbol{f}_{i}^{ext})/m_{i}$ . The two expressions are equal:

(6) \begin{align}\boldsymbol{h}=\left(\boldsymbol{v}_{i}^{n+1}-\boldsymbol{v}_{i}^{n}\right)/{\Delta} t-\left(\boldsymbol{f}_{i}^{{int}^{n+1}}+\boldsymbol{f}_{i}^{ext}\right)/m_{i} \end{align}

The roots which satisfy $\boldsymbol{h}=\mathbf{0}$ is going to be calculated. By solving this equation, the grid velocity is then updated. In Eq. 6, there are two variables of next time step $n+1$ , $\boldsymbol{v}_{i}^{n+1}$ and $\boldsymbol{f}_{i}^{{int}^{n+1}}$ , but the internal force is dependent on grid velocity. To represent the internal force in terms of grid velocity, we begin with the total strain energy:

(7) \begin{align}e\left(\boldsymbol{x}\right)=\int _{\Omega _{0}}\Psi (\boldsymbol{F}(\boldsymbol{X}))d\boldsymbol{X}=\sum _{p}V_{p}^{0}\Psi \left(\boldsymbol{F}_{p}\right) \end{align}

The updated intermediate grid position is $\hat{\boldsymbol{x}}_{i}=\boldsymbol{x}_{i}+{\Delta} t\boldsymbol{v}_{i}$ , and the stress-based force on grid $\boldsymbol{f}_{i}^{int}(\hat{\boldsymbol{x}}_{i})$ is

(8) \begin{align}-\frac{\partial e}{\partial \hat{\boldsymbol{x}}_{i}}=-\partial \sum _{p}V_{p}^{0}\Psi \left(\hat{\boldsymbol{F}}_{p}\left(\hat{\boldsymbol{x}}_{i}\left(\boldsymbol{v}_{i}\right)\right)\right)/\partial \hat{\boldsymbol{x}}_{i} \end{align}

The internal force can be written in terms of grid velocity by taking a second derivative:

(9) \begin{align}-\frac{\partial e}{\partial \hat{\boldsymbol{x}}_{i}}=-\frac{\partial e}{\partial \boldsymbol{v}_{i}}\frac{\partial \boldsymbol{v}_{i}}{\partial \hat{\boldsymbol{x}}_{i}}=-\frac{1}{{\Delta} t}\frac{\partial e}{\partial \boldsymbol{v}_{i}} \end{align}

In Eq. 8, the intermediate deformation gradient $\hat{\boldsymbol{F}}_{p}$ can be calculated by

(10) \begin{align}\hat{\boldsymbol{F}}_{p}\left(\hat{\boldsymbol{x}}_{i}\right)=\left(\boldsymbol{I}+{\Delta} t\boldsymbol{C}\left(\boldsymbol{v}_{i}\right)\right)\boldsymbol{F}_{p} \end{align}

For implicit time integration, by taking $\boldsymbol{v}_{i}=\boldsymbol{v}_{i}^{n+1}$ and Eq. 6 for $\boldsymbol{h}(\boldsymbol{v}_{i}^{n+1})$ becomes

(11) \begin{align}m_{i}\left(\boldsymbol{v}_{i}^{n+1}-\boldsymbol{v}_{i}^{n}\right)/{\Delta} t-\boldsymbol{f}_{i}^{ext}+1/{\Delta} t\left(\partial e/\partial \boldsymbol{v}_{i}^{n+1}\right) \end{align}

Eq. 11 has only one unknown variable $\boldsymbol{v}_{i}^{n+1}$ . Also, at $\boldsymbol{h}=\mathbf{0}$ it is a large set of nonlinear equations. In order to simplify the process of finding the roots, we perform an integration of Eq. 11 over $\boldsymbol{v}_{i}^{n+1}$ , to construct an objective function for optimization $E(\boldsymbol{v}_{i}^{n+1})$ .

(12) \begin{align}\sum _{i}\frac{1}{2{\Delta} t}m_{i}\left\| \boldsymbol{v}_{i}^{n+1}-\boldsymbol{v}_{i}^{n}-{\Delta} t\frac{\boldsymbol{f}_{i}^{ext}}{m_{i}}\right\| ^{2}+\frac{1}{{\Delta} t}e\left(\boldsymbol{v}_{i}^{n+1}\right) \end{align}

By using automatic differentiation, the derivative $\partial E(\boldsymbol{v}_{i}^{n+1})/\partial \boldsymbol{v}_{i}^{n+1}$ can be easily obtained. In our implementation, we use chain rule to decouple contribution from different grids, which fastens the computation.

3. Numerical implementation for simulation of soft robot

There are three stages to be taken in a time step of dynamic analysis. For the first stage, the physical properties of particles are transferred to the adjacent grids (spatial points at which the shape functions are centered). Mass, momentum, and external forces are transferred from particles to grids by weight factor. For weight $\alpha _{ip}$ , we choose quadratic B-spline function [41] as it takes less time to compute.

(13) \begin{align} N\left(x\right)= \left\{\begin{array}{l@{\quad}l} 3/4-\left| x\right| ^{2}, & 0\leq \left| x\right| \lt 1/2\\ 1/2\left(3/2-\left| x\right| \right)^{2}, & 1/2\leq \left| x\right| \lt 3/2\\ 0, & 3/2\leq |x| \end{array}\right. \end{align}

The external force is applied on grids, indicating that the driven sources can be simulated by analyzing and modeling the mechanics of driven sources and then computing their direction and magnitude. There are many driven sources, for example, pneumatic pressure, hydraulic pressure, and heat driven. We use pneumatic soft robots [24, 25] to demonstrate how to apply driven sources to soft robot simulation. Air drives the soft robots using pressure. The driven pressure is required to simulate the behavior of pneumatic soft robots. For every particle located at the driven loading surface, there exists a force proportional to the input pressure and parallel to the normal direction of local surface. Thus, by estimating the particle’s normal direction and area, the driven pressure can be applied in simulation. Nanson’s relation [26] can be used to compute. The formula is given as follows:

(14) \begin{align}\boldsymbol{n}\mathrm{d}a=J\boldsymbol{F}^{-T}\boldsymbol{N}\mathrm{d}A \end{align}

where $\mathrm{d}a$ and $\mathrm{d}A$ are the area of the particle located at driven loading surface in current and reference configuration, respectively, and $\boldsymbol{n}$ , $\boldsymbol{N}$ are the normal direction pointing to the interior of the material of $\mathrm{d}a$ and $\mathrm{d}A$ , respectively. When generating the geometric model, particles on the driven loading surface are marked, and normal directions $\boldsymbol{N}$ are estimated. The area $A$ is obtained when meshing the geometric model. For any time step after starting the simulation, the forces on particles located at the driven loading surface can be computed.

(15) \begin{align}\boldsymbol{f}_{p}^{ext}=p^{load}\boldsymbol{n}\mathrm{d}a \end{align}

where $p^{load}$ is the magnitude of the input pressure. Then, the forces on particles are interpolated to grids:

(16) \begin{align}\boldsymbol{f}_{i}^{ext}=\sum _{p}\alpha _{ip}p^{load}\boldsymbol{n}\mathrm{d}a \end{align}

By using Eq. 16, we can incorporate the pneumatic pressure into our simulation without having to compute the normal direction and magnitude of local force for each time step separately. For the second stage, the internal force $\boldsymbol{f}_{i}^{{int}^{n}}$ is computed using automatic differentiation, and intermediate grid variables are updated using implicit time integration. To utilize automatic differentiation, JAX [42], a programing tool developed by Google, is used. It supports both forward and backward mode of automatic differentiation. As simulation requires extensive computation, so we used JAX which provides parallel computation on GPU. In the last stage, the intermediate variables on grids are transferred back to update the particles as $\boldsymbol{v}_{p}^{n+1}=\sum _{i}\alpha _{ip}\boldsymbol{v}_{i}^{n+1}$ and $\boldsymbol{x}_{p}^{n+1}=\boldsymbol{x}_{p}^{n}+{\Delta} t\boldsymbol{v}_{p}^{n+1}$ . Algorithm 1 summarizes the implemented numerical simulation approach to clarify shooting interloop (Table 2).

Table 2. The implicit framework with automatic differentiation.

4. Validation and results

4.1. Setup of experiments

Soft finger is the fundamental component of soft robotic systems, as they are responsible for moving or controlling wearables, haptic, and medical devices [43]. We chose soft finger to validate our simulation because it possesses a high degree of flexibility and has been applied to a variety of tasks. We use DragonSkin-30 silicone rubber to construct the soft finger’s body and RTV silicone glue as the adhesive. The material parameters for DragonSkin-30 silicone rubber in Mooney-Rivlin is set as $C_{1}=1.19\mathrm{kPa}$ , $C_{2}=23.028\mathrm{kPa}$ , $1/D_{1}=k/2=1.75\mathrm{GPa}/2$ , where $k$ is the initial bulk modulus. For fluid simulation, its bulk modulus $\lambda$ is set to less than $2.1\mathrm{GPa}$ .

Our geometric model is generated using two steps shown in Figure 4. Firstly, we use Gmsh [44], an opensource software for FEA meshing, to generate tetrahedral mesh. Then, we make the center and volume of each element as the position and volume of the particle, respectively. In order to simulate the pressure acting on the inside chambers of the soft finger, a set of particles on chamber walls are specially marked out and assigned two extra properties. One is the area $\mathrm{d}a$ on the wall surface, and other is the normal direction of area $\boldsymbol{n}$ . These settings make it possible to apply the Nanson’s relation to compute the input pressure.

Figure 4. Particle based geometry model from the meshed geometry model.

Figure 5. Soft finger experiment and simulation results at different input pressures.

Figure 6. Comparison of fingertip Euclidean distance under varying input pressure between our method and CEL simulation.

Figure 7. Comparative analysis of deformation of soft finger in fluid environment between simulation and experimental results.

4.2. Validation of simulation’s accuracy

In our simulation’s object, the free boundary condition contains two cases: specifically nonlinear deformation and nonlinear interactions with nonsolid materials. Therefore, in order to validate the algorithm’s accuracy, we test two cases, respectively, which are soft finger simulation with nonlinear deformation and interaction. The figures and the error rate of experiments show that our framework can model the soft finger effectively to predict its behavior in complex environment.

Experiment I: Validation of simulation with nonlinear deformation

In this experiment, we validated the accuracy of our algorithm based on implicit time integration by using nonlinear deformation in air. We use von-Mises stress to demonstrate the concentration of stresses, as shown in Figure 5 (Visualization using Open3D [45]). The simulation results are consistent with experimental results. Moreover, our computation indicates that the stress concentration occurs simultaneously with severe local deformation at key locations. The results indicate that our method can accurately predict the deformation and stress distribution of soft robots.

When the deformation is nonlinear, the boundary coupling becomes complex. In FEA, handling such complex contact of free boundary requires extra computation. However, our approach solves such a problem autonomously.

Experiment II: Validation of simulation with nonlinear interactions

We choose a fluid environment as the representative of nonlinear interactions, and we compare our results with those obtained using the Coupled Eulerian–Lagrangian (CEL) in Abaqus with C3D10M elements for dynamic explicit analysis. The results are shown in Figure 6 and Figure 7. In Figure 6, the y-axis indicates the Euclidean distance of the fingertip from its original position, and the x-axis is the input pressure. In Figure 7, the tracked point on the soft finger is marked as red dot, and the comparison is made between simulation and experiment. The results show that our approach can recover the deformation effects from the environment of the soft robot. Furthermore, the comparison reveals that our proposed method is more accurate, demonstrating the superiority of our method in simulation accuracy.

In Table 3, we computed the error ratio $\gamma /\beta$ , where $\gamma$ is euclidean distance between simulation and results of experiment in water, while $\beta$ is the euclidean distance between simulation and results of experiment in air. During the initial pressure at 0 (kPa), we achieved an error rate of 0.6%. Similarly, for 5, 10, 15, and 20 (kPa), we achieved error rate of 13.7%, 16.7%, 15.1%, and 8.7%, respectively, which indicates that the computation results are nearly 10 times closer in fluid than in air. This indicates that our model meets the requirement of setting up a virtual environment where the interactions are fully exhibited.

5. Discussion

We proposed a novel method for soft robot simulation using MPM by employing automatic differentiation and implicit integration to address the challenges associated with free boundary interactions. We also designed experiments to validate our approach and performed comparisons to demonstrate the advantages of our proposed method in solving free boundary condition problems.

Moreover, there are two major aspects that will significantly increase the applicability of the proposed approach.

In order to emphasize the advantages of free contact simulation and accurate deformation prediction in complex environments, many possibilities for extending the approach to multiphysical phenomena are simplified, including Rayleigh damping and Coulomb friction (both of which dissipate energy in a nonconservation system), turbulent fluid models (which are effective in facilitating bidirectional interactions between fluids and solids), and magnetic and thermal effects in soft robots. In addition, applying different types of soft robot driven sources to simulation is an important topic.

This approach eliminates the need to explicitly determine contact surfaces within the soft robot itself or with other objects in its environment, so it can be used in exploring the adaptability of different soft robot systems to their working environments, which serves as the backbone for the design of optimal topological structures for soft robots.

Although we use a single example for a single type of soft robot, our approach has the potential to be applied to different kinds of soft robots. In practice, there are several complex soft robots with complex chambers, complex structures (even composed of multiple materials), and more diverse functions (not limited to bending or gripping). However, their simulation methods do not differ in principle from ours. The reason is that these soft robots will deform largely and have complex interactions with the environment, which causes the free boundary coupling and makes the simulation difficult. At the same time, our approach can handle these soft robots’ nonlinear interactions and deformations.

Table 3. Error ratio of simulation.

Many driven sources other than pneumatic pressure need to be modeled. For pneumatic pressure, we find out that the force obeys Nanson’s relation and then integrate it into our simulation. Nevertheless, the forces applied to soft robots are different for other driven sources. As the framework expressed, the external force is applied on grids, which means the driven sources can be simulated by analyzing and modeling the mechanics behind other driven sources, then computing their direction and magnitude.

This article primarily uses pneumatic bending soft robots as an example for simulation and experimentation. In future research, we will conduct more comprehensive simulations and experiments on other types of soft robots to analyze the impact of various parameters, including but not limited to pneumatic pressure, on the performance of soft robots.

6. Conclusion

In this article, we propose an approach for soft robot simulation using MPM, in order to overcome the free boundary interactions problem, particularly nonlinear deformation and interactions, so that a variety of soft robot scenarios (i.e., subsea bio-sampling, in-vivo surgery) of multiple phases and objects can be simulated. In the proposed approach, a vivid virtual environment can be set up easily to accurately predict the behavior of soft robot, which fastens and promotes the design of soft robots. Also, we explored the application of automatic differentiation and implicit integration, which achieve efficient and stable simulation, respectively. The error rate of experiments shows that our framework can simulate the soft robot effectively to predict its behavior in complex environment. In the future, we will extend our work based on soft finger for more complex morphological soft robots.