Journal of

Marine Science
and Engineering

Article
Hydrodynamic Analysis and Motions of Ship with Forward
Speed via a Three-Dimensional Time-Domain Panel Method
Peng Zhang 1 , Teng Zhang 2, * and Xin Wang 3

1 School of Ocean Science and Technology, Dalian University of Technology, Dalian 116024, China;
zp31803017@[Link]
2 School of Naval Architecture and Ocean Engineering, Dalian University of Technology, Dalian 116024, China
3 College of Navigation, Dalian Maritime University, Dalian 116026, China; [Link]@[Link]
* Correspondence: dllgbsh10@[Link]

Abstract: A new three-dimensional (3D) time-domain panel method is developed to solve the ship
hydrodynamic problem and motions. For an advancing ship with a constant forward speed in
regular waves, the ship’s hull can be discretized and processed into a number of quadrilateral panels.
Based on Green’s theorem, an analytical expression for Froude–Krylov (F–K) forces evaluation
on the quadrilateral panels is derived without accuracy loss. Within the linear potential theory,
the transient free surface Green function (TFSGF) is applied to solve the boundary value problem.
To improve the efficiency and numerical stability of TFSGF evaluation, a precise integration method
with variable parameters setting for extended identity matrix is developed to compute the TFSGF
in the computation domain. Then, radiation and diffraction forces can be evaluated by means of
the impulse response function method. The Wigley I hull form is taken as a study case, and the
computed hydrodynamic coefficients, wave exciting forces, and motions by the present method are
compared with previous literature experimental data and prior published results. It manifests that



the three-dimensional time-domain panel method proposed in this paper has good accuracy.



Citation: Zhang, P.; Zhang, T.; Wang, Keywords: time-domain panel method; hydrodynamic problem; Froude–Krylov forces; transient
X. Hydrodynamic Analysis and free surface Green function; precise integration method
Motions of Ship with Forward Speed
via a Three-Dimensional
Time-Domain Panel Method. J. Mar.
Sci. Eng. 2021, 9, 87. https:// 1. Introduction
[Link]/10.3390/jmse9010087
For the initial stages of ship design, accurate and reliable predictions of ship hydrody-
namic analysis and motions in waves are essential. Various numerical methods are required
Received: 27 November 2020
Accepted: 12 January 2021
to be developed for ship seakeeping analysis. For the wave-ship interaction problem with
Published: 15 January 2021
large size, the potential flow theory is much more efficient than RANS (Reynolds-averaged
Navier-Stokes) simulation [1], which is widely applied to a practical engineering problem.
Publisher’s Note: MDPI stays neutral
In the early researches, ship hydrodynamic analysis is developed based on two-
with regard to jurisdictional claims in
dimensional strip theories. Ogilvie and Tuck [2], Tasai [3], and Salvesen et al. [4] proposed
published maps and institutional affil- a new strip theory, rational strip theory, and STF, method respectively. The STF strip
iations. method is most widely used in ship motion calculation and structure design. Fonseca
and Guedes Soares [5,6] formulated the ship hydrodynamic analysis in the time domain.
Tavakoli et al. [7] investigated unsteady planning motion in waves using towing tank
tests, Computational Fluid Dynamics (CFD), and the 2D+t model. These two-dimensional
Copyright: © 2021 by the authors.
methods have been applied in the ship hydrodynamic and motion analysis for a long time.
Licensee MDPI, Basel, Switzerland.
However, for the strip theories, the flow is assumed to be constrained in two-dimensional
This article is an open access article
sections. Accurate hydrodynamic analysis can only be carried out for slender ships,
distributed under the terms and and high frequency and low-speed assumptions are also required.
conditions of the Creative Commons The shortcomings of strip theory can be overcome by 3D panel methods. Nakos [8]
Attribution (CC BY) license (https:// used the Rankine panel method to study the ship seakeeping problem in the frequency
[Link]/licenses/by/ domain. Kring [9] and Chen [10] applied the Rankine panel method to ship hydrodynamic
4.0/). analysis in the time domain. The Rankine panel methods employ the Green function

J. Mar. Sci. Eng. 2021, 9, 87. [Link] [Link]

J. Mar. Sci. Eng. 2021, 9, 87 2 of 19

1/r as a fundamental solution of the Laplace equation; the evaluation of 1/r integration
over the boundaries is easy to be carried out. However, the Rankine Green function 1/r
does not satisfy any boundary conditions, and many more panels are required for mesh
discretization of the free surface, which would greatly reduce computational efficiency.
The problems caused by the Rankine panel method can be avoided by using the 3D
free surface Green function method, which employs the 3D free surface Green function
as the fundamental solution of the Laplace equation. The 3D free surface Green function
method only requires the discretization of the ship wetted hull surface, but the evaluation
of the 3D free surface Green function is quite complex. Wehausen and Laitone [11] deduced
expressions for 3D free surface Green function in arbitrary water depth, which laid the
foundation of solving the hydrodynamic problem by 3D free surface Green function
method. Blandeau and Francois [12] solved the radiation and diffraction problem for
FPSOs by using the HydroSTAR software, in which the original codes were developed by
the frequency domain Green function method. Furthermore, Wu and Eatock Taylor [13]
solved the hydrodynamic problem for practical ships by using the frequency domain
Green function with speed. When the frequency domain Green function method is used
to solve the ship hydrodynamic problems with speed, the numerical calculation of the
frequency domain Green function is much more complex and time-consuming. Although
the numerical results are close to the experimental values, there are still some problems
in the treatment of waterline integral terms, which hinders its wide application. It’s easy
to formulate and solve ship hydrodynamic problems and motion problems in the three-
dimensional time domain. Liapis [14] applied the transient free surface Green function
(TFSGF) to solve the linear radiation problem for a ship with constant forward speed.
King [15] further extended to the linear diffraction problem for ships, and the Froude–
Krylov (F–K) forces were evaluated by Gaussian quadrature. However, the F–K forces near
the mean free surface may not be accurately evaluated due to wave volatility. Rodrigues
and Guedes Soares [16] evaluated the Froude–Krylov forces by analytical exact pressure
integration expressions, allowing for considerably coarse meshes with no loss of accuracy.
However, radiation and diffraction forces are kept linear by the indirect time-domain
method, which is difficult to calculate the hydrodynamic coefficients in the high-frequency
range. Zhang et al. [17] studied the influences of the water line integral terms on the
wave diffraction force of a moving floating body and pointed out that the influences of
the water line integral terms on the first-order force could be ignored. Sun [18] developed
a 3D time-domain program based on the TFSGF method, which can be used to solve the
ship hydrodynamic problems in waves. Singh and Sen [19] studied seakeeping problems
under different nonlinear levels, and computations were carried out for a Wigley hull
and an S175 hull in waves. Datta et al. [20] carried out modifications for three fishing
vessels and presented a variety of calculated motion results for different wave angles.
Lin [21,22] proposed a three-dimensional time-domain approach to study the ship’s large-
amplitude motions in a seaway and developed the software LAMP (Large Amplitude
Motion Program) for ship hydrodynamic analysis and wave loads calculation.
The accurate and efficient evaluation of TFSGF is essential for the ship hydrodynamic
analysis in the time domain. According to the oscillating properties of the wave part
of TFSGF, King [15] and Shan [23] divided the time computation domain into different
regions to evaluate the TFSGF and its derivatives, and series expansions and asymptotic
expansions were applied to the small time computation domain and large time computation
domain, respectively. Clement [24] firstly found that the TFSGF and its derivatives are the
solutions of the ordinary differential equations (ODEs). Shen et al. [25] solved the ODEs
by using the fourth-order Runge–Kutta method (RK44), in which the numerical instability
would occur after a long time simulation even with a very small time step size. Based
on the precise integration method (PIM) proposed by Zhong [26], Li et al. [27] solved the
ODEs, which could greatly improve the numerical stability even with a quite large time
step size. However, it can be very time-consuming due to the quite high order of the
coefficient matrix.
 ODEs by using the fourth-order Runge–Kutta method (RK44), in which the numerical in-
stability would occur after a long time simulation even with a very small time step size.
Based on the precise integration method (PIM) proposed by Zhong [26], Li et al. [27]
solved the ODEs, which could greatly improve the numerical stability even with a quite
J. Mar. Sci. Eng. 2021, 9, 87 large time step size. However, it can be very time-consuming due to the quite high 3order of 19
of the coefficient matrix.
Within the linear potential theory, the main objective of this paper is to develop a
three-dimensional time-domain panel method to study the ship’s hydrodynamic analysis
Within the linear potential theory, the main objective of this paper is to develop a
and motions. Analytical
three-dimensional integration
time-domain expressions
panel method to forstudy
F-K force formulations
the ship’s over theanaly-
hydrodynamic quad-
rilateral panels can be derived by using Green’s theorem, which
sis and motions. Analytical integration expressions for F–K force formulations over the can avoid computational
errors by numerical
quadrilateral panels integration
can be derived methods.
by usingTo improve the accuracy
Green’s theorem, andcan
which numerical stability
avoid computa-
of TFSGF evaluation, a numerical method is developed to solve
tional errors by numerical integration methods. To improve the accuracy and numerical the TFSGF by using a
precise integration
stability method with
of TFSGF evaluation, varying parameter
a numerical method is settings.
developed Based on the
to solve theimpulse
TFSGF by re-
sponse function method, the ship radiation problem and diffraction
using a precise integration method with varying parameter settings. Based on the impulse problem are solved
by using function
response TFSGF. The Wigley
method, theI ship
hull radiation
is taken as a studyand
problem case; convergence
diffraction studies
problem arefor ship
solved
hydrodynamic analysis and motions are conducted with respect
by using TFSGF. The Wigley I hull is taken as a study case; convergence studies for ship to time step size and hull
discretization. The computed hydrodynamic coefficients, wave
hydrodynamic analysis and motions are conducted with respect to time step size and hull exciting forces, and mo-
tions for the ship
discretization. Theare comparedhydrodynamic
computed to other solutions, such as Magee’s
coefficients, method
wave exciting [28], previous
forces, and mo-
literature
tions for the experimental data [29],
ship are compared and so
to other on. Thus,
solutions, such theasthree-dimensional
Magee’s method [28], time-domain
previous
panel method
literature proposeddata
experimental in this paper
[29], andissovalidated.
on. Thus, the three-dimensional time-domain
panel method proposed in this paper is validated.
2. Materials
2. Materials
2.1. Coordinate Systems
2.1. Coordinate Systems
For the present linear ship hydrodynamic analysis and motion problem, as shown in
For1,the
Figure present
a freely linearship
floating shipishydrodynamic analysiswith
considered to advance and constant
motion problem, as shown
forward speed U in
in Figure 1, a freely floating ship is considered to advance with constant
the presence of a linear incident wave field. The o − xyz is a reference, right-handed, Car- forward speed
U in the presence of a linear incident wave field. The o-xyz is a reference, right-handed,
tesian coordinate system with its origin o located amidship and travels along with the ship
Cartesian coordinate system with its origin o located amidship and travels along with the
at the same speed U, the o-xy plane is coincident with the mean free surface z = 0, the
ship at the same speed U, the o-xy plane is coincident with the mean free surface z = 0,
positive x-axis is pointing upstream, the positive y-axis is pointing portside, and the pos-
the positive x-axis is pointing upstream, the positive y-axis is pointing portside, and the
itive z-axis is pointing vertically upwards. G − xb yb zb is a body-fixed, right-handed, Car-
positive z-axis is pointing vertically upwards. G-x b yb zb is a body-fixed, right-handed,
tesian coordinate
Cartesian coordinate system,
system, andandthe theorigin
originGGisislocated
locatedatatthe thegravity
gravity ofof the ship. At
the ship. At initial
initial
time t = 0, the space fixed coordinate system O − XYZ coincides with
time t = 0, the space fixed coordinate system O-XYZ coincides with the reference coordinate the reference coor-
dinate system o − xyz , and
system o-xyz, and Gzb axis is aligned Gz b axis is aligned with ox
with ox axis. The fluid domain Ω is enclosed byΩtheis
axis. The fluid domain
ship hull surface
enclosed SB , free
by the ship hullsurface SF ,Sand
surface surface
B , free surface SFinfinity.
S∞ at n is theSunit
, and surface ∞ at normal n is
[Link]
pointing inward the ship hull surface.
the unit normal vector pointing inward the ship hull surface.

z y
Wave
zb yb U
∇ o x
xb
n G
S∞
SF SB Ω

Figure 1.
Figure 1. The
The coordinate
coordinate systems
systems and
andfluid
fluiddomain.
domain.

2.2.
2.2. Boundary
Boundary Value
Value Problem
Problem
Based on the linear potentialflow
Based on the linear potential flowtheory,
theory,thethefluid
fluidis is assumed
assumed to to
be be irrotational,
irrotational, in-
inviscid, andincompressible,
viscid, and incompressible,andandthere
thereisisnono fluid
fluid separation
separation or or lifting
lifting effect
effect [30].
[30]. In
In the
the
reference
reference coordinate system o-xyz,
coordinate system o − xyzthe total
, the velocity
total potentialΦ(Φr;(tr);tin
velocitypotential ) the infinite
in the water
infinite wa-
can be written as
ter can be written as
7
Φ(r; t) = −Ux + ΦS (r) + Φ0 (r; t) + ∑ Φk (r; t) (1)
k =1

where r = ( x, y, z) is position vector on the ship hull surface, t is the time instant,
−Ux + ΦS (r) is the steady wave velocity potential due to constant forward speed U, Φ0
is the incident wave velocity potential, Φk (k = 1, 2, · · · , 6) is radiation velocity potential,
and Φ7 is the diffraction velocity potential.
J. Mar. Sci. Eng. 2021, 9, 87 4 of 19

The perturbation potential Φk (k = 1, 2, · · · , 7) satisfies the following governing

equations and initial-boundary value conditions:

∇ Φk (r; t) = 0, wherer ∈ Ω
 2



∂ Φk
2

+ ∂Φ
 k
∂z = 0, onz = 0




 ∂t2

 .
∂Φk
∂n = mk ζ k + nk ζ k , k = 1, 2, · · · , 6, onSB






∂Φ7 ∂Φ0
∂n = − ∂n , onSB (2)

 ∇Φk → 0, as x2 + y2 → ∞ or z → −∞

 p




Φk and ∂Φ

∂t → 0, at t = 0, ( k = 1, 2, · · · , 6)

 k




 Φ7 and ∂Φ ∂t → 0, at t = − ∞.
7



where ζ k is unsteady motion in the kth mode; n = (n1 , n2 , n3 ); r × n = (n4 , n5 , n6 );

(m1 , m2 , m3 ) = (0, 0, 0); (m4 , m5 , m6 ) = (0, Un3 , −Un2 ); g is the gravity acceleration [31].
The TFSGF is adopted to solve the boundary value problem and written as [14]

G ( P, Q; t − τ ) = G ( P; Q)δ(t − τ ) + G
e( P, Q; t − τ ) H (t − τ ) (3)

where δ(·) is the Dirac function; H (·) is Heaviside unit step function; P(x, y, z) is field point;
Q(ξ, η, ζ) is source point; τ is the retard time. The Rankine part G ( P; Q) and memory part
e( P, Q; t − τ ) of G ( P, Q; t − τ ) can be given as
G

G ( P, Q) = 1/r − 1/r 0
(
(4)
e( P, Q; t − τ ) = 2 ∞ √ gυeυ(z+ζ ) sin √ gυ(t − τ ) J0 (υR)dυ
R  
G 0

where r = | P − Q|; r 0 = | P − Q0 |; Q0 (ξ, η, −ζ ) is the image point of Q about the mean free
surface; R is the horizontal distance between the field point P and source point Q; J0 is a
Bessel function of order zero; υ is wave number.
The integrations for G and its derivative over a quadrilateral panel can be computed
by Hess–Smith method [32]. G e and its derivatives can be solved numerically in Section 3.2.
The boundary integral equation for the perturbation potential Φk (k = 1, 2, · · · , 7)
can be written as
s s ∂Φk ( Q;t)
2πΦk ( P; t) + SB Φk ( Q; t) ∂n∂G
dSQ = SB G ∂n dS
Q Q
Rt s ∂Φk ( Q;τ ) Rt s e( P,Q;t−τ )
∂G
+ t0 dτ SB G ( P, Q; t − τ ) ∂n Q dS − t0
e dτ SB Φk ( Q, τ ) ∂n Q dS
2 e( P,Q;t−τ ) (5)
U2 t
Rt
+ Ug e( P, Q; t − τ ) ∂Φk (Q;τ ) dη − Φk ( Q; τ )
∂G
H R H
t0 dτ Γ0 G ∂ξ g t0 dτ Γ0 ∂ξ dη

2U t e( P,Q;t−τ )
∂G
Γ0 Φk ( Q; τ )
R H
+ g t0 dτ ∂τ dη

where Γ0 is the mean waterline; t0 is the initial time; nQ is the unit normal vector of source
point Q pointing inward the hull surface.
To solve the perturbation potential Φk , the trapezoidal rule is adopted for convolution
integration, and the mean wet hull surface SB in Equation (5) can be discretized by using
the constant panel method [31].

2.3. The Hydrodynamic Problem Formulation

Based on the impulse response function method, the radiation potential Φk (k = 1, 2, · · · , 6)
can be written as Z t .
Φk ( P; t) = Φ̂k ( P; t − τ )ζ k (τ )dτ (6)
0
J. Mar. Sci. Eng. 2021, 9, 87 5 of 19

where Φ̂k ( P; t) is the impulse response function of radiation potential Φk in the kth mode.
The Φ̂k ( P; t) can be decomposed into the following form [14]

Φ̂k ( P; t) = ψ1k ( P)δ(t) + ψ2k ( P) H (t) + χk ( P; t) (7)

where ψ1k and ψ2k are the impulsive potentials, and χk is the transient potential.
The expression of radiation forces Fjk can be written as

.. . Z t .
Fjk (t) = − a jk ζ k (t) − b jk ζ k (t) − c jk ζ k (t) − dτK jk (t − τ )ζ k (τ ) (8)
0

where Fjk is the radiation force in jth mode due to kth mode motion; K jk (t) =
s  k  s s
ρ SB ∂χ n
∂t j − χ k j dS; a jk
m = −ρ SB ψ1k n j dS; c jk = −ρ SB ψ2k m j dS; b jk =
s

−ρ SB ψ2k n j − ψ1k m j dS.
The added mass A jk (ω ) and damping coefficients Bjk (ω ) can be obtained via Fourier
transform (j = 1, 2, · · · , 6; k = 1, 2, · · · , 6) [15]:

 A jk (ω ) = a jk − c jk2 − 1 ∞ K jk (τ ) sin(ωτ )dτ

 R
ω ω 0
(9)
 B (ω ) = b + R ∞ K (τ ) cos(ωτ )dτ
jk jk 0 jk

In the infinite water depth, the analytical expression of linear incident wave potential
ΦI ( P; t) in the reference coordinate system o-xyz can be given as [31]

iη0 g υ[z−i( x cos α+y sin α)] iω. t

ΦI ( P; t) = e e e (10)
ω
where η0 is the amplitude of incident wave; α is wave propagation angle (α = π is head
waves); ω is the absolute frequency; υ = ω2 /g; ωe = ω − υU cos α is the encounter frequency.
The elevation ηI (t) of the incident wave at the origin o is given as

ηI (t) = η0 eiωe t (11)

Based on the impulse response function method [10], the incident wave velocity
∇ΦI ( P; t), diffraction velocity potential Φ7 ( P; t), and the diffraction force Fj7 (t) in the jth
mode (j = 1, 2, · · · , 6), K̂ ( P; t), Φ̂7 ( P; t), and K j7 (t) can be expressed in the following as
R∞
∇ΦI ( P; t) = −∞ K̂ ( P; t − τ )ηI (τ )dτ



 R∞
Φ7 ( P; t) = −∞ Φ̂7 ( P; t − τ )ηI (τ )dτ

(12)
 Fj7 (t) = −∞∞ K j7 (t − τ )ηI (τ )dτ

 R


where K̂ ( P; t) is the impulse response function of incident wave velocity ∇ΦI ( P; t), Φ̂7 ( P; t)
is the impulse response function of Φ7 ( P; t), and K j7 (t) is the impulse response function
of Fj7 (t).
In combination with Equations (2), (5) and (12), the K j7 (t) is given by

x 
∂Φ̂7 ( P; t)

K j7 (t) = ρ Φ̂7 ( P; t)m j − n j dS (13)
SB ∂t

2.4. Solving the Ship Motion Equations

Using Newton’s second law, the six degree of freedom motions of the rigid body in
space fixed coordinate system are determined by
n .. o
Mij ζ j (t) = Fi (t) (14)
J. Mar. Sci. Eng. 2021, 9, 87 6 of 19

6
where Fi (t) = FiI (t) + FiH (t) + Fi7 (t) + ∑ Fij (t), FiI (t) is F–K forces in the ith mode, FiH (t)
j =1
is hydrostatic forces in the ith mode, and Mij is the element of the general mass matrix.
In the body coordinate system G-xb yb zb , the six-degrees of motion equation can also
be given as
. .
mu + mΛ × u = Fb ; IΛ + Λ × IΛ = Mb (15)
where m is the mass
matrix, I is inertial moment matrix, u =
(u, v, w) is velocity vector,
Λ = Λ xb , Λyb ,
Λzb is the angular velocity, Fb = Fxb , Fyb , Fzb is force vector, and Mb =
Mxb , Myb , Mzb is moment vector.
The vectors between the space fixed coordinate system and the body coordinate
system can be transformed by a transformation matrix [33]. For linear problems, the trans-
formation matrix is the identity matrix.
The ship motion equations can be solved by various numerical methods, such
J. Mar. Sci. Eng. 2021, 9, x FOR PEER REVIEW 7 of as
22
Runge–Kutta method, predictor-correctors method, and so on. To avoid initial numerical
instability, a ramp function [34] is applied to solve the ship motion equations.
3. Numerical Methods
3. Numerical Methods
3.1. Analytical
3.1. Analytical Expression
Expression for
for F–K
F-K Forces Evaluation
Forces Evaluation
Consider two
Consider two coordinate
coordinatesystems
systemsillustrated
illustratedinin Figure
Figure 2: the reference
2: the coordinate
reference sys-
coordinate
tem o-xyz and the local panel coordinate system o ′ −0 x ′y
0 ′z0 ′ 0 defined
system o-xyz and the local panel coordinate system o -x y z defined by the vertices 1 to by the vertices 1 to 44
in the
in thecounterclockwise
counterclockwisedirection, namedP1 ,PP
direction,named 1 2 ,PP
, 2 3, ,Pand
3
, and P4 , respectively.
P4 , respectively. The The o 0 z0
positive
positive
o ′z ′ points
axis to the to
axis points exterior of the of
the exterior fluid.
the fluid.

z′
P4 y′
z P3
o′
P1 x′
P2
o
o y
x
0 0 y0 z0 .
Figure 2. The
Figure 2. The reference
reference coordinate systemo-xyz
coordinate system o-xyzand
andlocal
localpanel
panelcoordinate
coordinatesystem o′ − x′y′z′ .
systemo -x

The
The transformation
transformation matrix
matrixbetween
betweenthe
theposition vectorr in
positionvector r the reference
in the coordinate
reference coordi-
system oxyz and the position vector r0 in the local panel coordinate system o 0 -x 0 y0 z0 is
nate system oxyz and the position vector r′ in the local panel coordinate system
given as
o ′ − x ′y ′z ′ is given as r = r0 + Tr0 (16)
where ro = ( xo , yo , zo ) is the position vector r = r0 in
+ Tr ′
o-xyz, which is the origin of the panel (16)
0 0 0 0
coordinate o -x y z system, T is the unit transformation cosine-director matrix between
where or0o-x=0 y( x0 zo 0, yand
system ) is theo-xyz
o , z osystem position
and givenvectorby in o − xyz , which is the origin of the panel
coordinate o′ − x′y′z ′ system,  T Dis the E unit
D 0 transformation
E D 0 E  cosine-director matrix be-
0
tween system o ′ − x ′y ′z ′ and system o − xyz and given
x , x y , x z , byx
 D E D 0 E D 0 E 
0
T =  x , y′ y ,y z ,y  (17)
 
 D 0 xE, x D 0 y ′,Ex D z0 ′, xE  
 zz ′,,zy 
T =x , zx ′, y y ,yz′, y  (17)
 x ′, z y ′, z z0 ′, z0  0
where x, y, and z are unit base vectors in system o-xyz; x , y , and z are unit base vectors in
x ,0 y0yz,0 ;and
where o-x
system z are unit
h, i denotes the internal
base vectors product betweeno base
in system ; x ′ , y ′ , and z ′ are unit
− xyzvectors.
base vectors in system o − x ′y ′z ′ ; , denotes the internal product between base vectors.
Within linear dynamic conditions, the elevation η I ( t ) ( z ≤ 0 ) is given as

 pI ( t ) = ρ gη0 e cos ωe t − υ ( x cos α + y sin α ) 

υz

 (18)
 pH = − ρ gz

The pressure p IH can be given by

J. Mar. Sci. Eng. 2021, 9, 87 7 of 19

Within linear dynamic conditions, the elevation ηI (t) (z ≤ 0) is given as


pI (t) = ρgη0 eυz cos[ωe t − υ( x cos α + y sin α)]
(18)
pH = −ρgz

The pressure pIH can be given by

pIH = ρgηI (t)eυz − ρgz (19)

The resultant forces of F–K forces and hydrostatic forces acting on the mean wetted
hull surface in the jth mode can be given by
x
FIHj (t) = pIH (t)n j dS (20)
SB

The mean wetted hull surface can be discretized by N quadrilateral panels. For the ith
quadrilateral panel under the mean free surface, the F–K force FIi can be written as
x
FIi = ni ρgη0 eυz cos[ωe t − υ( x cos α + y sin α)]dS (21)
Si

where ni is the unit normal vector pointing outward the fluid, and Si is the area of the
quadrilateral panel.
The Equation (21) can be represented in the ith local panel coordinate system o 0 -x 0 y0 z0
x 0 0
eB( x ,y ) cos A x 0 , y0 dS


FIi = ni ρgη0 (22)
Si
h D 0 E D 0 E   D 0 E D 0 E  i
where A(x0, y0 ) = ωe t − υ xo + x , x x0 + y , x y0 cos α + yo + x , y x0 + y , y y0 sin α ,
 D 0 E D 0 E 
and B(x0, y0 ) = υ zo + x , z x0 + y , z y0 .
The jth ( j = 1, 2, 3, 4) edge for the ith quadrilateral panel can be parametrized by
   
gij (κ ) = x 0 0,ij , y0 0,ij + ∆x 0 ij , ∆y0 ij κ, κ ∈ [0, 1] (23)

where ∆x 0 ij = x 0 1,ij − x 0 0,ij , ∆y0 ij = y0 1,ij − y0 0,ij , x 0 0,ij = x 0 ij (0); x 0 1,ij = x 0 ij (1), y0 0,ij = y0 ij (0),
and y0 1,ij = y0 ij (1).
A, B, and their derivatives on the jth side of the ith panel can be expressed with respect
to parameter κ as
  h D 0 E
D 0 E i 

  xo + x , x x 0 0,ij + ∆x 0 ij κ + y , x y0 0,ij + ∆y0 ij κ cos α 
Aij (κ ) = ωe t − υ  h

 D 0 E D 0 E i 

 + y + x , y x 0 + ∆x 0 κ + y , x y0 + ∆y0 κ sin α 


o 0,ij ij 0,ij ij (24)

 h D 0 E  i
Bij (κ ) = υ zo + x , z x 0 0,ij + ∆x 0 ij κ + hy∗ , zi y0 0,ij + ∆y0 ij κ





Let Ψκ = A2κ + B2κ , and

(
eB (Bκ cos A+Aκ sin A)
Ψκ , Ψκ 6 = 0
CE = (25)
B
e κ cos A, Ψκ = 0
(
eB (Bκ sin A−Aκ cos A)
Ψκ , Ψκ 6 = 0
SE = (26)
eB κ sin A, Ψκ = 0
J. Mar. Sci. Eng. 2021, 9, 87 8 of 19

With respect to parameter κ, A and B’s derivatives on the ith panel are written as
 D 0 E D 0 E 

 A x 0 ,i = − υ x , x cos α + x , y sin α

 D 0 E D 0 E 

Ay0 ,i = − υ y , x cos α + y , y sin α (27)

 D 0 E D 0 E

 Bx0 ,i = υ x , z , By0 ,i = υ y , z


Let Ψ x0 = A2x0 + B2x0 and Ψy0 = A2y0 + B2y0 . The Green’s theorem with parameterization
in Equation (23) is applied, and integration by parts is adopted to solve Equation (22).
If Ψ x0 6= 0, Equation (22) is expressed as

ni ρgη0 4
Z 1

Ψ x0 j∑
∆yij0 eBij (κ ) Bx0 ,i cos Aij (κ ) + Ax0 ,i sin Aij (κ ) dκ
 
FIi = (28)
=1 0

∆CEij can be understood as CEij (1) − CEij (0) , and ∆SEij can be understood as
 

SEij (1) − SEij (0) . FIi is written as

ni ρgη0 4
∑ ∆yij∗ Bx0 ,i ∆CEij + Ax0 ,i ∆SEij


FIi = (29)
Ψ x 0 j =1

If Ψy0 6= 0, FIi is directly solved by

−ni ρgη0 4  
FIi = ∑
Ψ y 0 j =1
∆xij0 By0 ,i ∆CEij + Ay0 ,i ∆SEij (30)

Note that the analytical integration expressions for F–K moments, hydrostatic forces,
and hydrostatic moments can be solved in a similar approach.
With the established formulation for F–K forces evaluation, it is worth comparing the
present method with the method proposed by Rodrigues and Guedes Soares (2017) [16]
(see Table 1).

Table 1. Differences between the present method and Rodrigues’s method.

Present Method Rodrigues’s Method [16]

Froude–Krylov forces calculated in an
Froude–Krylov forces calculated and motions
Equation of motion inertial frame, and motions solved in a
solved in a moving system fixed to the ship.
moving system fixed to the ship.
Ship velocity arbitrary zero
Diffraction force direct time-domain method indirect time-domain method
Radiation force direct time-domain method indirect time-domain method
Expression of incident
cosine form sinusoidal form
wave elevation

3.2. Numerical Method for TFSGF Evaluation

The memory part of TFSGF in Equation (4) can be written in the non-dimensional form
q
e( P, Q; t − τ ) = 2
G g/r0 3 F (µ, β) (31)

where ν = υr 0 , µ = −(z + ζ )/r 0 , and β = g/r 0 (t − τ ). F (µ, β) is given by

p

Z ∞√ √

 q 
F (µ, β) = ν sin β ν e−νµ J0 ν 1 − µ2 dν (32)
0
J. Mar. Sci. Eng. 2021, 9, 87 9 of 19

F (µ, β) is the solution to the following fourth-order ODE

β2 + 16µ ∂2 F (µ, β) 7β ∂F (µ, β) 9F (µ, β)

∂4 F (µ, β) ∂3 F (µ, β)
+ µβ + + + =0 (33)
∂β4 ∂β3 4 ∂β2 4 ∂β 4

The fourth-order ODE Equation (33) can be given by

.
X
e ( β) = A
e ( β )X
e ( β) (34)

e ( β) = F, ∂F/∂β, ∂2 F/∂β2 , ∂3 F/∂β3 T , and

 
where X
 
0 1 0 0
 0 0 1 0 
A
e ( β) =  (35)


 0 0 0 1 

β2
 
− 49 − 47 β − 4µ + 4 −µβ

Equation (34) can be solved once the initial conditions are given.
β− β
For β ∈ [ β k , β k+1 ], let s = β −kβ (s ∈ [0, 1]); the relationship between the unit
k +1 k
time-variant system and time-variant system in Equation (34) can be written as

e ( β k + s·( β k+1 − β k )) = e
X x( s )
e ( β k + s·( β k+1 − β k )) = A(s) (36)
A

where A(s) = ∑2i=0 Ai si (Ai is the time-invariant coefficient matrix). The transformation
relationship is given by
dex/ds = ( β k+1 − β k )A(s)x(s) (37)
where the initial condition is x| s=0 = x(0).
The following equations can be obtained by using Equation (37)

x ( s ) = x( s )

 0

xi (s) = si x0 (s) (i = 1, 2, ..., m + 2; m > 2) (38)
 T
X(s) = xT0 xT1 xTm xTm+1 xTm+2
 
...

where m is an integer variable, and X(s) is 4 (m + 3) dimensional column vector.

The constant coefficient matrix M can be obtained by using the unit linear time-varying
coefficient matrix M(s) [18]. The elements of M on the principal diagonal are 1, and nonzero
elements of the principal diagonal of M are high power exponent of s. As m increases,
the power exponent of s can be high enough, and the coefficient matrix M tends to be
the unit constant coefficient matrix. The specific values of ‘m’ can be selected based on
convergence analysis for TFSGF evaluation, it’s fairly long, and it can refer to the paper by
Li [27]. Thus, it can be obtained from Equation (38) as follows
.
X(s) = MX(s) (39)

For the linear time-invariant differential equation (Equation (39)), it can be solved
by the 2 N algorithm in reference [26]. Once Equation (39) is solved, the TFSGF and its
derivatives can be easily obtained.
When µ = 0, the amplification of the oscillatory behavior of TFSGF should be noticed.
The analytical expression for the TFSGF is written as

πβ3
  2  2  2  2 
β β β β
F (0, β) = √ J1 ·J− 1 + J3 ·J− 3 (40)
16 2 4 8 4 8 4 8 4 8
 J. Mar.
J. Mar. Sci. Sci.
Eng. Eng. 2021,
2021, 9, 9,
87x FOR PEER REVIEW 11 of 22 10 of 19
J. Mar. Sci. Eng. 2021, 9, x FOR PEER REVIEW 11 of 22

4. Results and Discussions

4. Results
4. Resultsand
andDiscussions
Discussions
4.1. The Numerical Results of TFSGF
4.1. The
4.1. TheNumerical
Numerical Results of TFSGF
Results of TFSGF
In Figure 3, the “analytical method” denotes the solutions obtained by the analytical
In
InFigure
Figure 3,3,the “analytical
the method” denotes
“analytical the solutions obtained by the analytical
expression for TFSGF when μ = 0 method”
(presenteddenotes the solutions
in Equation obtained
(40)). Figure 3 shows by
expression for TFSGF when μ = 0 (presented in Equation (40)). Figure 3 shows the solu-
the the analytical
solu-
expression for TFSGF when = 0 (presented in Equation (40)).
tions obtained by the “PIM method”, with constant m = 50 for 0 ≤ β ≤ 100 [27].
µ Figure 3 shows the
tions obtained
solutions by theby
obtained “PIM
themethod”, with constant
“PIM method”, with m = 50 for m0 =
constant ≤β50≤for
1000 [27].
≤ β ≤ 100 [27].
100
100
Analytical method
Analytical method
PIM method
PIM method
50
50

0
0

-50
-50

-100
-100 0 20 40 60 80 100
0 20 40 60 80 100

Figure 3. The F ( 0, β ) computed by “precise integration method (PIM) method” with Δβ = 0.02 .
Figure The FF((0,0,ββ)) computed
Figure [Link] computedbyby“precise
“preciseintegration method
integration (PIM)(PIM)
method method” with Δ
method” β = 0.02
with ∆β .= 0.02.
In Figure 4, “present method” denotes solutions of TFSGF solved by PIM method,
In Figure
In Figure4,4,“present
“present method”
method” denotes
denotes solutions
solutions of TFSGF
of TFSGF solved solved
by PIM by PIM method,
method,
with variable m for 0 ≤ β ≤ 100 ; in the present study, m = 30 for 0 ≤ β ≤ 60 , m = 40 for
with variable
with variablemmfor β ≤β100
for00≤ ≤ ≤ ;100;
in the
in present study,study,
the present m = 30mfor
= 30 β ≤0 60
0 ≤for ≤, βm≤= 40
60,for
m = 40 for
60 < β ≤ 80 , and m = 50 for 80 < β ≤ 100 , respectively.
60 <
60 < ββ≤≤8080, and
, and mm = 50
= 50 for for < β<≤ 100
80 80 β ≤, respectively.
100, respectively.

100
100
Analytical method
Analytical method
Present method
Present method
50
50

0
0

-50
-50

-100
-100 0 20 40 60 80 100
0 20 40 60 80 100

FF((0,
Figure 4. The
Figure The F
Figure [Link] ( 0,0,βββ))) computed
computed by
computedby “present
by“present
“presentmethod” with Δβ ∆β
method”
method” withwith
= 0.02
= ..0.02.
Δβ = 0.02

Figure55shows
Figure shows absolute
absolute errors
errors computed
computed by “PIM
by “PIM method”
method” and “present
and “present method”,method”,
Figure 5 shows absolute errors computed by “PIM method” and “present method”,
with variablemmfor
with variable
J. Mar. Sci. Eng. 2021, 9, x FOR PEER REVIEW β
for00≤≤ ≤β100
≤ ;100; in
in the the present
present m = 30mfor
study,study, = 30 β
0 ≤for ≤0 60
≤ , βm≤ 60,
= 40
12
with variable m for 0 ≤ β ≤ 100 ; in the present study, m = 30 for 0 ≤ β ≤ 60 , m = 40 for
22m = 40 for
of for
60 < ββ≤≤8080,
60 < and
, and mm = 50
= 50 for for < β<≤ 100
80 80 β ≤, respectively.
100, respectively.
60 < β ≤ 80 , and m = 50 for 80 < β ≤ 100 , respectively.

10-10
1.5
PIM method
Present method

0.5

0
0 20 40 60 80 100

Figure
Figure5.5.
The absolute
The error
absolute σ σ
error computed by “present
computed method”
by “present and “PIM
method” method”
and “PIM with with ∆β = 0.02.
method”
Δβ = 0.02 .

In the present study, all computations are carried out on the platform Intel(R) Core
(TM) i7700HQ CPU 2.80 GHz. From Figures 3–5, the “PIM method” takes about 236.38 s
to generate TFSGF at 1× 5000 sets of ( μ , β ) for 0 ≤ β ≤ 100 at μ = 0 , while the “pre-
sent method” only takes about 85.8 s. The absolute errors of both “present method” and
−9
 In the present study, all computations are carried out on the platf
(TM) i7700HQ CPU 2.80 GHz. From Figures 3–5, the “PIM method” ta
to generate TFSGF at 1× 5000 sets of ( μ , β ) for 0 ≤ β ≤ 100 at μ = 0
J. Mar. Sci. Eng. 2021, 9, 87 11 of 19
sent method” only takes about 85.8 s. The absolute errors of both “pres
“PIM method” are within 10 −9 . Thus, the “present method” shows bett
ciency than the
In the “PIM
present method”.
study, all computations are carried out on the platform Intel(R) Core
(TM) i7700HQ CPU 2.80 GHz. From Figures 3–5, the “PIM method” takes about 236.38 s
to generate TFSGF at 1 × 5000 sets of (µ, β) for 0 ≤ β ≤ 100 at µ = 0, while the “present
4.2. Themethod”
StudyonlyCase
takesand
aboutParameter Setting
85.8 s. The absolute errors of both “present method” and “PIM
method” are within 10−9 . Thus, the “present method” shows better evaluation efficiency
The
thanhull form
the “PIM of Wigley I hull can be defined as in the reference
method”. [2
4.2. The Study Case and Parameter Setting
The y  of Wigley
hull form 2 x  I hull
 can bez defined
2
   as in the"reference
2
2 x   [29]# z 
2 2
  z 4  
=z 12 −   2x 12−   z 2  1 +
 z0.2
4   2x+2 4  1 −    1
"  2 #  " #
B1 −2  1 +L0.2   + D  1−  1L−  L  D   D   
     
y 2x
B/2
= 1−
L D L D D
 (41)
L is the length of the ship; B is the width of the ship; D is the draft of the ship.
wherewhere
L The ismain
the particulars
length of the ship; B is the width of the ship; D is th
of the Wigley I hull are given in Table 2; k yy denotes the pitch
The main
inertia radiusparticulars
of the ship, and of the Wigley
∇ denotes I hullvolume
the displacement are given in Table 2; k yy
of the ship.

inertiaTable
radius
2. Mainof the ship,
particulars ∇ denotes the displacement volume of
and I hull.
of the Wigley

The meshes
Ship of the
L/m ship’s B/m hull are illustrated
D/m in
kyyFigure 6.
∇ /m3
Wigley I hull 3.0 0.3 0.1875 0.25 L 0.0946

The meshes of the ship’s hull are illustrated in Figure 6.

Figure 6. The panels’ distribution on the Wigley I hull.

Non-dimensional forms of hydrodynamic coefficient, memory function, wave exciting

forces, and motion response are given below.
The non-dimensional added mass coefficients are defined as A33 0 = A / ( ρ ∇) and
33
Figure A6.0 The
= A panels’
55 / ρ ∇ L 2

distribution
, respectively. on the Wigley I hull.
55
0 =
p
The non-dimensional damping coefficients are defined as B33 L/gB33 /(ρ∇) and
0 = L/gB55 / ρ∇ L2 , respectively.
p

B55
Table 2. Main particulars of the Wigley I hull.
0 and pitch 0 are defined as K 0 =
The heave memory function K33 memory function K55 33
0
(K33 L)/(ρg∇) and K55 = K55 /(ρg∇), respectively.
ShipThe non-dimensional L /mwave exciting forcesB /m 0 D=/m
are defined as F3w F3w L/(ρgη0 ∇) and k yy
0
F = F5w /(ρgη0 ∇), respectively.
Wigley5wI The hullnon-dimensional 3.0 heave and pitch 0.3motion are defined0.1875
as ζ 30 = ζ 3 /η0 and ζ 50 0.25L
=
(ζ 5 L)/(2πη0 ), respectively.
The non-dimensional frequency ω 0 is defined as ω 0 = ω L/g, the non-dimensional
p
Non-dimensional
wave number υ0 is defined forms of and
as υ0 = υL, hydrodynamic
the non-dimensional coefficient, memory
frequency t0 is defined as fun
0
t = t/Te (Te is encounter period).
ing forces, and motion response are given below.
4.3. The Convergence Study on Panel Number N and Time Step ∆t
The panel number N and the time step ∆t have influences on the numerical results of
ship hydrodynamic analysis and ship motion. The convergence of the two parameters is
verified by taking the Wigley I hull as a study case. The error of 0.1% is adopted to
ptest
convergence. Wigley I hull is studied at Fn = 0.2 (Fn denotes Froude number, Fn = U/ gL)
in head seas (α = π, η0 = 0.036 m).
 e e

4.3. The Convergence Study on Panel Number N and Time Step Δt

The panel number N and the time step Δt have influences on the numerical results
of ship hydrodynamic analysis and ship motion. The convergence of the two parameters
J. Mar. Sci. Eng. 2021, 9, 87 12 of 19
is verified by taking the Wigley I hull as a study case. The error of 0.1% is adopted to test
convergence. Wigley I hull is studied at Fn = 0.2 (Fn denotes Froude number,
Fn = U gL ) in head seas ( α = π , η0 = 0.036m ).
The 22×
The × 22 Gaussian
Gaussian quadrature [15] can
quadrature [15] canbe
beused
usedtotocompute
compute thethe F–K
F-K forces
forces andand
hy-
hydrodynamic
drodynamic forces acting on the panel. As shown in Figure 7, the square panel can be
forces acting on the panel. As shown in Figure 7, the square panel can be
progressively
progressivelysubdivided
subdivideduntil
untilthe
theprescribed
prescribedprecision
precisionisisreached.
reached.

Figure 7.
Figure 7. The
The panel
panel subdivision
subdivision diagram.
diagram.

The
The four different
different panel
panelnumbers
numbersofofthe
the Wigley
Wigley I hull
I hull surface
surface discretization
discretization withwith
con-
constant time
stant time stepstep ∆tTe=40
Δt = Te /40
are are carried
carried out out (Table
(Table 3). 3).

[Link]
Table Panelnumbers
numbersfor
forthe
theconvergence
convergenceanalysis
analysisofofWigley
WigleyI Ihull
hulldiscretization.
discretization.

Case
Case 11 22 33 44
Wigley
Wigley I I NN== 20
20×× 44 NN== 40
40×× 44 NN==60 × 44
60× NN==80
80×× 44
J. Mar. Sci. Eng. 2021, 9, x FOR PEER REVIEW 14 of 22
J. Mar. Sci. Eng. 2021, 9, x FOR PEER REVIEW 14 of 22
Figure 88 illustrates
Figure illustrates the
the time
timehistory
historyof
ofthe
theheave
heavememory
memoryfunction and
function pitch
and memory
pitch mem-
function.
ory function.
4
4 N = 80 0.4 N = 80
N = 160
80 0.4 N = 160
80
3 0.3
3 160
N = 240 160
N = 240
N = 320
240 0.3 N = 320
240
2 0.2
N = 320 N = 320
2 0.2
0.1
1 0.1
1 0
0 0
0 -0.1
-0.1
-1 -0.2
-1 -0.2
0 2 4 6 8 0 2 4 6 8
0 2 t4 6 8 0 2 t4 6 8
(a)t (b)t
(a) (b)
Figure 8. Non-dimensional memory function of Wigley I hull: (a) heave; (b) pitch.
[Link]-dimensional
Figure Non-dimensionalmemory
memoryfunction
functionofofWigley
Wigley I hull:(a)(a)heave;
I hull: heave;(b)
(b)pitch.
pitch.
Figure 9 shows the time history of the non-dimensional F-K forces. The Wigley I hull
Figure
Figure
advances 9 9shows
in shows
head thetime
the
regulartime history
history
waves ofofthe
at Fn =the non-dimensional
non-dimensional
0.2, F-Kforces.
F–K
and the wavelength forces.
to The
shipThe Wigley
Wigley
length ILhull
is λI hull
=1
advances
advancesininhead
headregular
regularwaves
waves FnFn
atat = =0.2,
0.2,and the
and wavelength
the toto
wavelength ship
shiplength is λ =
is λ/L
length L=1.1
. F-K forces can be obtained by the Gaussian quadrature method.
F–K forces can be obtained by the Gaussian quadrature method.
. F-K forces can be obtained by the Gaussian quadrature method.

4
4 N = 80 2 N = 80
3.82 N = 160
80 2 N = 160
80
2.14
2 3.82 160
N = 240 160
N = 240
3.81 2.14
2 N = 320
240 1 2.12 N = 320
240
3.81 1
3.8 N = 320 2.128.2 8.25 8.3 N = 320
0 3.8 7 0 8.2 8.25 8.3
0 7 0
-1
-2 -1
-2
-2
-4 -2
-40 5 10 15 0 5 10 15
0 5 t 10 15 0 5 t 10 15
(a)t (b)t
(a) (b)
[Link]-dimensional
Figure F-Kforces
Non-dimensionalF–K forcesofofWigley
WigleyI Ihull ( λ L==1:
hull(λ/L 1 ):(a)(a)heave;
heave;(b)
(b)pitch.
pitch.
Figure 9. Non-dimensional F-K forces of Wigley I hull ( λ L = 1 ): (a) heave; (b) pitch.
As shown in Figures 8 and 9, the numerical results tend to be convergent when
As shown
N ≥ 240 in relative
, and the Figures error
8 andbetween
9, the numerical
N = 240 results
and N tend
= 320toisbewithin
convergent when
0.1%. Thus,
N ≥ 240
panel , and N
number the= 320
relative error between
is selected N = 240 and analysis
for ship hydrodynamic N = 320 and
is within 0.1%. in
ship motion Thus,
the
J. Mar. Sci. Eng. 2021, 9, 87 13 of 19

As shown in Figures 8 and 9, the numerical results tend to be convergent when

N ≥ 240, and the relative error between N = 240 and N = 320 is within 0.1%. Thus,
panel number N = 320 is selected for ship hydrodynamic analysis and ship motion in the
present study.
In the system o-xyz, a square panel is adopted as a study case advancing in head regu-
lar waves, whose four vertices are P1 (−0.5, 0, −1), P2 (0.5, 0, −1), P3 (0.5, 0, 0),
and P4 (−0.5, 0, 0), respectively. The square panel advances in head regular waves at
Fn = 0.2, and the wavelength to ship length is 1.5, and the wave steepness is 0.024.
At t = 0, the F–K force acting on the square panel by the analytical integration method
is 179.6193
J. Mar. Sci. Eng. 2021, 9, x FOR PEER REVIEW N (the calculation accuracy can be set to 0.0001 N). 15 of 22
From Table 4, as the subdivision time is 4, the square panel can be subdivided into
256 smaller subpanels, as shown in Figure 7. However, the computational results of the
Gauss integration method and the analytical integration method are the same. For simple
Table 4. Values of F-K forces acting on the panel at t = 0.
geometry of the ship’s hull, the magnitude of panel number for hull mesh generation
can be as low as O(10), such as barge vessel [16]. The analytical integration method can
Subdivision
improveTimethe F–K forces0evaluation accuracy
1 with no 2loss of accuracy.3 4

F-K forces
Table 4. Values of F–K179.3024
forces acting on179.5996
the panel at t = 179.6181
0. 179.6192 179.6193
(unit: N)
Non-dimen-
Subdivision Time 0 1 2 3 4
sional F-K 15.7080
F–K forces (unit: N)
15.7340
179.3024
15.7357
179.5996 179.6181
15.7358
179.6192
15.7358
179.6193
forces
Non-dimensional F–K forces 15.7080 15.7340 15.7357 15.7358 15.7358

Figure 10 shows the time history of the non-dimensional F-K forces in head regular
Figure 10 shows the time history of the non-dimensional F–K forces in head regular
waves at Fn = 0.2 obtained by “numerical method” and “present method”, respectively.
waves at Fn = 0.2 obtained by “numerical method” and “present method”, respectively.
The Wigley I hull surface is discretized and processed into N = 320 quadrilateral panels.
The Wigley I hull surface is discretized and processed into N = 320 quadrilateral panels.
“numerical method” denotes the numerical results of F-K forces obtained by 2 × 2 Gauss-
“numerical method” denotes the numerical results of F–K forces obtained by 2 × 2 Gaussian
ian quadrature, and “present method” denotes the results of F-K forces obtained by the
quadrature, and “present method” denotes the results of F–K forces obtained by the
analytical method.
analytical method.

4 2.5
Present method 2 Present method
Numerical method Numerical method
1.5
2
1
0.5
0 0
-0.5
-1
-2
-1.5
-2
-4 -2.5
0 5 10 15 0 5 10 15
t t
(a) (b)
Figure10.
Figure [Link]
Timehistory
historyofofnon-dimensional F-Kforces
non-dimensionalF–K forcesofofWigley
WigleyI Ihull
hullobtained
obtainedbyby“present
“presentmethod”
method”and
and“numerical
“numerical
method” ( λ L = 1 ): (a) heave; (b)
method” (λ/L = 1: (a) heave; (b) [Link].

From
FromFigure
Figure 10,10,
the the
F–K F-K
forcesforces
obtained by “numerical
obtained method”method”
by “numerical and “presentand method”
“present
are almost the same, and the relative error between the “numerical method”
method” are almost the same, and the relative error between the “numerical method” and “present
and
method”
“presentismethod”
0.1%. Theis CPU0.1%.time
Theconsumption for the “numerical
CPU time consumption for themethod” is about
“numerical 35.4 s,is
method”
while
aboutthe CPU
35.4 time consumption
s, while the CPU timefor the “present
consumption formethod” is about
the “present 27.8 s. isThe
method” “present
about 27.8 s.
method”
The “present method” is more efficient than the “numerical method”. The efficiency the
is more efficient than the “numerical method”. The efficiency and accuracy of and
“present
accuracy method” can be validated.
of the “present method” can be validated.
Figure
Figure11 11illustrates
illustratesthethetime
timehistory
historyofofheave
heavemotion
motionandandpitch
pitchmotion
motionofofWigley
WigleyI I
hull with four different time steps, which are Te /10, T /20, T /30, and T /40, respectively.
hull with four different time steps, which are Te 10e, Te 20e , Te 30 , and e Te 40 , respec-
tively.
 J. Mar.
J. Mar. Sci.
Sci. Eng.
Eng. 2021,
2021, 9, 9,
87x FOR PEER REVIEW 16ofof1922
14

1.5 1
1.008 t = Te /10 t = Te /10
1.006 0.638
1 1.004 t = Te /20 t = Te /20
0.636
1.002 0.5
1 t = Te /30 0.634 t = Te /30
0.5 0.632
0.998
t = Te /40 t = Te /40
7.25 7.3 10.5 10.55
0 0

-0.5
-0.5
-1

-1.5 -1
0 5 10 15 0 5 10 15
t t
(a) (b)
Figure
Figure [Link]
Time history
history ofof motions
motions ofof Wigley
Wigley I hull
I hull (λ =
(λ/L L =2,2N, N = 320):
= 320): (a)(a) heave;
heave; (b)(b) pitch.
pitch.

InInFigure
Figure11,11,
asas
time step
time ∆t decreases,
step a convergent
Δt decreases, trendtrend
a convergent can be
canobtained, and the
be obtained, andrel-
the
ative error between ∆t = Te /30 and ∆t = T /40 is within 0.1%. The time step ∆t = T /40
relative error between Δt = Te 30 and Δt = Te 40 is within 0.1%. The timee step
e
is selected for the subsequent ship hydrodynamic analysis and ship motions.
Δt = Te 40 is selected for the subsequent ship hydrodynamic analysis and ship motions.
4.4. Added Mass and Damping Coefficient
[Link]
Addedthe
Mass and Dampingcoefficients
hydrodynamic Coefficient on the main diagonal have a significant influence
on shipSince the hydrodynamic coefficients on the main diagonal 0 , A0 , influ-
have a significant
A33 0
hydrodynamic analysis and motion, the non-dimensional coefficients 55 B33 ,
and
ence 0
B55onare selected
ship for detailed
hydrodynamic numerical
analysis motion, the non-dimensional coefficients A33
andanalysis. ′ ,
A55′ , B33′ , and B55′ are selected for detailed numerical analysis.
For the Wigley I hull at Fn = 0.2, Figure 12 presents the comparisons of the non-
dimensional heave-heave added masses and damping coefficients, and Figure 13 presents
For the Wigley I hull at Fn = 0.2, Figure 12 presents the comparisons of the non-di-
the comparisons of the non-dimensional pitch-pitch added masses and damping coeffi-
mensional heave-heave added masses and damping coefficients, and Figure 13 presents
cients. In the legend of Figures 12 and 13, “present method” denotes numerical results
the comparisons of the non-dimensional pitch-pitch added masses and damping coeffi-
obtained by the TFSGF method without waterline terms in Equation (5). “experiment”
cients. In the legend of Figures 12 and 13, “present method” denotes numerical results
denotes previous literature experimental data [29] obtained by Journée, and the experi-
obtained by the TFSGF method without waterline terms in Equation (5). “experiment”
ments were carried out in the Shiphydromechanics Laboratory of the Delft University of
denotes previous literature experimental data [29] obtained by Journée, and the experi-
Technology. The experimental data on hydrodynamic coefficients, wave loads, and added
ments were carried out in the Shiphydromechanics Laboratory of the Delft University of
resistance for heave and pitch motions in head waves of Wigley I hull are sufficient.
Technology.
“Magee’s The experimental
method” data on
denotes numerical hydrodynamic
results obtained bycoefficients,
the TFSGF wave
method loads,
[28] and added
in which
resistance for heave and pitch motions in head waves of Wigley I hull
the waterline integral terms of boundary integral equations are included. “LAMP method” are sufficient.
“Magee’s
denotes
J. Mar. Sci. Eng. 2021, 9, x FOR PEER REVIEWnumericalmethod” denotes
results numerical
obtained by LAMP results obtained
software [21],by
in the
whichTFSGF method
the body [28]
nonlinear
17 of in
22
which the waterline integral terms of boundary integral equations are
method was adopted where the perturbation potential was computed on the instantaneous included. “LAMP
method”
wetted hulldenotes
under thenumerical
mean freeresults obtained by LAMP software [21], in which the body
surface.
nonlinear method was adopted where the perturbation potential was computed on the
2 instantaneous wetted hull under the mean free surface.
Present method Present method
Experiment 1.5 Experiment
Magee’s method Magee’s method
1.5 LAMP method LAMP method
1

1
0.5

0.5
0
0 2 4 6 8 0 2 4 6 8

(a) (b)
Figure12.
[Link]-dimensional
Non-dimensionalheave-heave
heave-heaveadded
addedmass
massand
anddamping
dampingcoefficients: ′ ; (b) 0 B′ .
coefficients:(a)(a)A0 A;33
Figure 33 (b) B33 . 33

0.25 0.1
Present method Present method
Experiment Experiment
0.2 0.08
Magee’s method Magee’s method
LAMP method LAMP method
0.15 0.06
 0.5
0
0 2 4 6 8 0 2 4 6 8

J. Mar. Sci. Eng. 2021, 9, 87

(a) (b) 15 of 19
′ ; (b) B33
Figure 12. Non-dimensional heave-heave added mass and damping coefficients: (a) A33 ′ .

0.25 0.1
Present method Present method
Experiment Experiment
0.2 0.08
Magee’s method Magee’s method
LAMP method LAMP method
0.15 0.06

0.1 0.04

0.05 0.02

0 0
0 2 4 6 8 0 2 4 6 8

(a) (b)
Figure13.
[Link]-dimensional
Non-dimensionalpitch-pitch
pitch-pitchadded
addedmass
massand
anddamping
dampingcoefficients:
coefficients:(a) 0A′; (b)
(a)A55 0 B ′ .
Figure 55 ; (b)
B55 . 55

Tables
Tables5–8 present
5–8 absolute
present relative
absolute errors
relative of hydrodynamic
errors of hydrodynamic coefficients for Wigley
coefficients I hull I
for Wigley
byhull
“present method” and “Magee’s method”, which are compared with previous
by “present method” and “Magee’s method”, which are compared with previous literature

ex-
lit-
perimental data [29]. The absolute relative error can be written
erature experimental data [29]. The absolute relative error can as ε = N − N
num be exp /N
writtenexp as

( )
(ε denotes absolute relative error, Nnum denotes numerical results obtained by various models,
ε = N denotes
and Nexpnum
− Nexp a value ε previous
Nexp (for denotes absolute
literature relative
Nnum
error,data
experimental denotes numerical re-
[29]).
sults obtained by various models, and N exp denotes a value for previous literature exper-
Table 5. Absolute
imental relative errors of heave-heave added mass A33 for Wigley I hull by various methods.
data [29]).
0
ω 2.2 2.8 3.3 3.9 4.4 5.5
Table 5. Absolute relative errors of heave-heave added mass A33 for Wigley I hull by various
Present method
methods. 31.7% 3.7% 2.8% 1.4% 0.2% 2.5%
Magee’s method 40.6% 12.2% 11.0% 2.4% 0.2% 5.6%
ω′ 2.2 2.8 3.3 3.9 4.4 5.5
Present method 31.7% 3.7% 2.8% 1.4% 0.2% 2.5%
Table 6. Absolute
Magee’s methodrelative
40.6%errors of12.2%
heave-heave 11.0% 2.4% B33 for0.2%
damping coefficients Wigley I hull
5.6%by
various methods.
Table 6. Absolute
0 relative errors of heave-heave damping coefficients B33 for Wigley I hull by
ω 2.2 2.8 3.3 3.9 4.4 5.5
various methods.
Present method 0.1% 3.9% 4.2% 1.9% 10.6% 33.5%
ω′ 2.2 2.8 3.3 3.9 4.4 5.5
Magee’s method 11.5% 5.1% 1.0% 4.8% 4.2% 23.9%
Present method 0.1% 3.9% 4.2% 1.9% 10.6% 33.5%
Magee’s method 11.5% 5.1% 1.0% 4.8% 4.2% 23.9%
Table 7. Absolute relative errors of pitch-pitch added mass A55 for Wigley I hull by various methods.
0
ω 2.2 2.8 3.3 3.9 4.4 5.5
Present method 9.9% 2.3% 3.3% 17.1% 23.3% 28.5%
Magee’s method 57.0% 34.5% 32.1% 35.0% 40.0% 43.7%

Table 8. Absolute relative errors of pitch-pitch damping coefficients B55 for Wigley I hull by vari-
ous methods.
0
ω 2.2 2.8 3.3 3.9 4.4 5.5
Present method 8.2% 9.7% 4.1% 2.6% 2.1% 14.0%
Magee’s method 14.6% 14.3% 14.7% 16.0% 11.0% 18.07%

From Figures 12 and 13, resonance occurs around the non-dimensional frequency
ω 0 = 1.40, and there is a larger deviation between numerical results and previous literature
experimental data [29]. As the non-dimensional frequency increases, hydrodynamic coeffi-
cients obtained from “present method”, “Magee’s method”, and “LAMP method” show
 ious methods.

ω′ 2.2 2.8 3.3 3.9 4.4 5.5

Present method 8.2% 9.7% 4.1% 2.6% 2.1% 14.0%
Magee’s method 14.6% 14.3% 14.7% 16.0% 11.0% 18.07%
J. Mar. Sci. Eng. 2021, 9, 87 16 of 19
From Figures 12 and 13, resonance occurs around the non-dimensional frequency
ω ′ = 1.40 , and there is a larger deviation between numerical results and previous litera-
ture experimental data [29]. As the non-dimensional frequency increases, hydrodynamic
coefficients
a similar change obtained
trend. from
From“present method”,
Tables 5–8, “Magee’s method”,
the numerical and “LAMP
results obtained method”
by the “present
show a similar change trend. From Tables 5–8, the numerical results
method” are in better agreement with the previous literature experimental data obtained by the “pre-
[29] than
sent method” are in better agreement with the previous literature experimental data [29]
those obtained by “Magee’s method”. Based on the potential flow assumption, the viscosity
than those obtained by “Magee’s method”. Based on the potential flow assumption, the
of a fluid is ignored. When both the field point and source point are on the mean free
viscosity of a fluid is ignored. When both the field point and source point are on the mean
surface, the amplitude
free surface, and oscillating
the amplitude frequency
and oscillating frequencyof TFSGF
of TFSGFincrease
increaserapidly
rapidly with time pa-
with time
rameter β, and there
parameter β , andwill be awill
there larger
be anumerical error iferror
larger numerical waterline integral
if waterline terms terms
integral of boundary
of
integral
boundary integral equations are included. The accuracy for hydrodynamic coefficient can
equations are included. The accuracy for hydrodynamic coefficient evaluation
be improved
evaluationby canthe
be present
improved method with waterline
by the present method with integral terms
waterline excluded.
integral terms excluded.

4.5. The Wave

4.5. The Exciting
Wave Forces
Exciting Forces
The The
Wigley I hull advances
Wigley I hull advancesin head regular
in head waves
regular at Fn
waves = 0.2.
at Fn TheThe
= 0.2. wave
waveexciting
excitingforces
are composed of F–K forces and diffraction forces in the frequency domain,
forces are composed of F-K forces and diffraction forces in the frequency domain, andand diffraction
diffraction
forces forces in the
in the frequency frequency
domain candomain can beby
be obtained obtained
Fourier bytransformation
Fourier transformation [15]. 14
[15]. Figure
Figure 14 shows the numerical results of amplitudes of the non-dimensional
shows the numerical results of amplitudes of the non-dimensional wave exciting forces wave excit-
ing forces
obtained obtainedmethods.
by various by various methods.
”SAMP ”SAMP denotes
method” method” numerical
denotes numerical
results results
obtainedob- by a
tained by a transient free surface Green function method in the linear
transient free surface Green function method in the linear time domain [21]. time domain [21].

20 2.5
Present method Present method
Experiment Experiment
2
15 Magee’s method Magee’s method
SAMP method SAMP method
1.5
10
1

5
0.5

0 0
0 5 10 15 0 5 10 15

(a) (b)
J. Mar. Sci. Eng. 2021, 9, x FOR PEER REVIEW 19 of 22

FigureFigure 14. Non-dimensional

14. Non-dimensional wavewave exciting
exciting forceamplitude
force amplitude at
at Fn
Fn ==0.2:
0.2:(a)
(a)heave;
heave;(b)(b)
pitch.
pitch.

From Figure
From 14,14,the
Figure themaximum relative
maximum relative error
error between
between “present
“present method”method” and “ex-
and “exper-
periment”
iment” within 10.0%, and the non-dimensional wave exciting forces obtained from “pre-from
is within 10.0%, and the non-dimensional wave exciting forces obtained
“present
sent method”, “Magee’smethod”
method”, “Magee’s method” and
and ”SAMP
”SAMP method”
method” showshow a similar
a similar change
change trend.
trend.
Thus,
Thus, the reliability
the reliability of of analyticalexpression
analytical expression for
for F-K
F–Kforces
forcesevaluation
evaluationis verified.
is verified.

4.6. The
4.6. The Motion
Motion Responses
Responses
Figure
Figure 15 shows
15 shows thethe time
time history
history ofofmotions
motionsofofWigley
WigleyIIhull
hull at
at Fn
Fn == 0.2
0.2in
inhead
headreg-
regular
ular
waves = π).( α = π ).
(αwaves

1.2 1.2

0.8 0.8

0.4 0.4

0 0

-0.4 -0.4

-0.8 -0.8

-1.2 -1.2
0 5 10 15 0 5 10 15
t t
(a) (b)

FigureFigure 15. Non-dimensional

15. Non-dimensional timetime history
history of of motionsofofWigley
motions λ L ==1.5
Wigley II hull ((λ/L ): (a)
1.5): (a)heave;
heave;(b)(b)
pitch.
pitch.

From Figure 15, the instantaneous influences of the initial disturbance on the time
history of motion disappear after about 3~5 wave periods. Thus, the three-dimensional
linear time-domain simulation program developed in the present study is numerically
stable.
Figure 16 shows heave response amplitude operators (Raos) and pitch response am-
 0 5 10 15 0 5 10 15
t t
(a) (b)
Figure 15. Non-dimensional time history of motions of Wigley I hull ( λ L = 1.5 ): (a) heave; (b) pitch.
J. Mar. Sci. Eng. 2021, 9, 87 17 of 19

From Figure 15, the instantaneous influences of the initial disturbance on the time
history of motion disappear after about 3~5 wave periods. Thus, the three-dimensional
linear
Fromtime-domain simulation
Figure 15, the program
instantaneous developed
influences of theininitial
the present studyonisthe
disturbance numerically
time his-
stable.
tory of motion disappear after about 3~5 wave periods. Thus, the three-dimensional linear
Figure 16
time-domain shows heave
simulation response
program amplitude
developed in theoperators (Raos)isand
present study pitch response
numerically [Link]-
plitude
Figureoperators for ship
16 shows heave motions in head
response regularoperators
amplitude waves at Fn = [Link]
(Raos) The pitch
response ampli-
response
tude of heave
amplitude motion
operators forcan be motions
ship defined as ζ 30 η0regular
in head , and ζwaves
30 is the amplitude
at Fn = 0.2. The of ship
responseheave
motions. The
amplitude response
of heave amplitude
motion can be of pitch motion
defined can
as ζ 30 /η be defined
0 , and (ζ 50 L ) ( 2πηof0 )ship
as amplitude
ζ 30 is the , and
ζ 50 is the amplitude of ship pitch motions. In Figure 16, “Kim’s method” denotes 0the
heave motions. The response amplitude of pitch motion can be defined as ( ζ 50 L ) / ( 2πη ),
and ζ 50 is the amplitude of ship pitch motions. In Figure 16, “Kim’s method” denotes the
numerical results obtained by a three-dimensional Rankine panel method in the linear
numerical results obtained by a three-dimensional Rankine panel method in the linear time
time domain [35].
domain [35].

1.2 2
Experiment Experiment
Present method Present method
Kim’s method 1.6 Kim’s method
1

1.2

0.8
0.8

0.6 0.4
1.2 1.7 2.2 1.2 1.7 2.2
/L /L
(a) (b)
Figure16.
Figure [Link]
Responseamplitude
amplitudeoperators
operatorsofof motion
motion for
for Wigley
Wigley I hull
I hull atat
FnFn
= =0.2
0.2ininhead
headregular
regularwaves:
waves:(a)(a)heave;
heave;(b)
(b)pitch.
pitch.

From Figure 16, for ship heave motion responses, both “present method” and “Kim’s
method” can be in good agreement with previous literature experimental data [29]; for ship
pitch motion responses, the numerical results obtained by “present method” are closer
to previous literature experimental data [29] than those obtained by “Kim’s method”.
When λ/L approaches 1.80, the numerical results of the pitch motion responses show
non-ignorable errors compared with the previous literature experimental data [29]. As
λ/L approaches 1.80, the wave encounter frequency is nearly equal to the natural fre-
quency of pitch motion. Thus, the resonance can take place, and the previous literature
experimental data [29] is quite larger than the numerical results. In the present study, the
TFSGF method is adopted to solve ship motion problems. The TFSGF can automatically
satisfy the free surface boundary condition, and the numerical errors involved in radiation
boundary conditions can be reduced. Thus, the numerical results obtained by the “present
method” can be in better agreement with previous literature experimental data [29] than
“Kim’s method”.

5. Conclusions
In the present study, a three-dimensional time-domain panel method is developed
to study the ship’s hydrodynamic analysis and motions in regular waves. The Wigley I
hull is taken as a study case. The following conclusions can be made based on numerical
simulation and investigations:
(1) The precise integration method with variable parameter m is adopted for TFSGF
evaluation, which can improve the efficiency and numerical stability. It can provide a
reliable solver for a ship’s hydrodynamic analysis.
(2) Based on the TFSGF method, the boundary integral equation without waterline terms
is established to solve the perturbation velocity potential. When µ = 0, the violent
oscillation and amplitude amplification characteristics of TFSGF could lead to worse
numerical calculation results. The numerical results of hydrodynamic coefficients
obtained by the present method can be in good agreement with previous literature
J. Mar. Sci. Eng. 2021, 9, 87 18 of 19

experimental data [29]. In comparison with the TFSGF method, including waterline
terms, the present method shows higher accuracy.
(3) The derived analytical integration expressions for F–K forces evaluation over quadri-
lateral panels have been proved to provide exact results. For a much simple hull
shape, like a barge vessel, only about ten quadrilateral panels are required to dis-
cretize the hull body, which needs much fewer mesh grids than the Gauss integration
method. The wave exciting forces of the Wigley I hull in regular head waves are in
good agreement with both previous literature experimental data [29] and numerical
results by other published results [21,28]. Thus, the algorithm developed for F–K
forces can be validated.
(4) Since TFSGF can automatically satisfy the free surface condition, the ship pitch
motion response results obtained by the present method are in better agreement with
previous literature experimental data [29] than the three-dimensional Rankine panel
method [35].
(5) Based on the present research work, different levels of nonlinearity can be considered
in future study work. The boundary integral equation can be built on the instanta-
neous wetted hull surface instead of the mean wetted hull surface, in which the body
nonlinearity can be incorporated more fully. Moreover, the second-order drift forces
can also be derived from the present boundary value problem formulation, which
can be used to study the ship maneuvering problem.

Author Contributions: Methodology, P.Z.; formal analysis, P.Z. and T.Z.; Investigation, P.Z.; writing—
original draft preparation, P.Z.; writing—review and editing, P.Z., T.Z., and X.W. All authors have
read and agreed to the published version of the manuscript.
Funding: This research is supported in part by the National Natural Science Foundation of China
(Grant. 51909022, 61976033), the Natural Science Foundation of Liaoning Provence (Grant. 2019-BS-
024), and the Fundamental Research Funds for the Central Universities (Grant. 3132019347).
Data Availability Statement: The data used to support the findings of this study are available from
the corresponding author upon request.
Acknowledgments: The authors gratefully acknowledge the financial support from the National
Natural Science Foundation of Chin.
Conflicts of Interest: The authors declare no conflict of interest.

References
1. Subramanian, R. A Time Domain Strip Theory Approach to Predict Maneuvering in a Seaway. Ph.D. Thesis, The University of
Michigan, Ann Arbor, MI, USA, 2012.
2. Ogilvie, T.E.; Tuck, E.O. A Rational Strip Theory for Ship Motions; Report; University of Michigan: Ann Arbor, MI, USA, 1969.
3. Tasai, F. On the swaying, yawing and rolling motions of ships in oblique waves. Int. Shipbuild. Prog. 1967, 14, 1–20. [CrossRef]
4. Salvesen, N.; Tuck, E.O.; Faltinsen, O. Ship motions and sea loads. Trans. Soc. Naval Archit. Mar. Eng. 1970, 78, 250–287.
5. Fonseca, N.; Guedes Soares, C. Time domain analysis of large amplitude vertical motions and wave loads. J. Ship Res. 1998, 42,
100–113. [CrossRef]
6. Fonseca, N.; Soares, C.G. Comparison of numerical and experimental results of non-linear wave induced vertical ship motions
and loads. J. Mar. Sci. Technol. 2002, 6, 193–204. [CrossRef]
7. Tavakoli, S.; Niazmand, R.; Mancini, S.; De Luca, F.; Dashtimanesh, A. Dynamic of a planing hull in regular waves: Comparison
of experimental, numerical and mathematical methods. Ocean Eng. 2020, 217, 107959. [CrossRef]
8. Nakos, D.E. Ship Wave Patterns and Motions by a Three-Dimensional Rankine Panel Method. Ph.D. Thesis, Massachusetts
Institute of Technology, Cambridge, MA, USA, 1990.
9. Kring, D.C. Time Domain Ship Motions by a Three-Dimensional Rankine Panel Method. Ph.D. Thesis, Massachusetts Institute of
Technology, Cambridge, MA, USA, 1994.
10. Chen, J.P.; Zhu, D.X. Numerical simulations of wave-induced ship motions in time domain by a Rankine panel method.
J. Hydrodyn. 2010, 22, 373–380. [CrossRef]
11. Wehausen, J.V.; Laitone, E.V. Surface Waves. Encyclopedia of Physics, Vol. IX/Fluid Dynamics III; Springer: Berlin, Germany, 1960.
12. Blandeau, F.; Francois, M. Linear and non-linear wave loads on FPSOs. In Proceedings of the ASME 9th International Conference
on Offshore Mechanics and Arctic Engineering, Brest, France, 30 May–4 June 1999; Available online: [Link]
ISOPEIOPEC/proceedings-abstract/ISOPE99/All-ISOPE99/ISOPE-I-99-039/24645 (accessed on 7 January 2021).
J. Mar. Sci. Eng. 2021, 9, 87 19 of 19

13. Wu, G.X.; Eatock Taylor, R. A Green’s function form for ship motion at forward speed. Int. Shipbuild. Prog. 1987, 34, 189–196. [CrossRef]
14. Liapis, S.J. Time Domain Analysis of Ship Motions. Ph.D. Thesis, The University of Michigan, Ann Arbor, MI, USA, 1986.
15. King, B.K. Time Domain Analysis of Wave Exciting Forces on Ships and Bodies. Ph.D. Thesis, The University of Michigan,
Ann Arbor, MI, USA, 1987.
16. Rodrigues, J.M.; Guedes Soares, C. Froude-krylov forces from exact pressure integrations on adaptive panel meshes in a time
domain partially nonlinear model for ship motions. Ocean Eng. 2017, 139, 169–183. [CrossRef]
17. Zhang, L.; Li, Y.B.; Huang, D.B. The Effect of the Water Line Term on Wave Diffraction by a Floating Body with Forward Speed.
J. Harbin Eng. Univ. 1998, 19, 1–7.
18. Sun, W.; Ren, H.L. Ship motions with forward speed by time domain Green function method. Chin. J. Hydrodyn. 2018, 33, 216–222.
19. Singh, S.P.; Sen, D. A comparative linear and nonlinear ship motion study using 3-D time domain methods. Ocean Eng. 2007, 34,
1863–1881. [CrossRef]
20. Datta, R.; Rodrigues, J.M.; Soares, C.G. Study of the motions of fishing vessels by a time domain panel method. Ocean Eng. 2011,
38, 782–792. [CrossRef]
21. Lin, W.M. Numerical Solutions for Large-Amplitude Ship Motions in the Time Domain. In Proceedings of the 18th Symposium
on Naval Hydrodynamics, Ann Arbor, MI, USA, 19–24 August 1990.
22. Lin, W.M.; Zhang, S.; Weems, K.; Yue, D.K. A mixed source formulation for nonlinear ship motions and wave-induced loads.
In Proceedings of the 7th International Conference on Numerical Ship Hydrodynamics, Nantes, France, 19–22 July 1999.
23. Shan, P.; Wang, Y.; Wang, F.; Wu, J.; Zhu, R. An efficient algorithm with new residual functions for the transient free-surface green
function in infinite depth. Ocean Eng. 2019, 178, 435–441. [CrossRef]
24. Clement, A.H. An ordinary differential equation for the green function of time-domain free-surface hydrodynamics. J. Eng. Math.
1998, 33, 201–217. [CrossRef]
25. Shen, L.; Zhu, R.C.; Miao, G.P.; Liu, Y. A practical numerical method for deep water time domain in Green function. Chin. J.
Hydrodyn. 2007, 3, 380–386.
26. Zhong, W.X. On precise integration method. J. Comput. Appl. Math. 2004, 163, 59–78.
27. Li, Z.F.; Ren, H.L.; Tong, X.W.; Li, H. A precise computation method of transient free surface Green function. Ocean Eng. 2015,
105, 318–326. [CrossRef]
28. Magee, A.R.; Beck, R.F. Compendium of ship Motion Calculations Using Linear Time-Domain Analysis; (Report No. 310); Department
Naval Architects Marine Engineering, University of Michigan: Ann Arbor, MI, USA, 1988.
29. Journée, J.M.J. Experiments and Calculations on Four Wigley Hull Form; (Report 0909); Faculty of Mechanical Engineering and
Marine Technology, Delft University of Technology: Delft, The Netherlands, 1992.
30. Kara, F. Time Domain Hydrodynamic and Hydroelastic Analysis of Floating Bodies with Forward Speed. Ph.D. Thesis, University
of Strathclyde, Glasgow, UK, 2000.
31. Zhang, T.; Ren, J.S.; Zhang, X.F. Mathematical model of ship motion in regular waves based on three-dimensional time-domain
Green function method. J. Traffic Transp. Eng. 2019, 19, 110–121.
32. Hess, J.L.; Smith, A.M.O. Calculation of non-lifting potential flow about arbitrary three-dimensional bodies. J. Ship Res. 1964, 8,
22–44. [CrossRef]
33. Kukkanen, T. Numerical and Experimental Studies of Nonlinear Wave Loads of Ships. Ph.D. Thesis, Vtt Technical Research
Centre of Finland, Espoo, Finland, 2012.
34. Sun, L. Study of Ship-Generated Waves and Its Effects on Structures. Ph.D. Thesis, Dalian University of Technology, Dalian,
China, 2009.
35. Kim, K.H.; Kim, Y. Comparative study on ship hydrodynamics based on Neumann-Kelvin and double-body linearizations in
time-domain analysis. Int. J. Offshore Polar Eng. 2010, 10, 265–274.