Analysis of Robot–Environment Interaction Modes in Anguilliform Locomotion of a New Soft Eel Robot

Abstract

Bio-inspired robots with elongated anatomy, like eels, are studied to discover anguilliform swimming principles and improve the robots’ locomotion accordingly. Soft continuum robots replicate animal–environment physics better than noncompliant, rigid, multi-body eel robots. In this study, a slender soft robot was designed and tested in an actual swimming experiment in a still-water tank. The robot employs soft pneumatic muscles laterally connected to a flexible backbone and activated with a rhythmic input. The position of seven markers mounted on the robot’s backbone was recorded using QualiSys^® Tracking Manager (QTM) 1.6.0.1. The system was modeled as a fully coupled fluid–solid interaction (FSI) system using COMSOL Multiphysics^® 6.1. Further data postprocessing and analysis were conducted, proposing a new mode decomposition algorithm using simulation data. Experiments show the success of swimming with a velocity of 28 mm/s and at a frequency of 0.9 Hz. The mode analysis allowed the modeling and explanation of the fluctuation. Results disclose the presence of traveling waves related to anguilliform waves obtained by the superposition of two main modes. The similarities of the results with natural anguilliform locomotion are discussed. It is concluded that soft robot undulation is ruled by dynamic modes induced by robot–environment interaction.

1. Introduction

Bionic slender swimming robots take inspiration from natural anguilliform locomotion, primarily observed in elongated animals such as eels, lampreys, and some amphibians. This type of locomotion falls under the category of undulatory swimming, which involves wave-like movements of the body to propel the swimmer forward [1,2]. The effectiveness of animal aquatic locomotion and an energetically efficient swimming mode is conceived as an ideal paragon in designing swimming robots that replicate animal gaits instead of using conventional marine machines with rotary propulsion [3,4,5,6]. On the other hand, robotic devices with mechanical properties optimized according to the principles observed in biological paragons are developed to study the undulatory locomotion of fish [7]. Additionally, from a biological point of view, robots that resemble biological systems in structure and motion can be studied to understand neuro-dynamics and the evolution of animals [8,9]. Hypotheses regarding the locomotion and functionality of animal nervous systems can be tested using bioinspired robots. The models derived and developed biomimetically form the basis for discovering technical design principles.

From an engineering point of view, the effectiveness of biological systems in aquatic locomotion is represented in different terms, such as the underlying control system producing the relevant locomotion patterns and the agility of the animals, which far surpasses that of concurrent robotic systems. Many fish, amphibians, and even reptiles propagate lateral undulations of their bodies and tails to generate propulsion in water [9]. This undulatory movement creates a wave-like motion that efficiently propels them through their aquatic environment. Salamanders, with their terrestrial locomotion pattern coupling and thus restricting their rhythmic horizontal bending, exemplify iconic models for “elementary” swift swimming dynamics control [10]. In its anguilliform swimming, the entire body of the salamander participates in creating the undulatory waves, which travel from head to tail. Researchers have developed sophisticated robots [11,12] replicating the anguilliform locomotion of salamanders. In constructing such robots, geared servomotors and rigid-body robotics play the main roles. However, a rigid-body system is fundamentally different from biological systems, leading to different physics. The difference is measured in terms such as actuation force (muscle force vs. motor torque) and material density or inertia. Soft robotics, on the other hand, benefits from similarity to biological systems in those terms.

Soft robots, generally continuum structures without the rotary joints of conventional systems, are concurrently developed to mimic the swimming gaits of aquatic animals [13,14,15,16,17]. Most soft robots are made from elastomeric materials actuated with different principles, mainly pneumatic pressurization [14]. Soft robots may consist of one single (often trunk-like [18,19]) soft body or may contain multiple trunks, like an octopus or starfish [13,20]. In particular, soft pneumatic actuators (SPAs) combine multiple resemblance factors, including flexibility, strength, and density, compared to natural organisms. Another similarity, referring to McKibben actuators, is that artificial muscles typically have their maximum force at the nominal length, which decays when the muscle contracts [21]. This aligns with the biological discovery that proprioceptive sensory receptors in lampreys reduce muscle activity on the concave side of the bend [22].

In eels and salamanders, the exhibited undulation relies on rhythmic contractions of the axial muscles along their highly flexible body [9]. The nervous system coordinates this periodic muscle activity through central pattern generators (CPGs) embedded in the spinal cord [23]. These CPGs generate rhythmic motor patterns that create the sinusoidal waves necessary for directed swimming. Electromyographic (EMG) recordings from swimming salamanders, such as the tiger salamander (Ambystoma tigrinum), show periodic bursts of activity in the axial muscles [24].

Rhythmic actuation implemented in rigid swimming robots has shown promising results. Rhythmic, in this context, refers to periodically repeating patterns applied as inputs or setpoints to the actuators of the robots. Such rhythms can be mathematically expressed, e.g., by sinusoidal or exponential spatiotemporal functions. A biology-inspired hypothesis is that rhythmic stimulation produces locomotion patterns, resulting in the swimmer’s forward motion. In rigid robots employing geared electromotors, the rhythmic actuation can be implemented as periodic desired angle inputs to the motors. In that case, the joint angles will follow the input instructions, with the negligible counter-influence of external loads on the motor shaft angle. In fact, as a general technological limitation of electric motors, the motors must be equipped with gearboxes to push water. However, in biological swimmers and the soft robots imitating them, the physics is governed by a two-way fluid-body interaction. Technically speaking, due to compliance, external forces applied by the fluid considerably influence the dynamics and deformation of the soft body. For this reason, a soft robot is a better reference for bioinspired locomotion. Additionally, metal components and motors make swimming robots heavy and bulky. Using silicon-based soft pneumatic actuators (SPAs), soft robots become low-weight systems that can be balanced in water with minimal additional weights or floats.

In recent years, many researchers have developed eel-like soft robots employing multiple (multi-segment, in other words) SPAs [25,26,27,28,29,30,31,32]. The eel robots belong to the most common types of soft robots, SPAs, but examples also include cable-driven soft robots, as seen in [33,34]. The soft eel robots are constructed of serially connected segments, which individually can bend laterally. The trunk-like body of soft eel robots allows for the implementation of pumps inside the robots, as seen in [25,26,29,32], though some of the robots, as seen in [27,29,30,31], depend on external pneumatic pressure and transmission piping. The studies generally concentrate on soft robotics and its challenges. The fabrication method and procedure may be the first, but also major, challenge and contribution within the soft robotics context. Being accented in aquatic applications, leakage is a general problem in SPAs. Research on the improvement of SPAs remains state-of-the-art [35]. In [31], a fabrication method for prototyping the soft eel robot is presented. One interesting aspect of the design proposed in [31] is that a smooth outside surface and fish-like anatomy are proposed, despite using chambers (teeth of pneunet structure, in contrast to other prototypes such as in [26,29]). The next soft robotics contribution in developing eel robots is related to control and robot locomotion studies. This part can be common with multi-body eel-like robots, e.g., using CPG-based controllers [27], deep reinforcement learning [33], etc. Due to the high density of water resulting in high inertial hydrodynamic forces besides viscous resistances, making soft robots swim is complex, as they are essentially for low interaction forces. The literature on soft eel robots is not abundant and deals with different aspects of robotics.

As a side effect of their infinite degrees of freedom (DoFs), soft robots are under-actuated. The technical simplification that the local bending of a soft robot results in a traveling wave remains an assumption that needs investigation. Furthermore, the existence of a smooth traveling wave as a testimonial of the anguilliform locomotion ability of the soft robot might not be visually apparent. This paper introduces a new eel robot that employs McKibben SPAs on a flexible beam to provide a two-dimensional (2D) system (from a hydrodynamics point of view). We increased the robot segments to six (compared to previous literature) to provide more flexibility. Experimentally, the actuators were excited rhythmically for steady swimming, and the robot showed successful swimming capability. Note that diving into acceleration (in contrast to steady swimming), obtained by arrhythmic and a different kinematic [36], is not within the scope and interest of this study.

This paper is organized as follows. In Section 2, first, the soft robot structure and its working principles are described in Section 2.1, then, a finite element (FE) model of the system developed using COMSOL Multiphysics^® 6.1 (COMSOL AB, Stockholm, Sweden) is introduced in Section 2.2. The physically realized system and the experimental setup are described in Section 2.3. Section 2.4 describes the underlying equations and relations governing the 2D laminar flow of incompressible fluid, moving mesh (fluid domain), and fluid–structure interaction (FSI). The FE simulation results are further analyzed using MATLAB^® R2022b (MathWorks^®, Natick, MA, USA) with an algorithm proposed to obtain the dynamic mode shapes, developed in Section 2.5. The experimental results are presented in Section 3.1. Section 3.2 is devoted to numerical simulations and extracting the mode shapes attributed to the traveling or stationary waves. Section 4 discusses how the proposed modal decomposition describes the backward traveling wave and the adversarial (seemingly noisy) behavior, based on the robot–environment modes.

2. Materials and Methods

2.1. Bio-Inspiration and Conceptual Design

In the animal locomotion context, four main modes of fish swimming are commonly reported: subcarangiform, carangiform, thunniform, and anguilliform. The anguilliform mode, named after the eel Anguilla anguilla, is more distinct from the others due to the undulation of the whole body in the form of backward-moving waves instead of a remarkable contribution of fins or limbs. Performers of anguilliform locomotion in nature encompass flexible-bodied and elongated organisms with body sizes ranging from microscopic nematodes to long snakes, moving in highly viscous environments such as open water and mud, or even terrestrial environments. Some fish larvae, before developing their mature swimming mode, perform anguilliform motion. Anguilliform swimming is known to be energy-efficient, i.e., low energy (or fat) consumption relative to the traveled distance.

Nevertheless, the reasons for the efficiency of this mode of motion and its optimality are not understood. Researchers develop laboratory settings and robotic models for empirical studies on natural swimmers and robots. Commercial software is widely used for numerical modeling and simulations of fluid dynamics. Nevertheless, intensified multidisciplinary studies are required to understand the underlying hydrodynamics and its interaction with swimmers’ bodies [37]. The study of swimming robots and fluid–robot interactions helps illuminate the underlying physics. Such studies on natural organisms are subject to behavioral aspects that can obscure underlying effects on swimming kinematics. The transfer of concepts from biology to engineering should begin early to minimize potential confounding factors present in natural models. This allows for a clearer focus on the fundamental physics that can serve as the foundation for engineering applications. Additionally, the duration of iterative interaction cycles between biology and engineering may be reduced.

The conceptual design of the robot can be described as follows. The soft robot consists of a highly flexible continuum beam, partitioned into six independently bending modules or segments. Each segment has lateral pneumatic actuators that enable planar bending. Schematically, in Figure 1, the blue curve represents the backbone beam, and the R_iR_j and L_iL_j (i and j = 1...7) lines represent the pneumatic actuators connected to the right and left sides of the backbone via rigid connectors A_iR_i and A_iL_i, which are also the floats for the buoyancy of the robot. On the neutral axis, A₁A₇, the A_i points are distributed equally. Due to segmental bending, the R_i and L_i points are dislocated to a new posture. To achieve this, the muscles contract when the segment is concaved towards the muscle, and antagonistic muscles are inactive. A waveform similar to the salamander swimming gait, which is a traveling sinusoidal wave with uniform amplitude and wavelength, is supposed along the swimming axis. As the exaggerated illustration in Figure 1b suggests, the muscles are thought of as linear actuators. Natural muscles generate force through neuromuscular activation, leading to the active contraction of muscle fibers. Additionally, when relaxed or subjected to external forces, the muscle can undergo passive elongation, allowing the muscle fibers to extend inactively. This dual functionality is essential for facilitating movement and maintaining structural integrity in biological systems. In contrast, artificial muscles, known as McKibben actuators, contract by pressurization but do not exhibit passive extension. To make the functionality more similar to the biological system, in this work, the actuator nominal length is taken as

$$l=(1+0.5{\lambda }_n)\left|\overline{A_i{A}_{i+1}}\right|,$$

(1)

where ${\lambda }_n=0.2$ is the actuator’s nominal contraction rate, i.e., the ratio of displacement to the initial actuator length. Therefore, muscle clearance allows for inactive extension. When pressure is applied to an actuator, the clearance is quickly resolved, and the muscle performs the contraction task similar to its biological counterpart. Note that, in this manner, the actual contraction rate is supposed to be $\lambda =0.1$.

2.2. Design and Model of the Robot

The proposed structure is given in Figure 2a. The overall dimensions of the 2D model (in mm) are shown in Figure 2b. The environment in this study is a shallow water tank with a width of 200 mm and a depth of 60 mm, on whose surface the robot swims. An FE discretization, as shown in Figure 2c, is used for simulation. Note that the SPAs are not modeled with FE elements but, as described in the next subsection, are modeled as external forces. The 2D model considering the shallow tank approximation was developed to simulate the robot in the environment using COMSOL Multiphysics^® 6.1. The fluid domain is modeled assuming laminar flow and incompressible fluid (water). The method of moving mesh (Yeoh model with stiffening factor C₂ = 100) was used to model large deformations of the deforming domain. No-slip boundary conditions (BCs) were adopted for the domain boundaries, except for the inlet and outlet. The robot backbone, floats, and tail are modeled with a linear elastic material. The forces of the muscles are modeled as external forces applied to the solid. A fully coupled FSI was considered on the solid BCs. Time-dependent simulations were performed with 0.001 s time steps. For the discretization, a user-controlled free mesh with triangular elements is used for the fluid and solid domains, along with two boundary layers and corner refinement (element size scaling factor 0.25) implemented on all boundaries except the inlet and outlet.

In two-dimensional models, triangular elements can easily conform to irregular shapes and boundaries, providing better accuracy in complex geometry and higher resolution in regions with high gradients. User-controlled meshing allows customization of the mesh density and distribution according to the specific needs of the simulation. Many physical phenomena, such as velocity gradients in fluid flow, occur near boundaries. Boundary layers in the mesh help capture these gradients accurately. Including boundary layers ensures that the mesh can resolve steep gradients near walls, leading to more precise simulation results. By refining the mesh near boundaries, where high gradients are expected, the mesh can be coarser elsewhere, optimizing computational resources. Corner refinement (element size scaling factor 0.25) means the size of the elements in these refined regions is 25% of the size of elements in the remaining domain. Corners often experience highly sharp changes in flow characteristics, and increasing mesh density near corners ensures these effects are accurately captured. The refinement is controlled and localized by scaling the element size, preventing an unnecessary increase in the overall number of elements. At the inlet and outlet, the flow is usually assumed to be more uniform and less prone to sharp gradients or complex interactions that require dense meshing. These techniques allow for using a coarse mesh, preventing overburdening computational resources.

Note that fixed meshes lead to inaccuracies when large deformations occur because they do not adapt to the changing shape. The moving mesh method ensures that the mesh conforms to the deformed shape, maintaining the accuracy of the simulation. The moving mesh technique maintains the continuity of the mesh, avoiding issues such as element distortion or mesh tangling. Using a material model like Yeoh ensures that the material behavior under large deformations is realistically captured, reflecting how the beam is influenced by the flow-induced forces. In the experiment, the tank is 2 m long, and the robot swims in stationary water. The FE model contains 1 m of the water tank, while the robot’s heading side (the right side, referring to Figure 2) is the inlet of the incompressible flow, and the tail side is the outlet. The BC at the inlet is a normal inflow velocity obtained from the experiment, as will be described in subsequent sections, and at the outlet is a fully developed flow with the given average velocity.

2.3. Hardware and Its Modeling Assumptions

The robot prototype used in the experiments is shown in Figure 3a. To achieve a 2D in-plane bending, at each side of each segment, typically shown in Figure 3b, a pair of connected thin McKibben actuators are implemented (the paired actuators are pressurized together with a single pressure pipe). The robot in the test environment is shown in Figure 3c. The filled water depth is approximately 20 mm more than the backbone width to make up for the corner chamfers. A QualiSys^® Tracking Manager (QTM) version 1.6.0.1 (QualiSys AB, Göteborg, Sweden) motion capture system is used to track the robot’s motion using five cameras. Seven markers were located on the robot’s floats. The floats are at the equally distanced points denoted as A₁ to A₇.

After system calibration, the robot was placed in the water tank and the electro-pneumatic controller was activated. The robot swam from a stationary state to the end of the tank, and the position of the markers was logged. The coordinate system is chosen so that X is along the motion and XY is the horizontal surface. The (x, y, z) values of markers were saved with 200 Hz incidences. The data, saved as ‘.tsv’ files, were processed offline using Python 3.8.19 (Utilizing Spyder 5.5.1 as the integrated development environment). The electro-pneumatic actuation system (depicted for one side of a typical segment for simplicity) is shown in Figure 4. In pneumatics, the actuator can be expressed as a single-acting pull-type pneumatic cylinder with a spring return, as in Figure 4a. The force can be modeled as a constant pneumatic force f_p, which is subject to a linear spring (elasticity), as shown in Figure 4b. In other words, in the FSI model, the McKibben actuators are modeled as ‘constant forces’ applied to the structure, while the ‘elasticity’ of the structure in the FSI model provides the spring effect, as shown in Figure 4c.

Note that soft robots generally have small elasticity moduli (if modeled as elastic materials). However, with some static solid mechanics FE simulation trials, the elasticity modulus required to stop the applied force (40 N from the actuators) at a point where the maximum contraction ($\lambda =0.1$) is reached was obtained by selecting E = 10 kPa. In the actual system, the dominant forces are hydrodynamic forces and those from McKibben actuators (the actual springiness of the robot backbone is negligible). However, in the FSI model, an elasticity modulus (E = 10 kPa) is considered to model the linear decay of the actuator forces. The increased elasticity allowed for modeling the soft robot with COMSOL Multiphysics^® 6.1, and helped achieve converging solutions, according to our experience. FSI problems are generally complex, and achieving converging solutions demands several considerations, including FE modeling tips and engineering intuition.

Another modeling approximation in the FSI simulation study is related to the mass effects. The robot’s buoyancy was achieved experimentally by mounting small metal weights and 3D-printed floats. This allows us to model the robot’s structure with a homogeneous material with water density. In computational methods such as the FE FSI, engineering intuition and approximations are essential for creating simplified models of a mechanical system. As in the case of this study, a 2D model can capture the dominant dynamics of an actual system, while the computation of a 3D model is considerably more expensive. Therefore, the 2D model (Figure 2) was designed with the overall sizes of the actual robot (Figure 3). The mass and volume of the small-sized actuators (1.9 mm un-actuated diameter, 4.1 mm actuated) are neglected in the 2D model. Another approximation in the FSI modeling is related to the body thickness.

The robot backbone thickness (approximately 1 mm in the actual robot) is scaled (to 10 mm) in the FSI model. This approximation is based on the fluid dynamic fact that, at the BCs, the fluid moves with the solid (in FSI, an approximation method known as the added mass is based on this phenomenon). As the robot’s body density is equal to that of the fluid, as discussed earlier, it is reasonable to implement such an assumption. This is, in practice, a regularization of the FE mesh, which makes the underlying matrix equations simpler to solve. The thinner the beam is modeled, the finer the mesh (triangular solid elements) is required, making in turn the matrix equations larger and causing convergence problems for the FSI solver. With the FE model dimensions and the SPA characteristics, a simple static problem (plus a linear regression) can be solved for the elasticity modulus that is compatible with the maximum contraction of the SPAs.

We leverage a fact in the theory of elasticity that the elasticity module and external forces are linearly related in small deflections, thus, can be scaled by the same factor. For example, a steel I-beam subjected to a load can show the same curvature as an aluminum beam with the same size and geometry but with a properly scaled load. To focus on steady swimming, a simple open-loop control was implemented to avoid acceleration and the effects of feedback controller dynamics on the robot.

The microcontroller simply provides a monotonous activation for the actuation of each segment, as shown in Figure 4a. The activation signal is simply a square wave that turns the solenoid pneumatic valves on and off. The segmental actuation pattern is shown in Figure 5. A microcontroller system provides the activation logic given the actuation period, T. Note that optimization of the swimming performance is not within the scope of this study. The methodology followed in this research was to achieve a working setup empirically, and then investigate and analyze the system using computational and modeling methods. On the other hand, the scope is not exclusively focused on soft robotics or the development of a soft robot, though a novel robot is introduced.

2.4. Robot–Environment Interaction

To study a system in which a fluid and a deformable body affect each other, various FSI methods have been developed to model the fluid and the solid domains (structures) with predefined conditions at the fluid–solid boundaries. The Arbitrary Lagrangian–Eulerian (ALE) formulation can effectively handle the dynamic interactions between fluid and solid domains. An ALE method incorporates the geometrical changes of the fluid domain, allowing the mesh to move and deform to accommodate the motion of the solid structure without the mesh being tied to the material points (as in a fully Lagrangian method). The solid domain is handled using a Lagrangian approach (the mesh moves with the material), while the fluid domain is treated using a Eulerian approach (the mesh deforms but does not move with the material).

In the swimming robot, the interaction between the fluid and the robot results in significant deformation of the fluid domain, demanding a mesh that can adapt to these changes. Using the moving mesh method, the fluid mesh deforms to match the movement of the robot boundary but remains independent of the fluid particles. The mesh follows the robot’s structure, capturing its motion accurately. The ALE method’s ability to manage large deformations makes it ideal for FSI problems involving moving boundaries. In this method, structural mechanics is used to compute the deformation of the solid domain based on the applied forces and the boundary conditions. The moving mesh model uses the solid’s displacement to update the mesh in the fluid domain, ensuring that the mesh conforms to the deformed solid boundaries. Note that the moving mesh, as shown in Figure 6, is the fluid domain in this study. The fluid flow physics solves the modified Navier–Stokes equations in the updated mesh to obtain the fluid velocity and pressure fields. Satisfying the boundary conditions at the fluid–solid interface ensures the continuity of velocity and stress. This process iterates over each time step, with the fluid and solid solvers exchanging information at the interface to maintain the coupled interaction, and the moving mesh module continuously adjusts the fluid mesh to follow the evolving geometry. This iterative coupling continues until convergence is achieved for each time step before proceeding to the next.

In the context of COMSOL Multiphysics^® 6.1, the modified Navier–Stokes equation is given by

$$\underset{Inertia\hspace{0.17em}term}{\underbrace{\rho {\displaystyle \frac{\partial \mathrm{u}}{\partial t}}+\rho (\mathrm{u}.\nabla )\mathrm{u}}}=\underset{\begin{array}{l}Stress\hspace{0.17em}divergence\\ term\end{array}}{\underbrace{\nabla \cdot \left[-p\mathrm{I}+\mathrm{K}\right]}}\underset{Shallow\hspace{0.17em}channel\hspace{0.17em}approximation}{\underbrace{-12\mu {\displaystyle \frac{\mathrm{u}}{d_z^2}}}}$$

(2)

where the parameter $\rho$ is the fluid density, $\mathrm{u}$ is the fluid velocity vector, and $\mu$ is the dynamic viscosity. The inertia term accounts for the changes in fluid momentum due to both local acceleration and the convective effects of fluid motion. The stress divergence term includes contributions from both the fluid pressure p (where I is identity matrix), and additional stresses represented by the stiffness matrix $\mathrm{K}=\mu \nabla \mathrm{u}+\mu {(\nabla \mathrm{u})}^T$. In COMSOL Multiphysics^® 6.1, the stiffness matrix K is computed from the coupling between the fluid’s stress tensor, the pressure field, and the structural response using FE and the integration of the effect on the interface. The viscous dissipation term is based on the shallow channel approximation, which is adopted when the depth is much smaller than the length of the channel to simplify the equation. The continuity equation for incompressible flow is written as follows:

$$\nabla \cdot \mathrm{u}=0 \mathrm{or} {\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial v}{\partial y}}=0$$

(3)

Knowing the lateral velocity v (conventionally $\mathrm{u}=u\widehat{i}+v\widehat{j}$), Equation (3) can be solved for u. Therefore, the component will be used to represent the state vector. The continuity equation and the modified Navier–Stokes equations are solved in the updated mesh for the fluid velocity and pressure fields. For the solid domain, the momentum equation implies

$$\underset{Inertia\hspace{0.17em}term}{\underbrace{\rho \left({\displaystyle \frac{{\partial }^2{\mathrm{u}}_{solid}}{\partial t^2}}\right)}}=\underset{\begin{array}{l}Internal\hspace{0.17em}\\ forces\end{array}}{\underbrace{\nabla \cdot {\left[\mathrm{FS}\right]}^T}}\underset{\begin{array}{l}External\\ forces\end{array}}{\underbrace{+f_a}}$$

(4)

where ${\mathrm{u}}_{solid}$ represents the displacement, F is the deformation gradient tensor given as $\mathrm{F}=\mathrm{I}+\nabla {\mathrm{u}}_{solid}$, and ${\left[\mathrm{FS}\right]}^T$ is the first Piola–Kirchhoff stress tensor transformed to the current configuration. As the robot is balanced using floats and added masses to buoyance on the water surface, the solid density is assumed to equal that of the fluid. For the specific problem, a rigid motion suppression condition is read ${\displaystyle \sum {\mathrm{u}}_{solid}=}{\displaystyle \sum \mathrm{X}\times {\mathrm{u}}_{solid}=}0$ where X represents the material coordinate. Additionally, the plane–strain assumption was adopted for the solid domain. The interface equations that make the FSI coupling between the flow equations and the solid mechanics are

$$f_a=[-p\mathrm{I}+\mathrm{K}]\cdot n \mathrm{and} {\mathrm{u}}_{\mathrm{tr}}={\displaystyle \frac{\partial {\mathrm{u}}_{\mathrm{solid}}}{\partial t}}$$

(5)

where f_a is the force vector applied to the structure from the fluid, and n is the normal vector to the robot surface. The first equation projects the force per unit area exerted by the fluid on the robot. The second equation ensures that the ALE mesh velocity, ${\mathrm{u}}_{\mathrm{tr}}$, is equal to the solid velocity. Note that the wall BCs imply that $\mathrm{u}={\mathrm{u}}_{\mathrm{tr}}$. In the case that fluid does not influence solid deformation, the first term of (5) is supposed to be zero. Considering both terms means solving a fully coupled FSI problem.

For each time step, solving the equations is followed by a mech update, which involves solving an equilibrium equation using the Yeoh model written as follows:

$$\nabla \cdot {\sigma }_{mesh}\hspace{0.17em}=0 , {\sigma }_{mesh}\hspace{0.17em}={\displaystyle \frac{\partial W_{mesh}}{\partial {\mathrm{F}}_{mesh}}}$$

(6)

where $W_{mesh}=C_2{(I_1-3)}^2$ is the strain energy, with $I_1=trace({{\mathrm{F}}_{mesh}}^T{\mathrm{F}}_{mesh})$ being the first invariant of the mesh deformation tensor, and ${\mathrm{F}}_{mesh}=\mathrm{I}+\nabla {\mathrm{u}}_{mesh}$. By specifying C₂ = 100 in the Yeoh model within COMSOL Multiphysics^® 6.1, we focus on the quadratic term in the strain energy density function, employing the ‘nonlinear’ stiffening effect as the mesh deforms. The external boundaries of the fluid domain are fixed (${\mathrm{u}}_{mesh}=0$). This setup is handy for simulating large deformations.

2.5. Proposed Modal Decomposition

To comprehend the robot–environment interaction and reveal the underlying patterns or modes of the interaction, the FE model is used as the digital twin, producing synthetic data for the mode extraction method described as follows. The FE simulation results, i.e., the time history of the system state vectors $\overrightarrow{v}\in {ℝ}^{m\times 1}$, are stored as a data matrix for postprocessing using MATLAB^® R2022b. The state vector for the fluid corresponds to the vertical velocity field written as a vector at a given time. For points of the solid domain, i.e., the spatial points where the solid is present, the state is filled arbitrarily with a value smaller than the minimum velocity (because COMSOL Multiphysics^® 6.1 returns NAN values for v in points where the solid exists). To eliminate the transient part of the FE solutions, the initial periods’ simulation results were discarded (the simulation was performed for 11 cycles, but the results of the first 5 cycles were discarded). The data matrix, $D\in {ℝ}^{m\times n}$, contains $\overrightarrow{v}$ as its columns.

$$D=\left[\begin{array}{ccc}& \left.\right|\hspace{0.33em}& \\ \cdots & \overrightarrow{v}& \cdots \\ & \left.\right|\hspace{0.33em}& \end{array}\right]$$

(7)

The next step of the proposed mode extraction involves the decomposition of the data matrix based on singular value decomposition (SVD) and order reduction, expressed as follows:

$$D=\Psi W\hspace{1em},\hspace{1em}W=\sum V^T$$

(8)

where $\Psi \in {ℝ}^{m\times r}$, $\sum \in {ℝ}^{r\times r}$, and $V\in {ℝ}^{n\times r}$ are obtained from the reduced-order SVD of order r. Note that the superscript, ^T, represents the transpose of the matrix. The matrix, $\Psi$, contains the spatial modes as column vectors, written as

$$\Psi =\left[\begin{array}{ccc}\left.\right|\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}& & \left.\right|\hspace{0.33em}\hspace{0.33em}\\ {\psi }_1& \cdots & {\psi }_i\\ \left.\right|\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}& & \left.\right|\hspace{0.33em}\end{array}\begin{array}{ccc}& \hspace{0.33em}& \left.\right|\hspace{0.33em}\hspace{0.33em}\\ \cdots & & {\psi }_r\\ & \hspace{0.33em}& \left.\right|\hspace{0.33em}\end{array}\right]$$

(9)

Reshaping the ${\psi }_i$ vectors in the reverse manner of the vectorization in the construction of the D matrix exhibits the FSI mode shapes of the system.

With s number of sensors (in this work, tracking markers) and n time snapshots, the experimental data is presented as $Y\in {ℝ}^{s\times n}$. Next, we seek for $\varphi$ so that $\varphi W=Y$. Here, $\varphi$ can be thought of as a new coordinate system that, due to the invertibility of $W$, can be obtained as follows:

$$\varphi =Y{W}^{-1}\hspace{1em},\hspace{1em}{W}^{-1}=V{\sum }^{-1}$$

(10)

The matrices $\varphi$ and $W^T$ are represented with their column vectors as:

$$\varphi =\left[\left.{\phi }_1\right|\left.\cdots \right|{\phi }_r\right], W^T=\left[\left.w_1^T\right|\left.\cdots \right|w_r^T\right]$$

(11)

Then, the reconstructed $Y$ with the superposition of ith to jth terms is defined as follows:

$${}_i{}^j\tilde{Y}={\displaystyle {\sum }_{k=i}^j{\phi }_k{w}_k}$$

(12)

In this equation, ${\phi }_k$ can be interpreted as mode shape of the kth mode, with temporal weights of $w_k$. In this context, instead of a reconstruction like $Y\approx {}_1{}^r\tilde{Y}$, we are interested in investigating the effects of specified or individual terms. By specifying i and j in (12), the contribution of modes i to j on the measured signal is calculated.

3. Results

3.1. Experimental Results

The experimental setup, explained in Section 2.3, was employed for the actual swimming tests and for collecting the experimental data. A snapshot of the robot during swimming along the tank is shown in Figure 7. Seven gray spherical tracking markers are located on the green floats. Four static markers are used to indicate the water tank borders. Figure 8 shows an observation snapshot taken by the QTM 1.6.0.1 system, highlighting the coordinate system.

The recorded data were further preprocessed and visualized using Python 3.8.19, as shown in Figure 9. The X component of the recorded positions, shown in Figure 9a, suggests a uniform average velocity of the robot’s motion. The Y components, shown in Figure 9b, correspond to the lateral displacement of the robot backbone at points A₁ to A₇. At the starting phase of the motion, in approximately 10 s, the graphs show transients. In biological swimmers, this phase is known as acceleration and is associated with different kinematics. This study focuses on steady swimming, so the data regarding the transient response are discarded. The steady swimming phase observed between t = 15 s and t = 25 s shows an average forward velocity of 28 mm/s and a frequency of 0.9 Hz (the period is approximately 1.1 s). Further analysis and explanation of the experimental results require more knowledge about the environment and the interaction phenomena. For this reason, simulation studies are proposed in the next section.

3.2. FSI and Mode Analysis

The system was modeled according to the details provided in the Methods section. The flow inlet velocity was equal to the robot’s velocity obtained experimentally. After the simulation, COMSOL Multiphysics^® 6.1 provides extensive postprocessing tools for visualizing the results. Streamlines help visualize the flow paths, while arrow plots show the direction and magnitude of the velocity field.

Volume and surface plots display scalar fields over volumes or surfaces within the domain, providing a comprehensive view of the distribution of variables. Simulation data can be exported in various formats for further analysis or presentation, including text files, Microsoft^® Excel^® spreadsheets, images, and videos. Figure 10 shows some results in ten snapshots. The simulation plots present important information about the system and fluid dynamics. The plot shows the evolution of the traveling waves along the robot’s backbone. Similar to the experimental observations, some imperfectness (in contrast to the perfectness of a sinusoidal wave) is observed in the shape of the backbone. The exotic fluctuation is represented in Figure 11, which shows the lateral displacement for the midline points (A₁ to A₇, plus the distal end at the tail). The simulation data suggest that the behavior is due to the physical interactions (in contrast to noise or disturbance).

These visualizations can provide valuable insights for researchers studying anguilliform locomotion. Nevertheless, raw data from a numerical method can be used to discover more interpretable phenomena through postprocessing based on machine learning techniques. Mode decomposition techniques allow us to interpret the dataset and reveal patterns that reflect the underlying system dynamics. The system mode shapes and their time evolution were measured using the simulation result data and the method mentioned in Section 2.5. Figure 12 represents the schematic representation of the results for the data matrix, D, and its relation to the mode shapes. Note that in the data matrix, the time snapshots are reshaped into vertical columns, and for the vectorization, the horizontal rows are read left to right (transposed). The calculated FSI mode shapes, ${\psi }_i$, for the first modes are shown in Figure 13a.

Using the experimental data, $Y$, the corresponding spatial modes, namely ${\phi }_k$, were measured and illustrated in Figure 13b. The robot’s head is shown as A₁, and the last backbone marker (not the tail) is A₇. The corresponding temporal modes are shown in Figure 13c. The first mode can be interpreted as the overall lateral motion of the robot’s backbone. The plots imply that mode 2 and mode 3 contribute to the traveling wave related to the anguilliform locomotion of the robot. Other higher modes are interpreted as fluctuations in the robot–environment system. For better presentation and comparison, the backbone mode shapes are shown in Figure 14, for modes 2 and 3, and in Figure 15, for the higher modes. Both the mode shapes shown in Figure 14, and the relevant temporal evolution, Figure 13c, represent out-of-phase waves. The remaining modes, shown in Figure 15, appear as fluctuations that cannot produce a wave traveling along the robot’s backbone. The results show that the fluctuations have smaller amplitudes and higher frequencies (than the main modes) but are not negligible.

Thus, it will be interesting to reconstruct the robot’s backbone using the modes responsible for the traveling wave. In other words, the decoded data are encoded, keeping modes 2 and 3 and eliminating other modes. For this purpose, ${}_2{}^3\tilde{Y}$ was calculated to represent the pattern of anguilliform motion of the robot. Figure 16 shows the result for one period. The backward traveling wave is visually trackable, confirming the existence of an anguilliform regime in the locomotion. The wavelength-to-length ratio is calculated from the data to be $0.63$ (without considering the length of the inactive tail).

The values, including the length of the tail, would be 0.49. Anguilliform swimmers undulate their bodies at wavelengths between 0.5 and 1L, and in contrast, for subcarangiform and carangiform swimmers, the value is equal to or greater than L [37]. In [38], values between 0.65 and 0.80 for optimal hydrodynamic performance are estimated based on numerical simulations. Superposition of the higher modes, ${}_4{}^7\tilde{Y}$, can exhibit the pure effect of the higher modes. The calculated result is presented in Figure 17, which shows that the higher modes appear as stationary (concerning the robot) bending, causing additional fluctuations along the midline. The modes with the belly and nodes do not contribute to a traveling wave but, in reality, are transferred to the robot and create distortions. They are also evidence of the softness of the robot in the environment. In other words, the robot is compliant, as it is influenced by external forces.

4. Discussion

The efficient locomotion of aquatic and amphibious animals has attracted robotics research, leading to the development of bioinspired robots that mimic the topology and kinematics of these natural creatures. Rigid robots with servomotors as their actuators cannot replicate the compliance of such animals. In this study, a soft eel-like robot was designed and tested empirically for studying anguilliform locomotion. The robot contains six segments, each bending laterally using pneumatic actuation. The hypothesis that a rhythmic pattern of actuation and local bending causes anguilliform locomotion in the robot was implemented and investigated experimentally and numerically. Experiments showed the robot’s capability of moving forward in a still-water tank along the axial direction. The robot’s midline was digitized using seven tracing marks set equidistantly along the robot’s midline. The experimental data presented exotic fluctuations that could not be explained by a single sinusoidal mechanical wave moving along the robot’s midline.

A simulation study was performed employing COMSOL Multiphysics^® 6.1 to simulate the robot in the fluid as the environment and, in particular, to confirm the responsibility of robot–environment interaction effects on the seemingly noisy and non-sinusoidal backward-moving wave. The results convinced us that, due to compliance, the soft robot’s behavior is highly influenced by the robot–environment interaction. Further postprocessing of the data based on the proposed method revealed the anguilliform mechanical wave, consisting of two dominant modes. Furthermore, higher modes were detected as stationary modes responsible for lateral fluctuations along the robot’s backbone.

This conclusion supports the opinion in biology, which states that body deformation is the result of both internal musculoskeletal forces and external fluid dynamic pressures. Additionally, similar to this robot, natural anguilliform swimmers like eels are known to present both traveling and stationary body waves, which, based on the results of this paper, can be attributed to the influence of additional modes (not the muscle activations).

A controversial topic in the context of anguilliform swimmers is how the CPG and local mechanosensory systems interact. Employing the McKibben actuator’s force, which decreases with contraction, the robot is actuated similarly without needing additional feedback control. With the six segments and actuator characteristics, the robot replicates a swimmer controlled by a CPG and local (segmental) mechanosensory systems. This study suggests that global (i.e., over the whole body) deformations can be due to fluid–body interaction, without interconnection of the local system. The contradicting hypothesis about the lampreys is that global body regulation might be based on the communication of the local curvature regulation systems (using neural interconnections, as reported in [39] for the zebrafish). The results of this research can serve as a foundation for further studies in robotics and biology. In robotics, topics such as sensing (or filtering), motion optimization, and robot control, and in biology, the study of the fluid–body interaction are of interest. Our future work will focus on modeling or identifying the temporal behavior based on the modal decomposition proposed in this study.

5. Conclusions

In this study, a soft eel-like swimming robot was designed and tested empirically. The robot contains six segments, which are actuated by McKibben pneumatic actuators. Each segment is bent laterally by two individual pneumatic control valves. The robot was tested in a still-water tank, and seven markers on the midline were recorded using a tracking system. The system was modeled employing the fluid–structure interaction (FSI) module of COMSOL Multiphysics^® 6.1, and the system modes were extracted and analyzed. The results revealed rearward traveling mechanical waves attributed to the anguilliform undulation, along with low-amplitude stationary modes, in the robot’s movement. It is concluded that the robot–environment interaction resulted in the anguilliform waves similar to natural swimmers but not following the excitation rhythm in terms of wavelength and amplitude.