Hydrodynamic Analysis of Self-Propulsion Performance of Wave-Driven Catamaran

Abstract

The wave-driven catamaran is a small surface vehicle driven by ocean waves. It consists of a hull and hydrofoils, and has a multi-body dynamic structure. The process of moving from static state to autonomous navigation driven by ocean waves is called “self-propulsion”, and reflects the ability of the wave-driven catamaran to absorb oceanic wave energy. Considering the importance of the design of the wave-driven catamaran, its self-propulsion performance should be comprehensively analysed. However, the wave-driven catamaran’s multi-body dynamic structure, unpredictable dynamic and kinematic responses driven by waves make it difficult to analyse its self-propulsion performance. In this paper, firstly, a multi-body dynamic model is established for wave-driven catamaran. Secondly, a two-phase numerical flow field containing water and air is established. Thirdly, a numerical simulation method for the self-propulsion process of the wave-driven catamaran is proposed by combining the multi-body dynamic model with a numerical flow field. Through numerical simulation, the hydrodynamic response, including the thrust of the hydrofoils, the resistance of the hull and the sailing velocity of the wave-driven catamaran are identified and comprehensively analysed. Lastly, the accuracy of the numerical simulation results is verified through a self-propulsion test in a towing tank. In contrast with previous research, this method combines multi-body dynamics with computational fluid dynamics (CFD) to avoid errors caused by artificially setting the motion mode of the catamaran, and calculates the real velocity of the catamaran.

1. Introduction

The wave-driven catamaran is a small surface vehicle driven by ocean waves. It is composed of a catamaran hull and wave oscillating hydrofoils, and has a multi-body dynamic structure. Because ocean waves exist widely and are ever-present, the wave-driven catamaran is consistently powered by ocean waves, and has strong endurance. Compared with conventional marine exploration equipment, the wave-driven vehicle is free of the constraints of conventional power. As a result of its strong endurance, the wave-driven vehicle has now been applied to tasks such as marine meteorological monitoring and prediction, marine biological tracking, seawater quality monitoring, seabed terrain detection, fishery protection, drug and other maritime trafficking surveillance, and even mine clearance [1,2,3].

The development of wave-driven vessels has a long history. A typical representative of this industry is the “Autonaut” company in the United Kingdom, which has developed a series of wave-driven, unmanned vehicles named “Autonaut”. With strong endurance, the “Autonaut” wave-driven unmanned vehicles are highly suitable for fishery protection, marine environment monitoring and other tasks [4]. In Japan, great advances have also been made in the research of wave-driven vessels. In 2008, a wave-driven catamaran named “Suntory Mermaid II” was made, and on March 16 of the same year, “Suntory Mermaid II” set off from Honolulu port in Hawaii, driven by the famous navigator Kenichi Horie. It relied only on ocean waves, wind and solar energy to support its progress. It took 110 days over 7800 kilometres to reach Japan, setting a historical record for human beings using green energy to cross the ocean [5,6]. In China, South China University of technology developed the “wave-driven unmanned ship” in 2017 and successfully tested it in the sea area near Hebao Island in the city of Zhuhai [7]. Harbin Engineering University developed the wave-driven unmanned catamaran “Wave Controller” in 2018, and successfully tested it in the sea near the city of Dalian in the next year [8,9].

The process of the wave-driven catamaran’s movement from static state to autonomous navigation driven by ocean waves is called “self-propulsion”. The self-propulsion performance, which involves the hydrodynamic and kinematic characteristics of catamaran, reflects the ability of wave-driven catamaran to absorb oceanic wave energy [10]. Analysis of its self-propulsion performance is crucial for the further development of the wave-driven catamaran. In this regard, in 2014, Li et al. [11] compiled a calculation program based on the quasi-steady hydrofoil theory to calculate the thrust-generation process of the oscillating hydrofoil in an oscillating cycle. The resistance increment in ship motion and waves is calculated by three-dimensional potential flow theory, and the interaction between hydrofoil and ship is considered. Firstly, the curve of the thrust varying with the assumed velocity of the hull is calculated, and the curve intersects with the resistance curve in the wave to estimate the achievable velocity. The research method is based on the potential flow theory, which is far from the actual viscous fluid environment.

In 2014, based on the research results on the ship wave-resistance theory and the propulsion in the hydrofoil wave, Feng [12] established a frequency-domain coupled hydrodynamic model based on the potential flow theory, and predicted the heave and pitch motion of the wave-energy-recovery-hydrofoil vessel in regular waves. On this basis, the prediction method of resistance of the vessel in regular waves is studied, and the radiation resistance, reflection resistance and thrust generated by hydrofoil are predicted.

In 2015, Zheng [13] used the CFD method to calculate the hydrodynamic and motion response of the ship hull in waves. Assuming that the ship was sailing at a constant speed, the heave-motion amplitude of the connection point between the oscillating hydrofoil and the ship hull was equivalent to that of the oscillating hydrofoil, and then the thrust generated by the vertical heave motion of the oscillating hydrofoil was calculated separately. However, the effect of hydrodynamic forces on the pitching motion of the hull was ignored.

In 2018, based on the bionic propulsion principle, Liu [14] calculated the hydrodynamic performance of the two-dimensional oscillating hydrofoil with tandem arrangement by solving the RANS equations, and analysed the influence of key design parameters, such as oscillating angle, wing spacing and foil shape on the propulsion performance, as well as the flow field disturbance between wings under different working conditions. However, it was assumed that the hull was fully responsive to the waves, which was not consistent with the actual situation.

In 2019, Y. Amini [15] studied the turbulent flow of a pitching hydrofoil near the free surface by using the Volume of Fluid method (VOF), and studied the influence of the depth underwater and oscillating frequency on the resistance coefficient and lift coefficient of the hydrofoil. The results showed that the resistance coefficient was positively correlated with the water depth and the oscillating frequency, and the lift coefficient was positively correlated with the oscillating frequency, which first increased and then decreased with the increase in water depth.

In 2019, Deng et al. [16] established the fluid-structure coupling dynamic model of a wave-driven ship and oscillating hydrofoil in a still water environment, and analysed the hydrodynamic performance of a wave-driven ship by programming the multi-body dynamic equation into CFD software FLUENT. The heave and pitch motions of the hull with different amplitudes and periods were simulated by the software, the propulsion force and torque generated by the oscillating hydrofoil were calculated to solve the multi-body dynamic equation, and the thrust generated by the hydrofoil under different amplitudes and periods was simulated and analysed.

In 2020, Chang et al. [17] established the multi-body dynamic model of the wave-driven ship based on Kane method, and derived the motion equation including the wave excitation force and the water force. By solving the motion equation, the response of the system was obtained. Based on the analysis of the components of the equation, the contribution and composition of the forces driving the ship motion were studied. It was found that the power of the wave-driven ship mainly came from the heave motion of the hull, and the thrust generated by the swing of the hydrofoil changes periodically. In the whole motion cycle, the thrust and the hull resistance cancel each other. However, the pitching motion of the ship in waves is ignored in the calculation, which is not consistent with the actual situation.

The above scholars provide many useful ideas for the study of the self-propulsion performance of wave-driven vessels, but there are some defects. Some scholars regard the hull and the oscillating hydrofoil of the wave-driven ship as independent elements, ignoring the interaction between them, and their dynamic models are therefore incomplete. When calculating the self-propulsion performance of wave-driven vessels, some scholars assumed that the vessel was moving at constant and uniform velocity, so that the real velocity of the wave-driven vessel was not included in their calculations. Some scholars assumed an ideal fluid environment for the wave-driven vessel, and their solutions based on potential flow theory led to significant errors in the numerical solutions of the hydrodynamics.

In order to obtain an accurate self-propulsion performance analysis to ensure optimal design for the wave-driven catamaran, the defects mentioned above must be overcome. Considering that the wave-driven catamaran has a multi-body dynamic structure, and its dynamic and kinematic responses driven by ocean waves are both unpredictable, a numerical method was proposed for simulating the self-propulsion process of the wave-driven catamaran to calculate its hydrodynamic response and sailing velocity in real time. In this method, the multi-body dynamics and computational fluid dynamics (CFD)were combined. Therefore, firstly, the multi-body dynamic model of the wave-driven catamaran was established, which took the hull and hydrofoils as a whole. Secondly, a two-phase numerical flow field which contained water and air was established. Lastly, the self-propulsion process of the wave-driven catamaran was simulated using an overlapping grid. In the simulation, the motion of the catamaran was a completely passive navigation driven by wave, which meant there was no artificial assumption with regard to its movement.

Compared with previous research, this method combined multi-body dynamics with CFD, and calculations were based on viscous fluid and the velocity of the wave-driven catamaran when propelled by wave energy.

2. Methods

2.1. Multi-Body Dynamic Model of the Wave-Driven Catamaran

In this paper, the hydrodynamic verification model of the “Wave Controller” wave-driven catamaran developed by Harbin Engineering University was taken as the research object. Its basic structure is shown in Figure 1.

The three-dimensional model of the hull and oscillating hydrofoils of the wave-driven catamaran were established. The catamaran had two pairs of hydrofoils, one pair at the bow and one pair at the stern. The pair of hydrofoils at the bow was set as hydrofoil #1, and that at the stern was set as hydrofoil #2. The hydrofoils and the hull were connected by rotating hinge. The model of the catamaran was simplified: the superstructure and the connecting plate between the hull and oscillating hydrofoils were ignored during the simulation. The simplified model is shown in Figure 2. The basic parameters of wave-driven catamaran are shown in

Table 1. Basic parameters of wave-driven catamaran.

Parameters	Numerical Value
Length of demihull	4 m
Breadth of demihull	0.5 m
Height of demihull	0.65 m
Distance between two demihull	2.5 m
Distance between two demihull	2.5 m
Draft of demihull	0.3 m
Number of oscillating hydrofoils	4
Section type of hydrofoil	NACA0012
Wingspan of hydrofoil	1.1 m
Chord length of hydrofoil	0.48 m
Longitudinal position of hydrofoil #1 shaft	0.5 m from the bow of demihull
Longitudinal position of hydrofoil #2 shaft	0.5 m from the stern of demihull
Vertical position of hydrofoil	0.75 m below waterline
Displacement	500 kg

.

The self-propulsion process of the wave-driven catamaran is shown in Figure 3. Here, the relative velocity between the fluid and the hydrofoil is defined as $V_f$ and the angle between $V_f$ and hydrofoil as $\theta$. When the hull is on the upper-wave-surface (Position A in Figure 3a), the hull moves up and pitches up, driving up the hydrofoils. At this time, $V_f$ points to the hydrofoil from above, as is shown in Figure 3b. Under the action of $V_f$, the hydrofoil generates a lift force $F_L$ that is perpendicular to $V_f$, and a drag force $F_D$ that is in the same direction as $V_f$. At this time, the hydrofoils rotate downward around its rotation axis to a certain angle $\phi$ under the action of $F_L$ and $F_D$. It can be easily derived that the thrust generated by each hydrofoil is $T=F_L\cdot \sin (\theta +\phi )-F_D\cdot \cos (\theta +\phi )$. Under the action of this thrust, the catamaran is propelled to sail.

When the hull is on the down-wave-surface (Position B in Figure 3a), the hull moves down and pitches down, driving down the two pairs of oscillating hydrofoils. At this time, the hydrofoil again rotates upward to a certain angle, and the thrust in the forward direction is again generated to propel the catamaran to sail [18,19,20].

The force analysis of the wave-driven catamaran overall is shown in Figure 4. The forward direction of the wave-driven catamaran is defined as direction x, and the direction perpendicular to the horizontal plane is direction z. The direction to the port side of the catamaran is defined as direction y, and the positive direction of the pitching angle of the hull and the angle of each hydrofoil is specified to be positive when rotating about the positive direction of y-axis. In this paper, the self-propulsion performance of the wave-driven catamaran was mainly direct navigation; only the motion on the xoz plane was studied. As shown in Figure 4, $G_{hull}$ is the gravity of the hull, $G_{foil}$ is the gravity of each hydrofoil, $R_{hull}$ is the sailing resistance of the hull, $F_{zh}$ represents the vertical hydrodynamic forces of hull, $T_1$ is the thrust generated by hydrofoil #1, $T_2$ is the thrust generated by hydrofoil #2, $F_{z1}$ represents the vertical hydrodynamic forces of hydrofoil #1, $F_{z2}$ represents the vertical hydrodynamic forces of hydrofoil #2, $M_{hull}$ is the pitching movement of the hull, $M_1$ is the pitching movement of hydrofoil #1, and $M_2$ is the pitching movement of hydrofoil #2. Among them, $G_{hull}$ and $G_{foil}$ are known quantities, and $R_{hull}$, $F_{z1}$, $T_1$, $T_2$, $M_1$, $M_2$ and $M_{hull}$ are solved by the CFD method.

The general equations describing the dynamic behaviour of the multi-body system are as follows [21,22,23]:

$$\begin{array}{l}M\ddot{q}+C_q{}^T\lambda =f\\ C_q(q,t)=0\end{array}$$

(1)

where $q$ is the generalized coordinate vector and $f$ is generalized force, including gravity and hydrodynamic force. $M$ is the square diagonal matrix of the inertia matrix of the rigid body and $C(q,t)=0$ is the constraint equation. If there are S constraints, then:

$$C(q,t)=\left[\begin{array}{c}{C}_1(q,t)\\ C_2(q,t)\\ ⋮\\ C_S(q,t)\end{array}\right]$$

(2)

By deriving the above equation, the velocity constraint equation is obtained:

$$C_q\dot{q}+C_t=0$$

(3)

where $C_q$ is the Jacobian matrix of the constraint equation. Then, the acceleration equation is obtained:

$$C_q\ddot{q}-\gamma =0$$

(4)

where

$$\gamma =-{(C_q\dot{q})}_q\dot{q}-2C_{qt}\dot{q}-C_{tt}$$

(5)

The constraint equation at any time $t$ can be expanded as follows:

$$C(q_i+\Delta q_i,t)=C(q_i,t)+C_q\Delta q_i+o(\Delta q_i{}^2)=0$$

(6)

It can be derived that:

$$C_q\Delta q_i=-C(q_i,t)$$

(7)

$\Delta q_i$ can be obtained from the formula, and the vector correction of generalized coordinates is obtained:

$$q_{i+1}=q_i+\Delta q_i$$

(8)

The new generalized coordinates of the multi-body system are thus obtained. The solution is completed when the convergence criteria of the following formula are satisfied:

$$(\left|\Delta q_i\right|\le {\epsilon }_1)\cap (\left|C(q_i,t)\right|\le {\epsilon }_2)$$

(9)

2.2. Numerical Flow Field

Considering that the research content of this paper involves flow of a high Reynolds number, and the literature [24,25,26] uses an RNG $k-\epsilon$ turbulence model to simulate the flow field around a hydrofoil, the simulation effect was effective. The hydrodynamic force of a hydrofoil was calculated accurately, as verified by previous test data. Therefore, in this paper, the RANS equations including the continuity and momentum equations were used to model the unsteady and incompressible wave field. The continuity equation is:

$$\nabla \cdot \overrightarrow{V}=0$$

(10)

where $\overrightarrow{V}=u\overrightarrow{i}+v\overrightarrow{j}+w\overrightarrow{k}$. $u$, $v$, $w$ are the components of fluid velocity in the $x$, $y$ and $z$ directions, respectively. Expand Equation (9) to obtain:

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0$$

(11)

The momentum equation is:

$$\begin{array}{c}\rho \left[\frac{\partial u}{\partial t}+\nabla \cdot \left(u\overrightarrow{V}\right)\right]=-\frac{\partial p}{\partial x}+\mu {\nabla }^{2}u+\nabla \cdot \left(-\rho \overline{u^{\prime }{\overrightarrow{V}}^{\prime }}\right)+\rho f_x\\ \rho \left[\frac{\partial v}{\partial t}+\nabla \cdot \left(v\overrightarrow{V}\right)\right]=-\frac{\partial p}{\partial y}+\mu {\nabla }^{2}v+\nabla \cdot \left(-\rho \overline{v^{\prime }{\overrightarrow{V}}^{\prime }}\right)+\rho f_y\\ \rho \left[\frac{\partial w}{\partial t}+\nabla \cdot \left(w\overrightarrow{V}\right)\right]=-\frac{\partial p}{\partial z}+\mu {\nabla }^{2}w+\nabla \cdot \left(-\rho \overline{w^{\prime }{\overrightarrow{V}}^{\prime }}\right)+\rho f_z\end{array}$$

(12)

where $\rho$ is the fluid density, $p$ is pressure, $\mu$ is the dynamic viscosity of fluid, While $f_x$, $f_y$, and $f_z$ are the component in $x$, $y$, $z$ directions, and $\overline{{u^{\prime }}_i{u^{\prime }}_j}$ is the Reynolds stress tensor.

The RNG $k-\epsilon$ model was used to model the wave turbulence, and the transport equations of the turbulent kinetic energy $k$ and the turbulent dissipation rate $\epsilon$ of are depicted as follows:

$$\begin{array}{cc}\hfill \frac{\partial \left(\rho k\right)}{\partial t}+\frac{\partial \left(\rho u_{j}k\right)}{\partial x_j}& =\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]\hfill \\ & +{\mu }_t\frac{\partial u_i}{\partial x_j}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)-\rho \epsilon \hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill \frac{\partial \left(\rho \epsilon \right)}{\partial t}+\frac{\partial \left(\rho u_j\epsilon \right)}{\partial x_j}& =\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right]-C_{\epsilon 2}\rho \frac{{\epsilon }^2}{k}\hfill \\ & +C_{\epsilon 1}\frac{\epsilon }{k}{\mu }_t\frac{\partial u_i}{\partial x_j}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\hfill \end{array}$$

(14)

The VOF method was used to simulate the free surface between water and air. The basic principle of the VOF method is to determine the interface between the two flows according to the proportion of the volume of fluid in a grid cell. The governing equation is:

$$\frac{\partial \gamma }{\partial t}+u\frac{\partial \gamma }{\partial x}+v\frac{\partial \gamma }{\partial y}+w\frac{\partial \gamma }{\partial z}=0$$

(15)

where $\gamma$ is the volume fraction of water in a cell. When $\gamma =0$, it means that the fluid in the region is air. When $\gamma =1$, it means that the fluid in the region is water. When $0<\gamma <1$, this refers to its free surface, that is the wave surface.

The cuboid computational domain was established in CFD solver STAR-CCM+ as flow field. In order to ensure that the catamaran had a long-enough sailing distance in the x direction, and that the bottom side and the left and right sides of the cuboid flow field did not interfere with the flow field around the catamaran, the cuboid flow field needed to be large. Therefore, the distance from the bottom and the left and right sides of the cuboid flow field to the centre of gravity of the catamaran needed to be greater than 2 times the length of the catamaran, and the distance from the front side to the centre of gravity of the catamaran was temporarily set as 10 times the length of the catamaran. The dimensions of the cuboid flow field were: 50 m long, 30 m wide and 10 m high. The front side of the flow field was 40 m away from the centre of gravity of the catamaran, the left and right sides were 15 m away from the centre of gravity of the catamaran, and the bottom side was 10 m away from the centre of gravity of the catamaran.

In the numerical simulation, if the size of cuboid flow field met the requirements of the sailing distance of the catamaran, the current size was adopted. If not, the cuboid flow field was lengthened until its size was adequate. The boundary conditions of the cuboid flow field were as follows: the front side, the top and bottom side were set as velocity inlets, the left and right sides were set as symmetrical planes, the back side was set as a pressure outlet, and the catamaran surface was set as a non-slip wall. The water surface was 10 m from the bottom of the cuboid flow field. The fluid media above and below the water surface were water and air, respectively. In the numerical simulation, in order to prevent wave reflection near the pressure outlet, it was necessary to set a damping wave elimination zone on the pressure outlet, in which a resistance term in the z direction was added to the momentum equation of the fluid. According to the method proposed in the literature [27], the resistance term was calculated according to the following formula:

$$S_z^d=-\rho \left(f_1+f_2\left|w\right|\right)\frac{e^k-1}{e-1}w$$

(16)

where

$$k={\left(\frac{x-x_s}{x_e-x_s}\right)}^n$$

(17)

where $f_1$ and $f_2$ are constants, $n$ is the attenuation coefficient of which the value is 2, $w$ is the velocity component of the fluid along the z direction, $x_s$ is the starting point of wave damping and $x_e$ is the end point. In this paper, $x_e-x_s$, the length of the wave elimination zone is 1.5 times of the wavelength.

In the flow field, the overlapping grid was used to simulate the passive navigation of the wave-driven catamaran. A cuboid was built outside each demihull and hydrofoil as an overlapping grid region. In order to ensure calculation accuracy, the grid size of the background region needed to be consistent with that of the surface of the overlapping region within the motion range of the catamaran. Therefore, for each overlapping region, when setting its size, it was necessary to ensure not only that the quality of the grid on the surface of each demihull and each hydrofoil was good, but that the grid on the 6 surfaces of the cuboid overlapping region was of the same size. Hence, the size of the overlapping region of each demihull was temporarily set to 5 m long, 1 m wide and 1.2 m high. The size of the overlapping region of each hydrofoil was temporarily set to 0.64 m long, 1.25 m wide and 0.3 m high. When generating the grid, for each demihull and each hydrofoil, if the grid on the 6 surfaces of the cuboid overlapping region were of the same size, the current size was adopted. If not, the size of overlapping region was enlarged until the grid size on 6 surfaces of one overlapping region was the same.

Boundary conditions and grid plan are shown in Figure 5.

The multi-body dynamic model and numerical flow field were combined together, by which the self-propulsion process of the wave-driven catamaran driven by regular waves were simulated [28]. The hydrodynamic performance and sailing velocity of the catamaran was calculated. The simulation process is shown in Figure 6. It should be noted here that in the literature [28] previously published by the authors, the research object was a flexible, connected wave-driven robot that absorbed wave energy by the heave motion of the float, while the research object of this paper was a wave-driven catamaran of which the hull and hydrofoil were rigidly connected through a rotating hinge, which absorbed more wave energy by the pitch motion of the hull. The dynamic model was different from the robot in literature [28].

The numerical uncertainty of the CFD numerical simulation was analysed, which process was divided into two parts: verification and validation. The verification part included grid-convergence verification and time-step convergence verification. When generating the grid, the height of the 1st grid on the catamaran surface needed to meet the requirement that the value of Y+ was approximately 30. When carrying out grid-convergence verification, 4 groups of different grids were selected, of which the minimum grid size was 0.012 m. The refinement ratio of 4 sets of grids was $r_G=\sqrt[6]{2}$. In the numerical simulation, the 5th order Stokes regular wave was adopted and the wave height was set as 0.25 m and the wavelength as 8 m. The encounter angle of the wave was 0°, which meant the catamaran sailed into the head wave. In order to ensure the calculation accuracy of each time step, the Courant number needed to meet the following requirement [29]:

$$C=V\Delta t/\Delta x\le 1$$

(18)

where $\Delta x$ is the minimum grid size, $V$ is the velocity of fluid, and $\Delta t$ is the time step of the program. For a wave in which the water depth is infinite, $V=1.25\sqrt{\lambda }$ and $\lambda$ is wavelength.

Because $\Delta x=0.012\mathrm{m}$, and $V=3.54\mathrm{m}/\mathrm{s}$ when $\lambda =8\mathrm{m}$, $\Delta t\le 0.0034\mathrm{s}$. Therefore, during grid-convergence verification, it was decided that the time step $\Delta t=0.0025\mathrm{s}$. At the initial time when $t=0\mathrm{s}$, the starting position of the wave was located 0.5m in front of the stem of the catamaran. At this time, the position of the wave crest should have been 5.5 m away from the front side of the flow field, of which the amplitude was defined as ${\zeta }_a$.

The average value of $T_S$, the total thrust of the hydrofoils and $R_{hull}$ (the sailing resistance of the hull) were calculated. Considering that the numerical simulation results of the catamaran’s motion in regular waves at the initial stage were not stable, the duration to calculate the average value of $T_S$ and $R_{hull}$ comprised 7 wave periods, from the 5th wave period to the 12nd wave period, after the catamaran had experienced the 1st four wave periods. During the grid-convergence verification, ${\zeta }_a$ was monitored. The results of numerical simulation are shown in

Table 2. Physical quantities under different grids.

Grid	Total Cells Number	${\mathit{T}}_{\mathit{S}}\mathbf{ }\mathbf{\left(}\mathbf{N}\mathbf{\right)}$	${\mathit{R}}_{\mathit{h}\mathit{u}\mathit{l}\mathit{l}}\mathbf{ }\mathbf{\left(}\mathbf{N}\mathbf{\right)}$	${\mathit{V}}_{\mathit{C}}\mathbf{ }\mathbf{(}\mathbf{m}\mathbf{/}\mathbf{s}\mathbf{)}$	${\mathit{\zeta }}_{\mathit{a}}\mathbf{ }\mathbf{\left(}\mathbf{m}\mathbf{\right)}$
Grid 1	4,055,474	54.36	35.02	0.773	0.112
Grid 2	6,826,876	60.87	39.22	0.727	0.121
Grid 3	10,768,546	65.95	42.49	0.692	0.126
Grid 4	15,987,062	65.48	42.19	0.694	0.126

.

It can be seen from Table 2 that with the increase in the number of grids, the values of each quantity tended to be stable, and the values of $T_S$, $R_{hull}$ and $V_C$ corresponding to grid 3 and grid 4 were similar. In addition, ${\zeta }_a$ gradually approached the theoretical value.

Since the numerical results corresponding to grid 3 and grid 4 were nearly equal, only grid 1, grid 2 and grid 3 were selected for grid convergence verification.

According to the literature [30], the calculation method of the convergence ratio is:

$$R_G=\frac{{\epsilon }_{G32}}{{\epsilon }_{G21}}$$

(19)

The order-of-accuracy is:

$$p_G=\frac{\ln \left({\epsilon }_{G32}/{\epsilon }_{G21}\right)}{\ln \left(r_G\right)}$$

(20)

The error is:

$${\delta }_{R{E}_G}^{*}=\frac{{\epsilon }_{G32}}{r_G^{p_G}-1}$$

(21)

The correction factor is:

$$C_G=\frac{r_G^{p_G}-1}{r_G^{p_{kest}}-1}$$

(22)

where $p_{kest}=2$

Error with correction factor ${\delta }_G^{*}$ is:

$${\delta }_G^{*}=C_G{\delta }_{R{E}_G}^{*}$$

(23)

The numerical uncertainty is:

$$U_G=\left\{\begin{array}{l}\left[9.6{\left(1-C_G\right)}^2+1.1\right]\left|{\delta }_{R{E}_G}^{*}\right|, \left|1-C_G\right|<0.125\\ \left(2\left|1-C_G\right|+1\right)\left|{\delta }_{R{E}_G}^{*}\right|, \left|1-C_G\right|\ge 0.125\end{array}\right.$$

(24)

$$U_{G_C}=\left\{\begin{array}{l}\left[2.4{\left(1-C_G\right)}^2+0.1\right]\left|{\delta }_{R{E}_G}^{*}\right|, \left|1-C_G\right|<0.25\\ \left|1-C_G\right|\left|{\delta }_{R{E}_G}^{*}\right|, \left|1-C_G\right|\ge 0.25\end{array}\right.$$

(25)

It should be noted here that in the above formula, subscript $G$ represents grid convergence verification.

Table 3. Grid convergence verification.

	${\mathit{R}}_{\mathit{G}}$	${\mathit{p}}_{\mathit{G}}$	${\mathit{C}}_{\mathit{G}}$	${\mathit{U}}_{\mathit{G}}$	${\mathit{\delta }}_{\mathit{G}}^{\mathbf{*}}$	${\mathit{U}}_{{\mathit{G}}_{\mathit{C}}}$
$T_S$	0.78	2.15	1.08	0.089	0.083	0.006
$R_{hull}$	0.78	2.15	1.08	0.317	0.295	0.022
$V_C$	0.76	2.38	1.22	0.230	-0.195	0.035
${\zeta }_a$	0.56	5.02	3.02	0.252	0.151	0.101

shows the results of grid-convergence uncertainty analysis.

It can be seen from Table 3 that the convergence factors of $T_S$, $R_{hull}$, $V_C$ and ${\zeta }_a$ all met the requirement of $0<R_G<1$, and the grid convergence was verified. Therefore grid 3 was finally adopted, and the total number of grids was 10,768,546.

Then, the time-step convergence was verified. Grid3 was selected when generating the grid, and $\Delta t$, the time steps were 0.005 s, 0.0025 s and 0.00125 s, respectively, which meant the refinement ratio of $\Delta t$ was $r_T=2$. The results of the numerical simulation are shown in

Table 4. Physical quantities under different time steps.

Time Step (s)	${\mathit{T}}_{\mathit{S}} \left(\mathbf{N}\right)$	${\mathit{R}}_{\mathit{h}\mathit{u}\mathit{l}\mathit{l}} \left(\mathbf{N}\right)$	${\mathit{V}}_{\mathit{C}} (\mathbf{m}/\mathbf{s})$	${\mathit{\zeta }}_{\mathit{a}} \left(\mathbf{m}\right)$
0.005	65.92	41.46	0.688	0.1251
0.0025	65.95	42.49	0.692	0.1256
0.00125	65.96	42.50	0.693	0.1257

.

It can be seen from Table 4 that with the decrease in time step, the values of each quantities tended to be stable. When $\Delta t$ were 0.0025 s and 0.00125 s, respectively, the values of $T_S$, $R_{hull}$, $V_C$ and ${\zeta }_a$ were very similar.

Table 5. Time step convergence verification.

	${\mathit{R}}_{\mathit{T}}$	${\mathit{p}}_{\mathit{T}}$	${\mathit{C}}_{\mathit{T}}$	${\mathit{U}}_{\mathit{T}}$	${\mathit{\delta }}_{\mathit{T}}^{*}$	${\mathit{U}}_{{\mathit{T}}_{\mathit{C}}}$
$T_S$	0.33	1.58	0.66	5.57 × 10−3	3.32 × 10−3	1.13 × 10−3
$R_{hull}$	0.33	1.58	0.66	5.57 × 10−3	3.32 × 10−3	1.13 × 10−3
$V_C$	0.25	2.00	1	3.67 × 10−4	3.33 × 10−4	3.33 × 10−5
${\zeta }_a$	0.20	2.32	1.33	5.53 × 10−5	3.33 × 10−5	1.10 × 10−5

shows the results of time-step convergence uncertainty analysis, and the subscript $T$ represents the time-step convergence verification.

It can be seen from Table 5 that the convergence factors of $T_S$, $R_{hull}$, $V_C$ and ${\zeta }_a$ all met the requirement of $0<R_G<1$, and the time step convergence was verified. Therefore, it was ultimately decided that the time step $\Delta t=0.0025\mathrm{s}$.

After convergence verification, the validation was carried out. Since only $V_C$ (the sailing velocity of catamaran) was measured in the test in the towing tank, only the validation of the numerical uncertainty of $V_C$ was carried out here. For details of the test in towing tank, please refer to Section 3.3. When the wave height is 0.25 m and wavelength is 8 m, $V_C=0.704\mathrm{m}/\mathrm{s}$. The error $E$ is:

$$E=D-S$$

(26)

where $D$ is $V_C$ in the test, which is 0.704 m/s, and $S$ is $V_C$ in the simulation, which is 0.692 m/s. Therefore, $E=0.012\mathrm{m}/\mathrm{s}$. $U_V$ is:

$$U_V^2=U_D^2+U_T^2+U_G^2$$

(27)

where $U_D$ is the uncertainty in the test, which is 2.5% in this paper. The validation of $V_C$ is shown in

Table 6. Validation of $V_C$.

	$\mathit{E}$	${\mathit{U}}_{\mathit{V}}$	$\mathit{E}(\mathbf{\%}\mathit{D})$	${\mathit{U}}_{\mathit{V}}(\mathit{\%}\mathit{D})$
$V_C$	0.012	0.231	1.7	32.6

.

.

It can be concluded from Table 6 that $\left|E\right|<U_V$, and validation is achieved at the $U_V$ level.

3. Results

The self-propulsion process of the wave-driven catamaran was simulated under regular wave and typical wave parameters, including wave height and wavelength, as shown in Table 3, in which each wave height corresponded to four wavelengths. It should be noted here that the wave parameters were selected according to the values in

Table 7. Wave parameters.

Wave Height (m)	Wavelength (m)
0.15	6, 8, 10, 12
0.2	6, 8, 10, 12
0.25	6, 8, 10, 12
0.3	6, 8, 10, 12

for comparison with the data of the test in the towing tank, because the towing tank where the authors work can produce waves with a wave height of 0.15~0.3 m and a wavelength of 6~12 m.

The dynamic and kinematic responses of the wave-driven catamaran were calculated during the simulation. The head of the regular wave was directly in front of the bow of the demihull at the initial moment when t = 0 s. The wave encounter angle of the catamaran was 0°, which meant the catamaran was sailing into the head wave.

For a wave in which the water depth is infinite, the wave period is calculated as follows [31,32]:

$$\tau =0.8\sqrt{\lambda }$$

(28)

where $\tau$ is the wave period and $\lambda$ is wavelength. When the wavelengths are 6 m, 8 m, 10 m and 12 m, the wave periods are 1.96 s, 2.26 s, 2.53 s and 2.77 s, respectively.

3.1. Analysis of the Hydrodynamics of the Wave-Driven Catamaran

Firstly, the hydrodynamic response of a wave-driven catamaran in the self-propulsion process was analysed, including the thrust generated by the oscillating hydrofoils and the sailing resistance of the hull. A self-propulsion diagram of a catamaran driven by wave energy as calculated using CFD is shown in Figure 7. The corresponding wave parameters in Figure 7 show a wave height of 0.25 m and a wavelength of 8 m. The animations of the whole self-propulsion process of the catamaran under different wave parameters are shown in supplementary materials, in which there are 4 videos. Each video shows the process of catamaran from static state to self-propulsion at time-varying velocity.

3.1.1. Analysis of the Thrust of Hydrofoils

Firstly, taking the wave parameters as being 0.25 m wave height and 8m wavelength as an example, the curves of the rotating angle of oscillating hydrofoils #1 and #2 and their respective thrust are shown in Figure 8a,b. It can be seen from Figure 8a that the hydrofoils #1 and #2 oscillated periodically, but the amplitude of the rotating angle of hydrofoil #1 was much larger than that of hydrofoil #2, and the oscillation of hydrofoil #2 was slightly behind that of hydrofoil #1. In the initial stage, and hydrofoils #1 and #2 both rotated downward and upward to the maximum angle. The maximum downward and upward angles of hydrofoil #1 were −26.6° and 34.8°, and those of hydrofoil #2 were −28.4° and 15°, respectively. As the self-propulsion process of a catamaran driven by waves tends to be stable, the angle of hydrofoil #1 was stable within −15° to 25°, and the angle of hydrofoil #2 was stable within −10° to 10°. The average values of $T_1$ and $T_2$, the thrusts generated by hydrofoils 1# and 2#, are shown in

Table 8. Average values of $T_1$ and $T_2$ when wave height is 0.25 m and wavelength is 8 m.

Hydrofoil	#1	#2
Thrust (N)	61.27	4.69

.

It can be seen from Figure 8a and Table 8 that the rotation amplitude of hydrofoil #1 was much larger than that of hydrofoil #2, so $T_1$ was much larger than $T_2$ numerically. The reason for this result is that the angle between the incident direction of the wave and the course of the catamaran was 0°, which meant the catamaran was sailing against the wave. Assuming that the vertical displacement at the centre of gravity of the ship is $z_o$, the pitching angle of hull is $\theta$, the pure and vertical distances between the rotation axis of each hydrofoil and the centre of gravity of the hull are $L$ and $z_v$, respectively, as is shown in Figure 8c, then $z_1$, the vertical displacement of hydrofoil #1, is:

$$z_1=z_o+\left[z_v-L\cdot \cos \left(\theta +\arccos \frac{z_v}{L}\right)\right]$$

(29)

Similarly, $z_2$, the vertical displacement of hydrofoil #2 is:

$$z_2=z_o-\left[z_v-L\cdot \cos \left(\theta +\arccos \frac{z_v}{L}\right)\right]$$

(30)

Because $z_v-L\cdot \cos \left(\theta +\arccos \frac{z_v}{L}\right)$ is numerically nonnegative, so as long as the catamaran is sailing against the waves, the numerical value of $z_1$ must be greater than that of $z_2$. As a result, the hydrodynamic moment of hydrofoil #1 was greater than that of hydrofoil #2, and the rotation amplitude of hydrofoil #1 was greater than that of hydrofoil #2, which made $T_1$ much larger than $T_2$ numerically. Similarly, when the catamaran sailed along the wave, which meant the angle between the incident direction of the wave and the course of the catamaran was 180°, $T_2$ was larger than $T_1$ numerically.

Moreover, when the wavelength was 8m, the wave period was 2.26 s. It can be seen from Figure 8b that in a wave period, there were two peaks in the curves of $T_1$ and $T_2$. This was because in each wave period, when the hull was on the upper-wave surface (Position A in Figure 3) and the lower-wave surface (Position B in Figure 3), the hydrofoils rotated downward and upward, respectively, which generated thrust twice in the forward direction.

Assume that the sum of the thrust forces generated by hydrofoil #1 and hydrofoil #2 is $T_S$, and $T_S=T_1+T_2$. The influence of wave parameters on $T_S$ was studied and the average values of $T_S$ under different wave parameters are shown in

Table 9. Average values of $T_S$ (Unit is N).

	Wavelength	6 m	8 m	10 m	12 m
Wave Height
0.15 m		43.81	35.65	21.48	14.83
0.20 m		54.21	44.10	26.58	18.35
0.25 m		81.06	65.95	39.75	27.44
0.30 m		92.62	75.36	45.42	31.36

.

The curve of the average value of $T_S$ changing with wave height and wavelength is shown in Figure 9a, and the curves of $T_S$ varying with time under different wave parameters are shown in Figure 9b,c. Among them, Figure 9b shows the variation curve of $T_S$ under the wave height of 0.25 m and different wavelengths, and Figure 9c shows the variation curve of $T_S$ under the wavelength of 8m and different wave heights.

It can be seen from Figure 9a–c that $T_S$ increases numerically with the increase in wave height and decreases numerically with the increase in wavelength. In addition, the period of the curve of $T_S$ is greatly affected by the wavelength and less affected by the wave height. In Figure 9b, the period of the curve of $T_S$ increases with the increase in wavelength. Under different wavelengths, the curve of $T_S$ has obvious phase differences. In Figure 9c, under different wave heights, the phase difference of $T_S$ curve is not obvious, and when $t\ge 20\mathrm{s}$, the phase difference between curves becomes obvious gradually. The reason for the phase difference in the curve of $T_S$ under different wave heights is that $V_C$, the sailing velocity of the catamaran, was different, which made the time taken for the catamaran to sail over a single wavelength distance also different. Therefore, there is a phase difference in the curve period. However, compared with the phase difference in Figure 9b, the phase difference in the curves of $T_S$ under different wave heights is still very small.

3.1.2. Analysis of Sailing Resistance of Hull

The influence of wave parameters on $R_{hull}$ was studied and the average values of $R_{hull}$ under different wave parameters are shown in

Table 10. Average values of $R_{hull}$ (Unit is N).

	Wavelength	6 m	8 m	10 m	12 m
Wave Height
0.15 m		30.96	22.97	13.24	5.18
0.20 m		34.34	25.47	14.68	5.75
0.25 m		57.29	42.49	24.50	9.59
0.30 m		65.46	48.55	27.99	10.96

.

The curve of the average value of $T_S$ changing with wave height and wavelength is shown in Figure 10a. It can be seen from Figure 10a that $R_{hull}$ increased numerically with the increase in wave height and decreased with the increase in wavelength, similarly to $T_S$. The variation laws of $R_{hull}$ under different wave parameters were compared respectively. The variation curves of $R_{hull}$ under different wavelengths when wave height was 0.25 m are shown in Figure 10b, and those under different wave heights when wavelength was 8 m are shown in Figure 10c.

It can be seen from Table 10 that $R_{hull}$ also varies periodically, as with $T_S$, and decreases numerically with the increase in wavelength, and the period of the curve of $R_{hull}$ increases with wavelength when the wave height is fixed. Moreover, it can be seen from Figure 10c that when the wavelength is fixed, $R_{hull}$ increases numerically with the increase in wave height, and the period of the curves of $R_{hull}$ are basically the same under the same wavelength. This is because when the wave height was fixed, with the decrease in the wavelength, $T_S$, the thrust generated by the hydrofoils, increased correspondingly, resulting in the increase in the velocity of catamaran, which made $R_{hull}$, the sailing resistance of the hull, also increase. Moreover, similarly to $T_S$, the period of the curve of $R_{hull}$ was greatly affected by the wavelength and less affected by the wave height. In Figure 10c, the phase difference between the curves of $R_{hull}$ under different wave heights is not obvious compared with Figure 10b.

In addition, taking the wave parameters of 0.25 m wave height and 8m wavelength as an example, the curves of $T_S$ and $R_{hull}$ are shown in Figure 10d. It can be seen from Figure 10d that the period of the curve of $T_S$ is almost half of that of $R_{hull}$. This is because the period of $R_{hull}$ was greatly affected by the wave period, and was almost consistent with the wave period. Moreover, comparing Figure 10d, Table 9 and Table 10, it can be seen that although the peak value of the curve of $R_{hull}$ is greater than that of $T_S$ in a single wave period, and the average value of $T_S$ is always greater than that of $R_{hull}$ in general.

3.2. Analysis of the Kinematic Response of the Wave-Driven Catamaran

3.2.1. Heave and Pitch Motion of the Wave-Driven Catamaran

The responses of the heave motion of the hull under different wave parameters were compared. The curves of heave motion under different wavelengths when the wave height was 0.25 m are shown in Figure 11a, and those under different wave heights when the wavelength was 8 m are shown in Figure 11b.

It can be seen from Figure 11a that when the wave height was fixed, the amplitude and period of heave motion increased with the increase in wavelength. This is because the vertical motion of hydrofoils #1 and #2 in the flow field made them subject to the force produced by the surrounding fluid, which could be divided into the longitudinal force along the x direction and the vertical force along the z direction. The longitudinal forces of each hydrofoil were $T_1$ and $T_2$, and the vertical forces, $F_{z1}$ and $F_{z2}$, had a damping effect on the heave motion of the hull, and increased with the increase in $T_1$ and $T_2$. Specifically, when the hull was on the upper-wave surface (Position A in Figure 3), the flow field exerted a vertical downward force on the hydrofoil to prevent the hull from moving up; when the hull was on the lower-wave surface (Position B in Figure 3), the flow field exerted a vertical upward force on the hydrofoil to prevent the hull from moving down. When the wave height was fixed, with the decrease in the wavelength, $T_1$ and $T_2$, the thrust generated by each hydrofoil, increased, and $F_{z1}$ and $F_{z2}$, the vertical force of each hydrofoil, also increased. Therefore, the larger the $T_S$, the greater the damping effect by $F_{z1}$ and $F_{z2}$ on the heave motion of the hull. When the wavelength was 6 m, $T_S$ was at its maximum, but the heave amplitude of the hull was at a minimum.

Still, when the wavelength was fixed, the influence of wave height on the heave amplitude of the hull was quite obvious, as shown in Figure 11b. It can be seen that when the wave height is 0.3 m, $T_S$ is the maximum, but the heave amplitude increases with the increase in wave height, which indicates that the effect of $T_1$ and $T_2$ on the heave amplitude was negligible at this time.

The responses in pitch motion of the hull under different wave parameters were compared. The curves of pitch motion under different wavelengths when the wave height was 0.25 m are shown in Figure 11c and those under different wave heights when the wavelength was 8 m are shown in Figure 11d.

It can be seen from Figure 11c that when the wave height is fixed, the amplitude and period of pitch motion increase with the increase in wavelength, as with heave motion. This is also because $F_{z1}$ and $F_{z2}$, the vertical forces of each hydrofoil produced by the surrounding fluid, had a damping effect on the pitch motion of the hull. When the hull was on the upper-wave-surface (Position A in Figure 3), the flow field produced a vertical downward force on the hydrofoil to prevent the hull from pitching up; when the hull was on the down-wave-surface (Position B in Figure 3), the flow field produced a vertical upward force on the hydrofoil to prevent the hull from pitching down. Therefore, the larger the $T_S$, the greater the damping effect by $F_{z1}$ and $F_{z2}$ on the pitch motion of the hull. When the wavelength was 6 m, $T_S$ was at a maximum, but the pitch amplitude of the hull was at a minimum.

When the wavelength was fixed, the effect of $F_{z1}$ and $F_{z2}$ on the pitch amplitude of the hull was far less than that of the wave height, as shown in Figure 11d. It can be seen that when the wave height is 0.3 m, $T_S$ is the maximum, but the heave amplitude still increases with the increase in wave height.

Moreover, similarly to $T_S$ and $R_{hull}$, the period of heave and pitch motion is greatly affected by wavelength and less affected by wave height. In Figure 11b,d, the phase difference between the curves of heave and pitch motion under different wave heights is not significant compared with Figure 11a,c.

3.2.2. Analysis of the Velocity of the Wave-Driven Catamaran

The variation law of $V_C$, the self-propelled velocity of the wave-driven catamaran, was studied. The average values of $V_C$ under different wave parameters are shown in

Table 11. Average values of $V_C$ (Unit is m/s).

	Wavelength	6 m	8 m	10 m	12 m
Wave Height
0.15 m		0.22	0.22	0.17	0.19
0.20 m		0.53	0.52	0.40	0.31
0.25 m		0.70	0.69	0.53	0.41
0.30 m		0.83	0.82	0.63	0.49

. It should be noted here that when calculating the average value of $V_C$, the previous acceleration stage is omitted; only the average value of $V_C$ when the curve is stable is calculated.

The curve of the average value of $V_C$ changing with wave height and wavelength is shown in Figure 12a, and the curves of $V_C$ varying with time under different wave parameters are shown in Figure 12b,c. Among them, Figure 12b shows the variation curves of $V_C$ under different wavelengths when the wave height is 0.25 m. Figure 12c shows the variation curves of $V_C$ under different wave heights when the wavelength is 8 m.

It can be seen from Figure 12a that $V_C$ increases numerically with the increase in wave height, and decreases numerically with the increase in wavelength. Moreover, with the decrease in wave height and the increase in wavelength, the decrease in the average value of $V_C$ is accelerated.

It can be seen from Figure 12b,c that when the catamaran was propelled by wave energy, $V_C$, the velocity of the catamaran increased rapidly in the initial stage of the self-propulsion process, then slowed down, and finally tended to be stable. The variation curve of $V_C$ shows a periodic fluctuation of small amplitude, and with the increase in wavelength, the fluctuation period of the curve of $V_C$ also increases. When the wavelength was fixed, $V_C$ increased numerically with the increase in wave height, and as with $T_S$ and $R_{hull}$, the period of the curve of $V_C$ was greatly affected by wavelength and less affected by wave height. In Figure 12c, under different wave heights, the phase difference between the curves of $V_C$ is not obvious.

In addition, by comparing Figure 12b with Figure 12c, it can be seen that the numerical changing amplitude of $V_C$ caused by the change in wavelength in Figure 12b is obviously smaller than that caused by the change in wave height in Figure 12c, which indicates that for the wave-driven catamaran in which the hull and hydrofoils were connected by hinges in this paper, the influence of wave height on the velocity of the catamaran driven by wave energy was obviously stronger than that of wavelength.

3.3. Self-Propulsion Test in Towing Tank

In order to verify the accuracy of the numerical results, it was necessary to carry out a self-propulsion test on the wave-driven catamaran prototype in regular waves. The test was carried out in the comprehensive test pool at Harbin Engineering University, and the test prototype was the hydrodynamic verification model of the wave-driven catamaran named “Wave Controller” developed by Harbin Engineering University. The towing tank was 50 m long, 30 m wide and 10m deep, as shown in Figure 13a. It is equipped with an X-Y carriage structure for the tow-testing of ship models. Wave-making equipment is installed at one side of the length direction of the tank. The wave parameters that can be generated by the wave-making equipment are a wave height of 0.15 m~0.3 m and a wavelength of 6 m~12 m. Current-making equipment is installed at the other side of the tank, and the maximum velocity of the water current is 0.6 m/s. In this paper, only the wave-making equipment was used; the current-making equipment and X-Y carriage structure were not used. The wave-making equipment is shown in Figure 13b.

Because the data regarding sailing velocity comprise the most important physical quantity for the self-propulsion performance of wave-driven catamarans, and are also the most intuitive expression of the self-propulsion process, the numerical simulation results were verified directly by measuring the data of the sailing velocity of the prototype in the test. The prototype was equipped with indoor positioning equipment, which could measure the position and velocity of the catamaran in real time, and transmit the velocity data to the monitoring terminal beside the tank every 0.5 s. The indoor positioning equipment consisted of 4 base stations beside the tank and 1 receiver installed on the catamaran. The 4 base stations were respectively arranged at the 4 vertices of the cube towing pool, as is shown in Figure 14a. The validity of the numerical method was verified by comparing the $V_C$ in the test and that in the simulation. Figure 14b,c show the test scene in top and side view, respectively.

Taking the working condition of 8m wavelength as example, the comparison between variation curves of $V_C$ in the test and that in the simulation under the wave heights of 0.15 m, 0.2 m and 0.25 m are shown in Figure 15a–c. The average values of $V_C$ in the test and the simulation are shown in

Table 12. Average values of $V_C$ in test and simulation (Wavelength is 8 m).

Wave Height	${\mathit{V}}_{\mathit{C}} \mathbf{in} \mathbf{Test} (\mathbf{m}/\mathbf{s})$	${\mathit{V}}_{\mathit{C}} \mathbf{in} \mathbf{Simulation} (\mathbf{m}/\mathbf{s})$
0.15m	0.201	0.216
0.2m	0.498	0.519
0.25m	0.704	0.692

. Similarly to $V_C$ in the numerical simulation, when calculating the average value of $V_C$, the previous acceleration stage is omitted; only the average value of $V_C$ when the curve is stable is calculated.

It can be seen from Figure 15 and Table 12 that the in the test was reasonably close to that $V_C$ in the simulation, and the average value of $V_C$ in the test under the wave heights of 0.15 m and 0.2 m was less than that in simulation. When the wave height reached 0.25 m, the average value of $V_C$ in the test exceeded that in the numerical simulation.

Taking the working condition of 0.25 m wavelength as an example, the comparison between the variation curves of $V_C$ in the test and those in simulation under the wavelengths of 6m, 10m and 12m are shown in Figure 16a–c. The average values of $V_C$ in the test and in the simulation are shown in

Table 13. Average values of $V_C$ in test and simulation (Wave height is 0.25 m).

Wavelength	${\mathit{V}}_{\mathit{C}} \mathbf{in} \mathbf{Test} (\mathbf{m}/\mathbf{s})$	${\mathit{V}}_{\mathit{C}} \mathbf{in} \mathbf{Simulation} (\mathbf{m}/\mathbf{s})$
6m	0.709	0.699
8m	0.704	0.692
10m	0.544	0.532
12m	0.386	0.411

.

It can be seen from Figure 16 and Table 13 that the $V_C$ in the test was reasonably close to that in the simulation, and the average value of $V_C$ in the test under the wavelength of 6 m, 8 m and 10 m was larger than that in the simulation. When the wavelength reached 12 m, the average value of $V_C$ in the test was less than that in numerical simulation.

In addition, the variation trend of the curve of $V_C$ measured in the test was similar to that in the numerical simulation, which increased from the initial time and finally tended to be stable, and the curve of $V_C$ in the test increased from 0 m/s to the final stable stage.

However, the curve of $V_C$ in the test did not show the obvious periodic oscillation seen in the numerical simulation. The main reason for this phenomenon is that there is a filtering algorithm in the processor of the indoor positioning equipment, which greatly reduced the oscillation of the curve of $V_C$ during the test in towing tank. Moreover, the control system of the prototype itself transmitted the velocity data to the monitoring terminal beside the tank every 0.5 s, while in the numerical simulation, the time step of the program was 0.0025 s, which meant the velocity data were recorded every 0.0025 s. The significant difference between these two time steps led to the data points in the numerical simulation being much denser, showing the periodic oscillation characteristics of the curve of $V_C$. The time step in the pool test was significantly greater, making it difficult to show the detailed characteristics of the curve of $V_C$. As a result, the periodic oscillation did not appear in the curve of $V_C$ in the test.

In addition, it should be noted here that in the CFD numerical simulation, the 3D model of the catamaran only retained the hull and hydrofoil; the connecting rod between them was ignored. In the self-propulsion test in the towing tank, the hull of the prototype was connected with the hydrofoil through a cylindrical stainless-steel rod of which the diameter was 20 mm. Through comparison of the data in Table 12 and Table 13, it can be concluded that the difference between the results in the test and those in the simulation is small, and therefore the resistance generated by the cylindrical connecting rod can be ignored.

4. Conclusions

In this paper, based on the hydrodynamic model “Wave Controller” wave-driven catamaran developed by Harbin Engineering University, a numerical calculation method of the self-propulsion process of the wave-driven catamaran coupled with multi-body dynamics and CFD was developed. Firstly, the multi-body dynamic model of a wave-driven catamaran was established. Secondly, the numerical flow field was established. Finally, the self-propelled motion of the wave-driven catamaran in regular waves was realized using overlapping grids. The self-propulsion process of the wave-driven catamaran was numerically simulated, and the hydrodynamic response and motion response of the wave-driven catamaran were calculated. In particular, the thrust generated by oscillating hydrofoils, the sailing resistance, the heave and pitch motion of the hull and the velocity of the catamaran were comprehensively analysed. In addition, the accuracy of the numerical results was verified by a self-propulsion test in a towing tank. The conclusions are as follows:

(1)

In this paper, a numerical method for calculating the self-propulsion process of a wave-driven catamaran was proposed. In contrast with previous research, the self-propelled velocity of the catamaran was calculated by this method. By analysing the thrust generated by the oscillating hydrofoils, it can be concluded that the thrust of the hydrofoils at the bow of the catamaran was much greater than that at the stern, because the heave amplitude of the front hydrofoil was greater than that of the rear hydrofoil. In addition, the thrust of the hydrofoils and the sailing resistance of the hull were both positively correlated with wave height and negatively correlated with wavelength.

(2)

By analysing the heave and pitch motions of the hull under typical wave parameters, it can be concluded that when the wave height was fixed, the heave and pitch amplitudes of the hull increased with the increased of wavelength. This is because the vertical force on the oscillating hydrofoils had a damping effect on the heave and pitch motions of the hull, and the larger the thrust generated by the hydrofoils, the larger the vertical force on the hydrofoils. As a result, the heave and pitch motion of the hull was weakened. This also reflected the strong interaction between each monomer in the multi-body system.

(3)

By analysing the velocity of the catamaran under typical wave parameters, it can be concluded that velocity of the catamaran was positively correlated with wave height and negatively correlated with wavelength. When the wave height was 0.3m and the wavelength was 6m, the velocity of the catamaran reached 0.83 m/s. As for the wave-driven catamaran studied in this paper, the influence of wave height on its velocity was obviously stronger than that of wavelength. In addition, compared with the curves of velocity of the wave-driven robot in the literature [28], the wave-driven catamaran on which the hull and hydrofoil were rigidly connected through a rotating hinge was able to sail at a more stable velocity, the fluctuation of the curve of sailing velocity changing with time being much smaller.

On this basis, further research on the influence of structural parameters on the self-propulsion performance of the wave-driven catamaran, such as the depth of hydrofoils, the longitudinal position of hydrofoils and the distance between two demihulls, as well as research on self-propulsion performance in different wave directions and irregular waves, will be carried out in the next stage.