Finite Element Modeling of Soft Fluidic Actuators: Overview and Recent Developments

Abstract

Many soft robots are composed of soft fluidic actuators that are fabricated from silicone rubbers and use hydraulic or pneumatic actuation. The strong nonlinearities and complex geometries of soft actuators hinder the development of analytical models to describe their motion. Finite element modeling provides an effective solution to this issue and allows the user to predict performance and optimize soft actuator designs. Herein, the literature on a finite element analysis of soft actuators is reviewed. First, the required nonlinear elasticity concepts are introduced with a focus on the relevant models for soft robotics. In particular, the procedure for determining material constants for the hyperelastic models from material testing and curve fitting is explored. Then, a comprehensive review of constitutive model parameters for the most widely used silicone rubbers in the literature is provided. An overview of the procedure is provided for three commercially available software packages (Abaqus, Ansys, and COMSOL). The combination of modeling procedures, material properties, and design guidelines presented in this article can be used as a starting point for soft robotic actuator design.

Introduction

Soft robotics has experienced an exponential growth in publications in the last two decades.^[1] Unlike rigid conventional manipulators,^[2,3] soft robots based on hydrogels,^[4,5] electroactive polymers,^[6,7] and elastomers^[7–9] are physically resilient and can adapt to delicate objects and environments due to their conformal deformation.^[10,11] They also show increased safety and dexterity can be lightweight and used within constrained environments with restricted access.^[12,13] Many soft robots have a biologically inspired design coming from snakes,^[14–17] worms,^[18–20] fishes,^[21–24] manta rays,^[25,26] and tentacles.^[27–29] The scope of applications includes minimally invasive surgery,^[30,31] rehabilitation,^[32,33] elderly assistance,^[34] safe human–robot interaction,^[35,36] and handling of fragile materials.^[37,38] Important features of soft robotics design, fabrication, modeling, and control are covered in the soft robotics toolkit.^[39,40]

The building blocks of soft robots are the soft actuators. The most popular category of soft actuator is the soft fluidic actuator (SFA), where actuation is achieved using hydraulics or pneumatics.^[8,41] These actuators are usually fabricated with silicone rubbers following a 3D molding process,^[42] although directly 3D printing the soft actuators is also possible.^[43,44] Silicone rubber is a highly flexible/extensible elastomer with high-temperature resistance, low-temperature flexibility, and good biocompatibility.^[45] Elastomers can withstand very large strains over 500% with no permanent deformation or fracture.^[46] For relatively small strains, simple linear stress–strain relationships can be used, and two of the following parameters can be used to describe the elastic properties: bulk compressibility, shear modulus, tensile modulus (Young's modulus of elasticity), or Poisson's ratio.^[45] For large deformations, nonlinear solid mechanic models using hyperelasticity should be considered.^{[8,32,47–50]} Due to the strong nonlinearities in SFAs and their complex geometries, analytical modeling is challenging.^[51] A brief review of the analytical methods for modeling of soft robotic structures is provided in the following.

Analytical Modeling of Soft Actuators

The majority of soft/continuum robots with bending motion can be approximated as a series of mutually tangent constant curvature sections, i.e., piecewise constant curvature.^[52] This is a result of the fact that the internal potential energy is uniformly distributed along each section for pressure-driven robots.^[53] This approach has also been validated using Hamilton's principle by Gravagne et al.^[54] As discussed by Webster and Jones,^[52] the kinematics of continuum robots can be separated into robot-specific and robot-independent components in this approach. The robot specific mapping transforms the input pressures P or actuator space q to the configuration space $\kappa ,\varphi ,l$, and the robot-independent mapping goes from the configuration space to the task space x. The actuator space contains the length of tubes or bellows. The configuration space consists of the curvature κ, the angle of the plane containing the arc ϕ (also called azimuthal angle or angle of curvature) and arc length l (or $s\in [\mathrm{0,}l]$). Alternatively, κ or l can be replaced by the bending angle θ or the radius of curvature r through the relations $\theta =\kappa s$ and $\kappa =1/r$.

The robot-independent mapping can be obtained with arc geometry,^[55,56] Denavit–Hartenberg parameters,^[3,57,58] differential geometry (Serret–Frenet frame),^[55,59,60] integral representation,^[59,61,62] or exponential coordinates.^[63] For the robot-specific transformation, one must account for the transformation from input pressures to actuator lengths or input pressures directly into configuration space. In the latter case, Suzumori et al.^[64,65] obtained the robot-specific transformation by linear analysis based on the theory of infinitesimal elastic deformation and the constant curvature assumption. Under the constant curvature assumption, Drotman et al.^[66] developed an analytical expression for the blocked force of a parallel bellows actuator using Castigliano's method.

The piecewise constant curvature approach is practical when inertia effects are negligible due to low speed and mass.^[28] To account for the mass of the actuator and external loading, models using the Euler–Bernoulli equation^[8,54,67] and the theory of Cosserat rods^[68–70] have been considered.

The Euler–Bernoulli beam theory is used to model bending actuators as ideal beams loaded with a pure moment at its edges, which originates from the axial component of the loading pressure and is counteracted by the bending stiffness EI of the beam, where E is Young's modulus, and I is the second moment of area.^[8] The Euler–Bernoulli equation has been applied to the modeling of cylindrical actuators with eccentric void,^[67,71] multimaterial asymmetry,^[72] and omnidirectional actuators.^[65] However, this method is not applicable to large bending deformations or corrugated cross sections. In addition, Young's modulus does not describe the complex stress–strain behavior of most materials used in soft robotics.^[8]

Many soft robots have a slender structure with one dimension much larger than the other two; hence, they can be modeled using the theory of Cosserat rods.^[68–70] In the study given by Trivedi et al.,^[73] a geometrically exact model is presented that accounts for the large deformations and loading using concepts from nonlinear elasticity and Cosserat rod theory for the manipulator dynamics. This modeling approach was shown to be ten times more accurate than the constant curvature model for the OctArm manipulator when gravitational loading is considered.

An alternative is to use an empirical approach with energy-based inspiration, where the dynamic behavior of a bending actuator can be approximated by a lumped second-order dynamic equation.^[14,74] The model parameters can be determined by least-squares curve fitting^[74,75] or system identification with a periodic input signal.^[76,77] Furthermore, a purely data-driven approach can be used to predict the bending angle of soft actuators.^[78]

Finite Element Method in Soft Robotics

Finite element modeling (FEM) provides an effective solution to handling nonlinearities in soft robotics and avoids the need for an explicit analytical model. Furthermore, the FEM is not limited to the slender structures required in the theory of Cosserat rods or the modeling of bending motion using the piecewise constant curvature assumption or Euler–Bernoulli beam theory. FEM has found many applications in soft robotics due to the following reasons:^[51,79–81] 1) FEM can cope with the large deformations and material nonlinearities during the pressurization of SFAs. 2) FEM can be used to predict performance and evaluate the capabilities and limitations of soft actuator designs under various inputs. This rapid and efficient design framework reduces cost and development time, because the fabrication of SFAs is very time-consuming. 3) FEM can improve our understanding of the stress concentration and strain distribution in soft actuators. This feature leads to a better understanding of the influence of local strain on global actuator performance and can be used, for example, to determine potential locations of fatigue. 4) FEM can handle contact nonlinearities associated with surfaces that come into contact upon deformation.

Contributions and Outline of the Article

Recent review papers have covered different aspects of soft robotics design,^[8,41] fabrication,^[9,82] biological inspiration,^[83,84] locomotion,^[85] actuation methods,^[86,87] sensing,^[88,89] kinematic modeling,^[52] control,^[90,91] stiffening techniques,^[92] and biomedical applications.^[93,94] Despite the number of articles which use FEM, there is very little comparison of methods or guidelines for successful modeling. This means that each new study often requires significant trial and error to identify a compatible modeling and simulation procedure. Furthermore, these difficulties are compounded by the lack of a comprehensive database for nonlinear material properties of common silicone rubbers.

To address the first issue, different strategies and recommendations for analysis and design of SFAs with FEM are presented using Abaqus (Daussault Systèmes SE), Ansys (Ansys, Inc.), and COMSOL Multiphysics (COMSOL Inc.). These recommendations draw from procedures described in the literature and our group's experience with FEM.

A step toward resolving the second issue is the July 2020 article of Marechal et al.,^[95] where the authors present the mechanical characterization and the resulting material properties for a wide set of commercially available silicone rubbers. Here, the material constants presented in the literature are compiled and serve as a starting point for soft roboticists in the modeling and design of soft robots.

This article first reviews the different categories of SFAs in Section 2. The continuum mechanics background and hyperelastic models that underlie the modeling of soft actuators are introduced in Section 3. The strategies for characterizing the elastomeric materials used in the fabrication of soft actuators are discussed in Section 4. Section 5 provides an extensive review of the hyperelastic model constants in the literature for the most popular silicone rubbers within the soft robotics community. Using these model parameters, the FEM procedure can then be introduced. Section 6 presents an overview of the FEM procedure using commercially available software. In Section 7, a comprehensive set of design guidelines for soft actuators is presented based on results drawn from the literature and experience. Finally, the challenges and opportunities of FEM for SFAs are discussed in Section 8.

Soft Fluidic Actuators

SFAs can be classified according to their motion into four categories: extending, contracting, bending, and twisting.^[8,41] These actuators have specifically engineered anisotropic flexible structures to achieve each of the four motions.^[8,96] For bending actuators, the anisotropic behavior can be achieved through the use of multi-material, corrugated membrane (multi-chamber) or eccentric void asymmetries (Figure 1a).^[8] In the first case, the actuator has different rubber compositions or a strain limiting layer.^[97] In the second case, folds (fins) on one side of the actuator expand under pressure. In the third case, the inflatable void of the actuator is placed off-center, leading to different layer thicknesses in the actuator.^[71] Combinations of the three scenarios are also possible; Sun et al.,^[97] for example, evaluated multi-chambered actuators with fabric as a strain limiting layer.

FIGURE: SFAs designs. a) Bending actuators. Reproduced with permission.^[8] Copyright 2017, John Wiley & Sons, Inc. b) Omnidirectional actuators. Reproduced with permission.^[80] Copyright 2014, Mary Ann Liebert, Inc. c) Multi-chambered actuators. Reproduced with permission.^[47] Copyright 2014, John Wiley & Sons, Inc. d) Fiber-reinforced actuators. Reproduced with permission.^[102] Copyright 2017, PNAS. e) Bidirectional actuator. Reproduced with permission.^[75] Copyright 2017, IEEE. f) One degree of freedom revolute joint. Reproduced under the terms of the CC-BY 4.0 License.^[109] Copyright 2018. The authors, published by PLOS. g) Parallel bellows actuator. Reproduced with permission.^[49] Copyright 2017, IEEE.

The simplest design for SFAs consists of a single chamber. However, unless fibers are wrapped around the chamber, the actuator behaves like a balloon with high radial expansion.^[48,98] One of the most investigated designs in the literature is the multi-chambered actuator. A popular example is the PneuNets,^[32,47,99] which consists of an extensible top layer and an inextensible but flexible bottom layer. Slow PneuNets (sPNs) use a block of silicone rubber with embedded air chambers,^[97,99] whereas the fast PneuNets (fPNs) developed by Mosadegh et al.^[47] contain gaps between the inside walls of each chamber (Figure 1c). Other designs for multi-chambered actuators include the use of trapezoidal^[48] and triangular^[100] bellows.

Another common design is the fiber-reinforced actuator,^[81,101,102] which draws inspiration from McKibben actuators, also referred to as pneumatic muscles. For these actuators, fibers can be arranged along the length of the actuator to achieve different motions (Figure 1d). Double wrapping of fiber is used in extending, contracting, and bending actuators, whereas single wrapping of fiber can be used to introduce twisting. Moreover, an actuator can be fabricated with different motion types in different segments, thus combining multiple behaviors. Fiber-reinforced actuators are also referred to as PneuFlex actuators.^[38,103,104]

Soft actuators with multiple degrees of freedom have also been developed. Bidirectional actuators^{[14–16,105]} can be achieved using soft actuators with two chambers or by joining two bending actuators via the bottom layer. Omnidirectional actuators were proposed by Suzumori et al.^[64,65] and further explored by Yan et al.,^[106] Elsayed et al.,^[80] and Khan and Li.^[76] The omnidirectional actuator has with three or more internal chambers (Figure 1b). When three chambers are actuated with the same pressure, the actuator stretches. In contrast, when only one or two chambers are actuated, the actuator bends in the direction opposite to the pressurized chambers. These actuators have three degrees of freedom, which are pitch, yaw, and stretch. Actuators with three degrees of freedom can also be fabricated using three (fiber-reinforced) extending actuators between two plates forming an equilateral triangle^[107,108] (Figure 1g) in a design inspired by the parallel bellows actuators in pneumatic continuum robots.^[27,35] Alternatively, cylindrical fiber-reinforced actuators can be combined in smart structures to achieve bending (Figure 1f) as demonstrated in the one degree of freedom revolute joint mechanism in the study given by Skorina et al.^[74,109]

Background

Hyperelastic models are considered that assume silicone rubber is isotropic and incompressible. Furthermore, all inelastic phenomena such as viscoelasticity and stress-softening are neglected.

Fundamentals of Continuum Mechanics

The deformation of a solid can be described by the relationship between the spatial coordinate frame x (current configuration) and the material coordinate frame X (original configuration) as[Image Omitted. See PDF]where u is the displacement vector.^[110–112]

The deformation gradient is given by[Image Omitted. See PDF]

F is the Jacobian matrix of the transformation from X to x; hence, $J=\text{det}(\mathbf{F})$ is the local volume scale factor. For an incompressible material, $J=1$.^{[110,111,113]}

The polar decomposition theorem states that any second-order tensor can be decomposed into a product of pure rotation (R) and a symmetric deformation tensor (U the right stretch tensor or V the left stretch tensor)[Image Omitted. See PDF]

Thus[Image Omitted. See PDF][Image Omitted. See PDF]where C is the right Cauchy–Green deformation tensor, and B is the left Cauchy–Green deformation tensor.^[110–112] C and B admit the same three principal invariants $I_1$, $I_2$, and $I_3$ given by^{[110,111,113]}[Image Omitted. See PDF]

The stretch ratios ${\lambda }_i$ represent the deformation of a differential cubic volume element along the principal axes of a Cartesian coordinate system.^[112,114] They are given by the ratio of deformed length ($l_i$) to undeformed length ($L_i$), that is[Image Omitted. See PDF]

The three eigenvalues of U (or V), i.e., ${\lambda }_1$, ${\lambda }_2$, and ${\lambda }_3$, are called principal stretches. The corresponding eigenvectors of U give three orthogonal directions in the material frame. Using the principal stretches, the principal invariants reduce to[Image Omitted. See PDF]

A hyperelastic material model relies on the definition of a strain energy function W (or strain energy density), which is the amount of energy stored elastically in a unit volume of material under the state of stretch specified by ${\lambda }_1$, ${\lambda }_2$, and ${\lambda }_3$.^{[110,111,113]} For an incompressible material, $I_3=1$ and $W=W(I_1,I_2)$. Material models using these functions are discussed in Section 3.2.

For an isotropic incompressible hyperelastic material, stress–stretch relations are obtained from the strain energy function by virtual work considerations[Image Omitted. See PDF]where σ is the principal Cauchy stress, and p is the hydrostatic pressure (indeterminate Lagrange multiplier), which can be determined from the equilibrium equations and boundary conditions. Note that p does not represent the actual pressure on the chosen surfaces.^[45,110,115] For a uniaxial tensile experiment, ${\lambda }_1=\lambda$, ${\lambda }_2={\lambda }_3=1/\sqrt{\lambda }$, and those boundary conditions are ${\sigma }_2={\sigma }_3=0$. Then, the stress as a function of the stretch is^{[45,110,113,114]}[Image Omitted. See PDF]

Hyperelastic Material Models

The main hyperelastic models used in soft robotics are reviewed in this section. A typical stress–strain curve for hyperelastic materials is shown in Figure 2, alongside the fitting of hyperelastic models (2)–(5). The curve shows the highly extensible and elastic behavior of hyperelastic materials where large recoverable strains are produced at low stress levels.^[45,116] Note that the Cauchy strain e is defined as the ratio of total deformation to the initial dimension of the body; hence, $e=\lambda -1$. The models described here assume incompressible rubber, i.e., $I_3=1$ and $D_i=0$, where $D_i=0$ are material constants describing the bulk compressibility.^[115] For other constitutive models, the reader is referred to the studies by Marckmann and Verron,^[113] Martins et al.,^[115] Aidy Ali and Hosseini,^[46] and Mihai and Goriely.^[117]

FIGURE: Typical stress–strain curve for a hyperelastic material and fitting of different hyperelastic models with uniaxial tensile testing data. Simple models such as Neo-Hookean and Mooney–Rivlin are only accurate for small strains. In contrast, the Yeoh and Ogden models have the ability to match experimental data at small and large strain values. Reproduced under the terms of the CC-BY 4.0 License.^[118] Copyright 2015. The authors, published by ABM, ABC, and ABPol.

Polynomial Model (Generalized Rivlin Model)

[Image Omitted. See PDF]This model is often truncated in terms of the second and third order.

Mooney–Rivlin Model

[Image Omitted. See PDF]This model is used for moderate deformation, i.e., lower than 200%.

Ogden Model

[Image Omitted. See PDF]Six parameter model (N = 3) is the most used for large strain problems, i.e., at or above 400%.

Experimental results show that $\frac{\partial W}{\partial I_1}$ is approximately constant, and that $\frac{\partial W}{\partial I_2}$ decreases with the amount of strain. Therefore, one simplification is to ignore terms in $I_2$.^[45]

Yeoh Model

[Image Omitted. See PDF]Also, one of the most used models for large strain problems, i.e., at or above 400%.

Neo-Hookean Model

[Image Omitted. See PDF]where μ is the shear modulus. Good agreement for small strains, i.e., lower than 50%.

Material Testing and Curve Fitting Material Testing

The coefficients in the hyperelastic models should be determined using uniaxial, biaxial, and shear test data.^[46,113] However, most studies have simply used uniaxial testing, because this is easy to perform with readily available material testing equipment.^[118] An outline for the procedure in the fabrication and uniaxial testing of silicone rubber samples is shown in Figure 3.

FIGURE: Fabrication and testing of elastomer sample: 1) silicone rubber is mixed with recommended proportions from manufacturer, 2) mixture is degassed and 3) poured into 3D printed molds, 4) samples are cured for recommended time, 5) dimensions are measured, and 6) quality of specimen is verified. Finally, uniaxial tensile testing or other material characterization procedures are performed, and non-contact strain measurements are obtained using video or laser extensometer. Reproduced with permission.^[95] Copyright 2020, Mary Ann Liebert, Inc.

According to Marckmann and Verron,^[113] Shahzad et al.,^[118] Rackl,^[114] and Meier et al.,^[119] a unique experiment is not sufficient to characterize the behavior of elastomers, especially if the elastomer experiences multiaxial stress states, which is the particular case with the pressurization of soft actuators. Under the incompressibility assumption, all deformations are determined by two independent stretch ratios. Consequently, a series of biaxial tests is suggested to completely characterize the multiaxial stress states in the hyperelastic constitutive models.^[113]

For models (4) and (5), Aidy Ali and Hosseini^[46] and Gent^[45] argue that uniaxial tensile testing is enough to characterize the material response. In particular, the Yeoh model describes the elastic behavior with reasonable success over quite large ranges of strain and is able to predict stress–stretch behavior in different deformation modes from data gained with simple uniaxial testing. On the other hand, attention should be paid not to use the Ogden model with limited testing data, e.g., just uniaxial tension.^[118]

Testing of elastomers has not yet been clearly defined by international standards.^[120] For uniaxial tensile testing of elastomers, Marechal et al.^[95] have recommended the use of the ASTM D142 standard.^[121] In addition, the following guidelines are suggested by Miller:^[120] 1) In uniaxial tensile tests, the specimen must be at least ten times longer than the width or thickness. 2) In pure shear tests, the specimen must be at least ten times wider than the length in the stretching direction. 3) Biaxial tensile tests are recommended for a full description of the deformation modes of the elastomer. 4) Use non-contacting strain measuring devices such as video or laser extensometer and measure strains away from the clamp. 5) Due to stress softening, consider 3–20 repetitions of the tests to obtain stable levels. 6) Tests should be performed at strain levels relevant for the application.

Curve Fitting

In the determination of the hyperelastic material constants, a particular expression is obtained, substituting each of the models discussed in the previous section into Equation (10). The respective stress–stretch equations are shown in Table 1.^[95,114,115] These relations are derived as a function of stress σ and stretch λ, while tensile testing data are force F and change in length $\delta L$, respectively. Therefore[Image Omitted. See PDF]where A is the initial cross-sectional area of the reduced section of the sample, and $L_0$ is the undeformed length.

Table 1 Stress–stretch equations for curve fitting with uniaxial tensile testing data

Model	Stress–stretch equation
Mooney–Rivlin	$\sigma =2({\lambda }^2-{\lambda }^{-1})(C_1+C_2{\lambda }^{-1})$
Ogden	$\sigma ={\sum }_{p=1}^n{\mu }_p({\lambda }^{{\alpha }_p-1}-{\lambda }^{-({\alpha }_p/2+1)})$
Yeoh	$\sigma =2({\lambda }^2-{\lambda }^{-1}){\sum }_{i=1}^{n}i{C}_i{({\lambda }^2+2{\lambda }^{-1}-3)}^{i-1}$
Neo-Hookean	$\sigma =2({\lambda }^2-{\lambda }^{-1})C_1$

To determine the coefficients $C_i$, ${\alpha }_p$, and ${\mu }_p$, a nonlinear least squares optimization method is used to minimize the error with respect to the parameters of the model. This can be achieved using, for example, Scilab,^[114] MATLAB, or Python^[95] and also within the standard FEM software discussed here. It is important to note that the models can lead to error if they are used out of the deformation range in which the parameters were identified.^[113] Furthermore, the calculation of the coefficients in the Ogden model is more computationally intensive than in the other models and can be more accurate in fitting experimental results when data are available from multiple experimental tests. Therefore, models with few parameters are preferred for the purpose of computational efficiency.^[46]

The nonlinear least squares method involves minimizing the sum of squared errors between the theoretical and experimental stresses as a function of stretch. In MATLAB, for example, this is achieved using the function fit with a nonlinear least squares method in fitoptions and using one of the equations in Table 1 in the definition of fittype. As another option, Marechal et al.^[95] have provided a set of Python functions for the implementation of the fitting algorithm.

Hyperelastic Models for Soft Actuator Materials

The most widely used silicone rubbers within the soft robotics community are Ecoflex (Shore Hardness 00-10 to 00-50),^[122] DragonSkin (Shore Hardness 10–30 A),^[123] Smooth-Sil (Shore Hardness 36–60 A),^[124] and Elastosil M4601 (Shore Hardness 28A).^[125] A hardness scale is provided in Figure 4, and uniaxial tensile testing results are shown in Figure 5. Ecoflex is softer than the other elastomers, which results in lower pressure levels being required for desired bending or extension strokes of the soft actuators. Many hyperelastic models for these elastomers have been proposed in the literature, as reviewed in Table 2 (N/A stands for not available). The parameters for additional materials used in the fabrication of soft actuators are shown in Table 3, where E denotes Young's modulus, ρ is the density, and ν is Poisson's ratio. The reader is also referred to Table 3 in the study by Marechal et al.^[95] for further hyperelastic material constants.

FIGURE: Shore hardness scale for silicone rubbers in soft robotics.

FIGURE: Experimental uniaxial tensile stress–stretch pull to failure responses for a selection of silicone rubbers. The mean curves and 95% confidence bands are presented from five samples for each material. Reproduced with permission.^[95] Copyright 2020, Mary Ann Liebert, Inc.

Table 2 Constitutive model parameters for silicone rubbers in soft robotics

[Table omitted. See PDF]

Table 3 Additional materials used in the modeling of soft actuators with FEM

[Table omitted. See PDF]

Two approaches can be followed to model the strain limiting layer. One alternative is the inclusion of a thin surface layer within the computer-aided design (CAD) models and then assigning the corresponding material properties to such layer. Alternatively, to improve convergence and computational efficiency, the behavior of the strain limiting layer can be incorporated into the silicone rubber parameters using a constant that considers the increased stiffness of the combined layer, as performed in the study by Polygerinos et al.;^[81] see $C_{\text{combined}}$ in Table 3.

As discussed in the study by Moseley et al.^[79] and Martins et al.,^[115] the Neo-Hookean and Mooney–Rivlin models are based on linear approximations of the strain invariants and are accurate at strains below 50% and 200%, respectively. The Ogden, Yeoh, or higher order polynomial models are recommended for higher deformations.^[113,126] In the FEM of SFAs, the Mooney–Rivlin and Neo-Hookean models can be safely used for silicone rubbers of higher shore hardness, such as Elastosil M4601 and Smooth-Sil, and perhaps DragonSkin. However, for softer silicone rubbers such as Ecoflex, the Yeoh and Ogden models are preferred, as suggested in Table 2. In addition, as shown in Figure 6, fiber-reinforced or multi-chambered actuators show reduced volume change and ballooning, i.e., reduced radial expansion; hence, lower order models can be used for a wider range of silicone rubbers.

FIGURE: Hyperelastic materials and models used in SFAs. The Neo-Hookean and Mooney–Rivlin models are limited to low-strain regimes, whereas the Ogden and Yeoh models are also accurate at high strains. SFAs composed of softer materials, such as Ecoflex and DragonSkin, are expected to experience large levels of deformation; consequently, the Yeoh and Ogden models are recommended. However, fiber-reinforced and multi-chambered actuators show reduced ballooning, which has permitted a number of authors to utilize lower order models such as Neo-Hookean or Mooney–Rivlin in the FEM of these actuators.^[81,127,147]

Overview of FEM Procedure with Commercially Available Software General Procedure for FEM of SFAs

The three prevalent software packages for FEM of SFAs are Abaqus, Ansys, and COMSOL, as outlined in Table 2 and other works discussed in this article. The wide adoption of Abaqus for modeling soft actuators can be associated with the highly cited works of Connolly et al.,^[127] Polygerinos et al.,^[81] and Mosadegh et al.^[47] and the availability of FEM tutorials using Abaqus in the soft robotics toolkit.^[128,129] An overview of the procedure is shown in Figure 7.

FIGURE: Overview of the FEM procedure for SFAs: 1) drawing the soft actuator geometry in CAD software, 2) assignment of material properties (refer to Table 2 and 3 or Marechal et al.^[95]), 3) meshing, 4) modeling of pressurization and mechanical fixture using boundary conditions and loads, and 5) analysis of results.

All three software packages offer capabilities for drawing, but models can also be imported from, for example: Solidworks (Daussault Systèmes SE), Autodesk Inventor (Autodesk Inc.), Fusion 360 (Autodesk Inc.), and Creo (PTC Inc.). While drawing the SFA geometry, the user should consider avoiding small features that require a fine mesh to provide adequate convergence. Air flow inlets or connections between chambers, for example, may be omitted without effecting the simulation accuracy.

For multi-material actuators, a split tool in the CAD software can be used to generate two bodies that are assigned different materials within the FEM software. As an alternative, an assembly file can be imported. Note that, for bending actuators, an extra layer of the same SFA material should be added before using the split tool to select only the bottom cavity rather than the whole bottom layer.

Independent descriptions of the procedure for FEM of SFAs using Abaqus 6.11, Ansys 2020, and COMSOL 5.5 are presented in the Appendix. Recommended meshing strategies based on each software's capabilities are discussed, and further recommendations for convergence are presented in Section 6.2. It is important to note that reliable results should be relatively mesh-independent. Specific settings for boundary conditions, strain limiting layer, contact boundaries, and fiber reinforcements in each of three software packages are introduced in the Appendix. In particular, contact boundaries are attached to external walls in multi-chambered actuators to simulate the contact of these walls under pressurization.

Ansys 2020 and COMSOL 5.5 simulation files for a variety of soft actuator designs are provided in the Supporting Information. More specifically, the following files are included: 1) Semi-circular bending actuator with two materials in Ansys (Supporting Information S1). 2) Fiber-reinforced bending actuator in Ansys (Supporting Information S2). 3) Multi-chambered bending actuator (fPN) in COMSOL (Supporting Information S3). 4) Multi-chambered manta fin prototype (sPN) in COMSOL (Supporting Information S4).

Note that Abaqus files can be found in the soft robotics toolkit.^[128–130]

Recommendations for Improved Convergence

1) According to Tawk and Alici,^[51] as SFAs can undergo large deformations, a very fine mesh is not recommended, and a relatively coarse mesh instead is desired for convergence. However, a coarse mesh might result in lower deformations for the same pressure value. In the authors’ experience, while coarse meshes have facilitated convergence, they have failed to capture strong ballooning effects during the pressurization of soft actuators. 2) As a rule of thumb, reduced mesh sizes are only required in areas of large deformations and stresses. For soft actuators of typical sizes, i.e., 50–200 mm length and 10–25 mm width/diameter,^{[47,79–81,102,127]} a mesh size of 1–3 mm is a good starting point. 3) A 3D problem is always more complicated to analyze than a 2D one, especially in nonlinear FEM. Whenever possible, take advantage of symmetry in the model to reduce the total number of elements in the model.^[45] 4) If the mesh distorts badly during the simulation, the mesh density needs to be changed between load steps.^[45] 5) During meshing, it is preferable to use a mesh with quadratic order elements and a hex dominant configuration, because hexahedron elements are less distorted when capturing curved geometries. 6) A reduction in the size of time steps can also improve convergence. As another option, a larger number of substeps can be used. This is an important parameter to consider, because it prevents the load from being applied in one step, which usually results in convergence issues.^[51] 7) Increase the number of iterations for simulations that are close to converging. In Ansys, for example, this can be achieved using the command NEQIT,50 to increase to 50 iterations (standard of 26).

Validation of FEM Results

Once the design is complete and relevant geometries and material properties have been optimized, the final geometry can be fabricated for validation of the results. This can be achieved by directly 3D printing the geometry^[51,131] or 3D printing molds for the desired shape.^[9,82]

Experimental characterization of SFAs can be performed using compressed air systems,^[81,132,133] syringe pumps,^[134–136] fluidic drive cylinders,^[137,138] or simply by manually actuating a syringe.^[139] Bending, extension, or twisting levels are post-processed from video recordings and/or obtained from motion trackers. These measurements are then combined with the corresponding pressure values for comparison with FEM results.

FEM of SFAs: Design Guidelines

In this section, a collection of design guidelines for the different types of SFAs is presented. Most of these results are drawn from FEM studies; nonetheless, relevant concepts from experimental characterization are also included. While these conclusions have been organized into four sections, note that the conclusions for single chamber actuators also apply for fiber-reinforced actuators. In addition, bidirectional and omnidirectional actuators fabricated using multiple fiber-reinforced actuators inherit the characteristics of single fiber-reinforced actuators.

In the design of SFAs, two universal parameters have a strong influence in the pressurization of soft actuators: wall thickness and silicone rubber hardness. Reduced wall thickness and silicone hardness lead to higher deformation at lower pressure values. For fiber-reinforced actuators, the wrapping density and angle have high impact on the performance of the soft actuator. For omnidirectional actuators, increased void area gives rise to higher bending. For multi-chambered actuators, thinner internal chamber walls, higher number of chambers, and the corrugated design from fPNs result in improved strokes. Further design guidelines are presented in the following.

Single Chamber Actuators

1) Silicone rubbers with lower durometer (hardness) exhibit very high levels of bending at low pressures and lower von-Mises stresses (Figure 8a).^[80,101,107] The low pressure level is a desired feature for safe application to biomedical devices. The disadvantage, however, is lower force at the tip, because they support lower pressures. 2) Bending is maximized when one of the layers is two to three times more rigid or thick than the other.^[8,9,71] 3) SFAs with lower wall thickness show higher bending angles and extensions (Figure 9a).^[81,140] 4) Actuation in two directions can be obtained with the actuator deflecting in opposite directions below or above a threshold.^[71] 5) Optimal force is obtained when the ratio of length to width of the inflatable void is ≈10.^[8] 6) Semi-circular actuators are optimal for bending, and cylindrical actuators are optimal for extension.^[81,141,142] 7) For semi-circular actuators, higher bending angles are obtained with larger actuator radius and larger actuator length.^[81] 8) Cylindrical actuators have a more linear response and show reduced wear and tear.^[142]

FIGURE: FEM of SFAs. a) Omnidirectional actuators. Reproduced with permission.^[80] Copyright 2014, Mary Ann Liebert, Inc. b) Fiber-reinforced bending actuators. Reproduced with permission.^[81] Copyright 2015, IEEE. c) Multi-chambered actuators, sPN is the slow PneuNet, and fPN is the fast PneuNet. Reproduced with permission.^[47] Copyright 2014, IEEE. Three important results are outlined: 1) silicone rubbers with lower hardness exhibit higher bending at low pressures and lower von Mises stresses,^[80,101,107] 2) fiber reinforcements limit the radial extension and enhance bending of SFAs,^[81,140] and 3) fPNs require lower pressures and volumes of input fluid for actuation in comparison with sPNs.^[47]

FIGURE: FEM of fiber-reinforced cylindrical SFA. a) Higher internal radius r implies lower wall thickness and, consequently, better axial extension. b) Fiber reinforcements minimize radial expansion and enhance the axial extension performance of the SFA. Reproduced under the terms of the CC-BY 4.0 License.^[140] Copyright 2019. The authors, published by IEEE.

Fiber-Reinforced Actuators

1) Actuators with no fiber wrapping or low fiber wrapping density exhibit a ballooning effect, i.e., high radial expansion.^{[81,101,140,143]} 2) Fiber-reinforced actuators show enhanced extension or bending strokes, require lower amounts of input flow, and minimize the energy lost in radial expansion of the rubber (Figure 9b).^[140,141] 3) If the fiber density is too large, the actuator requires higher input air pressure.^[140] 4) Dense reinforcements improve linearity, reliability, and durability of SFAs.^[140,142] 5) A higher density of wrapping results in higher bending levels^[141] and reduces the nonlinear behavior resulting from the ballooning effect. This leads to reduced wear of the flexible material.^[142] 6) With single fiber wrapping, maximum axial extension occurs for a fiber angle at 0°, maximum radial expansion and no axial extension for 90°, and maximum twist is obtained for fiber angles around 30°.^[127] 7) Symmetric double fiber wrapping results in no twisting.^[127] 8) For McKibben actuators, actuators with braided angle greater than 54.7° extend, and otherwise contract, when pressurized. Fiber-reinforced SFAs with greater difference between braided angles on either side present higher bending angles and force at the same pressure level.^[144] Combining the fiber angles of 70° and 35° in a single actuator produces the highest bending angle.^[145]

Bidirectional and Omnidirectional Actuators

1) For omnidirectional actuators, higher bending is achieved, in order of importance, with lower wall thickness, larger length, larger chamber diameter, and lower central diameter.^[106,146] 2) For three-chambed omnidirectional actuators, bending ability of a triangular cross section is superior to that of a circular shape.^[146] 3) For omnidirectional actuators, chambers with semi-circular cross section have the least amount of ballooning, and chambers with a ring-sector cross section show the highest bending.^[80] 4) Constrained parallel bellows actuators show higher bending angles in comparison with their unconstrained counterpart. In addition, lower thickness of bellows and larger inner chamber diameter result in higher bending and blocked forces.^[66]

Multi-Chambered Actuators

1) Multi-chambered actuators have a more linear response compared with single inflatable void actuators, because they offer resistance to radial extension.^[47,48] 2) sPNs require approximately three times higher pressure and eight times higher volume to fully bend compared with fPNs (Figure 8c). fPNs also show improved speed and force by the factors of 25 and 1.4, respectively.^[47] 3) The most significant factors that influence the bending angle are in order: bottom layer thickness, wall thickness, and gap size. In particular, smaller gaps result in higher bending but might cause damage to the channels.^[147] 4) For rapid actuation and low radial expansion, the internal walls should be thinner and with larger surface area than the other exterior walls.^[32,47] 5) For a fixed length, more chambers enable greater bending at lower pressures (Figure 10).^[47,48,148] 6) For a fixed length, thicker chamber walls result in lower bending and lower output force/torque.^{[32,47,97,148]} 7) Increased chamber height leads to increased force output.^[32] 8) Free bottom actuators have the outer sides of the actuator bonded to the bottom layer, whereas the intervals between chambers are free. They show ≈20% higher bending and 40% higher force compared with conventional multi-chambered actuators.^[148] 9) Helical multi-chambered actuators have higher mechanical blocking force than the typical bending actuator. Increased chamber angle results in lower bending and higher twisting. In addition, the length of the helical multi-chambered actuator does not influence the winding radius and pitch length but only the number of loops that are created.^[23]

FIGURE: Deflection of multi-chambered actuators at 55 kPa. For a fixed length, more chambers enable greater bending at lower pressures and result in a more linear response to applied pressure. Red and blue indicate 200 kPa and 0 Pa von-Mises stresses, respectively. Reproduced with permission.^[48] Copyright 2015, Elsevier.

Challenges and Opportunities in the FEM of SFAs Challenges and Limitations

This article surveys the use of FEM in modeling SFAs. Using FEM, the soft roboticist can evaluate the effect of material properties and geometry without time-consuming fabrication. Nevertheless, this approach has several limitations that support the need for experimental characterization to validate the FEM results.

First, it is clear from Table 2 that some authors have incorrectly used material constants from a different type of silicone rubber in their simulations. In addition, significant variations in the material coefficients can be seen even within the same hyperelastic model for a specific type of silicone rubber. As discussed in Section 5, the material model has substantial influence in the stroke of the soft actuator. Accordingly, the large differences for material constants presented in the literature can hinder the evaluation of realistic bending, extension, and twisting levels. Moreover, simplified models, such as Mooney–Rivlin and Neo-Hookean, are not recommended for actuators with large deformations due to ballooning. Moseley et al.^[79] also suggest that the hyperelastic model needs to be calibrated across a large range of realistic strains appropriate for the particular application.

Second, care must be taken during meshing of soft robotic structures due to the strong nonlinearities and high deformations they exhibit. Recommendations for meshing have been provided, and it is important to emphasize that different meshes may result in different bending, extension, or twisting levels that are not representative of the motion for the physical device. Note that a coarse mesh may facilitate convergence, but it may also underestimate the actuator's deformation.

Third, the majority of studies use quasi-static simulations in the FEM of soft actuators. However, dynamic effects might need to be included in the simulations at high pressures or for fast actuation, wherein the quasi-static assumption does not hold, and vibrations can be observed.^[76] In this context, transient structural studies are required. Alternatively, to improve the convergence of the model at high pressures, a small amount of Rayleigh damping can be added, which keeps kinetic effects to a minimum and ensures quasi-static conditions.^[79,127]

Fourth, for fiber-reinforced actuators, convergence issues arise due to the small dimension of the fiber windings in comparison with the body of the soft actuator. The quadratic beam elements used in the modeling of the fiber windings introduce bending stiffness; hence, the diameter of fiber reinforcements needs to be reduced by a factor of 2 in the FEM simulations to decrease the modeled stiffness while still reducing radial expansion.^[81,139] Furthermore, although the physical fiber reinforcements do not have compressive stiffness, the beam elements of the FEM model add some.^[81] In addition, the extra silicone layer added to fix the fiber reinforcements and not modeled in FEM result in overestimation of the deformation in the simulations. For a cylindrical fiber-reinforced SFA, for example, Decroly et al.^[139] suggest that simulation and experimental results almost perfectly match for a shift in the input pressure of 25%. From a fabrication perspective, this problem can be solved by creating the reinforcement first rather than applying the fiber winding on the surface of the SFA.^[42,142]

Many factors have an influence in the accuracy of the FEM procedures discussed here. From a materials perspective, constants obtained from uniaxial testing data might not be representative of the load conditions and multi-axial stress–strain relationships in the pressurization of SFAs. The properties of hyperelastic materials are also affected by curing temperature, mixing ratio, and degassing. Moreover, compressibility, viscoelasticity, stress softening, and the Mullins effect are usually ignored in FEM but also impact the performance of SFAs.

From a fabrication perspective, discrepancies are associated with 3D-printing artifacts or manual fabrication of soft actuators, especially for fiber-reinforced SFAs with manual wrapping of fiber. From a measurement perspective, image post-processing techniques, camera resolution, and performance of trackers also affect the accuracy of results. From an actuation perspective, the precision of pressure measurement devices needs to be considered.

Another reason for inconsistencies is the effect of gravity on simulations, as shown in Figure 11. This effect, referred to as “instability phenomenon” in the study by Peng et al.,^[149] has a dramatic influence on multi-section soft manipulators, mainly because of the compliant mechanism of soft materials. In particular, results show large differences between models with the inclusion of gravity, because this phenomenon prevents the soft manipulator from reaching a larger workspace.^[149]

FIGURE: Comparison of FEM results for a soft manipulator with the inclusion of gravity. Reproduced with permission.^[149] Copyright 2019, IEEE.

Opportunities

Despite the advantages of FEM, a number of difficulties remain. A commonly encountered issue is the lack of a comprehensive material database for the properties of silicone elastomers suited to soft robotic applications. For example, it would be ideal if the database contained raw material testing results and parameters for common model structures. This work aims to partially address this issue by consolidating material parameters that have been used in the research literature. This work also provides an introduction to the use of these parameters for the challenging task of finite element analysis of soft robotic structures.

Further investigation is required to optimize the simulations and reduce computational time while still maintaining accuracy. In this regard, the behavior of the soft actuator under pressurization can be simplified using quarter-symmetry for linear actuators and half-symmetry for bending actuators.^[79] The slow computational speed of FEM software packages discussed in this work implies that they are only usable for off-line simulation of soft robots. In contrast, Simulation Open Framework Architecture (SOFA), an open-source toolkit geared toward interactive medical simulation,^[150] allows for fast real-time control using a direct/inverse FEM solver.^[151,152] The keypoint of this approach is that the same software is used for interactive simulation and control of soft robots. The simulation framework also allows modeling of the robot's environment and interactions. The SoftRobot plugin^[153] can be used in the modeling of soft robots with tendon-driven or pneumatic actuation. The deformable model is a non-linear geometric model under linear elasticity assumption; i.e., the implementation does not consider hyperelasticity. Further information is provided in a tutorial introduction in the soft robotics toolkit.^[154] Thieffry et al.^[155] have derived a linearized model with reduced order to control the dynamics of soft robots based on the nonlinear model from SOFA. Alternatively, machine learning methods^[156,157] can be applied to learn the nonlinear kinematics of SFAs from FEM results; i.e., the finite element simulation is treated as a data generator mechanism that yields the required training data sets for artificial neural networks.^[50]

Appendix A. FEM with Abaqus

1) Start a fluid-structure interaction problem treated with a standard/implicit FEM. 2) Draw geometry or import CAD geometry. Abaqus also offers an assembly option, so parts can be individually imported and constrained within the software. 3) Create a homogeneous section for each material and assign materials for each geometry. 4) Create a new surface for the internal cavities. Use tools → view cut manager to select hidden surfaces. 5) Activate Nlgeom to include nonlinear effects of large displacements. 6) Insert boundary conditions and loads: Steps → static general → pressure. Select surfaces with internal cavities and edit the magnitude with uniform distribution and ramp amplitude. According to the soft robotics toolkit, under the incrementation tab, set 1000 as the maximum number of increments, 0.01 as the initial increment size, 1E-5 as the minimum increment size, and 0.01 as the maximum increment size. Boundary conditions → select desired fixed face → check the boxes U1, U2, and U3, and set values of 0 for fixed support. 7) Mesh the parts. Tetrahedral elements with quadratic order and hybrid formulation are recommended. The approximate global size might need to be reduced to 1–2 for the soft actuator in particular. 8) Submit job and view results.

Strain limiting layer: bottom layer → surfaces. Create a shell section and then a skin for this surface. Contact boundaries: self-contact (standard).^[158] Fiber reinforcement: mesh with quadratic beam elements. Tie constraints added after meshing. The soft robotics toolkit also provides scripts for FEM of fiber-reinforced actuators, which can be used for more advanced users of Abaqus. In addition, for unstable quasi-static simulations at higher pressures, dynamic simulations can be performed with a small damping factor.^[127]

B. FEM with Ansys

1) Create a static structural study. 2) Insert material properties into engineering data. 3) Draw geometry in Design Modeler or SpaceClaim. Alternatively, a geometry generated using standard CAD software can be imported into Design Modeler. In this case, generate geometry before proceeding. 4) Assign materials for different body parts. 5) Named selections: include a selection with all cavities in the soft actuator. This can be achieved using a section plane and the hide faces option to select hidden cavities inside the actuator. 6) Perform meshing. Recommended settings: nonlinear mechanical physics, quadratic elements, and aggressive mechanical error limits. A mesh size of 1–5 mm can be used as a starting point. 7) Set up analysis settings. Auto time stepping on → initial time step and minimum time step in the range 0.001–0.01, maximum in the range 0.01–0.05. Large deflection on. 8) Insert boundary conditions and loads: Static structural → pressure → named selection → cavities. Normal to surface, ramped. Fixed support: applied to one of the faces of the soft actuator. 9) Select desired quantities to be evaluated such as total deformation and equivalent von-Mises stress. 10) Solve and evaluate results. During the solution, it might be worth evaluating force convergence using the solution information tab. Furthermore, in the results tab, set the results to true scale (1.0).

Strain limiting layer: Concept → surface from faces → assign user defined thickness. Contact boundaries: frictional with symmetric option to minimize penetration, which results in more accurate results and realistic behavior.^[51] Fiber reinforcement: under model, select the fiber and assign a beam model type. Under connections, assign a bonded contact region between the actuator and the fiber reinforcement.

C. FEM with COMSOL

1) Create a 3D (model wizard) study → Physics: structural mechanics, select solid mechanics → stationary study. 2) Draw or import geometry. COMSOL offers a range of 3D primitives that integrated with Boolean operations allow for intuitive 3D drawing of geometries. 3) Assign material properties for each domain: solid mechanics → material models → hyperelastic material. Select user defined model parameters and insert constants. 4) Under definitions, create an explicit selection and select the internal faces of the cavity using boundary as the geometric entity level. 5) Under global definitions, create a parameter (“ramp”), so that the pressure is increased linearly during simulation. 6) Insert boundary conditions and loads: Fixed constraint: applied to one end of the soft actuator. Boundary load: pressure load type is applied to every internal cavity. Define p to be the value of final pressure multiplied by the parameter in global definitions (“ramp”). 7) Perform meshing, see note 2 as follows. 8) Set up the stationary study. Under study extensions, enable auxiliary sweep and define the parameter value list to be range(0.025,0.025,1). 9) Solve and evaluate results. Strain limiting layer: modeled by combining the behavior of the strain limiting layer into the silicone rubber parameters for the bottom layer. Contact boundaries: contact pairs. Fiber reinforcement: prescribed displacement fixed in radial directions but free in the longitudinal direction, simulates the behavior of fiber wrapping in extending actuators.

Note 1: FEM of SFAs with COMSOL Multiphysics requires a license for the nonlinear structural materials module.

Note 2: In the authors’ experience with COMSOL, a physics-controlled mesh with a “finer” element size enables convergence for a substantial range of pressures, whereas further refinement of the meshes or the use of mesh adaptation techniques showed little improvement.

Acknowledgements

The authors would like to thank Ben Routley for his assistance with simulations in COMSOL Multiphysics and Octavio Reyes for his initial work on FEM of fiber-reinforced soft actuators in Ansys.

Conflict of Interest

The authors declare no conflict of interest.