Numerical study of a novel small waterplane area USV advancing in calm water and in waves using the higher-order Rankine source method

Abstract

The wave loads and motion responses of an Unmanned Surface Vehicle (USV) in water will directly influence the powering requirement, energy supply design and function of the installed sensors. In the present article, to investigate the wave loads and motions response characteristics of a USV advancing on the free surface, a Higher-Order Rankine Source (HORS) method is proposed. During the discretization of the boundary elements, bi-quadratic B-splines are applied to distribute the velocity potentials on the body surface and free surface. Based on the proposed HORS method, steady state simulation and time-domain simulation approaches are both used for predicting the wave loads and motion responses of a novel Small Waterplane Area (SWA) USV in calm water and regular waves, respectively. The predicted wave loads in calm water and the predicted heave and pitch in heading wave are compared with corresponding model test data. The good agreement found indicates the validity of the proposed HORS method. Finally, the heave and pitch motion responses and the wave components of the USV including diffraction and radiation forces and Froude−Krylov forces are investigated when the SWA USV is running against heading wave under different conditions of wavelength, wave steepness and advancing velocity.

KEYWORDS:

• Wave loads
• quadratic B-splines
• small waterplane area USV
• higher-order Rankine source method
• regular waves
• model tests

1. Introduction

Currently, autonomous marine robots like Unmanned Surface Vehicle (USV) have a wide range of applications that can realize marine environment monitoring, marine resource sampling, ocean rescuing, etc. (Li et al., 2022). However, during the implementation of time-consuming tasks under high sea states, USV performance will decrease and the endurance is greatly shortened because, with limited energy stored in the battery device, more power is required to overcome the acting wave loads. Also, its speed may be greatly reduced because, as the wave resistance becomes large, the speed will have to be reduced to cater for the limited power available. Consequently, the USV may not be able to return after using up all its energy or may not reach the required speed to complete its tasks. Furthermore, too large motion responses under wave excitations may cause the malfunctioning of USVs’ sensors, greatly reducing the task capability of these marine robots. There are mainly two aspects to solving this problem: to improve the hull-form design of USVs and to predict accurately the performance of USVs’ advance in the water.

Numerous innovations have been introduced into USVs such as refinements to hull-form components such as stern flaps, alternative propulsors such as waterjets, ducted propellers and tunnel drives (Gong et al., 2022; Park et al., 2021; Song et al., 2021; Sun et al., 2016; Xia et al., 2020) and, most importantly, improvements in the design of the hull form itself with the continuous development of new hull types such as trimarans, slender high-speed catamarans, Small-Waterplane-Area Twin Hull ships, etc.

Because of the limits and restrictions on hull form or its component options to the very limited amount of existing data, these new hull forms or technology have no existing data for reference. The only way available to designers of novel hulls has been to resort to more complex boutique design tools peculiar to each type of hull. An efficient numerical tool capable of defining hull geometry fast and evaluating the resistance and seakeeping properties of a variety of non-conventional hulls would greatly enhance early-stage designer creativity and efficiency, i.e. the designer’s ability to evolve effective ship concepts (Xiang and Xiang, 2021).

Initially, low fidelity two-dimensional (2D) methods were proposed for resistance and seakeeping simulations (Stern et al., 2006). Later, more generalized high-fidelity codes, such as the commercial codes ANSYS^® Fluent^®, Star CCM+™, open source codes, Open-FOAM^® and the inhouse codes CFDShip-Iowa have been all widely validated for predicting the seakeeping performance of floating vessels in Begovic et al. (2015), Feng et al. (2019), Poundra et al. (2017), Xiang and Soares (2020, 2022), Zhai et al. (2021) and Deng et al. (2022).

However, as Stern et al. (2008) mentioned, hydrodynamic analysis tools and ship motion simulators based on three-dimensional (3D) potential-flow can potentially offer a more computationally efficient solver compared with the viscous flow CFD codes mentioned above. The potential theory is generally implemented by the Boundary Element Method (BEM) approach, entailing the solution of a boundary integral formulation on the body and free surface boundaries for the governing Laplace equation. Meanwhile, both linear and nonlinear formulations for the wave flow and body motions can be implemented. Hong et al. (2016) calculated the hydrodynamic forces and motions of the Wigley III hull and the S-175 containership advancing in waves with various forward speeds by using 3D frequency domain potential flow theories. The 3D translating and pulsating (3DTP) source panel method is employed to solve the radiation and diffraction potential problems, in which a semi-analytical scheme is established to get more robust results from the integration of Green’s function over panels. The proposed method is found to be of satisfactory accuracy and efficiency. Yu, Zhu, Xu, Huang, and Hong (2021) solved the large amplitude motion problems in the time domain when the ship is advancing with constant forward speed in waves by combining a transient Green’s function method and a Rankine source method. A good agreement can be obtained by comparing simulated results to those of the frequency domain 3D panel method (NEWDRIFT) of NTUA-SDL and available experimental data. Lee et al. (2021) presents a comparative study of the motion responses and added resistance of four representative types of ship−LNG carrier, tanker, containership and bulk carrier−by several numerical methods (asymptotic formulae, 2D strip theory, the 3D panel method and CFD). The comparison results validated the accuracy and reliability of each analysis technique, and related characteristics and limitations were investigated. Zhang et al. (2020) numerically simulated the S-175 containership advancing in calm water and in regular oblique waves using a double-body potential and trailing vortex sheets combined method. A systematic analysis demonstrated that the transverse wave drift force attained considerable amplitudes, which increased rapidly at larger drift angles in short beam waves.

In the present article, to increase the discretization accuracy further, a bi-quadratic B-splines based Higher-Order Rankine Source (HORS) method is proposed for numerically studying the wave loads on a novel Small Waterplane Area (SWA) USV in regular waves and calm water, respectively. The predicted wave loads on the SWA USV in calm water under different Froude numbers and its motion response, including heave and pitch in heading wave, is compared with corresponding experimental data to check the validity of the proposed HORS method. Finally, the wave radiation and diffraction force, Froude−Krylov (F-K) force and the resultant total wave loads are investigated in the time domain when the SWA USV is running at different advancing velocities against a fixed regular wave and running against different waves at a constant advancing velocity, respectively.

2. Mathematical formulations

2.1. Governing equations

The 3D Rankine source method is based on the ideal flow assumption that the fluid is inviscid, irrotational and incompressible. Assume a USV with a body fixed system $xyz$ is moving in the inertial system $X_0{Y}_0{Z}_0$ with a constant advancing velocity of $\overrightarrow{U}\left(\overrightarrow{x}\right)$ as shown in Figure 1; there exists a total disturbance velocity potential $\mathrm{\Psi }\left(\overrightarrow{x}\right)$ meeting the requirement of the Laplace equation, (1) ${\mathrm{\nabla }}^2\mathrm{\Psi }=0$(1) By applying Green’s theorem to the Laplace equation in Equation (1), the Boundary Integral Equation (BIE) is obtained: (2) $\begin{array}{rl}\alpha \mathrm{\Psi }\left(\overrightarrow{x}\right)& =\int \int {}_S\int \int \frac{\mathrm{\partial }\mathrm{\Psi }\left(\overrightarrow{x}\right)}{\mathrm{\partial }n}G\left({\overrightarrow{x}}^{\mathrm{\prime }},\overrightarrow{x}\right)\text{d}S\left(x^{\mathrm{\prime }}\right)\\ & -\int \int \int \int {}_S\frac{\mathrm{\partial }G\left({\overrightarrow{x}}^{\mathrm{\prime }},\overrightarrow{x}\right)}{\mathrm{\partial }n}\mathrm{\Psi }\left(\overrightarrow{x}\right)\text{d}S\left(x^{\mathrm{\prime }}\right)\end{array}$(2) where $G\left({\overrightarrow{x}}^{\mathrm{\prime }},\overrightarrow{x}\right)$ is the Rankine source potential, $G\left({\overrightarrow{x}}^{\mathrm{\prime }},\overrightarrow{x}\right)=1/\left|{\overrightarrow{x}}^{\mathrm{\prime }}-\overrightarrow{x}\right|$, $S$ denotes the boundary including the free surface boundary, $S_f$, and the ship surface boundary, $S_b$, $\mathrm{\Psi }\left(\overrightarrow{x}\right)$ is the total velocity potential on the defined nonlinear free surface profile of $z=\eta$, and $\alpha$ is the solid angle of a field point. In Equation (2), two critical unknowns of normal velocity $\mathrm{\partial }\mathrm{\Psi }\left(\overrightarrow{x}\right)/\mathrm{\partial }n$ on $S_f$ and $\mathrm{\Psi }\left(\overrightarrow{x}\right)$ on $S_b$ need to be solved for the evaluation of $\mathrm{\Psi }\left(\overrightarrow{x}\right)$ on the whole boundary $S$. The two integrals in Equation (2) can be interpreted as the influence of a source distribution of strength, $\mathrm{\partial }\mathrm{\Psi }\left(\overrightarrow{x}\right)/\mathrm{\partial }n$, and a normal dipole distribution of strength $\mathrm{\Psi }\left(\overrightarrow{x}\right)$ on $S$, respectively. It is clear from Equation (2) that $\mathrm{\Psi }\left(\overrightarrow{x}\right)$ inside $S$ is completely determined provided $\mathrm{\partial }\mathrm{\Psi }\left(\overrightarrow{x}\right)/\mathrm{\partial }n$ on $S_f$ and $\mathrm{\Psi }\left(\overrightarrow{x}\right)$ on $S_b$ are known.

Figure 1. Two coordinate systems.

In general, the BIE needs to be solved numerically. To do that, two boundaries − $S_f$ and $S_b$ − are subdivided into their own sets of small boundary elements: (3) $S=S_f+S_b=\sum _{i=1}^{N_f}\text{B}\text{E}{\text{F}}_i+\sum _{i=1}^{N_b}\text{B}\text{E}{\text{B}}_i$(3) where $\text{BEF}$ and $\text{BEB}$ represent $N_f$ boundary elements on $S_f$ and $N_b$ boundary elements on $S_b$, respectively. The discretized integral equation can be written in the form (4) $\begin{array}{rl}\alpha \mathrm{\Psi }\left(\overrightarrow{x}\right)& =\sum _{i=1}^{N_b}\{\int {\int }_{\text{BE}{\text{B}}_i}\frac{\mathrm{\partial }\mathrm{\Psi }\left(\overrightarrow{x}\right)}{\mathrm{\partial }n}G({\overrightarrow{x}}^{\mathrm{\prime }},\overrightarrow{x}\left)\text{d}S\right(x^{\mathrm{\prime }})\\ & -\int {\int }_{\text{BE}{\text{B}}_i}\frac{\mathrm{\partial }G\left({\overrightarrow{x}}^{\mathrm{\prime }},\overrightarrow{x}\right)}{\mathrm{\partial }n}\mathrm{\Psi }\left(\overrightarrow{x}\right)\text{d}S\left(x^{\mathrm{\prime }}\right)\}\\ & +\sum _{i=1}^{N_f}\{\int {\int }_{\text{BE}{\text{F}}_i}\frac{\mathrm{\partial }\mathrm{\Psi }\left(\overrightarrow{x}\right)}{\mathrm{\partial }n}G({\overrightarrow{x}}^{\mathrm{\prime }},\overrightarrow{x}\left)\text{d}S\right(x^{\mathrm{\prime }})\\ & -\int {\int }_{\text{BE}{\text{F}}_i}\frac{\mathrm{\partial }G\left({\overrightarrow{x}}^{\mathrm{\prime }},\overrightarrow{x}\right)}{\mathrm{\partial }n}\mathrm{\Psi }\left(\overrightarrow{x}\right)\text{d}S\left(x^{\mathrm{\prime }}\right)\}\end{array}$(4) In general, the linearization is necessary to solve the total velocity potential, $\mathrm{\Psi }\left(\overrightarrow{x}\right)$. In present article, the aspiration model proposed by Kring (1994) is used for linearization. The aspiration model not only allows a specified normal flux through the hull and imposes symmetry about the $z=0$ plane, as in the double-body model, but also produces a particular transom treatment for ships with complex sterns that cannot be modelled well using double-body linearization or Neumann−Kelvin linearization. By linearization, the total disturbance velocity potential $\mathrm{\Psi }\left(\overrightarrow{x}\right)$ is composed of four components: $\mathrm{\Phi }\left(\overrightarrow{x}\right)$, the basis flow component; ${\mathrm{\Phi }}_l\left(\overrightarrow{x}\right)$, the local flow component; ${\mathrm{\Phi }}_m\left(\overrightarrow{x}\right)$, the memory flow component; and the incident wave flow component, ${\mathrm{\Phi }}_i\left(\overrightarrow{x}\right)$. These are expressed by Equation (5): (5) $\mathrm{\Psi }\left(\overrightarrow{x}\right)=\mathrm{\Phi }\left(\overrightarrow{x}\right)+{\mathrm{\Phi }}_l\left(\overrightarrow{x}\right)+{\mathrm{\Phi }}_m\left(\overrightarrow{x}\right)+{\mathrm{\Phi }}_i\left(\overrightarrow{x}\right)$(5) The basis velocity potential $\mathrm{\Phi }\left(\overrightarrow{x}\right)$ is assumed to be the dominant one, which provides the linearization. The local perturbation potential ${\mathrm{\Phi }}_l\left(\overrightarrow{x}\right)$ indicates the instantaneous fluid response, which provides the force for unsteady motion. The memory perturbation potential ${\mathrm{\Phi }}_m\left(\overrightarrow{x}\right)$ represents the solution for memory effects that arise from the propagation of wave energy away from the ship such as the steady, radiated and scattered wave patterns. Therefore, the solutions for the sum of the local and memory flows will be the wave radiation and scattering forces. The incident wave potential, ${\mathrm{\Phi }}_i\left(\overrightarrow{x}\right)$, is expressed by Equation (6) for the deep-water condition: (6) $\begin{array}{rl}{\mathrm{\Phi }}_i\left(\overrightarrow{x}\right)& =A\frac{g}{{\omega }_i}{e}^{k_{i}z}\sin \left\{k_i\right(x\cos \text{β}+y\sin \text{β})\\ & -\left[{\omega }_i-k_i(u\cos \text{β}+v\sin \text{β})\right]t\}\end{array}$(6) where ${\omega }_i$ is the frequency of the incident wave, $\text{β}$ is the wave encounter angle between the incident wave and the advancing ship, $A$ is the wave amplitude and $k_i$ is the wave number of the incident wave for the deep-water condition via $k_i={\omega }_i^2/g$.

2.2. Free surface boundary condition

The basis velocity potential $\mathrm{\Phi }\left(\overrightarrow{x}\right)$ also needs to fit the combined kinetic and dynamic boundary conditions on the free surface, which are represented by Equations (7) and (8): (7) $\begin{array}{rl}& \frac{\mathrm{\partial }\eta }{\mathrm{\partial }\text{t}}-\left(\overrightarrow{U}\left(\overrightarrow{x}\right)-\mathrm{\nabla }\mathrm{\Phi }\left(\overrightarrow{x}\right)\right)\cdot \mathrm{\nabla }\left(\eta \right)\\ & =\frac{\mathrm{\partial }}{\mathrm{\partial }z}\left(\mathrm{\Phi }\left(\overrightarrow{x}\right)+{\mathrm{\Phi }}_l\left(\overrightarrow{x}\right)+{\mathrm{\Phi }}_m\left(\overrightarrow{x}\right)\right)\\ & \text{on }z=\eta (x,y,t)\end{array}$(7) (8) $\begin{array}{rl}& \left(\frac{\mathrm{\partial }}{\mathrm{\partial }\text{t}}-\left(\overrightarrow{U}\right(\overrightarrow{x})-\mathrm{\nabla }\mathrm{\Phi }(\overrightarrow{x}\left)\right)\cdot \mathrm{\nabla }\right)\left({\mathrm{\Phi }}_l\left(\overrightarrow{x}\right)+{\mathrm{\Phi }}_m\left(\overrightarrow{x}\right)\right)\\ & =-g\eta +\left(\overrightarrow{U}\left(\overrightarrow{x}\right)\cdot \mathrm{\nabla }\mathrm{\Phi }\left(\overrightarrow{x}\right)-\frac{1}{2}\mathrm{\nabla }\mathrm{\Phi }\left(\overrightarrow{x}\right)\cdot \mathrm{\nabla }\mathrm{\Phi }\left(\overrightarrow{x}\right)\right)\\ & \text{on }z=\eta (x,y,t)\end{array}$(8) By applying a Taylor expansion for small disturbance about the $z=0$ plane, noting that $\mathrm{\partial }\mathrm{\Phi }\left(\overrightarrow{x}\right)/\mathrm{\partial }z=0$ and retaining linear terms, the free-surface conditions now become Equations (9) and (10): (9) $\begin{array}{rl}& \frac{\mathrm{\partial }\eta }{\mathrm{\partial }\text{t}}-\left(\overrightarrow{U}\right(\overrightarrow{x})-\mathrm{\nabla }\mathrm{\Phi }(\overrightarrow{x}\left)\right)\cdot \mathrm{\nabla }\eta \\ & =\frac{{\mathrm{\partial }}^2\mathrm{\Phi }\left(\overrightarrow{x}\right)}{\mathrm{\partial }{z}^2}\eta +\frac{\mathrm{\partial }{\mathrm{\Phi }}_l\left(\overrightarrow{x}\right)}{\mathrm{\partial }z}+\frac{\theta {\mathrm{\Phi }}_m\left(\overrightarrow{x}\right)}{\mathrm{\partial }z}\text{on }z=0\end{array}$(9) (10) $\begin{array}{rl}& \left(\frac{\mathrm{\partial }}{\mathrm{\partial }\text{t}}-\left(\overrightarrow{U}\right(\overrightarrow{x})-\mathrm{\nabla }\mathrm{\Phi }(\overrightarrow{x}\left)\right)\cdot \mathrm{\nabla }\right)\left({\mathrm{\Phi }}_l\left(\overrightarrow{x}\right)+{\mathrm{\Phi }}_m\left(\overrightarrow{x}\right)\right)\\ & =-g\eta +\left(\overrightarrow{U}\left(\overrightarrow{x}\right)\cdot \mathrm{\nabla }\mathrm{\Phi }\left(\overrightarrow{x}\right)-\frac{1}{2}\mathrm{\nabla }\mathrm{\Phi }\left(\overrightarrow{x}\right)\cdot \mathrm{\nabla }\mathrm{\Phi }\left(\overrightarrow{x}\right)\right)\\ & \text{on}z=0\end{array}$(10) To deal with the treatment of a transom, Kutta conditions are used: (11) $\begin{array}{rl}& \frac{\mathrm{\partial }{\mathrm{\Phi }}_m\left(\overrightarrow{x}\right)}{\mathrm{\partial }\text{z}}=-\frac{\mathrm{\partial }{\eta }_T}{\mathrm{\partial }\text{t}}+\left(\overrightarrow{U}\right(\overrightarrow{x})-\mathrm{\nabla }\mathrm{\Phi }(\overrightarrow{x}\left)\right)\cdot \mathrm{\nabla }{\eta }_T\\ & +\frac{{\mathrm{\partial }}^2\mathrm{\Phi }\left(\overrightarrow{x}\right)}{\mathrm{\partial }{z}^2}\left({\eta }_T\right)+\frac{\mathrm{\partial }{\mathrm{\Phi }}_l\left(\overrightarrow{x}\right)}{\mathrm{\partial }\text{z}}\\ & \text{on }(x,y,z)=\left(x_T,y_T,0\right)\end{array}$(11) (12) $\begin{array}{rl}& \left(\frac{\mathrm{\partial }}{\mathrm{\partial }\text{t}}-\left(\overrightarrow{U}\right(\overrightarrow{x})-\mathrm{\nabla }\mathrm{\Phi }(\overrightarrow{x}\left)\right)\cdot \mathrm{\nabla }\right)\left({\mathrm{\Phi }}_l\left(\overrightarrow{x}\right)+{\mathrm{\Phi }}_m\left(\overrightarrow{x}\right)\right)\\ & =-g{\eta }_T+\left(\overrightarrow{U}\left(\overrightarrow{x}\right)\cdot \mathrm{\nabla }\mathrm{\Phi }\left(\overrightarrow{x}\right)-\frac{1}{2}\mathrm{\nabla }\mathrm{\Phi }\left(\overrightarrow{x}\right)\cdot \mathrm{\nabla }\mathrm{\Phi }\left(\overrightarrow{x}\right)\right)\\ & \text{ on }(x,y,z)=\left(x_T,y_T,0\right)\end{array}$(12) For the local flow, the linearized free surface condition is (13) ${\mathrm{\Phi }}_l\left(\overrightarrow{x}\right)=0\text{on}z=0$(13) The local potential can be further decomposed into the following: (14) ${\mathrm{\Phi }}_{lk}\left(\overrightarrow{x}\right)=N_k\left(\overrightarrow{x}\right){\dot{\eta }}_k\left(t\right)+M_k\left(\overrightarrow{x}\right){\eta }_k\left(t\right)\text{for}k=1,\dots ,6$(14) where ${\eta }_T$ is the wave elevation at the transom, and $N_k\left(\overrightarrow{x}\right)$ and $M_k\left(\overrightarrow{x}\right)$ are canonical potentials to satisfy six rigid body motion modes. For the memory flow, the corresponding free surface conditions with an existing incident wave are as follows: (15) $\begin{array}{rl}& \frac{\mathrm{\partial }\eta }{\mathrm{\partial }\text{t}}-\left(\overrightarrow{U}\right(\overrightarrow{x})-\mathrm{\nabla }\mathrm{\Phi }(\overrightarrow{x}\left)\right)\cdot \mathrm{\nabla }\eta =\frac{{\mathrm{\partial }}^2\mathrm{\Phi }\left(\overrightarrow{x}\right)}{\mathrm{\partial }{z}^2}\left(\eta +{\eta }_i\right)\\ & +\frac{\mathrm{\partial }{\mathrm{\Phi }}_l\left(\overrightarrow{x}\right)}{\mathrm{\partial }z}+\frac{\mathrm{\partial }{\mathrm{\Phi }}_m\left(\overrightarrow{x}\right)}{\mathrm{\partial }z}\\ & -\mathrm{\nabla }\mathrm{\Phi }\left(\overrightarrow{x}\right)\cdot \mathrm{\nabla }{\eta }_i\text{on }z=0\end{array}$(15) (16) $\begin{array}{rl}& \frac{\mathrm{\partial }{\mathrm{\Phi }}_m\left(\overrightarrow{x}\right)}{\mathrm{\partial }t}-\left(\overrightarrow{U}\right(\overrightarrow{x})-\mathrm{\nabla }\mathrm{\Phi }(\overrightarrow{x}\left)\right)\cdot \mathrm{\nabla }{\mathrm{\Phi }}_m\left(\overrightarrow{x}\right)\\ & =-g\eta -\mathrm{\nabla }\mathrm{\Phi }\left(\overrightarrow{x}\right)\cdot \mathrm{\nabla }{\mathrm{\Phi }}_i\left(\overrightarrow{x}\right)\text{on }z=0\end{array}$(16) where ${\eta }_i(x,y,t)$ represents the incident wave elevation.

2.3. Solid boundary conditions

The solid boundary conditions mainly refer to the bottom boundary condition and the body boundary condition. The bottom boundary condition meets the no-flux or no-penetration conditions. Under the assumption that, by the basis flow and the memory flow on the mean body position induced by the steady forcing on the body, the instantaneous body boundary condition can be defined following Equation (17): (17) $\begin{array}{rl}& \frac{\mathrm{\partial }}{\mathrm{\partial }z}\left(\mathrm{\Phi }\left(\overrightarrow{x}\right)+{\mathrm{\Phi }}_l\left(\overrightarrow{x}\right)+{\mathrm{\Phi }}_m\left(\overrightarrow{x}\right)\right)\\ & =\left(\overrightarrow{U}\left(\overrightarrow{x}\right)+\frac{\mathrm{\partial }\eta }{\mathrm{\partial }t}\right)\cdot \overrightarrow{n}\text{on}{S}_b\end{array}$(17) For the local flow boundary value problem, the boundary condition is defined by Equation (18): (18) $\frac{\mathrm{\partial }{\mathrm{\Phi }}_l\left(\overrightarrow{x}\right)}{\mathrm{\partial }n}=\sum _{i=1}^6\frac{\mathrm{\partial }{\eta }_i}{\mathrm{\partial }t}+{\eta }_i{m}_i\text{on}{S}_b$(18) with (19) $\begin{array}{rl}& \left(n_1,n_2,n_3\right)=\overrightarrow{n}\\ & \left(n_4,n_5,n_6\right)=\overrightarrow{x}\times \overrightarrow{n}\\ & \left(m_1,m_2,m_3\right)=(\overrightarrow{n}\cdot \mathrm{\nabla })\left(\overrightarrow{U}\right(\overrightarrow{x})-\mathrm{\nabla }\mathrm{\Phi }(\overrightarrow{x}\left)\right)\\ & \left(m_4,m_5,m_6\right)=(\overrightarrow{n}\cdot \mathrm{\nabla })(\overrightarrow{x}\times (\overrightarrow{U}\left(\overrightarrow{x}\right)-\mathrm{\nabla }\mathrm{\Phi }\left(\overrightarrow{x}\right)\left)\right)\end{array}$(19) The m-terms, $m_n$, provide a coupling between the basis flow and the existing unsteady wave solution. These terms tend to be the largest at the ship and can be significant for full form ships. For memory flow with an existing incident wave, the body boundary condition becomes (20) $\frac{\mathrm{\partial }{\mathrm{\Phi }}_m\left(\overrightarrow{x}\right)}{\mathrm{\partial }n}=\left(\overrightarrow{U}\right(\overrightarrow{x})-\mathrm{\nabla }\mathrm{\Phi }(\overrightarrow{x}\left)\right)\cdot \overrightarrow{n}-\frac{\mathrm{\partial }{\mathrm{\Phi }}_i\left(\overrightarrow{x}\right)}{\mathrm{\partial }n}\text{on}{S}_b$(20)

2.4. Radiation condition

In order to satisfy the radiation condition on a truncated free surface, a numerical beach is designed to absorb waves generated by the body. A Newtonian term is applied to the kinematic free surface condition that damps all wavelengths less than about twice the extent of the numerical beach. The effect of the Newtonian cooling physically corresponding to a mass flux through the free surface is illustrated by the combination of the following free surface conditions represented by the simplest ‘Kelvin’ form together with the related wave dispersion relation: (21) $\begin{array}{rl}& \{\begin{array}{l}{\displaystyle \frac{\text{d}{\mathrm{\Phi }}_m\left(\overrightarrow{x}\right)}{\text{d}t}}=-g\eta \\ {\displaystyle \frac{\text{d}\eta }{\text{d}t}}={\mathrm{\Phi }}_{m_Z}-2\upsilon \eta +{\displaystyle \frac{v^2}{g}}{\mathrm{\Phi }}_m\end{array}\\ & \omega =i\upsilon \pm {\left(g\sqrt{k_x^2+k_y^2}\right)}^{0.5}\end{array}$(21) where $k_x$ and $k_y$ denote the wavenumbers in the $x\text{ - }$ and $y\text{-}$directions, respectively, and $\omega$ and $\upsilon$ represent the wave frequency and the uniformly distributed strength of the Newtonian cooling, respectively. Here, $\upsilon$ can be interpreted as the Rayleigh viscosity. For example, when $\upsilon >0$, the wave becomes damped with the wave frequency $\omega$ deviating from the real axis of the complex plane.

2.5. Pressure calculation

The calculations for the total linear pressure within the fluid domain can be expressed according to the Bernoulli equation: (22) $\begin{array}{rl}p& =-\rho \left(\frac{\mathrm{\partial }}{\mathrm{\partial }\text{t}}-\left(\overrightarrow{U}\right(\overrightarrow{x})-\mathrm{\nabla }\mathrm{\Phi }(\overrightarrow{x}\left)\right)\cdot \mathrm{\nabla }\right)\left({\mathrm{\Phi }}_l\left(\overrightarrow{x}\right)+{\mathrm{\Phi }}_m\left(\overrightarrow{x}\right)\right)\\ & -\rho \left(\overrightarrow{U}\left(\overrightarrow{x}\right)\cdot \mathrm{\nabla }\mathrm{\Phi }\left(\overrightarrow{x}\right)-\frac{1}{2}\mathrm{\nabla }\mathrm{\Phi }\left(\overrightarrow{x}\right)\cdot \mathrm{\nabla }\mathrm{\Phi }\left(\overrightarrow{x}\right)+gz\right)\end{array}$(22) According to the linearization, the linear pressure can also be decomposed into the memory pressure, $p_m$, the local pressure, $p_l$, and the static pressure, $p_s$, where (23) $\begin{array}{rl}{p}_m& =-\rho \left(\frac{\mathrm{\partial }}{\mathrm{\partial }t}-\left(\overrightarrow{U}\right(\overrightarrow{x})-\mathrm{\nabla }\mathrm{\Phi }(\overrightarrow{x}\left)\right)\cdot \mathrm{\nabla }\right){\mathrm{\Phi }}_m\left(\overrightarrow{x}\right)\\ & -\rho \left(\overrightarrow{U}\left(\overrightarrow{x}\right)\cdot \mathrm{\nabla }\mathrm{\Phi }\left(\overrightarrow{x}\right)-\frac{1}{2}\mathrm{\nabla }\mathrm{\Phi }\left(\overrightarrow{x}\right)\cdot \mathrm{\nabla }\mathrm{\Phi }\left(\overrightarrow{x}\right)\right)\end{array}$(23) (24) $\begin{array}{rl}{p}_l& =-\rho \left(\frac{\mathrm{\partial }}{\mathrm{\partial }\text{t}}-\left(\overrightarrow{U}\right(\overrightarrow{x})-\mathrm{\nabla }\mathrm{\Phi }(\overrightarrow{x}\left)\right)\cdot \mathrm{\nabla }\right){\mathrm{\Phi }}_l\left(\overrightarrow{x}\right)\end{array}$(24)

2.6. Higher-order discretization method

For the evaluation of the integrals in Equation (4), approximations are usually introduced for the geometry of the boundary elements as well as the physical equations. Depending on the scheme and order of the approximations, a variety of boundary element methods has been proposed. As described by Nakos (1990), as one important component, the use of the high-order representation of the solution on the free surface realizes the implementation of non-dissipative free surface conditions. As a result, more robust and accurate computations of the body generated waves over the entire computational domain can be realized.

The BIE and boundary conditions are discretized using a higher-order method based on tensor product B-splines and bi-quadratic B-splines (Kring, 1994; Kring et al., 2000, 2004). The equations are integrated on the exact geometry defined by the imported CAD geometry. Typically, bi-quadratic B-splines are adopted to represent the velocity potential distributions on the body surface and free surface: (25) $\begin{array}{rl}2\pi \mathrm{\Phi }\left(\overrightarrow{x}\right)& =\sum _{i=1}^{N_b}\sum _{j=1}^9\{{\left(\frac{\mathrm{\partial }\mathrm{\Phi }\left(p,q\right)}{\mathrm{\partial }n}\right)}_{ij}\\ & \times \int {\int }_{\text{B}\text{E}{\text{B}}_i}G\left({\overrightarrow{x}}^{\mathrm{\prime }},p,q\right)L_j\left(p,q\right)J\left(p,q\right)\text{d}p\text{d}q\\ & -(\mathrm{\Phi }\left(p,q\right){)}_{ij}\int {\int }_{\text{B}\text{E}{\text{B}}_i}\frac{\mathrm{\partial }G\left({\overrightarrow{x}}^{\mathrm{\prime }},p,q\right)}{\mathrm{\partial }n}\\ & \times L_j\left(p,q\right)J\left(p,q\right)\text{d}p\text{d}q\\ & +{\left(\frac{\mathrm{\partial }\mathrm{\Phi }\left(p,q\right)}{\mathrm{\partial }n}\right)}_{ij}\int {\int }_{\text{B}\text{E}{\text{F}}_i}G\left({\overrightarrow{x}}^{\mathrm{\prime }},p,q\right)\\ & \times L_j\left(p,q\right)J\left(p,q\right)\text{d}p\text{d}q\\ & -(\mathrm{\Phi }\left(p,q\right){)}_{ij}\int {\int }_{\text{B}\text{E}{\text{F}}_i}\frac{\mathrm{\partial }G\left({\overrightarrow{x}}^{\mathrm{\prime }},p,q\right)}{\mathrm{\partial }n}\\ & \times L_j\left(p,q\right)J\left(p,q\right)\text{d}p\text{d}q\}\end{array}$(25) The value of $G$ at any location of an element can be defined in terms of the nine nodal values and Lagrangian interpolation functions: (26) $\begin{array}{rl}{L}_j(p,q)& =\frac{1}{4}p\left(p+p_j\right)q\left(q+q_j\right),j=1,3,5,7\\ L_j(p,q)& =\frac{1}{2}\left(1-p^2{q}_j^2-p_j^2{q}^2\right)\\ & \times \left[p{p}_j\left(1+p{p}_j\right)+q{q}_j\left(1+q{q}_j\right)\right],\\ & j=2,4,6,8\\ L_j(p,q)& =\left(1-p^2\right)\left(1-q^2\right),j=9\end{array}$(26) where $p_j$ and $q_j$ represent the coordinates of the $j\text{ - th}$ node in the parametric space.

3. Experimental tests

3.1. Model description

The SWA USV is composed of a pontoon at the bottom to provide the majority of the buoyancy force and a foil strut vertical to the pontoon. According to this special design, the whole hull of the USV has a very small waterplane area to reduce wave excitation and hence motions as shown in Figure 2. Additionally, the pontoon and the foil strut both contribute to a very narrow transom to delay the flow separation as shown in Figure 1. Compared with a similar displacement monohull, the USV will benefit from the design of a very slender waterplane and delayed flow separations at the stern to achieve a much smaller residual resistance. The geometric details including both model and prototype are listed in Table 1:

Table 1. Main parameters of the SWA USV.

Parameter	Model	Prototype	Units
$\mathrm{\Lambda }$	1	25
$\mathrm{\Delta }$	13.200	206,115.917	kg
LCB	1.018	25.457	m (aft)
VCB	0.052	1.311	m
KML	2.025	50.640	m
KMT	0.150	3.750	m
LCF	1.093	27.317	m (aft)
LWL (main hull)	2.019	50.488	m
WL beam	0.060	1.494	m
Wetted surface area	0.554	346.500	m

Note: Λ - scale ratio; Δ - displacement; LCB - Longitudinal Center of Buoyancy; VCB - Longitudinal Meta Height; KML - Transverse Meta Height; LCF - Longitudinal Center of Floating Point; LWL - Wet Length(main hull); WL beam - Wet Length Beam.

Figure 2. Construction diagram of the SWA USV.

3.2. SWA USV model resistance and seakeeping tests

The towing tests for the USV, including the calm water performance tests and the regular wave seakeeping tests, are conducted in ${20}^{\circ }$C fresh water at the Davidson Laboratory high-speed towing tank, Tank 3 at Stevens Institute of Technology, which is 95 m long, 5 m wide and 2 m deep. The towing carriage rides on a monorail situated above the centerline of the tank. The monorail is 1.14 m above the water surface. A steel cable driven by an electric motor at the far end of the tank drives the carriage. Figure 3 shows the SWA USV, which is ballasted to a specified loading condition with 0.1 m draft and $0^{\circ }$ trim by placing ballast weights on weight pegs located in the model. First, the SWA USV is towed at a set speed varying in the range Fr = 0.15–0.33. During the tests, two video cameras, Cam 1 and Cam 2, which are mounted at the front and back of the model, respectively, via a carriage can clearly capture the scene of the SWA USV running against the water as shown in Figure 4. The SWA USV is free to heave and pitch but fixed in yaw, roll, surge and sway. The vertical motion of the tow-point is measured using a motion transducer attached to the free-to-heave apparatus. The pitch of the model mounting pad relative to the horizon is measured using an inclinometer mounted inside the model. The resistance is measured using a drag balance located directly above a pivot box (a towing fixture with a hinge permitting pitch motion). These heave, pitch and drag sensors can accurately measure the parameters in time and output the data through three channels to the core computer as electronic signals. The calibration of these sensors is completed before the model tests start.

Figure 3. The novel SWA USV set-up in the towing tank.

Figure 4. Scene captured by frontal and stern cameras during a seakeeping test.

The wave force $F_w$ in the model tests is actually measured in the direction of surge motions and obtained in an indirect way, following the relation (27) $F_w=F_t-F_f$(27) where the components of the wave force $F_w$ mainly consist of the diffraction and radiation force and the Froude−Krylov force; $F_t$ is the measured total force acting on the SWA USV; $F_f$ is the friction force induced by the viscosity of the water and can be estimated following the empirical formulae (28) $\begin{array}{rl}{F}_f& =C_f\frac{1}{2}\rho S{V}^2\end{array}$(28) (29) $\begin{array}{rl}{C}_f& =\frac{0.075}{{\left({\log }_{10}\text{Re}-2\right)}^2}\end{array}$(29) where $C_f$ denotes the frictional resistance coefficient in ITTC 1957 Model-Ship Correlation Line.

4. Numerical simulations

4.1. Numerical model

Considering the geometric symmetry of the SWA USV, only the half body of the USV on the port side is modelled to reduce the calculating time. The geometric model of the half body of the USV is built in Rhino^®, a Non-Uniform Rational B-Spline (NURBS) based CAD tool, as shown in Figure 5. Afterwards, the CAD representations are imported into a numerical tool called Aegir, developed based on the proposed HORS method. With automating discretization of the surfaces on the body and free surface, the Rankine source elements can be assigned in Aegir.

Figure 5. Geometric model of the novel SWA USV.

As shown in Figure 6, the free surface and body surface of the SWA USV are discretized into small quadrilateral elements, each of which are assigned nine calculation nodes. There are six mesh assignments with the coarsest mesh, No. 1, ranging to the finest mesh, No. 6. In Table 2, the corresponding minimum mesh sizes for free surfaces and body surfaces are presented, respectively. The size of the free surface mesh is actually close to that of the body surface mesh to avoid divergence induced by too large adjacent size differences. The minimum body surface mesh size is a little smaller than the minimum free surface mesh size. The convergence study is shown in Figures 7–9. With the total number of elements increasing, the results vary very little, which indicates the convergence of the simulations. Mesh No. 5 is selected for simulations in consideration of its accuracy and efficiency. As per mesh No. 5, the computational domain for the free surface is a rectangle area of $9\text{m}\times 3\text{m}$ which can be divided into $180\times 60$ elements. The body surface is divided into 972 elements.

Figure 6. Grid assignment on the body surface and free surface (the length scale in the longitudinal direction, X, is twice that in the vertical direction, Z).

Figure 7. Convergence study of the total wave force on the number of mesh elements.

Figure 8. Convergence study of the heave motion on the number of mesh elements.

Figure 9. Convergence study of the pitch motion on the number of mesh elements.

Table 2. Mesh refinements.

Mesh no.	Number of free surface meshes (Y-X)	Number of body surface meshes (Y-X)	Minimum free surface mesh size (m) dY-dX	Minimum body surface mesh size (m) dY-dX
1	60–20	162	0.15–0.15	0.09755–0.0435
2	90–30	243	0.1–0.1	0.065–0.029
3	120–40	432	0.075–0.075	0.04875–0.02175
4	150–50	675	0.06–0.06	0.039–0.0174
5	180–60	972	0.05–0.05	0.0325–0.0145
6	210–70	1323	0.043–0.043	0.0279–0.0124

4.2. Validation study of the SWA USV advancing in calm water

The calm performance of the SWA USV is studied in a validation study. The SWA USV is simulated running in calm water at different advancing velocities of 0.67, 0.93, 1.11, 1.20 and 1.47 m/s, corresponding to Fr = 0.15–0.33 via a steady state solver in Aegir. During steady state simulations, all the degrees of freedom of the SWA USV are fixed, so the effects of dynamic responses like pitch and heave are ignored. But it is noted that this assumption is only reasonable if the dynamic responses of the SWA USV are not large.

By comparing the predicted total wave forces on the SWA USV using steady state simulation to the corresponding model test, the predicted total load is found to be always smaller than that obtained from the model tests. The discrepancy gradually increases with the Froude number, becoming larger especially when Fr exceeds 0.25. This is expected, because when the SWA USV is running fast, the lift force becomes larger, which will lift the SWA USV and the excited wave elevations around the SWA USV become large. As a result, the wet area of the SWA USV will become large. But in the estimations of frictional force, a constant wet area is used, which will give rise to a smaller friction resistance force according to Equation (10) and consequently a larger measured wave force according to the Equation (9). Another explanation for why the discrepancy exists is that all the six degrees of freedom are fixed in the simulations, which is different from the model tests with pitch and heave degrees of freedom released. It is found that, when the Fr number is larger than a certain value, e.g. 0.25 for the Wigley, KCS and S60 hulls, the wave making force of a fixed ship will be smaller than that of a pitching and heaving ship.

But by examining the total wave force−Fr curve, the trends of both curves show very good agreement. With the Froude number increasing from 0.15 to 0.25, there is a peak wave force at an Fr of 0.25. Once the Froude number exceeds 0.25, the wave force suddenly drops down to a trough at an Fr of 0.27. In general, this trough point happens for Fr numbers in the range 0.3–0.35 for the classical Wigley and S60 hulls. It is found that the wave force can be reduced at a slower advancing speed for our SWA USV. When the Froude number continues to increase from 0.27 to a larger value, the total wave force keeps increasing.

Another validation study for the USV advancing in calm water is based on transient state simulation. The simulated results in the time domain are compared with corresponding model test data and steady state simulation as well as a lower-order method − Michlet theory simulated data provided in ACCeSS (Royce & Mcdonald, 2012) as shown in Figure 10 to illustrate the accuracy and improvement of the proposed time domain HORS method. By comparing the resistance and Fr curve predicted by steady state simulations and transient state simulations as shown in Figure 10, it can be observed that transient state simulations can match the model test data better. This is first because, in transient simulations, the pitch and heave motions are taken into account, which completely matches the real model test setting. Secondly, the transient simulations solve the resistance in the time domain, so more unsteady flow details and motion responses happening in the model test can be captured.

Figure 10. Total load on the advancing novel SWA USV with two side hulls at different velocities.

4.3. Validation study of the SWA USV advancing in heading wave

During the following simulations, the equations of motion in the vertical plane of the SWA USV are solved, which indicates that the pitch and heave freedoms are released. Also, the HORS method is used to solve the wave loads acting on the SWA USV in the time domain. According to the dispersion relation, wavelengths relating to the wave frequency are the critical factor that may influence the dynamics of the SWA USV; thus, the effect of wavelength on the motions’ wave loads and motion responses of the SWA USV is investigated in detail. The SWA USV is simulated advancing at a velocity of 1.335 m/s, Fr = 0.3, in a regular heading wave with a constant wave steepness, $k$, of 0.022; here, the wave steepness k = 2A$\pi$/L, but with different wavelengths, $L$, of 8.74, 6.12 and 3.70 m, corresponding to wave amplitudes of A = 0.0305, 0.0213 and 0.0107 m.

Figures 11–13 present a comparison of the simulated heave and pitch responses of the SWA USV in a heading wave with wavelengths of 8.74, 6.12 and 3.70 m with the corresponding model test data. It is found that the predicted motion responses show very good agreement with those from model tests for both long waves like $L=8.74$ m and short waves like L = 3.70 m. The capability of the proposed higher-order boundary element method is validated for solving the ship seakeeping problem in regular waves. By comparing Figures 11–13, with a decrease in wavelength the heave and pitch motions become more frequent. The amplitudes of the heave and pitch motions start to increase as well. With wavelengths decreasing from 8.74 to 3.7 m, namely the wave encounter frequency increasing from 3.47 to 6.35 rad/s and ${\omega }_e/L$ increasing from 0.397 to 1.716, the heave and pitch response amplitude of operator (RAO)s can be clearly observed to increase in simulations and experiments as shown in Figures 14 and 15. This may be because the natural frequencies of the heave and pitch motions are 6.85 and 7.04 rad/s, respectively. With wavelengths decreasing from 8.74 to 3.7 m, i.e. the wave encounters frequencies increasing from 3.47 to 6.35 rad/s, the wave frequency becomes closer to the natural frequency; as a result, larger heave and pitch responses are observed in both numerical simulations and model tests. Ultimately, it is noted that some minor instabilities can be observed even if the simulation becomes converged. Because the weight of the USV’s hull is slightly larger than the required displacement, an unloader spring is applied to give some lift force for the USV to meet the required draft. As a consequence, the spring brings some damping effect and reduces the sensitivity of the model’s heave and pitch motions.

Figure 11. Heave and pitch responses at L = 8.74 m.

Figure 12. Heave and pitch responses at L = 6.12 m.

Figure 13. Heave and pitch responses at L = 3.7 m.

Figure 14. Heave RAOs at different wave encounter frequencies.

Figure 15. Pitch RAOs at different wave encounter frequencies.

4.4. Wave loads on the SWA USV in regular waves

4.4.1. The effect of wavelength

The SWA USV is simulated advancing at a velocity of 1.335 m/s with Fr = 0.3 in a regular heading wave. By using the validated numerical method, systematic simulations are carried out to examine the effect of wavelength on the wave loads acting on the SWA USV advancing in regular heading waves with different wavelengths, $L$, of 8.74, 6.12, 3.99, 3.70 and 2.77 m, corresponding to ${\omega }_e/L$ = 0.397, 0.734, 1.502, 1.716 and 2.805 rad/s/m and wave amplitudes of A = 0.0305, 0.0213, 0.014, 0.0107 and 0.0097 m.

The total wave loads on an SWA USV moving into waves mainly consist of the components of radiation and diffraction force and the F-K force. Figures 16–18 display the time history of the radiation and diffraction force, F-K force and total wave force, respectively. In Figure 16, the wave radiation and diffraction force is oscillating irregularly initially. After reaching a certain moment, the oscillations become regular, which means the results become converged. With ${\omega }_e/L$ increasing from ${\omega }_e/L$ = 0.397 to 2.805 rad/s/m, the radiation and diffraction force decreases. When ${\omega }_e/L$ = 2.805 rad/s/m, the radiation and diffraction force starts to increase. In Figure 17, the F-K force becomes reduced continuously with the increase of ${\omega }_e/L$ from 0.397 to 2.805 rad/s/m. This is because the F-K force is related to the time derivative of the velocity potential at the surface of the ship, namely the changing rate of the potential distributions. In Figure 18, the overall curve trend and the force value of the total wave load are very similar to the F-K force shown in Figure 12. Also, by comparing the maximum F-K force with the maximum radiation and diffraction force, the former is much larger than the latter. This means that the F-K force, namely the loads acting on the ship induced by incident waves without any interaction effect, is the dominant wave load component acting on the SWA USV.

Figure 16. Time history of radiation and diffraction wave force at different wave encounter frequency to wavelength ratios.

Figure 17. Time history of F-K force at different wave encounter frequency to wavelength ratios.

Figure 18. Time history of total wave load at different wave encounter frequency to wavelength ratios.

4.4.2. The effect of wave steepness

The SWA USV is simulated advancing at a velocity of 1.335 m/s and Fr = 0.3 in a regular heading wave. In the simulations, $L=8.74,6.12$ and $3.99\text{m}$ but, for each wavelength, three wave heights are applied to explore the effect of wave steepnesses of 0.024, 0.022 and 0.014 on the wave loads acting on the SWA USV. For each fixed wavelength case, the trends of the heave and pitch curves for various k values are nearly the same. But, with increasing wave steepness, the heave and pitch motion responses also increase, as shown in Figures 19–24.

Figure 19. Time history of heave response for different wave steepnesses at L = 8.74 m.

Figure 20. Time history of pitch responses for different wave steepnesses at L = 8.74 m.

Figure 21. Time history of heave response for different wave steepnesses at L = 6.12 m.

Figure 22. Time history of pitch response for different wave steepnesses at L = 6.12 m.

Figure 23. Time history of heave response for different wave steepnesses at L = 3.99 m.

Figure 24. Time history of pitch response for different wave steepnesses at L = 3.99 m.

Figures 25–33 show the time histories of the values of the diffraction and radiation force, the F-K force and the total wave loads for three wave steepnesses, i.e. $k=0.024$, 0.022 and 0.014. It is easy to observe that, for each fixed wavelength case, with an increase in wave steepness, all the wave load components, including the diffraction and radiation force and the F-K force, increase as well. After taking the average of the converged results presented in Figures 25–33, the wave radiation and diffraction force can be obtained. The percentage ratio of the diffraction and radiation force to the total wave load is nearly a constant value at 7.28% for various wave steepnesses at $L=8.74$ m.

Figure 25. Time history of radiation and diffraction wave force for different wave steepnesses at L = 8.74 m.

Figure 26. Time history of F-K force for different wave steepnesses at L = 8.74 m.

Figure 27. Time history of total wave load for different wave steepnesses at L = 8.74 m.

Figure 28. Time history of radiation and diffraction wave force for different wave steepnesses at L = 6.12 m.

Figure 29. Time history of F-K force for different wave steepnesses at L = 6.12 m.

Figure 30. Time history of total wave load for different wave steepnesses at L = 6.12 m.

Figure 31. Time history of radiation and diffraction wave force for different wave steepnesses at L = 3.99 m.

Figure 32. Time history of F-K force for different wave steepnesses at L = 3.99 m.

Figure 33. Time history of total wave load for different wave steepnesses at L = 3.99 m.

4.4.3. The effect of advancing velocity

The SWA USV is simulated to move against regular heading wave at a series of advancing velocities of 0.67, 0.76, 0.93, 1.20, 1.47 and 1.74 m/s, which correspond to Fr = 0.15−0.39. The wave condition is kept unchanged during all the simulations with a wave height of 0.308 m, a wavelength of 3.51 m and an incoming angle of 180°. The wave force components and motion responses are obtained in the time domain at Fr = 0.15−0.39, as shown in Figures 34–38, respectively. In Figure 34, with increasing velocity and Froude number, increasing, the radiation and diffraction force amplitudes are increasing. Meanwhile, the oscillation phases also become very different. But in Figure 35, the F-K forces stay almost the same for different Froude numbers. This may be because the F-K forces are mainly induced by the incident waves without considering the interaction effect from the SWA USV. Here, the wave conditions are not changed so the F-K force is also a constant. As a result, in Figure 36, the amplitude of the oscillations of total wave load becomes different with increasing Froude number, but this difference is very small because the percentage or contribution of the radiation and diffraction forces are very small compared with the total wave load. In Figures 37 and 38, with the velocity and Froude number increasing, the heave and pitch responses increase slightly and reach a peak value of 2 and 2.1 at Fr = 0.33 for heave and pitch, respectively; afterwards, the heave and pitch responses drop to 1.7 and 0.4, respectively, at Fr = 0.39.

Figure 34. Time history of radiation and diffraction wave force for different advancing velocities.

Figure 35. Time history of F-K force for different advancing velocities.

Figure 36. Time history of total wave load for different advancing velocities.

Figure 37. Heave response of the USV for different advancing velocities.

Figure 38. Pitch response of the USV for different advancing velocities.

As shown in Figure 39(a), with increasing Froude number from 0.15 to 0.27, the radiation and diffraction forces reach peak values. With Froude numbers larger than 0.27, the radiation and diffraction forces firstly go to a trough value and then go straight up. Figure 39(b) shows the total wave force acting on the SWA USV with the Froude number increasing from 0.15 to 0.27. Comparing with Figure 39(a), a similar trend can be seen. Also compared with the total wave force−Fr curve in calm water, the total wave force−Fr curve in regular waves displays similar trends.

Figure 39. Variations of radiation and diffraction wave force and total wave load with advancing velocities.

5. Conclusions

In the present article, a higher-order Rankine source (HORS) method is proposed for the numerical study of the motion responses and wave loads of a proposed novel Small Waterplane Area (SWA) USV advancing in calm water and in regular waves, respectively. The predicted heave and pitch responses to wave loads by the SWA USV are compared with model test data. The main findings are as follows.

1. For a USV advancing in calm water and heading wave, time-domain simulations using the HORS method show very good agreement with model test data, which indicates the validity of the proposed higher-order Rankine source method (HORS). Steady state simulations can only be used for advancing in calm water. Although minor discrepancies exist, good agreement with model test data can be achieved.

2. With the Froude number of the SWA USV increasing from 0.15 to 0.25, a peak wave force is observed at an Fr of 0.25. Once the Froude number exceeds 0.25, the wave force suddenly drops down to a trough at an Fr of 0.27, which is slightly smaller than the 0.3–0.35 for the classical Wigley and S60 hulls. After the trough, the wave force rises again with increasing Froude number.

3. The F-K force is the dominant wave load and behaves more regularly than the wave radiation and diffraction force. Both the F-K force and the radiation and diffraction force are greatly affected by wave steepness and wavelength. The ship’s advancing velocity affects both the amplitude and phase of the induced wave disturbance effect, including the radiation and diffraction force.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work is supported by the Hubei Provincial Natural Science Foundation, China for Innovation Groups [No. 2021CFA026]; the Open Foundation of the State Key Laboratory of Fluid Power and Mechatronic Systems [Grant GZKF-202001]; and in part by the Fundamental Research Funds for the Central Universities [Grant 2021XJJS016].