Acknowledgements k ω Abstract

(1)

IN

DEGREE PROJECT ENGINEERING PHYSICS, SECOND CYCLE, 30 CREDITS

,

STOCKHOLM SWEDEN 2020

Evaluation of drag

estimation methods for ship

hulls

(2)

(3)

Abstract

This study aims to evaluate to which degree CFD can complement traditional towing-tank experiments, and to develop a computationally cheap and robust methodology for these type of simulations. Two radically different surface ship hulls were chosen for the trials: a capesize bulk carrier and a fast displacement hull. A bare submarine hull was also used to benchmark the software in the early stages of the study. All simulations were Reynolds–Averaged Navier–Stokes (RANS) simulations using the

k

-

ω

-SST turbulence model. The chosen software was OpenFOAM 5.x and foam-extend 4.1 coupled with the commercial expansion Naval Hydro Pack, which is supposed to handle high Courant numbers well.

Producing a high-quality mesh which is able to capture both the free surface and the boundary layers proved to be of utmost importance and the meshing procedure is thus discussed in detail. A majority of the surface ship simulations suffered from a phenomenon known as numerical ventilation. The effect seemed to be much worse for the fast-displacement hull and the drag estimates for this hull deviated around 16.1% from experimental data. The bulk carrier was less affected and the drag estimates for this hull only deviated around 6.3% from experimental data. In order to reduce the computational cost, all surface ship simulations were run with a maximum Courant number of 25 and some consequences of this are also discussed.

Acknowledgements

(4)

(5)

1 Introduction

Traditional towing-tank tests are notoriously complex to perform. First an accurate scale model of the ship hull needs to be manufactured, where parameters such as displacement, weight distribution and surface geometry will be sources of error. The next step is to design a measurement campaign where the model is tested under the right conditions. Here the measurement devices, and possibly the testing facilities will introduce errors to the measurements. Finally the data needs to be scaled up to give estimates of the performance of the full scale ship and this process will both introduce new errors and possibly scale up the previously mentioned errors. For these reasons, towing tank tests can easily become very expensive and are usually only performed near the end of the design campaign.

A numerical towing tank based on Computational Fluid Dynamics (CFD) could instead allow evaluations of designs multiple times during the design process. See [1] for a relatively recent summary of state of the art numerical ship hydrodynamics. Simulations of full-scale ships would most likely be too computationally expensive due to the high Reynolds numbers and large meshes, but it would be possible to run simulations of models of different scales to further investigate how the final results can be scaled up as accurately as possible. When developing the simulation methodology emphasis was put on keeping down computational costs to allow entire series of hull forms to be tested within a reasonable time frame. In order to ensure that the methodology would work for many different types of ships three vastly different hull types were tested, one high speed displacement hull, one bulk carrier hull belonging to the largest class of bulk carriers and one bare submarine hull. The software chosen as basis for the numerical towing tank was the open source finite volume based software OpenFOAM.

1.1 Submarines

The idea of submarines has been around for hundreds of years but it wasn’t until the end of the 18th century that working prototypes were built, the first being Bushnell’s Turtle followed by Fulton’s Nautilus. These early submarines were man-powered military vessels designed for the purpose of planting explosive charges on the underside of warships. Although the vessels were technologically advanced for their time, they had little to no success.

(7)

Figure 1: Top: German Type VII-C U-boat (1930s) Bottom: German Type XII U-boat (1940s) Towards the end of the Second World War German designers started to realize the improvements in underwater performance that could be gained by streamlining the hull and removing appendages. The difference in design philosophy can be seen by comparing the earlier Type VII with the later Type XII [2] (see Figure 1).

After the war both the British and American navies obtained Type XII’s for evaluation and were stunned by the technological advances that the Germans had made. Their subsequent designs, the British Porpoise-class and American Tang-class, were heavily influenced by the German submarine [3] [4]. Modern day submarines are even more streamlined and are usually based around an axisymmetric base hull, such as the Joubert Bare Hull which is discussed in Chapter 2.

1.1.1 Resistance

The resistance acting on a body moving through a 1-phase fluid (constant density) is usually divided into two components, pressure resistance (or form resistance), R_p, and frictional resistance, R_f.

Pressure resistance: Since all fluids are viscous to some degree there will be a boundary layer on all surfaces of the body in contact with the fluid. This disrupts the fore-aft symmetry of the pressure acting on the body, resulting in a net drag force known as the pressure resistance.

Frictional resistance: The no-slip condition at the surface of the submarine means that the fluid counteracts the movement of the body by exerting tangential shear forces on the surface as it moves forward. The resulting drag coming from these shear forces is referred to as the frictional resistance

These resistances are often scaled with the dynamic pressure acting on the body to produce the non-dimensional coefficients

C

p and

C

f. The scaling is done according to

C

p

=

R

p 1 2

ρu

2∞

S

(1)

and similarly for

C

f, where

ρ

is the density of the fluid,

u

∞ is the inflow velocity and

S

is a characteristic area which in this study will be taken as the wet area of the submarine. We want to emphasize that

C

p, in this study, refers to the streamwise (

x

) component of the pressure coefficient. Furthermore, when the coefficients are written with capital letters such as

C

pthis refers to the integral of the local coefficient values over the entire body, according to

C

p

=

Z

Sh

(8)

1.2 Surface Ships

Unlike submarines, surface ships have been around for thousands of years and even the earliest designs showed some understanding of ship hydrodynamics. Up until the invention of propeller driven ships it was however not as important to estimate the hydrodynamic drag on the hull. In order to properly dimension the propulsion system of a ship it is crucial to have a good estimate of its resistance at design speed, so when the first steamships were developed in the mid 1800s this immediately became a problem and many of the early propeller driven ships thus had badly dimensioned propulsive systems.

In the 1870s William Froude started developing methods for estimating the resistance of full scale ships. He experimented with scale models and performed towed resistance tests using planks of different length and surface texture. A few of Froude’s ideas are discussed in more detail in Section 1.4

1.2.1 Resistance

For bodies moving close to or on the free surface of the fluid another type of resistance becomes important, namely the resistance associated with the production of waves. With this in mind it is possible to make another division of the total resistance, R_T, into viscous resistance, R_v, and residual resistance, R_R. Using equation (1) this can be written in terms of non-dimensional coefficients as

C

T

= C

v

+ C

R

.

(3)

For surface ships the area under water will be used in the scaling of the coefficients since the drag from the water is vastly dominating. This is discussed further in Chapter 4.

Viscous resistance: The component of the total resistance here referred to as viscous resistance is the sum of the frictional resistance, Rf, and the pressure resistance Rp. The reason for this grouping is that the pressure resistance can be argued to be the result of viscous effects, much like the frictional resistance. The pressure resistance is sometimes also referred to as viscous pressure resistance. This viscosity dependence means that the viscous resistance will scale with the Reynolds number of the flow, defined as

Re

=

ρuL

µ

=

uL

ν

=

Inertial forces Viscous forces

,

(4)

where

ρ

is the density of water, L is a characteristic length which in most naval applications is taken as the length of the ship,

µ

is the dynamic viscosity of water and

ν

is the kinematic viscosity of water.

Residual resistance: The residual resistance is what is left of the total resistance once the viscous resistance has been subtracted. This component is usually dominated by the so called wave-making resistance and is even denoted as such in some references dealing with large displacement ships such as bulk carriers or container ships. When considering high speed crafts, other factors such as wave-breaking and spray can also become important. All these components being related to surface waves, which are governed by gravity, means that the residual resistance will scale with the Froude number which is defined as

Fn

=

√

u

L

wl

g

=

Inertial forces

Gravitational forces

,

(5)

(9)

1.3 Ship terminology

For clarification a few of the terms specific to the ship building community are here explained since they are used both throughout this report as well as in many of the references.

Figure 2: Some ship measurements on a bulk carrier

Length between perpendiculars (L_pp) is the distance between the forward and aft perpendiculars. On

wooden ships the stem was usually taken as the forward perpendicular and the stern post as the rear perpendicular, however on many modern ships it is not as simple. For instance on the bulk carrier in Fig (2) the forward perpendicular is right above the bulb on the bow and the aft perpendicular is the planned axis of rotation of the rudder.

Waterline length (L_wl) is the length of the ship measured at the waterline, at design load. Beam at waterline (Bwl) is the widest width of the ship at the waterline, at design load. Draft (T) is the distance from the lowest point of the hull to the surface.

Wet area (S) is the area of the hull which is under water.

Displacement (

∇

) is the volume of water being displaced by the presence of the hull. It is directly proportional to the weight of the ship and the density of the water.

Block coefficient (C_B) is the ratio of the displacement

∇

to a box with the dimensions L_wl x B_wl x T. Amidships is the midpoint between the two perpendiculars.

Sinkage is the difference in draft experienced when the ship is moving relative to when it is stationary. Where along the hull the sinkage is measured can vary, but in this study it is defined as the vertical movement of the center of mass.

Trim is the difference in pitch experienced when the ship is moving relative to when it is stationary. It is sometimes also expressed as the difference in sinkage between the fore and aft perpendiculars.

Center of Gravity (CoG) is the center of mass of the ship.

(10)

1.4 Semi-Empirical Methods

Semi-empirical methods are a type of performance prediction methods used extensively by the ship building community. As the name implies they are based to some degree on empirical data, but also contain estimations and modeling. A common approach is to start with experimental data from scale models which can then both be extrapolated to full scale and interpolated to make estimates about similar hull forms.

One major problem when scaling up resistance data is the fact that C_vand C_Rscale with different non-dimensional numbers. Since it is practically impossible to keep both the Reynolds number and Froude number constant between a model test and a full scale ship, the total resistance cannot be scaled directly. The first systematic approach to this problem was suggested by Froude in the 1870s and consisted of running model trials at the correct Froude number, measuring the total resistance and then subtracting an estimate of the frictional resistance, leaving only the residual resistance.

C

R

= C

T

− C

f (6)

Note that Froude used the frictional resistance instead of the viscous resistance as in Section 1.2.1, simply because no theory explaining viscous boundary layers existed at that time.

This method was adopted by many ship builders, and at the International Towing Tank Conference (ITTC) in 1957 [6] a standardized formula for calculating the frictional resistance was accepted in favor of Froude’s experimental plank-towing data. The formula was originally derived from the skin friction of a flat plate but has been corrected to also take into account some three-dimensional effects. It can be written as

C

f

=

0.075 log(Re) − 2

2

.

(7)

This procedure was later developed further and now also includes a form factor (k) to take the pressure (form) resistance into account as well (see ITTC 1978 [7]). The expression for the viscous resistance then becomes

C

v

= (1 + k)C

f

,

(8)

and the expression for the total resistance becomes

C

T

= (1 + k)C

f

+ C

R (9)

(11)

1.5 Fluid Mechanics

The ground pillar of modern fluid mechanics is undoubtedly the Navier–Stokes equations. They are a set of coupled differential equations describing the conservation of mass, momentum and energy of a moving fluid. For the applications discussed in this report, both air and water can be considered to be incompressible fluids. Under these conditions the energy equation uncouples from the mass and momentum equations and it will therefore not be considered. Conservation of mass, often referred to as the continuity equation, can be written as

∂ρ

∂t

+

∂(ρu

i

)

∂x

j

= 0

(10)

and since water and air are Newtonian fluids the conservation of momentum can be written as

∂(ρu

i

)

∂t

+

∂(ρu

j

u

i

)

∂x

j

= ν

∂

2

_u

i

∂x

2 j

−

1 ρ

∂p

∂x

i

+ f

i

,

i = 1, 2, 3 ,

(11)

where

ρ

is the density of the fluid,

u

is the flow velocity,

p

is the pressure and

f

is the sum of all body forces. The subscripts

i

and

j

denotes directions in the cartesian coordinate system. Solving this set of equations yields the pressure and velocity fields of the flow domain which can then be used to calculate other flow quantities, such as the skin friction and pressure drag.

1.6 Pressure Gradient Flows

As soon as there is curvature in the geometry pressure gradients will form which will either decelerate or accelerate the flow, depending on the curvature. Pressure gradients which decelerates the flow are called adverse pressure gradients (APG,

∂p/∂x

> 0) while pressure gradients accelerating the flow are called favorable pressure gradients (FPG,

∂p/∂x

< 0).

Adverse pressure gradients will often increase the complexity of fluid mechanics problems and can in some cases lead to flow separation. This happens when the flow closest to the wall is decelerated to a standstill and starts flowing opposite to the outer flow (see Figure 3). Adverse pressure gradients can be expected to arise near the converging stern sections on both submarines and surface ships. The importance of pressure gradient effects on turbulent boundary layers is discussed in detail in two articles by Vinuesa et al. [8] [9].

Figure 3: Boundary layer with an APG

1.7 Turbulence

When the inertial forces of a flow grow large enough relative to the dampening viscous forces, the flow will turn turbulent. A turbulent flow is characterized by chaotic fluctuations in both flow velocity and pressure. These fluctuations give rise to eddies in a range of different sizes, from the scale of the geometry or domain (depending on the problem) down to the smallest viscous length scales. The presence of these eddies greatly increases the mixing rate of both momentum and energy, meaning larger wall shear stresses in wall bounded flows.

(12)

2 Cases

In this chapter the three test cases of the study are discussed and the chosen sources for reference data are summarized. The test cases were one submarine hull without appendages, one bulk carrier and one high speed displacement hull.

2.1 Case 1: Joubert Bare Hull

The DSTO Joubert Bare Hull (JBH) comes from a series of generic submarine geometries based on the work done by P.N. Joubert in the early 2000s [3] [4]. The series was developed to be representative of many modern submarines without having any full scale equivalent, meaning that the geometries can remain unclassified and be discussed in the open literature. A number of articles documenting both experimental testing and CFD of the bare hull configuration have been published and some of them are summarized below.

Figure 4: Joubert Bare Hull

Since this study mainly focuses on capturing the drag and the magnitude of its components using CFD, a number of articles presenting relevant experimental data have been chosen as references. The table below summarizes what kind of data can be found in these articles as well as some information about the experiments.

References Re Exp. / Sim Data obtained Cf Cp CT

Manovski et al. 2018 [10] 3.8x106 Exp. Integrated 3.0490

Clarke et al. 2016 [11] 12x106 Exp. (Raw) Integrated 2.8819 0.82246 3.7044

12x106 Exp. (Corr) Integrated 2.5954 0.63529 3.2307

Jones et al. 2013 [12] 3.5x106 Exp. Integrated 3.0941 1.6876 4.7817

5.4x106 Exp. Integrated 2.8757 1.6876 4.5632

Anderson et al. 2012 [13] 5.4x106 Exp. Measured 4.7581

5.4x106 RANS Measured 3.4821

Quick et al. 2012 [14] 5.4x106 Exp. (Trip) Measured 4.76

5.4x106 Exp. (No trip) Measured 4.03

Table 1: Summary of data available from the references for the JBH

(13)

Quick et al. 2012

The article by Quick et al. [14] documents the first phase of a larger study aiming to build a body of information about a generic submarine shape using both experimental testing and CFD. The report documents a series of wind-tunnel tests performed in the DSTO Low Speed Wind Tunnel with the purpose of gathering steady-state aerodynamic force data, as well as asses the flow characteristics around the body. In this early phase of the study only the bare hull of the submarine was considered. A model was machined from aluminum, measuring 1.35 m long with a fineness ratio (length over max diameter) of 7.3. This resulted in an estimated blockage ratio of 2.1% relative to the wind tunnel test-section. The surface roughness of the model was measured to be less than 0.8

µ

m. Force and moment coefficients are presented for a number of different angles of attack and sideslip. The authors mention that tests were done to check the influence that the support pylon and pitch-arm might have had on the flow angularities and the results showed that the influence was larger than expected. However, no attempts were made to correct the force and moment data since it will be part of an "incremental database" and subsequent tests of fully appended submarines will use the same mounting. The wake effect from the model support can clearly be seen in the asymmetry in the Cx (CT in

x

-direction) vs. angle of attack data as well, indicating that the presented total drag coefficient might very well be too high.

Anderson et al. 2012

Anderson et al. [13] documents a joint study between Swedish FOI and Australian DSTO. Both experimental and computational results are presented. The experiments were performed in the DSTO Low Speed Wind Tunnel located in Melbourne, Australia. Two different models were used, one of the bare hull and one of the fully appended configuration. Both models were 1.35 m long and the bare hull version was the same model as the one used by Quick et al. [14]. The bare hull configuration, which is the main interest of the current study, was run at a free stream velocity of 60 m/s, meaning a Reynolds number of 5.4 million. The skin-friction coefficient, cf, was measured along the model using the Preston tube method, described in detail in the original article by Preston [15]. The authors claim that the total drag coefficient, C_T was scaled using the surface area of the submarine, but when studying their data it seems more likely that the length,

L

, squared was used. Since the experimental setup was the same as that used by Quick et al., the previously mentioned effect of the support pylon can be expected to be visible in these results as well. The authors mention that the effect can be seen in the PIV measurements in that the turbulence levels are higher on the lower side of the model, and the downwards divergence of the wake. It is also mentioned that the pressure coefficient is slightly asymmetric in the lower half of the near wake, which according to the authors is most likely due to the support pylon being misaligned.

Jones et al. 2013

The article by Jones et al. [12] documents a study investigating different kinds of tripping devices to be used with the Joubert generic submarine geometries. The study focuses on the bare hull configuration and the effect that the tripping devices have on the skin friction and pressure coefficients along the model. All measurements were done in the DSTO Low Speed Wind Tunnel and the bare hull model used was the same 1.35 m model used by Quick et al. [14]. The skin friction coefficient was measured using the Preston tube method [15]. The tripping devices used were wires with diameters of 0.1 mm, 0.2 mm and 0.5 mm, and also a 3 mm wide tape on which grit of an average particle diameter of 190

µ

m had been distributed. Curves of the skin friction coefficient, cf, and pressure coefficient, cp, along the body are presented for a range of velocities for all the different tripping devices as well as for the case without tripping. For the 0.1 mm tripping wire, measurements were only taken at 3 measurement locations since it was found that this tripping wire merely moved the transition point a bit upstream rather than all the way to the wire.

(14)

Clarke et al. 2016

Clarke et al. [11] is an experimental study of the JBH in the cavitation tunnel at the Australian Maritime Collage. The model used was the same 1.35 m long model described by Quick et al. [14], using a different more streamlined mounting solution. The blockage was estimated to be around 8.5% and thus a blockage correction based on CFD-results was applied to the measured data. Skin friction measurements were done at Reynolds numbers of 6x106 and 12x106 using the Preston tube method [15] and this data was presented in terms of c_f vs

x

/L curves. The characteristics of the cavitation tunnel, including free stream turbulence intensity, are discussed in another article by D. Butler et al. [16].

Manovski et al. 2018

Manovski et al. [10] describes an experimental study where a long-distance Particle Image Velocimetry-system (PIV) was benchmarked for the purpose of measuring mean as well as fluctuating velocity components in turbulent boundary layers. The system was tested using a 2 m long carbon-fiber composite model of the JBH in the DSTO Low Speed Wind Tunnel. All tests were run at zero angle of incidence with a free stream velocity of 28.8 m/s, meaning a Reynolds number of roughly 3.8

×

106.

The PIV data was used to produce estimates of boundary layer profiles and skin friction along the model. As a reference the skin friction coefficient, c_f, was also measured along the model using the Preston tube method [15].

Conclusions

After reading through all the references, the articles by Anderson et al. [13] and Clarke et al. [11] were chosen as the primary references. Even though the wind-tunnel experiments by Anderson et al. seem to have been affected by the mounting pylons to some degree, c_f was measured on the top side of the model where the effect is expected to be small. Furthermore the experiments, including facilities, models and tripping devices, are very well documented which is of great value. The study by Clarke et al. was chosen since it uses the same model as Anderson et al. but a different testing facility and mounting solution. It is also good that the authors estimate the blockage and makes an attempt to correct it.

Integration of cf and cp

One problem with references using Preston tube measurements for c_f is that the method relies on the boundary layer being turbulent. For this reason, c_f data is only available from, at best, right after the tripping location. This means that the integration of c_f misses out a significant part of the total friction. The contribution from these neglected parts was estimated for the data from Anderson et al. and Clarke et al. This was done by cutting the c_f curves from the corresponding simulations and then comparing them with the uncut curves.

In the experiments by Anderson et al. c_f was sampled from x/

L

= 0.05 to x/

L

= 0.79 (measured from

the nose of the model). In the experiments by Clarke et al., c_f was sampled from x/

L

= 0.13 to x/

L

= 0.96 (measured from the nose of the model). The results from the comparisons are presented in the tables below.

C_f (Full) 3.1265

Cf (Cut) 2.6605

C_f (Full) 2.7726

Cf (Cut) 2.3999

Table 2: Left: Anderson et al. Right: Clarke et al.

(15)

2.2 Case 2: Japanese Bulk Carrier

The Japanese Bulk Carrier (JBC) is a generic capesize bulk carrier designed through a joint program involving the National Maritime Research Institute (NMRI), Yokohama National University (YNU) and the Ship Building Research Centre of Japan (SRC), all from Japan. Capesize meaning that it is too large to transit any of the canals, and thus has to travel around the capes. Figure 5 shows the JBC model from the side and Table 3 presents a few properties of the hull.

Figure 5: Japanese Bulk Carrier

symbol Model Full Scale

Length (between perpendiculars) LP P 7.0 m 280 m

Length (Waterline) L_wl 7.125 m 285 m

Beam at water line B_wl 1.125 m 45 m

Draft T 0.4125 m 16.5 m

Wet area S 12.22 m2 19556.1 m2

Displacement volume

∇

2.787 m3 178369.9 m3

Block coefficient C_B 0.8580 0.8580

Form factor k 0.314 0.314

Long. pos. of CoB fwd amidship L_CoB 2.5475 % of L_{P P} 2.5475 % of L_{P P}

Slenderness L_wl /

∇

1/3 5.0630 5.0630

(16)

Hino et al. 2016

Hino et al. [17] was a study with the purpose of producing benchmark data for a ship equipped with an energy saving device (ESD), to be used as validation for future CFD-analyses. The chosen ship was the JBC of which two scale models were constructed, one 7 m long at NMRI and one 3.2 m long at Osaka University.

Both towed tests and self propulsion tests were performed. Stereo-PIV was used in both cases to capture the flow field between the ESD and the propeller, and resistance data was collected for the 7 m model. The towed resistance measurements were performed in the large towing tank at NMRI, measuring 400 m long, 18 meters wide and 8 m deep. Data was collected for Froude numbers ranging from 0.12 to 0.16 meaning Reynolds numbers ranging from 6.28x106 to 8.38x106. Both the Froude numbers and Reynolds numbers were calculated using the length between perpendiculars, L_pp.

The data was presented in the form of the residual resistance and through private communications with one of the authors it was concluded that they had used equation (8) with the formfactor from Table 3 to subtract an estimate of the viscous resistance. The data is plotted in Figure 6 together with the viscous resistance estimate, as well as the sum of the two which should be the measured total resistance. This data will be used as reference data for all the simulations of the JBC.

0.12 0.13 0.14 0.15 0.16 Fr 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 C x10 3 C R C f C v 0.12 0.13 0.14 0.15 0.16 Fr 4.3 4.31 4.32 4.33 4.34 4.35 4.36 4.37 4.38 CT x10 3

Figure 6: Left CR(Cw), Cf and Cv vs Fr Right CT vs Fr

Also available in the article was the trim and sinkage at the design Froude number. These are presented in Table 4. Notice that the trim is here expressed as the difference in sinkage between the perpendiculars. The trim and sinkage will also be used to check the validity of the results, since capturing the ships position in the water is essential to get a good estimate of the resistance.

Sinkage (% of L_pp) -0.086

Trim (% of Lpp) -0.180

(17)

2.3 Case 3: MARIN Systematic Series FDS-5

The MARIN Systematic Series is a series of high-speed displacement ship hulls which were developed in the 1980s through a joint industry program involving Delft University and the Maritime Research Institute Netherlands (MARIN) as well as the navies of the Netherlands, Australia and the US. The purpose of the program was to evaluate the seakeeping and endurance of high speed displacing mono-hull vessels in waves. The basis of the series is the parent hull called the “fast displacement ship 5” (FDS-5) which was chosen through a pre-study. The other hulls in the series are systematic iterations of the FDS-5 where the parameters being varied are length over beam ratio, beam over draft ratio and block coefficient. The entire program, including all results, was summarized in a book by Kaspenberg et al. [18] published by MARIN in 2014. Figure 7 shows the FDS-5 model used and Table 5 presents a few properties of the hull.

Figure 7: MARIN Systematic Series FDS-5

symbol Model

Length (between perpendiculars) LP P 5.0 m

Length (Waterline) L_wl 5.0 m

Beam at water line Bwl 0.625 m

Draft T 0.15625 m

Wet area S 2.9288 m2

Displacement volume

∇

0.195 m3

Block coefficient C_B 0.396

Long. pos. of cB fwd amidship LcB -5.11% of LP P

Slenderness L_wl /

∇

1/3 8.62

(18)

Kaspenberg et al. 2014

As mentioned, the book by Kaspenberg et al. [18] summarizes the entire joint program including, of interest to this study, a measurement campaign where models of all the hull forms were manufactured and tested. The measurements were performed in the towing tank at MARIN in the Netherlands and consisted of both seakeeping tests in waves and calm water resistance tests. In the book there is also information about the design of the parent hull as well as instruction on how the data can be interpolated to other similar hull forms.

Tabulated data from the resistance tests is presented for all hull forms and a range of velocities. The data includes the total drag as well as the sinkage at the fore and aft perpendiculars. The data for the FDS-5 is plotted in Figure 8, where C_f was calculated using the ITTC-57 friction line and C_R is simply C_T minus C_f. The sinkage at the fore and aft perpendiculars was recalculated to the trim and sinkage at the CoG and this is also presented below.

0.2 0.4 0.6 0.8 1 1.2 Fn 1 1.5 2 2.5 3 3.5 C x 10 3 C_f C_R 0.2 0.4 0.6 0.8 1 1.2 Fn 3.5 4 4.5 5 5.5 C T x 10 3 0.2 0.4 0.6 0.8 1 1.2 Fn -20 -15 -10 -5 0 5 Sinkage [mm] 0.2 0.4 0.6 0.8 1 1.2 Fn -2 -1.5 -1 -0.5 0 Trim [ ° ]

Figure 8: Left C_Rand C_f vs Fn Right C_T vs Fn

(19)

3 Method - Computational Fluid Mechanics

Developing a theory for solving the Navier–Stokes equations with turbulence analytically has proven to be notoriously difficult, so difficult in fact that a large part of modern day fluid mechanics research is instead conducted within the field of Computational Fluid Dynamics (CFD). Here the aim is to discretize the governing equations and solve them numerically using the ever increasing power of computers.

The process of discretizing the governing equations and then solving them numerically without any further simplifications is called Direct Numerical Simulation (DNS). Such simulations resolve the entire spectra of turbulent fluctuations and thus require very fine computational grids and are often extremely computationally expensive. The number of floating-point operations required also increases as Re3, meaning that this type of simulation only is feasible for low to moderate Reynolds numbers. Direct numerical simulations are almost exclusively used for research and the industry instead relies on different kinds of modeling to reduce the complexity of the problem.

The approach used in this study is known as Reynolds–averaged Navier–Stokes (RANS) and is discussed in more detail in the section below. Other crucial elements to surface ship simulations are the ability to model two phases (air and water) at the same time, as well as the ability to model the movement of the ship hull with respect to the free surface. The software used throughout this study was the open source software OpenFOAM, both version 5.X and the fork foam-extend coupled with the commercial extension Naval Hydro Pack. This is discussed in more detail in Chapter 4.

3.1 RANS

As discussed resolving the entire spectra of turbulent scales is often extremely computationally expensive and in many practical applications the bulk flow is of much greater significance than the turbulent fluctuations so studying the averaged flow quantities can often be sufficient. This is done by first decomposing the flow into an averaged part and a fluctuating part in a process known as Reynolds decomposition. Using this decomposition the flow quantities can be written as

u

i

= U

i

+ u

0i

p = P + p

0 (12)

where

U

i and

P

are the averages and

u

0i and

p

0are the fluctuations with zero mean. Substituting (12) into the momentum equations (11) and then averaging results in the Reynolds–averaged Navier–Stokes equations,

∂U

i

∂t

+ U

j

∂U

i

∂x

j

=

∂

∂x

j

h

ν

∂U

i

∂x

j

−

1 ρ

P δ

ij

− hu

0 i

u

0 j

i

,

i = 1, 2, 3 ,

(13)

where

h

...

i

denotes average. These equations are almost identical to the regular Navier–Stokes momentum equations (11) with the exception of the new term

hu

0_i

u

0_j

i

, which is a tensor containing the covariances of the velocity fluctuations. This term is usually referred to as the Reynolds stress tensor even though it is technically

ρhu

0_i

u

0_j

i

which has the dimension of a stress. The introduction of this stress tensor means that there are now 6 additional unknowns resulting in a total of 10 unknowns and 4 equations. This leads to the so called closure problem of RANS. The two most common methods for dealing with this problem are to use either eddy viscosity models (EVMs) or Reynolds-stress models (RSMs).

Eddy viscosity models model the Reynolds stresses using the Boussinesq eddy–viscosity hypothesis which states that

hu

0_i

u

0_j

i =

2

3 kδ

ij

− 2ν

T

hS

ij

i ,

(14)

(20)

In Reynolds-stress models individual transport equations are instead solved for each Reynolds-stress component, thus removing the need for any additional modeling. This does however also make RSMs relatively computationally expensive for being RANS models. This study focuses on the use of eddy viscosity models which are explained in more detail below.

3.2 Eddy viscosity models

The purpose of eddy viscosity models is simply to determine an expression for the eddy viscosity and this expression usually consists of a velocity scale

u

∗

(x

i

, t)

and a length scale

l

∗

(x

i

, t)

. The simplest EVMs, algebraic models, relates these scales to the geometry of the flow and thus require the user to know details of the flow beforehand. The much more commonly used two-equation models instead relate

u

∗

(x

i

, t)

and

l

∗

(x

i

, t)

to turbulent quantities of the flow and then solve separate transport equations for these quantities. The workings of a few such models is discussed below [19] [20].

3.2.1

k

-

model

The

k

-

model was first developed by B.E Launder and W.P Jones in the 1970s [21] and has since evolved and improved significantly. The model was one of the first examples of a two-equation model, a model where the eddy viscosity is expressed in terms of two turbulence quantities which are solved for using two model transport equations [20]. The

k

-

model uses the turbulent kinetic energy,

k

, and the turbulence dissipation rate,

. The

k

-transport equation can be written as

∂k

∂t

+ U

j

∂k

∂x

j

=

∂

∂x

i

"

ν

T

σ

k

∂k

∂x

i

#

+ hu

i

u

j

i

∂u

i

∂x

j

−

(15)

and the

-transport equation for incompressible flow can be written as

∂

∂t

+ U

j

∂

∂x

j

=

∂

∂x

i

"

ν

T

σ

∂

∂x

i

#

+ C

1

k

hu

i

u

j

i

∂u

i

∂x

j

− C

2

k

,

(16)

where

σ

k,

σ

,

C

1 and

C

2 are model constants. The eddy viscosity is then estimated as

ν

T

=

c

ν

ρk

2

,

(17)

where

c

ν is a constant which has been empirically determined to be approximately 0.09. The model is widely used today, but is known to struggle with adverse pressure gradient flows for which it predicts excessively high shear stresses. [22]

3.2.2

k

-

ω

model

The

k

-

ω

model was presented by David C. Wilcox in 1988 [23] and uses the turbulent kinetic energy,

k

, and the specific rate of dissipation of

k

,

ω

, as its two model parameters. Wilcox wrote that

ω

can be thought of as“... the ratio of the turbulent dissipation rate

to the turbulent mixing energy

k

”. It is defined as

ω =

β · k

,

(18)

where

β

is a closure coefficient equal to 3/40. The transport equation for

k

is exactly the same as in the

k

-

model and the transport equation for

ω

can be written as

(21)

where

σ

ω,

C

ω1and

C

ω2are model constants. Unlike the the

k

-

model, the

k

-

ω

model has proven to be very accurate when modeling boundary layers with adverse pressure gradients. It does however come with one major setback: It requires a non-zero inflow condition for

ω

and is very sensitive to this value. According to F. R. Menter [22] the effect of the inlet

ω

value on the magnitude of the modeled eddy viscosity can be as large as 100%. The model is also know to under-predict the spreading rate for free shear layers.

3.2.3

k

-

ω

-SST model

The

k

-

ω

shear stress transport model,

k

-

ω

–SST in short, was published by F. R. Menter in 1994 [22] and can be thought of as a merger between the the

k

-

model and the

k

-

ω

model. The model makes use of a blending function which allows it to use the

k

-

ω

formulation in the vicinity of a wall and then gradually change to the

k

-

formulation as the distance to the wall grows larger. The blending function takes on values between 0 and 1 and uses the wall normal distance as a switching variable.

The

k

-

ω

-SST model is able to handle strong adverse pressure gradients without the heavy reliance on free stream conditions. For this reason it was chosen as the primary turbulence model to be used in this study.

3.2.4 Turbulence at the free surface

Traditional two-equation turbulence models, like the ones used in this study, are known to struggle near the free surface. The abrupt change in fluid density at the surface often makes the models over-predict the production of turbulence. In reality the surface restricts the turbulent length scales, and stretches the eddies parallel to the surface which leads to increased dissipation in the direction normal to the surface. [24]

The effect of the free surface on the friction force is however assumed to be negligible and since the focus here is on drag, the inability of traditional two-equation models to capture turbulence near the free surface is considered to be outside the scope of this study.

3.3 Solving algorithm

Due to the coupled nature of the Navier-Stokes equations the pressure field needs to be known in order to calculate the velocity field and the velocity field needs to be known to calculate the pressure field. One way around this problem is to use a predictor-corrector algorithm, i.e. guess one of the fields, calculate the other field, go back and correct the initial guess and then repeat until convergence. Two very common algorithms of this type are the PISO and SIMPLE algorithms, which are discussed briefly below.

3.3.1 SIMPLE

The Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) is an iterative algorithm that often requires multiple iterations per time step before convergence. An overview of the algorithm-steps is presented below. [25] 1) Guess the pressure field,

p

, or if available take the pressure field from the previous time step.

2) Calculate the velocity field,

u

∗, from the pressure field using a linearized version of the momentum equations. 3) Calculate a pressure correction,

p

corr.

4) Use the correction to obtain a corrected pressure field,

p

∗, and a corrected velocity field,

u

∗∗, which satisfies continuity.

5) Solve any additional transport equations.

(22)

3.3.2 PISO

The PISO (Pressure-Implicit with Splitting of Operators) algorithm is based on the SIMPLE algorithm but splits the solution process into an implicit predictor step and one to several explicit corrector steps where operations on the pressure field are uncoupled from the velocity. This means that the momentum equation is only solved once per time step which in turn means that any errors in the initial pressure field remain unaffected by the following iterations. [26] An overview of the steps of the PISO algorithm is presented below.

Predictor step) This step is exactly the same as step 1 and 2 in the SIMPLE algorithm although the predictor step is only performed once per time step. A first estimate of the velocity field, u∗, is calculated using the linearized momentum equation and the pressure field, p, which is either taken from the previous time step or guessed if not available.

Corrector step) The corrector step is the same as steps 3 and 4 in the SIMPLE algorithm and this step can be repeated multiple times per time step. The velocity field calculated from the momentum equation will not be divergence free and will thus not satisfy the continuity equation. A pressure correction is therefore used to obtain a corrected pressure, p∗. This pressure is then used to calculate the corrected velocity, u∗∗. As stated earlier this corrector step can be repeated multiple times by calculating a new pressure correction from the corrected fields p∗ and u∗∗. Although convergence is often reached within just a few iterations and the algorithm is even considered by some to be non-iterative. [25]

3.3.3 PIMPLE

The PIMPLE algorithm is a hybrid between SIMPLE and PISO. It basically solves a predetermined number of SIMPLE loops each time step, with a predetermined number of PISO pressure corrections within each loop. The number of SIMPLE loops to perform is ususally refered to as the number of outer correctors (nOuterCorrectors) and the number of PISO loops to perform is called simply the number of correctors (nCorrectors). The workings of the PIMPLE algorithm as well as its implementation in OpenFOAM is explained in detail in the book by T. Holzmann [27].

3.3.4 Convergence

OpenFoam allows 3 solver convergence criteria to be specified and if one of them is reached the solver will stop. The criterias are maxIter which is the maximum allowed number of iterations allowed, tolerance which is the residual value considered to be sufficiently low and relTol which is the value of the current residual divided by the initial residual at which to stop.

3.3.5 Hydrostatic Pressure Effects

For the two-phase simulations an alternative pressure, pd, is solved for instead of the total pressure. This alternative pressure is defined as

pd

= p − ρgh .

(20)

According to the authors of the OpenFOAM Users Guide [28] this formulation is numerically convenient for cases where the hydrostatic pressure is important.

3.3.6 Unsteady RANS

Even though RANS solutions are averaged it is possible to capture the time dependence of the mean flow if small enough time steps are used. Such simulations are referred to as Unsteady RANS (URANS). In OpenFOAM the PISO and PIMPLE solvers (PisoFoam and PimpleFoam) are transient solvers which can be used for URANS. They are thus somewhat sensitive to the size of the time step.

(23)

3.4 Wall functions

In order to properly predict the friction forces on any surface it is of great importance to capture the viscous and turbulent mechanisms of the boundary layer. The inner regions of the boundary layer adds significant complexion to the process of turbulence modeling. Here the wall both dampens wall normal components of the turbulence, rendering it anisotropic, and increases the production of turbulence through shear forces. Moreover, the profiles of

U

and

are very steep in this region meaning that a very fine mesh is required to capture the physics. The result is that a substantial fraction of the computational cost of wall-resolved RANS simulations comes from this near wall region, and many two-equation models even fail to properly model the near wall physics. An often much computationally cheaper approach is to use wall functions. These functions add boundary conditions a small distance out from the wall as well, thus reducing the need to solve the RANS equations in this region. The boundary conditions are generated based on analytical models describing the near wall flow.

When measuring distances close to the wall one often uses a non-dimensional wall–normal distance (

y

+) scaled with the viscous length scale. It is defined as

y

+

≡

u

τ

y

ν

(21)

where

u

τ is the friction velocity and

y

is the wall normal distance. The friction velocity is a velocity scale based on the viscous stresses at the wall and is defined as

u

τ

≡

pτ

w

/ρ

where

τ

w is the wall shear stress. When using wall functions it is desirable to have the first grid point in the lower log-law region, i.e. a

∆y

+ w somewhere between 30 and 300. Wall functions have proven to work well for parallel flows with weak or ideally zero pressure gradients making them a viable choice for many ship hull geometries.

3.5 Two-phase flows

Simulations of surface ships adds one extra degree of difficulty compared to simulations of vessels moving through a single medium (such as submarines or aircraft). In order to properly resolve the domain around a surface ship both of the two phases (air and water), as well as the free surface between them needs to be modeled accurately. The most common way to do this is to consider both phases as part of one continuous medium and model the free surface as a jump in density and viscosity, meaning that only one set of governing equations is required to describe the motion of the entire domain.

Since marine hydrodynamics problems can be considered incompressible, the density of both air and water are assumed to be constant. The difference in density across the free surface can be expressed as

[ρ] = ρ

−

− ρ

+

_,

₍₂₂₎

where the superscripts + and - denotes an infinitesimal distance from the surface on the water and air side respectively. The velocity field is is also defined as being continuous according to

[u] = u

−

− u

+

= 0 ,

(23)

which basically says that the velocity of the air in the immediate vicinity of the surface must be the same as the velocity of the water just under the surface. Finally the surface tension is neglected which results in a continuous pressure field,

[p] = 0 .

(24)

Equations 22, 23 and 24 together make up the jump conditions which needs to be satisfied across the free surface. How the viscosity is handled is described in the next section.

(24)

3.5.1 Volume of fluids

The Volume of Fluids method was first described by Hirt and Nichols in 1981 [29] and is based on an indicator function which represents the volume fraction of water (V_w) in each cell to the total volume (V) of that cell. This indicator function can be written as

α =

V

w

V

,

(25)

where

α = 1

means water and

α = 0

means air. The effective kinematic viscosity is modeled as being

continuous across the surface. Using the indicator function

α

the effective viscosity can be written as

ν

e

= αν

e,w

+ (1 − α)ν

e,a

,

(26)

where

ν

e,w is the kinematic viscosity of water and

ν

e,a is the kinematic viscosity of air. This allows the use of conventional eddy-viscosity models to model the turbulence.

The indicator function is treated as a passive scalar which is convected by the velocity field. Its motion is described by the transport equation

∂α

∂t

+

∂(u

i

α)

∂x

i

+

∂

∂x

i

(u

r,i

α(1 − α)) = 0 ,

(27)

where the last term is a compression term to prevent excessive smearing of the surface and u_r is the relative velocity field normal to the surface. Equation 27, together with equations 11 and 10 makes up the set of governing equations for the two-phase flow cases. The compression term is described in more detail in the PhD thesis by Rusche [30], and its implementation in foam-extend is explained in the article by Vukčević et al. [31].

3.5.2 Ghost Fluid Method

One feature included in the Naval Hydro Pack extension is the ghost fluid method for polyhedral meshes. The Ghost Fluid Method is a method which implicitly enforces the jump conditions over the surface. It creates ghost cells at all points in the vicinity of the free surface, meaning that all these points contains the mass, momentum and energy of both the fluid in that cell and the ghost fluid which is the fluid on the other side of the interface (i.e. the ghost fluid for "air cells" is water and vice versa). The value of the VoF-indicator function (

α

) is then used to determine which of the fluids to consider in each cell. The Ghost Fluid Method is described in more detail in the articles [32] and [33] by Fedkiw et al. The implementation of the Ghost Fluid Method in foam-extend is described in the article by Vukčević et al. [34]. The method was used in all simulations of surface ships in this study.

3.6 Rigid Body Motion

Capturing the movement of the ship with respect to the surface is crucial when trying to estimate the drag. As the ship rotates and translates the wetted surface changes which affects both the skin friction and the procution of waves. For a ship under normal operating conditions all 6 degrees of freedom (DoF) needs to be considered, but during calm water resistance test the problem can be simplified to two DoFs. These two are translation along the z-axis (sinkage) and rotation about the y-axis (trim).

In nearly all marine hydrodynamics problems the motion of the vessel is almost entirely governed by pressure forces, greatly overshadowing the viscous forces. This means that the process of solving the pressure equations

and the process of solving the rigid body equations needs to be tightly coupled as well. The 6DoF code

(25)

4 Computational Setup

In this chapter the chosen software is briefly discussed followed by a detailed discussion of the computational setup of the three cases, including meshing, boundary- and initial conditions, solvers and schemes.

4.1 OpenFOAM

OpenFOAM is an open source software for solving continuum mechanics problems. Its CFD solvers are based on the Finite Volume Method meaning that the problem domain is divided into a finite number of control volumes over which the governing equations are then discretized and solved.

Naval Hydro Pack is a commercial expansion to OpenFOAM built on the fork foam-extend. The expansion offers an under-relaxed PIMPLE algorithm which is supposed to be very robust even at high Courant numbers. As previously discussed, it also offers the implementation of the ghost fluid method to better capture the free surface jump conditions. For a description of Naval Hydro Pack refer to the homepage of the company developing the expansion, Wikki Ltd.

4.2 Meshing

All meshes were generated using functions included in either OpenFOAM 5.x or foam-extend 4.1.

The first step in the construction of all meshes was to make a base-mesh using blockMesh. These base-meshes were then refined and snapped to the geometry using snappyHexMesh. For the two surface ships refineMesh was also used immediately after blockMesh to make sure that the cells in the vicinity of the hull were as "cubic" as possible as this tends improve the performance of snappyHexMesh. This was needed since the base-meshes for the final surface ship meshes had vertical grid refinement near the surface. This is discussed in more detail in Section 4.4 and 4.5

As mentioned, using wall functions to model the inner part of the viscous boundary layer puts some requirements on the mesh. The first grid point should ideally be between a

y

+ _{of 30 and 300, and this is best achieved by} generating a series of prism cell layers on the surface. All meshes in this study had 6 layers of rectangular prism cells closest to the surface. An expansion ratio of 1.2 was used and the thickness of the final cell layer was half of that of the first cell outside the layer. This is usually done to ensure a good transition between the prism layer and the outer mesh.

This meshing methodology was chosen since it is very flexible when changing between similar geometries and still produces quality meshes relatively fast. Another positive with snappyHexMesh is that it can be run in parallel with relative ease, thus reducing the time required to generate a mesh.

(26)

4.3 Case 1: Joubert Bare Hull

A computational domain was constructed for the JBH spanning 8

L

in the x-direction and 4

L

in the y- and z-directions, where

L

is the length of the model. This resulted in a blockage ratio of roughly 0.1%. The blockage ratio is calculated as the projected frontal area of the model divided by the cross sectional area of the domain.

4.3.1 Meshing

Due to the axisymmetric shape of the bare hull and the 0 degree inflow angle of the case, it would be possible to run half, or even a quarter of the body and then mirror the results. This was however not done since the methods developed might be used for dynamic maneuverability simulations in the future.

Initially two meshes were developed, one coarse and one finer, as a first step in a grid dependence study. A few mesh parameters from the two meshes are presented in the tables below.

Cells (Total) 830.964

Max Aspect Ratio 7.17

Non orthogonality (Max) 64.7

Non orthogonality (Avg) 6.60

Max skewness 1.04

(a) Mesh 1

Cells (Total) 1.511.314

Max skewness 0.79

(b) Mesh 2

Table 6: Mesh properties for meshes 1 and 2

Close ups of the two meshes can be found in Figures 9 and 10 and additional screenshots of both meshes can be found in Appendix A.

(27)

Figure 10: Close up of Mesh 2

4.3.2 Simulation

The steady-state solver simpleFoam, based on the SIMPLE-algorithm, was used to run the JBH cases. All simulations were run for 1000 iterations which proved to be more than enough for convergence (see Appendix B for residuals). As discussed in Chapter 3 the turbulence model used was

k

-

ω

-SST. Both meshes were run at Re = 5.4

×

106 and mesh 2 was also run at Re = 12

×

106. The boundary conditions are summarized in Table 7 and the initial conditions are presented in Table 8 (a) and (b).

k

ν

t

ω

p u

Inlet fixedValue calculated fixedValue zeroGradient fixedValue

Outlet inletOutlet calculated inletOutlet fixedValue inletOutlet

Sides Slip calculated slip slip slip

JBH kqWallFunction nutkWallFunction omegaWallFunction zeroGradient noSlip

Table 7: Boundary conditions for the JBH

(28)

The initial value of

k

was calculated as

k

0

=

3

2 (u

∞

I)

2 ₍₂₈₎

and the initial value for

ω

was calculated as

ω

0

= C

−1 4 µ

√

k

L

(29)

where

C

µ is a constant equal to 0.09 and

L

is the turbulent length scale of the problem, in this case taken as the length of the model (1.35 m). The turbulent intensities were taken from the two reference papers and these were 0.67% in the wind tunnel used by Jones et al. [12] and 0.5% in the cavitation tunnel used by Clarke et al. [16]. The resulting initial condition are presented in Table 8.

u

∞ 60 [m/s]

p

0 0

k

0 0.24 [m2/s2]

ω

0 0.66 [1/s]

Re 5.4

×

106

(a) Flow conditions calculated from Jones et al. [12]

u

∞ 134 [m/s]

p

0 0

k

0 0.67 [m2/s2]

ω

0 1.11 [1/s]

Re 12

×

106

(b) Flow conditions calculated from Clarke et al. [11]

Table 8: Flow conditions from the two chosen JBH references

A linear interpolation scheme with gauss discretization was used for all gradient calculations. The divergences were also discretized using a Gauss scheme and the interpolation schemes used are summarized in Table 9.

u

k

ω

Scheme linearUpwind Upwind Upwind

Table 9: The divergence schemes used in all JBH simulations.

(29)

4.4 Case 2: Japanese Bulk Carrier

The computational domain for the JBC spanned 6.5

L

in the x-direction, 3

L

in the y-direction and 4

L

in the z-direction resulting in a blockage ratio of 0.11% .

4.4.1 Meshing

Much like for the JBH, two meshes were also made for the Japanese bulk carrier as a first step of a grid dependence study. Some properties of the meshes are presented in the tables below.

Max skewness 2.09

∆

z at the free surface [mm] 22

(a) Mesh 1

Max skewness 8.04

∆

z at the free surface [mm] 11

(b) Mesh 2

These first two meshes were built in a way much like the meshes for the JBH, with the exception of a set of refinement boxes along the free surface to allow it to be properly resolved. Emphasis was put on having good resolution on the sides of the hull, expected to be in contact with water, and a coarser mesh was allowed along the upper surface of the model since the contribution from the air skin friction was considered to be negligible. For a bulk carrier of similar proportions, with all superstructures, the total air drag is estimated to be around 2.5% [5]. In these simulations towing tank models are used, which have no appended structures, meaning that the aerodynamic drag is considerably lower than 2.5%. Figures 11 and 12 show close ups of the two initial meshes.

(30)

Figure 12: Close up of Mesh 2

A set of simulations were performed to evaluate the quality of these meshes, including one stationary test where the ship was "dropped" with the waterline at the specified design waterline, and a dynamic test where the ship was run at the design Froude number. The results are presented in Figure 13 where "Diff" indicates the difference between the values from the dynamic and the static tests. Also included in Figure 13 is reference data from Hino et al..

0 20 40 60 80 100 Time (s) -0.02 -0.018 -0.016 -0.014 -0.012 -0.01 -0.008 -0.006 -0.004 -0.002 0 Sinkage [m] Sinkage 0 20 40 60 80 100 Time (s) -0.02 0 0.02 0.04 0.06 0.08 0.1 0.12 Trim [ °] Trim Fn = 0 - Mesh1 Fn = 0 - Mesh2 Fn = 0.142 - Mesh1 Fn = 0.142 - Mesh2 Diff - Mesh1 Diff - Mesh2 Fn = 0.142 - Exp.

Figure 13: Comparison of trim and sinkage for Mesh 1 and Mesh 2

(31)

Max skewness 2.49

∆

z at the free surface [mm] 9.9

(a) Mesh 3

Max skewness 2.52

∆

(b) Mesh 4

The new meshes were refined using directional refinement in z around the free surface already on the blockMesh stage. A block 1 m thick in z, centered on the free surface, spanning the entire domain in x and y was constructed and then split in the z-direction, 20 times for mesh 3 and 40 times for mesh 4. In order to create a smooth transition between this free surface mesh and the outer coarser mesh, grading was added on both sides. In order for snappyHexMesh to work properly it is important to keep the cells which are to be further refined as cubical as possible. For this reason the refineMesh utility was used to refine the meshes in the xy-plane within a series of trapezoidal prism around the model until the cells in the vicinity of the hull were of roughly the same refinement level in all directions. The trapezoidal refinement regions had an angle of 19◦ to properly capture the wave patterns. The regions can be seen in Figure 15. Finally snappyHexMesh was run with refinement and layer generation enabled on the hull. The main difference between mesh 3 and 4 is the level of directional refinement in the surface region added in the initial blockMesh stage. Figures 14 and 15 show Mesh 3 from two different angles. More figures of all meshes can be found in Appendix A. Discussion around the trim and sinkage predictions of mesh 3 and 4 can be found in Chapter 5.

(32)

Figure 15: Mesh 3 from above, cut along the free surface at rest.

4.4.2 Simulation

The transient solver navalFoam, based on the PIMPLE-algorithm, was used to run all JBC cases. All simulations used an initial

∆

t of 0.01 s and variable time stepping and meshes 1, 2 and 3 were all run for 100 s. Mesh 4 showed convergence problems in the drag force for some of the velocities so the simulations with this mesh were run for at least 300 s. This is discussed more in Chapter 5. As discussed in Chapter 3 the turbulence model used was

k

-

ω

-SST. Choosing free stream values of

k

and

ω

is not as straightforward for the surface ships as for the JBH since all experiments were performed in calm water conditions, i.e. no free stream turbulence. For this reason a few of the Naval Hydro Pack validation cases using similar geometries were studied and based on this the value of the free stream turbulence intensity was chosen to 5%. The initial conditions for

k

and

ω

were then calculated using equations (28) and (29).

Focus was put on getting good results from the case where Fn = 0.142 since this is the specified design Froude number and since this point had the most available reference data. A few more Froude numbers were also tested and the entire range is presented in the table below.

Fn 0.12 0.13 0.142 0.15 0.16

u 0.9944 1.077 1.179 1.2430 1.3259

Re 6.40x106 6.93x106 7.59x106 8.00x106 8.53x106

(33)

The boundary conditions used in all JBC simulations are summarized in Table 13.

α k ω pd u

Inlet waveAlpha fixedValue fixedValue zeroGradient waveVelocity

Outlet waveAlpha inletOutlet inletOutlet zeroGradient waveVelocity

Sides zeroGradient inletOutlet inletOutlet zeroGradient zeroGradient

Top inletOutlet inletOutlet inletOutlet totalPressure pressureInletOutletVelocity

Seabed inletOutlet kqRWallFunction omegaWallFunction zeroGradient fixedValue

JBC zeroGradient kqRWallFunction omegaWallFunction zeroGradient movingWallVelocity

Table 13: Boundary conditions for the JBC

The boundary condition for

u

on the hull had to be set to movingWallVelocity since an ordinary no-slip

boundary condition wouldn’t work with a moving geometry. waveVelocity is a boundary condition included in the Naval Hydro Pack which simply means that the velocity is set using the initWaveField command. Since no incoming waves are used it is essentially the same as using fixedValue.

As mentioned earlier a variable time-step was used, and this time step was adjusted according to a Courant number limiter. Since a part of this study consisted of evaluating the performance of Naval Hydro Pack, a quite high maximal Courant number of 25 was chosen. This was the value recommended by Wikki in their validation cases. Allowing a higher Courant number practically means allowing larger time steps, which can reduce the computational cost of the simulation. This is in line with the other goal of the study which was to develop a robust yet cheap methodology for performing this type of simulation. A test run with a maximum Courant number of 0.5 was also run using mesh 4 and the results from this simulation are discussed in Chapter 6. Dynamic meshing and 6DoF was enabled with the model free to move in pitch (rotation about y) and heave (translation in z), and locked in the other degrees of freedom. This allows the solver to dynamically find the correct trim and sinkage.

The implicit Euler scheme was used for all time derivatives in the JBC simulations and the gradients were

calculated using a least-squares scheme. The divergences were discretized using a Gauss scheme and the

interpolation schemes used are summarized in Table 14.

u

α

k

ω

Scheme linearUpwind vanLeer01DC Upwind Upwind

Table 14: The divergence schemes used in all JBC simulations. vanLeer01DC is a combination of upwind and central differencing that uses a flux limiter to ensure that alpha is strictly bounded between 0 and 1. For more information about the schemes refer to the OpenFOAM Users Guide [28].

(34)

4.5 Case 3: MARIN Systematic Series FDS-5

The computational domain for the FDS-5 spanned 7.5

L

in the x-direction, 3

L

in the y-direction and 4

L

in the z-direction resulting in a blockage ratio of roughly 0.07%.

4.5.1 Meshing

Two initial meshes were also made for the FDS-5 as a first step of a grid dependence study. Some properties of the meshes are presented in the tables below.

Max skewness 3.07

∆

(a) Mesh 1

Max skewness 3.53

∆

(b) Mesh 2

With the meshing process from the JBC in mind, having sufficient z-resolution around the free surface was a priority from the beginning when meshing the FDS-5 and a similar blockMesh design was therefore used. Much like for the JBC, emphasis was also put on having a good resolution around the hull and a coarser mesh was allowed along the top of the model. Trapezoidal refinement regions were also used for the FDS-5 to allow the wave pattern to be properly captured while still keeping down the mesh size and these can be seen in Figure 18. Due to the shape of the hull the layer generation was failing in some areas of mesh 1, especially around the sharp part of the bow (see Figure 16). For this reason additional refinement was added in the near hull region of mesh 2 and this lead to successful prism layer generation along most of the geometry. Additional figures of both meshes are available in Appendix A.

(35)

Figure 17: Close up of mesh 1

Acknowledgements k ω Abstract

Evaluation of drag

estimation methods for ship

hulls

Abstract

k

ω

Acknowledgements

Contents

1

Introduction

1.1

Submarines

C

C

C

=

R

ρu

S

C

ρ

u

S

C

x

C

C

=

Z

1.2

Surface Ships

C

= C

+ C

.

=

ρuL

µ

=

uL

ν

=

,

ρ

µ

ν

=

√

u

L

g

=

,

1.3

Ship terminology

∇

∇

1.4

Semi-Empirical Methods

C

= C

− C

C

=

0.075

log(Re) − 2

.

C

= (1 + k)C

,

C

= (1 + k)C

+ C

1.5

Fluid Mechanics

∂ρ

∂t

+

∂(ρu

_u