Finite Element Modeling in the Design Process of 3D Printed Pneumatic Soft Actuators and Sensors

The modeling of soft structures, actuators, and sensors is challenging, primarily due to the high nonlinearities involved in such soft robotic systems. Finite element modeling (FEM) is an effective technique to represent soft and deformable robotic systems containing geometric nonlinearities due to large mechanical deformations, material nonlinearities due to the inherent nonlinear behavior of the materials (i.e., stress-strain behavior) involved in such systems, and contact nonlinearities due to the surfaces that come into contact upon deformation. Prior to the fabrication of such soft robotic systems, FEM can be used to predict their behavior efficiently and accurately under various inputs and optimize their performance and topology to meet certain design and performance requirements. In this article, we present the implementation of FEM in the design process of directly three-dimensional (3D) printed pneumatic soft actuators and sensors to accurately predict their behavior and optimize their performance and topology. We present numerical and experimental results to show that this approach is very effective to rapidly and efficiently design the soft actuators and sensors to meet certain design requirements and to save time, modeling, design, and fabrication resources.


Introduction
Soft robots are ideal to interact safely with humans and operate in highly dynamic environments since they are made of excessively deformable and stretchable materials [1,2]. The development of these soft systems is inspired by soft biological structures present in nature, such as an elephant trunk, octopus arm, squid tentacles, and worms that are made mostly of compliant materials and liquids [3]. The structures, actuators, sensors, electronics, and power sources of a soft robotic system should primarily be made of soft materials [4]. However, the realization of completely soft robots remains a great challenge for scientists and engineers.
The soft actuators and sensors involved in every soft robotic system are the major and critical components that make such systems functional and useful. Soft robots demand dexterous soft actuators that can generate relatively low and high forces (i.e., controllable forces) to facilitate soft and adaptive interaction with their environments [5]. In addition, soft robots require robust flexible, stretchable, or compressive soft sensors that can sustain large deformations repeatedly over sustained periods, while providing consistent output signals [6]. Such reliable and stable sensors are essential for soft robots to develop reliable feedback control systems [7,8].

Design Process
The design process starts with modeling the 3D computer-aided design (CAD) models of the soft actuators and sensors, as shown in Figure 1. The structures are designed so that they can be 3D printed using affordable and accessible fused deposition modeling (FDM) 3D printers without using support materials and postprocessing. Additionally, the soft material used to 3D print the soft actuators and sensors is characterized to extract its stress-strain data that can be used to develop a hyper-elastic material model for use in FE simulations. Afterwards, the designed CAD models are imported to the FEM software along with the material model developed to simulate the soft actuators and sensors to predict their behavior and optimize their performance. Based on the FEM results, the CAD models are modified, and their simulation is repeated accordingly until the desired performance is achieved before their fabrication. Once the desired or optimized performance is obtained, the 3D CAD models are sliced in 3D printing commercial slicing software where the printing parameters are optimized to fabricate airtight and functional soft actuators and sensors. Finally, the fabricated prototypes are characterized to assess their experimental performance and to compare it with its FEM counterpart to verify and validate the simulation results.

Design Process
The design process starts with modeling the 3D computer-aided design (CAD) models of the soft actuators and sensors, as shown in Figure 1. The structures are designed so that they can be 3D printed using affordable and accessible fused deposition modeling (FDM) 3D printers without using support materials and postprocessing. Additionally, the soft material used to 3D print the soft actuators and sensors is characterized to extract its stress-strain data that can be used to develop a hyper-elastic material model for use in FE simulations. Afterwards, the designed CAD models are imported to the FEM software along with the material model developed to simulate the soft actuators and sensors to predict their behavior and optimize their performance. Based on the FEM results, the CAD models are modified, and their simulation is repeated accordingly until the desired performance is achieved before their fabrication. Once the desired or optimized performance is obtained, the 3D CAD models are sliced in 3D printing commercial slicing software where the printing parameters are optimized to fabricate airtight and functional soft actuators and sensors. Finally, the fabricated prototypes are characterized to assess their experimental performance and to compare it with its FEM counterpart to verify and validate the simulation results.

TPU Stress-Strain Data
The stress-strain relationship of the TPU used was extracted to assess its behavior using an experimental uniaxial tensile test. Eight TPU samples were printed using a longitudinal infill pattern and another eight samples were printed using a crosswise infill pattern to assess the effect of different

TPU Stress-Strain Data
The stress-strain relationship of the TPU used was extracted to assess its behavior using an experimental uniaxial tensile test. Eight TPU samples were printed using a longitudinal infill pattern and another eight samples were printed using a crosswise infill pattern to assess the effect of different infill patterns on the behavior of the TPU. The two infill patterns had an insignificant effect on the average stress-strain data of the TPU, as shown in Figure 2.
Robotics 2020, 9, x FOR PEER REVIEW  4 of 14 infill patterns on the behavior of the TPU. The two infill patterns had an insignificant effect on the average stress-strain data of the TPU, as shown in Figure 2.
The tests were conducted on the TPU samples, according to the ISO 37 standard, where all the samples were stretched by 800% at a rate of 100mm/s using an electromechanical Instron Universal Testing machine (Instron8801, Instron, Norwood, MA, U.S.A.).

Bending Soft Vacuum Actuators (SOVA)
Bioinspired bending soft vacuum actuators are developed for use in diverse soft robotic applications, including locomotion robots, adaptive grippers, artificial muscles, and modular soft robotic pneumatic hinges [32]. A single SOVA is shown in Figure 3. The SOVA actuators are activated using negative pressure (i.e., 90% vacuum). When negative pressure is applied to the actuator, consisting of five pneumatic bending hinges, it generates a bending motion, as shown in Figure 3. The objective of the FE simulations was to design a single chamber that could achieve a bending angle higher than 80° when 90% vacuum was applied. The geometry of a single chamber was modified and optimized to achieve this bending angle when 90% vacuum was applied.

Linear Soft Vacuum Actuators (LSOVA)
Similarly, linear soft vacuum actuators are developed for use in diverse soft robotic applications, including locomotion robots, adaptive grippers, artificial muscles, parallel manipulators, and soft prosthetic fingers [36,37]. A single LSOVA actuator, consisting of five pneumatic chambers, is shown in Figure 4. The objective of the FE simulations was to design a single chamber that could generate a The tests were conducted on the TPU samples, according to the ISO 37 standard, where all the samples were stretched by 800% at a rate of 100mm/s using an electromechanical Instron Universal Testing machine (Instron8801, Instron, Norwood, MA, U.S.A.).

Bending Soft Vacuum Actuators (SOVA)
Bioinspired bending soft vacuum actuators are developed for use in diverse soft robotic applications, including locomotion robots, adaptive grippers, artificial muscles, and modular soft robotic pneumatic hinges [32]. A single SOVA is shown in Figure 3. The SOVA actuators are activated using negative pressure (i.e., 90% vacuum). When negative pressure is applied to the actuator, consisting of five pneumatic bending hinges, it generates a bending motion, as shown in Figure 3. The objective of the FE simulations was to design a single chamber that could achieve a bending angle higher than 80 • when 90% vacuum was applied. The geometry of a single chamber was modified and optimized to achieve this bending angle when 90% vacuum was applied.

Linear Soft Vacuum Actuators (LSOVA)
Similarly, linear soft vacuum actuators are developed for use in diverse soft robotic applications, including locomotion robots, adaptive grippers, artificial muscles, parallel manipulators, and soft prosthetic fingers [36,37]. A single LSOVA actuator, consisting of five pneumatic chambers, is shown in Figure 4. The objective of the FE simulations was to design a single chamber that could generate a maximum linear displacement under an applied negative pressure (i.e., 95.7% vacuum). The thickness of the thin walls and the topology of a single chamber were optimized to achieve the desired performance. The wall thickness played an important role in the design process since it is directly related to the airtightness of the 3D printed LSOVA actuators. robotic pneumatic hinges [32]. A single SOVA is shown in Figure 3. The SOVA actuators are activated using negative pressure (i.e., 90% vacuum). When negative pressure is applied to the actuator, consisting of five pneumatic bending hinges, it generates a bending motion, as shown in Figure 3. The objective of the FE simulations was to design a single chamber that could achieve a bending angle higher than 80° when 90% vacuum was applied. The geometry of a single chamber was modified and optimized to achieve this bending angle when 90% vacuum was applied.

Linear Soft Vacuum Actuators (LSOVA)
Similarly, linear soft vacuum actuators are developed for use in diverse soft robotic applications, including locomotion robots, adaptive grippers, artificial muscles, parallel manipulators, and soft prosthetic fingers [36,37]. A single LSOVA actuator, consisting of five pneumatic chambers, is shown in Figure 4. The objective of the FE simulations was to design a single chamber that could generate a maximum linear displacement under an applied negative pressure (i.e., 95.7% vacuum). The thickness of the thin walls and the topology of a single chamber were optimized to achieve the desired performance. The wall thickness played an important role in the design process since it is directly related to the airtightness of the 3D printed LSOVA actuators.

Soft Pneumatic Sensing Chambers (SPSC)
Soft pneumatic sensing chambers (SPSCs), that are sensitive to the four main mechanical input modalities of compression, bending, torsion, and rectilinear displacement, are developed for use in diverse robotic applications, such as soft touch human-machine interfaces, soft wearable gloves, soft controllers and throttles, and soft bending sensors. The distinct SPSCs are shown in Figure 5. The main objective of the FE simulations was to modify and optimize the topology of the various SPSCs to obtain a linear relationship between the applied input mechanical loads and the output change in their internal volume [38]. This linearity means that the sensors can be directly 3D printed and used. The relationship between the input displacement and output pressure can be obtained by using two experimental data points and, consequently, can be used consistently since the SPSCs are stable over time, reliable, and repeatable. Therefore, there is no need for an empirical formula that requires an experimental evaluation using a specific experimental setup to obtain and describe the relationship between the input displacement and the output pressure for each 3D printed SPSC. Linearity is one of the desired performance metrics for actuators and sensors.

Soft Pneumatic Sensing Chambers (SPSC)
Soft pneumatic sensing chambers (SPSCs), that are sensitive to the four main mechanical input modalities of compression, bending, torsion, and rectilinear displacement, are developed for use in diverse robotic applications, such as soft touch human-machine interfaces, soft wearable gloves, soft controllers and throttles, and soft bending sensors. The distinct SPSCs are shown in Figure 5. The main objective of the FE simulations was to modify and optimize the topology of the various SPSCs to obtain a linear relationship between the applied input mechanical loads and the output change in their internal volume [38]. This linearity means that the sensors can be directly 3D printed and used. The relationship between the input displacement and output pressure can be obtained by using two experimental data points and, consequently, can be used consistently since the SPSCs are stable over time, reliable, and repeatable. Therefore, there is no need for an empirical formula that requires an experimental evaluation using a specific experimental setup to obtain and describe the relationship between the input displacement and the output pressure for each 3D printed SPSC. Linearity is one of the desired performance metrics for actuators and sensors.
The relationship between the input displacement and output pressure can be obtained by using two experimental data points and, consequently, can be used consistently since the SPSCs are stable over time, reliable, and repeatable. Therefore, there is no need for an empirical formula that requires an experimental evaluation using a specific experimental setup to obtain and describe the relationship between the input displacement and the output pressure for each 3D printed SPSC. Linearity is one of the desired performance metrics for actuators and sensors.

Finite Element Modeling
The FE diagram that shows and describes all the steps involved in a single simulation of a soft pneumatic actuator or sensor is shown in Figure 6.

Finite Element Modeling
The FE diagram that shows and describes all the steps involved in a single simulation of a soft pneumatic actuator or sensor is shown in Figure 6.

FEM Software
The FE simulations were performed in ANSYS Workbench (ANSYS Inc.) using a "Static Structural Analysis". The 3D CAD models of the soft actuators and sensors (i.e., geometries) were directly imported to ANSYS "Design Modeler". Here, ANSYS was used to perform the FE simulations of the soft actuators and sensors presented, since it offers various hyper-elastic materials models and it is ideal for performing static structural simulations using hyper-elastic materials.

Material Model
The average experimental stress-strain data of the TPU were imported and used in a manually defined material (i.e., NinjaFlex) in the "Engineering Data". The data were imported as "Uniaxial Test Data" for the "Hyperelastic Experimental Data". The TPU was modeled using a five-parameter Mooney-Rivlin hyper-elastic material model that fitted its average stress-strain data. The fiveparameter Mooney-Rivlin model was chosen since it best fitted the TPU experimental data compared to other available models, such as Neo-Hookean, Yeoh, and Ogden. The parameters of the material model are listed in Table 1. The model was fitted to the experimental stress-strain data using the available curve fitting tools in ANSYS.

FEM Software
The FE simulations were performed in ANSYS Workbench (ANSYS Inc.) using a "Static Structural Analysis". The 3D CAD models of the soft actuators and sensors (i.e., geometries) were directly imported to ANSYS "Design Modeler". Here, ANSYS was used to perform the FE simulations of the soft actuators and sensors presented, since it offers various hyper-elastic materials models and it is ideal for performing static structural simulations using hyper-elastic materials.

Material Model
The average experimental stress-strain data of the TPU were imported and used in a manually defined material (i.e., NinjaFlex) in the "Engineering Data". The data were imported as "Uniaxial Test Data" for the "Hyperelastic Experimental Data". The TPU was modeled using a five-parameter Mooney-Rivlin hyper-elastic material model that fitted its average stress-strain data. The five-parameter Robotics 2020, 9, 52 7 of 14 Mooney-Rivlin model was chosen since it best fitted the TPU experimental data compared to other available models, such as Neo-Hookean, Yeoh, and Ogden. The parameters of the material model are listed in Table 1. The model was fitted to the experimental stress-strain data using the available curve fitting tools in ANSYS.

Meshing
The CAD models of the soft actuators and sensors are meshed using an adaptive mesh with higher-order tetrahedral elements. A sizing function is applied when necessary to obtain a mesh with a specific element size. The mesh used in all cases is suitable for hyper-elastic materials. Since such materials can undergo large deformations, a very fine mesh is not recommended. A relatively coarse mesh is desired for convergence. The mesh size is carefully studied and selected to ensure that the solution is not dependent on it.

Analysis Settings
The "Large Deformation" option is activated in all the performed simulations to account for the large deflections exhibited by the soft structures upon activation. Additionally, a specific number of "Substeps" of 20 is imposed in each simulation to apply the loads gradually. This is an important parameter to consider since it prevents the load from being applied in one step which, in most cases, results in convergence issues, especially when dealing with soft materials that have a relatively low Young's modulus.

Contacts
Frictional contact pairs are defined between all the surfaces that come into contact upon deformation. The behavior of the contacts is set to "Symmetric" to minimize penetration, which results in more accurate results and realistic behavior. The "Augmented Lagrange" formulation is used since it is suitable for frictional and self-contacts.

Boundary Conditions
In addition to the contact pairs defined in all the structures, for SOVA, a "Fixed Support" boundary condition is applied to their base to fix them and a negative pressure (up to 90% vacuum) is applied to the internal surfaces of their hollow chambers. Similarly, for LSOVA, a "Fixed Support" is applied to their base and a negative pressure (up to 95.7% vacuum) is applied to the internal surfaces of their hollow chambers. The FEM and experimental boundary conditions applied to the various SPSCs are shown in Figure 7.
A "Fixed Support" is defined on one side of each structure and an appropriate "Displacement Support" (i.e., a linear or rotational displacement) is imposed on their opposite ends to simulate the mechanical deformations applied for each mode of deformation, as shown in Figure 7.
In addition to the contact pairs defined in all the structures, for SOVA, a "Fixed Support" boundary condition is applied to their base to fix them and a negative pressure (up to 90% vacuum) is applied to the internal surfaces of their hollow chambers. Similarly, for LSOVA, a "Fixed Support" is applied to their base and a negative pressure (up to 95.7% vacuum) is applied to the internal surfaces of their hollow chambers. The FEM and experimental boundary conditions applied to the various SPSCs are shown in Figure 7. The linear sensor is attached to a motor that generates a rectilinear stroke. (c) The bending sensor is attached to a soft flexure joint that generates a bending angle between 0° and 90° when a tendon is pulled using a linear motor. (d) The torsional sensor is attached to a servo motor that generates an angular displacement between 0° and 90°.
A "Fixed Support" is defined on one side of each structure and an appropriate "Displacement Support" (i.e., a linear or rotational displacement) is imposed on their opposite ends to simulate the mechanical deformations applied for each mode of deformation, as shown in Figure 7.

FEM Results
Various important results can be effectively used to assess and interpret the simulation of a soft body, as shown in Figure 6. The total deformation is the result of the load applied to a soft body and it shows how the body is deformed at different stages of the load applied, as shown in Figure 8. This deformation can be directly compared with its experimental counterpart to verify it. The volume change of the various SPSC was assessed at different stages of the applied load to obtain a relationship between this change and the applied mechanical load, as shown in Figure 9. As stated earlier, the aim was to obtain a linear relationship between the applied mechanical loads and the change in their internal volume. Ideally, a relationship exists between the change in the internal volume (V 1 /V 2 ) of the soft chambers and the experimental pressure change (P 1 V 1 = P 2 V 2 ), obtained due to the mechanical deformation applied.

FEM Results
Various important results can be effectively used to assess and interpret the simulation of a soft body, as shown in Figure 6. The total deformation is the result of the load applied to a soft body and it shows how the body is deformed at different stages of the load applied, as shown in Figure 8. This deformation can be directly compared with its experimental counterpart to verify it. The volume change of the various SPSC was assessed at different stages of the applied load to obtain a relationship between this change and the applied mechanical load, as shown in Figure 9. As stated earlier, the aim was to obtain a linear relationship between the applied mechanical loads and the change in their internal volume. Ideally, a relationship exists between the change in the internal volume (V1/V2) of the soft chambers and the experimental pressure change (P1V1 = P2V2), obtained due to the mechanical deformation applied.
Additionally, the blocked force, which is the reaction force obtained when the actuators are activated but restricted from moving at both sides, can be computed and compared to its experimental counterpart to verify it, as shown in Figure 10. This reaction force is a measure of the output force that a soft actuator can generate.
Finally, the stresses and strain resulting from the load applied can be assessed as presented in Table 2. These two parameters are very important since they should be minimized as much as possible to enhance the lifetime of the soft actuators and sensors developed.    Figure 9. FEM deformed shape for (a) the push button when 24.5° remote rotating displacement is applied on the crank which not shown for visual purposes, (b) linear sensor when a 10 mm displacement is applied, (c) bending sensor when a 90° bending angle is applied by pulling the tendon, and (d) torsional sensor when a 90° remote rotating displacement is applied on its end. Additionally, the blocked force, which is the reaction force obtained when the actuators are activated but restricted from moving at both sides, can be computed and compared to its experimental counterpart to verify it, as shown in Figure 10. This reaction force is a measure of the output force that a soft actuator can generate.

Verification of FEM Results
The FE results, in terms of total deformation and output blocked force for SOVA and LSOVA, were compared to their experimental counterparts to assess their validity. In addition, the relationship between the input mechanical loads and the change in the air pressure was obtained experimentally to verify that this relationship is linear, as the FE simulations of the SPSCs predicted. Finally, the stresses and strain resulting from the load applied can be assessed as presented in Table 2. These two parameters are very important since they should be minimized as much as possible to enhance the lifetime of the soft actuators and sensors developed.

Verification of FEM Results
The FE results, in terms of total deformation and output blocked force for SOVA and LSOVA, were compared to their experimental counterparts to assess their validity. In addition, the relationship between the input mechanical loads and the change in the air pressure was obtained experimentally to verify that this relationship is linear, as the FE simulations of the SPSCs predicted.

Total Deformation
The experimental total deformation of SOVA and LSOVA actuators were obtained and compared with their FEM counterparts. The difference between the experimental bending angle of SOVA of 265.50 • and its FEM bending angle of 269.63 • is 1.56%. Similarly, the difference between the experimental linear displacement of LSOVA of 35.03 mm and its FEM linear displacement of 39.47 mm is 12.67% (this difference is higher due to the printing artifacts in LSOVA that interfered with its linear displacement [36]). These differences verify that the FEM results can accurately predict the deformation of the soft pneumatic bending and linear actuators.

Blocked Force
Likewise, the experimental blocked force of SOVA and LSOVA actuators were obtained and compared with their FEM counterparts. The difference between the experimental blocked force of SOVA of 10.72N and the FEM blocked force of 12.33N is 15.02%. This difference is higher since, at higher pressures, the force sensor moved slightly, which reduced the experimental blocked force of SOVA and the two actuators did not perfectly align, as shown in Figure 10b in the experimental setup, compared to the simulation [32]. Although ultimate attention was provided to fabricate two identical SOVA actuators, there were differences in their mechanical output under the same pressure. Similarly, the difference between the experimental blocked force of LSOVA of 27.66N and its FEM blocked force of 28.66N is 3.62%. These differences verify, again, that the FEM results can accurately predict the blocked force of the soft pneumatic bending and linear actuators.

Volume Change Linearity
The change in the internal volume of the various SPSCs changed linearly with the applied load [38]. The experimental load was ramped up and down by a step of 10 • , as shown in Figure 11, where at each step enough time (i.e., 5 s) was given for the various SPSCs to ensure that the pressure value recorded was the final steady-state value. This change was verified experimentally by obtaining a linear relationship between the internal pressure change in the various SPSCs and the corresponding applied load, as shown in Figure 7, since a relationship exists between the change in the internal volume (V 1 /V 2 ) and the experimental pressure change (P 1 V 1 = P 2 V 2 ). V 1 is the initial volume of single SPSC and V 2 is the internal volume at a corresponding input displacement. V 1 is computed from the initial CAD models of the SPSCs and V 2 is computed from the deformed CAD models, resulting from the FE simulations. Figure 11a shows the experimental linear relationship between the input displacement and the output pressure for the bending sensor and its corresponding FEM linear relationship between the same input displacement and the output volume change. The linear fit of both relationships resulted in a coefficient of determination of R 2 = 0.995 for the experimental data and R 2 = 0.988 for the FEM data. Similarly, Figure 11b shows the experimental linear relationship between the input displacement and the output pressure for the torsional sensor and its corresponding FEM linear relationship between the same input displacement and the output volume change. The linear fit of both relationships resulted in a coefficient of determination of R 2 = 0.996 for the experimental data and R 2 = 0.997 for the FEM data. These results show that the FEM predicted the behavior of the physical bending and torsional sensors accurately.
applied load, as shown in Figure 7, since a relationship exists between the change in the internal volume (V1/V2) and the experimental pressure change (P1V1 = P2V2). V1 is the initial volume of single SPSC and V2 is the internal volume at a corresponding input displacement. V1 is computed from the initial CAD models of the SPSCs and V2 is computed from the deformed CAD models, resulting from the FE simulations. Figure 11a shows the experimental linear relationship between the input displacement and the output pressure for the bending sensor and its corresponding FEM linear relationship between the same input displacement and the output volume change. The linear fit of both relationships resulted in a coefficient of determination of R 2 = 0.995 for the experimental data and R 2 = 0.988 for the FEM data. Similarly, Figure 11b shows the experimental linear relationship between the input displacement and the output pressure for the torsional sensor and its corresponding FEM linear relationship between the same input displacement and the output volume change. The linear fit of both relationships resulted in a coefficient of determination of R 2 = 0.996 for the experimental data and R 2 = 0.997 for the FEM data. These results show that the FEM predicted the behavior of the physical bending and torsional sensors accurately.

Stress-Strain Results
The stress and strain results were obtained to assess the regions of each structure that experience high stresses and large strains, as presented in Table 2. High stress and strain values led to a significantly shorter lifetime due to the failure of the material (i.e., separation of the printed layers). However, it is important to note that, although in some regions the stress values might be high, they might not have a significant effect on the lifetime or performance of the structure, since they do not lead to the failure of the actuator or sensor. For instance, for the push button, the high stress leads to a shorter lifetime compared with the lifetime of the remaining SPSCs [38]. The push button sustained 60,000 activation cycles, while the bending sensor, which exhibits higher stresses and strains, sustained more than 150,000 activation cycles, since these high stresses and strains occurred on the top of the sensor in the groove area, which is far from the thin walls. However, for the LSOVA, since the stress concentration occurs on the corners of the thin walls, their lifetime was significantly shorter (~20,000 activation cycles) compared to the SOVA, which sustained more than 120,000 cycles. The stress and strain data were taken into consideration in the design process to minimize them and to ensure that they minimally affected the lifetime of the soft actuators and sensors.

Challenges and Recommendations
One challenge encountered is the distortion of some elements in the FE model due to the large mechanical deformations. However, this issue was alleviated by incorporating a coarser mesh, when needed, that was suitable for hyper-elastic materials. The mesh used in each case was selected to verify that the results were accurate and mesh independent.
Another possible challenge was the difficulty to achieve convergence (i.e., solution) due to the multiple contacts involved in each structure. The reason is that the "Symmetric" frictional contacts imposed imply that no penetration is allowed between the surfaces that come into contact. Although in the contact settings, the "Penetration Tolerance" is set to "Program Controlled", in some cases it is recommended to put a reasonable manual value to relax this constraint to achieve convergence. This tolerance should be chosen such that the simulations must be realistic where no excess penetration is allowed.
Another recommendation for simulating pneumatic soft bodies successfully is that there is no need to include all the details involved in the real geometry in the simulation. However, the simulated geometry must be simplified as much as possible as long as it can be still considered as the real body and captures the most significant features. This simplification will eliminate several problems that might be encountered during meshing and solving to save considerable amounts of computing power and time. For instance, for SOVA and LSOVA, the holes connecting the different chambers were ignored in the simulated body, since they have an insignificant effect on the final solution. Similarly, the input holes in the different SPSCs are ignored in the simulated bodies, as shown in Figures 5 and 9.

An Effective Solution to a Soft Robotic Challenge
Developing soft pneumatic actuators and sensors for soft robotic systems using various manufacturing techniques and technologies is both time and material consuming. As stated earlier, although the modeling of soft actuators and sensors is an indispensable part of the soft robotics field, it is a very challenging problem, due to nonlinearities. However, in this study, we showed that FEM can be effectively used to handle the nonlinearities involved in modeling soft robotic actuators and sensors to optimize and predict their behavior accurately. Therefore, this study has effectively presented a roadmap for efficiently designing and modeling soft pneumatic actuators and sensors.

Future Work
In future work, various materials can be considered to fabricate various soft actuators and sensors using different manufacturing technologies with different geometries where their appropriate material models can be developed for use in simulations. Additionally, coupled simulations that account for fluid-structure interaction can be used for both actuators, and sensors and dynamic simulations can be performed to assess their behavior under dynamic conditions where inertial forces are significant. Finally, FEM is not limited to soft actuators and sensors that are made of TPU or silicone materials that are not responsive (i.e., cannot be controlled) to external stimuli. Smart materials (e.g., shape memory alloys, dielectric, ionic polymer-metal composites, coiled polymer fibers, hydrogels, humidity-responsive materials, magnetic responsive structures, etc.) that can be excited by external stimuli can be used as actuators and sensors in some cases and, therefore, can be also simulated.

Conclusions
The modeling of soft and highly deformable robotic structures and systems is challenging, due to nonlinearities, such as geometric, material, and contact nonlinearities. Although the modeling of such soft robotic systems is challenging, it should be integrated into the design process of these systems before their fabrication to efficiently and accurately predict their behavior and to optimize their performance and topology to meet certain design and performance requirements. Soft pneumatic actuators and sensors that can be directly 3D printed using open source and low-cost FDM 3D printers can be optimized in terms of performance and topology using FEM that predicts their behavior and performance accurately. This means that tremendous amounts of time and resources can be saved, and that FEM can be efficiently used to handle the nonlinearities involved in modeling soft robotic actuators and sensors, which are essential elements of soft robotic systems.