Energy Consumption Modeling for Underwater Gliders Considering Ocean Currents and Seawater Density Variation

Abstract

Energy management is a critical and challenging factor required for efficient and safe operation of underwater gliders (UGs), and the energy consumption model (ECM) is indispensable. In this paper, a more complete ECM of UGs is established, which considers ocean currents, seawater density variation, deformation of the pressure hull, and asymmetry of gliding motion during descending and ascending. Sea trial data are used to make a comparison between ECMs with and without the consideration of ocean currents, and the results prove that the ECM that considers the currents has a significantly higher accuracy. Then, the relationship between energy consumption and multiple parameters, including gliding velocity relative to the current, absolute gliding angle, and diving depth, is revealed. Finally, a simple example is considered to illustrate the effects of the depth-averaged current on the energy consumption.

1. Introduction

After about 30 years of development, underwater gliders (UGs) have become an efficient ocean observation instrument, which are widely applied due to their high endurance, low energy consumption, and low cost without expensive supporting vessels [1].

Generally, UGs are designed to be very compact to facilitate deployment and concealment, which limits the space for the battery onboard, making it difficult for a UG to carry an energy intensive sensor, which is, however, necessary in some situations. For example, a UG cannot obtain its precise position without using a high-energy consumption navigation sensor, such as the Doppler velocity logger, which will result in route deviation and increased energy consumption [2]. Therefore, improving energy efficiency without adding an extra battery is critical to improving the performance indicators of UGs, such as motion accuracy, endurance and duration [3]. In addition, as battery power failures have been the second leading cause of glider mission failures [4], accurately predicting the energy consumption of the UG and then implementing energy management are very important and challenging factors for operating the glider efficiently and safely.

To solve the above problem, energy consumption models (ECMs) have been explored and established in many studies, which can reflect the relationship between energy consumption and various parameters of UGs. In the design and perfection stage of the glider, the ECM is usually used to evaluate the rationality of the design parameters, such as hydrodynamic coefficients. During the operation of the glider, the ECM is usually employed to optimize mission deployment, mainly involving path planning, motion parameter optimization, and sensor arrangement.

Sun et al. optimized the shape of a Blended-Wing-Body glider by an ECM to improve the hydrodynamic performance and increase the gliding range [5]. Zhang et al. designed a miniaturized underwater glider, and a simple ECM was used to understand the trade-off between the glide performance and energy cost [6]. Zhou et al. established an ECM for a deep-sea glider to evaluate the benefits of a hybrid buoyancy-regulating system they proposed [7]. Yang et al. used the ECM of the Petrel-L glider to verify the effectiveness of their shape optimization method [8]. Cao et al. optimized the 3D Dubin curve path of UGs by employing a highly nonlinear and coupled ECM as the heuristic function [9]. Using the same ECM in [9], Cao et al. obtained the rendezvous trajectories with minimal energy consumption for multiple gliders [10]. To address the problem of energy-optimal path planning in the presence of flow fields for UGs with complex dynamics, an ECM considering the change of the heading angle was presented by Lee et al. [11]. Woithe and Kermer investigated various sensor control techniques of a Slocum glider with the sensor energy models [12]. Yang et al. optimized the motion parameters of the Petrel-L glider by combining an ECM and the inner penalty-function method [13]. Song et al. established an ECM of the Petrel-II glider to identify the key parameters that affect the energy consumption and gliding range, with which the optimal motion parameters of the Petrel-II glider were determined [14]. Based on an ECM, Wu et al. proposed a multi-objective optimization method to make a trade-off between the energy utilization rate and motion accuracy [15]. Li et al. [16] established an ECM and sailing-range equation suitable for the electric, thermal and hybrid propulsion underwater gliders, and then performed analysis to investigate the effects of configuration parameters and navigation parameters on the endurance of three types of underwater gliders.

The above ECMs are effective in revealing the basic relationship between energy consumption and various parameters of UGs, which is very valuable in the shape design, motion parameter optimization and mission planning of UGs. The existing ECMs, however, may be inaccurate or inconsistent in a real situation in actual marine environments, with ocean currents and gradient of seawater density which have a significant impact on dynamic behaviour and the energy consumption of UGs [17,18,19]. The ocean currents will cause changes in the velocity relative to the current and thus the variation of hydrodynamics. The change of seawater density with depth will lead to buoyancy loss (BL) [1,7]. The dynamic behaviour of the UG when influenced by the ocean current and density variation has been fully investigated by Fan and Woolsey [18] and Yang et al. [1], respectively. However, the energy consumption of UGs under the above two factors has not been studied.

In the present work, we first establish a dynamic model of UGs, considering the depth-averaged current (DAC), and the BL caused by seawater density variation and deformation of the glider’s pressure hull. Then, a more complete and realistic ECM is presented based on open-loop control, which means the control parameters are fixed during the ascending or descending phase of one work cycle and will not be affected by the motion state and the ambiance of the glider. In addition, asymmetry of the glider motion during ascending and descending is taken into account. The established ECM is validated to be effective by sea trials. Finally, several simulations are performed to investigate the relationship between the energy consumption and multiple parameters, as well as the impact of DAC on the energy efficiency of the UGs.

The main contributions of this work can be concluded as follows:

• A more complete and realistic ECM of the underwater gliders is developed, which considers the effects of horizontal DAC, BL and asymmetry of the glider motion during ascending and descending. The ECM is validated by a comparison with sea trial data.
• The relationship between the energy consumption and diving depth, gliding velocity relative to the current and gliding angle is revealed by simulations.
• The effects of horizontal DAC on energy consumption are analyzed by a simple example.

The paper is organized as follows. Section 2 presents the dynamic model of UGs, which considers the ocean currents and BL. In Section 3, the ECM is established and validated by sea trial data of a Petrel-II underwater glider. In Section 4, simulations are carried out to study the relationship between energy consumption and various parameters, and the impact of DAC on the energy efficiency of UGs. Section 5 ends the paper by presenting the main conclusions drawn from this work.

2. Dynamic Model

According to the work of Thomasson and Woolsey [20] and Fan and Woolsey [18], we use Lagrange’s equations to establish a dynamic model of the Petrel-II underwater glider, whose structure and workflow can be found in our previous work [14].

For simplicity, we refer to [21] to give the following assumptions and simplifications:

• The UG has two planes of symmetry (see Section 2.5 for details).
• The dynamic equations of the UG do not include the disturbance forces of waves and wind.

2.1. The Reference Frames

As shown in Figure 1, an inertial frame, a body frame, and a flow frame are established to describe the motion of an underwater glider. The origin of the body frame B₀ (b₁, b₂, b₃) coincides with the centre of buoyancy B₀, and the inertial frame I₀ (i₁, i₂, i₃) is established at a point on the sea surface. i₁, i₂ and i₃ axes point North, East and Down, respectively. The origin of the flow frame V₀ (c₁, c₂, c₃) also coincides with B₀. This flow frame can make the calculation of hydrodynamics more convenient.

Let X = [x, y, z]^T and Θ = [φ, θ, ψ]^T denote the position and attitude of the glider in the inertial frame, respectively. Vectors v = [u, v, w]^T and ω = [p, q, r]^T are defined as the translational velocity and the angular velocity of the glider in the body frame, respectively. In the flow fame, α and β are attack angle and slip angle, respectively.

Here, rotation matrixes are defined to transform the motion parameters from one frame to another.

$${\mathit{R}}_{\mathrm{BI}}=\left[\begin{array}{ccc}\cos \theta \cos \psi & \sin \phi \sin \theta \cos \psi -\cos \phi \sin \psi & \cos \phi \sin \theta \cos \psi +\sin \phi \sin \psi \\ \cos \theta \sin \psi & \cos \phi \cos \psi +\sin \phi \sin \theta \sin \psi & -\sin \phi \cos \psi +\cos \phi \sin \theta \sin \psi \\ -\sin \theta & \sin \phi \cos \theta & \cos \phi \cos \theta \end{array}\right]$$

(1)

$${\mathit{R}}_{\mathrm{IB}}={\mathit{R}}_{\mathrm{BI}}{}^{-1}={\mathit{R}}_{\mathrm{BI}}{}^{\mathrm{T}}$$

(2)

$${\mathit{\Omega }}_{\mathrm{BI}}=\left[\begin{array}{ccc}1& \sin \phi \tan \theta & \cos \phi \tan \theta \\ 0& \cos \phi & -\sin \phi \\ 0& \sin \phi \sec \theta & \cos \phi \sec \theta \end{array}\right]$$

(3)

$${\mathit{R}}_{\mathrm{VB}}=\left[\begin{array}{ccc}\cos \alpha \cos \beta & -\cos \alpha \sin \beta & -\sin \alpha \\ \sin \beta & \cos \beta & 0\\ \sin \alpha \cos \beta & -\sin \alpha \sin \beta & \cos \alpha \end{array}\right]$$

(4)

where R_*# represents the linear velocity transformation matrix from frame * to #, and Ω_BI maps the angular velocity in the body frame to the rate of change of Θ in the inertial frame.

2.2. Mass Distribution and Model of Mechanics

As shown in Figure 2, the total mass of a glider m_v is divided into three parts: stationary mass, moving internal point mass, and variable ballast point mass, which are described in more detail below.

2.2.1. Stationary Mass

The stationary mass m_s is the total mass of the pressure hull, fixed wings, ballasts, and other objects relatively stationary to the body. The stationary mass can be regarded as a rigid body. Let ${\mathit{r}}_{\mathrm{s}}={\left[r_{\mathrm{s}1}, r_{\mathrm{s}2}, r_{\mathrm{s}3}\right]}^{ \mathrm{T}}$ denote the vector from B₀ to the centre of stationary mass, and ${\mathit{l}}_{\mathrm{s}}={\left[l_{\mathrm{s}1}, l_{\mathrm{s}2}, l_{\mathrm{s}3}\right]}^{ \mathrm{T}}$ represent the vector from I₀ to the centre of stationary mass.

2.2.2. Moving Internal Point Mass

The moving internal point mass m_p is the mass of the attitude regulation unit. An internal battery pack with eccentric offset R_p is used as a mass block to adjust the attitude of the glider, and the principle is illustrated in Figure 3. The pitch angle of the glider is adjusted by moving the battery pack back and forth to keep the glider moving forward, and the roll angle is controlled by rotating the battery pack to achieve the turning motion. Let ${\mathit{r}}_{\mathrm{p}}={\left[r_{\mathrm{p}1}, r_{\mathrm{p}2}, r_{\mathrm{p}3}\right]}^{ \mathrm{T}}$ denote the vector from B₀ to the centre of the moving internal point mass, and ${\mathit{l}}_{\mathrm{p}}={\left[l_{\mathrm{p}1}, l_{\mathrm{p}2}, l_{\mathrm{p}3}\right]}^{ \mathrm{T}}$ represent the vector from I₀ to the centre of the moving internal point mass.

If $r_{\mathrm{p}1}{}^0$ and $\Delta r_{\mathrm{p}1}$ are the initial position and displacement of the m_p along the b₁ axis, respectively, we know that

$$r_{\mathrm{p}1}=r_{\mathrm{p}1}{}^0+\Delta r_{\mathrm{p}1}$$

(5)

Let $\xi (-\pi /2\le \xi \le \pi /2)$ denote the rotation angle of battery pack, and then the position of m_p in the body frame can be expressed as

$$\{\begin{array}{l}{\mathit{r}}_{\mathrm{p}}={\mathit{R}}_{\mathrm{BP}}{\overline{\mathit{r}}}_{\mathrm{p}}\\ {\mathit{R}}_{\mathrm{BP}}=\left[\begin{array}{ccc}1& 0& 0\\ 0& -\sin \xi & -\cos \xi \\ 0& \cos \xi & -\sin \xi \end{array}\right]\\ {\overline{\mathit{r}}}_{\mathrm{p}}={\left[r_{\mathrm{p}1}, R_{\mathrm{p}}, 0\right]}^{\mathrm{T}}\end{array}$$

(6)

2.2.3. Variable Ballast Point Mass

The variable ballast point mass m_b is the net buoyancy of the glider, which can be adjusted by the buoyancy engine with the ability to transport the hydraulic oil between an external bladder and an internal bladder (tank). Assuming that the buoyancy engine is symmetrically arranged around the b₁ axis, and the position change of the hydraulic oil can be negligible, we define ${\mathit{r}}_{\mathrm{b}}={\left[r_{\mathrm{b}1}, 0, 0\right]}^{ \mathrm{T}}$ as the vector from B₀ to the centre of the variable ballast point mass, and ${\mathit{l}}_{\mathrm{b}}={\left[l_{\mathrm{b}1}, l_{\mathrm{b}2}, l_{\mathrm{b}3}\right]}^{ \mathrm{T}}$ as the vector from I₀ to the centre of the variable ballast point mass.

Then, the following equations can be obtained:

$$m_{\mathrm{v}}=m_{\mathrm{s}}+m_{\mathrm{p}}+m_{\mathrm{b}}$$

(7)

$$\Delta m=m_{\mathrm{v}}-\overline{m}$$

(8)

where m_v is the total mass of the glider, $\overline{m}$ is the mass of water displaced by the glider, and $\Delta m$ is the net mass of the glider.

Note that before deployment, the glider should be in suspension (completely submerged in seawater), when m_b is set to zero. According to Equation (7), we have $\overline{m}=m_{\mathrm{s}}+m_{\mathrm{p}}$ on the water surface. Substituting this to Equation (8), we know that the net mass of the glider is determined by the variable ballast point mass m_b, i.e., $\Delta m=m_{\mathrm{b}}$.

2.3. Kinematics

It takes the glider several hours to complete one work cycle to cover a gliding range of several kilometers, during which the ocean currents are less likely to change drastically. Therefore, DAC is taken as the ocean current in one work cycle. In the inertial frame, the DAC vector is provided as North/ East/ Down vector components at each direction (i₁, i₂, i₃) in the ocean, which can be expressed as

$${\mathit{V}}_{\mathrm{f}}={\left[V_{\mathrm{f}1}, V_{\mathrm{f}2}, V_{\mathrm{f}3}\right]}^{ \mathrm{T}}$$

(9)

Then, the DAC in the body frame can be obtained using Equation (2).

$${\mathit{v}}_{\mathrm{f}}={\left[u_{\mathrm{f}}, v_{\mathrm{f}}, w_{\mathrm{f}}\right]}^{ \mathrm{T}}={\mathit{R}}_{\mathrm{IB}}{\mathit{V}}_{\mathrm{f}}={\mathit{R}}_{\mathrm{BI}}{}^{\mathrm{T}}{\mathit{V}}_{\mathrm{f}}$$

(10)

The linear velocity and rotation velocity relative to the current (flow) in the body frame are

$${\mathit{v}}_{\mathrm{r}}=\mathit{v}-{\mathit{v}}_{\mathrm{f}}={\left[u_{\mathrm{r}}, v_{\mathrm{r}}, w_{\mathrm{r}}\right]}^{\mathrm{T}}$$

(11)

$${\mathit{\omega }}_{\mathrm{r}}=\mathit{\omega }-{\mathit{\omega }}_{\mathrm{f}}={\left[p_{\mathrm{r}}, q_{\mathrm{r}}, r_{\mathrm{r}}\right]}^{\mathrm{T}}$$

(12)

where ${\mathit{\omega }}_{\mathrm{f}}$ is the rotation velocity of the current in the body frame.

In this study, we assume that the ocean current is irrotational, indicating that ${\mathit{\omega }}_{\mathrm{f}}={\left[0, 0, 0\right]}^{\mathrm{T}}$ and ${\mathit{\omega }}_{\mathrm{r}}=\mathit{\omega }$. Then, the kinematics equations of a glider in the ocean currents are

$$\dot{\mathit{X}}={\mathit{R}}_{\mathrm{BI}}\mathit{v}$$

(13)

$${\dot{\mathit{R}}}_{\mathrm{BI}}={\mathit{R}}_{\mathrm{BI}}\widehat{\mathit{\omega }}$$

(14)

For vectors $\mathit{a}={\left[a_1, a_2, a_3\right]}^{\mathrm{T}}$ and $\mathit{b}={\left[b_1, b_2, b_3\right]}^{\mathrm{T}}$, the rule for the symbol $\widehat{·}$ is

$$\widehat{\mathit{a}}\mathit{b}=\mathit{a}\times \mathit{b}=\left[\begin{array}{ccc}0& -a_3& a_2\\ a_3& 0& -a_1\\ -a_2& a_1& 0\end{array}\right]\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\end{array}\right]$$

2.4. Dynamics

For simplicity, this section only studies the dynamics of the stationary mass (rigid body part), whose dynamic behaviour can represent the macroscopic behaviour of the whole glider in the ocean. The forces of the moving internal point mass and variable ballast point mass are regarded as external forces.

We first define the generalized velocity of the glider in the body frame as

$${\mathit{v}}_{\mathrm{g}}=\left[\begin{array}{c}\mathit{v}\\ \mathit{\omega }\end{array}\right]={\left[u, v, w, p, q, r\right]}^{\mathrm{T}}$$

(15)

The generalized velocity of the DAC in the body frame is

$${\mathit{v}}_{\mathrm{f}}^{\mathrm{g}}=\left[\begin{array}{c}{\mathit{v}}_{\mathrm{f}}\\ \mathbf{0}\end{array}\right]={\left[u_{\mathrm{f}}, v_{\mathrm{f}}, w_{\mathrm{f}}, 0, 0, 0\right]}^{\mathrm{T}}$$

(16)

Then, the generalized velocity relative to the current (flow) in the body frame is

$${\mathit{v}}_{\mathrm{r}}^{\mathrm{g}}={\mathit{v}}_{\mathrm{g}}-{\mathit{v}}_{\mathrm{f}}^{\mathrm{g}}=\left[\begin{array}{c}{\mathit{v}}_{\mathrm{r}}\\ \mathit{\omega }\end{array}\right]$$

(17)

The kinetic energy of the combined fluid and rigid body system can be expressed as [20,22]

$$T=\frac{1}{2}{\left({\mathit{v}}_r^g\right)}^T({\mathbb{M}}_f+\overline{\mathbb{M}}){\mathit{v}}_r^g+\frac{1}{2}{\mathit{v}}_g^T{\mathbb{M}}_s{\mathit{v}}_g-\frac{1}{2}{\mathit{v}}_g^T\overline{\mathbb{M}}{\mathit{v}}_g$$

(18)

where ${\mathbb{M}}_f$ is the generalized added inertia matrix [23], $\overline{\mathbb{M}}$ is the generalized inertia matrix of the fluid replaced by the glider, ${\mathbb{M}}_s$ is the generalized inertia matrix of the rigid body of the glider. Their specific expressions are

$${\mathbb{M}}_f=\left[\begin{array}{cc}{\mathit{M}}_f& {\mathit{C}}_f^T\\ {\mathit{C}}_f& {\mathit{J}}_f\end{array}\right]$$

(19)

where ${\mathit{M}}_f$ is the added mass matrix, ${\mathit{J}}_f$ is the added inertia matrix, and ${\mathit{C}}_f$ is the cross term.

$$\overline{\mathbb{M}}=\left[\begin{array}{cc}\overline{m}\mathit{I}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}\end{array}\right]$$

(20)

where I is a 3 × 3 identity matrix, and 0 is 3 × 3 matrix of zero.

$${\mathbb{M}}_s=\left[\begin{array}{cc}{m}_s\mathit{I}& -m_s\widehat{{\mathit{r}}_s}\\ m_s\widehat{{\mathit{r}}_s}& {\mathit{J}}_s-m_s\widehat{{\mathit{r}}_s}\widehat{{\mathit{r}}_s}\end{array}\right]$$

(21)

where ${\mathit{J}}_{\mathrm{s}}$ is the inertia matrix of m_s.

Let

$$\mathit{ℒ}\left(\mathit{q}, \dot{\mathit{q}} , t\right)=T\left(\mathit{q}, \dot{\mathit{q}} , t\right)-V\left(\mathit{q}\right)$$

(22)

be the Lagrangian function for a mechanical system with generalized coordinate q. The Lagrange’s equation is

$$\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial ℒ}{\partial \dot{\mathit{q}} }\right)-\frac{\partial ℒ}{\partial \mathit{q}}=\mathit{Q}$$

(23)

where Q denotes the generalized exogenous forces.

For a UG, the conventional generalized coordinates q can be expressed as

$$\mathit{q}=\left[\begin{array}{c}\mathit{X}\\ \mathbf{\Theta }\end{array}\right]$$

(24)

Here, we use the alternative variables $\left(\mathit{q}, {\mathit{q}}_{\mathrm{w}}\right)$ rather than state elements $\left(\mathit{q}, \dot{\mathit{q}} \right)$ to express Lagrangian [18,20,22]. If ${\mathit{q}}_{\mathrm{w}}={\mathit{v}}_{\mathrm{g}}$, we have

$$\dot{\mathit{q}}=\mathit{D}\left(\mathit{q}\right){\mathit{q}}_{\mathrm{w}}$$

(25)

where $\mathit{D}\left(\mathit{q}\right)=\mathrm{diag}\left[{\mathit{R}}_{\mathrm{BI}}, {\mathit{\Omega }}_{\mathrm{BI}}\right]$.

Then, a new Lagrangian can be expressed as

$$\overline{ℒ}\left(\mathit{q}, {\mathit{q}}_{\mathrm{w}}, t\right)=ℒ\left(\mathit{q}, \mathit{D}\left(\mathit{q}\right){\mathit{q}}_{\mathrm{w}}, t\right)$$

(26)

Equation (23) can be rewritten as

$$\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial \overline{ℒ}}{\partial {\mathit{q}}_{\mathrm{w}}}\right)+G\left(\mathit{q}, {\mathit{q}}_{\mathrm{w}}\right)\frac{\partial \overline{ℒ}}{\partial {\mathit{q}}_{\mathrm{w}}}-\mathit{D}{\left(\mathit{q}\right)}^{\mathrm{T}}\frac{\partial \overline{ℒ}}{\partial \mathit{q}}=\mathit{D}{\left(\mathit{q}\right)}^{\mathrm{T}}\mathit{Q}$$

(27)

where the calculation of the elements of the matrix $G\left(\mathit{q}, {\mathit{q}}_{\mathrm{w}}\right)$ can be found in [20], and the final result is

$$G\left(\mathit{q}, {\mathit{q}}_{\mathrm{w}}\right)=\left[\begin{array}{cc}\widehat{\mathit{\omega } }& 0\\ \widehat{\mathit{v} }& \widehat{\mathit{\omega } }\end{array}\right]$$

(28)

According to Equation (18), the new Lagrangian of the combined fluid and rigid body system can be obtained by

$$\overline{\mathit{ℒ}}\left(\mathit{q}, {\mathit{q}}_w, t\right)=\frac{1}{2}{\left({\mathit{v}}_g-{\mathit{v}}_f^g\right)}^T\left({\mathbb{M}}_f+\overline{\mathbb{M}}\right)\left({\mathit{v}}_g-{\mathit{v}}_f^g\right)+\frac{1}{2}{\mathit{v}}_g^T{\mathbb{M}}_s{\mathit{v}}_g-\frac{1}{2}{\mathit{v}}_g^T\overline{\mathbb{M}}{\mathit{v}}_g$$

(29)

Substituting Equation (29) into Equation (27) and rearranging the terms, we can obtain the dynamic equations of the glider in the body frame.

$$\begin{gathered} \left({\mathbb{M}}_f+{\mathbb{M}}_s\right){\dot{\mathit{v}}}_g= \\ -\left[\begin{array}{cc}\hat{\mathit{\omega }}& \mathbf{0}\\ \hat{\mathit{v}} -{\hat{\mathit{v}}}_{f }& \hat{\mathit{\omega }}\end{array}\right]\left({\mathbb{M}}_f+\overline{\mathbb{M}}\right)\left({\mathit{v}}_g-{\mathit{v}}_f^g\right)+\left[\begin{array}{c}\mathit{f}\\ \mathit{\tau }\end{array}\right] \\ -\left[\begin{array}{cc}\hat{\mathit{\omega }}& \mathbf{0}\\ \hat{\mathit{v}}& \hat{\mathit{\omega }}\end{array}\right]\left({\mathbb{M}}_s-\overline{\mathbb{M}}\right){\mathit{v}}_g+\left({\mathbb{M}}_f+\overline{\mathbb{M}}\right)\left[\begin{array}{c}{\mathit{v}}_f\times \mathit{\omega }\\ \mathbf{0}\end{array}\right] \end{gathered}$$

(30)

where $\mathit{f}$ and τ are the external force and moment acting on the rigid body of the glider in the body frame.

To obtain the flow-relative dynamic equations, subtracting

$$\left({\mathbb{M}}_f+{\mathbb{M}}_s\right){\dot{\mathit{v}}}_f^g =\left({\mathbb{M}}_f+{\mathbb{M}}_s\right)\left[\begin{array}{c}{\mathit{v}}_f\times \mathit{\omega }\\ \mathbf{0}\end{array}\right]$$

(31)

from Equation (30) gives

$$\begin{gathered} \left({\mathbb{M}}_f+{\mathbb{M}}_s\right){\dot{\mathit{v}}}_r^g = \\ -\left[\begin{array}{cc}\hat{\mathit{\omega }}& \mathbf{0}\\ {\hat{\mathit{v}}}_r& \hat{\mathit{\omega }}\end{array}\right]\left({\mathbb{M}}_f+{\mathbb{M}}_s\right){\mathit{v}}_r^g+\left[\begin{array}{c}\mathit{f}\\ \mathit{\tau }\end{array}\right] \\ -\left[\begin{array}{cc}\hat{\mathit{\omega }}& \mathbf{0}\\ {\hat{\mathit{v}}}_r+{\hat{\mathit{v}}}_f& \hat{\mathit{\omega }}\end{array}\right]\left({\mathbb{M}}_s-\overline{\mathbb{M}}\right){\mathit{v}}_f^g-\left[\begin{array}{cc}\mathbf{0}& \mathbf{0}\\ {\hat{\mathit{v}}}_f & \mathbf{0}\end{array}\right]\left({\mathbb{M}}_s-\overline{\mathbb{M}}\right){\mathit{v}}_r^g \\ +\left(\overline{\mathbb{M}}-{\mathbb{M}}_s\right)\left[\begin{array}{c}{\mathit{v}}_f\times \mathit{\omega }\\ \mathbf{0}\end{array}\right] \end{gathered}$$

(32)

2.5. Forces and Moments Acting on the Rigid Body of the Glider

External forces acting on the rigid body include net buoyancy and hydrodynamic force, and external moments acting on the rigid body include the moments from the stationary mass, moving internal point mass and variable ballast point mass, and the hydrodynamic moments.

At depth z, the net mass of a glider in the inertial frame is

$$\Delta m\left(z\right)=m_{\mathrm{b}}+\overline{m}\left(0\right)-\overline{m}\left(z\right)$$

(33)

According to Section 2.2.3, we know that

$$\overline{m}\left(0\right)=m_{\mathrm{s}}+m_{\mathrm{p}}$$

(34)

The mass of water displaced by the glider at depth z can be calculated by

$$\overline{m}\left(z\right)=\rho \left(z\right)V_{\mathrm{h}}\left(z\right)$$

(35)

where $\rho \left(z\right)$ is the seawater density at depth z, which can be fitted by a polynomial [1]

$$\rho \left(z\right)=p_1{z}^4+p_2{z}^3+p_3{z}^2+p_{4}z+p_5$$

(36)

where $V_{\mathrm{h}}\left(z\right)$ is the volume of the glider at depth z, which can be expressed as

$$V_{\mathrm{h}}\left(z\right)=V_{\mathrm{h}0}-\Delta V_{\mathrm{h}}\left(z\right)$$

(37)

where $V_{\mathrm{h}0}$ and $\Delta V_{\mathrm{h}}\left(z\right)$ are the volume of the glider on the water surface and the volume reduction of the glider at depth z caused by the water pressure, respectively. They can be calculated by

$$V_{\mathrm{h}0}=\frac{\overline{m}\left(0\right)}{\rho \left(0\right)}$$

(38)

$$\Delta {\mathit{V}}_{\mathrm{h}}\left(z\right)=K_{\mathrm{vh}}z$$

(39)

where $K_{\mathrm{vh}}$ is the compressibility of the pressure of the glider.

Substituting Equations (34)–(39) into Equation (33) results in

$$\Delta m\left(z\right)=m_{\mathrm{b}}+m_{\mathrm{s}}+m_{\mathrm{p}}-\rho \left(z\right)\left(\frac{m_{\mathrm{s}}+m_{\mathrm{p}}}{p_5}-K_{\mathrm{vh}}z\right)$$

(40)

Then, the net buoyancy in the body frame can be represented as follows

$${\mathit{f}}_{\mathrm{nb}}={\mathit{R}}_{\mathrm{BI}}^{\mathrm{T}}\left[\begin{array}{c}0\\ 0\\ \Delta m\left(z\right)\mathrm{g}\end{array}\right]$$

(41)

Let ${\mathit{\tau }}_{\mathrm{s}}$, ${\mathit{\tau }}_{\mathrm{p}}$ and ${\mathit{\tau }}_{\mathrm{b}}$ denote the moments imposed by the stationary mass, moving internal point mass and variable ballast point mass, respectively. Then, we have

$${\mathit{\tau }}_{\mathrm{s}}={\mathit{r}}_{\mathrm{s}}\times {\mathit{R}}_{\mathrm{BI}}^{\mathrm{T}}{\left[0, 0, m_{\mathrm{s}}\mathrm{g}\right]}^{\mathrm{T}}$$

(42)

$${\mathit{\tau }}_{\mathrm{p}}={\mathit{r}}_{\mathrm{p}}\times {\mathit{R}}_{\mathrm{BI}}^{\mathrm{T}}{\left[0, 0, m_{\mathrm{p}}\mathrm{g}\right]}^{\mathrm{T}}$$

(43)

$${\mathit{\tau }}_{\mathrm{b}}={\mathit{r}}_{\mathrm{b}}\times {\mathit{R}}_{\mathrm{BI}}^{\mathrm{T}}{\left[0, 0, m_{\mathrm{b}}\mathrm{g}\right]}^{\mathrm{T}}$$

(44)

For a glider with two symmetry planes of ${\mathit{b}}_1{\mathrm{B}}_0{\mathit{b}}_3$ and ${\mathit{b}}_1{\mathrm{B}}_0{\mathit{b}}_2$, the inertial hydrodynamic force ${\mathit{f}}_{\mathrm{I}}$ and moment ${\mathit{\tau }}_{\mathrm{I}}$ can be expressed as

$$\left[\begin{array}{c}{\mathit{f}}_{\mathrm{I}}\\ {\mathit{\tau }}_{\mathrm{I}}\end{array}\right]={\mathbb{M}}_f{[{\dot{\mathit{v}}}_r\hspace{1em}{\dot{\mathit{\omega }}}_r]}^T=\left[\begin{array}{c}{f}_{I1}\\ f_{I2}\\ f_{I3}\\ {\tau }_{I1}\\ {\tau }_{I2}\\ {\tau }_{I3}\end{array}\right]=\left[\begin{array}{cccccc}{\lambda }_{11}& 0& 0& 0& 0& 0\\ 0& {\lambda }_{22}& 0& 0& 0& {\lambda }_{26}\\ 0& 0& {\lambda }_{33}& 0& {\lambda }_{35}& 0\\ 0& 0& 0& {\lambda }_{44}& 0& 0\\ 0& 0& {\lambda }_{53}& 0& {\lambda }_{55}& 0\\ 0& {\lambda }_{62}& 0& 0& 0& {\lambda }_{66}\end{array}\right]\left[\begin{array}{c}{\dot{u}}_r\\ {\dot{v}}_r\\ {\dot{w}}_r\\ {\dot{p}}_r\\ {\dot{q}}_r\\ {\dot{r}}_r\end{array}\right]$$

(45)

where ${\lambda }_{ij}$ is the inertial hydrodynamic coefficient.

In the flow frame, the viscous hydrodynamic force ${\mathit{f}}_{{}_{\mathrm{S}}}^V$ and moment ${\mathit{\tau }}_{{}_{\mathrm{S}}}^{{V}}$ are usually expressed as

$$\left[\begin{array}{c}{\mathit{f}}_{{}_{\mathrm{S}}}^{\mathrm{V}}\\ {\mathit{\tau }}_{{}_{\mathrm{S}}}^{\mathrm{V}}\end{array}\right]=\left[\begin{array}{c}-\mathit{D}\\ SF\\ -L\\ T_{\mathrm{DL}1}\\ T_{\mathrm{DL}2}\\ T_{\mathrm{DL}3}\end{array}\right]=V_{\mathrm{r}}^2\left[\begin{array}{c}-K_{\mathrm{D}0}-K_{\mathrm{D}}{\alpha }^2\\ K_{\beta }\beta \\ -K_{\mathrm{L}0}-K_{\mathrm{L}}\alpha \\ K_{\mathrm{MR}}\beta +K_{\mathrm{p}}{p}_{\mathrm{r}}\\ K_{\mathrm{M}0}+K_{\mathrm{M}}\alpha +K_{\mathit{q}}{\mathit{q}}_{\mathrm{r}}\\ K_{\mathrm{MY}}\beta +K_{\mathrm{r}}{r}_{\mathrm{r}}\end{array}\right]$$

(46)

where D, SF and L are the drag, side force, and lift, respectively, T_DL1, T_DL2, T_DL3 are the hydrodynamic moments about the axes of the flow frame c₁, c₂ and c₃, respectively, K_* and K_** are coefficients of viscous hydrodynamic forces and moments obtained by the computational fluid dynamic (CFD) method [24], and V_r is the value of the velocity relative to the current, which can be calculated by

$$V_{\mathrm{r}}=\Vert {\mathit{v}}_{\mathrm{r}}\Vert =\sqrt{u_{\mathrm{r}}{}^2+{\mathit{v}}_{\mathrm{r}}{}^2+w_{\mathrm{r}}{}^2}$$

(47)

The attack angle α and slip angle β can be expressed as

$$\alpha {=\tan }^{-1}\left(\frac{w_{\mathrm{r}}}{u_{\mathrm{r}}}\right)$$

(48)

$$\beta ={\sin }^{-1}\left(\frac{v_{\mathrm{r}}}{V_{\mathrm{r}}}\right)$$

(49)

Then, we map the viscous hydrodynamic forces and moments to the body frame as

$$\left[\begin{array}{c}{\mathit{f}}_{\mathrm{S}}\\ {\mathit{\tau }}_{\mathrm{S}}\end{array}\right]={\mathit{R}}_{\mathrm{VB}}\left[\begin{array}{c}{\mathit{f}}_{{}_{\mathrm{S}}}^{\mathrm{V}}\\ {\mathit{\tau }}_{{}_{\mathrm{S}}}^{{V}}\end{array}\right]$$

(50)

where ${\mathit{f}}_{\mathrm{S}}$ and ${\mathit{\tau }}_{\mathrm{S}}$ are the viscous hydrodynamic forces and moments in the body frame, respectively.

Finally, we can obtain the total external forces and moments by

$$\left[\begin{array}{c}\mathit{f}\\ \mathit{\tau }\end{array}\right]=\left[\begin{array}{c}{\mathit{f}}_{\mathrm{nb}}+{\mathit{f}}_{\mathrm{I}}+{\mathit{f}}_{\mathrm{S}}\\ {\mathit{\tau }}_{\mathrm{s}}+{\mathit{\tau }}_{\mathrm{p}}+{\mathit{\tau }}_{\mathrm{b}}+{\mathit{\tau }}_{\mathrm{I}}+{\mathit{\tau }}_{\mathrm{S}}\end{array}\right]$$

(51)

3. Energy Consumption Model

3.1. Steady-State Motion

For a glider operating in DAC based on open-loop control, the following vectors can be used to describe the steady-state motion of the glider at a certain depth z in the body frame.

$${\mathit{v}}_{\mathrm{r}}={\left[u_{\mathrm{r}}, 0, w_{\mathrm{r}}\right]}^{\mathrm{T}}$$

(52)

$${\mathit{\omega }}_{\mathrm{r}}=\mathbf{0}$$

(53)

$$\Theta ={\left[0, \theta , \psi \right]}^{\mathrm{T}}$$

(54)

$${\dot{\mathit{v}}}_{\mathrm{r}}^{\mathrm{g} }=0$$

(55)

$${\mathit{r}}_{\mathrm{p}}={\left[r_{\mathrm{p}1}, 0, R_{\mathrm{p}}\right]}^{\mathrm{T}}$$

(56)

Substituting Equations (52)–(56) into the dynamic model Equation (32) results in the steady-state dynamic equations of the glider at depth z as

$$0=-D\cos \alpha +L\sin \alpha -\Delta m\left(z\right)\mathrm{gsin}\theta$$

(57)

$$0=-D\sin \alpha -L\cos \alpha +\Delta m\left(z\right)\mathrm{gcos}\theta$$

(58)

$$\begin{array}{l}0=T_{\mathrm{DL}2}-m_{\mathrm{s}}\mathrm{g}\left(r_{\mathrm{s}1}\cos \theta +r_{\mathrm{s}3}\sin \theta \right)-m_{\mathrm{b}}\mathrm{g}{r}_{\mathrm{b}1}\cos \theta \\ -m_{\mathrm{p}}\mathrm{g}\left(r_{\mathrm{p}1}\cos \theta +R_{\mathrm{p}}\sin \theta \right)+\left({\lambda }_{33}-{\lambda }_{11}\right)u_{\mathrm{r}}{w}_{\mathrm{r}}\end{array}$$

(59)

Let γ denote the gliding angle, which represents the angle between the projection of v onto ${\mathit{b}}_1{\mathrm{B}}_0{\mathit{b}}_3$ and the horizontal plane of the inertial frame, as

$$\gamma =\theta -\alpha$$

(60)

Then, Equations (57) and (58) can be rewritten as [25]

$$\begin{gathered} \left[\begin{array}{c}0\\ \Delta m\left(z\right)\mathrm{g}\end{array}\right]=\left[\begin{array}{cc}\cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{array}\right]\left[\begin{array}{cc}\cos \alpha & -\sin \alpha \\ \sin \alpha & \cos \alpha \end{array}\right]\left[\begin{array}{c}L\\ D\end{array}\right] \\ =\left[\begin{array}{cc}\cos \gamma & \sin \gamma \\ -\sin \gamma & \cos \gamma \end{array}\right]\left[\begin{array}{c}{K}_{\mathrm{D}0}+K_{\mathrm{D}}{\alpha }^2\\ K_{\mathrm{L}0}+K_{\mathrm{L}}\alpha \end{array}\right]{\mathit{V}}_{\mathrm{r}}^2 \end{gathered}$$

(61)

Taking the first line of Equation (61), together with the gliding velocity ${\mathit{V}}_{\mathrm{r}}\ne 0$ and gliding angle $\gamma \ne 0$, we have

$${\alpha }^2+\frac{K_{\mathrm{L}}}{K_{\mathrm{D}}}\tan \gamma \cdot \alpha +\frac{1}{K_{\mathrm{D}}}\left(K_{\mathrm{D}0}+K_{\mathrm{L}0}\tan \gamma \right)=0$$

(62)

To solve α in Equation (62), we have

$$\alpha \left(\gamma \right)=\frac{1}{2}\frac{K_{\mathrm{L}}}{K_{\mathrm{D}}}\tan \gamma \left(-1+\sqrt{1-4\frac{K_{\mathrm{D}}}{K_{\mathrm{L}}^2}\cot \gamma \left(K_{\mathrm{D}0}\cot \gamma +K_{\mathrm{L}0}\right)}\right)$$

(63)

Then, velocities $u_{\mathrm{r}}$ and $w_{\mathrm{r}}$ can be expressed as

$$u_{\mathrm{r}}={\mathit{V}}_{\mathrm{r}}\cos \alpha$$

(64)

$$w_{\mathrm{r}}={\mathit{V}}_{\mathrm{r}}\sin \alpha$$

(65)

3.2. Energy Consumption of the Buoyancy Engine

To avoid confusion, hereinafter, subscripts 1 and 2 are used to represent the descending and ascending of the glider, respectively.

In one work cycle, the solenoid valve will open on the water surface to allow the hydraulic oil in the external bladder to flow into the internal tank, and the pump will work at the inflection point under the water pressure to transfer the hydraulic oil from internal tank to the external bladder.

According to the second line of Equation (61), if the glider reaches a given velocity ${\mathit{V}}_{\mathrm{r}}$ and a gliding angle γ at depth z, the net mass supplied by the buoyancy engine should be

$$\Delta m\left(z\right)=\left(-\sin \gamma \left(K_{\mathrm{D}0}+K_{\mathrm{D}}{\alpha }^2\right)+\cos \gamma \left(K_{\mathrm{L}0}+K_{\mathrm{L}}\alpha \right)\right){\mathit{V}}_{\mathrm{r}}{}^2/\mathrm{g}$$

(66)

Substituting Equation (40) into Equation (66) and rearranging the terms, we have

$$\begin{array}{l}{m}_{\mathrm{b}}=\frac{\left(-\sin \gamma \left(K_{\mathrm{D}0}+K_{\mathrm{D}}{\alpha }^2\right)+\cos \gamma \left(K_{\mathrm{L}0}+K_{\mathrm{L}}\alpha \right)\right){\mathit{V}}_{\mathrm{r}}{}^2}{\mathrm{g}}\\ +\rho \left(z\right)\left(\frac{m_{\mathrm{s}}+m_{\mathrm{p}}}{p_5}-K_{\mathrm{vh}}z\right)-\left(m_{\mathrm{s}}+m_{\mathrm{p}}\right)\end{array}$$

(67)

subject to

$$\left|m_{\mathrm{b}}\right|\le m_{\mathrm{b} \max }$$

The volume of the transported oil in one work cycle is

$$\Delta {\mathit{V}}_{\mathrm{oil}}=\left|\frac{m_{\mathrm{b}1}}{\rho \left(z_1\right)}\right|+\left|\frac{m_{\mathrm{b}2}}{\rho \left(z_2\right)}\right|$$

(68)

where $z_1$ and $z_2$ are the final depth during descending and ascending, respectively. In general, $z_2$ is zero, excluding for some special gliding models, such as the hidden gliding strategy model [13]. In this study, $z_2$ is set to zero, and then $z_1$ can be replaced by z.

The energy consumption of the buoyancy engine in one work cycle can be then calculated by

$$E_{\mathrm{b}}=\Delta {\mathit{V}}_{\mathrm{oil}}\left(\frac{P_{\mathrm{pump}}\left(z\right)}{{\mathit{q}}_{\mathrm{pump}}\left(z\right)}+\frac{P_{\mathrm{v}}}{{\mathit{q}}_{\mathrm{v}}\left(0\right)}\right)$$

(69)

where $P_{\mathrm{pump}}\left(z\right)$ and ${\mathit{q}}_{\mathrm{pump}}\left(z\right)$ are the power and flow rate of the pump at depth z, respectively, $P_{\mathrm{v}}$ and ${\mathit{q}}_{\mathrm{v}}\left(0\right)$ are the power and flow rate of the solenoid on the water surface, respectively, which can be found in [14].

Substituting Equations (67) and (68) into Equation (69) results in

$$\begin{array}{l}{E}_{\mathrm{b}}=\left(\begin{array}{l}\left|\frac{\left(-\sin {\gamma }_1\left(K_{\mathrm{D}0}+K_{\mathrm{D}}{\alpha }_1{}^2\right)+\cos {\gamma }_1\left(K_{\mathrm{L}0}+K_{\mathrm{L}}{\alpha }_1\right)\right)}{\rho \left(z\right)\mathrm{g}}{\mathit{V}}_{\mathrm{r}1}^2+\frac{m_{\mathrm{s}}+m_{\mathrm{p}}}{p_5}-\frac{m_{\mathrm{s}}+m_{\mathrm{p}}}{\rho \left(z\right)}-K_{\mathrm{vh}}z\right|\\ +\left|\frac{\left(-\sin {\gamma }_2\left(K_{\mathrm{D}0}+K_{\mathrm{D}}{\alpha }_2{}^2\right)+\cos {\gamma }_2\left(K_{\mathrm{L}0}+K_{\mathrm{L}}{\alpha }_2\right)\right)}{\rho \left(0\right)\mathrm{g}}{\mathit{V}}_{\mathrm{r}2}^2+\frac{m_{\mathrm{s}}+m_{\mathrm{p}}}{p_5}-\frac{m_{\mathrm{s}}+m_{\mathrm{p}}}{\rho \left(0\right)}\right|\end{array}\right)\\ \times \left(\frac{P_{\mathrm{pump}}\left(z\right)}{{\mathit{q}}_{\mathrm{pump}}\left(z\right)}+\frac{P_{\mathrm{v}}}{{\mathit{q}}_{\mathrm{v}}\left(0\right)}\right)\end{array}$$

(70)

3.3. Energy Consumption of the Attitude Regulation Unit

Before deployment, the initial position of the moving internal point mass in the body frame is

$$r_{\mathrm{p}1}^0=\frac{-m_{\mathrm{s}}{r}_{\mathrm{s}1}}{m_{\mathrm{p}}}$$

(71)

According to Equation (59), if the glider reaches a given velocity ${\mathit{v}}_{\mathrm{r}}$ and a gliding angle γ at depth z, the position of the moving internal point mass is

$$\begin{array}{l}{r}_{\mathrm{p}1}=\frac{T_{\mathrm{DL}2}-m_{\mathrm{s}}\mathrm{g}\left(r_{\mathrm{s}1}\cos \left(\alpha +\gamma \right)+r_{\mathrm{s}3}\sin \left(\alpha +\gamma \right)\right)-m_{\mathrm{b}}\mathrm{g}{r}_{\mathrm{b}1}\cos \left(\alpha +\gamma \right)+\left({\lambda }_{33}-{\lambda }_{11}\right)u_{\mathrm{r}}{w}_{\mathrm{r}}}{m_{\mathrm{p}}\mathrm{gcos}\left(\alpha +\gamma \right)}\\ -R_{\mathrm{p}}\tan \left(\alpha +\gamma \right)\end{array}$$

(72)

subject to

$$r_{\mathrm{p}1 \min }\le r_{\mathrm{p}1}\le r_{\mathrm{p}1 \max }$$

The displacement of the moving internal point mass is

$$\Delta r_{\mathrm{p}1}=\left|r_{\mathrm{p}1}-r_{\mathrm{p}1}^0\right|$$

(73)

As shown in Figure 4, for a glider based on open-loop control, the total displacement of the moving internal point mass in one work cycle is $2\left(\Delta r_{\mathrm{p}11}+\Delta r_{\mathrm{p}12}\right)$.

Then, the energy consumed by the attitude regulation unit can be expressed as

$$E_{\mathrm{p}}=\frac{2\Delta r_{\mathrm{p}11}}{{\mathit{v}}_{\mathrm{p}}}{P}_{\mathrm{p}}\left(\left|{\gamma }_1\right|\right)+\frac{2\Delta r_{\mathrm{p}12}}{{\mathit{v}}_{\mathrm{p}}}{P}_{\mathrm{p}}\left(\left|{\gamma }_2\right|\right)$$

(74)

where $\Delta r_{\mathrm{p}11}$ and $\Delta r_{\mathrm{p}12}$ are the displacement of the moving internal point mass during descending and ascending, respectively, ${\mathit{v}}_{\mathrm{p}}$ is the velocity of the moving internal point mass along the ${\mathit{b}}_1$ axis in the body frame, and $P_{\mathrm{p}}\left(\left|\gamma \right|\right)$ is the input power of the motor.

3.4. Energy Consumption of the Control Unit

Since the control unit keeps working throughout the operation of the glider, the energy it consumed is determined by the elapsed time of the glider. To be more accurate and practical, we divide the elapsed time into navigation time t, which represents the time taken to complete the ascent and descent movements, and additional floating time t_add, during which the glider floats on the surface to perform instructions and data transfers. In order to reduce the number of variables, hereinafter, the average velocity of the glider is replaced by V_r.

For a glider with the diving depth z, the elapsed time of the control unit can be expressed by the following equations:

$$t_{\mathrm{c}}=t+t_{\mathrm{add}}$$

(75)

$$t=\frac{z}{\left|V_{\mathrm{r}1}\sin {\gamma }_1+V_{\mathrm{f}3}\right|}+\frac{z}{\left|V_{\mathrm{r}2}\sin {\gamma }_2+V_{\mathrm{f}3}\right|}$$

(76)

Let $P_{\mathrm{c}}$ denote the average power of the control unit, then the energy consumption in one work cycle is

$$E_{\mathrm{c}}=t_{\mathrm{c}}{P}_{\mathrm{c}}=\left(\frac{z}{\left|V_{\mathrm{r}1}\sin {\gamma }_1+V_{\mathrm{f}3}\right|}+\frac{z}{\left|V_{\mathrm{r}2}\sin {\gamma }_2+V_{\mathrm{f}3}\right|}+t_{\mathrm{add}}\right)P_{\mathrm{c}}$$

(77)

3.5. Energy Consumption of the Detection Unit

The detection unit includes mission sensors. Let $\Delta z_i$, $\Delta t_i$ and $P_{\mathrm{s}i}$ denote the single sampling time, the sampling depth interval and the power of the ith sensor, respectively, and the energy consumption of the detection unit in one cycle can be expressed as

$$E_{\mathrm{d}}={\displaystyle \sum _{i}t{P}_{\mathrm{s}i}}={\displaystyle \sum _j\left(\frac{z\Delta t_i}{\Delta z_i+\Delta t_i\left|V_{\mathrm{r}1}\sin {\gamma }_1+V_{\mathrm{f}3}\right|}+\frac{z\Delta t_i}{\Delta z_i+\Delta t_i\left|V_{\mathrm{r}2}\sin {\gamma }_2+V_{\mathrm{f}3}\right|}\right)}{P}_{\mathrm{s}i}$$

(78)

For the uninterrupted sensors which keep running throughout the whole operation of the glider, $\Delta z_i$ in Equation (78) should be set to zero.

3.6. Energy Consumption of the Communication Unit

Communication occurs only when the glider is floating on the surface. The energy consumption for communication E_com is determined by the elapsed time $t_{\mathrm{com}}$ and the average power $P_{\mathrm{com}}$ of the communication module, such as Iridium. Therefore, we have

$$E_{\mathrm{com}}=t_{\mathrm{com}}{P}_{\mathrm{com}}$$

(79)

3.7. Energy Consumption Model and Its Characteristics

Summing the energy consumption of each unit, we can obtain the energy consumption model of a UG in ocean currents based on open-loop control as

$$E=E_{\mathrm{b}}+E_{\mathrm{p}}+E_{\mathrm{c}}+E_{\mathrm{d}}+E_{\mathrm{com}}$$

(80)

Compared with the previous ECMs, this model has the following characteristics:

• The velocity relative to the current is used to reflect the effects of ocean currents;
• The BL caused by seawater density variation and deformation of the pressure hull is considered;
• The asymmetry in velocity and gliding angle during ascending and descending is taken into account.

3.8. Model Validation

A Petrel-II glider is used to validate the proposed ECM, which completed 188 profiles in the South China Sea. This glider is equipped with an uninterrupted conductivity–temperature–depth (CTD) sensor and an interrupted altimeter. The parameters and coefficients of the glider (

Table A1. Mass and geometric characteristic parameters.

Parameter	Value	Parameter	Value
$m_{\mathrm{s}}$	50 kg	${\mathit{q}}_{\mathit{v}}$(0)	7.5 × 10−6 m3/s
${\mathit{r}}_{\mathrm{s}}$	[−0.00271, 0, 0.0199]T m	${\mathit{v}}_{\mathrm{p}}$	0.001 m/s
$m_{\mathrm{p}}$	18.9 kg	$P_{\mathrm{com}}$	3 W
$R_{\mathrm{p}}$	0.016 m	$P_{\mathrm{s}1}$	0.2 W
$r_{\mathrm{b}1}$	0.9 m	$P_{\mathrm{s}2}$	1.92 W
$K_{\mathrm{vh}}$	0.27172 mL/m	$P_{\mathit{v}}$	5 W
${\mathit{r}}_{\mathrm{w}}$	[−0.5, 0, 0]T m	$P_{\mathrm{c}}$	2 W
$r_{\mathrm{p}1 \min }$	−0.016 m	$\Delta z_2$	30 m
$r_{\mathrm{p}1 \max }$	0.054 m	$\Delta t_2$	3 s
$m_{\mathrm{b} \max }$	0.7 kg

,

Table A2. Inertial hydrodynamic coefficients.

Coefficient	Value	Coefficient	Value
${\lambda }_{11}$	1.30 kg	${\lambda }_{55}$	23.56 kg·m2
${\lambda }_{22}$	60.53 kg	${\lambda }_{66}$	21.83 kg·m2
${\lambda }_{33}$	79.66 kg	${\lambda }_{26}={\lambda }_{62}$	−19.35 kg·m
${\lambda }_{44}$	0 kg·m2	${\lambda }_{35}={\lambda }_{53}$	12.37 kg·m

and

Table A3. Viscous hydrodynamic coefficients.

Coefficient	Value	Coefficient	Value
$K_{\mathrm{D}0}$	7.65 kg/m	$K_{\mathrm{p}}$	−18.98 kg·s/rad
$K_{\mathrm{D}}$	357.97 kg/m/rad2	$K_{\mathrm{M}0}$	0.3 kg
$K_{\beta }$	−100.13 kg/m/rad	$K_{\mathrm{M}}$	−61.02 kg/rad
$K_{\mathrm{L}0}$	−0.5 kg/m	$K_{\mathit{q}}$	−196.78 kg·s/rad
$K_{\mathrm{L}}$	381.73 kg/m/rad	$K_{\mathrm{MY}}$	32.63 kg/rad
$K_{\mathrm{MR}}$	−55.81 kg/rad	$K_{\mathrm{r}}$	−372.54 kg·s/rad

) and environment (

Table A4. Environmental parameters.

Parameter	Value	Parameter	Value
g	9.8 m/s2	$p_3$	−2.75 × 10−5 kg/m5
$p_1$	−5.083 × 10−12 kg/m7	$p_4$	0.02248 kg/m4
$p_2$	1.95 × 10−8 kg/m6	$p_5$	1022.7 kg/m3

) are in Appendix A. As the vertical components of the currents are generally negligible in comparison to the lateral currents, herein, we assume that V_f3 = 0 and the velocity of the DAC in the horizontal plane is calculated as shown in Figure 5.

• Calculate the velocity of the glider along i₃ axis using

  $$V_{\mathrm{act}\_\mathrm{i}3}\left(t_i\right)=\frac{z_{\mathrm{act}}\left(t_i\right)-z_{\mathrm{act}}\left(t_{i-1}\right)}{t_i-t_{i-1}},i=1,2,\dots , n$$

  (81)
• Calculate the velocity to the North and East without the currents by

  $$V_{\mathrm{cal}\_\mathrm{i}1}\left(t_i\right)=\frac{V_{\mathrm{act}\_\mathrm{i}3}\left(t_i\right)}{\left|\tan \left({\gamma }_{\mathrm{act}}\left(t_i\right)\right)\right|}\cos \left({\psi }_{\mathrm{act}}\left(t_i\right)\right),i=1,2,\dots , n$$

  (82)

  $$V_{\mathrm{cal}\_\mathrm{i}2}\left(t_i\right)=\frac{V_{\mathrm{act}\_\mathrm{i}3}\left(t_i\right)}{\left|\tan \left({\gamma }_{\mathrm{act}}\left(t_i\right)\right)\right|}\sin \left({\psi }_{\mathrm{act}}\left(t_i\right)\right),i=1,2,\dots , n$$

  (83)

  where ${\gamma }_{\mathrm{act}}$ and ${\psi }_{\mathrm{act}}$ denote the actual gliding angle and heading, respectively.
• Calculate the displacement to the North (x_cal) and East (y_cal) without the currents by

  $$x_{\mathrm{cal}}={\displaystyle \underset{0}{\overset{t_n}{\int }}{V}_{\mathrm{cal}\_\mathrm{i}1}\left(t\right)\mathrm{d}}t$$

  (84)

  $$y_{\mathrm{cal}}={\displaystyle \underset{0}{\overset{t_n}{\int }}{V}_{\mathrm{cal}\_\mathrm{i}2}\left(t\right)\mathrm{d}}t$$

  (85)
• Calculate the DAC velocity to the North and East according to the actual displacement of the glider on the surface $\left(x_{\mathrm{act}}, y_{\mathrm{act}}\right)$ obtained by the global positioning system (GPS).

$$V_{\mathrm{f}1}=\frac{x_{\mathrm{act}}-x_{\mathrm{cal}}}{t}$$

(86)

$$V_{\mathrm{f}2}=\frac{y_{\mathrm{act}}-y_{\mathrm{cal}}}{t}$$

(87)

The actual energy consumption is calculated using the voltage, current, and time recorded by the onboard controller. The energy consumed by the motor for roll regulation, which is random and difficult to calculate in the ECM, is subtracted from actual energy consumption [14]. Figure 6 shows the actual energy consumption and that calculated by the ECM without considering the DAC of the 188 profiles. Figure 7 shows the actual energy consumption and that calculated by the ECM considering the DAC of the 188 profiles.

The mean absolute percentage error (MAPE) is employed to investigate the accuracy of the ECM. The MAPE of ECM without considering DAC is 47.58%, while that of ECM considering DAC is 14.45%, showing a difference of 33.13%, which indicates that the ocean currents have significant impact on the accuracy of the ECM. From Figure 6 and Figure 7, we can find that the energy consumption calculated by the ECM is generally greater than the actual value. The main reason is that the final velocity, which is smaller than the average velocity, is used to calculate the navigation time t, leading to the greater energy consumption of the control unit. Other factors of the error mainly include seawater density error, which may mainly cause the inaccurate calculation of m_b, and other marine environmental factors, such as biofouling, variation in temperature, and biological attacks.

4. Simulations and Discussions

Simulations are performed to investigate the relationship between energy consumption and different parameters. First of all, the factors that may have the most important impact on energy consumption [14] are first analysed, including the gliding velocity relative to the current V_r, gliding angle γ, and diving depth z. In the simulations, we assume that the analysed parameters are consistent during the ascending and descending process. The range of V_r is from 0.198 m/s to 0.766 m/s, the range of $\left|\gamma \right|$ is from 15 degrees to 45 degrees, while the diving depth is set to five typical values, including 200 m, 400 m, 600 m, 800 m and 1000 m. The results are shown in Figure 8. Because of the seawater density variation and the limitations of $m_{\mathrm{b} \max }$ and r_{p1 max/min}, the range of achievable V_r decreases as the diving depth increases and absolute gliding angle decreases. The energy consumption increases with the increasing depth and decreases with the increasing absolute gliding angle, which is consistent with the facts that the power of the buoyancy engine and the navigation time of the glider will increase with diving depth, and smaller absolute gliding angles will also result in longer navigation time. With fixed $\left|\gamma \right|$ and z, the energy consumption decreases first and then increases with the growth of V_r, which is the result of a trade-off between the energy consumption of the buoyancy engine and the navigation time of the glider, and this phenomenon becomes more apparent as diving depth increases.

To illustrate the effects of the DAC, an example of a UG moving in an ocean current is considered. For simplicity, we assume that a glider is sailing at an absolute northerly velocity V_i1, a certain gliding angle and a depth in the inertial frame without a roll. There is only a current due north V_f1. Therefore, we have ψ = 0, φ = 0, V_f2 = V_f3 = 0. Then, we take V_i1 = 0.2 m/s as an example to calculate the energy consumption with different absolute gliding angles, diving depths, and due north current. To ensure that reasonable range of the gliding velocity relative to the current V_r can be covered, the range of V_f1 is set to −0.5 m/s to 0.1 m/s. The range of $\left|\gamma \right|$ is from 15 degrees to 45 degrees, while the diving depth is set to five typical values including 200 m, 400 m, 600 m, 800 m, and 1000 m. The energy consumption obtained by the ECM is shown in Figure 9, where the limits of m_b and r_p1 are taken into account, and the parts that cannot be realized in practice have been removed. The results demonstrate that at diving depths from 200 m to 1000 m, the energy consumption decreases first and then increases with the growth of V_f1, which has an important effect on the velocity relative to the current V_r. Since V_i1 is fixed, V_r will decrease as V_f1 grows. Therefore, the results are consistent with those in Figure 8. To maintain the desired V_i1, the buoyancy engine has to transfer more hydraulic fluid to increase the driving force at a higher current velocity in the direction opposite to V_i1 (negative V_f1), and then the navigation time of the glider to complete a work cycle will be reduced, resulting in higher energy consumption of the buoyancy engine and lower energy consumption of control unit and detection unit. For a higher current velocity in the same direction as V_i1 (positive V_f1), the change in energy consumption is reversed. Therefore, the result of energy consumption is a trade-off between the energy consumption of the buoyancy engine and the navigation time of the glider. The above phenomenon becomes less obvious as the depth decreases. This can be explained by the fact that for an underwater glider operating in a shallow area, the change in energy consumption caused by navigation time is smaller, and the influence of the buoyancy engine on energy consumption is more significant.

5. Conclusions

This work presents an ECM of UG based on open-loop control considering the DAC and the BL caused by seawater density variation and deformation of the pressure hull. In the ECM, the asymmetry of gliding motion during ascending and descending is also taken into account. The ECMs with and without the consideration of ocean currents are compared using sea trial data, and the results show that the accuracy of the ECM considering ocean currents is significantly higher. Several simulations are performed to reveal the relationship between energy consumption and multiple parameters. The results demonstrate that the seawater density variation, the maximum capacity of the buoyancy engine, and the moving range of the moving internal point mass will limit the gliding velocity relative to the current in the deep sea. The energy consumption increases with the increase of diving depth and the decrease of absolute gliding angle, while with the growth of gliding velocity relative to the current, the energy consumption decreases first and then increases. In addition, the effects of DAC on energy consumption are illustrated by a simple example, and the specific reasons are explained.

Future work will focus on the establishment and validation of the energy consumption model of the underwater glider based on closed-loop control in time-varying and space-varying currents, and the development of the energy consumption control strategies.