Development of a Numerical Model of a Bio-Inspired Sea Lion Robot

1. Introduction

There has been a significant rise in the development of underwater robots, with the goal of improving the ability of these robots to swim and maneuver effectively in dynamic, high-energy flow environments [1,2,3]. This growth is driven by the increasing demand for robotic systems to assist with offshore tasks such as scientific research, critical underwater missions, and ocean exploration, needs recognized by organizations like the Office of Naval Research (ONR) and the National Oceanic and Atmospheric Administration (NOAA) [1,2,4,5]. These applications often require underwater robots to navigate tight spaces, avoid obstacles, operate at high speeds, and execute complex maneuvers with agility in energetic flows and occluded environments. To operate effectively in these conditions, underwater robots need to swim with great agility and maneuverability to take advantage of and adapt to multi-directional flows.

One promising approach to improve the performance of underwater robots is to develop robots that are able to understand and adapt the locomotion strategies of biological swimmers [6,7,8]. Many aquatic animals utilize their bodies and propulsors to swim and maneuver effectively in complex, unsteady flow environments. As a result, engineers want to adapt swimming techniques employed by biological systems onto underwater robots to improve swimming performance [9,10]. Current biologically inspired underwater robots often feature multiple body segments, fins, or flippers of various shapes, size, and degrees of freedom (DoF). For example, fish robots utilize multiple fins to produce and shape propulsive forces [11,12,13,14], sea turtle robots use soft morphing flippers to swim and traverse on land [15,16], and sea snakes and eel robots have multiple articulated bodies [17,18,19] to navigate confined places. As these robotic platforms become more advanced, it is increasingly important to understand how the actuation of various body segments influences overall swimming and maneuvering performance.

Numerical modeling is a powerful tool for studying the biomechanics of multi-body underwater robots; however existing approaches often rely on simplifying the model that limit their applicability for complex bio-inspired systems. The limited multi-body models that do exist primarily focus on snake-like robots composed of uniform, symmetric links [20,21,22]. Although deriving the equations of motion in closed form solutions provides significant advantages for stability analysis [23] and control design [20,22,24], such formulations have only been demonstrated for planar or similar linked systems rather than for freely swimming three-dimensional platforms with non-uniform segments. Unlike traditional underwater robots that rely on propellers for propulsion [25,26], bio-inspired unmanned underwater vehicles (UUV) utilize propulsors such as fins and flippers that generate unsteady, three-dimensional, time-dependent forces that vary in direction throughout each stroke cycle [13,27]. Moreover, each articulated body segment is not necessarily uniform in shape, size and mass. These differences affect distribution of inertia, mass and hydrodynamic forces for the entire system. This makes modeling particularly difficult, as each segment’s motion not only influences the forces and torques acting on itself but also affects adjacent segments, leading to strong interdependencies between the degrees of freedom. Another major challenge lies in the analytical computation of the hydrodynamic forces, which depends on accurate approximations of hydrodynamic coefficients [28]. Experimental methods, such as tow-tank experiments, have proven to be useful for computing these coefficients [29,30,31]. However, these methods require large, dedicated facilities, making them costly and often impractical. To address these challenges, a numerical modeling framework that integrates a closed-form dynamical model, accurately estimates hydrodynamic forces produced by multiple body segments, and provides effective kinematic visualization is essential.

To that end, the objective of this work was to develop and validate a numerical model for a biologically inspired sea lion robot to simulate, analyze and visualize the body motions of the robot in water (Figure 1). A bio-robotic sea lion was chosen for this study because its complex, multi-segmented body and three-dimensional flipper movements provided an ideal test case for validating the proposed methodology for the development of the numerical model. Additionally, as the robot was developed in-house, the authors had full access to its geometric and inertial properties, which was crucial for accurate model development and validation. While previous numerical models have addressed simplified or planar systems, this study, to the best of authors knowledge, is the first to derive the equations of motion in closed form solution for a freely swimming, multi-body unmanned underwater vehicle (UUV) with non-uniform body segments. The critical difference between our approach and existing models lies in our framework which overcomes the inherent challenges of modeling such systems. The equations of motion for the entire model were solved symbolically using the Euler-Poincare formulation [24], which take the generalized form of the Lagrange’s equations. Unlike, the Lagrange’s equations of motion that rely on the time derivative of the generalized coordinates, the Euler-Poincare model uses quasi-velocities which are velocity components defined in the local body-fixed frame [24]. This formulation made it easier to estimate the hydrodynamic forces produced by the articulation of different body segments. To accurately compute these forces, we developed a hybrid approach combining computational fluid dynamics (CFD) with analytical strip theory. Another significant contribution to this work is the application of a genetic algorithm to optimize hydrodynamic coefficients and other key parameters, which ensured a close alignment between our numerical model and the physical robot's performance. This optimization step is vital; it addresses a major challenge in aquatic robotics by effectively reducing the sim-to-real gap and providing an accurate, predictive tool for motion analysis and control. To visualize the kinematics of the vehicle, the Euler-Poincare model was integrated with MATLAB Simscape (The MathWorks, Inc., Portola Valley, CA, USA) to show how the bio-robotic platform moves over time. The joints and their respective kinematics in the numerical model were controlled using trajectory tracking controllers. These controllers were developed using extensions of feedback linearization methods to accurately follow the desired motion path.

Numerical modeling has proven effective to analyze the performance of multi-body robotic systems in terrestrial and aerial domains; however, its application in the underwater domain remains limited due to the challenges of capturing the dynamics of multi-body articulation of non-uniform body segments in a fluid environment. In the terrestrial domain, dynamic models of legged robots have been developed with accurate representations of ground contact forces, which enabled precise control and analysis of multi-limbed locomotion [32,33,34]. In the aerial domain, numerical models have been developed to capture the dynamics of both biologically inspired and engineered systems. For example, a numerical model of the common guillemot was developed to capture the kinematic and inertial effects of its flapping wings during aerial locomotion [35], and demonstrated how articulated control surfaces could be effectively modeled. Similarly, [36] developed a numerical model of a bird scale aerial vehicle that incorporated the effects of inertia change in wings and its passive deformation effects to analyze and demonstrate flight stability. In the underwater domain, early research primarily focused on torpedo and box shaped vehicles, with studies such as [25,26] using classical numerical modeling techniques based on [28] approach. Building on these foundations, only a few numerical models have been developed for bio-inspired underwater robots. One such study [37] , modeled an octopus-inspired robot capable of underwater legged locomotion (ULL). This work uniquely combined hydrodynamic sculling forces with ground contact forces to simulate benthic propulsion, demonstrating how body morphology, buoyancy distribution, and multi-modal actuation can influence underwater locomotion in shallow or near bottom environments. Most existing numerical models for multi-body underwater systems have been developed for snake-like robots which possess uniform body segments using both closed-form and non-closed form approaches [20,38,39,40]. Another important aspect of underwater robot modeling lies in accurately estimating the hydrodynamic forces. While many have computed coefficients using tow-tank experiments, more recently, computational fluid dynamic (CFD) methods [26,41,42] have been increasingly used demonstrating effectiveness and offering results comparable to experimental approaches. Thus, there remains a clear need for a physics-based, closed-form numerical framework capable of accurately capturing the hydrodynamics and kinematics of a freely swimming, multi-segmented, non-uniform underwater vehicle.

The remainder of this paper will discuss: (1) derivation of the closed form dynamical equations of motion for a bio-robotic sea lion platform using the Euler-Poincare formulation; (2) determination of hydrodynamic coefficients using CFD and strip theory; (3) experimental results of the bio-robotic sea lion executing flipper and body motions with numerical model validation.

2. Materials and Methods

2.1. Overview

A numerical model was developed to simulate, visualize and analyze the body kinematics for a biologically inspired sea lion robot during different swimming and maneuvering trials. The equations of motion for the bio-robotic swimmer were derived in closed form using the Euler-Poincare approach [24] (Section 2.2). To simulate the articulation of individual body segments in the UUV, independent trajectory tracking controllers were developed that followed the prescribed kinematics using extensions of feedback linearization methods [24] (Section 2.3). Hydrodynamic forces that were produced by the articulation of various body segments of the robotic swimmer were analytically computed and applied to the center of mass of each body segment in the numerical model (Section 2.4). The hydrodynamic coefficients, which were essential for hydrodynamic forces, were computed using computational fluid dynamic simulations and strip theory. A genetic algorithm was then used to refine these coefficients, that reduced the sim-to-real gap in terms of performance between the robotic swimmer and its numerical model. For visualization of the kinematics of the numerical model, the Euler-Poincare model was integrated into the MATLAB Simscape environment (Section 2.5). Performance metrics were then defined that allowed comparison and validation of the numerical model against the experimental platform (Section 2.6).

The bio robotic sea lion, referred to as the Stroke Experimentation and Maneuver Optimizing Underwater Robot (SEAMOUR) (Figure 2a), was developed to investigate how sea lions utilize flippers and body to swim and maneuver. The robotic system comprised of seven major body segments: a main body, head (Figure 2b), a pair of foreflippers (Figure 2c), pelvis (Figure 2d), and a pair of hind flippers. The head and a pelvis section each with two degrees of freedom in pitch and yaw relative to its body fixed frame (Figure 3a), were both connected to the central main body. A pair of rigid foreflippers were attached to the main body with three degrees of freedom (roll, pitch, and yaw) each, while a pair of rigid hind flippers (Figure 2d) was connected to the pelvis with two degrees of freedom (roll and yaw). SEAMOUR was an open system that allowed for water ingress and egress during operation. An external antenna was mounted on a small surface float during experiments. SEAMOUR was set to be neutrally buoyant with horizontally level orientation. A comprehensive description of SEAMOUR’s design, fabrication, and component specifications has previously been described in [43,44].

2.2. Derivation of Closed-Form Equations of Motion (EoM)

Numerical model of the bio-robotic sea lion was developed using the Euler-Poincare formulation, which takes the generalized form of the Euler-Lagrange equations [24]. In this approach, the equations of motion are expressed using the quasi-velocity variables, which define the linear and angular velocities in the local body-fixed frame (Figure 3a). This differs from the traditional lagrangian method, where velocities are defined as the time derivatives of the generalized coordinates [24]. The equations of motion were symbolically derived in Mathematica (Wolfram Research, Inc., Mathematica, Champaign, IL (2024)) using ProPac [45], a toolbox that integrates into Mathematica as an add on package. ProPac provides the necessary kinematic and dynamic functions for modeling multi-body vehicles. The entire kinematic and dynamic formulation of this sea lion robot was carried out with Mathematica using these tools. ProPac also includes functions to automatically generate C code, which can be compiled into Simulink compatible S-functions for simulation purposes [24]. A schematic of the overall numerical model structure is shown in Figure 4.

Detailed instructions on downloading and installing ProPac, setting it up with Mathematica, and running simulations in the Simulink environment can be found in [24].

SEAMOUR was modeled as a mechanical tree structure, where the system is represented as a chain of interconnected bodies [24,46]. Each body contained a distinct node, which served either as a connection point to another body segment or as the location where generalized forces were applied. The node placements for the various body segments of SEAMOUR are illustrated in Figure 3b. The main rigid body had 6 degrees of freedom (3 translation and 3 rotation). The head and pelvis were both connected to the main body with 2 degrees of freedom. A pair of hind flippers were connected to pelvis with 2 degrees of freedom each. In total, the model encompassed 6 degrees of freedom (DoF) in cartesian space and 14 degrees of freedom in individual joint space. To simulate the kinematics of these body segments, individual trajectory tracking controllers were developed; further details are provided in Section 2.3. The position and orientation variables were defined using the SNAME notation [40]. The generalized coordinates associated with this entire model are represented by a vector $q$of length 20 composed of the elements:

$$q_{mainbody}={\left[{\phi }_1,{\theta }_1,{\psi }_1,x,y,z\right]}^T,$$

(1)

$$q_{foreflipper}={\left[{\phi }_i,{\theta }_i,{\psi }_i\right]}^T$$

(2)

where i=2,3 representing left and right fore flipper,

$$q_{hindflipper}={\left[{\phi }_i,{\psi }_i\right]}^T$$

(3)

where i=4,5 representing left and right hind flipper,

$$q_{head}={\left[{\theta }_6,{\psi }_6\right]}^T,$$

(4)

$$q_{pelvis}={\left[{\theta }_7,{\psi }_7\right]}^T,$$

(5)

where $\phi ,\theta ,\psi$ (i=1, 2,…,7) represents the rotations about X, Y and Z axes for each body part. Each vector$q$specifies a configuration of 7 bodies. The set of all configurations is a smooth manifold called the configuration space, M. At each point $q\in M$ the time derivative of$q$is the generalized velocity $\dot{q}$.

The time derivative of the generalized coordinates, i.e., the generalized velocity, $\dot{q}$, is related to the quasi-velocity, $p$, using the following equations:

$$p_{mainbody}={\left[p_1,q_1,r_1,u,v,w\right]}^T,$$

(6)

$$p_{foreflipper}={\left[p_i,q_i,r_i\right]}^T$$

(7)

where i=2,3 representing left and right fore flipper,

$$p_{hindflipper}={\left[p_i,q_i,r_i\right]}^T$$

(8)

where i=4,5 representing left and right hind flipper,

$$p_{head}={\left[q_6,r_6\right]}^T,$$

(9)

$$p_{pelvis}={\left[q_7,r_7\right]}^T.$$

(10)

The notation $p_i,q_i,r_i$ (i=1, 2,…,7) represents the time angular velocities about the joint X, Y and Z axis for each body part. The kinematics of this multi-body model are represented by (11):

$$\dot{q}=V\left(q\right)p.$$

(11)

$$q=\left[\begin{array}{c}{q}_{mainbody};q_{leftfore};q_{rightfore};\\ q_{lefthind};q_{righthind};q_{head};q_{pelvis}\end{array}\right],$$

(12)

$$p=\left[\begin{array}{c}{p}_{mainbody};p_{leftfore};p_{rightfore};\\ p_{lefthind};p_{righthind};p_{head};p_{pelvis}\end{array}\right]$$

(13)

where q is the generalized coordinates vector and p is the quasi-velocity vector, each of length 20. The matrix $V\left(q\right)$ is the velocity transformation matrix, a 20x20 square matrix of full rank. The columns of $V\left(q\right)$, $v_i\left(q\right),i=1,\dots ,20$ form a set of independent vector fields on the configuration manifold M. The construction of these vector fields is based on joint configurations.

Similar to the standard form of dynamical equations used for underwater vehicles [28], the dynamics (Figure 4) of SEAMOUR are represented by (14):

$$M\left(q\right)\dot{p}+C\left(p,q\right)p+F\left(p,q\right)=Q$$

(14)

where M(q) is the mass matrix, the matrix C(p,q) contains centripetal and Coriolis components, F(p,q) is the force vector internal to the system, and Q is the vector of externally applied generalized forces. The complete mass matrix for the entire model, including details of how each body segment is connected to the main body, is provided in the supplementary material. The principal and product moments of inertia were calculated using a computer aided design (CAD) model of SEAMOUR on CREO [Creo Parametric, PTC, Boston, MA, USA]. In addition to these inertial parameters, the total mass of the main body was adjusted to include the additional water mass from the UUV’s flooded section. The CAD model was used to determine the geometrical properties of the main body and its various segments, while each body segment was individually weighed to determine its dry mass. The mass and geometric parameters used in the numerical model of SEAMOUR are summarized in Table I.

2.3. Trajectory Tracking Controllers

Individual tracking controllers were developed to control the fore and hind flippers, head and pelvis with desired trajectories (Figure 5). A two-level strategy was used that first implemented the feedback linearizing control and then linear feedback to regulate the linearized system [24]. For each degree of freedom (roll, pitch and yaw) in the foreflippers, individual trajectory tracking controllers were developed. For each of the head, pelvis and the hind flippers, a 2DoF tracking controller was developed. A detailed description on model setup and the use of ProPac toolbox is described in [24].

2.4. Hydrodynamic Forces

Hydrodynamic forces such as drag, lift, lateral and added mass forces were incorporated into the numerical model to simulate SEAMOUR swimming and maneuvering in water. As shown in (14), which defines the equations of motion, vector Q represents the total external hydrodynamic forces acting on the system. These forces include the effects of drag, lift, and lateral forces, applied to each body segment. Drag forces in the main body were modelled as resistive forces acting opposite to the X-axis (axial direction) as shown in Figure 3a. Lateral forces are applied along the Y axis and lift forces are applied along the Z axis. Using the Generalized Force function within the ProPac toolbox, hydrodynamic forces were computed and applied at the center of mass nodes for the foreflippers, hind flippers, head, and pelvis. For the main body, drag and added mass forces in the X direction were applied to its center of mass. To account for changes to relative velocities across the body in Y and Z directions, the entire model was subdivided into eight segments, including the head, pelvis, and six sections along the main body. Forces were applied to the center of mass of each segment and scaled according to their proportional surface area in the relevant direction. The force coefficients for drag, lift, and lateral components were obtained through computational fluid dynamics (CFD) simulations, while the added mass coefficients were computed analytically using strip theory, applied individually to each body segment.

Some assumptions were taken to simplify the computation of the hydrodynamic coefficients. The main body was considered a prolate spheroid and the foreflippers were considered rectangular flat plates. Water was assumed to be inviscid, incompressible and irrotational-an assumption commonly adopted in prior studies, such as in [28]. It was also assumed that the model has three planes of symmetry with the model always completely submerged in water. In addition, the center of mass and center of buoyancy for each flipper as well as the head and pelvis sections are coincident. All body segments were assumed to be individually neutrally buoyant.

Computational fluid dynamic (CFD) simulations were developed to determine the force coefficients in all directions for each body segment of SEAMOUR (Figure 6). The CAD geometry of the main body was directly imported into COMSOL Multiphysics, where a user-controlled mesh was utilized, with the element size calibrated for fluid dynamics (Figure 6a). A k–ε turbulence model was used, with an inflow velocity of 10 m/s applied across the model domain. To calculate the force coefficients, the line integral of the total stress along the boundaries of each geometry was evaluated, yielding force coefficients along the X, Y, and Z axes (Figure 6b). For the X-direction force coefficient of the main body, the full CAD assembly was imported into COMSOL. For force coefficients in the y and z directions, each body segment was imported individually (Figure 6c), and CFD simulations were run separately to obtain the segment-specific coefficients. Similarly, for both the foreflippers and hind flippers, each flipper’s geometry was analyzed independently to compute force coefficients in all three directions. Drag forces acting on the main body, head, pelvis, and flippers were calculated using (15). Additionally, the velocities at the center of mass nodes were sensed and subsequently used to compute forces at the center of mass nodes of each flipper as well as body segment. NodeVelocity function in ProPac was used to sense the velocities at every node of the numerical model [24]. More information on this can be found in the attached supplementary material. The drag force equation is

$$F_D={\displaystyle \frac{C_D\rho Av\left|v\right|}{2}}$$

(15)

where F_D is the drag force exerted, C_D is the coefficient of drag, ρ is the density of water, A is the reference area and v is the velocity of the object relative to the fluid.

To compute the added mass forces, it was important to compute the added mass coefficients which are used in (22) [28]. Due to the symmetrical shape of the model, the added mass coefficients related to the product of two axes were ignored. To calculate these coefficients for the main body, [37] provided formulas that utilize the geometrical properties of the body. The eccentricity of the prolate spheroid, and the added mass coefficients for a prolate spheroid is given below:

$$e=\sqrt{1-{\left({\displaystyle \frac{b}{a}}\right)}^2},$$

(16)

$${\alpha }_0={\displaystyle \frac{2\left(1-e^2\right)}{e^3}}\left({\displaystyle \frac{1}{2}}ln{\displaystyle \frac{1+e}{1-e}}-e\right),$$

(17)

$${\beta }_0={\displaystyle \frac{1}{e^2}}-{\displaystyle \frac{1-e^2}{2e^3}}ln{\displaystyle \frac{1+e}{1-e}},$$

(18)

$$X_{\dot{u}}=-{\displaystyle \frac{{\alpha }_0}{2-{\alpha }_0}}m,$$

(19)

$$Y_{\dot{v}}=Z_{\dot{w}}=-{\displaystyle \frac{{\beta }_0}{2-{\beta }_0}}m.$$

(20)

To compute the hydrodynamic added mass coefficients for the fore and hind flippers, strip theory was employed. This approach involved dividing the submerged shape into strips, computing the two-dimensional cross-sectional coefficients, and integrating them along the shape’s length to determine the three-dimensional coefficients (Figure 6c). These added mass coefficients were computed based on methodologies presented by [47], and (21) was used:

$$M_a=k_{trans}{\displaystyle \frac{\pi \rho c^{2}b}{4}}$$

(21)

where K_trans is the coefficient of additional mass as a function of aspect ratio of rectangular flat plate in [Figure 3, 31], where $\rho$ is the density of water, c is the chord and b is the span.

$$Addedmass=M_a.a$$

(22)

where a is the acceleration sensed at the same node. In the Simscape multibody environment, a transform sensor block was placed between the body frame and the world frame to sense accelerations in all axes at various nodes. Coefficients of force for the main body as well as, the front face of the foreflippers is presented in Table 1.

A genetic algorithm was used to obtain more accurate estimates of the hydrodynamic coefficients and the volume of water ingress within SEAMOUR. The calculation of hydrodynamic derivatives using empirical formulas is particularly challenging for robots with complex shapes, as existing formulas in literature are typically based on simplified geometries. The optimized parameters included the volume of water inside SEAMOUR, the added mass coefficient in the z direction, the buoyancy offset in the main body, and the drag coefficient in the z direction. To address these challenges, an optimization process focused on minimizing an error function based on Dynamic Time Warping (DTW) using a multi-objective approach was used. The two primary objectives were ensuring accurate alignment between the robot’s position and the numerical model, and consistent pitch angle between the model and the robot for the complete duration during different strokes. The algorithm was configured with population size of 80 and the maximum generations of 50 with a crossover fraction of 0.8 and the pareto function was 0.4. The genetic algorithm was applied to two distinct trajectories: (1) the characteristic stroke, a gait commonly observed in sea lions, and (2) the characteristic stroke with an actuated pelvis. Following the optimization process, the resulting coefficients were fine-tuned manually to identify values that worked best for both trajectories. These refined parameters were then used for the analysis of swimming and maneuvering tests.

2.5. Simscape Visualization

A multibody assembly of the bio-robotic sea lion was developed using MATLAB’s Simscape environment to visualize the 6 DoF body movements (Figure 4). This assembly served as a dynamic representation of the robot’s kinematics. Each anatomical section of the bio robotic sea lion – such as the main body, head, pelvis, fore and hind flippers was modeled in CREO and imported into the Simscape environment. Each body segment was connected to the main body with revolute joints equivalent to its degrees of freedom. The main body of the model was connected to the simulation environment through a 6DoF joint to visualize translation and rotation of the entire assembly. Outputs from the Poincare model were directly integrated into the Simscape model for visualization.

2.6. Model Validation Experiments and Data Processing

Experiments were conducted in a pool to validate the numerical model against SEAMOUR. First, to ensure that various body segments in the numerical model followed the prescribed trajectories, the performance of the individual trajectory tracking controllers was evaluated by analyzing the differences between the input and output trajectories. Following this, a characteristic stroke of a sea lion, as previously described in [28], was implemented on both SEAMOUR and the numerical model. The same stroke was then tested with additional pelvis actuation to compare the response between the two systems. Lastly, simultaneous actuation of the head and pelvis sections in both pitch and yaw was conducted to validate the numerical model against SEAMOUR during maneuvers.

To quantify the kinematic validation of the center of mass of the numerical model against the robot, several metrics were defined. These included: total distance travelled in surge and heave directions, maximum and mean resultant velocity magnitude for every cycle within a stroke, final pitch angle after very cycle, and maximum and mean angular velocity during each cycle of a stroke. For the maneuver tests, final orientation as well as, maximum and mean angular velocity at the end of four seconds was recorded. For the evaluation of the performance of the tracking controller, rise time and settling time is used to describe the differences between the input and output trajectories.

SEAMOUR was tested in a pool where videos were recorded to capture its position and orientation across different trials. To ensure a consistent starting position and orientation, the robot was magnetically attached at its midsection to an aluminum docking structure mounted against one wall of the pool. The docking system included a magnetic end effector with prongs to secure the robot prior to each trial. All experiments were conducted at a depth of approximately 1.22 meters. Video footage was captured using a GoPro Hero 13 (GoPro, San Mateo, CA, USA) mounted on the adjacent pool wall, recording at 30 frames per second in 4K resolution. A reference scale with a known length was used in the field of view to convert pixel measurements to meters. Processed clips were then imported into MATLAB, where the Image Processing and Computer Vision Toolboxes were used to track the robot’s motion. Two visual markers (Figure 7a), placed on the head and pelvis shells, were used to extract position and pitch angle data over the duration of each trial. During maneuver tests, orientation data were also recorded by an onboard inertial measurement unit (BNO085, PN 4754, Adafruit, USA), at a sampling rate of 20 Hz. Each maneuver trial began with three paddle strokes to bring SEAMOUR upto speed, followed by simultaneous actuation of the head and pelvis sections. The foreflippers were kept in a streamlined configuration and were not actuated during maneuvers. Although the foreflippers used for these tests were flexible, and similar to those described in [43], they were treated as rigid for the purposes of this analysis. Since they were not actuated during the trials, their influence on overall performance was considered negligible. Separate video recordings were made for different types of maneuvers. For pitch trials (Figure 7b), the camera was mounted on the side wall, and for yaw trials (Figure 7c), a camera was placed on the bottom of the pool. To ensure fair comparison with numerical simulations, each simulation was initialized with the same pre-actuation velocity measured from experiments (approximately 0.21 m/s). Three trials were conducted for each test condition, and the results were averaged to quantify the robot’s position and orientation.

The numerical model was solved within the MATLAB Simscape environment. A variable-step ODE45 (Dormand-Prince) solver was used to solve the differential equations. The maximum step size was taken to be 0.1 with an automatically determined minimum step size depending on the trial. The relative tolerance was set at 1e-3. A detailed description of the modeling methodology can be found in [48].

3. Results

3.1. Effectiveness of Trajectory Tracking Controller

Across different input trajectories, the tracking controllers consistently generated outputs that closely aligned with the desired input signals. The flipper roll trajectory was tracked near perfectly (Figure 8a). When the pelvis pitch trajectory was implemented, the tracking controller achieved a 2% settling time of 0.58 seconds (Figure 8b). Furthermore, for a steady-state input of zero with an initial condition of 0.6 radians, the controller exhibited a 2% settling time of just 0.017 seconds (Figure 8c), demonstrating strong agreement between input and output.

3.2. Translation for Characteristic Stroke

The numerical model closely matched the bio-robotic system's surge translation but exhibited a larger deviation in heave. When the characteristic stroke was implemented on the bio-robotic system, it translated in surge to 1.79m. When implemented on the numerical model, it translated to 1.80m which was only 0.56 % higher than the robotic system (Figure 9a). In heave, the variation was higher. The bio-robotic sea lion translated to 0.73m at the end of three cycles. The numerical model translated 0.57m in heave which was approximately 22% lower than the robot.

3.3. Velocity Magnitude for Characteristic Stroke

The mean resultant velocity magnitude increased progressively for each cycle for both the bio-robotic system and the numerical model (Figure 9b). For example, when the characteristic stroke was implemented on the numerical model, the mean of the resultant velocity magnitude for the first cycle was 0.09m/s. This increased by about 167% to 0.24m/s in cycle 2 and 256% to 0.32m/s when compared to cycle 1 (Table II). Similarly, for the robotic system, the mean for the first cycle was 0.15m/s. For the second and third cycle the mean was 0.26 m/s and 0.27 m/s which was about 73% and 80% higher respectively when compared to cycle 1. When comparing the two systems, for cycle 1, the robotic system had almost 67% higher mean than the numerical model. For cycle 2, the mean for both systems was much closer with the robotic system being only approximately 8% higher. For the third cycle the mean for the numerical model was approximately 16% higher when compared to the robotic system.

The maximum resultant velocity magnitude increased for each cycle for the numerical model, but for the robotic system, it increased from cycle 1 to cycle 2 and slightly reduced in cycle3. Specifically in the numerical model, when characteristic stroke was applied, the maximum resultant velocity magnitude for the first cycle was 0.21m/s. This increased by approximately 52% to 0.32m/s in cycle 2 (Table II). For cycle 3, the maximum was 0.39m/s which was about 86% higher compared to cycle 1. For the robotic system, the maximum was 0.33m/s in cycle 1. In cycle 2 and 3, the maximum increased by approximately 12% and 9% to 0.37 m/s and 0.36m/s respectively when compared to cycle 1. Additionally, when the numerical model is compared to bio robotic system, for cycle 1, the robotic system had higher maximum velocity by approximately 57%. However, for cycle 2, the max velocity for the robot was much closer to the numerical model and only 11% higher. Like cycle 2, the maximum velocities were close for both systems with the numerical model only 8% higher than the robot.

Table II. MAXIMUM AND MEAN RESULTANT VELOCITY MAGNITUDE FOR EACH CYCLE FOR ROBOT AND SIMULATION DURING CHARACTERISTIC STROKE (CS) AND CS WITH PELVIS ACTUATION.

3.4. Orientation for characteristic stroke

Final pitch angle increased progressively after each cycle for characteristic stroke for both the bio robotic system as well as the numerical model (Figure 10a). Specifically, for the numerical model, final pitch after cycle 1 was 12°. The pitch angle increased by approximately 158% to 31° at the end of cycle 2. At the end of cycle 3, the pitch angle further increased by 11% to 35° (Table III). Similarly, for the robot, the final pitch angle after cycle 1 was 17°. The pitch angle increased by approximately 71% to 29° at the end of cycle 2. The pitch angle further increased by approximately 45% to 42° at the end of cycle 3. When comparing the robot and the model, at the end of cycle 1, the robot’s final pitch angle was approximately 42 % higher than the model. However, at the end of cycle 2, the final pitch angle for both systems was much closer and only about 7% higher in the numerical model. At the end of cycle 3, the robot’s final pitch angle was about 35% higher than the model.

Table III. ANGULAR POSITION AND ANGULAR VELOCITY FOR EACH CYCLE FOR THE ROBOT AND THE SIMULATION DURING CHARACTERISTIC STROKE (CS) AND CS WITH PELVIS ACTUATION.

3.5. Angular Velocity for Characteristic Stroke

The mean angular velocity followed an inconsistent trend between the numerical model and the robotic system across each cycle when the characteristic stroke was executed (Figure 10b). Specifically, in the numerical model, the mean angular velocity increased from cycle 1 to cycle 2, then decreased in cycle 3, while the robotic system showed the opposite pattern. The mean angular velocity during cycle 1 in the numerical model was 5°/sec. During cycle 2, the mean angular velocity increased by 20% to 6°/sec (Table III). In cycle 3, the mean angular velocity decreased to 2°/sec. For the robotic system, the mean angular velocity was 6°/sec. During cycle 2, this decreased by about 33% to 4°/sec. In cycle 3, the mean again increased by 25% to 5°/sec. When comparing both systems, in cycle 1, the mean was relatively similar with the robot being 20% higher. In cycle 2, the mean angular velocity of the model was 50% higher. During the third cycle, the mean angular velocity of the robot was significantly higher than the numerical model, almost by 150%.

The maximum angular velocity for each cycle decreased for the numerical model, but for the robotic system, it decreased from cycle 1 to cycle 2 and slightly increased in cycle 3. For the numerical model, the maximum angular velocity for cycle 1 was approximately 11 °/sec. The max angular velocity decreased to approximately 9°/sec, which was approximately 18% lower. The max velocity further decreased by about 56% to 4°/sec. For the robotic system, the max angular velocity in cycle 1 was 23°/sec. The max angular velocity decreased in cycle 2 and cycle 3 to about 15°/sec. For each cycle, the maximum angular velocity was consistently higher in the robot compared to the numerical simulation.

3.6. Translation for Characteristic Stroke with Pelvis Actuation

Similar to the results of the characteristic stroke, when the pelvis was actuated with the characteristic stroke on both systems, the numerical model closely matched translation in surge with a larger deviation in heave. The numerical model translated 2.10m in surge at the end of three cycles (Figure 11a). When the same stroke was applied on the robotic system, it translated 2.09m in surge which was approximately only 0.47%lower. In contrast, the numerical model traversed 0.25m in heave compared to no translation in heave for the robot at the end of three cycles.

3.7. Velocity Magnitude for Characteristic Stroke with Pelvis Actuation

The mean resultant velocity magnitude increased for every cycle when the characteristic stroke with pelvis actuation was implemented on both the numerical model as well as the bio-robotic system. Specifically, for the numerical model, the mean resultant velocity in the first cycle was 0.11m/s. For the second cycle, the mean increased by approximately 146% to 0.27m/s. This mean velocity further increased by 33% to 0.36m/s (Table II). Similarly, for the robotic system, the mean resultant velocity in cycle 1 was 0.16m/s. The mean resultant velocity increased by approximately 69% to 0.27m/s in cycle 2. This further increased by only 7% to 0.29m/s in cycle 3. When comparing the robotic system to numerical model, the mean resultant velocity for robotic system in cycle 1 was approximately 46% higher. For cycle 2, both systems had similar mean values. However, in cycle 3, the mean resultant velocity in numerical model was only approximately 24% higher than robotic system.

Like the mean velocity magnitude, the maximum resultant velocity magnitude increased progressively for every cycle when characteristic stroke with pelvis was implemented on both bio robotic system as well as the numerical model. For instance, in the numerical model, the maximum velocity in cycle 1 was 0.23m/s (Table II). The maximum velocity increased by approximately 52% to 0.35m/s. This further increased by 20% to 0.42m/s. Similarly, for the robotic system, the maximum velocity in cycle 1 was 0.35m/s. The mean value increased by approximately 9% to 0.38ms/s in cycle 2 and cycle 3. When comparing the robotic system to the numerical model, the maximum resultant velocity in cycle 1 was 52% higher. However, the maximum velocities were much closer in cycle 2 for both systems, with maximum velocity of robot being approximately 8% higher. During cycle 3, the maximum velocity in the simulation was approximately 11% higher than the robotic system.

3.8. Orientation for Characteristic Stroke with Pelvis Actuation

The final pitch angle decreased progressively after each cycle when the characteristic stroke with pelvis actuation was implemented on the bio-robotic system but first increased and then decreased for the numerical model. In the numerical model, the final pitch angle after cycle 1 was 2° (Figure 12a). At the end of cycle 2, the final pitch angle increased by 50% to 3°. The final pitch angle decreased by about 166.70% going from cycle 2 to cycle 3 to approximately -2° (Table III). For the robotic system, at the end of cycle 1, the final pitch angle was 8°. The final pitch angle decreased by approximately 25% and 88% in cycle 2 and cycle 3 respectively. When comparing both systems, at the end of cycle 1,2, and 3, the numerical model was about 75%, 50%, and 300% lower respectively.

3.9. Angular Velocity for Characteristic Stroke with Pelvis Actuation

The mean angular velocity progressively decreased for every cycle for the robotic system, but for the numerical model, it slightly increased from cycle 1 to cycle 2 and then decreased from cycle 2 to cycle 3 (Figure 12b). For the numerical model, the mean angular velocity was approximately 0.8°/sec for cycle 1 and 0.9°/sec for cycle 2 which was approximately 13% higher. However, going from cycle 2 to cycle 3, the mean angular velocity decreased by about 255% to -1.4°/sec. For the robotic system, during cycle 1, the mean angular velocity was approximately 2.7°/sec. During cycle 2, the mean angular velocity decreased by about 133% to -0.9°/sec. This further decreased by about 77% to -1.6°/sec. For cycle 1, the mean angular velocity of the numerical model was 66% lower. During cycle 2, the mean angular velocity for the numerical model was about 200% higher (Table III). The mean angular velocity for the numerical model and the robotic system was much closer with the numerical model only about 13% lower. For both the numerical model and the robotic system, the maximum angular velocity consistently decreased from cycle 1 to cycle 2 and then increased from cycle 2 to cycle 3. In the numerical model, the maximum angular velocities across each cycle were within 8% of each other. Similarly, in the robotic system, the maximum angular velocities were within 6% of each other. During cycle 1 for the numerical model, the maximum angular velocity was 24°/sec. The max angular velocity decreased by 8% to 22°/sec in cycle 2 and then increased back to 24°/sec in cycle 3. Similarly, in the robotic system, the maximum angular velocity in cycle 1 was 18°/sec. The max angular velocity decreased by 6% to 6°/sec during cycle 2 and then increased back to 18°/sec in cycle 3. When compared to the numerical model, the maximum angular rate for the robotic system across each cycle was about 24% lower.

3.10. Head and Pelvis Pitch Maneuver

The numerical model’s prediction for the final pitch angle, as well as, maximum and mean angular velocity across the four seconds during a pitch maneuver very closely matched the bio robotic system. The numerical model predicted a final pitch angle of 44.27° at the end of four seconds. On the other hand, the bio robotic system produced a total pitch angle of 42.40° which was only approximately 4% lower (Figure 13a). In addition, the maximum and mean pitch velocity across the four second duration for the numerical model was 41.17°/sec and 11.27°/sec (Figure 13b). When compared to the numerical model, the maximum pitch velocity for the robot was only about 4% higher. Additionally, the mean pitch velocity for the robot was only approximately 6% lower when compared to the numerical model.

3.11. Head and Pelvis Yaw Maneuver

Like the pitch maneuver, when head and pelvis were actuated for a yaw maneuver, the final yaw angle, and the maximum and mean angular velocity across the four second duration, showed very close alignment between the numerical model and the robotic system. The final yaw angle for the numerical model at the end of four seconds was 35.84° (Figure 14a). For the robotic system, the final yaw angle was 34.80° which was just 3% lower. The maximum and mean angular velocity across four seconds in the numerical model was 24.67°/sec and 8.10°/sec respectively (Figure 14b). Similarly, the maximum angular velocity across the total duration I the robotic system was 24.63°/sec which was less than 2% lower. The mean angular velocity or the robotic system was 8.70°/sec which corresponded to an approximate increase of about 7% when compared to the numerical model.

4. Discussion

This study successfully developed and validated a numerical model of a bio-robotic sea lion using closed-form nonlinear equations of motion. The model achieved over 99% accuracy in the prediction in surge and approximately 85% accuracy in mean angular pitch rates during swimming tests, with larger deviations observed in heave. During maneuver tests in pitch and yaw directions, where the head and pelvis were actuated simultaneously, the model closely aligned with experimental results, achieving approximately 90% accuracy in predicting rotational rates and final orientation. Crucially, the equations of motion for the model incorporated all major body segments ensuring a comprehensive representation of the system's dynamics. Unlike prior numerical models that are developed for uniform segmented vehicle swimming in planar motion, this numerical model captured the three-dimensional hydrodynamic forces produced by non-uniform body segments in a UUV.

Initial efforts to align the numerical model with the bio-robotic system used traditional strip theory methods to estimate hydrodynamic added mass coefficients. However, this approach led to limited agreement between the simulation and experimental results. The robot’s main body was initially approximated as a prolate spheroid with equal minor axes to simplify the computation of added mass coefficients. The cross-sectional profile of the body was elongated along the Y-axis relative to the Z-axis, introducing asymmetries that strip theory failed to model accurately.

To address these limitations, a genetic algorithm was employed to optimize key hydrostatic and hydrodynamic parameters. This approach enabled refinement of the hydrodynamic coefficients—particularly along the Z-axis—to better capture the true geometry and dynamic response of the robot. In addition, because SEAMOUR was built for water going in, the internal water volume was treated as a variable in the optimization process. Analytical modeling of water ingress and egress in SEAMOUR poses a significant challenge, as the amount of water entering the system varies with the robot’s speed and orientation. The genetic algorithm proved effective in tuning these complex parameters, resulting in a significantly improved alignment between the numerical model and the experimental observations.

While the numerical model demonstrated high predictive accuracy for surge and rotational maneuvers, the primary discrepancies were observed in heave and pitch dynamics. These deviations likely resulted from differences in buoyancy distribution between the simulation and the physical robot, internal fluid movement within SEAMOUR during trials, and the influence of cable dynamics, which were not included in the simulation. In the numerical model, each body segment except the main body was assumed to be neutrally buoyant, with coincident centers of mass and buoyancy. In contrast, SEAMOUR’s head and pelvic sections, although neutrally buoyant, had non-coincident hydrostatic centers, which affected the restoring forces acting on the system. Over time, the physical robot also tended to become negatively buoyant as internal components gradually absorbed water. Furthermore, for numerical simplicity, the flippers and motors were modeled as neutrally buoyant, whereas in the physical system, they were negatively buoyant, contributing to shifts in the overall center of mass during different experimental trials. The simulation also excluded the internal movement of water, which likely varied with SEAMOUR’s speed and orientation and influenced its dynamic response. Lastly, the physical robot included a trailing antenna cable with a small surface float that was not represented in the model; however, their effect on the overall dynamics was likely marginal compared to the more significant contributions of buoyancy distribution and internal water movement.

5. Conclusions

In this study, a numerical model of a bio-robotic sea lion was successfully developed and validated. The kinematic validation of the center of mass of the numerical model against the experimental platform demonstrated high accuracy in the prediction of position and orientation over time during free-swimming and maneuver trials. The modeling framework successfully derived closed-form equations of motion for the multi-body system, enabled visualization of full-body kinematics, and applied estimated fluidic forces at the center of mass of each segment. While the current study did not include a comparison between simulated joint torques and actuator torques measured from the physical robot, incorporating this form of dynamic validation, along with fluid–structure interaction effects, represents an important next step to further improve model fidelity and dynamic accuracy.

The derivation of equations of motion for a multi-body UUV in closed-form has potential benefits beyond the scope of this study, particularly in control theory, stability analysis, and reinforcement learning applications [23,49,50,51,52].The closed-form equations of motion provide explicit expressions for the mass matrix, coriolis and centripetal terms, and generalized force matrices, that offer deep insight into the system’s inherent dynamics [24,53]. Such a formulation enables the application of advanced model-based control strategies, including feedback linearization, adaptive control, and sliding mode control [24]. Moreover, these closed form dynamical equations facilitate the construction of Lyapunov functions which are critical for formal stability analysis [24].Furthermore, the model serves as an effective training environment for reinforcement learning algorithms and has previously been used to generate a novel gait for straight swimming [43].

The hydrodynamic coefficients, which were critical for estimating fluid forces in the numerical model, were initially underestimated when using strip theory calculations and computational fluid dynamics (CFD) simulations. This discrepancy was likely due to the complex geometry of the robot and the influence of internal water movement. To address this, optimization technique such as a genetic algorithm was employed to refine the coefficients, which resulted in an improved alignment between the numerical model and the experimental platform.

The closed-form equations of motion developed for the multi-body UUV in this study provide a strong foundation for future work in both analysis and control of underwater robotic systems. Future research will focus on conducting stability analysis using techniques such as Lyapunov methods and eigenvalue analysis of linearized dynamics. The model also enables the development of advanced model-based controllers, as well as reinforcement learning approaches that can be trained in simulation to optimize swimming gaits and maneuvering strategies. To further enhance model fidelity, fluid–structure interaction effects and the impact of ingress and egress should be investigated. Static stability tests along the X and Y axes should be conducted on the bio-robotic swimmer, and the hydrostatic centers in the numerical model should be adjusted accordingly. This would help align the buoyancy and mass distribution between the numerical model and the physical system, leading to improved accuracy in the numerical simulation.

Further experimental testing of the robotic platform with various combinations of actuated control surfaces, coupled with iterative refinements to the simulation, will enhance the model’s ability to accurately predict system performance. The modeling approach presented in this study offers a versatile framework that can be extended to other multi-body underwater robotic systems.