*Article* **On the Comparative Seakeeping Analysis of the Full Scale KCS by Several Hydrodynamic Approaches**

#### **Florin Pacuraru, Leonard Domnisoru and Sandita Pacuraru \***

Department of Naval Architecture, "Dunarea de Jos" University of Galati, 800008 Galati, Romania; florin.pacuraru@ugal.ro (F.P.); leonard.domnisoru@ugal.ro (L.D.)

**\*** Correspondence: sorina.pacuraru@ugal.ro

Received: 17 October 2020; Accepted: 19 November 2020; Published: 25 November 2020

**Abstract:** The main transport channel of the global economy is represented by shipping. Engineers and hull designers are more preoccupied in ensuring fleet safety, the proper operation of the ships, and, more recently, compliance with International Maritime Organization (IMO) regulatory incentives. Considerable efforts have been devoted to in-depth understanding of the hydrodynamics mechanism and prediction of ship behavior in waves. Prediction of seakeeping performances with a certain degree of accuracy is a demanding task for naval architects and researchers. In this paper, a fully numerical approach of the seakeeping performance of a KRISO (Korea Research Institute of Ships and Ocean Engineering, Daejeon, South Korea) container ship (KCS) container vessel is presented. Several hydrodynamic methods have been employed in order to obtain accurate results of ship hydrodynamic response in regular waves. First, an in-house code DYN (Dynamic Ship Analysis, "Dunarea de Jos" University of Galati, Romania), based on linear strip theory (ST) was used. Then, a 3D fully nonlinear time-domain Boundary Element Method (BEM) was implemented, using the commercial code SHIPFLOW (FLOWTECH International AB, Gothenburg, Sweden). Finally, the commercial software NUMECA (NUMECA International, Brussels, Belgium) was used in order to solve the incompressible unsteady Reynolds-averaged Navier–Stokes equation (RANSE) flow at ship motions in head waves. The results obtained using these methods are represented and discussed, in order to establish a methodology for estimating the ship response in regular waves with accurate results and the sensitivity of hydrodynamical models.

**Keywords:** strip theory; BEM; RANS; regular wave; seakeeping

#### **1. Introduction**

Shipping, often considered as the main transport channel of the global economy, is responsible for approximately 80 percent of world trade. An increasing focus on more environmentally friendly shipping pushes hull designers to further ensure fleet safety, proper operation of the ships, and, more recently, compliance with International Maritime Organization (IMO) regulatory incentives regarding the energy efficiency operational indicator (EEOI) and emissions reduction. Considerable efforts have been devoted to an in-depth understanding of the hydrodynamics mechanism and the prediction of ship behavior in waves. In addition to safety, efficiency and operability, the waves' loads may determine structural failure. Prediction of seakeeping performances with a certain degree of accuracy is a demanding task for naval architects and of great practical interest for shipbuilders, owners, operators, as it affects both the ships' design and operation [1]. More recently, the shipping industry has embraced more and more digitalization. Digital twin concept couples physical–numerical modelling and simulation to evolve solutions that improve vessel predictability, behavior control and response.

There are several aspects concerning seakeeping that make it one of the most challenging problems in ship hydrodynamics. The nonlinearities induced by the fluid viscosity and hydrodynamic pressure, boundary conditions at free-surface formulation, interference between external waves and ship's body, specific geometry of the hull shape, all require advanced iterative time-domain procedures for the ship's oscillations response analysis, involving significant computational resources [2].

The problem of a moving body interacting with waves has been pursued since the very first Froude and Michel studies. Since then, different approaches, such as experimental fluid dynamics, potential flow and recently computational fluid dynamics, have been developed to estimate the seakeeping performances. Traditionally, the hydrodynamic performances of a full-scale ship are determined by extrapolating the results of model-scale towing tank tests. Another approach commonly adopted in ship hydrodynamics is the employment of numerical techniques to solve the equations' system of ship motions. The prediction of non-linear phenomena specific to ship motions in waves is difficult to be solved numerically, not only due to the complexity of the problem, but also the results are highly dependent on the details of the hull form and the incident wave condition. Due to the development of numerical modelling methods and computing power increase, direct prediction of full-scale ship performance became a practical approach. Most of the available techniques used to predict ship motions rely on assumptions from potential flow theory. Even if the natural trend in ship hydrodynamics is to move from frequency-domain to time-domain approach, from linear 2D strip theory type to fully 3D nonlinear techniques, and from potential to viscous computation, potential flow methods are still highly used in ship design, since they provide robust and quite accurate results in low to moderate sea states [3].

The frequency-domain approach is carried out mainly for linear or weakly nonlinear wave theories, where time dependence can be removed by assuming that solution is harmonic in time, in consequence only steady solutions are solved. Assuming that the ship body is slender, strip theory simplifies the 3D flow problem into a 2D formulation, modelling the ship hull as a set of multiple 2D ship stations [4]. The independent boundary value problem for each station can be solved analytically (Lewis form method) or by boundary element method. For the 2D hydrodynamic formulation, since the 1950s, various strip theories have been developed for the seakeeping problem [5]. The pioneer work of Korvin-Kroukovsky [6] set the principal feature of the strip theory for calculating ship motions based on the slender body assumption. This was the first suitable theory for numerical computations of ship motions that had adequate accuracy for engineering applications [7]. A modified strip theory approach of Gerritsma and Beukelman [8] was shown to obtain a good agreement with experimental tests for head wave case. Ogilvie and Tuck [9] developed a mathematically consistent approach based on short wave-length approximation, conducting a systematic analysis for the slender body problem to determine the added mass and damping for heave and pitch motions. Most of today's strip methods are variations of the approach proposed by Salvesen et al. [7], which is one of the most complete versions, solving five degrees of freedom in ship's motion equations, as the surge component was neglected, being not relevant for a slender body oscillation. Since then, more comprehensive strip theories have been developed for ship design, such as [10–14]. Few comprehensive reviews of strip theory variations have been reported by [4,15,16]. Most recently, a combination of two-dimensional strip theory with the two-dimensional Green function based on the potential theory to solve boundary values and motion responses of a semi-planing craft have been reported by [17]. Despite theoretical shortcomings, strip theory has the advantage of being fast, cheap and sufficiently accurate for a range of hull forms and moderate speeds. However, for higher speed vessels, highly flared hull forms, wave loads or extreme motions, other approaches have to be used.

The effort devoted to overcoming the weakness of the strip theory methods and to modelling the non-linear phenomena leads to the development of other potential flow alternatives based on numerical 3D methods that enhance the modelling technique of the physical domain of interest. By the late 1970s, a new treatment of the boundary conditions led to the Neumann–Kelvin approach [2], which supposes that the body boundary condition is applied to the mean position of the exact body surface and linearized free-surface boundary condition [2]. The most efficient way to solve the Neumann–Kelvin problem is to employ boundary integral methods where the solution is formulated in

terms of singularities, as sources or dipoles, over the hull and free surface [2]. Two different categories of methods have been developed based on the type of singularities used in integral equations, one relies on Green function [18], which automatically satisfies both the radiation and the linearized free surface condition, and the other one relies on the distribution of the simple Rankine sources [18], not only on the hull surface, but also on the free surface. Some of the weaknesses of strip theory and Green function approach can be overcome by using a 3D time-domain panel method based on Rankine sources, including the fully 3D effects of the flow and forward speeds effects.

Hess and Smith [19] proposed in 1964 a method for evaluating velocity and pressure fields around a fully immersed arbitrary three-dimensional body. The method, also called the panel method or the boundary element method allowed, for the first time, to calculate the flow around an arbitrary three-dimensional body imposing the boundary condition exactly on the body surface. Firstly Gadd [20], then Dawson [21], applied the panel method for steady free surface flow distributing Rankine sources on hull and undisturbed free surface. A wide variety of different approaches for solving the nonlinear steady forward motion problem, based on Rankine sources distribution, has been reported by [22–25]. Due to the development of computing capabilities, the panel method has been continuously improved to solve the problems of ship motions in waves considering different modelling of nonlinearities of the hull and free surface boundary conditions. Nakos et al. [26] used a Rankine source method to solve transient wave–body interactions. This method has linearized the solutions about the double body flow. Large ship motion analysis performed by [27] was based on the application of the desingularised source method. Söding et al. [28] applied patch method instead of the collocation method to satisfy boundary conditions on the solid body surface, considering for the free surface a nonlinear condition. Recently, an improved desingularised Rankine method has been proposed by Mei et al. [29] to solve 3D diffraction and radiation problems. Dai and Wu [30] combined a method based on a Rankine panel method for the near-field and a transient Green function method for the far-field to predict large amplitude ship motions.

Potential flow solvers for seakeeping prediction have been continuously improved in terms of efficiency and accuracy, but still, the effect of viscosity, turbulence, wave dispersion, wave breaking, green water and slamming impact loading, deck green sea cannot be captured properly using the method based on potential flow assumption [31,32]. Considering the steadily increasing computational power and parallel computing, computational fluid dynamics methods based on the Reynolds-averaged Navier–Stokes (RANS) approach is more frequently applied to solve unsteady seakeeping problems. Most of the commercial solvers adopt a combination of unsteady RANS method, based on finite volume method and multi-phase flow approach, using the volume of fluid method for the treatment of free surface, which is rapidly gaining popularity for ship's motions applications. The advantages of the nonlinear computation techniques without using analytical formulas for added resistance or empirical values for viscous effect are significant for the seakeeping analysis accuracy, but they are considerably more time consuming than potential flow approaches. Therefore, a good balance between accuracy and computational speed is required, especially for the ship design process.

The majority of RANS seakeeping simulations have been performed at model scale. The first attempt to solve ship motions in waves was presented in [33]. The results show some problems regarding the accuracy of free-surface due to the limited grid quality. Later, Simonsen et al. [34] carried out motions analysis in heave and pitch motions in regular head waves for KRISO container ship (KCS) hull using the CFDSHIP-IOWA code. Recently, Lungu and Bekhit [35,36] have assessed the seakeeping performances of the KCS and KVLCC (KRISO Very Large Crude Carrier) ship models for regular head wave conditions using NUMECA/FineMarine code. Their results show good agreement with the experiment for model scale. On the other hand, the scale effect seems to be important, as Hochkirch and Mallol [37] presented significant scale effects due to the differences between model-scale flows and full-scale flows. Tezdogan et al. [38] investigate seakeeping behavior and performance of the KCS model at a slow forward speed using Star-CCM+ package (Siemens PLM Software, Plano, TX, USA).

Considering the context described above, the present work is focused on a systematic comparison study concerning the analysis of ship's pitch, heave and roll motions in regular waves, for a full-scale ship model. The study was made using three different numerical approaches: the linear strip theory based on the Lewis form shape parameterization developed by Domnisoru [39], as in-house code DYN (module OSC - Oscillations), the non-linear potential flow based on Rankine source distribution using the commercial code SHIPFLOW and, for a limited number of cases, due to it being a highly time consuming solution, the viscous RANS method using NUMECA/FineMarine commercial code. Several series of computations have been performed for the regular wave height of 2 m, heading angle in the range of 0 to 180 deg, with a step of 45 deg. For the BEM model also, a wave height of 8 m has been considered. The computations have been carried out for the full-scale KCS hull numerical model, considering two forward speeds, 12 and, respectively, 24 knots.

The authors of this article aim to integrate all the presented methods in the same study in order to achieve a methodology that serves the rapid decision-making process in the case of further studies concerning the choice of the most accurate estimation method for ship motions on waves, depending on the specific requirements of the study and taking into account costs and calculation time. The proposed methodology meant approaching the mentioned methods, globally, with in-house code and commercial software, and the major advantage consisted of a better management of the input data, working conditions, calculation hypotheses, analyzed cases. It is worth mentioning that in the present article, full-scaled ship computations are considered, while most of the literature references reported analyses on model-scaled ship.

#### **2. Methods for Ship Motions Prediction**

As most of the available techniques to predict ship motions for design purposes rely on assumptions of potential flow theory, two different methods from this category have been chosen for comparison. The third approach considered for the present study is based on RANS solver, but, despite the parallel computing, this method does not prove yet to be able to solve the seakeeping problem in reasonable time to be consider in the practical design procedures. Taking into account that the methods involved in the present study rely on well-known mathematical formulation presented in seakeeping references [4,7,22,40–42], in this chapter only a brief presentation for each method is included.

The first method is based on a potential flow hydrodynamic linear 2D strip theory, with Lewis ship's stations parameterization and regular wave Airy model excitation, with a frequency-domain solution of the coupled heave–pitch motion equations and uncoupled roll motion equation. This method is implemented by Domnisoru [39]—as in-house code DYN, module OSC—for the linear oscillations response amplitude operator (RAO) response amplitude operator computation. The mathematical formulation of the method used is described in [39]. The DYN code, module OSC, has been validated by experimental tests on scaled models of a fishing vessel [43] and a survey vessel [44], on follow, head, beam and quartering regular wave conditions.

For the second set of computations, a 3D fully nonlinear time-domain boundary element method, implemented into the commercial code SHIPFLOW [22,40] has been employed for ship motion estimation, considering only potential flow formulation being suitable for practical use. The method solves the missed boundary value problem for Laplace equation by distributing sources on the hull and free surface. Integrating the kinematic and dynamic boundary conditions with a fourth order Adam–Bashford–Moulton method, the elevation of free surface is obtained at each time step based on a Mixed Euler–Lagrange method. A blending zone based on an analytical solution is introduced in order to avoid the reflection from domain boundaries and in order to generate a certain incident wave. Once the velocity potential is obtained on the hull, using Unsteady Bernoulli formulation, the pressure distribution on the hull can be obtained. Finally, after integrating pressure on the hull, forces, moment and then motions are determined. A full description of the mathematical model and the numerical method may be found in [40,45]. Validation of the above-described method has been

reported by Larsson et al. [46] for added resistance, pitch and heave, and by Coslovich et al. [47] mainly for roll.

The commercial software NUMECA/FineMarine [41,42] has been used in this study to solve the incompressible unsteady RANSE flow at ship motions in head waves. The solver relies on the finite volume method to build the spatial discretization for the governing equation. Closure to the turbulence is achieved by making use of the k-ω SST model with wall function formulation [41]. Pressure–velocity coupling is enforced through a SIMPLE (Semi-Implicit Method for Pressure Linked Equations)-like approach, where the velocity updates come from the momentum equation and the pressure is extracted from the mass conservation constraints transformed into pressure equation. Convection and diffusion terms in the RANSE are discretized using a second-order up-wind scheme and a central difference scheme, respectively. Free-surface capturing strategy is based on multi-phase flow approach using a volume of fluid method with high-resolution interface schemes [42]. Ship's 6 DOF (Degrees Of Freedom) can be solved by the solver, but also some degree of freedom can be restrained. Mathematical model and the numerical framework are nicely described in [41]. Several validation studies of the solver have been reported by [48] for ship resistance and self-propulsion, and by [31,35,36,49] for ship motions in head waves. Their results show good agreement with the experiment for the scaled model.

In the preliminary ship design process, time is an important aspect to be considered due to the fact that naval architect usually cannot afford too long computational time. Considering the present methods comparison for one case (*v* = 12 Kn, ω = 0.4, *hw* = 2 m and head wave), by far strip theory is the less time-consuming method, the time spend for a single calculation case was 3 min using a single processor. The nonlinear BEM computations have been performed on 4 processor machines and the physical computational time for the same case was 19 h and 5 min for a number of about 68,000 panels. The RANS computations have been carried out on a high-performance computing (HPC) machine using 120 processors for a grid of about 20 million cells, which led to a 41 h and 47 min computational time. On the other hand, time consuming should be also correlated with the accuracy and with the amount of information received from a specific flow post-processing. The RAO functions for heave, pitch and roll based on the linear strip theory implemented in DYN code have been validated in references [43,44] for heading angles, 0, 90, 180 degrees, proving good results for practical investigations, revealing average differences from 12% to 18%. A validation BEM method in terms of heave and pitch was performed by the authors in [50] and the computed results showed good agreement of 3% to 10% with towing tank measurements. As mentioned before, extensive validation studies may be found in [45,46]. In addition, grid convergence and validations studies based on FineMarine RANS solver have been reported by Lungu and Bekhit [35,36]. Their results revealed good agreement with the towing tank test measurements for model scale, 1–7%. Taking all the above mentioned into consideration, a trade-off between accuracy and computational costs must be counted in order to obtain feasible results for the ship design practical use.

In order to compare directly the computational fluid dynamics (CFD) results with those obtained by using linear strip theory, for all numerical analyses, the regular Airy model is considered as excitation source. Both CFD methods are hydrodynamic nonlinear, so that even if the excitation has one harmonic, the time-domain response has several harmonics. For the two CFD methods, heave, pitch and roll time records are spectrally analyzed by using Direct Fourier Transformation Method. Regarding the response amplitude spectra (see Section 6), only the motions amplitudes with encountering wave frequency harmonic are considered. The response amplitude operators (RAO s) for ship motions obtained by the two CFD methods represent the ratio between heave, pitch and roll amplitudes on wave's harmonic and the regular wave amplitude.

#### **3. Hull Geometry and Conditions**

The objective of the present study is to investigate the seakeeping performances of the KRISO container ship (KCS) while sailing in regular wave at two speeds, 12 and 24 Kn, corresponding to Froude numbers of 0.13, and 0.26, respectively. The KCS hull is a benchmark test case for ship hydrodynamics widely used in the marine hydrodynamic scientific community [34]. The hull is a modern commercial container vessel with bulbous bow and flare fore above the waterline. The stern is a typical pram type stern with transom. Bare hull, without rudder, has been considered for the seakeeping calculations. KCS hull geometry is presented in Figure 1 and the main dimensions and conditions used for numerical simulations are presented in Table 1.

**Figure 1.** KRISO container ship (KCS) hull geometry.



In Table 2, the test cases for seakeeping analysis are presented. The ship motions are approached considering the ship on regular Airy wave, having the wave height of 2 and 8 m, the wave frequency between 0 and 2 rad/s and heading angle in range of 0 to 180 deg, with a step of 45 deg. Even though the ship response in a regular wave at zero speed was computed, in the following paragraphs the ship motions for only two speed values of 12 and 24 Kn, respectively, are presented. Although, the study was carried out for three values of vertical center of gravity, *KG* = 7.3, 10.8 and 14.3 m, in the paper, the results for only *KG* = 10.8 m, which corresponds to a medium loading case of the ship, are presented.


**Table 2.** Test cases.

#### **4. Numerical Setup**

Several details about the computational conditions will be briefly described in the following. For the discretization of the hull used for the linear 2D strip theory calculations, a model of 526 stations has been used.

For the BEM method, the boundary value problem modelling the field equation and boundary conditions have been solved, making use of the commercial code SHIPFLOW. The computations carried out with BEM is a fully nonlinear unsteady three-dimensional potential flow method, that takes into account the geometric non-linearities of ship shape and the induced hydrodynamic non-linearities, although the external excitation is kept linear. Wave reflection from domain's boundaries is a pure numerical problem that affects all the simulations with a restricted domain. In order to avoid wave reflection from the free surface truncation boundaries, the flow on the intersection between the outer and inner domain has to be the same. To achieve this condition, a numerical damping zone is located close to the domain's boundaries and such a zone only has to dampen the difference between the velocity potential in the inner and outer domain. This difference in the solution between the inner and the outer domain is mainly due to the presence of the body with radiation and diffraction effects but can also derive from a numerical error. This damping is achieved adding a term in the free surface boundary condition for all those panels that belong to such a zone [51,52]. The fluid domain assumes that the flow is known a priori in the outer domain and the current method uses analytically described waves to represent the flow field in the outer domain. The hull and free surface have been discretized by panels. The calculations have been performed for a surface domain of one ship length upstream, two ship length downstream, the width of the free surface being two ship length. For BEM method, the ship hull is discretized both sides without the request of symmetry condition, so that any heading angle condition can be computed. The number of panels distributed on hull and free surface varied between 60,000 and 90,000 panels function of the wave frequency case studied, which assures a minimum number of 40 panels per wavelength. The computations have been performed on local desktop machines with four cores of 3.1 GHz.

For the unsteady RANS simulation of incompressible flow of ship motions in head waves, the commercial solver NUMECA/FineMarine was used. A reference length has been considered *Lref* = max (*LBP*, λ). The dimensions of the computation domain have been chosen according to the ship and wavelength, considering 2.0 *Lref* upstream, 4.0 *Lref* downstream, 2.0 *LBP* on the side, 4.0 *LBP* underneath, and 2.0 *LBP* above the undisturbed free-surface level. The boundary conditions imposed on the solid wall and domain boundaries, but also the dimensions considered for the computational domain, are depicted in Figure 2. A mirror condition is applied on the centerline plane of the ship and on the side boundary of the domain, to avoid wave reflections at lateral boundary. Considering the head wave case computed by the RANS method, the corresponding wave is generated at the upstream boundary and the inlet boundary is chosen to ensure at least two full waves to be generated before encountering the ship. The free surface is refined for the entire domain using three different refinement boxes, as can be observed in Figure 2. The first starts from the inlet boundary and is extended until it reaches 1.0 *Lref* behind the ship. The cell size is chosen for this refinement zone in *x*- and *y*-direction to

provide minimum 60 cells per wavelength, such that Δ*x*1 = Δ*y*1 = λ/60, while in the *z*-direction the refinement criterion is selected as Δ*z* = *hw*/16, where *hw* represents the wave height. The refinement depth in *z*-direction is extended for 3.0 *hw* equally distributed above and beneath the non-disturbed free-surface level, which is set at the design draft *T* = 10.8 m from the ship base line. The second and third boxes are coarsened gradually in *x*- and *y*-directions to generate sufficient numerical damping zones to prevent any reflections from the exit boundary. The second box is extended for 1.0 *Lref* and it is used as an initial damping relaxation zone, in which the cell size is coarsened in *x*- and *y*-directions by a factor of 4 (i.e., Δ*x*2 = Δ*y*2 = 4.0 × Δ*x*1). The third refinement zone is used as a final damping zone, having the cell size coarsened in *x*- and *y*-directions by a factor of 8.0 (i.e., Δ*x*3 = Δ*y*3 = 8.0 Δ*x*1), while the refinement in the *z*-direction for the entire domain is maintained unchanged. The effect of damping zone on the RANS solution is clearly seen in Figure 3.

**Figure 2.** RANS Computational domain, dimensions and boundary conditions.

**Figure 3.** RANS damping zone effect.

Body-fitted full hexahedral unstructured meshes of about 22 million cells have been generated **(**Figure 4). The time step Δ*t* has been chosen in order to fulfil the condition of 300-time steps per wave period, with fourth order convergence criteria. All the computations have been performed on a HPC machine with 120 cores of 3.3 GHz.

**Figure 4.** Computational grid. (**a**) Domain discretization; (**b**) mesh details in bow area.

In general, methodology used in the present paper for the seakeeping performance prediction is based on the practical guidelines recommended by the International Towing Tank Conference (ITTC), by software providers and different sensitivity and validation studies [35,36,49–52], but for some particular cases the recommendations could not be accomplished due to the practical reasons.

#### **5. Grid Convergence Test**

A grid study for SHIPFLOW calculations was conducted for three grids based on the Richardson extrapolation [53] in order to investigate the numerical simulation error and uncertainty. Two calculation cases have been considered for the grid test: ship on regular waves, head waves (μ = 180 deg) and beam sea (μ = 180 deg). For both cases ship speed was *v* = 12 Kn. The results of grid convergence study for RAO heave, pitch and roll are summarized in Table 3. The convergence ratio is defined as *RG* = ε21/ε32, where ε<sup>21</sup> = *S2* − *S1* and ε<sup>32</sup> = *S3* − *S2* stands for simulation error, *S1* represents the solution calculated for fine grid (68,716 panels), *S2* for medium grid (32,692 panels) and *S3* for coarse grid (16,346 panels), according to Richardson extrapolation approach [53], with the expressions (1) ÷ (4). The grid corresponding to each level of refinement is depicted in Figure 5 for hull and free surface.


**Table 3.** Grid convergence test results.

**Figure 5.** Grid refinement levels used for the grid convergence test. (**a1**) Coarse grid on hull surface; (**a2**) Coarse grid on free-surface; (**b1**) Medium grid on hull surface; (**b2**) Medium grid on free-surface; (**c1**) Fine grid on hull surface; (**c2**) Fine grid on free-surface.

Considering the computed values of *RG* for all variables of grid convergence study, one can observe that monotonic convergence is achieved for all the simulation variable. The order of accuracy *pG* can be calculated based on the refinement ration *rG*, as follows:

$$p\_G = \frac{\ln\left(\frac{\varepsilon\_{32}}{\varepsilon\_{21}}\right)}{\ln(r\_G)}.\tag{1}$$

The calculated order of accuracy is used further to determine the correction factor *CG* and the error δ*G*, as follows:

$$c\_G = \frac{r\_G^{pc} - 1}{r\_G^{p\_{th}} - 1},\tag{2}$$

$$
\delta\_G = \frac{\pounds\_{21}}{r\_G^{p\_G}} \tag{3}
$$

where *pth* is the theoretical grid convergence order usually considered *pth* = 2. Finally, the uncertainty is computed by the expression:

$$\mathcal{L}I\_{\mathcal{G}} = |\mathcal{C}\_{\mathcal{G}} \times \delta\_{\mathcal{G}}| + |(1 - \mathcal{C}\_{\mathcal{G}}) \times \delta\_{\mathcal{G}}| \tag{4}$$

Based on this study, bearing in mind that the percentage estimated error with respect to the fine grid is very small, i.e., ε21%*S1* << 1, which concludes that the solution computed by SHIPFLOW is grid independent.

#### **6. Results and Discussion**

In this section, the ship hydrodynamic response in regular waves is presented and discussed. The results are represented in the frequency domain, in terms of response amplitude operators (RAO) for heave, pitch and roll motion, respectively. The charts contain the results obtained by linear strip theory method (DYN code) and the boundary element method (SHIPFLOW software), also including a few results obtained using the RANS method (NUMECA/FineMarine).

In Figure 6, the ship on regular waves, for the case of top speed *v* = 12 Kn, μ = 90 deg, ω = 0.4 rad/s and *hw* = 2 m, is presented. The images represent the hull position at five different fractions of ship motion periods (Δ*T*) related to free surface topologies. Hydrodynamic forces and moments are determined by integrating pressure over the ship hull based on the quadratic pressure distribution. Pictures on the left, Figure 6a1–e1, present the pressure distribution over the ship hull. One may see that the flow solution for ω = 0.4 rad/s reveals that the wave is fully developed, and artificial damping condition imposed on the boundaries works fine. Moreover, the interference of the wave system generated by the ship and the incident wave is well captured by the boundary element method and some nonlinear effects are expected (see the pictures on the right, Figure 6a2–e2).

**Figure 6.** *Cont.*

**Figure 6.** Pressure distribution and free surface topology computed for the case μ = 90 deg. (**a1**) Pressure distribution on the hull, ΔT = 0; (**a2**) Free-surface topology, ΔT = 0; (**b1**) Pressure distribution on the hull, ΔT = 0.125; (**b2**) Free-surface topology, ΔT = 0.125; (**c1**) Pressure distribution on the hull, ΔT = 0.250; (**c2**) Free-surface topology, ΔT = 0.250; (**d1**) Pressure distribution on the hull, ΔT = 0.375; (**d2**) Fre-surface topology, ΔT = 0.375; (**e1**) Pressure distribution on the hull, ΔT = 0.500; (**e2**) Free-surface topology, ΔT = 0.500.

In Figures 7–11, the response amplitude operators for heave, pitch and roll motions of KCS ship advancing in regular waves are represented. Heave RAO (m/m), Pitch RAO (rad/m) and Roll RAO (rad/m), respectively, for the speed case of ship speed 12 Kn and the wave height 2 m (amplitude *aw* = 1 m). In Figure 7, the heave response amplitude operator, Heave RAO (m/m), for the speed case of 12 Kn, obtained using a 2D strip and BEM methods is represented. As a general remark, for both methods, the ship response is significant in the lower frequency domain, while, as the frequency increases, the ship response decreases considerably. The maximum values are recorded for the heading angle μ = 90 deg, at ω = 0.756 rad/s. Once the circular frequency of the regular wave increases, one may observe that both methods are in good agreement.

**Figure 7.** Heave response amplitude operator (RAO) (m/m) for *v* = 12 Kn, *hw* = 2 m. *(***a**) μ = 0 deg; (**b**) μ = 45 deg; (**c**) μ = 90 deg; (**d**) μ = 135 deg; (**e**) μ = 180 deg.

**Figure 8.** Pitch RAO (rad/m), for *v* = 12 Kn, *hw* = 2 m. *(***a**) μ = 0 deg; (**b**) μ = 45 deg; (**c**) μ = 90 deg; (**d**) μ = 135 deg; (**e**) μ = 180 deg.

**Figure 9.** Heave RAO (m/m), DYN–SHIPFLOW–NUMECA comparison, *v* = 12 Kn, μ = 180 deg, *hw* = 2 m.

**Figure 10.** Pitch RAO (rad/m), DYN–SHIPFLOW–NUMECA comparison, *v* = 12 Kn, μ = 180 deg, *hw* = 2 m.

**Figure 11.** Roll RAO (rad/m), *v* = 12 Kn, *KG* = 10.8 m, *hw* = 2 m. (**a**) μ = 45 deg; (**b**) μ = 90 deg; (**c**) μ = 135 deg.

In Figure 8, the response amplitude operator for pitch motion, Pitch RAO (rad/m) is depicted, for ship speed 12 Kn, heading angles from 0 to 180 deg, and wave height of 2 m. For pitch motion at beam sea condition, the response amplitude is the smallest one, recording scattered values between the 2D strip and the BEM hydrodynamic models. In particular, at 90 deg, the variation of BEM results has two peaks, the second one at similar frequency as the DYN results; however, the values obtained by the boundary element method are higher than the strip model results at beam sea. The maximum value is recorded at μ = 135 deg for ω = 0.57 rad/s. A good agreement between both hydrodynamic models is found at a heading angle of 45 and 135 deg.

In Figures 9 and 10, the following are presented: numerical results for Heave RAO (m/m) and Pitch RAO (rad/m) obtained with linear 2D strip method (DYN), BEM boundary element method (SHIPFLOW) and RANS method (NUMECA/FineMarine), considering the ship on head regular wave, μ = 180 deg, with 2 m height and a ship speed of 12 Kn.

In Figure 11, the following are given: the results for the roll motion in terms of response amplitude operator, Roll RAO (rad/m). The results are represented for the heading angle of 45, 90 and 135 deg, respectively. The maximum value for Roll RAO is recorded for μ = 90 deg and, as long as it moves away from the beam waves, the roll motion of the ship attenuates.

Next, a comparison between the results computed with linear 2D strip and BEM hydrodynamic models, for a regular wave height *hw* = 8m(*aw* = 4 m) for BEM model, is depicted in Figures 12 and 13. First, one may see the Heave RAO (m/m) values at five different heading angles (Figure 12). The trend

of the differences between the two methods remains the same as for the case of wave height *hw* = 2 m. Very reduced differences occur in the case of the BEM method results at 2 and 8 m wave height, due to the reduced hydrodynamic nonlinearities. Higher differences appear at beam wave condition between 0.4 and 0.6 rad/s.

**Figure 12.** Heave RAO (m/m), *v* = 12 Kn, *hw* = 8 m for BEM model. (**a**) μ = 0 deg; (**b**) μ = 45 deg; (**c**) μ = 90 deg; (**d**) μ = 135 deg; (**e**) μ = 180 deg.

**Figure 13.** Pitch RAO (rad/m), *v* = 12 Kn, *hw* = 8 m for BEM model. (**a**) μ = 0 deg; (**b**) μ = 45 deg; (**c**) μ = 90 deg; (**d**) μ = 135 deg; (**e**) μ = 180 deg.

In Figure 13, the variation curves of Pitch RAO (rad/m) for wave height *hw* = 8 m for BEM model, with the peak value for the heading angle of μ = 135 deg around ω = 0.5 rad/s, are depicted. At beam waves, μ = 90 deg, differences between the hydrodynamic approaches are scattered, where the values of Pitch RAO are very reduced.

The next results are computed for the second ship speed case, *v* = 24 Kn, for both wave heights of 2 and 8 m, respectively. There were also considered the same medium loading conditions of the ship. In Figure 14, the RAO heave (m/m) is presented. The differences between the linear 2D strip and BEM hydrodynamic approaches have slightly increased compared to 12 Kn ship speed, especially for wave circular frequency lower than ω = 0.4 rad/s.

**Figure 14.** Heave RAO (m/m), *v* = 24 Kn, *hw* = 2 m for BEM model. (**a**) μ = 0 deg; (**b**) μ = 45 deg; (**c**) μ = 90 deg; (**d**) μ = 135 deg; (**e**) μ = 180 deg.

In Figure 15, the Pitch RAO (rad/m) for 24 Kn ship speed and 2 m wave height are presented. The Pitch RAO has the maximum peak values for the head wave case, around ω = 0.46 rad/s. At the

same time, one may observe that the hydrodynamic models provide better agreement for the beam wave case μ = 90 deg.

**Figure 15.** Pitch RAO (rad/m), *v* = 24 Kn, *hw* = 2 m. (**a**) μ = 0 deg; (**b**) μ = 45 deg; (**c**) μ = 90 deg; (**d**) μ = 135 deg; (**e**) μ = 180 deg.

In Figure 16, the Roll RAO (rad/m) for 24 Kn ship speed and 2 m wave height are presented. The trend of results is maintained for the roll motion as for the 12 Kn ship speed case, the boundary element method being in agreement with the strip theory, for all the heading waves considered. The maximum values appear at beam wave, 0.16 rad/m around ω = 0.47 rad/s. At the same time, as long as the ship is in oblique waves, the values of roll motion decrease.

**Figure 16.** Roll RAO (rad/m), *v* = 24 Kn, *KG* = 10.8 m, *hw* = 2 m. (**a**) μ = 45 deg; (**b**) μ = 90 deg; (**c**) μ = 135 deg.

At 24 Kn ship speed and 8 m wave height for the BEM model, the Heave RAO (Figure 17) records visible differences to the BEM solution at the same speed, but with 2 m wave height, due to hydrodynamic nonlinearities for the wave's circular frequency up to 0.6 rad/s. Figure 18 presents the pitch RAOs for 24 Kn ship speed and 8 m wave height. The differences between the BEM model, with 8 m wave height, and the linear 2D strip model become more significant in comparison to the case with wave height 2 m, due to the recorded hydrodynamic nonlinearities.

**Figure 17.** Heave RAO (m/m), *v* = 24 Kn, *hw* = 8 m for the BEM model. (**a**) μ = 0 deg; (**b**) μ = 45 deg; (**c**) μ = 90 deg; (**d**) μ = 135 deg; (**e**) μ = 180 deg.

**Figure 18.** Pitch RAO (rad/m), *v* = 24 Kn, *hw* = 8 m for the BEM model. (**a**) μ = 0 deg; (**b**) μ = 45 deg; (**c**) μ = 90 deg; (**d**) μ = 135 deg; (**e**) μ = 180 deg.

The Heave RAO peak is 1.071 m/m around ω = 0.5 rad/s for μ = 135 deg (as can be seen in Figure 17) while the maximum value of the Pitch RAO is 0.02 rad/m around ω = 0.45 rad/s (shown in Figure 18), when the ship is in head waves μ = 180 deg.

In Figure 19, a comparison between heave amplitude spectrum, by DFT (Direct Fourier Transformation) method implemented in DYN code, for the two speed cases (12 and 24 Kn) and wave height *hw* = 2 m, heading angle μ = 180 deg, wave circular frequency ω = 0.4 rad/s (with ship–wave *fe* encountering frequencies 0.080 and 0.096 Hz), is made. The same comparison for pitch amplitude spectrum, by DFT method is made. Both diagrams represent the results obtained using the BEM boundary element method by SHIPFLOW. In Figure 20, the same comparison is made for wave height of *hw* = 8 m.

**Figure 19.** SHIPFLOW—Amplitude Spectrum by DFT method, *v* = 12 and 24 Kn, *hw* = 2 m, μ = 180 deg, ω = 0.4 rad/s (*fe* = 0.080 and 0.096 Hz). (**a**) Heave Amplitude Spectrum (m); (**b**) Pitch Amplitude Spectrum (rad).

**Figure 20.** SHIPFLOW— Amplitude Spectrum by DFT method, *v* = 12 and 24 Kn, *hw* = 8 m, μ = 180 deg, ω = 0.4 rad/s (*fe* = 0.080 and 0.096 Hz). (**a**) Heave Amplitude Spectrum (m); (**b**) Pitch Amplitude Spectrum (rad).

Both heave and pitch amplitude spectra by DFT method, in the case of higher ship speed of 24 Kn, have narrowed frequency bands. Additionally, at *hw* = 8 m wave height, the heave and pitch amplitude spectra present a series of higher harmonics, with small amplitudes, due to the hydrodynamic nonlinearities.

In Figures 21 and 22, a comparison between BEM and RANS hydrodynamic models is performed. A comparison of Heave Time Record and Pitch Time Record computed for the KCS vessel at *v* = 12 Kn at heading angle μ = 180 deg and circular wave frequency ω = 0.4 rad/s (*fe* = 0.080 Hz) is presented in Figure 21.


**Figure 21.** NUMECA–SHIPFLOW Time Records comparison *v* = 12 Kn, *hw* = 2 m, μ = 180 deg, ω = 0.4 rad/s (*fe* = 0.080). (**a**) Heave Time Record (m) (**b**) Pitch Time Record (rad).

**Figure 22.** Comparison of NUMECA–SHIPFLOW Amplitude Spectrum by DFT method, *v* = 12 Kn, *hw* = 2 m, μ = 180 deg, ω = 0.4 rad/s (*fe* = 0.080). (**a**) Heave Amplitude Spectrum [m]; (**b**) Pitch Amplitude Spectrum (rad).

A comparison of Heave Amplitude Spectrum and Pitch Amplitude Spectrum computed for the KCS vessel at *v* = 12 Kn, at heading angle μ = 180 deg and circular wave frequency ω = 0.4 rad/s (*fe* = 0.080 Hz) is presented in Figure 22. While for the Heave Amplitude Spectrum, SHIPFLOW and NUMECA show good agreement, for the Pitch Amplitude Spectrum some difference occurs between the two models.

#### **7. Conclusions**

In this paper, for the KCS full-scale model (Section 3), several numerical seakeeping analyses have been performed, in order to compare three hydrodynamic methods for the prediction of the dynamic response on heave, pitch and roll oscillations, in terms of response amplitude operators. The selected hydrodynamic methods are the most used approaches for ships motions, linear 2D strip theory, BEM and RANS methods, involving in-house codes and commercial codes.

A grid convergence study has been conducted to determine the accuracy of the grid used to compute the flow based on the potential BEM approach (Section 5). The study revealed that solution computed by SHIPFLOW is grid independent and the grid error for the grid used is less than 1%.

The BEM and RANS methods, with time-domain solutions, were compared to the linear strip (ST) method, with frequency-domain solution; although the external excitation is a regular wave, the ship's response shows higher harmonics (Figures 19 and 20), due to the hydrodynamic nonlinearities. Therefore, in order to ensure the same reference for all three methods, the nonlinear response amplitude operators (RAOs), obtained using BEM and RANS methods, were computed considering only the main spectral amplitude for the ship-wave encountering frequency, being the only result delivered by the linear 2D strip method. Nevertheless, the BEM and RANS studies have revealed, for the KCS model, reduced nonlinear dynamic response components, even if the regular wave height was up to 8 m and the speed of the ship was 24 Kn (Figures 19, 20 and 22).

For a detailed benchmark between the three hydrodynamic models (Section 2), this parametric study developed for the KCS model (Sections 3 and 4), at one loading case, has been focused on the influence of ships speed (12 and 24 Kn), wave heading angle (0, 45, 90, 135, 180 deg) and height (2 and 8 m), on the ship's dynamic response in regular wave, using deterministic analyses (Section 6).

Comparing the Heave RAO functions between the three hydrodynamic models' results (Figures 7, 9, 12, 14 and 17) (Section 6), for ship speed of 24 Kn versus 12 Kn and 8 m wave height versus 2 m, the BEM method delivers visible differences due to the hydrodynamic nonlinearities, especially for the wave circular frequencies up to 0.6 rad/s. In the case of beam waves, for wave circular frequency within the range of 0.4 to 0.6 rad/s, the BEM values are quite scattered around the 2D strip model results. In addition, visible differences, for head waves, have been seen between the RANS method and the other two (Figure 9).

Comparing the Pitch RAO functions between the three hydrodynamic models (Figures 8, 10, 13, 15 and 18) (Section 6), it turned out that the maximum values are obtained for head wave conditions and are very small for beam wave conditions. As the ship's speed and wave height increase, the hydrodynamic nonlinearities induce visible differences between the results, especially close to the peak values of Pitch RAO. Except for 0.3 rad/s, on head waves, the differences between the RANS and the other two methods are quite scattered for the analyzed wave's frequency domain (Figure 10).

Comparing the Roll RAO functions between the BEM and 2D strip hydrodynamic results (Figures 11 and 16) (Section 6), it turned out that the maximum values are obtained for beam waves condition, being a narrow frequency band response amplitude operator. The differences are noticeable at a beam wave condition close to the RAO peak value, especially due to the different hydrodynamic damping formulation of the two methods.

As an overall analysis, we can conclude that the RAO's values predicted by the linear 2D strip theory can be considered as average values compared to the BEM and RANS results, being suitable for practical seakeeping analyses of mono-hull slender body ships.

Further studies will extend this parametric analysis to more numerical results based on the RANS method considering other ship types, being focused on the sensitivity of the hydrodynamic methods used for ship's motions prediction. Additionally, the next step consists of short-term seakeeping analysis in irregular waves, where, besides the BEM and RANS methods, also a nonlinear strip theory with time-domain solution is to be applied.

**Author Contributions:** Conceptualization, F.P. and L.D.; methodology, F.P., L.D. and S.P.; software, F.P., L.D. and S.P.; validation, F.P. and L.D.; formal analysis, S.P.; investigation, F.P., L.D. and S.P.; resources, F.P., L.D. and S.P.; data curation, F.P.; writing—original draft preparation, F.P., L.D. and S.P.; writing—review and editing, S.P., F.P. and L.D.; visualization, S.P. and F.P.; supervision, L.D.; project administration, S.P.; funding acquisition, S.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The research study presented in this article was developed in the frame of Dunarea de Jos University of Galati, Naval Architecture Research Centre.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **Abbreviations**


#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Journal of Marine Science and Engineering* Editorial Office E-mail: jmse@mdpi.com www.mdpi.com/journal/jmse

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18