Modeling, Trajectory Analysis and Waypoint Guidance System of a Biomimetic Underwater Vehicle Based on the Flapping Performance of Its Propulsion System

Abstract

The performance of biomimetic underwater vehicles directly depends on the correct design of their propulsion system and its control. These vehicles can attain highly efficient motion, hovering and thrust by properly moving part(s) of their bodies. In this article, a mathematical modeling and waypoint guidance system for a biomimetic autonomous underwater vehicle (BAUV) is proposed. The BAUV achieves sideways and dorsoventral thunniform motion by flapping its caudal fin through a parallel mechanism. Also, an analysis of the vehicle’s design is presented. A thrust analysis was performed based on the novel propulsion system. Furthermore, the vehicle’s kinematics and dynamic models were derived, where hydrodynamic equations were obtained as well. Computed models were validated using simulations where thrust and moment analysis was employed to visualize the vehicle’s performance while swimming. For the path tracking scheme, a waypoint guidance system was designed to correct the vehicle’s direction toward several positions in space. To accurately obtain waypoints, correction over the propeller’s flapping frequency and bias was employed to achieve proper thrust and orientation of the vehicle. The results from numerical simulations showed how by incorporating this novel propulsion strategy, the BAUV improved its performance when diving and maneuvering based on the dorsoventral and/or sideways configuration of its swimming mode. Furthermore, by designing proper strategies to regulate the flapping performance of its caudal fin, the BAUV followed the desired trajectories. The efficiency for the designed strategy was obtained by comparing the vehicle’s traveled distance and ideal scenarios of straight-line trajectories between targets. During simulations, the designed guidance system presented an efficiency of above 80% for navigation tasks.

1. Introduction

Nowadays, biomimetic autonomous underwater vehicles (BAUVs) represent an alternative strategy to accomplish navigation tasks without disrupting the natural cohesion in aquatic ecosystems. These vehicles deliberately imitate a diverse set of locomotion skills found on fish to enhance their underwater performance. One important reason for this is that fish excel in propelling their body using only their fins [1]. Also, compared to classic nonbiomimetic systems, BAUVs have higher levels of maneuverability and efficiency, outperforming traditional propulsion systems. The capacity to generate thrust and moment without the use of noisy and energy-expensive turbines is one of the many advantages of these mechanisms [2].

Like fish, BAUVs move a part (or parts) of their bodies to generate a pressure difference in water in order to achieve momentum. However, to attain undulatory or oscillatory propulsion, swimming styles may vary according to species. In general, two types of swimming methods can be found in nature: one where medium and/or paired fins (MPF) are employed to propel, and the other where the body and/or a caudal fin (BCF) is moved to generate thrust [3]. BCF swimmers comprise approximately 85% of fish species, including many fast swimmers such as sailfish, tuna, and pike [4]. For the BCF swimming style, a more detailed classification may be made, dividing roughly into anguilliform, subcarangiform, carangiform, and thunniform locomotion. While other fish undulate their whole (anguilliform) or most of their bodies (carangiform and subcarangiform), thunniform swimmers generate thrust by only flapping their caudal fins, making this the best swimming strategy in terms of efficiency [5].

The performance of BAUVs with BCF locomotion will directly depend on the mechanics employed inside their propellers. Specifically, their efficiency, speed, thrust generation and maneuverability will be limited by the way their propulsion system is configured. Many novel designs have been developed to properly generate oscillations to the last section of the bodies of these vehicles. Some of them are based on multijoint robotic propellers, where several rotary actuators are serially linked to imitate fish-like vertebrae and carangiform locomotion. By generating rhythmic movements from joints, sideways or dorsoventral flapping might be achieved inside the propulsion system [6,7]. Also, to improve turning without taking speed into account, some larger configurations might be implemented. Among these cases, configurations of four or more actuators may be employed to attain a larger range of bias from the moving part of the body and/or caudal fin while flapping [8]. On the other hand, and to reduce energy consumption, some designs focus on only producing motion on the caudal fin. This way, by limiting the number of actuators and with the help of transmission systems to enhance the flapping frequency from the vehicle’s tail, faster robotic fish may be designed [9,10]. Then, based on how the propulsion system was configured and the achievable swimming fashion, some tradeoffs in the performance of an underwater vehicle should be considered. On one hand, high energy consumption mechanisms may attain lower cruising speeds but excel at maneuverability, while highly energy efficient mechanisms may produce faster speeds but present relatively poor maneuverability [11].

Regardless of the locomotion achieved, a BAUV position and orientation over time will depend on how fast (frequency), how large (amplitude), and how biased the propeller’s flapping can become, and how its motion is regulated. Also, the trajectory described by the vehicle will be conditioned by the forces and moments exerted by the propulsion system’s motion, the added mass effect produced, and the hydrodynamics presented by the vehicle [9]. Hence, to correctly implement guidance strategies in the path tracking task of these vehicles, it is necessary to understand how all these factors will affect performance. A proper derivation of the kinematics and dynamic modeling of the vehicle, a force and moment analysis of the propulsion system, and open-loop validations are required to identify which strategies may overcome unwanted reactions.

Several studies have focused on the derivation of the kinematics and dynamic models of biomimetic underwater vehicles. These studies are conducted to identify which behaviors the vehicle will tend to present while swimming. In [12], Ravichandran et al. performed numerical simulations based on the dynamic modeling of a REMUS autonomous underwater vehicle to which a caudal fin was fixed as the main propeller. By changing its orientation, the caudal fin was able to reproduce dorsoventral and sideways flapping. However, no information of the mechanical design required for the propulsion system to attain such configurations was presented, and a force and moment analysis was not conducted. They found that based on the swimming fashion, the vehicle would tend to follow a circular path while flapping sideways, and straight lines while flapping dorsoventrally. Szymak [13] developed a mathematical model for a BAUV and a forces and moments analysis produced by a sideways undulating propeller. In his design, the vehicle’s movable tail was considered for thrust generation and two independent pectoral fins were employed to change the cruising orientation. The vehicle achieved a proper thrust once the tail flapping frequency reached 2 Hz at a total amplitude of 20$°$. During the open-loop simulations, the trajectories presented circular shapes, and the turning direction could be changed when pectoral fins were activated. Majeed and Ali [14] detailed a mathematical model of a carangiform BAUV. For the vehicle’s dynamics, hydrodynamic added mass and damping terms were computed by approximating the vehicle’s shape to a prolate ellipsoid. Also, to estimate the thrust exerted by the caudal fin, the authors used Lighthill’s Large Amplitude Elongated Body Theory [15], and the heave and pitch from the vehicle were assumed to be sinusoid functions. The vehicle presented circular paths during open-loop simulations with a tail-sideways-flapping frequency of 1 Hz at a total amplitude of 42$°$. The authors also demonstrated that biasing the caudal fin’s flapping amplitude to one side made it possible to change the forward direction of the vehicle.

Guiding a BAUV toward a goal is a high priority task. Therefore, a guidance system is required to process the navigation and trajectory presented by the vehicle. This system considers the BAUV’s velocity and attitude, and information regarding the desired goal. Hence, it defines the correction that must be performed to keep on the correct track. Additionally, a control scheme is used to generate signals on actuators to produce the required correction that will help the vehicle swim toward set points. Based on its biomimetic features, the vehicle’s swimming performance will mainly depend on the motion of the propeller and the ways in which forces and moments are exerted.

In the waypoint tracking strategy, the vehicle should reach a preestablished set of coordinates in space (waypoints). The line of sight (LOS) algorithm is one common guidance strategy, where the underwater vehicle maneuvers itself to correct its direction and arrive at predetermined intermediate targets. This is usually done by using coordinate transformations where the desired and actual position of the vehicle are iteratively compared to reduce the error distance, and to compute the desired orientation [16]. Then, for the complementary control scheme, strategies such as dynamic and robust sliding mode methodologies may be employed to produce desired corrections [16,17]. Guo [18] designed a path tracking controller for a carangiform robotic fish. Inside the control loop, the vehicle’s forward speed and average heading error were fed back to correct the flapping frequency and the bias from the robotic fish tail. By offsetting the tail undulation, turning was achieved and the orientation of the vehicle was iteratively adjusted. Kopman et al. [19] implemented a waypoint tracking controller based on the heading control of a small two-link robotic fish. The vehicle followed a predetermined path approximated by waypoints positioned close to each other. Some other authors have successfully implemented this strategy by combining control over pectoral fins attached to the vehicle. Wang et al. [20] controlled the path followed by their robotic ray by producing forward thrust with a tail fin and correcting its direction using its two pectoral fins, accounting for the vehicle’s speed, depth, and course control.

In the present research, a waypoint guidance system for a BAUV is developed and simulated. The BAUV presents a vectored thruster with a novel propulsion system designed by Aparicio et al. [21]. The vehicle follows thunniform locomotion and is driven by a caudal fin. The tail’s flapping frequency, amplitude and offset are regulated by a three-degrees-of-freedom parallel mechanism [22]. This novel propulsion system allows the caudal fin to change its orientation via sideways and dorsoventral flapping configurations. Then, by only controlling linear actuators, the tail can oscillate in different dispositions, providing excellent vectored thrust and allowing the BAUV to dive and maneuver while swimming.

The present work develops a method of analysis of force and moment produced by such a propulsion system based on the performance of the designed parallel mechanism and its workspace. Also, the vehicle’s kinematics and dynamic modelling are detailed. Further studies on hydrodynamics based on the hull design of the BAUV are described. Then, by considering such models and the novel feature of switching among swimming modes, the underwater performance of the vehicle is simulated in open-loop trajectories. Finally, based on an analysis of trajectories determined by the vehicle’s biomimetic features, a waypoint guidance strategy is developed. The results from open-loop simulations showed how the vehicle tended to describe defined trajectories according to the swimming style and flapping performance from the fin. However, by maintaining vehicle’s forward speed, and by properly implementing a strategy to correct the fish’s heading angle, the BAUV successfully reached predetermined waypoints on the horizontal plane. Hence, by the adequate definition of flapping parameters such as frequency and bias, thrust and moment are regulated to correct unwanted BAUV behaviors. To the best of the authors’ knowledge, even when parallel mechanisms are incorporated to produce vectored thrust on nonbiomimetic underwater vehicles (with turbine-based-propellers) [23,24,25], their incorporation to change biomimetic features (i.e., swimming mode) on BAUVs has not previously been developed. The means by which this novel feature enhances the cruising trajectories of BAUVs is explained in this paper.

This article is divided as follows: Section 2 details the methods considered during the development of the research. The BAUV design, and mathematical definitions of its hull, kinematics, dynamics, and hydrodynamics are developed. Further, the waypoint guidance strategy and all considerations for the implementation of simulations are detailed. In Section 3, the results from numerical simulations are presented. A discussion of the results obtained for thrust and moment generation, open-loop trajectories and the incorporated guidance strategy based on the flapping performance of the propeller is presented. Finally, Section 4 includes final remarks and the conclusion of this paper, as well as explaining some future work that this investigation may lead to.

2. Materials and Methods

2.1. Vehicle Design

The design of the proposed biomimetic underwater vehicle is presented on Figure 1. Two rigid pectoral fins are positioned at the vehicle’s midbody section for stabilization purposes. In the vehicle’s last section, the propulsion system, consisting of a parallel mechanism, is found. A lunate-shaped caudal fin is attached on the peduncle, serving as the main system propeller. The vehicle’s hull shape, presented in Figure 2, can be divided into four parts:

• The nose section, which presents the form of a Myring hull, is responsible for reducing fluid resistance while the vehicle is swimming [26]. Based on nose geometries from Myring equations [27], the nose will present an initial offset $a_{off}$, since this section will house a camera. Mathematically, radius $r$ of the nose section hull, measured normal to the vehicle’s centerline, is represented by:

$$r\left(x_n\right)=\frac{1}{2}d{\left[1-{\left(\frac{x_n+a_{off}-a}{a}\right)}^2\right]}^{\frac{1}{n}};$$

(1)

where $a$ is the nose full length, $x_n$ is the axial position that varies from 0 to $a$, $d$ is the maximum body diameter, and $n$ is the exponential parameter that defines the nose hull shape.

• The midbody section, which contains all the electronics, and where the vehicle’s rigid pectoral fins start, is formed by an acrylic cylinder of diameter $d.$ The total length of the midbody section is defined by $b$ on Figure 2.
• The propulsion system’s hull, which is formed by a silicone cover, houses the parallel mechanism employed to move the caudal fin. It is mathematically defined by the equation:

$$r\left(x_t\right)=\frac{3919}{355}{x}_t^3-\frac{3255}{842}{x}_t^2+\frac{428}{1705}{x}_t+\frac{d}{2};$$

(2)

where $r$is the hull radius, $x_t$ is the axial position along the centerline ranging from the end of the midbody section $b$ until the peduncle for a total length of $c$, and $d$ is the vehicle’s body diameter.

• The peduncle section, which attaches the caudal fin to the parallel mechanism, presents a semicircular shape. Like the nose section, the peduncle presents an offset $p_{off}$ on its last part to correctly connect with the propeller (caudal fin). The peduncle’s hull radius $r$ is defined by:

$$r\left(x_p\right)=p_{off}+\sqrt{c_o^2-{\left(x_p-ac\right)}^2};$$

(3)

where $c_o$ is the inner radius of the peduncle’s shape, $ac$ is the initial length of the section ($a-a_{off}+b+c)$, and $x_p$ is the axial position ranging from $ac$ until the vehicle’s total length $l.$

Table 1. BAUV hull parameters, weight, and buoyancy.

Symbol	Value	Units	Parameter
$a$	$59.00\times {10}^{-3}$	m	Nose section full axial length.
$a_{off}$	$2.20\times {10}^{-3}$	m	Nose section axial offset.
$b$	$318.43\times {10}^{-3}$	m	Midbody section full axial length.
$c$	$230.82\times {10}^{-3}$	m	Parallel mechanism’s hull axial length.
$c_o$	$38.50\times {10}^{-3}$	m	Peduncle’s inner radius.
$d$	$114.3\times {10}^{-3}$	m	Midbody diameter.
$e$	$38.50\times {10}^{-3}$	m	Peduncle’s full axial length.
$p_{off}$	$6.30\times {10}^{-3}$	m	Peduncle’s aft offset.
$l$	$644.6\times {10}^{-3}$	m	Vehicle’s total axial length.
$n$	$2$	-	Nose shape coefficient.
$W$	$39.24$	N	Vehicle’s weight.
$B$	$44.54$	N	Vehicle’s buoyancy.
$\left(x_g,y_g,z_g\right)$	$\left(241.95, 0, 50\right)\times {10}^{-3}$	m	Vehicle’s center of gravity 1.
$\left(x_b,y_b,z_b\right)$	$\left(241.95, 0, 0\right)\times {10}^{-3}$	m	Vehicle’s center of buoyancy 1.

¹ Measured with respect to nose tip.

shows the parameters that were used to mathematically define the vehicle’s hull shape. The vehicle’s weight $W$ may vary based on the components included in its middle section. The distribution of components inside the vehicle may cause variations in its center of gravity $\left(x_g,y_g,z_g\right)$. On the other hand, the vehicle’s center of buoyancy $\left(x_b,y_b,z_b\right)$was computed based on the CAD design of the exterior hull, as shown in Figure 1. Also, both centers are placed on the same axis separated by a metacentric distance to improve stability while swimming [10]. For simplicity purposes, buoyancy, weight, center of buoyancy and center of gravity are considered as constants throughout all analyses presented in this article; see Table 1. Both center of gravity and buoyancy are measured from the vehicle’s nose tip.

2.2. Vehicle Kinematics and Coordinate Systems

To correctly track the vehicle’s position in space, two reference frames were employed: an inertial Earth-fixed frame $\left(x,y,z\right)$ and a body-fixed frame $\left(x_b,y_b,z_b\right)$ [28]. The moving body-fixed frame $\left\{O_B\right\}$ is located at the BAUV’s center of buoyancy, as shown on Figure 3. The vehicle’s position and orientation are referenced to the inertial frame $\left\{O_I\right\}$, while linear and angular velocities, as well as forces and moments, are referenced to the body-fixed frame coordinates.

Moreover, six parameters were used to represent the vehicle’s six degrees of freedom (DOF): parameters $\left(u,v,w\right)$ represent the translational velocities, while $\left(p,q,r\right)$ represent the angular velocities measured on the vehicle respecting to $\left(x_b,y_b,z_b\right)$ coordinates.

Table 2. Nomenclature employed for underwater motion.

DOF	Position 1	Velocity 2	Force/Moment 2
$x$ direction (surge)	x	u	X
$y$direction (sway)	y	v	Y
$z$direction (heave)	z	w	Z
$x$ axis rotation (roll)	ϕ	p	K
$y$ axis rotation (pitch)	θ	q	M
$z$ axis rotation (yaw)	ψ	r	N

¹ Measured with respect to earth-fixed frame. ² Measured with respect to body-fixed frame.

shows the notation and nomenclature based on the SNAME convention for the position, velocity, forces, and moments used for underwater vehicles [29]. The same nomenclature is displayed in Figure 3.

Hence, the general motion of the vehicle with respect to inertial frame, its translational and rotational velocity with respect to its body-fixed frame, and the total forces and moments applied to it will be defined by vectors $\mathit{\eta },\mathit{\upsilon },$ and $\mathit{\tau }$, respectively:

$$\mathit{\eta }={\left[x,y,z,\varphi ,\theta ,\psi \right]}^T;$$

(4)

$$\mathit{\upsilon }={\left[u,v,w,p,q,r\right]}^T;$$

(5)

$$\mathit{\tau }={\left[X,Y,Z,K,M,N\right]}^T.$$

(6)

Transformation matrices ${\mathit{J}}_1\left({\mathit{\eta }}_1\right)$ and ${\mathit{J}}_2\left({\mathit{\eta }}_2\right)$ were used to convert linear and angular velocities on the body-fixed coordinate system to velocities $\dot{\mathit{\eta }}$ on the inertial frame system:

$$\dot{{\mathit{\eta }}_1}={\left[\dot{x},\dot{y},\dot{z}\right]}^T={\mathit{J}}_1\left({\mathit{\eta }}_1\right){\left[u,v,w\right]}^T, \dot{{\mathit{\eta }}_2}={\left[\dot{\varphi },\dot{\theta },\dot{\psi }\right]}^T={\mathit{J}}_2\left({\mathit{\eta }}_2\right)\left[p,q,r\right];$$

(7)

where using Euler convention for rotation angles results in:

$${\mathit{J}}_1\left({\mathit{\eta }}_1\right)=\left[\begin{array}{ccc}c\theta c\psi & s\varphi s\theta c\psi -c\varphi s\psi & c\varphi s\theta c\psi +s\varphi s\psi \\ c\theta s\psi & s\varphi s\theta s\psi +c\varphi c\psi & c\varphi s\theta s\psi -s\varphi c\psi \\ -s\theta & s\varphi c\theta & c\varphi c\theta \end{array}\right];$$

(8)

$${\mathit{J}}_2\left({\mathit{\eta }}_2\right)=\left[\begin{array}{ccc}1& s\varphi t\theta & c\varphi t\theta \\ 0& c\varphi & -s\varphi \\ 0& s\varphi /c\theta & c\varphi /c\theta \end{array}\right];$$

(9)

where $s\varphi$ is short notation for sin$\left(\varphi \right)$, $t\theta$ for tan$\left(\theta \right)$, $c\psi$ for cos$\left(\psi \right)$, and so on. In the vehicle kinematics Equations (8) and (9), it was determined that the BAUV’s pitch angle $\theta$ was small enough and would not reach values of $\pm \frac{\mathsf{\pi }}{2}.$ Further kinematics analyses might be required in cases whereby the vehicle’s motion approaches singularities; nonetheless, the vehicle is designed to avoid extreme pitch angles.

2.3. Propulsion System

A caudal fin, driven by a 3 DOF parallel mechanism, is the BAUV’s main propeller (Figure 1b). These mechanisms are configured with two platforms, i.e., one fixed and one moved by several actuators or limbs. By changing the position and orientation of the moving platform, an end effector can achieve different poses. With proper control over multiple combinations of dispositions of the limbs, the platform presents a defined workspace, as well as its end effector. Generally, these mechanisms employ linear actuators connected to both platforms through passive spherical and/or universal joints, providing structural stiffness, stability, and better support [30]. Hence, this mechanism is employed to increase the position accuracy of the caudal fin during high-frequency flapping.

The designed propulsion system produces vectored thrust by biasing the flapping direction from the fin. Using a 3UCU-1S (three universal-cylindrical-universal and one spherical joint) configuration, the moving platform to which the caudal fin is attached is set to only produce oscillatory motion. This means that the parallel robotic system consists of a series of three robotic arms with the same configuration linking both platforms. Each arm is attached to the fixed and moving base by nonactuated universal joints, while motion is produced using linear (prismatic) actuators. These linear actuators vary in lengths $d_i$, thereby exerting force and moment on the moving platform and the caudal fin. Then, by varying the lengths of the limbs, the caudal fin is set to flap, producing thrust and moment. Moreover, the fixed and the moving platforms are linked at their centers by a fourth restrictive limb. Hence, the fin cannot attain translational motion but rotational only, allowing the vehicle to achieve thunniform motion. The technical specifications of the employed linear actuators may be found in Appendix A.

2.3.1. Propeller Workspace and Kinematics

Figure 4 shows a frontal view of the parallel robot, where two actuated arms and the restrictive limb can be seen. To properly comprehend how the limbs are attached to the platforms, passive joints are presented. The angular displacement of the universal joints is defined by ${\delta }_{1,i }$ and ${\delta }_{2,i },$ where $i=1, 2, 3$ represents the number of the limb that each passive joint is attaching to the fixed platform. Figure 4b shows a free body diagram of the parallel mechanism with two added coordinate systems: $\left\{O_F\right\},$ where the origin is placed at the center of the fixed platform, and $\left\{O_M\right\}$, placed at the center of the moving platform. The base platform presents a radius $R_A$ of 0.05 m, while the moving platform is smaller, with a radius $R_B$ of 0.025 m.

The rotational motion of the upper platform is defined by the roll-pitch-yaw Euler angles convention $\left(\alpha ,\beta ,\gamma \right)$, and is measured with respect to frame $\left\{O_F\right\}.$ Flapping is produced by moving two limbs a certain distance, $d_i$, and using the third limb as a pivot. The limbs are connected to the fixed platform at attaching points $A_i,$ while their other ends are connected to the moving platform by joints $B_i.$ The moving platform can be set to oscillate over the $y_M$ axis at a total range for $\beta$ from $\pm 30°$. This is what allows the caudal fin to flap at rates of 0.5 to 5 Hz. Moreover, the parallel mechanism may be controlled to present yaw motion $\gamma$ over the $z_M$ axis, taking values of 0 for sideways flapping or $90°$ for dorsoventral flapping.

For thrust and moment generation, the parallel mechanism is set to follow a predetermined trajectory (orientation) based on its moving platform workspace, i.e., by setting functions to pitch $\left(\beta \right)$ and yaw $\left(\gamma \right)$ angles over time (measured with respect to $\left\{O_F\right\})$. Flapping is obtained by setting the upper platform to oscillate over $y_M$; hence, the angular position of the caudal fin may be defined by a sine function as follows:

$$\beta \left(t\right)=A_M\sin \left(2\pi ft+\phi \right)+b.$$

(10)

Further, the flapping performance of the propeller will depend on the definition of various parameters in Equation (10). $A_M$ will define the flapping amplitude from the tail, and its value will depend on whether or not oscillation is biased. Moreover, $\phi$ and $b$ are the shifting and bias parameters that will set the platform to always start at a zero position and reach a maximum of $\beta =\pm 30°.$ Finally, flapping frequency $f$ determines how fast the propeller oscillates. The thrust and moment exerted by the fin are expected to increase when flapping frequency increases; thus, the proper definition of frequency during biased flapping is of paramount importance in terms of pushing the vehicle forward. Figure 5 shows the two swimming modes. Figure 5a shows sideways flapping, while Figure 5b presents dorsoventral flapping, both of which have a full flapping range of $\beta =\pm 30°$ at $\gamma =0$ and $\gamma =90°,$ respectively.

Table 3. Flapping parameters and expected performance of the caudal fin.

Case	Parameter Values	$\mathbf{Range} \mathbf{for} \mathit{\beta }$	Flapping Direction during Sideways Swimming	Flapping Direction during Dorsoventral Swimming
1.	$A_M=30;$ $\phi =b=0.$	$-30°\le \beta \le 30°$	Port-starboard full range.	Top-bottom full range.
2.	$A_M=b=15;$ $\phi =-0.5\pi .$	$0°\le \beta \le 30°$	Center to starboard range.	Center to bottom range.
3.	$A_M=15;$ $b=-15;$ $\phi =0.5\pi .$	$-30°\le \beta \le 0°$	Center to port range.	Center to top range.

details the expected behavior of the vehicle’s tail based on the defined flapping parameters. Then, the tail can follow either a full-range trajectory, i.e., $\beta =\pm 30°,$ or a biased trajectory, i.e., $\beta =\left[0,30°\right]$ or $\beta =\left[-30°,0\right].$

Since the orientation of the caudal fin is referenced (locally) to the fixed platform coordinate system $\left\{O_F\right\}$, homogeneous transformations are employed to convert such information and reference it to different systems, i.e., $\left\{O_B\right\}$ and $\left\{O_I\right\}$. Hence, to properly define the position and orientation of the several coordinate frames inside different parts of the vehicle, all information is referenced to the inertial earth-fixed frame. By computing a $\left(4\times 4\right)$ homogeneous matrix $T$, it is possible to determine the relationship between two consecutive frames, i.e., $\left\{O_{j-1}\right\}$ to $\left\{O_j\right\}$ by ${}_{j-1}{}^{j}T_{}^{},$ where the inertial frame would be defined by:

$$\left\{O_I\right\}={}_0{}^{0}T_{}^{}=\left[\begin{array}{cccc}{n}_1& o_1& a_1& x\\ n_2& o_2& a_2& y\\ n_3& o_3& a_3& z\\ 0& 0& 0& 1\end{array}\right];$$

(11)

where the $\left[n o a\right]$ submatrix is the identity matrix, and the center of the coordinate system ${\left[x,y,z\right]}^T$ is positioned in an arbitrary location, set as $\left(0,0,0\right).$Then, the relation between the last coordinate frame inside the vehicle $\left\{O_M\right\}$ and $\left\{O_I\right\}$ will be defined by:

$${}_0{}^{M}T_{}^{}={}_0{}^{B}T_{}^{}·{}_B{}^{F}T_{}^{}·{}_F{}^{M}T_{}^{} ,$$

(12)

where ${}_0{}^{B}T_{}^{}$ is the pose of the BAUV’s centroid with respect to an earth-fixed frame, ${}_B{}^{F}T_{}^{}$ is the position and orientation of the fixed platform, and ${}_F{}^{M}T_{}^{}$ is the spatial position and orientation of the propeller.

To build the required transformation matrices ${}_{j-1}{}^{j}T_{}^{}$, the Denavit-Hartenberg $\left(DH\right)$ convention [31] was used throughout this analysis. This formulation was used to translate between the coordinate systems, since they are all part of the vehicle’s body (as found in kinematic chains). This was done by following four basic transformations (two translations and two rotations):

$$DH\left({\mathsf{\Theta }}_i,{\partial }_i,{\epsilon }_i,{\propto }_i\right)=Rot\left(z_{j-1}, {\mathsf{\Theta }}_i\right)·Tran\left(z_{j-1}, {\partial }_i\right)·Tran\left(x_j,{\epsilon }_i\right)·Rot\left(x_j,{\propto }_i\right)=\left[\begin{array}{cccc}c{\mathsf{\Theta }}_i& -s{\mathsf{\Theta }}_{i}c{\propto }_i& s{\mathsf{\Theta }}_{i}s{\propto }_i& {\epsilon }_{i}c{\mathsf{\Theta }}_i\\ s{\mathsf{\Theta }}_i& c{\mathsf{\Theta }}_{i}c{\propto }_i& -s{\propto }_{i}c{\mathsf{\Theta }}_i& {\epsilon }_{i}s{\mathsf{\Theta }}_i\\ 0& s{\propto }_i& c{\propto }_i& {\partial }_i\\ 0& 0& 0& 1\end{array}\right];$$

(13)

where ${\mathsf{\Theta }}_i$ is angular displacement along axis $z_{j-1}$ that makes $x_{j-1}$ and $x_j$ match, ${\partial }_i$ is the distance translated along axis $z_{j-1}$ to match the system origins $\left\{O_j\right\}$ and $\left\{O_{j-1}\right\}$, ${\epsilon }_i$ is the separation distance along axis $x_j$ between both origins, and ${\propto }_i$ is the angular displacement that matches $z_{j-1}$ and $z_j.$ For this analysis, $DH$ units are considered in meters and degrees for linear and angular displacement, respectively. In this way, every coordinate system inside the vehicle was computed, allowing us to determine the actual position and orientation of its main parts, i.e., the propeller’s pose over time.

The parallel mechanism is responsible for moving the caudal fin and propelling the vehicle’s body $\left\{O_B\right\}$. The kinematic analysis of the platforms was rearranged so that systems $\left\{O_F\right\}$ and $\left\{O_M\right\}$ were referenced to the vehicle’s centroid. This means that the position and orientation of both the fixed and moving platforms will be measured based on where the vehicle’s center of buoyancy is placed. From the design of the propulsion system, some considerations must be stated before the coordinate transformations can be computed:

• For the current analysis, it was assumed that matrix ${}_0{}^{B}T_{}^{}$ for system $\left\{O_B\right\}$ was formed from Equations (7)–(9), and therefore, that vector $\mathit{\eta }$ was known. In this way, ${}_0{}^{B}T_{}^{}$ will contain the rotation matrix formed by $\left[\varphi ,\theta ,\psi \right]$and its position information $\left[x,y,z\right]$ with respect to the inertial frame.
• By design, the distance $\mu$ between the vehicle’s center of buoyancy $\left\{O_B\right\}$ and system $\left\{O_F\right\}$ remains constant, i.e., approximately 13 cm, as measured along $x_b.$
• $\mathrm{The} y_F$ axis from the fixed platform’s coordinate system $\left\{O_F\right\}$ was arbitrarily set to always point toward attaching point $A_1$ from the first limb. Also, the three actuated limbs were uniformly distributed along the platform and separated from the origin by 5 cm. The relative position $\left({}_F{}^{A_i}P\right)$ for joints $A_i$, measured from $\left\{O_F\right\}$, was defined by:

$${}_F{}^{A_1}P_{}^{}=\left[\begin{array}{c}{R}_{A}c\left(90°\right)\\ R_{A}s\left(90°\right)\\ 0\end{array}\right], {}_F{}^{A_2}P_{}^{}=\left[\begin{array}{c}{R}_{A}c\left(210°\right)\\ R_{A}s\left(210°\right)\\ 0\end{array}\right], {}_F{}^{A_3}P_{}^{}=\left[\begin{array}{c}{R}_{A}c\left(330°\right)\\ R_{A}s\left(330°\right)\\ 0\end{array}\right] ;$$

(14)

Then, transformation matrix ${}_0{}^{F}T_{}^{}$ was obtained as follows:

$${}_0{}^{F}T_{}^{}={}_0{}^{B}T_{}^{}·{}_B{}^{F}T_{}^{};$$

(15)

where:

$${}_B{}^{F}T_{}^{}=DH\left(90°,0,0,-90°\right)·DH\left(0,\mu ,0,0\right)·DH\left(180°,0,0,0\right)$$

(16)

• The distance between $\left\{O_M\right\}$ and $\left\{O_F\right\}$ will always be restricted by the fourth limb and will stay constant. Moreover, $z_F$ is set to point toward the origin of $\left\{O_M\right\}$. By design, the restrictive limb presents a total height $\ell$ of 19.2 cm, measured along $z_F.$
• The moving platform is in constant motion throughout the navigation task. When the platform oscillates (and therefore, the caudal fin), the spatial positions of attaching joints $B_i$ change iteratively. However, when measured from $\left\{O_M\right\}$, $B_i$ relative positions $\left(PM{B}_i\right)$ will remain constant as follows:

$${}_M{}^{B_1}P_{}^{}=\left[\begin{array}{c}{R}_{B}c\left(90°\right)\\ R_{B}s\left(90°\right)\\ 0\end{array}\right], {}_M{}^{B_2}P_{}^{}=\left[\begin{array}{c}{R}_{B}c\left(198°\right)\\ R_{B}s\left(198°\right)\\ 0\end{array}\right], {}_M{}^{B_3}P_{}^{}=\left[\begin{array}{c}{R}_{B}c\left(342°\right)\\ R_{B}s\left(342°\right)\\ 0\end{array}\right]$$

(17)

When measuring with respect to $\left\{O_F\right\}$, attaching points $B_i$ are defined based on the desired orientation of the moving platform $\left[\alpha ,\beta ,\gamma \right]$, while the rotation matrix $\left({}_F{}^{M}R_{}^{}\right)$ is obtained using Euler angles, as in Equation (8). The orientation and position of the moving platform defined by ${}_0{}^{M}T_{}^{}$ are finally obtained by:

$${}_0{}^{M}T_{}^{}={}_0{}^{F}T_{}^{}·{}_F{}^{M}T_{}^{} ;$$

(18)

where:

$${}_F{}^{M}T_{}^{} =DH\left(0,\ell ,0,0\right)·\left[\begin{array}{cc}{}_F{}^{M}R_{}^{}{}_{\left(3\times 3\right)}& 0_{\left(3\times 1\right)}\\ 0_{\left(1\times 3\right)}& 1\end{array}\right].$$

(19)

Figure 6 shows a schematic of the coordinate systems in space when the vehicle is positioned at an arbitrary location with $\mathit{\eta }={\left[1,0,0,0,0,-90°\right]}^T.$ In Figure 6a, the moving platform roll-pitch-yaw angles are set to $0°$, showing the sideways flapping feature. On the other hand, Figure 6b shows the same schematic but with values for $\left[\alpha ,\beta ,\gamma \right]$ of $\left[0,0,90°\right],$ where dorsoventral swimming mode can be achieved. For brevity purposes, the analysis of inner kinematic chains from limbs inside the parallel mechanism is skipped in this article, but full details are provided in [22].

2.3.2. Thrust and Moment Analysis

Thrust and moment is generated by water displacement, resulting from the motion of the platform and fin. Moreover, the propulsion system applies forces through each linking limb. Each force is the result of the control action over each limb to make the propeller flap at a certain amplitude, bias, and frequency. This means that the three linear actuators are constantly exerting forces, pushing the vehicle forward and laterally.

To compute the forces of the limbs that make the vehicle’s propeller move as desired, a dynamic equation for the designed parallel mechanism was required. The dynamic modeling of such systems is not a straightforward process; its complexity mainly depends on the configuration of the inner closed-loop structures. Based on the distribution of passive and active joints displayed in Figure 4, and by applying reduced inverse dynamic modeling using Lagrange formalism, force ${\mathcal{F}}_i$ for each linear actuator was obtained by solving a system of equations:

$$\frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial \dot{q_j}}\right)-\left(\frac{\partial \mathcal{L}}{\partial q_j}\right)+{\displaystyle \sum }_{k=1}^6{\lambda }_k \frac{\partial f_k}{\partial q_j}={\mathbb{Q}}_j; for j=1,..,9,$$

(20)

where:

$q_j=\left\{\begin{array}{cc}{d}_i& j=i=1,2,3\\ {\delta }_{1,i}& j=4,5,6;i=1,2,3\\ {\delta }_{2,i}& j=7,8,9;i=1,2,3\end{array}\right.$ and ${\mathbb{Q}}_j=\left\{\begin{array}{cc}{\mathcal{F}}_i& j=i=1,2,3\\ 0& j=4,\dots ,9\end{array}\right.$.

From Equation (20), actuating forces ${\mathcal{F}}_1, {\mathcal{F}}_2, {\mathcal{F}}_3$ applied by the limbs were found based on the derivation of the Lagrange equation of the system, where the kinetic and potential energies of the moving platform are computed in [22]. Also, Equation (20) defines a system of nine equations in which normalized limb lengths $d_i$ are the independent generalized coordinates and revolute ${\delta }_{1,i}$ and ${\delta }_{2,i}$ from universal passive joints are dependent generalized coordinates. The six Lagrange multipliers ${\lambda }_k$ were obtained by iteratively solving six simultaneous equations when $j=4,\dots ,9$ using kinematics constraint equations $f_k$ (where $k=1,\dots ,6$). To solve the control problem of the parallel mechanism, it was assumed that platform’s orientation, velocity, and acceleration were known. Then, based on the desired flapping performance on the propeller, the forces required by each limb were determined.

Since each limb exerts acting forces on the propeller to produce motion, it is necessary to comprehend how such forces add to the generation of thrust and moment measured by the vehicle’s coordinate system $\left\{O_B\right\}.$ First, it was assumed that the force acting on each limb was vectorized along its axis of action (Figure 4b). Hence, the direction of the force applied per limb could be determined by computing a unit vector ${\widehat{n}}_{AB}$ that went from $A_i$ to $B_i,$ or $\stackrel{\rightharpoonup }{A_i{B}_i},$ as follows:

$$\begin{array}{c}\stackrel{\rightharpoonup }{A_i{B}_i}=\left(B_i{}_x-A_i{}_x\right)i+\left(B_i{}_y-A_i{}_y\right)j+\left(B_i{}_z-A_i{}_z\right)k \therefore \\ {\widehat{n}}_{AB}=\frac{\overline{A_i{B}_i }}{\left|\overline{A_i{B}_i }\right|}. \end{array}$$

(21)

Moreover, the spatial positions of attachment points $A_i’s$ and $B_i’s$ in Equation (21) should be iteratively computed and referenced to $\left\{O_B\right\}$. Also, $\left\{O_B\right\}$ will be considered as the inertial coordinate system. This will define the direction of the forces exerted by each limb while the platform is flapping, and their effects on the vehicle’s centroid. Then, to compute the attaching positions referenced to $\left\{O_B\right\}$, transformation matrices may be employed:

$${\left[\begin{array}{c}{A}_i\\ 1\end{array}\right]}_{x_b{y}_b{z}_b}={}_B{}^{F}T_{}^{}·\left[\begin{array}{c}{}_F{}^{A_i}P_{}^{}\\ 1\end{array}\right];$$

(22)

and

$${\left[\begin{array}{c}{B}_i\\ 1\end{array}\right]}_{x_b{y}_b{z}_b}={}_B{}^{F}T_{}^{}·{}_F{}^{M}T_{}^{}·\left[\begin{array}{c}{}_M{}^{B_i}P_{}^{}\\ 1\end{array}\right].$$

(23)

The vectorized force ${\overrightarrow{\mathcal{F}}}_i$ applied per limb is then computed:

$$\begin{array}{l}{\overrightarrow{\mathcal{F}}}_i={\mathcal{F}}_i·{\widehat{n}}_{AB}\in {ℝ}^3\therefore \\ {\overrightarrow{\mathcal{F}}}_i=\left[{\mathcal{F}}_i{}_x,{\mathcal{F}}_i{}_y,{\mathcal{F}}_i{}_z\right].\end{array}$$

(24)

Then, total force ${\mathit{\tau }}_1$ produced by the propulsion system to the BAUV will be defined by:

$${\mathit{\tau }}_1={\left[F_X,F_Y,F_Z\right]}^T={\left[{\displaystyle \sum }_{i=1}^3{F}_i{}_x, {\displaystyle \sum }_{i=1}^3{F}_i{}_y,{\displaystyle \sum }_{i=1}^3{F}_i{}_z \right]}^T.$$

(25)

For the moment generated by the designed propeller, the vectorized forces per limb as determined using Equation (24) were required. Also, the moment was measured around the vehicle’s centroid, which, by design, was aligned to the center of the moving platform $\left\{O_M\right\}$. Hence, the torque produced could be determined dynamically by computing the radius of action ${\overrightarrow{RB}}_i$ from each universal joint $B_i$ to the center of the upper platform and the vectorized force ${\overrightarrow{\mathcal{F}}}_i$ applied to that union. ${\overrightarrow{RB}}_i$ was obtained using the first three elements from the ${\left[\begin{array}{c}{B}_i\\ 1\end{array}\right]}_{x_b{y}_b{z}_b}$ matrix, computed in Equation (23). The vectorized moment ${\overrightarrow{\mathcal{M}}}_i$ produced per limb was the cross product between the radius of action and the applied vectorized force as:

$${\overrightarrow{\mathcal{M}}}_i={\overrightarrow{RB}}_i\times {\overrightarrow{\mathcal{F}}}_i\in {ℝ}^3$$

(26)

Finally, the net torque and moment ${\mathit{\tau }}_2$ applied to the BAUV by the propulsion system was defined by:

$${\mathit{\tau }}_2={\left[K_X,M_Y,N_Z\right]}^T={\left[{\displaystyle \sum }_{i=1}^3{\mathcal{M}}_i{}_x, {\displaystyle \sum }_{i=1}^3{\mathcal{M}}_i{}_y,{\displaystyle \sum }_{i=1}^3{\mathcal{M}}_i{}_z \right]}^T.$$

(27)

2.4. Parallel Mechanism Controller

The flapping performance of the caudal fin depends on the linear displacement of the actuators inside the parallel mechanism. Then, once the upper platform is set to follow a desired trajectory, the limbs should reach certain distances $d_i$ to attain the desired roll-pitch-yaw $\left(\alpha ,\beta ,\gamma \right)$ rotational motion. To properly regulate the performance of the limbs, a feedforward plus feedback PD controller was employed. The feedforward term from the control scheme used the dynamics model of the parallel mechanism to proactively compute the forces required to drive the caudal fin at certain speeds and positions. In compliment, the feedback PD controller computed the error between the desired limb distance $d_i{}_D$ and the real limbs positions $d_i$. The mathematical analysis of the dynamics of the parallel mechanism used inside the BAUV propulsion system, and the implementation of the feedforward plus feedback controller are duly explained in [22]. Figure 7 presents a block diagram of the closed-loop control scheme incorporated to regulate limb positions based on the desired platform performance over time.

In the control scheme in Figure 7, the dispositions of the limbs depend on the swimming mode, defined by the value assigned to $\gamma .$ Moreover, the positions of the limbs are iteratively changing based on the preestablished trajectory of $\beta .$ Roll motion $\alpha$ in the platform is always set to zero. Then, the control law of the PD controller based on the computation of error is defined by:

$$f_i\left(t\right)=kp·e_i\left(t\right)+kv·{\dot{e}}_i\left(t\right), i=1,2,3;$$

(28)

where:

$$\begin{array}{l}{e}_i\left(t\right)=d_i{}_D\left(t\right)-d_i\left(t\right);\\ {\dot{e}}_i\left(t\right)={\dot{d}}_i{}_D\left(t\right)-{\dot{d}}_i\left(t\right). \end{array}$$

(29)

The proportional and derivative gains $\left(kp,kv\right)$ from the feedback controller were set to 2500 and 100 respectively. The relation between gains was established to achieve a critically damped response. Elements $M_P\left(\dot{d}\right), C_P\left(d,\delta ,\dot{\delta }\right), g_P\left(\delta \right)$ and ${\lambda }_k$ inside the block diagram are the mass, Coriolis, gravitational forces, and Lagrange constraint matrices respectively. Such matrices define the dynamics model from the parallel mechanism. Hence, the feedforward component complements the control action taken by the PD by anticipating the dynamic behavior of the parallel robotic system. Then, the forces $F_i$ required to actuate each limb are defined by the actions computed from both strategies as:

$$F_i\left(t\right)=M_P\left(\dot{d}\right)·f_i^{\prime }\left(t\right)+H, i=1,2,3;$$

(30)

where:

$$\begin{array}{c}{f}_i^{\prime }\left(t\right)={\ddot{d}}_i{}_D+f_i\left(t\right);\\ H=C_P\left(d,\delta ,\dot{\delta }\right)+g_P\left(\delta \right)+{\lambda }_k.\end{array}$$

(31)

The implementation of this control strategy makes it possible to adequately correct the limbs in the event of sudden changes. The flapping frequency and bias of the caudal fin should be updated according to the thrust and moment required to steer the BAUV. Then, by implementing a feedforward plus feedback controller into the designed propulsion system, a smooth response with low error rates may be attained [22].

2.5. Vehicle Dynamics

The dynamic modeling for biomimetic vehicles followed the same criteria as the methodologies employed to model conventional underwater vehicles [12]. Nonetheless, the biomimetic aspects of the prototype had to be considered to adequately estimate some hydrodynamic terms. As in the previous section, where the thrust and moment produced by the biomimetic propeller required its own analysis, some dynamic terms such as drag, added mass, lifting, and/or crossflow coefficients will depend on the shape of the vehicle and will require their own analysis.

2.5.1. Assumptions

Before computing the vehicle’s dynamics model, some assumptions were made for simplification purposes:

• The study was done using a 6 DOF BAUV, capable of achieving surge-heave-sway translational motions and roll-pitch-yaw rotational motions.
• The BAUV was assumed to be a rigid body of constant mass and shape throughout the mission.
• The vehicle was not considered to be neutrally buoyant. The origin of the body-fixed reference frame was assigned to the vehicle’s center of buoyancy. The center of gravity and buoyancy were considered to be aligned along the $z_b$ axis, separated by a metacentric height of $z_g$ m.
• The vehicle was assumed to be symmetric over the $\left(x,y\right)$ and $\left(x,z\right)$ planes.
• No environmental disturbances such as wind, waves, current, vehicle wave-making loads, and other external forces apart from those produced by the propeller were taken into account in this study.
• The vehicle’s source of thrust and moment was only obtained by the motion of the vehicle’s caudal fin and its orientation. These forces were measured with respect to the body-fixed frame. The pectoral fins at the midbody were considered to be rigid. In other words, their purpose was for stability, and they were considered in this light in the hydrodynamics analysis.

2.5.2. Rigid-Body Dynamics

The computation of the dynamics of the designed BAUV was based on the standard six DOF underwater vehicle model, considering the aforementioned assumptions. Using the Newton-Euler methodology [28], the dynamic equation for the BAUV was defined by:

$$\mathit{M}\dot{\mathit{v}}+\mathit{C}\left(\mathit{v}\right)\mathit{v}+\mathit{D}\left(\mathit{v}\right)\mathit{v}+\mathit{g}\left(\mathit{\eta }\right)=\mathit{\tau }$$

(32)

where $\mathit{M}\in {ℝ}^{6\times 6}$ is the inertial and added mass matrix, $\mathit{C}\left(\mathit{v}\right)\in {ℝ}^{6\times 6}$ is the Coriolis, centripetal forces and added mass terms matrix, $\mathit{D}\left(\mathit{v}\right)\in {ℝ}^{6\times 6}$ is the vehicle’s damping matrix with drag terms, $\mathit{g}\left(\mathit{\eta }\right)\in {ℝ}^{6\times 1}$ is the vector of hydrostatic forces, $\mathit{\tau }\in {ℝ}^{6\times 1}$ is the input vector of forces and moments produced by the propeller, and $v\in {ℝ}^{6\times 1}$ are the linear and angular velocities measured from the body-fixed frame $\left\{O_B\right\}.$

The equations of motion from the vehicle were defined based on the body-fixed frame origin located at the BAUV’s center of buoyancy. For the present study, the BAUV was considered as a rigid body with constant mass $m$ and a diagonal inertia tensor matrix $I=diag\left[I_{xx},I_{yy},I_{zz}\right]$ about $\left(x_b,y_b,z_b\right).$ Then, the rigid body set of equations for each degree of freedom (surge, sway, heave, roll, pitch, yaw) from the underwater vehicle was defined as follows [28]:

$$\begin{array}{c}X=m\left[\dot{u}-vr+wq-x_g\left(q^2+r^2\right)+y_g\left(pq-\dot{r}\right)+z_g\left(pr+\dot{q}\right)\right];\\ Y=m\left[\dot{v}-wp+ur-y_g\left(p^2+r^2\right)+z_g\left(qr-\dot{p}\right)+x_g\left(qp+\dot{r}\right)\right];\\ Z=m\left[\dot{w}-uq+vp-z_g\left(p^2+q^2\right)+x_g\left(rp-\dot{q}\right)+y_g\left(rq+\dot{p}\right)\right];\\ K=I_{xx}\dot{p}+\left(I_{zz}-I_{yy}\right)qr+m\left[y_g\left(\dot{w}-uq+vp\right)-z_g\left(\dot{v}-wp+ur\right)\right];\\ M=I_{yy}\dot{q}+\left(I_{xx}-I_{zz}\right)pr+m\left[z_g\left(\dot{u}-vr+wq\right)-x_g\left(\dot{w}-uq+vp\right)\right];\\ N=I_{zz}\dot{r}+\left(I_{yy}-I_{xx}\right)pq+m\left[x_g\left(\dot{v}-wp+ur\right)-y_g\left(\dot{u}-vr+wq\right)\right].\end{array}$$

(33)

To achieve static stability while swimming, the centers of buoyancy and gravity from the vehicle were aligned along the $x_b$ and $y_b$ axes to reduce roll and pitch rotations [10]. Then, the center of gravity of the vehicle was set as low as possible with a metacentric height $z_g.$

2.5.3. Hydrodynamic Added Mass Terms

To complete the terms from the $\mathit{M}$ and $\mathit{C}\left(\mathit{v}\right)$ matrices, added mass coefficients from the vehicle were computed. The added mass terms considered the volume of surrounding fluid that the vehicle displaces while accelerating, and their analytical estimation was done based on the hull shape of the proposed design [32]. To compute hydrodynamic added mass terms, a methodology like the one developed by Prestero [33] was employed. Usually, underwater vehicle hull designs are torpedo-shaped, even when some biomimetic aspects are introduced. Then, the computation of added mass coefficients is done by approximating such shapes to geometric ones, i.e., a cylinder or a prolate ellipsoid [14,34]. For the current analysis, shape Equations (1)–(3) were employed for the estimation of added mass coefficients. Furthermore, a top-bottom and port-starboard hull symmetry was assumed. Forces and moments due to the effect of added mass were defined by the following equations [33]:

$$\begin{array}{l}{X}_A= X_{\dot{u}}\dot{u}+Z_{\dot{w}}wq+Z_{\dot{q}}{q}^2-Y_{\dot{v}}vr-Y_{\dot{r}}{r}^2;\\ Y_A=Y_{\dot{v}}\dot{v}+Y_{\dot{r}}\dot{r}+X_{\dot{u}}ur-Z_{\dot{w}}wp-Z_{\dot{q}}pq;\\ Z_A= Z_{\dot{w}}\dot{w}+Z_{\dot{q}}\dot{q}-X_{\dot{u}}uq+Y_{\dot{v}}vp-Y_{\dot{r}}rp;\\ K_A= K_{\dot{p}}\dot{p};\\ M_A= M_{\dot{w}}\dot{w}+M_{\dot{q}}\dot{q}-(Z_{\dot{w}}-X_{\dot{u}})uw-Y_{\dot{r}}vp+\left(K_{\dot{p}}-N_{\dot{r}}\right)rp-Z_{\dot{q}}uq;\\ N_A= N_{\dot{v}}\dot{v}+N_{\dot{r}}\dot{r}-(X_{\dot{u}}-Y_{\dot{v}})uv+Z_{\dot{q}}wp-\left(K_{\dot{p}}-M_{\dot{q}}\right)pq+Y_{\dot{r}}ur.\end{array}$$

(34)

From Equation (34), axial added mass $X_{\dot{u}}$ and rolling added mass $K_{\dot{p}}$ were estimated from the shape relationship defined by Blevins [35], where the vehicle’s total axial length $l$ and diameter $d$ had to be considered. Then, these coefficients were obtained by the following computations:

$$\begin{array}{l}{X}_{\dot{u}}=-\frac{4{\beta }_B\pi }{3}{\left(\frac{d}{2}\right)}^3,\\ K_{\dot{p}}=-{{\displaystyle \int }}_{xpf}^{xp}\left(0.5\pi ·\rho ·h{f}^4\right)dx.\end{array}$$

(35)

The Blevins’s empirical parameter ${\beta }_B$ obtained from the relationship $\left(l/d\right)$ for the BAUV was estimated to be $0.2780.$ For the estimation of $K_{\dot{p}},$ $hf$ is the fin height over the centerline, considered to be of 13.90 cm, $\rho$ is the density of seawater, i.e., $1.03\times {10}^{-3} \mathrm{kg}/{\mathrm{m}}^3$, and the range of integration limits from the aft end of the peduncle section $xpf$ to its forward end $xp$ was measured with respect to $\left\{O_B\right\}.$

To estimate the crossflow added mass terms from Equation (34), strip theory was used, and the vehicle was considered to comprise a slender body over its centerline [36]. This means that the vehicle’s added mass terms were estimated based on the integration of the two-dimensional (2D) cross section added mass coefficients along its body. The geometrical shape of the BAUV hull was approximated by cylindrical slices along its centerline, and the 2D added mass $a_m$was computed by using the formulation proposed by Newman [37]:

$$a_m\left(x\right)=\pi ·\rho ·r{\left(x\right)}^2;$$

(36)

where $\rho$ is the surrounding fluid density, and $r\left(x\right)$ is the hull radius as a function of axial position, as defined in Equations (1)–(3). Moreover, in the sections where the fins are placed, i.e., midbody and peduncle, the 2D added mass $a_{mf}$ was computed by approximating the vehicle’s cross section to a circle with fins using the following equation proposed by Newman:

$$a_{mf}\left(x\right)=\pi \rho \left[r{\left(x\right)}^2+\frac{{\left(r{\left(x\right)}^2-h{f}^2\right)}^2}{h{f}^2}\right];$$

(37)

where $hf$ is the maximum distance reached from the fins, measured with respect to the centerline. Then, the crossflow terms for the BAUV were defined by the formulations presented in the SNAME Principles of Naval Architecture [38]:

$$\begin{array}{l}{Y}_{\dot{v}}=-{{\displaystyle \int }}_L{A}_m\left(x\right) dx;\\ M_{\dot{w}}=-{{\displaystyle \int }}_{L}x{A}_m\left(x\right) dx;\\ N_{\dot{r}}=-{{\displaystyle \int }}_L{x}^2{A}_m\left(x\right) dx;\end{array}$$

(38)

where $Z_{\dot{w}}=Y_{\dot{v}}, N_{\dot{v}}=-M_{\dot{w}}, Y_{\dot{r}}=N_{\dot{v}}, Z_{\dot{q}}=M_{\dot{w}}, M_{\dot{q}}=N_{\dot{r}}.$From Equation (38), $L$ means that the integration is done over the vehicle’s total length. Furthermore, $A_m\left(x\right)$ is the two-dimensional added mass substituted by Equations (36) and (37) based on the presence of fins. Then, matrix $\mathit{M}$ may be complemented by the added mass terms computed using Equation (38):

$$\mathit{M}=\left[\begin{array}{cccccc}m-X_{\dot{u}}& 0& 0& 0& z_{g}m& -y_{g}m\\ 0& m-Y_{\dot{v}}& 0& -z_{g}m& 0& x_{g}m-Y_{\dot{r}}\\ 0& 0& m-Z_{\dot{w}}& y_{g}m& -x_{g}m-Z_{\dot{q}}& 0\\ 0& -z_{g}m& y_{g}m& I_{xx}-K_{\dot{p}}& 0& 0\\ z_{g}m& 0& -x_{g}m-M_{\dot{w}}& 0& I_{yy}-M_{\dot{q}} & 0\\ -y_{g}m& x_{g}m-N_{\dot{v}}& 0& 0& 0& I_{zz}-N_{\dot{r}}\end{array}\right].$$

(39)

The Coriolis and centripetal forces matrix $C\left(v\right)$ may be completed with the remaining cross-terms that result from added mass coupling, presenting its final form as:

$$\mathit{C}\left(\mathit{v}\right)=\left[\begin{array}{cccccc}0& -\left(m+X_{vr}\right)r& \left(m-X_{wq}\right)q& (y_{g}q+z_{g}r)m& -\left(x_{g}m+X_{qq}\right)q& -\left(x_{g}m+X_{rr}\right)r\\ \left(m-Y_{ur}\right)r& -Y_{uv}u& -\left(m+Y_{wp}\right)p& (x_{g}m-Y_{pq})q& z_{g}mr& 0\\ -\left(m+Z_{uq}\right)q& \left(m-Z_{vp}\right)p& -Z_{uw}u& -z_{g}mp& -z_{g}mq& \left(x_{g}m-Z_{rp}\right)p\\ -\left(z_{g}r+y_{g}q\right)m& y_{g}mp& z_{g}mp& 0& 0& \left(I_{zz}-I_{yy}\right)q\\ \left(x_{g}m-M_{uq}\right)q& -\left(x_{g}m+M_{vp}\right)p& z_{g}mq-M_{uw}u& \left(\left(I_{xx}-I_{zz}\right)-M_{rp}\right)r& 0& -z_{g}mv\\ -\left(x_{g}m+N_{ur}\right)r& y_{g}mr-N_{uv}u& \left(x_{g}m-N_{wp}\right)p& 0& \left(\left(I_{yy}-I_{xx}\right)-N_{pq}\right)p& 0\end{array}\right];$$

(40)

where most of the Coriolis matrix cross terms were evaluated from the previously derived added mass coefficients, as follows [33]:

$$\begin{array}{l}{X}_{vr}=-Y_{\dot{v}}, X_{wq}=Z_{\dot{w}, } X_{qq}=Z_{\dot{q}}, X_{rr}=-Y_{\dot{r}};\\ Y_{ur}=X_{\dot{u}}+Y_{urf}, Y_{uv}=Y_{uvl}+Y_{uvf}, Y_{wp}=-Z_{\dot{w}}, Y_{pq}=-Z_{\dot{q}};\\ Z_{uq}=-X_{\dot{u}}+Z_{uqf}, Z_{vp}=Y_{\dot{v}}, Z_{uw}=Z_{uwl}+Z_{uwf}, Z_{rp}=Y_{\dot{r}};\\ M_{uq}=-Z_{\dot{q}}+M_{uqf}, M_{vp}=-Y_{\dot{r}}, M_{uw}=X_{\dot{u}}-Z_{\dot{w}}+M_{uwl}+M_{uwf}; \\ M_{rp}=K_{\dot{p}}-N_{\dot{r}}, N_{ur}=Y_{\dot{r}}+N_{urf}, N_{uv}=Y_{\dot{v}}-X_{\dot{u}}+N_{uvl}+N_{uvf};\\ N_{wp}=Z_{\dot{q}}, N_{pq}=M_{\dot{q}}-K_{\dot{p}}.\end{array}$$

(41)

From Equation (41), some added mass cross-terms had to be completed by considering vehicle’s body lift forces and moments, i.e., $Y_{ur}, Y_{uv}, Z_{uq}, Z_{uw}, M_{uq}, M_{uw}, N_{ur},$ and $N_{uv}.$ The BAUV presents pectoral fins on its middle body, and the lifting forces and moments that they produce were also considered. Body lift force and moment coefficients $Y_{uvl}, Z_{uwl}, M_{uwl}$ and $N_{uvl}$ were estimated by Hoerner’s formulations [39]:

$$\begin{array}{l}{Y}_{uvl}=Z_{uwl}=-\frac{1}{2}\rho d^2{c}_{\beta },\\ M_{uwl}=-N_{uvl}=-\frac{1}{2}\rho d^2{c}_{\beta }{x}_{pr};\end{array}$$

(42)

where $\rho$ is the density of surrounding fluid, $d$ is the vehicle’s hull diameter, $c_{\beta }$ is the Hoerner’s lift coefficient (estimated to be of $0.9692),$ and $x_{pr}=-0.1770 \mathrm{m}$ is the pressure point where viscous forces center on the vehicle, located at approximately 65% of vehicle’s total length from its nose tip and measured with respect to $\left\{O_B\right\}.$Finally, the fin lift force and moment coefficients were estimated based on the characteristics of the pectoral fins, as follows:

$$\begin{array}{l}{Y}_{uvf}=Z_{uwf}=-\frac{1}{2}\rho c_L{\delta }_{fin},\\ Y_{urf}=-Z_{uqf}=-\frac{1}{2}\rho c_L{\delta }_{fin}{x}_{fin},\\ M_{uwf}=-N_{uvf}=\frac{1}{2}\rho c_L{\delta }_{fin}{x}_{fin},\\ M_{uqf}=N_{urf}=-\frac{1}{2}\rho c_L{\delta }_{fin}{x}_{fin}^2;\end{array}$$

(43)

where ${\delta }_{fin}=9.0322\times {10}^{-3} {\mathrm{m}}^2$ is the fin’s planform area, $x_{fin}=0.16 \mathrm{m}$ is the fin’s axial position with respect to $\left\{O_B\right\}$, and $c_L$ is the rate of change of the lift coefficient, estimated with respect to the fin’s aspect ratio $AR$as:

$$c_L={\left[\frac{1}{1.8\pi }+\frac{1}{\pi \left(AR\right)}\right]}^{-1}.$$

(44)

From Equation (44), the effective fin aspect ratio $\left(AR\right)$ was estimated by the relationship between the fin’s planform area ${\delta }_{fin}$ and span $S_{fin}:$

$$AR=2\left(\frac{S_{fin}^2}{{\delta }_{fin}}\right).$$

(45)

2.5.4. Hydrodynamic Damping Coefficients

To continue with the hydrodynamics analysis of the vehicle, drag forces have to be computed. These forces are considered to act in the opposite direction of the BAUV’s motion due to the viscosity of water [40]. As such, vehicle’s drag will directly depend on the surrounding fluid density $\rho ,$the vehicle’s projected area $A,$ drag coefficient $C_d,$ and body velocity $V$, defined by Morison equation as:

$$f\left(V\right)=\frac{1}{2}\rho C_{d}A\left|V\right|V.$$

(46)

Equation (46) is used to define the drag term due to velocity in direction $V$. For axial drag $X_{u\left|u\right|}$, coefficient $C_d\cong 0.0582$ was estimated based on Triantafyllou empirical formulations [41], using the vehicle’s ratio between its diameter and length and by setting the frontal area of the vehicle $A_f$ to $1.68\times {10}^{-2} {\mathrm{m}}^2.$ The damping matrix $\mathit{D}\left(\mathit{v}\right)$ is formed by drag coefficients on each of the vehicle’s DOF, $X_{u\left|u\right|}, Y_{v\left|v\right|}, Z_{w\left|w\right|}, K_{p\left|p\right|}, M_{q\left|q\right|}$ and $N_{r\left|r\right|}.$ Furthermore, drag components were added considering pectoral fin drag forces. Then, the final form of the built damping matrix will be:

$$\mathit{D}\left(\mathit{v}\right)=\left[\begin{array}{cccccc}{X}_{u\left|u\right|}\left|u\right|& 0& 0& 0& 0& 0\\ 0& Y_{v\left|v\right|}\left|v\right|& 0& 0& 0& 0\\ 0& 0& Z_{w\left|w\right|}\left|w\right|& 0& 0& 0\\ 0& 0& 0& K_{p\left|p\right|}\left|p\right|& 0& 0\\ 0& 0& 0& 0& M_{q\left|q\right|}\left|q\right|& 0\\ 0& 0& 0& 0& 0& N_{r\left|r\right|}\left|r\right|\end{array}\right].$$

(47)

To solve the drag terms from Equation (47), strip theory was used as for the added mass terms estimated in the previous section. Then, the vehicle’s hull drag was approximated by adding the drag estimations of its two-dimensional cross-sections. This methodology may yield some inaccurate estimations for drag coefficients, and some corrections may be required through experimental data, as explained in [33]. However, for the present report, the methodology was used to complete the equations of vehicle motion and to obtain proper simulations of the BAUV’s performance. Then, by using the strip theory over the BAUV’s centerline $x$, and by considering the features of the pectoral fins, drag coefficients could be approximated by the following equations:

$$\begin{array}{l}{Y}_{v\left|v\right|}=Z_{w\left|w\right|}=-\frac{1}{2}\rho c_{dc}{{\displaystyle \int }}_{L}2r\left(x\right) dx-\frac{1}{2}\rho c_{df}\left(2{\delta }_{fin}\right),\\ M_{q\left|q\right|}=N_{r\left|r\right|}=-\frac{1}{2}\rho c_{dc}{{\displaystyle \int }}_{L}2x^{3}r\left(x\right) dx-\frac{1}{2}\rho c_{df}\left(2x_{fin}^3{\delta }_{fin}\right).\end{array}$$

(48)

On the right-hand side of Equation (48), drag effects from the vehicle’s hull shape and pectoral fins are added and defined as negative, since they oppose the vehicle’s motion. Moreover, the hull drag coefficient $c_{dc}$ was set to 1.1 based on the Hoerner drag estimation for a cylinder, while that of the pectoral fins, $c_{df},$ was set to 0.68 based on their taper ratio between height and span. Furthermore, ${\delta }_{fin}$ and $x_{fin}$ are, respectively, the fin’s planform area and position measured along the centerline with respect to $\left\{O_B\right\}$.

The final values of estimated hydrodynamic added mass and damping coefficients are provided in Appendix B.

2.5.5. Hydrostatic Forces and Moments

Finally, matrix $\mathit{g}\left(\mathit{\eta }\right)$ takes into consideration the restoring forces and moments acting on the vehicle. Then, based on the vehicle’s attitude, its weight ($W)$ and buoyancy $\left(B\right)$ forces acting through centers of gravity $\left(x_g,y_g,z_g\right)$ and buoyancy $\left(x_b,y_b,z_b\right)$ (defined on Table 1), the $\mathit{g}\left(\mathit{\eta }\right)$ matrix is defined as follows:

$$\mathit{g}\left(\mathit{\eta }\right)=\left[\begin{array}{c}\left(W-B\right)sin\left(\theta \right)\\ -\left(W-B\right)cos\left(\theta \right)sin\left(\varphi \right)\\ -\left(W-B\right)cos\left(\theta \right)cos\left(\varphi \right)\\ \left(z_{g}W-z_{b}B\right)\cos \left(\theta \right)\sin \left(\varphi \right)+\left(y_{g}W-y_{b}B\right)cos\left(\theta \right)cos\left(\varphi \right)\\ \left(z_{g}W-z_{b}B\right)\sin \left(\theta \right)+\left(x_{g}W-x_{b}B\right)\cos \left(\theta \right)cos\left(\varphi \right)\\ \left(x_{g}W-x_{b}B\right)\cos \left(\theta \right)\sin \left(\varphi \right)+\left(y_{g}W-y_{b}B\right)\sin \left(\theta \right)\end{array}\right].$$

(49)

From Equation (49), the orientation of the vehicle is referenced to inertial frame $\left\{O_I\right\},$ while the positions of the centers of gravity and buoyancy are referenced to body-fixed coordinates $\left\{O_B\right\}.$ $\left(\theta ,\varphi ,\psi \right)$ angular positions are iteratively computed using Euler angles convention, as defined on the kinematics analysis.

2.6. Open-Loop Trajectories and Waypoint Guidance System

Open-loop simulations were implemented to validate the vehicle’s dynamic model. Hence, by doing some algebraic manipulation of Equation (32), the plant equation was defined as:

$$\dot{v}={\mathit{M}}^{-1}\left[\mathit{\tau }-\mathit{C}\left(\mathit{v}\right)\mathit{v}-\mathit{D}\left(\mathit{v}\right)\mathit{v}-\mathit{g}\left(\mathit{\eta }\right)\right].$$

(50)

Then, the vehicle’s behavior will depend not only on the modeled hydrodynamic effects, but also on forces vector $\mathit{\tau }$ entering the plant and produced by the propulsion system. From plant’s Equation (50), linear and angular accelerations of the BAUV were obtained as the system output. The vehicle’s pose with respect to earth-fixed frame was then computed by integrating $\dot{v}$ and using the coordinate transformations from Equation (7). Furthermore, simulations were conducted to determine the required flapping frequency that produced forward thrust on the vehicle considering full range and biased oscillations.

2.6.1. Waypoint Guidance System

To navigate the designed BAUV, a waypoint guidance system was implemented. The vehicle was guided toward specific coordinates on the horizontal $\left(x,y\right)$ plane via sideways flapping. The strategy consisted of correcting the vehicle’s orientation based on its heading angle toward a specific target, i.e., by updating the flapping function $\beta$ of the fin. Then, heading angle correction was achieved by biasing the flapping performance from the vehicle’s caudal fin.

Figure 8 shows a schematic of the coordinates between the position of the vehicle and an arbitrary waypoint, all measured with respect to $\left\{0_I\right\}$. The desired heading yaw angle ${\psi }_d$ was computed based on the trigonometric relationship between the desired target position or waypoint $\left(x_d,y_d\right)$ and the vehicle’s real position $\left(x,y\right)$. Then, the error signal was determined as the difference between the desired ${\psi }_d$ and real yaw angle $\psi .$

Based on the information displayed in Figure 8, the desired heading angle ${\psi }_d$ can be obtained by:

$${\psi }_d=atan2\left(\left(y_d-y\right),\left(x_d-x\right)\right) \left(\mathrm{rad}\right);$$

(51)

where $atan2\left(∎Y,∎X\right)$ is the two-argument arctangent and $2\mathsf{\pi }$ variant of ${\tan }^{-1}\left(\frac{∎Y}{∎X}\right).$ The desired heading angle will range between $-\mathsf{\pi }\le {\psi }_d\le \mathsf{\pi }$ radians, as well as the vehicle’s real orientation $\psi .$

Since vehicle’s spatial position and orientation change over time, ${\psi }_d$ should be iteratively computed, as well as the error signal, where:

$$e={\psi }_d-\psi \left(\mathrm{rad}\right).$$

(52)

Furthermore, the Euclidean distance between a waypoint and the vehicle’s position is computed as a convergence criterion. This is helpful when the vehicle should reach several positions in a given trajectory. By defining a waypoint area, the guidance system can determine how far the vehicle is from the current goal or if it has reached it. The Euclidean distance $E_d$ is computed by:

$$E_d=\sqrt{{\left(x_d-x\right)}^2+{\left(y_d-y\right)}^2 } \left(\mathrm{m}\right).$$

(53)

If the vehicle has entered the waypoint region $(E_d\le r)$, then the convergence criterion has been met and the BAUV should head toward a new target. If all targets have already been reached, the mission has been completed.

Figure 9 shows a closed-loop block diagram of the waypoint guidance system. First, the coordinates of the desired waypoints are introduced into the system. The vehicle’s real attitude is iteratively compared to the desired heading orientation, and the error signal is processed. The flapping parameters of the tail are updated to produce the desired performance over $\beta$. Inside the propeller’s dynamics and control block, the feedforward plus feedback controller regulates the limb performance inside the parallel mechanism to achieve the desired flapping behavior. Then, forces produced by the fin’s motion input the BAUV’s dynamic model block and the vehicle is propelled. Finally, a change of coordinates is computed.

Finally, corrective actions were determined based on the error in Equation (52); these consisted of the adequate parameter selection for the caudal fin’s flapping performance $\beta$ (Table 3). Corrections were developed as if–then rules based on the computed value for $e$, where turning was achieved by biasing the motion from vehicle’s tail. For instance, if the vehicle has to correct its orientation, i.e., to the left or right, the control action will bias flapping toward either the port or starboard, respectively. The guidance system was set to wait for a period of 2 s to update new parameters for motion $\beta .$ This was strategically determined, since the minimum increase for flapping frequency is 0.5 Hz, allowing the platform to reach a neutral position before receiving new instructions. A parameter-selection flowchart is presented in Figure 10.

2.6.2. Simulations

The BAUV’s swimming performance was simulated using Simulink from MATLAB^®. During the simulations, the vehicle dynamics, as developed in the present article, were employed. To properly visualize how the vehicle’s motion was affected by the performance of the propulsion system, three sets of simulations were conducted:

• Set 1: The caudal fin was set to move at different flapping frequencies and dispositions. Based on the defined motion, thrust was estimated. This set of simulations was conducted to determine the required flapping frequencies to achieve forward thrust $\left(F_X\right)$ when the caudal fin presented full- and mid- range flapping amplitudes. Moreover, the moment produced when biasing flapping was also estimated.
• Set 2: Open-loop simulations were conducted. The differences among the swimming performance when the BAUV presented sideways and dorsoventral flapping were analyzed. The main attributes from each swimming style are described below.
• Set 3: The waypoint guidance closed-loop system to regulate the BAUV’s navigation was simulated. During these simulations, several waypoints were introduced. The vehicle’s heading angle correction, its translational motion, and the propeller’s behavior were iteratively tracked.

To test the accuracy of the designed guidance strategy, the Euclidean distance between desired targets and the vehicle’s final position was computed. Further, to compute its efficiency $E_{rel}$, the BAUV’s real trajectory (and total distance traveled) was compared to the length of a straight line between waypoints. This parameter helped in measuring the relation between the attained performance and an ideal scenario. The ideal trajectory was considered as a straight line between the vehicle’s departure position $\left(x_i,y_i,z_i\right)$ and the target $\left(x_d,y_d,z_d\right)$. The ideal distance $I_D$ was measured as the Euclidean distance between such points, while the real total distance $R_D$ traveled by the BAUV was measured according to the trajectory $T\left(t\right)$ described during the simulation.

The vehicle’s real trajectory $T\left(t\right)$ was described by the $\left(x,y,z\right)$ coordinates obtained over time, and the total distance traveled $R_D$ was computed by integrating the BAUV speed as follows:

$$R_D={{\displaystyle \int }}_{t_i}^{t_f}\sqrt{{\left(x^{\prime }\left(t\right)\right)}^2+{\left(y^{\prime }\left(t\right)\right)}^2+{\left(z^{\prime }\left(t\right)\right)}^2 }dt \left(\mathrm{m}\right);$$

(54)

where $x^{\prime }\left(t\right), y^{\prime }\left(t\right)$, and $z^{\prime }\left(t\right)$ are the time derivatives from vehicle’s position, $t_i$ is the departure time and $t_f$ is the time when the BAUV reached the desired goal.

Finally, the efficiency of the designed guidance strategy was obtained as follows:

$$E_{rel}=\frac{I_D}{R_D}\times 100 \left(\%\right);$$

(55)

where:

$$I_D=\sqrt{{\left(x_d-x_i\right)}^2+{\left(y_d-y_i\right)}^2+{\left(z_d-z_i\right)}^2 } \left(\mathrm{m}\right).$$

(56)

3. Results and Discussion

Figure 11 presents graphs of the force produced when the propeller started oscillating. The thrust levels along the vehicle’s centerline $x_b$ increased once the flapping frequency increased. This means that to push the vehicle and achieve forward momentum, it was necessary to make the caudal fin oscillate more rapidly. During the simulations, the frequency increased by steps of 0.5 Hz every four seconds. As shown in Figure 11, an estimation of force $F_X$ was conducted based on the total amplitude oscillation of the caudal fin $\left(-30°\le \beta \le 30°\right)$, with a flapping frequency ranging from 0.5 to 2 Hz. Figure 11a describes the caudal fin’s angular position over time, and Figure 11b shows the force produced from such motion.

Figure 12 shows the results from the same simulated experiment but only considering biased flapping (i.e., a half fin workspace). Compared to full-range flapping, biased flapping required higher frequencies to push the vehicle forward. Then, based on these results, the propeller was set to flap at 2 Hz when moving at its full range and set to 3.5 Hz when only moving half span.

Since the vehicle’s path is controlled based on the correction of its heading angle, the generation of turning moment is of paramount importance. Figure 13 shows moment $N_Z$ over the $z_b$ axis based on the flapping performance of the propeller. In Figure 13a, three different arrangements for flapping disposition are shown. First, the caudal fin followed a full-range trajectory, and was then biased to each side. As determined in the force analysis, the oscillating frequency of the fin was set to 2 and 3.5 Hz for full- and half-range flapping, respectively. Figure 13b shows how the moment direction changed when biased flapping was induced in the vehicle’s tail. Hence, biasing the caudal fin oscillations changes the vehicle’s orientation, and thus, correction of the BAUV’s heading angle toward desired targets on the horizontal plane may be achieved.

For the second set of simulations, Figure 14 shows the open loop trajectory of the BAUV when swimming sideways. For this analysis, the propeller was set to flap at full port-starboard range at 2 Hz. The starting position and orientation of the BAUV was set at an arbitrary value, i.e., $\mathit{\eta }=\left[1,0,0,0,0,-\frac{\mathsf{\pi }}{2}\right]$, as established in the examples in Figure 6. In Figure 14a, the helicoidal trajectory described by the vehicle is shown from an isometric view, while Figure 14b presents the top view with a circular path on the horizontal plane. Once the propeller starts oscillating, thrust makes the vehicle move forward and dive. Also, during sideways flapping, a lateral force is produced, coupling force components in the surge and sway directions. This event causes the vehicle to start turning over yaw angle $\psi$ in the direction of positive surge and sway. Since the propeller kept the same flapping parameters throughout the simulation, the vehicle kept turning, describing a circular path over the $\left(x,y\right)$ plane. The turning behavior of the BAUV could be corrected by changing the moment produced over the yaw axis to make it rotate in different directions. The linear displacement of the BAUV is shown in Figure 14c–e, while its angular displacement is shown in Figure 14f–h.

In Figure 15a, flapping performance $\beta$ of the fin during the sideways swimming open-loop simulation is described. Figure 15b–d describe the linear force that was exerted in sideways swimming mode. In Figure 15e–g, the rotational moment over $\left(x_b,y_b,z_b\right)$ exerted by fin’s motion is shown. The thrust and moment exerted by the propeller presented consistent behavior throughout the whole simulation. To avoid the loss of detail in graphs of forces and moments, only a 10 s span is shown. During sideways flapping, linear forces $F_X$ and $F_Y$ couple, causing motion in directions $x_b$ and $y_b$. Also, this swimming mode tends to produce moment $N_Z$ over $z_b$. Then, based on the way in which thrust and moment were produced by the designed propulsion system during sideways flapping, enhanced steering maneuvers over the $\left(x,y\right)$ plane could be achieved.

Figure 16 shows the results of setting the BAUV on an open loop mission while swimming dorsoventrally. For a 60 s simulation, the propeller was set to flap at $-30°\le \beta \le 0$ during the first 30 s, and then at $0\le \beta \le 30°$. This allowed us to change thrust direction in order to visualize how the vehicle reacted. When setting the BAUV to swim dorsoventrally, the vehicle also achieved forward thrust, and better control over diving was obtained when compared to sideways flapping (Figure 16a). Figure 16b shows that by applying this swimming method, the BAUV tends to describe straight trajectories. Moreover, the vehicle started diving when the propeller was set to move on the upper side of its workspace (center to top). While diving, buoyant forces acting on the vehicle affected its motion, driving it back over the $z$ axis, as seen on Figure 16e. When the flapping direction changed to lower-side workspace (center to bottom) of the fin, the thrust exerted on the vehicle allowed it to reach the surface sooner. Also, the BAUV’s linear and angular displacements (Figure 16c–h) presented more stability than those obtained during lateral swimming. This was because of the up-and-downward motion of the caudal fin during dorsoventral flapping, where lateral forces produced over the vehicle’s body were diminished, and vertical forces could be properly directed.

Figure 17a describes the flapping behavior of the propeller during the dorsoventral swimming simulation. This time, the caudal fin was biased to each side. Figure 17b–d shows the linear force that was exerted by the fin’s motion with both configurations. Forward thrust $F_X$ was not affected once the flapping direction changed. However, dorsoventral flapping exerted forces along the $z_b$ axis. According to the biased flapping produced (center to top/center to bottom), $F_Z$ forces could be properly directed. This caused the vehicle to achieve motion along surge and heave. Finally, dorsoventral swimming mode tended to produce moment over the vehicle’s $y_b$ axis. Figure 17e–g shows the moment computed in the simulation.

As seen in Figure 14 and Figure 16, a BAUV with thunniform motion will be capable of producing thrust efficiently, thereby attaining excellent forward speed. Nonetheless, a lack of efficient maneuvering will occur when turning and diving. Since the caudal fin’s motion is the only parameter responsible for changing moment and force directions, the vehicle takes longer to react. Furthermore, since changing between flapping modes while swimming is not an energy-efficient strategy, the designed propulsion system only follows one swimming style while performing a task. Then, if the mission comprises following straight trajectories while gaining depth, the best swimming option will be the dorsoventral. Otherwise, sideways swimming would be the best option.

To overcome the low degree of maneuverability of the BAUV, a strategy to regulate the flapping performance of the vehicle’s tail should be developed. Since the vehicle moves forward thanks to the thrust produced by its caudal fin, correcting the direction of such motion becomes a control task. Then, if the flapping performance of the propeller is correctly updated, unwanted behavior may be reduced.

Figure 18 shows the results obtained from the waypoint guidance strategy in the third set of simulations. Figure 18a illustrates how the BAUV was able to correct its erratic behavior over the horizontal plane while performing sideways swimming. The heading angle correction was achieved by biasing the flapping motion from the caudal fin, as shown in Figure 18b. Then, by maintaining forward speed and properly updating the flapping parameters after each control cycle, it was possible to correct the orientation of the vehicle relative to the desired goal. In Figure 18d, tracking errors between ${\psi }_d$ and $\psi$ are plotted. The Euclidean distance was iteratively reduced (Figure 18c), and once the vehicle reached the waypoint area, the propeller was turned off. Figure 18e–g shows the linear coordinates of the BAUV throughout the mission. Once the propeller was shut down, the BAUV floated back to surface (as seen in Figure 18g).

Figure 19 shows the ways in which forces and moments changed during the simulated trial described in Figure 18. The caudal fin was able to produce a constant forward thrust, $F_X$, throughout the mission. Figure 19a shows how thrust decayed to zero once the propeller was shut down at 36 s. Figure 19b,f shows the changing behavior of $F_Y$ and $N_Z$ over time.

Table 4. Accuracy and efficiency obtained by the proposed guidance system during the first simulation.

Initial-Final Coordinates	$\mathbf{Ideal} \mathbf{Distance} {\mathit{I}}_{\mathit{D}}$	$\mathbf{Real} \mathbf{Distance} \mathbf{Traveled} {\mathit{R}}_{\mathit{D}}$	$\mathbf{Efficiency} {\mathit{E}}_{\mathit{r}\mathit{e}\mathit{l}}$	$\mathbf{Accuracy} {\mathit{E}}_{\mathit{d}}$
$\left(1,0,0\right)-\left(-5,15,0\right)$	16.1555 m	19.2415 m	83.96%	0.1150 m

states the accuracy and efficiency obtained from the waypoint guidance strategy applied in this scenario.

Finally, Figure 20a shows the BAUV performance when several waypoints were introduced into the guidance system. The vehicle corrected its course (Figure 20d) toward each of the established targets, ultimately reaching waypoint regions with a proximity $E_d$ of less than 0.2 $\mathrm{m}$ (Figure 20c). Based on these results, the vehicle demonstrated its ability to reach different positions on the horizontal plane; however, its accuracy and performance were limited by its thunniform locomotion. Nonetheless, the proposed propulsion system can overcome such limitations. By using a parallel mechanism, the vehicle’s tail flapping may achieve different speeds and present different behaviors (Figure 20b). This design for motion generation provides efficient and well-directed thrust. For maneuvering, having a system that can expeditiously bias the propeller is as important as thrust generation. Then, by the development of proper control strategies over the flapping performance of the robotic fish tail, the BAUV will achieve the desired outcomes. Figure 20e–g shows the vehicle’s behavior while en route to the desired coordinates. In

Table 5. The accuracy and efficiency achieved by the proposed guidance system during the second simulation.

Initial-Final Coordinates	$\mathbf{Ideal} \mathbf{Distance} {\mathit{I}}_{\mathit{D}}$	$\mathbf{Real} \mathbf{Distance} \mathbf{Traveled} {\mathit{R}}_{\mathit{D}}$	$\mathbf{Efficiency} {\mathit{E}}_{\mathit{r}\mathit{e}\mathit{l}}$	$\mathbf{Accuracy} {\mathit{E}}_{\mathit{d}}$
$\left(1,0,0\right)-\left(-5,15,0\right)$ m	16.1555 m	19.2415 m	83.96%	0.1150 m
$\left(-5,15,0\right)-\left(0,30,0\right)$ m	15.8114 m	19.2798 m	82.01%	0.1580 m
$\left(0,30,0\right)-\left(5,15,0\right)$ m	15.8114 m	17.5905 m	89.89%	0.1850 m
$\left(5,15,0\right)-\left(0,0,0\right)$ m	15.8114 m	18.7031 m	84.54%	0.0730 m

, the accuracy and efficiency of the designed strategy based on the distance traveled by the BAUV toward each waypoint are detailed.

As seen in Set 3 of the simulated scenarios, the designed strategy did not consider either the vehicle’s orientation upon arrival or the reduction of the path-tracking error when compared to straight-line trajectories. Some novel controllers may be able to reduce the path-tracking error by controlling the production of forces and moments from the vehicle’s propeller. This might be achieved by the incorporation of robust controllers such as dynamic sliding mode methods or adaptive and/or intelligent controllers. However, in this case, thrust and moment exertion depends on the vehicle’s biomimetic features, which are the result of the proposed propulsion system. In other words, the vehicle’s thrust and maneuvering capabilities are limited by the flapping performance of its caudal fin and the way in which it is driven by the designed parallel mechanism. Even if a control strategy could properly define the required flapping frequency, amplitude, and bias to reduce path-tracking error, the navigation task will depend exclusively on the ability of the propulsion system to produce such movement at the desired moments. From the final simulated scenario, the BAUV showed an average travel-efficiency of about 85% in its navigation toward several waypoints. Also, the average accuracy when arriving at targets was 0.1327 m. These results were achieved thanks to the flapping parameters which were applied to the guidance system, and the control which was applied to the parallel robotic system inside the propeller.

4. Conclusions

In this article, a method for the mathematical modeling of a 6 DOF BAUV with thunniform locomotion was developed. The kinematics and dynamics of the robotic fish and its inner components were studied to comprehend how moments and forces were produced to move the vehicle. For the trajectory analysis, the vehicle’s hydrodynamics were determined. Further, a waypoint guidance strategy based on the flapping performance of the BAUV’s caudal fin was implemented. The vehicle was able to advance and maneuver thanks to its propulsion system, which was responsible for driving the vehicle’s tail at different flapping frequencies and with different dispositions.

The use of a parallel mechanism inside the propulsion system allowed the vehicle’s caudal fin to oscillate at different frequencies and in different ways. Due to this novel feature, the thunniform locomotion of the designed vehicle was enhanced, due to its ability to switch between two swimming modes, i.e., lateral and dorsoventral. From the obtained results, it was seen that the swimming mode directly impacted the vehicle’s performance and its trajectory. While sideways flapping allowed motion on the horizontal plane, dorsoventral flapping presented the advantage of producing stable motion along the vertical plane. Depth and planar motions were achieved by selecting the appropriate swimming mode. Hence, the BAUV’s motion was determined by controlling its tail, since this was the only parameter responsible for generating thrust and moment. Finally, based on the limitations of this locomotion, strategies were developed to correct the course of the vehicle according to the desired goals. With proper control over the flapping frequency and bias of the fin, the designed BAUV achieved forward thrust and moment, allowing it to orient itself relative to several established targets. The waypoint guidance strategy was found to be the optimal means by which to overcome the unwanted performance characteristics that the vehicle tended to present.

In future work, a faster, more robust, and adaptive smart controller for BAUV path tracking needs to be incorporated. More complex techniques such as deep learning and/or neuro-fuzzy algorithms may enhance the process of selecting the flapping parameters. Improved control over parameters such as frequency, bias and amplitude will reduce cruising cost during missions, thereby improving travel efficiency. Hence, complementary controllers to enhance the performance of the parallel mechanism inside the propulsion system should be analyzed. The presence of external disturbances based on different aquatic ecosystems should be incorporated in future analyses to define the expected performance of the proposed design. Further, future work should focus on path planning strategies to allow the vehicle to not only reach preestablished goals, but to properly define which trajectories are the most convenient in the search of such targets.