Effect of wave periods on added resistance and motions of a ship in head sea simulations

Ocean Engineering 137 (2017) 309–327

Contents lists available at ScienceDirect

Ocean Engineering journal homepage: www.elsevier.com/locate/oceaneng

Effect of wave periods on added resistance and motions of a ship in head sea simulations

MARK

⁎

Seonguk Seoa, Sunho Parka, , BonYong Koob a b

Department of Ocean Engineering, Korea Maritime and Ocean University, Busan, Republic of Korea School of Mechanical Convergence System Engineering, Kunsan National University, Gunsan, Jeollabuk-do, Republic of Korea

A R T I C L E I N F O

A BS T RAC T

Keywords: Added resistance Computational fluid dynamics Head wave OpenFOAM Ship motion Uncertainty assessment

Given that the International Maritime Organization (IMO) recently introduced the Energy Efficiency Design Index (EEDI) for new ships and the Energy Efficiency Operational Indicator (EEOI) for ship operations, accurate prediction of the added resistance of ships advancing in waves has become necessary. In the present study, the open-source computational fluid dynamics library, OpenFOAM, was used to predict the added resistance with the motions of a ship. Unstructured grid using a hanging-node and cut-cell method was used to generate fine grid around a free-surface and ship. Five wavelengths from a short wavelength (λ/LPP=0.65) to a long wavelength (λ/LPP=1.95) were considered as head waves. The grid uncertainty was evaluated with respect to three grid systems for a wavelength of λ/LPP=1.15. The predicted added resistance and heave and pitch motions with various head waves were compared with the data from model experiments. The wave pattern around the ship and the motions of the ship were analyzed with respect to the wave encounter period. When the wavelength was similar to the ship length, the ship moved against the waves, and thus the added resistance was greater compared with other wavelengths.

1. Introduction With the current international interest in reducing greenhouse gas emissions, efforts to enhance the eco-friendliness of ships have recently been drawing more attention. To regulate the quantity of greenhouse gases emitted to the atmosphere from ships, the International Maritime Organization (IMO) introduced the Energy Efficiency Design Index (EEDI) for new ship designs and the Energy Efficiency Operational Indicator (EEOI) for ships already in operation. To accurately estimate fuel consumption, the added resistance due to waves and winds is considered in the basic design stage to meet the requirements of the EEDI and EEOI. Generally, the total resistance of ships in waves is increased by about 15–30% in comparison with the resistance of ships in calm water (Arribas, 2007). Therefore, an accurate prediction of the added resistance has been issued. Many studies on the added resistance of ships have been dependent on experiments. Strom-Tejsen et al. (1973) conducted an added resistance experiment with the Series 60 hull form, while both Fujii and Takahashi (1975) and Nakamura and Naito (1977) carried out experiments using the S175 container ships. In addition, Journee (1992) conducted systematic added resistance experiments with four Wigley hull forms. Park et al. (2015) suggested a systematic uncertainty

⁎

Corresponding author. E-mail address: [email protected] (S. Park).

http://dx.doi.org/10.1016/j.oceaneng.2017.04.009 Received 26 October 2016; Received in revised form 21 March 2017; Accepted 6 April 2017 0029-8018/ © 2017 Elsevier Ltd. All rights reserved.

study that enabled the reliability of the experimental results of the added resistance to be investigated through repeated model experiments under various conditions. Numerical studies based on the potential flow have been conducted to estimate the added resistance. Studies have been continuously conducted using the far-field method (Gerritsma and Beukelman, 1972), the near-field method (Faltinsen et al., 1980), and the Rankine panel method (Joncquez et al., 2008; Kim and Kim, 2011). Seo et al. (2013) compared three different numerical methods (strip method, Rankine panel method, and Cartesian grid method) of calculating the added resistance of ships in waves. As numerical studies based on computational fluid dynamics (CFD), which consider the viscosity directly, have made rapid progress, extending their applications, CFD-based studies are actively being conducted with regard to the added resistance of ships. Sato et al. (1999) investigated heave and pitch motions of the Wigley hull forms and the Series 60 hull forms in waves by applying the marker-densityfunction method. Orihara and Miyata (2003) employed WISDAM-X software and an overset mesh system to analyze the added resistance of the SR-108 ships. Simonsen et al. (2013) compared the CFD results obtained by the URANS codes and the model experiment results with respect to the KCS. Similarly, Sadat-Hosseini et al. (2013) compared

Ocean Engineering 137 (2017) 309–327

S. Seo et al.

UGC Corrected uncertainty for grid size t Time Te Wave encounter period up Moving velocity of mesh point u, v, w Velocity in each direction u̇,v̇,ẇ Acceleration in each direction Vi The volume of i-th cell v Velocity vr Artificial velocity xG, yG, zG Distance between center of rotation and gravity in each direction xp Mesh point at finally deformed position over time xp0 Mesh point at previous position α Volume fraction γ Diffusion coefficient γI Initial phase of wave δb Vertical moving distance of bow with reference to design draft δs Vertical moving distance of stern with reference to design draft δ* Correction numerical error ζln Incident wave amplitude corresponding to n-th harmonic amplitude λ Wavelength μ Viscosity coefficient ρ Density τ Turbulent stress tensor

Nomenclature BWL C Cad d F fe Fr Fx,calm Fx,wave g I k LPP M p PG PG,est p, q, r ṗ,q̇,r ̇ Ri ri S Sˆn Uship UG

Maximum breadth at waterline Correction factor Anti-diffusion coefficient Distance between a cell center of a moving mesh and a center of an arbitrary cell Force Wave encounter frequency Froude number Force in calm water in x-direction Force in waves in x-direction Gravity Inertial mass Wave number Ship length between perpendiculars Moment Static pressure Order of accuracy Limiting order of accuracy for the input parameter Angular velocity in each direction Angular acceleration in each direction Convergence ratio Refinement ratio for the input parameter Source term Numerical solution for n-th grid sizes Ship speed Uncertainty for grid size

naoe-FOAM-SJTU and predicted added resistance at various wave stiffnesses. Shen et al. (2015) implemented dynamic overset grid technique into OpenFOAM and applied this to a propeller rotation for KCS self-propulsion and maneuvering. Seo et al. (2016) reviewed numerical uncertainty for grid size, time step, and iteration number and identified the OpenFOAM uncertainty for a ship's resistance performance. Vukcevic et al. (2016a, 2016b) proposed a fluid dynamic model by combining non-viscous and viscous analysis codes. In the present study, the open-source CFD libraries were used to predict the added resistance and motions of the KCS with regular head waves. To predict free-surface flow and ship motions, interDyMFoam, one of the basic application solvers of OpenFOAM, was used. Then waves2Foam was incorporated into this solver to simulate waves (Jacobsen et al., 2012). Since the encounter frequency for a wave generation did not include the linear velocity (Uship) term, as shown in Eq. (1), the term was added to waves2Foam. For a numerical validation of the result, the uncertainty of the grids was evaluated with respect to a medium wavelength (λ/LPP=1.15). To analyze the added resistance of a ship, the resistance performance of the ship in calm water was first analyzed, and then the added resistance with heave and pitch motion performances were studied with respect to five wavelengths, from a

and verified the results from a CFD method and a model experiment data with regard to the KVLCC2 ships in the cases where pitching was restricted and not restricted. Tezdogan et al. (2015) compared the CFD, the potential flow-based method, and the model experiment results with respect to full-scale KCS in low-speed sailing, and then analyzed the differences. Yet despite many numerical studies on added resistance, the prediction of added resistance is still difficult. At the 2010 Gothenburg Workshop (Larsson et al., 2013) and the 2015 Tokyo Workshop (Larsson et al., 2015), multiple researchers calculated the added resistance of a specified ship, and their numerical methods, grids, and results were compared. Recently, open-source CFD libraries, called OpenFOAM, have been used in many shipbuilding and maritime industries. In particular, as an example of a ship resistance study using OpenFOAM, Park et al. (2013) developed SNUFOAM, which applies open source libraries to shipbuilding and marine engineering. Shen and Wan (2013) developed Table 1 Principal particulars of model scale KCS. Principal particulars

Symbol

Scale

Scale factor Length between perpendiculars Maximum beam of waterline Depth Design draft Displacement volume Ship Wetted surface with rudder Longitudinal center of buoyancy (fwd+) Longitudinal center of gravity (from aft) Vertical center of gravity (from keel) Moment of inertia

λ (-) LPP (m) BWL (m) D (m) T (m) ∇ (m3) S (m2) LCB (%LPP) LCG (m) KG (m) Kxx/BWL (-) Kyy/LPP (-) Kzz/LPP (-) Re (-) Fr (-)

37.9 6.0702 0.8498 0.5015 0.2850 0.9571 6.6978 −1.48 2.945 0.192 0.40 0.25 0.25 1.07×107 0.261

Reynolds number Froude number

Table 2 Simulation conditions.

310

Case number

Froude number

Wave steepness (H/λ)

Wavelength (λ/LPP)

Encounter frequency (Hz)

Encounter period (s)

0 (calm water) 1 2 3 4 5

0.261

–

–

–

–

0.261 0.261 0.261 0.261 0.261

1/60 1/60 1/60 1/60 1/60

0.65 0.85 1.15 1.37 1.95

1.136 0.940 0.762 0.676 0.533

0.878 1.064 1.312 1.480 1.874

Ocean Engineering 137 (2017) 309–327

S. Seo et al.

Fig. 1. Domain extent, boundary conditions, and grid.

short wavelength (λ/LPP=0.65) to a long wavelength (λ/LPP=1.95). Aided by the latest CFD method, the presented paper visualized vortices structures around the stern for various wavelengths and ship motions. The vortices structures can be used for stern shape, propeller, and rudder design. In particular, the variation of the vortices structure due to the pitch motion under the transom can be used as a valuable reference for understanding the stern slamming. This paper is organized as follows. The next section describes the problem and discusses the simulation conditions. Next, the computational methods are described, followed by the analysis method. Finally, some conclusions are derived.

fe =

g /(2πλ ) +Uship /λ

(1)

Te=1/ fe

(2)

Here, g is the acceleration of the gravity, and Uship is the ship speed. The x-directional resistance coefficient acting on the ship in the calm water was calculated as

CT =

Fx, calm 2 0. 5ρAUship

(3)

where Fx,calm is the force acting on the ship in calm water in xdirection. ρ and A are the density and wetted surface area of the ship, respectively. The added resistance of the ship caused by the waves was calculated by subtracting the resistance in calm water from the resistance in the waves. The added resistance coefficient of the ship in the waves was expressed as

2. Model description and simulation conditions The subject ship of the present study was a 3600 TEU KRISO container ship (KCS) designed by the Korea Research Institute of Ship and Ocean Engineering (KRISO). The length was established as 6.07 m, which was the same as that of the KCS used in the model experiment conducted by the FORCE Technology (Sadat-Hosseini et al., 2015). The model of the KCS included the rudder. Table 1 lists the principal particulars of the KCS. As listed in Table 2, the simulation conditions include calm water (Case 0), short wavelength (Cases 1 and 2), medium wavelength (Case 3), and long wavelength (Cases 4 and 5). The heave and pitch motions of the ship were free, and other motions were restricted. The wave steepness (H/λ) was fixed at 1/60 in all wavelength cases, and the water depth was assumed to be infinite. The encounter frequency (fe) and encounter period (Te) that were generated due to the advancing velocity of the ship (Uship) and regular waves were calculated as

σaw =

(Fx, wave − Fx, calm ) 2 ρgζl21 BWL / LPP

(4)

where Fx,wave is the x-directional force in waves corresponding to the 0-th harmonic amplitude, ζI1 is the wave amplitude corresponding to the first harmonic amplitude, and BWL is the maximum breadth at the water line. To compare the motion responses with experimental data, the heave and pitch motions are expressed, respectively, as

Nondimensionalized heave motion=

z A

(5)

Nondimensionalized pitch motion=

Θ kA

(6)

Fig. 2. Simulation domain for wave generating and damping zones.

311

Ocean Engineering 137 (2017) 309–327

S. Seo et al.

Fig. 3. y+ contours on ship surface for various wavelengths.

312

Ocean Engineering 137 (2017) 309–327

S. Seo et al.

Table 3 Total resistance coefficient and heave and pitch motions in calm water.

Present (Case 0) EFD (Sadat-Hosseini et al., 2015) Difference (% of EFD)

Table 4 1st harmonic wave amplitude for various wavelengths.

CT (×10−3)

sinkage/ LPP (×10−3)

trim (deg.)

4.018 3.835

−1.989 −2.074

−0.1676 −0.1646

−4.757

4.216

−1.845

Wave amplitude (ζ1/LPP×103) Case λ/LPP Present EFD (Sadat-Hosseini et al., 2015) Difference (% of EFD)

where z is the heave displacement, Θ is the pitch angle and k is the wave number.

2 0.85 6.226 6.409 2.849

3 1.15 10.033 10.152 1.167

4 1.37 12.169 12.300 1.068

5 1.95 16.111 16.114 0.044

Table 5 Grid sizes for uncertainty analysis.

3. Computational methods 3.1. Governing equations The mass and momentum conservation equations were considered to obtain the velocity and pressure fields. The mass and momentum conservation equations can be written, respectively, as

∇∙v⎯→ m =0

1 0.65 4.883 5.132 4.846

Grid system

Number of background of mesh

Number of mesh

Coarse (N3) Medium (N2) Fine (N1)

128×16×48 160×20×60 200×25×75

2,026,386 3,700,684 6,877,501

indicates the mixture phase. The mixture properties of the density and viscosity are calculated as

(7)

∂ρm v⎯→ m → +∇∙(ρm v⎯→⎯ m vm )=−∇p +∇∙(τ )+ S (8) ∂t → ⎯ where v , p, τ and S are the velocity vector, static pressure, turbulent stress tensor and the source term, respectively. The subscript m

ρm = αa ρa +α w ρw

(9)

μm = αa μa +α w μw

(10)

where ρ, μ, and α are the density, viscosity coefficient, and volume fraction, respectively. Subscripts a and w indicate the air and water

Fig. 4. Time history of wave elevation for various wavelengths (EFD: symbol, Present: line).

313

Ocean Engineering 137 (2017) 309–327

S. Seo et al.

Fig. 5. One encounter period of wave amplitude, added resistance coefficient, heave and pitch motion for grid convergence study (Case 3, λ/LPP=1.15).

Table 6 Grid uncertainty analysis for Case 3. Solution

RG a

PG

C

δ*

UG

UGC

−0.29×10−3 (−0.48%) 1.97 (20.15%) −2.05×10−2 (−2.25%) −3.94×10−2 (−5.38%)

0.48×10−3 (0.78%) 3.13 (31.89%) 2.67×10−2 (2.92%) 5.84×10−2 (7.97%)

0.18×10−3 (0.30%) 1.15 (11.74%) 0.62×10−2 (0.68%) 1.90×10−2 (2.60%)

S3

S2

S1

EFD

ζI1

0.0605

0.0609

0.0611

0.0616

MC

4.22

2.72

σaw

12.188

10.807

9.797

9.911

MC

1.36

0.63

z/A

0.882

0.901

0.911

0.899

MC

2.66

1.43

Θ /Ak

0.713

0.724

0.733

0.748

MC

0.70

0.67

a

MC: Monotonic Convergence

the y-axis, and the upward direction was the z-axis. Fig. 1 shows the domain extent, boundary conditions, and grid. The length measured from the bow to the inlet boundary was set as 2LPP, and the length measured from the stern to the outlet boundary was set as 5LPP. The width measured from the ship to the side boundary was set as 1LPP. The height above the water level was set as 0.5LPP, and the depth below the water level was set as 3LPP (Park et al., 2017). Here, LPP was the ship length. On the inlet, side, bottom and top boundaries, the Dirichlet condition was applied for the velocity, turbulence properties and volume fraction, while the Neumann condition was applied for the pressure. On the outlet boundary, the Neumann condition was applied for the velocity, turbulence properties and volume fraction, while the Dirichlet condition was applied for the pressure, which generated the hydrostatic pressure (Park et al., 2015). The symmetry condition was applied to the plane, y=0, because the model was physically symmetric to the plane. In general, the no-slip condition was applied to the ship surface, but because the ship showed motions due to the waves, the

phases, respectively. The volume fraction transport equation was considered to get the free-surface wave flow, which can be expressed as

∂ (αρm )+∇∙(αρm v⎯→ m) = ∂t

−Cad ∇∙(α (1 − α )→ v⎯r )

(11)

where the volume fraction is 1 in the water phase, and 0 in the air phase. The term on the right-hand side is an anti-diffusion source term, which reduced a solution smearing (Lee and Rhee, 2015; Park et al., 2016). Cad was an anti-diffusion coefficient. Here, Cad of 1 was used. → v⎯r was the artificial velocity to reduce interface smearing near the interface. 3.2. Numerical descriptions The Cartesian coordinates system was used in the computations. With reference to the ship's center of gravity, the direction toward the stern of the ship was set as the x-axis, the direction toward the port was 314

Ocean Engineering 137 (2017) 309–327

S. Seo et al.

mean velocity, turbulence kinetic energy, and turbulence dissipation rate at wall boundaries. The time step size of 0.001 s and 0.002 s showed similar results in terms of the wave propagation, while the time step size of 0.004 s causes different results. Thus, the time step size of 0.002 s was selected. For the selected time step size, the maximum and minimum Courant numbers for any meshes in the computational domain were 1.870 and 0.002, respectively, for case 3. In the present study, 439, 656, and 937 time steps per encounter period were used to predict ship responses by incident regular waves for cases 1, 3, and 5, respectively. These time steps per encounter period satisfy the ITTC recommendation, which was 100 time steps per encounter period (ITTC, 2011).

moving speed of the ship was applied to the ship surface as the boundary condition. Waves were generated in a wave-generating zone and dissipated in a damping zone to prevent reflection at the boundaries. The established wave-generating and damping zones are shown in Fig. 2, in which the waves moved from the left to the right. The x-directional size of the wave generating and damping zones was set to be 2λ for the short wavelength and 1λ for the medium and long wavelengths. The incidence wave was generated using the second-order Stokes’ wave theory (Stokes, 1847). As the initial condition, the waves were set in the whole domain and the wave crest was located at the bow of the ship (x/ λ=0). The center of gravity of the ship was at x/λ=0.513. For grid generation, automatic grid generation utilities provided by OpenFOAM called blockMesh and snappyHexMesh were used. The blockMesh utility was used to generate structured grids for the background mesh on the whole domain. Then the snappyHexMesh utility was used to generate unstructured grids, including hanging node and cut cell, around the ship and free-surface. The grids were refined around the ship and free-surface, while the density of the grids was coarse in the regions closer to the outlet, top, and bottom boundaries. The total number of grids was about 3.7×106, and about 3×104 grids were used for the surface of the ship. This number is a reference number for the medium grids in the grid dependency test that will be discussed later. Fig. 3 shows the y+ contours of the ship surface. The mean y+ values of cases 1, 2, 3, 4 and 5 are 171, 172, 173, 178 and 171 when the wave crest passes the bow. Because these y+ values are included in the overlap layer, the selected log-law based wall function (Park et al., 2013) is valid to calculate the solution variables including

3.3. Numerical methods Time derivative terms were discretized using the first-order accurate backward implicit scheme, which has proved sufficient for engineering accuracy. The solution gradients at the cell centers were evaluated using the Euler method. The convection terms were discretized using a total variation diminishing (TVD) scheme with the vanLeer limiter (van Leer, 1979), and for the diffusion terms, a central differencing scheme was used. The velocity-pressure coupling and overall solution procedure were based on the PIMPLE algorithm (Vukcevic et al., 2016a), which combines the SIMPLE (Patankar and Spalding, 1972) and PISO (Issa, 1985) algorithms. SST k-ω turbulence model (Menter, 1993) with the wall function (Park et al., 2013) was used for the turbulence closure. A higher-order interface capturing scheme, which combined the upwind and downwind schemes, was

Fig. 6. Time history for CT, heave and pitch motions in waves (EFD: symbol, Present: line).

315

Ocean Engineering 137 (2017) 309–327

S. Seo et al.

Fig. 6. (continued)

and the sum of the force and moment acting in each direction was calculated as

used to capture the free-surface behavior. The discretized algebraic equations were solved using a pointwise Gauss-Seidel iterative algorithm, while an algebraic multi-grid method (Weiss et al., 1999) was employed to accelerate solution convergence.

∑ Fx = m [u−̇ vr +wq−xG (q 2+r 2)+yG (pg−r )̇ +zG (pr +q )̇ ]

3.4. Mesh deformation

∑ Fy = m [v−̇ wp+ur−yG (r 2+p2 )+zG (qr−p ̇)+xG (qp+r )̇ ]

Meshes deformed to express motions of a ship resulting from wave loads. In the whole domain, the meshes were deformed by solving the Laplace equation expressed as

∇∙(γ ∇u⎯→ p )=0

∑ Fz = m [ẇ −uq+vp−zG ( p2+q 2)+xG (rp−q )̇ +yG (rq+p ̇) ]

(12)

∑ Mx = Ix p ̇ +(Iz−Iy) qr +m [ yG (ẇ −uq+vp)−zG (v−̇ wp+ur )]

where γ denotes the diffusion coefficient. The diffusion coefficient was calculated from γ = 1/ d 2 , which was a second-order accurate method. Here, d denotes the distance between a cell center of a moving mesh and a center of an arbitrary cell, and u⎯→ p denotes the moving velocity of the mesh. The moving velocity of the mesh may be adjusted according to the magnitude of the diffusion coefficient. After obtaining the moving velocity of the mesh, the internal mesh were deformed using the following equation

∑ My = Iy q + ̇ (Ix −Iz ) rp + m [z G (u − ̇ vr +wq ) −xG (ẇ −uq +vp)] ∑ Mz = Iz r + ̇ (Iy−Ix ) pq+m [xG (v− ̇ wp +ur )−yG (u− ̇ vr +wq )]

⎯x→ = ⎯x⎯⎯→+u⎯→∆t p p0 p

(14)

(13) ⎯ ⎯⎯ → ⎯ → where xp0 denotes the mesh at the previous position, and xp denotes the mesh at the deformed position.

where F , M, and I are the force, moment, and inertial mass, respectively. In addition, u, v, and w, and p, q, and r are the velocity and angular velocity in the directions of x, y, and z, respectively. u̇, v̇, and ẇ , and ṗ, q̇, and r ̇ are the acceleration and angular acceleration in the directions of x, y, and z, respectively. xG, yG, and zG are the distance between the center of the rotation and the center of the gravity in the directions of x, y, and z, respectively.

3.5. Six Degrees of Freedom Motion To analyze the six degrees of freedom motion of a ship by waves, the ship was considered as a rigid body in the Cartesian coordinate system, 316

Ocean Engineering 137 (2017) 309–327

S. Seo et al.

Fig. 6. (continued)

4. Results and discussion

4.2. Wave generation

4.1. Calm water

The waves were generated using the second-order Stokes’ wave theory (Stokes, 1847). Accurate and consistent generation of all waves was an important factor for the analysis of the added resistance. To validate the accurate and consistent generation of the waves in the domain, the simulated wave propagations over time at x/λ=0 were plotted and compared with the experimental data (http://www.t2015. nmri.go.jp). The wave profile was captured using the volume fraction value of 0.5. The wave propagation studies were carried out without the ship. These studies were carried out with same mesh, and thus the mesh count was not changed with the wavelength. Thus, 104 and 312 grids for one wave period for cases 1 and 5 were used, which was more than the ITTC-recommended number of 80 (ITTC, 2011). To verify the general tendency of wave propagation, the two wave periods are considered and compared with the experimental results, as shown in Fig. 4. The overall wave propagation was similar to that of the experimental data. The first harmonic amplitude for the generated wave is compared with the experimental data listed in Table 4. Since

Table 3 lists the total resistance acting on the ship, sinkage, and trim in the calm water (Case 0) with the Froude number of 0.261. A negative value for the sinkage and trim motions indicated that the ship descended and the bow of the ship inclined downward, respectively. The result showed that, when the KCS was moving forward in the calm water, the ship descended and the bow was inclined downward. Here the trim angle was calculated by using the vertical moving distance of the bow (δb) and the stern (δs) with reference to the design draft, which was tan−1 (2 δb+δs / LBP ). The comparison with the experimental data showed that the difference between the total resistance and sinkage was less than 5% and the difference for the trim was less than 2%, indicating that the computational results were close to the experimental results (Sadat-Hosseini et al., 2015).

317

Ocean Engineering 137 (2017) 309–327

S. Seo et al.

Fig. 7. Time history for total resistance, heave and pitch motions over 3 wave periods of Case 3 (λ/LPP=1.15). Table 7 CT and σaw for various wavelengths. σaw

CT (×103) Case λ/LPP Present EFD Difference (% of EFD)

1 0.65 4.38 4.13 −6.12

2 0.85 5.04 4.63 −8.94

3 1.15 7.45 7.08 −5.27

4 1.37 7.56 6.98 −8.34

5 1.95 5.51 5.42 −1.73

1 0.65 4.80 3.39 −41.85

2 0.85 8.32 6.10 −36.35

3 1.15 10.81 9.91 −9.04

4 1.37 7.58 6.51 −16.34

5 1.95 1.99 1.91 −4.04

4.3. Uncertainty assessment of added resistance

the same mesh was used for all cases, more meshes were included in one wavelength as the wavelength was increased, and thus the difference was decreased. Overall, the differences in comparison with the experimental data were less than 5%, indicating that the boundary conditions, domain extent, and numerical methods for wave generation were adequately selected.

Numerical uncertainty of the added resistance was assessed for three types of grids. As listed in Table 5, three types of grid size, i.e., coarse, medium, and fine grids, were considered. The background structured grid size was increased or decreased in all directions by 1.25

318

Ocean Engineering 137 (2017) 309–327

S. Seo et al.

Table 8 Heave and pitch motions for various wavelengths. z/A Case λ/LPP Present EFD Difference (% of EFD)

Θ /Ak

1 0.65 0.12 0.13 9.94

2 0.85 0.24 0.24 −0.45

3 1.15 0.90 0.90 −0.24

4 1.37 0.90 0.88 −2.76

5 1.95 0.89 0.93 4.62

Table 9 Phase difference of added resistance, heave and pitch motions for various wavelengths. Phase difference (%)

Case 1

Case 2

Case 3

Case 4

Case 5

Added resistance Heave motion pitch motion

17.04 195.21 72.87

62.44 205.16 154.14

– 101.60 123.24

5.33 130.40 180.95

2.52 186.70 13.23

N

i =1

⎠

1 ln (Sˆ3−Sˆ2 )/(Sˆ2−Sˆ1) +q (PG ) ln (r21)

⎛ r PG − sign ((Sˆ3 − Sˆ2 )/(Sˆ2 − Sˆ1)) ⎞ ⎟⎟ q (PG ) = ln ⎜⎜ 32 P ⎝ r21G − sign ((Sˆ3 − Sˆ2 )/(Sˆ2 − Sˆ1)) ⎠

(15)

(16)

(17)

where Sˆ1, Sˆ2 , and Sˆ3 are numerical solutions for the three grid sizes, respectively. The value of 0 < (Sˆ3−Sˆ2 )/(Sˆ2−Sˆ1) < 1 indicated monotonic convergence. In addition, the uncertainty (UG), corrected numerical error (δ*), and corrected uncertainty (UGC) can be defined by applying the correction factor (C) as follows

C=

PG r21 −1 P

r21G, est −1

4 1.37 0.97 0.97 −0.58

5 1.95 1.07 1.12 4.65

⎛ˆ ˆ⎞ * = C ⎜⎜ S2 − S1 ⎟⎟ δ* = CδRE P ⎝ r21G−1 ⎠

(19)

* + (1 − C ) δRE * UG = CδRE

(20)

* UGC = (1 − C ) δRE

(21)

Fig. 6 shows the total resistance coefficient (CT), and heave (z/A) and pitch (Θ/Ak) motion responses in the waves for one wave period compared with the existing experimental data. Fig. 7 shows the three encounter periods after the solutions have been converged. Fig. 6 is the selected last encounter period from Fig. 7. The experimental data were obtained from the model experiment with 6.07 m KCS provided by the FORCE Technology (http://www.t2015.nmri.go.jp). In the short wavelength (Cases 1 and 2), the pitch motion response was slightly different from the experimental data. In particular, in case 2, the total resistance coefficient showed a different pattern from that seen in the experimental results. Some numerical results presented in the 2015 Tokyo CFD Workshop also showed different results compared with the experimental data (Stern, 2015). Potential reasons for the difference might be the model test uncertainty for the generating short waves, different initial conditions between the model experiment and the present study, and reflected waves. The heave and pitch motion responses in the medium and long wavelength cases (Cases 3, 4, and 5) were relatively consistent with the experimental data. The total resistance of the model experiment for case 3 was shown as an average value over the period because the total resistance had great variation. Similarly, the present total resistance also included irregular variation. The total resistance for case 4 showed a similar trend to the result of the model experiment, although the some oscillation was shown. The resistance for case 5 was similar to the experimental result. Tables 7, 8 provide quantitative comparisons of the total resistance coefficient, added resistance coefficient (σaw), and heave (z/A) and pitch (Θ / kA) motion responses for each wavelength with the experimental results. The overall difference was greater in the added resistance and pitch motion than in the total resistance and heave motion. In addition, the difference was greater in the short wavelength cases (Cases 1 and 2) than in other wavelength cases. The difference in the total resistance calculated for each wavelength was less than 9% in

where ΔVi is the volume of the i-th cell. The refinement ratio (ri ) was expressed as the ratio between the input parameters (r21 = ∆x2 /∆x1, r32=∆x3 /∆x2 ). For r21 ≠ r32, the order of accuracy (PG) is calculated as

PG =

3 1.15 0.72 0.75 3.29

4.4. Added resistance and motion response in waves

⎞1/3

∑ (∆Vi ) ⎟⎟

2 0.85 0.23 0.15 −53.57

Table 6 lists the uncertainty assessment results of the wave amplitude, added resistance coefficient, and heave and pitch motions for the three types of grid for case 3. The limiting order of accuracy (PG,est) of 2 was used. All of the results showed monotonic convergence. The uncertainty of the wave amplitude was 0.78% and the corrected uncertainty was 0.30%, indicating that the wave amplitude was less sensitive to the change of the grid. The uncertainty of the heave motion was 2.92% and the corrected uncertainty was 0.68%, while the uncertainty of the pitch motion was 7.97% and the corrected uncertainty was 2.60%. The uncertainty of the added resistance was 31.89% and the corrected uncertainty was 11.74%, which were higher than those of the motion responses. The corrected and uncorrected uncertainty of the added resistance was higher than that obtained through the experiment, which was 7.88% (Simonsen et al., 2014).

times each. The refinement ratio of 2 satisfied the recommendation of the refinement ratio of 1.3 (Celik et al., 2008). Here, a 1.25 times increase in each direction is equivalent to about 2 times volumetric increase in the 3-dimensions. In addition, the aspect ratio of the three grid systems was set to be the same. An unstructured grid around the ship and free-surface was generated using hanging nodes and cut-cell meshes. Despite being unstructured, the total grid number increased with consistency (Seo et al., 2016). The finally generated total grid number was 2.03×106 (N3) for the coarse grid, 3.70×106 (N2) for the medium grid, and 6.88×106 (N1) for the fine grid. To carry out the grid uncertainty assessment, case 3 was selected. Fig. 5 shows one encounter period of the wave amplitude, added resistance coefficient, and motion responses for three grids. Through the uncertainty analysis for three types of grid, it was found that the wave amplitude, added resistance coefficient, and motion responses were dependent on the grid change. In particular, the ratio of the change in the added resistance coefficient was excessive. Magnification of the plot showed that the change in the magnitude was decreased as the grid size was decreased from the coarse grid to the fine grid. This tendency was shown as an indication of a monotonic convergence, and the uncertainty could be analyzed by the Richardson extrapolation. The representative cell size (Δx) was calculated as follows (Celik et al., 2008)

⎛1 ∆x = ⎜⎜ ⎝N

1 0.65 0.01 0.02 42.00

(18) 319

Ocean Engineering 137 (2017) 309–327

S. Seo et al.

Fig. 8. Added resistance coefficient, heave and pitch motions in head waves.

320

Ocean Engineering 137 (2017) 309–327

S. Seo et al.

Fig. 9. One encounter period of added resistance coefficient, heave and pitch motions for various wavelengths.

the experimental results in the short wavelength cases. Since the pitch motion had very small values in the short wavelength cases, a slight numerical deviation from the experimental result caused a great difference. Table 9 lists the phase difference of the added resistance, heave and pitch motions for various wavelengths. In terms of the added resistance, case 2 showed a large phase difference. For the heave and pitch motions, the difference was relatively large because the absolute values were sufficiently small.

comparison with the experimental result. However, the difference in the added resistance was greater because the difference of total resistance in the calm water and the difference of the first harmonic amplitude were added to the difference of the total resistance. The difference (41.85%) of the added resistance between the present finding and the experiment for the case 1 was similar to the sum of the grid uncertainty (31.89%) and the experimental uncertainty (7.88%). The difference in the pitch motion was greatly different from 321

Ocean Engineering 137 (2017) 309–327

S. Seo et al.

Fig. 10. Wave pattern (top view) and motions (side view) over one encounter period of Case 1 (λ/LPP=0.65).

Fig. 11. Wave pattern (top view) and motions (side view) over one encounter period of Case 2 (λ/LPP=0.85).

322

Ocean Engineering 137 (2017) 309–327

S. Seo et al.

the trough of the wave passed the bow, and the bow reached the highest position. In other words, the bow motion was opposite to the wave amplitude, and thus the added resistance coefficient increased. Thus, for case 3, the added resistance coefficient is higher than that of other cases, as shown in Fig. 8. In case 5, when t/Te was 0, the crest of the wave passed the bow, and the pitch had a positive peak, which was in contrast to case 3. In other words, when the ship length was longer than the wavelength, the ship's motion was adjusted to the wave elevation. Thus, for case 5, the added resistance was lower than that of other cases. The maximum level of the added resistance coefficient for various wavelengths was similar, while the minimum level of the added resistance coefficient was different. In case 3, the minimum level of the added resistance coefficient was higher than in other cases. The wave pattern around the ship and the heave and pitch motions of the ship for various wavelengths are shown in Figs. 10–14 at t/Te =0, 1/4, 2/4, and 3/4, respectively. Here, when a crest of the incident wave passed the bow of the ship, t/Te was set as 0. The ship motions were magnified as the wavelength increased. In case 3, the ship moved in opposition to the wave's elevation, while in case 5, the ship moved with the wave's elevation. Figs. 15–17 show vortices structures generated around the stern for cases 1, 3, and 5. The vortices structures are shown when the Qcriterion value is 1, and are expressed as the velocity magnitude normalized by the ship speed. In the Figures, the vortices structures are generated around the ship bottom, transom, and rudder. It can be seen that the vortices structures caused by the rudder developed to the

Fig. 8 shows the magnitude of the added resistance and the motion responses for various wavelengths. In addition to the model experiment results with the 6.07 m KCS (Sadat-Hosseini et al., 2015), which were the same as in the present study, the results obtained from the model experiment with the 4.37 m KCS (Simonsen et al., 2013) and the 2.10 m KCS (Sadat-Hosseini et al., 2015) were also plotted. The present added resistance was greater than the result from the model experiment conducted with the same length of the ship, but within the range of the results from the model experiments conducted with different ship lengths. The heave and pitch motions were almost the same as the results from the model experiment conducted with the same size of model. The results in the short wavelength cases, shown in Fig. 6, are different from the experimental results due to the phase difference, but the magnitude at each wavelength was similar to the experimental result. As such, an accurate analysis of the added resistance should include not only the magnitude at each wavelength but also the change in the magnitude over time. Fig. 9 shows the added resistance coefficient, and heave and pitch motions in one wave encounter period (Te). When the wave crest passed the bow and stern, t/Te =0 and 0.5, respectively. The pitch and heave motions increased with an increase of the wavelength. In case 1, slight heave and pitch motions were found at positions lower than the center of the ship. In case 3, the wavelength and wave amplitude had a great effect on the added resistance and motion response. When t/Te was 0, the pitch had a negative peak, and thus the bow was inclined downward, making it greatly affected by the waves. When t/Te was 0.5,

Fig. 12. Wave pattern (top view) and motions (side view) over one encounter period of Case 3 (λ/LPP=1.15).

323

Ocean Engineering 137 (2017) 309–327

S. Seo et al.

Fig. 13. Wave pattern (top view) and motions (side view) over one encounter period of Case 4 (λ/LPP=1.37).

effect on the analysis of the added resistance. A comparison of the present and model experiment results for one encounter period showed that the total resistance and motion responses had some discrepancy with the short waves, while showing a reasonable agreement with the long waves. Through the Fourier series analysis, it was found that the total resistance and motion responses were similar to the experimental data. The maximum difference of the added resistance was 41.85%, which was for case 1, and the minimum difference was 4.04%, which was for case 5. The heave motion was similar to the experimental data, while the pitch motion showed different behavior from the experimental data with the short wavelength. The added resistance of the ship was small with the short wavelengths because the absolute wave amplitude was small. For long wavelengths, the ship showed motion in accordance with the wave elevation; thus, the added resistance was small. However, for the medium waves, in which the wavelength was similar to the ship length, the ship moved against the waves, and thus the added resistance was greater compared with other wavelengths. Large vortices structures occurred under the transom and ship bottom due to the large ship motions. It is believed that the presented numerical modelling can predict the added resistance and motions of a ship. Future studies will be needed to assess the added resistance and motion responses of a ship with oblique waves.

wake and fell off. In case 1, the size of the vortices structures is relatively small compared to cases 3 and 5 because the influence of the waves is small. Since the heave and pitch motions of case 1 were small, it can be seen that the size of the vortices structures remained almost unchanged, even though the encounter period changed. In cases 3 and 5, the vortices structures varied due to the motions of the ship. When the stern was higher than the design draft (t/ Te =0.5 in case 3, t/ Te =0.25, 0.5 in case 5), the vortices structures did not occur under the transom and large amounts of vortices structures occurred on the ship bottom. On the other hand, when the stern was lower than the design draft, large vortices structures occurred under the transom, which was larger than that of case 1. Such vortices structures can be used as an understanding of the stern slamming. 5. Conclusions In the present study, numerical simulations for the added resistance and motions of the ship in the regular waves were performed using an open-source library called OpenFOAM. Five incident waves with different wavelengths were considered. The generated waves were compared with existing experimental data, and it was found that a longer wave period enabled easier prediction. The grid uncertainty was assessed. As the grid size increased, the wave amplitude, total resistance and motion responses showed a monotonic convergence. The uncertainty of the added resistance was 31.89% and the corrected uncertainty was 11.74%, indicating that the grid refinement had a great

324

Ocean Engineering 137 (2017) 309–327

S. Seo et al.

Fig. 14. Wave pattern (top view) and motions (side view) over one encounter period of Case 5 (λ/LPP=1.95).

Fig. 15. Vortices structures around stern for one encounter period of case 1 (λ/LPP=0.65).

325

Ocean Engineering 137 (2017) 309–327

S. Seo et al.

Fig. 16. Vortices structures around stern for one encounter period of case 3 (λ/LPP=1.15).

Fig. 17. Vortices structures around stern for one encounter period of case 5 (λ/LPP=1.95). resistance and propulsion of a ship in a seaway. In: Proceedings of the 13th Symposium on Naval Hydrodynamics. Tokyo, Japan, 6-10 October. Fujii, H., Takahashi, T., 1975. Experimental study on the resistance increase of a ship in regular oblique waves. In: Proceedings of the 14th ITTC, Ontario, Canada, 2-11 September. Gerritsma,, J., Beukelman,, W., 1972. Analysis of the resistance increase in waves of a fast cargo ship. Int. Shipbuild. Progress. 19 (217), 285–293. Issa,, R.I., 1985. Solution of the Implicitly discretised fluid flow equations by operatorsplitting. J. Comput. Phys. 62, 40–65. International Towing Tank Conference (ITTC), 2011. Practical guidelines for ship CFD applications. In: Proceedings of the 26th International Towing Tank Conference, Rio de Janeiro, Brazil, 28 August – 3 September. Jacobsen,, N.G., Fuhrman,, D.R., Fredsøe,, J., 2012. A wave generation toolbox for the open-source CFD library: OpenFOAM. Int. J. Numer. Methods Fluids 70 (9), 1073–1088. Joncquez, S.A.G., Bingham, H.B., Andersen, P. and Kring, D., 2008. Validation of Added

Acknowledgement This work was supported by the National Research Foundation (2015037577) of the Korea government. References Arribas,, F.P., 2007. Some methods to obtain the added resistance of a ship advancing in waves. Ocean Eng. 34, 946–955. Celik,, I.B., Ghia,, U., Roache,, P.J., Freitas,, C.J., Coleman,, H., Raad,, P.E., 2008. Procedure for estimation and reporting of uncertainty due to discretization in CFD applications. J. Fluids Eng. 130 (7), 078001. Faltinsen, O.M., Minsaas, K.J., Liapis, N. and Skjørdal, S.O., 1980. Prediction of

326

Ocean Engineering 137 (2017) 309–327

S. Seo et al.

Sadat-Hosseini,, H., Wu,, P.C., Carrica,, P.M., Kim,, H., Toda,, Y., Stern,, F., 2013. CFD verification and validation of added resistance and motions of KVLCC2 with fixed and free surge in short and long head waves. Ocean Eng. 59, 240–273. Sadat-Hosseini, H., Toxopeus, S., Kim, D.H., Castiglione, T., Sanada, Y., Stocker, M., Simonsen, C., Otzen, J.F., Toda, Y. and Stern, F., 2015. Experiments and computations for KCS added resistance for variable heading. World Maritime Technology Conference, Providence, RI, 1-5 November. Sato,, Y., Miyata,, H., Sato,, T., 1999. CFD simulation of 3-dimensional motion of a ship in waves: application to an advancing ship in regular heading waves. J. Mar. Sci. Technol. 4, 108–116. Seo,, M.-G., Park,, D.-M., Yang,, K.-K., Kim,, Y., 2013. Comparative study on computation of ship added resistance in waves. Ocean Eng. 73, 1–15. Seo,, S., Song,, S., Park,, S., 2016. A study on CFD uncertainty analysis and its application to ship resistance performance using open source libraries. J. Soc. Nav. Archit. Korea 53 (4), 338–344. Shen,, Z., Wan,, D., 2013. RANS computations of added resistance and motions of a ship in head waves. Int. J. Offshore Polar Eng. 23 (4), 263–271. Shen,, Z., Wan,, D., Carrica,, P.M., 2015. Dynamic overset grids in OpenFOAM with application to KCS self-propulsion and maneuvering. Ocean Eng. 108, 287–306. Simonsen,, C.D., Otzen,, J.F., Joncquez,, S., Stern,, F., 2013. EFD and CFD for KCS heaving and pitching in regular head wave. J. Mar. Sci. Technol. 18, 435–459. Simonsen, C.D., Otzen, J.F., Nielsen, C. and Stern, F., 2014. CFD prediction of added resistance of the KCS in regular head and oblique waves. In: Proceedings of the 30th Symposium on Naval Hydrodynamics, Hobart, Austrailia, 2-7 November 2014. Stokes,, G.G., 1847. On the theory of oscillatory waves. Trans. Camb. Philos. Soc. 8, 441–455. Storm-Tejsen,, J., Yeh,, H.Y.H., Moran,, D.D., 1973. Added resistance in waves. Soc. Nav. Archit. Mar. Eng. Trans. 81, 250–279. Tezdogan,, T., Demirel,, Y.K., Kellett,, P., Khorasanchi,, M., Incecik,, A., Turan,, O., 2015. Full-scale unsteady RANS CFD simulations of ship behaviour and performance in head seas due to slow steaming. Ocean Eng. 97, 186–206. van Leer,, B., 1979. Towards the ultimate conservative difference scheme. J. Comput. Phys. 32 (1), 101–136. Vukcevic,, V., Jasak,, H., Malenica, S., 2016a. Decomposition model for naval hydrodynamic applications, Part I: Computational method. Ocean Eng. 121, 37–46. Vukcevic,, V., Jasak,, H., Malenica, S., 2016b. Decomposition model for naval hydrodynamic applications, Part II: verification and validation. Ocean Eng. 121, 76–88. Weiss,, J.M., Maruszewski,, J.P., Smith,, W.A., 1999. Implicit solution of preconditioned Navier–Stokes equations using algebraic multigrid. AIAA J. 37 (1), 29–36.

Resistance Computations by a Potential-Flow Boundary-Element Method. In: Proceedings of the 27th Symposium on Naval Hydrodynamics. Seoul, Korea, 5-10 October 2008. Journee,, J.M.J., 1992. Experiments and Calculations on 4 Wigley Hull Forms in Head Waves. Delft University of Technology, (Report No. 0909). Stern, F., 2015. Case2.10: KCS calm water and regular head waves – time series. Tokyo 2015 Workshop on CFD in Ship Hydrodynamics. Tokyo, Japan, 2–5 December. Kim,, K.H., Kim,, Y., 2011. Numerical study on added resistance of ships by using a timedomain Rankine panel method. Ocean Eng. 38, 1357–1367. Larsson,, L., Stern,, F., Visonneau,, M., 2013. Numerical Ship Hydrodynamics - An Assessment of the Gothenburg 2010 Workshop. Springer, Dordrecht. Larsson, L., Stern, F., Visonneau, M., Hirata, N., Hino, T. and Kim, J., 2015. Tokyo 2015 Workshop on CFD in Ship Hydrodynamics. [Online] (Updated 1 August 2016) Available at: 〈http://www.t2015.nmri.go.jp〉 (accessed 18 October 2016). Lee,, H., Rhee,, S.H., 2015. A dynamic interface compression method for VOF simulations of high-speed planing watercraft. J. Mech. Sci. Technol. 29 (5), 1849–1857. Menter, F.R., 1993. Zonal two equation k-w turbulence models for aerodynamic flows. 24th Fluid Dynamics Conference, AIAA Paper 93-2906, Orlando, Florida, 6-9 July 1993. Nakamura,, S., Naito,, S., 1977. Propulsive performance of a container ship in waves. Soc. Nav. Archit. Jpn. 15, 24–48. Orihara,, H., Miyata,, H., 2003. Evaluation of added resistance in regular incident waves by computational fluid dynamics motion simulation using an overlapping grid system. J. Mar. Sci. Technol. 8, 47–60. Park,, D.M., Lee,, J., Kim,, Y., 2015. Uncertainty analysis for added resistance experiment of KVLCC2 ship. Ocean Eng. 95, 143–156. Park,, S., Oh,, G., Rhee,, S.H., Koo,, B.-Y., Lee,, H., 2015. Full scale wake prediction of an energy saving device by using computational fluid dynamics. Ocean Eng. 101, 254–263. Park,, S., Park,, S.W., Rhee,, S.H., Lee,, S.B., Choi,, J.-E., Kang,, S.H., 2013. Investigation on the wall function implementation for the prediction of ship resistance. Int. J. Nav. Archit. Ocean Eng. 5, 33–46. Park,, S., Lee,, H., Rhee,, S.H., 2016. Numerical investigation of anti-diffusions term for free-surface wave flow. J. Adv. Res. Ocean Eng. 2 (2), 48–60. Park,, S., Yang,, J., Rhee,, S.H., 2017. Parametric study on ship's exhaust-gas behavior using computational fluid dynamics. Eng. Appl. Comput. Fluid Mech. 11, 159–171. Patankar,, S.V., Spalding,, D.B., 1972. A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows. Int. J. Heat. Mass Transf. 15, 1787–1806.

327