Schedule for Presentation Meeting of Master-Degree Thesis in 2018
International Course of Maritime & Urban Engineering
Date :
Venue :
No.

Time

1

09:30 – 09:55

PRADHAN Sujan

SanadaIijima2

09:55 – 10:20

ROYCHOUDHURY
Seemontini

Seismic Performance Evaluation of RC Beam Considering Buckling of Longitudinal
Reinforcing Bar

Toda

CFD Study on KCS Maneuvering in Calm Water

Iijima3

10:20 – 10:45

KHOIRI Putika Ashfar

Nishida/Irie

Parameter Estimation of a Distributed Hydrological Model in the Ibo River Basin Using
the Polynomial Chaos Expansion Method

Iijima4

10:45 – 11:10

Thant Zin Htun

Suzuki

Motion Analysis of ROV Coupled With Tether

Iijima5

11:10 – 11:35

Myo Zin Aung

Umeda

6

13:00 – 13:25

ADILAH Ahmad

IijimaStructural Analysis of a Floating Vertical Axis Wind Turbine

Suzuki7

13:25 – 13:50

ZHOU Dali

Kashiwagi

Numerical Calculation of Ship Hydroelasticity Using Enhanced Unified Theory

Suzuki8

13:50 – 14:15

Aung Kyaw Htet

Fujikubo

9

14:30 – 14:55

AMRA Rajuli

Aoki

Adaptive Mitigation Approach for Tsunami Ascending to Meandering River

Irie

10

14:55 – 15:20

Toda

SPIV (Stereo-Particle Image Velocimetry) Measurement for a Self-propelled Ship inRegular Head Waves with Different

Irie

11

15:20 – 15:45

Osawa

Verification and Validation of a Cellular Automaton Simulating Under-Film Corrosion

Irie

16:00 –

Name

August 1st (Wednesday) starting from 09:30
S1-312 Lecture Room

MWANGI
Benson Oyunge
KRISTOV
Ivan William
at S1-323

Supervisor

Title of thesis

Minimum Propulsion Power Prediction of a Ship under Adverse Weather Conditions
with Dynamics of Diesel Engine and Turbocharger Taken into AccountLunch Break

Ultimate Strength of a Continuous Stiffened Panel under Combined Transverse Thrustand In-plane Shear
BreakMeeting for Evaluating Mater-thesis Presentations, and Steering Committee Meeting

Meeting for Evaluation: From 16:00 at S1-323 Meeting Room
Members are Steering committee members and the supervisors

Chairman

IijimaSuzukiInternational Course of Maritime & Urban Engineering / Abstract of Master Thesis

July 27, 2018

Seismic Performance Evaluation of RC Beam Considering Buckling ofLongitudinal Reinforcing Bar
Sujan Pradhan
Concrete Structure Laboratory, Department of Architecture Engineering
Key Words: Full scale test, Performance based design, Reinforced concrete, Safety limit, Seismic assessment

1. Introduction
The strong column-weak beam type frame structures possess
better energy dissipation capacity during the earthquake. Hence,
based on the performance based seismic design, strong
column-weak beam type RC moment resistant frame structures
are effectively used. The safety limit of RC beam members
applied to the seismic performance evaluation of RC buildings 1)
is defined as the smallest of drift ratios at which (1) the flexural
resistance decreases to 80% of its maximum strength, (2) shear
failure occurs after flexural yielding and (3) tensile rupture ofbeam longitudinal reinforcement bar occurs.
However, in a series of static loading experiment of strong
column-weak beam type RC frames 2), damage conditions defined
as above (from (1) to (3)) were not observed even under large
drifts, but the buckling phenomena of the beam longitudinal bars
were visually observed after the cover concrete was peeled off at
the beam ends. This experiment showed that the reduction instrength due to the buckling was related to the safety limit ofstrong column-weak beam type RC moment-resisting buildings.
Many studies have focused on the buckling of longitudinal
reinforcing bars in columns under high axial forces. However,limited studies focused that in beams. Hence, the method toestimate the buckling occurrence point for beam longitudinal
reinforcing bars had been proposed in past research 3). In thecurrent research, the buckling of beam longitudinal reinforcing
bars observed in the previous experiment 2) was experimentally
re-examined and the applicability of the proposed evaluation
method 3) was verified using full scale beam specimens, as theprevious research 2) tested small scaled frame specimens with
common beam structural details and the evaluation method 3) wasapplied to only one common type of beam.

2. Experiments

on

Full

Scale

Beam Specimens

Focusing on Buckling of Longitudinal Reinforcing
Bars
2.1 Specimens and loading system
Three full scale beam specimens were designed based on the
beams in the scaled specimen, each with height of 1540 mm andcross-section of 450 X 700 mm. They had different diameters oflongitudinal reinforcing bars and spacing of shear reinforcing bars,
as shown in Table 1, which were likely to affect bucking of beam
longitudinal bars.
The specimens were tested using the static loading system at
Osaka University. Every specimen was subjected to a constantaxial load of N/N0 = 0.07 (where N is the axial load, and N0 is thecompressive strength of the specimen). The axial load wasapplied to the beam specimens, as the previous experiment 2)reported that the axial load of 7% to the beam compressive

strength was induced due to the restriction to beam elongation by
the columns. Equation (1) gives the axial forces for the south (Ns)and north (NN) vertical jacks during the experiment to maintain abending moment diagram representing that in the previous beam
member 2), as shown in Fig. 1.
Table 1 Specifications of the specimens
Descriptions
Specimen1
Longitudinal
6-D25
bar
Shear
reinforcing
D13@200bars
NS = N/2 + Q. a/lSpecimen 2

Specimen 3

10-D196-D25

D13@200

D10@100

NN = N/2 - Q. a/l(1)

Fig.1 Target bending moment
diagram.
2.2 Experimental results
Fig. 2 shows the comparison of beam shear force vs.
experimental drift ratio relationships. Flexural cracks occurred ateach beam end during the first cycle to R = 0.55% rad, and every
specimen formed a yield mechanism after flexural yielding at thecritical section. Thus, the stiffness degraded in the following
cycles. However, the shear resistance still increased with anincrease of drift before the buckling of longitudinal reinforcing
bars.
The vertical compression cracks along the longitudinal
reinforcing bars on each specimen appeared during the negative
cycle to R = 2.18% rad, as well as the flexural resistancedecreased during the next cycle to R = 3.30% rad, which indicatedthat the buckling of longitudinal reinforcing bars in eachspecimen was initiated from the cycle to R = 2.18% rad.
Comparing specimens 1 and 2, the shear reinforcing bar
spacing was the same but the diameter of longitudinal reinforcing
bars was different. The shear resistance was significantly
decreased for specimen 2 with the smaller diameter of
longitudinal reinforcing bar, whereas it was a little moderate for
specimen 1 after the initiation of longitudinal reinforcing bar
buckling.

Flexural cracking;

Yielding of longitudinal reinforcing bars;

Maximum strength;

Ultimate state with buckling of longitudinal reinforcing bar

300

100

300

0
-100

200

-200

Specimen100
1

-300
-4

-2

0

0
4

2

Drift Ratio (% rad)300

Shear force(kN)

200Shear Force (kN)Shear Force(kN)
Shear Force(kN)

300

200
100
0
-100
-200

Specimen 2

-300

200
100
0
-100
-200

Specimen 3

-300
-4

-2

0

2

4

-4

Drift Ratio (% rad)

-2

0

2

4

Drift Ratio (% rad)calQy= 204
Fig. 2 Beam shear force vs. drift ratio relationships
Q
=
206
y
-200
Qmax= 242
Comparing specimens 1 and 3, the-300
diameter of longitudinal
investigatedstress-strain curve of longitudinal reinforcing bar andreinforcing bars was the same but the shear
bar 0 then compared-4 reinforcing
-2
2
4
with
respective working stress. The buckling
Drift Ratio (% rad)spacing was different. The shear resistance was dropped

significantly for specimen 1 with larger shear reinforcing bar
spacing, whereas it was little moderate for specimen 3 after the
buckling of longitudinal reinforcing bar.
During the static loading experiments, the tensile force pulled
out the longitudinal reinforcing bars from the stub, which affected
the experimentally measured drift. To obtain the pull-out length,
strains of longitudinal reinforcing bars at different depth weremonitored by strain gauges. However, the strain gauge data wasfound to be reliable only for specimen 1. Hence, in this research,
the flexural drift ratio was evaluated excluding the pull-out
behavior and comparison of experimental and analytical drift at
the buckling of the longitudinal reinforcing bar was performed
only for specimen 1. After modifying the flexural drift without
pull-out behavior (Rn), the buckling of longitudinal reinforcing
bars of specimen 1 was found to be initiated during the cycle toRn = 1.52% rad., as shown in Fig. 3.

3. Analytical Evaluation for Drift at Buckling of BeamLongitudinal Reinforcing Bars

resistance strength (σE) became less than the working stressduring the cycle to Rn = 1.52% rad. Hence, the beam longitudinal
reinforcing bars buckling drift for specimen 1 was estimated to be
Rn = 1.52% rad. Fig. 3 shows the experimental and analytical
shear force and drift ratio (without pull-out behavior) relationshipwith experimental and analytical buckling drifts. This shows thatthe drift at beam longitudinal reinforcing buckling from theanalysis had a good agreement with that from the experiment.

Fromexperiment;
cyclic cross-section
Fromexperiment;
From From
cyclic cross-section
analysis;analysis;
Buckling
experiment; Buckling
Buckling
driftfrom
cycle
from analysis
Buckling
driftdrift
fromfrom
experiment;
drift cycle
analysis
200
200

(kN)
force (kN)Shear force
Shear

-100

100
100

00
-100
-100
-200
-200
-3
-3

-2-2

3.1 Cyclic cross-section analysis
To investigate the stress-strain behavior of longitudinal
reinforcing bars to evaluate the drift at its buckling during thecyclic loading, the cyclic cross-section analysis was performed byfiber approach. Trilinear model and Menegotto Pinto steel model
were used for concrete and reinforcing bars, respectively. The
relation between the drift ratio and the curvature was evaluatedassuming deflection characteristics of the beam specimen to
perform the analysis with the curvature equivalent to theexperimental drift without pull-out behavior and the stress-strain
behaviors of longitudinal reinforcing bars were investigated.
3.2 Estimation for longitudinal reinforcing bars buckling
The estimation of longitudinal reinforcing bars buckling drift
was based on the previous evaluation method 3). The buckling
resistance strength ( σ E) of longitudinal reinforcing bars was
obtained by using the equation (2).

E 

4   2.  I
 Er
L2B  Ar(2)

where, I and Ar are the moment of inertia and thecross-sectional area of longitudinal reinforcing bar; LB is thecenter to center spacing of shear reinforcing bars and Er is
equivalent stiffness coefficient.
The value of equivalent stiffness coefficient was evaluatedbased on the method used by Akiyama et. al.4). The buckling
resistance strength (σE) was calculated at every point of the

-1 -1

0

0

1

Drift
ratio
(% rad)Driftratio(% rad)

1 2

23

3

Fig.3 Comparison of the drifts at bucking from the experimentand analysis

4. Conclusions
1. The buckling of beam longitudinal reinforcing bars observedin the previous study was experimentally verified in full scalebeam specimens.
2. Buckling initiation drift was found to be independent to thediameter of longitudinal reinforcing bars and the spacing ofshear reinforcing bars within this series of tests.
3. The diameter of longitudinal reinforcing bars and the spacingof shear reinforcing bars affected the drop in shear resistanceof RC beam after the buckling of longitudinal reinforcing bar.
4. The evaluation method for beam longitudinal reinforcing bar
buckling proposed in the past research was found to be
applicable to full scale beam members.

References
1) AIJ: Guidelines for Performance Evaluation of Earthquake
Resistant Reinforced Concrete Buildings (draft), 2004.
2) Yoon, R., Sanada, Y. and Akahori, T.: Journal of Advanced
Concrete Technology, Vol.15, pp.5 44-557,2017.
3) Akahori, T.: Master Thesis, Department of ArchitecturalEngineering, Osaka University, 2016.
4) Akiyama, M., Naito, H. and Suzuki, M.: Journal of JSCE,
Vol.58, No.725, pp.113-129,2003.

International Course of Maritime & Urban Engineering / Abstract of Master Thesis

July 27, 2018

CFD Study on KCS Maneuvering in Calm Water
Seemontini RoyChoudhury
Hull Form Design Sub-Area, Department of Naval Architecture and Ocean Engineering
Key Words: CFD, Ship Motions, Neutral Rudder angle, Flow field visualizations, Free Maneuver

1. Introduction
As assessed in the reports by the Maneuvering Committee of theInternational Towing Tank Conference1), modern prediction tools
for ship maneuverability mainly include empirical formulae, fullscale trials, model tests (also known as Experimental Fluid
Dynamics, EFD), system identification and numerical simulation.
The conventional prediction mechanism for ship maneuvering ismodel testing, which provides direct and reliable assessment of the
maneuvering characteristics but is costly and time consuming andmost significantly the flow field is seldom measured. Alternatively,
Computational fluid dynamics (CFD) provides answers to thepitfalls of model testing by providing flow field visualization at any
given point and time of a maneuvering ship. The visualizations are
an aid in understanding the mechanism of the ship maneuver. In
recent years CFD simulations have progressed from captive to freerunning six degrees of freedom (6DOF) conditions. Zigzag
maneuvers of KRISO Container Ship (KCS) implementing an
actual propeller model using Delayed Detached Turbulence Model
has been investigated2). The approach provided satisfactory results
but the computational cost was very high. Reynolds Averaged
Navier Stokes (RANS) maneuvering simulation coupled withvarious propeller modelling for KVLCC2 was studied3). Thoughcomputationally not as exhaustive, the drawback with theaxisymmetric body force propeller model is its inability to provide
propeller side and vertical forces.
In this study, we propose an approach of using RANS coupledwith Osaka University (OU) propeller model for KCS
maneuvering studies. The advantage of this method is that the hullpropeller-rudder interaction is determined without requiring thedetailed modelling of the propeller, thereby reducing
computational cost. The model also provides the propeller side andvertical forces which is an added advantage to axisymmetric bodyforce model.

2. Methodology
CFD has been used to predict the self-propulsion and freemaneuvering characteristics of the KCS model for Froude number
0.26, with six degrees of freedom (DOF) conditions, in calm water,as recommended by CFD workshop 20154). The length between
perpendiculars (LPP) of the model ship is 2.7m for whichmeasurements have also been carried out at IIHR Wave Basin,
University of Iowa.
The grid comprising of 11 blocks has been generated using
Gridgen. Thereafter SUGGAR library has been used for the oversetof grids. CFDShip-IOWA V4.5 (which is an overset, block
structured, and incompressible URANS code coupling with motion
solver) has been used for numerical prediction. The k-ω SST
turbulence model was employed using no wall function. A single
phase level-set method has been used to model the free surface and
6DOF ship motions are considered. The computational domain andboundary condition is shown in Fig. 1

Fig. 1 Computational domain and boundary conditions.

3. Results
3.1 Self propulsion
KCS is a single screw single rudder ship fitted fitted with righthanded propeller, thus moves portside when the rudder angle is set
to zero. It is therefore necessary to rotate the rudder on the
starboard side to maintain straight ahead condition. This angle is
known as neutral rudder angle. The self-propulsion in calm water
was performed with both rudder and propeller equipped with PID
controller. The obtained self-propulsion factors are placed at Table
1.
Table 1 Self Propulsion factors.
Parametern'δ (deg)u/UThrust (N)Sinkage
Trim

CFD
33.61223
0.774161.001176.393540.00207-0.13093

EFD
37.72000
0.75000
1.00000
6.57000
0.00220-0.11000% Deviation
10.890-3.221-0.1172.686
6.062-19.024

The local flow field were compared with OU PIV measurement
data at a location behind the rudder (x/Lpp =1.025) as shown in Fig.
2. The grid used is able to resolve the two vortices near the hub.
The vortex induced vertically lower portion of the propeller alongwith the vortex induced near starboard upward propeller is also
predicted. Thus the same grid was further used in maneuvering
simulations.
3.2 Zigzag 20°/20° Maneuvering
The comparison between CFD and EFD results of zigzag
maneuvers has been undertaken for ±10°/10°,±20°/20° Resultsfor only 20°/20° have been shown in Fig. 3. Fig. 4 shows the freesurface elevation corresponding to point A (maximum yaw), B(zero yaw starboard to port), C (zero yaw port to starboard) and D(minimum yaw) marked on Fig.3. At point B and point D, ship andpropeller wakes are aligned and follows the trajectory. At point A,
the rudder is deflected to port side thus the ship wake near the freesurface is directed to port. The propeller is still at starboard withpropeller wake following the trajectory. Thus the wake is split in
two directions.

(a)

(b)
(b)

A

B

C

D

Fig. 2 Local flow comparison at x/Lpp =1.025 (a) OU PIV(b) OU CFD.

Fig. 4 Free surface elevation for zigzag 20/20 at maximum
yaw(A), Zero yaw (B and D) and Minimum Yaw (C).

A
D

B

C
Fig.3 Starboard side zigzag 20°/20°
3.3 Turning Circle
The comparison between CFD and EFD results of turning
circle maneuvers has been undertaken for ±35°. Results for only
-35° have been shown in Fig. 5. The trajectory is under predictedwith 35° rudder angle. This could be attributed to the simplified
rudder geometry used for this study leading to higher rudder force
prediction and higher speed drops at higher angles of attack. To
analyzed the stated reason, the turning circle maneuver was
performed for 31° and 30° rudder angle. The predicted trajectorywith 30° rudder angle matched well with the EFD data thus
establishing the need to accurately model the rudder geometry.

4. Conclusions
The predicted self-propulsion factors show good agreement with
the EFD data. The propeller revolution is under predicted as theopen water chart for IIHR was not used for computations. In
addition, the CFD wake field at the analyzed plane shows very
good agreement with PIV measurements.
The turning circle maneuver is underpredicted for 35° rudder
angle due to overprediction of the rudder force. With lesser rudder
angle, the prediction was improved.

Fig.5 Turning Circle portside maneuver.
The zigzag maneuver shows good agreement with the EFD data.
The maneuvering flow field was analyzed. CFD is undeniably
attractive to provide analysis and visualization of the flow.
For future work, the rudder geometry will be improved toincorporate the real geometry. Simulations for course keeping inwaves will be performed and analyzed in detail.

References
1) The Manoeuvring Committee.: Final Report andRecommendations to the 25th ITTC, Vol. I, Proc. of 25th
ITTC, 2008.
2) Mofidi, T. and Carrica, P.M.: J of Computers and Fluids,
Vol.96, pp.191-203, 20143) Hosseini, H.S, Wu P.C. and Stern F.: CFD simulations of
KVLCC2 maneuvering with different propeller modeling,
Proc. of Workshop on Verification and Validation of Ship
Maneuvering Simulation Methods (Lyngby, Denmark), 2014.
4) http://www.t2015.nmri.go.jp/Instructions_KCS/Case_2.10/Ca
se_2-10.html

International Course of Maritime & Urban Engineering / Abstract of Master Thesis

July 27, 2018

Parameter Estimation of a Distributed Hydrological Model in the Ibo
River Basin Using the Polynomial Chaos Expansion Method
Putika Ashfar Khoiri
Water Engineering Laboratory, Department of Civil Engineering
Key Words: Hydrologic prediction, Distributed Hydrological Model, Parameter estimation, Polynomial Chaos Expansion

1.

Introduction

3.

Water discharge prediction is important in order to improve
flood management. Modeling actual hydrology systems in riverwatersheds is difficult because accurate modeling needs torepresent land use and other properties of the watersheds. A
distributed hydrological model (DHM) is an adequate tool torepresent spatial information for each sub-catchment of thewatersheds.1) This prediction is sometimes poor due to various
uncertainties in the model and observations. One of their sources is
inaccurate parameter estimation, resulting large misfit arises
between observed and simulated discharges. In order to reduce themisfit, in this research, polynomial chaos expansion (PCE) isutilized for the parameter estimation. In previously research, PCE
techniques have been applied to estimate time-dependent values
for two-biological parameters in an ocean model 2) and otherengineering applications. However, the PCE application forhydrological model is rarely used. The objectives of this study areto verify the applicability of PCE in DHM by revealing theparameter combinations and its statistical significance to improve
simulation results.

2.

Parameter estimation method

In this study, an emulator of the polynomial chaos expansion(PCE) method was adopted to estimate model parameters in DHM
in order to improve the accuracy of simulated discharge. Parameter(θ) with uniform distribution were adapted to a Legendre-Chaos
polynomial. Only two parameters can be optimized in one
parameter space.3) General expression of the Legendre polynomialused in this study is Eq. (1)𝑘𝑚𝑎𝑥𝑓(𝜃, 𝑥𝑡 , 𝑢𝑡 ) = ∑𝑘=0
𝑎𝑘 (𝑥𝑡 , 𝑢𝑡 )∅𝑘 (𝜃) +∈𝑡𝑟𝑢𝑛𝑐 (𝜃)(1)where, 𝑓(𝜃, 𝑥𝑡 , 𝑢𝑡 ) = distributed hydrological model output,𝑎𝑘 (𝑥𝑡 , 𝑢𝑡 ) = expansion coefficients of uncertain input θ in theterms of distance and velocity within time, ∅𝑘 (𝜃) = polynomial oforder k in the parameter space defined by θ, 𝑘𝑚𝑎𝑥 = the maximum
order of polynomial, ∈𝑡𝑟𝑢𝑛𝑐 (𝜃) = truncation error, for the 𝑘𝑚𝑎𝑥 =∞ , the approximation is exact and ∈𝑡𝑟𝑢𝑛𝑐 (𝜃) = 0. 𝑁𝑘 =normalization factor specific to the kth polynomial and independentof θ. Here, maximum order of Legendre polynomial in Eq. (1) with
(𝑘𝑚𝑎𝑥 = 6) used for both parameters in each simulation. Then, 49
model runs of two parameters combination are required. The
results obtained is the global minimum RMSE from two parameter
simulation results.
Due to the orthogonally of the polynomials, conditional
expectation and variance are given by Eq. (2), (3)𝐸(𝑓(𝑥, 𝑢, 𝜃)|𝜃) = 𝑎0 𝑡

(2)𝑣𝑎𝑟(𝑓(𝑥, 𝑢, 𝜃)|𝜃) = ∑𝑛𝑘=1 𝑎𝑘 2 (𝑥, 𝑢)𝑁𝑘

(3)

Calculation methods

3.1 Distributed Hydrological Model

Fig. 1 The schematically diagram of distributed hydrologicalmodel used in this study. 4)
Figure 1 shows the structure of Distributed HydrologicalModel used in this study4) This model consists of 3 sub-model: 1)
forest blocking and evapotranspiration sub-model, 2) penetrationand transpiration sub-model, and 3) slope and river channel
outflow sub-model. In the forest sub-model, the arrival rainfall is
considered as a forest layer. Rainfall reached the crown of theforest canopies and moves to the ground surface after filling tank I(the crown part) and tank II. After water in tank II reached its
maximum storage, the water outflows over ground surfaces. When
the water in the tank II exceeds the limit, the water moves to tank
III next. In the penetration and transpiration sub-model, the rainwater migrates to soil, and stored into tank III and IV. The storage
functions methods in Eq. (4) is used for the calculation.
𝑑𝑆𝑖
𝑑𝑡

𝑃

= 𝑟𝑖 − 𝑄𝑖 , 𝑆𝑖 = 𝐾𝑖 ∙ 𝑄𝑖 𝑖

(4)where, Si : storage of tank i (mm), t : period (h), ri : effective amountof rainfall (mm/h), Qi = runoff height (mm/h), Ki : storage constantof tank i, Pi = storage power constant of tank i, i : tank number （IIIor IV）.
In the slope and river outflow sub-model, the kinematic wave
model given in Eq. (5), (6) is used.
𝜕ℎ 𝜕𝑞
+
= 𝑅𝑒
𝜕𝑡 𝜕𝑥
𝑞 = 𝛽 ∙ ℎ𝑚

(5)
(6)

where, h = water depth (mm), q : unit width discharge (mm2/ sec),x = distance from the top edge of the grid (mm), β, m = model
parameters.

3.2

Study area, Data and calculation setting
The study area is the Ibo River basin in Hyogo Prefecture,
Japan. The area is divided into 197 grids. The size of the grids is 2
km x 2 km. The discharge observation stations are located at
Tatsuno, Yamazaki and Kamigawara. Hourly observed water
discharge in 2015 was used.5) The input data are hourly
precipitation data from XRAIN with 250 m x 250 m resolution
from January 1st, 2015 to December 31st, 2015. The period of
calculations is separated into three parts, which is in the period oflarge flooding, medium flooding, and the period before large
flooding occurs. Sensitivity analysis was conducted before theparameters is put for estimation and the effect of parameter
optimization on simulation is evaluated. Parameter estimation was
performed in 4 cases with the target parameters which were storage
constant of tank III and IV (K3, K4), maximum storage in tank I(S1), layer A-thickness (D), and effective soil layer thickness of
16% and 50% runoff (D16, D50).

4.

Results and discussion

Fig. 2. The results of simulated discharge estimated by PCE

All parameters have not been estimated here, while only
parameters that have high sensitivity for discharge simulations are
determined. The estimated parameters vary in space and time. For
example, at Kamigawara station, based on results from the seriesof sensitivity analyses, parameters D16 and K4 were selected as atarget to be estimated (see Table 1). The parameters D16 and K4 are
put into PCE calculations and the optimum values of these twoparameters were obtained.
The simulated discharge is better reproduced when using the
optimum parameter values. During the period of large discharge
from June 20th, 2015 to August 15th, 2015, at Yamazaki station, thecombination of maximum storage in tank I (S 1) and storage
constant in tank III (K3) give a largest effect on simulated discharge
by reduce the RMSE by 29 % from the initial simulation.
The amount of discharge increased during the large discharge
period in the case with estimated parameters in Kamigawara station
and Tatsuno station, the peak discharge on July 18th, 2015 is higher
than initial simulation (see Fig.2).
Table 1. List of optimum parameter value during large
discharge period.
Kamigawara
Parameter name
Effective soil layer
thickness, of 16% runoff(D16) (mm)Storage constant of tank
IV (K4) (mm23/25h2/25)
TatsunoEffective soil layer
thickness, of 50% runoff(D50) (mm)Storage constant of tank
III (K3) (hour)
Yamazaki
Maximum storage in tank
I (S1) (mm)Storage constant of tank
III (K3) (hour)

Parameter
value

RMSE
Initial
PCE

11
18.66

14.2

24.08

20.56

21.7249

15.25

90

40
6

2.16
4

Because the PCE method can estimate the limited number ofparameter, further development is required for better optimizationof all parameters in DHM. In addition, further evaluation to know
the effect of rainfall data distributions should be carried out toobtain more accurate discharge simulation.

5.

Conclusions

The parameter estimation method by PCE is capable to reducethe misfit between observed and simulated discharges in the Ibo
River basin at 3 observation stations. The storage function of tank
III and IV (K3, K4) with other parameters is important andsubstantial to acquire the better simulation. The PCE method canminimize the amount of exhaustive sampling to obtain optimal
setting of parameters.

References
1) Fukushima, Y.：A model of river flow forecasting for a small
forested mountain catchment ， Hydrological Processes ，Vol.2，pp.167-185，1988
2) Mattern, Jann Paul., Fennel, Katja., Dowd, Michael.:
Estimating time-dependent parameters for a biological ocean
model using an emulator approach. Journal of Marine
Systems. 96-97. pp. 32-47, 2012.
3) Xiu, D.B., Karniadakis, G.E.: The Wiener–Askey polynomialchaos for stochastic differential equations. SIAM J. Sci.
Comput. 24 (2), pp. 619–644, 2002.
4) Oda, K., Irie, M., Toi, H., Ishizuka, M., Tanaka, K.: Parameter
estimation of a distributed hydrological model using an
adjoint method, Journal of JSCE B1 (Hydraulic Engineering)(submitted).
5) Ministry of Land, Infrastructure, Transport and Tourism：
Hydrological
and
water
quality
database,
http://www1.river.go.jp，accessed on 2018/3/31.

International Course of Maritime & Urban Engineering / Abstract of Master Thesis        July 27, 2018



Motion Analysis of ROV Coupled with Tether
Thant Zin Htun
Marine Mechanical System Engineering Laboratory
Naval Architecture and Ocean Engineering Department, Osaka University
Key Words: ROV, tethered cable, ANCF, Boundary condition, Numerical simulation

1. Introduction
Continuous explorations by humans to extend theiraccessibilities in the deep oceans are demanding a drastic increasein the deployment of underwater vehicles such as AUVs
(Autonomous Underwater Vehicles) and ROVs (Remotely
Operated Underwater Vehicles) in the deep oceans. Thanks to thereal-time telepresence to give the operator a sense of being in theplace of the machine, and continuous real-time data transmission,
the tether-operated underwater systems have been the favored
ones among existing means in deep sea’s explorations. ROV canbe upgraded by mounting additional equipment to expand thevehicle’s capabilities. The tether in low tension state, meanwhile,
significantly affects the ROV’s motions by creating bending andtorsional moments combined with hydrodynamic forces acting on
it. The interactions between ROV and tether in real operations
complicate the job of ROV’s operators. Therefore, an accurate
motion simulator is prerequisite for training the ROV’s pilots andof predictive assessments.
In this study, a numerical motion analysis of ROV coupledwith tether is developed by using 24 dof ANCF (Absolute Nodal
Coordinate Formulation)1) beam elements for kinematics of thetether, which can also account for the torsion and shear. The
mathematical modelling of ROV’s motion is developed by usinghydrodynamic coefficients, which were obtained by the numericalforced oscillation method presented in the previous research2).
The elastic forces developed in the tether, thrusters’ forces, andhydrodynamic forces are treated as the external forces. In thepreceding works3),4), the authors presented initial curl effect of thetether on ROV’s motion and defined the boundary condition
between the tether and ROV by applying gradient deficient cablebeam elements which cannot describe torsion and shear. The
anisotropic bending stiffness of the tether was also investigatedby conducting cantilever bending test for four different directions.
The focus of this research is to modify the previous researchand develop fully-coupled motion analysis between ROV andtether. In this research, the authors presented the effect of
torsional moment restored in tether on ROV’s yaw motion, and to
modify initial conditions of the tether and boundary condition
between the tether and ROV. Internal viscous damping was
implemented into the motion equation of cable and its effect on
cable motion is also discussed. The internal damping forces andtorsional rigidity of the tether are evaluated, developing atorsional pendulum simulation and the results are validated byconducting the torsional pendulum experiment of the cable. The
torsional effect of the tether on ROV’s yaw angle is also
discussed by performing turning test of ROV in simulation.

2. CU-ROV and Tether
2.1 Overview of CU-ROV
CU-ROV, which was developed by JAMSTEC is used as the
sample ROV in this research. It is equipped on the Deep Sea
Scientific Drilling Vessel “CHIKYU” and used for observation of
the drilling points and drilling holes on reinstallation of thedrilling pipes. Two LED lights and a wide-angled camera aremounted at the front of the ROV to provide visuals for the pilot,

and it is propelled by a total of five propellers, two of which areinstalled on the upper face while the other threes, at the stern. The
principle dimensions and characteristics of CU-ROV are shownin Table 1 and Fig. 1.
Table 1

Characteristics of CU-ROV

Weight in air
Length×Breadth×Height
Speed
Operational depth
Table 234.65 kg
0.6 m×0.365m×0.325m
02.0 knot
7000 mCharacteristics of Tethered Cable

Weight per unit length

0.1342 kg/m

Diameter
Cross section area
Geometric Moment of inertia

0.01 m
7.854×10-5 m2
4.909×10-10 m4

Fig.1 Overview of CU-ROV (developed by JAMSTEC)
2.2 Modelling Motion Equation of ROV
The dynamics of the ROV is described by the motion in thebody-fixed coordinate system which is defined at the ROV’s
geometric center.

!" = $(1)! = #, %, &, ', (, ) * , ! = #$ , #& , #' , ($ , (& , (')(2)

where M accounts for the total mass matrix of ROV, ! is theacceleration vector and ! represents the generalized external
force vector applied on ROV.
2.3 Modelling Motion Equation of the Tether
A non-incremental absolute nodal coordinate formulation(ANCF) for three-dimensional beam elements is used to analyze
the kinematics of the tethered cable. In this formulation, theposition of each node is described by absolute nodal coordinatesand gradients of the positions with respect to material coordinates
instead of infinitesimal angles. A distinguishing feature in ANCF
is the use of slope vectors to describe the orientation of thecross-section of the beam, thus, continuity of slope is achieved.
This method leads to a constant mass matrix.
& &
& &
&
&
&
! = $%& $%,($%,) $%,*$%+, $%+,,($%+,,)$%+,,*&(3)

! = # $, &, ' ((4)! = !# $ !& $ !' $ !( $ !) $ !* $ !+ $ !, $(5)

where S is the element shape function for 24 DOFs beam, andeach node is defined by the position vector ! and the three
gradient vectors ! ,$ , !,% , ! ,& .

! + !# $ + %&' $+() $ + (* $-$, + (- $-$, = /0

(6)

where
and !" represents the element mass matrix and added
mass matrix respectively. !" represents the initial coordinate
vector of the cable, and ( ) is the differentiation of global
coordinate vector w.r.t time. !"# is the internal viscous damping
factor and !" , !" , !" represent the stiffness matrices
corresponding to the axial deformation, bending and torsion
respectively. !" represents the generalized external force vector.

Fig. 3 Trajectories of ROV’s Turnings (Showing the Effect ofTorsional Moment Restored in Tether on ROV’s Turning Motion
Propelled by Port /Starboard Side Thruster with 6 rev/s for Each
Direction in Simulation)

2.4 Boundary Condition between the Tether and ROV
The boundary condition is defined such that the cable’s centerline at the sling point, which is exactly same as the tangential
vector of the cross-section of the Euler-Bernoulli’s beam,
coincides with the negative normal vector of the ROV’s uppersurface. The coupling terms are, thereby, established between
ROV and the tether at the attached node.
!' = $ ! =&'()&'(*(+,.(+,)&'(*-&'(.(+,*&'(.(+,)&'(* + (+,.(+,*+,-.
!"_$%"&',) = --/01
0

-/01
+,-.
0

&'()(+,*(+,.(+,)(+,* + &'(.&'(*&'(.(+,)(+,*-(+,.&'(*-(+,)
0(+,.&'() 0&'(.&'() −100 ∙ -0' ∙ |!"_$%"&',) |1(7)

where E (φ, θ, ψ) is the Euler’s rotation matrix, n' represents thenegative normal vector of the ROV’s upper face, and !"_$%"&',) isthe gradient vector along the center line of the cable at the slingpoint.

3.2 Future Endeavors
In the present work, Berzeri’s T1 model (Berzeri andShabana)5) was used to develop the bending forces of the finite
elements for better simulation execution. In this model, theconstant curvature assumption is used when the longitudinal
deformation is not quite large. However, it has been observed thatthe program converges slowly when increasing the number of
elements for thin beam models compared to other beam models.
That convergence problem is due to constant bending strains overthe length of the element. Meanwhile, the inclusion of torsion
effect to the tether shows considerable effect on ROV’s yaw
motion. In the current work, the assumed shape function is
obtained by applying linear interpolation in x and y directions andcubic interpolation in z direction, where the z direction is initiallyalong the central axis of the beam. To improve the accuracy ofbending and torsional effects, higher order elements will be usedby introducing additional shape functions6). The depth varying
ocean current velocity, and the movement of mother ship andcage will also be taken into account in the cable dynamics.
References

Fig. 2 Boundary Condition between ROV and Tether
3. Conclusions and Future Endeavors
3.1 Conclusions
In this research, the stiffness matrices and constant mass
matrix that can also account for rotary inertia in three dimensional
24dof beam were explicitly presented. The cable is considered asa very slender beam, i.e., the shear of the tether was neglected inthis work. The author also presented the coupling interaction
between the tethered cable and ROV’s six degrees of freedom
motion. The boundary condition between the tether and ROV at
the sling point was clearly expounded. The implementation of
torsion effect was presented by carrying out the simulation of theturning maneuver in both clockwise and counter-clockwise
directions.

1) Shabana, A. A., Hussien, H. A., and Escalona, J.
L.:Application of the Absolute Nodal Coordinate
Formulation to Large Rotation and Large Deformation
Problems, ASME J. Mech. Des., Vol.120, pp.188-195, 1998.
2) Suzuki, H., Yoshida, H., Inoue, T.:Evaluation of Methods toEstimate Hydrodynamic Force Coefficients of Underwater
Vehicle Based on CFD, Proc. of the 9th IFAC Conference2013(CAMS2013).
3) Suzuki, H., Tomobe, H., Kuwano, A., Takasu, K., Inoue, T.,
Ishiwata, J.:Numerical Analysis of the Motion of ROV
Applying ANC Method to the Motion of Tethered Cable,
Proc. of the ISOPE Conference, pp.447-453, 2017.
4) Suzuki, H., Tomobe, H., Kuwano, A., Takasu, K., Thant Zin
Htun :Numerical Motion Analysis of ROV applying ANCF
to Tether Cable Considering its Mechanical Property, Proc.
of the ISOPE Conference, pp.365-372, 2018.
5) Berzeri, M., Shabana, A. A.:Development of Simple Models
for the Elastic Forces in the Absolute Nodal Coordinate
Formulation, Journal of Sound and Vibration, Vol.235, No.4,
pp.539-565, 2000.
6) J. Gerstmayr, A.A. Shabana :Analysis of Thin Beams andCables Using the Absolute Nodal Coordinate Formulation,
Nonlinear Dynamics, Vol.45, pp.109, 2006.

International Course of Maritime & Urban Engineering / Abstract of Master Thesis

August 1, 2018

Minimum Propulsion Power Prediction of a Ship under Adverse
Weather Conditions with Dynamics of Diesel Engine and Turbocharger
Taken into Account
Myo Zin Aung
Ship Design Subarea, Department of Naval Architecture and Ocean Engineering
Key Words: Minimum Propulsion Power Requirements, Mean Value Engine Model, Turbocharger, Bulk carrier

1. Introduction(𝑚 + 𝑚𝑦 )𝑣̇ + (𝑚 + 𝑚𝑥 )𝑢𝑟 − 𝑚𝑦 𝑙𝑦 𝑝̇

To reduce the amount of greenhouse gas emissions frominternational shipping and to encourage increasing the efficiencyof the new ships at the design stage, IMO adopted the EEDI
regulations in Annex VI of MARPOL convention. EEDI means
Energy Efficiency Design Index and it is measured in CO2 g/tonnmile. Therefore, higher EEDI value means more CO2 emission
from a ship. EEDI regulation limits the amount of CO2 emission
from ships and this regulation will be gradually tighten phase byphase every 5 years starting from 2015. Since EEDI is proportional
to the ship’s propulsion power and the power is also roughly
proportional to the cube of ship’s speed, ship design speed
reduction is the easiest response to the EEDI requirement.
Although the main purpose of EEDI is to encourage for thedevelopment of more efficient ships, it becomes a limit for
maximum power when expensive new technologies and low
emission fuels are not used. In the near future, the ship owners anddesigners may reduce the installed propulsion power just to comply
with EEDI regulations. Simple reduction of engine power could
result in operational danger under adverse weather conditions.
Therefore, IMO also adopted a guideline for minimum propulsion
power.
To improve this guideline, manoeuvring simulations in adverse
weather conditions are essential for the assessment of navigationalsafety due to propulsion power deficiency. EU project SHOPERA
and JASNAOE projects (Ohiwa and Umeda, 2016) contributed to
discussion on the guideline at the IMO. Better numerical prediction
models with precise modelling of wave drift forces, wind forces
and manoeuvring forces were also developed. However, theirengine models were based on the diagram of engine speed andpower without explicit consideration of engine-turbocharger andso on. Since the minimum propulsion requirements are directly
related to dynamics of diesel engine and turbocharger, more precise
modelling of diesel engine and turbocharger is desirable. Thus, inthis work, the author proposes a new prediction model with a more
precise engine-turbocharger model in the time domain andexecutes manoeuvring simulations of a Handymax bulk carrier
under adverse weather conditions to access its safe navigational
limit.

= 𝑌𝐻 + 𝑌𝑅 + 𝑌𝐴 + 𝑌𝑊

(2)

(𝐼𝑥 + 𝐽𝑥 )𝑝̇ − 𝑚𝑦 𝑙𝑦 𝑣̇ − 𝑚𝑥 𝑙𝑥 𝑢𝑟 + 2𝑑1 𝜙̇ + 𝑑2 𝜙̇|𝜙̇|̅̅̅̅̅𝜙 = 𝐾𝐻 + 𝐾𝑅 + 𝐾𝐴
+ 𝑚𝑔𝐺𝑀(𝐼𝑧 + 𝐽𝑧 )𝑟̇ = 𝑁𝐻 + 𝑁𝑅 + 𝑁𝐴 + 𝑁𝑊

(3)
(4)

2.2 Environmental Forces and Moments
The wind forces and moments are estimated with the Fujiwara’s
empirical formula (Fujiwara 1998). The added wave resistance fordifferent ship speeds are estimated by Maruo’s theory (Maruo
1963). The wave-induced steady sway force, and yaw moment arecalculated by a linear 3D panel method using the Green function
for zero forward speed (Yasukawa et al., 2018). These waveinduced steady forces and moments in irregular waves arecalculated with the JONSWAP wave spectrum by assuming they
are proportional to the square of wave height.

3. Mean Value Engine Model for Propulsion System
The main contribution of this thesis is to include the engineturbocharger and governor model in the manoeuvring simulation.
For the engine-turbocharger modelling, cycle mean value approachis used as shown in Fig. 1 whereas the governor is designed using
PI controller with the integral wind-up action.

2. Manoeuvring Model
Fig. 1 Engine Components and Operating Parameters.

2.1 Equations of Motion
For the maneuvering simulations, 4 DOF equations of motion
are represented as follow using the MMG model and the relevant
standard symbols (Ohiwa and Umeda, 2016):(𝑚 + 𝑚𝑥 )𝑢̇ − (𝑚 + 𝑚𝑦 )𝑣𝑟= 𝑋𝐻 + 𝑋𝑃 + 𝑋𝑅 + 𝑋𝐴 + 𝑋𝑊

(1)The engine-turbocharger model was split into different
subsystems as shown in Fig. 1; the compressor, turbine, intake
manifold, exhaust manifold and engine cylinders. Each part wasmodelled using physical principles, and information was
exchanged between them with the means of flow of mass andenergy.

3.1 Dynamics of Engine Assembly
The dynamics at the crankshaft, turbocharger, intake(scavenging) manifold, exhaust manifold are represented as follow.
Using the Newton’s law, crankshaft and turbocharger dynamics aremodelled as follows (Karlsen 2012, Xiro 2002):
𝜔̇ 𝐸 =𝑄𝐸 −𝑄𝑃 −𝑄𝑓 (𝜔𝐸 )𝜔̇ 𝑇𝐶 =

𝐽𝑡𝑜𝑡𝑎𝑙𝑄𝑇 −𝑄𝐶

(5)

𝐽𝑇𝐶

(6)
Fig. 4 Trajectories of Ship.

Using the first law of thermodynamics, adiabatic formulationof intake and exhaust manifold dynamics are as follows (Karlsen
2012, Guzzella, Onder 2010):𝑝̇𝑖𝑚 =
𝑇̇𝑖𝑚 =
𝑝̇𝑒𝑚 =
̇ =
𝑇𝑒𝑚

𝜅𝑎 𝑅𝑎
𝑉𝑖𝑚

(𝑚̇𝑐 𝑇𝑖 − 𝑚̇𝑎 𝑇𝑖𝑚 )𝑅𝑎 𝑇𝑖𝑚𝑝𝑖𝑚 𝑉𝑖𝑚
𝜅𝑒 𝑅𝑒
𝑉𝑒𝑚

(7)

(𝜅𝑎 𝑚̇𝑐 𝑇𝑖 − 𝜅𝑎 𝑚̇𝑎 𝑇𝑖𝑚 − 𝑇𝑖𝑚 (𝑚̇𝑐 − 𝑚̇𝑎 ))(𝑚̇𝑒 𝑇𝑒 − 𝑚̇𝑡 𝑇𝑒𝑚 )𝑅𝑒 𝑇𝑒𝑚𝑝𝑒𝑚 𝑉𝑒𝑚(8)(9)(𝜅𝑒 𝑚̇𝑒 𝑇𝑒 − 𝜅𝑒 𝑚̇𝑡 𝑇𝑒𝑚 − 𝑇𝑒𝑚 (𝑚̇𝑒 − 𝑚̇𝑡 ))(10)

3.2 Validation with Manufacturer’s Data
The mean value engine model is validated with the steady-state
performance data provided by the engine manufacturer (MAN
Diesel & Turbo 2016) using the 1 DOF manoeuvring model ofHandymax bulk carrier. The engine specified MCR point (design
operating point) is 7930 kW at 137 rpm, which corresponds to theheavily running condition. The validation results show the good
agreement with the manufacturer’s data.

Fig. 2 Interaction of Engine Validation Model Subsystems
Fig. 5 Time Series Data.
9000

5. Conclusions

Engine Power (kW)

8000
7000
6000
5000
4000
3000
2000
1000
0
90

95

100

105

110

115

120

Engine Speed (RPM)
Manufacturer

125

130

135

140

Simulation

Fig. 3 Heavily Running Engine Design Curve (Propeller Curve)

4. Manoeuvring Simulation in Adverse Weather
The ship is initially assumed to be in following sea condition andit tries to turn to head sea condition using PID autopilot. In this
scenario, the ship is able to make turning to the head sea up to
Beaufort Scale (BF) No. 10 condition. The final steady-state ship
speeds in the head sea are 8.5 knots, 3.0 knots and 0.7 knots for BF
No. 8, BF No. 9 and BF No. 10, respectively. Therefore, we may
assume that the ship with this particular installed power is safe up
to the BF No. 10 weather condition for keeping head sea conditions.
The ship trajectories and time series data are shown in Fig.4 andFig.5, respectively.

A new simulation model for the prediction of minimum
propulsion power is proposed using the mean value engine
modelling approach for the propulsion system. The proposedengine model is validated with the steady-state performance data
provided by the engine’s manufacturer. Manoeuvring simulations
are executed for a Handymax bulk carrier to observe the dynamics
of propulsion engine and to access the safe navigational limit under
different weather conditions. The proposed simulation model canscope realistic responses of engine power and rpm which areimportant when considering the power limit lines defined by thesevariables. Moreover, the new simulation model is applicable forwider type of manoeuvring scenario since the engine or propeller
rpm can be controlled using the governor.

References
1) Ohiwa, S., Umeda N., “Probabilistic aspect on minimum
propulsion power requirement issue under adverse weather
conditions,” Proc. of DfS 2016, Hamburg, pp. 127-131, 2016.
2) Karlsen, T. A. “On modelling of a ship propulsion system forcontrol purpose,” Master Thesis, NTNU, Norway, 2012.
3) MAN Diesel & Turbo, Copenhagen, Denmark, “MAN B&W
S40ME-B9.5-TII Project Guide,” July, 2016.

International Course of Maritime & Urban Engineering / Abstract of Master Thesis

July 27, 2018

Structural Analysis of a Floating Vertical Axis Wind Turbine
Ahmad Adilah
4th Laboratory, Department of Naval Architecture and Ocean Engineering
Key Words: Floating vertical axis wind turbine, Gyroscopic effect, Vertical bending moment

1. Introduction

2.2 Equations of motion

In recent days floating Vertical Axis Wind Turbine (VAWT)
are attracting researchers’ attention. This is because it is moresuitable for floating offshore wind turbine than the horizontal ones.
It has a low center of gravity, simplicity of the installation andinsensitivity to yaw and pitch1).
There are a few researches that analyze the motion of floating
VAWT. Meanwhile, almost no researchers pay attention to thestructural analysis of floating VAWT. This analysis becomes more
important considering the rapid development of floating verticalaxis wind turbine with larger structure in the future.

The equations of translational and rotational motion of thestructure are given as follows.

2. Research methodology
2.1 Two-step method
A practical so-called two-step method is used in which a rigid
body model is used to obtain the load distribution, and the
distribution of the load is given to a structural model adopted in thesecond step for stress/deformation analysis.

𝐹𝑚𝑜𝑟𝑖𝑠𝑜𝑛 + 𝐹𝑎 + 𝐹𝑏 + 𝐹𝑔 + 𝐹𝑚𝑎𝑔𝑛𝑢𝑠 + 𝐹𝑚𝑜𝑜𝑟𝑖𝑛𝑔 = 𝑚𝑑𝑑𝑡(𝐼𝜔) + 𝜔 × (𝐼𝜔) = ∑ 𝑀𝑖

𝑑𝑣
𝑑𝑡

(2)

(3)

where, m is the total mass of the structure, v is the relative velocity
vector of the structure, I is the inertial moment matrix and ω is theangular velocity vector of the structure in body fixed coordinate
system. Fa is the aerodynamic force that acting on the structure, Fb
and Fg is the buoyancy and gravity force of the structure. The
rotational moment caused by the total forces are written as Mi.
Quaternion is used to rotationally transforming the coordinate
system.
2.3 Model of analysis
In this research, a spar structure is used as a floating structureof the turbine to accommodate 1 MW vertical axis wind turbine.
The SPAR structure consisting of various diameter. Table 1 shows
the structural configuration of spar.

Fig. 1 Block diagram of two-step method
A rigid body numerical simulation procedure for the floating
VAWT is developed and compared against the data fromexperiment result2). The simulation can account for thehydrodynamic loads, aerodynamic loads, rigid body motions,
gyroscopic effects, Magnus effect, mooring loads, hydrostatic andgravity loads. While evaluate the hydrodynamic load, SPAR
structure is divided into short cylinders for numerical evaluation.
Magnus effect is considered when the SPAR rotates in marine
current and/or wave. The direction of the Magnus force isperpendicular to the flow and the rotational axis. The magnitude of
the force per unit length is,
𝐹𝑚𝑎𝑔𝑛𝑢𝑠 = 𝜌𝑤 𝑢𝑟 Γ

Table 1 Structural configuration of spar
Draft
48.3 mUpper diameter
3.3 m
Bottom diameter
5.0 m
Shell thickness
24.2 mm
For evaluating the structural behavior, the spar structure isdivided into beam elements along the axis. Figure 2b shows thesubdivision with beam elements for the spar structure.

(1)

where ur is the relative flow speed component perpendicular tothe cylinder axis. Γ is the circulation around the cylinder, Γ =(πd2/4) · 2ωrot. Here, ωrot is the rotational speed of cylinder and dis the diameter of cylinder.
Further, the structural analysis is performed by using the load
distribution that obtained from the rigid body analysis. The
external force and inertial force of the structure should be balanced
while perform the structural analysis. Flexibility of the structure
should be considered in this analysis. The load distribution is
integrated to obtain shear force, vertical bending moment andfinally bending stress distribution of the structure.

(a)

(b)

Fig. 2 (a) Floating axis wind turbine and (b) beam element
subdivision of the spar structure

3. Result and Discussion
3.1 Rigid body motion
The rigid body motion is simulated in regular wave condition
with wave height H is 0.2 m and wave period T varying between 5s to 12 s.

Fig. 5 Bending stress distribution along the spar axis. T = 12 s
From the result we can know that the maximum bending stressoccurs in element near still-water level, this because of the small
cross-section area in the transition part compared to other parts, andlarge vertical bending moment there.

Fig. 3 Trajectory in xy-plane of spar axis in the regular wave.
T = 12 s
Figure 3 shows the trajectory in xy-plane of the spar axis. The
simulation results show that Magnus effect and gyroscopic effects
occur in the motion of the structure. The Magnus effect cause theturbine experience the drift motion and the gyroscopic effects
contributes to the stabilization of the turbine.
3.1 Structural analysis
Structural analysis is performed under the static (aerodynamic)
and dynamic (wave) load condition. Based on the load distribution,
we can know the shear force and vertical bending moment
distribution along the spar axis. Figures 4 and 5 show the
distribution of shear force and vertical bending moment for regularwave with period T = 12 s.

Fig. 6 RAO stress at beam element transition part
By plotting the response amplitude operator (RAO) stress at
the element near still-water level for all wave frequency range we
can know that the maximum bending stress is about 28.2 MPa perunit wave amplitude. For comparison, the pitch acceleration RAO
is presented in Fig. 7. The two curves in Figs. 6 and 7 resemble
each other because the inertial load is causing the structural
response, which results in the larger stress in higher frequency.

Fig. 7 RAO of angular acceleration of spar structure
Fig. 4 Shear force distribution along the spar axis. T = 12 s
For static aerodynamic load case, the maximum bending
stress in the transition part is measured to be 123 MPa. It is almost
exceeding the yield stress which is 175 MPa, therefore some designmodification is necessary.

4. Conclusions

Fig. 5 Vertical bending moment distribution along the spar axis.
T = 12 s

From rigid body simulation, it is found that the turbine
experiences the Magnus effect and gyroscopic effect. However
gyroscopic effect on turbine enhances the stability of the structure.
The maximum bending stress are found about 28.2 MPa for
dynamic load case and 123 MPa for static load case.

References
The zero point in x-axis represent the top of the structures andthe bottom position of the structure is 96.5 m. After verticalbending moment is known, the bending stress can be obtained by
dividing the vertical bending moment with sectional modulus of
the spar.

1) H. J. Sutherland et al.: A Retrospective of VAWT Technology,
SAND2012-0304, 2015.
2) H. Akimoto.: Floating vertical axis wind turbine simulator
SimFAWT (not published)International Course of Maritime & Urban Engineering / Abstract of Master Thesis

July 28, 2017

Numerical Calculation of Ship Hydroelasticity
Using Enhanced Unified Theory
ZHOU Dali

Ocean-Space Development Laboratory, Dept. of Naval Architecture & Ocean Engineering
Key Words: Ship hydroelasticity, Enhanced unified theory, 3D effect, Mode function, Green function method

1. Introduction
As the main means of transportation, the ship is still keeping
popularity in the trading industry over centuries. Since the needfor cargo transportation capability is gradually increasing, thedimension of ships becomes larger subsequently. Thus, ships
tend to be flexible with increasing the ship length. The
deformation of a body, the effect of structure stiffness, and thefluid structure interaction are of much importance in theseakeeping analysis.
Various analysis methods for the global hydroelastic behavior
have been developed by many researchers. To name a few
among them, we may write as follows. Bishop and Price1)
presented a fundamental concept of the hydroelastic problem
using the strip theory. Jensen and Dogliani calculated nonlinear
springing load using a nonlinear strip theory. Malenica et al.
proposed a method combining the frequency-domain 3D
boundary integral equation method (BIEM) with the finite
element method (FEM) for a non-uniform beam, results of whichwere compared successfully with experimental data.
Considering these developments so far and a practicalrequirement that the global hydroelastic response must be
computed with reasonable accuracy and less computation time,
we need to establish a practical computer code for shiphydroelasticity problems. In this paper, Enhanced Unified
Theory (EUT) 2) is applied to the hydrodynamic analysis. For the
structural analysis, it is common to define a set of natural modeshapes, which correspond to actual elastic deflections of a bodyin a specified physical context. The structural deflection isrepresented by a superposition of these natural mode shapes.
Compared to the analysis for a rigid body, an importantdifference in the hydroelastic analysis exists in the bodyboundary condition which represents the normal velocity on the
body surface induced by the structural deflection. Computed
results for hydrodynamic forces and wave-induced motions arevalidated for the zero-speed case through comparison with moreexact results by the 3D Green function method, and then the
results for the forward-speed case are presented and discussed.

2.

Structural Analysis

2.1 Structural Problem
In the structural analysis, it is common to define a set oforthogonal mode functions3), which correspond to the actual
elastic deflections of a body in a specified physical context.
Generally this is difficult when the mode shapes are affected by
the hydrodynamic pressure field and cannot be specified inadvance. This difficulty can be avoided if the structuraldeflection is represented instead by a superposition of simpler
mathematical mode shapes which are sufficiently general andcomplete to represent the physical motion.

2.2 Mode function
In the modal superposition method, the total displacement of aship H ( x , t ) is expressed by superposition of dry natural mode
functions as follows:N

H ( ξ , t ) = ∑ ξ j (t )h j ( ξ )

(1)

j =1

where x = ( x, y, z ) is the position vector, index j denotes themode number, and ξ j (t) is the time dependent amplitude of thej-th mode of motion. The mode number j=1~6 and j=7~Ncorrespond to the rigid- and elastic-body motions, respectively,
and N is the total number of modes to be considered.
As an example, in the case of vertical bending the components
of normalized mode-function vector h j ( x ) can be written as=
h1j 0,=
h2j 0,=
h3j w j(2)where h1j has a nonzero component, but due to order estimation
with slender ship assumption, it can be regarded as negligibly
small. On the other hand, w j (x) is the vertical displacement in
the j-th mode of motion. If the free-free beam modes of a
uniform beam are adopted for w j (x), they are given analyticallyand their shapes can be given as shown in Fig. 1.

Fig. 1 Vertical displacement in the j-th mode of motion

3.

Hydrodynamic Analysis

3.1 Computation method
Let us consider a ship advancing with constant speed U andundergoing small-amplitude harmonic motions with angularr
frequency ω in deep water. A steady reference frame is chosenwith the x-axis pointing to the direction of forward motion andthe z-axis pointing downward. The plane progressive wave isincident upon the ship with incidence angle χ , in whichχ = 180 deg is defined as head wave.
Assuming the flow inviscid with irrotational motion, thevelocity potential is expressed as=φΦ=−Ux + Re[φ ( x, y, z )eiωt ](3)ga{φ0 ( x, y, z ) + φ7 ( x, y, z )} + iω ∑ X jφ j ( x, y, z )
iω0
j =1

(4)

6

φ0 = exp[−k0 z − ik0 ( xcoscc+ ysin )]ω=ω0 − k0Ucosc ,k0 =ω0 / g2Here φ0 denotes the incident-wave potential, φ7 the scatteredpotential, and φ j the radiation potential of the j-the mode withcomplex amplitude X j . Also a , ω0 , k0 denote the amplitude,
circular frequency, wave number of the incident wave; g theacceleration due to gravity. Besides, the disturbance due to thesteady forward motion is neglected as higher order.

As an example, computed results of the wave exciting force atzero speed for elastic motion (mode number equal to 7, 8 and 9)
are shown in Fig. 2, from which we can see god agreement withdifferent methods. In Fig. 3, the results of the amplitude ofwave-induced ship motions for both rigid and elastic bodies arepresented. Good agreement also could be found except for theelastic mode j = 8 , for which the result by the frequency-domain
3D Green function method must be something wrong, but thereason of the problem could not be detected up to now.

3.2 Hydrodynamics force
Integrating the pressure over the mean wetted surface of a ship,
and using Tuck's theorem, the hydrodynamic force in theradiation problem2) acting in the i-th direction can be expressedasFi =−(iω ) 2 ∑[ Aij + Bij / iω ] X j

(5)

U ~
U
M i ){j j + jj }ds
iω
iω
~
U ~
− ρ ∫ dxC j ( x) ∫ ( N i − M i ){j3 − j3*}ds
L
SH
iω

(6)

~

− ρ ∫ dx ∫ ( N i −Aij + Bij / iω =L

SH

Fig. 2 Wave exciting forces for elastic modes, Fn=0.0.

where Aij and Bij are the added-mass and damping
coefficients in the i-th direction due to the j-th mode of motion,and SH denotes the sectional contour below the free surface atstation x. The first line on the right-hand side of Eq. (6) isessentially the same as that in the strip theory and the second line
represents the 3D-effect which can be computed in terms ofhomogeneous component in EUT.
In the diffraction problem2), the hydrodynamic pressure isgiven from the sum of the incident-wave and scattered potentials.
The symmetric part of the diffraction force can be expressed as~

E j = − ∫∫ Pd N j dSSHρ ga ∫ dxC7s ( x)eix ∫
L

SH

~

N j {y 7s + e − k0 z cos(k0 ysinc )}eix ds(7)Fig. 3 Ship motion amplitudes in rigid and elastic modes.

4. Numerical calculation

5. Conclusions

Numerical computations are performed for a 2-meter
modified Wigley model. The principal dimensions are shown in
Table 1.

In this research, the hydroelastic forces and the elastic
deformation of a ship are computed by using the enhanced
unified theory (EUT) combined with the modal superposition
method for representing the elastic deformation of a body, andthe normal velocity on the body surface due to elastic
deformation is computed from the mode functions obtained by
the uniform-beam approximation for a ship.
Through the comparison between EUT, Green-function
method and Rankine panel method, it was confirmed that theresults by EUT can provide reasonable accuracy and agreement
and hence can be a practical tool. By using EUT, plenty of time
could be saved due to very fast computation while enoughaccuracy still can be kept.

Table 1 Principal dimensions of modified Wigley model.
Length L
Breadth B
Draft dDisplacement volume ∇Height of the center of gravity OG
Gyrational radius in pitch κ yy

2.000 m0.300 m0.125 m0.04205 m3
0.0404 m
0.456 m

4.1 Result
The numerical computations were launched for the zero speedcase for validation and then for the forward speed case.
Computed results by EUT for the radiation force, wave exciting
and ship motions are compared with the frequency-domain 3D
Green function method for the zero-speed case, the time-domain
Green function method, and the Rankine panel method for theforward speed case.

References
1)

2)
3)

Bishop. R.E.D. and Price. W.G.: Hydroelasticity of Ships,
Cambridge University Press, 1979
Kashiwagi, M.: Transactions of West-Japan Society ofNaval Architects, No.89, pp.77-89, 1995
Newman. J.N.: Applied Ocean Research, Vol. 16, pp. 47-59,
1994.

International Course of Maritime & Urban Engineering

July 27, 2018

Ultimate Strength of a Continuous Stiffened Panel under Combined
Transverse Thrust and Inplane Shear
Aung Kyaw Htet
Ship Structural Integrity Subarea, Dept. of Naval Architecture and Ocean Engineering
Key Words: Stiffened panel, Ultimate strength, Buckling, Transverse thrust, Shear

1. Introduction
A ship structure basically consists of thin plates on which a
number of stiffeners are provided to increase their strength andstiffness. The side shell plating of bulk carrier, Fig.1, is subjectedto the thrust in the shorter edge direction due to horizontal hull
girder bending, and to the in-plane shear due to torsional moment.
It has been recognized that larger plate and side-frame scantlings
are required for the side shell plating ever since the strength
requirements for the horizontal bending and torsional moment
were introduced to the Common Structural Rule for Bulk-Carrier
and Oil-Taker (CSR-BC&OT) issued by IACS (International
Association of Classification Societies) in 2006.

With regard to the in-plane displacements, all edges are
assumed to be straight during whole collapse process consideringinteraction of adjacent panels. The periodic boundary condition is
imposed on the in-plane displacements along the edges. In theprevious studies of buckling behavior under in-plane shear,uniformly distributed shear stresses have been assumed. But inthe proposed model, the shear stress distribution can change along
the edges as the buckling deformation develops. This enables us
to simulate the capacity reduction after shear buckling. The two
representative nodal displacement UR and VR are defined at thecorner node A. The transverse thrust is applied by VR and thein-plane shear by UR. Combined thrust and shear can be applied
under any stress ratio and in either load or displacement control.
D

A

VR

y

UR

B

C

Fig. 1 Side shell plating of Bulk Carrier under combined loads.
Fig. 2 Model for analysis using periodic boundary condition.
In this study, the buckling/ultimate strength of a continuous
stiffened panel subjected to combined transverse thrust andin-plane shear is investigated using nonlinear FEM. To consider
the effect of continuity of the plating, the periodic boundary
condition developed in the laboratory is used. The background of
the plate-capacity and stiffener-capacity formulas given in the
IACS/CSR is clarified. The simplified ultimate-strength formuladeveloped by Fujikubo et al (OU formula)1) is also considered.
The ultimate strength is compared between FEM and the CSR
and OU formulas and applicability of the formulas are examined.

The initial deflection of local buckling mode is expressed asmwo   A0 sini 0When buckling/plastic collapse analysis is performed on acontinuous stiffened panel, the minimum required region that canaccount for all possible bucking modes of plate and stiffener
should be properly selected using relevant boundary conditions.
Here we employ the periodic boundary conditions as shown in
Fig. 2. When the rectangular plate of length a and width b isanalyzed, the region 2a x 2b including two stiffeners is taken as amodel for analysis. The condition that all displacements are sameat the longitudinal edges AD and BC, and also at the transverse
edges AB and DC, is imposed, assuming the same out-of-plane
deformation is exerted periodically. This periodic boundary
condition allows the buckling mode of either even or odd numberof half waves to take place.

(1)

where the first term is the buckling mode under thrust and thesecond term the buckling mode under in-plane shear. The lattermode is calculated by FEM analysis. The initial deflection of
stiffener is expressed by the tripping modes and the flexural modeasvo =Co2. FEM Collapse Analysisi x  ysin
 A0 f s ( x, y )
a
b

x
z
x ,
woall  Bosin(
)
sin(
)
a
hs
a

(2)

where Bo=Co=0.001×a. The amplitude A0 is determined so thatthe maximum initial deflection is equal to the specified value. The
dimensions and material properties are summarized in Table 1.
Table 1 Dimensions and material properties.
Plate width (b)Aspect ratio (a/b)
Plate thickness (tp)Angle-bar stiffener
Young’s modulus (E)Poisson’s ratio
Yield stress (Y)

800 mm2, 4, 69, 16, 22, 28 mm
h×bf×tw/tf = 300×90×10/16
205,800 MPa
0.3
313.6MPa

3.

Simplified Ultimate Strength Formula

3.1 CSR formula
In the CSR, the plate capacity, σup, and the stiffener capacity,
σus, are calculated by different formulas, and the ultimate strengthof a stiffened panel is determined by the smaller value of σup andσus.The plate capacity, σup, is given by the plastic correction of theelastic buckling strength of a bare plate. The stiffener capacity,
σus, under the action of transverse thrust and in-plane shear is
calculated by a concept of an equivalent lateral loads as follows:(a)in-plane shear. The stiffener failure, as assumed in the CSR
formula, has not been observed for the stiffened panels underconsideration. As shown in Fig. 4, the OU formula gives good
agreement with the FEM results slightly on the safe side. On theother hand, the CSR formula gives reasonable estimate when thein-plane shear is dominating, but it becomes over-conservative
when the transverse thrust is dominating especially for large
aspect ratio and large plate thickness. This is because the stiffener
capacity equation gives smaller capacity than the plate capacity as
the aspect ratio, and hence the span increases. Such stiffener
failure is unreasonable referring to the FEM results. Some
improvement is needed for the CSR formula. Although the results
are omitted due to page limitation, the OU formula has also aroom for improvement especially for very thin plate.
a/b=6 ,tp =28mm

0.7

(b)

FEM result

0.6

CSR formulaFig. 3 Equivalent lateral load concept in CSR-BC&OT.
For instance, when transverse thrust acts, transverse bending
moment is exerted as shown in Fig. 3(a). The equivalent lateralload qy in the Fig. 3(b) is defined as a load per unit length thatcauses the same transverse bending moment. The stiffener
capacity is evaluated by the longitudinal bending stress of thestiffener caused by the qy distributing along the stiffener span.
The applicability of this strength model is examined through thecomparison of the collapse mode and strength with the FEM
results.

τ/σY

0.5

OU formula

0.4
0.3
0.2
0.1
0
0

0.1

0.2

0.3

0.4

0.5

σ/σY

0.6

0.7

0.8

0.9

1

Fig. 4 Ultimate Strength interaction relationship of continuous
stiffened panels under combined transverse and shear.

3.2 OU formula
Fujikubo et al1) derived the ultimate strength formula of acontinuous unstiffened panel under transverse thrust as u b  2.4 1.4   b   0.06 0.6 


 1
 2
 y a    2   a   
 

(3)

where is the slenderness parameter given by



b y
t E

(4)

The parameter  has a physical meaning of a square root of theyield strength to elastic buckling strength ratio. In the stiffened
panel, the elastic buckling strength of local plate is increased by
the torsional stiffness of stiffener. Using the increasing rate y,derived analytically1), the equivalent slenderness ratio with theeffect of stiffener βs is defined as

s 

b y
t s E(5)

Substituting the s into Eq. (3) instead of  , the ultimate strengthof stiffened panel under transverse thrust is calculated.

4. Result and Discussion
The ultimate strength of stiffened panel obtained by FEM iscompared with that predicted by the CSR and OU formulas. Fig.
4 shows the ultimate strength interaction curve obtained for theplate with the aspect ratio of 6. Fig.5 shows the deflection modeat the Post-ultimate behaviour for two calculation cases.
According to the FEM, the buckling collapse occurs in thelocal plate bucking mode under combined transverse thrust and(a) a/b=2; tp=16mm; σ :τ=1:2(b) a/b=4; tp=16mm; σ :τ=1:0
Fig. 5 Buckling deformation under combined transverse and shear

5. Conclusions
For larger aspect ratio, especially for thicker plate, theultimate strength obtained by CSR formula is remarkably lower
than that of FEM results. This appears due to improper concept ofequivalent lateral loads. This error may significantly influence thedesign scantling of the side shell plating of Bulk Carrier.
OU formula considering the effect of the stiffener gives
reasonable estimate of the ultimate strength. But for the strength
assessment of very thin plate, such as severely corroded plate,
improvement should be made. This remains as a future work.

Reference
1) Fujikubo, M., Yao, T., Mohammad R K, Harada, M., andYanagihara, D.: Estimation of Ultimate Strength ofContinuous Stiffened Panel under Combined Transverse
Thrust and Lateral Pressure Part2: Continuous Stiffened Panel,
Marine Structure, pp411-427, 2005.

International Course of Maritime & Urban Engineering / Abstract of Master Thesis

July 12, 2018

Adaptive Mitigation Approach for Tsunami Ascending to
Meandering River
Rajuli Amra
Coastal Engineering Laboratory, Department of Civil Engineering
Key Words: Intruding tsunami, meandering river, adaptive mitigation

1. Introduction
Previous tsunami disasters revealed that tsunami could penetrate
up far upstream of the river. It brought serious damage to river
structures and flooded the residential area1). Tsunami propagating
river exhibited very different actions from those observed on an
open coast. Unfortunately, the existing tsunami protections wereless effective in mitigating this typical tsunami 1).
Previous research highlighted that morphological features, tideand river discharge were some dominant factors in influencing
propagation behavior 1) 2). However, the effect of river planform
were not significantly addressed, although some research notedthat the meandering part of the river might play the significant
role1). Therefore, this study focuses on investigating the role ofriver geometry to tsunami characteristic. This research also
proposes the new mitigation approach with respect to propagation
behaviors and limitation of existing protection concept.
The proposed mitigation approach carried the concept of “extra
room for the river”. The result showed that the new approach was
effective in countering the overflow, attenuating the wave height
distribution along river, delaying tsunami arrival and reducing thecapital cost.

2. Methodology
2.1 Scope and limitations
The investigation was conducted numerically by using 2D
shallow water equation of COMCOT model, in which the tsunami
was generated by 1m, 2m and 3m amplitude of solitary wave3). In
order to clearly understand the effect of river planform, then theinfluence of other components, e.g. river discharge, tide andmorphological features, were not considered. The sedimenttransport process was also not discussed. In addition, to isolate theoverland flow from the shore, the very high wall was establishedand the wall structures both at the river banks and the shore wereassumed stiff.
The tsunami behaviors were addressed by analyzing the
maximum tsunami attributes, i.e. wave height and resultant
velocity, along and lateral-river, tsunami arrival, overflow patterns
and range of inundation. After that, the comparison was drawn
between the two conditions, before and after proposed mitigation
approach was applied.
2.2 River planform
Sine-generated curve function was employed in modellingsymmetric meandering geometry, while the river bed elevation wasmodelled by empirical function of the transverse bed. Threeplanforms configuration were used for 10 km river-length torepresent the complexity of natural river, which are low (θ0=300),
moderate (θ0=600) and high-sinuosity meandering planform(θ0=900).

( )=

(2

⁄ )

( )=∫

( )

(2)

( )=∫

( )

(3)

(1)

Fig. 1 Sketch of simulation area and notation for equations
2.3 Simulation Scenario
To understand the propagation behavior, three scenarios wereset which were CS1, CS2 and CMS. Both CS2 and CMS scenarios
were also used for investigating the effectivity of proposed
mitigation approach, along with PM1 and PM2 scenarios. Each ofcases was simulated for 60 minutes condition. The observation
points were stationed per 1 km along river and each 20 m lateralriver at the bend apex. The river slope was assumed flat with 10
km length and 400 m width.
Table 1 Simulation Scenario
No
Scen.
Description
1
CS1
No river embankment installed
2
CS2
River wall height: 1.5 m uniform
3
CMS
River wall height: max. wave height of CS2
4
PM1
River wall height: varied along river
5
PM2
PM1 combined with “extra room” concept

3. Result and Discussion
3.1 Propagation characteristics
The simulation results described that the tsunami height along
river exponentially decreased to the upstream distance. This pattern
was qualitatively similar to all meandering level, even though there
were a slight difference particularly at the inflection point. In
addition, propagated wave height comparison among low,
moderate and high-sinuosity planform also showed a small
difference. For 1m, 2m and 3m initial tsunami wave height, it was
known that the moderate-sinuosity planform dissipated tsunami
better, compared to two other planforms. Yet, the difference wasonly by about 2%. The wave transmission ratio were correlated

well with the preserved energy. For 1m initial tsunami case, whenoverflow was not occurred, wave transmission ratio for all
sinuosity level at the most upstream point was bigger than 85%,
while it was only about 20% when overflow occurred. Surprisingly,
the dissipation rate at the most upstream point increased
significantly for case 2m and 3m initial tsunami height.
At the bend apex, the lateral distribution of tsunami heightmatched the previous study finding, in which the wave height at
the outer bank tended to be higher, compared to the inner bank. In
the contrary, the higher velocity distributions were at the innerbank. The interesting patterns were found at the case of 2m and 3m
initial wave height for no overflow condition, in which the
maximum lateral velocity distribution indicated the flow separation
might occurred at the bend apex.
3.2 Overflow and inundation pattern
The overflow pattern was depended on the initial tsunami
amplitude and correlated well with propagation characteristics, inwhich the propagating tsunami tended to easily overtop theembankment at the downstream river, as shown in Fig. 2. The CS2
results showed that the northern part of the river were inundated
greater, while longest overflow occurred at the southern part of theriver.
The result also described that the higher inundation were more
dominant at the concave banks. For high sinuosity planform, thisbehavior was more significant, compared to other two otherplanform. This pattern might relate to centrifugal force that appliedin curved path.

about 57% of the initial construction cost. Still, both uniformly or
gradually elevated wall will increase the wave height at the most
upstream point and decrease dissipation rate by 5%. Indeed, forPM1 scenario, the small part of area that was classified as densely
inhabited zone was inundated.
Another alternative was increasing the river capacity byrelocating the embankment at the downstream river. In the PM2
scenario, both side of embankment up to 1.5 river km, which crosssectional to the centerline, was relocated about 350m away tocreate 100ha “extra room”. Combination of “extra room” concept
and gradual embankment positively decreased wave distribution in
the river, compared to previously mentioned approaches. As shownin Fig.3, the most drastic wave reduction was at 1km point fromthe river mouth, in which the wave height has decreased by about
40%. While it was only about 20%, 10% and 9% for CS2, CM andPM1 cases consecutively. The arrival time of tsunami at the most
upstream point was also decelerated. In addition, the riverside area
was not inundated.

Fig. 3 Comparison on wave height and arrival time for θ0=600,
lines indicated wave height, markers showed arrival time
Several variation of PM2 scenario revealed that the effectivityof the proposed concept was constrained by “extra room” capacity,enlargement orientation, i.e. lateral or along river expanding, andthe starting point of river embankment relocation. Since thepropagation behavior was, almost, similar in all sinuosity level,
then the proposed concept is assumed applicable to any
meandering river planform.

Fig. 2 Sample of overflow distribution for CS2 scenarios
3.3 Effectivity of proposed mitigation approach
Elevating the river embankment as CMS scenario proved
effective in countering the overflow along the river. However,preserving the propagated wave in the river caused the increasingof wave height in the river. The increasing wave height followedby the velocity and ultimately raised the force distribution. It could
endanger the in-river structure, e.g. bridge and jetty, and indeed theriver embankment itself, particularly at concave banks. The
elevated wall also accelerated tsunami arrival time at the most
upstream point. The typical characteristic applied to all sinuosity
level of the river.
For the case of 2m initial tsunami height, raising theembankment from 1.5m to 3.25m uniformly along river boosted
the initial cost about 117%. The cost could be reduced if theembankment was elevated gradually along the river, based on
maximum wave height distribution. This method only increased

4. Conclusions
Generally, it could be concluded that different meanderinglevel of river respond in similar mechanism to the intruding
tsunami wave. However, the higher sinuosity planform tends to beinundated greater at the bend apex section. In addition, propagation
patterns tends to follow the river shape. Finally, the proposed
mitigation concept proved more effective in attenuating wave
height distribution and decelerating tsunami arrival time.

References
1) Tanaka H: Study on the relation of river morphology andtsunami propagation in rivers, Ocean Dynamic, Springer, 2014.
2) Elena T: Tsunami Penetration in Tidal Rivers, with
Observations of the Chile 2015 Tsunami in Rivers in Japan,
Pure and Applied Geophysics, 2016.
3) Philip L.F. Liu, Harry Yeh and Synolakis C: Advanced
Numerical Models for Simulating Tsunami Waves and Runup,
W Scientific, 2008.

International Course of Maritime & Urban Engineering / Abstract of Master Thesis

July 27, 2018

S-PIV (Stereo-Particle Image Velocimetry) Measurement for a
Self-propelled Ship in Regular Head Waves with Different Amplitudes
Mwangi Benson Oyunge
Hull Form Design Sub-area, Department of Naval Architecture and Ocean Engineering

Key Words: KVLCC2, EFD, Particle Image Velocimetry, Flow Field Measurement, Wave Amplitude Dependency

1. Introduction
The self-propulsion performance in regular head waves of
KVLCC2 model including ship motion were investigated by Wu
et al1). The CFD computation results were compared withexperiment. The CFD could predict the detailed thrust fluctuationin one period including higher harmonics. But the computationand experiments were done for only one wave height. The flow
field around self-propelled ship model including nominal velocity
field for without propeller condition were measured andcompared with CFD computation. In the Study of Wu, therelation between thrust time history and flow field were discussed
and the vertical relative motion between propeller disk and thebilge vortices are very important to explain the thrust fluctuation.
So, it seems that the wave amplitude has a lot of effect onthrust fluctuation and averaged wake. In this study, the flow fieldaround stern is measured using the SPIV as well as thrust, torque,towing force and 3DOF motion to investigate how the velocity
distribution changes with the wave amplitude and its effect on theamount of thrust delivered by the propeller.

Fig. 1 KVLCC2 Tanker model and coordinate System.

2. Experiment Overview and Test Conditions
2.1 Experiment Overview
The experiments were conducted at the Osaka University
towing tank using a fully-loaded 3.2m KVLCC2 tanker model inlong regular head waves with different wave amplitudes at adesign Froude number, Fr=0.142 (0.795m/s). The motions andforce measurement experiment was done first for with andwithout propeller condition while PIV experiment was done later
for without propeller condition. The flow field at AP section wasmeasured because the laser reflection effect was large forpropeller plane and the flow field is similar as propeller plane toconsider the thrust fluctuation from boundary layer assumption.
For the self-propulsion experiment, 4 blades, 0.986m diameter
MOERI propeller model was used. The ship model particulars are
as shown in Table1 and Fig.1 shows the ship model and thecoordinate system used in this experiment. The propeller model
used is as shown in Fig 2.
Table 1 KVLCC2 Tanker model particulars.

Fig. 2 MOERI Propeller model.

2.2 Test Conditions
For this experiment, three waves of λ/L=1.1, 1.6 and 2.0 werechosen because the large relative motion between propeller disk
and low velocity region due to bilge vortices was observed inprevious study for long waves. Various wave amplitudes i.e.
A=0.75, 1.5, 3, 4, 4.3cm were used for each wavelength. Ship
motion in this experiment is 3DOF (degree of freedom) in which
the ship was free to heave, pitch and surge. These motions were
measured by potentiometers. In self-propulsion experiment, thrust
and torque were measured by a motor-driven dynamometer
connected to the propeller by a shaft. Force in x-direction wasmeasured by a dynamometer which was attached to the heavingrod. The incident wave data was measured at the wave gauge thatwas installed 3.52m (1st experiment) and 3.314m (2nd experiment)
upstream of FP.
The flow field at AP was measured by the S-PIV system for thesame conditions as above. For larger wave amplitudes the PIV
measurement was done in two vertical positions with an overlap

so as to capture the whole flow field. The average number ofimages collected per phase is 180 so as to get sufficient data forimage analyses to get the phase averaged flow field.

wave amplitude whereas the inflow velocity into the propeller
plane increases with increase in wave amplitude as shown in Fig.
6.

3. Results and Analysis
The motions and force measurement results were analyzed by
Fourier series expansion and the data was reconstructed byignoring the higher harmonic than 7th harmonic to understand thefluctuation shape in one period. So, for thrust and torque a lot ofnoise due to vibration of connection joint between motor andself-propulsion dynamometer was eliminated. It was noted thatthe 1st, 2nd and 3rd harmonic components are responsible for thethrust oscillations with a quickly increasing upward slope and aslowly decreasing downward slope. Similar observation wasmade by Emel et al2) with a conclusion that the cause could be therelative motions between bilge vortices and ship motions andpressure gradient between inner and outer boundary layers.
From the Fourier analysis, the harmonic components
responsible for the original trends of motions and forces were
determined. For motions, it was found that 1st harmonics isenough to achieve the original motion trend whereas for thrust
and torque we need up to 3rd harmonic components. From theheave time history in Fig. 3, its amplitudes of oscillation increases
with increase in wave amplitude and Fig. 4 shows that the 1st
harmonic amplitudes of heave increase linearly with wave
amplitude. The same trends are shown by pitch motion.

Fig. 3 Heave time history for wave peak at propeller plane.

Fig. 5Thrust time history with mean values in λ/L=1.1 fordifferent wave amplitudes.

Fig. 6

Velocity distribution for phase 90 at AP for different
wave amplitudes.

4. Conclusions
The velocity distribution at AP was analyzed and discussedwith respect to the thrust fluctuation in waves. The trend can be
explained although the velocity field was measured for without
propeller condition. The motions and thrust time histories inwaves of different amplitudes were computed by CFD and the
CFD results show fairly good prediction as in the presentation.

References

Fig. 4 Heave 1st Harmonic amplitudes.
Fig. 5 shows that mean thrust decreases slightly with increase in

1) Ping-Chen et al: Proc. of 30th Symposium on Naval
Hydrodynamics (Tasmania, Australia), 2014
2) Emel et al: JASNAOE, pp.255-258,2015

International Course of Maritime & Urban Engineering / Abstract of Master Thesis

August 1, 2018

Verification & Validation of a Cellular Automaton Simulating UnderFilm Corrosion
Ivan Kristov,Ocean Material Engineering Laboratory,
Department of Naval Architecture & Ocean Engineering.
Key Words: Under-Film Corrosion, Random Fields, Cellular Automaton, Verification & Validation process, Covariance Structure

1. Introduction
To understand the influences that under-film corrosion has on
metal selection is currently insufficient in providing adequate
performance predictions for new Corrosion Resistance Steels(CRSs), without long-term experiments. There a few methods thatcan model the incubation and extension of epoxy coating failures.
This has prompted the need for a 2D cellular automaton (CA) isbased on a modified version of IACS CSR-H’s 3-phases
probabilistic model.
As part of the CA development, the model needs to undergo a
Verification and Validation process (V&V) consisting of asensitivity analysis on cell size and optimization examination on
the number of internal runs for a simulation. From thesemodifications, a re-examination will be performed on the results ofa single scribe defect reported by Osawa et al [1]; and the double
scribe defect reported by Oda et al [2].

2. Verification and Validation process
The purpose of this V&V is to improve the accuracy of the CA
while maintaining (or improving) the simulations run-time. This
will be achieved through two methods: the evaluation of the results
acceptable error based on the number of runs within the simulation;and the influence that cell size has to the CA.
2.1 Optimization study of Internal Computational Runs PerSimulationinternal runs for the maximum error vs each cell position from
-25mm to 25mm. The graph also shows the duration for eachsimulation indicating that 30 internal runs provides the best
reduction of overall error within an acceptable processing time.
2.2 Cell Size Sensitivity Analysis
After a lengthy investigation, a correlation was proposed that theprobability of a pit occurring within the cell (Ppit) was influencing
the depth. As the area of a cell was reduced/increased by changing
the cell size, Ppit was unchanged resulting in a higher probabilityof a pit occurring in the same area. This would cause an earlier
initiation of a pit within a cell resulting in a deeper pit due to alonger duration. This relationship is shown in equation (2),
𝐶𝑜𝑟𝑔
𝐶𝑛𝑒𝑤

𝑃𝑝𝑖𝑡

= √
𝑃

𝑜𝑟𝑔

(2)

𝑝𝑖𝑡 𝑛𝑒𝑤

where, Corg is the original cell size, Cnew is the new cell size. Ppitorg
is the original probability of a pit occurring in the cell and Ppitnew
is the new probability of a pit occurring in the cell.
Results for both CRS and conventional steel where comparedwith the original 2.5mm CA for 5mm, 1.25mm and 0.625mm. Fig
2 shows the original simulations with 2.5mm cell size and theoriginal Ppit of 0.036 for both CRS and conventional steel. The
figure also includes the new simulations with a cell size of 1.25mm
and new Ppit of 0.009 for comparison. All simulations arepresented with a minimum and maximum error profile.

The confidence interval formula, equation (1), shown below isused to determine the balance between accuracy and the number ofruns for the CA.
𝑥𝑏𝑎𝑟 = ±𝑧

𝑠
√𝑛

(1)

where, xbar is the sample mean, z is the confidence level set to 99%(= 2.575), and s is the standard deviation and n is the size of thesample. This equation is used to measure the standard error of 10
simulations with changing number of internal runs (n) for each
averaged cell position value over 10 simulations (s) that forms theuniformed surface profiles (corroding valley).

Fig 2. The original CRS and convention steel simulations
compared against new simulations with changed Ppit (0.0009) andcell size (1.25mm).
Fig 2 confirms that changing Ppit will allow the CA to be scaled.
However, there is an increase in divergence when the cell size is
reduced further to 0.625mm from 2.5mm, which may be caused by
the cellular automaton neighbor rules, e.g. coating life and coating
degradation, which controls the corroding valley’s slope.

Fig 1. Standard error for each cell position for simulations with 10,30, and 100 internal runs against processing time.
Fig 1 shows the original 10-run simulation against 30 and 100

3.

Re-examination of Scenarios from Previous Work

From the V & V, a re-examination is now required of past work
to ensure that the simulations produced by the CA are the best
depiction of the experiments possible.
This will only require slight changes to a select number of input
parameters. This investigation will repeat the process in simulating
both the single and double scribed analysis reported by Osawa etal [1] and Oda et al [2].

The generation of Random Fields for a and A0 was repeatedwith a range of 20mm to provide covariance structure for theseparameters, with the results shown in Fig 4.

3.1 Re-examination of Single Scribe Defect [1]The V & V found that the results of the CA had maintained an
accurate correlation with the experiments blister width for a
1.25mm cell size simulation. However, the depth and inclinationof the corroding valley could be improved upon. This was achievedby repeating the simulation refinement process presented inreference [1] from step 3 by decreasing the percentage of the meanand standard deviation of a, which controls the depth of thecorroding valley. Followed by step 4 in altering the following
parameter As, which controls the slope of the corroding valley.

Fig4. Double scribe defect simulations with and without
covariance structure for a and A0.
The probability of occurrence was then re-examined with minor
alterations to the methodology in reference [2]. The method starts
with quantifying the relationship between the two corroding
valleys’ depths to a single entity (), which is achieved withequation [2],𝛼 =Fig3. 2D corrosion wastage profiles comparison of CRS andconventional steel experiments against simulation created by theCA.
These results have indicated that the accuracy of a CA’s
simulation can be improved upon by changing the scale of themodel’s geometry with only slight corrections required of key
parameters, a and As. With Conversional steel decreasing only a by
85% and a and As for CRS by 65%. With only As requiring attention
when simulating corrosion wastage of CRS as the shielding effecthas a larger influence to the mathematical model then whenmodeling conventional steel.
3.2 Re-examination of double scribe defect [2]The double scribe simulation begins by defining the input
parameters from a single scribe experiment, which are then used in
the double scribe simulations. This methodology is used to allow
more complex scenarios to be simulated once the coating andmetal’s physical characteristics have been defined from the simpler
scenario under the same environmental conditions.
The conditions changed for this re-examination of the single
scribe defect from reference [2] was the cell size (1.25mm) andonly a and A0, by 65%. However, the corroding valley for thisexperiment showed two slopes create by increased corrosion from
15mm the scribe. The previous examination took the average but
for this examination, the inner slope is selected. These parameters
were then used for the double scribe defect simulation, which
concurred with previous results showing symmetry.

|𝐿𝑚𝑎𝑥 |−|𝑅𝑚𝑎𝑥 ||𝑑𝑎𝑣𝑔 |

[2]

where, Lmax and Rmax are the maximum depth of the left half andright half of the averaged surface profile, and davg is the average of
Lmax and Rmax.
The statistical population data of is generated by the CA of
one internal run per 1000 simulations, which differed, to theoriginal ten runs per 300 simulations. Table 3 shows the percentage
of probability of occurrence for the CA against each experiment
including their average. The results presented have confirmed theneed for using random fields to provide certain parameters in the
CA with covariance structure to resolve the non-symmetry foundin corrosion.
Table 1. Comparison of the probability of occurrence for  withand without covariance structure.
Case
Experiment 1
Experiment 2
Experiment 3
Experiment Avg.

No Covariancestructure included
4.37%
0.47%
0.02%
0.62%

Covariancestructure included
21.63%
11.18%
5.09%
12.07%

4. Conclusions
The V&V process has meet its objectives of providing amethod to refine the CA’s accuracy against a scenario and allowing
the CA to be scalable. This has been confirmed from the reexamination process showing improvements in simulating bothscenarios.

5. References
1) Osawa, et al (2016) ”Development of Under-Film Corrosion
Simulation Method based on Cellular Automaton”,
Proc.OMAE2016,paper OMAE2016-54508.
2) Oda, et al (2018) "Study on Under-Film Corrosion Simulation
Considering Covariance Structure of the Random
Parameters",Proc.JSNAOE2018, paper JASNAOE2018-0323.