Effect of Ship Motions on Propeller-Hull

Interaction

by
c Md Ashim Ali

A thesis submitted to the School of Graduate Studies in partial fulfillment of the

requirements for the degree of
Doctor of Philosophy

Department of Ocean and Naval Architectural Engineering

Faculty of Engineering and Applied Science

Memorial University of Newfoundland

February 2020

St. John’s Newfoundland and Labrador, Canada

Abstract

It is important to consider the detailed hydrodynamics of the propeller, the associated

flow conditions and the hull interaction with effective wake distribution for developing
the best ship propulsion system design. Model experiments are conventionally used
to assess propeller flow and performance but are subject to scale effect and difficulties
in flow measurement. Computational Fluid Dynamics (CFD), using a viscous flow
solver, provides an alternate way to simulate propeller flow conditions. Although
simplified propeller flow models are currently available, this thesis extends the CFD
approach to develop a more comprehensive model that includes detailed propeller
geometry, ship motions and wave action.

Numerical self-propulsion simulations with propeller action and ship motions are
performed using a newly developed dynamic motion class based on the sliding mesh
method in OpenFOAM. Validation and verification are carried out in a step by
step manner. As a first step, a model incorporating bare hull resistance, wave
elevation on the hull surface, free surface wave contour, and axial velocity distribution
at several transverse planes for a fixed ship condition and dynamic ship condition
was developed and validated using available experimental data, and the comparison
between the simulations and the experiments showed good agreement. Propeller open
water simulations with detailed propeller geometry and using the body-force method
were developed and validated with experimental data in the second step. Because of

ii
the highly non-linear and rotational nature of the flow near to the propeller region,
curvature correction was implemented to the standard k − ω SST turbulence model
to improve the numerical predictions.

Prior to the third stage self-propulsion simulation, which combined the stage 1 and
2 simulations, a numerical wave tank was developed and validated. Numerical
self-propulsion simulations were performed for a ship at an even keel condition
with constant propeller rotational speed. The numerical calculation of the
skin friction correction, thrust coefficient, torque coefficient, and axial velocity
distribution downstream of the propeller, confirmed good accuracy against the
available experimental data. Then the developed dynamic motion class was extended
to models with ship sinkage and trim. Comparison with available data showed a
excellent agreement for all self-propulsion parameters. An expected increase in total
resistance and decrease in propeller thrust was observed due to ship motion.

As an alternative of detailed propeller geometry modelling, two body-force models

were developed for self-propulsion simulation when detailed flow features are not
required. The effective local velocity body-force model showed better numerical
predictions than nominal average velocity body-force method. The self-propulsion
simulation was performed in waves with forward speed using the developed dynamic
motion class and it was able to incorporate the rotating propeller with large ship
motions. The thrust generated by a propeller varied with encounter wave period and
was reduced compared to the even keel condition.

iii
Acknowledgements

Firstly, I would like to express my sincere gratitude to my supervisor, Dr. Heather

Peng, for her continuous support throughout my PhD studies and research. I am very
grateful for her patience and kindness. Not only her guidance consistently helped
me throughout conducting my research, but also her encouragement has helped me
through the difficult times I have encountered.

Besides my supervisor, I would like to thank the rest of my thesis committee, Dr.
Wei Qiu and Dr. David Molyneux, for their insightful comments. They have helped
me to better understand the field of marine hydrodynamics. Without their precious
support, it would not have been possible to conduct this research. I would also like
to show my gratitude to my colleagues and friends at Memorial University, for the
helpful discussions and suggestions in all the stages of my research.

I am very grateful to the support by NSERC through its CREATE program. Without
this support, I would have not been able to complete this research.

Last but not the least, I would like to thank my parents and my family. I would like
to thank my wife, my son and my daughter, who sacrifices their time throughout my
PhD studies. I also owe my deepest gratitude to my parents for their support and
understanding all the time.

iv
Table of Contents

Abstract ii

Acknowledgments iv

Table of Contents xi

List of Tables xiv

List of Figures xxiii

1 Introduction 1
1.1 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.1.1 Experimental Studies . . . . . . . . . . . . . . . . . . . . . . . 4
1.1.2 Numerical Studies . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.1.3 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . 14
1.1.4 Research Aims and Objectives . . . . . . . . . . . . . . . . . . 14
1.1.5 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2 Mathematical Formulation 19
2.1 Governing Equation for Multiphase Flow . . . . . . . . . . . . . . . . 19
2.1.1 Continuity and Momentum Equations . . . . . . . . . . . . . 20
2.1.1.1 Pressure Term . . . . . . . . . . . . . . . . . . . . . 21

v
 2.1.1.2 Viscous Term . . . . . . . . . . . . . . . . . . . . . . 21
2.1.2 Volume of Fluid . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.1.3 Interface Capturing Scheme . . . . . . . . . . . . . . . . . . . 24
2.1.3.1 interfaceCompression Scheme . . . . . . . . . . . . . 25
2.1.3.2 HRIC Scheme . . . . . . . . . . . . . . . . . . . . . . 26
2.2 Discretization of Governing Equations . . . . . . . . . . . . . . . . . 27
2.2.1 Time Derivative . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.2.2 Convective Term . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.2.3 Diffusion Term . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.2.4 Source Term . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.2.5 Temporal Discretization . . . . . . . . . . . . . . . . . . . . . 30
2.3 Pressure Velocity Coupling . . . . . . . . . . . . . . . . . . . . . . . . 31
2.4 Turbulence Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.4.1 k − ω SST Turbulence Model . . . . . . . . . . . . . . . . . . 33
2.4.2 LRR Turbulence Model . . . . . . . . . . . . . . . . . . . . . 36
2.5 Near-Wall Treatment . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.6 Sliding Mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
2.7 Dynamic Motion Solver . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.8 Body-Force Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
2.8.1 Fixed Velocity Body-Force Model (FV-BFM) . . . . . . . . . 45
2.8.2 Local Velocity Body-Force Model (LV-BFM) . . . . . . . . . . 46
2.9 Convergence Criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
2.9.1 Residual Convergence . . . . . . . . . . . . . . . . . . . . . . . 47
2.9.2 Convergence Error . . . . . . . . . . . . . . . . . . . . . . . . 48
2.10 Numerical Solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
2.11 Verification and Validation Methodology . . . . . . . . . . . . . . . . 50

vi
 2.11.1 Verification . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
2.11.1.1 Least Squares Root (LSR) Method . . . . . . . . . . 50
2.11.2 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
2.12 Numerical Wave Tank . . . . . . . . . . . . . . . . . . . . . . . . . . 53
2.12.1 Wave Generation . . . . . . . . . . . . . . . . . . . . . . . . . 54
2.12.2 Wave Absorption . . . . . . . . . . . . . . . . . . . . . . . . . 55
2.13 Vortex Identification Criteria . . . . . . . . . . . . . . . . . . . . . . . 56
2.13.1 λ2 -criterion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
2.13.2 Q-criterion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
2.14 Self-propulsion Modelling . . . . . . . . . . . . . . . . . . . . . . . . . 58

3 Resistance and Flow Validation 63

3.1 Bare Hull . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.1.1 KCS Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
3.1.1.1 Geometry and Test Conditions . . . . . . . . . . . . 66
3.1.1.2 Computational Domain and Grid Generation . . . . 66
3.1.1.3 Resistance Validation . . . . . . . . . . . . . . . . . 70
3.1.1.4 Free Surface Wave Contour . . . . . . . . . . . . . . 74
3.1.1.5 Flow Velocity . . . . . . . . . . . . . . . . . . . . . . 75
3.1.1.6 Effect of Turbulence Model . . . . . . . . . . . . . . 80
3.1.2 Japan Bulk Carrier . . . . . . . . . . . . . . . . . . . . . . . . 86
3.1.2.1 Geometry and Test Conditions . . . . . . . . . . . . 86
3.1.2.2 Computational Domain and Grid Generation . . . . 86
3.1.2.3 Resistance Validation . . . . . . . . . . . . . . . . . 88
3.1.2.4 Free Surface Wave Contour . . . . . . . . . . . . . . 89
3.1.2.5 Flow Velocity . . . . . . . . . . . . . . . . . . . . . . 96
3.1.3 Fishing Vessel . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

vii
 3.1.3.1 Geometry and Test Conditions . . . . . . . . . . . . 98
3.1.3.2 Computational Domain and Grid Generation . . . . 99
3.1.3.3 Resistance Validation . . . . . . . . . . . . . . . . . 101
3.1.3.4 Free Surface Wave Contour . . . . . . . . . . . . . . 105
3.1.3.5 Flow Velocity . . . . . . . . . . . . . . . . . . . . . . 106
3.2 Hull with Appendages . . . . . . . . . . . . . . . . . . . . . . . . . . 108
3.2.1 KCS Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
3.2.1.1 Geometry and Test Conditions . . . . . . . . . . . . 108
3.2.1.2 Computational Domain and Grid Generation . . . . 108
3.2.1.3 Resistance, Sinkage, and Trim Validation . . . . . . 111
3.3 Scale Effect on Resistance, Sinkage and Trim . . . . . . . . . . . . . . 115
3.3.1 Computational Domain and Grid Generation . . . . . . . . . . 115
3.3.2 Resistance, Sinkage and Trim . . . . . . . . . . . . . . . . . . 117
3.3.3 Free Surface Wave Contour . . . . . . . . . . . . . . . . . . . 119
3.3.4 Flow Velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

4 Propeller Open Water Validation 123

4.1 KP505 Test Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
4.1.1 Geometry and Test Conditions . . . . . . . . . . . . . . . . . 124
4.1.2 Actual Propeller Geometry . . . . . . . . . . . . . . . . . . . . 124
4.1.2.1 Computational Domain and Grid Generation . . . . 125
4.1.2.2 Numerical Results . . . . . . . . . . . . . . . . . . . 127
4.1.2.3 Steady Simulation . . . . . . . . . . . . . . . . . . . 127
4.1.2.4 Unsteady Simulation . . . . . . . . . . . . . . . . . . 130
4.1.3 Body-Force Method . . . . . . . . . . . . . . . . . . . . . . . . 131
4.1.3.1 Computational Domain and Grid Generation . . . . 131
4.1.3.2 Grid Dependency Study . . . . . . . . . . . . . . . . 132

viii
 4.1.3.3 Sampling Plane Distance . . . . . . . . . . . . . . . . 133
4.2 MP687 Test Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
4.2.1 Propeller Geometry . . . . . . . . . . . . . . . . . . . . . . . . 136
4.2.1.1 Computational Domain and Grid Generation . . . . 136
4.2.1.2 Steady Simulation . . . . . . . . . . . . . . . . . . . 137
4.2.1.3 Unsteady Simulation . . . . . . . . . . . . . . . . . . 139
4.3 NRC-IOT Propeller Test Case . . . . . . . . . . . . . . . . . . . . . . 143
4.3.1 Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
4.3.2 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . 144

5 Numerical Wave Tank Validation 147

5.1 Computational Domain . . . . . . . . . . . . . . . . . . . . . . . . . . 148
5.2 Grid Dependency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149
5.3 Time Step Size Dependency . . . . . . . . . . . . . . . . . . . . . . . 150
5.4 Wave Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152

6 Self-Propulsion Simulation 153

6.1 KCS Test Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153
6.1.1 Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155
6.1.2 Computational Domain and Boundary Conditions . . . . . . . 155
6.1.3 Self-Propulsion at Restrained Even Keel Condition . . . . . . 156
6.1.3.1 Pressure Distribution . . . . . . . . . . . . . . . . . . 159
6.1.3.2 Flow Velocity . . . . . . . . . . . . . . . . . . . . . . 159
6.1.4 Self-Propulsion with Sinkage and Trim . . . . . . . . . . . . . 161
6.1.4.1 Pressure Distribution . . . . . . . . . . . . . . . . . . 163
6.1.4.2 Flow Velocity . . . . . . . . . . . . . . . . . . . . . . 165
6.1.5 Seakeeping Simulation . . . . . . . . . . . . . . . . . . . . . . 169

ix
 6.1.6 Self-Propulsion in Waves . . . . . . . . . . . . . . . . . . . . . 171
6.1.6.1 Ship Response . . . . . . . . . . . . . . . . . . . . . 174
6.1.6.2 Free Surface Wave . . . . . . . . . . . . . . . . . . . 174
6.1.6.3 Pressure Distribution . . . . . . . . . . . . . . . . . . 177
6.2 JBC Test Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179
6.2.1 Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179
6.2.2 Computational Domain and Grid Generation . . . . . . . . . . 180
6.2.3 Results Validation . . . . . . . . . . . . . . . . . . . . . . . . 180
6.2.4 Flow Velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . 183
6.2.5 Vortex Visualization . . . . . . . . . . . . . . . . . . . . . . . 190
6.3 Fishing Vessel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191
6.3.1 Wave Elevations . . . . . . . . . . . . . . . . . . . . . . . . . . 193
6.3.2 Velocities on Propeller Plane . . . . . . . . . . . . . . . . . . . 194
6.3.3 Thrust, Torque, Tow Force and Self-Propulsion Point . . . . . 197
6.3.4 Self-Propulsion Results . . . . . . . . . . . . . . . . . . . . . . 198
6.3.5 Pressure Distribution . . . . . . . . . . . . . . . . . . . . . . . 199
6.4 HPC and Computational Time . . . . . . . . . . . . . . . . . . . . . 202

7 Conclusions and Recommendations 203

7.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203
7.2 Best practice for propeller-hull interaction. . . . . . . . . . . . . . . . 209
7.2.1 Selection of computational domain. . . . . . . . . . . . . . . . 209
7.2.2 Mesh generation. . . . . . . . . . . . . . . . . . . . . . . . . . 210
7.2.3 Selection of propeller modelling . . . . . . . . . . . . . . . . . 211
7.2.4 Solver selection . . . . . . . . . . . . . . . . . . . . . . . . . . 211
7.2.5 Boundary conditions set-up . . . . . . . . . . . . . . . . . . . 211
7.2.6 Selection of turbulence model . . . . . . . . . . . . . . . . . . 212

x
 7.2.7 Selection of numerical schemes . . . . . . . . . . . . . . . . . . 212
7.2.8 Time step size selection . . . . . . . . . . . . . . . . . . . . . 213
7.2.9 Setting convergence criteria . . . . . . . . . . . . . . . . . . . 213
7.3 Recommendations for Future Work . . . . . . . . . . . . . . . . . . . 213

Bibliography 215

Appendices 227

A Propeller Tip Vortex Simulation 228

A.1 Geometry and Test Conditions . . . . . . . . . . . . . . . . . . . . . . 228
A.2 Computational Domain and Grid generation . . . . . . . . . . . . . . 229
A.3 Numerical Results Validation . . . . . . . . . . . . . . . . . . . . . . 231

B Wave Run-up Simulation 247

B.1 Model Experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . 247
B.2 Computational Domain and Grid Generation . . . . . . . . . . . . . . 248
B.3 Numerical Set-up . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 250
B.4 Wave Generation Modules . . . . . . . . . . . . . . . . . . . . . . . . 250
B.5 Grid Resolution Study . . . . . . . . . . . . . . . . . . . . . . . . . . 252
B.6 Time Step Dependency . . . . . . . . . . . . . . . . . . . . . . . . . . 252
B.7 Effect of Divergence Scheme . . . . . . . . . . . . . . . . . . . . . . . 253
B.8 Effect of Turbulence Model . . . . . . . . . . . . . . . . . . . . . . . . 257
B.9 Wave Run-up Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 257

xi
List of Tables

3.1 Principal particulars of the KCS bare hull. . . . . . . . . . . . . . . . 66

3.2 Test case for bare hull resistance . . . . . . . . . . . . . . . . . . . . 67
3.3 Boundary conditions for RANS simulations . . . . . . . . . . . . . . . 69
3.4 Computational grid for bare hull simulation. . . . . . . . . . . . . . . 69
3.5 Steady resistance and uncertainty of the KCS bare hull at F n = 0.26. 69
3.6 Effect of turbulence model on total resistance for GRID-C . . . . . . 80
3.7 RMS error of the velocity components for different turbulence models
for GRID-C. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
3.8 Principal particulars of the JBC hull. . . . . . . . . . . . . . . . . . . 87
3.9 JBC test cases considered for validation of numerical results. . . . . . 87
3.10 Computational grid for JBC bare hull simulation. . . . . . . . . . . . 88
3.11 Grid dependency of CT , sinkage and trim at F n = 0.142. . . . . . . . 92
3.12 Comparison of | USN %S | for CT , sinkage and trim at F n = 0.142. . . 92
3.13 Principal particulars of the ship model and the propeller. . . . . . . . 99
3.14 Test conditions for fishing vessel model. . . . . . . . . . . . . . . . . . 99
3.15 Computational grids for simulations of fishing vessel. . . . . . . . . . 100
3.16 Comparison of total resistance coefficient, CT . . . . . . . . . . . . . . 101
3.17 Comparison of |E%D| for CT , sinkage and trim at F n = 0.34 . . . . 103
3.18 Principal particulars of rudder for KCS hull. . . . . . . . . . . . . . . 109

xii
3.19 Test conditions for appended KCS model. . . . . . . . . . . . . . . . 109
3.20 Computational grid for KCS with rudder. . . . . . . . . . . . . . . . 111
3.21 Grid dependency of CT , sinkage and trim at F n = 0.26 . . . . . . . . 111
3.22 Comparison of | USN %S | for CT , sinkage and trim at F n = 0.26. . . 111
3.23 Results of unsteady simulation of KCS hull with rudder using GRID-C. 112
3.24 Computational grid for scale effect study. . . . . . . . . . . . . . . . . 117

4.1 Principal particulars of KP505 propeller. . . . . . . . . . . . . . . . . 125

4.2 Computational grid for KP505 propeller open water simulation . . . . 127
4.3 Grid dependency of open water performance at J = 0.7 using steady
solver. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
4.4 Grid dependency of open water performance at J = 0.7 using unsteady
solver. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
4.5 Principal particulars of MP687 propeller. . . . . . . . . . . . . . . . . 137
4.6 Computational grids for MP687 propeller open water simulation . . . 138
4.7 Grid dependency of open water performance at J = 0.7 using steady
solver. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
4.8 Grid dependency of open water performance at J = 0.7 using unsteady
solver. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
4.9 Principal particulars of the NRC-IOT propeller. . . . . . . . . . . . . 143

5.1 Wave conditions for numerical tank validation . . . . . . . . . . . . . 148

6.1 Test cases for the validation of self-propulsion simulation of KCS model.154
6.2 Numerical results of self-propulsion simulation. . . . . . . . . . . . . . 158
6.3 Numerical results of self-propulsion test with ship motion. . . . . . . 161
6.4 Numerical results of self-propulsion test with ship motion in wave C3. 171
6.5 Test conditions for JBC self-propulsion simulation . . . . . . . . . . . 179

xiii
6.6 Numerical results of self-propulsion of JBC model. . . . . . . . . . . . 181
6.7 Comparison of |E%D| for FT ow , sinkage and trim at F n = 0.34 . . . 195
6.8 Self-propulsion simulation results using the load varying method by
FV-BFM at F n = 0.34 . . . . . . . . . . . . . . . . . . . . . . . . . . 196
6.9 Self-propulsion simulation results using the load varying method by
LV-BFM at F n = 0.34 . . . . . . . . . . . . . . . . . . . . . . . . . . 196
6.10 Self-propulsion parameters of fishing vessel at F n = 0.34 . . . . . . . 198
6.11 Computational time in skipper and openfoam server. . . . . . . . . . 202

A.1 Propeller particulars . . . . . . . . . . . . . . . . . . . . . . . . . . . 228

A.2 Computational grids for DTM5168 propeller. . . . . . . . . . . . . . . 232
A.3 Boundary conditions for tip vortex simulations. . . . . . . . . . . . . 232
A.4 Comparison of KT and 10KQ . . . . . . . . . . . . . . . . . . . . . . 233
A.5 RMS error of Vx , Vt and Vr at x/R = 0.1756 . . . . . . . . . . . . . . 241
A.6 RMS error of Vx , Vt and Vr at x/R = 0.2386 . . . . . . . . . . . . . . 241

B.1 Wave probe locations. . . . . . . . . . . . . . . . . . . . . . . . . . . 248

B.2 Wave matrix for wave run-up simulation. . . . . . . . . . . . . . . . 249
B.3 Boundary conditions for numerical wave generation. . . . . . . . . . . 250

xiv
List of Figures

1.1 Contour of phase-averaged axial velocity in the longitudinal plane at

the phase angle φ = 00 . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Axial velocity comparison in front of CP and CRP . . . . . . . . . . . 7
1.3 The non-dimensional axial velocity at x/R = −0.514 . . . . . . . . . 8
1.4 The non-dimensional axial velocity at r/R = 0.7. . . . . . . . . . . . 8
1.5 Numerical wake contours behind the ship hull of a LPG carrier at
F n = 0.242 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.1 Calculation of cell face volume fraction. . . . . . . . . . . . . . . . . . 25

2.2 Parameters for finite control volume discretization (Rusche, 2002). . . 27
2.3 PIMPLE algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.4 Law of the wall . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.5 Illustration of the sliding mesh method. . . . . . . . . . . . . . . . . . 40
2.6 Dynamic update of ship and propeller domains. . . . . . . . . . . . . 42
2.7 Calculation procedure using dynamic motion class. . . . . . . . . . . 43
2.8 Example of iterative convergence for variable x. . . . . . . . . . . . . 48
2.9 Schematic of numerical wave tank . . . . . . . . . . . . . . . . . . . . 55
2.10 Variation of γR (χR ) for damping zone . . . . . . . . . . . . . . . . . 56

3.1 Body plan and profile elevation of KCS. . . . . . . . . . . . . . . . . 65

xv
3.2 KCS bare hull geometry. . . . . . . . . . . . . . . . . . . . . . . . . . 66
3.3 Extent of the computational domain for bare hull simulation. . . . . . 67
3.4 Grid distribution around the ship hull. . . . . . . . . . . . . . . . . . 68
3.5 Grid refinement of Kelvin wedge on the free surface. . . . . . . . . . . 68
3.6 Convergence of total resistance coefficient, CT for GRID-C. . . . . . . 70
3.7 Residual convergence for GRID-C. . . . . . . . . . . . . . . . . . . . . 71
3.8 Convergence error for total resistance coefficient for GRID-C at F n =
0.26. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
3.9 Numerical uncertainty of CT at F n = 0.26. . . . . . . . . . . . . . . . 72
3.10 Wave elevation along the ship hull at F n = 0.26. . . . . . . . . . . . . 72
3.11 Free surface wave pattern comparison at F n = 0.26. . . . . . . . . . . 73
3.12 Wave elevation at y/LP P = 0.1024 . . . . . . . . . . . . . . . . . . . . 74
3.13 Wave elevation at y/LP P = 0.1509 . . . . . . . . . . . . . . . . . . . . 75
3.14 Wave elevation at y/LP P = 0.30 . . . . . . . . . . . . . . . . . . . . . 76
3.15 Wave elevation at y/LP P = 0.40 . . . . . . . . . . . . . . . . . . . . . 76
3.16 Axial flow velocity contour on the propeller plane at F n = 0.26. . . . 77
3.17 Axial flow velocity on the propeller plane at z/LP P = −0.03. . . . . . 78
3.18 Transverse flow velocity on the propeller plane at z/LP P = −0.03. . . 78
3.19 Vertical flow velocity on the propeller plane at z/LP P = −0.03. . . . 79
3.20 Effect of turbulence models on wave elevation along the hull surface for
GRID-C. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
3.21 Effect of turbulence models on free surface wave pattern for GRID-C. 82
3.22 Axial flow velocity contour on the propeller plane using different
turbulence models for GRID-C. . . . . . . . . . . . . . . . . . . . . . 83
3.23 Effect of turbulence model on axial flow velocity on the propeller plane
at z/LP P = −0.03 for GRID-C. . . . . . . . . . . . . . . . . . . . . . 84

xvi
3.24 Effect of turbulence model on transverse flow velocity on the propeller
plane at z/LP P = −0.03 for GRID-C. . . . . . . . . . . . . . . . . . . 84
3.25 Effect of turbulence model on vertical flow velocity on the propeller
plane at z/LP P = −0.03 for GRID-C. . . . . . . . . . . . . . . . . . . 85
3.26 JBC bare hull geometry. . . . . . . . . . . . . . . . . . . . . . . . . . 86
3.27 Computational domain of JBC hull. . . . . . . . . . . . . . . . . . . . 89
3.28 Grid distribution around the JBC hull. . . . . . . . . . . . . . . . . . 89
3.29 Convergence of CT for GRID-C at F n = 0.142. . . . . . . . . . . . . . 90
3.30 Time history of residuals convergence for GRID-C. . . . . . . . . . . 90
3.31 Convergence error of CT for GRID-C at F n = 0.142. . . . . . . . . . 91
3.32 Grid uncertainty for CT at F n = 0.142. . . . . . . . . . . . . . . . . . 91
3.33 Wave elevation along the ship hull. . . . . . . . . . . . . . . . . . . . 92
3.34 Free surface wave pattern comparison at F n = 0.142. . . . . . . . . . 93
3.35 Wave elevation at y/LP P = 0.1043 at F n = 0.142. . . . . . . . . . . . 94
3.36 Wave elevation at y/LP P = 0.1900 at F n = 0.142. . . . . . . . . . . . 94
3.37 Axial flow velocity contour on x/LP P = 0.9625 for F n = 0.142. . . . . 95
3.38 Axial flow velocity contour on x/LP P = 0.9843 for F n = 0.142. . . . . 96
3.39 Axial flow velocity contour on x/LP P = 1.0 for F n = 0.142. . . . . . . 97
3.40 Ship model of optimized fishing vessel. . . . . . . . . . . . . . . . . . 98
3.41 Computational domain for Fishing vessel. . . . . . . . . . . . . . . . . 100
3.42 Grid distribution around the bulbous bow. . . . . . . . . . . . . . . . 101
3.43 Total resistance coefficient for F n = 0.34. . . . . . . . . . . . . . . . . 102
3.44 Total resistance coefficient of Model-H using GRID-C. . . . . . . . . . 103
3.45 Free surface wave contours for Model-H at F n = 0.34 . . . . . . . . . 104
3.46 Bow waves of Model-H at F n = 0.34 . . . . . . . . . . . . . . . . . . 105
3.47 Comparison of wave elevation along the hull of Model-H at F n = 0.34 106

xvii
3.48 Axial velocity contour on the propeller plane at F n = 0.34 for GRID-C. 107
3.49 KCS model with rudder. . . . . . . . . . . . . . . . . . . . . . . . . . 108
3.50 Semi-balanced horn rudder for KCS. . . . . . . . . . . . . . . . . . . 110
3.51 Grid distribution for sinkage and trim computation. . . . . . . . . . . 110
3.52 Comparison of total resistance coefficient using GRID-C. . . . . . . . 112
3.53 Comparison of sinkage of KCS model using GRID-C. . . . . . . . . . 113
3.54 Comparison of trim of KCS model GRID-C. . . . . . . . . . . . . . . 113
3.55 Computational domain for the study of scale effect. . . . . . . . . . . 115
3.56 Grid distribution at bow and stern of the model. . . . . . . . . . . . . 116
3.57 Effect of scale factor on frictional resistance. . . . . . . . . . . . . . . 118
3.58 Effect of scale factor on pressure resistance. . . . . . . . . . . . . . . 118
3.59 Effect of scale factor on sinkage. . . . . . . . . . . . . . . . . . . . . . 119
3.60 Effect of scale factor on trim. . . . . . . . . . . . . . . . . . . . . . . 120
3.61 Effect of scale factor of free surface wave contour. . . . . . . . . . . . 121
3.62 Effect of scale factor on axial velocity distribution on propeller plane. 122

4.1 Geometry of KP505 propeller. . . . . . . . . . . . . . . . . . . . . . . 125

4.2 Computational domain of KP505 propeller open water simulation. . . 126
4.3 Grid distribution on KP505 propeller blade. . . . . . . . . . . . . . . 126
4.4 Convergence of KT , 10KQ for J = 0.7 using steady solver for GRID-B. 128
4.5 Residual convergence using steady solver for GRID-B at J = 0.7. . . . 128
4.6 KP505 propeller open water simulation using steady solver for GRID-B. 129
4.7 Convergence of KT , 10KQ for J = 0.7 using transient solver for
GRID-B. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
4.8 Residual convergence using transient solver for GRID-B at J = 0.7. . 131
4.9 Open water simulation of KP505 using transient solver for GRID-B. . 132

xviii
4.10 Computational domain for KP505 open water simulation using
body-force method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
4.11 Effect of number of cells in radial direction. . . . . . . . . . . . . . . 133
4.12 Effect of velocity sampling plane distance. . . . . . . . . . . . . . . . 134
4.13 Propeller open water characteristics of KP505 using body-force method. 134
4.14 MP687 propeller geometry. . . . . . . . . . . . . . . . . . . . . . . . . 136
4.15 Grid distribution on the propeller blade and hub of MP687. . . . . . 137
4.16 Convergence history of KT , KQ and η0 using steady solver for GRID-B
at J = 0.70. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
4.17 Residuals convergence for advance coefficient J = 0.70 for GRID-B. . 139
4.18 Open water simulation of MP687 propeller using steady solver for
GRID-B. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
4.19 Convergence history of KT , KQ and η0 using transient solver for
GRID-B at J = 0.70. . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
4.20 Residuals convergence for GRID-B at advance coefficient J = 0.70. . . 141
4.21 Open water simulation of MP687 propeller using transient solver for
GRID-B. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
4.22 Effect of the cell numbers on propeller parameters at J = 0.70 using
transient solver. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
4.23 Effect of the velocity sampling plane distance on propeller parameters
at J = 0.70 using transient solver. . . . . . . . . . . . . . . . . . . . . 144
4.24 Comparison of open water data of NRC-IOT propeller using transient
solver. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145

5.1 Computational domain for numerical wave tank validation. . . . . . . 148

5.2 Grid distribution in the computational domain. . . . . . . . . . . . . 149
5.3 Comparison of wave elevation at WP-02 using TS03 of wave C1. . . . 150

xix
5.4 Comparison of wave elevation at WP-02 using GRID-C of wave C1. . 151
5.5 Comparison of wave elevation at WP-02 using GRID-C of wave C3. . 151

6.1 Geometry for KCS self-propulsion. . . . . . . . . . . . . . . . . . . . 155

6.2 Computational grid for self-propulsion of KCS model. . . . . . . . . . 155
6.3 Dynamic pressure distribution on the stern of KCS for self-propulsion. 157
6.4 Axial velocity distribution at 0.25D downstream of the propeller plane. 158
6.5 Velocity components comparison at x/LP P = 0.9910 (0.25D
downstream of the propeller plane) and z/LP P = -0.0291. . . . . . . . 160
6.6 Residual convergence for self-propulsion with sinkage and trim. . . . 162
6.7 Total resistance convergence error for self-propulsion with sinkage and
trim. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162
6.8 Propeller thrust convergence error for self-propulsion with sinkage and
trim. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163
6.9 Pressure distribution on the hull and the propeller. . . . . . . . . . . 164
6.10 Virtual pressure gauge locations. . . . . . . . . . . . . . . . . . . . . 165
6.11 Numerical pressure measurement at the pressure gauges. . . . . . . . 166
6.12 Axial velocity distribution on the centreplane. . . . . . . . . . . . . . 167
6.13 Axial velocity distribution at x/LP P = 0.9910 (0.25D downstream of
the propeller plane). . . . . . . . . . . . . . . . . . . . . . . . . . . . 168
6.14 Comparison of total resistance coefficient in wave C3. . . . . . . . . . 169
6.15 Motion response of KCS bare hull in wave C3. . . . . . . . . . . . . 170
6.16 Residual convergence for self-propulsion in wave C3. . . . . . . . . . . 172
6.17 Total resistance coefficient for self-propulsion in wave C3. . . . . . . . 172
6.18 Thrust coefficient over one encounter wave period for self-propulsion in
wave C3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173

xx
6.19 Torque coefficient over one encounter wave period for self-propulsion
in wave C3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173
6.20 Motion response of KCS hull during self-propulsion in wave C3. . . . 175
6.21 Free surface wave contour over one encounter wave period, Te for
self-propulsion in wave C3. . . . . . . . . . . . . . . . . . . . . . . . . 176
6.22 Pressure distribution on hull surface over one encounter wave period,
Te for self-propulsion in wave C3. . . . . . . . . . . . . . . . . . . . . 177
6.23 Geometry for JBC model for self-propulsion simulation. . . . . . . . . 180
6.24 Computational grid for self-propulsion of JBC model. . . . . . . . . . 181
6.25 Residual convergence of self-propulsion simulation of JBC model. . . 182
6.26 Convergence of total resistance of JBC model. . . . . . . . . . . . . . 182
6.27 Blade angle definition. . . . . . . . . . . . . . . . . . . . . . . . . . . 183
6.28 Axial velocity distribution at x/LP P = 0.9625 . . . . . . . . . . . . . 185
6.29 Axial velocity distribution at x/LP P = 0.9843 for blade angle 00 . . . 186
6.30 Axial velocity distribution at x/LP P = 0.9843 for blade angle 480 . . 187
6.31 Axial velocity distribution at x/LP P = 1.00 for blade angle 00 . . . . . 188
6.32 Axial velocity distribution at x/LP P = 1.00 for blade angle 480 . . . . 189
6.33 Vortex structures based on Q = 500 for the self-propulsion simulation
with sinage and trim at F n = 0.142. . . . . . . . . . . . . . . . . . . 190
6.34 Computational grid for self-propulsion of Fishing vessel. . . . . . . . . 191
6.35 Comparison of free surface wave contour for GRID-C at F n = 0.34 in
calm water. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192
6.36 Comparison of wave elevation along the hull for GRID-C at F n = 0.34. 193
6.37 Input velocity distribution on the propeller plane at F n = 0.34. . . . 194
6.38 Axial velocity distribution on the propeller plane at F n = 0.34. . . . 195
6.39 Effective wake distribution on the propeller plane at F n = 0.34. . . . 196

xxi
6.40 Propeller suction effect at x/Lpp = 0.10D at F n = 0.34. . . . . . . . . 197
6.41 Comparison of tow forces at F n = 0.34. . . . . . . . . . . . . . . . . . 199
6.42 Pressure distribution on the hull at F n = 0.34. . . . . . . . . . . . . . 200
6.43 Streamlines around the propeller plane for F n = 0.34 . . . . . . . . . 201

A.1 Geometry and computation domain for grid generation . . . . . . . . 230

A.2 Type of refinement. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 230
A.3 Refinement level. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 231
A.4 Comparison of KT and 10KQ with experimental data (Chesnakas and
Jessup (1998)) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233
A.5 Axial velocity Vx contour at x/R = 0.2386 . . . . . . . . . . . . . . . 234
A.6 Tangential velocity Vt contour at x/R = 0.2386 . . . . . . . . . . . . 235
A.7 Radial velocity Vr contour at x/R = 0.2386 . . . . . . . . . . . . . . . 236
A.8 Vx , Vt and Vr across vortex core at x/R = 0.1756 . . . . . . . . . . . 237
A.9 Vx , Vt and Vr across vortex core at x/R = 0.1756 . . . . . . . . . . . 238
A.10 Vx , Vt and Vr across vortex core at x/R = 0.2386 . . . . . . . . . . . 239
A.11 Vx , Vt and Vr across vortex core at x/R = 0.2386 . . . . . . . . . . . 240
A.12 Vx , Vt and Vr across vortex core at x/R = 0.1756 . . . . . . . . . . . 242
A.13 Vx , Vt and Vr across vortex core at x/R = 0.2386 . . . . . . . . . . . 243
A.14 Pressure coefficient on suction side. . . . . . . . . . . . . . . . . . . . 244
A.15 Pressure coefficient of pressure side. . . . . . . . . . . . . . . . . . . . 245
A.16 Vortex trajectory for GRID-C3. . . . . . . . . . . . . . . . . . . . . . 246

B.1 Wave probes locations . . . . . . . . . . . . . . . . . . . . . . . . . . 248

B.2 Schematic diagram of computational domain for wave tank validation. 249
B.3 The computational grid for wave run-up simulation. . . . . . . . . . . 249

xxii
B.4 Comparison of wave elevation at x/λ = 0 for T15S110 using Tw /∆t =
2000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 251
B.5 Comparison of wave elevation at x/λ = 0 for T15S110 using Tw /∆t =
2000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 251
B.6 Comparison of wave elevation at x/λ = 0 for T15S110 using Tw /∆t =
2000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253
B.7 Comparison of wave elevation x/λ = 0 for T15S116 using Tw /∆t = 2000.254
B.8 Comparison of wave elevation at x/λ = 0 for T15S110 using GR-03. . 254
B.9 Comparison of wave elevation at x/λ = 0 for T15S116 using GR-03. . 255
B.10 Comparison of wave elevation at x/λ = 0 for T12S110 using GR-03. . 255
B.11 Comparison of wave elevation at x/λ = 0 for T12S116 using GR-03. . 256
B.12 Effect of divergence scheme on wave generation. . . . . . . . . . . . . 256
B.13 Effect of turbulence model on wave generation. . . . . . . . . . . . . . 257
B.14 Time history of wave force on the cylinder for wave T12S110. . . . . . 258
B.15 Time history of wave force on the cylinder for wave T15S110. . . . . . 259
B.16 Wave elevation comparison for wave T12S110 at A2 and A4. . . . . . 260
B.17 Wave elevation comparison for wave T12S110 at B2 and B4. . . . . . 260
B.18 Wave elevation comparison for wave T12S110 at C2 and C4. . . . . . 261
B.19 Instantaneous free surface for wave T12S110. . . . . . . . . . . . . . . 262
B.20 Wave elevation comparison for wave T12S116 at A2. . . . . . . . . . . 263
B.21 Wave elevation comparison for wave T12S116 at A4. . . . . . . . . . . 263
B.22 Wave elevation comparison for wave T12S116 at B2. . . . . . . . . . 264
B.23 Wave elevation comparison for wave T12S116 at B4. . . . . . . . . . 264
B.24 Wave elevation comparison for wave T12S116 at C2. . . . . . . . . . 265
B.25 Wave elevation comparison for wave T12S116 at C4. . . . . . . . . . 265
B.26 Turbulence viscosity distribution on the centre plane. . . . . . . . . . 266

xxiii
Chapter 1

Introduction

High propulsive efficiency, minimum noise, and low vibrations are desirable features
in any ship design. All these features are associated with the propeller and the flow
characteristics around the propeller. For merchant ships, the propulsive efficiency
is most important, whereas, for naval vessels, minimum vibration and noise level
are more important. It is essential to consider the detailed hydrodynamics of the
propeller, the associated flow conditions, and the propeller-hull interaction for the
best propulsion system design. Model experiments are conventionally used to assess
propeller flow and performance but these are subject to scale effect. The flow pattern
in front of the propeller also differs from the model scale to the full scale due to
the different Reynolds Numbers. The current development of Computational Fluid
Dynamics (CFD) enables researchers to simulate complex flow conditions in model
scale as well as in full scale. The propeller-hull interaction can be predicted by
numerical simulation using a viscous flow solver considering the free surface and
detailed propeller geometry. Although a simplified propeller model can be used to
model the propeller action, this neglects the flow alteration by the propeller itself.
With the increasing computational power, it is becoming possible to use detailed

1
propeller geometry in numerical propeller-hull interaction studies to better model a
very complex flow.

In the numerical simulation of ships moving in calm water or in waves, the propeller
effect can be modelled using either detailed propeller geometry or a simplified propeller
model. The most accurate method to simulate propeller-hull interaction in waves
is to model the detailed propeller geometry and hull in viscous flow with the free
surface. For a simplified propeller model, the propulsion effect can be simulated using
the body-force method. The propeller effect can also be obtained using a potential
flow solver. Although the wake distribution behind the ship hull can be modelled
with a certain level of accuracy, the potential flow model is unable to predict the
effective wake due to flow alteration by the propeller itself. Furthermore, when a
ship experiences six degree of freedom motions, the interaction and flow between the
propeller and hull is expected to be further influenced due to the load variation on
the propeller caused by generated waves and ship motions. Also, for ships in waves,
the effective wake varies depending on the position of the ship in the waves, which
also influences the propeller performance. Research on the propeller-hull interaction
in waves is essential to better ship design.

In realistic sea states, the motions of the hull create unfavourable wake and separated
flow into the propeller, which reduces the propeller performance. The propeller
loading conditions also change when the ship moves through waves due to the change
of propeller submergence with respect to the free surface and the large motion of
the vessel itself. In a propulsion system design, knowledge of the effects of these
important factors would improve design for real-life operations. In the numerical
simulation, the rotation of a propeller coupled to ship motions can be simulated using
either the overset grid method or the sliding mesh method. In terms of accuracy and

2
computational time, the sliding mesh method is a good option for detailed propeller
simulation coupled to ship motion.

This thesis details the development of a full numerical model of a working propeller
on a ship undergoing forward motion in a wave field with the resultant six degree of
freedom motions of the ship, using the open source viscous flow solver OpenFOAM.
This development provides a full flow model of a working propeller on a moving ship
and allows the realistic flow conditions of such a propeller to be modelled and assessed.

Due to the very complex nature of the final simulation, a step-wise development and
validation process is followed where the model is built up in stages, starting with the
model of a hull undergoing forward motion, and then adding in the propeller action,
wave action, and finally ship with motions. Validation is performed at each stage using
experimental data available in the open literature. According to the International
Towing Tank Conference (ITTC) procedures (ITTC, 1978), three sets of model tests
are required for power performance prediction behind the hull, including bare-hull
resistance tests, propeller open water tests, and self-propulsion tests. In reality, it
is very difficult to conduct the experimental flow measurements around a working
propeller on a fully moving ship model. Whereas, CFD provides the opportunity for
numerical simulation of self-propulsion tests covering the complex flow conditions of
the full simulation.

1.1 Literature Review

The following review of the literature covers the current extent of available
experimental and numerical works on propulsion systems and pays particular
attention to the studies of propeller-hull interaction.

3
1.1.1 Experimental Studies

Extensive experimental results from benchmark ship models have been developed over
the years for the validation of numerical methods for ship resistance, wave elevation on
the hull surface, free surface wave contour, and wave elevation at different transverse
and longitudinal positions. Self-propulsion data were also developed for different flow
parameters to facilitate the validation of numerical self-propulsion simulations. The
effect of propeller suction on the effective wake was studied by Nagamatsu (1975) using
flow visualization by an air injection method. A considerable change was observed
in the effective wake due to the propeller suction effect for a tanker model and a
cargo liner model. Fujisawa et al. (2000) studied the resistance, wave elevation, and
mean velocity for the KCS model, with and without a rotating propeller, at the even
keel condition, in the towing tank of the National Maritime Research Institute in
Japan (NMRI) following ITTC guidelines. A similar ship model was used by Kim
et al. (2001) for the study of the bare hull resistance, mean flow velocity, and the
free surface elevation of KRISO Container Ship (KCS) and two other vessels KVLCC
and KVLCC2 experimentally in the towing tank of the Korea Research Institute for
Ships and Ocean Engineering (KRISO). The unsteady wave pattern on the free surface
behind the transom was reported at the design speed. A difference between the NMRI
and KRISO results was reported due to a small difference in the ship model geometry.

Propeller-hull interaction was studied by Felli and Felice (2004) using LDV phase
sampling techniques for a twin screw vessel. The wake was measured on the upstream
section of the propeller disk and another section behind the rudder. From the
experimental results, the propeller induced suction was found significant in the region
of x/R = 0.44 to 0.68. This interference led to the peak velocity at the entrance of the
propeller. The effect of propeller induced suction was negligible in the downstream

4
Figure 1.1: Contour of phase-averaged axial velocity in the longitudinal plane at the
phase angle φ = 00 .

direction for distances over x/R = 0.68. The viscous wake generated by appendages
such as shaft brackets also influences the incoming flow to the propeller. The
asymmetric wake distribution was observed which leads to an unbalanced loading
on the propeller, which causes vibrations, noise, and fatigue stress on the rotating
shaft. Velocity reduction in the flow was also observed due to appendage vortices in
the two stern sections. The deformation of the wake due to the combined effect of
the propeller tip and vortices was observed on the downstream section. The vortices
were found stretched on one surface of the rudder, which suffered vibration arising
from interaction with the rudder tip vortex.

Paik et al. (2004, 2007) also investigated the effect of a rotating propeller behind the
KCS model using a two-frame particle image velocimetry (PIV) technique. In the
experiment, the origin was positioned at the center of the propeller boss. The Z−axis
was set along the propeller shaft, pointing toward the stern of the ship; the X−axis was

5
set horizontally toward the starboard, and the Y −axis points upward. The difference
between the incoming flow above the propeller axis and below the propeller axis is
shown in Figure 1.1. The non-dimensionalized axial velocity component (W/W0 ) is
shown in the longitudinal plane for phase angle φ = 00 . The phase angle, (φ) is the
angle (deg) measured from the vertical axis which is pointed in the upward direction.
The distances are normalized by the propeller diameter D in this figure. It was
found that the axial and tangential velocities in the propeller plane were modified
by the hull wake and bilge vortices. The hull wake and free surface influenced the
axial velocity distribution in the propeller plane, which also affected the tip vortex
propagation in the downstream direction. The blade-to-blade interaction led to the
trailing vortex rotating faster than the tip vortex in the longitudinal direction Z/D
up to 0.5. An asymmetric inflow structure was observed due to the bilge vortex
and hull wake interaction, which could be a source of cavitation, noise, and pressure
fluctuations on the hull.

Inukai and Ochi (2009) studied the propeller suction effect on the effective wake
for conventional propeller (CP) and contra-rotating propeller (CRP) on the different
sections upstream of the propeller plane. Figure 1.2 shows a comparison of
experimental results and calculated results using the infinite blade theory. There
is a considerable deviation observed between the numerical results and experimental
measurements.

Pecoraro et al. (2013) studied the inflow characteristics of a single-screw chemical

tanker model due to the large flow separation in the stern region using LDV in the
Large Circulating Water Channel of CNR-INSEAN. A 15% increase in flow rate was
observed on the measurements plane, and non-uniform loading was observed on the
propeller due to the propeller suction effect. The suction effect of the propeller on

6
 Figure 1.2: Axial velocity comparison in front of CP and CRP

incoming flow is shown in Figure 1.3 as a map for the change in axial velocity in the
upstream plane at x/R = −0.514. The effect of propeller suction is more evident in
Figure 1.4 which shows the same data as the change in axial velocity r/R = 0.7. The
non-dimensional axial velocity is plotted against angular position (φ) from the vertical
axis which is pointed in the upward direction. The curve of the non-dimensional
velocity with propeller moves upward due to the propeller suction effect.

Most of the experimental studies have been carried out for the ship model at the
even keel condition with a more or less incoming flow. Because of the cost, time, and
difficulties involved in data acquisition for tests of the self-propulsion with ship motion

7
 (a) Without propeller (b) With propeller

Figure 1.3: The non-dimensional axial velocity at x/R = −0.514

Figure 1.4: The non-dimensional axial velocity at r/R = 0.7.

8
and with incoming waves, limited information is available (Sigmund and el Moctar,
2017). A knowledge gap is still present in this research area. Recent advancements
of computational fluid dynamics help to study the propeller-hull interaction involving
different levels of complexities. The following review of numerical studies explains the
current status of the propeller-hull interaction study.

1.1.2 Numerical Studies

Numerical studies of flow around ships are a more recent development. The use of
numerical methods has grown with increasing computer power and the concurrent
development of more capable software packages. The propeller-hull interaction
was studied by Stern et al. (1988) using the RANS method. In the numerical
simulation, the propeller was modelled using a body-force distribution model, which
was calculated using the vortex-lattice lifting-surface method developed by Kerwin
and Lee (1978). To consider the influence of the effective inflow velocity on the
body-force distribution, the calculation was carried out iteratively. A similar approach
was adopted by Zhang et al. (1992). An actuator disk with distributed body-force
approach was used to model the propeller effect. The body forces were calculated
using lifting-line propeller analysis program and incorporated in viscous flow solver.
Watanabe et al. (1994) used a pressure jump condition to simulate the propeller
effect. The value of uniform pressure jump was calculated from the measured thrust of
self-propulsion experiment. The interaction between the propeller and stern flow was
not considered in the calculation. Propeller-hull interaction was studied by Kawamura
et al. (1997) using a finite volume viscous flow solver WISDAM-V combined with a
potential flow solver developed by Nakatake (1989). The velocity distribution on
the propeller disk was used from the converged RANS simulation(without propeller
modelling) for body-force calculation based on the simplified propeller model. Then

9
the body-forces were introduced into the RANS simulation. The RANS simulation was
performed repeatedly until the flow was converged for a specified thrust controlling the
propeller rate of revolution. For the turbulence modelling, a hybrid turbulence model
of Baldwin-Lomax model with an SGS sub-grid model was used to avoid an unrealistic
velocity profile in the thin boundary layer region by sub-grid scale (SGS) model. An
11% deviation was observed in the thrust deduction factor between numerical and
experimental results. In the numerical calculation, unsteady behaviour of the flow was
neglected as an infinite number of blades were considered in the propeller model. Tocu
and Amoraritei (2008) studied the propeller-hull interaction of an LPG carrier. The
numerical computation was done using Shipflow, and an in-house code for propeller
design. For turbulence modelling, the k − ω SST model was used, which acts as
the k − ω model in the boundary layer region and k −  outside the boundary layer
region. The computations were carried out interactively with the propeller effect,
which was obtained using an actuator disk or the lifting line method. The effective
wake behind the ship hull was compared with experimental results and showed only
a small deviation. The bilge vortices, which are usually generated by large block
coefficient ships, can lead to wake hooks, as shown in Figure 1.5.

Banks et al. (2010) studied the resistance components of KCS using ANSYS CFX
at model scale. The effect of the turbulence model on the numerical results was
investigated using two turbulence models available in the software package, the k −
ω SST and the BSL Reynolds stress turbulence models. The BSL Reynolds stress
model performed better than the k − ω SST model. Dhinesh et al. (2010) studied the
propeller-hull interaction at model scale using a RANS solver with the realizable k − 
turbulence model using STAR CCM+ for a self-propulsion condition. The propeller
mesh motion was accounted for using rigid body motion. A sliding mesh interface
was used to connect the rotating portion and the stationary portion of the mesh. The

10
Figure 1.5: Numerical wake contours behind the ship hull of a LPG carrier at F n =
0.242

numerical results for resistance and thrust were compared with experimental results.
The flow developed by the propeller behind the ship was over-predicted, although
the effect of propeller action was clearly observed on the pressure distribution on the
ship hull and in the free surface elevation. Zhang (2010) used a sliding mesh method
to study the propeller-hull interaction. Rotation of the propeller was modeled using
rigid body motion and a sliding mesh interface. A simple body-force method was also
used to model the propeller effect. The numerical simulations were carried out for the
even keel condition only.

Win et al. (2013) studied the propeller-hull interaction for the Series 60 ship with
CB =0.6, using a combination of viscous flow solver and potential flow solver. The
body-force distribution for the propeller was obtained using the quasi-steady blade
element theory. The numerical results were validated with experimental results behind
the ship model. The self-propulsion performance of a ship in waves was studied by
Winden et al. (2013). In their simulation, the propeller was modelled using the
Boundary Element Method with a RANS solver for the ship simulation. Body-force

11
modelling based on BEM was chosen over a dynamic mesh due to lower computational
time.

The interaction between the hull and propeller was studied by Castro et al. (2011)
at full scale and model scale for the KCS ship using an overset mesh technique. The
overset interpolation was performed using SUGGAR. The numerical computation
was done using CFDShip-Iowa, which is a RANS solver with a blended k − ω/k − 
turbulence model. For the free surface simulation, single-phase level set method was
adopted in the CFDShip-Iowa code. Differences in the torque coefficient KQ and the
nominal wake factor 1 − wn for model scale and full scale simulations were observed.
This was because the boundary layer thickness for the model is relatively thicker
than the full scale, and the Reynolds number in each case is different. The incoming
velocity to the propeller and outgoing velocity were found to be higher, and more
uniform, for full scale than for model scale.

The overset grid method was also used by Shen et al. (2015), who studied the zigzag
maneuvering of the KCS model in OpenFOAM. Actual propeller geometry was used
for self-propulsion simulations. Carrica et al. (2016) studied the zigzag maneuvering
of the KCS model in shallow water using CFDShip-Iowa with an overset method. The
numerical results were compared with the experimental data.

Pereira et al. (2017) used an in-house code ReFRESCO, for the KVLCC2 at
model scale and full scale Reynolds numbers. The effect of turbulence models was
investigated without a wall function using a double body method. The y + values were
set to less than 1.0 for all computations. Results for ship resistance and flow field at

12
the propeller plane showed that the modelling error was strongly dependent on the
turbulence model.

Hull pressure fluctuation due to the rotating propeller behind the hull is a significant
source of hull vibration Hylarides (1978). The non-uniform wake flow into the
propeller introduces unsteady propeller blade loading and pressure, which excites
vibration due to the propeller rotation. The effective wake field due to propeller-hull
interaction is key to the fluctuating pressure excitations. The wake is usually
determined in model basin experiments. In practice, the unsteady propeller forces
for the given wake are determined in numerical simulations Bertram (2012). It is
essential to keep the excitation force within acceptable limits for sophisticated ships.
Lee and Chen (2005) investigated the influence of the propeller-hull interaction effect
on induced hull pressure for a cavitating propeller. The vortex lattice method was used
for propeller simulation with RANS and Euler solvers. The propeller induced pressure
was under-predicted using the Euler solver, which did not include propeller-hull
interaction. Paik et al. (2013) used FLUENT to study the hull pressure fluctuation
for a rotating propeller behind a ship. The importance of accurate wake prediction
was identified for propeller-induced pressure fluctuation. Gaggero et al. (2016) studied
propeller generated noise and hull pressure fluctuation for a single screw research vessel
numerically and experimentally in model scale and full scale. The Boundary Element
Method (BEM) method was used to model the propeller effect. Accurate prediction of
the effective wake was emphasized for propeller radiated noise and pressure fluctuation
on the ship hull. The shortcoming in the prediction of the propeller suction effect was
mentioned. The BEM method tends to over-estimate the cavitation effect, which in
turn influences the pressure prediction. Kim et al. (2018) used potential flow method
to study the propeller-induced hull pressure fluctuation. The Boundary Element

13
Method (BEM) was coupled with the RANS method to predict the induced pressure.
It was concluded that accurate effective wake prediction is crucial near the ship hull.

1.1.3 Problem Statement

From the literature review, it is found that very few attempts have been made to
study the effective wake involving free surface using detailed propeller geometry. It is
also well known that ship motions introduce unsteady flow patterns leading to load
variation on the propeller and change the effective wake, which reduces the propulsive
performance of a ship. The unsteady incoming flow and loading conditions lead to
propeller vibration. The propeller generated vortices make the ship noisy and easily
detectable. From an environmental point of view, propeller noise has a significant
impact, particularly on marine mammals. A better understanding of the flow features
can help to mitigate these unwanted effects.

This thesis uses CFD to get a better understanding of the detailed flow features
in the interaction between the ship hull and the propeller under realistic operating
conditions. This scenario presents considerable challenge for experimentation and
thus CFD presents a viable option for design and evaluation.

1.1.4 Research Aims and Objectives

The aim of this research is to develop an efficient method to consider the propeller
suction effect while considering the hull wake and ship motions using the open source
Computational Fluid Dynamics (CFD) code OpenFOAM. The approach uses the
sliding mesh interface method rather than the overset method for low computational
cost. In addition to a model that uses the detailed propeller geometry, two body-force
models are developed for faster self-propulsion simulation that still considers the

14
motion of the ship as well. The development of the full method requires the following
intermediate objectives to be achieved in support of the full simulation. Five of these
steps involve new software, or modifications of existing software, that were developed
as part of this research.

1. In modelling the water surface, the volume fraction, which is discontinuous

on the water surface, needs to be solved. To resolve the flow features, isotropic
refinement of the grids in the free surface region is very costly due to the division
of each cell in eight for a hexahedral cell. Thus, it becomes impossible to use
isotropic meshing because of the large grids. However, the local flow features on
the free surface can be solved by refining the cells in the normal direction of the
free surface using anisotropic meshing, which greatly reduces the total number
of cells required to solve the free surface accurately. In OpenFOAM, the native
meshing module snappyHexMesh does not support anisotropic meshing. To
overcome this limitation, an anisotropic meshing capability is developed based
on snappyHexMesh for the numerical simulation involving the water surface.
The developed anisotropic meshing method is used in the grid generation process
for simulation of the propeller-hull interaction.

2. In recent versions of OpenFOAM, propeller rotation can be applied for a fixed

condition of the ship. However, the motion of the ship changes the incoming
flow to the propeller. A dynamic motion class is developed, which can be used
for numerical simulation of ship motions combined with rotating propellers for
the detailed study of propeller-hull interaction using actual propeller geometry.

3. The incoming flow to the propeller is highly non-linear due to the curved
geometry of a ship hull and the propeller rotation itself. In the standard
k − ω SST turbulence model, the effect of rotation is not considered. To

15
 consider the effect of rotation and curvature of flow, a curvature correction
is implemented in the standard k − ω SST turbulence model in OpenFOAM to
improve the numerical prediction of propeller-hull interaction.

4. The wave generation and wave damping module in OpenFOAM is adapted for
numerical wave simulation. The effect of ship motion in waves is simulated
using the implemented wave module to investigate the effect of ship motion on
propeller-hull interaction.

5. Two different velocity distribution based body-force methods are developed

using propeller open water data for faster simulation of propeller effects in
OpenFOAM. Although this body-force method is unable to predict flow features
with the same accuracy as the detailed geometry simulation, the propulsion
parameters can be predicted with reasonable accuracy in much shorter time.
Depending on the time and resource availability, the implemented body-force
method can be used for the numerical simulation of ships with or without ship
motions.

6. There is no experimental data available in the open literature to validate a full

CFD simulation that includes propeller action, ship forward motion and wave
effects and wave induced motions. However, data is available for subsets of this
total problem. Thus, the final simulation has been built up incrementally in
stages consisting of:

(a) Ship resistance prediction, mean flow velocity, and pressure distribution
on the hull surface with incoming current. Validation data available.

(b) Propeller open water performance, and upstream and downstream flow.
Validation data available.

16
 (c) Water wave generation. Validated based on analytical solution.

(d) A complete self-propulsion simulation that combines hull modelling,

propeller modelling, wave modelling and wave induced motions. This
simulation provides the final simulation objective of predicting the effective
wave considering the viscous flow with free surface, and pressures around a
working propeller under realistic ship operating conditions. No validation
data are available.

1.1.5 Thesis Outline

The results of this research work are explained in this thesis. This thesis is organized
in seven chapters as follows:

• Chapter 1 presents the current status of the propeller-hull interaction study

using experimental and numerical approaches. The conclusions are summarized
from the related research publications, and a knowledge gap is identified. The
scope of this research is explained at the end of the chapter.

• In Chapter 2, the CFD methodologies, which are used in this thesis,

are presented. Mathematical derivation of governing equations with the
discretization procedures is described. The implemented curvature correction
for the turbulence model, developed dynamic motion class, wave generation,
and absorption modules are included in this chapter.

• The main focus of Chapter 3 is the validation and verification of numerical

results. The detailed validation of the bare hull resistance, free surface wave
contour, mean flow velocity, and pressure distribution on the hull surface of
the KCS, JBC and Fishing vessel with and without ship motions using the
experimental data are presented.

17
• The numerical propeller open water results using detailed propeller geometry
and the body-force methods are validated against the experimental data and
presented in Chapter 4. The objective of this chapter is to find the optimum
settings for propeller modelling for self-propulsion simulation using a viscous
flow solver. Three propeller models KP505, MP687 and a stock propeller, are
used for this part of the study.

• Chapter 5 concentrates on the water wave validation. It explores different

aspects of numerical wave generation, detailed grid dependency, and the
time step size dependency for different wave steepness. The numerical wave
generation for the self-propulsion is also validated.

• The objective of Chapter 6 is to demonstrate full implementation of this

research. At first, the self-propulsion simulations are carried out for the
fixed condition with detailed propeller geometry and body-force method. The
predicted numerical results are validated with the experimental data. Then
the self-propulsion simulations with the ship motions were performed, and
numerical results are discussed. The seakeeping simulation is also validated with
experimental data before the self-propulsion simulation in waves with forward
speed. In the end, the previous simulation modules are integrated into a full
self-propulsion in waves with ship motions.

• The main findings and conclusions are summarized in Chapter 7. The strengths
and limitations of the current development with recommendations for further
research are discussed.

18
Chapter 2

Mathematical Formulation

The following sections detail the mathematical basis for the general flow simulation
and the developed simulation modules incorporated into OpenFOAM as part of this
research.

2.1 Governing Equation for Multiphase Flow

A single-field formulation for two immiscible, incompressible fluids of constant

viscosity is presented here. Neglecting heat and mass transfer only surface tension is
considered in the formulation. In this representation, the fluids are identified by a
step function α, which is 1 for one particular fluid and 0 for the other. Hence, the
density and other material properties can be written in terms of their constant values
for the two fluids as well as for interface.

ρ = αρ1 + (1 − α)ρ2 (2.1)

19
 µ = αµ1 + (1 − α)µ2 (2.2)

where, ρ1 is the density for fluid 1, and ρ2 is the density for fluid 2.

2.1.1 Continuity and Momentum Equations

Assuming the fluid motion is governed by the Navier-Stokes equations, in the

conservation equations for a single-field formulation (Tryggvason et al., 2001; Rusche,
2002), the differences of the material properties and the surface tension force acting
at the interface are considered. The mass continuity and momentum equations are as
follows:

∇•U =0 (2.3)

∂ Z
(ρU ) + ∇ • (ρU U ) = −∇p + ∇ • τ + ρf + σκ0 n0 δ β (x − x0 )dS (2.4)
∂t S(t)

where U is the velocity field, p is the pressure, τ is the stress tensor, f is the
acceleration due to body-forces which is gravitational force, g, δ β is the three
dimensional δ function, σ is the surface tension, and κ is twice the mean curvature
for 3D flows.
Using the continuum surface force model, Brackbill et al. (1992) has shown that the
last term of equation 2.4 is the continuous volumetric force acting within the transition
region from one fluid to another fluid. The last term of equation 2.4 can be written
as

20
 Z
σκ0 n0 δ β (x − x0 )dS ≈ σκ∇α (2.5)
S(t)

2.1.1.1 Pressure Term

The specification of the pressure boundary conditions is simplified if the modified

pressure is used as a dependent variable. It is defined as:

p = p∗ + ρg • x (2.6)

where, p∗ is the modified pressure and ρg • x is the hydrostatic pressure. It helps to

specify the pressure at the boundaries of the space domain. Taking the gradient of
equation 2.6 gives

∇p = ∇p∗ + ∇(ρg • x)

= ∇p∗ + ρg + g • x∇ρ (2.7)

2.1.1.2 Viscous Term

For the Newtonian fluid the stress tensor can be written as.

 
τ = µ ∇U + (∇U )T (2.8)

Using the equation 2.8, the viscous stress term can be simplified for numerical
calculation as follows:

21
   
∇ • τ = ∇ • µ ∇U + (∇U )T
 
= ∇ • (µ∇U ) + ∇ • µ (∇U )T

= ∇ • (µ∇U ) + (∇U ) • ∇µ + µ∇ (∇ • U )

= ∇ • (µ∇U ) + (∇U ) • ∇µ (2.9)

Substituting all the above derived terms in the momentum equation gives

∂
(ρU ) + ∇ • (ρU U ) = −∇p∗ + ∇ • (µ∇U ) + (∇U ) • ∇µ − g • x∇ρ + σk∇α
∂t
(2.10)

2.1.2 Volume of Fluid

The volume averaged flow equations are used to derive the governing equation for the
volume of fluid method (Hirt and Nichols, 1981). The interface is tracked using the
previously defined phase indicator function α (volume fraction) as follows:


1, If the control volume is filled with phase 1







α= 0, If the control volume is filled with phase 2



0 < α < 1, If interface is present




The volume fraction function α is further defined as

H
Vi α(x)dVi
α= (2.11)
|Vi |

22
where, Vi is the volume ith cell
Assuming the velocity is continuous across the interface, the governing equations can
be written as

∂ρ
+ ∇ • (ρu) = 0 (2.12)
∂t

∂
(ρu) + ∇ • (ρuu) = ∇ • T − ρg + Fs (2.13)
∂t

where, Fs is the surface tension force and computed as

Fs = σκ(x)nδΓ + (∇σ)δΓ (2.14)

For constant surface tension force,

Fs = σκ(x)nδΓ (2.15)

∇α 5α
where, δΓ = |∇α| and the unit vector n = and κ (x) = −∇ • |5α| is the curvature
|∇α|
at the interface. The surface tension force can be written as:

Fs = σκ(x)∇α (2.16)

The volume fraction α can be computed from the general transport equation:

23
 ∂α
+ ∇ • (αU ) = 0 (2.17)
∂t

In OpenFOAM, an additional artificial compression term is introduced to achieve the

necessary compression of the free surface, which is only active in the interface region.
The transport equation for the VOF method can be given as follows:

∂α
+ ∇ • (αU ) + ∇ • (α(1 − α)Ur ) = 0 (2.18)
∂t

where Ur is a velocity field suitable to compress the interface, which is calculated

using the following equation:

!!
∇α ∇α
Ur = min Cα n • u, maxf • u (2.19)
|∇α| |∇α|

with Cα the interface compression coefficient, which is typically set as one in numerical
computations.

In the present study, interf aceCompression and High Resolution Interface

Compression (HRIC) (Muzaferija and Peric, 1997) scheme are used to capture the
free surface. As a single velocity U is considered in the whole domain, Ur or the cell
face flux φf is approximated.

2.1.3 Interface Capturing Scheme

The interface capturing scheme discretizes the convective term of volume fraction to
prevent the smearing of the free surface due to the numerical diffusion and maintain

24
 Figure 2.1: Calculation of cell face volume fraction.

the monotonic distribution of the volume fraction. The interface capturing scheme
calculates the cell face value (φf ) of the volume fraction using the neighbouring
cell values (Figure 2.1). The interf aceCompression scheme and HRIC scheme is
explained in this subsection.

2.1.3.1 interfaceCompression Scheme

In interf aceCompression scheme the cell face flux (φf ) is obtained from the following
formulation

φf = λ (φC − φD ) + φD (2.20)

where φC is the value of volume fraction at the current cell, and φD is the value of
volume fraction in the neighbour cell. In interf aceCompression scheme the value of
blending function is given by

 h n o i 
λ = min max 1 − max (1 − 4φC (1 − φC ))2 , (1 − 4φD (1 − φD ))2 , 0 , 1 (2.21)

25
Although, this scheme is does not follow the Total Variation Diminishing (TVD) or
Normalized Variable Diagram (NVD), it provides better boundedness and convergence
(Weller, 2008).

2.1.3.2 HRIC Scheme

The High Resolution Interface Capturing (HRIC) scheme is based on the blending of
the bounded downwind and upwind differencing schemes. This method is independent
of the explicit Courant number. However, it depends on the local face Courant number
(Cf ). The cell face flux (φf ) is calculated using the following formulations

φC − φU
φfC = (2.22)
φA − φU



φfC if φfC < 0 or φfC > 1








φff =  2φfC if 0 < φfC < 0.5 (2.23)





1 if 0.5 ≤ φfC ≤ 1.0





φff if Cf < 0.3







∗    0.7 − C
φff = f (2.24)
φfC + φff − φfC if 0.3 ≤ Cf ≤ 0.7



 0.7 − 0.3


φf if 0.7 ≤ Cf

 C

∆φC
where, θf = cos−1 (|dn|), n = |∆φC |
and d = CD
|CD|
.

∗∗ ∗√ √
φff = φff cosθ + φfC (1 − cosθ) (2.25)

26
 Figure 2.2: Parameters for finite control volume discretization (Rusche, 2002).

∗∗
φf = φff (φD − φU ) + φU (2.26)

2.2 Discretization of Governing Equations

In Finite Volume Method, the RANS equations are solved in a fluid continuum,
which is considered as the composition of discrete finite control volumes (CVs).
The domain discretization can be divided into spatial discretization and temporal
discretization. Many authors have described the discretization of the Finite Volume
Method (Villiers, 2006; Rusche, 2002; Ferziger and Peric, 2012; Patankar, 1980;
Versteeg and Malalasekera, 2007). The similar discretization approaches are used
in this thesis. The spatial discretization represents the individual discrete control
volume as shown in Figure 2.2. The face area vector S is normal to the face and its
magnitude is equal to the area of the face, and the normal vector on the face f is
S
n = |S|
. It points out of the cell of interest P into the neighbouring cell N . The
centres of cell P and N are connected with the vector d. The standard transport

27
equation will be used to describe the discretization process. The general transport
equation in integral form can be express as

∂ Z Z Z Z
ρφdV + ∇ (ρU φ)dV =
• ∇ (Γ∇φ)dV + Sφ (φ)dV
• (2.27)
∂t V V V V

Here, φ is the transport quantity such as velocity, U is the velocity, Γ is the diffusivity
and Sφ is the force term. The volume integrals are reduced to surface integrals by
using Gauss’s theorem. where S is the surface bounding the volume V and dS is an
infinitesimal surface element with outward pointing normal on the surface S.

Z Z
∇φdV = φdS (2.28)
V S

2.2.1 Time Derivative

∂ R
The time derivative ρφdV is discretized by integrating it over the control volume
∂t V
CV . For incompressible flow, assuming the linear variation of φ over the time step
given by

∂ Z φn − φoP
ρφdV ≈ ρ P VP (2.29)
∂t V ∆t

where φn ≡ φ(t + ∆t) is the new value at the current time step and φ0 ≡ φ(t) the old
values from previous time step.

28
2.2.2 Convective Term

The convective term ∇ • (ρU φ) is discretized in the following form using Gauss’s
theorem

Z Z X
∇ (ρU φ)dV =
• (ρU φ) • dS = S • (ρU )f φf (F,S) (2.30)
V S f

Here S • (ρU )f is the mass flux through the face f . The value of φf can be calculated
using Central Differencing, Upwind Differencing and Blended Differencing schemes.

2.2.3 Diffusion Term

The diffusion term ∇ • (Γ∇φ) is discretized using the similar procedure as the
convective term. After applying Gauss’s theorem, the diffusion term can be written
as

Z Z X
∇ • (Γ∇φ)dV = (Γ∇φ) • dS ≈ Γf (S • ∇f φ) (2.31)
V S f

The face gradient of φ can be calculated as

φN − φP
S • (∇φ)f = |Sd | + S∆ • (∇φ)f (2.32)
|d|

Here the vector Sd is the component parallel to the d and S∆ is the remainder which
needs to calculated by the non-orthogonality treatment to preserve the second order
accuracy. The over-relaxed approach proposed by Jasak (1996) is used in this study.

29
2.2.4 Source Term

The source term Sφ φ can be expressed as

Sφ φ = Sc + Sp φ (2.33)

where Sc and Sp can also depend on φ. The volume integral of the source term can
be written as

Z
Sφ (φ)dV = Sc Vp + Sp Vp φp (2.34)
V

2.2.5 Temporal Discretization

The temporal discretization of equation 2.27 can be written as

" #
Z t+∆t ∂ Z Z
ρφdV + ∇ • (ρU φ)dV dt
t ∂t V V
Z t+∆t Z Z 
= ∇ (Γ∇φ)dV +
• Sφ (φ)dV dt (2.35)
t V V

Using the derived discretized term of the transport equation the above equation can
be rewritten as follows

30
  
Z t+∆t φnP − φoP X
ρ VP + S • (ρU )f φf (F,S)  dt
t ∆ f
 
Z t+∆t X
=  Γf (S • ∇f φ) + Sc Vp + Sp Vp φp  dt (2.36)
t f

2.3 Pressure Velocity Coupling

Since there is no explicit equation for pressure calculation, an equation is derived

from the continuity equation for pressure calculation using momentum conservation
equation. There are several pressure-velocity coupling algorithms, such as SIMPLE,
SIMPLEC, PISO, and PIMPLE available in OpenFOAM. The PIMPLE algorithm is
the combination of SIMPLE and PISO algorithm. The PIMPLE algorithm provides
the advantages of both SIMPLE and PISO algorithm.

SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) was originally

developed by Patankar and Spalding (1972). In the SIMPLE algorithm, a initial
guess pressure field is set to solve the momentum equation for velocities. Since the
velocities are calculated using the guessed pressure, the calculated velocities do not
satisfy the continuity equation. Then the pressure and velocities are modified using
pressure correction equations in order for modified velocities to satisfy the continuity
equation. The difference between the SIMPLEC (SIMPLE Consistent) (Doormaal
and Raithby, 1984) and SIMPLE is the how the pressure is corrected. SIMPLEC
shows better performance over the SIMPLE algorithm.

The PISO (Pressure Implicit with Splitting Operator) was developed by Issa (1986)
for pressure velocity calculation procedure. Oliveira and Issa (2001) presented the

31
 Set boundary conditions

Solve momentum equation

Calculate U /

Compute face flux ϕf

No
Solve pressure equation

PIMPLE
Converged Correct ϕf

Yes Correct U /

Solve other discretized equations

Update boundary conditions

Yes PISO No
Converged
Move to the next time step

Figure 2.3: PIMPLE algorithm

improved version of PISO algorithm. In contrast to the SIMPLE algorithm where

pressure correction is employed once for each time step, PISO algorithm involves two
pressure correction steps. A small time step size and Courant number less than 1.0 are
the major drawback of PISO algorithm. In PIMPLE algorithm, the whole solution
procedure is repeated for a number of times until the convergence achieved. This
allows it to use for a larger Courant number. The flow chart of PIMPLE algorithm
is shown in Figure 2.3.

32
2.4 Turbulence Modeling

2.4.1 k − ω SST Turbulence Model

The standard k − ω SST model is a blended turbulence model that uses the k − ω
model for near wall flow and transitions to k− model away from the wall as suggested
by Menter et al. (2003) to overcome the individual weaknesses of the k − ω and k − 
models. The governing equations of the standard model are defined as follows:

" #
∂(ρk) ∂(ρuj k) ∂ ∂k
+ = P + β ∗ ρωk + (µ + σk µt ) (2.37)
∂t ∂xj ∂xj ∂xj
" #
∂(ρω) ∂(ρuj ω) γ ∂ ∂σ ρσω2 ∂k ∂ω
+ = P − βρω 2 + (µ + σω µt ) + 2 (1 − F1 )
∂t ∂xj νt ∂xj ∂xj ω ∂xj ∂xj
(2.38)
The coefficients are blended forms of the two baseline models. The blending function
F1 is defined as
 
F1 = tanh arg14 (2.39)

where
" √ ! #
k 500ν 4ρσω2 k
arg1 = min max , 2 , (2.40)
Cµ ωy y ω CDkω y 2

Here y is the normal distance to the wall and CDkω is the positive portion of the
cross-diffusion term

!
1 ∂k ∂ω
CDkω = max 2ρσω2 , 10−20 (2.41)
ω ∂xi ∂xi

The coefficients σk ,σω , γ and β in above equation are computed with the general form

φ = F1 φ1 + (1 − F1 ) φ2 (2.42)

33
where φ1 and φ2 are the coefficients of the k − ω and k −  turbulence models. The
eddy viscosity is calculated from

a1 k
νt = (2.43)
max (a1 ω, ΩF2 )

where F2 is given by

 
F2 = tanh arg22 (2.44)

and
√ !
2 k 500ν
arg2 = max , (2.45)
Cµ ω y 2 ω

The flow system rotation and streamline curvature influence the turbulence generation
on convex and concave surfaces. The standard k − ω SST turbulence model is unable
to capture this effect. Smirnov and Menter (2009) proposed a rotation and curvature
correction term, fr1 for the turbulence production term of the standard model. The
production term P is defined as:

∂ui
P = fr1 τij (2.46)
∂xj

!
2 ∂uk 2
τij = µt 2Sij − δij − ρkδij (2.47)
3 ∂xk 3

where
fr1 = max [min (frotation , 1.25) , 0.0] (2.48)

2r∗ h −1
i
frotation = (1 + cr1 ) 1 − c r3 tan (c r2 r̂) − cr1 (2.49)
1 + r∗

34
where cr1 , cr2 and cr3 are set as 1.0, 2.0 and 1.0. Einstein notation is adopted in the
equations. The other terms are defined as:

r∗ = S/Ω (2.50)

2Ωik Sjk DSij

 
r̂ = + (εimn Sjn + εjmn Sin )Ωrot
m (2.51)
ΩD3 Dt

!
1 ∂ui ∂uj
Sij = + (2.52)
2 ∂xj ∂xi

" ! #
1 ∂ui ∂uj
Ωij = − + 2εmji Ωrot
m (2.53)
2 ∂xj ∂xi

S 2 = 2Sij Sij (2.54)

Ω2 = 2Ωij Ωij (2.55)

35
  
D2 = max S 2 , 0.09ω 2 (2.56)

where Ωrot is the rotation of the reference frame and the term DSij /Dt represents the
Lagrangian derivative of the strain rate tensor. Sij is the strain rate tensor and Ωij
is the rotation rate tensor. The anisotropic effects were considered in the numerical
simulation by solving the turbulence equation.

2.4.2 LRR Turbulence Model

The effects of streamline curvature, changes in pressure strain and secondary flow
are considered in a Reynolds Stress Turbulence model developed by Launder et al.
(1975). The components of Reynolds stress and the resulting dissipation are calculated
by solving the transport equations. The equations of the LRR model are:

Du0i u0j
" #
∂ ūi ∂ ūi 2
= − u0j u0k + u0i u0k − δij  (2.57)
Dt ∂xk ∂xk 3
 2
 
− C1 u0i u0j − δij k + (φij + φji )2 + (φij + φji )w
k 3
∂u0j u0k 0 0 ∂u0j u0j
" !#
∂ k 0 0 0 0 ∂uk ui 0 0
− Cs ui ul + uj ul + uk ul
∂xk  ∂xl ∂xl ∂xl

The influence of the mean velocity gradient on the pressure strain correlation is
expressed by

36
 !
(C2 + 8) 2 (30C2 − 2) ∂ ūi ∂ u¯j
 
(φij + φji )2 = − Pij − P δij − k + (2.58)
11 3 55 ∂xj ∂xi
(8C2 − 2) 2
 
− Dij − P δij
11 3

where

∂u0 ∂u0
!
Pij ≡ − u0i u0k j + u0j u0k i (2.59)
∂xk ∂xk

∂u0 ∂u0
!
Dij ≡ − u0i u0k k + u0j u0k k (2.60)
∂xj ∂xi

The transport equation for  is given by

u0i u0k ∂ ūi

!
D ∂ k 0 0 ∂ 2
= C u u − C1 − C2 (2.61)
Dt ∂xk  k l ∂xl k ∂xk k

The constants for LRR model are C1 = 1.8, C2 = 0.6, CS = 0.25, C1 = 1.44,
C2 = 1.92, and C = 0.15

2.5 Near-Wall Treatment

To resolve the boundary layer, it is very important to know the near-wall grid
resolutions. A non-dimensional parameter, is used to introduce the near-wall region,
y + is defined as

37
 Figure 2.4: Law of the wall

uτ y
y+ = (2.62)
ν

s
τw
uτ = (2.63)
ρ

where uτ is the friction velocity, τw is the wall shear stress, y is the distance to the
wall, ν is the kinematic viscosity, and ρ is the density. The non-dimensional velocity
is given by

38
 u
u+ = (2.64)
uτ

Experimental measurement shows that near-wall region consists of three layers,

viscous sub-layer, buffer layer and log-law layer as shown in Figure 2.4 (OpenFOAM,
2019)

• The viscous sub-layer exist approximately upto y + = 5. In the viscous sub-layer

viscous stress is dominant. The shear stress can be assumed the wall shear stress
and the velocity profile is linear, which is given by

u+ = y + (2.65)

• The logarithmic layer exist between y + = 5 to y + = 200. In logarithmic layer,

the turbulence stress dominate and control the flow. The velocity profile follows
the logarithmic function given by

1
u+ = × log(Ey + ) (2.66)
κ

where κ is the von Karman constant which is equal to 0.41, and E = 9.8

• In the buffer layer, Reynolds stress and viscous stress are both important, this
layer lies between y + = 5 to y + = 30. The velocity profile does not follow any
specific function.

39
 Figure 2.5: Illustration of the sliding mesh method.

In OpenFOAM, viscous shear stress can be calculated either by solving the viscous
sub-layer directly or using the empirical wall functions. The approach depends on
the non-dimensional parameters y + . A number of wall functions are available in
OpenFOAM for modelling the shear stress. The blended wall functions are also
available in OpenFOAM, which switches from viscous sub-layer to fully developed
turbulent formulation depending on the local y + . In this thesis both approaches
are used. Several turbulence models are also available for high Reynolds number
turbulence modelling and low Reynolds number turbulence modelling.

2.6 Sliding Mesh

The Arbitrary Mesh Interface (AMI) is a sliding mesh technique which is implemented
using the local Galerkin projection as described in Farrell and Maddison (2011)
for unsteady flow simulation in turbo-machinery applications. In sliding mesh
method two or more cell zones are used. At each time step, the mesh of the
rotating/translating zones move relative to each other along the mesh interface. An

40
illustration of a rotating mesh zone is shown in Figure 2.5. At every time step, the
cells of rotating zone updated using a sliding mesh method, and the values lying on
the interface are interpolated to the updated mesh. This enables simulation across
disconnected, non-conformal patches but adjacent mesh domains. In OpenFOAM,
the implementation is fully parallelized, using constrained decomposition. Bensow
(2013) studied propeller performance including the transient wake field and cavitation
behaviour using the sliding mesh interface in OpenFOAM for a 7000 DWT chemical
tanker in model scale. The simulations were performed using a double body model
with a rotating propeller. The vorticity components in planes across the propeller disc
were well predicted using the AMI method without any dissipation or generation of
spurious solution. Chandar and Gopalan (2016) studied the comparative performance
of the AMI, the Generalized Grid Interface (GGI) and the overset method, in terms
of predicted force coefficients, mass conservative properties and parallel scalability for
flow past a spinning cylinder. The global mass conservation for the AMI and overset
were 5.83 × 10−11 and 5.04 × 10−10 respectively. The AMI method has an order of
magnitude lower error than the overset method. Less computational time for the AMI
method is also reported relative to the overset method using parallel computation.

2.7 Dynamic Motion Solver

A new dynamic motion class is implemented in OpenFOAM as part of this research

to study propeller-hull interaction with both the motion of the ship and the rotating
propeller. In the motion class, the propeller rotation is accomplished using the solid
body motion and the cell zone approach. The propeller domain is updated with
respect to the rotation center of the ship. Then the ship motions are calculated using
a six degree of freedom motion solver. The rotating propeller domain is connected to

41
 Figure 2.6: Dynamic update of ship and propeller domains.

the moving ship domain using the Arbitrary Mesh Interface (AMI) to transfer data
from the ship domain to the propeller domain. In the motion solver, the propeller
first rotates at the given rotational speed and the mesh points of the propeller domain
are updated. The resultant forces and moments acting on the whole system are then
calculated by integrating the fluid pressure, p over entire surface area S.

ZZ
Fp = pn̂dS (2.67)
S

ZZ
Mp = rCS × n̂dS (2.68)
S

where rCS is the position vector of each surface relative to the center of rotation of the
body. The velocity and acceleration of the ship are calculated by solving the equation
of motion using the motion solver. Then the mesh is moved to the new position using
the mesh morphing based on spherical linear interpolation method. Figure 2.6 shows

42
 Update propeller domain using Calculate force and moments
given rotation

Solve 6DoF motion and update

Start simulation position

Solve alpha equation

Solve momentum equation

Solve pressure equation

No

Update velocity

Yes No
PIMPLE PISO
Converged? Converged?

Yes

Solve turbulence

Update boundary conditions Move to next time step

Figure 2.7: Calculation procedure using dynamic motion class.

a schematic diagram of the domain update considering the propeller rotation. The
calculation steps are presented in Figure 2.7.

43
2.8 Body-Force Method

As an alternative approach to the previous one, which is relatively computationally

intensive, the effect of the propeller was also modelled using the actuator disk based
on the load distribution proposed by Hough and Ordway (1964). The momentum due
to the propeller action is added as a body-force in the RANS solver. The body-force
varies in the radial direction and the resulting total thrust and torque imposed on the
fluid is the same as the prescribed thrust and torque. The body-force is computed
in such a way that the force distribution approximately follows Goldstein’s optimum
(Goldstein, 1929) distribution. The non-dimensional body-forces in the axial and
tangential directions are calculated using the following equations:

√
fbx = Ax r∗ 1 − r∗ (2.69)

√
r∗ 1 − r∗
fbθ = Aθ ∗ (2.70)
r (1 − rh0 ) + rh0

r0 − rh0 0 r RH
where r∗ = 0
,r = , rh0 = , r is the radius of a propeller section, RH is
1 − rh RP RP
the radius of propeller hub and RP is the radius of propeller tip. Since the body-forces
must result in the prescribed thrust, T , and the torque, Q, coefficients, Ax and Aθ ,
can be computed from:

CT 105
Ax = (2.71)
4 16 (4 + 3rh0 ) (1 − rh0 )

44
 KQ 105
Aθ = (2.72)
4 J 16 (4 + 3rh0 ) (1 − rh0 )
2

T 8KT
CT = 1 = (2.73)
2
ρVa2 πRP2 πJ 2

Va
J= (2.74)
nD

where 4 is the mean chord length projected into the xz plane.

The advance coefficient, J is calculated based on the advance velocity which can
be obtained depending on the the body-force model. Then thrust coefficient KT
and torque coefficient KQ are calculated for the respective J from the open-water
propeller curve. The body-force for each cell in the actuator disk is calculated using
the propeller geometric information and open-water data.

Two different body-force models were developed based on the above formulation. One
is the fixed velocity body-force model (FV-BFM) and the other one is the local velocity
body-force model (LV-BFM). These two models are explained in more detail in the
following subsections.

2.8.1 Fixed Velocity Body-Force Model (FV-BFM)

In the FV-BFM, the average velocity, calculated based on the computed axial velocity
distribution on the propeller plane without involving propeller modelling, is applied

45
on the propeller disk as input. The major drawback of this approach is, the average
velocity does not change with time and the local velocity effect is not considered in
this model. This is important for the ship hull interaction case where the flow differs
across the propeller plane due to the presence of the hull surface as well as the presence
of hull appendages. The self-propulsion simulation was carried out using this method
to provide a comparison with the numerical predictions with the detailed propeller
geometry model.

2.8.2 Local Velocity Body-Force Model (LV-BFM)

As for LV-BFM, the input velocity for each cell is obtained from the velocity sampling
plane upstream of the propeller disk and updated at every time step. Note that
correction is made to the input velocity based on the distance of the velocity sampling
plane from the propeller plane, rotational speed of the propeller, and the propeller
diameter.

2.9 Convergence Criteria

In computational fluid dynamics, a continuous problem domain is replaced with a

discrete domain using a grid. In the continuous domain, each flow variable is defined
at every point in the domain. This discrete system is a large set of coupled, algebraic
equations in the discrete variables. These large systems of linear equations are solved
using an iterative procedure and it is important to monitor the convergence of the
numerical solution. Overall convergence of numerical results is based on the residual
convergence and the convergence error.

46
2.9.1 Residual Convergence

The descretized governing equations are solved as large system of linear equations in
iterative way. The problem can be expressed in matrix form

Ax = b (2.75)

and the residual is defined as:

r = b − Ax (2.76)

The residuals are then normalized using a normalizing factor:

X
n= (|Ax − Ax| + |b − Ax|) (2.77)

where x is the solution vector and x is the average of solution vector.

Then the normalized residuals are expressed as:

1X
r= | b - Ax | (2.78)
n

Since the CFD solution is iterative, it is not possible to achieve exactly zero residuals.
However, the lower the residual value is, the more numerically accurate the solution.
The solution is considered converged when the residuals do not change significantly
over the simulation time. In most of the cases, pressure residuals were obtained at
10−2 or less, and for other parameters the initial residuals were 10−6 .

47
 Figure 2.8: Example of iterative convergence for variable x.

2.9.2 Convergence Error

The convergence error, E%C is monitored for integrated quantities such as resistance,
and propeller thrust is used to judge the convergence of the numerical simulation. The
convergence error E%C for the quantity of interest x is illustrated in Figure 2.8

1 NX data
xi = xj (2.79)
N data j=1

xi+1 − xi
E%C = × 100 (2.80)
xi+1

The numerical simulation can be considered converged when it meets both of the
residual convergence and the convergence error criteria.

48
2.10 Numerical Solvers

In OpenFOAM, a large numbers of solvers are available for steady and unsteady
simulation. In this thesis simpleF oam, pimpleDyM F oam, interF oam, and
interDyM F oam solvers were used for steady and unsteady simulations. A short
description of the mentioned solvers are given below.

simpleF oam is a steady-state solver for incompressible, turbulent flow, single phase
flow. The pressure-velocity coupling is obtained using the SIMPLE (Semi-Implicit
Method for Pressure Linked Equations) algorithm. In the newer releases, it also
includes an option to use the SIMPLEC (Semi-Implicit Method for Pressure Linked
Equations Consistent) algorithm. The multiple reference frames (MRF) models can
be used for the steady-state simulation of turbo-machinery problems.

pimpleDyM F oam- is a transient solver for incompressible, turbulent flow of

Newtonian fluids on a moving mesh for single phase flow. The pressure-velocity
coupling is obtained using PISO and PIMPLE algorithm in OpenFOAM. The time
step size in PISO algorithm is restricted to the Courant number 1.0 for stability
reason, whereas the PIMPLE algorithm allows to use higher time step size.

interF oam and interDyM F oam - are solvers for two incompressible, isothermal
immiscible fluids using a VOF (volume of fluid) phase-fraction based interface
capturing approach, with optional mesh motion and mesh topology changes including
adaptive re-meshing. Both the PISO and PIMPLE algorithm can be used for
pressure-velocity coupling. For the steady-state simulation, Local Time Stepping
(LT S) method is used for time scheme.

49
2.11 Verification and Validation Methodology

2.11.1 Verification

Several methods are available for verification of numerical results. The first one
is the grid-triplet method, which only considers the first term of the Richardson
extrapolation and assumes that the solutions are in the asymptotic range. The
Grid Convergence Index (GCI) derived by Roache (1998) can be used to estimate
the uncertainties due to spatial and temporal discretizations and is widely used and
recommended by ASME and AIAA. The other two common verification methods
in ship hydrodynamics are the Factor of Safety methods (FS) developed by Xing
and Stern (2010) and the Least Squares Root (LSR) method introduced by Eça and
Hoekstra (2014). The LSR method is explained in the following subsections

2.11.1.1 Least Squares Root (LSR) Method

In ship hydrodynamics applications, it is challenging to keep the variation of the grid

exactly systematic (Zou and Larsson, 2014). Since the LSR method considers the
scatter of numerical solutions, it was used in this work to assess the uncertainties
in the numerical solutions due to spatial discretization. The discretization error was
obtained by the Richardson extrapolation as below:

δRE = φi − φ0 = βhPi (2.81)

" N
# 13
1 X
h= (4Vi ) (2.82)
N i=1

50
where φi is the numerical solution of any parameter on a given grid, φ0 is the estimated
exact solution, β is a constant, h is the cell size, P is the observed order of accuracy,
∆Vi is the volume of the ith cell, and N is the numbers of cells. If the numerical
solutions are available on more than three grids, the parameters, φ0 , β and P , can be
obtained in the least square sense with or without weight by minimizing the following
function (Eça and Hoekstra, 2006):

v
u ng
uX 2
S (φ0 , β, P ) = t w i (φi − (φ0 + βhPi )) (2.83)
i=1

where ng is the number of grid sets. For the non-weighted approach, wi = 1, and for
1/hi
the weighted approach, wi = ng .
P
(1/hi )
i=1
The convergence condition is then decided based on the P , following the rules below:

1. Monotonic divergence: P < 0

2. Monotonic convergence: P > 0

3. Oscillator convergence: nch ≥ IN T (ng /3), where nch is the number of triplets
with (φi+1 − φi )(φi − φi−1 ) < 0

4. Otherwise, anomalous behavior

The following three additional error estimator are introduced including the general
RE error estimator to consider the scatter in the numerical solutions (Eça et al., 2010;
Larsson et al., 2014):

02
δRE = φi − φ0 = β02 h2 (2.84)

12
δRE = φi − φ0 = β11 h + β12 h2 (2.85)

51
 max(|φi − φj |)
δ∆M = 
hng
 (2.86)
h1
−1

where, 1 ≤ i, j ≤ ng
Based on the convergence condition, the numerical uncertainty is defined as follows:

1. Monotonic convergence:

(a) 0.95 ≤ P ≤ 2.05 : USN = 1.25δRE + σSD

12 12
(b) P ≤ 0.95 : USN = min (1.25δRE + σSD , 3δRE + σSD )

02 02
(c) P ≥ 2.05 : USN = max (1.25δRE + σSD , 3δRE + σSD )

2. Oscillatory convergence: USN = 3δ∆M

12 12
3. Anomalous behavior: USN = min (3δ∆M , 3δRE + σSD )

02 12
where σSD , σSD , σSD are the standard deviations of the curve fit. The standard
deviation of the least square fit is given by
v
u ng   2
− φ0 + βhPi
uP
u i=1 ng wi φi
u
(2.87)
σSD = t
ng − 3

2.11.2 Validation

The validation of the numerical results was carried out using the simplified version
of the ASME V&V 20 Committee Standard (ASME, 2009) for modelling deficiencies.
The validation error is defined as follows

E =S−D (2.88)

52
 S−D
|E%D| = × 100 (2.89)
D

where S is the numerical solution and D is experimental data. The uncertainty of

validation error E is given by

UE2 = UD2 + USN

2 2
+ USM (2.90)

where UD is the experimental uncertainty, USN is the numerical uncertainty and

USM is the simulation modelling uncertainty. According to Stern et al. (1999), it is
impossible to calculate USM . So the validation uncertainty UV is introduced instead
of UE . The validation error UV is is given by

UV2 = UD2 + USN

2
(2.91)

If |E| ≤ UV , it is assumed that the validation is achieved at UV level, and if |E| ≥ UV ,

the simulation modelling needs to be improved.

2.12 Numerical Wave Tank

Numerical wave generation is an important step in the study of wave structure

interaction, propeller-hull interaction in waves, seakeeping problems etc. For wave
generation in the present work, the wave generation and absorption modules are
studied from IHFOAM (Higuera et al., 2013) and waves2Foam (Jacobsen et al.,

53
2012). In IHFOAM, waves are generated at the inlet boundary, and absorbed right
on the outlet boundary. The wave generation right on the boundary reduces the
computational grid by one wave length at least. The wave absorption formulation in
IHFOAM is based on the shallow water wave absorption. The wave reflection occurs
from the outlet boundary, particularly when the waves are generated for deep water
cases and have high wave steepness (Higuera, 2020). In waves2Foam, relaxation zone
techniques are used for wave generation and wave absorption. The waves are generated
using an analytical expression. The wave velocity, and wave height are imposed
using a blending function in the inlet relaxation zone. Then the generated waves
are propagated through the computational domain. This wave generation technique
increases the computational domain by one wave length which is computationally
expensive for most of the cases. For wave absorption, the wave damping functions are
introduced in the wave absorption zone. Wave reflection is minimal from the outlet
boundary due to the wave damping zone.

In this study, the wave generation approaches are adopted from both modules but for
wave absorption, only the relaxation zone approach is used. For the study of wave
generation effectiveness, the wave generation on the boundary is named Module-01
and wave generation using inlet relaxation is named Module-02. A schematic of a
numerical wave tank is shown in Figure 2.9.

2.12.1 Wave Generation

The wave generation module from IHFOAM will be discussed in this section. Several
waves theories are available for wave generation. Details can be found in IHFOAM
user guide (Higuera et al., 2013). The wave height is calculated based on the wave
theory based on the input wave type in wave dictionary. This boundary condition

54
includes wet, dry and partially wet cells depending on the calculated wave elevation.
Then the wave velocity and volume fraction are calculated and imposed on the inlet
boundary face by face approach. The pressure at the inlet boundary is calculated
using fixed flux pressure boundary condition, which works as a zero gradient boundary
condition.

2.12.2 Wave Absorption

The wave damping method is implemented same as in waves2Foam for wave damping
at the outlet boundary of the computational domain. The damping technique is
applied from the work of Mayer et al. (1998). A damping function:

exp (χ3.5
R )−1
γR (χR ) = 1 −
exp (1) − 1 (2.92)
f orχR ∈ [0; 1]

is applied inside the damping zone in the following way:

φ = γR φcomputed + (1 − γR ) φtarget (2.93)

Figure 2.9: Schematic of numerical wave tank

55
 Figure 2.10: Variation of γR (χR ) for damping zone

where φ is either U or γ. The variation of γR is the same as in Figure. 2.10. The

definition of χR is such that γR is always 1.0 at the interface between the non-damping
part of the computational domain and the damping zone is shown in Figure. 2.9.

2.13 Vortex Identification Criteria

In OpenFOAM, two vortex identification criteria are available based on the velocity
gradient. The velocity gradient D can be decomposed into symmetric tensor and
skew-symmetric components which are strain rate tensor and vorticity tensor.

D = Sij − Ωij (2.94)

1  ∂ui ∂uj
 1  ∂ui ∂uj

where Sij = + and Ωij = − . The characteristic of velocity
2 ∂xj ∂xi
2 ∂xj ∂xi

gradient is given by

56
 λ3 + P λ2 + Qλ + R = 0 (2.95)

where, P , Q and R are the three invariants of the velocity gradient tensor. The
invariants can be decomposed into symmetric and anti-symmetric parts

P = −tr(D) (2.96)

1h i 1
Q= tr(D)2 − tr(D2 ) = ||Ω||2 − ||S||2 (2.97)
2 2

R = −det(D) (2.98)

2.13.1 λ2 -criterion

According to Holmen (2012), λ2 criterion is based on the local minimum pressure

criteria which entails the vortex. The gradient of the acceleration aij can be obtained
from the gradient of the Navier-stokes equations as

1
aij = − pij + νui,jkk (2.99)
ρ

57
The transport equation of vorticity can be obtained using the symmetric and
anti-symmetric parts of the acceleration gradient

DSij 1
− νSij,kk + Ωik Ωkj + Sik Skj = − pij (2.100)
Dt ρ

The negative eigenvalues of S 2 + Ω2 defines the vortex core (Chong et al., 1990). The
−λ2 iso-surfaces present the vortex.

2.13.2 Q-criterion

The fluid region connected with the positive second invariant of the velocity gradient
presents the vortex. This vortex criterion imposes a lower pressure requirement in the
vortex core from ambient pressure. The positive Q value presents the region where
the rotation rate tensor is higher than the strain rate tensor.

2.14 Self-propulsion Modelling

The self-propulsion point for a ship model can be obtained by the load varying method,
the constant loading method, or the mixed loading method (ITTC, 2008). In the load
varying (constant speed) method, the model speed remains constant and the propeller
load is changed by adjusting the revolution rate of the propeller. A set of numerical
simulations are carried out for different propeller revolutions, and the thrust T , the
torque Q, and the tow force FT ow are calculated. The propeller revolution rate at the
model self-propulsion point is then determined by interpolation of the tow force curve
versus the propeller revolution rate. Note that the tow force is referred to the force
on the tow post in physical model tests, which is determined from FT ow = RT (SP ) − T ,
where RT (SP ) is the total resistance of the vessel including the propeller. The tow

58
force at the ship self-propulsion point, also known as the skin friction correction SF C
is calculated using

1
SF C = ρSWm Vm (CTs − CTm ) (2.101)
2

Here, m subscript refers to the model scale and s refers to the full scale data. SW is
wetted surface, V is the velocity, and CT is the total resistance coefficient.

After running the self-propulsion simulation for different rates of propeller revolution,
the tow force is plotted against the rate of revolutions. The ship self-propulsion
point can be found using the calculated SF C. The rate of revolution nm , thrust
Tm and torque Qm are all calculated at model scale. The calculated data at model
scale is converted to full scale using the thrust identity method. The self-propulsion
parameters at full scale are calculated following procedures in Lewis (1988)
The propeller rotational speed at full scale

nm
ns = √ (2.102)
λ

The thrust and Torque is given by

Ts = Tm λ3 (2.103)

Qs = Qm λ4 (2.104)

The thrust deduction factor at self-propulsion is calculated

59
 t = ts = tm = 1.0 + (F D − RTm )/Tm (2.105)

The open water propeller efficiency η0 , and advance coefficient, J0 are calculated using
the propeller open water curve. The wake fraction at model scale is given by

Jo
wm = 1.0 − (2.106)
JSP

The wake fraction from model scale to full scale is calculated using Lewis (1988),

!
CF CFs
ws = wm . s + (t + 0.04) 1.0 − (2.107)
CFm CFm

The relative rotative efficiency is given by

Q0
ηR = (2.108)
QSP

The hull efficiency is given by

1 − ts
ηH = (2.109)
1 − ws

The Quasi Propulsive Efficiency(QPC)

60
 PE
QP C = (2.110)
PD

where, effective power, PE = RTs Vs and delivered power, PD = 2πQs ns .

61
This page intentionally left blank.
Chapter 3

Resistance and Flow Validation

A comprehensive flow simulation that included ship resistance, sinkage, trim, wave
elevation along the hull surface and the wave elevation at different longitudinal cuts
was conducted using the steady and unsteady solvers in OpenFOAM and validated
using experimental data available from the CFD workshops T2005 (2005); G2010
(2010); T2015 (2015). The following sections present the numerical results for the
KRISO Container Ship (KCS), Japan Bulk Carrier (JBC), and a multi-species Fishing
vessel.

The pressure resistance, frictional resistance, and total resistance coefficients were
validated against the experimental data. The resistance coefficients are defined as:

RT
CT = 1 (3.1)
2
ρSW U 2

RF
CF = 1 (3.2)
2
ρSW U 2

63
 RP
CP = 1 (3.3)
2
ρSW U 2

where RT is the total resistance, RF is the frictional resistance, RP is the pressure

resistance, ρ is the water density, SW is the wetted surface and U is the ship velocity,
CT is the total resistance coefficient, CF is the frictional resistance coefficient, and CP
is the pressure resistance coefficient. For transient simulation, the time step size is
calculated based on the Courant number. The Courant number is defined as:

U ∆t
Co = (3.4)
∆x

where U is the ship velocity, ∆t is the time step size and ∆x is the grid size. The
average cell size for the computational grid is calculated

" N
# 13
1 X
h= (∆Vi ) (3.5)
N i=1

where N is the number of cells, and ∆Vi is the volume of the ith cell

3.1 Bare Hull

For verification and validation of the numerical results, the bare hull numerical results
from the KCS test case 2.1 of the Gothenburg 2010 CFD workshop are considered.
The test conditions are listed in Table 3.2. The bare hull resistances of the KCS
hull at model scale were computed using the Local Time Stepping (LT S) method in

64
OpenFOAM. The LT S method reduces the computational time to obtain a steady
simulation significantly. Steady simulation results can also be used as the initial
conditions for an unsteady simulation, which helps to get a steady start-up of the
unsteady simulation and speed up the convergence. Since the JBC hull was free
to heave and pitch in the model tests, the bare hull resistance was computed using
the unsteady solver. At first the steady simulations were carried out using the LT S
method, then the steady solution was used as the initial condition of the unsteady
simulation.

3.1.1 KCS Model

A well-known KRISO Container Ship (KCS) model, which was built at a scale
1:31.5994, was tested at the Maritime and Ocean Engineering Research Institute
(MOERI), Korea, and National Maritime Research Institute (NMRI) in Japan at
different Froude numbers. Resistance, sinkage, trim, and local flow measurement on
the propeller plane were recorded. This data is in public domain (T2005; G2010;
T2015) were used for the validation of the current numerical results.

Figure 3.1: Body plan and profile elevation of KCS.

65
 Figure 3.2: KCS bare hull geometry.

Table 3.1: Principal particulars of the KCS bare hull.

Parameters Full scale Model Scale Unit

Length, LP P 230.00 7.2786 m
Breadth, B 32.20 1.0190 m
Depth,H 19.00 0.6013 m
Draft, T 10.80 0.3418 m
Displacement,∇ 52030.00 1.6489 m3
Wetted surface, S 9530.00 9.5437 m2
Design speed, V 12.35 2.1960 m/s

3.1.1.1 Geometry and Test Conditions

The body plan and profile of the KCS model are shown in Figure 3.1 and Figure 3.2.
The surface file for the geometry of the KCS model was obtained from the T2015
CFD workshop website. The principal particulars of the KCS model are summarized
in Table 3.1. The bare hull test conditions are listed in Table 3.2.

3.1.1.2 Computational Domain and Grid Generation

The computational domain for the numerical simulation of the KCS bare hull was
created following the ITTC guidelines (ITTC, 2011) for CFD applications. A
schematic diagram of the computational domain is shown in Figure 3.3. The inlet
boundary was set 1.0LP P in front of the model, and the outlet boundary was 3.0LP P

66
 Table 3.2: Test case for bare hull resistance

Case 1.1
Model Model 1
Wave Calm
Condition Towing
LP P 7.2786
Fn 0.26
7
Rn × 10 1.26
Rudder without
Propeller without
Attitude No sinkage & trim
Validation variables Resistance
wave elevation
Mean Flow velocity
EFD provider NMRI*

Figure 3.3: Extent of the computational domain for bare hull simulation.

behind the model. The side boundaries were set 1.5LP P from the central plane to avoid
the wall effect on the numerical simulation. The bottom boundary was set 1.5LP P
from the free surface, and the top boundary was set 1.0LP P from the free surface. For

67
 Figure 3.4: Grid distribution around the ship hull.

Figure 3.5: Grid refinement of Kelvin wedge on the free surface.

the numerical simulation, four computational grids were generated using the mesh
generation module navalSnappyHexMesh in OpenFOAM. The navalSnappyHexMesh
mesh generation module was developed for this research, based on snappyHexMesh,
in order to generate the anisotropic cells in the free surface region. The generation of
the anisotropic cells in the free surface region greatly reduces the number of cells in
the computational domain. The grid distribution around the hull is shown in Figure
3.4. Local grid refinements were employed in the bow and stern region of the model
to capture the flow changes. A cylindrical refinement domain is used to capture the

68
 Table 3.3: Boundary conditions for RANS simulations
Parameters Inlet Outlet top bottom/sides hull
U fixedValue outletPhaseMeanVelocity pressureinletOutletVelocity symmetry/slip noSlip
p_rgh fixedFluxPressure zeroGradient totalPressure symmetry/slip fixedFluxPressure
alpha fixedvalue variableHeightFlowRate inletOutlet symmetry/slip zeroGradient
k fixedValue inletOutlet inletOutlet symmetry/slip kLowReWallFunction
ω fixedValue inletOutlet inletOutlet symmetry/slip omegaWallFunction
νt fixedValue zeroGradient zeroGradient symmetry/slip nutUSpaldingWallFunction

Table 3.4: Computational grid for bare hull simulation.

GRID Base cell size (m) h hi /h1 No. of cells Average, y +

GRID-A 1.00 0.1768 1.84 1649858 52.01
GRID-B 0.80 0.1458 1.52 2943851 53.27
GRID-C 0.67 0.1293 1.35 4219043 52.64
GRID-D 0.50 0.0960 1.00 10315239 46.59

Table 3.5: Steady resistance and uncertainty of the KCS bare hull at F n = 0.26.

GRID CP CF CT |E%D| |USN %S|

GRID-A 0.000807 0.002700 0.003507 0.76 8.05
GRID-B 0.000824 0.002730 0.003554 0.57 5.74
GRID-C 0.000788 0.002723 0.003511 0.66 4.80
GRID-D 0.000758 0.002707 0.003464 1.97 3.16
Experimental 0.003534

wake in the propeller plane as well as behind the propeller plane. A Kelvin wedge grid
refinement was used near to the ship hull to capture the generated wave systems more
accurately. The Kelvin wedge refinement for the KCS bare hull simulation is shown in
Figure 3.5. For the bare hull simulation, a symmetry plane boundary condition was
used on the center plane to reduce the computational cost and time. The boundary
conditions used in the numerical simulation are listed in Table 3.3. The summary of
the computational grids for the bare hull resistance validation are presented in Table
3.4.

69
 0.01
Numerical
Experimental

0.008

Total resistance coeff. CT

0.006

0.004

0.002

0
0 2000 4000 6000 8000 10000 12000 14000
Number of Iteration

Figure 3.6: Convergence of total resistance coefficient, CT for GRID-C.

3.1.1.3 Resistance Validation

The resistance coefficients for the bare hull were calculated for the design Froude
number, F n = 0.26 using the four generated sets of grid and steady solver. The
model was restrained from any motion. The numerical results of the force calculation
are shown in Table 3.5. The convergence history of the total resistance coefficient for
GRID-C is shown in Figure 3.6. The simulation was considered converged when the
residual did not change with the number of iterations and the physical quantity (such
as resistance) was stable or followed the same pattern of oscillation. The convergence
of residuals for GRID-C at F n = 0.26 are presented in Figure 3.7. The convergence
error E%C was calculated from averages of 500 iterations for the total resistance
coefficient. The convergence error for GRID-C is shown in Figure 3.8. It is seen that
the convergence error for CT is less than 0.01% after 12000 iterations.

70
 1
p_rgh
alpha.water
0.1 k
omega

0.01
Initial residuals

0.001

0.0001

1e-05

1e-06

1e-07
0 2000 4000 6000 8000 10000 12000 14000
Number of Iteration

Figure 3.7: Residual convergence for GRID-C.

1.0
Avg. over 500 data

0.5
E%C

0.0

−0.5

−1.0
0 2000 4000 6000 8000 10000 12000 14000
Number of Iteration

Figure 3.8: Convergence error for total resistance coefficient for GRID-C at F n = 0.26.

71
 3.7
Numerical

Total resistance coefficient, CT x 103

Fitted Curve
3.65 Experimental

3.6

3.55

3.5

3.45

3.4
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Grid refinement ratio, h i/h 1

Figure 3.9: Numerical uncertainty of CT at F n = 0.26.

GRID-A
0.025 GRID-B
GRID-C
Experimental
0.020
Wave Elevation, ζ/LPP

0.015

0.010

0.005

0.000

-0.005

-0.010

-0.015

0.40 0.20 0.00 -0.20 -0.40

x/LPP

Figure 3.10: Wave elevation along the ship hull at F n = 0.26.

72
 (a) GRID-A

(b) GRID-B

(c) GRID-C

(d) Experimental (NMRI)

Figure 3.11: Free surface wave pattern comparison at F n = 0.26.

73
The validation error for resistance coefficients from the steady simulation are presented
in Table 3.5. The validation error against the experimental data for the total resistance
coefficient CT is less than 2.0%. The uncertainty of the numerical results for the bare
hull resistance is calculated using Least Squares Root (LSR) method. The fitted
curve using LSR method is shown in Figure 3.9. The numerical uncertainties of the
resistance components are also shown in Table 3.5. The numerical result obtained
using GRID-C shows only 5% numerical uncertainty.

3.1.1.4 Free Surface Wave Contour

The numerically predicted wave elevations along the hull surface using three grids are
compared with the experimental data for F n = 0.26 in Figure 3.10. Although there
is a small difference in the prediction for different grids, the overall numerical wave

0.015
GRID-A
GRID-B
GRID-C
Experimental
0.010
Wave Elevation, ζ/LPP

0.005

0.000

-0.005

-0.010

-0.015
0.40 0.20 0.00 -0.20 -0.40 -0.60 -0.80 -1.00 -1.20
x/L PP

Figure 3.12: Wave elevation at y/LP P = 0.1024

74
 0.015
GRID-A
GRID-B
GRID-C
Experimental
0.010

Wave Elevation, ζ/LPP 0.005

0.000

-0.005

-0.010

-0.015
0.40 0.20 0.00 -0.20 -0.40 -0.60 -0.80 -1.00 -1.20
x/LPP

Figure 3.13: Wave elevation at y/LP P = 0.1509

elevation agrees very well with the experimental data. The generated wave pattern
on the free surface is also computed for different grids and is compared with the
experimental results in Figure 3.11. From the figure, it is seen that small cell size, same
as the minimum cell size on the hull surface is necessary to capture the free surface
deformation behind the stern. The wave elevations at the four different transverse
positions y/LP P = 0.10, 0.15, 0.30, and 0.40 are compared with the experimental data
in Figures 3.12 to 3.15. The numerical wave prediction using GRID-C shows good
agreement with the experimental results.

3.1.1.5 Flow Velocity

The axial flow velocity at the propeller plane is compared for the three grids with
the experimental results in Figure 3.16. The wake fraction on the propeller plane is
also compared with the experimental results. From the wake fraction comparison, it

75
 0.015
GRID-A
GRID-B
GRID-C
Experimental
0.010

Wave Elevation, ζ/LPP 0.005

0.000

-0.005

-0.010

-0.015
0.40 0.20 0.00 -0.20 -0.40 -0.60 -0.80 -1.00 -1.20
x/LPP

Figure 3.14: Wave elevation at y/LP P = 0.30

0.015
GRID-A
GRID-B
GRID-C
Experimental
0.010
Wave Elevation, ζ/LPP

0.005

0.000

-0.005

-0.010

-0.015
0.40 0.20 0.00 -0.20 -0.40 -0.60 -0.80 -1.00 -1.20
x/LPP

Figure 3.15: Wave elevation at y/LP P = 0.40

76
 (a) GRID-A, (1 − w) = 0.705 (b) GRID-B, (1 − w) = 0.711

(c) GRID-C, (1 − w) = 0.717 (d) Experimental, (1 − w) = 0.719

Figure 3.16: Axial flow velocity contour on the propeller plane at F n = 0.26.

is seen that GRID-C predicts the wake fraction with more accuracy than the other
grids. The numerical simulation shows two strong opposite vortices just below the
stern bulb, which can not be confirmed because experimental data near the bulb are
not available. The velocity components (u, v, w) are also compared on the propeller
plane at z/LP P = −0.030 in Figures 3.17, 3.18, and 3.19 respectively using all
three computational grids. The velocity components predicted by GRID-C show best
agreement with the experimental data.

77
 1.00

u/U 0.50

0.00

-0.50 GRID-A
GRID-B
GRID-C
Experimental

-0.05 -0.04 -0.03 -0.02 -0.01 0.00

y/LPP

Figure 3.17: Axial flow velocity on the propeller plane at z/LP P = −0.03.

1.00

0.50
v/U

0.00

-0.50 GRID-A
GRID-B
GRID-C
Experimental

-0.05 -0.04 -0.03 -0.02 -0.01 0.00

y/LPP

Figure 3.18: Transverse flow velocity on the propeller plane at z/LP P = −0.03.

78
 1.00

w/U 0.50

0.00

-0.50 GRID-A
GRID-B
GRID-C
Experimental

-0.05 -0.04 -0.03 -0.02 -0.01 0.00

y/LPP

Figure 3.19: Vertical flow velocity on the propeller plane at z/LP P = −0.03.

79
3.1.1.6 Effect of Turbulence Model

The effect of turbulence model on the resistance computation was investigated using
k−ω SST, k−ω SST with curvature correction (k−ω SSTCC), and the Reynolds Stress
Turbulence model LRR for GRID-C. The total resistance coefficient, free surface wave
contour, wave elevation along the hull surface, and nominal wake were computed using
GRID-C with different turbulence models. The numerical results were compared with
the experimental results.

The comparison of total resistance coefficient for different turbulence models is

shown in Table 3.6. The LRR turbulence model shows a larger deviation from the
experimental results than k − ω SSTCC and k − ω SST turbulence models because
the frictional resistance is over-predicted using RSM model. The effect of turbulence
models on the wave elevation along the hull surface is compared in Figure 3.20.
Different turbulence models predict the wave elevation with similar level of accuracy.
The free surface wave contours are also compared with the experimental measurements
in Figure 3.21. No significant difference in wave contours is observed for different
turbulence models.

The axial velocity distribution for the different turbulence models is compared in
Figure 3.22. From the comparison, it seems that the Reynolds Stress Turbulence
model predicts the vortex more accurately, but without experimental data, it is

Table 3.6: Effect of turbulence model on total resistance for GRID-C

GRID CP CF CT |E%D|
k − ω SST 0.000788 0.002723 0.003511 0.66
k − ω SSTCC 0.000822 0.002823 0.003644 3.12
LRR 0.000779 0.003009 0.003788 7.19
Experimental 0.003534

80
 k-ω SST
0.025 k-ω SSTCC
LRR
Experimental
0.020

Wave Elevation, ζ/LPP

0.015

0.010

0.005

0.000

-0.005

-0.010

-0.015

0.40 0.20 0.00 -0.20 -0.40

x/LPP

Figure 3.20: Effect of turbulence models on wave elevation along the hull surface for
GRID-C.

Table 3.7: RMS error of the velocity components for different turbulence models for
GRID-C.

Turbulence Model u/U v/U w/U

k − ω SST 0.0386 0.0108 0.0160
k − ω SSTCC 0.0311 0.0147 0.0184
LRR 0.0357 0.0190 0.0183

difficult to draw specific conclusions because the Reynolds Stress model under predicts
the nominal wake fraction relative to the k − ω SSTCC model. Figures 3.23, 3.24 and
3.25 show the comparison of the velocity components (u, v, w) on the propeller plane
at z/LP P = −0.030. The RM S errors are compared for predicted velocity components
using different turbulence models in Table 3.7. The k − ω SSTCC turbulence model
shows minimum error for axial velocity component, whereas k − ω SST turbulence
model presents minimum error for both transverse and vertical velocity components.

81
 (a) k − ω SST Model.

(b) k − ω SSTCC Model.

(c) LRR Model.

(d) Experimental (NMRI).

Figure 3.21: Effect of turbulence models on free surface wave pattern for GRID-C.

82
 (a) k − ω SST, (1 − w) = 0.720. (b) k − ω SSTCC, (1 − w) = 0.717.

(c) LRR, (1 − w) = 0.711. (d) Experimental, (1 − w) = 0.719.

Figure 3.22: Axial flow velocity contour on the propeller plane using different
turbulence models for GRID-C.

83
 1.00

0.50
u/U

0.00

-0.50 k-ω SST

k-ω SSTCC
LRR
Experimental

-0.05 -0.04 -0.03 -0.02 -0.01 0.00

y/LPP

Figure 3.23: Effect of turbulence model on axial flow velocity on the propeller plane
at z/LP P = −0.03 for GRID-C.

1.00

0.50
v/U

0.00

-0.50 k-ω SST

k-ω SSTCC
LRR
Experimental

-0.05 -0.04 -0.03 -0.02 -0.01 0.00

y/LPP

Figure 3.24: Effect of turbulence model on transverse flow velocity on the propeller
plane at z/LP P = −0.03 for GRID-C.

84
 1.00

0.50
w/U

0.00

-0.50 k-ω SST

k-ω SSTCC
LRR
Experimental

-0.05 -0.04 -0.03 -0.02 -0.01 0.00

y/LPP

Figure 3.25: Effect of turbulence model on vertical flow velocity on the propeller plane
at z/LP P = −0.03 for GRID-C.

85
3.1.2 Japan Bulk Carrier

3.1.2.1 Geometry and Test Conditions

The geometry file for the JBC hull was obtained from the T2015 workshop website
as an IGES file. Figure 3.26 shows the hull surface for the resistance validation of
the JBC. The principal particulars of the JBC hull at full scale and model scale are
shown in Table 3.8.
The CFD simulation conditions were set to be the same as the test conditions
determined by the T2015 organizing committee. The model tests were carried out
for resistance, sinkage and trim for the JBC hull at Froude number, F n = 0.142. All
tests were carried out without rudder. Since the propeller has an overlapping area
with the Energy Saving Device (ESD), the test cases without ESD were considered
for this comparison. The test conditions for resistance, sinkage, trim, and mean flow
are listed in Table 3.9.

3.1.2.2 Computational Domain and Grid Generation

The computational domain for the numerical simulation of the JBC hull was also
created using the ITTC guidelines (ITTC, 2011) for CFD simulation of bare hull

Figure 3.26: JBC bare hull geometry.

86
 Table 3.8: Principal particulars of the JBC hull.

Parameters Full scale Model Scale Unit

Length, LP P 280.00 7.000 m

Length, LW L 285.00 7.125 m
Breadth, B 45.00 1.125 m
Depth,D 25.00 0.625 m
Draft, T 16.50 0.4125 m
Displacement,∇ 178369.90 2.787 m3
Wetted surface w/o ESD, SW 19556.10 12.223 m2
Wetted surface of ESD, SESD 745.20 0.466 m2
Block coefficient, CB 0.8580 0.8580
Midship section coefficient, CM 0.9981 0.9981
LCB (%LP P ).f wd+ 2.5475 2.5475
V CG 13.29 0.3323 m
Design speed, F n 0.142 0.142

Table 3.9: JBC test cases considered for validation of numerical results.

Case 2.1
Loading condition Design full
Wave Calm
Condition Towing
ESD W/O
Rudder without
Propeller without
Attitude sinkage & trim
Validation variables Resistance, Mean velocity
sinkage and trim
wave pattern, turbulence
EFD provider NMRI

87
 Table 3.10: Computational grid for JBC bare hull simulation.

GRID Base cell size (m) h hi /h1 No. of cells Average, y +

GRID-A 1.00 0.1511 2.50 648657 42.53

GRID-B 0.67 0.1094 1.81 1707903 43.13
GRID-C 0.50 0.0811 1.34 4203020 54.77
GRID-D 0.33 0.0604 1.00 10132950 57.71

resistance. A schematic diagram of the computational domain is shown in Figure

3.27. The boundaries were placed in the same way as the KCS model computational
domain. The grids were also generated using navalSnappyHexMesh in OpenFOAM
with local grid refinement. The grid distribution around the JBC hull is shown in
Figure 3.28. Four sets of grids were generated to study the bare hull resistance of the
JBC hull using the unsteady solver. The grid parameters are listed in Table 3.10.

3.1.2.3 Resistance Validation

The bare hull resistance, sinkage and trim of the JBC hull were calculated using the
unsteady multiphase solver interDyM F oam. The model was allowed for sinkage and
trim only. The calculated resistance coefficients, sinkage and trim are compared with
the experimental data in Table 3.11 for Froude number, F n = 0.142. The validation
error for total resistance coefficient, and sinkage for GRID-C is smaller than the other
grids. For trim, the validation error for GRID-C is higher than the other grids. The
overall prediction using GRID-C shows good agreement. The numerical uncertainty
for CT , sinkage and trim is presented in Table 3.12. The uncertainty for CT and
trim obtained using GRID-C is around 10%. Whereas, for trim the uncertainty is
just less than 3% of the numerical results. The uncertainty calculation shows good
convergence for CT and trim. In the following subsection the numerical results for the
first three grids are presented.

88
 Figure 3.27: Computational domain of JBC hull.

Figure 3.28: Grid distribution around the JBC hull.

3.1.2.4 Free Surface Wave Contour

The wave elevation along the hull surface was computed using three grids at F n =
0.142. The numerical results are compared with the experimental values in Figure
3.33. The numerical results show good agreement with the experimental data. The
free surface wave contours for the three grids are also compared with the experimental

89
 0.01
Numerical
Experimental

0.008

Total resistance coeff. CT

0.006

0.004

0.002

0
0 2 4 6 8 10
Time, t (s)

Figure 3.29: Convergence of CT for GRID-C at F n = 0.142.

Residuals
1
alpha.water
k
omega
0.1 p_rgh
pcorr

0.01

0.001
Initial residual

0.0001

1e-05

1e-06

1e-07

1e-08
0 1 2 3 4 5 6 7 8 9 10
Time [s]

Figure 3.30: Time history of residuals convergence for GRID-C.

90
 1.0
Avg. over 250 data

0.5
E%C

0.0

−0.5

−1.0
1 2 3 4 5 6 7 8 9 10
Time, (s)

Figure 3.31: Convergence error of CT for GRID-C at F n = 0.142.

10-3
5.5
Numerical
Fitted Curve
Experimental
Total resistance coefficient, CT

4.5

3.5
0 0.5 1 1.5 2 2.5 3
Grid refinement ratio, h i/h 1

Figure 3.32: Grid uncertainty for CT at F n = 0.142.

91
 GRID-A
0.025 GRID-B
GRID-C
Experimental
0.020

0.015
Wave Height, ζ/LPP

0.010

0.005

0.000

-0.005

-0.010

-0.015

0.40 0.20 0.00 -0.20 -0.40

x/LPP

Figure 3.33: Wave elevation along the ship hull.

Table 3.11: Grid dependency of CT , sinkage and trim at F n = 0.142.

GRID CT |E%D| Sinkage (m) |E%D| Trim (deg) |E%D |

GRID-A 0.00457 6.50 -0.00548 8.98 0.10807 4.78
GRID-B 0.00448 4.45 -0.00563 6.53 0.10742 6.66
GRID-C 0.00421 1.86 -0.00651 7.97 0.11130 7.92
GRID-D 0.00437 1.93 -0.00586 2.65 0.10907 5.75
Exp. 0.00429 -0.00602 0.10313

Table 3.12: Comparison of | USN %S | for CT , sinkage and trim at F n = 0.142.

GRID hi /h1 CT Sinkage Trim

GRID-A 2.50 23.32 3.09 0.14
GRID-B 1.81 13.76 8.47 1.87
GRID-C 1.34 9.34 9.10 2.17
GRID-D 1.00 6.24 10.51 2.32

92
 (a) GRID-A

(b) GRID-B

(c) GRID-C

(d) Experimental (NMRI)

Figure 3.34: Free surface wave pattern comparison at F n = 0.142.

93
 0.015
GRID-A
GRID-B
GRID-C
Experimental
0.010

Wave Elevation, ζ/LPP 0.005

0.000

-0.005

-0.010

-0.015
1.00 0.50 0.00 -0.50 -1.00
x/LPP

Figure 3.35: Wave elevation at y/LP P = 0.1043 at F n = 0.142.

0.015
GRID-A
GRID-B
GRID-C
Experimental
0.010
Wave Elevation, ζ/LPP

0.005

0.000

-0.005

-0.010

-0.015
1.00 0.50 0.00 -0.50 -1.00
x/LPP

Figure 3.36: Wave elevation at y/LP P = 0.1900 at F n = 0.142.

94
 (a) GRID-A (b) GRID-B

(c) GRID-C (d) Experimental (NMRI, 2015)

Figure 3.37: Axial flow velocity contour on x/LP P = 0.9625 for F n = 0.142.

data in Figure 3.34. It is observed that GRID-C predicts the wave contour with more
accuracy than the other two grids.

The numerical wave elevations at two longitudinal wave cut locations y/LP P = 0.1043
and y/LP P = 0.19 are compared in Figure 3.35 and 3.36 for different grid densities.
Although the overall comparison shows good agreement for GRID-C, a small difference
is observed at y/LP P = 0.1043 near x/LP P = 0.5. However a better prediction is
observed at y/LP P = 0.19.

95
 (a) GRID-A (b) GRID-B

(c) GRID-C
(d) Experimental (NMRI, 2015)

Figure 3.38: Axial flow velocity contour on x/LP P = 0.9843 for F n = 0.142.

3.1.2.5 Flow Velocity

The normalized axial velocity at three different longitudinal positions x/LP P = 0.9625,
0.9843 and 1.00 were computed using k − ω SSTCC turbulence model for F n = 0.142.
Comparison is made with the experimental data from NMRI. Figure 3.37 shows the
comparison at x/LP P = 0.9625 for the different grid resolutions. The hull vortex is
captured quite well with GRID-C, but the vortex is stretched-out from the center
plane. A similar phenomenon is observed for the axial velocity contour at x/LP P =

96
 (a) GRID-A (b) GRID-B

(c) GRID-C
(d) Experimental (NMRI, 2015)

Figure 3.39: Axial flow velocity contour on x/LP P = 1.0 for F n = 0.142.

0.9843 and x/LP P = 1.00 using GRID-C. It can be concluded that GRID-C has
sufficient grid density to capture the local features of the flow.

97
3.1.3 Fishing Vessel

A multi-species fishing vessel model with an optimized bulbous bow was also used
for validation of the resistance calculation and numerical self-propulsion using the
body-force method. The detailed grid dependency study can be found in (Ali et al.,
2018) and (Ali et al., 2019a).

3.1.3.1 Geometry and Test Conditions

The hull form of the fishing vessel with bulbous bow is shown in Figure 3.40. The
vessel has a small length-beam ratio of 3.67. The principal particulars of the ship
model and the NRC-IOT stock propeller are given in Table 3.13.

Figure 3.40: Ship model of optimized fishing vessel.

The bare-hull resistance model tests were carried out for the fishing vessel in calm
water and in head seas at the Memorial University towing tank at various speeds,
including the full scale design speed of 12 knots (Friis et al., 2007). The vessel was
free to heave and pitch during the model tests. The test conditions for the fishing
vessel are listed in Table 3.14.

98
 Table 3.13: Principal particulars of the ship model and the propeller.

Parameters Model scale Full scale Units

Length, Lpp 1.833 33.603 m

Breadth, B 0.499 9.148 m
Lpp /B 3.67 3.67
Draft, d 0.224 4.106 m
Wetted Surface Area, Sw 1.110 372.90 m2
Volume of Displacement, ∇ 0.072 442.30 m3
LCG 0.878 16.096 m
V CG 0.259 4.748 m
Propeller Diameter, D 0.1205 2.209 m
P/D 1.27 1.27
AE /A0 0.906 0.906
Number of Blades, Z 4 4

Table 3.14: Test conditions for fishing vessel model.

Case 3.1
Loading condition Design full
Wave Calm
Condition Towing
Rudder without
Propeller without
Attitude sinakge & trim
Validation variables Resistance,

3.1.3.2 Computational Domain and Grid Generation

The computational domain for the numerical simulation of the fishing vessel is shown
in Figure 3.41. The origin o, of the coordinate system, oxyz, was set at the intersection
point of the calm water surface, the midship section, and the ship’s centreplane. The
positive oz points up vertically, and the positive ox points from the stern to the bow.

99
 Figure 3.41: Computational domain for Fishing vessel.

In the computational domain, the inlet and outlet boundaries are 1.5 LP P and 2.5 LP P
from the midship section, respectively. The side boundaries are 1.0 LP P from the
centreplane. The bottom is 1.0 LP P below the free surface, and the top boundary is
placed 0.5 LP P above the free surface. The symmetry boundary condition is imposed
on the centreplane to reduce the number of grids and computational time.

Grids were generated using blockMesh and navalSnappyHexMesh in OpenFOAM.

Using the STL surface of the vessel model, isotropic cells were generated near the
ship hull. To reduce the number of cells in the free surface region, anisotropic cells
were generated. Fine cells were also generated in the propeller domain so that a
sufficient number of cells can be distributed inside the actuator disk region. Different
grid refinement ratios were used in the grid generation process because of difficulties

Table 3.15: Computational grids for simulations of fishing vessel.

LP P /∆x No. of cells (million) hi /h1

GRID-A 470 0.70 2.36
GRID-B 700 1.48 1.84
GRID-C 940 4.16 1.31
GRID-D 1400 9.24 1.00

100
 Figure 3.42: Grid distribution around the bulbous bow.

Table 3.16: Comparison of total resistance coefficient, CT .

Grid Cp CF CT |E%D| |USN %S|

GRID-A 0.006431 0.002885 0.010424 8.82 31.30
GRID-B 0.006166 0.002779 0.010011 4.51 19.64
GRID-C 0.005874 0.002732 0.009655 0.79 10.46
GRID-D 0.005914 0.002626 0.009547 0.33 6.28
Experimental 0.009578

in maintaining a constant grid refinement ratio using navalSnappyHexMesh. The

generated grids for bare hull resistance simulation are summarized in Table 3.15.

3.1.3.3 Resistance Validation

The bare hull resistance of the fishing vessel was calculated using the steady solver in
OpenFOAM. The model was fixed at the experimental sinkage and trim condition.
Table 3.16 presents |E%D| of total resistance coefficient for different grids. It can
be seen that |E%D| gradually decreases with the increased number of cells in the
computational domain. The |E%D| values for the predicted total resistance coefficient
for both models are less than 1.0 % using GRID-C and GRID-D.

101
 11.5
Numerical

Total resistance coefficient, CT x 103

Fitted Curve
Experimental
11

10.5

10

9.5

9
0 0.5 1 1.5 2 2.5
Grid refinement ratio, h i/h 1

Figure 3.43: Total resistance coefficient for F n = 0.34.

The convergence of the numerical solution to the grid was assessed using the LSR
method for the total resistance coefficients. As shown in Table 3.16, the uncertainty
for each grid is given in percentage term |USN %S|, where USN is the uncertainty of
the solution, S. GRID-C and GRID-D have uncertainties around 10.0% and 5.0%.
The fitted curve using the weighted LSR method for the total resistance coefficient
is plotted in Figure 3.43 with the experimental data. A satisfactory grid convergence
was achieved with an order of accuracy P = 3.09 for the model.

GRID-C was then used to predict the total resistance coefficients of the model for
another five Froude numbers (F n). Figure 3.44 shows comparison of the numerical
results and experimental data. The numerical results show a good agreement for all
Froude numbers.

102
 Table 3.17: Comparison of |E%D| for CT , sinkage and trim at F n = 0.34

GRID CT |E%D| Sinkage (m) |E%D| Trim (deg) |E%D|

GRID-A 0.0106 12.80 0.0143 6.26 0.5103 6.99
GRID-B 0.0104 10.70 0.0142 5.69 0.5194 5.34
GRID-C 0.0096 1.79 0.0137 1.82 0.5265 4.04
Experimental 0.0094 - 0.0135 - 0.5487 -

0.025
Experimental
Numerical
Total resistance coefficient, CT

0.02

0.015

0.01

0.005
0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45
Froude Number, Fn

Figure 3.44: Total resistance coefficient of Model-H using GRID-C.

The bare-hull resistance, sinkage, and trim were computed using the unsteady solver,
interDyM F oam at the design Froude number, F n = 0.34. The predicted total
resistance coefficient, sinkage, and trim with the three grids are compared with the
experimental data in Table 3.17. The relative errors of predicted total resistance,
sinkage and trim are less than 5.0 % when GRID-C is used.

103
 (a) GRID-A

(b) GRID-B

(c) GRID-C

(d) GRID-D

Figure 3.45: Free surface wave contours for Model-H at F n = 0.34

104
 (a) GRID-A

(b) GRID-D

Figure 3.46: Bow waves of Model-H at F n = 0.34

3.1.3.4 Free Surface Wave Contour

The non-dimensional wave patterns generated by the fishing vessel model are
compared in Figure 3.45 for the four computational grids. The simulated wave
contours show differences in the bow region. Wave elevations in the bow region by

105
 0.05
GRID-A
GRID-B
GRID-C
0.04 GRID-D

0.03

Wave Elevation, ζ/LPP 0.02

0.01

-0.01

-0.02

-0.03
-0.4 -0.2 0 0.2 0.4
x/LPP

Figure 3.47: Comparison of wave elevation along the hull of Model-H at F n = 0.34

GRID-A are less than those by GRID-D. To investigate the generated wave in the
bow region, snapshots of the bow wave using two grids are compared in Figure 3.46.
From the comparison, it is seen that GRID-D was able to capture the overturning of
the bow wave better than GRID-A. It is evident that fine cells are required to capture
the highly non-linear waves in the bow region.

Figure 3.47 presents the predicted wave elevations along the hull using different grid
resolutions. Some differences in the bow shoulder waves (x/LP P = 0.15 ∼ 0.4) can
be observed.

3.1.3.5 Flow Velocity

The axial flow velocity for the bare hull was computed on the propeller plane using
the unsteady solver and is shown in Figure 3.48 for F n = 0.34. A strong boundary
layer effect due to the skeg is visible.

106
Figure 3.48: Axial velocity contour on the propeller plane at F n = 0.34 for GRID-C.

107
3.2 Hull with Appendages

3.2.1 KCS Model

The resistance, sinkage, and trim of the KCS model with rudder were numerically
computed using the unsteady multiphase solver in OpenFOAM. The geometry of the
KCS model and rudder is shown in Figure 3.49.

Figure 3.49: KCS model with rudder.

3.2.1.1 Geometry and Test Conditions

For the test case 2.1 from Tokyo 2015 CFD workshop, the KCS model was equipped
with a semi-balanced horn rudder. The geometry of the rudder is shown in Figure
3.50. Table 3.18 gives the particulars of the KCS rudder at full scale and model scale.
The test conditions for the appended KCS model are shown in Table 3.19

3.2.1.2 Computational Domain and Grid Generation

A right-handed Cartesian coordinate system was used for numerical simulation set-up.
The origin O was set at the intersection of midship section, transverse plane, and
calm water surface. The x-axis is positive pointing to the bow from the stern, the
positive y is in the starboard direction, and the positive z is in the upward direction.
The same computational domain from the bare hull resistance computation was used
for transient simulation. The same grid resolution was also used for the dynamic

108
 Table 3.18: Principal particulars of rudder for KCS hull.

Parameters Full scale Model Scale Unit

Profile type NACA 0018 NACA 0018

Rudder area AR 54.45 1.7232 m2
Percentage of lateral area AR /LT 2.21 2.21 %
Movable rudder area AM 45.3 1.4336 m2
Percentage of lateral area AM /LT 1.83 1.83 %
Balancing rudder area AB 13.62 0.4310 m2
Rudder height b 9.9 0.3133 m
Mean chord length cm 5.5 0.1741 m
Mean thickness tm 0.99 0.0313 m
Aspect ratio, geometric 1.8 1.8
Thickness ratio tm /cm 0.18 0.18
Degree of twist 0 0 deg

Table 3.19: Test conditions for appended KCS model.

Case 1.2
Model Model 1
Wave Calm
Condition Towing
LP P 7.2786
Fn 0.108,0.152, 0.195,
0.227, 0.26, 0.282
Rn × 106 5.23, 7.33, 9.42
11.0, 12.6, 13.6
Rudder with
Propeller without
Attitude sinkage & trim
Validation variables Resistance
sinkage and trim
EFD provider NMRI*

109
 Figure 3.50: Semi-balanced horn rudder for KCS.

Figure 3.51: Grid distribution for sinkage and trim computation.

simulation case. The computational grids for resistance, sinkage, and trim calculation
are summarized in Table 3.20. The grid distribution along the ship hull with local
refinement is shown in Figure 3.51.

110
 Table 3.20: Computational grid for KCS with rudder.

GRID Base cell size (m) h hi /h1 No. of cells Average, y +

GRID-A 1.00 0.1716 1.79 1804097 52.01
GRID-B 0.80 0.1409 1.47 3257069 53.27
GRID-C 0.67 0.1269 1.32 4461541 52.64
GRID-D 0.50 0.0960 1.00 10344597 41.58

Table 3.21: Grid dependency of CT , sinkage and trim at F n = 0.26

GRID CT ×103 |E%D| Sinkage (m) |E%D| Trim (deg) |E%D|

GRID-A 3.577 3.6014 -0.01440 3.2999 0.1616 4.3787
GRID-B 3.568 3.8520 -0.01459 4.6628 0.1605 5.0296
GRID-C 3.594 3.1449 -0.01441 3.3716 0.1650 2.3669
GRID-D 3.600 2.9957 -0.01501 7.6758 0.1638 3.0769
Experimental 3.711 -0.01394 0.1690

Table 3.22: Comparison of | USN %S | for CT , sinkage and trim at F n = 0.26.

GRID hi /h1 CT Sinkage Trim

GRID-A 1.79 0.09 0.90 0.55
GRID-B 1.47 0.45 1.91 0.85
GRID-C 1.32 0.52 2.45 1.18
GRID-D 1.00 0.59 2.77 1.59

3.2.1.3 Resistance, Sinkage, and Trim Validation

The total resistance coefficient (CT ), sinkage, and trim for the three grids were
computed for F n = 0.26. The model was allowed only heave and pitch. The numerical
results are validated against the experimental data in Table 3.21. The numerical
uncertainty for resistance, sinkage and trim are calculated using LSR method and
summarized in Table 3.22 for the four computational grids. The numerical uncertainty
for total resistance coefficient of GRID-C and GRID-D is less than 1%. For sinkage
and trim, the numerical uncertainty of GRID-C and GRID-D are also less than 3%.

111
Table 3.23: Results of unsteady simulation of KCS hull with rudder using GRID-C.

Experimental Numerical
F n CT × 103 Sinkage, σ/Lpp Trim τ CT × 103 Sinkage, σ/Lpp Trim τ
0.108 3.796 -0.00012 -0.017 3.626 -0.00027 -0.021
0.152 3.641 -0.00038 -0.053 3.461 -0.00056 -0.047
0.195 3.475 -0.00082 -0.097 3.350 -0.00105 -0.089
0.227 3.467 -0.00130 -0.127 3.390 -0.00141 -0.123
0.260 3.711 -0.00192 -0.169 3.594 -0.00198 -0.165
0.282 4.501 -0.00234 -0.159 4.410 -0.00255 -0.143

5
Numerical
Total resistance coefficient, CT x 103

Experimental

4.5

3.5

3
0.1 0.15 0.2 0.25 0.3
Froude number, Fn

Figure 3.52: Comparison of total resistance coefficient using GRID-C.

GRID-C was used for the computation of total resistance coefficient, sinkage and
trim for five other Froude numbers. The results are summarized in Table 3.23.
For the low-speed case, it is found that the error in numerical prediction is quite
large. The total resistance coefficient, CT is compared with the experimental data in
Figure 3.52. A large deviation from the experimental data is observed for low speeds.
The discrepancy is reduced with increasing Froude numbers. A similar trend is also
observed in sinkage results in Figure 3.53. However, numerical results for trim show

112
 0
Numerical
Experimental
-0.0005

Sinkage, σ/Lpp -0.001

-0.0015

-0.002

-0.0025

-0.003
0.1 0.15 0.2 0.25 0.3
Froude number, Fn

Figure 3.53: Comparison of sinkage of KCS model using GRID-C.

0
Numerical
Experimental

-0.05
Trim, τ

-0.1

-0.15

-0.2
0.1 0.15 0.2 0.25 0.3
Froude number, Fn

Figure 3.54: Comparison of trim of KCS model GRID-C.

113
good agreement with the experimental data at low speed in Figure 3.54. At the low
speed, the Reynolds number is lower than the higher speed. It increases the laminar
boundary layer and turbulent boundary layer thickness. To keep the same y + for
low speed simulation, the near-wall thickness should be adjusted for similar level of
accuracy as for high speed. The grids were designed for Froude number F n = 0.26
and it was necessary to modify the grids according to the Froude number.

114
3.3 Scale Effect on Resistance, Sinkage and Trim

The scale effect on resistance coefficient, sinkage, trim, free surface wave contour,
axial flow velocity was studied for the design speed (F n = 0.26) using the similar
set-up for the computational domain as for bare hull resistance. Four different scales
of the KCS model were used for this study.

3.3.1 Computational Domain and Grid Generation

The computational domain for grid generation to study the scale effect is shown
in Figure 3.55. The domain shape and size with respect to the ship length was
same for four scale factors. The grid density along the ship length, in the free
surface region, and the in Kelvin wedge refine was similar. To reduce the number
of computational grids for the different scale, the near-wall cell height was increased
and the number of prism layer also was changed. It is quite impractical to generate
grids with the same y + for different model scale. In this case, the total number of

Figure 3.55: Computational domain for the study of scale effect.

115
 (a) Scale factor = 31.599

(b) Scale factor = 16.

(c) Scale factor = 8.

(d) Full scale.

Figure 3.56: Grid distribution at bow and stern of the model.

116
 Table 3.24: Computational grid for scale effect study.

Scale factor No. of cells Average, y +

1 7437954 715.42
8 6430654 293.48
16 6430654 109.71
31.599 4461541 44.21

grids was very high. In multiphase simulations, high aspect cells on the hull surface
leads to the numerical ventilation. This unwanted phenomenon lead to the wrong
numerical results and eventually the simulation crash. Special attention was paid
during the mesh generation step to maintain the cell aspect ratio as low as possible.
The comparison of the distribution on the bow for four scale factor is shown in Figure
3.56. The grid information is summarized in Table 3.24.

3.3.2 Resistance, Sinkage and Trim

The frictional resistance and pressure resistance were computed using the multiphase
dynamic motion solver interDyM F oam. The ship was appended with semi-balanced
rudder. The numerical frictional resistance coefficients are compared with ITTC-57
for different scale factors in Figure 3.57. The difference between the numerical results
and ITTC data is higher for the full scale ship than at model scale. The higher
Reynolds number reduces the relative boundary layer thickness for the full scale case
which increases the frictional resistance. The pressure resistance coefficients are also
compared for the different scale factors in Figure 3.58. A large deviation is observed
for the full scale simulation.

The scale effect on sinkage and trim was investigated, and the comparisons
are presented in Figures 3.59 and 3.60 for sinkage and trim respectively. The

117
 5.00
ITTC
Numerical
4.00

3.00
CF x 103

2.00

1.00

0.00
0 5 10 15 20 25 30 35
Scale factor, λs

Figure 3.57: Effect of scale factor on frictional resistance.

2.00
Experimental
Numerical

1.50
CP x 103

1.00

0.50

0.00
0 5 10 15 20 25 30 35
Scale factor, λs

Figure 3.58: Effect of scale factor on pressure resistance.

118
 0.000
Numerical
Experimental
-0.001

Sinkage, σ/Lpp -0.002

-0.003

-0.004

-0.005
0 5 10 15 20 25 30 35
Scale factor, λs

Figure 3.59: Effect of scale factor on sinkage.

non-dimensional sinkage shows similar values for all scale factors except 16. The
comparison for trim shows that trim is higher at model scale.

3.3.3 Free Surface Wave Contour

The effect of scale factor on the free surface was examined for the cases given in
Table 3.24. Figure 3.61 shows the comparison for four scale factors. The difference in
the bow region is not significant. However, the difference at the stern is significant.
Higher wave height is observed for full scale simulation than model scale behind the
stern of the vessel.

3.3.4 Flow Velocity

The axial velocity distribution in the propeller plane is compared in Figure 3.62 for
different scale factors. The vortex around the stern bulb is stronger in model scale than

119
 0.00
Numerical
Experimental
-0.10

-0.20
Trim, τ

-0.30

-0.40

-0.50
0 5 10 15 20 25 30 35
Scale factor, λs

Figure 3.60: Effect of scale factor on trim.

in full scale. The boundary layer is also smaller in the full scale than at model scale
which causes higher incoming velocity to the propeller. This difference in incoming
velocity leads to a difference in propulsive performance from model scale to full scale.
The challenges involve in the full scale numerical simulation is the grid generation.
The wall function approach predicts the pressure force on the hull surface when y +
is less than 300. The insufficient number of grid points near to the wall lead to the
large deviation in pressure force prediction. From the comparison of force coefficients,
it is recommended to increase the grid resolution in the log law region. Items for
summary covers the first stage of the simulation development and shows that the ship
resistance, flow and pressures can be reasonably predicted.

120
 (a) Scale Factor 31.599

(b) Scale Factor 16

(c) Scale Factor 8

(d) Scale Factor 1

Figure 3.61: Effect of scale factor of free surface wave contour.

121
 (a) Scale factor, λs = 31.599 (b) Scale factor, λs = 16

(c) Scale factor, λs = 8 (d) Scale factor, λs = 1.0

Figure 3.62: Effect of scale factor on axial velocity distribution on propeller plane.

122
Chapter 4

Propeller Open Water Validation

This chapter presents the validation of propeller open water simulations of three
different propellers KP505, MP687, and NRC IOT stock propeller using the detailed
propeller geometry or/and body-force method. Detailed propeller geometry was
used for the KP505, and MP687 propellers and the body-force method, described in
section 2.8, was used for the KP505 and NRC IOT stock propeller. The open water
hydrodynamic characteristics of the propellers are compared with the experimental
data using thrust coefficient KT , torque coefficient KQ , and efficiency ηo as defined
by the following equations:

T
KT = (4.1)
ρn2 D4

Q
KQ = (4.2)
ρn2 D5

123
 J KT
ηo = (4.3)
2π KQ

where T is the thrust, Q is the torque, ρ is the fluid density, n is the rate of revolution,
D is the propeller diameter, and J is the advance coefficient which is defined as:

Va
J= (4.4)
nD

where Va is the advance velocity. Ali et al. (2017b) carried out a detailed investigation
of tip vortex prediction of the DTMB5168 propeller can be found in Appendix-A.
Based on the experience of tip vortex study, propeller open water data validation
were performed for the KP505, MP687 and NRC-IOT propellers.

4.1 KP505 Test Case

4.1.1 Geometry and Test Conditions

4.1.2 Actual Propeller Geometry

The self-propulsion test case 2.7 from the T2015 workshop was performed using
a five bladed right-hand fixed pitch propeller designated KP505. The open water
experiments were carried out at NMRI, Japan using a propeller provided by KRISO.
The surface file for the KP505 propeller was obtained from the T2015 workshop
website. The geometry of the KP505 propeller is shown in Figure 4.1. The main
particulars of the KP505 propeller are listed in Table 4.1.

124
 Figure 4.1: Geometry of KP505 propeller.

Table 4.1: Principal particulars of KP505 propeller.

Item Full scale Model scale

Type FP FP
No. of blades 5 5
D (m) 7.9 0.25
P/D(0.7R) 0.997 0.997
Ae/A0 0.8 0.8
Rotation Right hand Right hand
Hub ratio 0.18 0.18

4.1.2.1 Computational Domain and Grid Generation

Cylindrical computational grids were generated for the propeller open water
simulation of KP505. The inlet boundary was placed 2.0D from the propeller center,
the outlet boundary was set at 4.0D behind the propeller plane, and the outer
boundary was set 2.0D from the propeller axis of rotation as shown in Figure 4.2.
The computational domain consists of two regions, one is stator which is stationary
about the global coordinate system and the other one is the rotor which rotates

125
 Figure 4.2: Computational domain of KP505 propeller open water simulation.

Figure 4.3: Grid distribution on KP505 propeller blade.

about the propeller axis of rotation in a transient simulation. The rotating region
covers the propeller diameter, which is connected with the static region using the
sliding mesh interface. Three computational grids were generated using NUMECA
Hexpress, cfMesh (Juretic, 2017) and navalSnappyHexMesh in OpenFOAM. The grid
distribution on the propeller blade is shown in Figure 4.3. The computational grids
for propeller open water simulations are summarized in Table 4.2.

126
 Table 4.2: Computational grid for KP505 propeller open water simulation

Grid Stator Rotor Average y +

GRID-A 2182803 1347576 45.51
GRID-B 2182803 2905635 22.72
GRID-C 2182803 4250712 15.23
GRID-D 2182803 5868703 4.23

4.1.2.2 Numerical Results

The propeller open water simulations of the KP505 propeller were performed using
the steady and unsteady solvers for a wide range of advance coefficients. For steady
simulations, the simpleF oam solver was used, and for transient simulations the
pimpleDyM F oam solver was used.

4.1.2.3 Steady Simulation

The steady simulations were performed using a moving reference frame (MRF)
method. In the MRF method, a Coriolis force was applied to the cells inside the
rotating region to mimic the rotational effect of the propeller. Physically the rotational
part of the domain stays fixed. In OpenFOAM, SIM P LE algorithm was used for
pressure-velocity coupling. The grid uncertainty is studied using the above mentioned
grids using LSR method for J = 0.7. Table 4.3 presents the validation error with the
grid uncertainty for each grid for steady solution. The order of accuracy, P for KT
is 1.09 and for 10KQ is 2.01. The numerical uncertainty for KT varies from 25% to
30% and for KQ the numerical uncertainty varies from 0.5% to 1.0%. The validation
error for both KT and KQ is less than 10%. The convergence of KT , 10KQ and ηo
for J = 0.7 is shown in Figure 4.4. The parameters converged after 50 iterations,
whereas the residuals did not converge until 1200 iterations as shown in Figure 4.5.
Simulations were stopped after the residuals converged.

127
 4.0
KT
10KQ
ηo
3.0

2.0
KT, 10, KQ , ηo/100

1.0

0.0

-1.0

-2.0
0 200 400 600 800 1000 1200 1400 1600
Number of Iteration

Figure 4.4: Convergence of KT , 10KQ for J = 0.7 using steady solver for GRID-B.

1
Ux
Uy
0.1 Uz
p
k
omega
0.01
Initial residuals

0.001

0.0001

1e-05

1e-06

1e-07
0 200 400 600 800 1000 1200 1400 1600
Iterations

Figure 4.5: Residual convergence using steady solver for GRID-B at J = 0.7.

128
Table 4.3: Grid dependency of open water performance at J = 0.7 using steady solver.

GRID hi /h1 KT |E%D| |USN %S| 10KQ |E%D| |USN %S|

P = 1.09 P =2.01
GRID-A 1.32 0.1915 3.50 31.12 0.3204 3.03 0.32
GRID-B 1.17 0.1832 1.00 28.61 0.3286 5.66 0.67
GRID-C 1.08 0.1814 1.93 26.62 0.3206 3.07 0.87
GRID-D 1.00 0.1793 3.07 24.90 0.3243 4.27 1.02

1.4 100
KT - Experimental
KT - Numerical
10KQ - Experimental
1.2 10KQ - Numerical
ηo - Experimental 80
ηo - Numerical
1.0
KT, 10KQ

0.8 60

ηo
0.6 40
0.4
20
0.2

0.0 0
0 0.2 0.4 0.6 0.8 1
Advance Coefficient, J

Figure 4.6: KP505 propeller open water simulation using steady solver for GRID-B.

The results obtained using the MRF method using GRID-B are compared with the
experimental results from Fujisawa et al. (2000) in Figure 4.6. The predicted KT and
10KQ showed good agreement with the experimental measurements. However, the
open water efficiency was under predicted for J = 0.7 and higher values.

129
Table 4.4: Grid dependency of open water performance at J = 0.7 using unsteady
solver.

GRID hi /h1 KT |E%D| |USN %S| 10KQ |E%D| |USN %S|

P = 2.35 P =2.86
GRID-A 1.32 0.1868 0.99 4.13 0.3118 0.24 2.54
GRID-B 1.17 0.1789 3.32 4.29 0.3187 2.47 2.37
GRID-C 1.08 0.1808 2.29 4.22 0.3181 2.30 2.37
GRID-D 1.00 0.1876 1.41 4.06 0.3103 0.22 2.43

4.1.2.4 Unsteady Simulation

Unsteady numerical simulations of propeller open water tests were also carried out for
the same advance coefficients as the steady simulations. In the transient simulations,
the rotating part of the computational grid rotates at a given rotational speed. The
data from the static part are transferred to the rotating part using the sliding mesh
method. The grid uncertainty was also calculated for J = 0.7 for four grids. The

1.0
KT
10KQ
ηo

0.8
KT, 10KQ , ηo/100

0.6

0.4

0.2

0.0
0 0.05 0.1 0.15 0.2 0.25
Time, (s)

Figure 4.7: Convergence of KT , 10KQ for J = 0.7 using transient solver for GRID-B.

130
 1
Ux
Uy
Uz
0.1 p
k
omega

Initial residuals 0.01

0.001

0.0001

1e-05

1e-06
0 0.05 0.1 0.15 0.2 0.25
Time, (s)

Figure 4.8: Residual convergence using transient solver for GRID-B at J = 0.7.

validation error and the numerical uncertainty are presented in Table 4.4. It is
observed that the numerical uncertainties for KT are less than 5%, and for KQ are
less than 3%. Even for the fine grid resolution, the grid uncertainty is close to the
other grid resolutions. The validation error for KT and KQ is less than 4.0%. The
convergences of KT , 10KQ and ηo for J = 0.7 are shown in Figure 4.7. The residual
convergences are shown in Figure 4.8. The comparison of open water prediction with
the measured data performance is shown in Figure 4.9. Good agreement was found for
KT , 10KQ , whereas the efficiency η0 shows deviation for higher advance coefficients.

4.1.3 Body-Force Method

4.1.3.1 Computational Domain and Grid Generation

For the body-force modelling of the KP505 propeller, a rectangular computational

domain was created for open water simulation. The grid is shown in Figure 4.10. The

131
 1.4 100
KT - Experimental
KT - Numerical
10KQ - Experimental
1.2 10KQ - Numerical
ηo - Experimental 80
ηo - Numerical
KT, 10KQ 1.0

0.8 60

ηo
0.6 40
0.4
20
0.2

0.0 0
0 0.2 0.4 0.6 0.8 1
Advance Coefficient, J

Figure 4.9: Open water simulation of KP505 using transient solver for GRID-B.

inlet boundary was set at 4D and the outlet boundary was set at 8D from the propeller
plane. The side boundaries were set at 2D from the centre plane. A cylindrical cell
zone was defined for the propeller disk where the body-force was applied to the cells.

4.1.3.2 Grid Dependency Study

The effect of the number of cells (MR ) along the radial direction of the propeller disk
was studied using the transient solver with local velocity body-force model, LV-BFM
for advance coefficient J = 0.70. The sampling was set d = 0.3D upstream from
the propeller plane. In numerical simulation, MR = 8, 16, 32, and 64 were used for
the grid dependency study. Figure 4.11 shows the effect of the number of cells on
numerical prediction accuracy with respect to the experimental data. It is found that
MR = 16 and higher shows good convergence with the experimental results.

132
Figure 4.10: Computational domain for KP505 open water simulation using
body-force method.

0.8
KT - Experimental
KT - Numerical
0.7 10KQ - Experimental
10KQ - Numerical

0.6

0.5
KT, 10KQ

0.4

0.3

0.2

0.1

0.0
0 10 20 30 40 50 60 70 80
Number of cells in radial direction, MR

Figure 4.11: Effect of number of cells in radial direction.

4.1.3.3 Sampling Plane Distance

The effect of the velocity sampling plane distance on the propeller performance
predicted by the body-force method was studied from d/D = 0.10 to d/D = 0.40,

133
 0.8
KT - Experimental
KT - Numerical
0.7 10KQ - Experimental
10KQ - Numerical

0.6

0.5
KT, 10KQ

0.4

0.3

0.2

0.1

0.0
0 0.1 0.2 0.3 0.4 0.5
Velocity sampling plane distance, d/D
Figure 4.12: Effect of velocity sampling plane distance.

1.4 100
KT - Experimental
KT - Numerical
1.2 10KQ - Experimental
10KQ - Numerical
η - Experimental 80
η - Numerical
1.0
KT, 10KQ

0.8 60

0.6 40 ηo
0.4
20
0.2

0.0 0
0 0.2 0.4 0.6 0.8 1
Advance Coefficient, J

Figure 4.13: Propeller open water characteristics of KP505 using body-force method.

134
where d is the offset distance of the sampling plane from the propeller plane, and
D is the diameter of the propeller. The LV-BFM was used for body-force modelling
without induced velocity correction. Comparison between the numerical data and
experimental values is shown in Figure 4.12 for J = 0.70. The influence of the
propeller disk position on the incoming velocity is higher for lower offset distances.
The influence effect is not significant for d/D = 0.30 and higher.

Then MR = 32 was used for different advance coefficients for the validation of the
open water data. Comparison between the numerical results and experimental data
shows good agreement in Figure 4.13.

135
4.2 MP687 Test Case

4.2.1 Propeller Geometry

The propeller geometry for the self-propulsion simulation of the JBC hull is shown
in Figure 4.14. The propeller is designated MP687. The principal particulars of the
MP687 propeller are listed in Table 4.5. The surface file for the propeller geometry
and open water data were obtained from the T2015 workshop website. The open
water performance tests were carried out at NMRI, Japan.

Figure 4.14: MP687 propeller geometry.

4.2.1.1 Computational Domain and Grid Generation

A similar cylindrical computational domain to that used for the KP505 propeller was
also used for the grid generation of MP687 propeller for open water simulation. The

136
 Figure 4.15: Grid distribution on the propeller blade and hub of MP687.

grid distribution on the propeller blade is shown in Figure 4.15. The computational
grids for open water simulations of the MP687 propeller are summarized in Table 4.6.

4.2.1.2 Steady Simulation

The MRF method was used to perform the steady simulation of the MP687 propeller
using the steady solver. The grid uncertainty for the numerical results was studied

Table 4.5: Principal particulars of MP687 propeller.

Item Full scale Model scale

Type FP FP
No. of blades 5 5
D (m) 8.12 0.203
P/D(0.7R) 0.75 0.75
Ae/A0 0.5 0.5
Rotation Right hand Right hand
Hub ratio 0.18 0.18

137
 Table 4.6: Computational grids for MP687 propeller open water simulation

Grid Stator Rotor Average y +

GRID-A 797083 547454 20.03
GRID-B 797083 1632402 20.57
GRID-C 797083 2736403 20.71
GRID-D 797083 3113553 15.20

Table 4.7: Grid dependency of open water performance at J = 0.7 using steady solver.

GRID hi /h1 KT |E%D| |USN %S| 10KQ |E%D| |USN %S|

P = 1.79 P =1.88
GRID-A 1.43 0.0864 0.29 0.63 0.1669 10.62 6.72
GRID-B 1.17 0.0870 0.38 0.08 0.1640 8.70 4.97
GRID-C 1.03 0.0862 0.55 0.18 0.1623 7.54 4.13
GRID-D 1.00 0.0876 1.01 0.23 0.1637 8.50 3.90

0.8
KT
10KQ
ηo
0.7

0.6
KT, 10KQ , ηo/100

0.5

0.4

0.3

0.2

0.1

0.0
0 200 400 600 800 1000 1200 1400 1600 1800 2000
Number of Iteration

Figure 4.16: Convergence history of KT , KQ and η0 using steady solver for GRID-B
at J = 0.70.

138
 1
Ux
Uy
0.1 Uz
p
k
omega
0.01
Initial residuals

0.001

0.0001

1e-05

1e-06

1e-07
0 200 400 600 800 1000 1200 1400 1600 1800 2000
Iterations

Figure 4.17: Residuals convergence for advance coefficient J = 0.70 for GRID-B.

using the four grids in Table 4.6 for J = 0.70. The validation error and the uncertainty
of the numerical results is shown in Table 4.7. The validation error of KT for the four
grids is around 1% and for KQ is around 8% except for the coarse grid. The numerical
uncertainty for the fine grids are similar. The convergences of KT , KQ and η0 for
GRID-B are shown in Figure 4.16. The residuals are also plotted in Figure 4.17. The
open water performance was simulated using the steady solver for different advance
coefficients and the results are compared with experimental data in Figure 4.18.

4.2.1.3 Unsteady Simulation

The unsteady numerical open water simulations for the MP687 propeller were also
carried out using transient solver pimpleDyM F oam for the four grids at J = 0.70.
The validation error and numerical uncertainty are presented in Table 4.8. The
convergences of force coefficients are shown in Figure 4.19. The residuals are also

139
 1.0 100
KT - Experimental
KT - Numerical
10KQ - Experimental
10KQ - Numerical
0.8 ηo - Experimental 80
ηo - Numerical

0.6 60
KT, 10KQ

ηo
0.4 40

0.2 20

0.0 0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
Advance Coefficient, J
Figure 4.18: Open water simulation of MP687 propeller using steady solver for
GRID-B.

shown in Figure 4.20. The change in residuals after 0.2 is not significant and the
solution can be considered converged. Figure 4.21 shows the comparison of KT , 10KQ
and η with the experimental data for a range of advance coefficients for GRID-B, and
the comparison shows better agreement than steady simulation.

Table 4.8: Grid dependency of open water performance at J = 0.7 using unsteady
solver.

GRID hi /h1 KT |E%D| |USN %S| 10KQ |E%D| |USN %S|

P = 1.38 P =1.92
GRID-A 1.43 0.0863 0.41 7.92 0.1654 9.64 3.96
GRID-B 1.17 0.0889 2.56 5.51 0.1648 9.21 2.88
GRID-C 1.03 0.0880 1.45 4.45 0.1630 8.00 2.41
GRID-D 1.00 0.0895 3.20 4.11 0.1635 8.33 2.28

140
 0.8
KT
10KQ
ηo
0.7

0.6
KT, 10KQ , ηo/100

0.5

0.4

0.3

0.2

0.1

0.0
0 0.05 0.1 0.15 0.2 0.25 0.3
Time (sec)

Figure 4.19: Convergence history of KT , KQ and η0 using transient solver for GRID-B
at J = 0.70.

1
Ux
Uy
0.1 Uz
p
k
omega
0.01
Initial residuals

0.001

0.0001

1e-05

1e-06

1e-07
0 0.05 0.1 0.15 0.2 0.25 0.3
Time (sec)

Figure 4.20: Residuals convergence for GRID-B at advance coefficient J = 0.70.

141
 1.0 100
KT - Experimental
KT - Numerical
10KQ - Experimental
10KQ - Numerical
0.8 ηo - Experimental 80
ηo - Numerical

0.6 60
KT, 10KQ

ηo
0.4 40

0.2 20

0.0 0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
Advance Coefficient, J
Figure 4.21: Open water simulation of MP687 propeller using transient solver for
GRID-B.

142
4.3 NRC-IOT Propeller Test Case

4.3.1 Geometry

A rectangular computational domain was created similar to the KP505 propeller

domain to study the open water test of the NRC-IOT propeller using the body-force
method. The principal particulars are listed in Table 4.9. A sampling plane for the
velocity inlet was set 1.0D upstream from the propeller plane.

Table 4.9: Principal particulars of the NRC-IOT propeller.

Item Model scale

Type FP
No. of blades 4
D (m) 0.1205
P/D(0.7R) 1.27
Ae/A0 0.906

1.0
KT -Experimental
KT -Numerical
10KQ -Experimetal
10KQ -Numerical
0.8

0.6
KT, 10KQ

0.4

0.2

0.0
0 10 20 30 40 50 60 70 80
Number of cells in radial direction, MR

Figure 4.22: Effect of the cell numbers on propeller parameters at J = 0.70 using
transient solver.

143
 1.0
KT -Experimental
KT -Numerical
10KQ -Experimental
10KQ -Numerical
0.8

KT, 10KQ 0.6

0.4

0.2

0.0
0 0.1 0.2 0.3 0.4 0.5
Velocity sampling plane distance, d/D

Figure 4.23: Effect of the velocity sampling plane distance on propeller parameters
at J = 0.70 using transient solver.

4.3.2 Numerical Results

The effect of number of cells in the radial direction of the propeller disk was studied
for four resolutions, MR = 8, 16, 32 and 64. In these cases, the velocity sampling plane
was set at 1.0D upstream from the propeller plane. The simulations were performed
for advance coefficient, J = 0.70. For grid resolution 16 and more along the radial
direction, the numerical results show good agreement with the experimental data as
seen in Figure 4.22. The effect of velocity sampling plane distance from the propeller
plane was also investigated. The open water performance was calculated for d/D =
0.1, 0.2, 0.3 and 0.5 using the LV-BFM without induced velocity correction. The
numerical results are compared with the experimental data in Figure 4.23. The effect
of sampling plane distance greatly influences the numerical prediction when it is close
to the propeller disk due to the higher input velocity to the propeller disk. Then the
open water simulations for the other advance coefficients were carried out using d/D =

144
 1.5 100
KT -Experimental
KT -Numerical
10KQ -Experimental
10KQ -Numerical
1.2 η - Experimental 80
η - Numerical

KT, 10KQ 0.9 60

ηo
0.6 40

0.3 20

0.0 0
0.2 0.4 0.6 0.8 1
Advance Coefficient, J

Figure 4.24: Comparison of open water data of NRC-IOT propeller using transient
solver.

1.0 and MR = 32. The comparison of numerical results and experimental results is
shown in Figure 4.24. The numerical results agreed well with the experimental data.

145
This page intentionally left blank.
Chapter 5

Numerical Wave Tank Validation

The numerical wave tank validation is an important step before application in

wave-body interaction simulations. In this chapter, the available generation modules
are studied for a steep wave. Depending on the numerical results obtained by different
modules the most efficient module is adopted for the grid dependency and time
dependency studies. The effect of turbulence model on numerical wave generation
was also studied. Based on the findings of the grid dependency and time step size
dependency analysis, the optimum settings are used for wave run-up study of a
truncated square cylinder (Ali et al., 2017a). The detailed numerical results validation
can be found in Appendix-B. Since the grid resolution and time step size depends
on the wave steepness and wave length, the wave validation also was studied for the
waves prior to the self-propulsion simulation in waves.

To study the effect of ship motion on propeller-hull interaction in waves, the numerical
wave tank without ship model was validated for the waves listed in Table 5.1 for the
KCS model. The numerical wave elevations at the three different locations as shown
in Figure 5.1 were compared with the analytical wave elevation. The details of the
validation process are described in the following subsections.

147
5.1 Computational Domain

The schematic diagram of the numerical wave tank is shown in Figure 5.1. The
computational domain with grid distribution for the wave validation with the virtual
ship model is shown in Figure 5.2. The inlet boundary was set 1.0Lpp in front of the
ship model. The outlet boundary was set 3.0Lpp behind the ship model, and the side
boundaries were set 1.5Lpp from the center plane. The wave was generated at the
inlet boundary and the numerical wave damping zone was used for wave damping in
front of the outlet. The bottom boundary was set 1.5Lpp from the calm water surface.

Figure 5.1: Computational domain for numerical wave tank validation.

Table 5.1: Wave conditions for numerical tank validation

No C1 C3
Wave length, λ (m) 4.7351 8.3683
Wave height,H (m) 0.0743 0.1475

148
 Figure 5.2: Grid distribution in the computational domain.

The same grid resolutions as the grid for the bare hull resistance study were used
for wave validation. The computational grids were generated using blockMesh and
navalSnappyHexMesh. For incoming wave simulation, a computational domain with
a symmetry plane was generated to reduce the computational cost. The symmetry
plane was set at the centreplane.

5.2 Grid Dependency

The sensitivity of grid resolution on the numerical wave elevation was studied for
the wave C1 with three grids resolutions having λ/∆x = 40, 50 and 60 and H/∆z
= 5, 6, and 7, where is H is the wave height. The normalized wave elevation with
respect to the wave amplitude A1 = H/2.0 are compared at longitudinal position
of the forward perpendicular (WP-01), midship section (WP-02), and at the aft
perpendicular (WP-03). The numerical simulations for grid dependency were carried
out using time step size TS-03 (Tw /∆t = 500). Figure 5.3 shows the comparison

149
 2.0
GRID-A
GRID-B
GRID-C
1.5 Analytical

1.0
Wave Height, η / A1

0.5

0.0

-0.5

-1.0

-1.5
15 16 17 18 19 20
Time, t/Tw

Figure 5.3: Comparison of wave elevation at WP-02 using TS03 of wave C1.

between the analytical and numerical prediction. It is observed that the GRID-C
prediction shows good agreement with the analytical values.

5.3 Time Step Size Dependency

For the study of the influence of time step size on the numerical wave generation three
time step sizes, TS01(Tw /∆t = 125), TS02(Tw /∆t = 250), TS03(Tw /∆t = 500), for
GRID-C were used in numerical simulations. The comparison is shown in Figure 5.4
for wave probe WP-02. It is seen that TS-03 predicts the wave elevation accurately.
From the comparison, it was found that both TS-01 and TS-02 under predict the
wave elevation.

150
 2.0
TS-01
TS-02
TS-03
1.5 Analytical

Wave Height, η / A1 1.0

0.5

0.0

-0.5

-1.0

-1.5
15 16 17 18 19 20
Time, t/Tw

Figure 5.4: Comparison of wave elevation at WP-02 using GRID-C of wave C1.

2.0
Numerical
Analytical

1.5

1.0
Wave Height, η / A1

0.5

0.0

-0.5

-1.0

-1.5
15 16 17 18 19 20
Time, t/Tw

Figure 5.5: Comparison of wave elevation at WP-02 using GRID-C of wave C3.

151
5.4 Wave Validation

For the self-propulsion simulation in waves, the C3 wave was validated against the
analytical wave elevation. The grid resolution for C3 is λ/∆x = 100 and H/∆z = 14
with the time step size Tw /∆t = 500 for the numerical wave simulation. The numerical
wave elevation at the midship section is compared with the analytical wave elevation
in Figure 5.5. The wave crests are slightly over predicted compared with analytical
wave elevation.

152
Chapter 6

Self-Propulsion Simulation

This chapter presents the validation of numerical self-propulsion simulation of the

selected ship models. The numerical self-propulsion simulations for the KCS and JBC
model were carried out using detailed propeller geometry. The developed body-force
models were also used for self-propulsion simulations for the KCS model and the
Fishing vessel. For the self-propulsion simulation with detailed propeller geometry,
the developed dynamic motion class was used. Since the code was developed for the
constant rotational speed of the propeller, fixed rpm was used for propeller rotation.

6.1 KCS Test Case

The self-propulsion simulations for the KCS model were performed for the even keel
condition, with sinkage and trim, and in head wave conditions. The test conditions
for the self-propulsion simulations are presented in Table 6.1.

153
 Table 6.1: Test cases for the validation of self-propulsion simulation of KCS model.

Case 6.1 6.2 6.3

Model KCS KCS KCS
Wave Calm Calm Head sea
Condition Self-prop (ship) Self-prop (ship) Self-prop (ship)
LP P 7.2786 7.2786 7.2786
Fn 0.26 0.26 0.26
Rn × 107 1.26 1.26 1.26
Rudder without without without
Propeller with with with
∗
Attitude F R0 F Rzθ F Rzθ
154

Validation variables Thrust, resistance,

wave elevation, prop factor,
mean flow velocity
EFD provider NMRI

∗
Hull Attitude Definitions:
F R0 = Zero sinkage and zero trim at the design speed.
F Rzθ = Free to heave and pitch.
6.1.1 Geometry

The self-propulsion for the KCS ship model was carried out for fixed rpm of the
propeller. The ship model was equipped with the KP505 propeller as shown in Figure
6.1. The rudder was not included with the ship model during the model test.

6.1.2 Computational Domain and Boundary Conditions

From the bare hull resistance validation and propeller open water test validation,
the computation grid for numerical self-propulsion simulation, the optimum bare hull
grid was combined with the propeller grid. For the self-propulsion simulation of the
KCS model, the bare hull grid GRID-C was combined with propeller grid GRID-B.

Figure 6.1: Geometry for KCS self-propulsion.

Figure 6.2: Computational grid for self-propulsion of KCS model.

155
Due to the non-symmetric condition of the propeller wake, the whole computational
domain was used in this numerical simulation. The same coordinate system as the
bare hull case was used for the self-propulsion simulation. The grid distribution using
the actual propeller geometry in shown Figure 6.2. For the self-propulsion simulation
using the body-force method, the bare hull grid GRID-C was used. To reduce the
computational time, a steady simulation was carried out initially, and then the steady
solution was used as the initial condition for the unsteady simulation.

The self-propulsion simulation using the body-force model, as descried in section 2.8
was carried out for n = 9.5 at the design speed F n =0.26. The incoming velocity to
the propeller domain was the average of the local velocity on the sampling plane which
was set in the upstream of the propeller domain. The local velocity on the sampling
plane was also used as input velocity in the propeller body-force calculation to improve
the velocity distribution and the propulsion factor prediction. The diameter of the
actuator disk was set equal to the diameter of the propeller and the longitudinal
extent is the width of the projected chord length on the propeller plane.

6.1.3 Self-Propulsion at Restrained Even Keel Condition

The self-propulsion simulations were performed using the detailed propeller geometry,
the body-force model based on fixed input velocity, and the body-force model using
local velocity. The numerical results from the self-propulsion simulation are compared
in Table 6.2. For turbulence modelling, k − ω SSTCC and k − ω SST were used
for detailed propeller cases. The Skin Friction Correction (SF C = RT (SP ) − T )
using k − ω SSTCC turbulence model with detailed propeller geometry shows less
than 1.0 % error against the experimental values. However, the thrust and torque
were over predicted by a small amount. From the comparison it can be concluded

156
 (a) Experimental

(b) Detailed propeller geometry

(c) FV-BFM

(d) LV-BFM

Figure 6.3: Dynamic pressure distribution on the stern of KCS for self-propulsion.

157
 (b) Detailed propeller geometry.
(a) Experimental

(c) FV-BFM (d) LV-BFM

Figure 6.4: Axial velocity distribution at 0.25D downstream of the propeller plane.

Table 6.2: Numerical results of self-propulsion simulation.

Method RT (SP ) (N) CT Thrust, T (N) Torque, Q (N-m) RT (SP ) − T (N) |E%D| KT KQ VA (m/s) w
k − ω SST 91.05 4.00E-03 62.79 2.7314 28.26 1.990 0.1781 0.0310 1.695 0.2281
k − ω SSTCC 92.88 4.08E-03 63.49 2.7503 29.39 0.852 0.1801 0.0312 1.689 0.2309
LV-BFM 89.00 3.91E-03 57.04 2.4880 31.96 1.712 0.1620 0.0280 1.776 0.1910
FV-BFM 88.55 3.89E-03 71.66 2.9391 17.09 13.359 0.2033 0.0334 1.570 0.2851
Experimental 3.97E-03 30.25 0.1700 0.0288

that self-propulsion with detailed propeller geometry and k − ω SSTCC shows better
accuracy than the other approaches although the computation time is longer than that
for the simple body-force method. The simulation time can be reduced considerably
by using the steady solution as the initial condition for the unsteady simulation.

158
6.1.3.1 Pressure Distribution

Figure 6.3 illustrates the pressure distribution on the stern region of the KCS model
using the detailed propeller geometry and two body-force models. A difference is
observed on the lower region of the stern bulb due to the incoming flow alteration by
the propeller action itself. The overall pressure distribution on the hull surface shows
a good agreement with experimental results.

6.1.3.2 Flow Velocity

The axial velocity distribution at 0.25D downstream of the propeller plane (x/LP P
= 0.9910 from the forward perpendicular) with the working propeller is compared
in Figure 6.4 with the experimental data. The comparison shows that, although
the body-force method can predict the KT , KQ , and effective wake fraction w with
reasonable accuracy, it failed to predict the velocity distribution.
The velocity components obtained for self-propulsion simulation using k − ω SSTCC
with working propeller are compared at 0.25D behind the propeller plane and
z/LP P = −0.0291 with the experimental data in Figure 6.5. The comparison shows
higher deviation near the hub center. The results indicate the hub vortex strength is
under predicted.

159
 1

0.5
u/U,v/U,w/U

0
u/U-Actual propeller
v/U-Actual propeller
w/U-Actual propeller
u/U-EFD
v/U-EFD
-0.5 w/U-EFD
u/U-Body-Force
v/U-Body-Force
w/U-Body-Force

-0.03 -0.02 -0.01 0 0.01 0.02 0.03

y/LPP

Figure 6.5: Velocity components comparison at x/LP P = 0.9910 (0.25D downstream

of the propeller plane) and z/LP P = -0.0291.

160
6.1.4 Self-Propulsion with Sinkage and Trim

The effect of ship motion on self-propulsion parameters was investigated for the KCS
hull with the KP505 propeller. The computational grid used for self-propulsion for
the even keel condition was used for numerical simulation of self-propulsion with ship
motions. The ship model was free to heave and pitch only. Mesh morphing method
was adopted to consider the ship motion in numerical simulations. The simulation
was carried out using the k − ω SSTCC turbulence model. The results from the
numerical simulation are presented in Table 6.3. For the dynamic case, the skin
friction correction SF C is about 21% higher than the even keel condition. The wake
fraction w, for the dynamic condition is also reduced at the same ship speed and
propeller rotational speed. In the case of propeller design, these changes should be
included in the optimum propulsive system design considering the ship motion.

Figure 6.6 presents the time history of residual convergence for the numerical
simulation when the ship model is allowed heave and pitch motions. The sudden
change in residual parameters is due to the restart of the case from the previous
saved data or changes in time step size. The residuals for alpha.water, k and omega
converged to 10−6 . For pressure the residual obtained was 10−4 . The convergence
error, E%C, for total resistance, RT (SP ) is shown in Figure 6.7. The error was
calculated based on the average of 250 data points, is less than 1.0%. Figure 6.8
shows the convergence error for propeller thrust, which is much lower than 1.0% for
the average of 250 data points.

Table 6.3: Numerical results of self-propulsion test with ship motion.

Condition RT (SP ) (N) CT Thrust, T (N) Torque, Q (N-m) RT (SP ) − T (N) VA (m/s) w
Even keel 92.88 4.08E-03 63.4847 2.7503 28.8331 1.689 0.2309
Dynamic 96.83 4.26E-03 61.8382 2.2684 34.9949 1.711 0.2209

161
 1
p_rgh
0.1 alpha.water
k
omega
0.01

0.001
Initial residuals
0.0001

1e-05

1e-06

1e-07

1e-08

1e-09
0 2 4 6 8 10 12 14
Time t

Figure 6.6: Residual convergence for self-propulsion with sinkage and trim.

Avg. over 250 data

4

2
E%C

−2

−4

0 2 4 6 8 10 12 14
Time, (s)

Figure 6.7: Total resistance convergence error for self-propulsion with sinkage and
trim.

162
 Avg. over 250 data
4

2
E%C

−2

−4

0 2 4 6 8 10 12 14
Time, (s)

Figure 6.8: Propeller thrust convergence error for self-propulsion with sinkage and
trim.

6.1.4.1 Pressure Distribution

The effect of ship motion on pressure distribution on the hull was investigated for
the self-propulsion simulation with ship motion. The comparison of dynamic pressure
coefficient (Cp ) is shown in Figure 6.9. For the dynamic condition, the pressure
distribution is found more asymmetric than even keel condition on the hull surface
above the propeller, which will lead to vibration in the hull. The detailed hull pressure
fluctuation due to the propeller rotation for the even keel condition and the dynamic
condition was investigated by measuring the pressure on the hull surface using pressure
probes as shown in Figure 6.10. The pressure gauges P 1, P 2, P 3 were placed on the
center line, P 4 on starboard, and P 5 on the port side. Pressure gauge P 2 was placed
at the intersection of the propeller plane and the ship center plane. The pressure
fluctuations due to the propeller rotation are compared in Figure 6.11 for all the

163
 (a) Sinkage and trim constrained to zero.

(b) Free sinkage and trim.

Figure 6.9: Pressure distribution on the hull and the propeller.

164
 (a) Plan view (b) 3D view

Figure 6.10: Virtual pressure gauge locations.

probes. A larger difference in pressure coefficient is observed in the downstream

direction than that for the sideways direction. No difference in pressure coefficient is
observed upstream of the propeller due to the ship motion.

6.1.4.2 Flow Velocity

The normalized axial velocities (with respect to the ship velocity) on the center plane
are compared in Figure 6.12. The axial velocity distribution above the propeller tip
shows larger negative value than the even keel condition due to the ship motions. The
vortex core trajectory is also stretched towards the propeller rotational axis due to
the ship motions. The axial velocity distribution is also compared at x/LP P = 0.9910
looking upstream in Figure 6.13. The axial velocity distribution downstream of the
propeller plane shows a small difference.

165
 0.26
Probe-P1 Even keel
0.24 Dynamic

0.22

Cp
0.20
0.18
0.16
0 50 100 150 200 250 300 350
Blade angle, (deg)

0.26
Probe-P2 Even keel
0.24 Dynamic

0.22
Cp

0.20
0.18
0.16
0 50 100 150 200 250 300 350
Blade angle, (deg)

0.26
Probe-P3 Even keel
0.24 Dynamic

0.22
Cp

0.20
0.18
0.16
0 50 100 150 200 250 300 350
Blade angle, (deg)

0.26
Probe-P4 Even keel
0.24 Dynamic

0.22
Cp

0.20
0.18
0.16
0 50 100 150 200 250 300 350
Blade angle, (deg)

0.26
Probe-P5 Even keel
0.24 Dynamic

0.22
Cp

0.20
0.18
0.16
0 50 100 150 200 250 300 350
Blade angle, (deg)

Figure 6.11: Numerical pressure measurement at the pressure gauges.

166
 (a) Sinkage and trim constrained to zero.

(b) Free sinkage and trim.

Figure 6.12: Axial velocity distribution on the centreplane.

167
 (a) Sinkage and trim constrained to zero.

(b) Free sinkage and trim.

Figure 6.13: Axial velocity distribution at x/LP P = 0.9910 (0.25D downstream of the
propeller plane).

168
6.1.5 Seakeeping Simulation

Before proceed to the self-propulsion simulation in wave, the bare hull simulation was
carried out in wave with forward speed corresponding to F n = 0.26 for one wave
condition. A similar computational domain as the bare hull resistance simulation was
used for the seakeeping simulation. The wave generation and and wave absorption
conditions were set up for the incoming wave and wave damping. Seakeeping
simulation was performed for 10 wave encounter periods for wave λ/Lpp = 1.15.
Figure 6.14 shows the comparison of the total resistance comparison in the head sea
condition and the even keel condition. In compare to the even keel condition, the
total resistance coefficient was increased by about 75 %. The motion response are
compared with the experimental data (T2015, 2015) in Figure 6.15. The numerical
simulation motion response shows a good agreement with the experimental data.

0.020
CT
CT - Mean
0.015 CT - Even Keel (No wave)
CT - FORCE
Total resistance coeff, CT

CT - FORCE-Mean
0.010

0.005

0.000

-0.005

-0.010
5 5.5 6 6.5 7 7.5 8
t/Te

Figure 6.14: Comparison of total resistance coefficient in wave C3.

169
 1.50
Experimetal
Numerical
1.00

0.50
Heave, z/ ζs

0.00

-0.50

-1.00

-1.50
0 1 2 3 4 5
t/Te

(a) Heave motion.

1.50
Experimetal
Numerical
1.00

0.50
Pitch, θ/ k ζs

0.00

-0.50

-1.00

-1.50
0 1 2 3 4 5
t/Te

(b) Pitch motion.

Figure 6.15: Motion response of KCS bare hull in wave C3.

170
6.1.6 Self-Propulsion in Waves

The self-propulsion simulation of the KCS model in wave was performed with the
detailed propeller geometry. A regular incoming wave of λ/LP P = 1.15 and HS /λ =
1/56.74 with current corresponding to the design speed was imposed on the inlet
boundary. The numerical wave tank validation of the selected wave was studied
prior to application in the self-propulsion simulation. The effect of the wave on the
propeller-hull interaction was investigated using the developed dynamic class based
on the sliding mesh approach. The numerical results are presented in Table 6.4. The
results show an increase in total resistance coefficient and decrease in wake fraction
compared to the even keel condition. The residual convergence history is shown in
Figure 6.16. The residuals converged to 10−2 ∼10−6 . Figure 6.17 presents the time
history of total resistance coefficient, CT , for last two wave encounter periods. For the
last two wave encounter periods, the resistance curve shows a cyclic pattern, which
can be considered as a converged simulation. The variation of thrust coefficient KT
and torque coefficient KQ over one wave encounter period is shown in Figure 6.18 and
6.19 respectively. The average of KT and KQ is lower than the even keel condition for
the selected wave condition. The generated thrust from the propeller is reduced by
20%, and the torque is also reduced by 13% which causes reduction of the propulsion
efficiency.

Table 6.4: Numerical results of self-propulsion test with ship motion in wave C3.

Wave crest RT (N) CT Thrust, T (N) Torque, Q (N-m) RT − T (N) VA (m/s) w

At FP 166.53 7.325E-3 45.0574 2.3376 121.47 1.938 0.1175
Midship 158.12 6.955E-3 54.0788 2.4777 104.04 1.812 0.1758
Even keel 92.88 4.08E-03 63.4847 2.7503 28.83 1.689 0.2309

171
 1

0.01

0.0001
Initial residuals

1e-06

1e-08

1e-10

1e-12 p_rgh
alpha.water
k
omega
1e-14
0 2 4 6 8 10 12
Time t

Figure 6.16: Residual convergence for self-propulsion in wave C3.

0.020
Experimental:: Towed
Numerical:: Towed
0.015 Numerical:: Self-propelled
Total resistance coeff, CT

0.010

0.005

0.000

-0.005

-0.010
6 6.5 7 7.5 8
t/Te

Figure 6.17: Total resistance coefficient for self-propulsion in wave C3.

172
 0.25
Numerical :: In waves
Numerical :: In waves (mean)
Experimental :: Even keel (No wave)
0.20

Thrust coeff, KT
0.15

0.10

0.05

0.00
7 7.2 7.4 7.6 7.8 8
t/Te

Figure 6.18: Thrust coefficient over one encounter wave period for self-propulsion in
wave C3.

0.50
Numerical :: In waves
Numerical :: In waves (mean)
Experimental :: Even keel (No wave)
0.40
Torque coeff, 10KQ

0.30

0.20

0.10

0.00
7 7.2 7.4 7.6 7.8 8
t/Te

Figure 6.19: Torque coefficient over one encounter wave period for self-propulsion in
wave C3.

173
6.1.6.1 Ship Response

The response of the model with a rotating propeller in a wave with forward speed
is shown in Figure 6.20. The simulation was carried out for eight encounter wave
periods. A regular response of the model was observed for heave and pitch motion.

6.1.6.2 Free Surface Wave

Figure 6.21 illustrates the generated periodic wave pattern over one encounter wave
period which is the combination of an incoming wave and the wave generated by the
model. Since λ/LP P = 1.15, when the bow is in the middle of the wave crest, the
stern region is also in a wave crest. At t/Te = 0, the bow is deeply submerged in
the wave and the bow flare generates more waves which results in more wave making
resistance. As the wave propagates, the bow comes out from the water and generates
less wave. At t/Te = 1/2, the model is basically in the hogging condition which means
less frictional resistance and less wave making resistance. As the wave propagation
continues, the bow again goes into wave crest and generates a similar wave pattern
as t/Te = 0.

174
 1.50
Numerical

1.00

Heave, z/ ζs 0.50

0.00

-0.50

-1.00

-1.50
0 1 2 3 4 5 6 7 8 9
t/Te

(a) Heave motion.

1.50
Numerical

1.00

0.50
Pitch, θ/ k ζs

0.00

-0.50

-1.00

-1.50
0 1 2 3 4 5 6 7 8 9
t/Te

(b) Pitch motion.

Figure 6.20: Motion response of KCS hull during self-propulsion in wave C3.

175
 (a) t/Te = 0.0

(b) t/Te = 1/4

(c) t/Te = 1/2

(d) t/Te = 3/4

Figure 6.21: Free surface wave contour over one encounter wave period, Te for
self-propulsion in wave C3.

176
 (a) t/Te = 0.0

(b) t/Te = 1/4

(c) t/Te = 1/2

(d) t/Te = 3/4

Figure 6.22: Pressure distribution on hull surface over one encounter wave period, Te
for self-propulsion in wave C3.

6.1.6.3 Pressure Distribution

The pressure coefficient, CP , on the stern part of the model is compared in Figure
6.22 over one encounter wave period. Since the bow is deeply submerged in the wave

177
crest at t/Te = 0.0, the bow experiences high pressure. The stern part of the model
also experiences high pressure due to the wave crest present in the stern region. The
bow starts to experience less pressure since it starts to emerge from the wave crest
which can be seen at t/Te = 1/4. As the wave crest moves to the midship section,
the bow and stern region exerts low pressure at t/Te = 1/2. At t/Te = 3/4, the bow
region and stern region start to experience high pressure again due to the incoming
wave crest. A similar pressure distribution was observed for the next incoming wave.

178
6.2 JBC Test Case

The JBC test case allowed comparison between the full numerical self propulsion
simulation and a set of experimental data that included measurements of flow velocity
and vorticity around the propeller. The test conditions for the self-propulsion of JBC
model are listed in Table 6.5. In numerical simulation, fixed rpm was used for propeller
rotation.

6.2.1 Geometry

The self-propulsion simulation for the JBC ship model was also performed with a fixed
rpm for the propeller. The ship model was equipped with the MP687 propeller as
shown in Figure 6.23. The rudder and energy saving device (ESD) were not included
with the model during the experiments, and also not included in the numerical
simulation.

Table 6.5: Test conditions for JBC self-propulsion simulation

Case 6.4
Wave Calm
Condition Self-prop (ship)
LP P 7.00 m
Fn 0.142
7
Rn × 10 0.746
Rudder without
Propeller with
Attitude sinkage & trim
Validation variables prop factor, wave elevation
mean flow velocity
EFD provider NMRI

179
6.2.2 Computational Domain and Grid Generation

As a first step in the self-propulsion simulation, the bare hull resistance and
propeller open water tests were validated. Based on the numerical results analysis,
the computation grid for the numerical self-propulsion simulation was generated
combining the optimum bare hull grid with the propeller grid. Bare hull grid, GRID-C
was combined with propeller grid GRID-B for the computational grid in the JBC
self-propulsion simulation. Due to the non-symmetric condition of the propeller wake,
the whole computational domain was used in the numerical simulation. The grid
distribution of the self-propulsion simulation using the actual propeller geometry is
shown in Figure 6.24. To reduce the computational time, a steady simulation was
carried out initially, and then the steady solution was used as the initial condition for
the unsteady simulation.

6.2.3 Results Validation

The model was free to heave and pitch during the model test as well as in the numerical
simulations. The numerical simulation was started with the larger time step size.
When the ship motion was stable, the time step size was reduced to the corresponding
1.00 rotation of the propeller for every time step. The effect of the change in time step

Figure 6.23: Geometry for JBC model for self-propulsion simulation.

180
 Figure 6.24: Computational grid for self-propulsion of JBC model.

size is clearly seen in the residual convergence which is shown in Figure 6.25. The time
step size was changed after about 8.40 sec and the residuals dropped another order
of magnitude. The convergence error for total resistance is calculated and presented
in Figure 6.26. It was found that the convergence error at the lower time step was
much lower than 1.0%. The numerical results for the self-propulsion simulation are
compared with the experimental data in Table 6.6. The comparison of CT against the
experimental data shows just above 5.0% error. The error in SF C was about 7.0%.
The overall comparison shows good agreement with the experimental results.

Table 6.6: Numerical results of self-propulsion of JBC model.

RT (SP ) (N) CT Thrust, T (N) Torque, Q (N-m) RT (SP ) − T (N) VA (m/s) w

Numerical 43.30 5.10E-03 26.45 0.6073 16.85 0.4909 0.5836
Experimental 40.82 4.81E-03 22.40 0.5846 18.20 0.5624 0.5230
|E%D| 6.08 7.39 11.60

181
 0.01
p_rgh
alpha.water
k
0.001 omega

Initial residuals 0.0001

1e-05

1e-06

1e-07

1e-08
0 1 2 3 4 5 6 7 8 9 10
Time t

Figure 6.25: Residual convergence of self-propulsion simulation of JBC model.

1.0
Avg. over 100 data

0.5
E%C

0.0

−0.5

−1.0
0 2 4 6 8 10
Time, (s)

Figure 6.26: Convergence of total resistance of JBC model.

182
 Figure 6.27: Blade angle definition.

6.2.4 Flow Velocity

The non-dimensional axial velocity distribution is validated with the experimental

data for three longitudinal locations at S2 (x/LP P = 0.9625), S4 (x/LP P = 0.9843),
and S7 (x/LP P = 1.00). For S2 and S7 sections, the numerical results are compared
for the blade angles 00 and 480 . The definition of blade angle is explained in Figure
6.27. Figure 6.28 shows the axial velocity distribution at x/LP P = 0.9625. The bilge
vortices in the numerical simulation separated earlier and moved further outwards
than the experimental vortices. However, the overall velocity distribution shows good
agreement with the experimental results. Figure 6.29 and 6.30 present the comparison
of experimental and numerical results at x/LP P = 0.9843 for blade angles 00 and 480 .
This section is just upstream of the propeller. It was expected to see the suction
effect of propeller on the velocity distribution and the numerical simulation was also
able to capture the effect. The comparison shows fairly good agreement with the

183
experimental results. The axial velocity was also compared just downstream of the
propeller (x/LP P = 1.0). The implemented dynamic motion class turbulence model
was able to capture the hub vortex more accurately and shows good agreement with
the experimental data for both blade angles.

184
 (a) Experimental.

(b) Numerical

Figure 6.28: Axial velocity distribution at x/LP P = 0.9625

185
 (a) Experimental.

(b) Numerical

Figure 6.29: Axial velocity distribution at x/LP P = 0.9843 for blade angle 00

186
 (a) Experimental.

(b) Numerical

Figure 6.30: Axial velocity distribution at x/LP P = 0.9843 for blade angle 480

187
 (a) Experimental.

(b) Numerical

Figure 6.31: Axial velocity distribution at x/LP P = 1.00 for blade angle 00 .

188
 (a) Experimental.

(b) Numerical

Figure 6.32: Axial velocity distribution at x/LP P = 1.00 for blade angle 480 .

189
Figure 6.33: Vortex structures based on Q = 500 for the self-propulsion simulation
with sinage and trim at F n = 0.142.

6.2.5 Vortex Visualization

The Q−criterion was used to identify the bilge vortex as well as vortices generated by
the propeller and hub. The iso-surface of vortices shown in Figure 6.33 was generated
using Q = 500. It is clearly seen that the bilge vortex was at the aft-shoulder and
it detached from the hull downstream of the propeller plane. In the absence of ESD,
the bilge vortices were contracted in the propeller disk. The vortices were generated
by the tips of the propeller blades, and the hub were well produced in the numerical
simulation. The vortex structures are coloured based on the pressure coefficient.

190
6.3 Fishing Vessel

The numerical self-propulsion simulations of the fishing vessel were performed using
the body-force method. The preliminary study is published by Ali et al. (2019b).
A whole computational domain was used for self-propulsion simulation due to the
asymmetric propeller force distribution. GRID-C was selected for self-propulsion
based on the uncertainty analysis and validation error in the numerical prediction. A
cylindrical refined domain was used as a propeller disk, which ensures the required grid
resolution based on open water data validation. During the numerical simulation, the
ship model was free to heave and pitch. The grid distribution with refined propeller
domain is shown in Figure 6.34. The numerical simulations were carried out using
the FV-BFM and LV-BFM methods.

Figure 6.34: Computational grid for self-propulsion of Fishing vessel.

191
 (a) Bare hull

(b) FV-BFM with n = 15.92

(c) LV-BFM with n = 15.82

Figure 6.35: Comparison of free surface wave contour for GRID-C at F n = 0.34 in
calm water.

192
 0.05
Bare Hull
FV-BFM
0.04 LV-BFM

Wave Elevation, ζ/Lpp 0.03

0.02

0.01

-0.01

-0.02

-0.03
-0.4 -0.2 0 0.2 0.4
x/Lpp

Figure 6.36: Comparison of wave elevation along the hull for GRID-C at F n = 0.34.

6.3.1 Wave Elevations

The free surface wave elevation behind the stern of the fishing vessel for the bare-hull
and the self-propulsion simulations are compared in Figure 6.35. It is observed that
the flow behind the vessel was accelerated due to the propeller action which increased
the wave elevation behind the transom. Differences were also observed in the flow
predicted by the two body-force models.

The non-dimensional wave elevations along the hull surface are compared in Figure
6.36 for the bare hull and the self-propulsion cases. The predicted wave elevation for
the self-propulsion case is slightly higher than that of the bare hull case near the bow
and is lower near the middle of the model. The difference in the wave elevation is due
to different sinkage and trim for the bare hull and the self-propelled hull.

193
 (a) FV-BFM (b) LV-BFM

Figure 6.37: Input velocity distribution on the propeller plane at F n = 0.34.

6.3.2 Velocities on Propeller Plane

The difference in the distribution of input velocity for the two body-force models is
shown in Figure 6.37 on the propeller plane. In case of FV-BFM, the input velocity
is the same as that on the propeller plane, which neglects the local incoming flow.
The input velocity for LV-BFM method varies over the propeller disk as it samples
the incoming velocity on a plane upstream of the propeller disk.

The axial velocities with and without propeller modelling were computed on the
propeller plane using the unsteady solver and compared in Figure 6.38 for F n =
0.34. As shown in the figure, a strong boundary layer effect due to the skeg is visible
in the axial velocity distribution. In the case of the axial velocity with the FV-BFM
method, the average flow is accelerated due to the propeller action. The axial velocity
distribution is almost symmetrical about the center plane. For the axial velocity
distribution based on the LV-BFM method, the effect of local flow can be clearly
observed. A stronger velocity is seen on one side than the other side. The effect of

194
 (a) Without propeller modeling

(b) FV-BFM (c) LV-BFM

Figure 6.38: Axial velocity distribution on the propeller plane at F n = 0.34.

Table 6.7: Comparison of |E%D| for FT ow , sinkage and trim at F n = 0.34

RT (SP ) (N) T (N) Q (N-m) n (rps) FT ow (N) |E%D| Sinkage (m) Trim (deg) |E%D|
k − ω SST 14.6048 15.9 0.438 15.9218 -1.2952 61.99 -0.0156 0.2810 18.76
k − ω SSTCC 14.4751 15.9 0.438 15.9218 -1.3682 59.82 -0.0156 0.2936 15.12
Experimental 15.902 0.464 15.8255 -3.4071 - 0.3459 -

rotational flow on the propeller plane can also be observed. Strong vortices are seen
in the axial velocity contours due to the propeller action.

195
 (a) FV-BFM (b) LV-BFM

Figure 6.39: Effective wake distribution on the propeller plane at F n = 0.34.

Table 6.8: Self-propulsion simulation results using the load varying method by
FV-BFM at F n = 0.34

Numerical Experimental
n (rps) FT ow (N) Trim (deg) n (rps) FT ow (N) Trim (deg)
CASE01 13.8751 4.1924 0.4318 13.7760 2.1189 0.4327
CASE02 14.7391 1.7548 0.3635 14.6566 -0.3246 0.4095
CASE03 15.4899 -0.2972 0.3200 15.3975 -2.4043 0.3431
CASE04 15.9218 -1.3682 0.2936 15.8255 -3.4072 0.3459

Table 6.9: Self-propulsion simulation results using the load varying method by
LV-BFM at F n = 0.34

Numerical Experimental
n (rps) FT ow (N) Trim (deg) n (rps) FT ow (N) Trim (deg)
CASE01 13.7760 4.4160 0.4363 13.7760 2.1189 0.4327
CASE02 14.6566 1.3620 0.4786 14.6566 -0.3246 0.4095
CASE03 15.3975 0.4074 0.3854 15.3975 -2.4043 0.3431
CASE04 15.8255 -0.1583 0.3107 15.8255 -3.4072 0.3459

The effective axial velocity distributions on the propeller plane by two different
body-force models are compared in Figure 6.39. Note that the effective axial velocity,
u, is defined as the difference between the induced axial velocities obtained from open

196
 (a) FV-BFM (b) LV-BFM

Figure 6.40: Propeller suction effect at x/Lpp = 0.10D at F n = 0.34.

water and self-propulsion simulations, and it is normalized using the ship speed. The
effective axial velocity by FV-BFM is in general greater than that by LV-BFM method.
In other words, the flow was more accelerated in the simulations by the FV-BFM
method. The predicted propeller suction effects by the two methods at d/D = 0.10
upstream of the propeller plane are compared in Figure 6.40. The predicted suction
effect is more asymmetric about the ship center plane by the LV-BFM method.

6.3.3 Thrust, Torque, Tow Force and Self-Propulsion Point

In the self-propulsion simulations, the effect of turbulence model on the numerical

solution was investigated using two turbulence models, k − ω SST, and k − ω SSTCC
for F n = 0.34 with the FV-BFM method. The predicted tow force, FT ow , thrust,
T , torque, Q, propeller revolution rate, (n), and trim and sinkage of the model are
presented in Table 6.7. The relative errors, |E%D|, for FT ow and trim, are calculated
using the experimental data.

197
 Table 6.10: Self-propulsion parameters of fishing vessel at F n = 0.34

ns ηR ηH ηo PE (kws) PD (kws) QP C
FV-BFM 3.4134 1.0019 0.9714 0.4954 364.90 836.68 48.22 %
LV-BFM 3.3697 0.9968 1.0040 0.5022 364.90 775.13 50.26 %
Experimental 3.2287 0.9004 1.0786 0.5012 377.05 759.75 48.67 %

To determine the self-propulsion point, additional cases were simulated using different
propeller rates of revolution with the k − ω SSTCC turbulence model and using
FV-BFM and LV-BFM. The numerical results are summarized in Tables 6.8 and 6.9.
The tow force, FT ow , KT , and 10KQ are plotted against the rate of revolution, n,
and compared with the experimental data in Figure 6.41. A larger deviation between
the numerical results and the experimental data can be observed in the results by
both body-force models. At higher n, both FV-BFM and LV-BFM led to a large
discrepancy from the experimental data. The ship self-propulsion point obtained
from the experiment is at n = 13.82 rps, while n = 14.61 rps and n = 14.42 rps were
determined from FV-BFM and LV-BFM, respectively.

6.3.4 Self-Propulsion Results

Results of self-propulsion simulations in full scale for the fishing vessel at 12.28 knots
are presented in Table 6.10 in terms of efficiency and power. At this vessel speed,
the skin friction correction SF C is 2.0873 N. Note that the results were extrapolated
from model scale to full scale. From the table, it can be seen that LV-BFM was able
to predict the self-propulsion results with a greater accuracy than FV-BFM although
it slightly over-predicted the QP C.

198
 8.0 1
KT -Exp 10KQ -Exp FTow -Exp
KT -FV-BFM 10KQ -FV-BFM FTow -FV-BFM
6.0 KT -LV-BFM 10KQ -LV-BFM FTow -LV-BFM
0.8
4.0

0.6

KT, 10KQ
FTow (N)

2.0

0.0
0.4

-2.0
0.2
-4.0

-6.0 0
13 13.5 14 14.5 15 15.5 16

Rate of revolution, n (rps)

Figure 6.41: Comparison of tow forces at F n = 0.34.

6.3.5 Pressure Distribution

The distributions of dynamic pressure coefficient CP on the hull are presented in

Figure 6.42 for the bare hull and the hull with the propeller. The pressure at the end
of skeg for the self-propulsion case is lower than that for the bare hull case due to the
flow acceleration in front of the propeller. As expected, the pressure distribution on
the other part of the hull generally remains the same in the three cases. The effect of
the body-force modelling on the streamlines in the downstream of the propeller plane
is shown in Figure 6.43. Slight differences can be observed.

199
 (a) Bare hull

(b) FV-BFM

(c) LV-BFM

Figure 6.42: Pressure distribution on the hull at F n = 0.34.

200
 (a) FV-BFM

(b) LV-BFM

Figure 6.43: Streamlines around the propeller plane for F n = 0.34

201
6.4 HPC and Computational Time

The numerical simulations for this thesis were carried in skipper and openfoam server.
Both of the server are located at Memorial University. Most of the time, 120 to 132
computing cores were used on 7 nodes, 196 core cluster. Each of nodes were equipped
with two Intel processors at 1.20 GHz. The cluster was running 64-bit Centos 6.5
with Rocks 6.1.1 (Sand Boa) Cluster Management. A single node openfoam server,
56 cores with two Intel(R) Xeon(R) CPU E5-2697 @ 2.60GHz processors was used
for small simulations. The openfoam server was running Red Hat Enterprise Linux 7.
The number of calculation cores were different for different type of simulations. Table
6.11 presents the computational time for different type of cases.

Table 6.11: Computational time in skipper and openfoam server.

Cases Number Number of Computing Time Required

of cells Processors per time step or
(millions) Iteration

Steady Resistance 4.22 72 33 sec 5,000

Steady Resistance 8.43 144 33 sec 3,000
Unsteady 10.31 120 66 sec
Resistance
Wave Simulation 1.32 48 38 sec 500 ∼ 1000
steps per wave
period
Self-propulsion 12.00 100 55 sec
(Actual propeller)
Self-propulsion 8.98 72 78 sec
(body-force)
Self-propulsion in 12.00 96 363 sec 1,000 steps per
wave wave period

202
Chapter 7

Conclusions and Recommendations

7.1 Conclusions

Detailed propeller and inflow modelling is an important step for the accurate
prediction of the effective wake due to propeller-hull interaction, including the
propeller suction effect. The ship’s motion affects the incoming flow to the propeller,
which affects the loading on the propeller, the propeller performance, and the propeller
induced noise and vibration. Consideration of all these effects allows greatly improved
predictions of the actual propeller loading, leading to more efficient propulsion system
designs that consider realistic operating conditions. In the present research, numerical
simulations of propeller-hull interaction including ship motions were performed using
a multiphase viscous flow solver in OpenFOAM with the k −ω SST and the k −ω SST
with curvature correction turbulence models. The propeller-hull interaction was
studied using both the propeller detailed geometry and two body-force models.

Due to the complexity of the final numerical model, the simulation models
were developed and validated in several steps. The steps in the propeller-hull
interaction study include (I) bare hull resistance, wave profile, and flow velocity

203
model development and validation, (II) propeller open water model development
and data validation, (III) numerical wave tank implementation and validation,
(IV) self-propulsion simulation implementation at an even keel condition, (V)
self-propulsion simulation with sinkage and trim and (VI) final full self-propulsion
simulation in waves with ship motions.

In the case of bare hull resistance, wave profile, and flow velocity model validation,
the three ship models, KCS, JBC and a Fishing vessel, were used to validate
the numerical simulations. An anisotropic meshing capability was developed and
implemented in OpenFOAM meshing module snappyHexMesh to overcome the
limitation of anisotropic mesh generation. The computational grids for the bare
hull simulations were generated using the newly developed meshing module, named
navalSnappyHexMesh. Anisotropic meshing in the free surface region greatly reduces
the overall grid numbers by providing the required grid resolution in the free surface
region, which is essential to capture the generated wave elevation, but allowing for
lower grid resolutions in less critical regions. For the KCS and JBC models, a Kelvin
wedge shape grid refinement was used on the free surface with further local grid
refinement in the bow and stern regions. Wall functions were used for near-wall
treatment in the numerical simulations.

Simulations for fixed condition cases were performed using a multiphase steady
solver. The steady solver was also used for dynamic cases, but in the dynamic
cases, the steady solutions were used as initial conditions to reduce the computational
time required for the unsteady simulations. Verification and validation included
uncertainty calculations. The resistance validation and verification show that a
minimum of about seven hundred cells over the ship’s length are necessary for
adequate numerical prediction. This conclusion is also supported by the comparison

204
between predicted and measured wave elevation along the hull surface. Predictions of
wave elevation at different longitudinal cuts are also confirmed with good accuracy.
Axial velocity distributions at different transverse locations, including the propeller
plane, were compared with experimental data using local grid refinement in the
same way as for the hull grid refinement. From comparison with experimental
measurements, it is concluded that numerical simulation using the anisotropic grid
refinement reduces the number of computational grids without compromising the
accuracy of the numerical results.

Propeller open water simulations were validated using data for four propeller models;
KP505, MP687, DTMB5168, and a stock propeller. The KP505, MP687, and
DTMB5168 propellers were used to validate simulation cases that used detailed
geometry modelling. The KP505 and the stock propeller were used for the simulation
cases that used the body-force models. Due to the poor quality mesh generated by
snappyHexMesh for the complex geometries, such as propellers, the propeller meshes
were generated using NUMECA Hexpress. Propeller open water simulations were
carried out using a single-phase viscous solver. Wall functions were also employed for
near-wall treatment of turbulence modelling. The thrust and torque generated by the
detailed propeller geometry simulations showed good agreement with the experimental
data. However, the numerical results for efficiency showed large deviation at higher
advance coefficients. The implemented effective local velocity body-force model agreed
well with the experimental data for torque, thrust and propeller efficiency.

For the numerical wave tank validation, different wave steepness was studied with
different grid resolutions and different time-step sizes. In case of the high wave
steepness, more cells were required per wavelength and wave height. A grid
dependency study showed that 160 cells per wavelength and 20 cells per wave height

205
are required. The required number of time steps per wave period was less for low
wave steepness than for high wave steepness. In the case of high wave steepness,
a minimum of 1000 time steps per wave period were required with a higher-order
divergence scheme for the convective term. If the second-order divergence scheme is
used, almost double the number of time steps is required for the same wave. For
low wave steepness cases, the required number of time steps was almost half that
of the high wave steepness cases. The effect of the turbulence models on numerical
wave generation was also investigated using the realizable k −  and k − ω SST
turbulence models. In all cases 0.10% turbulence intensity was used. It was found
that grid resolutions and time step sizes were adequate for numerical wave generation
with both turbulence models. The numerical wave tank was validated for low wave
steepness for self-propulsion simulation in waves. It was found that less number of
cells per wave length and less number of time steps are adequate for low steepness
waves than high steepness waves.

For the study of propeller-hull interaction, the self-propulsion simulations were carried
out at a fixed even keel condition for the KCS ship model. The optimum grids for
the hull and propeller from the bare hull study and the propeller open water study
were merged for the case of the rotating propeller behind the hull. The simulations
were also carried out using both the nominal average velocity body-force model and
the effective local velocity body-force model. Steady simulations were performed
until the resistance converged. Then the steady simulation results were used as the
initial condition for the unsteady simulations. This approach significantly reduces
the computational time. Rotational speed was applied to the propeller at a fixed
rate of revolution. The numerical results using the detailed propeller geometry agree
well with the available experimental data. Since the flow around the propeller is
highly non-linear and rotational, the k − ω SSTCC turbulence model was used to

206
better model the rotational effect of the flow. The axial velocity distribution at
0.25D downstream of the propeller plane was compared with experimental data. The
derived skin friction coefficient was also compared with that calculated as part of the
experimental data analysis. The effective local velocity body-force model and nominal
average velocity body-force model were also used for self-propulsion simulation. In
terms of propeller force coefficient, the effective local velocity body-force model showed
better performance than the nominal average velocity body-force model.

Self-propulsion simulations were carried out with ship motion in calm water for the
KCS and JBC models using the detailed propeller geometry models. Both body-force
models were used for the self-propulsion simulation of the fishing vessel. During
these numerical simulations, the ship models were limited to heave and pitch motion
only. For self-propulsion simulation with detailed propeller geometry, the developed
dynamic class with sliding mesh was used to allow the combination of ship motions
and a rotating propeller. The overall comparison of the numerical results for the JBC
hull simulation with the experimental data found to show good agreement. In the case
of the KCS model, the total resistance was higher than in the even keel condition. The
propeller thrust was decreased due to the changes in incoming flow. The motion of the
KCS ship model introduced an asymmetric pressure distribution on the hull above the
propeller. All these results were as expected based on the numerical modelling that
allowed these features to be captured. Ships are known to suffer increased resistance
when motions are introduced and propeller thrust is known to decrease. Asymmetric
pressure distributions arise from the propeller rotation and this can lead to vibration
on the hull. The ability to capture the phenomena shows that the developed dynamic
class can be successfully used for self-propulsion simulation with ship motions and a
rotating propeller.

207
A seakeeping simulation was carried out in advance of the self-propulsion simulation
in waves. The optimum numerical settings based on the previous validations
of the resistance, propeller open water, and wave modelling were used in the
seakeeping numerical simulation. Comparison of the simulated motion responses
shows reasonable agreement with the available experimental data. As a concluding
demonstration of the capability of the developed dynamic class for this research, a
self-propulsion simulation with ship forward speed and ship motion was conducted in a
regular wave. Due to a lack of experimental data, it was not possible to compare these
numerical results. The pitch and heave motion plots show that the self-propulsion
simulation can be successfully conducted with the fully assembled numerical model
when the ship experiences relatively large motions. This model provides a realistic
simulation of propeller performance in operating conditions that fully captures all
the effects of a ship operating in a seaway. This allows the variability in flow and
propeller performance inherent in realistic operations to be analyzed as part of a
complete propulsion system design process.

The step-by-step validation of numerical simulation results against the available

experimental data using the developments in this Ph.D. demonstrates the effectiveness
and importance of the continuing development of CFD models and modelling
techniques for flow conditions that ultimately may be too challenging for experimental
measurement. The effective wake and self-propulsion parameters were well predicted
for the ship in calm water and ships with motion using detailed propeller geometry
modelling. The developed dynamic motion class using the sliding mesh method
exhibits good performance without increasing computational time with respect to the
overset grid method. Simplified propeller body-force models that can be used to speed
up computations or for cases where detailed propeller geometry is not yet known were
also developed for self-propulsion simulation and validated with experimental data.

208
7.2 Best practice for propeller-hull interaction.

The following sections provide guidance for the development of numerical models
involving the interaction between propellers and ship hulls. This guidance is intended
to speed development of future models and is based on the experiences of the present
work. These tips are specific to the OpenFOAM meshing and solver software and to
the algorithms developed as part of this study, but it is expected that many items
would also be applicable to other viscous flow solver packages.

7.2.1 Selection of computational domain.

1. The inlet boundary should be at 1.0 ∼ 2.0 ship lengths in front of the model.
For the simulation with incoming wave, the inlet boundary should placed at
least 2.0 ship lengths in front of the model.

2. The outlet boundary should be 3.0 ∼ 4.0 ship lengths behind the model to avoid
wave reflection from the outlet boundary.

3. The side boundaries should be at least one ship length from the hull.

4. The bottom boundary should be one ship length below the hull to reduce the
blockage effect.

5. The top boundary can be set 0.5 ∼ 1.0 ship lengths above the calm water
surface.

6. The bottom boundary in wave simulation should not be placed less than water
depth of the experimental facility. When the experimental facility’s water depth
is unknown, the water depth should be set according to the deep water wave
condition.

209
7.2.2 Mesh generation.

1. The mesh quality and resolution is crucial in numerical simulations with

OpenFOAM.

2. The non-dimensional first layer thickness, y + ≥ 30 allows the use of wall

functions to reduce the total number of cells. To solve the detailed flow features,
such as tip vortex, y + should be less than 5.0.

3. The mesh resolution should be fine enough to capture the geometry accurately.

4. Local refinements are necessary for regions where flow changes. The cell size
should be of the same order as that of the cell near the body.

5. For numerical simulation of ship resistance with free surface, 700 to 1000 cells
per ship length should be distributed.

6. For numerical simulation of propeller-hull interaction, 160 cells per wave

simulation, and 20 cells per wave height is recommended for accurate input
wave simulation for high steepness waves. The required cell number can be
reduced for low steepness waves.

7. The ratio of near wall thickness to the undistorted cell size along the body
surface should not be less than 0.1. Otherwise numerical ventilation will occur
in numerical simulation. The mesh quality should be checked in OpenFOAM
before running the simulation. Attention should be paid to non-orthogonality
and skewness of cells.

210
7.2.3 Selection of propeller modelling

1. The choice of propeller modelling method depends on the required parameters

from the simulation as well as the data availability.

2. Detailed flow features, tip vortex propagation, pressure fluctuation due to the
propeller rotational effect and suction effect, and wake distribution can be
obtained using the detailed propeller model.

3. Because of two different time scales for the ship hull and propeller, the simulation
time is higher for the detailed propeller method than that for the body-force
method.

4. For fast self-propulsion simulation, the body-force method can be used, which
requires propeller open water curves.

7.2.4 Solver selection

1. The steady solver is adequate to predict the force in steady flow with reasonable
accuracy.

2. If the wave making resistance in a low speed region (0.1 < F n < 0.2), the
single-phase steady simulation would be a good choice. If not, unsteady
simulation with free surface is recommended.

3. The steady simulation is recommended for the initial condition of unsteady

simulation, which can reduce the computational time significantly.

7.2.5 Boundary conditions set-up

1. Velocity boundary conditions should be imposed on the inlet boundary.

211
 2. Neumann boundary conditions can be imposed on the outlet boundary since it
is far from the body.

3. For deep water simulations, side boundaries and the bottom boundary can be
set as symmetry planes. The top boundary can be set as pressure inletOutlet
boundary.

4. The body needs to be set as a no-slip wall boundary.

5. Depending on y + , the wall functions for turbulence parameters need to be set.

If y + ≤ 1.0, fixed value condition is recommended for turbulence parameters. In
case of y + ≥ 30, a fully turbulent wall function is recommended for turbulence
parameters.

7.2.6 Selection of turbulence model

1. Reynolds stress models predict the vortex more accurately with a greater
computing cost.

2. Turbulence models with curvature correction are recommended.

3. realizable k −  turbulence model is recommended for wave simulation since

over-production of the turbulence viscosity in the k − ω SST turbulence model
damps out the wave height as the simulation progresses.

7.2.7 Selection of numerical schemes

1. First-order schemes can be used to start the numerical simulation.

2. To improve accuracy, numerical schemes need to be changed to second order

after initial convergence. The change in numerical scheme can be done manually
or automatically using timeActivatedFileUpdate function.

212
7.2.8 Time step size selection

1. In OpenFOAM, the PIMPLE algorithm helps to increase the time step size, i.e.,
a higher Courant number. In this case, the number of outer iterations should
be increased. For better accuracy, the Courant number for free surface should
be less than 5.0.

2. For self-propulsion simulation, the time scales for the hull and propeller are
different. The Courant number should be less than 1.0. A greater Courant
number could lead to incorrect results.

3. With the sliding mesh technique, 360 to 720 time steps per revolution of the
propeller is recommended.

7.2.9 Setting convergence criteria

1. For convergence of numerical simulation, the physical quantity (such as

resistance for ship) should be monitored along with the residuals.

2. The initial residual for pressure should be set to 10−2 or less.

3. For other parameters, such as volume fraction, k and ω, the initial residual
should be set to 10−6 .

4. For resistance convergence, 0.01% in error is suggested as the convergence

criterion.

7.3 Recommendations for Future Work

Based on the experience of the present research, the following should be considered for
future work. This work can improve some aspects of the methods and codes developed

213
as part of this research or in some cases build additional capability using the present
work as a base.

1. The prism layer generation in snappyHexMesh is inadequate, particularly when

the cell has a high aspect ratio. The prism layer generation module can be
improved.

2. Only constant rotational speed is implemented for propeller rotation in the

present study. Time dependent motion can be developed for propeller rotation
and can also be applied for rudder motion.

3. The present simulation is limited to regular waves. Irregular wave trains could
be incorporated into the simulation to provide more realistic motion simulations.

4. Because of the different time scale for ship hull and propeller simulations, the
detailed propeller numerical simulation requires longer simulation time. When
using the time scale for the hull, the propeller force and moment predictions are
incorrect. A numerical method might be developed to achieve propeller force
and moment convergence at a higher time scale.

5. Ultimately ship maneuvering simulation can be implemented, based on the

present code, by adding sway and yaw motions which would allow the effects
of horizontal plane ship motions on the propeller inflow and loading to be
simulated.

214
Bibliography

Ali, M. A., Peng, H., and Qiu, W. 2017a. Benchmark Studies of Wave Run-Up and
Forces on a Truncated Square Cylinder. International Conference on Offshore
Mechanics and Arctic Engineering, Volume 1: Offshore Technology.

Ali, M. A., Peng, H., Qiu, W., and Bensow, R. 2017b. Prediction of Propeller Tip
Vortex Using OpenFOAM. International Conference on Offshore Mechanics and
Arctic Engineering, Volume 7A: Ocean Engineering.

Ali, M. A., Peng, H., and Qiu, W. 2018. Resistance computation of an optimized
fishing vessel. In: 26th annual conference of the computational fluid dynamics
society of canada.

Ali, M. A., Peng, H., and Qiu, W. 2019a. Resistance prediction of two fishing
vessel models based on RANS solutions. Physics and Chemistry of the Earth,
Parts A/B/C, 113, 115–122. Computational Environmental Flow and Transport
Phenomena.

Ali, M. A., Peng, H., and Qiu, W. 2019b. Simulations of self-propulsion model tests of
a fishing vessel using a body-force method coupled with a RANS solver. In: The
Joint Canadian Society for Mechanical Engineering and CFD Society of Canada
International Congress 2019.

215
ASME. 2009. Standard for verification and validation in computational fluid
dynamics and heat transfer-ASME V&V 20-2009. The American Society of
Mechanical Engineers (ASME).

Banks, J., Phillips, A. B., and Turnock, S. R. 2010. Free surface CFD prediction of
components of ship resistance for KCS. 13th Numerical Towing Tank Symposium,
Germany.

Bensow, R. E. 2013. Simulation of unsteady propeller blade loads using OpenFOAM.

16th Numerical Towing Tank Symposium, NuTTS’13.

Bertram, V. 2012. Practical Ship Hydrodynamics. Second edn. Oxford:

Butterworth-Heinemann.

Brackbill, J. U, Kothe, D. B., and Zemach, C. 1992. A continuum method for modeling
surface tension. Journal of Computational Physics, 100(2), 335–354.

Carrica, P. M., Mofidi, A., Eloot, K., and Delefortrie, G. 2016. Direct simulation
and experimental study of zigzag maneuver of KCS in shallow water. Ocean
Engineering, 112, 117–133.

Castro, A. M., Carrica, P. M., and Stern, F. 2011. Full scale self-propulsion
computations using discretized propeller for the KRISO container ship KCS.
Computers & Fluids, 51(1), 35–47.

Chandar, D. D, and Gopalan, H. 2016. Comparative analysis of the Arbitrary Mesh

Interface (AMI) and overset methods for dynamic body motions in OpenFOAM.
46th AIAA Fluid Dynamics Conference.

Chesnakas, C. J., and Jessup, S. D. 1998. Propeller tip vortex measurements using
3-component LDV. Proc 22nd Symp on Naval Hydrodynamics, 156–170.

216
Chong, M. S., Perry, A. E., and Cantwell, B. J. 1990. A General Classification
of Three-Dimensional Flow Fields. Physics of Fluids A: Fluid Dynamics, 2(5),
765–777.

Dhinesh, G, Murali, K, and Subramanian, A. V. 2010. Estimation of hull-propeller

interaction of a self-propelling model hull using RANSE solver. Ships and Offshore
Structures, 5(2), 125–139.

Doormaal, J. P., and Raithby, G. D. 1984. Enhancements of the SIMPLE method

for predicting incompressible fluid flows. Numerical Heat Transfer Part A -
Applications, 7(2), 147–163.

Eça, L., and Hoekstra, M. 2006. Discretization Uncertainty Estimation based on a

Least Squares version of the Grid Convergence Index. In: 2nd Workshop on CFD
Uncertainty Analysis, Lisbon. International Conference on Offshore Mechanics and
Arctic Engineering.

Eça, L., and Hoekstra, M. 2014. A procedure for the estimation of the numerical
uncertainty of CFD calculations based on grid refinement studies. Journal of
Computational Physics, 262, 104–130.

Eça, L., Vaz, G., and Hoekstra, M. 2010. Code Verification, Solution Verification and
Validation in RANS Solvers. In: 29th International Conference on Ocean, Offshore
and Arctic Engineering: Volume 6.

Farrell, P. E., and Maddison, J. R. 2011. Conservative interpolation between volume

meshes by local Galerkin projection. Computer Methods in Applied Mechanics
and Engineering, 200(1), 89–100.

Felli, M, and Felice, F. D. 2004. Analysis of the propeller-hull interaction by LDV

phase sampling techniques. Journal of Visualization, 7(1), 77–84.

217
Ferziger, J. H., and Peric, M. 2012. Computational methods for fluid dynamics.
Springer.

Friis, D., Bass, D., Gradner, A., and Edward, L. S. 2007. Design of two multi-species
fishing vessel. Phase 2. Industrial Outreach Group, Faculty of Engineering and
Applied Science. MUN.

Fujisawa, J., Ukon, Y., Kume, K., and Takeshi, H. 2000. Local velocity field
measurements around the KCS model in the SRI 400 m towing tank. Tech. rept.
Ship Performance Division Report No. 00-003-02, The Ship Research Institute of
Japan.

G2010. 2010. Gothenburg 2010 A Workshop on CFD in Ship Hydrodynamics,

Gothenburg, Sweden.

Gaggero, S., Gaggero, T., Rizzuto, E., Tani, G., Villa, D., and Viviani, M. 2016. Ship
propeller side effects: Pressure pulses and radiated noise. Noise Mapping, 3(01).

Goldstein, S. 1929. On the vortex theory of screw propellers. Proceedings of the Royal
Society of London. Series A, Containing Papers of a Mathematical and Physical
Character, 123(792), 440–465.

Higuera, P. 2020. Enhancing active wave absorption in RANS models. Applied Ocean
Research, 94, 102000.

Higuera, P, Lara, J. V., and Losada, I. J. 2013. Realistic wave generation and active
wave absorption for NavierStokes models: Application to OpenFOAM. Coastal
Engineering, 71, 102–118.

Hirt, C. W., and Nichols, B. D. 1981. Volume of fluid (VOF) method for the dynamics
of free boundaries. Journal of Computational Physics, 39(1), 201–225.

218
Holmen, V. 2012. Methods for vortex identification. Masters thesis, Lund University.

Hough, G. R., and Ordway, D. E. 1964. The generalized actuator disk. Tech. rept.
Therm Advanced Research, Inc.

Hsiao, C., and Pauley, L. L. 1999. Numerical computation of the tip vortex flow
generated by a marine propeller. ASME Journal of Fluids Engineering, 121(3),
638–645.

Hylarides, S. 1978. Some hydrodynamic consideration of propeller induced ship

vibrations. In: The Ship Vibration Symposium. Society of Naval Architects
and Marine Engineers.

Inukai, Y., and Ochi, F. 2009. A study on the characteristics of self-propulsion factors
for a ship equipped with contra-rotating propeller. First International Symposium
on Marine Propulsors, SMP09.

Issa, R. I. 1986. Solution of the implicitly discretised fluid flow equations by

operator-splitting. Journal of Computational Physics, 62(1), 40–65.

ITTC. 1978. ITTC recommended procedures and guidelines. ITTC performance

prediction method, 7.5-02-03-01.4.

ITTC. 2008. ITTC recommended procedures and guidelines. ITTC testing and
extrapolation methods propulsion, performance propulsion test, 7.5-02-03-01.1.

ITTC. 2011. ITTC recommended procedures and guidelines. Practical Guidelines for
Ship CFD Applications, 7.5-03-02-03.

Jacobsen, N., Fuhrman, D. R., and Fredsøe, J. 2012. A wave generation toolbox for
the opensource CFD library: OpenFOAM. International Journal for Numerical
Methods in Fluids.

219
Jasak, H. 1996. Error analysis and estimation for the finite volume method with
applications to fluid flows. Ph.D. thesis, Imperial College of Science, Technology
and Medicine, Department of Mechanical Engineering.

Juretic, F. 2017. cfMesh: Advanced Meshing Tool.

Kawamura, T., Miyata, H., and Mashimo, K. 1997. Numerical simulation of the flow
about self-propelling tanker models. Journal of Marine Science and Technology, 2,
245–256.

Kerwin, J. E., and Lee, C. S. 1978. Prediction of steady and unsteady marine propeller
performance by numerical lifting-surface theory. The Society of Naval Architects
and Marine Engineers, Transaction., 86:21853.

Kim, S., Su, Y., and Kinnas, S. A. 2018. Prediction of the Propeller-induced Hull
Pressure Fluctuation via a Potential-based Method: Study of the Influence of
Cavitation and Different Wake Alignment Schemes. In: Proceedings of the 10th
International Symposium on Cavitation (CAV2018). ASME Press.

Kim, W. J., Van, S. H., and Kim, D. H. 2001. Measurement of flows around modern
commercial ship models. Experiments in Fluids, 31(5), 567–578.

Larsson, L., Stern, F., and Visonneau, M. 2014. Numerical Ship Hydrodynamics An
assessment of the Gothenburg 2010 Workshop. Springer, Dordrecht.

Launder, B. E., Reece, G. J., and Rodi, W. 1975. Progress in the development of a
Reynolds-stress turbulence closure. Journal of Fluid Mechanics, 68(3), 537566.

Lee, S., and Chen, H. 2005. The influence of propeller/hull interaction on propeller
induced cavitating pressure. In: The Fifteenth International Offshore and Polar

220
 Engineering Conference, International Society of Offshore and Polar Engineers
ISOPE.

Lewis, E. V. (ed). 1988. Principles of Naval Architecture: Volume II: Resistance,

propulsion and vibration. Society of Naval Architects and Marine Engineers.

Mayer, S., Garapon, A., and Sørensen, L. S. 1998. A fractional step method
for unsteady freesurface flow with applications to nonlinear wave dynamics.
International Journal for Numerical Methods in Fluids, 28, 293–315.

Menter, F. R., Kuntz, M., and Langtry, R. B. 2003. Ten years of industrial experience
with the SST turbulence model. Turbulence, Heat and Mass Transfer, 625–632.

Muzaferija, S., and Peric, M. 1997. Computation of free-surface flows using the
finite-volume method and moving grids. Numerical Heat Transfer, Part B:
Fundamentals, 32(4), 369–384.

Nagamatsu, T. Sasajima, T. 1975. Effect of propeller suction on wake. Journal of

the Society of Naval Architects of Japan, 1975(137), 58–63.

Nakatake, K. 1989. On ship hullpropellerrudder interactions (in Japanese). In

Proceedings of the 3rd JSPC symposium on flows and forces of ships, 231259.

Oliveira, P. J., and Issa, R. I. 2001. An improved PISO algorithm for the computation
of buoyancy-driven flows. Numerical Heat Transfer, Part B: Fundamentals, 40(12),
473–493.

OpenFOAM. 2019 (December 23). OpenFOAM Users Guide-1912. OpenFOAM.

Paik, B. G, Lee, C. M., and Lee, S. J. 2004. PIV analysis of flow around a container
ship model with a rotating propller. Experiments in Fluids, 36, 833–846.

221
Paik, B. G., Kim, G. D., Lee, J. Y., and Lee, S. J. 2007. Visualization of the inflow
ahead of a rotating propeller attached to a container ship model. Journal of
Visualization, 10(1), 47–55.

Paik, K-J, Park, H-G, and Seo, J. 2013. RANS simulation of cavitation and
hull pressure fluctuation for marine propeller operating behind-hull condition.
International Journal of Naval Architecture and Ocean Engineering, 5(4), 502–512.

Patankar, S. V. 1980. Numerical heat transfer and fluid flow. Hemisphere Pub.

Patankar, S. V., and Spalding, D. B. 1972. A calculation procedure for heat, mass and
momentum transfer in three-dimensional parabolic flows. International Journal of
Heat and Mass Transfer, 15(10), 1787–1806.

Pecoraro, A., Felice, F. D., Felli, M., Salvatore, F., and Viviani, M. 2013.
Propeller-hull interaction in a single-screw vessel. Third International Symposium
on Marine Propulsors, 10.

Peng, H., Qiu, W., and Ni, S. 2013. Effect of turbulence models on RANS computation
of propeller vortex flow. Ocean Engineering, 72, 304–317.

Pereira, F. S., Eça, L., and Vaz, G. 2017. Verification and validation exercises for the
flow around the KVLCC2 tanker at model and full-scale Reynolds numbers. Ocean
Engineering, 129, 133–148.

Qiu, W., Peng, H., Ni, S., Liu, L., and Mintu, S. 2013. RANS computation of propeller
tip vortex. International Journal of Offshore and Polar Engineering, 23(1), 73–79.

Roache, P. J. 1998. Verification and validation in computational science and

engineering. Albuquerque, NM: Hermosa.

222
Rusche, H. 2002. Computational Fluid Dynamics of dispersed two-phase flows at
high phase fractions. Ph.D. thesis, Imperial College of Science, Technology and
Medicine.

Shen, Z., Wan, D., and Carrica, P. M. 2015. Dynamic overset grids in OpenFOAM
with application to KCS self-propulsion and maneuvering. Ocean Engineering,
108, 287–306.

Sigmund, S., and el Moctar, O. 2017. Numerical and experimental investigation of

propulsion in waves. Ocean Engineering, 144, 35–49.

Smirnov, P. E., and Menter, F. R. 2009. Sensitization of the SST turbulence model
to rotation and curvature by applying the Spalart - Shur correction term. Journal
of Turbomachinery, 4(131).

Stansberg, C. T., Baarholm, R., Kristiansen, T., Hansen, E. W. M., and Rortveit,
G. 2005. Extreme wave amplification and impact loads on offshore structures. In:
Offshore technology conference.

Stern, F., Kim, H. T., and Patel, V. C. 1988. A viscous flow approach to the
computation of propeller-hull interaction. Journal of Ship Research, 32, 246–262.

Stern, F., V. Wilson, R. V., Coleman, H. W., and Paterson, E. G. 1999. Verification
and validation of CFD simulations. Tech. rept. Iowa Institute of Hydraulic Research,
The University of Iowa.

T2005. 2005. CFD Workshop 2005 Tokyo, Japan.

T2015. 2015. Tokyo 2015 A Workshop on CFD in Ship Hydrodynamics, Japan.

Tocu, A. M., and Amoraritei, M. 2008. Numerical study of the flow field around a
ship hull including propeller effects. Journal of Maritime Research, 5(3), 67–78.

223
Tryggvason, G., Bunner, B., Esmaeeli, A., Juric, D., Al-Rawahi, N., Tauber, W., Han,
J., Nas, S., and Jan, Y.-J. 2001. A front-tracking method for the computations of
multiphase flow. Journal of Computational Physics, 169(2), 708–759.

Versteeg, H. K., and Malalasekera, W. 2007. An introduction to computational fluid

dynamics: the finite volume method. John Wiley.

Villiers, E. 2006. The potential of large eddy simulation for the modeling of wall
bounded flows. Ph.D. thesis, Imperial College of Science, Technology and Medicine,
Department of Mechanical Engineering.

Watanabe, O., Masuko, A., and Shiorose, Y. 1994. Numerical simulation of a

wave-viscous flow about a ship including propeller effects. Proceedings 20th
Symposium on Naval Hydrodynamics, 946–958.

Weller, H. 2008. A new approach to VOF-based interface capturing methods for

incompressible and compressible flow. Tech. rept. OpenCFD Ltd.

Win, Y. N., Tokgoz, E., Wu, P. C., Stern, F., and Toda, Y. 2013. Computation of
propeller-hull interaction using simple body-force distribution model around Series
60 CB = 0.6. Journal of the Japan Society of Naval Architects and Ocean Engineers,
18, 17–27.

Winden, B., Badoe, C., Turnock, S., Phillips, A., and Hudson, D. 2013. Self
propulsion in waves using a coupled RANS-BEMt model and active RPM control.
In: Proceedings of the 16th Numerical Towing Tank Symposium.

Xing, T., and Stern, F. 2010. Factors of safety for Richardson extrapolation. Journal
of Fluids Engineering, 132(9), 691–696.

224
Zhang, B., Lou, J., Kang, C-W., Wilson, A., Lundberg, J., Svennberg, U., and
Bensow, R. 2014. CFD modeling of propeller tip vortex over large distances. In:
International Journal of Offshore and Polar Engineering. ISOPE.

Zhang, D. H., Broberg, L., and Larsson, L. 1992. A method for computing stern flows
with an operating propeller. Transaction of Royal Institute of Naval Architect,
134, 245–259.

Zhang, Z. 2010. Verification and validation for RANS simulation of KCS container
ship without/with propeller. Journal of Hydrodynamics, Ser. B, 22(5, Supplement
1), 932–939.

Zou, L., and Larsson, L. 2014. CFD verification and validation in practice - A study
based on resistance submission to the Gothenburg 2010 Workshop on numerical
ship hydrodynamics. In: 30th Symposium on Naval Hydrodynamics.

225
This page intentionally left blank.
Appendices

227
Appendix A

Propeller Tip Vortex Simulation

A.1 Geometry and Test Conditions

A five blade fixed-pitch propeller model, DTMB5168 was used in the numerical
simulation for the propeller tip vortex study. The principal particulars of the
DTMB5168 propeller model are given in Table A.1. For geometric simplification,
an infinitely long constant radius cylinder was considered as the hub on which the
blades were mounted, and the root fillet was ignored.

Table A.1: Propeller particulars

Parameters Variable Value Unit

Diameter D 0.403 m
Inflow Velocity U∞ 10.70 m/s
Chord Length at 0.7R C 0.175 m
Advance Coefficient J 1.1
Rotation Speed N 24.167 rps
Combined Velocity at 0.7R V∞ 23.93 m/s
Reynolds Number Re 4.2×106

228
A.2 Computational Domain and Grid generation

A spiral-like computational domain was adopted similar to the one in the previous
work by Peng et al. (2013) so that the grid concentration can be achieved along the tip
vortex core to minimize the flow across the periodic boundaries. The computational
domain was created using the blade to blade passage of the suction side and the
pressure side of the blades with two periodic boundaries as shown in Figure A.1.
According to the domain-dependency study by Qiu et al. (2013), it was found that
the outlet boundary from D to 1.5D gave almost the same results. In this study,
the computational domain was created with the inlet boundary at 0.5D upstream,
the outlet boundary at 1.0D downstream, and the outer boundary at 1.0D from the
propeller axis. A small modification, extending the inlet and outlet boundary in the
axial direction of the computational domain was made to avoid the problematic cells
during the grid generation process.

In previous work, the grid was generated based on the work of Hsiao and Pauley
(1999). It was found that the grid with lower y + value and/or a large number of cells
within the tip vortex core led to convergence problems in ANSYS CFX due to the
large Jacobian values for the generated grid. Thus the smallest first grid spacing was
restricted to y + = 7.5. The local cell refinement was also at a constant radius from the
propeller axis rather than along the tip vortex trajectory, which was not sufficiently
accurate to capture the tip vortex in the far field.

In this work, an unstructured hexahedral mesh was generated using NUMECA

Hexpress to overcome the shortcoming in the previous grid generation method. At the
first step, the computational domain was divided into a number of small blocks and
the refinement level on the blade surface was obtained by dividing the small blocks

229
 (a) Blade passage (b) Computational domain.

Figure A.1: Geometry and computation domain for grid generation

(a) GRID-A (b) GRID-B (c) GRID-C/GRID-D

Figure A.2: Type of refinement.

to a certain level. In the second step, the local refinement was applied at the desired
location depending on the type of refinement method. Then the boundary vertices
were snapped on the boundaries. Finally, the prism layer cells were inserted in the
viscous region and the grid was optimized and smoothed. Following the above steps,
it is possible to keep the cell quality in terms of the skewness and non-orthogonality
in the acceptable range for the solver with a small y + .

It is well known that the numerical predictions of the tip vortex are largely dependent
on the grid concentration in the vortex core (Zhang et al., 2014). To obtain the
grid concentration along the vortex core, three types of grid refinement methods were
employed in the grid generation process. An annular cylindrical refinement region

230
 (a) Level 3 (b) Level 4 (c) Level 5

Figure A.3: Refinement level.

was used for GRID-A. The cells within the specified region were refined to a specified
refinement level. This refinement type was used to track the initial predicted tip vortex
trajectory. In the second approach, the cells around the tip vortex were refined to
improve the numerical prediction of the vortex strength which is shown in GRID-B.
In the third approach, the cells around the experimental tip vortex trajectory were
refined for the better prediction of tip vortex in GRID-C/ GRID-D. Figure A.2 shows
the three refinement types. The base cell size was set 0.01 m in grid generation. Then
the cells were refined in three directions repeatedly depending on the refinement level.
In this study, three refinement levels were used which are shown in Figure A.3. The
difference between GRID-C and GRID-D is the y + values on the blades surfaces. For
GRID-C, the average y + = 5.0 and for GRID-D the average y + on the blade is 2.0.
The grid information is summarized Table A.2. In the first-grid spacing dependency
study, two y + values, 2.0 and 5.0, were used in the grid generation.

A.3 Numerical Results Validation

Validation studies were carried out for the DTMB5168 propeller model at J = 1.1.
Since the two periodic boundaries are non-conformal, the rotational cyclic Arbitrary

231
 Table A.2: Computational grids for DTM5168 propeller.

Grid Number of cells (In million) y+

GRID-A 9.47 5
GRID-B 15.56 5
GRID-C1 10.74 5
GRID-C2 13.67 5
GRID-C3 15.84 5
GRID-D1 18.38 2
GRID-D2 21.31 2
GRID-D3 27.19 2

Mesh Interface (cyclicAMI) was applied. The boundary conditions are summarized
in Table A.3.
The computed KT and KQ based on different grids are presented and compared with
the experimental results by Chesnakas and Jessup (1998) and numerical results by
Peng et al. (2013) in Table A.4. Figure A.4 shows the comparison of computed KT
and KQ using GRID-C3, and GRID-D3 in comparison with the experimental data. It
can be seen that the predictions using the two grids are almost the same. Relatively
greater discrepancies are shown at the lower advance coefficients. From Table A.4, it
can also be seen that the accuracy of the predictions of KT and KQ depends on the
first-grid spacing. Accuracy increases with the decrease of the first-grid spacing and
the increase of the number of grids on the blade surface. However, a greater relative

Table A.3: Boundary conditions for tip vortex simulations.

Parameters Inlet Outlet Outer Hub Blade

U fixedValue inletOutlet slip fixedValue fixedValue
P zeroGradient fixedValue zeroGradient zeroGradient zeroGradient
nut zeroGradient zeroGradient zeroGradient nutUSpaldingWallFunction nutUSpaldingWallFunction
K fixedValue inletOutlet slip fixedValue fixedValue
Omega fixedValue inletOutlet slip fixedValue fixedValue

232
 Table A.4: Comparison of KT and 10KQ

Grid KT E%D 10KQ E%D

GRID-A 0.312 -0.319 0.833 6.386
GRID-B 0.310 -0.958 0.824 5.236
GRID-C1 0.317 1.278 0.841 7.407
GRID-C2 0.317 1.278 0.841 7.407
GRID-C3 0.312 -0.319 0.826 5.492
GRID-D1 0.316 0.958 0.834 6.513
GRID-D2 0.316 0.958 0.834 6.513
GRID-D3 0.316 0.958 0.834 6.513
Peng et al. (2013) 0.312 0.320 0.805 2.81
Experimental 0.313 0.783

1.4
Experimental
GRID-C3
GRID-D3
1.2

1
KT ,10 KQ

0.8

0.6

0.4

0.2

0
0.6 0.8 1 1.2 1.4 1.6
J
Figure A.4: Comparison of KT and 10KQ with experimental data (Chesnakas and
Jessup (1998))

error was found in comparison with the results of the previous study (Peng et al.,
2013).

Figures A.5, A.6 and A.7 present the computed velocity contours for Vx , Vt and Vr ,
respectively at x/R = 0.2386 downstream. The numerical results in general show a

233
 1

Vx
1.5700
1.4864
0.5 1.4029
1.3193
1.2357
1.1521
1.0686
0.9850
0

z
0.9014
0.8179
0.7343
0.6507
0.5671
0.4836
-0.5 0.4000

-1
-1 -0.5 0 0.5 1
y

(a) Experimental. (b) Numerical result by Qiu et al. (2013)

Y Y
1.00 0.50 0.00 -0.50 -1.00 Vx 1.00 0.50 0.00 -0.50 -1.00 Vx
1.00 1.00 1.5700 1.00 1.00 1.5700

1.4864 1.4864

1.4029 1.4029

1.3193 1.3193
0.50 0.50 0.50 0.50
1.2357 1.2357

1.1521 1.1521

1.0686 1.0686

Z 0.00 0.00 Z 0.9850 Z 0.00 0.00 Z 0.9850

0.9014 0.9014

0.8179 0.8179

0.7343 0.7343
-0.50 -0.50 -0.50 -0.50
0.6507 0.6507

0.5671 0.5671

0.4836 0.4836

-1.00 -1.00 0.4000 -1.00 -1.00 0.4000

1.00 0.50 0.00 -0.50 -1.00 1.00 0.50 0.00 -0.50 -1.00
Y Y

(c) GRID-C3 (d) GRID-D3

Figure A.5: Axial velocity Vx contour at x/R = 0.2386

good agreement with the experimental data. To measure the accuracy of the numerical
prediction quantitatively, the velocity components, Vx , Vt and Vr across the tip vortex
core, based on grids with y + = 2 and y + = 5 are compared with the experimental
data and the previously computed results by Peng et al. (2013) in Figures A.10 and
Figure A.11, respectively. It should be noted that the vortex core was defined as
the point on the plane with the minimum pressure and the corresponding value θ =
0. The grid sensitivity is analyzed using root square mean (RMS) error between the
experimental and numerical results. Table A.5 and A.6 show the RMS error for Vx , Vt
and Vr for x/R = 0.1756 and x/R = 0.2386. The numerical predictions of the velocity

234
 1

Vt
2.7502
2.5884
0.5 2.4267
2.2649
2.1032
1.9414
1.7797
1.6180
0

z
1.4562
1.2945
1.1327
0.9709
0.8092
0.6474
-0.5 0.4857

-1
-1 -0.5 0 0.5 1
y

(a) Experimental. (b) Numerical results by Qiu et al. (2013)

Y Vt Y
1.00 0.50 0.00 -0.50 -1.00 1.00 0.50 0.00 -0.50 -1.00 Vt
1.00 1.00 2.7502 1.00 1.00 2.7502
2.5880 2.5884
2.4263 2.4267
2.2645 2.2649
0.50 0.50 0.50 0.50
2.1028 2.1032
1.9410 1.9415
1.7792 1.7797

Z 0.00 0.00 Z 1.6175 Z 0.00 0.00 Z 1.6179

1.4558 1.4562
1.2940 1.2945
1.1322 1.1327
-0.50 -0.50 -0.50 -0.50
0.9705 0.9709

0.8088 0.8092

0.6470 0.6475

-1.00 -1.00 0.4857 -1.00 -1.00 0.4857

1.00 0.50 0.00 -0.50 -1.00 1.00 0.50 0.00 -0.50 -1.00
Y Y

(c) GRID-C3 (d) GRID-D3

Figure A.6: Tangential velocity Vt contour at x/R = 0.2386

components are improved by using the third grid refinement approach for Vx and Vr .
For the tangential velocity components Vt the present approach could not improve
the results. However, the numerical prediction of the peak value at the vortex core
was improved. The prediction better agrees with the experimental data when the
grid refinement around the vortex core is increased and the first-grid spacing, y + , is
decreased. Similar behaviour is also observed at x/R = 0.1756 downstream in Figure
A.8 and Figure A.9. The results were improved for Vx , and Vr only. Note that the
presented velocity components are normalized using the free stream velocity U∞ .

235
 1

vr
0.4126
0.3332
0.5 0.2539
0.1745
0.0951
0.0157
-0.0636
-0.1430
0

z
-0.2224
-0.3017
-0.3811
-0.4605
-0.5399
-0.6192
-0.5 -0.6986

-1
-1 -0.5 0 0.5 1
y

(a) Experimental
(b) Numerical results by Qiu et al. (2013)
Y Vr Y
1.00 0.50 0.00 -0.50 -1.00 1.00 0.50 0.00 -0.50 -1.00 Vr
1.00 1.00 0.4126 1.00 1.00 0.4126
0.3332 0.3332
0.2539 0.2539
0.1745 0.1745
0.50 0.50 0.50 0.50
0.0951 0.0951
0.0157 0.0157
-0.0636 -0.0636

Z 0.00 0.00 Z -0.1430 Z 0.00 0.00 Z -0.1430

-0.2224 -0.2224
-0.3017 -0.3017
-0.3811 -0.3811
-0.50 -0.50 -0.50 -0.50
-0.4605 -0.4605

-0.5399 -0.5399

-0.6192 -0.6192

-1.00 -1.00 -0.6986 -1.00 -1.00 -0.6986

1.00 0.50 0.00 -0.50 -1.00 1.00 0.50 0.00 -0.50 -1.00
Y Y

(c) GRID-C3 (d) GRID-D3

Figure A.7: Radial velocity Vr contour at x/R = 0.2386

The velocity components across the tip vortex core based on different types of grids
are compared in Figures A.12 and A.13. From these figures, it can be found that the
numerical results based on Grid Type III with refinement level 5 show good agreement
with the experimental data.

The pressure coefficients on the suction side and the pressure side of the blade
are compared with the numerical results obtained by Hsiao and Pauley (1999) and
previously computed results using ANSYS CFX in Figures A.14 and A.15. There are
discrepancies in pressure distribution on the suction side in comparison with those by

236
 2.00
Experimental
Peng et al. (SST)
Y+2-GRID-D1
Y+2-GRID-D2
Y+2-GRID-D3

1.50

Vx
1.00

0.50

0.00
-0.40 -0.20 0.00 0.20 0.40
θ
3.50
Experimental
Peng et al. (SST)
Y+2-GRID-D1
Y+2-GRID-D2
Y+2-GRID-D3

3.00
Vt

2.50

2.00

1.50
-0.40 -0.20 0.00 0.20 0.40
θ
1.50
Experimental
Peng et al. (SST)
Y+2-GRID-D1
Y+2-GRID-D2
Y+2-GRID-D3
1.00

0.50
Vr

0.00

-0.50

-1.00
-0.40 -0.20 0.00 0.20 0.40
θ

Figure A.8: Vx , Vt and Vr across vortex core at x/R = 0.1756

237
 2
Experimental
Peng et al. (SST)
Y+5-GRID-C1
Y+5-GRID-C2
Y+5-GRID-C3

1.5

Vx
1

0.5

0
-0.4 -0.2 0 0.2 0.4
θ
3.5
Experimental
Peng et al. (SST)
Y+5-GRID-C1
Y+5-GRID-C2
Y+5-GRID-C3

3
Vt

2.5

1.5
-0.4 -0.2 0 0.2 0.4
θ
1.5
Experimental
Peng et al. (SST)
Y+5-GRID-C1
Y+5-GRID-C2
Y+5-GRID-C3
1

0.5
Vr

-0.5

-1
-0.4 -0.2 0 0.2 0.4
θ

Figure A.9: Vx , Vt and Vr across vortex core at x/R = 0.1756

238
 2.00
Experimental
Peng et al. (SST)
Y+2-GRID-D1
Y+2-GRID-D2
Y+2-GRID-D3

1.50

Vx
1.00

0.50

0.00
-0.40 -0.20 0.00 0.20 0.40
θ
3.50
Experimental
Peng et al. (SST)
Y+2-GRID-D1
Y+2-GRID-D2
Y+2-GRID-D3

3.00
Vt

2.50

2.00

1.50
-0.40 -0.20 0.00 0.20 0.40
θ
1.50
Experimental
Peng et al. (SST)
Y+2-GRID-D1
Y+2-GRID-D2
Y+2-GRID-D3
1.00

0.50
Vr

0.00

-0.50

-1.00
-0.40 -0.20 0.00 0.20 0.40
θ

Figure A.10: Vx , Vt and Vr across vortex core at x/R = 0.2386

239
 2
Experimental
Peng et al. (SST)
Y+5-GRID-C1
Y+5-GRID-C2
Y+5-GRID-C3

1.5

Vx
1

0.5

0
-0.4 -0.2 0 0.2 0.4
θ
3.5
Experimental
Peng et al. (SST)
Y+5-GRID-C1
Y+5-GRID-C2
Y+5-GRID-C3

3
Vt

2.5

1.5
-0.4 -0.2 0 0.2 0.4
θ
1.5
Experimental
Peng et al. (SST)
Y+5-GRID-C1
Y+5-GRID-C2
Y+5-GRID-C3
1

0.5
Vr

-0.5

-1
-0.4 -0.2 0 0.2 0.4
θ

Figure A.11: Vx , Vt and Vr across vortex core at x/R = 0.2386

240
 Table A.5: RMS error of Vx , Vt and Vr at x/R = 0.1756

Vx (10−1 ) Vt (10−1 ) Vr (10−1 )

GRID-A 1.025 2.262 1.175
GRID-B 1.210 2.519 1.393
GRID-C1 0.885 1.929 1.179
GRID-C2 0.823 2.274 1.083
GRID-C3 0.889 2.584 1.163
GRID-D1 0.857 1.736 1.153
GRID-D2 0.775 2.046 1.042
GRID-D3 0.877 2.277 1.165
Peng et al. 1.574 0.403 1.592

Table A.6: RMS error of Vx , Vt and Vr at x/R = 0.2386

Vx (10−1 ) Vt (10−1 ) Vr (10−1 )

GRID-A 1.238 1.702 0.951
GRID-B 1.135 1.496 1.018
GRID-C1 1.156 1.258 1.184
GRID-C2 1.224 1.423 0.810
GRID-C3 1.876 2.055 0.589
GRID-D1 1.904 1.238 0.838
GRID-D2 1.179 1.256 0.803
GRID-D3 0.994 1.484 1.022
Peng et al. 1.652 0.619 1.747

Hsiao and Pauley (1999) and by Qiu et al. (2013). For the pressure side, high pressure
on the trailing edge as predicted by the present method and by Qiu et al. (2013) was
not shown in the numerical results by Hsiao and Pauley (1999). However the general
pressure distributions on the pressure side are similar.

241
 2
Experimental
Peng et al. (SST)
GRID-A-Y+5
GRID-B-Y+5
GRID-C3-Y+5
GRID-D3-Y+2
1.5

Vx
1

0.5

0
-0.4 -0.2 0 0.2 0.4
θ
3.5
Experimental
Peng et al. (SST)
GRID-A-Y+5
GRID-B-Y+5
GRID-C3-Y+5
GRID-D3-Y+2
3
Vt

2.5

1.5
-0.4 -0.2 0 0.2 0.4
θ
1.5
Experimental
Peng et al. (SST)
GRID-A-Y+5
GRID-B-Y+5
GRID-C3-Y+5
1 GRID-D3-Y+2

0.5
Vr

-0.5

-1
-0.4 -0.2 0 0.2 0.4
θ

Figure A.12: Vx , Vt and Vr across vortex core at x/R = 0.1756

242
 2
Experimental
Peng et al. (SST)
GRID-A-Y+5
GRID-B-Y+5
GRID-C3-Y+5
GRID-D3-Y+2
1.5

Vx
1

0.5

0
-0.4 -0.2 0 0.2 0.4
θ
3.5
Experimental
Peng et al. (SST)
GRID-A-Y+5
GRID-B-Y+5
GRID-C3-Y+5
GRID-D3-Y+2
3
Vt

2.5

1.5
-0.4 -0.2 0 0.2 0.4
θ
1.5
Experimental
Peng et al. (SST)
GRID-A-Y+5
GRID-B-Y+5
GRID-C3-Y+5
1 GRID-D3-Y+2

0.5
Vr

-0.5

-1
-0.4 -0.2 0 0.2 0.4
θ

Figure A.13: Vx , Vt and Vr across vortex core at x/R = 0.2386

243
(a) Numerical results by Hsiao and Pauley (b) Numerical result by Qiu et al. (2013)
(1999)

P
P 1.0001
1.0000 0.8072
0.8071
0.6144
0.6143
0.4215
0.4214
0.2286
0.2286
0.0357 0.0358

-0.1571 -0.1571
-0.3500 -0.3499
-0.5429 -0.5428
-0.7357 -0.7357
-0.9286
-0.9285
-1.1214
-1.1214
-1.3143
-1.3143
-1.5071
-1.5071
-1.7000
-1.7000

(d) GRID-D3
(c) GRID-C3

Figure A.14: Pressure coefficient on suction side.

244
(a) Nmerical results by Hsiao and Pauley
(b) Numerical results by Qiu et al. (2013)
(1999)

P
P 1.0000
1.0000
0.8571
0.8571
0.7143
0.7143
0.5714
0.5714
0.4286
0.4286
0.2857 0.2857

0.1429 0.1429
0.0000 0.0000
-0.1429 -0.1429
-0.2857 -0.2857
-0.4286
-0.4286
-0.5714
-0.5714
-0.7143
-0.7143
-0.8571
-0.8571
-1.0000
-1.0000

(d) GRID-D3
(c) GRID-C3

Figure A.15: Pressure coefficient of pressure side.

245
In Figure A.16, isosurfaces of the vortex structure based on GRID-C3 are plotted
using Q = 100, 000, which is the second invariant of the velocity gradient tensor. The
experimental tip vortex trajectory is presented using a black dotted line to show the
comparison. The strength of tip vortex core is presented using the normalized axial
velocity color label.

Figure A.16: Vortex trajectory for GRID-C3.

246
Appendix B

Wave Run-up Simulation

The wave run-up simulations were performed for a single truncated square cylinder
using the laminar model. For the simulation with the structure in the domain, a
complete computational domain was used.

B.1 Model Experiment

Stansberg et al. (2005) carried out an experimental study of wave run-up around
a vertical truncated square cylinder, which was used as an ITTC benchmark study
during the 27th Workshop, Nantes, France (2013). In the current work, a similar
model set-up was used for validation purposes. In full scale, the diameter of the
square cylinder d is 16 m and the draft Td is 24 m. A scale factor 1:48.93 was used
in the experiment for the truncated square cylinder. Wave elevations were measured
at twelve locations as shown in Figure B.1. The coordinates of the wave probes are
summarized in Table B.1. The wave matrix for this study is given in Table B.2.

247
 Table B.1: Wave probe locations.

Series x/d y/d

A (0.000, 0.000, 0.000, 0.000) (0.503, 0.592, 0.797, 1.000)
B (0.429, 0.492, 0.637, 0.780) (0.429, 0.492, 0.637, 0.780)
C (0.503, 0.592, 0.797, 1.000) (0.000, 0.000, 0.000, 0.000)

Figure B.1: Wave probes locations

B.2 Computational Domain and Grid Generation

For the wave tank validation, the length of the computational domain was set as
L=6λ, and the width B = 32d. Here λ is the wavelength and d is the cylinder
diameter. The schematic diagram of the computational domain with wave probes
arrangements is shown in Figure B.2. Figure B.3 presents the computational grid
for wave tank validation. The grids for the 3D wave tank validation were generated
using blockMesh, topoSet, and refineMesh utilities which are available in OpenFOAM.
The computational grid for a single truncated square cylinder was generated using
navalSnappyHexMesh.

248
Figure B.2: Schematic diagram of computational domain for wave tank validation.

Figure B.3: The computational grid for wave run-up simulation.

For the numerical wave tank, a symmetry plane was used at the centreplane of the
computational domain to reduce the number of grids. The symmetry plane was set
at y/d = 0. The wave generation boundary was applied to the inlet for velocity,
and waveInletAlpha was used for the volume fraction alpha, which was calculated
using a specified wave type in the wave dictionary. The inlet boundary was placed at

Table B.2: Wave matrix for wave run-up simulation.

λ/d = 14.04 λ/d = 21.94

H/λ = 1/10 T12S110 T15S110
H/λ = 1/16 T12S116 T15S116

249
longitudinal distance 2λ in front of the cylinder, the outlet boundary was placed at
x/λ = 4 in the downstream direction, and the side boundaries were placed at distance
16d from the center of the structure on both sides in case of the wave structure
interaction simulation. The structure was defined as a non-slip wall boundary. The
boundary conditions used in the numerical computations are listed in Table B.3. In
this work, Stokes second-order wave was used based on the wave theory applicability
depending on the wave parameters used during the model test.

B.3 Numerical Set-up

Three computational grids namely GR-01 which has 10 cells per wave height and
80 cells per wave length, GR-02 which has 14 cells per wave height and 120 cells per
wavelength, and GR-03 which has 20 cells per wave height and 160 cells per wavelength
were used in the numerical computation. Four time step sizes, TS04(Tw /∆t = 1000),
TS05(Tw /∆t = 2000), TS06(Tw /∆t = 2500), TS07(Tw /∆t = 3000) were used in the
wave tank validation study.

B.4 Wave Generation Modules

The validation studies for the two wave generation modules were performed for
the wave T15S110 shown in Table B.2 using time step TS05. For Module-02, the
computational domain was extended 2λ in front of the inlet boundary. For the

Table B.3: Boundary conditions for numerical wave generation.

Parameters inlet outlet bottom left/right top cylinder
U waveInletVelocity fixedValue noSlip slip pressureInletOutletVelocity noSlip
prgh fixedFluxPressure zeroGradient zeroGradient slip totalPressure fixedFluxPressure
νt fixedValue zeroGradient nutUSpaldingWallFunction slip zeroGradient nutUSpaldingWallFunction
k fixedValue inletOutlet fixedValue slip inletOutlet fixedValue
ω fixedValue inletOutlet fixedValue slip inletOutlet fixedValue
α waveInletAlpha zeroGradient zeroGradient slip inletOutlet zeroGradient

250
 2.0
Module-01
Module-02
Analytical
1.5

1.0
Wave Height, ζ/A1

0.5

0.0

-0.5

-1.0

-1.5
20 21 22 23 24 25
Time, t/Tw

Figure B.4: Comparison of wave elevation at x/λ = 0 for T15S110 using Tw /∆t =
2000.

2.0
Module-01
Module-02
Analytical
1.5

1.0
Wave Height, ζ/A1

0.5

0.0

-0.5

-1.0

-1.5
20 21 22 23 24 25
Time, t/Tw

Figure B.5: Comparison of wave elevation at x/λ = 0 for T15S110 using Tw /∆t =
2000.

251
analytical and numerical simulations, the Stokes second order wave was used. The
numerical wave elevation for the two modules are compared for wave probe 03 using
grid resolution GR-02 and GR-03. The difference in the two modules in not significant,
but there is a significant difference in the total computational grid numbers. On this
basis, wave generation Module-01 is used for inlet wave generation with the relaxation
based wave damping method for the purposes of the self-propulsion simulation.

B.5 Grid Resolution Study

For the grid dependency analysis, the numerical and analytical time history of wave
elevation using time step size TS05 for three grid resolutions for wave T15S110
at x/λ = 0.0 are compared in Figure B.6. From the comparison, it is seen that
the numerical time history of wave elevation shows a good agreement using GR-03.
Another comparison is made for wave T15S116 using time step TS05 in Figure B.7
at x/λ = 0. It is observed that a small time step is necessary for the wave T15S116
for better numerical prediction.

B.6 Time Step Dependency

The dependency on time step size was studied using GR-03 for wave T15S116. From
Figure B.8, it is seen that time step TS05 gives the numerical time history of wave
elevations with reasonable accuracy. However, the wave crest is over predicted. Figure
B.9 shows the comparison between numerical and analytical results for wave T15S116.
A better agreement was found for both wave crest and trough. The time step
dependency was also studied for T12S110 and T12S116 using GR-03 at x/λ = 0
and the comparisons are presented in Figures B.10 and B.11 respectively. From the
grid resolution and time step study, it is seen that GRID-03 with time step TS-05

252
 2.0
Analytical
GR-01
GR-02
1.5 GR-03

1.0
Wave Height, η / A1

0.5

0.0

-0.5

-1.0

-1.5
20 21 22 23 24 25
Time, t/Tw

Figure B.6: Comparison of wave elevation at x/λ = 0 for T15S110 using Tw /∆t =
2000.

gives acceptable numerical wave height compared to the analytical solution. Please
note that the limitedLinear (OpenFOAM, 2019) divergence scheme was used in the
grid dependency and time step size dependency study.

B.7 Effect of Divergence Scheme

The effect of divergence scheme was investigated for T12S116 using limitedLinear and
limitedCubic (OpenFOAM, 2019) using different time step sizes. The comparison is
shown in Figure B.12. The close view of the wave trough on the generated waves
shows that the same level of accuracy can be obtained with large time step size using
the limitedCubic divergence scheme. The difference between the numerical results
for TS-04 and TS-05 is not significant compared with the analytical wave elevation.

253
 2.0
Analytical
GR-02
GR-03
1.5

1.0
Wave Height, ζ/A1

0.5

0.0

-0.5

-1.0

-1.5
20 21 22 23 24 25
Time, t/Tw

Figure B.7: Comparison of wave elevation x/λ = 0 for T15S116 using Tw /∆t = 2000.

2.0
Analytical
TS-03
TS-04
1.5 TS-05

1.0
Wave Height, η / A1

0.5

0.0

-0.5

-1.0

-1.5
20 21 22 23 24 25
Time, t/Tw

Figure B.8: Comparison of wave elevation at x/λ = 0 for T15S110 using GR-03.

254
 2.0
Analytical
TS-03
TS-04
1.5 TS-05

1.0
Wave Height, ζ/A1

0.5

0.0

-0.5

-1.0

-1.5
20 21 22 23 24 25
Time, t/Tw

Figure B.9: Comparison of wave elevation at x/λ = 0 for T15S116 using GR-03.

2.0
Analytical
TS-05
TS-06
1.5

1.0
Wave Height, ζ/A1

0.5

0.0

-0.5

-1.0

-1.5
20 21 22 23 24 25
Time, t/Tw

Figure B.10: Comparison of wave elevation at x/λ = 0 for T12S110 using GR-03.

255
 2.0
Analytical
TS-05
TS-06
1.5

Wave Height, ζ/A1 1.0

0.5

0.0

-0.5

-1.0

-1.5
20 21 22 23 24 25
Time, t/Tw

Figure B.11: Comparison of wave elevation at x/λ = 0 for T12S116 using GR-03.

2.0
Analytical
limitedCubic-TS-03
limitedCubic-TS-04
1.5 limitedCubic-TS-05
limitedLinear-TS-05

1.0
Wave Height, ζ/A1

0.5

0.0

-0.5

-1.0

-1.5
20 21 22 23 24 25
Time, t/Tw

Figure B.12: Effect of divergence scheme on wave generation.

256
 2.0
Analytical
Laminar
realizable k-ε
1.5 k-ω SSTCC

1.0
Wave Height, ζ/A1

0.5

0.0

-0.5

-1.0

-1.5
20 21 22 23 24 25
Time, t/Tw

Figure B.13: Effect of turbulence model on wave generation.

B.8 Effect of Turbulence Model

The effect of turbulence model on the numerical wave generation was studied for
three conditions, laminar, realizable k − , and k − ω SSTCC turbulence model.
A 0.10% turbulence with eddy viscosity ratio µt /µ = 10 was used for turbulence
parameters calculation. The comparison shown in Figure B.13 between the analytical
wave elevation and the numerical wave elevation for the different turbulence model
shows almost similar identical results.

B.9 Wave Run-up Results

The non-dimensional horizontal and vertical forces acting on the square cylinder
for wave T12S110 are plotted in Figure B.14. Figure B.15 shows the forces acting
on the square cylinder for wave T15S110. From these figures, it is seen that the

257
 4.0

3.0

2.0

1.0

Fx/(ρgd2A1/4) 0.0

-1.0

-2.0

-3.0

-4.0

-5.0

-6.0
35 36 37 38 39 40
Time, t/Tw

(a) Horizontal force

10.0

9.5

9.0
Fz/(ρgd2A1/4)

8.5

8.0

7.5

7.0

6.5

6.0
35 36 37 38 39 40
Time, t/Tw

(b) Vertical force

Figure B.14: Time history of wave force on the cylinder for wave T12S110.

258
 3.0

2.0

1.0

Fx/(ρgd A1/4) 0.0

2

-1.0

-2.0

-3.0

-4.0
35 36 37 38 39 40
Time, t/Tw

(a) Horizontal force

7.5

7.0

6.5

6.0
Fz/(ρgd2A1/4)

5.5

5.0

4.5

4.0

3.5

3.0

2.5
35 36 37 38 39 40
Time, t/Tw

(b) Vertical force

Figure B.15: Time history of wave force on the cylinder for wave T15S110.

259
 3.0
Experimental-A2
Numerical-A2
2.5 Experimental-A4
Numerical-A4

2.0
Wave Height, η / A1
1.5

1.0

0.5

0.0

-0.5

-1.0

-1.5
35 36 37 38 39 40
Time, t/Tw

Figure B.16: Wave elevation comparison for wave T12S110 at A2 and A4.

3.0
Experimental-B2
Numerical-B2
2.5 Experimental-B4
Numerical-B4

2.0
Wave Height, η / A1

1.5

1.0

0.5

0.0

-0.5

-1.0

-1.5
35 36 37 38 39 40
Time, t/Tw

Figure B.17: Wave elevation comparison for wave T12S110 at B2 and B4.

260
 3.0
Experimental-C2
Numerical-C2
2.5 Experimental-C4
Numerical-C4

2.0
Wave Height, η / A1
1.5

1.0

0.5

0.0

-0.5

-1.0

-1.5
35 36 37 38 39 40
Time, t/Tw

Figure B.18: Wave elevation comparison for wave T12S110 at C2 and C4.

non-dimensional forces for the short wave are higher than those for the long wave
having the same wave steepness. The wave run-up is measured at 12 different
locations. The time series of the wave elevations are compared in Figures B.16, B.17,
and B.18 for wave T12S110. The predicted wave run-up shows good agreement with
the experimental results. The instantaneous free surface for the wave T12S110 with
the single truncated square cylinder at different times is presented in Figure B.19.

The effect of the turbulence model on wave run-up for a single cylinder is also
investigated for wave T12S116. The k − ω SSTCC, and realizable k −  turbulence
models were used in the numerical simulation. The turbulence parameters were
calculated using 0.10 % turbulence intensity. The wave run-up values are compared
at probe A2, A4, B2, B4, C2 and C4 in Figures B.20 to B.25. From the comparison,
it can be seen that the wave run-up using the k − ω SSTCC turbulence model is
under predicted. The wave run-up is damped out considerably. The turbulence eddy

261
(a) t = 30Tw + 0Tw (b) t = 30Tw + Tw /8

(c) t = 30Tw + Tw /4 (d) t = 30Tw + 3Tw /8

(e) t = 30Tw + Tw /2 (f) t = 30Tw + 5Tw /8

(g) t = 30Tw + 3Tw /4 (h) t = 30Tw + 7Tw /8

Figure B.19: Instantaneous free surface for wave T12S110.

262
 3.0
Experimental-A2
Laminar
2.5 k-ω SSTCC
realizable k-ε

2.0
Wave Height, η / A1
1.5

1.0

0.5

0.0

-0.5

-1.0

-1.5
35.0 36.0 37.0 38.0 39.0 40.0
Time, t/Tw

Figure B.20: Wave elevation comparison for wave T12S116 at A2.

3.0
Experimental-A4
Laminar
2.5 k-ω SSTCC
realizable k-ε

2.0
Wave Height, η / A1

1.5

1.0

0.5

0.0

-0.5

-1.0

-1.5
35.0 36.0 37.0 38.0 39.0 40.0
Time, t/Tw

Figure B.21: Wave elevation comparison for wave T12S116 at A4.

263
 3.0
Experimental-B2
Laminar
2.5 k-ω SSTCC
realizable k-ε

2.0
Wave Height, η / A1
1.5

1.0

0.5

0.0

-0.5

-1.0

-1.5
35.0 36.0 37.0 38.0 39.0 40.0
Time, t/Tw

Figure B.22: Wave elevation comparison for wave T12S116 at B2.

3.0
Experimental-B4
Laminar
2.5 k-ω SSTCC
realizable k-ε

2.0
Wave Height, η / A1

1.5

1.0

0.5

0.0

-0.5

-1.0

-1.5
35.0 36.0 37.0 38.0 39.0 40.0
Time, t/Tw

Figure B.23: Wave elevation comparison for wave T12S116 at B4.

264
 3.0
Experimental-C2
Laminar
2.5 k-ω SSTCC
realizable k-ε

2.0
Wave Height, η / A1
1.5

1.0

0.5

0.0

-0.5

-1.0

-1.5
35.0 36.0 37.0 38.0 39.0 40.0
Time, t/Tw

Figure B.24: Wave elevation comparison for wave T12S116 at C2.

3.0
Experimental-C4
Laminar
2.5 k-ω SSTCC
realizable k-ε

2.0
Wave Height, η / A1

1.5

1.0

0.5

0.0

-0.5

-1.0

-1.5
35.0 36.0 37.0 38.0 39.0 40.0
Time, t/Tw

Figure B.25: Wave elevation comparison for wave T12S116 at C4.

265
 (a) realizable k − 

(b) k − ω SSTCC

Figure B.26: Turbulence viscosity distribution on the centre plane.

viscosity distribution on the center plane is shown in Figure B.26. It is seen that for
the k ω SSTCC turbulence model, the turbulence viscosity is much higher than that
for the realizable k −  model. This is a common problem with eddy-viscosity based
turbulence models. A numerical treatment must be applied to limit the turbulence
viscosity.

266