Numerical study on ship-generated unsteady waves based on a Cartesian-grid method ＊

Kyung-Kyu Yang, Masashi Kashiwagi, Yonghwan Kim

1. Department of Naval Architecture and Ocean Engineering, Seoul National University, Seoul, Korea

2. Department of Naval Architecture and Ocean Engineering, Osaka University, Osaka, Japan

Abstract: The added resistance of a ship in waves can be related to ship-generated unsteady waves. In the present study, the unsteady wave-pattern analysis is applied to calculate the added resistance in waves for two modified Wigley models using a Cartesian-grid method. In the present numerical method, a first-order fractional-step method is applied to the velocity-pressure coupling in the fluid domain, and one of volume-of-fluid (VOF) methods is adopted to capture the fluid interface. A ship is embedded in a Cartesian grid, and the volume fraction of the ship inside the grid is calculated by identifying whether each grid is occupied by liquid, gas, and solid body. The sensitivity to the location of measuring position of unsteady waves as well as the number of solution grids is examined. The added resistance computed by direct pressure integration and wave pattern analysis is compared with experimental data. In addition, nonlinear characteristics of the added resistance in waves are investigated by detailed analyses of unsteady flow field and resulting wave pattern.

Key words: Unsteady wave pattern, added resistance in waves, Cartesian-grid method, grid convergence index

Introduction

When a ship is moving in calm water or incoming waves, ship-generated waves can be observed. They can be decomposed into steady and unsteady waves. The former waves are known as Kelvin waves. There have been many studies of the Kelvin wave by means of theoretical, experimental,and computational methods[1-3]and the Kelvin wave is closely related to the wave making resistance of a ship advancing in calm water

On the other hand, there can be an additional steady wave component due to the mean value of interaction between the Kelvin wave and unsteady wave in the presence of incoming waves. The unsteady waves are mainly composed of the radiation waves induced by ship’s oscillation and the diffraction waves induced by scattering of incident waves due to a fixed body. The unsteady waves appear in the presence of incident waves. Thus, the resistance of a ship in waves will increase in general due to shipgenerated unsteady waves compared with the calm water resistance. This is a main component of the added resistance of a ship in waves, which is of practical importance for designing the propulsion power of a ship. Furthermore, for reduction of greenhouse gas emission, measuring the energy efficiency levels, such as the energy efficiency design index (EEDI), is required for newly built ships and the added resistance in waves is a crucial factor to determine the EEDI[4]. For this practical requirement,many studies have been conducted to estimate the added resistance in waves along with investigating associated physical phenomena[4-9].

In addition to the practical purpose of predicting the added resistance of a ship in waves, the study on unsteady waves itself is required to understand underlying physics connected with ship-generated waves. The analyses of unsteady wave pattern have been performed based on the pioneering work of Maruo[10], who derived the relation between the added resistance and ship-generated waves at far-field based on the energy and momentum conservation principles.Ohkusu[11]proposed a method to estimate the added resistance in waves by measuring ship-generated unsteady waves and calculating the wave amplitude function, the so called Kochin function. An improved method to measure diffraction and radiation waves around a ship advancing in waves was also proposed by Ohkusu and Wen[12]and the added resistance was predicted using measured two-dimensional diffraction and radiation waves. Further thorough investigation of the ship-generated unsteady waves in terms of nonlinear effects and hydrodynamic relations to the added resistance were conducted by Kashiwagi[13].The results showed that the discrepancy of peak value of the added resistance between direct measurement by a dynamometer and unsteady wave analysis was due to nonlinearity in the unsteady wave. Thus, the reproduced unsteady waves by superposing linear wave components in the diffraction and radiation problems provided the added resistance with an acceptable accuracy. Further comparisons of unsteady wave patterns and profiles in the diffraction and radiation problems were performed with the computational results using a time-domain higherorder boundary element method in He and Kashiwagi[14]. Recently, Sadat-Hosseini et al.[15]showed the distribution of unsteady waves in terms of Fourier components for the KRISO’s very large crude carrier 2 (KVLCC2) using a code named CFD Ship-Iowa, which is based on Unsteady Reynolds Average Navier-Stokes (URANS) equations. They concluded that the added resistance is mainly caused by the generated radiation and diffraction waves initiating at fore-body shoulder and transom corner and diverging from the ship.

In the preset study, the wave-induced ship motions and added resistance of two different modified Wigley hulls are computed based on a Cartesian-grid method[16,9], which is recently named SNU-MHL-CFD ver. 2.0 (Seoul National University Marine Hydrodynamics Laboratory CFD). It can be applied to simulate highly nonlinear wave-ship interaction problems. In this code, the force acting on a ship is obtained by integrating the pressure on the ship surface and the added resistance is calculated by subtracting the calm water resistance from the total resistance of the ship in waves. The unsteady wave analysis is also applied to calculate the added resistance in waves using computed unsteady waves along a straight line parallel to the advancing direction of the ship. Before comparing the results with experimental data, sensitivity studies for the lateral distance of wave measuring and the grid spacing are performed. The uncertainty of grid resolution is quantified based on grid-convergence tests and computed results of ship motion responses, added resistance in waves, and unsteady wave-pattern are compared with experimental and/or computational results based on potential flows. Finally, different wave amplitudes are considered to examine the nonlinear characteristics in the flow field and resulting unsteady wave pattern.

1. Methodology

1.1 Cartesian-grid method

The present numerical code, SNU-MHL-CFD ver. 2.0 was developed by Yang et al.[16]and improved by Yang and Kim[9]. It was validated to calculate the ship motions[16]and the added resistance in waves[9]. Because details of the present numerical method are described in those papers, here we briefly explain the numerical procedure. The wave-ship interaction problem was considered as a multi-phase problem with water, air, and solid phases. A solid body was embedded in a Cartesian grid, and to identify the different phases in each grid, the volumefraction functions,mφwere defined for the liquid(m=1), gas (m=2), and solid body (m=3) as shown in Fig. 1. To capture the free surface, the tangent of hyperbola for interface capturing (THINC)scheme[17]which is one of volume of fluid (VOF)methods was used by a weighted line interface calculation (WLIC) method[18]. In the fluid domain,the governing equations for an incompressible and inviscid fluid are the continuity and Euler equations,which are written in conservative forms as follows:

whereΩindicates the control volume,Γrepresents the control surface enclosing the control volume,nis the outward unit normal vector onΓ,ρis the fluid density andpandudenote the pressure and velocity vectors, respectively. In addition,fbindicates the body-force vector. The velocity and pressure were coupled using the fractional-step method and the spatial discretization was carried out based on the finite-volume approach with staggered variable allocation.

Fig. 1 (Color online) Coordinate system in the present numerical code

A ship can be represented with triangular surface meshes and the volume fraction of a solid body embedded in a Cartesian-grid system was calculated using a level-set algorithm. That means the signed distance from the ship surface for each Cartesian-grid point is calculated using a simple transformation of a three-dimensional triangle into a two-dimensional unit right triangle. After obtaining the signed distance field from the triangular surface for each Cartesian-grid point, the volume-fraction function can be calculated using a smoothed Heaviside function[16].

The body boundary condition was imposed using a volume-weighted formula[19]

wherenis the surface normal,IandNindicate the identity tensor and normal dyad, respectively.Because the present numerical method identifies a ship in a Cartesian-grid using the volume-fraction function in each computational cell, the ship interface has finite thickness. The iso-surface of the volumefraction function of the body, which is equal to 0.5,can be considered as the ship surface. Thus, the present numerical method has a limitation on wavebody interaction problems with very sharp geometry.

The force and moment acting on the body are calculated as follows:

wherenFacedenotes the number of triangular surfaces,pl,nland ΔSldenote the interpolated pressure, normal vector, and area of the lth triangular surface, respectively.is the center coordinate of each triangular surface. The added resistance can be calculated by subtracting the calm-water resistance from the total resistance in the presence of the incident wave. In both cases, the same mesh around the ship should be used to obtain a consistent added resistance,and it should be noted that there may be a difference in the calm-water resistance with experimental data because the viscosity was ignored in the present numerical method. The wave elevation,ηat a certain point (xi,yj) was obtained by summing up the volume fraction of water,φ1as follows

whereNzindicates the total number of grid points in thez-direction, Δzkmeans the grid spacing ofk- th grid inz-direction. It should be noted that the depth of numerical tank should be subtracted to make the wave elevation becomes zero in calm water. The other details of numerical method can be found in Refs. [16, 9]

1.2 Unsteady wave analysis

Kashiwagi[13]investigated the effects of nonlinear ship-generated unsteady waves on the added resistance by use of the unsteady wave analysis method for blunt and slender modified Wigley hulls.In the present study, this method is also applied to calculate the added resistance in waves for both models using computed unsteady waves. The added resistance in head wavesRAWcan be obtained in terms of the Kochin functionC(k) as follows:

Here,k0=ω02/gis the wave number of incident wave,0ωis the circular frequency of incident wave andgis the gravitational acceleration.Uis the ship speed andωeis the circular encounter wave frequency. In the coordinate systemO-xymoving at constant speedUwith a ship, the measured unsteady waves along a straight line parallel to the advancing direction of ship aty=y0can be expressed in the following form:

The relation between Fourier transformed wave elevation and the Kochin function is

whereεk= sgn(ωe+kU) and “sgn ” is the sign function. Thus, the added resistance in waves finally can be obtained in terms of the Fourier transformed ship-generated unsteady wave as follows

Because the range of measuring the wave elevation was finite, the following approximation is applied to calculate the Fourier transformed wave elevation. If the variation of wave elevations between adjacent points is linear, then the Fourier transform of the measured unsteady wave can be obtained as follows:

Here,aandbindicate the actual range of wave measurements and the number of total measured points isM+1.

2. Results and discussions

2.1 Test case

In the present study, the two modified Wigley hulls, which are called “slender” and “blunt”Wigley[13], are considered and the hull shape is expressed mathematically as follows:

Table 1 Principal particular of both Wigley hulls

The Sample of triangular surface meshes of two Wigley models are shown in Fig. 2.

Fig. 2 Sample of triangular surface meshes of two Wigley models

2.2 Sensitivity studies

The added resistance in waves is usually sensitive to the resolution of the computational grid[20].To investigate the uncertainty of the grid size for wave-induced ship motions and added resistance, the discretization error for each quantity was calculated according to the procedure suggested by Celik et al.[21].Through the use of three different grids-subscript “1”is the finest grid, “3” is coarse grid, and “2” is in between them-the grid refinement factorsr21andr32and the numerical solutions in each gridQiwere calculated. In this study, the number of mesh points was varied from two million to five million.Based on the Richardson extrapolation, the apparent orderqof the method can be calculated as follows:

The grid convergence indices for ship motions and added resistance of slender and blunt Wigley hulls with two different wavelengths,λ/L=0.5, 1.1 are summarized from Table 2 to Table 5. Here,RAW′indicates the normalized added resistance value. The grid uncertainty in the added resistance at short wave(λ/L=0.5) is larger than that at longer wave case(λ/L=1.1) for both ships. It was shown in the previous study[20]that the added resistance of a ship in short wave was quite sensitive to the computational grid. Moreover, the value of added resistance in short waves is smaller than the peak value and thus it is difficult to measure accurately the added resistance in short waves even in experiment. The grid uncertainty of heave motion amplitude is larger than that of pitch motion amplitude because the heave added mass is sensitive to the grid resolution in the present method[22]. The results of grid uncertainty analysis are used to plot the error bar in the next section.

Table 2 Calculation results of discretization error (λ/L=0.5, slender Wigley)

Table 3 Calculation results of discretization error (λ/L=1.1, slender Wigley)

Table 4 Calculation results of discretization error (λ/L=0.5, blunt Wigley)

Table 5 Calculation results of discretization error (λ/L=1.1, blunt Wigley)

The sensitivity of the lateral distance of wave measuring position was also investigated and the results are shown in Fig. 3(b). The measuring locations are indicated in Fig. 3(a) forλ/L=0.5, 1.1 cases and the closest location is the same as that in the experiment[13]. The added resistance calculated by the unsteady wave analysis decreases as the lateral distance of measuring position increases. The nonuniform grid was applied in the present computation because it is practically impossible to use the uniform grid spacing for three-dimensional computations.Thus, the numerical dissipation inevitably occurs as the wave propagates far from the ship. Consequently,it attenuates the wave elevation and hence the added resistance from the unsteady wave analysis becomes smaller in the far field.

Fig. 3 (Color online) Sensitivity of lateral distance of wave measuring position

2.3 Motion responses and added resistance in waves

Based on the previous sensitivity studies,computations with fine grid and at the closest wave measuring position were used for subsequent study,and all computational results are compared with experiment and other computed results by enhanced unified theory (EUT) and rankine panel method (RPM)given in Kashiwagi[13]. The motion amplitudes and phase angles are summarized in Figs. 4 and 5 for slender and blunt Wigley hulls, respectively. Overall results are similar to those in the experiment and the grid uncertainty values calculated in the previous section are plotted as an error bar and they are almost negligible in the motion responses.

Fig. 4 (Color online) Vertical motions of slender Wigley model:Amplitude (upper) and phase angle (lower)

The added resistance of both hulls in waves is shown in Fig. 6. The same symbol and line style as shown in the motion responses are used for the values obtained by direct measurement and pressure integration. Additionally, the results of added resistance obtained by the unsteady wave analysis are plotted with inverse triangle for the measured wave and with closed triangle for the computed wave. Both computed values by the pressure integration and the wave analysis are similar to corresponding values in the experiment. The magnitude of added resistance by the unsteady wave analysis is smaller than that by the pressure integration near the ship resonant wavelength as shown in experiment[13], while both of pressure integration and unsteady wave analysis provide similar added resistance in short wavelength region.

Fig. 5 (Color online) Vertical motions of blunt Wigley model:Amplitude (upper) and phase angle (lower)

Fig. 6 (Color online) Added resistance of two Wigley models in waves

If the ship motion becomes large, the ship-generated unsteady wave is easily broken and separation flows are generated from the bottom and shoulder parts of a ship. Due to these breaking waves,the wave energy will be dissipated. Furthermore, it was pointed out in the previous study[13]that there exists an interaction between steady and unsteady waves. The added resistance by the unsteady wave analysis cannot reflect such nonlinear phenomena with much energy dissipation. Thus, the unsteady wave analysis provides underestimated added resistance in ship resonant wavelength. It should be noted that the grid uncertainty in the added resistance by pressure integration in short wavelength(λ/L=0.5) is noticeable, especially for the blunt Wigley hull.

Wave profiles of zeroth-, first-, and secondharmonic components are summarized in Figs. 7 and 8 for slender and blunt Wigley hulls, respectively. The zeroth-harmonic wave elevation is different from each other. In particular, the Kelvin wave with small wavelength is not correctly captured in the present computational results. One possible reason of this discrepancy is due to inaccurate treatment of body boundary condition, where the body boundary condition is satisfied in terms of cell volume rather than in the exact position of the ship surface. Because of the interaction between steady and unsteady waves,the zeroth-harmonic component inλ/L=1.1 case becomes different from the Kelvin wave, especially in the fore-front part (the first trough) in the measured results.

On the contrary, the overall tendency in the firstand second-harmonic components is similar between measured and computed results because both unsteady wave components are mainly affected by inertia effect rather than viscosity. This means the present numerical code is suitable to simulate inertia dominant problems, whereas it is limited to use for the steady problem. It was indicated in Kashiwagi[13]that the wave elevation in the range ofx/(L/2) from -1.0 to 1.0 has significant contribution on the added resistance in waves. Thus, the present numerical method can provide a similar value of the added resistance to that by the experimental unsteady wave analysis even though the accuracy of steady wave elevation is not enough. The second-harmonic component is ignorable in short wavelength range,while there are prominent nonlinear unsteady waves in resonance.

Fig. 7 (Color online) Wave profiles generated by slender Wigley at λ / L =0.5 (left) and λ / L =1.1 (right)

Fig. 8 (Color online) Wave profiles generated by blunt Wigley at λ / L =0.5 (left) and λ / L =1.1 (right)

The contour plots of the first- and secondharmonic unsteady waves are shown from Figs. 9 to 16 for different wavelengths (λ/L=0.5, 1.1 and 1.6).The unsteady wave contour from measured data in the experiment[13]is plotted together in Fig. 9. The measuring range in the experiment is smaller than that in the computation and the dashed box in the computational results indicates the range of wave measurement in the experiment. The pattern of wave propagation seems similar with each other, whereas the detailed wave elevation is slightly different. The unsteady waves do not propagate ahead of the ship because the parameterτis larger than 1/4 for all cases considered in this study. The overall non-dimensional value of the first-harmonic component atλ/L=1.1 is larger than that in shorter (λ/L=0.5) and longer(λ/L=1.6) wavelength cases. It can be clearly seen that the radiation waves due to heave and pitch motion propagate toward the side direction of the ship.Those components are main contribution to the added resistance in waves for resonance case[11]. The secondharmonic component is ignorable inλ/L=0.5, 1.6 cases, while the small amount of second-harmonic component can be observed in the case of ship motion resonance.

Fig. 9 (Color online) Unsteady wave contour generated by slender Wigley at λ / L =0.5, cosine (left) and sine (right) components

Fig. 10 (Color online) Unsteady wave contour generated by slender Wigley at λ / L =1.1, cosine (left) and sine (right) components

Fig. 11 (Color online) Unsteady wave contour generated by slender Wigley at λ / L =0.5, cosine (left) and sine (right) components

Fig. 12 (Color online) Unsteady wave contour generated by blunt Wigley at λ / L =0.5, cosine (left) and sine (right) components

Before studying the motion-free cases with different amplitudes of incident wave, numerical implementation was done for the diffraction problem with different amplitudes of incident wave and the radiation problem with different excitation amplitudes.The normalized wave profiles of diffraction and radiation wave components generated by the blunt Wigley are shown in Fig. 15. Both of scattering wave and pitch radiation wave show opposite phase at bow and stern, while the same phase can be found in the heave radiation wave. Most of the normalized zerothand first-harmonic components are almost identical even for different amplitudes. On the other hand, the small wavelength component is not visible as the excitation amplitude increases in forced pitch motion case. The motion responses and the added resistance at fixed wavelengthλ/L=1.1 but different values of wave steepness are shown in Figs. 16 and 17 for the slender and blunt Wigley hulls, respectively. The dashed line with inverse triangle symbol indicates the magnitude of pitch motion and the solid line with closed triangle symbol is the magnitude of heave motion. The normalized magnitude of pitch motion slightly decreases as the wave steepness increases and consequently the added resistance decreases. This indicates that the motion response is not exactly proportional to the incident-wave amplitude and the added resistance is not exactly proportional to the square of incident-wave amplitude.

Fig. 13 (Color online) Unsteady wave contour generated by blunt Wigley at λ / L =1.1, cosine (left) and sine (right) components

Fig. 14 (Color online) Unsteady wave contour generated by blunt Wigley at λ/L = 1.6, cosine (left) and sine (right) components

The wave profiles generated by both ships in motion-free case with different incident-wave amplitudes are shown in Figs. 18 and 19 as nondimensional values. The first-harmonic component of unsteady waves becomes smaller as the wave steepness increases, whereas no significant difference can be observed in the second-harmonic component.In the previous diffraction and radiation tests, the normalized first-harmonic component was almost the same even for different excitation amplitudes. This indicates that there may be a nonlinear interaction between scattered and radiation wave components in motion-free cases, which cannot be considered in the unsteady wave analysis.

The contour plots of the first-harmonic unsteady waves for the slender Wigley are shown in Fig. 20 forλ/L=1.1 and two different values of wave steepness.As shown in the previous wave profiles along a longitudinal line with fixed lateral distancey, the magnitude of first-harmonic component of normalized unsteady waves decreases as the wave steepness increases, whereas the wave pattern looks similar even for different values of wave steepness.

The sequential snapshots of free surface and velocity vector around bow region of the slender Wigley forλ/L=1.1 are shown in Figs. 21 and 22 for wave steepnessk0A=0.0286, 0.0714, respectively. The time interval in these four figures isTe/8(whereTeindicates the encounter period) and the ship’s bow is moving down in these time instants. The bottom of ship bow is not exposed to the air for smaller wave steepness, whereas the bottom of ship bow is almost out of water in larger wave steepness.In both cases, the dominant flow is generated from the ship bottom as the ship’s bow is moving down and overturning waves occur around the ship bow. The height of overturning waves atk0A= 0.0714 is obviously higher than that atk0A= 0.0286 because of larger ship motion. Occurrence of breaking waves affects the discrepancy of the added resistance between the direct pressure integration and the unsteady wave analysis. If breaking waves occur, the amplitude of unsteady waves will decrease, providing underestimated added resistance by the unsteady wave analysis.

Fig. 15 (Color online) Wave profiles generated by blunt Wigley: Zeroth-harmonic (left) and first-harmonic (right) components

Fig. 16 Motion responses and added resistance of slender Wigley for different wave steepness and λ /L=1.1

Fig. 17 Motion responses and added resistance of blunt Wigley for different wave steepness and λ /L=1.1

Fig. 18 (Color online)Wave profiles generated by slender Wigley at λ / L =1.1 for different wave amplitudes

Fig. 19 (Color online) Wave profiles generated by blunt Wigley at λ / L =1.1 for different wave amplitudes

There are two locations where the fluid velocity becomes smaller than the ship speed. One is near the ship stem and the other is a place where the overturning waves plunge into the ambient water.Because of the body boundary condition, the incoming flow near the ship stem is changing in its direction toward the shoulder of the ship and accelerated. However, the accelerated flow meets the overturning wave soon and the flow velocity decreases.This tendency can be observed more clearly for the blunt Wigley case from Figs. 23 and 24 shown respectively for different values of wave steepnessk0A= 0.0228 andk0A= 0.0685. We can observe that overturning waves are initiated from the bottom of ship bow as the ship plunges into the water and also we can see breaking waves even for small wave steepness. Reduction of the flow velocity near the ship stem can be clearly seen for the case of blunt Wigley model. Moreover, overturning waves are more prominent compared with the slender Wigley model and consequently the flow velocity is smaller near the re-entering position. Those regions of smaller flow velocity may have higher pressure than that in other regions. Thus, the added resistance will increase due to blunt bow shape and breaking waves. However,further systematic research may be required for understating more the relation between the local flow and dynamic pressure distribution for highly nonlinear free-surface flows.

Thin water films attached around ship’s bow can be seen in top images in Figs. 22 and 24, which are unrealistic and caused by numerical error. Because the body boundary condition is satisfied in terms of cell volume rather than in the exact position of the ship surface, a certain amount of water tends to be attached to the body surface and moves together with the ship.However, its effects on the hydrodynamic force can be ignored because the volume of attached water is very small compared with ship’s displacement and there is no velocity gradient across the ship surface and attached water.

3. Conclusions

Fig. 20 (Color online) Contour of the first-harmonic component of unsteady waves generated by slender Wigley at λ / L =1.1,cosine (left) and sine (right) components

Fig. 21 (Color online)Sequential snapshots of free-surface shape (left) and velocity vector (right): Slender Wigley at λ/ L =1.1, k0 A = 0.0286 with time interval T e/8

Fig. 22 (Color online)Sequential snapshots of free-surface shape (left) and velocity vector (right): Slender Wigley at λ/ L =1.1, k0 A = 0.0714 with time interval T e/8

Fig. 23 (Color online)Sequential snapshots of free-surface shape(left)and velocity vector(right):Blunt Wigley at λ / L =1.1, k0 A = 0.0228 with time interval T e/8

Fig. 24 (Color online)Sequential snapshots of free-surface shape(left) and velocity vector(right):Blunt Wigley at λ/ L =1.1, k0 A = 0.0685 with time interval T e/8

In the present study, comparisons between computational and experimental results of wave-induced ship motions and added resistance of both slender and blunt Wigley models have been made. Furthermore,the added resistance in waves was calculated by both pressure integration and unsteady wave analysis and the ship-generated unsteady waves and local flow characteristics were investigated for different values of wave steepness. Based on this study, the following conclusions can be drawn:

(1) The sensitivity studies of grid resolution showed that the grid uncertainties in the heave and pitch motions and the added resistance in motion resonance were acceptable, whereas the grid uncertainty in the added resistance in short waves was noticeable.

(2) The added resistance computed from the unsteady wave analysis using longitudinal wave cut decreased as the lateral distance of measuring position increased. The numerical dissipation exists inevitably in the non-uniform grid spacing. Consequently, it attenuated the wave elevation and the added resistance computed from the unsteady wave analysis becomes smaller in the far field.

(3) The added resistance by the unsteady wave analysis using computed waves showed similar tendency with the experiment; that is, underestimated added resistance was obtained in the resonance region because of nonlinear phenomena, while both results by the pressure integration and the unsteady wave analysis provided similar results in short wavelength case.

(4) The first-harmonic unsteady waves normalized with incident-wave amplitude for the motionfree case decreased as the wave steepness increased,whereas independent diffraction and radiation tests showed almost similar non-dimensional values of unsteady waves even for different amplitudes of incident wave or forced excitation. Thus, the nonlinear interaction among diffraction and radiation waves as well as steady waves should be taken into account to improve the unsteady wave analysis for calculating the added resistance in waves.

(5) Detailed analyses of computed flow around the ship bow revealed that the overturning waves started to occur around the bottom of ship bow when the ship was submerged into the water even for the case of small wave steepness considered in this study.Due to breaking waves, the accuracy of unsteady wave analysis was degenerated and the local flow velocity decreased. Further systematic research is required for a deeper understating of local flow characteristics of highly nonlinear free-surface flows and their interaction with a floating body.

Acknowledgements

This study was conducted as a part of the research in the promotion program for international collaboration supported by Osaka University. The first and third authors were partly supported by the Ministry of Trade, Industry and Energy (MOTIE),Korea, through the project “Technology Development to Improve Added Resistance and Ship Operational Efficiency for Hull Form Design” (Grant No.10062881) and the Lloyd’s Register Foundation(LRF)-Funded Research Center at Seoul National University. Administrative support was also received from RIMSE and ERI at the Seoul National University. Authors appreciate all of the support received.