Coordinate-free Lie-group-based modeling and simulation of a submersible vehicle

1.   Introduction

Submersible vehicles represent a remarkable product of combined efforts in engineering and applied science, allowing researchers to explore and interact with the world beneath the Earth's oceans, rivers, and lakes ^[3,21,35]. These machines, designed to operate underwater, have evolved from simple submersibles to sophisticated vessels that facilitate scientific research, transportation, and even underwater tourism.

The category of unmanned submersible vehicles encompasses a diverse range of autonomous or semi-autonomous underwater systems, with notable representatives being autonomous underwater vehicles (AUVs) and unmanned surface vessels (USVs). These innovative vehicles operate autonomously or under remote guidance. Autonomous underwater vehicles, in particular, are self-propelled, programmable machines designed to execute specific underwater tasks autonomously ^[15]. Such vehicles are deployed for oceanographic surveys, underwater mapping, and deep-sea exploration. AUVs are equipped with a suite of sensors, cameras, and advanced navigation systems, which enable them to carry out missions independently and with precision.

Submersible vehicles technology relies on mathematical models to navigate the complex and often harsh underwater environments ^[5]. These models are essential for predicting how submersible vehicles will behave, ensuring safe operations, enabling to predict consumption and costs, and achieving mission objectives. In fact, mathematical models constitute the backbone of submersible vehicle technology ^[13] and provide the analytic foundation for designing, operating and controlling such devices. The continuous refinement of these models, coupled with advances in sensor technology and computational power, ensures that submersible vehicles continue to evolve toward reliable and widely-usable models.

At the heart of submersible vehicle modeling lies hydrodynamics and fluid mechanics. These mathematical constructions describe how the vehicle interacts with water, accounting for factors such as buoyancy and drag ^[24]. Understanding these forces is crucial for describing the vehicle's motion, stability, and propulsion. Mathematical equations of mechanical dynamics and kinematics govern the motion of submersible vehicles. These equations take into account the forces acting on the vehicle, its mass distribution as well as the resulting accelerations and velocities, and help in predicting how the vehicle will respond to changes in control inputs and external forces.

Submersible vehicles operate in diverse underwater environments, each with its own set of challenges ^[1,35]. Mathematical models help accounting for environmental variables such as water temperature, salinity, and currents and help vehicles adapt to changing conditions and make informed decisions about navigation and data collection. When submersible vehicles are deployed for specific missions, mathematical models aid in mission planning and optimization. Model-based navigation algorithms consider mission objectives, vehicle capabilities, and constraints to determine the optimal path, waypoints, and data collection strategies. Prior to deployment, mathematical models are used in simulation environments to test a vehicle's performance under plausible scenarios as well as to build-up real-time simulators ^[14,25]. Such preventive activity enables researchers to refine a vehicle's design, control algorithms, and mission plans, reducing the risk of occurrence of unexpected issues during actual operations.

In the present paper, our goal is to lay down a mathematical model of a submersible vehicle formulated in a coordinate-free manner, expressed in the language of Lie groups. To construct one such mathematical model, the minimal-action (Hamilton) principle for a rigid body will be invoked, and adjusted by the d'Alembert virtual-work principle to account for external factors, specifically drag forces and control inputs. A review about mathematical modeling of a complex mechanical device on Lie groups may be accessed in ^[10]. Since a body submerged in a non-ideal fluid interacts with the fluid in a complex manner, the concept and the standard model of 'added mass' will be recalled and profited of, the dynamics of which influence the vehicle's motion. The effects explained by the added mass will be incorporated into the Lagrangian of the vehicle to build up a comprehensive model of motion.

The resulting equations of motion will be presented in the form of (non-pure) Euler-Poincaré relations defined over the Lie group of rigid motions. The present approach to modeling may, in fact, be regarded as belonging to the modern 'geometric framework' ^[20]. Euler-Poincaré equations stand as a powerful mathematical framework at the intersection of physics and mathematics, particularly in the study of dynamical systems and fluid mechanics. Euler-Poincaré equations constitute a cornerstone of modern applied mathematics as they are able to offer a profound understanding of the behavior of complex systems.

In the present paper, a six-degree-of-freedom vehicle model shall be considered. Such an approach constitutes a deviation from the typical assumption that the motion in the horizontal plane is decoupled from the remaining degrees of freedom ^[20]. In fact, coupling between degrees of freedoms may prove advantageous in situations where a vehicle is by design (or becomes by accident) under-actuated.

To simulate the behavior of such a submersible vehicle, a numerical simulation technique based on Lie-group integrators will be reviewed briefly and applied. Lie-group integration is a mathematical technique that holds profound importance in the realm of differential equations and the study of continuous symmetries. Rooted in the profound mathematical concepts of Lie groups and group actions, this method offers a unique approach to solving differential equations and to understanding the behavior of dynamic systems. Lie-group integration deals seamlessly with the symmetries that underlie physical phenomena.

The present paper is organized as follows. Section 2 introduces a prototypical submersible vehicle, recalls the fundamental assumptions behind the presented mathematical model and presents a derivation of the Lagrangian function associated to such kind of vehicle. Section 3 recalls the minimal-action principle to infer the equations of motion of the dynamics of a rigid body tempered with the d'Alebert principle of virtual work to account for non-conservative external forces; the resulting optimization principle is tackled by a variational approach, which gives rise to a system of two equations of motion plus two reconstruction equations. Section 4 presents detailed calculations to apply the devised equations of motion to the studied submersible vehicle. Section 5 illustrates the results of numerical simulations about elementary maneuvers such as 'station keeping'. Section 6 concludes the paper.

2.   Lagrangian function associated to a submersible vehicle

The present section deals with a submersible vehicle by defining its associated Lagrangian function which embodies information on the vehicle and on the surrounding viscous fluid. In particular, Section 2.1 presents a review of the mathematical notation and of the general properties of the motion of rigid bodies as well as fundamental hypotheses to model a submersible vehicle. Section 2.2 introduces a kinetic energy function for the vehicle as well as for the surrounding fluid mass. Section 2.3 introduces a potential energy function to model buoyancy effects. Section 2.4 gathers the above results into a single Lagrangian function to explain the dynamics of a vehicle-fluid complex.

2.1.   Notation, general properties and fundamental hypotheses

In the present study, the translation space $\mathbb{R}^3$ is equipped with the standard Euclidean metric, namely, $\langle u, v\rangle^{ \mathbb{R}^3} : = u^\top v$, for any $u, v\in \mathbb{R}^3$. The special orthogonal group of three-dimensional rotations is defined as the matrix group $\mathrm{SO}(3) : = \{R \in \mathbb{R}^{3\times 3} \mid R^\top R = R R^\top = I_3, \ \det(R) = 1\}$ and its associated Lie algebra is defined as $\mathfrak{so}(3) : = \{\Omega \in \mathbb{R}^{3\times 3} \mid \Omega + \Omega^\top = 0_3\}$, where the symbol $^\top$ denotes matrix transposition, the quantity $I_3$ represents a $3\times 3$ identity matrix, the symbol $\mathbb{R}^{3\times 3}$ represents the space of $3\times 3$ real-valued (unstructured) matrices, while the quantity $0_3$ represents a $3\times 3$ zero matrix. In addition, the Lie group $\mathrm{SO}(3)$ shall be equipped with its canonical metric $\langle \Omega, \Psi\rangle^{ \mathfrak{so}(3)} : = \mathrm{tr}(\Omega^\top\Psi)$ for every $\Omega, \Psi\in \mathfrak{so}(3)$ and the matrix space with the Euclidean metric $\langle M, N\rangle^{ \mathbb{R}^{3 \times 3}} : = \mathrm{tr}(M^\top N)$ for any pair $M, N\in \mathbb{R}^{3\times 3}$.

In order to lay down a mathematical model of the unrestricted motion of a rigid body, it is necessary to fix an inertial (or absolute) reference frame, that shall be denoted as $\mathcal{F}_\mathrm{E}$, as well as a vehicle-fixed reference frame that will be denoted as $\mathcal{F}_\mathrm{V}$. A schematic of a work-class submersible vehicle is shown in Figure 1. The position of the origin of the reference frame $\mathcal{F}_\mathrm{V}$ with respect to the origin of the reference frame $\mathcal{F}_\mathrm{E}$ will be denoted by the vector $p\in \mathbb{R}^3$, while the orientation of the vehicle-fixed reference frame with respect to the inertial reference frame will be denoted by a rotation matrix $R\in \mathrm{SO}(3)$. In the present framework, there will be no reference to the NED (north-east-down) inertial reference system routinely invoked in maritime navigation.

A distinguishing feature of the present formulation is that there will be no need to fix any specific basis nor coordinate system since the primary variables to describe the motion of the submersible vehicles make up the pair $(R, p)\in \mathrm{SO}(3)\times \mathbb{R}^3$, where the state space may be summarized as the special Euclidean group $\mathrm{SE}(3)$. The present mathematical model is based on a coordinate-free Lie-group representation, hence no individual components of matrices will be accessed or made reference to. For this reason, within the present paper there will be no reference to the 18-variable SNAME (1950) notation ^[28] for marine vessels.

The motion of a submersible vehicle will be modeled under a number of assumptions about the viscous fluid that it navigates. Ocean currents are considered stationary and irrotational, namely, it is assumed that the fluid surrounding the vehicle moves at a constant speed $v_\mathrm{c}\in \mathbb{R}^3$ with respect to the inertial reference frame and flows along straight streams (see, for example, ^[6]). Wave-induced motion is supposed to be negligible, namely the natural undulatory motion of the fluid surrounding the vehicle is assumed to be of negligible magnitude, as much as the inertial motions due to the rotation of the Earth.

To what concerns the vehicle itself, it is assumed that the vehicle-fixed reference frame is centered, namely that the origin of the reference frame $\mathcal{F}_\mathrm{V}$ coincides with the center of gravity (CoG) of the vehicle. However, it is hereafter assumed that the center of gravity of the vehicle does not necessarily coincide to its center of buoyancy (CoB) and that the vehicle is not necessarily neutrally buoyant. Such assumptions imply that the gravitational pull and the buoyancy push do no compensate and may yield a vertical (typically upward) thrust as well as a righting moment that acts on the vehicle's hull.

In the dynamics of submersible vehicles, one crucial concept that plays a significant role is that of "added mass" ^[12]. Added mass, often referred to as hydrodynamic or virtual mass, represents the apparent increase in inertia that a submersible vehicle experiences due to its motion through a viscous fluid medium, typically sea/ocean water. This phenomenon bears a profound impact on the vehicle's dynamics, stability, and control. Added mass arises from the displacement of fluid as a vehicle moves through it. When the vehicle accelerates or changes its velocity, it causes the surrounding fluid to be displaced, effectively creating a virtual mass that appears to be added to the vehicle's actual mass ^[8]. Such additional mass is dynamic and its effects depend on the vehicle's shape, speed, and the density of the surrounding fluid. Added mass affects the vehicle's inertia matrix, influencing its response to accelerations and changes in velocity. This impact is particularly significant in high-speed or agile submersible vehicles. Added mass plays a role in hydrodynamic interactions between multiple submersible vehicles or between a vehicle and other objects in the fluid, even influencing collision dynamics and cooperative missions. Neglecting added mass might lead to inaccurate predictions. Mathematically, added mass is often represented through additional terms in the equations of motion for submersible vehicles affecting the vehicle's accelerations and angular velocities. The numerical characteristics of added mass may be determined experimentally or through dedicated simulation software ^[19].

In computational simulations of submersible vehicle dynamics, added mass is considered alongside other hydrodynamic factors, such as drag and lift forces. Such factors are essential for accurately modeling a vehicle's behavior in different operating conditions and environments. To what concerns how to treat added mass in the modeling process, it is assumed that the added mass behaves as a rigid body that takes the same attitude of the vehicle and whose attitude is described by the same rotation matrix $R$. This is tantamount to assuming that the principal axes of the added mass coincide to the principal axes of the vehicle ^[20].

In the equations of rotational motion, a number of operators will appear, whose definitions are recalled here for the sake of self-containedness. The matrix commutator is denoted and defined as $[M, N] : = MN-NM$ for any pair $M, N\in \mathbb{R}^{3\times 3}$. The adjoint operator $\mathrm{ad}: \mathfrak{so}(3)\times \mathfrak{so}(3)\to \mathfrak{so}(3)$ is defined as $\mathrm{ad}_\Omega(\Psi) : = \Omega\Psi-\Psi\Omega = [\Omega, \Psi]$, for any $\Omega, \Psi\in \mathfrak{so}(3)$ (notice that the matrix commutator is more general than the adjoint operator since the commutator is defined for any square matrix pair). Its dual with respect to the chosen inner product in $\mathfrak{so}(3)$, denoted as $\mathrm{ad}^\star$, is equal to $- \mathrm{ad}$. A skew-symmetrization operator $\sigma: \mathbb{R}^{3\times 3}\to \mathfrak{so}(3)$ is defined as $\sigma(M) : = \frac{1}{2}(M - M^\top)$ for any matrix $M\in \mathbb{R}^{3\times 3}$. A 'diamond' product $\diamond: \mathbb{R}^3\times \mathbb{R}^3\to \mathfrak{so}(3)$ is defined as $u\diamond v : = -\sigma(uv^\top)$ for any $u, v\in \mathbb{R}^3$. Both the diamond product and the adjoint operator are counterparts of the cross product in the ordinary space, translated to the Lie algebra $\mathfrak{so}(3)$.

A noteworthy observation is that, in a certain circumstance, the anti-commutator of two matrices may be written in terms of the skew-symmetrization operator. In fact, let $S\in \mathbb{R}^{3\times 3}$ denote a symmetric matrix and $\Omega\in \mathfrak{so}(3)$ a skew-symmetric one: In this case, it is straightforward to prove the identity $S\Omega+\Omega S = 2\, \sigma(S\Omega)$. A further noteworthy observation is that, as long as skew-symmetric matrices are concerned, the commutator of two matrices may be written in terms of the skew-symmetrization operator $\sigma$. In fact, for any pair $\Omega, \Psi\in \mathfrak{so}(3)$, it holds that $[\Omega, \Psi] = 2\, \sigma(\Omega\Psi)$.

2.2.   Kinetic energy of the submersible vehicle and of the added mass

Let us denote by $\mathrm{d} m$ a small volume element of the vehicle's body and by $\xi\in \mathbb{R}^3$ its position with respect to the reference frame $\mathcal{F}_\mathrm{V}$. Its position with respect to the inertial reference frame is then given by $x = p + R \xi\in \mathbb{R}^3$. The kinetic energy of the vehicle with respect to the inertial reference frame may be written as:

where $\mathcal{V}$ denotes the volume of space occupied by the vehicle at a given time instant and an over-dot denotes time derivative.

The linear velocity of such a mass element is expressed as $\dot{x} = \dot{p} + \dot{R} \xi$ since $\dot \xi = 0$ in the frame $\mathcal{F}_\mathrm{V}$. From the fact that $\|\dot{x}\|^2 = \mathrm{tr}(\dot{x} \dot{x}^\top)$, it follows that:

where

The quantity $\hat{J}_\mathrm{V}$ represents the non-standard inertial tensor of the vehicle which is related to its standard inertial tensor by the relation $\hat{J}_\mathrm{V} = \frac{1}{2} \mathrm{tr}(J_\mathrm{V})I_3 - J_\mathrm{V}$. Notice that, while the standard inertia tensor is a symmetric, positive-definite matrix ($J_\mathrm{V} > 0$), the non-standard inertia is not necessarily positive-definite. The scalar $m_\mathrm{V} > 0$ denotes the vehicle's mass (thought of to be concentrated in its CoG) and the vector $c_\mathrm{V}\in \mathbb{R}^3$ denotes the position of the center of gravity of the vehicle in the reference frame $\mathcal{F}_\mathrm{V}$. By the assumptions made, it holds that $c_\mathrm{V} = 0$, therefore, the kinetic energy of the vehicle reads

where the first term represents a translational-energy contribution while the second term represents a rotational-energy contribution to the total kinetic energy of the vehicle.

By a strict analogy, building on prevailing literature on the subject ^[12,20], the kinetic energy of the added mass is written as:

where the vector $v_\mathrm{c}\in \mathbb{R}^3$ denotes the local velocity of the fluid current expressed in the frame $\mathcal{F}_\mathrm{E}$, the matrix $M_\mathrm{A}\in \mathbb{R}^{3\times 3}$ denotes an added-mass translational inertia tensor and the matrix $\hat{J}_\mathrm{A}\in \mathbb{R}^{3\times 3}$ denotes its (non-standard) rotational inertia tensor. The translational component of the kinetic energy was written by assuming that the added mass translates with a relative velocity quantified as $\dot{p} - v_\mathrm{c}$. The added mass is only defined in the vehicle-fixed frame $\mathcal{F}_\mathrm{V}$, henceforth, its expression in the reference frame $\mathcal{F}_\mathrm{E}$ reads $RM_\mathrm{A}R^\top$. On the basis of the expression (2.5), one may gather that a purely formal interpretation of added mass is the proportionality constant that relates the kinetic energy associated to the fluid surrounding the vehicle to the (square of the) vehicle's speed ^[4].

According to ^[12], it may be assumed that the added-mass matrix $M_\mathrm{A}$ is constant and strictly positive. Moreover, when modeling operations at great depth, it may be assumed that the added-mass matrix is symmetric. As it is underlined in ^[12], the kinetic-energy approach to model the effects of the added mass is compatible with Kirchhoff's equations for a fluid surrounding a submerged body.

Hence the kinetic energy of the vehicle/added-mass ensemble takes the form

where $\hat{J} : = \hat{J}_\mathrm{V} + \hat{J}_\mathrm{A}$ denotes its total rotational inertia. It is worth underlining that, within the present paper, there will be no invocation of the small-angle approximation for the attitude matrix, as opposed to, e.g., ^[14]. In fact, the invoked Lie-group formulation affords to tackle seamlessly with large-tilting-based maneuvers.

2.3.   Potential energy of the vehicle

Let $|\mathcal{V}|$ denote the volume of fluid displaced by the vehicle, $g$ the (local) acceleration of gravity and $\rho$ the fluid density (both regarded as constant parameters in the present endeavor). The magnitudes of the gravity pull and of buoyancy push acting on a submerged system are expressed as $m_\mathrm{V}g$ and $\rho g|\mathcal{V}|$ respectively. While the gravity pull acts on the CoG of the vehicle, the buoyancy push acts on its CoB. In addition, let us denote by $\mathrm{e}_z$ the direction of gravity in the inertial reference system (which is assumed to stay constant in the volume of space interested by a vehicle maneuver) however directed upward, and by $\ell\in \mathbb{R}^3$ the position of the CoB with respect to the vehicle-fixed reference frame. Normally, it holds that $\ell^\top \mathrm{e}_z > 0$, which characterizes a 'bottom-heavy' vehicle ^[20].

The potential energy function corresponding to the gravitational and buoyancy restoring forces, written in the inertial reference frame, reads

The first term represents the gravitational potential, while the second term represents Archimedes' potential. The above expression does not include any term associated to the added mass because the restoring forces acting on it are assumed to be in equilibrium.

A further potential-energy term arises from the thrust of the fluid current that acts on the hull of the vehicle and tends to align the body of the vehicle to the current speed ^[34]. Such term reads

and gives rise to a deflection moment that tends to align a direction vector $\eta\in \mathbb{R}^3$, fixed in the body reference frame, to the current speed $R^\top v_\mathrm{c}$ as perceived by the vehicle. The direction $\eta$ may, for instance, coincide to one of the principal axes of inertia of the vehicle and, most likely, it is such that $\eta^\top \mathrm{e}_z = 0$. Its amplitude and direction depend on the shape of the vehicle's hull as it quantifies the resultant of the pressure of the fluid current on the hull of the submerged vehicle.

The potential energy function associated to the submersible vehicle then casts in the expressive form

The first term accounts for the vertical thrust on the vehicle's body due to the mismatch between the gravity pull and the buoyancy push, while the second term accounts for the righting torque that tends to keep the sagittal ($x-y$) plane of the vehicle horizontal with respect to the inertial reference frame.

Whenever a vehicle is designed and operated to be neutrally buoyant, the buoyancy push balances the gravity pull, namely, $m_\mathrm{V} = \rho |\mathcal{V}|$, hence, the body of the vehicle is not subjected to any vertical thrust due to such restoring forces (namely, one such vehicle will neither sink nor rise spontaneously). Such formal assumption is often encountered in the literature. In the present paper, it is rather assumed that $m_\mathrm{V} < \rho |\mathcal{V}|$, which is beneficial to make sure a submersible vehicle may be rescued in case of malfunctioning thanks to a slight thrust directed upward.

2.4.   Lagrangian of the vehicle/added-mass system and reduced Lagrangian

The Lagrangian function $L$ associated to the vehicle/added-mass ensemble is defined as the sum of the kinetic energy components minus the potential energy, namely, $L : = K - U$. In the present case, the Lagrangian function takes the expression

The expression of the above Lagrangian function may be simplified by introducing two auxiliary variables, namely $\Omega : = R^\top\dot{R}\in \mathfrak{so}(3)$, a matrix that represents the instantaneous angular velocity in the vehicle-fixed reference frame, and $v : = R^\top\dot{p}\in \mathbb{R}^3$, a vector that represents the instantaneous translational velocity in the vehicle-fixed frame. On the basis of such new variables, the Lagrangian function, referred to as 'reduced Lagrangian', recasts as

By a slight abuse of notation, the same symbol has been used to denote the Lagrangian as well as its reduced version. It is interesting to notice that the terms that come from the kinetic energy only depend on the new variable-pair $(\Omega, v)$, which is referred to as 'symmetry' in these terms; the symmetry is however broken by gravity and by the fluid current, since the terms coming from a potential energy depend on the original variable-pair $(R, p)$.

In the body-fixed reference frame, the velocity of the ocean current is expressed by the term $R^\top v_\mathrm{c}$, which coincides to the expression given, for instance, in ^[27].

3.   Lagrange-d'Alembert principle and calculus of variation for a submersible vehicle

The aim of the present section is to recall the general Lagrange-d'Alembert principle (LdAp) of mechanics and apply it to a submersible vehicle to get the corresponding equations of motion of a so-called 'simple mechanical system' ^[9]. In particular, Section 3.1 recalls the Lagrange-d'Alembert principle to describe, in variational terms, the dynamics of a rigid body subjected to external solicitations. Sections 3.2 to 3.4 present the related equations of motion together with a short discussion on unmodeled interaction phenomena.

3.1.   The Lagrange-d'Alembert principle of mechanics applied to the submersible vehicle

The Lagrange-d'Alembert principle is a fundamental concept in classical mechanics that provides a powerful framework for describing the behavior of mechanical systems (see, e.g., ^[32] and references therein). It serves as an alternative approach to Newton's laws of motion and is particularly useful in formalizing complex problems involving constraints and generalized coordinates. The LdAp finds applications in various fields of physics and engineering, including robotics, aerospace engineering, and structural mechanics. The Lagrange-d'Alembert principle is based on the idea of virtual work. It states that for a mechanical system in equilibrium (i.e., experiencing no acceleration), the sum of the virtual work effected by all forces and constraint reactions is zero for any infinitesimal virtual displacement consistent with the constraints.

For a continuous system such as a vehicle submerged in a viscous fluid, the Lagrangian is a function of the trajectory of the system in its state-space bundle. Therefore, the LdAp is formulated as a variational problem, namely,

where the symbol $\delta$ denotes a variation corresponding to a shift from a trajectory to an infinitely close nearby trajectory, $f_\mathrm{P}\in \mathbb{R}^3$ and $f_\mathrm{F}\in \mathbb{R}^3$ denote the thrust from the propellers and the resistance due to friction, respectively, while the terms $\tau_\mathrm{P}\in \mathfrak{so}(3)$ and $\tau_\mathrm{F}\in \mathfrak{so}(3)$ denote the torques due to the propellers and the torque due to rotational friction, respectively. In the present endeavor, the thrust and torque exerted by the propellers are assumed to be control inputs and the time index is assumed to span the interval $\Theta : = [t_0, \ t_\mathrm{f}]$, where $t_\mathrm{f} > t_0 \geq 0$.

3.2.   Calculations of the involved variations

All variables introduced so far depend on the two fundamental quantities $p$ and $R$ which describe in an exhaustive way the motion of the vehicle as a rigid body. In order to calculate the variations of the reduced Lagrangian, let us introduce two families of curves, namely, a $\mathbb{R}^3$-valued family $p(t, s)$ and a $\mathrm{SO}(3)$-valued family $R(t, s)$, with $s\in[-\epsilon, \epsilon]$, where $\epsilon > 0$, and $t\in[t_0, \ t_\mathrm{f}]$. The time-variable $t$ parameterizes each curve belonging to such continuous family, a member of which is selected upon assigning a value to the index $s$. The curve corresponding to the index value $s = 0$ represents the actual trajectory followed by the mechanical system in the state-space bundle.

On the basis of the above two families of curves, the following quantities are defined:

The variational principle (3.1) is then cast explicitly as

The families of curves $R(t, s)$ and $p(t, s)$ are chosen so that the endpoints $R(t_0, s)$, $R(t_\mathrm{f}, s)$, $p(t_0, s)$ and $p(t_\mathrm{f}, s)$ are independent of the value of the variable $s$, which implies that the variations $\delta\Omega$ and $\delta v$ vanish at the boundary of the interval $[t_0, \ t_\mathrm{f}]$. The variations are further assumed to be smooth, although arbitrary within the time interval $\Theta$. The requirement of smoothness, in particular, includes the commutation of partial derivatives, namely,

The reduced Lagrangian appears ultimately as a function of four variables, namely, $L = L(R, p, \Omega, v)$, therefore, it holds that

where $\frac{\partial L}{\partial \Omega}\in \mathfrak{so}(3)$ is a skew-symmetric derivative, while the matrix $\frac{\partial L}{\partial R}\in \mathbb{R}^{3\times 3}$ denotes the Euclidean gradient of the Lagrangian with respect to the matrix-valued variable $R$. In the following, the above terms will be evaluated one by one with the aim of writing any inner product in terms of either the variation $\delta\Omega$ or the variation $\delta v$, whichever is appropriate.

From the definitions (3.2), it is immediate to see that $\frac{\partial p}{\partial s} = R\, \delta v$, therefore,

Likewise, since $\frac{\partial R}{\partial s} = R\, \delta \Omega$, it holds that

From the definitions (3.2), it follows that

As a result, the corresponding term in the variation of the reduced Lagrangian reads

Analogously, to what concerns the last term in the variation of the Lagrangian function, it is immediately seen that

The corresponding term in the variation of the Lagrangian reads

Setting $s = 0$, the sought variation of the Lagrangian around the actual trajectory was gleaned to be

The second and fourth terms in the above right-hand side sum give rise to two terms in the expression (3.1) that need to be rewritten as inner products involving explicitly the variations.

Invoking the properties of the trace operator and the formula for integration by parts leads to

In the last term, the following property was invoked: For any pair of matrices $M\in \mathbb{R}^{3\times 3}$ and $\Omega\in \mathfrak{so}(3)$ it holds that $\mathrm{tr}(\Omega^\top M) = \mathrm{tr}(\Omega^\top\sigma(M))$, namely, in one such inner product only the skew-symmetric component of a matrix $M$ effectively plays a role.

Analogously, usage of integration by parts and of the notion of adjoint operator with respect to the inner product leads to the expression

Gathering all terms that appear in the Lagrange-d'Alembert conservation Equation (3.1), and factoring out the independent variations, yields an extended energy/work conservation law for the submersible vehicle. Such conservation law embodies terms that describe the mechanical structure of the moving vehicle as well as terms that express what causes or impedes its motion. Expressed in integral form, the gleaned conservation law takes the form

Since each sum of terms on the left-hand side must vanish irrespective of the variations, which are arbitrary, the above relation yields the following system of differential equations, which represents the conservation principle (3.1) expressed in differential form

Such equations are the analogous of Newton's equations of dynamics for the rotational and for the translational components of motion. In addition, the following reconstruction equations

afford the determination of the positional-orientational variable-pair $(R, p)$ for a vehicle under study.

3.3.   Generalized-momentum-based representation

An alternative way to restate the equations of dynamics is by the introduction of a generalized translational momentum $q : = \frac{\partial L}{\partial v}$ and a generalized rotational (or angular) momentum $\Pi : = \frac{\partial L}{\partial \Omega}$ (see, for example, ^[9] and references therein). On the basis of such new variables, whenever the fluid surrounding the vehicle may be considered locally at rest (namely, $v_\mathrm{c}\equiv 0$), the equations of motion recast as

In one such representation, each term may be attributed a physical meaning. The term $\mathrm{ad}_\Omega^\star(\Pi)\in \mathfrak{so}(3)$ represents a centrifugal torque due to rotational inertia, while the term $q\diamond v\in \mathfrak{so}(3)$ represents a centrifugal force due to translational inertia. Both terms are generally due to an uneven distribution of mass within the vehicle/added-mass ensemble. In fact, in a simple body whose translational inertia would be represented by a scalar mass, it would hold that $q\parallel v$ and hence $q\diamond v = 0$; due to the added mass burdening the vehicle, such simplistic assumption could hardly hold true. Likewise, in a body whose mass distribution is even, such as in a uniform sphere, it would hold that $\Pi = \mbox{scalar}\times \Omega$, and hence $[\Omega, \Pi] = 0$; even such assumption would be out of place in the present development, hence the term $\mathrm{ad}_\Omega(\Pi)$ in the equation of rotational motion needs to be duly taken into account.

The term $-\sigma\left(R^\top\frac{\partial U}{\partial R}\right)\in \mathfrak{so}(3)$ represents the righting moment due to the buoyancy push. The term $\Omega^\top q\in \mathbb{R}^3$ denotes a Coriolis-type fictitious force due to the fact that the translational motion of the vehicle's body is observed from a rotating reference frame $(\mathcal{F}_\mathrm{V})$. In those instances where the submersible vehicle navigates slowly, the Coriolis/centripetal forces might be neglected, as in ^[31], although one such shortcut has not been taken in the present paper. Moreover, the term $- R^\top\frac{\partial U}{\partial p}\in \mathbb{R}^3$ represents the hydrostatic thrust arising by the combined action of gravity and Archimedes' thrust.

The first two equations in (3.18) are expressed in the vehicle-fixed reference frame, implying that the non-conservative forces $f_\mathrm{P}, \, f_\mathrm{F}$ and torques $\tau_\mathrm{P}, \, \tau_\mathrm{F}$ need to be expressed in the frame $\mathcal{F}_\mathrm{V}$ as well. This is certainly the case for the components due to friction, and it is made to be the case for the components due to the propellers since, in an autonomous vehicle, the control actions are based on readouts from sensors that are fixed to the body of the vehicle itself.

3.4.   Further interaction phenomena not included in the present study

The Lagrange-d'Alembert principle (3.1) was expressed in a fairly general way, yet it misses a number of interaction phenomena that may affect the motion of a marine vehicle. A number of such phenomena was modeled in the specialized literature and gives rise to analytic expressions to be included in the variational principle (3.1), while further interaction phenomena might not be formalized and may hence be included in the energy/work conservation law as disturbances, either random or periodic (to represent, e.g., waves of ocean currents) ^[31,36].

Examples of complex interaction phenomena studied in the specialized literature include the Froude-Krylov and diffraction force in sway, the dependence of the added mass and of damping terms from the frequency of the undulatory phenomena ^[14], as well as the wave excitation force acting on a vehicle's hull ^[27]. Such phenomena are of lesser importance in the operation of submerged vehicles and were henceforth not considered in the present mathematical model. They were considered, albeit quantitatively, for instance in the applied paper ^[36].

A further interesting instance of interaction phenomena between the body of the vehicle and the surrounding fluid occurs when underwater vehicles are endowed with 'fins' (such as the ones installed in the tail of a torpedo). Artificial fins play a crucial role in the hydrodynamics and control systems of underwater vehicles, contributing to their overall performance and functionality. Fins help in maintaining stability and balance by adjusting the vehicle's orientation, as they may provide stability in different axes and help prevent unintended rotations. Some underwater vehicles use control surfaces, including fins, to regulate their depth. By adjusting the lift generated by the fins, the vehicle can ascend or descend in the water. The lift thrust and torque generated by fins may be modeled, to some extent, by expressions involving the velocity $v$ of the vehicle and its rotation speed $\Omega$, as illustrated, for example, in ^[4,17,22].

4.   Detailed calculations and complete mathematical model

The aim of the present section is to convey the results of detailed calculations of the mechanical forces and torques intervening in the equations of motion and to summarize the obtained relations.

4.1.   Gradients, fiber derivatives and related terms in the equations of motion

On the basis of the functional expression of the Lagrangian (2.10) that summarizes the dynamics of the submersible vehicle according to the Lagrange-d'Alembert principle of mechanics, the four derivatives that appear in the equations of motion (3.16) may be calculated.

The gradient $\frac{\partial L}{\partial p}$ that quantifies the hydrostatic thrust takes the expression

hence the hydrostatic force itself takes the expression

It is worth underlining that the product $-R^\top \mathrm{e}_z$ denotes the direction of gravity as seen from the vehicle-fixed reference system.

The gradient $\frac{\partial L}{\partial R}$ that quantifies the righting moment, the drag moment of the added mass and the deflection moment of the fluid current on the vehicle's hull takes the expression

The combined action of the righting (or hydrostatic) torque, of the added-mass-drag torque, and of the deflection moment due to the fluid current, results in the mechanical torque

Whenever the current velocity is negligible, namely, $v_\mathrm{c}\approx 0$, the drag moment as well as the deflection moment are likewise negligible and the resulting torque may be approximated as just the righting moment due to flotation thrusts $\rho g |\mathcal{V}|\sigma(R^\top \mathrm{e}_z\ell^\top)$.

The fiber derivative $\frac{\partial L}{\partial \Omega}$ that quantifies the intrinsic angular momentum of the vehicle has expression

where the total non-standard matrix of inertia of the vehicle/added-mass ensemble was defined as $\hat{J} : = \hat{J}_\mathrm{V} + \hat{J}_\mathrm{A}$. This expression appears particularly simple (and the total inertia is constant, in fact) because the added mass was assumed to be co-rotating with the vehicle. As a consequence, it holds that $\frac{ \mathrm{d}}{ \mathrm{d} t}\frac{\partial L}{\partial \Omega} = \sigma(\hat{J}\dot\Omega)$. Likewise, the fiber derivative $\frac{\partial L}{\partial v}$ that quantifies the translational momentum takes the expression

where the total mass of the vehicle-fluid ensemble was defined as $M : = m_\mathrm{V}I_3 + M_\mathrm{A}$. As a consequence

It is interesting to notice that the momentum of the vehicle also depends on the momentum of the added mass, since they are assumed to form a sort of compound in which, however, each component translates with a different speed.

In terms of generalized momenta, in the presence of a still fluid (or quasi still, namely, with $\|v_\mathrm{c}\|$ negligible) it would hold that $q = Mv$ and $\Pi = \sigma(\hat{J}\dot\Omega)$. On a moving fluid, however, even in a state of rest the vehicle possesses a non-zero translational momentum, due to the drift velocity of the added mass, along a direction opposite to that of the current. In fact, the notion of added mass is far from trivial.

4.2.   Force and torque related to friction

The effects of friction due to a fluid, no matter if it is an ideal fluid such as air or a non-ideal fluid such as sea/ocean water, are generally not negligible. The translational as well as the rotational friction are generally modeled as second-order or third-order polynomials of the corresponding velocities. In the present paper, a 'quadratic polynomial' approximation shall be opted for, to be adapted to the present Lie-group type representation. The quadratic term has indeed a specific meaning within the mathematical model (3.16), as it is meant to describe hydrodynamic forces experienced by a submerged vehicle due to the normal pressure on its hull ^[9]. Cubic terms, on the other hand, may be justified by experimental evidence ^[9] and may help describe unmodeled frictional effects.

To what concerns the translational friction, the expression adopted here reads

where $L_\mathrm{f}, \, Q_\mathrm{f} > 0$ are $3\times 3$ diagonal damping matrices, $v$ denotes again the translational speed of the vehicle in the reference frame $\mathcal{F}_\mathrm{V}$, the symbol $|\cdot|$ denotes component-wise absolute value, and the symbol $\circ$ denotes component-wise (Hadamard) vector multiplication.

It is interesting to remark that, since the friction is due to the interaction between the body of the vehicle and the surrounding fluid and since the mechanical force of friction is (nonlinearly) proportional to the velocity of the body in relation to the velocity of the fluid, the above terms include the relative velocity $v - R^\top v_\mathrm{c}$. Such assumption carries a somewhat counterintuitive consequence. Assume, in fact, that at some moment in time the vehicle is at rest ($v = 0$) while the fluid current is not ($v_\mathrm{c} \neq 0$). In this case, the damping term (4.8) takes the expression $f_\mathrm{F} = L_\mathrm{f}\, R^\top v_\mathrm{c} + Q_\mathrm{f}(|R^\top v_\mathrm{c}|\circ (R^\top v_\mathrm{c}))$, which indeed appears as an active force pushing the vehicle along the direction of the current.

Likewise, the rotational friction torque takes the form

where $L_\tau, Q_\tau > 0$ are $3\times 3$ diagonal damping matrices, $\Omega$ denotes the angular speed of the vehicle in the reference frame $\mathcal{F}_\mathrm{V}$, the symbol $|\cdot|$ denotes again component-wise absolute value, while the symbol $\circ$ denotes component-wise (Hadamard) matrix multiplication. The above expression may be recast in terms of the operator $\sigma$ as

thanks to a property recalled in Subsection 2.1.

It pays to notice that the product $|\Omega|\circ\Omega$ preserves the structure of skew symmetry and provides a quadratic-type dependence on the rotation matrix that preserves the direction of rotation, so as to make sure that the friction torque acts in the opposite direction to that of rotation since the damping matrices are assumed to be positive-definite. The expressions (4.8) and (4.10) reflect the corresponding expressions suggested, e.g., in the paper ^[36].

4.3.   Propellers action

The body of the submersible vehicle is subjected to active thrust exerted by a number of propellers. Such propellers are generally positioned along the hull of the vehicle at fixed locations and orientations in such a way to warrant complete maneuverability. The resultants of the thrusts and of the corresponding torques were indicated in the Equation (3.1) as $f_\mathrm{P}$ and $\tau_\mathrm{P}$, respectively.

A control algorithm will calculate, at every moment, the value of the pair $(f_\mathrm{P}, \tau_\mathrm{P})$ necessary to attain a desired motion. A built-in hardware or software, specifically tailored to the vehicle's propeller system configuration, will convert such pair into a set of commands to be imparted to the propellers such that the resultant of their combined thrust matches the inputs from the control algorithm.

Since the total mechanical thrust and torque that the propellers exert on the body of the vehicle are generated from within the vehicle, the variables $f_\mathrm{P}$ and $\tau_\mathrm{P}$ are always expressed in the vehicle-fixed reference frame $\mathcal{F}_\mathrm{V}$. An example of a $(f_\mathrm{P}, \tau_\mathrm{P})$ pair for a torpedo-like vehicle, due to a set of propellers that exert a constant thrust along the direction of motion and torque about the axis of motion, while the other components of thrust and torque are modulated by commanding a stern angle and a rudder angle, is explained in the paper ^[17].

The propellers considered in the present mathematical model are only those fixed to the hull of the vehicle. A submersible vehicle may be endowed with internal actuators whose purpose would be to stabilize its motion. The effects of a reaction wheel coaxial to one of the main axes of the vehicle are surveyed in Appendix 8.

4.4.   Complete mathematical model

The complete set of equations that constitute the mathematical model of the submersible vehicle under consideration takes the expression:

For the sake of clarity, a short description of each equation has been included to recall its physical/mechanical meaning within the devised mathematical model. Notice that the equations of dynamics do not depend explicitly on the variable $p$, which may be however recovered to evaluate the location of the vehicle at any given time.

A special mention goes to the term $\sigma((v - R^\top v_\mathrm{c})(v - R^\top v_\mathrm{c})^\top M_\mathrm{A})$, which represents the dragging action of the fluid mass on the hull of the vehicle. In order to express such action in a compact way, the observation that $\sigma(vv^\top) = 0$ has been exploited. Whenever the velocity of the vehicle coincides to the velocity of the added mass considered in the vehicle's reference frame, namely whenever $v = R^\top v_\mathrm{c}$, such dragging effect vanishes. It is also interesting to notice that in the case that the added mass is isotropic, namely, the matrix $M_\mathrm{A}$ coincides to a scaled identity, such contribution is null.

The Table 1 recaps the variables used within the mathematical model and their physical/mechanical meaning, while Table 2 summarizes the constant parameters (either scalar, vectorial or matricial) that appear in the mathematical model (4.11).

In addition, Table 3 summarizes the definition of a few mathematical operators that were made use of to express the equations of motion in a concise way.

The Table 4 summarizes those terms that quantify the hydrostatic and the totality of hydrodynamic thrusts and torques considered in the present mathematical model.

A simplified version of the model (4.11) is given in the Appendix 9 under a number of simplificative assumptions.

5.   Numerical simulations

The present section aims at illustrating a few characteristics of the devised mathematical model through a number of basic numerical simulations. Such simulations concern station-keeping as well as further simplistic tasks. (No complex maneuvers are considered here since devising a strategy to control the present model is a task for a future research endeavor.)

A constant non-zero absolute fluid current velocity is assumed, which is accounted for by the velocity vector $v_\mathrm{c} = [0.4\ -0.3\ 0]^\top$ m/s (i.e., the fluid current is purely longitudinal since it bears no ascensional component). Such numerical values are compatible with a real-world setting as mentioned, for example, in ^[18].

As a vehicle type to experiment with, the 'Seaking I' considered in detail in the paper ^[36] was taken as a reference model. There exist a variety of submersible vehicles, designed for different purposes, exhibiting different shapes and features and, hence, described by different sets of parameter values ^[29].

5.1.   Parameters values

The parameter values, adapted from the research paper ^[36], are summarized in Table 5, where the operator $\mathrm{diag}(\cdot, \cdot, \cdot)$ returns a diagonal matrix whose in-diagonal elements are indicated as arguments.

Several values in the present table differ from the corresponding values in Table I of ^[36] because the present paper makes use of a different representation of the equations of motion; hence, parameters like the rotational inertia matrix and the rotational friction coefficients take a slightly different meaning (and, hence, appear with different values). The vehicle-fluid-compound mass, in particular, results to be $M = \mathrm{diag}(4640, 4136, 5500)$ (kg).

5.2.   Numerical recipes

The formulation of the model based on a Lagrange-d'Alembert approach is compatible with sophisticated numerical simulation methods, such as the 'variational integrator' described in ^[23] and with more general schemes well-suited to simulate mechanical systems ^[7]. Within the present paper, the simplest numerical simulation schemes will be invoked, which are generally suitable to provide consistent simulation results.

The mathematical model (4.11) is summarized as follows

where the $\mathbb{R}^3$-vector field $a$ denotes the translational acceleration of the vehicle due to the forces that act on its body and the $\mathfrak{so}(3)$-vector field $A$ denotes the rotational acceleration of the vehicle due to the combined effect of the torques that act on its body.

In order to solve numerically such system of four differential equations, it is worth observing that the first three equations are based on the underlying vector (linear) spaces $\mathbb{R}^3$ and $\mathfrak{so}(3)$; hence, such first-order differential equations are solved numerically by a Euler-type stepping method, namely, by the iterations

where $h > 0$ denotes a stepsize (expressed in seconds) and the arrow '$\leftarrow$' denotes assignment of a new value upon erasing the old value.

The fourth equation in the system (5.1) is solved iteratively by a Lie-group integration method based on matrix exponential 'exp', namely, by the iteration ^[11]

which guarantees the structure of the attitude matrix $R$ to be preserved over iterations. The following numerical simulation results are referred to a sampling period $h = 0.1$ seconds, which is assumed to be a suitable timescale given the size of the submerged vehicle and of the velocities at stake.

More complex and precise numerical methods were developed in the specialized literature. In the present endeavor, Euler-like stepping methods are deemed of sufficient precision to simulate the behavior of a submersible vehicle.

5.3.   Numerical simulation: Drifting with motors off

The first simulation refers to the case that the vehicle is initially moving with random translational as well as rotational speed and that all direct-current-motor-driven propellers are switched off. In this case, the submerged vehicle is subjected only to the gravity/buoyancy forces/moments, to inertial/virtual (Coriolis) forces/moments, and to the thrust of the fluid current on its hull.

On the basis of the stated conditions, it is expected that the vehicle will rise due to the prevalence of the upward buoyant thrust over the downward gravity pull, and that the submerged vehicle will drift sideways and its longitudinal axis will be deflected due to the thrust of the fluid current. In addition, it is expected that the vehicle's sagittal plane will become horizontal due to the righting moment of the buoyant thrust through a rolling/pitching-based settlement. Such expected effects are clearly visible in Figure 2.

Due to the combined action of the righting moment and of the current thrust on the hull, it may be predicted that that the pitch and roll angle will go to zero, while the yaw angle will take a specific, non-zero value. Since the vehicle is initially tilted, the thrust of the fluid current is unbalanced over the vehicle's body and produces a rotation of the vehicle with respect to its vertical ($z$) axis. When the sagittal plane of the vehicle becomes horizontal and the profile offered by the vehicle to the fluid current is just right, the current thrust is only able to push the vehicle but is unable to cause any further rotation along the vertical axis. The asymptotic rotation angle along the vertical axis may be quantified by imposing the yaw-equilibrium condition

for a counterclockwise rotation of an angle $\gamma$ along the vertical axis expressed by the matrix

The zero-deflection condition (5.4), applied to a current velocity $v_\mathrm{c} = [v_\mathrm{c}^x\ v_\mathrm{c}^y\ 0]^\top$, gives rise to the prediction

which, throwing in numbers, gives $\tan\gamma\approx 1.3333$, which is good agreement with the result of this simulation $\tan\, \gamma \approx 1.3336$.

A detailed analysis of the ascensional motion due to the buoyant thrust is provided in Appendix 7.

5.4.   Numerical simulation: Station keeping

Station keeping consists of taking and then keeping a submersible vehicle in a given location, described by the position vector $q_\mathrm{s}\in \mathbb{R}^3$, and in a specific attitude in which the frame $\mathcal{F}_\mathrm{V}$ coincides to the frame $\mathcal{F}_\mathrm{E}$. In order to do so, a basic control algorithm was designed, whose purpose is to cancel the vehicle's own dynamics and to win over all disturbances (current thrust on the hull, current dragging) and to take the vehicle in the desired location/attitude (namely, the desired attitude matrix is $R_\mathrm{s} = I_3$).

Station keeping as well as guidance requires the vehicle to activate the internal reaction thrusters. In order to make sure that the controller will be properly tuned so as not to exceed the mechanical capabilities of the propellers that the 'Seaking I' submersible vehicle is equipped with, it is necessary to know the maximum thrust and torque values provided by the propulsion system. Such values were adapted from Table I of the paper ^[36] and are shown in Table 6. Interestingly, the ascensional thrust is noticeably unbalanced in comparison to the other thrust values, and the yaw-component of the propulsion torque is noticeably larger than the other torque values.

The control law to achieve station keeping was designed by two basic principles: cancellation of part of the internal dynamics and deployment of a proportional-derivative control strategy. Such control law was not chosen to be particularly optimized since this would not be the main purpose of the present paper nor of the present section devoted to numerical simulations. The devised control law is expressed as follows:

In the above equations, the scalar coefficient $K_\mathrm{p} > 0$ denotes the gain of the proportional component of the position controller, while the scalar coefficient $K_\mathrm{d} > 0$ denotes the gain of the derivative component of the position controller. Likewise, the scalar coefficient $K_\mathrm{a} > 0$ denotes the gain of the proportional component of the attitude controller, while the scalar coefficient $K_\mathrm{b} > 0$ denotes the gain of the derivative component of the attitude controller. Also, $\mathrm{e}_x : = [1\ 0\ 0]^\top$ was set as the unit vector indicating the direction of the $x$ axis in the reference system $\mathcal{F}_\mathrm{E}$.

In the present simulation, it was set $q_\mathrm{s} = [0\ 0\ -5]$ (m). To what concerns the control law, it was set $K_\mathrm{p} = 50$, $K_\mathrm{v} = 1000$, $K_\mathrm{a} = 500$, $K_\mathrm{b} = 100$. The results of this simulation, displayed in terms of positional/attitudinal variables, are shown in Figure 3. Within the timeframe of the simulation, the desired location, as well as the desired attitude, are achieved and kept. The initial transients observed in the middle- and right-hand panels are due to the initial translational/rotational speeds and attitude set at random.

The values of the total propulsion thrust and mechanical torque calculated by the control law (5.7) are displayed in Figure 4. As it may be verified directly from this figure, the values of the propulsion thrusts and torques do not exceed the physical limits reported in Table 6. The numerical values look certainly large, though compatible with the large tonnage of the 'Seaking I' vehicle.

The values of the thrusts corresponding to the surge and sway stay nonzero to counter the action of the ocean currents, while the thrust component corresponding to heave stays nonzero to counter the buoyant thrust. The values resulting from the simulation are in excellent agreement with the values predicted from the first of Equation (5.7), obtained by setting $\Omega = 0$, $v = 0$ and $R = I_3$, namely:

As for the the values of the torque related to pitching and rolling, these vanish to zero since there are no intrinsic nor external torques to win along these axes, while the torque corresponding to the yawing axis needs to stay non-zero to win the deflecting action of the ocean current. Even these values in the results displayed in Figure 4 are in excellent agreement with the values predicted from the second of Equation (5.7), obtained by setting $\Omega = 0$, $v = 0$ and $R = I_3$, namely:

Side note: The rationale for the devised control law may be summarized as follows. The sum of terms $(\Omega^\top M_\mathrm{A} - M_\mathrm{A}\Omega^\top)R^\top v_\mathrm{c} + (m_\mathrm{V} - \rho|\mathcal{V}|)g R^\top \mathrm{e}_z + L_\mathrm{f}\, (v - R^\top v_\mathrm{c}) + Q_\mathrm{f}(|v - R^\top v_\mathrm{c}|\circ (v - R^\top v_\mathrm{c}))$ serves to cancel some components of the translational dynamics of the vehicle, notably those related to the current velocity. The term $K_\mathrm{p}\, R^\top(q_\mathrm{s} - q)$ in the expression of the propulsion thrust serves to take the vehicle over the station location, while the term $- K_\mathrm{v}\, v$ helps it to stay put. The sum of terms $- \sigma((v - R^\top v_\mathrm{c})(v - R^\top v_\mathrm{c})^\top M_\mathrm{A}) - \sigma(R^\top v_\mathrm{c}\eta^\top)$ serves to cancel some components of the rotational dynamics of the vehicle, and again an emphasis is put on those depending on the fluid current pressure. Conversely, the terms $K_\mathrm{a}\, \sigma(R^\top \mathrm{e}_x \mathrm{e}_x^\top)$ and $K_\mathrm{b}\, \sigma(\Omega^\top R^\top \mathrm{e}_x \mathrm{e}_x^\top)$ serve to take the yaw angle to zero and to ensure it also stays put. In particular, the proportional component of the control torque was designed to mimic the deflection component of the potential, $U_\mathrm{D}$, namely,

where the scalar $\mathrm{e}_x^\top R^\top \mathrm{e}_x$ coincides to the cosine of the angle between the $x$ axis of the vehicle and the absolute $x$ axis. The control torque only affects the yaw angle since the inclination of the sagittal plane is determined by the buoyancy thrust.

6.   Conclusions

The aim of the present paper was to propose a Lie-group view of mathematical modeling and simulation of a submersible vehicle. To construct this mathematical model, the minimal-action (Hamilton) principle for a rigid body was invoked, adjusted by the d'Alembert virtual-work principle to account for external factors such as drag forces and control inputs. Such an approach is also able to account for the interaction between a vehicle and the non-ideal fluid that it is submerged into. The resulting equations of motion were presented in the form of (non-pure) Euler-Poincaré relations defined over the Lie group of rigid motions. To simulate the behavior of such a submersible vehicle, a numerical simulation technique based on simple Euler and specialized Lie-group integrators was reviewed briefly and applied.

From manned submarines that probe the depths of the ocean to undersea remotely operated vehicles, submersible vehicles play a pivotal role in scientific missions. While it is true that submersible vehicles have expanded our access to deep-sea environments, certain extreme depths of the ocean remain challenging to reach even with advanced technology. For example, the Mariana Trench's Challenger Deep, the deepest known point in the ocean, poses a serious challenge to exploration due to its extreme depth and conditions. Looking ahead, submersible vehicles continue to evolve, with developments in autonomous underwater vehicles and advancements in deep-sea exploration technology. Their critical roles in scientific research underscore their enduring significance in understanding the aquatic world. Key to their deployment is the ability to model and simulate their behavior as accurately as possible.

Use of AI tools declaration

The author declares he has not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgment

The present paper was completed while the author was visiting the Keio University, Yagami Campus, Japan (Peng's Lab) during January-February 2024. Support from this organization is gratefully acknowledged.

Conflict of interest

The author declares no conflict of interest.

Appendix

A.   Model of pure ascensional motion

The motion of a submersible vehicle, albeit possibly complex, may be broken down into three 'elementary' motions, namely, pure vertical ascent (heave), pure horizontal advancement (surge) and pure rotation along its vertical axis (yaw) ^[2]. Under the hypothesis that the propellers are switched off and that the transient caused by the deflection action of the fluid current has petered out, it holds that $\Omega = 0$ and $R = R_\gamma$, namely, the vehicle presents a fixed orientation due to the pressure of the fluid current on its hull. The motion reduces to the only translational component, whose speed is governed by the equation

The term $R_\gamma v_\mathrm{c}$ is constant and represents the effect of the fluid current on the profile of the hull that is at a just-right angle to ensure that there will happen no further rotation around the $z$ axis. After the transient action of the deflection moment has petered out, the orientation of the sagittal ($x-y$) plane of the vehicle tends to adapt to the fluid current speed, namely, $v^x = (R_\gamma v_\mathrm{c})^x$ and $v^y = (R_\gamma v_\mathrm{c})^y$, therefore, it holds that $\dot v^x = \dot v^y = 0$. Only the component $v^z\neq 0$ differs from $(R_\gamma v_\mathrm{c})^z = 0$, where $v^z : = v^\top \mathrm{e}_z$. As a result, the velocity of the vehicle composes as

hence it holds that $Mv = MR_\gamma v_\mathrm{c} + M^{zz}v^z \mathrm{e}_z$. As a consequence, the relation (A.1) simplifies to

The above equation describes a pure ascensional motion of the vehicle due to the prevalence of buoyant thrust over the gravity pull, slowed down by the friction to the surrounding fluid.

B.   Effects of a reaction wheel

Let us assume that a cylindrical wheel of mass $m_\mathrm{W}$, height $a$ and radius $r$ is positioned right in the center of the vehicle, so as not to change its center of gravity, along the $x$ axis of the vehicle. The inertia tensor of a cylindrical object is $J_\mathrm{W} = \mathrm{diag}\left(\tfrac{1}{2}m_\mathrm{W}r^2, \tfrac{1}{12}m_\mathrm{W}(3r^2+a^2), \tfrac{1}{12}m_\mathrm{W}(3r^2+a^2)\right)$. A reaction wheel is generally made to spin at constant speed, hereafter denoted as $\Omega_\mathrm{W}\in \mathfrak{so}(3)$.

The considered reaction wheel possesses a kinetic energy and a potential energy, which are calculated as

where the variables $v$, $p$ and $\Omega$ take the same meaning as in the main text and $\hat{J}_\mathrm{W}$ denotes the non-standard inertia tensor of the wheel that exhibits the expression $\mathrm{diag}\left(\tfrac{1}{12}m_\mathrm{W}a^2, \tfrac{1}{4}m_\mathrm{W}r^2, \tfrac{1}{4}m_\mathrm{W}r^2\right)$.

The effects of a reaction wheel on the equations of motion may be evaluated by computing the derivatives

These relations reveal that the reaction wheel contributes the total mass of the vehicle as well as its total weight and total inertia (along every direction of rotation of the vehicle's body). The main point is however the presence of the constant term $\hat{J}_\mathrm{W}\Omega_\mathrm{W}$ to the total angular momentum of the vehicle, which is able to counter the deflection effects of the ocean current since it makes it difficult to change the direction of the $x$ axis of the vehicle ^[33].

C.   Simplified model under pitch/roll stabilization and motion along the local level plane

The devised 6-degree-of-freedom model may be simplified under the assumptions, which indeed hold in some cases of interest, that a submersible vehicle's autopilot embodies a roll and pitch stabilization algorithm and that the vehicle navigates a local level plane ^[26]. These assumptions imply that the attitude of the vehicle may be described by a pure rotation around the vertical axis, namely, $R = R_\gamma$ as defined in (5.5), and that the ascensional velocity is zero, namely, $v^\top \mathrm{e}_z = 0$. Further assumptions that are often made are that the fluid current velocity is negligible, namely, $v_\mathrm{c} = 0$ and that the vehicle is neutrally buoyant.

On the basis of the assumptions made, the mathematical model (4.11) simplifies noticeably since several terms turn out to be zero. In particular, the first two equations simplify to

where the control inputs must coherently be such that $f_\mathrm{P} \parallel \mathrm{e}_z$ and $\tau_\mathrm{P} \parallel \Omega$.

The complete mathematical model (4.11) may be simplified at will with the aim of focusing one's attention on particular aspects of navigation or phases of motion such as, for instance, diving (downward motion along the plane $x-z$) ^[16,30].