Numerical modeling of a 3D printed soft pneumatic actuator

Soft robotics represents a significant evolution in the field of robotics, characterized by systems that can adapt their shape and functionality to interact safely and effectively with complex environments and living organisms (Tawk and Alici, 2020; Kim et al., 2019). This paradigm shift draws inspiration from biological systems, allowing soft robots to exhibit high compliance and flexibility, which enable them to address challenges that traditional rigid robots cannot tackle (Mosadegh et al., 2014; Trivedi et al., 2008; Laschi et al., 2016; Tang et al., 2022). This feature is particularly useful in situations involving direct human-robot interaction (Schmitt et al., 2018). There are several types of actuation systems for soft robots, each tailored to the application’s needs and available energy source. Based on their working principle, soft actuators can be categorized as cable-driven (CDA), shape memory alloy (SMA) thermal-driven, electro-driven (EMA) and fluidic elastomer actuators (FEAs) (Whitesides, 2018). These latter actuators, operating using pressurized fluids, are the focus of the present work. A specific subcategory of these actuators is soft pneumatic actuators (SPAs), which use pressurized air as the working fluid. SPAs have a relatively simple structure, consisting of a series of chambers connected by a central duct through which compressed air is fed. When pressurized, the actuator’s walls expand, deforming the entire structure (Hines et al., 2017; Polygerinos et al., 2017). The type of action, such as bending, twisting or extending, can be programmed based on the actuator’s morphological characteristics. Among the various models, the present study takes into consideration a “planar bending” motion. This model generates an expansion of the aligned chambers due to increased pressure, leading to contact between them and resulting in pure bending of the entire actuator on a plane. The actuator’s response to applied pressure depends on both the materials used and the actuator’s geometry. Therefore, geometrical and material properties are considered critical design parameters. The influence of geometrical parameters such as the gap size between two consecutive bellows, wall thickness, bottom layer thickness, chamber width, number of chambers and the type of cross-section on the bending angle of SPAs was extensively studied in Stano et al. (2020), Gao et al. (2022), Sun et al. (2019), Hu et al. (2018). Typically, two methods are applied to characterize actuator responses. The first method involves creating and testing numerous specimens to gather sufficient data for a full factorial analysis, evaluating the effect of plausible design factors. This approach allows an accurate understanding of a specific actuator’s behavior under certain conditions. However, collecting enough data requires a significant investment of time and resources. In addition, if the conditions or the actuator itself changes, the tests must be repeated. The alternative method involves using Finite Element Method (FEM) simulations. Once certain boundary conditions are identified, these simulations can replicate the behavior of a specific element. This approach is more efficient as it eliminates the need for physical testing and can easily adapt to changes in conditions or the actuator itself. Nevertheless, predicting the behavior of soft actuators can be challenging due to their inherent flexibility and nonlinear behavior. Under applied pressure, SPAs exhibit complex deformation due to their geometry and the hyperelastic – and occasionally anisotropic – properties of the material used. Several papers emphasize the importance of FEM for predicting the behavior of these soft devices, given their inherent nonlinearity resulting from large deformations, material behavior and surface contacts (Tawk, and Alici, 2020; Lalegani Dezaki, 2023; Rad, 2022; Wang et al., 2017). FEM modeling allows engineers to simulate these deformations, providing insights into the actuator’s range of motion, force output and deformation (Maruthavanan et al., 2021). This information is critical for designing SPAs that precisely meet various application requirements. Soft systems can benefit significantly from FEM analysis for several reasons:

• By comparing the outcomes of numerous simulations, performance can be evaluated and the optimal actuator geometry can be chosen to meet design phase specifications (Keong and Hua, 2018), saving time and resources.

• Once that a satisfactory digital twin of the real actuator is produced, simulation settings can be adjusted repeatedly to align the real behavior with the simulated one. This is particularly useful when dealing with highly elastic materials, as their behavior can be significantly influenced by external variables such as humidity, temperature, load types and positions, constraints.

Indeed, one of the main challenges in developing SPAs FEM analyses is accurately replicating the material’s behavior. Despite the existence of numerous models in the literature that describe materials’ hyper-elastic behavior (Xavier et al., 2022; Batsuren and Yun, 2019; Steck et al., 2019), the manufacturing process and geometric characteristics of the specimen can complicate accurate modeling.

In this study, we developed a FEM model to simulate the behavior of a planar-bending SPA fabricated from commercial thermoplastic polyurethane (TPU) using fused filament fabrication (FFF) 3D printing with the ultimate goal of producing a valid design tool for the production of actuators that are optimized according to operating conditions.

The present study meticulously evaluated the influence of three key parameters: the number of bellows (N), the thickness of the bellow wall (T) and the operating pressure (P) on the performance of pure-bending SPAs. Specifically, performances were evaluated in terms of bending angle (α) and extreme force exerted (F). To comprehensively understand this relationship, we selected three different wall thickness values (1.6 mm, 2 mm and 2.4 mm), three operating pressures (1 bar, 2 bar and 3 bar) and three numbers of bellows (9, 11 and 13), as shown in Table 1. Such values were selected considering a generic planar bending SPA, plausible overall dimensions for integration into a hand gripper and according to literature values and preliminary results carried out in previous research (Torzini et al., 2024). This resulted in a total of 27 possible parameter combinations (considering also the operating pressure among these) which describe a plausible design range of the device. The overall methodology followed in the study is presented in Figure 1. To acquire ground truth values to support the development of a faithful model for predicting/simulating the behavior of the considered actuator, performance tests were conducted on real actuators, covering all configurations, to understand their actual behavior under various conditions. Results are reported in the bending behavior and blocked force sections. Experimental data was used to drive the optimization of a material model used in a static structural FEM analysis simulating the actual test within Ansys (2024) R3. In a preliminary phase, three material models were taken from the literature (Reppel and Weinberg, 2019; Tawk, 2022) and were tested to select the one best fitting the experimental data in a series of FEM simulations performed on a reduced number of configurations (from 27 to 9), obtained by applying a fractional factorial design to speed up the process. Experimental data was then used to optimize the parameters defining the selected material model in a series of FEM analyses, thanks to an optimization algorithm applied to the same fractional factorial design. Four material parameters defining a 2nd order Ogden hyperelastic model were tuned to best fit simulated bending data to experimental measures according to an error function estimating the position of the actuator’s tip. The resulting material model was validated with respect to experimental data on the full factorial design.

Finally, a different – more comprehensive – strategy was attempted for optimizing the material model, taking into account the full shape of the pressurized actuator and trying to match it to the experimental data. Moreover, this last approach simultaneously considered three actuating pressures to produce general results. This complex FEM-driven optimization approach achieved an overall 9,7% reduction in error, resulting in a more accurate simulation of the actuator’s behavior.

The motion path in soft robotics is typically embedded within the actuator design; geometry plays an active role in determining the motion law of an actuator: its correct identification is critically important to assure the production of an actuator fit for its task. In the present study, SolidWorks, (2024) 3D computer-aided design (CAD) environment was used to model the soft actuators. Considering production using FFF 3D printing process, the design of the main features of the actuator was driven by a reduction of the overhanging surfaces to avoid the need for supports during the printing process. Apart from the three design and operating parameters examined in this study, the choice of the cross-section is a critical factor in the mechanical performance. According to the literature, rectangular cross-sections perform better than other shapes when applying the same air pressure in terms of bending and force exerted (Schmitt et al., 2018). Consequently, only rectangular cross-sections were taken into account. The size and form of the SPA are intended to take inspiration from human fingers because these actuators are frequently utilized in applications involving mimicking the characteristics of the human hand or interacting with it (Bianchi, 2016). As shown in Figure 2, the overall length of the actuator is set at 115 mm, the height at 22 mm and the width at 16 mm.

All the other dimensions were chosen during a preliminary test campaign that allowed the identification of plausible values to exclude air leakages and allow realization by 3D printing. In particular, T = 1,6 mm was chosen as the minimum thickness to allow for an adequate number of shell layers; the T-values used, shown in Table 1, allow for the deposition of 4, 5 and 6 shell layers, respectively.

In this study, SPAs were produced with an FFF 3D printer (Original Prusa i3 MK3+). Because the actuators were manufactured without support, postprocessing was unnecessary. This was accomplished by printing the actuators sideways, with the bending plane parallel to the building platform; this removed the requirement for internal support and left the actuator with only its upper part overhanging. Author-conducted studies in the past have shown that this orientation produced deposition paths that are beneficial for the reduction of air leakage. The soft actuators were 3D printed using a commercial TPU with a shore hardness A of 85 called NinjaFlex®. PrusaSlicer V2.8.1 was used to read and process STL files. A nozzle with a diameter of 0.4 mm and a layer height of 0,15 mm has been used for printing; an extrusion temperature of 245°C was selected to avoid the formation of micro-holes generated by incomplete melting of the filament, which would lead to air leakage. To minimize errors related to the production process, three samples were produced for each of the nine geometrical combinations, resulting in a total of 27 3D-printed actuators.

The experiment for assessing the bending capabilities of the specimens’ set was carried out using a fixed rig holding the actuator in a controlled position and pressurizing it. A portable compressor (Hyundai™ KWU750-24L) supplies the pressurized air. The airflow is regulated by a pressure gauge at the compressor’s outlet, maintaining a pressure range of 1–3 bar for the experiments. An Arduino® Uno board (Arduino AG, Chiasso CH) combined with an electrical circuit allows control, via PWM signal, of an electro-pneumatic transducer, which regulates the pressure in the circuit. The actuators are mounted horizontally on a 3D-printed custom support; pressurized air is delivered through a 5 mm PVC pipe, as illustrated in Figure 3. A CCD Mono 480p camera (IDS UI 2310 M) captures images of the pressurized actuator for measurement purposes.

A separate setup was used to measure the force applied at the pressurized actuators’ tip: a digital scale with a 0.1 g sensitivity was used. As seen in Figure 4, the actuator was fixed at one end and was positioned 5 mm above the scale’s plate to allow it to bend and apply force to a 3D-printed square when it was pressured.

To evaluate whether the FEM models could accurately predict the behavior of the real models, the bending angle (α), as shown in Figure 5, was measured at maximum deflection.

This measurement of the bending angle (α) formed by the horizontal axis and the line connecting the support and tip of the actuator (L) was performed on the captured image using a custom computer vision routine coded in MATLAB®. After calculating the length of L using known objects captured in the image as a reference, the displacements with respect to Z and Y (UZ and UY) of the actuator tip relative to its unpressurized horizontal position were determined based on the known angle (α). The bending test was conducted with variations in the three parameters mentioned before, resulting in 27 unique combinations. For each combination, the test was performed on three identical specimens, with each specimen tested three times, yielding a total of 243 images. Once all the images were processed, the average values were calculated and the results are shown in Table 2.

Blocked force is a crucial performance metric for soft actuators, quantifying the force produced at the actuator’s tip. This metric reflects the actuator’s efficiency in converting input pressure into an output force. In these tests, one end of the soft actuators is securely fixed, leaving the other end free to bend. The force readings are initially obtained in grams-force (gf). These values are converted to Newtons (N) using the conversion factor of 1 gf ≈ 0.0098 N to ensure accurate characterization of the actuator’s force generation capabilities. As the input air pressure is incrementally increased in 1 bar steps up to a maximum of 3 bars, the resulting output force generated by the actuator correspondingly increases. Each combination of geometrical parameters was tested on three distinct actuators, with each test repeated three times; obtained values were averaged, the confidence intervals were calculated (Table 3).

To quantitatively evaluate the statistical significance of the investigated design and operating parameters – number of bellows (⁠$N$⁠), wall thickness (⁠$T$⁠) and operating pressure (⁠$P$⁠) – a full-factorial Analysis of Variance (ANOVA) including main effects and two-way interaction terms (⁠$N×T$⁠, $N×P$⁠, $T×P$⁠) was performed on both the bending angle (⁠$α$⁠) and the blocked force (⁠$F$⁠).

The ANOVA results for the bending angle (⁠$α$⁠) are summarized in Table 4. The analysis demonstrates that all main factors and two-way interactions exert a statistically significant influence (⁠$p<0.001$⁠).

Operating Pressure (⁠$P$⁠) is the dominant factor driving deflection (⁠$F=4766.78$⁠, $p=4.53×10-242$⁠) accounting 68% of the total sum of squares; wall thickness constitutes the second most influential parameter (⁠$F=1970.45$⁠, $p=8.84×10-183$⁠) whereas bellows number (⁠$F=280.70$⁠, $p=1.29×10-71$⁠) exhibits a smaller yet statistically significant contribution . Among the two-ways interaction terms, the coupled effect of wall thickness and pressure (⁠$T×P$⁠) displays the highest significance (⁠$F=104.68$⁠, $p=3.07×10-57$⁠), followed by $N×T$ (⁠$F=56.72$⁠, $p=2.75×10-36$⁠) and $N×P$ (⁠$F=7.44$⁠, $p=9.62×10-6$⁠).

Regarding the blocked force (⁠$F$⁠) (Table 5), the ANOVA table confirms that operating pressure plays the primary role in force generation (⁠$F=18611.31$⁠, $p=1.01×10-249$⁠), contributing 71% of the overall variance. Wall thickness also exhibits an extremely strong influence (⁠$F=5305.48$⁠, $p=2.12×10-189$⁠), while the number of bellows remains significant (⁠$F=674.79$⁠, $p=1.50×10-95$⁠). Furthermore, strong two-way interaction effects were observed, particularly between $T×P$(⁠$F=519.22$⁠, $p=5.06×10-112$⁠) and between $N×T$(⁠$F=183.97$⁠, $p=1.44×10-69$⁠), highlighting a complex structural-pneumatic coupling in force output.

The “Static Structural Analysis” module of ANSYS Workbench 2019 R3 (ANSYS Inc.) has been used to perform FE simulations. The soft actuators’ 3D CAD models (for n = 9, n = 11 and n = 13) were loaded directly into the “Design Modeler” program for ANSYS. ANSYS conveniently incorporates several hyperelastic material models ready to be applied.

To accurately replicate the significant deformations exhibited by the actuators when internal pressure was applied, the “large deflection” option was selected in the analysis settings. To improve the stability of the simulation, a gradual load application was used, ranging from 0.1’’ to 1.0’’ with a time step of 1e-4.’’ “Fixed Support” was used at the actuator end to hold it in place while positive pressure was applied to all internal surfaces of the actuator’s cavity to recreate the test scenarios. The simulations also took gravity into account. Pressure was applied to the internal surfaces of the actuators using the desired value in a ramped form. “Frictionless contacts” were considered between the actuator bellows. Nonlinear mechanical quadratic elements were used to mesh the CAD models of the soft actuators. The mesh used in all cases is suitable for hyper-elastic materials. An excessively fine mesh is not recommended because such materials can experience significant deformations (Whitesides, 2018); a sensitivity analysis, testing element sizes of up to 1e-4m, led to the selection of a 1e-3 m as valid resolution for the application. To further explore the possibility of adequately increasing and adapting the mesh to the pressurized actuator, a simulation using a mesh density of 1e-3 m and activating the “Nonlinear Adaptive Region” option was run. Despite the significantly longer runtime compared to a simulation without that option, it did not yield a significantly more accurate result. Therefore, it was decided not to use this option for subsequent runs. Then, two configurations were modeled differently based on the simulated test type: a bending condition, where the actuator is locked at one end and left free to flex at the other end, with planar deflections being calculated and a force configuration, where the free end of the actuator contacts an undeformable plate and the reaction force is calculated.

To monitor displacements in the Z–Y plane during bending simulations, probes were inserted to track the position of the tip point on the bending plane throughout all the simulations. The bending error (BE) function was defined accordingly to equation (1) using a “User Defined Function” tool of Ansys: this function allows to calculate the difference between the displacements measured in the tests compared to the simulated model values and effectively expresses the Euclidean distance error separating the estimated tip position of the actuator and the experimental one:

where $Ufz$ and $Ufy$ are the displacement values of the actuator tips, with respect to the Z and y axes, that were measured in the corresponding real experiment. The values of $UZ$ and $UY$ are the actuator tip displacements obtained during simulations in the Z and y axis, respectively.

A fixed flat plate, modeled as a square of 25 × 25 × 1 mm made of steel, was included in the model to simulate the digital scale plate. “Frictionless contacts” between the plate and the actuator were considered. A probe was placed on the “fixed support,” positioned on the plate, to measure the force exerted by the actuator’s tip.

TPU exhibits large deformations typical of rubber-like materials, which cannot be accurately modeled using linear elasticity. For the selection of the hyperelastic material model, an extensive analysis of the state of the art was conducted. Specifically, works dealing with FEM modeling of soft actuators made by FFF 3D printing were selected. However, this analysis highlighted the lack of studies and the challenges in modeling this type of material, especially when combined with the manufacturing process, which further influences the behavior of the resulting elements. Once the sorting has been done, two hyperelastic models have been considered due to their ease of modeling and suitability for the structural characteristics of 3D printed TPU. Tawk and Alici, (2020) proposed the Mooney–Rivlin 5 parameters model, which is an improvement of the Neo-Hookean model. It is ideal for elastomers and rubbers in the medium-to-large deformation range and is typically shown as a polynomial curve, as specified by equation (2):

where W is strain energy potential, I_n are the invariants of the deformation tensor, C_nn and D are the material parameters and J is the volumetric invariant of the deformation tensor. Tawk et al. then fit the model to the experimental stress-strain data using the available curve fitting tools in ANSYS and the parameters are listed in Table 6.

Experiments were also carried out by Reppel and Weinberg (Reppel and Weinberg, 2019) to study the elastic characteristics of printed NinjaFlex® TPU. They used a second-order Ogden model to match the experimental data. The Ogden model is a constitutive model that can be used on a variety of rubber-like materials, polymers and biological tissues (Keong and Hua, 2018). The Ogden hyperelastic model for incompressible isotropic material under uniaxial load is defined by equation (3):

where W  is the strain energy density function, λ is the stretch, N is model order and p and are material coefficients. In Reppel’s work, two different sets of material parameters were obtained based on a uniaxial tension test made on unified DIN EN ISO 527–2 type 1BA specimens, characterized by a different number of deposition perimeters parallel to the external surfaces of the object; the two parameters set are listed in Table 7.

The first set of simulations aimed to identify the material model and parameters that most closely matched the outcomes of the actual tests. Regarding simulations, compared to the 27 combinations of geometric and operational parameters seen in Table 2, only nine are needed to evaluate the effect of the three factors N, T and P. To determine the degree of resemblance between the simulated model and the real test, the BE function computed using equation (1) was used. Figure 6 reports the measured errors for all the material models across the selected configurations. The findings indicate that the “002” specimen material characteristics in the second-order Ogden material model best approximates the actual behavior of the printed actuators and was accordingly selected as the starting model for the further optimization phase planned ahead.

For optimization purposes, Ogden 002 was selected as the reference material model. Despite these efforts, there remains a significant BE between the experimental and simulated models. In fact, even considering a standard deviation in real angle measurements of 1.98°, the average error derived from simulations using the material parameters of Reppel and Weinberg, (2019) is 2.64°. To address this discrepancy, has been decided to implement a “direct optimization” procedure. Evidently, the standard approach for further improving the obtained results would involve performing tensile tests on specimens made of the same material used for the actuators, printed under similar conditions. This seemingly rigorous methodology has significant inherent limitations. The samples thus obtained would inevitably exhibit a nonisotropic material structure, with a distribution of mechanical properties that would diverge significantly from that of the final actuator. Moreover, the deposition path of the material during fabrication would structurally differ between the test samples and the actuator, introducing variability factors that would compromise the reliability of the experimental results. Consequently, after a thorough methodological evaluation, we opted for an alternative approach involving direct connections between input and output. This strategy allows us to optimize material parameters according to the performance observed on the actuator, providing more immediate and contextualized experiential feedback than conventional methods. This methodological choice represents a reasoned compromise between scientific rigor and practical characterization needs, aiming to provide a description of material properties that is not only accurate but more importantly, contextually relevant. In this optimization process, the BE [equation (1)] served as an objective function, allowing the process to find the combination of material parameters that would minimize the function. A key aspect of the optimization process concerns the choice of the best solving method. Therefore, available Ansys optimization algorithms such as adaptive single-objective (ASO), multiobjective genetic algorithm (MOGA) and Screening were preliminarily tested; the same optimization process was launched using the three different methods, measuring time and results achieved. In conclusion, ASO was set up as the best resolution method; for further information on the methods, see the Ansys knowledge source (Ansys, 2024). An initial set of 20 samples was selected, with a maximum of 80 assessments and a convergence tolerance of 1e-6m was established. Starting from the material parameters values of Ogden 002, the four material parameters (α1, μ1, α2 and μ2) were set as optimization variables, allowing for a possible variation of ± 10%, with the optimization objective being the minimization of BE – equation (1). As previously mentioned, the optimization process was applied to 9 SPAs configurations, independently.

Nine different sets of material parameters, each specific to a given combination of N, T and P, were obtained once the optimization procedure for the nine combinations was finished. As seen in Figure 7, optimized parameters at the optimum design points show a high degree of agreement between simulated and real behavior.

At the end of the optimization processes, we looked for some patterns in the material parameters that would allow to connect actuators with similar characteristics design features; however, no specific parameter set was found to match actuators with similar traits (Figure 8). Ultimately, obtained values are spread on the entire solution space granted to the optimization algorithm, with the only exception of μ2, which converged toward a significant reduction in value. However, considering equation (3), it is important to note the combined effect of parameters, which do not allow independent evaluations.

BEequation (1) was used at the end of the optimization phase to calculate the error between the position of the actual and simulated actuator tip. A comparison between the results of the simulations using Ogden 002 and those with our optimized parameters shows an improvement in the prediction of actuator behavior in the latter case. As shown in Table 8, there is a mean error reduction of 37% and a median error reduction of 26%, obtained using each time the optimized material model for that specific combination of design parameters.

To determine if the optimized values for a particular combination of geometrical parameters could faithfully simulate the behavior of the identical actuator operating at varying pressures, we performed 18 additional FEM simulations. Material parameters obtained by simulating an actuator with a specific combination of geometrical parameters at a given pressure were used to simulate the behavior of the same actuator at different pressures. We repeated this procedure using the literature material model” Ogden 002” from Reppel and Weinberg (2019), to ensure a fair comparison. The BE objective function [equation (1)] was again used as a comparison measure to evaluate the performance of the simulations against the corresponding real cases. After calculating the average and median of these distances, reported in Figure 9, it was observed that using the optimized parameters reduced the distance between the tips of the simulated and real actuators globally by approximately 12%.

After examining the data in Figure 9, it was found that both of the material models used had inaccuracies localized at particular configurations, many of them found for pressure values p = 3 bar. This evidently imposes some limitations on the possible extrapolation of material parameters to be used on SPAs pressurized at higher pressures. As the final term for evaluating the behavior of the simulated actuator compared with the real one, it was determined to compute through simulations the force applied at the tip and then compare the findings with actual results. To do this, the FEM simulation was set up to mimic the actual test circumstances. A probe was placed at the actuator’s end to measure the force applied at the desired pressure. The analysis settings used in these simulations were similar to previous ones. Both material models were used, and Figure 10 shows the outcomes. As can be seen, the use of optimized parameters did not result in an improvement over the reference ones in predicting the forces exerted by the simulations.

To quantitatively evaluate the predictive fidelity and trend consistency of the primary finite element model across the nine representative specimen configurations of N, T and P, a sample-by-sample comparative regression analysis was conducted for both bending angle and blocked force. To account for experimental variance across treatment groups, weighted least squares (WLS) regression – weighted by the inverse experimental variance (w_i = 1/σ_i²) – was fitted to the empirical data set. Conversely, given the deterministic nature of finite element simulations, which lack dispersion across repetitions, ordinary least squares (OLS) regression was applied to the ANSYS numerical predictions across the categorical specimen index.

As illustrated in Figure 11, the baseline numerical predictions demonstrate strong point-by-point agreement with empirical observations across all tested configurations, successfully tracking localized response variations:

• Bending Angle (α):

  • Experimental Trend (WLS): y  =  0.93x  +  31.35 (R² = 0.14).

  • ANSYS Simulation Trend (OLS): y  =  0.81x  +  29.84 (R² = 0.12).

• Blocked Force (F):

  • Experimental Trend (WLS): y  =  0.08x  +  1.09 (R² = 0.36)

  • ANSYS Simulation Trend (OLS): y  =  0.13x  +  1.09 (R² = 0.19).

A comparative evaluation of the two outputs reveals a clear distinction in predictive capability between free kinematic deflection and constrained force generation:

• For what concern the bending angle, the baseline FEA model exhibits solid predictive accuracy and strong physical consistency with experimental behavior. The regression slopes (beta = 0.93 for experimental WLS vs 0.81 for numerical OLS) and baseline intercepts (31.4 deg vs 29.8 deg) demonstrate close alignment across all tested configurations. The FEA framework reliably reproduces the unconstrained hyperelastic deformation driven by internal pressure, proving to be an accurate tool for evaluating the kinematic bending response of the actuators.

• In contrast, predicting the blocked force reveals noticeable difficulties and greater divergence from empirical measurements. Although the simulation captures the general stepwise increase in force with rising actuation pressure, it significantly overestimates the force output at higher pressure levels (specifically at 3 bar). This discrepancy highlights the inherent complexity of numerically modeling blocked force states, where ideal rigid boundary constraints, localized hyperelastic stiffening under large strains and contact mechanics introduce artificial numerical rigidity that is difficult to replicate experimentally.

Overall, while the baseline FEA formulation proves highly reliable for predicting free bending kinematics, accurately capturing blocked force output presents clear limitations, motivating the advanced optimization steps discussed in subsequent sections.

The optimization process previously presented resulted in various sets of material parameters, one for each input combination of geometrical and operating parameters, such that the objective function equation (1) was minimized. Though the result obtained led to an overall improvement in the accuracy of SPA behavior prediction, two problems can be identified with the applied procedure:

• The optimization considered pressure as one of the conditions defining the behavior of the actuator. Theoretically, however, the material model should not be affected by operating conditions of the device. The independent optimization of material models of identical actuators under different pressures led to significant differences when comparing the material parameters of the actuators model with the same geometric parameters but operating at different pressures. A more refined optimization strategy should aim at the achievement of global optimization, using at the same time data acquired at different operating pressures.

• The second problem is that the optimization process considered a single measured point (i.e. the actuator tip) as ground truth to run the optimization/fitting problem between the simulated vs experimental data. Although this choice was justified by the tip’s importance for the actuator performances, it did not provide enough information considering that the study performed the optimization of an entire hyperelastic model, influencing the deformation of the entire actuator, using a single point. This could also be considered as one of the causes of the mismatch observed in the force data. By adding more reference points, the hyperelastic model would more accurately reflect reality. In addition, incorporating more points reduces the impact of any measurement errors made during the experimental procedure, which, although averaged by taking multiple measurements, could still be significant.

Therefore, it was decided to carry out an additional optimization to address these issues and obtain a single set of optimized material parameters for a given actuator model [known (N) and (T)] that would reproduce the behavior of the corresponding real actuator actuated at pressures of 1, 2 and 3 bar. Probes were used in the FEM model to track the displacement of multiple points of the simulated actuator, correspondingly to the actual points depicted in Figure 12, and of each probe the value of BE was calculated according to equation (1). The experimental displacement of each point was measured with respect to the starting (horizontal) positions previously explained in Chapter 4.1.

To perform a parallel optimization of three operating conditions for the actuator, a simulation consisting of several interconnected structural modules was set up, as visible in Figure 13.

Three structural static analysis modules were used, setting all simulations similarly, except for the pressure exerted, which was set to 1, 2 and 3 bar. The engineering data of the three simulations were linked so that all three simulations shared the same material parameters. A new function, used to calculate the optimal material parameters, was created:

where$BEN$ is the single objective function inherent to the N-th probe located at the base of the actuator; N corresponds to the number of humps. The simulations were run on a 9 bellows actuator and at pressures of 1, 2 and 3 bar; the function was calculated over a total of 27 points. The direct optimization process then takes as input the 27 values of $BEN$from the three static analyses and tests various combinations of the material parameters until the minimum of the new objective function [equation (3)] is identified within a maximum of 250 evaluations. The simulation ended due to resource limit: a material parameter set was identified, and the average error among all points of the three different simulations was 2.32 mm.

Considering only the error made on the actuator’s tip, from the comparison between this new set of material parameters and that obtained by optimizing only the position of the actuator tip there is an improvement in the prediction; in fact, the average BE committed using the parameters obtained from the optimization considering only the position of the actuator tip is 4.096 mm, while using this new set of material parameters the error is 3.701 mm, leading to a reduction of the error by 9.7%. The obtained improvement can also be noticed graphically by comparing images of real actuators with those of actuators simulated using the two different parameter sets Figure 14.

Indeed, it can be seen from the image that the model simulated in the first row succeeds better in simulating the behavior of its real version. For this type of actuator, (characterized by n = 9 and T = 1.6 mm), the parameters derived from this optimization result are shown in Table 9.

The study presented a comprehensive analysis of SPAs fabricated through FFF 3D printing, focusing on their design, modeling and experimental validation. The research highlighted the significant advantages of soft robotics, particularly their ability to safely interact with complex environments and living organisms, making them critical for applications in fields such as rehabilitation and wearable technology. A FEM model was devised to simulate the bending behavior of 3D-printed SPAs made from TPU material. The study meticulously assessed the influence of key design parameters, including wall thickness, the number of bellows and operating pressure, on the actuators’ performance. By systematically varying these parameters, the research aimed to understand their impact on the bending angle and force exerted by the actuators, leading to a total of 27 different configurations analyzed. To quantitatively evaluate factor significance, a full-factorial ANOVA was performed on the empirical data. The analysis revealed that operating pressure (⁠$P$⁠) is the dominant factor driving both bending angle and blocked force, accounting for over 68% and 71% of total variance, respectively. Wall thickness and bellows count also demonstrated statistically significant contributions (⁠$p<0.001$⁠), alongside strong two-way interaction effects such as $T×P$⁠. The experimental results confirmed the effectiveness of the FEM simulations, demonstrating that the proposed model could accurately predict the behavior of the SPAs under determinate loading conditions. The validation of the FEM model was achieved by comparing simulated data with experimental results, which revealed a close match. This baseline model was quantitatively evaluated using WLS regression for experimental data and OLS for ANSYS predictions. The comparative regression demonstrated strong physical consistency in unconstrained bending kinematics (similar experimental and numerical slopes from linear regression), while revealing that the model tends to overestimate blocked force at higher operating pressures (3 bar) due to ideal rigid constraints and localized hyperelastic stiffening. Once the set of material parameters that best simulated the real behavior of SPAs was identified, a further optimization process was carried out to reduce the BE in prediction: through an iterative process, the material parameters were changed until this error was minimized in nine different combinations of geometric/operational parameters. However, the optimization process is inherently dependent on pressure, leading to significant differences when comparing the material parameters of actuators with the same geometric parameters but operating at different pressures. In addition, the optimization process initially considered a single measured point as the ground truth, which did not provide enough information to allow the final deformed actuator’s shape to match the experimental one. An additional set of optimizations was set to address these issues, resulting in a single set of optimized material parameters for a given actuator model that could reproduce the behavior of the corresponding real actuator at various pressures. To address the single measured point problem, multiple reference points have been considered, which more accurately reflected reality and reduced the impact of measurement errors. This approach led to a significant improvement in the accuracy of the simulated behavior, with a mean reduction in BE of 9.7% compared to the material parameter obtained from single-point optimization. In conclusion, the integration of advanced modeling techniques with practical experimentation has yielded a robust understanding of the behavior of 3D-printed SPAs. This work demonstrates the potential for further exploration of soft robotics, particularly in developing more sophisticated actuators tailored to specific applications. The ongoing exploration of material behaviors and actuator designs will enhance the versatility and functionality of soft robotic systems, ultimately contributing to their broader adoption in real-world applications.