Numerical Study on Hydrodynamics of Ships with Forward Speed Based on Nonlinear Steady Wave

Abstract

In this paper, an improved potential flow model is proposed for the hydrodynamic analysis of ships advancing in waves. A desingularized Rankine panel method, which has been improved with the added effect of nonlinear steady wave-making (NSWM) flow in frequency domain, is employed for 3D diffraction and radiation problems. Non-uniform rational B-splines (NURBS) are used to describe the body and free surfaces. The NSWM potential is computed by linear superposition of the first-order and second-order steady wave-making potentials which are determined by solving the corresponding boundary value problems (BVPs). The so-called m_j terms in the body boundary condition of the radiation problem are evaluated with nonlinear steady flow. The free surface boundary conditions in the diffraction and radiation problems are also derived by considering nonlinear steady flow. To verify the improved model and the numerical method adopted in the present study, the nonlinear wave-making problem of a submerged moving sphere is first studied, and the computed results are compared with the analytical results of linear steady flow. Subsequently, the diffraction and radiation problems of a submerged moving sphere and a modified Wigley hull are solved. The numerical results of the wave exciting forces, added masses, and damping coefficients are compared with those obtained by using Neumann–Kelvin (NK) flow and double-body (DB) flow. A comparison of the results indicates that the improved model using the NSWM flow can generally give results in better agreement with the test data and other published results than those by using NK and DB flows, especially for the hydrodynamic coefficients in relatively low frequency ranges.

1. Introduction

Over the last decades, the rapid development of computing power and the emergence of more sophisticated numerical methods have promoted the applications of numerical methods in ship hydrodynamics problems. Nevertheless, these problems still need to be simplified due to the complexity behind the physical models. It becomes even more complicated when different models need to be coupled, which is for example the case when ship maneuvering in waves is considered.

In the early stage, two-dimensional strip theory was developed as a practical way to evaluate ship hydrodynamic performances [1,2]. However, relying on the assumption that the ship is a slender body, strip theory is only suitable for low speed and high encounter wave frequency cases. In order to consider more realistic three-dimensional (3D) effects, it is not appropriate to assume that the ship is slender.

In the study of ship hydrodynamics problems, 3D potential flow theory has focused mainly on linear analysis [3]. The theory assumes that the disturbance due to the presence of a ship in waves is relatively small. When using Rankine panel methods, two linearization methods can be distinguished: the Neumann–Kelvin (NK) linearization and the double-body (DB) linearization. The former considers uniform flow as basic flow to linearize the free surface boundary conditions. The latter is essentially based on a slow-ship assumption which obtains the double-body velocity potential by treating the free surface as a rigid horizontal plane and takes the DB flow as basic flow to linearize the free surface boundary conditions. Numerous studies have been published using these two methods. For instance, Kim and Kim [3] presented a study on ship hydrodynamics comparing the NK and DB linearization methods. Similar researches discussing the advantages and disadvantages of these two methods can be found in Zhang et al. [4], Zhang and Beck [5], Zhang and Zou [6]. Attempts have also been made to include ship maneuvering in waves in the analysis, e.g., Seo and Kim [7], Zhang et al. [8]. In their studies, the mean second-order wave force was evaluated by Rankine panel method using NK or DB linearization, which was then treated as the input force in the equations for predicting maneuvering behavior. However, these two linearization methods, as described in [4], can be justified in the case of a slender ship, but they are not suitable for blunt bodies or ships moving at high speeds [9]. In light of the limitations of the NK and DB linearizations, works addressing ship hydrodynamics by using steady wave-making flow as basic flow for linearization were carried out; e.g., Gao and Zou [10] computed the linear steady wave-making flow beforehand and then applied the results to solve the diffraction and radiation problems. Recently, researchers have considered nonlinear steady flow to study the interactions between the linear periodic wave-induced flow and the nonlinear steady flow caused by the ship’s forward speed in calm water, such as the studies by Bunnik [11], Söding et al. [12] and Chillcce and el Moctar [13] in frequency domain. As for the time domain method, studies can be found in Riesner et al. [14], Riesner and el Moctar [15] and Chen et al. [9]. Though the transient effect of flow can be investigated in time domain, the boundary integral equation should be solved at each time step, which is more computationally expensive than that with frequency domain method.

In this study, a new model is proposed to compute the ship hydrodynamic forces in the frequency domain. In contrast to other methods, nonlinear steady flow is considered and the interaction between nonlinear steady flow and unsteady flow is considered not only in the body boundary condition, but also in the free surface boundary conditions in the corresponding diffraction and radiation problems. The main objective of this method is to capture the coupling factors as accurately as possible. The boundary value problems (BVPs) for the first-order and second-order steady wave-making potentials are first derived and solved, and the nonlinear steady wave-making (NSWM) potential is then approximated by linear superposition of the first-order and second-order steady wave-making potentials. Subsequently, the wave exciting forces and the radiation forces are evaluated based on the obtained NSWM flow.

A desingularized Rankine panel method with distributed sources at a small distance inside the body and above the free surface [16,17] is applied to numerically solve the problems. This method has the advantage over conventional boundary integral methods in that it separates the integration surface and the collocation surface, which in turn results in a boundary integral equation with non-singular kernels. In addition, the second-order or even higher-order derivatives of the velocity potential can be directly evaluated without complicated numerical treatments to eliminate the singularities in the integral equation; thus, the method is faster and easier to implement. In recent years, this method has been extended and applied in the analysis of 2D wave-body interaction problems, such as Feng et al. [18,19,20].

To verify the proposed model, non-uniform rational B-splines (NURBS) are used to generate the mesh both on the body surface and the free surface. The desingularized Rankine panel method is then employed to discretize and solve the boundary integral equation. The m_j terms in the body boundary condition are evaluated with nonlinear steady flow, and the free surface boundary conditions in the diffraction and radiation problems are also derived by taking the nonlinear steady flow into account. In order to identify the effects of steady flow on the unsteady flow, the present method is compared with those based on NK and DB linearizations. The computations are carried out for a submerged moving sphere and a modified Wigley hull advancing in head waves. The numerical results including the wave exciting forces, added masses and damping coefficients with the effects of different steady flows are presented.

2. Mathematical Formulations

Figure 1 shows the two coordinate systems that are used: an earth-fixed coordinate system $o_0-x_0{y}_0{z}_0$ and a coordinate system $o-xyz$ moving along with the ship at a constant speed $U$ with $ox$ positive to the bow, $oy$ positive to the port side and $oz$ directing upwards.

In $o_0-x_0{y}_0{z}_0$, based on the assumptions of ideal fluid and irrotational flow, the total velocity potential $\mathsf{\Psi }({\stackrel{\rightharpoonup }{x}}_0,t)$ should satisfy the following equations:

Laplace’s equation in fluid domain:

$${\nabla }^2\mathsf{\Psi }=0$$

(1)

The kinematic and dynamic boundary conditions on the free surface $S_F$:

$$\left(\frac{\partial }{\partial t}+\nabla \mathsf{\Psi }\cdot \nabla \right)\left[z_0-\eta \left(x_0,y_0,t\right)\right]=0$$

(2)

$$g\eta +{\mathsf{\Psi }}_t+\frac{1}{2}\nabla \mathsf{\Psi }\cdot \nabla \mathsf{\Psi }=0$$

(3)

where $\eta$ is the free surface elevation, $g$ is the gravitational acceleration. The subscripts (i.e., t, x₀, y₀) denote the derivatives with respect to the corresponding variables.

By combining Equation (2) and Equation (3), the following boundary condition on $S_F$ is derived:

$${\mathsf{\Psi }}_{tt}+\nabla \mathsf{\Psi }\cdot \nabla {\mathsf{\Psi }}_t+{\mathsf{\Psi }}_{x_0}{\mathsf{\Psi }}_{x_{0}t}+{\mathsf{\Psi }}_{y_0}{\mathsf{\Psi }}_{y_{0}t}+{\mathsf{\Psi }}_{x_0}\nabla \mathsf{\Psi }\cdot \nabla {\mathsf{\Psi }}_{x_0}+{\mathsf{\Psi }}_{y_0}\nabla \mathsf{\Psi }\cdot \nabla {\mathsf{\Psi }}_{y_0}+g{\mathsf{\Psi }}_{z_0}=0$$

(4)

The boundary condition on the body surface $S_B$:

$$\nabla \mathsf{\Psi }\cdot \stackrel{\rightharpoonup }{n}={\stackrel{\rightharpoonup }{V}}_B\cdot \stackrel{\rightharpoonup }{n}$$

(5)

where $\stackrel{\rightharpoonup }{n}$ is the unit normal vector directed inward of the body surface with $\left(n_1,n_2,n_3\right)=\stackrel{\rightharpoonup }{n},$$\left(n_4,n_5,n_6\right)=\stackrel{\rightharpoonup }{x}\times \stackrel{\rightharpoonup }{n}$; ${\stackrel{\rightharpoonup }{V}}_B$ is instantaneous velocity of the body surface $S_B$

Moreover, a radiation condition should be satisfied. The details for implementing the numerical radiation condition will be introduced in Section 3.

By using the Galilean transformation, the relation from $o_0-x_0{y}_0{z}_0$ to $o-xyz$ can be transformed as:

$$\frac{d}{dt}=\frac{\partial }{\partial t}-U\frac{\partial }{\partial x}$$

(6)

where $\frac{d}{dt}$ is the time derivative in coordinate system $o_0-x_0{y}_0{z}_0$ and $\frac{\partial }{\partial t}$ is the time derivative in the moving coordinate system $o-xyz$.

In $o-xyz$, assuming that the velocity potential $\mathsf{\Psi }(\stackrel{\rightharpoonup }{x},t)$ can be written as:

$$\mathsf{\Psi }(\stackrel{\rightharpoonup }{x},t)=-Ux+{\mathsf{\Phi }}^S(\stackrel{\rightharpoonup }{x})+\mathrm{Re}\left[A{\phi }^I(\stackrel{\rightharpoonup }{x})e^{i\omega t}+A{\phi }^D(\stackrel{\rightharpoonup }{x})e^{i\omega t}\right]+\mathrm{Re}\left\{{\displaystyle \sum _{j=1}^6\left[{\overline{\xi }}_j{\phi }_j^R(\stackrel{\rightharpoonup }{x})e^{i\omega t}\right]}\right\}$$

(7)

where $\left[-Ux+{\mathsf{\Phi }}^S(\stackrel{\rightharpoonup }{x})\right]$ is the steady velocity potential; ${\phi }^I(\stackrel{\rightharpoonup }{x})$, ${\phi }^D(\stackrel{\rightharpoonup }{x})$ and ${\phi }_j^R(\stackrel{\rightharpoonup }{x})$ $(j=1,2,\cdots ,6)$ are the spatial parts of the incident, diffraction and radiation velocity potentials, respectively; $A$ is the incoming wave amplitude, ${\overline{\xi }}_j (j=1,2,\cdots ,6)$ is the amplitude of $j$-th mode of oscillating motion, and $\omega$ is the encounter frequency.

2.1. Nonlinear Steady Wave-Making (NSWM) Problem

Substituting Equation (7) into Equations (1)–(5), using $\mathsf{\Psi }(x_0,y_0,z_0,t)=\mathsf{\Psi }(x+Ut,y,z,t)$ and Equation (6), and extracting the terms unrelated to time t, the BVP of the steady wave-making velocity potential ${\mathsf{\Phi }}^S(\stackrel{\rightharpoonup }{x})$ can be expressed in the moving coordinate system $o-xyz$ as:

Laplace’s equation in fluid domain:

$${\nabla }^2{\mathsf{\Phi }}^S=0$$

(8)

The boundary condition on the free surface $S_F$:

$$U^2{\mathsf{\Phi }}_{xx}^S-U\nabla {\mathsf{\Phi }}_{}^S\cdot \nabla {\mathsf{\Phi }}_x^S-U{\mathsf{\Phi }}_x^S\cdot {\mathsf{\Phi }}_{xx}^S-U{\mathsf{\Phi }}_y^S{\mathsf{\Phi }}_{xy}^S+{\mathsf{\Phi }}_x^S\nabla {\mathsf{\Phi }}_{}^S\cdot \nabla {\mathsf{\Phi }}_x^S+{\mathsf{\Phi }}_y^S\nabla {\mathsf{\Phi }}_{}^S\cdot \nabla {\mathsf{\Phi }}_y^S+g{\mathsf{\Phi }}_z^S=0$$

(9)

The boundary condition on the body surface $S_B$:

$$-U{n}_1+\stackrel{\rightharpoonup }{n}\cdot \nabla {\mathsf{\Phi }}^S=0$$

(10)

By using Equation (6), the steady hydrodynamic pressure can be obtained from Bernoulli’s equation:

$$p^S=-\rho \left(\frac{1}{2}\nabla {\mathsf{\Phi }}^S\cdot \nabla {\mathsf{\Phi }}^S-U{\mathsf{\Phi }}_x^S\right)$$

(11)

The steady force $F_i^S (i=1,2,\cdots ,6)$ can then be calculated by integrating the pressure over the wetted body surface:

$$F_i^S={\displaystyle \underset{S_B}{\iint }{p}^S}{n}_{i}ds, i=1,2,\cdots ,6$$

(12)

By using Equation (6), the steady free surface elevation can be obtained from Equation (3):

$${\eta }^S=\frac{U}{g}{\mathsf{\Phi }}_x^S-\frac{1}{2g}\nabla {\mathsf{\Phi }}_{}^S\cdot \nabla {\mathsf{\Phi }}_{}^S$$

(13)

The boundary condition Equation (9) is nonlinear. To solve the resulting nonlinear BVP, the velocity potential ${\mathsf{\Phi }}^S$ and the free surface elevation ${\eta }^S$ are expressed by perturbation expansion until second order as:

$$\begin{array}{l}{\mathsf{\Phi }}^S\approx {\mathsf{\Phi }}^{S(1)}+{\mathsf{\Phi }}^{S(2)}\\ {\eta }^S\approx {\eta }^{S(1)}+{\eta }^{S(2)}\end{array}$$

(14)

Substituting Equation (14) into Equations (8)–(10), the BVPs for the first- and second-order steady velocity potentials can be obtained by Taylor expansion on $z=0$ and about the mean wetted body surface ${\overline{S}}_b$. The BVP for the first-order steady velocity potential is given as:

$${\nabla }^2{\mathsf{\Phi }}^{S(1)}=0, \mathrm{in} \mathrm{fluid} \mathrm{domain}$$

(15)

$$\frac{U^2}{g}{\mathsf{\Phi }}_{xx}^{S(1)}+{\mathsf{\Phi }}_z^{S(1)}=0, \mathrm{on} z=0$$

(16)

$$\stackrel{\rightharpoonup }{n}\cdot \nabla {\mathsf{\Phi }}^{S(1)}=U{n}_1, \mathrm{on} \mathrm{the} \mathrm{mean} \mathrm{wetted} \mathrm{body} \mathrm{surface} {\overline{S}}_b$$

(17)

The BVP for the second-order steady velocity potential is given as:

$${\nabla }^2{\mathsf{\Phi }}^{S(2)}=0, \mathrm{in} \mathrm{fluid} \mathrm{domain}$$

(18)

$$\begin{array}{l}\frac{U^2}{g}{\mathsf{\Phi }}_{xx}^{S(2)}+{\mathsf{\Phi }}_z^{S(2)}=-\frac{\partial }{\partial z}\left(\frac{U^2}{g}{\mathsf{\Phi }}_{xx}^{S(1)}+{\mathsf{\Phi }}_z^{S(1)}\right){\eta }^{S(1)} +\frac{U}{g}\nabla {\mathsf{\Phi }}_x^{S(1)}\cdot \nabla {\mathsf{\Phi }}^{S(1)}\\ +\frac{U}{g}\left({\mathsf{\Phi }}_x^{S(1)}{\mathsf{\Phi }}_{xx}^{S(1)}+{\mathsf{\Phi }}_y^{S(1)}{\mathsf{\Phi }}_{xy}^{S(1)}\right),\end{array} \mathrm{on} z=0$$

(19)

$$\stackrel{\rightharpoonup }{n}\cdot \nabla {\mathsf{\Phi }}^{S(2)}=0, \mathrm{on} \mathrm{the} \mathrm{mean} \mathrm{wetted} \mathrm{body} \mathrm{surface} {\overline{S}}_b$$

(20)

2.2. Diffraction Problem

For the diffraction problem, the ship moves with a constant speed $U$ in waves without oscillations. In deep water, the incident wave velocity potential is given as:

$${\phi }^I(x,y,z)=\frac{ig}{{\omega }_0}{e}^{kz}\cdot e^{-ik(x\cos \beta +y\sin \beta )}$$

(21)

where ${\omega }_0$ is the incident wave frequency, $k={\omega }_0^2/g$ is the wave number, $\beta$ is the wave angle and $\beta =\pi$ represents head sea condition. The encounter frequency $\omega$ is defined as:

$$\omega ={\omega }_0-kU\cos \beta$$

(22)

Substituting Equation (7) into Equations (1)–(5), using $\mathsf{\Psi }(x_0,y_0,z_0,t)=\mathsf{\Psi }(x+Ut,y,z,t)$ and Equation (6), and extracting the terms related to time t and ${\phi }^D(x,y,z)$, the BVP of the diffraction potential ${\phi }^D(x,y,z)$ can be derived as:

Laplace’s equation in fluid domain:

$${\nabla }^2{\phi }^D=0$$

(23)

The free surface boundary condition on $z=0$:

$$\begin{array}{l}-{\omega }^2{\phi }^D-2i\omega U{\phi }_x^D+U^2{\phi }_{xx}^D+g{\phi }_z^D+i\omega \nabla {\mathsf{\Phi }}_{}^S\cdot \nabla {\phi }_{}^D+i\omega {\mathsf{\Phi }}_x^S{\phi }_x^D+i\omega {\mathsf{\Phi }}_y^S{\phi }_y^D-U\nabla {\mathsf{\Phi }}_{}^S\cdot \nabla {\phi }_x^D\\ -U\nabla {\phi }_{}^D\cdot \nabla {\mathsf{\Phi }}_x^S-U{\mathsf{\Phi }}_x^S{\phi }_{xx}^D+{\mathsf{\Phi }}_x^S\nabla {\mathsf{\Phi }}_{}^S\cdot \nabla {\phi }_x^D +{\mathsf{\Phi }}_x^S\nabla {\phi }_{}^D\cdot \nabla {\mathsf{\Phi }}_x^S-U{\phi }_x^D{\mathsf{\Phi }}_{xx}^S+{\phi }_x^D\nabla {\mathsf{\Phi }}_{}^S\cdot \nabla {\mathsf{\Phi }}_x^S\\ -U{\mathsf{\Phi }}_y^S{\phi }_{xy}^D-U{\phi }_y^D{\mathsf{\Phi }}_{xy}^S+{\mathsf{\Phi }}_y^S\nabla {\mathsf{\Phi }}_{}^S\cdot \nabla {\phi }_y^D+{\mathsf{\Phi }}_y^S\nabla {\phi }_{}^D\cdot \nabla {\mathsf{\Phi }}_y^S+{\phi }_y^D\nabla {\mathsf{\Phi }}_{}^S\cdot \nabla {\mathsf{\Phi }}_y^S=RHS\end{array}$$

(24)

where RHS $=$${\omega }^2{\phi }^I+2i\omega U{\phi }_x^I-U^2{\phi }_{xx}^I-g{\phi }_z^I-i\omega \nabla {\mathsf{\Phi }}_{}^S\cdot \nabla {\phi }_{}^I-i\omega {\mathsf{\Phi }}_x^S{\phi }_x^I$$-i\omega {\mathsf{\Phi }}_y^S{\phi }_y^I+U\nabla {\mathsf{\Phi }}_{}^S\cdot \nabla {\phi }_x^I$ $+U\nabla {\phi }_{}^I\cdot \nabla {\mathsf{\Phi }}_x^S+U{\mathsf{\Phi }}_x^S{\phi }_{xx}^I-{\mathsf{\Phi }}_x^S\nabla {\mathsf{\Phi }}_{}^S\cdot \nabla {\phi }_x^I$$-{\mathsf{\Phi }}_x^S\nabla {\phi }_{}^I\cdot \nabla {\mathsf{\Phi }}_x^S+U{\phi }_x^I{\mathsf{\Phi }}_{xx}^S-{\phi }_x^I\nabla {\mathsf{\Phi }}_{}^S\cdot \nabla {\mathsf{\Phi }}_x^S$$+U{\mathsf{\Phi }}_y^S{\phi }_{xy}^I$$+U{\phi }_y^I{\mathsf{\Phi }}_{xy}^S$ $-{\mathsf{\Phi }}_y^S\nabla {\mathsf{\Phi }}_{}^S\cdot \nabla {\phi }_y^I-{\mathsf{\Phi }}_y^S\nabla {\phi }_{}^I\cdot \nabla {\mathsf{\Phi }}_y^S-{\phi }_y^I\nabla {\mathsf{\Phi }}_{}^S\cdot \nabla {\mathsf{\Phi }}_y^S$

The boundary condition on the mean wetted body surface ${\overline{S}}_b$:

$$\stackrel{\rightharpoonup }{n}\cdot \nabla {\phi }^D=-\stackrel{\rightharpoonup }{n}\cdot \nabla {\phi }^I$$

(25)

It is worth noting that in Equation (24), the nonlinear steady potential ${\mathsf{\Phi }}_{}^S$ is also considered in the free surface boundary condition.

Once the diffraction potential ${\phi }^D$ is obtained, the wave exciting forces on the hull can be computed as:

$$F_j^{}=\mathrm{Re}\left(A{f}_j{e}^{i\omega t}\right), j=1,2,\cdots ,6$$

(26)

$$f_j=-\rho {\displaystyle \underset{{\overline{S}}_b}{\iint }\left[i\omega ({\phi }^I+{\phi }^D)-U({\phi }_x^I+{\phi }_x^D)+\nabla {\mathsf{\Phi }}_{}^S\cdot \nabla ({\phi }_{}^I+{\phi }_{}^D)\right]} n_{j}ds$$

(27)

2.3. Radiation Problem

For the radiation problem, it is assumed that the ship undergoes a harmonic oscillation. Similar to the BVP of diffraction potential, the radiation potential ${\phi }_j^R$ should satisfy the control equation and boundary conditions below:

Laplace’s equation in fluid domain:

$${\nabla }^2{\phi }_j^R=0,j=1,2,\cdots ,6$$

(28)

The free surface boundary condition on $z=0$:

$$\begin{array}{l}-{\omega }_{}^2{\phi }_j^R-2i\omega U{\phi }_{jx}^R+U^2{\phi }_{jxx}^R+g{\phi }_{jz}^R+i{\omega }_{}\nabla {\mathsf{\Phi }}_{}^S\cdot \nabla {\phi }_j^R+i{\omega }_{}{\mathsf{\Phi }}_x^S{\phi }_{jx}^R+i{\omega }_{}{\mathsf{\Phi }}_y^S{\phi }_{jy}^R\\ -U\nabla {\mathsf{\Phi }}_{}^S\cdot \nabla {\phi }_{jx}^R-U\nabla {\phi }_j^R\cdot \nabla {\mathsf{\Phi }}_x^S-U{\mathsf{\Phi }}_x^S{\phi }_{jxx}^R +{\mathsf{\Phi }}_x^S\nabla {\mathsf{\Phi }}_{}^S\cdot \nabla {\phi }_{jx}^R+{\mathsf{\Phi }}_x^S\nabla {\phi }_j^R\cdot \nabla {\mathsf{\Phi }}_x^S\\ -U{\phi }_{jx}^R{\mathsf{\Phi }}_{xx}^S+{\phi }_{jx}^R\nabla {\mathsf{\Phi }}_{}^S\cdot \nabla {\mathsf{\Phi }}_x^S-U{\mathsf{\Phi }}_y^S{\phi }_{jxy}^R-U{\phi }_{jy}^R{\mathsf{\Phi }}_{xy}^S+{\mathsf{\Phi }}_y^S\nabla {\mathsf{\Phi }}_{}^S\cdot \nabla {\phi }_{jy}^R\\ +{\mathsf{\Phi }}_y^S\nabla {\phi }_j^R\cdot \nabla {\mathsf{\Phi }}_y^S+{\phi }_{jy}^R\nabla {\mathsf{\Phi }}_{}^S\cdot \nabla {\mathsf{\Phi }}_y^S =0\end{array}$$

(29)

The boundary condition on the mean wetted body surface ${\overline{S}}_b$:

$$\stackrel{\rightharpoonup }{n}\cdot \nabla {\phi }_j^R=-i{\omega }_{}{n}_j+m_j$$

(30)

where the $m_j$ terms representing the coupling effect between the steady and unsteady flows are given as:

$$\begin{array}{l}\left(m_1,m_2,m_3\right)=\left(\stackrel{\rightharpoonup }{n}\cdot \nabla \right)\left(\stackrel{\rightharpoonup }{U}-\nabla {\mathsf{\Phi }}^S\right) \\ \left(m_4,m_5,m_6\right)=\left(\stackrel{\rightharpoonup }{n}\cdot \nabla \right)\left[\stackrel{\rightharpoonup }{x}\times \left(\stackrel{\rightharpoonup }{U}-\nabla {\mathsf{\Phi }}^S\right)\right]\end{array}$$

(31)

where $\stackrel{\rightharpoonup }{U}=\left(U, 0, 0\right) .$

It is worth noting that the effect of the nonlinear steady potential ${\mathsf{\Phi }}_{}^S$ occurs not only in the so-called $m_j$ terms in Equation (30), but also in the free surface boundary condition Equation (29).

Once the radiation potential ${\phi }_j^R$ is determined, the added mass $a_{kj}$ and damping coefficient $b_{kj} (k,j=1,2,\cdots ,6)$ can be obtained as:

$$\begin{array}{l}{a}_{kj}=\frac{-\rho }{{\omega }_{}^2}\mathrm{Re}{\displaystyle \underset{{\overline{S}}_b}{\iint }(i\omega {\phi }_j^R-U{\phi }_{jx}^R+\nabla {\mathsf{\Phi }}_{}^S\cdot \nabla {\phi }_j^R)n_{k}ds}\\ b_{kj}=\frac{-\rho }{{\omega }_{}}\mathrm{Im}{\displaystyle \underset{{\overline{S}}_b}{\iint }(i\omega {\phi }_j^R-U{\phi }_{jx}^R+\nabla {\mathsf{\Phi }}_{}^S\cdot \nabla {\phi }_j^R)n_{k}ds}\end{array}$$

(32)

3. Desingularized Rankine Panel Method

In this paper, a desingularized Rankine panel method is applied, where the Rankine sources are distributed inside the body and above the free surface at a distance $L_d$ according to the formula $L_d=l_d{\left(D_m\right)}^{\upsilon }$ proposed by Cao et al. [17], $l_d$ and $\upsilon$ are equal to 1.0 and 0.5 respectively and $D_m$ is the local mesh size (the square root of the local mesh area), as demonstrated in Figure 2.

A suitable radiation condition should be implemented to ensure a unique solution for the specific BVP when using the Rankine source method. Typically, the numerical techniques can be classified as follows:

• The upstream radiation condition [21]: imposing a boundary condition by difference method at the upstream boundary of the truncated free surface to ensure that no scattered waves propagate ahead of the vessel.
• The staggered method [22]: shifting the source points above the free surface a certain distance downstream.
• The fluid domain decomposition method [23]: the flow field is divided into an inner domain and an outer domain by vertical control surfaces, where the Rankine source is adopted in the inner area, while the Kelvin source is adopted in the outer domain, and the solutions are matched on the control surfaces.
• The modified Sommerfeld radiation condition [24,25]: the modified Sommerfeld radiation condition is adopted by taking account the Doppler shift of the scattered waves at the control surface that truncates the infinite fluid domain.

In this study, the radiation condition is satisfied by using the staggered method for its simple implementation and good stability. The raised source points are moved a distance $\Delta x$ toward downstream. The recommended parameter in this study is $\Delta x=\delta$, where $\delta$ denotes the average longitudinal value between two adjacent collocation points on the free surface. However, it should be noted that this numerical treatment is only valid for steady wave-making problem and the radiation problem with the Brard number $\tau =U\omega /g>0.25$.

By using NURBS, the points $(x,y,z)$ on the body and free surfaces can be described with parameter coordinate $(u,v)$ as:

$$\left[x(u,v),y(u,v),z(u,v)\right]=\left[{\displaystyle \sum _{i=0}^m{\displaystyle \sum _{j=0}^n{\omega }_{ij}{D}_{ij}{N}_{i,k}(u)N_{j,l}(v)}}\right]/\left[{\displaystyle \sum _{i=0}^m{\displaystyle \sum _{j=0}^n{\omega }_{ij}{N}_{i,k}(u)N_{j,l}(v)}}\right]$$

(33)

where $D_{ij}$ are the control points on the body and free surfaces; ${\omega }_{ij}$ is the weight; $N_{i,k}(u)$ and $N_{j,l}(v)$ are the B-spline basis functions of k(l)-th order for a given knot sequence $u=(u_0,u_1,\cdots ,u_{n+k+1})$, defined as:

$$\{\begin{array}{l}{N}_{i,0}(u)=\{\begin{array}{l}1, u_i\le u<u_{i+1}\\ 0, \mathrm{otherwise}\end{array}\\ N_{i,k}(u)=\frac{u-u_i}{u_{i+k}-u_i}{N}_{i,k-1}(u)+\frac{u_{i+k+1}-u}{u_{i+k+1}-u_{i+1}}{N}_{i+1,k-1}(u)\end{array}$$

(34)

According to Green’s theorem, the velocity potential $\phi (P)$ in the flow field can be expressed as:

$$\phi (P)={\displaystyle \underset{S}{\iint }\sigma }(Q)\frac{1}{r_{PQ}}dS={\displaystyle \underset{S}{\iint }\sigma }({\stackrel{\rightharpoonup }{x}}_s)\frac{1}{\left|{\stackrel{\rightharpoonup }{x}}_c-{\stackrel{\rightharpoonup }{x}}_s\right|}dS$$

(35)

where $P$ is the field point, $Q$ is the source point on the integration surface $S$, $r_{PQ}$ represents the distance between the field point and the source point; $\sigma (Q)$ is the source strength distribution over the surface $S.$ ${\stackrel{\rightharpoonup }{x}}_c$ and ${\stackrel{\rightharpoonup }{x}}_s$ represent the coordinates of collocation point and source point, respectively.

Applying the corresponding boundary conditions on the free surface ${\mathsf{\Gamma }}_f$ and body surface ${\mathsf{\Gamma }}_b$, the integral equations for the unknown source strengths can be established and solved. The velocity potential on ${\mathsf{\Gamma }}_f$ and the normal derivative of the velocity potential on ${\mathsf{\Gamma }}_b$ are calculated by:

$${\displaystyle \underset{S_f}{\iint }\sigma }({\stackrel{\rightharpoonup }{x}}_s^f)\frac{1}{\left|{\stackrel{\rightharpoonup }{x}}_c^f-{\stackrel{\rightharpoonup }{x}}_s^f\right|}ds+{\displaystyle \underset{S_b}{\iint }\sigma }({\stackrel{\rightharpoonup }{x}}_s^b)\frac{1}{\left|{\stackrel{\rightharpoonup }{x}}_c^f-{\stackrel{\rightharpoonup }{x}}_s^b\right|}ds={\phi }_0({\stackrel{\rightharpoonup }{x}}_c^f), {\stackrel{\rightharpoonup }{x}}_c^f\in {\mathsf{\Gamma }}_f$$

(36)

$${\displaystyle \underset{S_f}{\iint }\sigma }({\stackrel{\rightharpoonup }{x}}_s^f)\frac{\partial }{\partial n}\left(\frac{1}{\left|{\stackrel{\rightharpoonup }{x}}_c^b-{\stackrel{\rightharpoonup }{x}}_s^f\right|}\right)ds+{\displaystyle \underset{S_b}{\iint }\sigma }({\stackrel{\rightharpoonup }{x}}_s^b)\frac{\partial }{\partial n}\left(\frac{1}{\left|{\stackrel{\rightharpoonup }{x}}_c^b-{\stackrel{\rightharpoonup }{x}}_s^b\right|}\right)d{S}_b=\frac{\partial }{\partial n}{\phi }_0({\stackrel{\rightharpoonup }{x}}_c^b), {\stackrel{\rightharpoonup }{x}}_c^b\in {\mathsf{\Gamma }}_b$$

(37)

where ${\stackrel{\rightharpoonup }{x}}_c^f$ and ${\stackrel{\rightharpoonup }{x}}_c^b$ denote the collocation points on the free surface ${\mathsf{\Gamma }}_f$ and the body surface ${\mathsf{\Gamma }}_b$; ${\stackrel{\rightharpoonup }{x}}_s^f$ and ${\stackrel{\rightharpoonup }{x}}_s^b$ denote the source points on the integration surface.$S_f$ is the integration surface above the free surface ${\mathsf{\Gamma }}_f$, and $S_b$ is the integration surface inside the body surface ${\mathsf{\Gamma }}_b$. ${\phi }_0({\stackrel{\rightharpoonup }{x}}_c^f)$ is the given velocity potential at ${\stackrel{\rightharpoonup }{x}}_c^f$ and $\frac{\partial }{\partial n}{\phi }_0({\stackrel{\rightharpoonup }{x}}_c^b)$ is the given normal derivative of the velocity potential at ${\stackrel{\rightharpoonup }{x}}_c^b$.

Discretizing the body surface and the free surface into $N_b$ and $N_f$ quadrilateral panels respectively, a set of discrete equations can be obtained from the integral equations. From Equations (36) and (37) it follows:

$${\displaystyle \sum _{j=1}^{N_f}{\sigma }_j^f({\stackrel{\rightharpoonup }{x}}_s^f)\frac{1}{\left|{\stackrel{\rightharpoonup }{x}}_{ci}^f-{\stackrel{\rightharpoonup }{x}}_{sj}^f\right|}}+{\displaystyle \sum _{j=1}^{N_b}{\sigma }_j^b({\stackrel{\rightharpoonup }{x}}_s^b)\frac{1}{\left|{\stackrel{\rightharpoonup }{x}}_{ci}^f-{\stackrel{\rightharpoonup }{x}}_{sj}^b\right|}}={\phi }_0({\stackrel{\rightharpoonup }{x}}_{ci}^f) , i=1,2,\cdots ,N_f$$

(38)

$${\displaystyle \sum _{j=1}^{N_f}{\sigma }_j^f({\stackrel{\rightharpoonup }{x}}_s^f)\frac{\partial }{\partial n_i}\left(\frac{1}{\left|{\stackrel{\rightharpoonup }{x}}_{ci}^b-{\stackrel{\rightharpoonup }{x}}_{sj}^f\right|}\right)}+{\displaystyle \sum _{j=1}^{N_b}{\sigma }_j^b({\stackrel{\rightharpoonup }{x}}_s^b)\frac{\partial }{\partial n_i}\left(\frac{1}{\left|{\stackrel{\rightharpoonup }{x}}_{ci}^b-{\stackrel{\rightharpoonup }{x}}_{sj}^b\right|}\right)}=\frac{\partial }{\partial n_i}{\phi }_0({\stackrel{\rightharpoonup }{x}}_{ci}^b) , i=1,2,\cdots ,N_b$$

(39)

As can be seen from the discrete equations of Equations (38) and (39), the total number of equations is equal to the number of unknowns, i.e., $N=N_b+N_f$. Therefore, by satisfying the corresponding boundary conditions on the body surface and free surface at the collocation points, a set of linear equations for the unknown source strengths can be obtained. By solving these equations, the source strengths can be determined.

4. Numerical Results and Discussion

Two cases are studied: a sphere given by Equation (40), and a Wigley I ship [26] given by Equation (41):

$$x^2+y^2+{\left(z-h\right)}^2=r^2$$

(40)

$$y=\frac{B}{2}\left\{\left[1-{\left(\frac{z}{T}\right)}^2\right]\left[1-{\left(\frac{2x}{L}\right)}^2\right]\left[1+0.2{\left(\frac{2x}{L}\right)}^2\right]+{\left(\frac{z}{T}\right)}^2\left[1-{\left(\frac{z}{T}\right)}^8\right]{\left[1-{\left(\frac{2x}{L}\right)}^2\right]}^4\right\}$$

(41)

where $r$ is the radius of the sphere and $h$ is the submerged depth; L, B and T are the length, the beam and the draft of the hull respectively. The Wigley I ship has the length to beam ratio $L/B=10$ and the beam to draft ratio $B/T=1.6$.

Figure 3 shows the typical panel arrangements of the submerged sphere and the Wigley I ship. In addition, the panel arrangements on the raised plane (cyan) above the free surface are shown. In Figure 3a, the free surface of the computational domain extends to 5.0r upstream, 5.0r sideways and 10.0r downstream. The discretized panels of the sphere and the half width free surface are 21 × 21 and 50 × 16, respectively. In Figure 3b, the free surface of the computational domain extends to 1.0L upstream, 0.75L sideways and 1.5L downstream. The discretized panels of the half Wigley I ship and half width free surface are 30 × 10 and 76 × 19, respectively.

4.1. Results of the NSWM Problem

In order to compute the NSWM flow, a submerged moving sphere of $r=1.0 \mathrm{m}$ at three different submerged depths $\left(h =1.5r, 2.0r, 3.0r\right)$ is chosen as the study case, the Froude number is defined as $Fr=U/\sqrt{gh}$, so the results at the same Froude number will correspond to different forward speeds when the submerged depth is varied. The numerical results are compared with the analytical results of Wu and Taylor [27]. Figure 4 shows the dimensionless nonlinear wave-making resistance coefficient $C_w$ and lift force coefficient $C_L$ of the sphere, where the “linear” results are obtained by solving the BVP of the first-order steady velocity potential, the “nonlinear” results are obtained from the BVPs of the superposition of the first-order and second-order steady velocity potentials; $C_w=-F_1^S/(\rho g\pi r^3)$ and $C_L=F_3^S/(\rho g\pi r^3)$.

As can be seen from Figure 4, the present results are in good agreement with the analytical results in Wu and Taylor [27] for the linear solution. As can be seen in Figure 4a, the nonlinear results are larger than the linear results when the Froude number is less than a certain threshold value, whereas the situation reverses when it exceeds the threshold. However, a converse trend can be seen for the lift force in Figure 4b. A similar result can also be found in Kim [28]. This may be attributed to the “bow and stern wave-making effect”, i.e., when a nonlinear free surface boundary condition is considered, the pressures at the bow and the stern will be different from those when a linear free surface boundary condition is considered. In addition, the differences between the linear and nonlinear results decrease when the submerged depth increases, which demonstrates that the effect of the nonlinear boundary condition on the free surface can be ignored when the submerged depth exceeds a certain value, which is also the case in reality.

4.2. m_j-Terms

The difficulty in solving a radiation problem lies in the accurate calculation of m_j terms, which contain the second-order derivatives of the steady velocity potential in the body boundary condition [29]. In order to calculate the m_j-terms to verify the calculation accuracy of the derivatives of the velocity potential on the body surface, the desingularized method is applied to a sphere ($r = 1.0 \mathrm{m}$) moving at a speed $U=1.0 \mathrm{m}/\mathrm{s}$ in unbounded fluid.

The results of the first-order and second-order derivatives of the velocity potential are shown in Figure 5a–c, and the results of $m_1, m_2$ are shown in Figure 5d. It shows that the numerical results virtually coincide with the analytical solutions, which demonstrates that the present method is suitable for calculating the first- and second-order derivatives of the velocity potential on the body surface.

4.3. Results of the Diffraction Problem

Table 1. Surge wave exciting forces on the sphere (Fr = 0.4, h = 2.0r, kcr = 0.26808).

$\mathit{k}\mathit{r}$	${\mathit{f}}_1/(\mathit{\rho }\mathit{g}\mathit{\pi }\mathit{A}{\mathit{r}}^3{\mathit{k}}_{\mathit{e}})$ (Real Part)			${\mathit{f}}_1/(\mathit{\rho }\mathit{g}\mathit{\pi }\mathit{A}{\mathit{r}}^3{\mathit{k}}_{\mathit{e}})$ (Imaginary Part)
	Analytical Results [27]	Present (Nonlinear Steady Wave-Making, NSWM)	Present (Neumann–Kelvin, NK)	Analytical Results [27]	Present (NSWM)	Present (NK)
0.4	−0.0081	−0.00898	−0.00973	−0.5827	−0.58452	−0.61944
0.5	−0.0102	−0.01025	−0.01121	−0.4524	−0.45413	−0.47710
0.6	−0.0098	−0.00975	−0.00997	−0.3525	−0.35387	−0.35335
0.7	−0.0082	−0.00826	−0.00860	−0.2758	−0.27684	−0.27274
0.8	−0.0065	−0.00644	−0.00668	−0.2166	−0.21749	−0.21142
0.9	−0.0049	−0.00487	−0.00513	−0.1707	−0.17146	−0.16897
1.0	−0.0036	−0.00356	−0.00333	−0.135	−0.13561	−0.13653
1.2	−0.0018	−0.00182	−0.00199	−0.0851	−0.08553	−0.08533
1.4	−0.0009	−0.00089	−0.00100	−0.0541	−0.05436	−0.05635
1.6	−0.0004	−0.00042	−0.00035	−0.0346	−0.03476	−0.03859
1.8	−0.0002	−0.00020	−0.00020	−0.0222	−0.02233	−0.02638
2.0	−0.0001	−0.00012	−0.00012	−0.0143	−0.01443	−0.01533

and

Table 2. Heave wave exciting forces on the sphere (Fr = 0.4, h = 2.0r, k_cr = 0.26808).

$\mathit{k}\mathit{r}$	${\mathit{f}}_3/(\mathit{\rho }\mathit{g}\mathit{\pi }\mathit{A}{\mathit{r}}^3{\mathit{k}}_{\mathit{e}})$ (Real Part)			${\mathit{f}}_3/(\mathit{\rho }\mathit{g}\mathit{\pi }\mathit{A}{\mathit{r}}^3{\mathit{k}}_{\mathit{e}})$ (Imaginary Part)
	Analytical Results [27]	Present (NSWM)	Present (NK)	Analytical Results [27]	Present (NSWM)	Present (NK)
0.4	−0.5691	−0.56878	−0.70940	0.0310	0.02993	0.03530
0.5	−0.4380	−0.43853	−0.41167	0.0248	0.02518	0.02611
0.6	−0.3398	−0.33971	−0.34704	0.0187	0.01872	0.01980
0.7	−0.2653	−0.26541	−0.27041	0.0136	0.01344	0.01479
0.8	−0.2082	−0.20856	−0.21633	0.0097	0.00977	0.01051
0.9	−0.1642	−0.16431	−0.17955	0.0069	0.00698	0.00727
1.0	−0.1299	−0.12995	−0.15140	0.0048	0.00483	0.00506
1.2	−0.0820	−0.08205	−0.09850	0.0023	0.00224	0.00239
1.4	−0.0522	−0.05224	−0.05490	0.0011	0.00107	0.00158
1.6	−0.0334	−0.03340	−0.03710	0.0005	0.00050	0.00052
1.8	−0.0214	−0.02146	−0.02938	0.0002	0.00021	0.00021
2.0	−0.0138	−0.01383	−0.01844	0.0001	0.00010	0.00017

present the non-dimensional real and imaginary parts of the surge and heave wave exciting forces on the submerged sphere (h = 2.0r) moving at $Fr=U/\sqrt{gr}=0.4$ in head waves as function of the non-dimensional wave number obtained by the present method in comparison with the analytical results in [27], where $k_e={\omega }^2/g$. As can be seen in Table 1 and Table 2, in general the present results based on the NSWM flow are in better agreement with those in [27] than the results based on the NK flow. Some deviations are observed at low frequencies, especially for the results based on the NK flow. There are two explanations for this larger deviation at low frequencies: firstly, NK flow cannot deal accurately with the relatively high forward speed because the wave disturbance induced by the forward speed of the body is neglected. Secondly, there exists a critical frequency $k_{c}r$ at the Brard number $\tau =U\omega /g=0.25$, which is associated with the radiation condition [10]. Since the critical frequency is $k_{c}r=0.2608$ for this case, poor accuracy is resulted when the frequency is near the critical frequency.

Figure 6 shows the non−dimensional amplitudes and corresponding phase angles of heave and pitch wave exciting force/moment for the Wigley I ship advancing at $Fr=U/\sqrt{gL}=0.4$ in head waves, where $\nabla$ is the displacement volume. In order to investigate the influence of different steady flow models on wave exciting forces at various wave frequencies, the results based on the NK, DB and NSWM flows are compared in Figure 6. From Figure 6 one can observe that the present results obtained based on the three different steady flow models are in favourable agreement with the results obtained by Kara and Vassalos [30] using a 3D time domain method based on a transient free surface Green function, as well as with the experimental results by Journée [26]. In addition, one can also find that the results based on the NSWM flow and other two methods based on the NK and DB flows do not show evident differences, the reasons can be explained as follows: on one hand, though the effect of nonlinear steady flow is considered in the free surface boundary condition Equation (24), the interaction has no relation with the diffraction potential in the body surface boundary condition Equation (25); on the other hand, the small differences can be attributed to the predominant proportion of the Froude–Krylov force in the wave exciting force. Therefore, it can be concluded that the effect of the NSWM potential contributes unremarkably to the wave exciting force.

Figure 7 shows the real part of the diffraction wave contour for the submerged sphere ($Fr=0.4$, $kr=0.5$) and the Wigley I ship ($Fr=0.4$, $kL=2\pi$) based on the NSWM flow.

4.4. Results of the Radiation Problem

Table 3. Added masses and damping coefficients of the sphere in surge motion (Fr = 0.4, h = 2.0r, k_ecr = 0.3906).

${\mathit{k}}_{\mathit{e}}\mathit{r}$	A11			B11
	Analytical Results [27]	Numerical Results [10]	Present (NSWM)	Analytical Results [27]	Numerical Results [10]	Present (NSWM)
0.6	1.2378	1.2358	1.25532	0.0362	0.0365	0.03447
0.7	1.1615	1.1614	1.16928	0.0247	0.0250	0.02586
0.8	1.1021	1.1021	1.10996	0.0195	0.0197	0.01972
0.9	1.0545	1.0544	1.06164	0.0169	0.0170	0.01724
1.0	1.0154	1.0153	1.02189	0.0154	0.0155	0.01572
1.5	0.8934	0.8933	0.89818	0.0113	0.0113	0.01144
2.0	0.8310	0.8310	0.83495	0.0079	0.0079	0.00794
2.5	0.7941	0.7940	0.79751	0.0052	0.0052	0.00519
3.0	0.7699	0.7699	0.77304	0.0033	0.0033	0.00327
3.5	0.7529	0.7529	0.75584	0.0021	0.0021	0.00207
4.0	0.7403	0.7403	0.74306	0.0013	0.0013	0.00131
4.5	0.7306	0.7305	0.73318	0.0008	0.0008	0.00081
5.0	0.7228	0.7228	0.72529	0.0005	0.0005	0.00050

,

Table 4. Added masses and damping coefficients of the sphere in sway motion (Fr = 0.4, h = 2.0r, k_ecr = 0.3906).

$\mathit{k}\mathit{r}$	A22			B22
	Analytical Results [27]	Numerical Results [10]	Present (NSWM)	Analytical Results [27]	Numerical Results [10]	Present (NSWM)
0.6	1.1330	1.1357	1.13521	0.0771	0.0769	0.09103
0.7	1.0517	1.0528	1.05084	0.0650	0.0655	0.06482
0.8	0.9933	0.9946	0.99414	0.0544	0.0551	0.05481
0.9	0.9406	0.9506	0.95030	0.0454	0.0459	0.04605
1.0	0.9159	0.9167	0.91640	0.0380	0.0385	0.03851
1.5	0.8215	0.8219	0.82204	0.0164	0.0166	0.01656
2.0	0.7783	0.7787	0.77908	0.0076	0.0076	0.00765
2.5	0.7535	0.7538	0.75435	0.0037	0.0037	0.00368
3.0	0.7371	0.7374	0.73805	0.0019	0.0019	0.00184
3.5	0.7254	0.7256	0.72638	0.0010	0.0010	0.00095
4.0	0.7165	0.7168	0.71757	0.0005	0.0005	0.00051
4.5	0.7095	0.7098	0.71066	0.0003	0.0003	0.00027
5.0	0.7040	0.7042	0.70507	0.0002	0.0002	0.00015

and

Table 5. Added masses and damping coefficients of the sphere in heave motion (Fr = 0.4, h = 2.0r, k_ecr = 0.3906).

$\mathit{k}\mathit{r}$	A33			B33
	Analytical Results [27]	Numerical Results [10]	Present (NSWM)	Analytical Results [27]	Numerical Results [10]	Present (NSWM)
0.6	1.0569	1.0526	1.05860	0.1090	0.1119	0.09671
0.7	0.9928	0.9925	0.98991	0.0862	0.0874	0.08827
0.8	0.9449	0.9447	0.94204	0.0709	0.0712	0.07146
0.9	0.9076	0.9076	0.90508	0.0597	0.0606	0.06026
1.0	0.8780	0.8775	0.87549	0.0511	0.0606	0.05162
1.5	0.7912	0.7911	0.78990	0.0264	0.0261	0.02659
2.0	0.7504	0.7503	0.74954	0.0145	0.0146	0.01473
2.5	0.7272	0.7271	0.72632	0.0083	0.0083	0.00828
3.0	0.7123	0.7123	0.71135	0.0048	0.0048	0.00478
3.5	0.7020	0.7021	0.70125	0.0028	0.0029	0.00276
4.0	0.6942	0.6942	0.69367	0.0017	0.0017	0.00169
4.5	0.6882	0.6882	0.68768	0.0010	0.0010	0.00100
5.0	0.6833	0.6833	0.68286	0.0006	0.0006	0.00064

present the added masses and damping coefficients of the submerged sphere (h = 2.0r) moving at $Fr=U/\sqrt{gr}=0.4$ in surge, sway and heave motions respectively, where the added masses and damping coefficients are non-dimensionalized as $A_{ij}=a_{ij}/(\pi \rho r^3)$, $B_{ij}=b_{ij}/(\pi \rho \omega r^3), i,j=1,2,3$. The k_ecr corresponds to the critical frequency ${\omega }_c$ at the Brard number $\tau =0.25$. As it can be seen from these tables, the numerical results based on the NSWM flow agree well with the analytical results in [27] and the numerical results in [10]. It should be noted that the linear steady wave-making potential ${\mathsf{\Phi }}^S$ was used to evaluate $m_j$ terms in [10], without considering ${\mathsf{\Phi }}^S$ in the free surface boundary condition. From these results, it can be concluded that the effects of the free surface nonlinearities are very weak due to the submerged depth.

Table 6. Coupling added masses and damping coefficients of the sphere (Fr = 0.4, h = 2.0r, k_ecr = 0.3906).

$\mathit{k}\mathit{r}$	A13		A31		B13		B31
	Analy. [27]	Present (NSWM)	Analy. [27]	Present (NSWM)	Analy. [27]	Present (NSWM)	Analy. [27]	Present (NSWM)
0.6	−0.0269	−0.02001	0.0269	0.01927	−0.0985	−0.09120	0.0985	0.08952
0.7	−0.0075	−0.00918	0.0075	0.00840	−0.0791	−0.07988	0.0791	0.07865
0.8	0.0030	0.00331	−0.0030	−0.00356	−0.0642	−0.06466	0.0642	0.06351
0.9	0.0089	0.00930	−0.0089	−0.00935	−0.0527	−0.05271	0.0527	0.05193
1.0	0.0122	0.01235	−0.0122	−0.01253	−0.0437	−0.04345	0.0437	0.04291
1.5	0.0135	0.01329	−0.0135	−0.01360	−0.0187	−0.01880	0.0187	0.01818
2.0	0.0095	0.00958	−0.0095	−0.00960	−0.0093	−0.00954	0.0093	0.00885
2.5	0.0061	0.00624	−0.0061	−0.00602	−0.0056	−0.00567	0.0056	0.00533
3.0	0.0038	0.00377	−0.0038	−0.00377	−0.0042	−0.00412	0.0042	0.00385
3.5	0.0023	0.00226	−0.0023	−0.00229	−0.0036	−0.00365	0.0036	0.00341
4.0	0.0014	0.00141	−0.0014	−0.00145	−0.0034	−0.00349	0.0034	0.00321
4.5	0.0009	0.00089	−0.0009	−0.00085	−0.0034	−0.00341	0.0034	0.00315
5.0	0.0006	0.00055	−0.0006	−0.00054	−0.0033	−0.00334	0.0033	0.00315

presents the coupling added masses and coefficients. It can be seen that the present results almost show the reverse relations, i.e., $(A_{13}, B_{13})=(-A_{31}, -B_{31})$, which is consistent with the results in [27].

Figure 8 shows the heave and pitch hydrodynamic coefficients of the Wigley I ship at $Fr=U/\sqrt{gL}=0.4$, where the coupling hydrodynamic coefficients are non-dimensionalized as $A_{35(53)}=a_{35(53)}/(\rho \nabla L)$, $B_{35(53)}=b_{35(53)}/(\rho \nabla L\sqrt{g/L})$. As it can be seen in Figure 8, good agreement is achieved among the present numerical results and the results in [30] using NK flow and a transient free surface Green function method, and the experimental results by Journée [26]. The results based on the NSWM flow show in general better agreement with the experimental results than those obtained using DB and NK flows, especially in the low frequency ranges. However, a remarkable deviation can be observed for the heave damping coefficient in Figure 8b. This is because the uniform flow is taken as the basic flow in the method using NK flow, correspondingly the second-order derivatives of the steady potential ${\mathsf{\Phi }}^S$ are neglected in the $m_j$ terms. As a result, this treatment cannot accurately reflect the interaction between the steady flow and the unsteady flow in the body surface boundary condition. On the other hand, as explained in [9], the larger contribution of the steady velocity potential in both the body surface and free surface boundary conditions at low frequencies also leads to relatively larger deviations, but these effects are not fully considered in the method using DB flow. Therefore, the coupling effects between the nonlinear steady flow and the unsteady flow, which are reflected in both the free surface boundary condition and the $m_j$ terms in the body surface boundary condition, are quite important for the prediction of hydrodynamic coefficients, especially at low frequencies.

Figure 9 shows the real part of the heave radiation wave contour of the sphere ($Fr=U/\sqrt{gr}=0.4$, $k_{e}r=2.0$) and the Wigley I ship ($Fr=0.4, \omega /\sqrt{g/L}=3.0$) based on the NSWM flow.

5. Conclusions

In this paper, a desingularized Rankine panel method based on the NSWM flow is applied for analysis of the hydrodynamic problems of a ship advancing in waves. NURBS are used to describe the body surface and the free surface. The wave exciting forces and the hydrodynamic coefficients are computed by solving the diffraction problem and radiation problem, respectively. A numerical study is carried out for a submerged sphere and a modified Wigley hull advancing in head waves. The numerical results are compared with the analytical solutions as well as other numerical results and experimental results available in literature. The following conclusions can be drawn.

(1) The numerical results of the wave exciting forces, added masses and damping coefficients computed using the present numerical method show good agreement with the published numerical and experimental results, which verifies the reliability of the present method. A comparison among the results indicates that the method based on the NSWM flow can generally give better agreement with the experimental and other published results than those based on NK and DB flows, especially for the hydrodynamic coefficients in relatively low frequency ranges.

(2) The NSWM potential has an influence on the prediction of the wave exciting forces. However, differences among different steady flow models are not very remarkable due to the dominant proportion of the Froude–Krylov force for the considered cases. The coupling effects between the nonlinear steady flow and the linear unsteady flow are important for the prediction of hydrodynamic coefficients, particularly at low frequencies.

(3) Compared with the time domain method, considering the NSWM flow as basic flow can be used as a more practical and faster numerical tool for evaluating the hydrodynamic performances of a ship in the early design stage.

In the present study, the method based on the NSWM flow is only applied for a submerged sphere and a modified Wigley hull. For reliable verification and application of this numerical method, further study on various ship forms needs to be carried out. Besides, in the numerical study, the squat of the hull (i.e., the trim and sinkage) is neglected. For the cases at larger forward speed, the numerical accuracy could be further improved by taking the effects of trim and sinkage into account. In addition, the desingularized Rankine panel method is only applied for the cases of Brard number $\tau$ larger than 0.25, where the radiation condition is satisfied by the staggered method. For $\tau <0.25$, more robust methods for satisfying the radiation condition, such as the modified Sommerfeld radiation condition in [24,25], should be adopted. These will be the focuses of the future studies.