Hydrodynamic Analysis and Motions of Ship with Forward Speed via a Three-Dimensional Time-Domain Panel Method

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.


Introduction
For the initial stages of ship design, accurate and reliable predictions of ship hydrodynamic analysis and motions in waves are essential. Various numerical methods are required to be developed for ship seakeeping analysis. For the wave-ship interaction problem with 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.
In the early researches, ship hydrodynamic analysis is developed based on twodimensional strip theories. Ogilvie and Tuck [2], Tasai [3], and Salvesen et al. [4] proposed a new strip theory, rational strip theory, and STF, method respectively. The STF strip 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 methods have been applied in the ship hydrodynamic and motion analysis for a long time. However, for the strip theories, the flow is assumed to be constrained in two-dimensional sections. Accurate hydrodynamic analysis can only be carried out for slender ships, and high frequency and low-speed assumptions are also required.
The shortcomings of strip theory can be overcome by 3D panel methods. Nakos [8] used the Rankine panel method to study the ship seakeeping problem in the frequency domain. Kring [9] and Chen [10] applied the Rankine panel method to ship hydrodynamic analysis in the time domain. The Rankine panel methods employ the Green function 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 and motions. Analytical integration expressions for F-K force formulations over the quadrilateral panels can be derived by using Green's theorem, which can avoid computational errors by numerical integration methods. To improve the accuracy and numerical stability of TFSGF evaluation, a numerical method is developed to solve the TFSGF by using a precise integration method with varying parameter settings. Based on the impulse response function method, the ship radiation problem and diffraction problem are solved by using TFSGF. The Wigley I hull is taken as a study case; convergence studies for ship hydrodynamic analysis and motions are conducted with respect to time step size and hull discretization. The computed hydrodynamic coefficients, wave exciting forces, and motions for the ship are compared to other solutions, such as Magee's method [28], previous literature experimental data [29], and so on. Thus, the three-dimensional time-domain panel method proposed in this paper is validated.

Coordinate Systems
For the present linear ship hydrodynamic analysis and motion problem, as shown in Figure 1, a freely floating ship is considered to advance with constant forward speed U in the presence of a linear incident wave field. The o-xyz is a reference, right-handed, Cartesian coordinate system with its origin o located amidship and travels along with the ship at the same speed U, the o-xy plane is coincident with the mean free surface z = 0, the positive x-axis is pointing upstream, the positive y-axis is pointing portside, and the positive z-axis is pointing vertically upwards. G-x b y b z b is a body-fixed, right-handed, Cartesian coordinate system, and the origin G is located at the gravity of the ship. At initial time t = 0, the space fixed coordinate system O-XYZ coincides with the reference coordinate system o-xyz, and Gz b axis is aligned with ox axis. The fluid domain Ω is enclosed by the ship hull surface S B , free surface S F , and surface S ∞ at infinity. n is the unit normal vector pointing inward the ship hull surface.
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.
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 and motions. Analytical integration expressions for F-K force formulations over the quadrilateral panels can be derived by using Green's theorem, which can avoid computational errors by numerical integration methods. To improve the accuracy and numerical stability of TFSGF evaluation, a numerical method is developed to solve the TFSGF by using a precise integration method with varying parameter settings. Based on the impulse response function method, the ship radiation problem and diffraction problem are solved by using TFSGF. The Wigley I hull is taken as a study case; convergence studies for ship hydrodynamic analysis and motions are conducted with respect to time step size and hull discretization. The computed hydrodynamic coefficients, wave exciting forces, and motions for the ship are compared to other solutions, such as Magee's method [28], previous literature experimental data [29], and so on. Thus, the three-dimensional time-domain panel method proposed in this paper is validated.

Coordinate Systems
For the present linear ship hydrodynamic analysis and motion problem, as shown in Figure 1, a freely floating ship is considered to advance with constant forward speed U in the presence of a linear incident wave field. The o xyz − is a reference, right-handed, Cartesian coordinate system with its origin o located amidship and travels along with the ship at the same speed U, the o-xy plane is coincident with the mean free surface z = 0, the positive x-axis is pointing upstream, the positive y-axis is pointing portside, and the positive z-axis is pointing vertically upwards.
is a body-fixed, right-handed, Cartesian coordinate system, and the origin G is located at the gravity of the ship. At initial time t = 0, the space fixed coordinate system O XYZ − coincides with the reference coor-

Boundary Value Problem
Based on the linear potential flow theory, the fluid is assumed to be irrotational, inviscid, and incompressible, and there is no fluid separation or lifting effect [30]. In the reference coordinate system o xyz − , the total velocity potential ( ) in the infinite water can be written as

Boundary Value Problem
Based on the linear potential flow theory, the fluid is assumed to be irrotational, inviscid, and incompressible, and there is no fluid separation or lifting effect [30]. In the reference coordinate system o-xyz, the total velocity potential Φ(r; t) in the infinite water can be written as 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. The perturbation potential Φ k (k = 1, 2, · · · , 7) satisfies the following governing equations and initial-boundary value conditions: where ζ k is unsteady motion in the kth mode; n = (n 1 , n 2 , n 3 ); r × n = (n 4 , n 5 , n 6 ); (m 1 , m 2 , m 3 ) = (0, 0, 0); (m 4 , m 5 , m 6 ) = (0, Un 3 , −Un 2 ); g is the gravity acceleration [31]. The TFSGF is adopted to solve the boundary value problem and written as [14] G(P, 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 G(P, Q; t − τ) of G(P, Q; t − τ) can be given as where r = |P − Q|; r = |P − Q |; Q (ξ, η, −ζ) 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; J 0 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 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 where Γ 0 is the mean waterline; t 0 is the initial time; n Q 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 S B in Equation (5) can be discretized by using the constant panel method [31].
The expression of radiation forces F jk can be written as where F jk is the radiation force in jth mode due to kth mode motion; The added mass A jk (ω) and damping coefficients B jk (ω) can be obtained via Fourier transform (j = 1, 2, · · · , 6; k = 1, 2, · · · , 6) [15]: 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] 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 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 F j7 (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 (12) whereK(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 F j7 (t).

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 where is hydrostatic forces in the ith mode, and M ij is the element of the general mass matrix.
In the body coordinate system G-x b y b z b , the six-degrees of motion equation can also be given as m where m is the mass matrix, I is inertial moment matrix, u = (u, v, w) is velocity 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 transformation matrix is the identity matrix.
The ship motion equations can be solved by various numerical methods, such as 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.

Analytical Expression for F-K Forces Evaluation
Consider two coordinate systems illustrated in Figure 2: the reference coordinate system o-xyz and the local panel coordinate system o -x y z defined by the vertices 1 to 4 in the counterclockwise direction, named P 1 , P 2 , P 3 , and P 4 , respectively. The positive o z axis points to the exterior of the fluid.

Analytical Expression for F-K Forces Evaluation
Consider two coordinate systems illustrated in Figure x y y y z y x z y z z z (17) where x , y , and z are unit base vectors in system o xyz The pressure IH p can be given by The transformation matrix between the position vector r in the reference coordinate system oxyz and the position vector r in the local panel coordinate system o -x y z is given as where which is the origin of the panel coordinate o -x y z system, T is the unit transformation cosine-director matrix between system o -x y z and system o-xyz and given by where x, y, and z are unit base vectors in system o-xyz; x , y , and z are unit 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 The pressure p IH can be given by 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 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 F Ii can be written as where n i is the unit normal vector pointing outward the fluid, and S i is the area of the quadrilateral panel.
The Equation (21) can be represented in the ith local panel coordinate system o -x y z The jth (j = 1, 2, 3, 4) edge for the ith quadrilateral panel can be parametrized by where (1), y 0,ij = y ij (0), and y 1,ij = y ij (1). A, B, and their derivatives on the jth side of the ith panel can be expressed with respect to parameter κ as With respect to parameter κ, A and B's derivatives on the ith panel are written as Let The Green's theorem with parameterization in Equation (23) is applied, and integration by parts is adopted to solve Equation (22). (22) is expressed as If Ψ y = 0, F Ii is directly solved by 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).

Numerical Method for TFSGF Evaluation
The memory part of TFSGF in Equation (4) can be written in the non-dimensional form where ν = υr , µ = −(z + ζ)/r , and β = g/r (t − τ). F(µ, β) is given by F(µ, β) is the solution to the following fourth-order ODE The fourth-order ODE Equation (33) can be given by .
; the relationship between the unit time-variant system and time-variant system in Equation (34) can be written as where the initial condition is x| s=0 = x(0).
The following equations can be obtained by using Equation (37) 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)
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

The Numerical Results of TFSGF
In Figure 3, the "analytical method" denotes the solutions obtained by the analytical expression for TFSGF when µ = 0 (presented in Equation (40)). Figure 3 shows the solutions obtained by the "PIM method", with constant m = 50 for 0 ≤ β ≤ 100 [27].

The Numerical Results of TFSGF
In Figure 3, the "analytical method" denotes the solutions obtained by the analytical expression for TFSGF when 0 μ = (presented in Equation (40)). Figure 3 shows the solutions obtained by the "PIM method", with constant m = 50 for 0 100 β ≤ ≤ [27].      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 "present method" only takes about 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 than the "PIM method".

The Study Case and Parameter Setting
The hull form of Wigley I hull can be defined as in the reference [29] where L is the length of the ship; B is the width of the ship; D is the draft of the ship. The main particulars of the Wigley I hull are given in Table 2; k yy denotes the pitch inertia radius of the ship, and ∇ denotes the displacement volume of the ship. The meshes of the ship's hull are illustrated in Figure 6.
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 μ = sent method" only takes about 85.8 s. The absolute errors of both "pres "PIM method" are within 9 10 − . Thus, the "present method" shows bett ciency than the "PIM method".

The Study Case and Parameter Setting
The hull form of Wigley I hull can be defined as in the reference [2 where L is the length of the ship; B is the width of the ship; D is th The main particulars of the Wigley I hull are given in Table 2; yy k inertia radius of the ship, and ∇ denotes the displacement volume of The meshes of the ship's hull are illustrated in Figure 6.  Non-dimensional forms of hydrodynamic coefficient, memory fun ing forces, and motion response are given below. Non-dimensional forms of hydrodynamic coefficient, memory function, wave exciting forces, and motion response are given below.
The non-dimensional frequency ω is defined as ω = ω L/g, the non-dimensional wave number υ is defined as υ = υL, and the non-dimensional frequency t is defined as t = t/T e (T e is encounter period).

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 test convergence. Wigley I hull is studied at Fn = 0.2 (Fn denotes Froude number, Fn = U/ gL) in head seas (α = π, η 0 = 0.036 m). The 2 × 2 Gaussian quadrature [15] can be used to compute the F-K forces and hydrodynamic forces acting on the panel. As shown in Figure 7, the square panel can be progressively subdivided until the prescribed precision is reached. e e

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 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 2 2 × Gaussian quadrature [15] can be used to compute the F-K forces and hydrodynamic forces acting on the panel. As shown in Figure 7, the square panel can be progressively subdivided until the prescribed precision is reached.  (Table 3).    (Table 3).   As shown in Figures 8 and 9, the numerical results tend to be convergent when 240 N ≥ , and the relative error between 240 N = and 320 N = is within 0.1%. Thus,    As shown in Figures 8 and 9, the numerical results tend to be convergent when 240 N ≥ , and the relative error between 240 N = and 320 N = is within 0.1%. Thus, panel number N = 320 is selected for ship hydrodynamic analysis and ship motion in the 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.
At t = 0, the F-K force acting on the square panel by the analytical integration method is 179.6193 N (the calculation accuracy can be set to 0.0001 N).
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 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 improve the F-K forces evaluation accuracy with no loss of accuracy.  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. 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 Gaussian quadrature, and "present method" denotes the results of F-K forces obtained by the analytical method.   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. 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 × Gaussian quadrature, and "present method" denotes the results of F-K forces obtained by the analytical method. From Figure 10, the F-K forces obtained by "numerical method" and "present method" are almost the same, and the relative error between the "numerical method" and "present method" is 0.1%. The CPU time consumption for the "numerical method" is about 35.4 s, while the CPU time consumption for the "present method" is about 27.8 s. The "present method" is more efficient than the "numerical method". The efficiency and accuracy of the "present method" can be validated. Figure 11 illustrates  From Figure 10, the F-K forces obtained by "numerical method" and "present method" are almost the same, and the relative error between the "numerical method" and "present method" is 0.1%. The CPU time consumption for the "numerical method" is about 35.4 s, while the CPU time consumption for the "present method" is about 27.8 s. The "present method" is more efficient than the "numerical method". The efficiency and accuracy of the "present method" can be validated. Figure 11 illustrates the time history of heave motion and pitch motion of Wigley I hull with four different time steps, which are T e /10, T e /20, T e /30, and T e /40, respectively. In Figure 11, as time step t Δ decreases, a convergent trend can be obtained, and the relative error between

Added Mass and Damping Coefficient
Since the hydrodynamic coefficients on the main diagonal have a significant influence on ship hydrodynamic analysis and motion, the non-dimensional coefficients 33 A′ , 55 A′ , 33 B′ , and 55 B′ 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 the comparisons of the non-dimensional pitch-pitch added masses and damping coefficients. In the legend of Figures 12 and 13, "present method" denotes numerical results obtained by the TFSGF method without waterline terms in Equation (5). "experiment" denotes previous literature experimental data [29] obtained by Journée, and the experiments were carried out in the Shiphydromechanics Laboratory of the Delft University of Technology. The experimental data on hydrodynamic coefficients, wave loads, and added resistance for heave and pitch motions in head waves of Wigley I hull are sufficient. "Magee's method" denotes numerical results obtained by the TFSGF method [28] in which the waterline integral terms of boundary integral equations are included. "LAMP method" denotes numerical results obtained by LAMP software [21], in which the body nonlinear method was adopted where the perturbation potential was computed on the instantaneous wetted hull under the mean free surface. In Figure 11, as time step ∆t decreases, a convergent trend can be obtained, and the relative error between ∆t = T e /30 and ∆t = T e /40 is within 0.1%. The time step ∆t = T e /40 is selected for the subsequent ship hydrodynamic analysis and ship motions.

Added Mass and Damping Coefficient
Since the hydrodynamic coefficients on the main diagonal have a significant influence on ship hydrodynamic analysis and motion, the non-dimensional coefficients A 33 , A 55 , B 33 , and B 55 are selected for detailed numerical analysis.
For the Wigley I hull at Fn = 0.2, Figure 12 presents the comparisons of the nondimensional heave-heave added masses and damping coefficients, and Figure 13 presents the comparisons of the non-dimensional pitch-pitch added masses and damping coefficients. In the legend of Figures 12 and 13, "present method" denotes numerical results obtained by the TFSGF method without waterline terms in Equation (5). "experiment" denotes previous literature experimental data [29] obtained by Journée, and the experiments were carried out in the Shiphydromechanics Laboratory of the Delft University of Technology. The experimental data on hydrodynamic coefficients, wave loads, and added resistance for heave and pitch motions in head waves of Wigley I hull are sufficient. "Magee's method" denotes numerical results obtained by the TFSGF method [28] in which the waterline integral terms of boundary integral equations are included. "LAMP method" denotes numerical results obtained by LAMP software [21], in which the body nonlinear method was adopted where the perturbation potential was computed on the instantaneous wetted hull under the mean free surface.   Tables 5-8 present absolute relative errors of hydrodynamic coefficients for Wigley I hull by "present method" and "Magee's method", which are compared with previous literature experimental data [29]. The absolute relative error can be written as ( ε denotes absolute relative error, num N denotes numerical results obtained by various models, and exp N denotes a value for previous literature experimental data [29]).  Tables 5-8 present absolute relative errors of hydrodynamic coefficients for Wigley I hull by "present method" and "Magee's method", which are compared with previous literature experimental data [29]. The absolute relative error can be written as ε = N num − N exp /N exp (ε denotes absolute relative error, N num denotes numerical results obtained by various models, and N exp denotes a value for previous literature experimental data [29]). 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 literature experimental data [29]. As the non-dimensional frequency increases, hydrodynamic coefficients obtained from "present method", "Magee's method", and "LAMP method" show a similar change trend. From Tables 5-8, the numerical results obtained by the "present method" are in better agreement with the previous literature experimental data [29] than those obtained by "Magee's method". Based on the potential flow assumption, the viscosity of a fluid is ignored. When both the field point and source point are on the mean free surface, the amplitude and oscillating frequency of TFSGF increase rapidly with time parameter β, and there will be a larger numerical error if waterline integral terms of boundary integral equations are included. The accuracy for hydrodynamic coefficient evaluation can be improved by the present method with waterline integral terms excluded.

The Wave Exciting Forces
The Wigley I hull advances in head regular waves at Fn = 0.2. The wave exciting forces are composed of F-K forces and diffraction forces in the frequency domain, and diffraction forces in the frequency domain can be obtained by Fourier transformation [15]. Figure 14 shows the numerical results of amplitudes of the non-dimensional wave exciting forces obtained by various methods. "SAMP method" denotes numerical results obtained by a transient free surface Green function method in the linear time domain [21]. , and there is a larger deviation between numerical results and previous literature experimental data [29]. As the non-dimensional frequency increases, hydrodynamic coefficients obtained from "present method", "Magee's method", and "LAMP method" show a similar change trend. From Tables 5-8, the numerical results obtained by the "present method" are in better agreement with the previous literature experimental data [29] than those obtained by "Magee's method". Based on the potential flow assumption, the viscosity of a fluid is ignored. When both the field point and source point are on the mean free surface, the amplitude and oscillating frequency of TFSGF increase rapidly with time parameter β , and there will be a larger numerical error if waterline integral terms of boundary integral equations are included. The accuracy for hydrodynamic coefficient evaluation can be improved by the present method with waterline integral terms excluded.

The Wave Exciting Forces
The Wigley I hull advances in head regular waves at Fn = 0.2. The wave exciting forces are composed of F-K forces and diffraction forces in the frequency domain, and diffraction forces in the frequency domain can be obtained by Fourier transformation [15]. Figure 14 shows the numerical results of amplitudes of the non-dimensional wave exciting forces obtained by various methods. "SAMP method" denotes numerical results obtained by a transient free surface Green function method in the linear time domain [21].  From Figure 14, the maximum relative error between "present method" and "experiment" is within 10.0%, and the non-dimensional wave exciting forces obtained from "present method", "Magee's method" and "SAMP method" show a similar change trend. Thus, the reliability of analytical expression for F-K forces evaluation is verified. From Figure 14, the maximum relative error between "present method" and "experiment" is within 10.0%, and the non-dimensional wave exciting forces obtained from "present method", "Magee's method" and "SAMP method" show a similar change trend. Thus, the reliability of analytical expression for F-K forces evaluation is verified. 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- 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 amplitude operators for ship motions in head regular waves at Fn = 0.2. The response amplitude of heave motion can be defined as ζ 30 /η 0 , and ζ 30 is the amplitude of ship heave motions. The response amplitude of pitch motion can be defined as (ζ 50 L)/(2πη 0 ), 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 time domain [35]. 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 amplitude operators for ship motions in head regular waves at Fn = 0.2. The response amplitude of heave motion can be defined as 30 0 ζ η , and 30 ζ is the amplitude of ship heave motions. The response amplitude of pitch motion can be defined as ( ) ( ) 50 0 2π L ζ η , 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 time domain [35].    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 frequency 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".

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 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 quadrilateral 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 discretize 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 instantaneous 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.; writingoriginal 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.

Data Availability Statement:
The data used to support the findings of this study are available from the corresponding author upon request.