A simple, fish-like geometry was constructed representing the spine and fin of the animal, with a core plastic sheet (Plastic Shim Stock™  0.5 mm, Artus). A silicone fast-PneuNet actuator molded with uncured elastomer (Dragon Skin™ 20, Smooth-On) (Mosadegh et al., 2014), 10 cm in length and 2 cm in height, was attached with a cured silicone sealant (ELASTOSIL™ E43, Wacker) to both sides of the sheet, representing the muscles of the fish for lateral undulation motion. The segmented and chambered soft actuator allows for a smooth bending pattern, similar to the undulation motion of fish (Figure 1A; Figure 2A).

The fabrication process for the soft strain sensors is shown in Figure 1B, based on the work from (Park et al., 2010). For the first layer with the microchannel pattern, uncured silicone rubber (EcoFlex™ 0030, elastic modulus E = 125 kPa) is poured into a 3D printed (Objet30™  Objet Geometries Ltd.) mold. For the second layer, uncured liquid silicone is coated on a flat substrate with a thickness of 1 mm. Both layers are than cured for an hour at 60°C in a convection oven. The flat layer is then spin-coated (1,000 rpm for 60 s) with uncured silicone and partially cured at room temperature for 10 min. Then, the layer with the microchannel pattern is gently placed on the flat layer, and the combined structure is cured again.

After bonding, conductive liquid metal eutectic Gallium-Indium (eGaIn, Sigma-Aldrich) is injected into the microchannels using two syringes, one for injection and the other to release the air captured in the channels during fabrication. Two wires are connected to both ends of the microchannel via the syringe holes. Finally, a silicone glue (Sil-Poxy, SmoothOn™) is applied and cured around the holes into which the wires are inserted, firmly holding the wires and sealing the holes.

Changes in the cross-sectional region of the eGaIn-filled channels result in changes in electrical resistance, which can be used to measure the bending amplitude of the soft robot. The sensors have low temperature sensitivity (Vogt et al., 2013) and the liquid’s incompressibility makes them ideal for water experiments with change of ambient pressure (Hellebrekers et al., 2018; Schwab et al., 2021). For calibration testing, the sensor was bent in the undulation range of −50 to 50°, and during this angle sweep the voltage was sampled. In comparison to the ground truth kinematics obtained from videography, the sensor response was linear (R2 = 0.952, Figure 2B), with a relative error well described by Gaussian noise with a standard deviation of 0.4% (Figure 2C).

In this work, hook-and-loop fasteners instead of embedded fabric like in (Mengüç et al., 2014) were used, for quick adjustment and replacement of the sensors onto the soft actuators.

To characterize swimming robots, various models for undulatory locomotion have been proposed (Salumäe et al., 2012; Bliss et al., 2013). In this paper, a data-driven, lumped parameter model was used, which is a popular method when modeling compliant robots (Sadler and Sandor, 1973; Wang and Huston, 1994; Nishikawa et al., 2015). The soft silicone robot was treated as a chain of rigid elements, which are connected by hinges with constant stiffness and damiping coefficients. The parameter of the hinges were then optimized with genetic algorithm using experimental ground-truth kinematics (Lin et al., 2021). Compared to finite element simulations (Allard et al., 2007) or machine learning tools (Gillespie et al., 2018), this data-driven approach is computationally less demanding while still managing to capture the behavior of the soft robotic fish (Lin et al., 2021).

The model uses the pressure as input and takes the co-contraction effects as well as hydrodynamic forces into account, successfully predicting the angle and position of the soft robotic fish and simulating the ground-truth behavior. This enables a faster testing of a controller, without the need of constant experiments for verification. In the future, the model will hopefully accelerate more sophisticated control designs and guide the further development of soft swimming robots.

The amplitude control system is designed in Simulink (MATLAB™ R2020b) and a flow chart is shown in Figure 2E. A proportional-integral-derivative (PID) controller is used for pressure adjustments in response to the measured curvature. Fish can be observed to adjust their undulation frequency to different flow conditions, while their amplitude remains constant (Tytell and Lauder, 2004). The goal of the controller is therefore to maintain a specific amplitude across a range of undulation frequencies and flow speeds, by controlling the pressure. The verification of the controller performance was previously done in a static tank (Lin et al., 2021) (Figure 2A) and we now evaluate the output of such a controller in a flow tank with varying flow rates.

Figure 2D shows the components for the controller and the experimental setup. Through a digital pressure regulator (ITV0050-3BS, SMC), 2.5 bar compressed air is supplied to the system. The PneuNet actuators are controlled by two directional solenoid valves (SYJ7320-5LOU-01F-Q, SMC) respectively. A micro-controller (ATmega328 on Arduino Uno, Arduino) sends the control signals to the valves as well as to the pressure regulator. The robotic platform is mounted on a 150 mm long, 20 × 30 mm rectangular aluminum section with a 2 mm wall thickness. Currently, the controller’s response is limited by the use of peak amplitude, which suggests that the control input is changed only twice per oscillation cycle.

Particle Image Velocimetry (PIV) (Willert and Gharib, 1991) is an established method used to study flow phenomena (Sfakiotakis et al., 1999) and is well adapted for investigating the physics surrounding the soft robotic fish platform discussed in this paper. Experiments were conducted at Swiss Federal Laboratories for Material Science large scale water flume (Engineering Laboratory Design, Inc.). This facility has a cross section of 0.6 × 1 m² and a fully transparent, 6 m long test section (Figure 3A). The system is driven by a variable speed pump, allowing precise control of the water bulk speed from 0.02 m/s to 1.5 m/s. Controlled velocity profiles at the inlet and outlet of the test section are ensured by a 6:1 inlet contraction and outlet diffuser with guide vanes.

The tail motion and velocity fields were captured using a PIV system that consists of a pulsed Nd:YAG laser, double cavity, with 100 mJ per pulse at 532 nm and a maximum repetition rate of 100 Hz and two 5.5 Mpx SCMOS cameras. The laser beam was guided through a laser arm to the sheet optics, allowing to illuminate a cross section of the robotic fish and the PIV tracers on a horizontal plane at 1/3 the foil height from the bottom. The PIV cameras were positioned underneath the water flume floor and calibrated simultaneously with a 3D calibration target. The cameras were run in double frame mode at 24 fps and equipped with 50 mm focal-length lenses achieving a spatial resolution of 10 px/mm. The combined field of view of both cameras covered an area of 400, ×, 200 mm, aligned with the flow direction. The field of view of Camera 1 was centered on the robotic fish, whereas Camera 2 was used to retrieve the velocity field at the trailing edge of the flexible foil. The origin of the coordinate system was placed at the leading edge of the robotic fish.

The position of the tail of the robotic fish was retrieved by segmentation and thresholding of the recorded images with an algorithm developed in MATLAB. Thanks to short pulsed laser illumination (5 ns pulse duration) the images showed little or no blurring, allowing us to reconstruct the position of the tail with an accuracy of ± 1 mm.

For the PIV measurements, spherical, glass seeding particles of 10 microns were used, while the laser sheet thickness was about 2 mm, in view of the strong 3D structure of velocity field. The images were stitched side-by-side and processed with a multi-grid approach having final interrogation window size of 32 × 32 Px and 50% overlap, yielding to a physical resolution of 0.625 velocity vectors/mm. For each trial, 1,000 images have been recorded to enable the calculation of the mean and rms velocity field, as well as Dynamic Mode Decomposition (DMD) (Schmid, 2010) of the velocity field.

The closed loop controller of the soft robotic fish has been tested at different flow conditions. First the water flow speed has been adjusted step-wise from 0 c m s to 19.4 c m s , while the actuation frequency was kept constant at 0.55 Hz. Afterwards, the actuation frequency has be varied step-wise from 0.3 to 0.7 Hz, while the water flow speed has been left constant at 5.5 c m s . The exact test parameters of the water flow speed sweep and of the frequency sweep are given in Table 1 and 2, respectively. Before varying the flow speed and frequency a reference trial (No. 1 and 7) at 0 c m s and a frequency of 0.55 Hz have been performed.

Figure 4 shows the sensor reading performance for two exemplary trials: trial 8 was performed at a undulation frequency of 0.3 Hz and a constant flow velocity of 5 c m sec (Figures 4A,C), while trial 11 (Figures 4B,D) was conducted at a frequency of 0.7 Hz. By performing a Discrete Fourier transform of the measurements of the undulatory soft robot, we observed rather high second harmonics in at 0.3 Hz (Figure 4A) compared to 0.6 Hz (Figure 4B). Those second harmonics are hypothesized to be introduced by unwanted structural interactions between the robot’s center plastic sheet and the actuators.

We analyzed the average flow speed around the swimming fish by PIV experiments. In the time averaged flow field, a clear upside V-shape pattern can be observed (Figures 5A,B), with two distinct jets originating from the tail motion. The angle enclosed by the jets varies in the different trials.

By plotting the angle enclosed by the jets as a function of the different flow parameters, we can illustrate how they affect the time-average features of the flow field associated to the tail motion. Figure 5C depicts the flow angle and displacement as a function of the water flow speed in the recirculating flow tank. The closed loop control is able to adjust the power input to the water flow speed, maintaining a rather constant displacement up to a water flow speed of 15 c m s .

Figure 5D compares the performance of the robotic platform in feed forward and closed loop control over a range of frequencies. We can observe that controller keeps the amplitude of the soft fish approximately constant while in feed forward control, the amplitude drops significantly with increasing frequency. This illustrates the potential of the controller to mimic the swimming behavior of fish, adapting the parameters to keep a constant tail-tip amplitude in changing environments.

Although the actual flow field exhibits a 3D structure, an approximate estimate of the wake power (the net energy flux) can be obtained from the 2D PIV measurements. Previous experimental work using flapping foils (Triantafyllou et al., 2000) and swimming live fish (Fish and Lauder, 2005) demonstrates that the Strouhal number, a fundamental dimensionless parameter, dominates the hydrodynamic performance of fish locomotion, because it correlates with vortex shedding dynamics (Eloy, 2012; van Leeuwen et al., 2015). S t = f A t a i l v s w i m , where v _swim is the mean swimming speed during undulatory movements, A _tail is the peak-to-peak amplitude of the tail tip, and f the frequency.

Studies suggest, that the optimal number is in general optimal between 0.1 and 0.55 St, (0.25–0.35 for carangiform swimmers and 0.4–0.5 for anguilliform swimmers) (Nudds et al., 2014; Taylor et al., 2003; Floryan et al., 2018).

Assuming a steady-state, uniform velocity profile upstream and a 2D flow field of thickness equal to the height of the robot, the wake power was estimated by integrating the kinetic energy measured by PIV analysis (TytellSchultz and Webb, 2002; 2007). We computed the kinetic energy flux under 2D assumption from the PIV averaged flow fields on a control volume surrounding the swimmer. In this case the upstream input can be expressed as 1 2 ρ h w U 3 , where U is the flow speed, h the height of the robot, and w the width of the wake. Therefore, the net wake power can den be expressed by P wake ′ = 1 2 ρ U ∫ U + u 2 + v 2 + w 2 − U 2 d S where u, v and w are the axial, lateral and vertical fluid velocities produced by the robot, respectively. The power coefficient C _p is obtained by normalizing the wake power by 1 2 ρ 2 h L U 3 , where L is the length of the swimmer ((Schultz and Webb, 2002; Tytell and Lauder, 2004; Krueger, 2006)).

In the present case we measured a C _p close to zero at a water speed of roughly 15 c m s , which would indicate an estimate of the maximum robotic fish swimming speed. At 5.3 c m s , we measured an optimal actuation frequency of 0.6 Hz. In Figure 6, we additionally plotted the results from open-loop experiments with a force sensor attached to the robotic fish platform (Jusufi et al., 2017). We can observe that the results from the PIV experiments agree qualitatively with the results from previous experiments. The slight differences might be explained, besides the limitations of this approach, by the new controller implemented in this study, which leads to a slightly different optimal frequency (0.6 Hz compared to 0.8). Furthermore, the wake power model estimates an higher self propelled speed as the C _p equal zero at roughly 15 c m s , while the previously measured thrust equals zero at 13 c m s .

In Figure 6C, it can be observed that as the St number increases, the thrust as well as the C _p gradually increase from negative to positive values.

It should be stressed that this analysis of 2D data is based on simplistic assumptions and provides only a rough estimate of the wake power and the power coefficient, providing just an indication of the expected swimming speed and the optimal actuation frequency. More accurate estimates require measuring the full 3D flow field. Nevertheless, these results show that a thrust assessment is possible without attaching sensors. We envision to use this method with a range of swimmers like e.g., eels, bluegills, or trouts.

Dynamic Mode Decomposition (DMD) is an emerging purely data-driven technique that provides linearly reduced order models for high-dimensional, complex systems (SchmidRowley et al., 2009; 2010; Tu et al., 2014; Kutz et al., 2016). Coupled spatio-temporal coherent patterns or modes can be extracted from the observed data. DMD has applications in fluid dynamics, neuroscience, robotics, or disease modeling.

As the swimming locomotion is periodic, DMD analysis can be useful to extract the dominant space-time modes of the flow field. As an example, Figure 7A shows the real and imaginary parts of the eigenvalues associated to each spatial mode for an undulation frequency of 0.55 Hz and a free stream velocity of 9.3 c m s . Besides the zero-frequency mode 7, which represents the average velocity field, we can identify modes 1, 3, and 5 as well as theirs conjugates 2, 4, and 6 that dominate the instantaneous flow field as in view of their amplitude.

Mode 1 (Figure 7B) has a frequency equal to 1.7 Hz, roughly three times the undulation frequency. This mode is associated to the presence of a vortex centered at x = 0, y = − 260 mm and shed by the tail motion at a frequency equal to three times the swimmer undulation frequency. Mode 3 (Figure 7C) has a frequency equal to 1.1 Hz or twice the undulation frequency. As the vertical velocity contour map is symmetric, we expect this mode to be associated to an instantaneous thrust that oscillates periodically between negative and positive values, influencing the instantaneous propulsion. Mode 5 (Figure 7D) has a frequency equal to the undulation frequency (0.55 Hz), and it is directly related to the oscillatory motion of the tail. As the vertical velocity contour map is anti symmetric, we do not expect this mode to be directly linked to a positive or negative trust along the stream wise direction. On the contrary, the lateral component of the thrust vector should be non-zero. We observed these modes also in the other trials, performed with different operating parameters.