A Numerical Study on the Performance of the 66k DWT Bulk Carrier in Regular and Irregular Waves

Abstract

Accurate resistance performance evaluation is essential to predict the minimum propulsion horsepower required for ships in real sea. Therefore, in this paper, the added resistance performance of K-Supramax in various wave conditions was calculated using CFD (Computational Fluid Dynamics). First, to compare the resistance performance according to the DOF (Degree of Freedom), the added resistance performance in the head regular wave was investigated and compared with the experiment. The 2DOF and 6DOF results in the head regular wave had no significant differences. Second, to compare the resistance performance among various wave conditions, the effect of the heading angle was analyzed by comparing resistance and motion RAOs (Response Amplitude Operators) in bow quartering sea conditions using 6DOF. The added resistance showed the highest value near the resonance frequency of 170°. As the heading angle increased, the added resistance tended to decrease, regardless of the wavelength ratio. Also, the added resistance performance in irregular waves was investigated with reference to the adverse conditions. The spectral method, which linearly superimposes regular wave results, and the added resistance generated during irregular waves directly through CFD were compared. The resistance in irregular waves increased by approximately 92% for the spectral method and by approximately 72% for the direct irregular wave calculation compared to the resistance in calm water. In addition, the PSD (Power Spectral Density) of resistance and the distribution of motion RAOs obtained from direct irregular wave calculations were investigated.

1. Introduction

At the 59th Meeting of the MEPC (Marine Environment Protection Committee) in 2009, it was announced that carbon dioxide emitted by ships is significant and that regulation of this is necessary. Accordingly, the IMO (International Maritime Organization) introduced the EEOI (Energy Efficiency Operational Indicator) and the EEDI (Energy Efficiency Design Index) to reduce CO₂ emissions. One way to reduce CO₂ emissions is to operate at a low operation speed by reducing the MCR (Maximum Continuous Rating) of the ship. In this way, the problems related to ship stability occurred. The MEPC proposed guidelines for minimum propulsion power to maintain the maneuverability of ships under adverse conditions [1]. An accurate ship performance estimation is essential to enhance ship stability under adverse conditions. One of the ways to improve the operational efficiency of ships is to reduce the added resistance caused by waves, which account for a large proportion of the resistance increase. Prpić-Oršić and Faltinsen [2] predicted that CO₂ emissions would increase by more than 15 percent compared to calm water when S-175 container ships sail in the North Atlantic seaway. There may be a difference between the current ship performance estimation process and the actual seaway, and one of the causes of this is added resistance. Several studies have been conducted on added resistance, which can be divided largely into experimental and numerical studies.

At first, early studies on added resistance were conducted by experiments. Denis and Pierson [3] derived motion responses to ships from various waves, suggesting that a ship’s motion can superimpose. Lewis and Bennet [4] predicted the ship motion response in irregular waves through this superimposing principle of ship motion. They suggested that the vertical acceleration of the motion of the ship received by waves can be expressed by combining the acceleration for heave and pitch motion with vectors. In addition, the expected motion of a ship in irregular waves was compared through actual ship data. Gerritsma et al. [5] reported that the main cause of added resistance by waves were the radiation waves generated by the ship motion and the diffraction force bounced by the bow shape of the ship. Since then, Strom-Tejen et al. [6] showed, using Series 60 ships, that ship motion is mainly caused by heave and pitch motion, and suggested that damping caused by fluid viscosity is negligible. Fujii and Takahashi [7] directly observed and analyzed the effect of diffraction force bounced by the bow shape of the S175 container ship. Journee [8] published experimental data on the heave and pitch motion, wave load, and added resistance through the Wigley ship. The effects of various waves and the shape of ships have been studied experimentally, including short wavelengths [9], various wavelengths [10], and the effects of bow shape [11].

To predict the added resistance, various theoretical methods using potential theory have been proposed based on the experimental study results. Far-field and near-field methods are two representative methods used to calculate the added resistance. Maruo [12] proposed an energy or momentum conservation method known as the far-field method, which has the advantage of relatively simple calculation because it does not directly calculate the pressure. On the other hand, Gerritsma and Beukhelman [13] suggested that the far-field method could predict resistance when calculating the short wavelength. The pressure integrating method proposed by Boese [14], known as the near-field method, directly integrates and calculates the secondary fluid pressure. Hence, the calculation process is relatively complicated, but a physical interpretation is convenient. These methods were applied mainly to the 2D strip method. Since then, new methods, such as the 3D panel method, have been proposed.

Recently, to study the performance of ships in waves, numerous studies have been conducted through CFD (Computational Fluid Dynamics), which performs repetitive calculations by solving the Navier-Stokes equation considering viscosity. The CFD can identify various fluid phenomena, which are hard to observe by experiment. Sadat-Horseini et al. [15] computed the added resistance for KVLCC2 (KRISO Very Large Crude Oil Carrier) in regular waves using CFD, which compared the difference between the 2DOF (heave- and pitch-free) and 3DOF (surge-, heave-, and pitch-free). Park et al. [16] estimated the ship power required to achieve a specific speed by performing CFD in irregular wave conditions of Beaufort 8. This study considered the SRtP(Safe Return to Port) regulations presented by Lloyd’s register for the twin-axis car ferry, in which the propeller on the port is not working due to engine failure. Shin et al. [17] examined the nominal wake in regular head waves for KVLCC2. Eom et al. [18] investigated the effect of the immersion depth and tilt angle of a propeller on the propeller performance in waves. Lee et al. [19,20] assessed self-propulsion performance methods that change with wave conditions using KVLCC2. K-Supramax used in this study was analyzed in various features. For example, Lee et al. [21,22] examined the effects of the diffraction force and compared the effects of bow shape. Also, Lee et al. [23] presents a comparative study of wave-induced motion responses and the added resistance of four representative types of ships. The effects of wave direction, such as head sea and oblique sea conditions, are considered, and added resistance performance is compared [24].

This study examined the added resistance acting on K-Supramax (scale ratio: 24) advancing in regular and irregular wave conditions using numerical simulations. A free decay test was performed to determine the resonance period of the ship to analyze the resistance performance in waves, and the effect of DOF (Degree of Freedom) was compared with the theory. Further, 1, 2, and 6DOFs were considered, and 6DOF were implemented through a soft-spring system. To validate results, resistance performance in calm water and head regular waves were investigated and compared with the experiment. Also, the resistance performance according to the DOF was analyzed. Based on these results, the resistance performance among various wave conditions were investigated. The 6DOF computations were performed on the bow quartering sea, with 10° intervals up to 150°, as well as the head sea. Using the regular wave results in the head sea, the spectral method was applied to estimate the added resistance in irregular wave conditions. Also, direct CFD computation was conducted. To do this, irregular waves were created by regular wave superposition, and the parameters included in irregular wave generation were tested. In addition, the motion RAOs according to the wave frequency, as well as the added resistance in irregular waves, were compared with the regular wave results.

2. Objective Ship

In this study, computations were conducted with the K-Supramax (66,000 DWT), one of the mid-sized bulk carriers. Figure 1 and

Table 1. Main particulars of K-Supramax bulk carrier.

	Full-Scale	IU & KRISO	PNU [21]	SNU [23]
Scale ratio, $\lambda$	1	24	33.33	60, 40.42
Length between perpendiculars, $L_{pp}\left[\mathrm{m}\right]$	192	8	5.76	3.20, 4.75
Length on waterline, $L_{wl}\left[\mathrm{m}\right]$	196	8.16	5.88	3.27, 4.85
Breadth, $B_{wl}\left[\mathrm{m}\right]$	36	1.50	1.08	0.60, 0.89
Draft, $T\left[\mathrm{m}\right]$	11.2	0.47	0.34	0.19, 0.28
Wetted surface area, $S_{w/o}\left[{\mathrm{m}}^2\right]$	9808	17.03	8.83	2.72, 6.00
Displacement, $\nabla \left[{\mathrm{m}}^3\right]$	65,028	4.70	1.76	0.30, 0.98
Block coefficient, $C_B$	0.8400
Midship section coefficient, $C_M$	0.9973

present the geometry, plans, and main particulars of the K-Supramax, respectively. The computation was performed with a 1:24 scale model, which is the same scale as the experimental results of KRISO (Korea Research Institute of Ships and Engineering). In addition, model tests and numerical analysis results with different scales (1:33.33 scale (PNU) [21], 1:60 scale (captive model, SNU) and 1:40.42 scale (free running model, SNU) [23]) were added to compare and verify.

3. Numerical Method and Condition

3.1. Numerical Set-Up

CFD simulations were performed using STAR CCM+ 15.06 version, a commercial CFD software. The fluid was assumed to be an incompressible viscous fluid. The governing equations are continuous equations, as shown in Equation (1), and in the Navier-Stokes equation, as shown in Equation (2):

$$\frac{\partial {\overline{u}}_i}{\partial {\overline{x}}_j}=0$$

(1)

$$\frac{\partial {\overline{u}}_i}{\partial t}+{\overline{u}}_j\frac{\partial {\overline{u}}_i}{\partial {\overline{x}}_j}=-\frac{1}{\rho }\frac{\overline{p}}{{\overline{x}}_i}+\frac{1}{\rho }\frac{{\partial }^2}{\partial {\overline{x}}_j}\left(\mu \frac{\partial {\overline{u}}_i}{\partial {\overline{x}}_j}-\rho \overline{{u^{\prime }}_i{u^{\prime }}_j}\right)$$

(2)

where $\rho$ and $\nu$ represent the fluid density and kinematic coefficient, and $x_j\left(x,y,z\right)$ and $u_i\left(u,v,w\right)$ are the coordinate and the fluid velocity component in the direction of each coordinate axis. $p$ is pressure.

The implicit unsteady model with a second-order temporal scheme was used as the time discretization method. The semi-implicit method for the pressure-linked equations (SIMPLE) algorithm was used as the pressure-velocity coupling. To set the convective courant number below 1, a value corresponding to 1/350 of the encounter wave period was used as a timestep. As the turbulence model, a realizable k-epsilon model was used, which shows a superior ability to capture the mean flow of complex structures. The realizable k-epsilon model is substantially better than the standard k-epsilon model for many applications [25]. Typical methods for simulating free surface include a VOF (Volume of Fluid) method and a level-set method. Among them, the VOF method is subordinate to the Eulerian method, which simulates multi-phase flow and is suitable for simulating flow around the ship. So, a VOF method was implemented for the free surface and the waves. For the ship motion, the CFD simulations considered 1DOF (free heave or free pitch), 2DOF (heave- and pitch-free), and 6DOF. The DFBI (Dynamic Fluid Body Interaction) model was used to implement the movement of the hull. The overset mesh deals with various wavelengths, configuring the internal and background grids independently. The overset mesh constructs background mesh and separates component mesh. This mesh interpolated cell data in an area where mesh overlap for two or more mesh systems. So, if the movement of the hull increases due to the wave, the overset grid has the advantage of following the wave by fixing the background grid system.

3.2. Grid System and Boundary Condition

As shown in Figure 2a, the symmetrical computations, such as calm water, regular wave, and irregular wave computations in head sea conditions, which consider the 2DOF, were performed by generating a half-domain grid. On the other hand, as the regular wave computations consider the 6DOF, the computation was performed in the grid, as shown in Figure 2b. In the computation domain, the distance to the inlet boundary surface, outlet boundary surface, left/right boundary surface, and top/bottom boundary surface based on the hull is expressed as a multiple of the hull length ($L_{pp}$), as shown in Figure 2. In regular wave and irregular wave computations, the velocity inlet was set at the inlet boundary and outlet boundary as boundary conditions. In addition, the velocity inlet condition was set at the left/right boundary and the top/bottom boundary. A symmetric condition was set at the symmetrical plane for a half-domain computation. This boundary conditions not only minimizes the dissipation of waves, but also minimizes reflection by waves. The number of grids used in the computation was approximately 1.9 million grids for resistance computation for calm water conditions. Various wavelengths as added resistance computation for regular wave computations were considered using 3.4 to 3.8 million grids for half-domain computation and 6.7 to 7.3 million grids for 6DOF computations. For irregular wave computation, 3.9 million grids were used with half domain.

3.3. DOF Simulation Method

In this study, CFD simulations were performed in 1DOF, 2DOF, and 6DOF. To apply the 6DOF, the model ship in the simulation was restrained with a soft spring-type mooring system consisting of four horizontal springs and a virtual towing carriage, as shown in Figure 3. The four virtual springs were attached to the ship bow and stern at the same height as the ship center of gravity position. This soft spring-type mooring system should not affect the motion response of the model ship significantly. The influence of the spring on the motion of the model ship was minimized by selecting a spring with more than six times the natural frequency by applying the spring coefficient. The virtual towing carriage was implemented by producing a coordinate system with the model ship to move at a constant x-direction speed.

4. Validation and Verification

4.1. Free Decay Test

Before analyzing the ship motion in a regular wave, a free heave/pitch decay test was performed to define the natural frequency of motions.

Table 2. Free decay test condition.

	Free Motion
$z/L_{pp}=-0.04$	① Heave ② Heave, Pitch ③ Surge, Sway, Heave, Roll, Pitch, Yaw
$\theta \left[\mathrm{d}\mathrm{e}\mathrm{g}\mathrm{r}\mathrm{e}\mathrm{e}\right]=5°$	① Pitch ② Heave, Pitch ③ Surge, Sway, Heave, Roll, Pitch, Yaw

lists the free decay test conditions. The heave/pitch natural frequency obtained through the free decay test was compared with the empirical formulas as shown in Equations (3) and (4):

$$f_{z,Heave}=\sqrt{\frac{g{C}_{Wp}}{8{\pi }^2{C}_{B}T}}$$

(3)

$$f_{z,Pitch}=\sqrt{\frac{C_{IT}g{B}^3}{96{\pi }^2{\widehat{I}}_y{L}_{pp}^4}}$$

(4)

where $C_{Wp}$ is the water-plane area coefficient, $C_B$ is the block coefficient, $C_{IT}$ is the coefficient of inertia of the water-plane area about the y-axis, $T$ is the draft, g is the gravitational acceleration, and ${\widehat{I}}_y$ is a non-dimensional mass moment of inertia about the y-axis.

Table 3. Natural frequency comparison.

	Empirical Formula	1DOF		2DOF		6DOF
Free condition		Heave	Pitch	Heave Pitch	Heave Pitch	All	All
Initial condition		$z/L_{pp}=-0.04$	$\theta =5°$	$z/L_{pp}=-0.04$	$\theta =5°$	$z/L_{pp}=-0.04$	$\theta =5°$
$f_{n,Heave}\left[\mathrm{H}\mathrm{z}\right]$	0.5438	0.5540		0.5208	0.5222	0.5236	0.5270
$f_{n,Pitch}\left[\mathrm{H}\mathrm{z}\right]$	0.5660		0.5475	0.5042	0.5051	0.5123	0.5097

lists the results of natural frequencies through empirical formulae and free decay tests. The natural frequency derived through the 1DOF is relatively high, and the natural frequency derived through the 2DOF and the 6DOF are similar. In the case of 1DOF, it is judged to be higher due to the restraint conditions, and it is judged that using at least 2DOF will obtain an accurate motion response. Figure 4 compares the encounter frequency by the wavelength and natural frequency derived through the free decay test at 2DOF. The derived natural frequency matches the encounter frequency of the wavelength ratio of 1.4. According to Lewis [26], the effects of the speed of the ship on the natural frequency are significantly small. Therefore, the natural frequency at the design speed is expected to be close to the encounter frequency of the wavelength ratio of 1.4. The validity of the restoring performance of a ship was reviewed by confirming that the theoretical natural frequency and the natural frequency obtained through computations are similar through the free decaying test.

4.2. Calm Water Simulation

The resistance, sinkage, and trim of the K-Supramax in calm water were verified by performing the resistance calculation on nine different ship speeds, including the design speed (14.5 knots, Fr = 0.172). The results on six different speeds of ship were compared with the KRISO experiments at the same scale ratio, as shown in Figure 5. The results were similar to the experiment, with an error of less than 2% for all resistances and less than 4% for all sinkage and trim.

In terms of fluid dynamics, the simulation compared the wake distribution with the experimental results [27], as shown in Figure 6. The computation result shows a clear hook-shaped vortex similar to the experiment. Figure 7 and Figure 8 were compared with the CFD analysis results [22] to compare the pressure distribution and wave pattern of the hull, respectively. As for the pressure distribution, Lee et al. [22] reported that the pressures at the fore, shoulder, and stern shoulder were similar to the CFD result. In addition, Figure 8 shows the distribution of waveforms around the hull because the vessel operates at a relatively slow speed, and the transverse wave is dominant. Lee et al. [22] reported similar results, confirming the reliability of the current CFD result.

5. Regular Wave

5.1. Regular Wave Simulation Condition

The regular wave simulations were performed at Fr = 0.172, which corresponds to the design speed (14.5 knots).

Table 4. Regular wave condition.

Wave Steepness	2DOF	6DOF
$$\begin{gathered} H/L_{pp}=0.00800 \\ {(A/L_{pp})}^{-1}=250 \end{gathered}$$	$\lambda /L_{pp}=0.3,0.4,0.5,0.6,$ $0.7,0.8,0.9,1.0$	$\lambda /L_{pp}=0.5,1.0$
$$\begin{gathered} H/L_{pp}=0.00667 \\ {(A/L_{pp})}^{-1}=300 \end{gathered}$$	$\lambda /L_{pp}=0.2,1.1,1.25,$ $1.4,1.6,2.0$	$\lambda /L_{pp}=1.25,1.4,1.6$

lists the wave conditions used to analyze added resistance for regular waves. The wave height was selected for $H/L_{pp}=0.00800$ and $H/L_{pp}=0.00667$. The total simulation conditions were 14 computations of 2DOF and five computations of 6DOF.

5.2. Regular Wave Test

For accurate wave generation, wave elevation was verified for three wavelength ratios $(\lambda /L_{pp}=0.5,1.0,\mathrm{a}\mathrm{n}\mathrm{d}1.6)$, as shown in Figure 9, which was tested without a ship. The verification compared the generated waves and theoretical waves at the bow and stern position of the ship in the domain. Various grid systems were used to compute various wave conditions accurately. As the grid system was used, approximately 80 to 120 cells were used, depending on the wavelength, and approximately 16 to 24 cells were used, depending on the wave height. The wave height of the generated wave is well generated within a 5% error with the theoretical wave.

5.3. Convergence Grid Test in Regular Wave

A convergence simulation was conducted at $\lambda /L_{pp}=1.0$ and $H/L_{pp}=0.008$ to confirm grid convergence in regular wave condition. Convergence was confirmed through the GCI values [28]. As a simulation condition, nine simulations were conducted with three grids and three timesteps, as shown in

Table 5. Comparison of grid convergence index (GCI).

	Case No.	Number of Grids	${\mathit{T}}_{\mathit{e}}/∆\mathit{t}$	${\mathbf{C}}_{\mathbf{A}\mathbf{W}}$	GCI (%)
Grid	G1T2	2,169,732	350	8.686	0.515
	G2T2	3,687,345	350	8.149
	G3T2	7,514,440	350	8.030
Timestep	G2T1	3,687,345	150	8.823	0.242
	G2T2	3,687,345	350	8.149
	G2T3	3,687,345	700	8.015

. The timestep was expressed based on the encounter wave period. When the GCI value was derived by considering the three timesteps, the convergence was valid at less than 1%. Therefore, the standard timestep was determined to be $T_e/∆t=350$. In addition, using this timestep, GCI values were derived through three grids. The values were valid at less than 1%. Finally, $T_e/∆t=350$ for the timestep and 3.6 million for the grid were used.

5.4. Resistance Performance in a Regular Wave

The computations were performed in the regular wave in the design speed (14.5 knots, Fr = 0.172). The added resistance coefficient and motion RAOs in the regular wave were verified with the experimental and numerical results, as shown in Figure 10. The equations for the added resistance are expressed as Equations (5) and (6):

$$R_{aw}=\overline{R_T}-R_{cw}$$

(5)

$$C_{aw}=\frac{R_{aw}}{\rho g{A}^2{B}^2/L_{pp}}$$

(6)

where $R_{aw}$ is the difference between the mean of the total resistance in waves ($\overline{R_T}$) and resistance in calm water ($R_{cw}$), and $C_{aw}$ is the added resistance coefficient using the wave height ($A$), ship breadth ($B$), and ship length ($L_{pp}$). $\rho$ and $g$ represent the fluid density and gravitational acceleration, respectively.

The added resistance coefficient for each encounter frequency was expressed as a dimensionless value through the natural heave frequency. In the case of the KRISO model, it is indicated by a black symbol in the graph. The model test and numerical analysis with different scales [22] are indicated by a blue line and symbol in the graph. In the case of this study, numerical analysis of 2DOF and 6DOF was shown and indicated by the red symbol in the graph. Figure 10 shows that the tendency of the added resistance coefficient in regular waves was similar to the results of the model test and numerical computations. As the wavelength ratio increased, the added resistance coefficient increased at $\lambda /L_{pp}<1.1$ and decreased after the peak near $\lambda /L_{pp}=1.1$. The highest point was judged to be the resonance wavelength ratio of the ship, and it is the part where the added resistance of the ship acts the most. This was larger than the dimensionless encounter frequency 1. On the other hand, the effect of the DOF was insignificant because the condition is the head sea.

Similar to the added resistance graph, the motion RAOs were compared. Figure 11 and Figure 12 show the heave coefficient and phase. Figure 13 and Figure 14 show the pitch coefficient and phase. The equations for the motion response are expressed as Equations (7) and (8). The CFD results of the heave motion and pitch motion showed high reliability compared with the experimental and numerical results.

$${{\xi }^{\prime }}_3=\frac{{\xi }_3}{A}$$

(7)

$${{\xi }^{\prime }}_5=\frac{{\xi }_5}{kA}$$

(8)

where ${\xi }_3$ is the heave amplitude, ${\xi }_5$ is the pitch amplitude, $A$ is the wave heigh, and $k$ is the wave number.

Figure 11 shows that the amplitude of the heave motion does not change significantly at $\lambda /L_{pp}<0.8$, which means the vertical motion of the hull does not affect the short wavelength significantly. The amplitude increases rapidly, showing a peak at $\lambda /L_{pp}=1.4$, and then decreases slowly at $\lambda /L_{pp}>1.4$. When the wavelength becomes similar to the length of the ship, the vertical movement of the hull increases, and the largest vertical movement is near $\lambda /L_{pp}=1.4$. In addition, the encounter frequency in the section representing the highest point is similar to the natural frequency obtained through the free decay test. In addition, the heave coefficient appears relatively low in the long wavelength.

Figure 12 shows that the heave phase indicates near 30° at $0.5<\lambda /L_{pp}<0.8$, whereas the heave phase change occurs significantly as the wavelength increases at a smaller short-wave section of $\lambda /L_{pp}<0.5$. The heave phase converges slowly to 30° when the encounter frequency increases more than the natural frequency. On the other hand, the large heave phase change occurs when the encounter frequency is more than twice due to nonlinear characteristics, such as the ship’s shape. The heave phase gradually decreases to approximately −90° at $0.8<\lambda /L_{pp}<1.0$, and then the heave phase is increased to converge approximately 0° at $\lambda /L_{pp}>1.0$. Furthermore, when similar to the encounter frequency and the natural frequency, the heave motion is the highest. Through this, it is thought that the vessel motion in the short wavelength region has a different phase than the incoming wave, and thus the motion is suppressed. Figure 13 shows that the amplitude of the pitch motion does not change significantly until $\lambda /L_{pp}=0.8$. The pitch motion also does not affect the movement of the hull significantly in the short wavelength. The amplitude pitch motion continues to increase, unlike the heave motion. The movement of the hull rapidly increases when the wavelength becomes similar to the length of the ship. In addition, Figure 14 shows the pitch phase according to the wavelength and the encounter frequency. The pitch phase changes significantly depending on the wavelength conditions, such as the heave phase. As the wavelength increases, the heave phase converges to 0°, but the pitch phase converges to 45°, so that the motion of the heave and the pitch has a difference of 45°. The same trend appears in other wavelengths. On the other hand, when the wavelength ratio is less than 0.5, the phase does not differ by 45°. As a result, the pitch motion has a phase difference of approximately 90° in the section where the wavelength ratio is smaller than 0.7. The phase changes significantly at the resonance frequency. Finally, a difference of approximately 45° can be seen at the long wavelength. The phase changes according to the wave introduced by the vessel shape, which greatly correlates with the change in motion RAOs. This was confirmed by the time series distribution reported by Lee et al. [29].

Next, the trend of motion and added resistance was analyzed by calculating the relative motion of the bow and stern of the hull. Figure 15 and Figure 16 show the relative motion of the bow (${RM}_{Bow}$) and stern (${RM}_{Stern})$, which can be estimated as Equations (9) and (10):

$${RM}_{Bow}=z_0+z_1\cos \left(2\pi f_{e}t+z_{\epsilon 1}\right)+x_G\sin \left({\theta }_0+{\theta }_1\cos \left(2\pi f_{e}t+{\theta }_{\epsilon 1}\right)\right)-A\cos \left(2\pi f_{e}t+{\gamma }_b\right)$$

(9)

$${RM}_{Stern}=z_0+z_1\cos \left(2\pi f_{e}t+z_{\epsilon 1}\right)-(L-x_G)\sin \left({\theta }_0+{\theta }_1\cos \left(2\pi f_{e}t+{\theta }_{\epsilon 1}\right)\right)-A\cos \left(2\pi f_{e}t+{\gamma }_s\right)$$

(10)

where $(z_0,z_1,z_2)$ and $\left({\theta }_0,{\theta }_1,{\theta }_2\right)$ are the mean value, the 1st amplitude, and the corresponding phase of heave and pitch. ${\gamma }_b$ and ${\gamma }_s$ are an incident wave phase at the bow and stern. $t=0\mathrm{s}$, respectively [15].

As shown in Figure 15, the relative motion at the bow does not change significantly in the short wavelength. The relative motion begins to increase near $\lambda /L_{pp}=0.8$ and decreases after the peak at $1.25<\lambda /L_{pp}<1.4$ near the natural frequency. The trend of the relative motion at the bow is similar to the added resistance, which considers that the relative motion of the bow influences the added resistance. On the other hand, in Figure 16, these relative motions at the stern are similar to the heave and pitch motion responses. Because the motion response shows the tendency of the ship’s center of gravity to move, the similarity of the relative motion at the stern means that the movement of the hull center and the relative motion at the stern is similar. As a result, it can be estimated that the added resistance is influenced by the bow movement, and the stern moves relatively along the motion response.

In this study, the 6DOF condition was adopted, which includes the surge motion. Figure 17 shows the surge amplitude time series in a soft spring-type mooring system. Because the spring exists, the surge amplitude affected by the spring appears, such as the blue line. The surge amplitude using FFT analysis is shown as a red line to predict the surge amplitude without a spring effect. Figure 18 shows the surge coefficient, calculated as Equation (11):

$${{\xi }^{\prime }}_1=\frac{{\xi }_1}{A}$$

(11)

where ${\xi }_1$ is the surge amplitude, and $A$ is the wave height.

The surge coefficients increase as the wavelength increases, and it is the largest when the encounter frequency is similar to the natural frequency, which is derived by the free decay test. In particular, the surge coefficients are relatively small, less than 0.3. In addition, Figure 10 shows that the difference is not large when comparing the added resistance coefficients of the 2DOF and 6DOF, which means surge motion is too small to affect the added resistance.

5.5. Resistance Performance in the Regular Wave with Various Heading Angles

For analysis of added resistance performance according to wave direction, the added resistance simulation was tested according to the bow quartering wave directions. The same wavelength as the 6DOF calculation in

Table 6. Irregular wave conditions in adverse conditions [1].

$\mathbf{Ship}\mathbf{Length},{\mathit{L}}_{\mathit{p}\mathit{p}}\left[\mathbf{m}\right]$	$\mathbf{Significant}\mathbf{Wave}\mathbf{Height},{\mathit{h}}_{\mathit{s}}\left[\mathbf{m}\right]$	$\mathbf{Peak}\mathbf{Wave}\mathbf{Period},{\mathit{T}}_{\mathit{p}}\left[\mathbf{s}\right]$	$\mathbf{Mean}\mathbf{Wind}\mathbf{Speed},{\mathit{V}}_{\mathit{w}}[\mathbf{m}/\mathbf{s}]$
Less than $200$	4.5	7.0 to 15.0	19.0
$200\le L_{pp}\le 250$	Parameters linearly interpolated, depending on the ship length
More than $250$	6.0	7.0 to 15.0	22.6

was simulated, and the wave direction was considered at 10° intervals from 180° to 150°. Added resistance performance and motion RAOs depending on the heading angles are compared in Figure 19.

In the case of added resistance, the added resistance coefficient value in the peak and long-wave region was high at heading angle 170°, but the added resistance in the short-wave region tended to be lower than that at 180°. The added resistance can be greater in the bow quartering sea than in the bow wave condition. ISO 15016:2015 [30] assumption of calculating the same added resistance up to 150° is used, and since the difference in added resistance occurs at the heading angle in this way, it is necessary to consider the difference according to the heading angle. After the added resistance showed a peak, it tended to decrease as the heading angle increased. It is judged that the encounter wave period and the speed felt by the ship change according to the wave direction, and the added resistance decreases due to the difference in resonance frequency and speed decrease for each heading angle. In the case of Heave RAOs, a tendency to rapidly increase after the condition of a wavelength ratio of 1.0 was shown in all wave directions, and a slight difference was shown depending on the wave direction at a wavelength ratio larger. In particular, the difference was cut under the condition of $\lambda /L_{PP}$ = 1.2, which is the resonance area, and it was the highest at the wave direction of 150°. In the case of Pitch RAOs, it was similar regardless of the wave direction, up to the wavelength ratio of 1.4, and then it appeared higher as the wave direction was closer to the head sea condition. Roll RAOs were insignificant in the forward bow condition considered in this study.

6. Irregular Wave

6.1. Resistance Performance in Irregular Waves Using the Spectral Method (Indirect Method)

Maruo [31] proposed a spectral method based on the theory that regular waves can superimpose linearly. The force (resistance, thrust, and torque) and motion response in irregular waves can be estimated using this linear superposition. In this method, the variables were assumed to be proportional to the square of the wave height for the superposition expressed as a quadratic transfer function. The added resistance in the regular wave can be determined by dividing the square of the wave amplitude to obtain the transfer function, as shown in Equation (12). In addition, a spectrum of waves is required because the sea state is considered. Equation (12) can be multiplied by the spectrum of waves to calculate the added resistance in irregular waves by applying the spectral method, as shown in Equation (13):

$$H\left(V_s,\omega \right)=\frac{R_{aw}(V_s,\omega )}{{\xi }_a^2}$$

(12)

$$R_{irregular}=2{\int }_0^{\infty }\frac{R_{aw}\left(V_s,\omega \right)}{{\xi }_a^2}S\left(\omega \right)d\omega$$

(13)

where $R_{aw}$ is the added resistance in regular waves; $V_s$ is the ship speed; $\omega$ is the frequency; ${\xi }_a$ is the wave amplitude; $S\left(\omega \right)$ is the wave spectrum.

In this study, the added resistance in irregular wave conditions was estimated in adverse conditions, which are presented in MEPC.1/Circ.850/Rev.3 [1], as shown in Table 6. The significant wave height was determined to be 4.5 m, and three peak wave periods were chosen: 10, 12, and 14 s.

Figure 20 shows the wave spectrum according to the three-peak wave period. In the wave spectrum, the peak shifts to a lower frequency as the peak wave period increases. Figure 21 shows the product process to estimate the added resistance in the irregular waves through the spectral method. The red line represents the JONSWAP spectrum, the blue symbol represents the result of added resistance divided by the square of the wave height, and the green line represents the product of the JONSWAP spectrum transfer function and the added resistance transfer function. Finally, the added resistance in irregular waves is estimated by integrating the spectrum of the green line.

Table 7. Estimated resistance in irregular waves using the spectral method.

Peak Wave Period ${\mathit{T}}_{\mathit{p}}\left[\mathbf{s}\right]$	Total Resistance ${\mathit{R}}_{\mathit{T},\mathit{i}\mathit{r}\mathit{r}\mathit{e}\mathit{g}\mathit{u}\mathit{l}\mathit{a}\mathit{r},\mathit{S}\mathit{M}}\left[\mathbf{N}\right]$
10	138.69
12	151.43
14	127.24

lists the derived resistance in irregular waves with a significant wave height ($H_s$) of 4.5 m and a ship speed ($V_s$) of 14.5 knots. The added resistance in irregular waves represents the highest point when the peak wave period is 12 s. Here, the added resistance is large because the peak frequency of the transfer function is similar to the peak frequency of the added resistance in the regular wave. The different peak periods of 10 and 14 s are due to the reduced multiplication ratio of the wave spectrum and transfer function. As the magnitude of the added resistance can vary according to the peak wave period of the sea area, it is necessary to select the seaway considering the natural frequency of the ship.

6.2. Resistance Performance in an Irregular Wave through Direct Computation

In this study, the wave test of irregular waves was conducted to generate a more accurate irregular wave, which is the same adverse condition considered in the spectral method (Section 6.1). The test condition chose the peak wave period of 12 s.

$$\eta =Asin(\omega t+\varphi )$$

(14)

Basically, the wave height ($A$), the angular frequency ($\omega$), and the phase angle ($\varphi$) are used as input parameters to generate the regular wave. Equation (14) is the theoretical expression of the regular wave. In this study, the method of superimposing regular waves described above was selected to generate irregular waves. Initially, after selecting the desired irregular wave spectrum, the wave height can be calculated using the PSD (Power Spectral Density) for each frequency to develop a regular wave condition. In the case of a phase angle, an arbitrary value should be set. A phase condition that can exclude the condition of instantaneous wave height increase was selected by comparing the time series wave height data produced according to the phase.

$$S(f_i)=\frac{A_n^2}{2}·∆f$$

(15)

$$A_n=\sqrt{S\left(f_i\right)·2·∆f}$$

(16)

$$k_n=\frac{{\omega }_n^2}{g}$$

(17)

$${\lambda }_n=2·\frac{\pi }{k_n}$$

(18)

The wavelength and wave height conditions are required for the superposition of regular waves, which can be calculated using the frequency of wave spectrum and PSD value. The frequency of wave spectrum equation was calculated, as shown in Equation (15). The wave height of the corresponding frequency can be obtained using Equation (16). Under the assumption of deep-water conditions, the wave number can be calculated as Equation (17), and the corresponding wavelength frequency can be calculated in the spectrum as Equation (18). In addition, to superimpose regular waves, it is necessary to select the range and interval of the frequency to be used. In this study, the frequency range was used for the JONSWAP spectrum, corresponding to a significant wave height of $4.5\mathrm{m}$ and a peak wave period of $12\mathrm{s}$. A range of frequencies was used (0.05–0.21 Hz) in which a spectrum form is manifested. When the frequency interval of the spectrum is equally divided, the frequency is divided at the same interval because the wave repeats after a certain time. In addition, when repeated waves are used, the amplitude obtained for each frequency can be estimated accurately during FFT analysis.

FFT analysis will experience a spectral leakage phenomenon if the start and end of the data do not match. In general, FFT analysis was performed using a window function that minimizes the influence of non-continuous data to minimize the spectral leakage phenomenon. On the other hand, a spectral leakage phenomenon occurs even if a window function is applied to obtain the exact amplitude required in this study, making it difficult to find the exact amplitude for each frequency. Even if the same data is used, a difference in amplitude occurs if the start and end are not continuous, as shown in Figure 22a. FFT analysis was performed to solve this problem by matching the start and end points as much as possible. As a result, it showed a good match with the target wave spectrum, as shown in Figure 22b.

Irregular waves are generated by superimposing regular waves. By selecting overlapping regular waves at the same frequency interval, the irregular waves generated are repeated for a certain period of time. The reason why irregular waves use the property of repeating at regular intervals is to minimize the spectral leakage phenomenon by performing FFT according to the repeated peaks. It allows to obtain accurate resistance and amplitude of motion for each frequency. Therefore, appropriate selection of the NOF (Number of Frequencies), which is the number of divided wave spectra, is needed to generate irregular waves. When generating an irregular wave, the cycle of the repeated irregular wave varies depending on the NOF, as shown in

Table 8. Wave period, according to the frequency number.

Frequency Number	Wave Period (s)
70	85.36
100	122.47
130	159.59
160	196.70

. In the current calculation, waves ranging from 70 to 160 were tested to check the wave repetition period different from the NOF, as shown in Figure 23. Among them, more than 100 were considered suitable to show the characteristics of the spectrum for each frequency. The NOF = 100 case sufficiently represents the shape of the wave spectrum, and the repetition period was judged to be an appropriate time to perform numerical analysis.

Figure 24 shows the generated wave elevations, which can be calculated theoretically (approximately 12 min on the full scale). The irregular wave data used in this study was repeated for a certain period, and the repeated waves could sufficiently represent the characteristics of the spectrum. The number of waves encountered within that time is approximately 80. The added resistance experiments in irregular waves recommended at least 50 waves be used, and the reliability of the number of waves encountered in this calculation was confirmed [31].

Case studies were conducted on time intervals, x-direction grid, and z-direction grid to review the generated irregular waves.

Table 9. Irregular wave test conditions.

	Z1X1	Z1X1T2	Z2X1	Z2X2
$T_{Pe}/∆t$	876	438	876	876
${\lambda }_{Pe}/∆x$	64	64	64	32
$H_S/∆z$	80	80	160	160

lists the variables used, and the timesteps and grid sizes were set based on the significant wave height and peak wave periods presenting as irregular. The wave spectrum has one peak, regardless of the ship speed, but the position where the resistance peak appears varies because of the ship’s speed. To apply this, the encounter peak wave period $\left(T_{Pe}\right)$ was calculated and used to select the timestep.

The timestep corresponding to 1/876-times (0.002 s) and 1/438-times (0.004 s) of the encounter peak wave period was reviewed. When the Z1X1 and Z1X1T2 conditions were used, the power spectrum density according to the conditions was compared, as shown in Figure 25. Z1X1T2, which has a relatively large time interval, tends not to simulate the spectrum peak well. Therefore, the 1/876 times of the encounter peak wave period (0.002 s) was applied to the calculation. As a result of comparing the z-direction grid corresponding to 1/80-times and 1/160-times the significant wave height to select the grid in the z direction, there was no difference at the peak. Z1X1 was reasonable, considering the calculation time, because a lower timestep level would be required if a denser grid was used. The encounter wavelength length $\left({\lambda }_{Pe}\right)$ can be obtained in the same manner, as calculated using Equation (18), through the encounter peak wave period. When comparing the x-direction grid corresponding to 1/64-times and 1/32-times the wavelength, the overall power spectrum density tends to decrease near the peak, such as in the case of Z2X2 using a coarse grid, and a grid size corresponding to 1/64-times was selected.

6.3. Resistance Performance in an Irregular Wave Using CFD Simulation

Based on the irregular waves test above, the added resistance computation was performed in irregular waves for the bollard condition (0 knots, Fr = 0) and the design speed (14.5 knots, Fr = 0.172). The analysis was performed for approximately 135 s.

The performance in irregular waves for the design speed can be obtained as shown in Figure 26. In terms of resistance, the effects of short and long waves were mixed. The PSD of the resistance time series data was derived to confirm whether the obtained time series resistance data reflected the influence of regular waves for each frequency. Figure 27 shows the power spectral density of the resistance performance. In the case of design speed, the PSD of the resistance is distributed at the same frequency as the frequency of the encounter JONSWAP spectrum, considering the ship speed. Through this, the resistance performance by frequency of irregular waves is well shown. The bollard condition showed a similar shape to the JONSWAP spectrum without the effect of the ship speed, and the PSD value of the resistance was relatively low.

The average value of the added resistance time series data in irregular waves was derived as shown in

Table 10. Comparison of the resistance in irregular waves.

Condition	Total Resistance [N]	${\mathit{R}}_{\mathit{T},\mathit{I}\mathit{r}\mathit{r}\mathit{e}\mathit{g}\mathit{u}\mathit{l}\mathit{a}\mathit{r},\mathit{i}}/{\mathit{R}}_{\mathit{T},\mathit{C}\mathit{a}\mathit{l}\mathit{m}}$
Calm water, $R_{T,Calm}$	78.68	-
Spectral method (Indirect method), $R_{T,Irregular,SM}$	154.43	1.96
Direct computation (Direct method), $R_{T,Irregular,DC}$	134.98	1.71

, and the results were compared with the result of the spectral method. Compared to the resistance in still water, the spectral method showed an increase in resistance of about 196% and direct random wave calculation of about 171%. According to the two estimation methods, the total resistance to perform CFD analysis directly was approximately 11% lower than the spectral method. In the case of the spectral method, the resistance was estimated by assuming that the resistance occurs linearly with wave height. The resistance was low because of the nonlinearity caused by the superposition of the wave that appears in the direct irregular wave calculation. This increase of resistance is a result of the proposed wave condition, and sufficient review through other irregular wave conditions and experiments is required.

Time series motion data were analyzed to check the motion performance in irregular waves, as shown in Figure 26b,c. The time series motion data were derived from the amplitude for each frequency using FFT analysis. The motion response coefficient was derived using the obtained amplitude for each frequency and the wave height data of the corresponding frequency, as shown in Figure 28. Except for the frequency near $f_e{/f}_{n,Heave}=1$, near the natural frequency, the motion response was similar to the motion RAOs resulting in the regular wave. In particular, the motion response was almost zero in the high-frequency section representing the short-wave region. The reason why the motion RAOs appeared large at the resonance frequency was because the wave steepness used for the irregular wave was small. Corresponding trends were also observed in the KVLCC2 of Lee et al. [20].

7. Conclusions

This study estimated the added resistance and motion RAOs of K-Supramax at the design speed (14.5 knots, Fr = 0.172) in various wave conditions using CFD. The reliability of the calculation was reviewed through the similarity of the added resistance coefficient and motion RAOs compared to the experimental results. At first, the difference in resistance performance according to the degree of freedom (2DOF and 6DOF) of the vessel in head sea condition was compared, and the added resistance performance during bow quartering sea was evaluated in 6DOF condition. Secondly, based on the results of the regular wave of the head sea, the added resistance performance in irregular waves was estimated using the spectral method. In addition, the irregular wave was generated through the superposition of regular waves and an added resistance computation was performed on the irregular waves directly.

In this study, heave and pitch were given free for 2DOF, and, in the case of 6DOF, a spring was used to maintain the ship’s heading. In the head sea condition, the difference between the results of 2DOF and 6DOF was insignificant. Also, there was no significant difference in the amplitude and phase of the motion of vessel. Surge vibration by the spring used to implement the 6DOF condition was observed, but it did not affect the value of the added resistance. On the other hand, the effect on ship performance in waves was observed through bow and stern relative motion calculations. The added resistance tendency is similar to the bow relative motion, and the change of heave and pitch RAOs was similar to the stern relative motion, which depends on the wavelength ratio.

Added resistance computations were performed in bow quartering sea conditions (10° intervals up to 150°), as well as head sea conditions, under 6DOF conditions. The added resistance showed the highest value near the resonance frequency of 170°, and then, as the heading angle increased, the added resistance tended to decrease, regardless of the wavelength ratio. Heave RAOs tended to shift the peak value to the short wavelength region due to the natural frequency change according to the heading angle. Pitch RAOs were all similar, except that they decreased as the heading angle increased in long waves. Roll RAOs were insufficient in bow quartering sea conditions up to 150°.

The added resistance in irregular waves was investigated using the spectral method, based on the theory that regular waves can superimpose linearly. As the peak wave period increased, the peak of the JONSWAP spectrum shifted to a lower frequency. When the peak wave period was 12 s, the added resistance was the highest in irregular waves, because the frequency of the peak of the transfer function is similar to the frequency of the peak of the added resistance in regular waves.

By generating irregular waves through CFD, the added resistance was directly estimated in irregular waves. Irregular waves were created through superimposing regular waves, and the optimal parameters were used by studying variables such as timestep, grid, and number of frequencies in the wave spectrum. As the computation results of direct irregular waves, the frequency of the encountered JONSWAP spectrum and the frequency of the power spectrum density of resistance were distributed equally. The resistance in irregular waves increased by approximately 92% for the spectral method and by approximately 72% for the direct irregular wave calculation compared to the resistance in calm water. In the case of the spectral method, the resistance was estimated by assuming that the resistance occurs linearly with the wave height. However, the added resistance among the direct irregular waves is thought to be estimated low because the nonlinearity of the incoming irregular waves is reflected, and additional research is needed on this.

In this study, 6DOF should be considered to investigate the resistance performance in the bow quartering sea. So, it was confirmed that there was no significant difference in resistance performance between 2DOF and 6DOF in head regular waves using a soft spring system, which is a 6DOF analysis method. The resistance performance in the bow quartering sea was investigated. Also, the optimal number of overlapping waves was selected to reflect the energy of the waves constituting the wave spectrum as much as possible to estimate the resistance performance in irregular waves, and economical calculations were performed by minimizing the calculation time through the selected NOF. As a result, it is judged that the resistance spectrum and the RAOs were reasonably obtained.

Future research plans to study not only the resistance but also the self-propulsion performance in waves and the effect of wind to consider more realistic seaway conditions.