Computational prediction of hydrodynamic coefficients for heave motion

Abstract

Keywords: Hydrodynamic coefficients; URANS; Wigley hull; Heave motion

Introduction

Prediction of ship motion characteristics of the ships is one of the most challenging research area in ship hydrodynamics. Several studies have been carried out in order to deal with this problem by the researchers and the studies are still proceeding. The very first contributions came with the Ursell’s milestone studies. He managed to solve the linear boundary value problem of an oscillating circular cylinder on the free surface in the frequency domain (Ursell, 1949a, 1949b). In his method, the velocity potential was represented as a summation of an infinite set of multipoles which satisfies the free surface boundary condition. His studies opened a wide perspective to the theoretical exploration of ship motion problems. The basic assumption of his theory was that the forward speed assumed to be zero. The other

assumption was that both incident wave and the body oscillate with the same frequency. Then, Tasai combined the conformal mapping transformation technique with Ursell’s first studies to get a solution for more realistic ship sections (Lewis, 1929; Tasai, 1959a, 1959b). His method was related to transforming the ship sections into a half circle using a scale factor and two or more mapping coefficients. Especially the ships with round bilge were represented well by this method. However, this method was not sufficient to obtain good results for the hulls like SWATH (Small Waterplane Area Twin Hull) or hulls with transom sterns and blunt or bulbous bows. Therefore, Frank worked on this problem and developed a method (Frank, 1967). His method, which is called Frank close-fit, was a better option in representing the hydrodynamic coefficients of the arbitrary ship sections because the solutions by this method could be used for ar-

*Corresponding author. *E-mail address: fcakici@yildiz.edu.tr Published by Yıldız Technical University Press, İstanbul, Turkey Copyright 2021, Yıldız Technical University. This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/by-nc/4.0/).

bitrary cross sections. The method adopts Green’s function which satisfies the linear free surface boundary condition. The density of the sources is an unknown function to be found from integral equations obtained by applying the body boundary condition (Frank, 1967). However, these aforementioned studies were limited with 1-DOF (degree-of-freedom) motion of a body in a fluid due to difficulties in coupling of hydrodynamic coefficients on each other. Salvesen et al. (1970) presented the original strip theory by using 2D hydrodynamic coefficients which were calculated by Frank’s method. As reported in the study of Beck and Reed, strip theory is used for seakeeping analyses intensively with its quickness and wide application area in the early of 2000s (Beck and Reed 2001; Tezdogan et.al 2015). However, this theory loses its effectiveness at particularly low frequency motions and high speeds. It means that the strip theory considers the forward speed correction in a simple manner. As another CFD (Computational Fluid Dynamics) approach, URANS and the additional turbulence equations, which are discretized by implementing finite volume method (FVM), have recently been used to calculate the vertical ship motions. The nonlinear URANS approach have been used by many researchers to find the hydrodynamic coefficients regarding the added mass as well as the damping. One of these studies was the work of Querard et al. (2009) where they dealt with the computations of added mass and damping of 2D sections. Calculations were done for a wide range of frequency. The main object of their work was to make a comparison with the results obtained by the potential theory and the experiments conducted by Vugts (1968). Bonfiglio et al. (2016), in their study, calculated radiation terms by performing URANS calculations, while excitation terms were found by the potential method. They obtained the radiation terms for 2-D ship sections. Then, they used the strip theory method to get 3-D radiation terms. In this study, the uncoupled hydrodynamic coefficients related to heave motion are obtained using URANS approach. Several CFD analyses are performed at Fn=0.3 and numerical results are then compared with the experimental data and strip theory outputs.

Figure 1. The sections of the Wigley hull. Table 1. Principal particulars of Wigley hull

MAIN Particulars And The Investigated Cases

The well-known Wigley hull is chosen for the hydrodynamic analyses. The numerical simulations are performed for the bare hull condition. The body plan of the Wigley hull is shown in Figure 1. The principal particulars of the model are given in Table 1. An Earth-fixed xyz-Cartesian coordinate system is selected for the solution domain. The xy-plane represents the free surface which is initially assumed as calm water at the draught of the model. z-axis is the vertical axis normal to the xy-plane. In the computation of radiation terms, the ship model is forced to oscillate in heave direction. Calcu-

Table 2. Definition of cases for strip theory and URANS calculations Case no

Methods

lations in the present method are performed at six different frequencies at Froude number 0.3. Dispersion relation for a regular wave in a deep water is given as follows: ω2=gk(1) Where g denotes the gravitational constant and k denotes the wave number which can be calculated by k=2π/λ. Here λ denotes the wave amplitude. Out of a large experimental database from the study of Journée (Journée, 2003), the cases shown in Table 2 are used. The amplitude of the oscillations (a) are selected as ak=0.025 on each case in the URANS calculations which are the same with the experiments. It should be noted that there is also one case in the experiments of Journée in which the oscillation frequency is 3.02 rad/s (case no 1). However, that case is excluded in this paper since the length of the wave created by the oscillation of the ship is

Mathematical Model

Strip Theory Calculations One can use the following strip theory formulations in order to find total heave added mass and damping (Salvesen et al, 1970). Transom correction terms do not include in the Equation (2) and Equation (3) since the Wigley hull does not have a transom stern.  

Figure 2. A square domain. pressure and the velocity. Volume of Fluid (VOF) method is used to represent the free water surface. The commercial CFD software Star-CCM+ is used to discretize the URANS equations by implementing FVM (CD-Adapco, 2014). Calculations of Hydrodynamic Coefficients The ship model is forced to make harmonic heave motion with several constant encounter frequencies and amplitudes while she has a forward speed (Fn=0.3). Then, A33 and B33 are calculated by applying Fourier analyses for each cases. The related formulations can be found in "Fuorier Series Expansion" section. Time Step Size Selection Time step size is selected to be 1/28 of T for the analyses. Here T denotes oscillation period which can be calculated by 2π/ω. The size of time steps in used this study are considered to be more accurate than the value recommended by ITTC (2011). Computational Domain The radiation forces and moments are calculated by using a square-shaped computational domain (Fig. 2). The sizes of the computational domain are chosen large enough in or-

Here a33 and b33 denote added mass and damping coefficients of 2D sections, respectively. For the Wigley hull, these coefficients are found by using a commercial software which uses three parameters Lewis Conformal Mapping technique. URANS Calculations An unsteady RANS solver based on a finite volume method is used in order to solve the flow around the ship hull. Reynolds stresses are calculated by using the standard k-ε with a near-wall treatment that is the most commonly used two-equation turbulence model in industrial applications (Querard et al., 2008). Segregated flow model is adopted for the solver and a second-order upwind scheme is applied to discretize the convection terms in the URANS equations while a first-order temporal scheme is used to discretize the unsteady terms. The SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithm is used to solve the continuity and momentum equations by coupling the

der to prevent the boundary effects. It should be noted that, since the problem has a symmetrical behavior, only half of the computational domain is used for the simulations in order to reduce the computational cost (Fig. 2, 3). The fluids (water and air) are considered to flow in the positive x-direction. The intersection of forward perpendicular and the free surface are set as the origin of the x-y-z Cartesian coordinate system (Fig. 3). To avoid the boundary effects, the faces in the negative x, negative y, positive z (top) and negative z (bottom) directions are considered as velocity inlets. Ship boundaries are defined as no-slip walls where the normal and tangential velocities are zero. Thus, kinematic boundary condition and no-slip condition are satisfied on the hull surfaces. The side face located at x-z plane (where y=0) is assigned as symmetrical boundary (Fig. 3). The computational domain dimensions are given in Table 3.

Grid Generation Finite volume method is applied to discretize the computational domain. Fully hexahedral elements are used in order to generate the computational grid. Local grid refinements are applied to refine the grid around the hull and free surface. Total cell numbers in the present study are about 1 million. Grid structure around the hull is given in Figure 4. To activate the near wall treatment in the turbulence model wall, y+ value is set to be in the range of 30−300 by creating structured hexahedral fine grid at very near of the hull. The black lines in Figure 4 represent the forward perpendicular and the free surface position. The cell sizes are gradually increased with a fixed ratio starting from the boundary layer of the hull to the domain boundaries. Overset grid type is used in the simulations.

Fourier Series Expansion An unsteady time histories of forces/moments, (t), can be represented as follows in Equation (4):

Urans

The detailed information about the overset mesh implementation can be found in related paper (Cakici et al. 2017; Tezdogan et al. 2015).

n=1,2,3... In Equation (4), 0, ωe and n represent zeroth harmonic of the signal, encounter frequency and nth harmonic amplitude of the signal, respectively. The first harmonics components of the signal a1 and b1 are calculated by using Equation (5) and Equation (6)  

In these equations, T denotes the sampling time. The first harmonic of the heave force or pitch moment can be calculated as following:  Figure 4. Grid structure at x/L=0.60.

In forced heave simulations, the coefficients are derived as follows:

Figure 5. Time histories of total force in heave direction. (8)

Here, β3 denotes the phase angle between ship heave motion and heave force. The coefficients A33, B33 and C33 stand for heave added mass, damping and restoring respectively. It also should be noted that C33 which is independent from the frequency, is computed 5883 (kg/s2). Since this value has to be subtracted from total measured force which is obtained by URANS, it is set to minus in the Equation (8).

spectively. Convergence condition (R) of the numerical study and apparent order of p can be obtained by Equation (11) and (12): (11)

The approximate relative error and extrapolated relative error are:

Verification And Validation Study

In the present study, Grid Convergence Index (GCI) method which is based on Richardson Extrapolation (Richardson, 1910), is implemented for the verification study. The GCI method was developed by Roache (1998) and it has been implemented on numerous studies in recent years. A methodology which is explained in detail by Celik et al. (2008) is applied to determine uncertainties of the grid spacing and time step size. The grid and time step size are refined systematically. The refinement factor (r) is selected as √2 because it is often used in CFD applications (Tezdogan et al. 2016; Sezen et al. 2018). The methodology which is put forward by Celik et al. (2008) can be summarized as follows: The difference between the solution scalars ( ε ) should be determined by Equation (10): 

Figure 6. Free surface deformations at Fn=0.3 (case no 3, t=6.02 s).

In these equations; ϕ1, ϕ2 and ϕ3 refer to the solution of fine, medium and coarse mesh grid or time step size, re-

Case no 3 is used for verification and validation study. The grid numbers are summarized in Table 4: In order to obtain nondimensionalized damping coefficients Equation (16) and Equation (17) are used: (16)

Where ρ denotes the water density. Nondimensionalized damping coefficients of forced heave motion obtained by

Figure 7. Comparison of added mass values (by Strip Theory and URANS) with experimental data.

Figure 8. Comparison of damping values (by Strip Theory and URANS) with experimental data. coarse, medium and fine time step size and grid spacing solutions are shown in Table 5. Please note that relative difference of the damping values between fine and medium time step size is approximately 1.1%. Therefore, medium time step size (Te/28) is selected in order to reduce the computational cost. The fine grid and medium time step size result (bold ones) is compared with the experimental results (Journée, 2003) for the validation. Table 6 shows the relative difference between the numerical result and experimental data.

Results And Discussion

Before presenting the added mass and damping values associated with forced heave motion, it may be useful to see the time histories of the total force (including the displacement force) obtained by URANS calculations (Fig. 5). In the same graph, the corresponding Fourier signal is also given. As it is known, a floating body moving with an oscillatory motion in the vertical plane with zero advance speed generates circular wave patterns which start at the center of motion and moves far away from the origin. The free sur-

face deformations are presented in Figure 6 for forced heave motion (case number 3) at 6.02 s. Since the ship has an advance speed, it is hard to observe the radiation waves in the dominant wave system as it is seen in Figure 6. The results obtained by using URANS, experimental data and strip theory are given with Figure 7 and Figure 8. It is obvious that the URANS approach underpredicts the added mass (A33) value at the lowest oscillation frequency. As the frequency increased, the URANS results are generally becoming close to the experimental results. As it is seen from Figure 7, one can easily see that strip theory outputs are able to catch the trend of experimental data but differences can be noticed from Table 7. Please note that ω'=ω/(g/L)0.5. For the whole frequency range, the URANS results connected with the damping values (B33) are in better agreement with experimental data (Fig. 8). As the low frequency region, strip theory results overpredict the value of B33. The discrepancies can be seen from Table 7. It is interesting to note that both calculation results are becoming very close as the frequency increases.

Conclusions

In the present work, the radiation coefficients associated with the forced heave motion are obtained by URANS approach. The Wigley hull model is chosen for the calculations of added mass and damping coefficients (A33, B33). Numerical simulations are performed for six different oscillation frequency while the ship has a forward speed at the Froude number 0.3. The first case (the lowest oscillation in the experiment) is excluded in CFD analyses in purpose due to requirement of a larger solution domain. The 3D ship model is forced with several frequencies and the hydrodynamic coefficients are found via Fourier series expansion. Numerical results are compared with those obtained by experimental data and strip theory. The results have revealed that URANS approach can be used for more accurate prediction of the radiation coefficients for heave motion. Since the theory takes into consideration the viscous effects, nonlinearities in free surface and ship geometry, it can be a good tool to solve some specific problems for which strip theory is insufficient.

Acknowledgments

The first author was supported by ASELSAN Graduate Scholarship for Turkish Academicians.

Remark From The Authors

This paper has been presented at the 3rd International Naval Architecture and Maritime Symposium, 24−25 April 2018, Yıldız Technical University, Istanbul, and published in its proceedings, pp. A21-A34, e-ISBN: 978-975-461548-7, http://www.int-nam.yildiz.edu.tr/2018/cd/INTNAM_2018_e_Proceedings.pdf.

Share and Cite