Modeling of soft fluidic actuators using fluid–structure interaction simulations with underwater applications

,


Introduction
Soft robots are made of highly deformable and compliant materials and usually follow a bioinspired design [1][2][3][4]. They show high dexterity, safety, physical resilience, and can adapt to delicate objects and environments due to their conformal deformation [2,5]. Their ability to elastically deform and adapt their shape to external constraints and obstacles make them ideal for applications such as gripping, bioinspired locomotion, minimally invasive surgery, wearable robotics, and implantable devices, where the environment is highly dynamic, sensitive to physical interaction, or constrained with restricted access [6][7][8][9].
The complex geometry of soft fluidic actuators and their hyperelastic material properties hinder the development of accurate mathematical models to describe their performance. Consequently, finite element modeling (FEM) has found many applications in soft robotics since it reduces cost and development time, and it can cope with the large deformations and material nonlinearities [16,45,46]. FEM has been used to evaluate soft actuator designs such as extending [47], fiber-reinforced [48,49], pneumatic network (PneuNet) [50][51][52][53], and omnidirectional actuators [54][55][56].
Although FEM studies have been widely employed for soft robotic modeling, further work is required to improve the modeling of environmental interactions, and to reduce the reality gap between simulations and experiments [17]. Custom physics-based and differential simulators with integrated gradient-based optimization algorithms have also been developed for soft robotic applications, such as ChainQueen [57], Elastica [58] and DiffPD [59], which allow for the co-design of robot geometry, materials, and closed-loop control strategies. Alternatively, mesh-free methods show strong potential to bridge this reality gap and have many advantages over mesh-based methods such as the FEM. Mesh-free methods include SPH (Smoothed Particle Hydrodynamics), MPS (Moving Particle Semi-implicit), MLPG (Meshless Local Petrov-Galerkin) and FVPM (Finite Volume Particle Method) [60,61]. SPH, for example, uses spatially distributed nodes, or particles, to represent matter (whether solid or fluid) [62,63]. However, unlike FEM, these nodes are not constrained by element connectivity, and these particles can flow and rearrange with the flow of material [64]. As a result, fluid interactions with deforming and moving solids can be handled naturally and without re-meshing [65], and complex geometries can be considered with little additional pre-processing.
Fluid-structure interaction (FSI) is the mutual interaction between a deformable solid body and an internal or surrounding fluid flow, where the flow has a strong impact on the structure, and vice versa. This fluid flow exerts hydrodynamic forces that deform the structure, and in turn modify the boundary conditions on the fluid domain, since the deforming structure imparts different velocity and physical constraints to the fluid [66]. FSI research has become a cornerstone in many different engineering applications such as aerodynamics, hydrodynamics, and structural analysis [67]. FSI analysis has also been shown to provide a method of evaluating the performance of flapping foils [68,69], and various flexible underwater and aerial robots [70][71][72]. In [73,74], FSI analysis with ANSYS was used to study the two-phase flow hydrodynamics and aquatic propulsion of a bionic flipper with a nonuniformly distributed stiffness along the surface [73,74]. FSI analysis of flexible caudal fins was performed in [75,76] by combining their in-house fluid solver with an FEM-based code CalculiX via preCICE, where flexible fins were shown to outperform their rigid counterparts in terms of thrust generation and propulsion efficiency. FSI analysis of flexible caudal fins was also performed in [77] using ADINA for twodimensional fish models with various combinations of rigid and flexible links. ADINA was also used to perform FSI simulations for the fin motion and thrust generation of a cownose ray undergoing oscillating and flapping behaviors in [78].
Due to their mesh-free Lagrangian nature, particle methods have been used to simulate free-surface fluid flows and their interactions in various ocean engineering applications, including fluid-structure interaction [79]. In [80], a benchmark setup consisting of a dam-break with an elastic plate was designed and experiments were conducted to investigate the phenomenon of free surface flow impacting on elastic structures. In [61], the accuracy, stability, convergence and conservation properties of a Hamiltonian MPS structure model was demonstrated using four benchmark tests. For a complete list of benchmarks tests and associated validation studies for hydroelastic FSI studies, the reader is referred to [60,81]. Additionally, SPH has been used for human movement simulation in sports such as swimming, diving, and kayaking in which robotics algorithms are used to represent skeletal motions [82][83][84].
Although the aforementioned works have used FSI to analyze the motion of benchmark test cases and flexible robots, they have not simulated the dynamics of fluid-driven actuation. In fact, very few works have employed FSI to model soft fluidic actuation. Instead, the vast majority of works using the FEM approach in soft robotics apply a uniform pressure boundary condition to the internal cavities of the soft actuator in quasi-static simulations. The assumption of uniform, steady-state internal pressures may limit the applicability of some quasi-static simulations. In reality, the pressurization of physical soft actuators is achieved by supplying flow into the inlet of the internal chambers from fluidic sources such as syringe pumps [85,86], pressure-regulated sources or flow sources [87,88]. In particular, pressure-regulated sources are a popular alternative since they provide larger system capacity, fast actuation and improved energy efficiency [89,90].
Two-way FSI analysis shows great potential for a more realistic method for modeling of the pressurization of soft actuators. For example, it allows the investigation of fluid flow and pressurization rate on the performance of soft actuators, which can assist with analysis and optimization of dynamic characteristics [17,91]. While one-way FSI analysis can be used for small deformation problems, the modeling of soft fluidic actuators requires two-way FSI simulations since both fluid and solid domains undergo large deformations. The meshless local Petrov-Galerkin method was used to perform two-way FSI analysis of a worm soft robot in [92]. In [91], COMSOL Multiphysics was used to perform two-way FSI simulations for PneuNet bending actuators using a time-dependent study under the assumption of incompressible and laminar flow. The FSI results were compared to static finite element simulations using Abaqus, where FSI simulations more accurately described the soft actuator motion for high pressurization rates.
In this work, fully-coupled, three-dimensional FSI simulations are applied to the modeling of soft fluidic actuators and robots. A mixed steady-transient approach is proposed, which significantly reduces the solution time. This approach is extended to model fluid flow both in the internal chambers and around the external walls of soft actuators in underwater applications. A number of case studies are presented including bending actuators, a soft robotic fish fin for propulsion, and experimental results of a bending actuator in a high-speed fluid. Two levels of complexity are considered: This article is organized as follows. Section 2 presents the governing equations for the mechanical and fluid domain and introduces the procedure for the FSI simulations. In Section 3, the proposed FSI approach is used to evaluate the internal pressure and flow dynamics of soft fluidic actuators. In Section 4, the FSI simulations are extended to underwater applications. In Section 5, the experimental platform is presented and the experimental results are compared to the FSI predictions. Finally, Section 6 concludes this article and outlines the future work.

Materials and methods
This section presents the background required for the proposed FSI simulations. Firstly, the hyperelastic material models describing the mechanical behavior of soft actuators fabricated using silicone rubbers and filaments are presented. Then, the equations for the fluid domain simulations are introduced. Finally, the methodology for the implementation of the FSI simulations is described, including boundary conditions and solver settings.

Hyperelastic material models
Soft fluidic actuators experience large deformations during pressurization. Therefore, the silicone rubbers used in their fabrication are described using hyperelastic material models [93][94][95]. In this work, the hyperelastic materials are considered to be isotropic and incompressible. Furthermore, all time-dependent properties of soft materials such as hysteresis, viscoelasticity, stress-softening and Mullins effect are neglected. These assumptions are reasonable for a wide class of soft fluidic actuators and largely employed in soft robotics [16,96,97], especially considering lower levels of volume change due to ballooning for fiber-reinforced actuators and fast pneumatic network actuators [98,99]. In addition, although for some soft materials the viscoelastic parameters play a larger role in the material behavior, for common silicone rubbers in soft robotics, the viscoelastic effects are minimal in comparison to the factors arising from the fabrication process and environment, including curing temperature, mixing ratio, and degassing pressure [46].
The stretch ratios represent the deformation of a differential cubic volume element along the principal axes of a Cartesian coordinate system [95,100]. They are given by the ratio of deformed length ( ) to undeformed length ( ), i.e., Using the principal stretches, the principal invariants are defined as A hyperelastic material model relies on the definition of a strain energy function (or strain energy density), which is the amount of energy stored elastically in a unit volume of material [93,94,101]. For an isotropic incompressible hyperelastic material, 3 = 1 and the stressstretch relations are obtained from the strain energy function by virtual work considerations: where is the principal Cauchy stress and is the hydrostatic pressure (indeterminate Lagrange multiplier), which can be determined from the equilibrium equations and boundary conditions [16,102,103].
The following hyperelastic models are considered in this work: (1) Yeoh model: (2) Polynomial model (generalized Rivlin model): The constant parameters 1 , 2 and are determined via material testing and curve fitting [16,104]. In this work, the silicone rubber Elastosil M4601 (28 A shore hardness and 1.13 g/cm 3 density) is used for the bending actuators. For this material, the Yeoh model is employed with 1 = 0.11, 2 = 0.02 and 3 = 0 MPa [49,50,105]. For the case study on a soft robotic fish, the TPU filament NinjaFlex (85 A shore hardness and 1.19 g/cm 3 density) is used. For this material, a fiveparameter Rivlin model is considered with 10

Governing equations for fluid domain
The equations governing the flow simulations in Fluent are the three-dimensional viscous incompressible Navier-Stokes equations [70,106] = 0 (6) where are the velocity components, is the pressure, and and are the fluid density and kinematic viscosity. The first equation is the continuity equation and the second equation describes the momentum balance.
At the boundary between the fluid and solid domains, the following equation applies [73] − where is the fluid effective shear stress tensor, is the Cauchy stress tensor and is the normal vector.

Fluid-structure interaction simulations
The soft fluidic actuators used in this work are modeled using Autodesk Inventor (Autodesk Inc.). In contrast to conventional finite element simulations for soft fluidic actuators, the actuator chambers in FSI simulations are filled with a fluid body which has the exact geometry of the cavity. The solid and fluid body geometries are then imported into ANSYS Workbench.
In the mechanical domain, a quasi-static assumption is made and simulations are performed using Static Structural Analysis in ANSYS Mechanical. Meshing is performed on the solid component by suppressing the geometries representing the fluid bodies and specifying an unstructured, tetrahedral mesh with a nonlinear mechanical physics preference and size of 2-5 mm. Analysis controls are specified as follows: auto time stepping is disabled, time steps are manually set to 0.01-0.05 s, and the large deflection setting is enabled. The time steps are selected according to the input inlet flow velocities and corresponding deformation levels, with smaller steps used for larger inlet flow velocities. This range is sufficiently refined to capture the soft actuator deformation and to provide converged results without the extreme computational overhead required for smaller time steps [16,46]. The boundary conditions of the FSI simulations are shown in Fig. 1a. A fixed support is applied to the face of the actuator inlet. The internal wall of the solid body is defined as the fluid-solid interface. An additional fluid-solid interface is introduced at the external walls of the solid body to simulate the actuator submerged under water (Fig. 1b).
A transient formulation is employed for the fluid solution and is applied using Fluent. The fluid body is meshed by first suppressing the solid body geometry and specifying an unstructured, tetrahedral mesh with a CFD physics preference and size of 2-5 mm. A viscous laminar model is assumed and the fluid material properties are assigned to represent water at ambient conditions. The boundary conditions are as follows (Fig. 1a): (1) the inlet is defined as a ''pressure-inlet'' to model a constant pressure inflow from pressure-regulated sources, and (2) the faces of the fluid inside the chambers (internal fluid walls) are defined as ''wall''. FSI between the solid and fluid domains is enabled by implementing a dynamic interface mesh that is created with the ''smoothing mesh'' method and the internal (fluid-side) walls are coupled to the fluid solution using the ''system coupling'' setting. The solution employs a pressure-based solver, pressure-velocity coupling, a second-order implicit transient formulation and hybrid initialization. The governing equations are discretized using a second-order interpolation scheme, which provides more accurate results for tetrahedral meshes and flows not aligned with mesh. Gradients are determined using the ''least squares cell-based'' method, where the solution is  Fig. 1b). In addition, the external fluid walls are set to system coupling.
The interaction between ANSYS Mechanical and Fluent is defined using System Coupling, a component system in ANSYS Workbench. A two-way data transfer is created between the internal cavity of the soft actuator (fluid solid interface 1) and the internal fluid walls. For underwater simulations, an additional two-way data transfer is created between the external walls of soft actuator (fluid solid interface 2) and the external fluid walls.
Through system coupling, the mechanical and fluid physics are solved separately and then coupled sequentially or simultaneously until equilibrium is reached. Forces and stresses from the fluid side of the interface are mapped to the solid side, and displacements and velocities from the solid side are mapped to the fluid side, as shown in Fig. 2. System coupling also ensures that time step settings are consistent across the Ansys Mechanical and Fluent solvers. The system coupling settings are as follows: (1) coupling end time: 1 s, (2) coupling step size: 0.01-0.05 s, and (3) maximum number of iterations per coupling step: 10.

Fully-coupled 3D FSI simulations
In this section, FSI simulations are presented considering internal FSI only. Here, the transient fluid simulations are used to represent the pressurization of the internal cavities of the soft actuator but no external fluid is considered. This section considers the pneumatic network As discussed in the Introduction, the vast majority of works on modeling of soft fluidic actuators using FEM employ a uniform pressure applied to the internal cavities [16,45], which does not represent the physical pressurization of the actuator using practical fluidic sources. In comparison to conventional FEM simulations, the FSI results display a nonuniform pressure during actuation, as shown in Figs. 3c and 3d. In particular, the internal pressure distribution displays larger pressures near the flow inlet, which may be used as a justification towards designing soft actuator with increased wall thickness around the pressure input [91].
A comparison between FSI and conventional FEM results for the bending angle of the pneumatic network actuator is presented in Fig. 4. The bending angle is measured between the initial and current vectors extending from the center of the base of the actuator and the center of the tip of the actuator [107]. In general, FEM results display slightly larger bending angles compared to FSI, which is likely to be associated with the nonuniform pressure predicted by the FSI simulations. This difference, however, is within tolerance associated with the uniaxial tensile test data, material processing, fabrication procedures and environment conditions.  In addition to the internal pressure distribution, the FSI simulations allow for the evaluation of the internal fluid flow in the chambers of the soft actuator. Figs. 3e and 3f show fluid streamlines colored by fluid speed. Higher fluid speeds are seen in the central channel as the fluid is able to accelerate under the applied inlet pressure without significant impedance from the internal walls of the actuator. In contrast, the fluid that enters the first few chambers is decelerated strongly by near perpendicular collision with the internal walls of the chamber. Fluid speed is predicted to be nonlinear, inversely proportional to calculated internal pressure and highly dependent on the shape of the actuator. The results from this analysis can be used to optimize the internal routing for improved deformation, response time and fatigue performance. The simulation time for an input pressure of 10 kPa and step size of 0.02 s was 22 minutes using 4 cores at 3.6 GHz.

Underwater simulations
In this section, the FSI simulations are extended to consider both internal and external fluid flow. The internal FSI component models the pressurization of the internal chambers of the soft actuator, while the external FSI component models the interaction between the external walls of the soft actuator and the underwater environment. Here, FSI simulations are used to describe a case study using a fish-inspired soft robot flapping in a fluid environment.
FSI simulations are not limited to modeling the internal flow used for pressurization of the soft actuator. Here, as an example application, a fish-inspired soft robot with a pneumatic network design is considered with 80 mm height, 28 mm width, 80 mm length. The passive fin has a 40 mm length and 4 mm thickness. The behavior of the soft fish in underwater applications is evaluated by inserting the soft robot model in a representation of a body of water which models a filled tank with an opening at the top. Preliminary simulations using Elastosil M4601 for the soft robotic fish revealed excessive deformation of the passive fin, which could be addressed by using a harder material, a multimaterial design, or a thick passive fin. Here, NinjaFlex was selected for the fish-inspired robot due to its larger hardness, in addition to its wide popularity in soft robotic applications. The use of a harder material results in increased fin stiffness and reduced deformation under external fluid flow, which consequently improves convergence of the underwater simulations.
The pressure and velocity contours for the water in the tank are shown in Fig. 5. As the internal pressure is increased within the chambers on either side of the soft robot, the bending deformation increases. This in turn increases the magnitude of the pressure and speed applied to the external fluid. Regions of high speed and pressure are created in the external fluid, which are focused at and beyond the fin. The reaction force from the fluid can be determined and used to predict propulsive action. These results can be further evaluated, for example, to design bioinspired underwater soft robots with improved propulsion performance and manoeuvrability, which is a topic for future research.

Experimental validation
In this section, a series of comparison experiments and FSI analysis are presented and discussed. First, the fabrication procedure for the soft fluidic actuator is presented, followed by a description of the custom experimental setup used to validate the FSI simulation results. Then, simulated and experimental deformations are compared for a soft bending actuator in an underwater environment where fluid flow is applied to the external walls of the actuator in both unpressurized and pressurized conditions.

Fabrication of the bending actuators
The actuator used in the experimental results is a fast pneumatic network bending actuator [50], which was fabricated using standard molding procedures [41,42]. Molds were designed in Autodesk Inventor and printed using a Flashforge Inventor (Flashforge Corporation, USA). Elastosil M4601 was mixed at the recommended ratio and degassed using a vacuum pump. The mixture was poured into the mold, and left to cure for 12 hours at room temperature. Then, the actuator was removed from the mold and a bottom layer of 3 mm thickness completed the actuator, which has a length of 100 mm and height of 23 mm. The geometrical features of the fabricated bending actuator is shown in Fig. 6c.

Experimental setup
The experimental setup is shown in Fig. 6a, which includes a water tank, sliding mechanism, hanger, pulley system and video-camera. The dimensions of the water tank are 121 cm (length) × 51 cm (width) × 50 cm (height). During the experiments, water was filled to a depth of 25 cm. The soft actuator is submerged at a depth of approximately 12 cm, ensuring that the surface and wall effects are small.
The locomotion of the actuator was restricted to the negativedirection using a platform consisting of two linear shafts connected to low-friction ball bearings on a 3D-printed slider. The linear shafts are connected to 3D-printed supports on both sides of the tank. A 3Dprinted hanger was used to secure the actuator under water at a fixed depth and is connected to the slider using threaded rods. The threaded rods are fixed using washers and nuts on both sides of the hanger and slider. All 3D-printed components were fabricated using a Flashforge Inventor printer.
The pulley system consists of a mass block supported by a string threaded through the pulley and fixed to the slider. To simulate water flow, 0.1 kg, 0.2 kg and 0.3 kg mass blocks were used. The average slider velocities, measured under four experiments, are 0.3 m∕s, 0.5 m∕s and 0.7 m∕s, respectively.
The internal pressure of the soft actuator is regulated using a 60 mL syringe pump [86], with an on/off controller implemented using an Arduino Due and pressure sensor (SEN0257). Image acquisition is performed at 60 frames per second using a webcam (Logitech C920) mounted on the slider.

Comparison between FSI and experimental results
In this section, FSI simulations are developed to evaluate the behavior of the soft actuator for underwater applications under the presence of fluid flow around the external walls of the actuator. In comparison to Fig. 1, the boundary conditions of the tank are modified as follows: (1) the left side wall of the tank becomes a ''velocity inlet'', (2) the right  side wall of the tank becomes a ''pressure outlet'', and (3) the top wall, which was previously a ''pressure outlet'', is assigned a ''wall'' boundary condition. The final boundary conditions are shown in Fig. 6b.
The simulated and experimental deformations are compared with constant flow when the actuator is either unpressurized or pressurized. The actuated pressure level was chosen to result in a bending angle of 45 • under static flow. Using FSI, this was predicted to require 35 kPa, and was experimentally found to be 33 kPa. The difference may be due to the uniaxial test data used to calculate the material constants of the silicone rubber, or dimensional tolerance of the fabrication process.
Instead of simulating the entire water tank, which would require a large number of mesh elements and significantly increase simulation time, a fluid body around the actuator is considered with dimensions 20 cm (length) × 18 cm (width) × 10 cm (height). This simplification is employed to reduce the computational effort and does not significantly affect the accuracy of the FSI simulations, as demonstrated in the following discussion.
In the FSI simulations, the inlet fluid velocity is matched to the experimental average slider velocities for each of the mass blocks, i.e., 0.3 m∕s, 0.5 m∕s and 0.7 m∕s. The experimental results and simulation predictions for pressurized actuators are shown in Fig. 7 and summarized in Fig. 8. Fig. 7a shows the actuator pressurized to a 45 • deformation without external fluid flow. When the fluid flow is increased from 0.3 m∕s to 0.7 m∕s in Figs. 7b to d, the fluid flow opposes the bending due to pressurization, eventually leading to backward bending and nonuniform curvature in Fig. 7d. As shown in Fig. 8, the experiments show good correlation between the simulated predictions and experimental results, which validates the reliability of the proposed FSI simulations. The small discrepancies may be associated with the uniaxial testing data used to characterize the hyperelastic model of the silicone rubber and the visual blockage of the base of the soft actuator in the experiments, which hinders the selection of an accurate measurement point for the bending angle. The evolution of the experimental deformations and velocity contours for selected cases are included in Supplementary Video S1.
The pressure and velocity contours at a cross-section of the water tank are shown in Fig. 9 to illustrate the flow solution around the actuator. The upstream velocity (left side) shows the inlet flow for each the three conditions. The velocity demonstrates three distinct flow characteristics of deceleration, acceleration, and diffusion. As the flow progresses in the -direction, it diffuses and recovers partially towards the free-stream condition near the exit boundary (on the right), as shown in Fig. 9a.
The pressure plot in Fig. 9 shows that regions of positive and negative pressure develop around the soft actuator. The velocity drops as the fluid approaches the soft actuator, reaching zero velocity and maximum pressure at the bottom layer of the soft actuator. Negative pressure regions can be seen near the fins and bottom of the actuator due to changes in the direction of the fluid flow as it travels around the actuator in the -direction. In addition, as the flow accelerates around the top and bottom edges of the actuator, the pressure drops. Fig. 10 shows the experimental and FSI simulation results for the unpressurized actuator with fluid velocities of 0.3 m∕s and 0.5 m∕s. In both cases, the fluid forces acting on the bending actuator cause it to bend backward. The pressure and flow dynamics in Fig. 10 are similar to the results in Fig. 9. Fluid deceleration, acceleration, and diffusion are observed.
Minor discrepancies between the experimental and simulation results can be observed in Fig. 10b. The differences are due to contact between the actuator chambers closest to the base, which occurs at velocities exceeding 0.4 m∕s. Contact mechanics are not presently accounted for in the simulation. Although contact walls can be added to the mechanical model of the bending actuator, this is not considered in the fluid domain. Resolving this issue is a topic of current research.

Conclusions
This article describes a computationally efficient method for threedimensional fluid-structure simulation of fluid-driven soft actuators and robots. The FSI simulations use a mixed steady-transient formulation, where the mechanical domain employs a quasi-static assumption and the fluid domain uses a transient formulation. Compared to conventional FEM simulations, where a uniform pressure is applied to the internal cavities of the soft actuator, the proposed FSI approach employs a pressure inlet which represents a practical scenario using pressure-regulated sources for soft robotic actuation. Further, the FSI simulations allow the evaluation of the internal pressure and flow dynamics.
The proposed approach is extended to simulate the behavior of soft actuators in underwater applications, where external fluid flow may exist. Although previous works have reported the use of FEM and FSI simulations to predict the deformation of soft actuators under internal pressurization, this work reports the first application of both internal and external fluid-structure interactions, which allow for the modeling of soft robots and actuators in underwater environments. Case studies include a bioinspired fish fin for propulsion, and a bending actuator exposed to significant external fluid flow. An experimental system was constructed to compare the simulated and measured deformation of bending actuators exposed to a range of external fluid velocities. The simulation results are observed to correlate closely with the experimental measurements.  The methods proposed in this work yield important insights into the behavior of soft actuators and robots. For example, the fluid pressure and velocity distributions can be used to design and optimize underwater soft robots, such as the soft robotic fish fin presented in the case study. Lift and drag forces can also be evaluated to determine the thrust performance. In future work, the proposed FSI approach will be extended to evaluate wing designs of a soft ray-inspired robot in order to obtain desired manoeuvrability properties. In addition, the impacts of fast pressurization cycles on the overall motion and pressure dynamics of soft fluidic actuators, such as vibration, will be studied.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data will be made available on request.

Appendix A. Supplementary data
Supplementary material related to this article can be found online at https://doi.org/10.1016/j.ijmecsci.2023.108437.