Numerical Study of Wave Drift Load and Turning Characteristics of KVLCC2 Ship in Regular Waves Based on TEBEM

: Maritime trafﬁc has increased considerably in recent years, making energy efﬁciency and navigation safety of ships more crucial than ever. Hence, a two-time scale model based on the Taylor expansion boundary element method (TEBEM) is proposed to predict ship turning trajectories in regular waves. The maneuvering motion is calculated using a three degrees of freedom MMG model that considers the wave drift loads. TEBEM overcomes the shortcomings of the constant panel method in solving tangential induced velocity at a non-smooth boundary and that of the high-order boundary element method in dealing with a high-order derivative of the velocity potential at the corner. This signiﬁcantly improves the calculation accuracy of the induced velocity and high-order derivative of velocity potential. Firstly, based on the TEBEM, the surge and sway wave drift forces and yaw moment of the KVLCC2 model with drift angle under full wave headings are calculated and compared with computational ﬂuid dynamics results, using which the calculation accuracy of TEBEM is veriﬁed. Subsequently, the two-time scale model is used to calculate the turning trajectories of the KVLCC2 model in regular waves with different wave headings, wave frequencies, and wave steepness. The numerical results show that the drift angle has a certain effect on the wave drift loads of the ship, and the proposed model can effectively predict the ship’s turning motion in regular waves.


Introduction
With the recent increase in maritime traffic, energy efficiency, and navigation safety of ships have become crucial, and ship coupled hydrodynamic performance analysis has become a research focus. In the early stages of the development of ship hydrodynamics, limited by theoretical models and computational power, some assumptions and simplifications must be made to conduct studies on ship hydrodynamics. Presently, the study of ship seakeeping mainly focuses on the wave drift loads at zero drift angle. However, the ship usually has a certain drift angle when moving under actual sea conditions, affecting the calculation of wave drift loads. Ship maneuverability prediction is typically conducted under calm water conditions, but during actual sea navigation, the seakeeping and the maneuvering motions occur simultaneously and affect each other. The International Towing Tank Conference (ITTC) set up a special committee on ship maneuvering in waves to review the advances in analyzing the coupled motion of ship maneuvering and seakeeping. At present, the popular coupling analysis methods for ship maneuvering and seakeeping can be divided into the following: experimental method, numerical simulation using computational fluid dynamics (CFD), numerical simulation based on the unified theory, and numerical simulation based on two-time-scale models.
The physical tank model test is the most reliable method to study the hydrodynamic performance of ships. The most representative experimental research is the work of Yasukawa's group at Hiroshima University [1][2][3][4], which investigated ship maneuverability in waves in a physical tank on SR108 and S175 container ship models and a KVLCC2 ship model. Xu et al. [5] and Kinoshita et al. [6] conducted planar motion mechanism (PMM) tests in waves to measure the hydrodynamic force and horizontal displacement of the ship, thereby analyzing the influence of the low-frequency wave drift forces on the ship maneuvering. Lee et al. [7] used a KVLCC model to determine the trajectory for the turning and zigzag maneuvers in regular waves with different wave heights and wavelengths. The test results show that second-order wave forces significantly impact ship maneuverability. Kim et al. [8] also conducted a tank-based study on ship maneuvering of a KVLCC2 ship in regular waves and analyzed the influence of wave heading, wave frequency, and wave steepness on the turning motion. These experiments also provide comparative analysis and validation criteria for numerical models of ship maneuvering prediction in waves.
With the recent increase in computational power, the CFD method has been adopted by the industry for the modeling of ship hydrodynamic performance. High accuracy can be achieved, with results being within 3% of the test results from a physical tank, especially in calm water. However, using the CFD method to simulate ship maneuvering directly in waves remains time-consuming. Kim et al. [9] used the Reynolds-averaged Navier-Stokes (RANS) solver to numerically simulate PMM maneuvering motion to obtain hydrodynamic derivatives of ship maneuvering. Shang et al. [10] used the overlapping grid technology and a propeller volumetric force model to realize numerical prediction of zigzag maneuvers in calm water conditions. Ma et al. [11] conducted numerical and physical tank studies on ship PMM movement in regular waves using the URANS (Unsteady RANS) viscous flow solution and overlapping grid technology.
Numerical simulation based on the two-time scale model and unified theory is more widely used than the CFD method. Both methods are based on the potential flow theory but differ in their treatment of the details. The method based on the unified theory adopts the 6 degrees of freedom (6-DOF) motion equation of the ship in calm water; the F-K, radiation, and diffraction forces caused by the waves are included in the equation as external forces, resulting in a set of unified 6-DOF rigid body motion equations to describe the ship's motion in the waves. Hamamoto and Kim [12] proposed the horizontal body-fixed coordinate system, thus simplifying the 6-DOF motion description of ships in waves. Sutulo and Soares [13] proposed an auxiliary state variable method to simplify the calculation of the convolution integral in the "memory effect." Fang et al. [14] calculated the additional mass and damping of the ship according to the instantaneous encounter frequency in the ship's turning process and proposed a nonlinear model to calculate the ship's turning trajectory in waves. Lin et al. [15] expanded LAMP, a three-dimensional seakeeping software, and conducted numerical prediction of the maneuverability of various ship models in waves. Subramanian and Beck [16] proposed a body-exact strip theory to simulate the maneuvering of a ship in a seaway, in which the wave drift load effects were considered via the square of the velocity potential gradient. Based on the model proposed by Hamamoto and Kim [12], Suzukie et al. [17] formed their own model by improving the external wave forces and the inertial and damping components of the hydrodynamic forces induced and simulated the turning motion of the KVLCC1 and KCS models in regular waves.
Although numerical simulation methods based on unified theory can consider the effect of wave drift forces, the calculation usually requires some approximation. Therefore, their calculation accuracy is often difficult to guarantee. In comparison, the two-time scale method has some advantages in calculating the wave drift forces. The two-time scale method divides the total ship motion into ship maneuvering and wave-induced motions, which are described by two different sets of motion equations. Seo and Kim [18], Zhang and Zou [19], and Xie et al. [20] built a prediction model of the ship's turning and maneuvering motion in waves based on the two-time scale method and three-dimensional potential flow theory. Yao et al. [21] used the two-time scale method to construct a prediction model of ship maneuvering in waves, but the difference is that a time-average model is proposed for the prediction model of ship seakeeping. Mei et al. [22] also built a coupled prediction model of ship motion using the two-time scale model. They used the RANS viscous flow CFD method to calculate the hydrodynamic derivative of the ship and the three-dimensional potential flow theory to predict wave drift loads. Jeon M. et al. [23] used EBM method to construct a fast prediction model of ship wave drift load, which was coupled with the 3-DOF MMG motion equation to form a ship handling motion prediction model in waves. Yasukawa et al. [24] calculated the wave-induced steady forces by using the zero-speed three-dimensional panel method and the wave-induced force formula based on the Cochin function of a hypothetical slender ship. Zhang et al. [25] calculated the 6-DOF motions and second-order wave drift force and moment of the KVLCC2 Tanker and S-175 Container ship in waves by using the near-field integral formula.
Presently, the seakeeping of ships is mainly calculated by the constant surface element method and the high-order boundary element method (HOBEM) based on the threedimensional potential flow theory. The constant surface element method is widely used because of its straightforward programming. However, it yields low accuracy in solving velocity potential and tangential induced velocity at the non-smooth boundary. Compared with the constant boundary element method, HOBEM improves the calculation accuracy under the same number of panel elements. Nevertheless, the calculation accuracy is often lower when the boundary has sharp angles and the higher-order derivative of the velocity potential is needed. Thus, in this study, the Taylor expansion boundary element method (TEBEM) is proposed to solve the velocity field and higher-order derivative field of the boundary value problem of the arbitrary flow field, thereby significantly improving the numerical accuracy of wave drift loads. Moreover, the two-time scale model based on TEBEM exhibits excellent numerical accuracy and computational efficiency in calculating the turning motion of ships in waves.
In this study, based on the Taylor expansion boundary element method (TEBEM) (Duan et al.) [26], the 6-DOF motions and wave drift loads of the KVLCC2 model with drift angle are calculated in oblique waves and compared with the CFD calculation results. The comparison results show that the two methods are in good agreement, thus verifying the accuracy of TEBEM for calculating wave drift loads. Then, based on the two-time scale model, the turning trajectory of a KVLCC2 ship is calculated in regular waves with different wave headings, wave frequencies, and wave steepness. The numerical simulation results agree with the experimental results reported by Kim et al. [8].

Two-Time Scale Model
In this study, the hydrodynamic problems of ships are addressed via potential flow theory. The total ship motion is divided into low-frequency ship maneuvering and highfrequency wave-induced motions, and the two are coupled with each other. The lowfrequency motion is affected by the second-order wave drift loads, which indirectly depend on the ship's wave-induced motion. Conversely, the high-frequency motion changes with the ship's speed and heading angle, which must be determined by solving the lowfrequency motion. Therefore, this study constructs a mathematical model of the ship's maneuvering motion based on the two-time scale model. In the numerical calculation, the following steps are adopted to deal with the coupling solution of the two sets of equations: (1) The time domain TEBEM is used to solve the boundary value problem of the ship at its initial velocity and position. Subsequently, the first-order hydrodynamic force is obtained by the pressure integration on the wet surface of the ship, and the highfrequency motion of the ship is obtained using the wave-induced 6-DOF motion equations. Finally, the wave drift forces and moment are directly calculated using the near-field integral method; (2) The boundary value problem of the ship is updated via the high-frequency motion obtained in step (1). Subsequently, the time domain TEBEM is used to solve the new boundary value problem of the ship at its initial velocity and position. Finally, the new high-frequency motion of the ship and wave drift forces and moment are obtained. Steps (1) and (2) are termed the seakeeping problem, and Nt steps are calculated circularly; (3) The obtained wave drift forces and moment are fitted via the least-squares method, and the result is incorporated into the MMG maneuvering motion. The ship velocity, heading angle, and position of the ship in the space-fixed coordinate system at the next moment are obtained by solving the MMG maneuvering motion equation; (4) The double-body (DB) flow velocity potential is recalculated according to the ship velocity, and the boundary value conditions satisfied by the disturbed potential are updated according to the calculated results; (5) The time-domain TEBEM is used to solve the boundary value problem of the updated disturbed potential. Subsequently, the wave frequency motion and the wave drift forces and moment are recalculated, and the N t seakeeping steps are calculated circularly; (6) The wave drift forces and moment obtained in step (5) are fitted via the least-squares method. They are then brought back to the maneuvering motion equation for the next step of the maneuvering numerical simulation. (7) The above six steps are repeated until the calculation time is met. Figure 1 shows the algorithm process of the two-time scale model.
near-field integral method; (2) The boundary value problem of the ship is updated via the high-frequency motion obtained in step (1). Subsequently, the time domain TEBEM is used to solve the new boundary value problem of the ship at its initial velocity and position. Finally, the new high-frequency motion of the ship and wave drift forces and moment are obtained. Steps (1) and (2) are termed the seakeeping problem, and Nt steps are calculated circularly; (3) The obtained wave drift forces and moment are fitted via the least-squares method, and the result is incorporated into the MMG maneuvering motion. The ship velocity, heading angle, and position of the ship in the space-fixed coordinate system at the next moment are obtained by solving the MMG maneuvering motion equation; (4) The double-body (DB) flow velocity potential is recalculated according to the ship velocity, and the boundary value conditions satisfied by the disturbed potential are updated according to the calculated results; (5) The time-domain TEBEM is used to solve the boundary value problem of the updated disturbed potential. Subsequently, the wave frequency motion and the wave drift forces and moment are recalculated, and the Nt seakeeping steps are calculated circularly; (6) The wave drift forces and moment obtained in step (5) are fitted via the least-squares method. They are then brought back to the maneuvering motion equation for the next step of the maneuvering numerical simulation. (7) The above six steps are repeated until the calculation time is met. Figure 1 shows the algorithm process of the two-time scale model.

Coordinate Systems
Two right-handed coordinate systems are considered, as shown in Figure 2. The first coordinate system, O-XYZ is a fixed coordinate system in space, with positive Z pointing downward. The second coordinate system, o-xyz, moves with the ship's maneuvering speed but does not oscillate with the ship. The oxy plane overlaps with the calm water surface, and its origin o is located at the intersection of the middle transverse section, middle longitudinal section, and the calm water surface. The ox, oy, and oz axes pointed to the bow, starboard, and vertically downward, respectively.

Coordinate Systems
Two right-handed coordinate systems are considered, as shown in Figure 2. The first coordinate system, O-XYZ is a fixed coordinate system in space, with positive Z pointing downward. The second coordinate system, o-xyz, moves with the ship's maneuvering speed but does not oscillate with the ship. The oxy plane overlaps with the calm water surface, and its origin o is located at the intersection of the middle transverse section, middle longitudinal section, and the calm water surface. The ox, oy, and oz axes pointed to the bow, starboard, and vertically downward, respectively. The space fixed coordinate system coincides with the ship moving coordinate system at the initial moment and have the following relationships: The space fixed coordinate system coincides with the ship moving coordinate system at the initial moment and have the following relationships: where Ψ is the heading angle, which is the included angle between the OX axis and the ox axis' positive direction, and (X 0 , Y 0 ) is the coordinate of the origin of the ship's moving coordinate system in the fixed space coordinate system, which is used to represent the ship's spatial position at different times. The incident wave angle χ is the angle between the wave heading and the X axis of the fixed coordinate system. When χ = 0 • , the incident wave propagates along the positive direction of the X axis and is positive clockwise. δ represents the rudder angle, which stipulates that the starboard is positive, and t denotes time.

Basic Equations of Ship Maneuvering Motion
The ship's motion equation is established in the ship's moving coordinate system. In the calculation of ship motion, it is assumed that the ship is a moving rigid body; its shape and internal mass distribution do not change with the motion. Its longitudinal velocity is u, transverse velocity is v, combined velocity is U = √ u 2 + v 2 , and yaw angular velocity is r. The basic equation of ship maneuvering motion is governed via the following equations: where m represents the ship's mass, I zz represents the moment of inertia around the Zaxis, x G is the longitudinal coordinate of the center of gravity, and F x , F y, and M z are the hydrodynamic forces and moment suffered by the ship, respectively. In maneuvering, the hydrodynamic force acting on the hull comprises inertial and non-inertial forces. In this study, the hydrodynamic forces and moment can be expressed by: where m x and m y are the corresponding added masses, and J zz is the added moment of inertia. X P , X R , and X H represent the longitudinal low-frequency hydrodynamic forces on the propeller, rudder, and hull, respectively. X W denotes the surge wave drift force. The lateral force caused by propeller rotation and its torque on the center of gravity are ignored. Combining Equations (2) and (3), a new equation of ship motion can be obtained:

Hydrodynamic Forces Acting on Ship Hull
The hydrodynamic forces acting on the ship hull are expressed as follows: where L PP is the ship length between perpendiculars. v and r denotes the non-dimensional lateral velocity and yaw rate defined via v/U and rL PP /U, respectively. X H , Y H , and N H are expressed as follows: where R' 0 represents the resistance coefficient, and X vv , X vr , X rr , X vvvv , Y v , Y r , N v , N r , etc., are termed the non-dimensional hydrodynamic derivatives. Compared with the other hydrodynamic derivatives, Y H and N H calculated via the 1st and 3rd order polynomial functions of v and r have better estimation accuracy.

Hydrodynamic Force Due to Propeller
The hydrodynamic force due to the propeller X p is expressed as: where t p is the thrust deduction fraction and is assumed to be constant. T is the propeller thrust and is expressed as follows: where ρ, K T , D p , and n p are the water density, thrust coefficient, propeller diameter, and propeller revolution, respectively. K T is expressed via the propeller advance ratio J P : J P is expressed as: where w p is the wake coefficient at the propeller, which is usually expressed via the wake coefficient w p0 in a straight motion. The Yasukawa et al. [27] model is adopted as follows: ( where C 1 and C 2 are the experimental constants representing the wave characteristics during maneuvering obtained via Yasukawa et al. [27]. β p is the geometrical inflow angle to the propeller in maneuvering motion and is calculated as follows: where x p is the non-dimensional longitudinal coordinate of the propeller by L PP and β is the drift angle.

Hydrodynamic Forces by Steering
The rudder forces and moment are expressed as: where F N , δ, and x R are the rudder normal force, rudder angle, and longitudinal position of the rudder, respectively. α H , t R , and x H are the rudder force increase factor, steering resistance deduction factor, and the position of an additional lateral force component, respectively.
Rudder normal force F N is expressed as: where A R and f α represent the profile area of the movable part of the rudder and rudder lift gradient coefficient, respectively. U R and α R represent the resultant rudder inflow velocity and the inflow angle, respectively, and are expressed as: Here, δ is the rudder angle. u R and v R represent the longitudinal and lateral inflow velocity components to the rudder, respectively.
Here, ε, l' R , κ, and γ R represent the ratio of wake fraction at the rudder position to that at the propeller position, the effective non-dimensional longitudinal coordinate of the rudder position, an experimental constant, and the flow straightening coefficient, respectively. η = D p /H R , where H R is the rudder span length. The relevant coefficient values are obtained via Yasukawa et al. [27].

Boundary Value Problem and Wave Drift Loads
Based on the potential flow theory, incompressible and inviscid water is assumed, and the flow is irrotational. Thus, the water motion can be described by the velocity potential Θ, which satisfies the Laplace equation and boundary value conditions. The boundary value problem can be linearized by decomposing the velocity potential Θ, including the base velocity potential Φ b , the steady wave flow velocity potential, Φ, and the unsteady wave flow velocity potential ϕ.
The base velocity potential Φ b can be divided into the incoming flow velocity potential and the DB flow velocity potential, Φ. The unsteady wave flow potential ϕ can be divided into the disturbed potential ϕ d and incident potential ϕ I . Wave elevation can also be divided into steady, disturbed, and incident terms (ζ S ,ζ d , and ζ I ), respectively.
Only the unsteady flow is discussed in this paper; the steady wave flow velocity potential and the steady wave elevation are ignored. Based on the DB assumption, the velocity potential Φ satisfies the following boundary value problem: where S H is the mean wet surface of the hull, and S F is the undisturbed free surface.
After the velocity potential Φ is solved, the m j terms can be expressed as follows: r denotes the distance between the fluid point and mass centroid of the ship.
The linearized incident velocity potential in deep water is known analytically, as follows: where A, k, g, and ω denote the wave amplitude, wave number, acceleration of gravity, and frequency of the incident wave, respectively. Ψ and χ are the heading and incident angles of the wave. The damping zone is used to deal with the distant radiation condition of the ship with a forward speed problem. Thus, the boundary value problem of the disturbed potential is: Here (n 1 , n 2 , n 3 ) = n j , (n 4 , n 5 , n 6 ) = → r × n j , and η j (j = 1 − 6) are the 6-DOF motion displacements of the ship. µ = 3µ 0 (l − l 0 ) 2 /L 2 − 2µ 0 (l − l 0 ) 3 /L 3 , f 1 = ∂ 2 Φ ∂z 2 − ∇Φ · ∇ ζ I and f 2 = −∇Φ · ∇ϕ I , where l 0 , L, and µ 0 are the initial edge, damping length, and the damping intensity of the damping zone, respectively.
When the unsteady disturbed potential is solved, the first-order wave force can be obtained via integrating the pressure on the wetted surface of the ship, as follows: We can also solve the restoring force via the restoring force coefficient matrix. Then, the motion of the ship can be obtained from the wave-induced 6-DOF motion equations, as follows: Here (ζ 1 , ζ 2 , ζ 3 ) and (ζ 4 , ζ 5 , ζ 6 ) represent the translation and rotational displacements of the ship, respectively. m ij represents the elements of the inertial matrix for the hull. F i denotes the first-order wave forces, including restoring, radiation, diffraction, and Froude-Kriloff forces. The first-order wave forces are also referred to as the wave frequency hydrodynamic forces.
The wave drift loads can be obtained via pressure integration on the body surface. We adopted the formulation of the pressure integration provided by Chen et al. [28] as follows: )ds (29) Here,

Taylor Expansion Boundary Element Method (TEBEM)
Given the low calculation accuracy of the constant surface element method for solving the sharp corner of the object surface and since the high-order boundary element method cannot directly calculate the high-order derivative of the velocity potential at the corner, the TEBEM is used to solve the disturbance velocity potential and its higher-order derivatives, which can effectively improve the calculation accuracy. The TEBEM is based on the source and dipole mixed distribution model. The potential functions ϕ and Green's function 1/r (p, q) are the basic solutions of the Laplace equation. r represents the distance between the field point p and the source point q. When the field point p is in the fluid domain, the following formula can be obtained using Green's third formula for ϕ and 1/r (p, q): Here, the unit normal vector n q points to the outer part of the fluid domain; S represents the boundary of the fluid domain and can be discretized into N small triangular or quadrilateral elements. The average value of each element node is taken as the center and the discrete format is as follows: When the field point P approaches the boundary of the fluid domain along the negative direction of the normal vector, the boundary integral can be obtained as follows: Here N represents the number of panels at the surface of the ship and the free surface.
For the dipole strength in Equation (32), the Taylor expansion is performed in the local coordinate system of each panel, and the first-order derivative term is retained as follows: where η and ξ are the local coordinates of the source point, as shown in Figure 3. For the dipole strength in Equation (32), the Taylor expansion is performed in the local coordinate system of each panel, and the first-order derivative term is retained as follows: where η and ξ are the local coordinates of the source point, as shown in Figure 3. The new boundary integral can be obtained via substituting Equation (33) into Equation (32) as: To solve Equation (34), we need to add a new boundary integral equation. The firstorder derivatives of the potential at the field point are considered, along with the two mutually orthogonal tangential directions. The supplementary equations are: The new boundary integral can be obtained via substituting Equation (33) into Equation (32) as: To solve Equation (34), we need to add a new boundary integral equation. The firstorder derivatives of the potential at the field point are considered, along with the two mutually orthogonal tangential directions. The supplementary equations are: 2π ∂ϕ(p,t) ∂y Finally, we can obtain the dipole strength and its first-order derivatives via solving Equations (34)-(36). All the influencing coefficients for the solution process can be found in Duan et al. [26].

Convergence Analysis
The KVLCC2 ship model is used to calculate the wave drift loads of ships with drift angles in oblique waves and the turning motion in regular waves. The main parameters of the KVLCC2 ship model are shown in Table 1. The convergence analysis is conducted to verify the reliability of the numerical results. Based on the TEBEM, the wave drift loads of the KVLCC2 model are calculated with wave height H = 0.02 L PP , wavelength λ = 1.0 L PP , wave heading χ = 120 • , drift angle β = 0 • , and U = 6.0 kn. Figure 4 is a grid diagram of the hull surface below the waterline. The bow and stern of the ship are triangular grids, and the rest are quadrilateral grids. Figure 5 is the grid diagram of the free surface. Since the damping region is around the computational domain, the grid becomes sparser the farther away from the ship's surface. Four soft springs are added in the horizontal direction of the ship to prevent the divergence of numerical simulation results, as shown in Figure 6. The spring's stiffness is chosen so that the natural period of planar motion is five times the highest wave period.
The bow and stern of the ship are triangular grids, and the rest are quadrilateral grids. Figure 5 is the grid diagram of the free surface. Since the damping region is around the computational domain, the grid becomes sparser the farther away from the ship's surface. Four soft springs are added in the horizontal direction of the ship to prevent the divergence of numerical simulation results, as shown in Figure 6. The spring's stiffness is chosen so that the natural period of planar motion is five times the highest wave period.    Figure 5 is the grid diagram of the free surface. Since the damping region is around the computational domain, the grid becomes sparser the farther away from the ship's surface. Four soft springs are added in the horizontal direction of the ship to prevent the divergence of numerical simulation results, as shown in Figure 6. The spring's stiffness is chosen so that the natural period of planar motion is five times the highest wave period.    Figure 9 shows the yaw wave drift moment of the KVLCC2 ship under H = 0.02 LPP, λ = 1.0 LPP, χ = 120°, β = 0°, and U = 6.0 kn. As shown in Figure 9, the yaw wave drift moment of the KVLCC2 ship converges well.   Figure 9 shows the yaw wave drift moment of the KVLCC2 ship under H = 0.02 L PP , λ = 1.0 L PP , χ = 120 • , β = 0 • , and U = 6.0 kn. As shown in Figure 9, the yaw wave drift moment of the KVLCC2 ship converges well. of the free surface. Finally, the convergence of the time-step scales (dt = T/45, dt and dt = T/65) is analyzed, where the panel numbers of the full-body surface and th numbers of the free surface are NB = 2966 and NF = 4584, respectively. Figure 9 sh yaw wave drift moment of the KVLCC2 ship under H = 0.02 LPP, λ = 1.0 LPP, χ = 1 0°, and U = 6.0 kn. As shown in Figure 9, the yaw wave drift moment of the KVLC converges well.  and dt = T/65) is analyzed, where the panel numbers of the full-body surface and th numbers of the free surface are NB = 2966 and NF = 4584, respectively. Figure 9 sh yaw wave drift moment of the KVLCC2 ship under H = 0.02 LPP, λ = 1.0 LPP, χ = 1 0°, and U = 6.0 kn. As shown in Figure 9, the yaw wave drift moment of the KVLC converges well.   Therefore, comprehensively considering the overall calculation accura efficiency, the panel numbers of the full-body and free surfaces are selected as 2 4584, respectively, and the time-step scale is T/55 for the subsequent calculation.
To verify the accuracy of the TEBEM, the Computational Fluid Dynamic method is also used to calculate the wave drift loads of the ship by the popular com program STAR-CCM+ (v15.04). The second-order wave drift load of the ship in wa be obtained by subtracting the surge, sway, and yaw forces and moment of the calm water from the surge, sway, and yaw forces and moment of the ship in Similarly, the convergence of the computational domain mesh, time-step scal damping zone length is analyzed for the surge, sway, and yaw forces and momen Therefore, comprehensively considering the overall calculation accuracy and efficiency, the panel numbers of the full-body and free surfaces are selected as 2966 and 4584, respectively, and the time-step scale is T/55 for the subsequent calculation.
To verify the accuracy of the TEBEM, the Computational Fluid Dynamic (CFD) method is also used to calculate the wave drift loads of the ship by the popular commercial program STAR-CCM+ (v15.04). The second-order wave drift load of the ship in waves can be obtained by subtracting the surge, sway, and yaw forces and moment of the ship in calm water from the surge, sway, and yaw forces and moment of the ship in waves. Similarly, the convergence of the computational domain mesh, time-step scales, and damping zone length is analyzed for the surge, sway, and yaw forces and moment of the ship in waves. Figures 10 and 11 are grid diagrams of the longitudinal section and damping zone of the computational domain, respectively. We retain only the first-order and second-order results of the calculated wave drift loads. Therefore, Equation (37) is used to carry out the least square fitting, where ω represents wave frequency. Figure 12 shows the least square fitting of the wave drift loads of the KVLCC2 ship under H = 0.02 L PP , λ = 1.0 L PP , χ = 120 • , β = 0 • , and U = 6.0 kn.
efficiency, the panel numbers of the full-body and free surfaces are selected as 2966 and 4584, respectively, and the time-step scale is T/55 for the subsequent calculation.
To verify the accuracy of the TEBEM, the Computational Fluid Dynamic (CFD) method is also used to calculate the wave drift loads of the ship by the popular commercial program STAR-CCM+ (v15.04). The second-order wave drift load of the ship in waves can be obtained by subtracting the surge, sway, and yaw forces and moment of the ship in calm water from the surge, sway, and yaw forces and moment of the ship in waves. Similarly, the convergence of the computational domain mesh, time-step scales, and damping zone length is analyzed for the surge, sway, and yaw forces and moment of the ship in waves. Figures 10 and 11 are grid diagrams of the longitudinal section and damping zone of the computational domain, respectively. We retain only the first-order and second-order results of the calculated wave drift loads. Therefore, Equation (37) is used to carry out the least square fitting, where ω represents wave frequency. Figure 12 shows the least square fitting of the wave drift loads of the KVLCC2 ship under H = 0.02 LPP, λ = 1.0 LPP, χ = 120°, β = 0°, and U = 6.0 kn.

Wave Drift Loads Analysis of Ships with Zero Drift Angle in Oblique Waves
Since the wave drift loads play an important role in the ship's turning motion, the wave drift loads and the 6-DOF motion of the KVLCC2 ship are calculated first via TEBEM Figure 13. The convergence analysis surge, sway, and yaw forces and moment of the KVLCC2 ship under different mesh numbers, time-step scales, and damping zone length, respectively.

Wave Drift Loads Analysis of Ships with Zero Drift Angle in Oblique Waves
Since the wave drift loads play an important role in the ship's turning motion, the wave drift loads and the 6-DOF motion of the KVLCC2 ship are calculated first via TEBEM and CFD method. Figure 14 shows the KVLCC2 ship's RAO of 6-DOF motion under H = 0.02 L PP , χ = 120 • , β = 0 • , and U = 6.0 kn. As shown in Figure 14, the computational results of CFD and TEBEM agree well with the experimental results of Seo et al. [29]. When λ/L PP = 1.7, the roll motion significantly increases, indicating that the ship is near the roll resonance. However, the amplitude of roll resonance is obviously smaller than the experimental value, because only linearized viscous roll damping is considered at approximately 2% of the critical roll damping, which is not accurate enough for cases with non-zero speed.  Figure 15 shows the surge, sway, and yaw wave drift forces and moment of the KVLCC2 ship under H = 0.02 L PP , χ = 120 • , β = 0 • , and U = 6.0 kn. The results of the TEBEM and CFD methods are also consistent with the trend of experimental values of Seo et al. [29]. The results calculated via TEBEM vary noticeably when λ/L PP ≥ 1.5, due to the roll resonance motion. Vertical displacement δ 3 is involved in Equations (28) and (29), which is the function of the roll motion displacement η 4 , δ 3 = η 3 + η 4 y − η 5 x. For surge, sway, and yaw wave drift forces and moment, the accuracy of TEBEM can meet the requirements compared with the experimental values of Seo et al. [29].

Wave Drift Loads Analysis of Ships with Different Drift Angles in Oblique Waves
Presently, the wave drift loads are mainly studied with zero drift angle. However, the existence of a drift angle will have a certain influence on the wave drift loads in the turning motion of the ship. Therefore, it is critical to study the wave drift loads of the ship with a drift angle in oblique waves. As described in this section, the RAO of 6-DOF motion and wave drift loads of the KVLCC2 ship are first calculated via TEBEM with H = 0.02 LPP, χ = 120°, β = 0°, 5.7°, 11.3°, 18.4°, and U = 6.0 kn. Then, the RAO of 6-DOF motion and wave drift loads of the KVLCC2 ship are calculated via the CFD method with H = 0.02 LPP, χ = 120°, β = 5.7°, 18.4°, and U = 6.0 kn, as shown in Figures 16 and 17. The comparison between the TEBEM and CFD methods shows that the 6-DOF motion, surge, and sway wave drift forces of the ship with β = 5.7° and 18.4° are in good agreement, whereas the yaw drift moments have large discrepancies in the short waves. This is because significant viscous effects due to flow separation occurred in the bow and stern areas, resulting in a small yaw moment calculated by the CFD method. Compared with the CFD results, the accuracy of the wave drift loads of the ship calculated based on TEBEM can meet the requirements. Figure 16 also shows that the existence of drift angle has little influence on the 6-DOF motion of the ship while having a noticeable effect on the ship's roll and yaw resonance regions; furthermore, the amplitude near the resonance point of roll and yaw decreases with the increase of drift angle. Figure 17 shows that the drift angle influences the calculation of wave drift loads. Especially when λ/LPP < 1.0, the surge and sway drift forces have considerable changes. The surge drift force decreases, while the sway drift force increases with the increase in the drift angle. This is because the disturbance of the flow field around the ship is more evident in the short-wave region. Figure 17 also illustrates that the existence of drift angle has little influence on the yaw moment.

Wave Drift Loads Analysis of Ships with Different Drift Angles in Oblique Waves
Presently, the wave drift loads are mainly studied with zero drift angle. However, the existence of a drift angle will have a certain influence on the wave drift loads in the turning motion of the ship. Therefore, it is critical to study the wave drift loads of the ship with a drift angle in oblique waves. As described in this section, the RAO of 6-DOF motion and wave drift loads of the KVLCC2 ship are first calculated via TEBEM with H = 0.02 L PP , χ = 120 • , β = 0 • , 5.7 • , 11.3 • , 18.4 • , and U = 6.0 kn. Then, the RAO of 6-DOF motion and wave drift loads of the KVLCC2 ship are calculated via the CFD method with H = 0.02 L PP , χ = 120 • , β = 5.7 • , 18.4 • , and U = 6.0 kn, as shown in Figures 16 and 17. The comparison between the TEBEM and CFD methods shows that the 6-DOF motion, surge, and sway wave drift forces of the ship with β = 5.7 • and 18.4 • are in good agreement, whereas the yaw drift moments have large discrepancies in the short waves. This is because significant viscous effects due to flow separation occurred in the bow and stern areas, resulting in a small yaw moment calculated by the CFD method. Compared with the CFD results, the accuracy of the wave drift loads of the ship calculated based on TEBEM can meet the requirements.

Wave Drift Loads Analysis of Ships with Different Drift Angles in Oblique Waves
Presently, the wave drift loads are mainly studied with zero drift angle. However, the existence of a drift angle will have a certain influence on the wave drift loads in the turning motion of the ship. Therefore, it is critical to study the wave drift loads of the ship with a drift angle in oblique waves. As described in this section, the RAO of 6-DOF motion and wave drift loads of the KVLCC2 ship are first calculated via TEBEM with H = 0.02 LPP, χ = 120°, β = 0°, 5.7°, 11.3°, 18.4°, and U = 6.0 kn. Then, the RAO of 6-DOF motion and wave drift loads of the KVLCC2 ship are calculated via the CFD method with H = 0.02 LPP, χ = 120°, β = 5.7°, 18.4°, and U = 6.0 kn, as shown in Figures 16 and 17. The comparison between the TEBEM and CFD methods shows that the 6-DOF motion, surge, and sway wave drift forces of the ship with β = 5.7° and 18.4° are in good agreement, whereas the yaw drift moments have large discrepancies in the short waves. This is because significant viscous effects due to flow separation occurred in the bow and stern areas, resulting in a small yaw moment calculated by the CFD method. Compared with the CFD results, the accuracy of the wave drift loads of the ship calculated based on TEBEM can meet the requirements. Figure 16 also shows that the existence of drift angle has little influence on the 6-DOF motion of the ship while having a noticeable effect on the ship's roll and yaw resonance regions; furthermore, the amplitude near the resonance point of roll and yaw decreases with the increase of drift angle. Figure 17 shows that the drift angle influences the calculation of wave drift loads. Especially when λ/LPP < 1.0, the surge and sway drift forces have considerable changes. The surge drift force decreases, while the sway drift force increases with the increase in the drift angle. This is because the disturbance of the flow field around the ship is more evident in the short-wave region. Figure 17 also illustrates that the existence of drift angle has little influence on the yaw moment.    Figure 16 also shows that the existence of drift angle has little influence on the 6-DOF motion of the ship while having a noticeable effect on the ship's roll and yaw resonance regions; furthermore, the amplitude near the resonance point of roll and yaw decreases with the increase of drift angle. Figure 17 shows that the drift angle influences the calculation of wave drift loads. Especially when λ/L PP < 1.0, the surge and sway drift forces have considerable changes. The surge drift force decreases, while the sway drift force increases with the increase in the drift angle. This is because the disturbance of the flow field around the ship is more evident in the short-wave region. Figure 17 also illustrates that the existence of drift angle has little influence on the yaw moment.
The wave drift forces and moment of the ship with drift angle at different wave incident angles are also calculated. Figure 18a- Figure 18c shows that the yaw wave drift moment is negative in stern waves while positive in bow waves.

Time History of Turning Motion in Regular Waves
After the wave drift loads of the ship are solved, the ship's turning motion in regular waves is calculated. The model ship speed is 1.179 m/s, corresponding to the full ship scale of 15.5 knots in calm water. The propeller speed of the model ship is fixed at 11.8 RPS in all runs. The rudder rate of the model ship is 15.8°/s, which corresponds to the full ship scale of 2.34°/s. The ship's speed decreases in waves due to the added resistance. The approach speed discussed in this section comes from the experimental value of Kim et al. [8]. The related coefficients and hydrodynamic derivatives used in the numerical simulation are shown in Table 2 and derived from Yasukawa et al. [27].
This section describes the simulation of the turning trajectory of the KVLCC2 model with δ = ±35° in calm water and its comparison with the experimental values of Kim et al. [8], as shown in Figure 19. The simulated trajectories generally agree with the experimental results, although the tactical diameter is slightly larger than the experimental value.

Time History of Turning Motion in Regular Waves
After the wave drift loads of the ship are solved, the ship's turning motion in regular waves is calculated. The model ship speed is 1.179 m/s, corresponding to the full ship scale of 15.5 knots in calm water. The propeller speed of the model ship is fixed at 11.8 RPS in all runs. The rudder rate of the model ship is 15.8 • /s, which corresponds to the full ship scale of 2.34 • /s. The ship's speed decreases in waves due to the added resistance. The approach speed discussed in this section comes from the experimental value of Kim et al. [8]. The related coefficients and hydrodynamic derivatives used in the numerical simulation are shown in Table 2 and derived from Yasukawa et al. [27].
This section describes the simulation of the turning trajectory of the KVLCC2 model with δ = ±35 • in calm water and its comparison with the experimental values of Kim et al. [8], as shown in Figure 19. The simulated trajectories generally agree with the experimental results, although the tactical diameter is slightly larger than the experimental value.   Figure 20 shows the time history of the heading angle, total, surge, sway, and yaw velocity compared with the experimental values of Kim et al. [8]. As shown in Figure 20a, the variation trend of the heading angle is consistent with the experimental value, but the turning period is slightly smaller than the experimental value. The surge, sway, yaw, and total velocity are periodical changes, as shown in Figure 20b-e. It can be found that the total velocity of the ship is almost consistent with the experimental values and the variation trend of the surge, sway, and yaw velocity are also consistent with the experimental values. The total velocity decreases obviously at 800-900 s and 1600-1700 s, indicating that the speed loss is more significant in bow waves, as shown in Figure 20a. The comparison results show that the numerical model used has sufficient reliability. Subsequently, the turning motion of the KVLCC2 ship is simulated with wave height H = 0.02 L pp , wavelength λ = 1.0 L pp , wave heading χ = 180 • , and rudder angle δ = +35 • . Corresponding to the full-scale, the numerical simulation time is 2000 s. Figure 20 shows the time history of the heading angle, total, surge, sway, and yaw velocity compared with the experimental values of Kim et al. [8]. As shown in Figure 20a, the variation trend of the heading angle is consistent with the experimental value, but the turning period is slightly smaller than the experimental value. The surge, sway, yaw, and total velocity are periodical changes, as shown in Figure 20b-e. It can be found that the total velocity of the ship is almost consistent with the experimental values and the variation trend of the surge, sway, and yaw velocity are also consistent with the experimental values. The total velocity decreases obviously at 800-900 s and 1600-1700 s, indicating that the speed loss is more significant in bow waves, as shown in Figure 20a. The comparison results show that the numerical model used has sufficient reliability. Figure 21 shows the time history of the ship's 6-DOF motion induced by first-order wave forces. It can be seen from Figure 21 that the 6-DOF motion of the ship changes periodically during each turning cycle and can be stably simulated for 2000 s. The roll resonance appears near 557 and 1308 s, and the ship is in beam waves at this time, as shown in Figure 21d. Figure 21b,c show that the sway and heave reach maximum values near the beam waves and minimum values in the heading and following waves. However, the minimum values of yaw motion are reached in the heading, beam, and following waves. Figure 22 shows the time history of the surge, sway, and yaw wave drift forces and moment of the KVLCC2 ship. The surge wave drift force is large near 560 and 1330 s and small near 400 and 1200 s; at that time, the ship is in bow and stern wave conditions, respectively. The sway wave drift force reaches its maximum value at 178, 940, and 1710 s, respectively, and the ship is in the beam waves. The yaw moment is small around 400 s and 800 s, and the ship is in the head and following waves, while it is larger in the bow and stern waves.
the time history of the heading angle, total, surge, sway, and yaw velocity compared with the experimental values of Kim et al. [8]. As shown in Figure 20a, the variation trend of the heading angle is consistent with the experimental value, but the turning period is slightly smaller than the experimental value. The surge, sway, yaw, and total velocity are periodical changes, as shown in Figure 20b-e. It can be found that the total velocity of the ship is almost consistent with the experimental values and the variation trend of the surge, sway, and yaw velocity are also consistent with the experimental values. The total velocity decreases obviously at 800-900 s and 1600-1700 s, indicating that the speed loss is more significant in bow waves, as shown in Figure 20a. The comparison results show that the numerical model used has sufficient reliability.  Figure 21 shows the time history of the ship's 6-DOF motion induced by first-order wave forces. It can be seen from Figure 21 that the 6-DOF motion of the ship changes periodically during each turning cycle and can be stably simulated for 2000 s. The roll resonance appears near 557 and 1308 s, and the ship is in beam waves at this time, as shown in Figure 21d. Figure 21b,c show that the sway and heave reach maximum values near the beam waves and minimum values in the heading and following waves. However, the minimum values of yaw motion are reached in the heading, beam, and following waves. Figure 22 shows the time history of the surge, sway, and yaw wave drift forces and moment of the KVLCC2 ship. The surge wave drift force is large near 560 and 1330 s and small near 400 and 1200 s; at that time, the ship is in bow and stern wave conditions, respectively. The sway wave drift force reaches its maximum value at 178, 940, and 1710 s, respectively, and the ship is in the beam waves. The yaw moment is small around 400 s waves. Figure 22 shows the time history of the surge, sway, and yaw wave drift forces and moment of the KVLCC2 ship. The surge wave drift force is large near 560 and 1330 s and small near 400 and 1200 s; at that time, the ship is in bow and stern wave conditions, respectively. The sway wave drift force reaches its maximum value at 178, 940, and 1710 s, respectively, and the ship is in the beam waves. The yaw moment is small around 400 s and 800 s, and the ship is in the head and following waves, while it is larger in the bow and stern waves.

Simulations of the Turning Trajectory in Regular Waves
The turning motion of the KVLCC2 ship in regular waves with different wavelengths, wave steepness, and wave headings is simulated, which is compared with the experimental values of Kim et al. [8]. The numerical simulation conditions are shown in Table 3.  24 compares the turning trajectory of +35° and −35° turns in head waves, respectively, with wave height H = 0.02 LPP and wavelength λ = 0.5, 0.7, 1.0, and 1.2. As described in this section, 150 wave periods are stably simulated. As shown in Figures 23 and 24, the turning trajectory of the ship drifts along the direction of wave propagation and perpendicular to the wave, and the trend of the turning trajectory is consistent with the experimental value of Kim et al. [8]. Figure 23 shows the trajectory drifting distance is the longest at λ/LPP = 0.5, and the drifting direction is almost parallel to the wave propagation direction. This is because the surge and sway wave drift forces increase rapidly with the decrease in wavelength. When λ/LPP = 1.0, the relation angle between the trajectory drift direction and the wave propagation direction is the largest. When λ/LPP = 1.2, the drift distance is not significant because of the small wave drift loads. Figure 24 shows the drifting distances of −35° turns are close to the +35° turns. However, when the wavelength is 0.5 LPP, the numerical drifting distances have a slight discrepancy in +35° and −35° turns.

Simulations of the Turning Trajectory in Regular Waves
The turning motion of the KVLCC2 ship in regular waves with different wavelengths, wave steepness, and wave headings is simulated, which is compared with the experimental values of Kim et al. [8]. The numerical simulation conditions are shown in Table 3. turns in head waves, respectively, with wave height H = 0.02 L PP and wavelength λ = 0.5, 0.7, 1.0, and 1.2. As described in this section, 150 wave periods are stably simulated. As shown in Figures 23 and 24, the turning trajectory of the ship drifts along the direction of wave propagation and perpendicular to the wave, and the trend of the turning trajectory is consistent with the experimental value of Kim et al. [8]. Figure 23 shows the trajectory drifting distance is the longest at λ/L PP = 0.5, and the drifting direction is almost parallel to the wave propagation direction. This is because the surge and sway wave drift forces increase rapidly with the decrease in wavelength. When λ/L PP = 1.0, the relation angle between the trajectory drift direction and the wave propagation direction is the largest. When λ/L PP = 1.2, the drift distance is not significant because of the small wave drift loads. Figure 24 shows the drifting distances of −35 • turns are close to the +35 • turns. However, when the wavelength is 0.5 L PP , the numerical drifting distances have a slight discrepancy in +35 • and −35 • turns.
to the wave propagation direction. This is because the surge and sway wave drift forces increase rapidly with the decrease in wavelength. When λ/LPP = 1.0, the relation angle between the trajectory drift direction and the wave propagation direction is the largest. When λ/LPP = 1.2, the drift distance is not significant because of the small wave drift loads. Figure 24 shows the drifting distances of −35° turns are close to the +35° turns. However, when the wavelength is 0.5 LPP, the numerical drifting distances have a slight discrepancy in +35° and −35° turns.   Figure 25 shows a noticeable discrepancy between the numerical simulation and experimental results when λ/LPP = 0.5. The main reason is that the model of the rudders and propellers in calm water are used in waves, ignoring the unsteady effect of waves on rudders and propellers. Especially in the short wavelength   Figure 25 shows a noticeable discrepancy between the numerical simulation and experimental results when λ/LPP = 0.5. The main reason is that the model of the rudders and propellers in calm water are used in waves, ignoring the unsteady effect of waves on rudders and propellers. Especially in the short wavelength region where the wave period is short; the unsteady effect of rudder and propeller is more  Figure 25 shows a noticeable discrepancy between the numerical simulation and experimental results when λ/L PP = 0.5. The main reason is that the model of the rudders and propellers in calm water are used in waves, ignoring the unsteady effect of waves on rudders and propellers. Especially in the short wavelength region where the wave period is short; the unsteady effect of rudder and propeller is more obvious. When λ/L PP = 0.7, 1.0, and 1.2, the trajectory drifting direction and the trajectory drifting distance of the ship are close to the experimental values. Figure 26 shows that the agreement between the numerical results and the experimental results is best when λ/L PP = 1.0. In summary, the drift tendency of turning trajectories can be obtained via numerical simulation, but some discrepancies remain in the predicted drift distance. According to the comparison results, when λ/L PP = 0.7, 1.0, and 1.2, the ship's turning trajectory predicted by the two-time scale model based on TEBEM in waves is acceptable. Figure 27 shows the turning trajectories of −35 • turns in head waves with wavelength λ = 1.0 L PP and wave height H = 0.01, 0.015, and 0.02 L PP . Figure 28 shows the time history of the surge, sway, and yaw wave drift forces and moment of the KVLCC2 under −35 • turns with a variation of the wave steepness in the head waves. The wave drift loads become larger with the increase in wave steepness, resulting in a significant increase in the drifting distances of the turning trajectory. Figure 27 also shows that the relative drifting angle between the trajectory drift and wave propagation directions also increases slightly. Figure 29 shows the turning trajectories of +35 • turns with wavelength λ = 1.0 L PP , wave height H = 0.02 L PP , and wave headings χ = 90 • , 120 • , 180 • , and 210 • . As shown in Figure 29, the turning trajectories of the ship are obviously different with the change of wave headings, and the trajectory drifting direction of the ship is in good agreement with the experimental values, but there are discrepancies between the trajectory drifting distances and the experimental values. Table 4 shows the advances and the tactical diameters of ships at different wave headings. It can be seen from Table 4 that the ship's advances and the tactical diameters differ as the direction of wave propagation changes.     Figure 28 shows the time history of the surge, sway, and yaw wave drift forces and moment of the KVLCC2 under −35° turns with a variation of the wave steepness in the head waves. The wave drift loads become larger with the increase in wave steepness, resulting in a significant increase in the drifting distances of the turning trajectory. Figure 27 also shows that the relative drifting angle between the trajectory drift and wave propagation directions also increases slightly.     Figure 29 shows the turning trajectories of +35° turns with wavelength λ = 1.0 LPP, wave height H = 0.02 LPP, and wave headings χ = 90°, 120°, 180°, and 210°. As shown in Figure 29, the turning trajectories of the ship are obviously different with the change of wave headings, and the trajectory drifting direction of the ship is in good agreement with the experimental values, but there are discrepancies between the trajectory drifting distances and the experimental values. Table 4 shows the advances and the tactical diameters of ships at different wave headings. It can be seen from Table 4 that the ship's advances and the tactical diameters differ as the direction of wave propagation changes.   Figure 29 shows the turning trajectories of +35° turns with wavelength λ = 1.0 LPP, wave height H = 0.02 LPP, and wave headings χ = 90°, 120°, 180°, and 210°. As shown in Figure 29, the turning trajectories of the ship are obviously different with the change of wave headings, and the trajectory drifting direction of the ship is in good agreement with the experimental values, but there are discrepancies between the trajectory drifting distances and the experimental values. Table 4 shows the advances and the tactical diameters of ships at different wave headings. It can be seen from Table 4 that the ship's advances and the tactical diameters differ as the direction of wave propagation changes.

Analysis of Trajectory Drifting Angle and Distances
The trajectory drifting angle and distances defined by Kim et al. [8] are used to analyze the data of the ship's turning trajectory, and the data are compared with the experimental values of Kim et al. [8]. The trajectory drifting distance Ddr 090-450 is defined as the vector magnitude between two points corresponding to the heading angles of 90 • and 450 • , and the trajectory drifting angle Adr 090-450 is the angle between the initial approach direction of the ship and that vector. Other definitions of trajectory drifting angles and distances are similar, where the subscript represents the heading angles at the start and end positions. Tables 5 and 6 compare the trajectory drifting angle and distances in head waves. When λ/L PP = 0.5 and 0.7, the discrepancies in trajectory drifting angle and distance between the numerical results and the experimental results are small, while the numerical discrepancies increase gradually in the long-wave region. Tables 7 and 8 compare the trajectory drifting angles and distances in beam waves. The drifting trajectory and drifting angle discrepancies in beam waves are similar to those in head waves. This is probably because of the large error of the numerical results of wave drift loads in the long-wave region compared with the experimental values. Furthermore, the load model of the propeller and rudder is based on the force model under calm water without considering the non-uniformity and unsteady characteristics of the flow field on the propeller surface caused by waves and the 6-DOF motion of the ship, and the influence of the unsteady flow field on the propeller load. Therefore, the propeller rudder load model in waves will be considered in the future. Tables 5-8 show that as the heading angle of ships increases, the error of the trajectory drifting angle and distance also increases. The reason is that the numerical discrepancies gradually accumulate and increase as the numerical simulation time increases. Although the discrepancies in trajectory drifting angles and distances are large at the long wavelength, the accuracy of numerical simulation is still acceptable compared with the experimental values.

Conclusions
This study applies the two-time scale model based on TEBEM to solve the ship's turning motion in regular waves. The total ship motion is broken down into high-frequency wave-induced motion and low-frequency maneuvering motion. The TEBEM solves the high-frequency seakeeping problem, and the MMG model calculates the low-frequency maneuvering motion. Finally, the wave drift loads are included in the low-frequency maneuvering motion. The major conclusions can be summarized as follows: 1.
The drift angle has a certain effect on the wave drift loads of the ship. In particular, when λ/L PP ≤ 1.0, the surge drift force gradually decreases while the sway drift force gradually increases with the increase of the drift angle.

2.
Although there are some discrepancies between the numerical and experimental results, the numerical simulation results are still acceptable for the turning trajectory of the KVLCC2 model. The simulated turning trajectories of the ship show that the incident wave significantly influences the ship's turning motion, and the drifting distance increases with the increase in wave steepness. When λ/L PP = 1.0, the relation angle between trajectory drifts and wave propagation directions is the largest. The numerical results also show that the turning trajectories of the ship change significantly with the change in wave heading.

3.
The two-time scale model based on TEBEM can effectively solve the turning trajectories of ships in regular waves using maneuvering coefficients and hydrodynamic derivatives in calm water.
In the future, we will further investigate the unsteady flow field distribution and load on the propeller surface to improve the accuracy of the ship's turning trajectory in waves.

Data Availability Statement:
The data presented in this study were available in this article (Tables  and Figures).