*Article* **Vertical Motions Prediction in Irregular Waves Using a Time Domain Approach for Hard Chine Displacement Hull**

#### **Ermina Begovic 1,\*, Carlo Bertorello 1, Ferdi Cakici 2, Emre Kahramanoglu <sup>2</sup> and Barbara Rinauro <sup>1</sup>**


Received: 23 March 2020; Accepted: 4 May 2020; Published: 9 May 2020

**Abstract:** In this paper, the validation of the hybrid frequency–time domain method for the assessment of hard chine displacement hull from vertical motions is presented. Excitation and hydrodynamic coefficients in regular waves are obtained from the 3D panel method by Hydrostar® software, while coupled heave and pitch motions are calculated in the time domain by applying the Cummins equations. Experiments using a 1:15 scale model of a "low-drag" small craft are performed in irregular head and following waves at Froude numbers Fr: 0.2, 0.4, and 0.6 at University of Naples Federico II, Italy. Results obtained by hybrid frequency–time domain simulations for heave, pitch, and vertical accelerations at center of gravity and bow are compared with experimental data and showed high accuracy.

**Keywords:** Cummins equations; vertical motions assessment; time domain simulations; experimental seakeeping; hard chine displacement hull form

#### **1. Introduction**

The assessment of ship behavior in an irregular seaway is one of the most difficult hydrodynamic problems. A large variety of different computational methods have been presented in the past three decades and are discussed in Hirdais et al. [1]. They proposed the subdivision into six levels, where each "level" introduces mathematical models closer to the physical models, generally moving from frequency domain calculations to time domain simulations, from linearized to nonlinear boundary conditions, and from small wave amplitudes to breaking waves, spray, and water flowing onto and off the ship's deck.

The simplest "Level 1" approach considers the potential flow linearized frequency domain methods, while "Level 6" deals with fully non-linear methods like Reynolds Averaged Navier Stokes (RANS) and Smooth Particle Hydrodynamics (SPH). Significant simplifications have been introduced in mathematical models considering separately vertical and horizontal motions and neglecting the viscous effects and the ship's transversal symmetry.

Two main approaches are known—the frequency and time domains. Frequency domain methods have been widely used as strip theories (2D) or as panel methods (3D) with different levels of nonlinearities considered in the mathematical model. They are valid under small wave amplitudes and small ship motions hypothesis, assuming linearized boundary conditions and the linear superposition principle. If the phenomenon involves non-linearities, rising from high wave amplitudes, high advancing speed, or hull forms with strong flare, then the time domain approach should be considered. Today, time domain simulations based on Cummins equation [2] are becoming standard

as they can accommodate nonlinearities due to nonlinear wetted surface and steeper waves and are very fast. The original Cummins' equation considers the fluid memory effects for radiation terms that are calculated using the damping and added mass values in the frequency domain. Although the direct calculation of convolution integrals in Cummins equation is possible (so called "full" time domain calculation), as stated in Perez and Fossen [3], this is very time consuming, and for analysis and control system design, the convolutions are not suitable. One of the possible solutions is the introduction of "parametric model identification" since the convolution is a dynamic linear operator and can be represented by a linear ordinary differential equation-state-space model or transfer functions in the Laplace domain [3]. Perez and Fossen [4] schematized works on parametric model identification developed during time as Time Domain Identification (TDI) and Frequency Domain Identification (FDI) and reported the major contributions by different authors. Armesto et al. [5] presented results of time domain simulations for the motion of the water inside an oscillating water column and a free decay test in the heave of a spar buoy by three techniques—the direct solution of convolution integral, an approximation of the convolution integral with state space model, and Prony's estimation of the convolution. The state space method and Prony's approximations are computationally cheaper, and their results are very close to those from the direct solution of convolution integrals. The authors reported different uncertainties seen in the identification of the coefficients in these methods and concluded that, with the increase of computational capabilities given by actual computers, they recommend the use of a direct integration method to compute the radiation term in Cummins' equation.

Some of the important works for the "Level 2" methodology validation and application cases presented in the last years are:


The present work is further contribution to the validations of the time domain simulations. The 2 DOF time domain simulation code, developed by the authors of this paper and explained in detail in Cakici et al. [8], has been validated for vertical motions in irregular head and following waves of the hard chine displacement boat and compared against experimental data obtained by the authors. The calculation of excitation and radiation terms in the frequency domain has been performed by Hydrostar® software based on the 3D panel methods. Direct computation of the convolution integrals has been used for fluid memory effects. Restoring terms are considered linear and are calculated from the ship main properties. The wave loads in the time domain are obtained by a realization technique from the Hydrostar® frequency domain calculations.

The considered ship is representative of actual small craft trends, oriented to the low drag simple hard chine hull form studied in Bertorello and Begovic [11]. Seakeeping characteristics of hard chine and warped hull at high speed are characterized by strong nonlinearities due to the dynamic trim and constant changing of wetted surface. The general approach for planing hulls vertical motions prediction follows time domain simulation by Zarnick's theory [12] based on the potential flow, full planing condition, and constant deadrise angle in regular waves. On the other hand, flow separation and spray formation due to the hard chine present a complex hydrodynamic problem, which can be correctly studied only by RANS methods and sophisticated mesh modelling techniques with high computational

efforts. At the moment, for a warped hard chine hull operating in displacement and semi-displacement regimes, there is no adequate numerical tool. Thorough validation of the applied method has been performed considering a very demanding hull form tested in irregular head and following seas at Froude numbers Fr = 0.2, 0.4, and 0.6. Simulations are performed with the time step equal to the frequency of sampling, and an identical analysis of the numerical and experimental time series is performed to obtain the fair comparison.

#### **2. Mathematical Model**

A brief overview of the fundamental equations used in the mathematical model for vertical motions of ship advancing in irregular waves is given.

#### *2.1. Cummins Equation for Coupled Heave and Pitch Motion*

Starting from the general Cummins' Equation (1),

$$\left(\overline{\mathbf{M}} + \mathbf{A}^{\otimes}\right)\bar{\eta}(t) + \int\_{0}^{t} \mathbf{K}(t-\tau)\eta(\tau)d\tau + \mathbf{C}\eta(t) = \mathbf{F}\_{\mathbf{W}}(t) \tag{1}$$

the following parameters can be defined:

A∞ is the added mass at infinite frequency,

K is the impulse response (retardation) function matrix, defined as

$$\mathbf{K}(t) = \frac{2}{\pi} \int\_0^\infty [\mathbf{B}(\omega)] \cos(\omega t) d\omega \tag{2}$$

B(ω) is the damping matrix in frequency domain.

C is the restoring matrix.

η(*t*) is the oscillatory response of the ship.

Fw(*t*) is the transient wave force vector that can be created by a linear superposition of frequency domain results for different wave spectra.

In this study, the Cummins' equation is solved for the coupled vertical motions of the ship advancing with the constant speed V, as written in Equations (3) and (4). Subscripts 3 and 5 refer to heave and pitch motions, respectively.

$$\begin{aligned} \left(\Delta + \mathcal{A}\_{33}^{\text{os}}\right)\ddot{\eta}\_{3}(t) + \mathcal{B}\_{33}(\infty)\eta\_{3}(t) + \int\_{0}^{t} \mathcal{K}\_{33}(t-\tau)\eta\_{3}(\tau)d\tau + \left[\mathcal{C}\_{33} + \mathcal{C}\_{33\gets}\right]\eta\_{3}(t) + \mathcal{A}\_{33}^{\text{os}}\ddot{\eta}\_{5}(t) + \\ \mathcal{B}\_{35}(\infty)\eta\_{5}(t) + \int\_{0}^{t} \mathcal{K}\_{35}(t-\tau)\eta\_{5}(\tau)d\tau + \left[\mathcal{C}\_{35} + \mathcal{C}\_{35\gets}\right]\eta\_{5}(t) = \mathcal{F}\_{3}(t) \end{aligned} \tag{3}$$

$$\begin{aligned} \left(\mathbf{I}\_5 + \mathbf{A}\_{55}^{\text{co}}\right)\ddot{\eta}\_5(t) + B\_{55}(\infty)\eta\_5(t) + \int\_0^t \mathbf{K}\_{55}(t-\tau)\eta\_5(\tau)d\tau + [\mathbf{C}\_{55} + \mathbf{C}\_{55\mathbb{C}}]\eta\_5(t) + \mathbf{A}\_{53}^{\text{co}}\ddot{\eta}\_3(t) \\ \mathbf{I}\_4 + B\_{53}(\infty)\eta\_3(t) + \int\_0^t \mathbf{K}\_{53}(t-\tau)\eta\_3(\tau)d\tau + [\mathbf{C}\mathbf{s} + \mathbf{C}\_{55\mathbb{C}}]\eta\_5(t) = \mathbf{F}\xi(t) \end{aligned} \tag{4}$$

Added mass values at infinite frequency are dependent on the ship geometry only and can be obtained as the convergence value from the frequency domain added mass graphs. Bij(V) term stands for the constant damping arises from forward speed of the ship. According to Riemann–Lesbesque, Lemma Bij(V) can be replaced with B(∞) [4,13,14]. Restoring coefficients C33, C35, C53, and C55 represent the constant values and they are calculated from ship main properties.

Convolution term for the restoring coefficients C33C, C35C, C53C, and C55C are also called radiation restoring terms, represent the correction to the hydrodynamic steady forces acting upon the ship surface and can be calculated as follows:

$$\mathbf{C\_{ijC}} = \boldsymbol{\omega\_e}^2 \left[ \mathbf{A\_{ij}^{\rm cs}} - \mathbf{A\_{ij}}(\boldsymbol{\omega\_e}) \right] - \boldsymbol{\omega\_e} \cdot \int\_0^\infty \left[ \mathbf{K\_{ij}}(t) \sin \boldsymbol{\omega\_e} t \right] dt, \text{ i. j.}\tag{5}$$

If regular waves are considered, as in Fonseca and Soares [15], then the CijC will be constant and independent of the encounter frequency. While restoring terms are depending on the ship geometry and mass distribution and are independent of wave frequencies, the radiation restoring coefficient is introduced to accommodate for a correction to the hydrostatic buoyancy force/moments due to the unsteady ship motions. In Fonseca and Soares [16–18] and in Vásquez et al. [19], the radiation restoring and the memory functions are obtained by relating the radiation forces in the time domain and in the frequency domain by means of Fourier analysis. Ma et al. [20] theoretically derived a formulation for radiation restoring coefficient by use of strip theory, and for the considered test cases, it seemed the more consistent formulation. Since the irregular waves are considered in the present study, the effects of radiation restoring are neglected.

#### *2.2. Calculation of Convolution Terms*

Kij(*t*) is the impulse response function, which can be calculated by Equation (6) for the coupled heave and pitch motions as:

$$\mathbf{K}\_{\ddot{\mathbb{I}}}(t) = \frac{2}{\pi} \int\_{0}^{\infty} \left[ \mathbf{B}\_{\ddagger \ddot{\mathbb{I}}}(\omega\_{\mathfrak{e}}) - \mathbf{B}\_{\ddot{\mathbb{I}}}(\infty) \right] \cos(\omega\_{\mathfrak{e}}t) d\omega\_{\mathfrak{e}} \tag{6}$$

where i, j = 3, 5.

The damping values Bij(ω*e*) are found using the 3D panel method implemented in Hydrostar® software by Bureau Veritas, France. Bij(∞) can be obtained as the convergence value from the graphs of damping coefficients as function of the wave frequency. In the present study, three forward speeds are considered, and for each speed, the encounter frequencies are calculated in the range of wave frequencies ω = 0.3 rad/s to 2.1 rad/s. The definite integrals for impulse response functions defined by Equation (6) are calculated for head waves since the encounter frequencies range are sufficiently large, and it covers the following wave frequencies as well. It is noted that once the impulse response functions are calculated for three forward speeds for head waves, these values can also be used for the following wave simulations.

For calculation of the convolution integrals in Equations (3) and (4), the following approximation is applied. If the convolution integrals are split into two additional parts, one can obtain the following statements:

$$\mathbf{X}\_{1} = \int\_{0}^{t} \mathbf{K}\_{33}(t-\tau)\eta\_{5}(\tau)d\tau + \int\_{0}^{t} \mathbf{K}\_{35}(t-\tau)\eta\_{5}(\tau)d\tau \tag{7a}$$

$$\mathbf{X}\_{2} = \int\_{0}^{t} \mathbf{K}\_{53}(t-\tau) \,\eta\_{3}(\tau) d\tau + \int\_{0}^{t} \mathbf{K}\_{55}(t-\tau) \,\eta\_{5}(\tau) d\tau \tag{7b}$$

As noted in Cakici et al. [8], Equations (7a) and (7b) can be approximated as [21–23]

$$\mathbf{X}\_{1} = \sum\_{n=0}^{N} \mathbf{K}\_{33}(n\Delta t) \,\eta\_{3}(t - n\Delta t)\Delta t + \sum\_{n=0}^{N} \mathbf{K}\_{35}(n\Delta t) \,\eta\_{5}(t - n\Delta t)\Delta t \tag{7c}$$

*J. Mar. Sci. Eng.* **2020**, *8*, 337

$$\mathbf{X}\_{2} = \sum\_{n=0}^{N} \mathbf{K}\_{\mathfrak{F}3}(n\,\Delta t)\eta\_{\mathfrak{F}}(t - n\,\Delta t)\Delta t + \sum\_{n=0}^{N} \mathbf{K}\_{\mathfrak{F}5}(n\,\Delta t)\eta\_{\mathfrak{F}}(t - n\,\Delta t)\Delta t \tag{7d}$$

where N = *<sup>t</sup>* 1 <sup>Δ</sup>*t*, Δ*t* is the time step size *t* <sup>1</sup> maximum time value,

The rule of thumb for choosing the time step and the maximum time value is given by the report of McTaggart [24], formulating

$$
\Delta t \approx 0.05 \sqrt{\frac{L}{\mathcal{S}}} \text{ and } t^1 \approx 5 \sqrt{\frac{L}{\mathcal{S}}}
$$

These formulations suggest that the time interval Δ*t* should be sufficiently small to capture the variation of retardation function, and the maximum time value *t* <sup>1</sup> should encompass the time when the retardation functions approach zero. In any case, the maximum time value should be determined according to retardation function plots.

#### *2.3. Heave Force and Pitch Moment Calculation*

To obtain randomized time record of heave force *F*3(*x*, *t*) and pitch moment *F*5(*x*, *t*), first the response spectra of heave force and pitch moment *SHF*/*PM*(ω*E*) are calculated, using the linear superposition principle proposed by St Dennis and Pierson [25] as in Equation (8):

$$\mathcal{S}\_{\rm HF/PM}(\omega\_E) = \mathcal{S}\_{\zeta}(\omega\_E) \times \left| TF\_{\rm HF/PM}(\omega\_{E'} \chi) \right|^2 \tag{8}$$

where

*S*ζ(ω*E*) is the encounter wave energy spectrum, in this work JONSWAP spectrum was used; *TFHF*/*PM*(ω*E*, χ) is the transfer function (RAO) of heave force or pitch moment;

χ denotes the encounter angle;

*SHF*/*PM*(ω*E*) is the heave force or pitch moment response spectrum.

According to linear random wave theory, unidirectional and long-crested waves can be expressed by the sum of finite regular wave components. The instantaneous heave force amplitude can be stated as follows:

$$F\_3(\mathbf{x}, t) = \sum\_{i}^{N} F\_{3i} \cos(\omega\_{i1} t + \varepsilon\_i + \arg[F\_{3i}(\omega\_{i1}, \chi)]) \tag{9}$$

*F*3*<sup>i</sup>* = %2*SHF*(ωei)Δωei represents the amplitude of i-th heave force component, ω*ei* is the i-th encounter wave frequency;

ε*<sup>i</sup>* is the i-th phase lag, randomly assigned between 0 and 2π. It is noted that once random phase lag is chosen, it will be the same for pitch moment;

arg[*F*3*i*(ω*ei*, χ)] denotes heave force phase angle.

#### *2.4. State Space Representation of Mathematical Model of Vertical Ship Motions*

To solve the system of equations, the state space representation is used as explained in the following Equation:

$$\frac{d}{dt}\eta(t) = \mathbf{A}\eta(t) + \mathbf{B}\_1[\mathbf{w}\text{ave}(t) - \text{convo}(t)]\tag{10}$$

where <sup>η</sup>(*t*) <sup>∈</sup><sup>n</sup> is the differentiable state vector, wave(*t*) <sup>∈</sup>*mW* is the wave load input vector, and convo(*t*) ∈ *mC* is the convolution vector. A, B1 are known appropriate state space matrices.

The system matrices, disturbance inputs and states of the system are defined as

$$\mathbf{A} = \begin{bmatrix} 0 & 0 & 1 & 1 \\ 0 & 0 & 1 & 1 \\ -\mathbf{M}^{-1}\mathbf{C} & -\mathbf{M}^{-1}\mathbf{B} \end{bmatrix}, \mathbf{B}\_1 = \begin{bmatrix} 0 & 0 \\ 0 & 0 \\ \mathbf{M}^{-1} \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \end{bmatrix}'$$

$$\text{wave}(t) = \begin{bmatrix} \mathbf{F}\_3 \\ \mathbf{F}\_5 \end{bmatrix}, \text{conv}(t) = \begin{bmatrix} \mathbf{X}\_1 \\ \mathbf{X}\_2 \end{bmatrix}, \ \eta(t) = \begin{bmatrix} \eta\_3 \\ \eta\_3 \\ \eta\_5 \end{bmatrix}$$

The statements in the state space matrices are defined as follows:

$$\mathbf{M} = \begin{bmatrix} \mathbf{A} + \mathbf{A}\mathbf{3}\mathbf{3}^{\mathrm{op}} & \mathbf{A}\mathbf{3}\mathbf{5}^{\mathrm{op}} \\ \mathbf{A}\_{53}\mathbf{3}^{\mathrm{op}} & \mathbf{I}\_{55} + \mathbf{A}\_{55}\mathbf{5}^{\mathrm{op}} \end{bmatrix}, \mathbf{C} = \begin{bmatrix} \mathbf{C}\_{33} & \mathbf{C}\_{35} \\ \mathbf{C}\_{53} & \mathbf{C}\_{55} \end{bmatrix}, \mathbf{B} = \begin{bmatrix} \mathbf{B}\mathbf{3}\mathbf{3}^{\mathrm{op}} & \mathbf{B}\mathbf{3}\mathbf{5}^{\mathrm{op}} \\ \mathbf{B}\mathbf{5}\mathbf{3}^{\mathrm{op}} & \mathbf{B}\mathbf{5}\mathbf{5}^{\mathrm{op}} \end{bmatrix}.$$

#### **3. Experimental Campaign for Hard Chine Displacement Hull Form**

#### *3.1. Experimental Setup and Model Description*

The experimental campaign was performed in the Towing Tank of University of Naples Federico II (UNINA), Italy. The towing tank dimensions are 135 × 9 × 4.2 m and it has a wave generator capable of generating waves from 0.20 /hz to 1.25 Hz and towing carriage with maximum speed of 8 m/s. The simulations are performed for the hard chine displacement hull form, studied in Bertorello and Begovic [11] and shown in Figure 1.

**Figure 1.** Hull form from Bertorello and Begovic [11].

A wooden model of 2.00 m LOA was tested at one displacement (342.17 N) in irregular waves, completing the previously performed experimental campaign in calm water and regular waves in Bertorello and Begovic [11]. The values of main characteristics for the tested model are given in Table 1.


**Table 1.** The main characteristics of the tested model.

The model has been ballasted again, and the center of gravity position and moments of inertia were measured by the model inertial balance shown in Figure 2.

**Figure 2.** The model on the inertial balance at the Towing Tank of University of Naples Federico II (UNINA).

The towing point was at x = 0.87 m from the stern, and the weight of the mechanical arm is the part of the model ballast. The model scale ratio is λ = 15, which corresponds to 30 m LOA.

Measurement of pitch and heave motions were performed by the mechanical arm R47, which holds the model restrained to surge, sway, roll, and yaw. Accelerations were measured at the CG and at the bow. Three ultrasonic wave gauges were used for the wave measurements. Two wave gauges were aligned in the front of the model at the distance of 1.93 m from the R47, and one was aligned with the R47 on the tank side. All data are sampled at a frequency of 500 Hz without filtering. The experimental set up is shown in Figure 3.

**Figure 3.** Experimental setup of low drag hard chine hull at UNINA.

#### *3.2. Experiments in Irregular Waves*

The experiments in irregular head and following waves were performed to validate the hybrid frequency–time domain numerical method. The standard JONSWAP theoretical spectra, defined as

$$S\_{\subset -\text{JONSWAP}}(\omega) = A\_{\mathcal{V}} S\_{\text{PM}}(\omega) \mathcal{V}^{\exp\left(\frac{-\mathbb{1}}{2} \left(\frac{\omega - \omega \cdot p}{\omega \cdot \omega p}\right)^{2}\right)}\tag{11}$$

where

SPM—Pierson-Moskowitz spectrum, as defined by DNV-RP-C205 [26];

γ—non-dimensional peak shape parameter;

σ—spectral width parameter, σ = 0.07 for ω ≤ ωP, σ = 0.09 for ω > ωP;

Aγ—normalizing factor, defined as *A*<sup>γ</sup> = 1 − 0.287 ln(γ);

were used for the representation of irregular waves. To obtain enough encounters, the runs were repeated 2–10 times, depending on the model speed and heading. For the head sea, the complete time series is 140 s, resulting in a minimum of 120 wave encounters. In the following sea, due to the very low (or negative) encounter frequency, the total time of experimental series is around 200 s, ending up with a minimum of 30 encounters. The JONSWAP spectrum parameters for model and ship scale are reported in Table 2. The example of the measured spectrum is given in Figure 4, together with the ideal wave and encounter spectra. Experimental set up and tested model velocities were identical to those in regular waves. In Table 3, the number of encounters, significant wave height, and m0 are given for the spectra measured at the following three Froude numbers: 0.2, 0.4, and 0.6 in head waves.

**Table 2.** Wave spectrum characteristics.


**Figure 4.** Ideal wave and encounter vs. measured JONSWAP spectra at Froude number Fr = 0.4.

**Table 3.** Measured spectra properties-head waves.


Examples of the time series of 140 s registration of the wave, heave, and pitch at Fr = 0.4 in the head and following seas are given in Figures 5a and 6. These show the appreciated different number of wave encounters, and consequently of heave and pitch oscillations for the same registration time in the head and following sea. In Figure 5b, measured accelerations at center of gravity CG and at the bow are given. In the following waves, the measured accelerations were very low, and the signals were noisy, so for the sake of clarity, these data have been neglected for the further comparisons. It can be noted from Figure 5b that the accelerations are normalized by g, so they oscillate around 1g value.

(**a**) **Figure 5.** *Cont.*

**Figure 5.** (**a**) Example of measured wave, heave, and pitch at Fr = 0.4 in head waves; (**b**) example of measured accelerations at Fr = 0.4 in head waves.

**Figure 6.** Example of measured wave, heave, and pitch at Fr = 0.4 in the following waves.

#### **4. Hydrodynamic Coe**ffi**cients and Exciting Forces Calculation**

To be able to calculate motions in the time domain, the hydrodynamic coefficients for added mass and damping and exciting forces are needed as input, and they were obtained by the hydrodynamic software HydroSTAR® developed by Bureau Veritas.

HydroStar® is a 3D diffraction/radiation potential theory software based on the Green function method for wave–body interactions, and it provides a complete solution of first- and second-order wave loads. Calculations have been performed in ship scale for the following three Froude number cases: 0.2, 0.4, and 0.6, which correspond to ship speed values of 3.385, 6.77, and 10.155 m/s. For each of the Froude number cases, the frequencies of forward incoming regular waves have been considered in the range from 0.3 to rad/s 2.1 rad/s, with the step of 0.05, corresponding to wave periods of 2 s to 20 s.

The panelized hull geometry representation of the model in HydroSTAR® is shown in Figure 7.

(

\$ NJ
P

**Figure 7.** Hull geometry in Hydrostar® software.

For the simulation in the time domain, the following hydrodynamic coefficients and transfer functions TF values have been calculated:


All diagrams show ship values as a function of wave frequencies ω (rad/s). The model values, used as input for the time domain simulations, are obtained by applying the scaling law as indicated in Table 4 and taking values at the highest wave frequencies.

( (

**Figure 8.** Added mass coefficients in the range of considered wave frequencies.

**Figure 9.** Damping coefficients in the range of considered wave frequencies.

**Figure 10.** Transfer function (TF) of heave and pitch exciting force.

#### **5. Time Domain Simulations**

#### *5.1. Numerical Results*

The simulations with the developed code are performed in model scale to be able to compare the time series directly with the measured data. The added mass and restoring coefficients from Hydrostar® are recalculated in model scale as reported in the Table 4. The values of the added mass coefficients at infinity used in further calculations for retardation functions are reported in the last column of the Table 4. While at the sufficiently high frequencies, added mass coefficients are converging to the same value for all Froude numbers, the values of the damping coefficients at infinity are dependent on the forward speed of the model. In Table 4, for the sake of compactness, only the values for Fr = 0.4 are reported.

Retardation function has been calculated for 5 s, with the time step equal to 0.002 s, and an example at Fr = 0.4 is reported in Figure 11. Details of the 30 s of time series of heave force and pitch moment acting on the hull surface calculated according to the realization technique, in head and following waves, are given in Figures 12 and 13.

**Table 4.** Scaling of hydrodynamic coefficients from Hydrostar® to model scale.

**Figure 11.** Impulse response functions for the Fr = 0.4.

**Figure 12.** Detail of wave excitations time series for the Fr = 0.4 in head sea.

**Figure 13.** Detail of wave excitation time series for the Fr = 0.4 in following sea.

For the solution, The fourth-order Runge Kutta Method is used, where the time step size is taken as 0.002, and the simulation time has been set equal to 140 s for head waves and 280 s for following waves to have a fair comparison with the experimental data. In Figures 14–19, the examples of comparison of 30 s of simulated and measured time histories of motions and accelerations are given. It has to be noted that this comparison has to be seen only as a qualitative control of simulated data. Both time series are random processes and are given for the same interval of time, without trying to overlap signals. Therefore, they can be compared only statistically.

**Figure 14.** Example of simulated and measured heave motion at Fr = 0.4 in head seas.

**Figure 15.** Example of simulated and measured pitch motion at Fr = 0.4.

**Figure 16.** Example of simulated and measured cg accelerations at Fr = 0.4 in head seas.

**Figure 17.** Example of simulated and measured bow accelerations at Fr = 0.4 in head seas.

**Figure 18.** Example of simulated and measured heave motion at Fr = 0.4 in following seas.

**Figure 19.** Example of simulated and measured pitch motion at Fr = 0.4 in following seas.

#### *5.2. Results Comparison*

Both numerical and experimental results have been analyzed in the time domain "peak to peak analysis." The complete set of performed calculations and experiments are summarized in Table 5. Tables 6 and 7 report the number of peaks in each time series, significant height, 1/10th of the highest value, maximum height value, and mean values. The differences among mean values have been expressed in the last column. It can be noted from Table 7 that calculated vertical accelerations in the following waves are very low, in the same order of magnitude as the measured ones. As the analysis of measured accelerations is strongly related to the filtering technique and cut off frequencies for the case of following waves, they have not been considered in the comparison.


**Table 5.** Summary of performed experiments and simulations.


**Table 6.** Head seas results comparison.

**Table 7.** Following seas results comparison.


In Tables 6 and 7 the comparison for head and following seas, respectively, are given.

As deduced from Table 6, the differences between numerical and experimental heave-pitch motions at head waves in terms of significant mean values remain in the range of 10% for all tested forward speeds, which indicates that the predicted motions have great level of accuracy when it is compared to experimental data. Only for the accelerations at CG at Fr = 0.2 and bow accelerations at Fr = 0.6 differences are 28% and 18%, respectively.

As shown in Table 7, in the case of following waves, the differences between numerical and experimental heave-pitch motions are around 10%, except for Fr = 0.4, where the difference is around 24%. It has to be noted that at this Fr, encounter frequencies are close to zero, which is well known situation where there are theoretical restrictions on damping coefficients. Since this Fr also corresponds to the case in which wave celerity and ship speed are the same, it is expected that linear tools might fail to represent this case. As stated previously, in following waves, both measured and calculated vertical accelerations are very low, in the order of magnitude 0.02g, and therefore, this comparison has not been considered.

#### **6. Conclusions**

In the present work, a hybrid frequency–time domain 2 DOF simulations code is presented. The prediction of the wave loads in the time domain was made by a realization technique from the linear frequency domain calculations by HydroStar® software. Direct computation of the convolution integrals was used for the fluid memory effects. The validation of the method was performed comparing the numerical results with the experimental data for a hard chine displacement hull at three model speeds in irregular head and following waves.

Results showed that the applied method for the prediction of vertical ship motions in the time domain is reliable for head and following wave cases. The accuracy of the proposed numerical implementation is dependent on the calculation of the impulse response function, which is, in turn, a function of damping values in the frequency domain. As underlined, a 3D panel method is used in this work. If strip theory is used for the evaluation of damping terms, attention must be taken because there are theoretical restrictions at the low-frequency region that directly affect the impulse response function signal. In this situation, filtering for the frequency dependent damping function may help despite the loss of accuracy. Very challenging validation of the developed method has been performed considering hard chine warped displacement hull, relative ship speed corresponding up to Fr = 0.6, and the following sea at the low and negative encounter frequencies. It has to be highlighted that the calculations in HydroStar® have been performed at the static trim, although at the Fr = 0.6 the dynamic trim is around 1.5 deg and this has been neglected in the calculation of the hydrodynamic coefficients. Obtained results in terms of percentage difference of the mean values is at max 15% in all cases, except in following sea at a model speed close to the wave celerity. In that case, the differences are around 25%. For the following sea, the shorter times of the experimental series and the obvious benefit in using the numerical tool, reasonably accurate, and fast to obtain vertical motions in the time domain in the design procedure should all be noted. Future work will consider the implementation and validation of six DOF methodologies.

**Author Contributions:** Conceptualization, data curation, investigation, methodology, validation, writing—original draft, writing—review and editing, supervision: E.B.; Conceptualization, Data Curation, Investigation, Writing—Original Draft: C.B.; Conceptualization, data curation, methodology, software, validation, writing—original draft, writing—review and editing: F.C.; Formal analysis, software, validation, writing—review and editing: E.K.; Data curation, investigation, writing—review and editing: B.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** Ferdi Çakici and Emre Kahramano ˘glu were supported by the Scientific and Technological Research Council of Turkey (TÜB˙ ITAK).

**Acknowledgments:** The calculation by HydroSTAR® software were possible thanks to the cooperation agreement between Bureau Veritas and the Department of Industrial Engineering of University of Naples Federico II (https://www.docenti.unina.it/webdocenti-be/allegati/materiale-didattico/576434).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 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/).
