Dynamics Modeling and Hydrodynamic Coefﬁcients Identiﬁcation of the Wave Glider

: The wave glider is an ocean-wave-propelled autonomous marine vehicle with unique dual-body architecture, which can converse the energy obtained from the ocean wave into the forward thrust. In this paper, the dynamic models of the submerged glider based on dynamic characteristics of tandem hydrofoils and the surface ﬂoat were separately established. The pitching angles of the hydrofoils and the submerged glider and the angle of attack between hydrofoils and relative current were considered for dynamic models and hydrodynamic coefﬁcients. The translational hydrodynamic coefﬁcient term for high-angle-of-attack passive motion of the submerged glider was calculated from static test simulations by using Computational Fluid Dynamics (CFD). Moreover, the rotational damping coefﬁcients and added mass coefﬁcients varying with the pitching angle of hydrofoils were analyzed by the simulation of the vertical planar motion mechanism (VPMM) tests. Furthermore, the numerical simulation of longitudinal motion with the computed hydrodynamic coefﬁcients was performed, and the simulation results were compared with the sea trial data. The analysis was performed, and conclusions were drawn, which would provide a theoretical reference for the design of the wave glider.


Introduction
In the 21st century, marine resources play an important role in supporting world economic prosperity. Numerous unmanned marine platforms are developing at an incredible pace. The wave glider invented by Robert Hine [1] is an ocean-wave-propelled autonomous marine vehicle which uses wave energy for propulsion and solar energy for onboard sensors and control and communication systems, etc. In comparison with other marine platforms, the wave glider has longer endurance and lower cost and has been widely applied to various marine scientific research in recent years.
The wave glider is composed of a submerged glider, a surface float and an umbilical cable, and its energy conversion system using wave energy is based on a dual-body design. When the surface float makes a heave motion in surface waves, the submerged glider also makes a heave motion by the umbilical cable. When the submerged glider rises, the hydrofoils rotate counterclockwise, and when the submerged glider falls, the hydrofoils rotate clockwise, as shown in Figure 1. The hydrodynamic forces generated by the relative movement between hydrofoils and seawater can be divided into vertical force and forward force. The forward hydrodynamic force will drive the wave glider forward. Some researchers analyzed the hydrodynamic performance of tandem hydr the propulsor of the wave glider considering some factors such as oscillation am and torsional spring stiffness [2][3][4][5]. Moreover, some researchers used nu simulation based on the mathematical model for performance evaluation of th glider. Kraus and Bingham proposed a two-dimensional model of the wave glid along-track and heave directions, which was simplified as a coupled rigid body [6 constructed a 6-DOF (Degree of Freedom) nonlinear dynamic model and id several key hydrodynamic coefficients by using analytical and experimental meth Qi et al. established a 3-DOF dynamic model in the longitudinal plane by the Kane [8]. Li et al., respectively, established the models of the surface float, the submerge and the connection umbilical based on multi-body dynamic theory and compl motion simulation by using ADAMS and MATLAB [9]. Rampersadh described body system of the wave glider using Denavit-Hartenberg parametrizati established the dynamic model for bodies by using the Lagrangian approach [10 et al. presented a 4-DOF model of the wave glider, which included the second-ord drift force on the horizontal plane and the first-order wave force in the vertical d [11]. Chen et al. described a nonlinear kinetic model in three reference coordinat with 6-DOF, analyzed the static stress and optimized the hydrodynamics of the h structure [12]. L. Wang et al. established a dynamic model of the wave glider base analysis of the state transfer problem of the umbilical [13].
The above-mentioned research works assume that the hydrodynamic perfo are almost not influenced by the pitching angles of the hydrofoils and the sub glider, and it is convenient to analyze the propulsive performance of the sub glider. However, for the hinge joint between the umbilical and the submerged gl moment of the submerged glider is unbalanced in the actual sailing. Under the a restoring the moment and hydrodynamic moment of the submerged glider, it periodical change of the pitching angle of the submerged glider in a certain range the hydrofoils pivoted on the submerged glider, the pitching angle of the sub glider is also an important factor for the hydrodynamic propulsion performanc wave glider. The aforementioned research works also assume that adde coefficients are constant and several coupled added mass coefficients of the sub glider are ignored. However, due to the dynamic change of the hydrofoil du motion of the wave glider, added mass coefficients are dependent on the body's g Some researchers analyzed the hydrodynamic performance of tandem hydrofoils as the propulsor of the wave glider considering some factors such as oscillation amplitude and torsional spring stiffness [2][3][4][5]. Moreover, some researchers used numerical simulation based on the mathematical model for performance evaluation of the wave glider. Kraus and Bingham proposed a two-dimensional model of the wave glider in the along-track and heave directions, which was simplified as a coupled rigid body [6]. Kraus constructed a 6-DOF (Degree of Freedom) nonlinear dynamic model and identified several key hydrodynamic coefficients by using analytical and experimental methods [7]. Qi et al. established a 3-DOF dynamic model in the longitudinal plane by the Kane method [8]. Li et al., respectively, established the models of the surface float, the submerged glider and the connection umbilical based on multi-body dynamic theory and completed the motion simulation by using ADAMS and MATLAB [9]. Rampersadh described a multi-body system of the wave glider using Denavit-Hartenberg parametrization and established the dynamic model for bodies by using the Lagrangian approach [10]. Wang et al. presented a 4-DOF model of the wave glider, which included the second-order wave drift force on the horizontal plane and the first-order wave force in the vertical direction [11]. Chen et al. described a nonlinear kinetic model in three reference coordinate frames with 6-DOF, analyzed the static stress and optimized the hydrodynamics of the hydrofoil structure [12]. L. Wang et al. established a dynamic model of the wave glider based on the analysis of the state transfer problem of the umbilical [13].
The above-mentioned research works assume that the hydrodynamic performances are almost not influenced by the pitching angles of the hydrofoils and the submerged glider, and it is convenient to analyze the propulsive performance of the submerged glider. However, for the hinge joint between the umbilical and the submerged glider, the moment of the submerged glider is unbalanced in the actual sailing. Under the action of restoring the moment and hydrodynamic moment of the submerged glider, it causes a periodical change of the pitching angle of the submerged glider in a certain range. Due to the hydrofoils pivoted on the submerged glider, the pitching angle of the submerged glider is also an important factor for the hydrodynamic propulsion performance of the wave glider. The aforementioned research works also assume that added mass coefficients are constant and several coupled added mass coefficients of the submerged glider are ignored. However, due to the dynamic change of the hydrofoil during the motion of the wave glider, added mass coefficients are dependent on the body's geometry and will change with the pitching angle of the hydrofoils.
The key step in numerical simulation is to determine the hydrodynamic coefficients applied in the mathematical model. With the rapid development of computing technology, it is feasible to analyze fluid-structure interaction problems based on CFD. Slone proposed a three-dimensional finite volume, unstructured mesh (FV-UM) method for dynamic fluid-structure interaction [14]. Westphalen, respectively, simulated the wave-structure interaction and wave loading by using the Finite Volume (FV) method, Control-Volume Finite Element (CV-FE) method, Smoothed Particle Hydrodynamics (SPH) and Cartesian Cut Cell method [15]. The Eulerian mesh method is limited by the grid size, but SPH is not limited, and it is better able to capture the jet. Ren proposed a two-dimensional SPH-DEM model to simulate the wave-structure interaction on a slope and predict the hydraulic stability of the 2D structure [16]. Xiang simulated the fluid impact forces on a bridge superstructure in the horizontal and vertical directions using the Finite Element-based Arbitrary Lagrangian-Eulerian Method [17]. Hasanpour presented a coupled SPH-FEM modeling approach that simulated the fluid with particles, and the flume, the debris and the structure with mesh-based finite elements [18]. Past studies have identified the importance and complexity of the dynamic fluid-structure interaction for different types of marine structures and this effect will be thoroughly investigated for wave gliders in our future research.
For hydrodynamic coefficients, CFD simulation is also more accurate compared to analytical and semi-empirical methods for complex shapes. Zhang et al. proposed a new method by using dynamic meshes in FLUENT to simulate hydrodynamic tests to obtain corresponding hydrodynamic coefficients of autonomous underwater vehicle (AUV) [19]. Malik and Guang obtained hydrodynamic coefficients for 6:1 prolate spheroid by simulating the planar motion mechanism (PMM) tests [20]. They maintained the domain mesh qualities during the movement of the rigid body by using three different mesh updating schemes and zoning of the computational field. Towing tests, PMM tests and rotating arm tests were simulated to obtain the hydrodynamic coefficients of AUV [21]. Singh et al. used CFD to estimate the hydrodynamic coefficients of an underwater glider, which was verified by comparing the experimental data of the gliding trajectory with the simulation results [22].
In this paper, the dynamic models of the surface float, the submerged glider with the hydrofoils and the umbilical were separately established in Section 2, which considered the influences of the dynamic variation of the pitching angles of the hydrofoils and the submerged glider. In Section 3, the hydrodynamic coefficients varying with the pitching angle of the hydrofoils were analyzed from static test simulations and dynamic PMM test simulations in FLUENT by using the Finite Volume Method (FVM). In Section 4, sea trial with the "Black Pearl" wave glider was conducted, and the numerical simulations under different sea conditions were carried out to predict the longitudinal velocity of the wave glider and the pitching angle of the submerged glider. The simulation results were reasonable by comparing with the sea trial data. Eventually, some conclusions were summarized in Section 5.

The Reference Frames
To facilitate modeling, three coordinate systems based on the right-handed rule were defined [10], as shown in Figure 2.
(1) Float frame (with the right-superscript F): the hinge point between the surface float and the umbilical is as the origin of the coordinate frame. x F , y F and z F are positive toward the bow, the starboard and the bottom of the surface float, respectively; (2) Glider frame (with the right-superscript G): the hinge point between the submerged glider and the umbilical is as the origin of the glider frame. x G , y G and z G positive toward the bow, the starboard and the bottom of the submerged glider, respectively; (3) Inertial reference frame (ξ − η − ζ): ξ, η and ζ point toward north, east and downwards, respectively.

Figure 2.
The coordinate frames of the "Black Pearl" wave glider.

Dynamic Model of the Surface Float
Xiang revealed a great dependence of the wave loads on the ratio of wavele structure width [23]. Rudman used the Smoothed Particle Hydrodynamics technique to simulate the influence on a semi-submersible platform at a range of d angles [24]. The interaction between the surface float and wave was simplified wave incident angle was assumed to be zero in this study. Moreover, due to the d of predicting numerically irregular waves in the experimental sea area, the surfac were assumed as regular waves. However, future studies will consider more conditions, such as irregular waves and breaking waves or bores. The latter type o is challenging because they can be dominated by significant turbulence effe entrapment and high-aleatory uncertainty [25][26][27] and can introduce totally d hydrodynamic loads compared to unbroken waves [28,29].
The heave motion of the surface float is assumed as a sinusoidal motion with waves. The longitudinal model of the surface float is formulated as follows [30]:   The wave force is simply considered as an excitation force, which can be de as: where f is the wave force per unit wave amplitude, A is the amplitude of wa

Dynamic Model of the Surface Float
Xiang revealed a great dependence of the wave loads on the ratio of wavelength to structure width [23]. Rudman used the Smoothed Particle Hydrodynamics (SPH) technique to simulate the influence on a semi-submersible platform at a range of different angles [24]. The interaction between the surface float and wave was simplified and the wave incident angle was assumed to be zero in this study. Moreover, due to the difficulty of predicting numerically irregular waves in the experimental sea area, the surface waves were assumed as regular waves. However, future studies will consider more realistic conditions, such as irregular waves and breaking waves or bores. The latter type of waves is challenging because they can be dominated by significant turbulence effects, air-entrapment and high-aleatory uncertainty [25][26][27] and can introduce totally different hydrodynamic loads compared to unbroken waves [28,29].
The heave motion of the surface float is assumed as a sinusoidal motion with surface waves. The longitudinal model of the surface float is formulated as follows [30]: where M F RB and M F A are the inertia matrix and the added mass matrix of the surface float in the rigid-body system, respectively. ν F is the velocity vector of the surface float and ν F r is the velocity vector of the surface float relative to ocean current. C F RB ν F is the rigidbody Coriolis and Centripetal matrix of the surface float and C F A ν F is the hydrodynamic Coriolis and Centripetal matrix of the surface float. D F ν F is the damping matrix of the surface float, τ 1_umb is the vector of umbilical force and τ F wave is the vector of wave force on the surface float.
The wave force is simply considered as an excitation force, which can be described as: where f w is the wave force per unit wave amplitude, A w is the amplitude of wave high, ω w is the wave excitation frequency, and ϕ w is the phase.

Dynamic Model of the Submerged Glider
The longitudinal 3-DOF (surge, heave and pitch) coupled model of the submerged glider can be expressed as: where M G RB and M G A are the inertia matrix and the added mass matrix of the submerged glider in the rigid-body system, respectively. ν G is the velocity vector of the submerged glider and ν G r is the velocity vector of the submerged glider relative to the ocean current. C G RB ν G is the rigid-body Coriolis and Centripetal matrix of the submerged glider and C G A ν G is the hydrodynamic Coriolis and Centripetal matrix of the submerged glider. D G ν G is the hydrodynamic damping matrix of the submerged glider and g(η) is the vector of restoring forces and moment of the submerged glider. τ 2_umb is the vector of umbilical force on the submerged glider and τ hf is the vector of hydrodynamic forces and moment of the hydrofoils. The submerged glider is symmetric about the xz-plane and longitudinal motion is only considered; thus M G RB , M G A , C G RB ν G and C G A ν G are, respectively, simplified as: and the coefficients of the hydrodynamic Coriolis and Centripetal matrix of the submerged glider (a G 1 and a G 3 ) are, respectively, written as: Hydrodynamic damping matrix D G ν G r includes the linear part D G and the nonlinear part D G n ν G r . D G n ν G r only includes additional nonlinear damping terms on the diagonal in surge and heave because other coupling terms are small, which are ignored in this model. Thus, D G and D G n ν G r are, respectively, expressed as: The vector of restoring force and moment g(η) is derived as: where W G and B G are the weight and the buoyancy force of the submerged glider, respectively.
The lift force (F L ) and drag force (F D ) of six pairs of tandem hydrofoils are obtained by the simulation. The vectors of drag and lift forces change with the vector of flow velocity, which is shown in Figure 3. During the falling phase of the glider, the drag force and the lift force of the hydrofoils in simulations can be resolved to the horizontal and vertical forces in the glider frame, which can be expressed as: The vector of restoring force and moment ( ) g η is derived as: where G W and G B are the weight and the buoyancy force of the submerged respectively.
The lift force ( L F ) and drag force ( D F ) of six pairs of tandem hydrofoils are ob by the simulation. The vectors of drag and lift forces change with the vector o velocity, which is shown in Figure 3. During the falling phase of the glider, the dra and the lift force of the hydrofoils in simulations can be resolved to the horizon vertical forces in the glider frame, which can be expressed as: sin cos cos sin During the rising phase of the glider, the drag force and the lift force of the hyd in simulations can be resolved to the horizontal and vertical forces in the glider which can be expressed as: The vector of hydrodynamic forces and moment of the hydrofoils can be writ M is the hydrodynamic moment of the hydrofoils acting on the subm glider. Figure 4 shows the force diagram of the hydrofoil, which mainly includ hydrostatics, hydrodynamic force and spring tension. The pitching acceleration hydrofoils can be expressed as:  During the rising phase of the glider, the drag force and the lift force of the hydrofoils in simulations can be resolved to the horizontal and vertical forces in the glider frame, which can be expressed as: The vector of hydrodynamic forces and moment of the hydrofoils can be written as: where M G h f is the hydrodynamic moment of the hydrofoils acting on the submerged glider. Figure 4 shows the force diagram of the hydrofoil, which mainly includes the hydrostatics, hydrodynamic force and spring tension. The pitching acceleration of the hydrofoils can be expressed as: where Approximating each pair of the hydrofoils as a series of flat strips aligned w hydrofoil chords, the pitching added mass of each pair of the hydrofoils ( hf q M  expressed as: where L hf denotes the chord-wise cross-sectional length as a function of distan the hydrofoil axis [31].

Umbilical Model
When the surface float does heave motion in surface waves, the umbilica always under tension. Comparing the distance between the surface float submerged glider ( FG d ) with the initial length of the umbilical ( umb L ), the stat umbilical can be judged. The umbilical is assumed as a spring and the tensio umbilical can be expressed as: where θ G h f is the pitching angle of the hydrofoils in the glider frame, θ n G the pitching angle of the submerged glider in the inertial reference frame, θ G s is the angle between the positive x-axis of the submerged glider and the spring, L f is the distance from the center of mass of a single pair of tandem hydrofoils to their pivot, L s is the distance from the pin attached to the spring to the pivot, k s is the spring stiffness coefficient, and ∆L sp is the elongation of the spring. θ G s can be expressed as: Approximating each pair of the hydrofoils as a series of flat strips aligned with the hydrofoil chords, the pitching added mass of each pair of the hydrofoils (M . q h f ) can be expressed as: where L h f denotes the chord-wise cross-sectional length as a function of distance along the hydrofoil axis [31].

Umbilical Model
When the surface float does heave motion in surface waves, the umbilical is not always under tension. Comparing the distance between the surface float and the submerged glider (d FG ) with the initial length of the umbilical (L umb ), the state of the umbilical can be judged. The umbilical is assumed as a spring and the tension of the umbilical can be expressed as: (22) where k umb is the stiffness coefficient of the umbilical. k umb can be expressed as: The umbilical force vector acting on the surface float can be expressed as: The angle between the umbilical and the z-direction of the surface float can be expressed as: where ξ F and ξ G are the longitudinal coordinates of the surface float and the submerged glider, respectively. ζ F and ζ G are the vertical coordinates of the surface float and the submerged glider, respectively. The umbilical force vector acting on the submerged glider can be expressed as: The angle between the umbilical and the x-direction of the submerged glider can be expressed as:

Simulation for Hydrodynamic Coefficients
Different from an underwater glider, the hydrofoils of the submerged glider can rotate in a certain angle range, which leads to the change of its shape. Therefore, the hydrodynamic coefficients of the submerged glider are related to the pitching angle of the hydrofoils. To reasonably reflect the propulsion performance of the wave glider, it is necessary to identify the unknown hydrodynamic coefficients in dynamic models. To define these coefficients, oblique motions in the vertical plane and the VPMM tests are simulated by using the computational fluid dynamics software FLUENT based on the RANS method. In simulations, the pitching angle of the hydrofoils ranges from −15 • to 30 • .

Turbulent Model
The Reynolds stress term in RANS equation has six independent elements and it is called as a closure problem of the RANS equation. Therefore, various turbulence models are used to find solutions to the Reynolds stress in terms of known quantities to allow closure of the RANS equations [32]. The realizable k − ε turbulence model is used for simulations in this study due to its robustness and ability to be applied to external flows with a complex geometrical shape.

Computational Domain Decomposition and Mesh Generation
To maintain the domain mesh qualities, the rectangular computational domain is divided into three zones, the inner, middle and outer zones, and data transfer between different zones is implemented by defining interface pairs. Considering the complex shape of the submerged glider, the unstructured mesh is used for the spherical inner domain containing the submerged glider and the mesh type is chosen as Tetra/Mixed type, and inflated prism layers are generated surrounding the surface of the submerged glider to improve the computational accuracy. The thickness of the first prism layer can be expressed as: where ρ and µ are, respectively, the fluid density and viscosity, the frictional velocity U τ = τ w /ρ, the wall shear stress τ w = 1 2 C f ρU 2 , U is the freestream velocity, and the skin friction coefficient of external flows is estimated by using the following equation: From the above equation, the thickness of the first prism layer is located such that is kept in a range of 30~300. In addition, structured mesh is used for the middle and outer zones. These meshes are produced by ANSYS ICEM meshing software. Moreover, the size of the whole computational domain is 20m × 12m × 12m. Figure 5 shows the division diagram of the computational domain where L, W and H are, respectively, the length, the width and the height of the submerged glider. Table 1 gives the average grid sizes of different zones. Moreover, the same mesh generation strategy is adopted for submerged glider models with different pitching angles of the hydrofoils.
domain containing the submerged glider and the mesh type is chosen as Tetra/Mixed type, and inflated prism layers are generated surrounding the surface of the submerged glider to improve the computational accuracy. The thickness of the first prism layer can be expressed as: where ρ and μ are, respectively, the fluid density and viscosity, the frictional velocity w U τ τ ρ = , the wall shear stress , U is the freestream velocity, and the skin friction coefficient of external flows is estimated by using the following equation: From the above equation, the thickness of the first prism layer is located such that is kept in a range of 30~300. In addition, structured mesh is used for the middle and outer zones. These meshes are produced by ANSYS ICEM meshing software. Moreover, the size of the whole computational domain is 20m 12m 12m × × . Figure 5 shows the division diagram of the computational domain where L, W and H are, respectively, the length, the width and the height of the submerged glider. Table 1 gives the average grid sizes o different zones. Moreover, the same mesh generation strategy is adopted for submerged glider models with different pitching angles of the hydrofoils.  Owing to the movement of the body, numerical simulations of the VPMM tests are regarded as transient; thus, the dynamic mesh capability and User Defined Functions (UDF) of FLUENT software are utilized to ensure the movement of the body to supply the demand of the VPMM tests. In dynamic simulations, the layering method of mesh up dating is adopted. The inner zone completely follows the motion of the submerged glider defined by the UDF, and the middle zone follows the heave motion of the submerged  Owing to the movement of the body, numerical simulations of the VPMM tests are regarded as transient; thus, the dynamic mesh capability and User Defined Functions (UDF) of FLUENT software are utilized to ensure the movement of the body to supply the demand of the VPMM tests. In dynamic simulations, the layering method of mesh up-dating is adopted. The inner zone completely follows the motion of the submerged glider defined by the UDF, and the middle zone follows the heave motion of the submerged glider. The outer zone remains fixed. The static test simulations are regarded as steady, and the computational domain can also be used in static test simulations.

Boundary Conditions
Numerical simulations for static tests are carried out at different angles of attack (angle of attack in this paper is defined as the angle between relative flow and the x-direction of the submerged glider). When it is 0 • or ±180 • , only one velocity-inlet and pressureoutlet boundary is used. When it is at other angles of attack, two velocity-inlet and two pressure-outlet boundaries are used. The velocity magnitude of velocity-inlet boundary is set at 1 m/s with flow turbulence of 5%. Components of flow velocity in surge and heave directions are dependent on the angle of attack. Moreover, the outlet boundary condition is with zero relative pressure.
In the VPMM test simulations, the inlet boundary condition is located at the +x boundary in front of the submerged glider and the outlet boundary condition with zero relative pressure is located at the −x boundary behind the submerged glider.
All simulations are performed by the FLUENT's double precision segregated pressurebased solver. The second order upwind scheme is used, and pressure-velocity coupling is adopted by SIMPLEC algorithm.

Grid Convergence
The grid convergence study is performed by following the procedure provided by Celik et al. To analyze the grid convergence, three different mesh resolutions were used for simulations maintaining a refinement ratio of 1.3 or more. Take the submerged glider model with 10 • pitching angle of the hydrofoils; for example, the used mesh resolutions were 2.37 million, 5.94 million and 14.36 million.
For the analysis, the order of accuracy p and Grid Convergence Index (GCI) were calculated by the following equations [33]: where The grid convergence analysis is shown in Table 2. The values of GCI are all small, it indicates that the computation is within the asymptotic range. In order to save computation time, a coarser mesh of 2.37 million was selected for the simulation.

Static Test Simulations
When the submerged glider rises under the tension of the umbilical or falls under the action of gravity, high positive angle of attack or high negative angle of attack will frequently occur. Therefore, in a steady-state flow field, numerical simulations for static tests under different angles of attack and pitching angles of the hydrofoils were performed with the 1m/s velocity magnitude to predict the hydrodynamic force and moment under any movement of the submerged glider.
In static test simulations, the submerged glider with six pairs of tandem hydrofoils was used and simulation data were the resultant force of the thrust and the drag force of the submerged glider. The translational hydrodynamic coefficient term τ G uw θ G h f , θ G ν was adopted to dispose the coupled large angle of attack of the submerged glider and the pitching angle of the hydrofoils, which includes the effects of the surge and heave velocities of the whole glider [34,35].
can be expressed as: where Therefore, the rotational damping coefficients D G R can be separated from hydrodynamic damping matrix D G ν G r and be identified by pure pitch test simulations, which can be expressed as: In order to use the acquired simulation data, some hydrodynamic terms in the dynamic model of the submerged glider are adjusted, thus Equation (3) can be rewritten as: The static tests for the pitching angle of the hydrofoils ranging from −15 • to 30 • and the angle of attack of the submerged glider ranging from −180 • to 180 • were conducted to identify the values of the translational hydrodynamic coefficient term. A high-order polynomial model was used for fitting of simulation data. Simulation cases are shown in Figure 6. The simulation and fitting results are shown in Figure 7. As shown the pitching angle of the hydrofoils is positive, and the angle of attac approximately 60 degrees to 120 degrees, or when the pitching angle The simulation and fitting results are shown in Figure 7. As shown in Figure 7a, when the pitching angle of the hydrofoils is positive, and the angle of attack is in the range of approximately 60 degrees to 120 degrees, or when the pitching angle of the hydrofoils is negative and the angle of attack is in the range of approximately −60 degrees to −120 degrees, the submerged glider can generate relatively large positive thrust. As shown in Figure 7b,c, the angle of attack has a greater effect than the pitching angle of the hydrofoils. The simulation and fitting results are shown in Figure 7. As shown in Figure 7a, when the pitching angle of the hydrofoils is positive, and the angle of attack is in the range of approximately 60 degrees to 120 degrees, or when the pitching angle of the hydrofoils is negative and the angle of attack is in the range of approximately −60 degrees to −120 degrees, the submerged glider can generate relatively large positive thrust. As shown in Figure 7b,c, the angle of attack has a greater effect than the pitching angle of the hydrofoils.

Simulation of the VPMM Test
Simulation of the VPMM test, including pure surge simulation, pure heave simulation and pure pitch simulation, are used to determine added mass coefficients and rotational damping coefficients. In simulations, oscillation amplitudes were set at 0.04 m, which could ensure reasonable accuracy of the added mass and damping coefficients [36]. Moreover, oscillation angular frequency was adjusted for 0.1 Hz, 0.125 Hz, 0.2 Hz, 0.25 Hz, 0.3125 Hz, 0.4 Hz and 0.5 Hz, respectively.

Pure Surge Simulation
In the pure surge test, the submerged glider is sinusoidally oscillated in the x-direction, which is shown in Figure 8.

Pure Surge Simulation
In the pure surge test, the submerged glider is sinusoidally oscillated in direction, which is shown in Figure 8. x is surge oscillation amplit Figure 9 shows the pressure contours of the submerged glider in pur simulation when pitching angle of the hydrofoils is 15 degrees. The pure surge test is simulated in FLUENT based on the dynamic mesh technique, with restrictions on other degrees of freedom, and the pure surge motion is defined as: where x, u and . u are, respectively, displacement, velocity and acceleration in surge direction, ω is oscillation angular frequency, and x 0 is surge oscillation amplitude. Figure 9 shows the pressure contours of the submerged glider in pure surge simulation when pitching angle of the hydrofoils is 15 degrees.

Pure Surge Simulation
In the pure surge test, the submerged glider is sinusoidally oscillated in the x direction, which is shown in Figure 8. x is surge oscillation amplitude. Figure 9 shows the pressure contours of the submerged glider in pure surg simulation when pitching angle of the hydrofoils is 15 degrees.  Equation (37) includes the effects of rigid-body forces induced by M G RB and C G RB ν G To determine the hydrodynamic forces, it is indispensable to remove rigid-body forces. Moreover, substituting Equation (38) into Equation (37), the hydrodynamic force acting on the submerged glider can be obtained as: Simulations with different oscillation angular frequencies and pitching angles of the hydrofoils were carried out. Take the submerged glider model with 15 • pitching angle of the hydrofoils as an example, the hydrodynamic forces and moment were fitted by Equation (39) firstly, and then fitting results were fitted by first-order polynomial as shown in Figure 10. The coefficients of the term −x 0 ω 2 are . u dependent added mass coefficients.
Simulations with different oscillation angular frequencies and pi hydrofoils were carried out. Take the submerged glider model with of the hydrofoils as an example, the hydrodynamic forces and mom Equation (39) firstly, and then fitting results were fitted by first-o shown in Figure 10. The coefficients of the term 2 0 x ω − are u  depe coefficients. Combining all fitting results of pure surge simulation data, u  mass coefficients under different pitching angles of the hydrofoils are Figure 11. The results of u  dependent added mass coefficients. Combining all fitting results of pure surge simulation data, . u dependent added mass coefficients under different pitching angles of the hydrofoils are shown in Figure 11. of the hydrofoils as an example, the hydrodynamic forces and mom Equation (39) firstly, and then fitting results were fitted by first-o shown in Figure 10. The coefficients of the term 2 0

Pure Heave Simulation
x ω − are u  dep coefficients. Combining all fitting results of pure surge simulation data, u  mass coefficients under different pitching angles of the hydrofoils are

Pure Heave Simulation
In the pure heave simulation, the submerged glider is moving forward at a specified velocity (u 0 = 1 m/s) and its x-axis always remains horizontal while the submerged glider is sinusoidally oscillated in the z-direction, which is shown in Figure 12.
In the pure heave simulation, the submerged glider is moving forward at a specified velocity ( 0 1m/s u = ) and its x-axis always remains horizontal while the submerged glider is sinusoidally oscillated in the z-direction, which is shown in Figure 12.
where z , w and w  are the displacement, velocity and acceleration in heave direction, respectively. 0 z is the heave oscillation amplitude. Figure 13 shows the pressure contours of the submerged glider in pure heave simulation when the pitching angle of the hydrofoils is 15 degrees. Substituting Equation (40) into Equation (37), the hydrodynamic forces and moment acting on the submerged glider can be obtained as: Combining all fitting results of pure heave simulation data by the same method, thus the w  dependent added mass coefficients under different pitching angles of the hydrofoils are shown in Figure 14.  The pure heave motion is defined as: where z, w and . w are the displacement, velocity and acceleration in heave direction, respectively. z 0 is the heave oscillation amplitude. Figure 13 shows the pressure contours of the submerged glider in pure heave simulation when the pitching angle of the hydrofoils is 15 degrees.

Pure Pitch Simulation
In the pure pitch simulation, the submerged glider moving forward at a s velocity ( 0 1m/s u = ) is sinusoidally oscillated in the z-direction and its x-axis is tangential to the sinusoidal path, which is shown in Figure 15. It is intended to obtain q  dependent added mass coefficients and q dep damping coefficients under different pitching angles of the hydrofoils. The pu motion is defined as: where θ is the pitching angle, q is the pitch angular velocity, q  is th acceleration, u is the surge velocity and u  is the surge acceleration. Figure 16 sh pressure contours of the submerged glider in pure pitch simulation when the p angle of the hydrofoils is 15 degrees.

Pure Pitch Simulation
In the pure pitch simulation, the submerged glider moving forward at a specified velocity (u 0 = 1 m/s) is sinusoidally oscillated in the z-direction and its x-axis is always tangential to the sinusoidal path, which is shown in Figure 15.

Pure Pitch Simulation
In the pure pitch simulation, the submerged glider moving forward at a specified velocity ( 0 1m/s u = ) is sinusoidally oscillated in the z-direction and its x-axis is always tangential to the sinusoidal path, which is shown in Figure 15. It is intended to obtain q  dependent added mass coefficients and q dependent damping coefficients under different pitching angles of the hydrofoils. The pure pitch motion is defined as: where θ is the pitching angle, q is the pitch angular velocity, q  is the pitch acceleration, u is the surge velocity and u  is the surge acceleration. Figure 16 shows the pressure contours of the submerged glider in pure pitch simulation when the pitching angle of the hydrofoils is 15 degrees. It is intended to obtain . q dependent added mass coefficients and q dependent damping coefficients under different pitching angles of the hydrofoils. The pure pitch motion is defined as: where θ is the pitching angle, q is the pitch angular velocity, . q is the pitch acceleration, u is the surge velocity and . u is the surge acceleration. Figure 16 shows the pressure contours of the submerged glider in pure pitch simulation when the pitching angle of the hydrofoils is 15 degrees. The hydrodynamic forces and moment acting on the submerged glider can be expressed as: Combining all fitting results of pure pitch simulation data by the same method, the q  dependent added mass coefficients under different pitching angles of the hydrofoils are shown in Figure 17a, and the q dependent added damping coefficient is shown in Figure 17b.  The hydrodynamic forces and moment acting on the submerged glider can be expressed as: Combining all fitting results of pure pitch simulation data by the same method, the . q dependent added mass coefficients under different pitching angles of the hydrofoils are shown in Figure 17a, and the q dependent added damping coefficient is shown in Figure 17b. The hydrodynamic forces and moment acting on the submerged glider can be expressed as: Combining all fitting results of pure pitch simulation data by the same method, the q  dependent added mass coefficients under different pitching angles of the hydrofoils are shown in Figure 17a, and the q dependent added damping coefficient is shown in Figure 17b.

Numerical Simulation Validation and Analysis
To verify the numerical simulation, simulation results obtained from solution of the dynamic model using the hydrodynamic coefficients are compared with the experimental data in the straight-line paths.
The non-dynamic parameters of the "Black Pearl" wave glider are used in both the sea trial and numerical simulations, which are given in Table 3. Table 3. Non-dynamic parameters of the "Black Pearl" wave glider.

Sea Trial
In December 2020, sea trial was carried out to test the sailing performance of the "Black Pearl" wave glider, as shown in Figure 18.

Sea Trial
In December 2020, sea trial was carried out to test the sailing perfo "Black Pearl" wave glider, as shown in Figure 18. Considering the degrees of freedom in the dynamic model, a desire path between point T1 and point T2 was set. The positions and the velocitie Pearl" wave glider can be obtained by using a GPS. Figure 19 gives the trajectory of the wave glider represented by the blue solid line and the straight-line path in sea trial. And the red box includes the section of sens work, these data were obtained in 2 1 T T  path because the difference betwe course and the wave direction was relatively stable in this path. The sailing T2 to T1 is about 5.5 km, and the sailing time is about 3 h. Considering the degrees of freedom in the dynamic model, a desired straight-line path between point T1 and point T2 was set. The positions and the velocities of the "Black Pearl" wave glider can be obtained by using a GPS. Figure 19 gives the actual sailing trajectory of the wave glider represented by the blue solid line and the red box is the straight-line path in sea trial. And the red box includes the section of sensor data. In this work, these data were obtained in → T 2 T 1 path because the difference between the desired course and the wave direction was relatively stable in this path. The sailing distance from T2 to T1 is about 5.5 km, and the sailing time is about 3 h. trajectory of the wave glider represented by the blue solid line and the red box straight-line path in sea trial. And the red box includes the section of sensor data. work, these data were obtained in 2 1 T T  path because the difference between the d course and the wave direction was relatively stable in this path. The sailing distanc T2 to T1 is about 5.5 km, and the sailing time is about 3 h. The sea conditions during the sea trial from T2 to T1 were obtained by using sensor; wave data are shown in Figure 20a. Moreover, the real-time values and the m average value of the pitching angle of the submerged glider are shown in Figure  can be seen that the average speed of the wave glider increases correspondingly w continuous increase of significant wave periods and wave heights. Moreover, the a pitching angle and the pitching angle amplitude of the submerged glider both in slightly. The sea conditions during the sea trial from T2 to T1 were obtained by using a wave sensor; wave data are shown in Figure 20a. Moreover, the real-time values and the moving average value of the pitching angle of the submerged glider are shown in Figure 20b. It can be seen that the average speed of the wave glider increases correspondingly with the continuous increase of significant wave periods and wave heights. Moreover, the average pitching angle and the pitching angle amplitude of the submerged glider both increase slightly.

Numerical Simulation with the Specific Sea State
Rampersadh calculated resultant propulsion effect on the float under wave height of 1m and wave period of 10 s [10]. In order to compare with Rampersadh's calculation results, the dynamic characteristics of the wave glider under the same sea state were examined.
In Figure 21a, the blue, red and blank solid lines, represent the trajectory of the surface float, the trajectory of the submerged glider and the umbilical, respectively. A comparison between Rampersadh's calculation result and this study's simulation result is shown in Figure 21b.

Numerical Simulation with the Specific Sea State
Rampersadh calculated resultant propulsion effect on the float under wave height of 1m and wave period of 10 s [10]. In order to compare with Rampersadh's calculation results, the dynamic characteristics of the wave glider under the same sea state were examined.
In Figure 21a, the blue, red and blank solid lines, represent the trajectory of the surface float, the trajectory of the submerged glider and the umbilical, respectively. A comparison between Rampersadh's calculation result and this study's simulation result is shown in Figure 21b. results, the dynamic characteristics of the wave glider under the same sea state were examined.
In Figure 21a, the blue, red and blank solid lines, represent the trajectory of the surface float, the trajectory of the submerged glider and the umbilical, respectively. A comparison between Rampersadh's calculation result and this study's simulation result is shown in Figure 21b. Under this specific sea state, the longitudinal average velocities in this study and Rampersadh's study were, respectively, about 0.41 m/s and 0.4 m/s.

Numerical Simulation Results with Different Current Velocities
The wave height of 0.8 m and the wave period of 4.5 s, which are the average values of experimental data in the sea trial area, are the typical sea state input values of this simulation to analyze the impact of current velocity [37]. Figure

Numerical Simulation Results with Different Current Velocities
The wave height of 0.8 m and the wave period of 4.5 s, which are the average values of experimental data in the sea trial area, are the typical sea state input values of this simulation to analyze the impact of current velocity [37]. Figure  As seen in Figure 22, the positive current velocity is conducive to improv navigation velocity of the wave glider, but the negative current velocity will red navigation performance of the wave glider, reflecting the bigger hydrodynamic da In addition, the greater the positive or negative current velocity, the greater the fa or hindering effect on the navigation performance of the wave glider.

Numerical Simulation Results with the Different Sea States
To study the effects of wave height and wave period on sailing performanc submerged glider, the simulated sea states through a two-by-two combination o height and wave period are given in Table 4. Due to the complexity and variabilit marine environment, the interaction between the surface float and wave was sim The initial speeds of the surface float and the submerged glider and the initial p As seen in Figure 22, the positive current velocity is conducive to improving the navigation velocity of the wave glider, but the negative current velocity will reduce the navigation performance of the wave glider, reflecting the bigger hydrodynamic damping. In addition, the greater the positive or negative current velocity, the greater the favorable or hindering effect on the navigation performance of the wave glider.

Numerical Simulation Results with the Different Sea States
To study the effects of wave height and wave period on sailing performance of the submerged glider, the simulated sea states through a two-by-two combination of wave height and wave period are given in Table 4. Due to the complexity and variability of the marine environment, the interaction between the surface float and wave was simplified.
The initial speeds of the surface float and the submerged glider and the initial pitching angle of the umbilical were all set to 0. The initial pitching angles of the submerged glider and the hydrofoils were set by calculating coupled hydrostatics of the submerged glider and the hydrofoils. The actual current velocity in the sea trial is unknown because the ADCP is not carried on the surface float. Thus, the current velocity in simulations is set to 0. degrees when the sea states were the same as the sea trial data. Moreover, the average pitching angle for the high wave height changes significantly in different wave periods.
In addition, it can be seen from Figure 23c that the high wave period is beneficial to reduce the pitch amplitude, especially for the high wave height. Moreover, the low wave height with the same wave period is conducive to decreasing the pitching angle amplitude of the submerged glider. Figure 23d presents the standard deviation of the pitching angle of the submerged glider under different sea states, which describes the deviation magnitude of the centerline of the pitching angle from its average value. The low wave period is the main factor leading to the increase of the standard deviation. Furthermore, the forward sailing performance of the wave glider is sensitive to wave period, especially in a high wave height with a low wave period. It can be seen from Figure  23b,c that under a high wave height with a low wave period, the average value, the amplitude and the standard deviation of the pitching angle suddenly change. Therefore, it means that the submerged glider is in an unstable state, which will lead to a significant reduction in the longitudinal average velocity.  The average pitching angle and the pitching angle amplitude of the submerged glider are shown in Figure 23b,c, respectively. It can be seen from Figure 23b that the average pitching angle of the submerged glider almost maintained between 15 degrees and 20 degrees when the sea states were the same as the sea trial data. Moreover, the average pitching angle for the high wave height changes significantly in different wave periods.

Numerical Simulation Validation
In addition, it can be seen from Figure 23c that the high wave period is beneficial to reduce the pitch amplitude, especially for the high wave height. Moreover, the low wave height with the same wave period is conducive to decreasing the pitching angle amplitude of the submerged glider. Figure 23d presents the standard deviation of the pitching angle of the submerged glider under different sea states, which describes the deviation magnitude of the centerline of the pitching angle from its average value. The low wave period is the main factor leading to the increase of the standard deviation.
Furthermore, the forward sailing performance of the wave glider is sensitive to wave period, especially in a high wave height with a low wave period. It can be seen from Figure 23b,c that under a high wave height with a low wave period, the average value, the amplitude and the standard deviation of the pitching angle suddenly change. Therefore, it means that the submerged glider is in an unstable state, which will lead to a significant reduction in the longitudinal average velocity.

Numerical Simulation Validation
Numerical simulations were performed to compare the average velocity of the wave glider, the amplitude and the moving average of the pitching angles of the submerged glider with sea trial data. To verify the effectiveness of the numerical simulation, sea states obtained every 15 minutes are used as input conditions for simulations.
The trends of simulated average pitching angles, average peaks, average troughs, and average velocities are shown in Figure 24a, which is similar with experimental data. The trends of simulated average pitching angles, average peaks, average troughs, and average velocities are shown in Figure 24a, which is similar with experimental data.
It can be seen from Figure 24a that simulated average peaks are consistent with the first half of experimental data, but it exists gradually, with increasing upward-offset compared with the other of experimental data. In the first half, the estimated average difference between simulated pitching average and the moving average of the pitching angle of experimental data is about one degree, which is small. However, in the latter part, it is approximately twice the average difference than that of the first half, and the maximum difference is about three degrees. In addition, the downward-offset of simulated average troughs relative to experimental data is about 5 degrees. It can be seen from Figure 24b that simulated longitudinal average velocity is in the upward trend with the same as experimental velocity, but the increasing trend is slower than that of experimental velocity. In the above comparison, it is acceptable in engineering practice to have a certain degree of error between the simulated data and experimental data. The reason for this is that the assumptions used in the numerical simulation are simplifications of the actual conditions. For example, the winds and the ocean currents are ignored, and actual wave conditions are constantly changing, but wave conditions in simulations are constant.

Conclusions
In this paper, to improve the propulsion performance of the wave glider, a 4-DOF dynamic model of the wave glider based on the dynamic characteristics of the tandem hydrofoils is built. The unknown hydrodynamic coefficients in the dynamic model are obtained. The "Black Pearl" wave glider is adopted for the sea trial, and straight-line path is set to verify the numerical simulation. The propulsion performance of the wave glider It can be seen from Figure 24a that simulated average peaks are consistent with the first half of experimental data, but it exists gradually, with increasing upward-offset compared with the other of experimental data. In the first half, the estimated average difference between simulated pitching average and the moving average of the pitching angle of experimental data is about one degree, which is small. However, in the latter part, it is approximately twice the average difference than that of the first half, and the maximum difference is about three degrees. In addition, the downward-offset of simulated average troughs relative to experimental data is about 5 degrees. It can be seen from Figure 24b that simulated longitudinal average velocity is in the upward trend with the same as experimental velocity, but the increasing trend is slower than that of experimental velocity.
In the above comparison, it is acceptable in engineering practice to have a certain degree of error between the simulated data and experimental data. The reason for this is that the assumptions used in the numerical simulation are simplifications of the actual conditions. For example, the winds and the ocean currents are ignored, and actual wave conditions are constantly changing, but wave conditions in simulations are constant.

Conclusions
In this paper, to improve the propulsion performance of the wave glider, a 4-DOF dynamic model of the wave glider based on the dynamic characteristics of the tandem hydrofoils is built. The unknown hydrodynamic coefficients in the dynamic model are obtained. The "Black Pearl" wave glider is adopted for the sea trial, and straight-line path is set to verify the numerical simulation. The propulsion performance of the wave glider is compared in simulation in terms of wave parameters, the velocities and some pitching angles of the wave glider. Some conclusions can be drawn as follows: (1) The proposed model and the numerical simulation using those hydrodynamic coefficients are found to be reasonable and acceptable in practice by comparison with the experimental data; (2) The numerical simulation can predict the longitudinal velocity of the wave glider and the pitching angle of the submerged glider and analyze the stability of the submerged glider from the mean value, the amplitude and the standard deviation of the pitching angle; (3) Under the same sea state, the positive current velocity is conducive to improving the navigation velocity of the wave glider, but the negative current velocity will reduce the navigation performance of the wave glider, reflecting the bigger hydrodynamic damping. In addition, the greater the positive or negative current velocity, the greater the favorable or hindering effect on the navigation performance of the wave glider; (4) The wave height and the wave period are the main factors affecting the straight-line sailing performance of the wave glider. Furthermore, the pitching stability of the submerged glider is also crucial for the wave glider to converse the energy obtained from ocean waves into the forward thrust.
In future research, ocean current velocity will be taken into account in dynamical models. Moreover, the influences of the complex ocean environment on wave gliders will also be considered to be closer to the realistic environment; for example, the wave loads and the motion responses of the surface float under irregular waves and different wave incident angles.

Acknowledgments:
The authors appreciate all staff members in research and development team of the "Black Pearl" wave gliders.

Conflicts of Interest:
The authors declare no conflict of interest.